Pulsational Pair-instability Supernovae. I. Pre-collapse Evolution and Pulsational Mass Ejection

The structure and evolution of massive stars depend on stellar mass, metallicity, and rotation (e.g., Nomoto & Hashimoto 1988; Arnett 1996; Heger et al. 2000; Heger & Woosley 2002; Nomoto et al. 2013; Hirschi 2017; Limongi 2017; Meynet & Maeder 2017). In stars with the zero-age main-sequence (ZAMS) mass of M ∼ 10 – 80 M_⊙, hydrostatic burning progresses from light elements to heavy elements in the sequence of H, He, C, O, Ne, and Si burning, and finally an Fe core forms and gravitationally collapses to form a compact object, such as a neutron star.

For very massive stars with M ≳ 80 M_⊙ (i.e., an He core of mass greater than 35 M_⊙), although exact correspondence is strongly dependent on metallicity (El Eid et al. 1983; Heger & Woosley 2002; Hirschi 2017; Woosley 2017), the effects of electron–positron pair production ($\gamma \to {e}^{-}+{e}^{+}$) on stellar structure and evolution are important when the O-rich core is formed (Fowler & Hoyle 1964; Fraley 1968). Pair production causes the dynamically unstable contraction of the O-rich core, which ignites explosive O-burning. For 140 M_⊙ ≲ M ≲ 300 M_⊙, the nuclear energy released is large enough to disrupt the whole star, so that the star explodes as pair-instability supernovae (PISN) (Barkat et al. 1967; Bond et al. 1984; Baraffe et al. 2001). Above ∼300 M_⊙, the star collapses to form a black hole (BH) again. As a result, no BH can be formed with a mass between ∼50 M_⊙ and ∼150 M_⊙ (Heger & Woosley 2002). Such a mass gap may provide distinctive features in the mass spectrum of BHs through the detection of a merger event of binary BHs.

For stars with $M=80\mbox{--}140\,{M}_{\odot }$ (He cores of 35 to ∼65 M_⊙, Woosley 2017), explosive O-burning does not disrupt the whole star, but creates strong pulsations (Barkat et al. 1967; Rakavy & Shaviv 1967), which is called PPI. These stars undergo distinctive evolution compared to more massive or less massive stars. PPI is strong enough to induce massive mass ejection as in PISNe, while the star further evolves to form an Fe core that collapses into a compact object later as a core-collapse supernova (CCSN). PPI supernovae (PPISNe) of 80−140 M_⊙ stars are thus the hybrid of a PISN and a CCSN.

The exact ZAMS mass range of PPISNe depends on the mass loss by stellar wind, thus on metallicity, and also on rotation. For PPISN progenitors, the wind mass loss during H- and He-burning phases could contribute to the loss of almost a half of the initial progenitor mass (see, e.g., Table 2 of Woosley 2017). Such mass loss can suppress the formation of a massive He core. The exact mass of the He core as a function of metallicity and ZAMS mass remains less understood because the mass loss processes in massive stars are not well constrained (Renzo et al. 2017). Rotation provides additional support by the centripetal force, which allows PPISNe to be formed at an even higher progenitor mass (Glatzel et al. 1985; Chatzopoulos & Wheeler 2012).

In order to pin down the mass range of PPISNe, a mass survey of main-sequence star models is done in Heger & Woosley (2002) and Ohkubo et al. (2009) with focus on zero-metallicity stars. Large surveys at other metallicities can be also found in, e.g., Heger & Woosley (2010) and Sukhbold et al. (2016). A large array of stellar models covering also PPISNe with rotation has been further explored in Yoon et al. (2012).

The evolution of a PPISN is very dynamical in the late phase. During the pulsation, the dynamical timescale can be comparable with the nuclear timescale when the hydrostatic approximation is no longer a good one. Also, when the star drastically expands after the energetic nuclear burning triggered by the contraction, the subsequent shock breakout near the surface is obviously a dynamical phenomenon. This suggests that during this dynamical but short phase, hydrodynamics instead of hydrostatics is required in order to follow the evolution consistently. Unlike hydrodynamical studies of PISNe (Barkat et al. 1967; Heger & Woosley 2002; Umeda & Nomoto 2002; Scammapieco et al. 2005; Chatzopoulos et al. 2013; Chen et al. 2014), a systematic hydrodynamical study of PPI has been conducted only recently (e.g., Woosley & Heger 2015; Woosley 2017).

The optical aspect of PPI and PPISNe might explain some superluminous supernovae (SLSNe), such as SN2006gy (Woosley et al. 2007; Kasen et al. 2011; Chen et al. 2014). Recent modeling of SLSN PTF12dam (Tolstov et al. 2017) has required an explosion of a 40 M_⊙ star with a 20–40 M_⊙ circumstellar medium (CSM) with a total of 6 M_⊙ ⁵⁶Ni in the explosion. The shape, rise time, and fall rate of the light curves of such supernovae provide constraints on the composition, density, and velocity of the ejecta, which give insights into the modeling of PPISNe. This demonstrates the importance of tracking the mass loss history of a star prior to its collapse. The rich mass ejection can explain the dense CSM observed in some supernovae, such as SN2006jc (Foley et al. 2007). Supernova models in the PPISN mass range are further applied to explain some unusual objects, including SN2007bi (Moriya et al. 2010; Yoshida et al. 2014) and iPTF14hls (Woosley 2018).

There is also a possible connection to the well observed Eta Carinae, which has demonstrated significant mass loss of about 30 M_⊙ (Smith et al. 2007; Smith 2008).

Furthermore, recent detections of gravitational waves (GWs) emitted by the merging of BHs (Abbott et al. 2016a, 2016b), such as GW150914 and GW170729, imply the existence of BHs of masses ∼30–50 M_⊙. In order to study the mass spectrum of BHs in this mass range, the evolutionary path of this class of objects becomes necessary. Such observations have led to interest in the evolutionary origin of massive BHs, including PPI phenomena (e.g., Belczynski et al. 2017; Woosley 2017; Marchant et al. 2019). Our calculations will update the lower end of the "mass gap" of massive BHs (not near the boundary between neutron star and BH).

Owing to the above importance of PPI, we re-examine it by using the open-source stellar evolution code MESA (v8118; Paxton et al. 2011, 2013, 2015, 2017).

We use this version because the recent update of the code (Paxton et al. 2015) has included an implicit energy-conserving (Grott et al. 2005) hydrodynamical scheme as one of its evolution options.

We study a series of evolutions of stars from the ZAMS for masses ranging from 80 to 140 M_⊙ and various metallicities. This corresponds to He-core masses from ∼40 to 65 M_⊙. Then we calculate the evolution of such He stars to study the hydrodynamical behavior of PPI including mass ejection.

In Section 2 we describe the code for preparing the initial models and the details of the one-dimensional implicit hydrodynamics code for the pulsation phase.

In Section 3 we examine the evolutionary path of PPISNe in the H- and He-burning phases and the influence of metallicity on the final He-core and CO-core masses.

Then in Section 4, we first present the pre-pulsation evolution of our models, which includes He- and C-burning phases. We study the dynamics of the pulsation and its effects on the shock-induced mass loss. After that, we present evolution models of He cores with 40–64 M_⊙. We examine their properties from four aspects—the thermodynamic, mass loss, energetics, and chemical properties.

In Section 5 we examine the connections of our models to SLSN progenitors.

In Section 6 we compare the final stellar mass of our PPISN models with the recently measured BH masses detected by GW signals.

In Section 7 we conclude our results.

We present in Appendix A a comparison of our numerical models with those in the literature. Appendix B and Appendix C examine the effects of some physical inputs in the numerical modeling, including convective mixing and artificial viscosity.

To prepare the pre-collapse model, we use the open-source code Modules for Experiments in Stellar Astrophysics (MESA) (v8118; Paxton et al. 2011, 2013, 2015, 2017). It is a one-dimensional stellar evolution code. Recent updates of this code have also included packages for stellar pulsation analysis and implicit hydrodynamics extension with artificial viscosity. We modify the package ccsn to build a model of an He-core or a main-sequence star directly and then we switch to the hydrodynamics formalism according to the global dynamical timescale of the star.

To understand the behavior of pulsation and runaway burning in the O core, we use the one-dimensional implicit hydrodynamics option. This option appears in the third instrument paper (Paxton et al. 2015). The energy-conserving scheme, coupled with the implicit mass-conserving property of the Lagrangian formalism, allows us to trace the evolution of the star consistently.

We refer readers to the instrument paper (Paxton et al. 2015, Section 4) where the detailed implementation of this mass- and energy-conserving implicit hydrodynamics scheme is documented. Here we briefly outline the specific points that are relevant to our calculation here.

The realization of this scheme relies on the use of artificial viscosity as a substitute for the exact Riemann solver. To capture the shock, the artificial viscosity takes the form

$$\begin{eqnarray}&&{Q}_{i}=-{C}_{{\rm{a}}}{\rho }_{i}\displaystyle \frac{4\pi {r}_{i}^{6}}{{{dm}}_{i}}\left(\displaystyle \frac{{v}_{i+1}}{{r}_{i+1}}-\displaystyle \frac{{v}_{i}}{{r}_{i}}\right),\end{eqnarray} \tag{ 1 }$$

which has the same unit as the pressure term and it enters the system of equations via $P\to (P+Q)$. v_i and r_i are the velocity and radius of mass shell i defined at the cell boundary. dm_i is the mass of the fluid element. We choose C_a = 0.002–0.02. However, the value of C_a has to be chosen by experience. Too large a value may dissipate the propagation of outgoing waves too early, which artificially suppresses the mass loss. Too small a value for C_a may create numerical difficulties when the shock becomes too strong for the hydrodynamics to handle, especially near the surface. We study the effects of the choice of C_a in Appendix C.

We define the physical quantities by convention as follows. Density, temperature, isotope mass fractions, and specific internal and related thermodynamics quantities are defined at the cell centers. Position, velocity, acceleration, and gravity source terms are defined at the cell boundaries. We impose the innermost boundary conditions as r₀ = 0.

The typical timescale during the pulsation is comparable to the dynamical timescale. However, after the pulsation phase, it is the Kelvin–Helmholtz (KH) timescale that dominates the contraction. Even with the implicit nature of the dynamics code, simply using the hydrodynamics formalism to evolve the whole pulsation phase is computationally challenging because the Courant–Friedrich–Levy condition limits the maximum possible timescale, despite the virtue of consistency in our calculation. We set conditions for the code to switch back to the hydrostatic approximation. When the star expands sufficiently after bounce that the evolutionary timescale is dominated by the thermal timescale, we increase the maximum time step at every 100 steps. When the star can evolve continuously with the maximum time step (10⁵ times the Courant time step), we change to the hydrostatic approximation to evolve the star until another pulsation starts. If the star appears to be non-static during the 100-step buffer, the buffer is extended until the star is fully relaxed. The convective mixing is also switched on only in the hydrostatic mode. In general, we find that dynamical treatment is necessary when the central temperature of the star exceeds 10^9.3 K.

In the pulsation phase, once the expansion of the star reaches the surface, it develops into a high-velocity outburst due to the density gradient near the surface. The fluid elements can have a velocity larger than the escape velocity. The ejected mass is dynamically irrelevant to the core evolution. We remove those mass elements that satisfy this condition and have a density below 10⁻⁶ g cm⁻³. We set a mass loss rate according to the velocity with which the outermost shell leaves our system. To avoid removing a mass shell at an unphysical rate due to interpolation, the mass loss rate is capped to a maximum value.

The code uses the Helmholtz equation of state (Timmes & Arnett 1999), which contains an electron gas with arbitrary relativistic and degeneracy levels, ions in the form of a classical ideal gas, a photon gas with a Planck distribution, and electron–positron pairs. To model the nuclear reactions, we use the "approx21_plus_co56.net" network. This includes the α-chain network (⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ³²S, ³⁶Ar, ⁴⁰Ca, ⁴⁴Ti, ⁴⁸Cr, ⁵²Fe, and ⁵⁶Ni), ¹H, ³He, and ¹⁴N for the hydrogen-burning and CNO cycle, and ⁵⁶Fe and ⁵⁶Co to trace the decay chain of ⁵⁶Ni. ⁵⁶Cr is included to mimic the neutron-rich isotopes formed after electron capture in nuclear statistical equilibrium (NSE).

The MESA equation of state is a blend of the OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), PTEH (Pols et al. 1995), HELM (Timmes & Swesty 2000), and PC (Potekhin & Chabrier 2010) equations of state.

Radiative opacities are primarily from OPAL (Iglesias & Rogers 1993, 1996), with low-temperature data from Ferguson et al. (2005) and the high-temperature regime dominated by Compton scattering from Buchler & Yueh (1976). Electron conduction opacities are from Cassisi et al. (2007).

Nuclear reaction rates are from JINA REACLIB (Cyburt et al. 2010) plus additional tabulated weak reaction rates (Fuller et al. 1985; Oda et al. 1994; Langanke & Martínez-Pinedo 2000). Screening is included via the prescription of Salpeter (1954), Dewitt et al. (1973), Alastuey & Jancovici (1978), and Itoh et al. (1979). Thermal neutrino loss rates are from Itoh et al. (1996).

As indicated in Woosley (2017), convective mixing is important because it redistributes the fuel and ash in the remnant core. This affects the subsequent nuclear burning when the star contracts again. We choose the mixing length theory (MLT) (Böhm-Vitense 1958; see, e.g., Cox & Giuli 1968 for a realization) to model the convective process with the Schwarzschild criterion. The MLT approximation is used in the main-sequence phase and also when the star enters the expansion phase.

We have attempted to couple the convective mixing in the dynamical phase but it results in numerical instabilities. We notice that the convective timescale during the dynamical phase is longer than the dynamical timescale. Furthermore, the more massive the star is, the shorter the burning timescale and its contraction timescale due to PPI. The mixing process becomes more inefficient than in lower-mass stars. We are interested in the mass loss process of PPISN, so convection can be neglected to a good approximation (also see the corresponding Kippenhahn diagrams in Section 4). Therefore, it becomes numerically manageable yet physically consistent to ignore convective mixing in the dynamical phase.

However, we also notice that in the lower mass regime (e.g., the 40 M_⊙ case), the contraction timescale is in fact long enough that convective mixing becomes more important. In the Appendix we examine the importance of mixing to the pulsation history for lower-mass PPISNe. We recall that during pulsations not only is convective mixing suppressed, but convective energy transport is also suppressed. This plays a major role in the weak pulsations found by other authors, such as the ones shown in models He36 and He40 in Figure 3 of Woosley (2017). In these cases there is only a mild collapse, and the energy produced through nuclear burning can be transported by convection without an eruption.

In Figures 1–3, we plot the Kippenhahn diagram of stars with M = 80,  100, and 120 M_⊙ at Z = Z_⊙. The red, green, and blue lines correspond to the He-, C-, and O-core mass coordinates respectively. Grey shaded regions are the convective zones inside the star. All models are run until the core reaches a central temperature of 10^9.4 K.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At solar metallicity the high metal content in the initial composition has largely increased the opacity, which allows strong mass loss during H-burning and He-burning due to the star's intrinsic high luminosity. It has an extremely large mass loss rate such that half of the matter is lost in the helium-burning phase for M = 80 and 100 M_⊙, and in the hydrogen-burning phase for 120 M_⊙. The whole H envelope is lost during He-burning, which occurs about 10⁵ yr before collapse. The initial He-core mass can reach about half of the initial total mass, but it gradually decreases due to the later mass loss. Also, in all three models, after the removal of the H envelope or He-burning, the C core quickly forms with a mass similar to that of the He core. Models of 80 and 100 M_⊙ have a C-core mass ∼20 M_⊙, which remains unchanged after it has been formed. The 120 M_⊙ model has a somewhat smaller one due to the previous drastic mass loss.

The convective pattern of the star is consistent with that of a typical massive star. In the H-burning phase, the core is mostly convective while the surface is radiative. In the He-burning phase, the core remains convective while some of the H envelope becomes convective. But this feature disappears when the mass loss sheds the H envelope. Once the C core has formed at about 10⁴ yr before pulsation, the star begins to contract rapidly. The core becomes radiative. But together with the C-core and burning of the C envelope, layers of convective shells appear. They gradually propagate and reach the surface of the C core. When the core starts O-burning (10⁻¹ yr from the onset of first pulsation), the strong energy generation triggers large-scale convection whereby the whole C envelope becomes convective. The inner core of the O-rich region also becomes convective.

We plot Figures 4–6 similarly to the previous three figures but at Z = 0.1 Z_⊙. Different from the models at solar metallicity, the low metallicity implies low opacity in the matter, which thus lowers the mass loss during the H- and He-burning phase. There is a clear signature of a massive He core from 30 to 50 M_⊙. The He-core mass remains constant after it has formed. Near the occurrence of first pulsation, a massive CO core ∼10 M_⊙ is also formed. A generally larger O-rich core is formed at the end of the simulation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Due to the preservation of the H envelope after He-burning, the star consists of a rich structure of convective activity before its collapse. The convective core has a similar structure to the case of higher metallicity, due to the extended H envelope that remains after He-burning. There is also a second convective zone, which gradually moves inward to the stellar core from its initial 60 M_⊙ to ∼40 M_⊙. Thereafter the He core is fully convective during He-burning and the convective zone extends into the H envelope. The surface is also convective. During C-burning, the convective layer propagates outward from the core to the surface of the C core. The outer layers of He and C envelopes are convective. During contraction before the onset of O-burning, the core returns to being mainly radiative. Similar to the case of high metallicity, the core becomes convective near the onset of pulsation.

In the left panel of Figure 7 we plot the H-R diagram for the main-sequence star models with M = 80, 100, 120, and 140 M_⊙ included. All models are fixed at Z = 0.002 Z_⊙. For numerical stability we do not include mass loss for the 140 M_⊙ model. In the pre-pulsation evolution, the models follow the typical H-R diagram of a main-sequence star. H-burning occurs after the star has contracted. After H is exhausted, the star develops into a red giant with He-burning, which greatly increases its luminosity. Depending on the mass loss, the effective temperature can reduce considerably. Also, the typical luminosity increases with mass.

Zoom In Zoom Out Reset image size

In the right panel we plot the evolution of the central temperature against the central density for the same set of models. We also draw the pair-instability zone (defined by the adiabatic index Γ < 4/3). There is no intersection among models, showing that the thermodynamic properties of the core before pulsation depend only on its mass. The contraction of the core is mostly adiabatic (with a slope −3) in the diagram.

To study the effects of metallicity and rotation, we run pre-pulsation stellar evolution models for different metallicities from Z = 10⁻³ Z_⊙ up to Z = 1 Z_⊙ for a non-rotating main-sequence star. In Table 1 we tabulate the pre-pulsation configurations of the main-sequence stars for their He- and CO-core masses when the core exhausts all H and He respectively. The CO-core masses are written in brackets.

Table 1.  The Pre-pulsation He-core Mass at the Exhaustion of H in the Core

Mass (M⊙)	Z = 10−3 Z⊙a	$Z={10}^{-2}\,{Z}_{\odot }$ a	$Z=0.1\,{Z}_{\odot }$	$Z=0.5\,{Z}_{\odot }$	$Z=0.75{Z}_{\odot }$	$Z=1\,{Z}_{\odot }$
80	34.05	37.40 (27.20)	33.80 (23.93)	30.10 (23.96)	23.60 (21.09)	22.70 (18.66)
100	44.51	49.44 (37.69)	47.16 (34.22)	33.00 (30.65)	31.70 (28.50)	30.30 (24.85)
120	54.87	64.71 (48.40)	59.95 (43.48)	57.10 (41.20)	37.40 (31.73)	15.50 (12.02)
140	65.87	nil	70.85 (56.67)	60.78 (50.36)	20.80 (16.90)	12.80 (9.60)
160	76.50	83.31 (82.40)	89.99 (89.12)	52.93 (46.46)	15.00 (11.63)	11.99 (8.90)

Notes. The numbers in brackets are the CO-core mass at the exhaustion of He in the core. All masses are in units of the solar mass.

^aThe models assume no mass loss.

Download table as:  ASCIITypeset image

He-core mass grows monotonically with M when Z < 0.625  Z_⊙. For star models with a higher metallicity, the mass loss rate, which is proportional to the metallicity, makes the He-core mass drop at the high-mass end. This transition starts at a lower mass for models with a higher metallicity. Notice that the change and the transition mass are not linearly proportional to Z due to the nonlinear dependence on mass loss rate. Also, the mass loss affects the gravity, which changes the equilibrium structure of the star even in the H-burning phase.

In Figure 8 we also plot He-core mass against progenitor mass for different metallicities. On the one hand, at low mass, the He-core mass is not so sensitive to metallicity, and it approaches its asymptotic value when Z ≤ 10⁻²  Z_⊙. On the other hand, at high mass, the He-core mass is very sensitive to metallicity, and from Z = 0.625 Z_⊙ to Z = Z_⊙ it can drop by 90% for the star model with M = 160 M_⊙, to about 15 M_⊙. At such low mass, the He core already leaves the regime of pulsation pair-instability and evolves as a normal CCSN. Furthermore, the maximum He-core mass for models at solar metallicity only barely reaches the transition mass of 40 M_⊙.

Zoom In Zoom Out Reset image size

For models that completely cover the PPISN mass range (He star of mass 40–64 M_⊙), we require a stellar model with a metallicity of at most 0.1 Z_⊙. This shows that the PPISN is very sensitive to the progenitor metallicity, while stars with solar metallicity are less likely to form PPISNe owing to their mass loss.

Then we examine the CO-core mass. Before the CO-core mass can be defined, the massive CO core has already started its contraction, which increases its mass. The effect of metallicity is similar to the trend for He-core mass. In all models, the CO-core mass in general increases with progenitor mass at the lower-mass branch, but it drops at the high-mass end. The CO-core mass also shows a monotonic decrease with metallicity for the same progenitor mass. The effect of mass loss in models with near solar metallicity is more significant because the CO-core mass is less than 10% of the stellar mass, while in models with a lower metallicity it can be about one third of the progenitor mass.

In Figure 9 we also plot the CO-core mass against progenitor mass at different metallicities. The significance of the metallicity for the mass loss rate can be seen. As models increase from 0.5 Z_⊙ to 0.75 Z_⊙, the CO-core mass can drop by 75% at M = 160 M_⊙. The CO-core mass shows a clearer relation than the He-core mass. The mass scales clearly with metallicity, except at M = 160 M_⊙, where the model without mass loss (Z = 10⁻² Z_⊙) has a lower mass than its counterpart with Z = 0.1  Z_⊙.

Zoom In Zoom Out Reset image size

In this part we examine the overall evolution of the PPISN from the He core until the onset of Fe-core collapse.

In Figures 10 and 11 we plot the Kippenhahn diagrams of models He40A, He50A, He60A, and He62A. The colored zone is again the convective zone while the lines (solid, dotted, dashed, long-dashed, dotted–dashed) are the He-, C-, O-, Si-, and Fe-core mass coordinates. The x-axis is the time counting backward from the star's collapse. We define the core boundary to be the inner boundary of the mass fraction for that corresponding element to drop below 10⁻². Therefore, since we start from an He core, the He-rich surface, which is also the total mass of the star, stands for the He core. Notice that for the cases with strong mass ejection, the whole He-rich surface can be shed. Here the time is defined by the time remaining until the onset of final collapse.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For model He40A, after the strong pulsation, the star expands and the outer part of the star above 18 M_⊙ becomes convective. Also, the star has established its O and Si cores at m(r) ∼ 5 M_⊙ and its Fe core at m(r) ∼ 2 M_⊙. Radiative transfer remains the major means of energy transport in the core. Thin layers of convection shells can be found in most parts of the He and C envelope. At about 10⁻¹ yr, the Si and O core can reach as far as ∼30 M_⊙. This is because during the propagation of the acoustic wave near the surface, the density gradient accelerates the wave into a shock, which heats up the matter around it. As a result, in such He-rich material, it facilitates He-burning and yields product including C, O, Ne, and Si. However, together with the extended convection during the expansion–contraction phase, the outer O- and Si-rich zones disappear and the values correspond to the inner layers, which come from previous hydrostatic burning.

For model He50A, after the pulsation, the C, O, and Si cores are produced simultaneously. But the O and Si cores quickly retreat from 30 M_⊙ to 5 and 10 M_⊙ respectively. The early formation of O and Si cores occurs because when the shock reaches the surface, the shock heating is capable of producing O- and Si-rich material around that region, but no significant O and Si production takes place away from the shock-heated zone. However, after its production, the mixing and mass loss caused by pulsation quickly remove this material. As a result, the O- and Si-core mass coordinates return to the corresponding inner values, where the real O and Si cores are located. At 10⁻³ yr before the final collapse, the contraction of the star allows the central density to be high enough for burning until NSE. The Fe and neutron-rich cores form almost simultaneously at ∼2 M_⊙. Unlike model He40A, the inner core is no longer convective, except after pulsation. There is an extended period of time at ∼10⁻² yr before its final collapse in which the star continues to hold fragmented convective layers.

For model He60A, there is also no inner convective core after its pulsation. Again, the shock heating creates a temporary outer surface to the O and Si cores, but they return to their inner ones after mass ejection and mixing, to ∼5 and 15 M_⊙. Unlike the previous two models, the expanded star after pulsation does not reach any convective state before its second pulse or final collapse. An outward propagating convective structure can be seen from ∼1 yr before collapse. It moves from m(r) = 20 M_⊙ to 40 M_⊙. The convective zone is small so it does not contribute in bringing the fuel from outer layers to the actively burning layer. Similar to model He50A, the Fe and neutron-rich cores appear at ∼2 M_⊙ at 10⁻³ yr before collapse.

Model He62A is different from the previous three models because of its extensive mass loss after pulsation. After the first pulse, the star reaches a very extended period of ∼10⁴ yr in a fully convective state. Again, the convection washes away the external C and O envelope. The first pulse creates final C and O cores, which are located at 20 and 5 M_⊙ respectively. In the second pulsation, the Fe core is also produced and has a mass ∼3 M_⊙. During its contraction at 1 yr before its final collapse, the core reaches the third convective state. During contraction, the outer extended convective zone also moves outward from 20 to 40 M_⊙. The convective structure is again fragmented. A 2 M_⊙ Fe core is formed only near 10⁻³ yr before the final collapse.

We first present the results for the pre-collapse profile based on the He main-sequence star in Table 2. The pre-pulsation evolution uses the hydrostatic approximation and it is continued until the central temperature reaches 10^9.4 K, where the dynamical timescale begins to be comparable with the O-burning timescale. Below 10^9.3 K the star evolves in a quasi-static manner but not assuming hydrostatic equilibrium. From the table we can see that the initial He-core mass affects the pre-pulsation C and O core. We choose the He-core models with a mass from 40 to 64 M_⊙, which produce CO cores from 30.82 to 50.42 M_⊙, with the remaining unburnt He in the envelope.

In Figure 12 we present the initial model profile and its composition. We find that most models are very similar to each other, so for demonstration we picked ${M}_{\mathrm{He}}=60\,{M}_{\odot }$ as an example. The star consists of three parts: a slowly varying core that extends up to 50 M_⊙, an envelope of rapidly falling density, and a surface with rapidly falling temperature.

Zoom In Zoom Out Reset image size

In the right panel we plot the chemical abundance profile for the same model. The model contains a flat core of mostly ¹⁶O up to ∼50 M_⊙. Then it becomes C-rich and then He-rich until the surface of the star.

We first study the time evolution of the pulsation. To do so, we examine the second pulse of model He60A, which is a strong pulse (with mass ejection) of mass ≈ 10 M_⊙. We choose this particular pulse because it is strong enough to create global change in the profile so that we can understand the changes during the contraction phase (before the maximum central temperature in the pulsation) and expansion phase (after the minimum in the pulsation).

The core is mostly supported by the radiation pressure. With the catastrophe in pair production, the supporting pressure suddenly drops, where the core softens with a corresponding adiabatic index of the equation of state of γ < 4/3 in the core. However, unlike stars with a mass of 10–80 M_⊙, which have rich Fe cores at the moment of their collapse, in PPISNe and PISNe the core is mostly made of ¹⁶O when contraction starts. The softened core allows a very strong contraction and the ¹⁶O-rich core can reach the explosive temperature, which releases a large amount of energy, sufficient to disrupt the star. ²⁸Si and ⁵⁶Ni can be produced during the contraction, when the central temperature can reach above 10^9.5 K. As a result, the star stops its contraction and expands. The rapid expansion causes strong compression of matter on the surface, which efficiently causes ejection of high-velocity matter from the surface and dissipates the energy. After that, the core becomes bound again. The pulsation restarts after it has lost most of its previously produced energy by radiation and neutrinos. The whole process repeats until the ⁵⁶Fe core, formerly ⁵⁶Ni, exceeds the Chandrasekhar mass and collapses under its own gravity before the compression heating can reach the further outgoing ¹⁶O-rich envelope.

In Figure 13 we plot the temperature, density, and velocity evolution at selected times in the top left, top right, and middle left panels, respectively. We pick the profiles when the core temperature reaches 10^9.3, 10^9.4, and 10^9.5 K before its peak during the pulse for profiles 1–3, at its peak temperature for profile 4, and after the core has reached its peak temperature for profiles 5–7 for the same central temperature interval. We plot the chemical abundance profiles for isotopes ¹⁶O, ²⁸Si, and ⁵⁶Ni in the middle right, bottom left, and bottom right panels respectively.

Zoom In Zoom Out Reset image size

First we study the hydrodynamics quantities. For the temperature, in the contraction (expansion) phase the star shows a global heating (cooling) due to the compression (expansion) of matter, and no temperature discontinuity can be observed. This shows that the whole star contracts adiabatically, without producing explosive burning in the star. By comparing the temperature profiles at the same central temperature (profiles 1 and 7 for T_c = 10^9.3 K, profiles 2 and 6 for T_c = 10^9.4 K, and profiles 3 and 5 for T_c = 10^9.5 K), the net effect of nuclear burning can be extracted. The part outside q ∼ 0.3 has a higher temperature after the pulse. A similar comparison can be carried out for the density profile. The inner core within q ∼ 0.3 is unchanged after pulsation, while the density in the outer part increases. The velocity profiles show more features during the pulse. Before the star reaches its maximally compressed state, the velocity everywhere is much less than 10⁸ cm s⁻¹. At the peak of the pulse, the envelope has the highest infall velocity of ∼2 × 10⁸ cm s⁻¹. After that, in profile 5, the core starts the homologous expansion phase, with a sharp discontinuous velocity peak near the surface between the outward-going core matter and the infalling envelope. Beyond profile 6, the discontinuity reaches the surface and creates a shock breakout. The surface matter can freely escape from the star.

For the chemical composition, the effects of the pulse become clear. By the time of the second pulse, part of the ¹⁶O core has already been consumed and converted to ²⁸Si in the first pulse. During the compression, before the core reaches its maximum temperature, ¹⁶O is significantly consumed and forms ²⁸Si. When the core reaches the peak temperature, the O within q ≈ 0.06 is completely burnt, which is where intermediate-mass elements, such as ²⁸Si, are produced. However, Fe-peak elements, such as ⁵⁶Ni, are not yet produced. On the other hand, during the expansion phase, most O-burning has ceased, making the ¹⁶O and ²⁸Si unchanged after the central temperature reaches 10^9.4 K, while advanced burning still proceeds slowly to form Fe-peak elements.

Here we study some representative models of an He core with a mass from 40 to 62 M_⊙. They show very different pulsation histories by their numbers of pulses and their corresponding strengths. In Table 3 we tabulate the stellar mass and the mass of elements in the star after each pulse.

For model He40A, most of the pulses are weak; however, following each pulse, the mass of ¹⁶O is gradually consumed and produces ²⁸Si. In late pulses, where the core exceeds 10⁷ g cm⁻³, NSE elements are also produced. In the last pulse, the core is sufficiently compressed that an Fe core beyond 1.4 M_⊙ is produced, which is followed by later mass loss. Most of the ejected mass is He.

For model He45A, most of the pulses are weak. As the number of pulses increases, not only Si but also ⁵⁶Ni is produced. The last pulse, which is the strongest overall, produces about 0.56 M_⊙ Ni, while the generated heat creates a shock to eject about 6 M_⊙ matter before the final collapse.

For model He50A, there are fewer pulses and again only the last pulse is a strong one that can eject mass. Compared to previous models, more ¹⁶O is consumed in each pulse, which produces Si. At the final strong pulse, less Ni is produced, while the accompanying mass loss ejects the He in the envelope. It should be noted that its lower mass ejection than model He45A arises because the O in model He45A is burnt in a much compressed state. This creates a much stronger pulsation when the expansion approaches the surface, which increases the mass loss.

Model He55A has two strong pulses, in contrast to lower-mass models that have only one. Its pulses are qualitatively similar to those in model He50A.

Model He60A has no weak pulse. The contraction always causes a significant mass of O to be burnt to produce the thermal pressure to support the softened core against its contraction. Due to the strong mass ejection, at the end of the simulation the star almost runs out of He. However, one difference between this model and the others is that it has a much lower Ni mass after pulsation. Most of the Fe, which leads to the collapse, is created during the contraction toward collapse.

Model He62A has a similar pulse pattern to model He60A but is stronger. Each pulse can consume about 3−4 M_⊙ of O. Unlike previous models, model He62A has an abundant amount of O even during its contraction toward collapse, and O continues to be consumed before it collapses.

Model He64A, which is a PISN instead of PPISN, has only one pulse before its total destruction. Due to its much lower density when large-scale O-burning occurs, even when about a few solar masses of O is burnt during the pulse, the energy is sufficient to eject all mass when the pulse reaches the surface.

In Figures 14–16 we plot the central density and temperature against time for models He40A, He45A, He50A, He55A, He60A, and He62A. To show that the rapid contraction comes from the PPI, we show in each plot the zones where electron–positron pair creation, the dynamical instability induced by photodisintegration of matter in NSE at Y_e = 0.5 (Ohkubo et al. 2009), and the dynamical instability induced by general relativistic effects (Osaki 1966) apply. The arrows in the figures show where the pulses take place. Here we define weak and strong pulses as the pulsation of the star without and with mass ejection. The strength of the pulse is further defined by how much the core expands and cools down.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For model He40A, at the beginning the central density is the highest among all six models. It thus has weaker pulses because the core is more compact and degenerate. It has five weak pulses and one strong pulse (indicated by arrows in the figure) where each of the small pulses only leads to a small drop in the central density and temperature. Then the core quickly resumes its contraction again. Only at the final pulse, when the core begins to reach the Fe photodisintegration zone, does the softened core lead to a fast contraction and reach a central temperature T_c = 10^9.8 K. This triggers large-scale O-burning in the outer core, which leads to a drastic drop in the central density and temperature, showing that the star is expanding, until T_c reaches 10^9.2 K. Then the core resumes its contraction. Since most O in the core is burnt, the Si-burning cannot produce adequate energy to create further pulsations. The star directly collapses.

Model He45A shows fewer pulses than model He40A. It has three weak pulses and one strong pulse. The initial path is closer to the PC instability zone. The last pulse is triggered at T_c = 10^9.7 K and has the lowest T_c of 10^8.9 K when it is fully expanded.

Model He50A has only two weak pulses and one strong pulse. Its evolution shows less structure than the previous two models because of the earlier trigger of large-scale O-burning in the core. The core starts the big pulse when T_c = 10^9.6 K and its expansion causes T_c to reach 10^8.8 K at minimum. Before its collapse, there is a small wiggle along its trajectory. We notice that in this phase the core has a small pulsation when it core becomes degenerate.

Model He55A has one weak pulse and two strong pulses. The two strong pulses start when T_c reaches 10^9.5 and 10^9.7 K respectively, with a minimum temperature after relaxation at 10^8.8 and 10^8.6 K.

For model He60A, there is no weak pulse and two strong pulses, where the stellar core intersects with the PC instability zone during its expansion. The two pulses start when T_c reaches 10^9.4 and 10^9.6 K. The core finishes its expansion when it reaches 10^9.0 and 10^8.3 K.

For model He62A, the star model becomes very close to the PC instability where the core enters the zone for a short period of time during its expansion. It is similar to model He60A in that there are two strong pulses. The two peaks start at 10^9.5 and 10^9.7 K and both pulses end at a minimum temperature of 10^8.4 K, showing that the two pulses are of similar strength. After that, the core starts collapsing similarly to the other five models.

By comparing all six models, we can observe the following trend for the pulse structure as a function of progenitor mass. First, when the progenitor mass increases, the number of small pulses decreases while the number of big pulses increases. Second, the strength of the big pulses increases with progenitor mass, which leads to a lower central temperature and density during its expansion. Third, the path during its early pulses becomes closer to the PC instability as mass increases. Fourth, the second strong pulse is stronger than the first.

In Figures 17–19 we plot the energy evolution for models He40A, He45A, He50A, He55A, He60A, and He62A, including the total energy E_total, internal energy E_int, gravitational energy E_grav, and kinetic energy E_kin. The energy is scaled in order to make the comparison easier.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In all six models, the energy evolution does not depend on the stellar mass strongly, except for the energy scale. In all these models, the small pulses do not make observable changes in the energy except for very small wiggles. The contraction before a pulse leads to a denser and hotter core, where neutrino emission continuously draws energy from the system. At a big pulse, the total energy shows a rapid jump, which increases close to zero, then the ejection of mass quickly removes the generated energy, making the star bound again. Similar jumps in E_int and E_grav show that the core is strongly heated due to contraction heating and nuclear reactions. After that, the star reaches a quiescent state with a mild increase in total energy due to the ⁵⁶Ni decay, followed by a quick drop when it contracts again.

In Figures 20–22 we plot the luminosity evolution for the six models similar to Figure 17. During the pulse, the extra energy from nuclear reactions allows the luminosity to grow by 3–4 orders of magnitude. For a short period of time (≤10⁻² yr), the star becomes dim suddenly. After that the star resumes its original luminosity quickly and remains unchanged until the next pulse or final collapse. We recall that the luminosity during shock breakout and shortly thereafter cannot be trusted because it requires in general non-equilibrium radiative transfer for an accurate treatment.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The neutrino luminosity is more sensitive to the structure of the star. The neutrino luminosity can also jump by 3–10 orders of magnitude from its typical value in the hydrostatic phase to the maximally compressed state. After the star has relaxed, the neutrino luminosity drops drastically. Depending on the strength of the pulse, neutrino cooling can become unimportant in the quiescent phase.

During the pulsation, when the bounce leads to explosive burning in the core and inner envelope, sufficient energy is produced to create an outgoing shock, where the outermost matter can gain sufficient energy to be ejected from the star. The ejected matter later cools down and becomes circumstellar matter (CSM). The existence of such H-free CSM is necessary in the circumstellar interaction models for Type I superluminous supernovae (SLSNe-I) (Sorokina et al. 2016; Tolstov et al. 2017). The chemical and hydrodynamic properties of the CSM thus become important, which influences the formation of the light curve of the explosion.

In Table 4, we tabulate the mass loss history of each model and its chemical composition. In Figures 23–26 we also plot the ejecta profiles of a representative pulsation taken from models He40A, He50A, He60A, and He62A respectively. Three patterns can be observed in the mass ejection. We choose these models because these examples characterize the typical ejecta features of strong pulses in the lower and higher mass regimes. We take the numerical values when the mass shells are ejected during the pulsation because that is the last moment at which the code keeps track of their evolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The first group is the strong pulse in the lower-mass branch. In models He40A and He45A the last pulse is the one that ejects mass. It displays wiggles in its density profiles, showing that the thermal expansion creates the first wave of mass ejection, followed by the shock as the velocity discontinuity approaches the surface, which creates the second wave of mass ejection. In both cases, only the He layer is affected, but as the He layer becomes thin, matter near the CO layer is ejected.

The second group is the weaker pulse of the more massive branch. In models He50A, He55A, He60A, and He62A the first strong pulse occurs after the core starts to consume O collectively. Since it burns much less O than other strong pulses, the ejection comes from the rapid expansion of the star, which includes matter in the He envelope.

The third group is the strong pulse of the more massive branch. In models He55A and He60A, the second pulse is stronger so the ejecta density gradually decreases. A continuous ejection of mass in terms of a smooth density profile is found. The mass ejection is sufficiently deep that, at the end of pulsation, traces of ¹²C and ¹⁶O can be found. We remark that the inclusion of massive elements (compared with H and He) will be important for future light-curve modeling because they contribute as the main source of opacity.

One of the pulsations needs to be discussed separately because of its very massive mass ejection, which involves a unique chemical composition in its ejecta. In model He62A (right plot), the second pulse becomes strong enough that, besides its decreasing density profile, the later ejected material contains a significant amount of heavier elements including C, O, Mg, and Si, showing that the He envelope is completely exhausted before the star is sufficiently relaxed.

In Figures 27–30 we plot the isotope profiles at different moments of models He40A, He50A, and He62A. We selected moments before and after each strong pulse to extract the nuclear burning history. The models are chosen to demonstrate how different strengths of pulsation and the convective mixing between pulses can create distinctive isotope abundances in the star. By comparing the isotope distribution, we can understand which part of burning contributes to the evolution of pulsation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In all models, it can be seen that the star is simply a pure O core with a minute amount of Si in the core or C in the envelope, covered by a pure He surface. However, their changes can be very different depending on the progenitor mass.

In model He40A, after the strong pulse, a range of elements are produced due to its previous weak pulses, which continue to burn matter in the core; these include ⁵²Fe, ⁵⁴Fe, ⁵⁶Fe, and also ⁵⁶Ni. There is a clear structure for each layer, which follows the order of ⁵⁶Fe, ⁵⁴Fe, ⁵⁶Ni, ⁴⁰Ca, ¹⁶O, and then ⁴He. After that, the core relaxes and becomes quiescent until it completely loses its thermal energy produced during the pulse, while at the same time convection redistributes the matter for a flat composition profile. Most convection occurs at q > 0.1, where the convective shells of different sizes make a staircase-like structure.

Models He45A and He50A share a similar nuclear reaction pattern. We choose model He50A as an example. The strong pulse provides the required temperature and density to make Ni in the center and Si in the outer zone. The Si-rich zone extends to q ≈ 0.2. During the quiescent phase, convection mixes the material not only in the envelope but also in the core, which is seen in the stepwise distribution of ⁵²Fe and ⁵⁴Fe.

Models He55A, He60A, and He62A also share a similar abundance pattern. We choose the evolution of He62A as an example. The first pulse makes the original O core into mostly Si and some Ca. Again, the convective mixing during the quiescent state redistributes the matter near the surface (q > 0.25). In the second strong pulse, the nuclear reaction is very similar to the late pulses of models He40A and He45A. Ni forms in the innermost part, with a small amount of Fe isotopes such as ⁵²Fe and ⁵⁴Fe. Then Si and Ca form the middle layer and finally the He envelope appears. During the quiescent phase, the convection occurs in a deeper layer that is absent in lower-mass models.

In this section we compare our results with some representative PPISN models in the literature.

In Yoshida et al. (2016), the PPISN models of mass from about 54 M_⊙ to 60 M_⊙ (corresponding to a progenitor mass from 140 M_⊙ to 250 M_⊙ with zero metallicity) are computed. In that work, the calculation is separated into two parts. During the quiescent and pre-pulsation phases, the hydrostatic stellar evolution code is used. During the pulsation phase, the star model is transferred to the dynamical code PPM, which follows the expansion of the star until the mass ejection has ended (∼10⁴ s). Then they map the results to the stellar evolution code again until the next pulsation.

Their 140 M_⊙ and 250 M_⊙ models have similar configurations to our models He55A and He60A. First, in their 140 M_⊙ model (250 M_⊙ model), they observe a total of six (three) pulses, which ejected 3.99 (7.87) M_⊙ of matter before collapse. Model He55A (He60A) exhibits three (two) pulses before collapse, which eject 6.78 (8.52) M_⊙ of matter. Our models show a smaller number of pulses, but give similar ejecta mass. This means that our models can capture the energetic pulse well, but not the smaller pulses.

Then we compare the ejection timescales. The 140 M_⊙ (250 M_⊙) model shows all pulses within a period of 0.92 (1434) yr, while model He55A (He60A) shows all pulses within a period of 1341 (2806) yr. There is a huge difference in the pulsation period in our model He55A and their 140 M_⊙ model. We notice that the difference comes from the strengths of the pulses. In particular, our second pulse leads to a transition at about 100 yr while ejecting 1.45 M_⊙. The most similar event in their model is the fourth pulse, but with a transition of only 0.279 yr.

At last we compare the final core composition. The 140 M_⊙ (250 M_⊙) model has an Fe (CO) core mass of 2.57 (43.51) M_⊙, while in our model, we have 2.49 (38.60) M_⊙ for the Fe (CO) mass. This shows that, despite the difference in the mass ejection history, our models can still capture the major mass ejection events, which results in a similar mass ejection and core composition. However, there is a strong pulse in our He55A model that is not seen in their 140 M_⊙ model.

Next, we compare our models with the models from Woosley (2017). We have chosen the PPISN close to that work; in particular, our models He40A, He50A, He60A, and He62A can be compared directly with the He40, He50, He60, and He62 models. In Woosley (2017), the Kepler code, which consists of both hydrostatic and hydrodynamic components, is used to follow the whole evolution of PPISN.

First we compare the mass ejection history. In Woosley (2017), there are 9, 6, 3, and 7 pulses with a total mass loss of 0.97, 6.31, 12.02, and 27.82 M_⊙ for models He40, He50, He60, and He62 respectively. In our models, we have 6, 3, 2, and 2 pulses with a total mass loss of 2.22, 2.61, 9.52, and 12.85 M_⊙ for models He40A, He50A, He60Am and He62A respectively. Again, our code tends to produce fewer pulses and the pulses in general eject less matter. One of the differences is how shock is treated. For a shock-capturing scheme with a larger dissipation, the kinetic energy will be partly dissipated into thermal energy, such that the star is globally thermalized instead of ejecting matter through kinetic pulses. Another origin of the differences can be related to the nature of the instability of PPISNe. The trigger of the explosive O-burning comes from the pair-instability, which is very sensitive to the initial condition (e.g., how we evolve the stellar evolution model before the pulsation and between pulses) and numerical treatment (e.g., how convection and mass ejection are treated). For example, a stronger contraction can lead to more O-burning in the core, which gives much stronger pulsation and hence more mass loss. In fact, such a dependence can also be seen in other fields: for example, in the propagation of flame, since it is unstable toward hydrodynamic instability (Glazyrin et al. 2013). The burning history can be highly irregular in the unstable regime.

Next we compare the timescale of the pulsation. In this work, the whole pulsation until collapse lasts for 0.38, 61.3, 2806, and 6610 yr for the four models, while in Woosley (2017) the durations are 2.48 × 10⁻³, 0.38, 2695, and 6976 yr. This shows that our results agree with their work for massive He cores, but there are large differences when the He core becomes less massive. In that case, our final pulse is always strong enough to re-expand the star again before the final collapse, which significantly lengthens the pulsation period.

Then we compare the Fe-core mass. In Woosley (2017) the core has 2.92, 2.76, 1.85, and 3.19 M_⊙ Fe. In our models, we have 3.42, 1.73, 1.64, and 2.66 M_⊙. There is a decreasing trend from model He40A to He60A, which corresponds to the trend that the pair-instability occurs at a lower density when the mass increases. On the other hand, near the pair-instability regime, the pulsation becomes sufficiently vigorous to enhance the burning in NSE.

Finally we compare the explosion energy. We compare model He62A, which has the largest explosion energy. In our model, the star has its total energy increased in the second big pulse by 2.0 × 10⁵¹ erg, while the maximum kinetic energy achieved is 2.8 × 10⁵¹ erg. This is very similar to the result in Woosley (2017), where the pulse is observed to have a kinetic energy of 2.8 × 10⁵¹ erg.

One major difference we notice is in the pair-instability limit: our model He64A shows a higher explosion energy. Across the strongest pulse, there is a change in total energy by 1.6 × 10⁵² erg, where the maximum kinetic energy of the system is ∼1.7 × 10⁵² erg. In Woosley (2017) the kinetic energy is reported to be 4 × 10⁵¹ erg. We observe that the difference comes from the number of pulsations, where our model He64A has two big pulses but there is only one in their work. The first pulse has incinerated the ¹⁶O in the core while ejecting mass from the surface. This means that the star has to reach a more compact state before it can explode. As a result, the amount of energy produced in the exploding pulse is much larger.

Our results show a systematically smaller number of pulses with slightly lower ejecta mass. The pulsation periods qualitatively agree with each other except for models with a final strong pulse, which may significantly lengthen our pulsation period. Also, in our explosion models, the system tends to store the energy in terms of internal energy instead of kinetic energy; as a result, the star tends to expand globally, where the excess energy and momentum of the star are transferred mostly to the surface. This ejects the low-density matter and leaves a bound and hot massive remnant. Despite the differences in the pulsation, globally the nucleosynthesis agrees in the two models because most of the heavy elements are produced by the strong pulses, where our results are consistent with those in the literature.

In Woosley (2019), the He star models are further evolved with mass loss. The solar metallicity of the Fe group is assumed in the mass loss rate for He stars. We note that this is not consistent with the mass loss history from the main sequence, because the He core of a star of solar metallicity becomes too small to undergo PPI (see our Figure 8). The final BH mass is lower than in his previous work. We notice that the PPISN models in this work have stronger pulsations in general. A 50 (60) M_⊙ He core ejects ∼7 (57) M_⊙, which is much higher than ∼6 (12) M_⊙ mass ejection in his previous work. Our results are closer to his previous work. This might depend on the thickness of the He layer, which is determined by the mass loss history.

This is a recent work that uses the same MESA code (Version 11123) (Paxton et al. 2017) to evolve the evolutionary path of PPISN. Their work has a similar setting to this work. Here we briefly compare their results with ours.

They have computed an array of single-star models from 40 to 240 M_⊙ with semi-convection, Riemann solver using the HLLC solver, and the approx21 nuclear reaction network. They treat the mass loss of the star by considering the average escape velocity.

Our work agrees qualitatively with theirs. For a model with a lower-mass He core, some distinctive differences can be seen. For example, multiple pulsations are observed in their models. They observe a total of four pulses for the 54 M_⊙ stars (corresponding to 39.73 M_⊙ at He depletion). On the other hand, our model with a 40 M_⊙ He star gives a total of six pulses. They observe in total 0.63 M_⊙ mass ejection before collapse while ours is about 2.2 M_⊙. The duration in their model is shorter (∼7 × 10⁻⁴ yr) but ours is longer (∼0.02 yr) after the onset of pulsation. For a model with a higher-mass He core such as ∼87 M_⊙ (corresponding to 60.04 M_⊙ at He depletion), they show only two pulses, which is the same as for our model of a 60 M_⊙ He star. The duration also agrees in the two models (theirs shows a duration of 7.5 × 10³ yr while ours is shorter at 3 × 10³ yr. A total of 4.6 M_⊙ mass loss is found in their model (the pre-He depletion mass loss is excluded) while ours is a higher value of ∼8.5 M_⊙.

We notice that their models and our models do not completely agree with each other. We notice that there are some critical differences in the implementation of this work from their work. First, they consider the evolution of H-free stars, with a metallicity at 0.1 Z_⊙. The He-core mass is therefore a function of the progenitor mass instead of a direct model parameter as controlled in our models. Furthermore, they use the Riemann solver (HLLC) in the newer version instead of the artificial viscosity scheme. How the pulse transfers into shock at the near-surface area can be different.