Collapse of Rotating Massive Stars Leading to Black Hole Formation and Energetic Supernovae

Numerical-relativity simulations are performed with our latest axisymmetric neutrino-radiation viscous-hydrodynamics code. Details of the code are described in Fujibayashi et al. (2017, 2020c). In this code, Einstein's equation is solved in the original version of the Baumgarte-Shapiro-Shibata-Nakamura formalism (Shibata & Nakamura 1995; Baumgarte & Shapiro 1998) with a constraint propagation prescription to make the constraint violation to propagate outward (Hilditch et al. 2013). A dynamical gauge condition described in Fujibayashi et al. (2017) is employed. To impose the axisymmetry for the geometrical variables, the so-called cartoon method (Alcubierre et al. 2001; Shibata 2000) with the fourth-order Lagrange interpolation is implemented.

The neutrino radiation transfer equations are approximately solved using a leakage scheme together with the truncated moment formalism (Fujibayashi et al. 2017; see also Sekiguchi 2010). In this formulation, the neutrino field is split into two components: trapped and free-streaming neutrinos. The trapped neutrinos are assumed to be tightly coupled with the fluid and have the same local temperature and velocity as those of the fluid. This component is treated as a part of the fluid and contributes to the internal energy and pressure. It becomes the free-streaming component with the generation rate controlled by the local diffusion rate of neutrinos. The free-streaming neutrinos are assumed to obey radiation transfer equations, which are solved by a truncated moment formalism with the M1 closure relation (Thorne 1981; Shibata et al. 2011). Following our previous work (Fujibayashi et al. 2020a, 2020b, 2020c), we solve the equations for the frequency-integrated energy and momentum density for the three neutrino radiation fields (electron, electron anti-, and other neutrinos).

The viscous-hydrodynamics equations are solved using the formulation described in Shibata et al. (2017), in which the energy-momentum tensor is written as

$$\begin{eqnarray}&&{T}_{\mu \nu }=\rho {{hu}}_{\mu }{u}_{\nu }+{{Pg}}_{\mu \nu }+\nu \rho h{\tau }_{\mu \nu }^{0},\end{eqnarray} \tag{ 3 }$$

where ρ = m_u n_b is the rest-mass density, h = c² + ε + P/ρ is the specific enthalpy with the specific internal energy ε and pressure P, ν is the kinematic viscous coefficient, and g_{μ ν} is the spacetime metric. The viscous tensor ${\tau }_{\mu \nu }^{0}$ is a symmetric tensor that satisfies ${\tau }_{\mu \nu }^{0}{u}^{\mu }=0$ (Israel & Stewart 1979), and it is determined by the following equation:

$$\begin{eqnarray}&&{{ \mathcal L }}_{u}{\tau }_{\mu \nu }^{0}=-\zeta ({\tau }_{\mu \nu }^{0}-{\sigma }_{\mu \nu }),\end{eqnarray} \tag{ 4 }$$

where ${{ \mathcal L }}_{u}$ denotes the Lie derivative with respect to u^μ , and ζ is a coefficient, for which we set ζ⁻¹ = O(10 μs). Assuming the form of the shear tensor as

$$\begin{eqnarray}&&{\sigma }_{\mu \nu }={{\rm{\nabla }}}_{\mu }{u}_{\nu }+{{\rm{\nabla }}}_{\nu }{u}_{\mu }={{ \mathcal L }}_{u}({g}_{\mu \nu }+{u}_{\mu }{u}_{\nu }),\end{eqnarray} \tag{ 5 }$$

where ∇_μ is the covariant derivative with respect to g_{μ ν} , we obtain the evolution equation for ${\tau }_{\mu \nu }:={\tau }_{\mu \nu }^{0}-\zeta ({g}_{\mu \nu }+{u}_{\mu }{u}_{\nu })$ as

$$\begin{eqnarray}&&{{ \mathcal L }}_{u}{\tau }_{\mu \nu }=-\zeta {\tau }_{\mu \nu }^{0}.\end{eqnarray} \tag{ 6 }$$

We only need to solve the spatial part of τ_{μ ν} because of the presence of the condition ${\tau }_{\mu \nu }^{0}{u}^{\mu }=0$. The spatial part obeys the following evolution equation (in Cartesian coordinates):

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}(\rho {u}^{t}\sqrt{-g}{\tau }_{\mathrm{ij}})+{\partial }_{k}(\rho {u}^{k}\sqrt{-g}{\tau }_{\mathrm{ij}})\\ +\rho {u}^{t}\sqrt{-g}({\tau }_{{ik}}{\partial }_{j}{v}^{k}+{\tau }_{{jk}}{\partial }_{i}{v}^{k})=-\rho \sqrt{-g}\zeta {\tau }_{\mathrm{ij}}^{0}.\end{array}\end{eqnarray} \tag{ 7 }$$

The effective viscosity in the disk is believed to arise as a result of the turbulence induced by magnetohydrodynmical instabilities such as MRI (Balbus & Hawley 1991, 1998) and Kelvin–Helmholtz instability (e.g., Obergaulinger et al. 2010). Following Shakura & Sunyaev (1973), we define the viscous coefficient by

$$\begin{eqnarray}&&\nu ={c}_{s}{{\ell }}_{\mathrm{tur}},\end{eqnarray} \tag{ 8 }$$

where ℓ_tur is the mixing length scale (or the largest eddy size) in the turbulence. In the α disk model, ℓ_tur is written as α_vis H (Shakura & Sunyaev 1973). In this study, we assume that ℓ_tur is proportional to the size of the black hole as

$$\begin{eqnarray}&&{{\ell }}_{\mathrm{tur}}=0.03\times \displaystyle \frac{2{{GM}}_{\mathrm{BH}}}{{c}^{2}},\end{eqnarray} \tag{ 9 }$$

where the black hole mass, M_BH, continuously increases with time due to the matter accretion. Thus, we employ a time-varying form of ℓ_tur. By the above definition, the mixing length scale becomes ℓ_tur ≈ 0.9 km(M_BH/10M_⊙). Since the disk scale height H is larger than 2GM_BH/c² for most parts of the disk around the black hole, Equation (9) implies that we assume α_vis ≤ 0.03, i.e., a conservative value of α_vis. Even for such a conservative value of α_vis, we will find a significant effect in Section 3. This spatially constant mixing length scale, Equation (9), may lead to a smaller kinetic viscous coefficient than that employed in Just et al. (2022a; see also Just et al. 2022b).

We note that the viscosity is incorporated just prior to the formation of the disk around the black hole. Specifically, we first perform simulations without viscosity until the disks are formed. We then go back to the time slice prior to the disk formation and rerun the simulation with viscosity. We check that there are any significant differences of the density and angular velocity profiles for infalling matter between the results of the simulations with and without viscosity. ⁵

The grid structure is the same as in the two-dimensional simulations recently performed with the same code (Fujibayashi et al. 2023), in which the cylindrical coordinates (R,z) are employed. In the inner cylindrical region of R ≤ 100Δx₀ and z ≤ 100Δx₀, a uniform grid with the grid spacing of Δx₀ is prepared, while in the outer region, a nonuniform grid with an increase rate of the grid spacing of 1 + δ is prepared. The value of δ, grid number N for each axis, and location of the outer boundaries along each axis (denoted by L) are listed in Table 1. We assume plane symmetry with respect to the z = 0 plane (the equatorial plane).

One of the theoretically accepted central engines of long-duration gamma-ray bursts is the system of a spinning black hole with a surrounding accretion disk (e.g., Woosley 1993). For the formation of a massive accretion disk around the black hole, the progenitor star has to be rapidly rotating. The progenitor models of Aguilera-Dena et al. (2020) may be promising for such a scenario. We employ, from their work, two of the rotating stars with zero-age main-sequence (ZAMS) masses of 9M_⊙ and 20M_⊙ (hereafter AD09 and AD20, respectively). Because of nearly chemically homogeneous evolution, the pre-collapse stars have very massive cores that are very compact. The compactnesses at the radius with its enclosed mass 2.5M_⊙ (referred to as ${r}_{2.5{M}_{\odot }})$, ${\xi }_{2.5}:=2.5/({r}_{2.5{M}_{\odot }}/1000\,\mathrm{km})$, are 0.68 and 0.66 for AD09 and AD20, respectively. This suggests that they are likely to form a black hole with no neutrino-driven explosion if it is nonrotating (see, e.g., O'Connor & Ott 2011, but see Burrows et al. 2019.)

It should be noted that the angular momentum transport via the convection, circulation, and magnetohydrodynamical interactions is taken into account only in an approximate way in current stellar evolution studies. The original rotation profile is applied for models AD09x1 and AD20x1, while for model AD20x2, the original angular velocity is doubled to investigate a possible case in which the rotation of the star is even faster.

Figure 1 shows the specific angular momentum along the equatorial direction as a function of enclosed mass for the models employed in this paper. The specific angular momentum at the innermost stable circular orbit for the black hole with the enclosed mass and angular momentum, j_ISCO, is also plotted in the dashed curves. For each progenitor model, the enclosed mass at which the solid and dashed curves cross approximately gives the black hole mass at which the infalling matter starts forming an accretion disk around the black hole. This figure indicates that disks are likely to be formed when the black holes grow to 4.7, 4.9, and 7.8M_⊙ for models AD20x2, AD09x1, and AD20x1, respectively. In Section 3, we will confirm that this prediction is approximately correct.

Zoom In Zoom Out Reset image size

At the time of the disk formation, the mass infalling rates for model AD20x2 are higher than those for AD09x1 and AD20x1. This makes a difference in the evolution process of the disk and the onset timing of the outflow from the disk as found in Sections 3.3 and 3.2.

In several previous studies, the mass ejection and nucleosynthesis calculation were carried out by simply modeling collapsar remnants as a black hole–disk system (Siegel et al. 2019; Miller et al. 2020). The evolution process of the black hole–disk system as a result of stellar collapse may be different from that of an isolated disk around a black hole. To clarify the possible differences, we also perform a simulation for a system of a spinning black hole and a massive disk (model BHdisk) as in our previous studies (Fujibayashi et al. 2020a, 2020b). The initial condition is constructed using code from Shibata (2007). The masses of the black hole and disk are 10 and 3M_⊙, respectively, and the dimensionless spin of the black hole is ≈0.6. The inner and outer edges of the disk are chosen to be 4GM_BH/c² ≈ 59 km and 400GM_BH/c² ≈ 5900 km, respectively. The constant entropy per baryon of s = 7k_B is simply assumed.

In this paper, we employ a tabulated equation of state (EOS) referred to as DD2 (Banik et al. 2014). We extended the table down to low density (ρ ≈ 0.17 g cm⁻³) and low temperature (k_B T = 10⁻³ MeV; see Hayashi et al. 2022) for the procedure.

To evolve black holes with good precision, the radius of the apparent horizon has to be resolved well (for the dependence of the black hole mass and spin on the grid resolution; see, e.g., Fujibayashi et al. 2020a). Employing a stiff EOS like DD2, which predicts a large maximum mass of neutron star, is advantageous to numerically resolve black holes with good accuracy (for a given computational resource), because the black hole mass (i.e., the radius of the apparent horizon) at its formation is larger for a stiffer EOS. Thus, for the present study, the DD2 EOS would be a better choice to save computational resources.

We define the mass-infall rate to the central region by

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{fall}}={\int }_{r={r}_{\mathrm{in}}}\sqrt{-g}\rho {u}^{k}{{ds}}_{k},\end{eqnarray} \tag{ 10 }$$

where ${{ds}}_{i}={r}^{2}{\delta }_{{ir}}d\cos \theta d\phi$ is the surface element of a sphere, $g=\det ({g}_{\mu \nu })$, and r_in is chosen in the following manner. In the early phase, it is the largest value among the radii of the surface of a standing accretion shock formed as a result of the core bounce; after this shock disappears due to the matter infall from the outer region, which induces the collapse of a proto-neutron star to a black hole, the maximum radius of the apparent horizon r_AH is used for r_in; after the formation of the disk, we again choose the largest value among the radii of the surface of a standing accretion shock, which is formed along the disk surface. In the late phase, we also define the mass accretion rate to the black hole by

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{BH}}={\int }_{r={r}_{\mathrm{AH}}}\sqrt{-g}\rho {u}^{k}{{ds}}_{k}.\end{eqnarray} \tag{ 11 }$$

The definitions of the unbound matter and explosion energy are the same as those in Fujibayashi et al. (2021). We first define the specific binding energy and binding energy flux density of the matter as

$$\begin{eqnarray}&&{e}_{\mathrm{bind}}:=\,-\displaystyle \frac{{T}_{t}{}^{t}}{\rho {u}^{t}}-(1+{\varepsilon }_{\min }),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{f}_{\mathrm{bind}}^{i}:=\,-{T}_{t}{}^{i}-\rho {u}^{i}(1+{\varepsilon }_{\min }),\end{eqnarray} \tag{ 13 }$$

where ${T}_{t}^{t}$ and ${T}_{t}^{i}$ are the time–time and time–space components of the energy-momentum tensor, and ${\varepsilon }_{\min }$ is the minimum specific internal energy of the employed EOS table (≈ −0.0013c² for DD2). We then define the condition for unbound matter as e_bind > 0. We find however that the binding energy of the outer layer of the star in our computational domain is not accurately calculated due to a numerical error accumulated in long-term simulations. Therefore, we decided to calculate the diagnostic explosion energy as well as the mass of the unbound matter as

$$\begin{eqnarray}&&\begin{array}{l}{E}_{\exp }:={\displaystyle \int }_{{e}_{\mathrm{bind}}\gt 0,r\lt {r}_{\mathrm{ext}}}{e}_{\mathrm{bind}}\rho {u}^{t}\sqrt{-g}{d}^{3}x\\ +{\displaystyle \int }^{t}{\displaystyle \int }_{{e}_{\mathrm{bind}}\gt 0,r={r}_{\mathrm{ext}}}{f}_{\mathrm{bind}}^{k}\,{{ds}}_{k}\,{dt},\\ {M}_{\mathrm{ej}}:={\displaystyle \int }_{{e}_{\mathrm{bind}}\gt 0,r\lt {r}_{\mathrm{ext}}}\rho {u}^{t}\sqrt{-g}{d}^{3}x\end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&+{\int }^{t}{\int }_{{e}_{\mathrm{bind}}\gt 0,r={r}_{\mathrm{ext}}}\rho {u}_{\mathrm{bind}}^{k}\,{{ds}}_{k}\,{dt},\end{eqnarray} \tag{ 15 }$$

where ds_k is the surface element of the sphere with radius r_ext. These quantities are defined by the volume integral for the matter of the positive binding energy inside an extraction radius r_ext plus the time integral for the components of the positive binding energy flux at the radius.

A part of the stellar matter is located outside the extraction radius and even outside the computational domain. Its mass and binding energy can contribute to the ejecta mass and explosion energy. To estimate their contribution, we first obtain the enclosed mass at the extraction radius when the shock wave reaches the extraction radius. We then estimate the binding energy of the matter above the radius with the same enclosed mass in the pre-collapse profile. We compare the explosion energies with different extraction radii in Section 3.4.

We note that the matter in the progenitor stars is composed mainly of ⁴He, ¹²C, ¹⁶O, and ²⁰Ne. When ¹⁶O burns into ⁵⁶Ni, the rest-mass energy of 0.67 MeV per baryon should be released into the internal energy, which can be an important energy source of the explosion of the star (if the internal energy is not carried away by the neutrino emission). However, this effect is absent in our simulation because of the assumption of nuclear statistical equilibrium (NSE) in constructing the EOSs, for which the low-temperature matter in the computational domain is assumed to be composed mostly of ⁵⁶Ni. In addition, ⁵⁶Ni is photodissociated into lighter nuclei if the temperature of the matter increases to ≳7 GK, consuming more internal energy than in the dissociation of ¹²C or ¹⁶O. These temperature conditions are found when the shock formed near the disk around the black hole has propagated sufficiently far out within the star and, hence, an excessively large energy of 8.6 MeV nucleon⁻¹, instead of ≈8.0 MeV nucleon⁻¹, is consumed. Thus, the explosion energy in our present simulations is underestimated. This point will be discussed in Section 3.6.

We also note that the inclusion of ${\varepsilon }_{\min }$ in Equations (12) and (13) is important to estimate the ejecta mass and explosion energy in a physically correct way. As mentioned above, the low-temperature matter in the outer region of the star is assumed to be composed of iron-group nuclei, which have smaller rest masses per baryon than the atomic mass unit (m_u ; the mass of ¹²C divided by 12). Due to the low temperature, such matter has a low specific internal energy, which leads to ε < 0 for the matter. If we do not include ${\varepsilon }_{\min }$ in Equations (12) and (13), such matter with ε < 0 is not recognized as ejecta, even though it gains sufficient energy to be unbound after being swept by the shock wave.

To perform the nucleosynthesis calculation, we apply our post-process tracer-particle method for the results of our simulations. The method is the same as that in Fujibayashi et al. (2023). The tracer particles are distributed for 128 polar angles in the range of θ = [0: π/2] on the arc with the radius of r_ext = 2 × 10⁴ km. The particles are continuously set with a time interval of Δt_set := r_extΔθ/〈v^r 〉, where Δθ = (π/2)/128 and 〈v^r 〉 is the average radial velocity of the ejecta at the extraction radius. The mass of each particle is determined based on the mass flux at the extraction radius as ${\rm{\Delta }}m={{r}_{\mathrm{ext}}}^{2}{\rm{\Delta }}{\rm{\Omega }}\rho {u}^{r}\sqrt{-g}{\rm{\Delta }}{t}_{\mathrm{set}}$, where ΔΩ is the solid angle element.

For all of the collapse models, a proto-neutron star is initially formed after the core bounce together with the standing accretion shock formation. Subsequently, due to the matter accretion from the outer region, the proto-neutron star collapses to a black hole. The black hole formation time is t_pb = 0.9, 0.9, and 2.3 s for models AD09x1, AD20x1, and AD20x2, respectively (denoted by t_BH in Table 2). Here, t_pb denotes the time after the core bounce. For model AD20x2, the formation of the black hole is delayed due to the significant centrifugal-effect associated with the rapid rotation. These results illustrate that the formation process of the black hole (and subsequent evolution process of the black hole and disk) depends on the profiles of the density and specific angular momentum of the progenitor stars.

For of the order of seconds after the black hole formation, the infalling matter does not have enough angular momentum to form disks, and thus, it is simply swallowed by the black hole. A geometrically thin disk starts forming at t_pb ≈ 7.0, 10.5, and 5.0 s at which the black hole mass is M_BH/M_⊙ ≈ 5.0, 7.0, and 4.5 for models AD09x1, AD20x1, and AD20x2, respectively (see, e.g., panel (a) of Figure 2 for model AD20x1). The above black hole mass is consistent with those inferred from the distributions of density and specific angular momentum of the progenitor stars (see Section 2.2).

The disks at their formation are geometrically thin because of the lower pressure inside the disk than the ram pressure of the infalling matter (Sekiguchi & Shibata 2011). The disks then become geometrically thick when the pressure in the outer part of the disk and ram pressure of the infalling matter become comparable. The vertical expansion of the disk occurs at t_pb ≈ 13, 21, and 10 s for models AD09x1, AD20x1, and AD20x2, respectively (see panel (b) of Figure 2 and panels (a) and (d) of Figure 3). In the presence of the viscosity, the disk expansion occurs at a time slightly later than that in the corresponding nonviscous model. The reason for this is that the viscous angular momentum transport accelerates the accretion of the disk matter onto the black hole, and the increase of the pressure in the disk is delayed.

After the vertical expansion of the disk, a shock surface between the disk and infalling matter expands (see, e.g., panel (b) of Figure 2). The geometrical cross section of the shock surface becomes large, which enhances the dissipation of the kinetic energy of the infalling matter at the shock more efficiently (Sekiguchi & Shibata 2011). As a consequence, the neutrino luminosity increases during this phase. The neutrino luminosity depends on the rate of mass supply to the disk (see Figure 4). For a rapidly rotating model AD20x2, the neutrino luminosity is higher because the mass-infall rate is higher.

After the vertical disk expansion, the models with the original rotational profiles AD09x1 and AD20x1, and a rapidly rotating model AD20x2 show different evolution process in terms of mass-infall rate and neutrino cooling efficiency. We describe the evolution processes in the following subsections separately.

For models AD09x1 and AD20x1, the rate of mass supply to the disk after the vertical expansion of the disk, ${\dot{M}}_{\mathrm{fall}}-{\dot{M}}_{\mathrm{BH}}$, is small (≲0.1M_⊙ s⁻¹; see the top panel of Figure 4). As a result, the disk temperature cannot be high enough for the efficient cooling by neutrino emission. Thus, the viscous heating dominates over the neutrino cooling in all phases after the disk expansion (the middle panel of Figure 4); a neutrino-dominated accretion disk is not formed in these models. The dominance of the viscous heating leads to an early launch of the outflow (several hundreds of milliseconds) after the disk expansion for these models. The outflow is launched mainly toward the equatorial direction (see panel (c) of Figure 2 and (b) of Figure 3). The final snapshots (see panel (d) of Figure 2 and (c) of Figure 3) also show the deformed profile of shock surfaces and the outflow from the disk toward the equatorial direction.

Zoom In Zoom Out Reset image size

The top panel of Figure 5 shows the mass histogram as a function of the electron fraction of the ejecta. For these models, the value of Y_e is not low (at lowest 0.47; see the top panel of Figure 5) because of their low disk density and low neutrino cooling efficiency throughout the disk evolution. This leads to weak electron degeneracy of the disk matter and thus keeps higher values of Y_e. Some components have an electron fraction even higher than 0.5. This reflects the fact that the positron capture proceeds in a shorter timescale than the electron capture under the condition of low electron degeneracy and mildly high temperature k_B T ≳ 1 MeV. In this condition, the positron capture n + e⁺ → p + ν_e is energetically more preferred than the electron capture because the mass difference between a free neutron plus an electron and a free proton (m_n + m_e − m_p )c² is not negligible compared with the kinetic energy of electrons ∼k_B T (see Just et al. 2022b; Arcones et al. 2010, and Beloborodov 2003). The results of these models are similar to those found in Just et al. (2022a; although our prescription of the viscous hydrodynamics is different from that of their study). ⁶ Because of the high electron fraction of the ejecta, we do not expect an r-process nucleosynthesis (see Section 4.2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The middle panel of Figure 5 shows the mass histogram as a function of entropy per baryon, s/k_B. These models have larger values of entropy per baryon s/k_B ≈ 30–50 and shorter expansion timescales r/v^r ∼ 50 ms than a rapidly rotating model AD20x2 (see Section 3.3). This result can be explained by the following three reasons. First, for moderately rotating models, the typical radius of the disk is smaller at the onset of the outflow because of the shorter time of the onset of the outflow after the disk formation. This requires a larger heating efficiency for the disk matter to be unbound and for the outflow matter to have the higher velocity. Second, since the outflow sets in before the disk settles into a quasi-steady phase, the shear of the poloidal velocity field in addition to that of the Keplerian motion leads to more efficient viscous heating. Third, the entropy generated by the viscous heating is not efficiently lost by the neutrino emission because of the lower neutrino cooling efficiency. These properties of moderately rotating models produce ejecta with a higher entropy and a shorter timescale.

This rapidly rotating model has appreciable differences from models AD09x1 and AD20x1 in two aspects. Specific angular momentum is much larger, and as a result, the mass infalling rate after the disk formation is much higher. In particular, the latter difference significantly modifies the evolution process of the disk from the model AD20x1.

For this model, a high-mass-infall phase continues prior to the disk outflow for a long timescale of ∼10 s. The reason for this is that during such a phase, the ram pressure from the infalling matter is strong and the temperature in the shocked region is high enough for enhancing the neutrino cooling, and hence, the thermal pressure generated by the viscous heating cannot be high enough for launching the disk outflow.

As a consequence, model AD20x2 exhibits a long-term quasi-steady phase of the disk in which the viscous heating and neutrino cooling rates are comparable for t_pb ≈ 10–16 s. In this quasi-steady phase, the neutrino cooling efficiency, defined by the neutrino luminosity L_ν divided by ${\dot{M}}_{\mathrm{BH}}{c}^{2}$, is several percent to 10% (see the bottom panel of Figure 4), which is sufficiently high for neutrinos to carry away the energy generated by viscous heating. That is, a neutrino-dominated accretion flow (NDAF) is established in this phase.

After the onset of the outflow from the disk, an expanding shock with a slightly oblate shape is formed (see panel (e) of Figure 3). During the shock expansion, the matter infall onto the black hole and disk still continues in the polar and equatorial directions because the outflow from the disk is launched along the surface of the geometrically thick disk (z ≈ 0.5R). The outflow toward the polar and equatorial directions is suppressed for different reasons. Because the outflow is launched mainly from the surface of the inner side of the disk, the outflow is prohibited in the equatorial direction due to the presence of the dense outer disk (bound matter). Near the polar axis, the infalling matter that passed through the expanding shock surface converges toward the polar region. As a result, the ram pressure near the polar region is enhanced and becomes larger than that in the other direction. This prevents the outflow toward the polar direction. The final snapshot for this model also shows that the outflow from the disk expands in the diagonal direction (see panels (f) and (i) of Figure 3; the matter with higher entropy, s/k_B ≳ 30, is the outflow component launched from the disk).

Although a neutron-rich region with Y_e < 0.2 is present in the inner mid-plane region of the disk reflecting the high density there (see panel (e) of Figure 3), the lowest value of Y_e of the ejecta is 0.40 for this model. The reason for this high value is that the outflow is launched from the disk surface region, in which the electron degeneracy is not as high as that in the disk mid-plane, and the value of Y_e is close to 0.5. The component with Y_e > 0.5 is present because of the same reason as for models AD09 and AD20x1 (see Section 3.2). From the ejecta with Y_e ≳ 0.5 and with T ≳ 5 GK, a substantial amount of ⁵⁶Ni can be synthesized. We discuss this topic in Section 3.6.

The typical value of s/k_B is found to be 10–20 for this rapidly rotating model, which is lower than those for models AD09x1 and AD20x1 (see the middle panel of Figure 5). This is because the outflow is developed well (∼10 s) after the formation of quasi-steady disks in a milder manner for the rapidly rotating model than for the moderately rotating models. This trend is also illustrated in the mass histogram of the expansion timescale (see the bottom panel of Figure 5). The expansion timescale for the rapidly rotating model is typically ≳0.1 s, which is longer than for the moderately rotating models. The entropy and expansion-timescale distributions for the rapidly rotating model are similar to those for model BHdisk, indicating that the formation of a dense disk in this model is a key to determining the properties of the ejecta.

The disk mass is ≈0.2M_⊙ for model AD20x2, and the cooling efficiency is low at the termination of the simulation. The electron fraction for the bulk of the disk matter is frozen to be Y_e = 0.4–0.5, because of a long weak interaction (electron/positron capture) timescale. Therefore, the electron fraction of the matter expected to be ejected in a longer timescale is likely to be 0.4–0.5 (see also Section 3.5).

As these results illustrate, the evolution of the disk, the timing for the onset of the explosion, and the electron fraction of the ejecta depend strongly on the distribution of the specific angular momentum of the progenitors and the resulting mass-infall rate on the disk. For the models with relatively low angular momentum (AD09x1 and AD20x1), the explosion occurs in a short timescale after the formation of a geometrically thick disk, at which the ram pressure of the infalling matter is low enough for launching an outflow. In this case, the electron fraction of the ejecta resulting from the explosion cannot be very low. By contrast, for the model with relatively high angular momentum (AD20x2), the explosion takes place at ∼10 s after the formation of the disk, and in this case, the electron fraction of the ejecta can be low with Y_e < 0.45, and thus, light trans-iron nuclei can be synthesized (see Section 4.2).

Figure 6 shows the mass of the ejecta and the diagnostic explosion energy as functions of time after the explosion, $t-{t}_{\exp }$. Here, the time of the explosion is defined as the time at which the explosion energy reaches 1 × 10⁵⁰ erg.

Zoom In Zoom Out Reset image size

For all of the viscous-hydrodynamics simulations, we find an explosion by the energy injection from the disk around the black hole. For models AD09x1 and AD20x1, for which the angular momentum of the progenitor stars is relatively low, the formation of the disk is delayed, and the onset of the explosion occurs in a late phase (see Table 2). As a consequence, only a small amount of mass remains outside the black hole at the explosion. The small energy budget at the formation of the disk leads to a slow increase in the explosion energy. For these models, the explosion energies amount to ≈(0.5–1) × 10⁵¹ erg at the termination of the simulations, which are comparable to or slightly smaller than the canonical value for core-collapse supernovae. The ejecta masses for these models are 0.1–0.3M_⊙ at the termination of the simulations, which are an order of magnitude smaller than that for canonical supernovae (but see below).

For the rapidly rotating model AD20x2, after the launch of the outflow, the explosion energy increases beyond 1 × 10⁵¹ erg with a timescale of ∼1 s. For this model, the explosion energy at the termination of the simulation reaches ≳3 × 10⁵¹ erg. At the termination of the simulations, disks with 0.2M_⊙ is present around the black holes for this model. In addition, matter infall still continues around the central region. Therefore, an energy budget to provide more explosion energy is still present. The result suggests that if the massive progenitors are rapidly rotating, a high-energy supernova-like explosion could occur from the disk outflow triggered by viscous heating. We discuss this point in Section 4.1.

For all of the collapse models, appreciable stellar matter is still present above the extraction radius. This can contribute additionally to the ejecta mass and explosion energy. The masses outside the extraction radius are larger than those of the ejecta, and thus, they may have an impact on the ejecta velocity. The terminal average velocity of the ejecta depends on the amount of the stellar matter that becomes ejecta in the later phase.

The mass and binding energy of the matter outside the extraction radii r_ext = 1 × 10⁹ cm and 2 × 10⁹ cm are listed in Table 2 (see the values in parentheses for columns of M_ej and E_exp; also see Figure 6). At the termination of the simulations, the absolute values of the binding energy are below the explosion energies for models AD09x1 and AD20x2. Thus, the explosion is likely successful for these models. For model AD20x1, on the other hand, the value is slightly above the explosion energy for r_ext = 2 × 10⁹ cm. Thus, it is not clear whether this model results in successful explosion. Note however that the bound matter outside the extraction radius at the termination of the simulation has a sufficiently large angular momentum to circularize around the black hole (see Figure 1). Thus, a possible increase may be expected in the explosion energy for the later phase if such matter falls into the central region to power an additional outflow. We plan to perform a longer-term simulation in our future work to address the possibility of the explosion for these models.

As already mentioned in Section 2.3, the composition of the computational domain is assumed to be that in NSE in our simulation. In reality, the outer region of the star is composed mainly of ¹²C, ¹⁶O, and ²⁰Ne, and hence, their nuclear burning in the stellar mantle swept by the shock wave can provide additional energy. In Section 3.6, we will investigate this more quantitatively.

For model BHdisk, which consists of a 3M_⊙ disk around a 10M_⊙ black hole with the dimensionless spin of 0.6, the disk matter accretes onto the black hole in a quasi-steady manner in the first 5 s. The outflow is then launched at t ≈ 5 s because the neutrino cooling rate has dropped far below the viscous heating rate (see the middle panel of Figure 4). The increased rate in the explosion energy for model BHdisk is much lower than for the collapsing star models (see bottom panel of Figure 6). The primarily reason for this is that, for the collapsing stars, a velocity shear is present not only in the accretion disk with the nearly Keplerian motion but also on the surface of the disk resulting from the infalling matter onto the disk, which is absent for model BHdisk. The strong shear on the disk surface significantly enhances the viscous heating rate, which results in the higher increase rate of the explosion energy.

The origin of the velocity shear associated with the infalling matter is different from the Keplerian motion of the disk. Here, for the latter we suppose that the MRI turbulence (Balbus & Hawley 1991) is the origin of the effective viscosity. For the surface region of the disk, on the other hand, we suppose that the shear region on the disk surface should induce the Kelvin–Helmholtz instability. In such regions, magnetic fields are supposed to be enhanced significantly leading to the development of turbulence and dissipating the kinetic energy of the infalling matter (e.g., Zhang et al. 2009; Obergaulinger et al. 2010; Rembiasz et al. 2016; Viganò et al. 2020). Thus, it is natural to consider that the effective viscosity on such a region is also high. ⁷ However, to clarify this process in the first-principle way, magnetohydrodynamics simulation is necessary. In future work, we plan to perform this to confirm that our assumption is indeed correct.

For model BHdisk, we stopped the simulation at t ≈ 21 s. For this model, mass ejection still continues with the ejection rate slightly higher than the mass accretion rate onto the black hole. This suggests that a large fraction of the disk matter will be eventually ejected from the system. At the termination of the simulation, the total ejecta mass is only ∼0.1M_⊙. However, the disk mass is still ≈2M_⊙. Extrapolating the mass ejection rate at the final time ∼10⁻² M_⊙ s⁻¹, we infer that the mass ejection continues for more than 100 s for this model.

The mass ejection for model BHdisk sets in after the neutrino cooling efficiency of the bulk of the disk drops (see the middle panel of Figure 4). The ejecta for model BHdisk has a low-Y_e component down to Y_e ≈ 0.4. The value is determined by the electron fraction in electron/positron capture equilibrium when the timescale of these reactions becomes comparable to that of the viscous expansion (see Fujibayashi et al. 2020c; Just et al. 2022b). Thus, the mass ejection mechanism for rapidly rotating model AD20x2 is qualitatively similar to that for model BHdisk.

Using the time evolution of the temperature and density along tracer particles, the post-process nucleosynthesis calculations are performed with the nuclear reaction network code rNET (Wanajo et al. 2018). The initial composition of the nucleosynthesis calculation depends on the thermal history of the tracer particles. If the maximum temperature along a particle is higher than 10 GK, we start the nucleosynthesis calculation at the time that the temperature decreases to T = 10 GK with the mass fraction of free protons and nucleons, Y_e and 1 − Y_e, respectively. For the tracer particles with the maximum temperature lower than 10 GK, we start the nucleosynthesis calculation at t = 0 of the simulation with the composition depending on the position of the particle in the progenitor star (mostly consisting of ¹⁶O and ²⁰Ne).

The resulting mass of ⁵⁶Ni is listed in Table 2. For models AD09x1 and AD20x1, the ⁵⁶Ni masses are smaller than 0.1M_⊙ (0.01 and 0.06M_⊙, respectively), because of their smaller ejecta masses. Thus, for moderately rotating progenitor models, the ⁵⁶Ni mass is likely to be comparable to or smaller than that for an ordinary supernova. These models predict the presence of moderately bright, but rapidly varying optical transients as found in Section 4.1.

It is found that for model AD20x2 the mass of ⁵⁶Ni amounts to 0.15M_⊙, reflecting the large mass of the ejecta that experiences high temperature ≳5 GK. Thus, the mass of ⁵⁶Ni found for this model could be high enough for explaining high-energy supernovae such as broad-lined type Ic (type Ic-BL) supernovae, considering that the ⁵⁶Ni mass inferred with the so-called "Arnett's rule" (Arnett 1982) is possibly overestimated (e.g., Meza & Anderson 2020, suggesting that the ⁵⁶Ni masses inferred from the radioactive tail luminosity for nearby two type Ic-BL supernovae, SN2009bb and SN2016coi, are 0.08 and 0.10M_⊙, respectively). These results suggest that massive and rapidly rotating stars leading to a black hole and massive disk are candidates for the progenitors of type Ic-BL supernovae.

We also note that the numerical simulation for this model underestimates the total ejecta mass because we stopped the simulations at a time when the ejecta mass and explosion energy are still increasing. Our results here indicate the lower bound for the ⁵⁶Ni mass.

It is important to note that, as can be found in Table 2, the produced ⁵⁶Ni mass fraction relative to the mass of ejecta exceeding 5 GK varies from 24% (AD20x2) to 43% (AD20x1) depending on the electron fraction, entropy, and expansion timescale of the outflowing matter for each model. The conditions of Y_e ≳ 0.5, low entropy, and slow expansion are favored for the efficient production of ⁵⁶Ni. Therefore, the mass of ejecta with ≥5 GK (as frequently used in the literature) only serves as a loose upper limit for the produced amount of ⁵⁶Ni.

Suppose that ⁵⁶Ni is synthesized from ¹⁶O, the rest-mass energy released into the internal energy due to the nuclear burning is 1.6 × 10⁴⁹ , 7.5 × 10⁴⁹ , and 1.9 × 10⁵⁰ erg for models AD09x1, AD20x1, and AD20x2, respectively. These contributions are several to tens percents of the explosion energy estimated in Section 2.3 and can have notable effects. Especially for model AD20x1, the explosion energy plus binding energy above the extraction radius r_ext = 2 × 10⁹ cm becomes negative, but the value is comparable to the energy generated by nuclear burning. Thus, to clarify whether such a marginal model explodes successfully, feedback of the nuclear reaction has to be taken into account in hydrodynamics simulations (see Bollig et al. 2021 and Navó et al. 2023 for a recent attempt).

Using the ejecta properties obtained in the present study, we analytically calculate bolometric luminosity models for photons following Arnett (1982). The thermalization efficiency for gamma-rays is estimated following Colgate et al. (1997) with the optical depth for nonthermal gamma-rays κ_γ = 0.03 cm² g⁻¹. For the optical depth of thermal photons, we simply set κ = 0.1 cm² g⁻¹. We note that the "Arnett" model tends to infer a larger ⁵⁶Ni mass by a factor of a few than that inferred by the "radioactive tail" luminosity of a supernova, the latter being less ambiguous (e.g., Meza & Anderson 2020; Afsariardchi et al. 2021; Rodríguez et al. 2022). This indicates that the luminosity predicted by the Arnett model for a given ⁵⁶Ni mass may be underestimated by a factor of a few (see Dessart et al. 2015, 2016 and Khatami & Kasen 2019). We also note that for the rapidly rotating model AD20x2, light trans-iron elements could be synthesized in the ejecta (see Section 4.2), and hence, the opacity for optical wavelengths may be higher than 0.1 cm² g⁻¹. For more quantitative study, we obviously need a radiation transfer simulation for photons taking into account a realistic opacity table.

Figure 7 shows the bolometric light curves for all of the models investigated in this paper. For deriving the solid curves of Figure 7, we take into account only the ejecta mass and explosion energy extracted at r_ext = 2 × 10⁹ cm. For models AD09x1 and AD20x1, the luminosity evolves rapidly. The rise times, defined by the time until the maximum luminosity is reached, are t_rise ≈ 3.3 and 4.4 days, respectively, for these models (see Table 2). Such fast transients may be discovered in the future high-cadence transient surveys. On the other hand, for model AD20x2, t_rise ≈ 10 days, which is consistent with that of type Ib/c supernovae (see, e.g., Taddia et al. 2018). This indicates that there could be a possible subclass of type Ib/c supernovae driven by the disk outflow in the black-hole-forming core collapses of rotating massive stars.

Zoom In Zoom Out Reset image size

To investigate the possible effects of the stellar matter outside the extraction radius, we also calculate the light curve for each model assuming that half of the mass and binding energy of the matter for r > r_ext = 2 × 10⁹ cm contribute to the ejecta mass and explosion energy (see the dashed curves in Figure 7). Because of the increased ejecta mass and decreased explosion energy, the peak luminosity and timescale for all of the light curves become smaller and longer, respectively.

For rapidly rotating model AD20x2, the timescale of the light curve is longer than that typically found for type Ib/c supernovae. Karamehmetoglu et al. (2022) reported such long-timescale type Ib/c supernovae recently. The long-timescale supernovae are reported as the explosions of massive (M_ZAMS ≳ 25M_⊙) stars with explosion energies comparable to typical type Ib/c supernovae, but they are inferred to have a larger amount (≳0.1M_⊙) of ⁵⁶Ni mass. These facts indicate that our rapidly rotating model can explain such a subclass of type Ib/c supernovae.

For the moderately rotating models AD09x1 and AD20x1, on the other hand, the timescales become ≈13 and 28 days, respectively, comparable to or longer than that of typical type Ib/c supernovae. The peak luminosity is, however, about 10 times dimmer than that of typical type Ib/c supernovae. This indicates that rotating massive stars exploded by outflows from black hole–disk systems may produce a variety of transients depending on the rotation profile of the progenitors and the presence of the stellar envelope, although the explosion mechanism is qualitatively universal.

As found in the comparison of the solid and dashed curves in Figure 7, the features of the bolometric light curves depend on the possible contribution of the matter outside the extraction radius. The quantitative prediction of the optical transients requires us to perform simulations for entire stars until the outer layer of the star is swept by the shock wave.

Another possible astrophysical transients can be powered by the interaction of the ejecta with a circumstellar medium that can result from the strong mass loss of their progenitor prior to the stellar core collapse. The progenitor models provided by Aguilera-Dena et al. (2018) are likely to be surrounded by a dense, massive (∼0.1–1M_⊙ within ∼10¹⁵ cm) circumstellar medium at the core collapse. Since the ejecta mass is of the order of 0.1–1M_⊙ for models AD09x1 and AD20x1, the ejecta will be significantly decelerated in the circumstellar medium, releasing a substantial fraction of its kinetic energy ≲10⁵¹ erg. The optical depth of the circumstellar medium is estimated as

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{CSM}}\sim \displaystyle \frac{3\kappa {M}_{\mathrm{CSM}}}{4\pi {{R}_{\mathrm{CSM}}}^{2}} & \approx & 170\left(\displaystyle \frac{\kappa }{0.35\,{\mathrm{cm}}^{2}/{\rm{g}}}\right)\\ & & \times \left(\displaystyle \frac{{M}_{\mathrm{CSM}}}{1{M}_{\odot }}\right){\left(\displaystyle \frac{{R}_{\mathrm{CSM}}}{{10}^{15}\,\mathrm{cm}}\right)}^{-2}.\end{array}\end{eqnarray} \tag{ 16 }$$

Here, the opacity for photons is assumed to be dominated by the Thomson scattering of fully ionized medium with Y_e = 0.875 (i.e., hydrogen and helium with mass fractions 0.75 and 0.25, respectively). Since the circumstellar medium is optically thick, the released energy diffuses out from the circumstellar medium with the diffusion time. The luminosity is then estimated as

$$\begin{eqnarray}\begin{array}{rcl}L & \sim & \displaystyle \frac{\epsilon {E}_{\exp }}{{t}_{\mathrm{diff}}}\cdot \displaystyle \frac{{R}_{\mathrm{CSM}}}{{v}_{\infty }{t}_{\mathrm{diff}}}\\ & \approx & 5.5\times {10}^{43}\mathrm{erg}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{\kappa }{0.35\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}}\right)}^{-1}\\ & & \times \left(\displaystyle \frac{\epsilon }{0.1}\right){\left(\displaystyle \frac{{v}_{\infty }}{1\times {10}^{9}\,\mathrm{cm}\,{{\rm{s}}}^{-1}}\right)}^{2}\left(\displaystyle \frac{{R}_{\mathrm{CSM}}}{{10}^{15}\,\mathrm{cm}}\right),\end{array}\end{eqnarray} \tag{ 17 }$$

where ${v}_{\infty }=\sqrt{2{E}_{\exp }/({M}_{\mathrm{ej}}+{M}_{\mathrm{CSM}})}$ is the terminal velocity of the ejecta plus the circumstellar medium, is the radiation efficiency, and

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{diff}} & \sim & \sqrt{\displaystyle \frac{\kappa ({M}_{\mathrm{ej}}+{M}_{\mathrm{CSM}})}{4\pi {v}_{\infty }c}}\\ & \approx & 26\,{\rm{d}}{\left(\displaystyle \frac{\kappa }{0.35\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}}\right)}^{1/2}\\ & & \times {\left(\displaystyle \frac{{M}_{\mathrm{ej}}+{M}_{\mathrm{CSM}}}{2{M}_{\odot }}\right)}^{3/4}{\left(\displaystyle \frac{{E}_{\exp }}{{10}^{51}\,\mathrm{erg}}\right)}^{-1/4}\end{array}\end{eqnarray} \tag{ 18 }$$

is the diffusion time of the expanding ejecta plus the circumstellar medium (see, e.g., Matsumoto & Metzger 2022 for a similar expression). The second factor of the first line in Equation (17) is the contribution of the adiabatic cooling. This can naturally lead to an optical transient like superluminous supernovae (see, e.g., Moriya et al. 2018 for a review). We note, however, that the properties of the transient depend not only on the mass and radius but also on the density profile of the circumstellar medium (see, e.g., Chevalier & Irwin 2011 and Suzuki et al. 2020 for the circumstellar medium, e.g., stationary wind).

Figure 8 shows isobaric mass fractions obtained by the nucleosynthesis calculations for our models. We find prominent peaks at α-nuclei (with A multiple of 4) and A = 56. The peaks at A = 12 and 16 reflect the initial composition of ¹²C and ¹⁶O, respectively, in progenitor stars. The bulk of these nuclei remains unprocessed owing to relatively low temperature achieved. On the other hand, the nuclei at peaks of A = 20–40 (corresponding to elements Ne, Mg, Si, S, Ar, and Ca) are synthesized from ¹²C and ¹⁶O in the ejecta that experience higher temperature but lower than that required for achieving NSE. The peak at A = 56 corresponds to ⁵⁶Ni, which is synthesized predominantly in NSE with Y_e ≳ 0.5. Interestingly, a certain amount of nuclei heavier than the iron group (A ≳ 60), up to A ≈ 90, is found to be synthesized for a rapidly rotating model AD20x2 (as well as BHdisk). It is known that in slightly neutron rich (Y_e ≳ 0.4) ejecta, such trans-iron nuclei are synthesized predominantly in quasi-nuclear statistical equilibrium (QSE; Meyer et al. 1998; Wanajo et al. 2018) under an α-rich condition. Since the neutron-richness is not very high, heavy r-process elements with A > 100 are not synthesized for any of the present models.

Zoom In Zoom Out Reset image size

The first peak nuclei of r-process, especially Zr and Y (synthesized in QSE here), are known to have opacities higher than those of iron-group elements (Kawaguchi et al. 2021). Therefore, if such elements are appreciably synthesized, the resulting optical transients may have longer timescales than those without such elements. In addition, the peak luminosity will be lower, and the spectrum could be redder. To quantify the light curve and spectrum, a radiation transfer simulation is needed in future work.

Figure 9 shows the masses (top) and dimensionless spins (bottom) of the black holes for our models as functions of post-bounce time. Here, the mass and dimensionless spin of the black holes are estimated from the equatorial and polar circumference radii of the apparent horizon (e.g., Shibata 2016).

Zoom In Zoom Out Reset image size

For models AD09x1 and AD20x1, in which we use the original angular momentum distribution of the stellar evolution simulations, the mass accretion and thus the spin-up of the black hole are suppressed because of the late-time disk formation and the quicker launch of the outflow from the disk. The dimensionless spins for these models are ∼0.4 and 0.7, respectively. For these models, the Blandford–Znajek mechanism (Blandford & Znajek 1977), which is one of the most promising mechanisms to power gamma-ray bursts (e.g., Gottlieb et al. 2022a), may provide only moderately large Poynting luminosity for launching intense electromagnetic waves because of its strong dependence on the black hole spin. Moreover, because of the low density of the disk for these models, strong magnetic fields may not be sustained in the vicinity of the black hole. This indicates that more rapidly rotating progenitors than those of AD09 and AD20 would be preferred for generating powerful relativistic jets.

For the rapidly rotating model AD20x2, the mass accretion from the disk onto the black hole continues for a long timescale (≳20 s) because of the presence of the long-term high-mass-accretion rate phase. This leads to a rapidly spinning black hole, the dimensionless spin of which is ∼0.8 at the termination of the simulations. During the long-term evolution of the black hole, the accretion disk is likely to be in a turbulent state associated with magnetohydrodynamical instabitlities such as the MRI (Balbus & Hawley 1991) in a realistic situation. By this, the magnetic-field strength should be enhanced, and the magnetic flux penetrating the black hole is increased as a result of the mass accretion (and magnetic flux accretion) onto the black hole. Such a highly spinning black hole penetrated by a strong magnetic field could be a promising central engine for relativistic jet via the Blandford–Znajek process (e.g., Christie et al. 2019; Hayashi et al. 2022 for a related topic).

It is known that at least a fraction of long-duration gamma-ray bursts are accompanied by type Ic-BL supernovae (e.g., Cano et al. 2017). As found in Section 3.6, rapidly rotating models synthesize a large amount of ⁵⁶Ni (>0.1M_⊙), which is consistent with the amount required to explain the light curves of type Ic-BL supernovae. In Section 3.4, it is also found that the rapidly rotating models show a large average velocity of the ejecta ≳2 × 10⁹ cm s⁻¹, which is necessary for the broad-line features for these supernovae. In addition, the large explosion energy ≳10⁵² erg observationally inferred from type Ic-BL supernovae is likely achieved if we consider longer evolution of the system than simulated in this study. Therefore, the rapidly rotating models in this study may reasonably represent the supernovae accompanying long-duration gamma-ray bursts.

We briefly discuss possible effects of a relativistic jet that may be launched in the polar direction by some mechanisms, accompanied with the formation of a rapidly rotating black hole. Here, we suppose that the jet is launched during the phase of a high mass accretion rate onto the black hole.

If a jet is powerful enough, it may partly prevent the infall of stellar matter (Tominaga 2009). If so, the ram pressure of the infalling matter decreases, and as a result, the outflow from the disk surface is launched earlier. The decrease of the matter infall also prevents the matter supply to the disk, and thus, the total energy for the outflow may be reduced.

Another possible effect of the relativistic jet is that the energy injection by the jet can be a source of ⁵⁶Ni production (see, e.g., Tominaga et al. 2007; Barnes et al. 2018; Leung et al. 2023). If the jet luminosity is high enough (≳10⁵³ erg s⁻¹; Tominaga et al. 2007), a significant amount of ⁵⁶Ni (≳0.1M_⊙) may be synthesized, and the optical transient may become more luminous. The energy injection by the jet may also modify the morphology of the ejecta, which can affect the features of the optical transient.