11institutetext:  1 Anton Pannekoek Institute for Astronomy and GRAPPA,
University of Amsterdam, NL-1090 GE Amsterdam, The Netherlands
2 Center for Computational Astrophysics, Flatiron Institute, New York, NY
10010, USA
3 School of Astronomy & Space Science, University of the Chinese Academy of
Sciences, Beijing 100012, China
4 National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100012, China
5 The observatories of the Carnegie institution for science, 813 Santa Barbara
St., Pasadena, CA 91101, USA
6 Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge,
MA 02138, USA
7 Geneva Observatory, University of Geneva, CH-1290 Sauverny, Switzerland
8 Department of Physics and Astronomy, Northwestern University, 2145 Sheridan
Road, Evanston, IL 60208, USA
9 Institute of Astrophysics, KU Leuven,Celestijnlaan 200D, 3001 Leuven,
Belgium
10 Steward Observatory, University of Arizona, 933 N. Cherry Ave., Tucson, AZ
85721, USA
11email<EMAIL_ADDRESS>
# Predictions for the hydrogen-free ejecta of pulsational pair-instability
supernovae

M. Renzo 1122 R. Farmer 1166 S. Justham 334411 Y. Götberg 55 S .E. de Mink
6611 E. Zapartas 77 P. Marchant 8899 N. Smith 1010

Present and upcoming time-domain astronomy efforts, in part driven by
gravitational-wave follow-up campaigns, will unveil a variety of rare
explosive transients in the sky. Here, we focus on pulsational pair-
instability evolution, which can result in signatures that are observable with
electromagnetic and gravitational waves. We simulated grids of bare helium
stars to characterize the resulting black hole (BH) masses together with the
ejecta composition, velocity, and thermal state. We find that the stars do not
react “elastically” to the thermonuclear ignition in the core: there is not a
one-to-one correspondence between pair-instability driven ignition and mass
ejections, which causes ambiguity as to what is an observable pulse. In
agreement with previous studies, we find that for initial helium core masses
of $37.5\,M_{\odot}\lesssim M_{\mathrm{He,init}}\lesssim 41\,M_{\odot}$,
corresponding to carbon-oxygen core masses $27.5\,M_{\odot}\lesssim
M_{\mathrm{CO}}\lesssim 30.1\,M_{\odot}$, the explosions are not strong enough
to affect the surface. With increasing initial helium core mass, they become
progressively stronger causing first large radial expansion
($41\,M_{\odot}\lesssim M_{\mathrm{He,init}}\lesssim 42\,M_{\odot}$,
corresponding to $30.1\,M_{\odot}\lesssim M_{\mathrm{CO}}\lesssim
30.8\,M_{\odot}$) and, finally, also mass ejection episodes (for
$M_{\mathrm{He,init}}\gtrsim 42\,M_{\odot}$, or $M_{\mathrm{CO}}\gtrsim
30.8\,M_{\odot}$). The lowest mass helium core to be fully disrupted in a
pair-instability supernova is $M_{\mathrm{He,init}}\simeq 80\,M_{\odot}$,
corresponding to $M_{\mathrm{CO}}\simeq 55\,M_{\odot}$. Models with
$M_{\mathrm{He,init}}\gtrsim 200\,M_{\odot}$ ($M_{\mathrm{CO}}\gtrsim
114\,M_{\odot}$) reach the photodisintegration regime, resulting in BHs with
masses of $M_{\mathrm{BH}}\gtrsim 125\,M_{\odot}$. Although this is currently
considered unlikely, if BHs from these models form via (weak) explosions, the
previously-ejected material might be hit by the blast wave and convert kinetic
energy into observable electromagnetic radiation. We characterize the
hydrogen-free circumstellar material from the pulsational pair-instability of
helium cores by simply assuming that the ejecta maintain a constant velocity
after ejection. We find that our models produce helium-rich ejecta with mass
of $10^{-3}\,M_{\odot}\lesssim M_{\mathrm{CSM}}\lesssim 40\,M_{\odot}$, the
larger values corresponding to the more massive progenitor stars. These ejecta
are typically launched at a few thousand ${\mathrm{km\ s^{-1}}}$ and reach
distances of $\sim$ $10^{12}-10^{15}\,\mathrm{cm}$ before the core-collapse of
the star. The delays between mass ejection events and the final collapse span
a wide and mass-dependent range (from subhour to $10^{4}$ years), and the
shells ejected can also collide with each other, powering supernova impostor
events before the final core-collapse. The range of properties we find
suggests a possible connection with (some) type Ibn supernovae.

###### Key Words.:

stars: massive, evolution, black holes, mass-loss — supernovae: general

##  1 Introduction

Figure 1: Evolution of a massive He core undergoing (pulsational) pair
instability evolution. Three final outcomes are possible: full disruption
without a compact remnant (4a.), formation of a BH because of the
photodisintegration instability (4c.), or episodic mass loss (4b.) and final
stabilization of the core, followed by a regular core-collapse event.

Massive stars can have diverse final fates depending on the structure of their
core at the end of their evolution. This diversity in terms of the kind of
collapse (“electron-capture”, iron-core collapse, pair-instability), whether
it triggers an explosion, albeit possibly weak, or not, and which remnant is
left behind is poorly understood as of yet. Typically, the core becomes
dynamically unstable and starts collapsing when no viable nuclear fuel is
left. However, for very massive radiation-pressure dominated stars, the core
becomes dynamically unstable while net energy generation by thermonuclear
reactions is still possible.

Stars that end their main sequence with a helium (He) core exceeding
$M_{\mathrm{He,init}}\gtrsim 80\,M_{\odot}$, which decreases to
$M_{\mathrm{He}}\simeq 60\,M_{\odot}$ accounting for the wind mass loss, are
predicted to end their evolution as pair-instability supernovae (PISN, Fowler
& Hoyle, 1964; Rakavy & Shaviv, 1967, see also Figure 1). They evolve in
hydrostatic equilibrium until they develop a carbon-oxygen (CO) core. Soon
after, the conversion of photons into electron-positron ($e^{\pm}$) pairs
(step 1 in Figure 1) causes a softening of the equation of state (EOS),
initiating the collapse of the star (step 2 in Figure 1). This increases the
inner temperature until explosive thermonuclear oxygen burning (step 3 in
Figure 1) reverts the collapse and fully disrupts the star (e.g., Barkat et
al., 1967; Fraley, 1968; Kasen et al., 2011; Yoshida et al., 2016; Woosley,
2017, 2019, step 4a in Figure 1). Such stars do not leave any compact remnant
at the end of their evolution.

For initial $M_{\mathrm{He,init}}\gtrsim 200\,M_{\odot}$, decreasing to
$M_{\mathrm{He}}\simeq 125\,M_{\odot}$ after wind mass loss, stars also
experience explosive thermonuclear oxygen burning but, owing to energy loss as
a result of the photo-disintegration of heavy nuclei, the explosion is not
energetic enough to reverse the collapse into an explosion and disrupt the
star, (e.g., Bond et al., 1984; Fryer et al., 2001; Heger et al., 2003, step
4c in Figure 1). In these cases, the final fate is core collapse (CC), forming
a massive black hole (BH). Therefore, if these stellar explosions do occur in
nature, a “PISN black hole mass gap” (also called “second mass gap”111The
“first gap” is the apparent lack of compact objects with masses between the
maximum neutron star mass, $\max\\{M_{\mathrm{NS}}\\}\simeq 2\,M_{\odot}$ and
the least massive BH known $\min\\{M_{\mathrm{BH}}\\}\simeq 5\,M_{\odot}$,
(e.g., Farr et al. 2011, but see also Wyrzykowski et al. 2016; Wyrzykowski &
Mandel 2019).) is expected between the most massive BH that can be formed
without encountering a PISN fate and the least massive BH formed because of
the photodisintegration instability.

The most massive BHs below the gap result from the evolution of He cores with
final masses just below $\sim$ $60\,M_{\odot}$ (e.g., Yoon et al., 2012;
Woosley, 2017; Farmer et al., 2019). In these stars, the explosive burning of
step 3 in Figure 1 releases less energy and thus is only able to eject a
fraction of the outer layers of the star. This produces a mass-loss pulse
(step 4b in Figure 1), without fully disrupting the star (Rakavy & Shaviv,
1967; Fraley, 1968; Woosley et al., 2002, 2007; Woosley, 2017, 2019). This
phenomenon is the lower-mass analog of a PISN, a pulsational pair-instability
(PPI). The star may undergo multiple such pulses until the combined effects of
pulsational mass loss, entropy loss to neutrinos (step 5 in Figure 1), and
fuel consumption stabilizes the core (Woosley, 2017; Marchant et al., 2019;
Farmer et al., 2019; Leung et al., 2019). Ultimately, this star is likely to
collapse to a BH, possibly with an associated supernova (SN), at step 7 in
Figure 1.

Given the impact on the distribution of BH masses (Belczynski et al., 2016;
Woosley, 2017; Marchant et al., 2019; Stevenson et al., 2019), the recent
direct detection of gravitational waves (Abbott et al., 2017, 2019) has
revived the interest in PPI evolution. Moreover, the follow-up of
gravitational wave merger events is driving large observational efforts in
time-domain astronomy, with new and upcoming facility such as the Zwicky
Transient Factory (Bellm, 2014), Large Synoptic Supernova Survey (LSST Science
Collaboration et al., 2009). When not following up gravitational wave events,
these instruments will perform surveys of various depth and cadence which will
soon unveil the variety of electromagnetic transients possible in the sky. The
James Webb Space Telescope will be able to probe the transients expected at
the death of the first stars in the Universe, increasing the chance of a
direct unambiguous detection of PISN or PPI (Whalen et al., 2013; Regos et
al., 2020). Another potential piece of indirect evidence for the occurrence of
PISNe is a peculiar distribution of isotopes in their yields, due to the
neutron-poor type of nucleosynthesis (so called odd-even effect, e.g., Woosley
et al., 2002). However, the detection of this effect in the surface mass
fractions of low metallicity stars remains debated (e.g., Aoki et al., 2014).

Therefore, it is important to characterize the observable characteristics of
PPI evolution, that is, address the question of what are the observable
signatures of a pulse. Of particular interest is the question of how much mass
do the pulses eject and at what velocity is it launched (Leung et al., 2019),
or in other words, what are the circumstellar material (CSM) structures that
this process can produce.

Previous studies from Chatzopoulos & Wheeler (2012a, b) investigated the fate
of (hydrogen-rich) stars with zero age main sequence (ZAMS) masses above
$40\,M_{\odot}$, with and without rotation during the pre-explosion
evolution222The effect of rotation on the explosion dynamics has been
investigated in Glatzel et al. (1985); Chatzopoulos et al. (2013)., and found
PPI evolution in the initial mass range $40-65\,M_{\odot}$ (for high rotation
rates) and $80-110\,M_{\odot}$ (without rotation). These ranges are also
sensitive to the details of the nuclear physics (e.g., Takahashi, 2018; Farmer
et al., 2019).

Woosley (2017, 2019), building up on previous work by Woosley et al. (2002,
2007), presented the first grids of stellar evolution calculations for a wide
mass range enclosing both PPI followed by a core collapse (PPI+CC) and PISN.
The light curves of the former are expected to show a series of brightening
events as the individual pulses collide with each other (Woosley, 2017), which
has been proposed to explain the extremely luminous light curve of SN2006gy
(Woosley et al., 2007). More recently, Arcavi et al. (2017); Woosley (2018)
also proposed PPI as a way to explain the peculiar photometric and
spectroscopic evolution of SN iPTF14hls, possibly coming from a merger
progenitor (Vigna-Gómez et al., 2019; Spera et al., 2019).

Although PISN need not be extremely luminous (Woosley, 2017), they are
routinely considered in the context of super-luminous supernovae (e.g., Gal-
Yam et al., 2009; Chatzopoulos et al., 2013). No unambiguous identification of
an astrophysical transient with a PISN is available as of yet, however
Kozyreva et al. (2018) proposed OGLE14-073 as a promising candidate. For
models evolving through PPI before their final collapse, recently claimed
observational candidates are PTF12dam, a fast rising type I super-luminous SNe
modeled by Tolstov et al. (2017) with a combination of CSM interaction and
radioactive decay; iPTF16eh, a type I super-luminous SN showing signs of a
shell of circumstellar material through the detection of a light echo (Lunnan
et al., 2018); and SN2016iet, for which a dense, hydrogen (H)- and He-free CSM
at $10^{15}\,\mathrm{cm}$ of the star can be invoked to explain the light
curve (Gomez et al., 2019). Other potential candidates are type Ibn SNe
showing relatively narrow He lines, such as SN2006jc, whose progenitor was
observed to experience an outburst two years before the final explosion
(Pastorello et al., 2007; Foley et al., 2007), and PS15dpn whose light curve
has also been modeled with a combination of CSM interaction and
${}^{56}\mathrm{Ni}$ decay by Wang & Li (2019).

Here we calculate the detailed evolution of massive He cores to characterize
jointly the effect that PPI evolution has on the final BH masses and on the
circumstellar material (CSM) structure. Our stellar evolution models can
provide input for the hydrodynamical evolution of the CSM, which can become
visible because of collisions between shells of CSM ejected at different times
(e.g., Woosley et al., 2007; Woosley, 2017), or possibly if the BH form with
an accompanied explosion. We also provide (_i_) a criterion to determine which
He cores encounter a _global_ instability, resulting in PPI-driven mass loss
and BH formation, and which He cores instead are fully disrupted in a PISN and
(_ii_) the bulk properties of the pulses and their distribution as a function
of mass.

In Section 2 we describe our calculations, before giving an overview of the
evolutionary outcome of our models in Section 3 and of the resulting BH masses
in Section 4. We focus on the PPI models in Section 5, where we describe three
physically motivated possible definitions of a “pulse”. While the basic ideas
on how the evolution of these models proceeds are well established from the
theoretical side, there is some ambiguity in the literature on what is called
a pulse. We discuss the CSM that our models can produce with a toy-model
assuming propagation of the ejecta at constant velocity in Section 6, and
provide input files for a more sophisticated modeling of the CSM at
https://doi.org/10.5281/zenodo.3406356. In Section 7 we discuss whether the
final core collapse after the PPI evolution would produce an associated SN
explosion, which would generate ejecta to interact with the previously ejected
stellar layers. We compare our results to a few observational transients that
have been interpreted as pulsational pair-instability events in Section 8. We
define a criterion to distinguish pulsational evolution from full disruption
without going through the hydrodynamic calculations in Section 9, before
highlighting the main limitations of this study. Section 10 summarizes our
main conclusions. Appendix A presents a resolution study of one of our models,
and Appendix B compares the evolution of a naked He core to a H rich star with
a similar He core mass.

##  2 Pulsational pair-instability evolution with MESA

We model the evolution of bare He cores because stars massive enough to
encounter the PPI are likely to have lost their H-rich envelop beforehand.
This could happen either because of the presence of a binary companion (e.g.,
Kippenhahn & Weigert, 1967), strong wind mass loss (e.g., Vink & de Koter,
2005), or because of rotational mixing preventing the formation of a core-
envelope structure (Maeder & Meynet, 2000; Yoon et al., 2006; de Mink et al.,
2009; Mandel & de Mink, 2016; Marchant et al., 2016).

Another way to form very massive stars which are expected to produce the most
massive (stellar mass) BHs is through runaway collisions in a dynamically
excited environment (e.g., van den Heuvel & Portegies Zwart, 2013), or binary
mergers (e.g., de Mink et al., 2014; Vigna-Gómez et al., 2019). Either might
result in the loss from the system of some H-rich material. Even if a binary
merges before the onset of pulsations and retains a significant amount of H
(e.g., Vigna-Gómez et al., 2019), the merger may well undergo
asteroseismologic (non-PPI) pulsations enhancing wind mass loss removing of
the remaining envelope (Moriya & Langer, 2015). Finally, any remaining H-rich
envelope is likely to be loosely bound and easily removed during the first PPI
pulse (Fraley, 1968; Leung et al., 2019, see also Appendix B).

We employ the open-source stellar evolution code MESA (release 11 701, Paxton
et al., 2011, 2013, 2015, 2018, 2019) to evolve a grid of He stars in the mass
range $35\,M_{\odot}\lesssim M_{\mathrm{He,init}}\lesssim 250\,M_{\odot}$.
Throughout this study, we define the He core mass as the total mass of our
models, and the CO core boundary as the outermost location where the mass
fraction of ${}^{4}\mathrm{He}$ drops below 0.01. We adopt an initial
metallicity of $Z=0.001$, and we also ran a limited sample of models with
$Z=0.00198$ (similar to the value quoted for SN2016iet, Gomez et al., 2019).
Both these values are below the upper limit for the occurrence of PISN of
$Z_{\odot}/3\simeq 0.006$ obtained from single star models (Langer et al.,
2007). Pair-instability evolution might even occur at higher metallicity
because of late stellar mergers in a binary (Vigna-Gómez et al., 2019), or if
magnetic fields funnel the mass lost to winds back to the star (Georgy et al.,
2017). We do not study here the impact of binarity, rotation, and magnetic
fields.

We include wind mass loss as in Marchant et al. (2019) and in the fiducial
model of Farmer et al. (2019), that is we use the rate from Hamann et al.
(1995); Hamann & Koesterke (1998) reduced by a factor of 10 to account for
wind clumpiness. For effective temperatures $T_{\mathrm{eff}}<23\,300$ K,
which can be achieved in between pulses due to the expansion of the star, we
employ the maximum between the Vink et al. (2000, 2001) and Nieuwenhuijzen &
de Jager (1990) wind mass loss rates. We turn off wind mass loss during the
(physically brief) dynamical phases of evolution: PPI-driven dynamical mass
ejections are the only source of mass loss in these phases. Uncertainties in
the wind mass loss rate can have an impact on the core structure (Renzo et
al., 2017). The main effect of varying the wind algorithm in our H-free models
is to change the mapping of the initial He core mass $M_{\mathrm{He,init}}$ to
the He core mass at the onset of the pulses. Stronger wind mass loss would
inevitably reduce the mass of PPI-produced CSM by removing the mass before the
instability. Moreover, the wind velocity might differ from the ejection
velocity in pulses, thus resulting in a different CSM density profile. In
Farmer et al. (2019) we explored the uncertainties in the wind mass loss
rates, varying their functional form (including empirically determined rates
from Nugis & Lamers 2000 and Tramper et al. 2016), and efficiency factors. We
found that these uncertainties do not significantly influence the range of
possible BH remnant masses and the PPI+CC properties when expressed as a
function of the carbon-oxygen core mass. For a recent comparison of wind mass
loss rates for Wolf-Rayet stars, we refer the interested reader to Yoon et al.
(2017) and to Woosley (2019) for a grid of models spanning the PPI-regime.

To follow the dynamical evolution of the pulses when they occur, we use MESA’s
Riemann HLLC solver (Toro et al., 1994; Paxton et al., 2018). We determine the
dynamical stability of the star based on the adiabatic index
$\Gamma_{1}=\partial\log(P)/\partial\log(\rho)|_{s}$. Regions of the star with
$\Gamma_{1}>4/3$ are formally stable (e.g., Kippenhahn et al., 2013). However,
this is a local quantity. To create a global metric descriptive of the entire
star, we follow Stothers 1999 in defining a volumetric pressure-weighted
average adiabatic index

$\langle\Gamma_{1}\rangle\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{\int\Gamma_{1}P\,d^{3}r}{\int
P\,d^{3}r}\equiv\frac{\int\Gamma_{1}\frac{P}{\rho}\,dm}{\int\frac{P}{\rho}\,dm}\
\ ,$ (1)

where $P$, $\rho$ are the local pressure and density, and we used the
continuity equation to transform the volumetric integral into an integral over
the mass domain. Weighting the _local_ $\Gamma_{1}$ with $P$ makes the average
$\langle\Gamma_{1}\rangle$ a dynamically relevant quantity, and guarantees
that the inner regions contribute more to the average. Whenever
$\langle\Gamma_{1}\rangle=4/3+0.01$, i.e., slightly before the stellar
structure becomes formally unstable, we switch to a hydrodynamical treatment
of the evolution and turn off the stellar winds (see also Marchant et al.
2019).

After a pulse, if the internal structure of the star meets the criteria
specified in Marchant et al. (2019) to conservatively ensure hydrostatic
equilibrium has been recovered, we excise the material moving faster than the
local escape velocity and create a new star with the entropy, chemical
composition, and mass of the layers remaining bound. Even for nonpulsating
models, we turn on the hydrodynamics to follow the onset of core-collapse,
when the core temperature rises above $T_{c}\gtrsim 10^{9.6}$ K.

We adopt a 22-isotope nuclear reaction network (approx21_plus_co56.net), which
is sufficient to trace the energy output during the relevant burning phases
but not the detailed nucleosynthesis (e.g., Farmer et al., 2016, 2019).

We assess convective stability using the Ledoux criterion, and adopt a mixing
length parameter of $\alpha_{\mathrm{MLT}}=2.0$. We consider semi-convective
mixing with an efficiency $\alpha_{s}=1.0$, whilst neglecting thermohaline
mixing. We assume an exponential undershooting and overshooting with
parameters333cf. Equation 2 in Paxton et al. 2011 and the MESA documentation
for the definition of f and f0. (f,f0)=(0.01, 0.005) for all convective
regions. We follow the approach of Marchant et al. (2019) based on Arnett
(1969) for the time-dependence of the convective velocity. This is required to
compute dynamical phases of the evolution with timesteps shorter than the
convective turnover timescale (Renzo et al., 2020).

We stop our evolution either at the onset of CC or at the onset of a PISN. We
define the former as when the infall velocity anywhere in the model exceeds
$1000\,\mathrm{km\ s^{-1}}$ (Woosley et al., 2002). For the latter we check
that the total energy (including the kinetic term) of the star is positive and
that the minimum radial velocity is non-negative. These conditions guarantee
that the star is unbound and there is an outflow of matter.

We define the CO core mass ($M_{\mathrm{CO}}$) of our models as the maximum
mass coordinate where the mass fraction of He is lower than $Y<0.01$ at He
core depletion, i.e. when the central mass fraction of He reaches
$X_{c}(^{4}\mathrm{He})<10^{-5}$. This allows us to characterize our models
with one single CO core mass, while the actual amount of CO-rich mass might
change because of the evolution. We define the iron core mass
($M_{\mathrm{Fe}}$) as the outermost mass location where the abundance of
${}^{28}\mathrm{Si}$ is lower than 0.01 and the abundance of elements with
mass number $A>46$ is greater than 0.1. We use discuss $M_{\mathrm{Fe}}$ only
at the onset of CC. We refer the reader to Appendix A for the description of
the numerical resolution and a quantitative assessment of its impact on our
results.

The input files (inlists) and customized routines added to the code
(run_star_extras.f) needed to reproduce our results are available at
https://doi.org/10.5281/zenodo.3406356. We also provide our numerical results
for the evolution and final structure of each of our models, including a
customized output file storing averaged information for each layer moving
beyond its local escape velocity during PPI-driven mass loss episodes. The
possibility of fallback is neglected in these files, even though our MESA
models allow for it. Such files can be used as inputs for hydrodynamic studies
of the CSM structure produced by these stellar models.

Figure 2: Final BH masses as a function of the initial He core mass. The scale
in the horizontal direction is logarithmic. The colors in the background
indicate the approximate range for each evolutionary path, see also Section 3.
The right panel shows the masses inferred from the first ten binary BH mergers
detected by LIGO/Virgo, with a red shade to emphasize the overlap between PPI
and CC, and green and blue hatches to indicate the fate of the progenitor in
different BH mass ranges.

##  3 Overview of the evolution of the progenitors

Figure 2 shows the BH masses resulting from our grid as a function of the
initial He core mass ($M_{\mathrm{He,init}}$, bottom axis) and approximate
maximum CO core mass reached during the evolution ($M_{\mathrm{CO}}$, top
axis). Both can decrease because of PPI mass-loss episodes toward the end of
the evolution. We estimate the BH mass as the mass coordinate where the
binding energy of the collapsing star reaches $10^{48}\,\mathrm{ergs}$, to
allow for the possibility of mass loss during the CC from, either a weak
explosion (Ott et al., 2018; Chan et al., 2020), or ejection of a fraction of
the envelope due to neutrino losses(Nadezhin, 1980; Lovegrove & Woosley,
2013). This estimate is typically within a few $0.01\,M_{\odot}$ of the total
final mass of the He star. We do not account for other energy loss terms
during the core collapse, such as neutrinos themselves which might carry away
(part of) the core binding energy. This effect is typically estimated to be
$\lesssim$ 10% of the precollapse core rest mass energy (e.g., O’Connor & Ott,
2011; Belczynski et al., 2016; Spera & Mapelli, 2017), and can shift our BH
mass estimates further down.

The colored background in the left panel of Figure 2 indicates approximately
the evolutionary path for the corresponding mass range. The four possibilities
are summarized as follows, in order of increasing initial He core mass:

#### CC:

Relatively low mass He cores end their lives in a core collapse (CC, blue on
the left of Figure 2) event without losing mass to pair-production driven
pulses. For these models, the layers which are unstable to pair production (if
any) are not massive enough to cause an episode of mass ejection. In this mass
range, the outcome of core-collapse is most likely BH formation, possibly
associated with a weak SN with large fallback (Ott et al., 2018; Kuroda et
al., 2018; Chan et al., 2018, 2020). We return on the “explodability” of our
grid of models in Section 7.

#### PPI+CC:

With increasing $M_{\mathrm{He,init}}$, the pair instability becomes
progressively more violent. The energy release by thermonuclear explosions
causes significant radial expansion. Increasing further in mass, models
experience one or more mass loss episodes, before the core is stabilized by
the consumption of fuel and entropy losses to neutrinos, and the stars finally
collapse (PPI+CC, green in Figure 2).

#### PISN:

For $80\,M_{\odot}\lesssim M_{\mathrm{He,init}}\lesssim 200\,M_{\odot}$, our
models are completely disrupted in a PISN, and produce no remnant (yellow
vertical area in Figure 2). Our lowest mass model going PISN and leaving no
remnant has $M_{\mathrm{He,init}}=80.75\,M_{\odot}$, corresponding to a
maximum CO core mass of $\sim$ $55\,M_{\odot}$ (see also Farmer et al., 2019).

#### CC:

For extremely massive cores, $M_{\mathrm{He,init}}\gtrsim 200\,M_{\odot}$, the
energy release by the explosive thermonuclear burning triggered by the pair
instability is insufficient to fully disrupt the star. This happens because
most of that energy is used to photodisintegrate the nuclear ashes and lost to
neutrinos, instead of becoming kinetic energy of the stellar gas (Bond et al.,
1984; Fryer et al., 2001). Therefore, models above a certain threshold reach
CC without any PPI-driven mass loss (blue area on the right of the panel of
Figure 2).

In Figure 2 the transition between the CC and PPI+CC is smooth, and we have
avoided quantifying the boundary between the low mass CC and PPI+CC because of
the subtleties in the definition of “pulse”. We outline three physically
motivated definitions, each one shifting the CC/PPI+CC boundary, in Section 5.

##  4 Resulting BH masses

The PISN BH mass gap is denoted by the hatched region in the left panel of
Figure 2. The lower and upper edge of the gap can be read from the y-axis.
With our numerical setup, we find a maximum BH mass below the PISN gap of
$\max\\{M_{\mathrm{BH}}\\}\simeq 45\,M_{\odot}$, in good agreement with the
lower boundary of the gap from previous studies Woosley (2017); Marchant et
al. (2019); Leung et al. (2019); Woosley (2019); Farmer et al. (2019). At the
upper-end, the PISN BH mass gap is closed by the photodisintegration
instability causing the direct collapse of an initially
$M_{\mathrm{He,init}}=200\,M_{\odot}$ He core, which builds up a CO core of
$M_{\mathrm{CO}}\simeq 114\,M_{\odot}$ and eventually forms to a BH of
$125\,M_{\odot}$. This value corresponds very closely to the final He core
mass of this model, after wind mass loss. This upper boundary too is in good
agreement with the results from Woosley et al. (2002); Woosley (2017),
although it is sensitive to the metallicity and uncertainties in the wind mass
loss rate.

We do not expect these boundaries would have varied if our models had a H-rich
envelope (e.g., Woosley, 2017), especially if the progenitor stars evolve in
close binaries which can remove the H-rich envelope long before the PPI. The
combination of the mass loss (Ott et al., 2018; Kuroda et al., 2018; Chan et
al., 2018, 2020) and energy loss to neutrinos at BH formation (e.g., Coughlin
et al., 2018) should be sufficient to unbind any the residual H envelope,
unless the progenitor is a blue super giant with a large binding energy of the
envelope (exceeding $\sim$ $10^{48}\,\mathrm{erg}$, Lovegrove & Woosley 2013).
Such blue supergiant pre-PPI structures might arise from binary mergers (e.g.,
Spera et al., 2019, for a population synthesis study). However, the stellar
structure calculations for the merger of two post-main-sequence stars from
Vigna-Gómez et al. (2019) show extended convective envelopes at the onset of
the instability, which support our expectation that the envelope would easily
be shed at the onset of the instability (see also Appendix B). Whether the
H-rich envelope can contribute to the BH mass or not deserves further
investigation.

The right panel of Figure 2 shows for comparison the individual BH masses of
the binary BH mergers detected to date by LIGO/Virgo444Other events have since
been reported by an independent analysis of the first two observing runs, see
Zackay et al. (2019) and references therein., with the 90% confidence level
uncertainty ranges. The masses of the two BHs in a merger event are not direct
observables, they are instead inferred from the chirp mass and total mass of
the binary. The color of the hatching in the right panel indicates the
possible progenitor evolution (see also Section 3): the red area emphasize the
range of BH masses that can be obtained by CC of a lower mass model, or by
severe PPI mass loss of the most massive PPI+CC models. Its extent to the
lower BH masses is somewhat dependent on the resolution in
$M_{\mathrm{He,init}}$ of our grid. However, in Marchant et al. (2019) we
showed that the minimum BH mass that can be obtained by PPI+CC evolution is
about 10 $M_{\odot}$ because of the production of radioactive material that
can unbind cores that recover hydrostatic equilibrium of lower masses. Most BH
progenitors for the gravitational wave mergers events detected to date are
compatible with encountering the PPI, although we do not expect most of them
to have gone through this evolution because progenitors with sufficient mass
are disfavored by the initial mass function.

##  5 The physics of pulses: cores, radii, and mass ejections

While the nuclear and thermal processes governing the evolution of a star
through pair instability are well understood, the characterization of the
observable properties of such events are not yet as clear. The main reason for
this is that stars do not react “elastically” to the pair instability: instead
the nuclear binding energy released by burning episodes at each pulse is
stored and re-distributed throughout the stellar structure, and there is not a
one-to-one correspondence between what happens in the core and what can be
observed at the surface.

To clarify the distinction between core behavior and observable properties
from the outermost layers of the star, in this section we describe the
physical processes that can be used to give three different physically-
motivated definitions of “pulse”, and illustrate them with an example
$M_{\mathrm{He,init}}=50\,M_{\odot}$ He core (for which we also present a
resolution study in Appendix A). These definitions do not cover all the
possible ways in which pulses can be defined and counted. For example, Woosley
et al. 2007; Woosley 2017 uses the core temperature $T_{c}$ while in Marchant
et al. 2019 we adopted a criterion based on the maximum velocity in the
stellar interior.

###  5.1 Thermonuclear ignition

Historically, studies on pair-instability evolution have focused on the core
of stars. Indeed, the region that becomes unstable because of the _runaway_
production of $e^{\pm}$ is typically deep in the star, and the subsequent
evolution is driven by the explosive burning of oxygen and heavier fuel (e.g.,
Barkat et al., 1967; Rakavy & Shaviv, 1967; Fraley, 1968; Woosley et al.,
2007; Woosley, 2017; Marchant et al., 2019; Leung et al., 2019). This allows
for a definition of a pulse based on the behavior of the deep interior of the
star.

Figure 3: Last year of evolution of the central temperature (blue, bottom
panel), nuclear (red) and neutrino (orange) luminosity (top panel) for a
$50\,M_{\odot}$ He core. The inset in the bottom panel shows the readjustment
of the core to hydrostatic equilibrium after the first pulse, which we resolve
even if this behavior is likely to be influenced by the imposed spherical
geometry.

Figure 3 shows the evolution of the central temperature ($T_{c}$, bottom
panel) and nuclear and neutrino luminosity ($L_{\mathrm{nuc}}$ and $L_{\nu}$
respectively, top panel) during the last year before CC for a
$M_{\mathrm{He,init}}=50\,M_{\odot}$ He core. The horizontal axis shows the
time to the onset of CC on a reversed logarithmic scale.

During the previous evolution (not shown), this model is thermally and
dynamically stable. At $\log_{10}\\{(t_{\mathrm{CC}}-t)/[\mathrm{yr}]\\}\simeq
2$, the core becomes thermally unstable because of the softening of the EOS,
causing the collapse of the core and rise in $T_{c}$ allowing for the
explosive ignition of fuel. The latter can be seen as spikes in the nuclear
luminosity $L_{\mathrm{nuc}}$. The thermonuclear release of energy expands the
core, cooling it adiabatically and causing a temperature drop. Eventually, the
core is stabilized by the loss of entropy to neutrinos and the burning of
nuclear fuel, and it ends its life steadily increasing its core temperature
until the onset of CC.

If we define PPI pulses based on core temperature spikes, or equivalently
spikes in nuclear and neutrino luminosity, then the lowest mass He core
showing hints of pulsational behavior is $M_{\mathrm{He,init}}\simeq
37.5\,M_{\odot}$, corresponding to a final $M_{\mathrm{CO}}\simeq
28\,M_{\odot}$. Even if the local adiabatic index $\Gamma_{1}<4/3$ somewhere
in this model, the volumetric pressure-weighted averaged adiabatic index is
always $\langle\Gamma_{1}\rangle>4/3$ for the entire evolution. The
thermonuclear ignition in the core never results in a _global_ instability of
the star. We find that $\langle\Gamma_{1}\rangle$ crosses the stability
threshold of 4/3 at some point in the evolution only for
$M_{\mathrm{He,init}}>40.5\,M_{\odot}$.

The inset plot in the bottom panel of Figure 3 shows that in a one-dimensional
spherical MESA model the core “bounces” off itself, which was already noted in
Paxton et al. (2018). These readjustments of the core cause secondary burning
episodes (Renzo et al., 2020) which can release further energy and introduce
complications in counting the $T_{c}$ spikes (see also Marchant et al. 2019).
While we resolve in time these bounces by taking timesteps shorter than the
dynamical timescale of the core, it is likely that multidimensional effects
and/or off-center energy release would affect them significantly. Even
counting the core oscillations as one individual pulse, our
$M_{\mathrm{He,init}}=50\,M_{\odot}$ model exhibits tens of thermonuclear-
ignition pulses.

This behavior of the core of very massive stars encountering the pair
instability is well established (e.g., Barkat et al., 1967; Woosley, 2017;
Marchant et al., 2019; Leung et al., 2019). However, since these processes
happens deep inside the optically thick layers of the star, their only direct
observable is the rapid variation of orders of magnitude of the neutrino
luminosity during each pulse (Fryer et al. 2001, and possibly during the post-
pulse bounces). However, because of the rarity of such massive stars in the
local Universe, such variations in the neutrino luminosity are unlikely to be
easily observed.

Figure 4: Radial evolution of a $M_{\mathrm{init}}=50\,M_{\odot}$ He core
through PPI pulses. The radius of the bound material (red, plotted only during
dynamical phases) oscillates because of the PPI. The photospheric radius
(orange) reaches $10^{4}\,R_{\odot}$ (dashed line) during the dynamical phase
of evolution, i.e., the location at optical depth 2/3 is in the material
already ejected.

###  5.2 Radial expansion

In Section 5.1 we discussed a definition of a PPI pulse based on the
thermonuclear behavior deep inside the core. However, for
$M_{\mathrm{He,init}}>41\,M_{\odot}$, the nuclear binding energy released by
the thermonuclear explosions deep down can have an observable impact on the
surface of the star. We emphasize that either because of other evolutionary
processes, or because of a previous PPI mass-loss episode, the surface can be
a He rich layer for these stars. Therefore, we can use our He core models to
define a PPI pulse based on surface properties, assuming the H-rich layers
have been lost before.

The thermonuclear burning injects energy into the core and drives a pulse
wave, which propagates down the decreasing density profile of the star and
eventually steepens into a shock. The core, post-explosion, readjusts and can
contribute to driving secondary shocks. There is a small range in mass,
$41\,M_{\odot}\lesssim M_{\mathrm{He,init}}\lesssim 42\,M_{\odot}$ in which
these shocks, which can often catch up with each other below the stellar
surface, are not energetic enough to dynamically unbind any significant amount
of matter (see also Section 5.3). Nevertheless, even in this mass range, they
produce a potentially observable radial expansion of the star. For models more
massive than $M_{\mathrm{He,init}}\gtrsim 42\,M_{\odot}$, the energy released
in the thermonuclear explosion also cause the ejection of material (see
Section 5.3).

Figure 4 shows the radial evolution of our $50\,M_{\odot}$ example. The orange
line shows the photospheric radius (defined as the location where the optical
depth is 2/3). For most of the evolution in hydrostatic equilibrium (until
$t_{\mathrm{CC}}-t\simeq 10^{-2}$ years), the stellar radius is on the order
of the solar radius ($R_{\odot}\simeq 6.9\cdot 10^{10}\,\mathrm{cm}$) or less.
As the star contracts and approaches the instability, we switch to the HLLC
solver at around $t_{\mathrm{CC}}-t\lesssim 1$ year, when the thicker red line
appears in Figure 4. This line shows the radius of the bound material
$R(v<v_{\mathrm{esc}})$ and is plotted only when the hydrodynamics is on.

We follow the thermal contraction of the star due to the pair instability, and
at $t_{\mathrm{CC}}-t\simeq 10^{-2}$ years a shock wave propagating from the
core causes a radial expansion by two orders of magnitude on a dynamical
timescale. The material remaining bound to the star extends to $\sim$
$10^{13}$ cm. For our $M_{\mathrm{He,init}}=50\,M_{\odot}$ model, the pulse
also ejects matter and the ejected layer extends beyond $10^{14}$ cm. For
numerical stability reasons we cap the radii at $10^{4}\,R_{\odot}\simeq
6.9\cdot 10^{14}\,\mathrm{cm}$ (horizontal dashed line in Figure 4), and treat
this limit as an open boundary555In none of our models is such an upper limit
in radius reached by the bound material. We discuss the ejected matter more
extensively in Section 5.3.

The ejected layer can obscure the bound surface of the star: the location
where the optical depth is 2/3 extends all the way to the outermost layers of
our Lagrangian mesh. However, we emphasize that the structure of the material
moving faster than the escape velocity should be recomputed accounting for
radiative losses for a better determination of the photosphere. It is possible
that such stars would exhibit large radius differences at different
wavelengths, with some that might even appear red during their maximal radial
expansion. For this particular model, the photospheric radius does not have
time to recover its prepulse value, since the radial expansion only started
days before the final core-collapse. More massive models have more violent
pulses that drive the inner core farther out of thermal equilibrium and for
which it takes longer to recover the condition for further (explosive or
stable) nuclear burning (see also Section 6): this can give time to the
photospheric radius to decrease again.

The bound radius instead experiences large oscillations between $10^{11}$ cm
and $10^{13}$ cm (the maximal expansion reached initially). In this case,
these might not be directly obervable since they are embedded within the
pseudo photosphere of the ejecta. Models more massive than our example might
have rather long lived phases with large radii, which might have implications
for binary interactions (e.g., Marchant et al., 2019) and wind mass loss
physics.

###  5.3 Ejection of material

Figure 5: Space-time diagram for propagation of the PPI ejecta of a
$50\,M_{\odot}$ He core. The color indicates the density, assuming radial
expansion at constant velocity. The top panel indicates the cumulative amout
of ejected at velocities larger than the escape velocity through pulses only
(i.e., excluding the wind mass loss). The cyan curve shows the radius of the
material instantaneously bound (cf. the red curve in Figure 4).

Although the two definitions of a pulse we introduced in the previous
subsections (based on the core explosive behavior and on the radial expansion,
respectively) might possibly give observable “pulses”, the processes they are
based on do not leave a direct imprint on the CSM structure, nor on the
remnant BH mass. Observational confirmation of the occurrence in nature of
PPI+CC evolution (and possibly PISN) is most likely to come from either
observations of transients which can probe the CSM around the exploding star
and/or the distribution of BH masses probed through gravitational waves. It is
therefore worth giving a definition of pulse based on the nonterminal ejection
of material from the stars: the ejecta carry away mass, decreasing the final
BH mass and shaping the CSM structure.

Our simulations produce output for the ejecta at each timestep. The top panel
of Figure 5 shows the cumulative mass lost to PPI-driven pulses for our
example $M_{\mathrm{He,init}}=50\,M_{\odot}$ He core. In this specific model,
the total (H-free) ejecta mass is about $1.2\,M_{\odot}$ by the end of the
evolution. Had our star retained an H-rich envelope until the onset of the
first pulse, the remaining H-rich envelope at the onset of the instability
would likely add to the amount of mass in the CSM (see also Appendix B).

The bottom panel of Figure 5 shows the density distribution around the star as
a function of distance from the star (y-axis) and time until the final CC
(x-axis). The cyan line shows the radius of the bound material (cf. the red
curve in Figure 4), which we assume to be the initial radius from which the
ejecta are launched. To compute the CSM density we assume propagation of the
ejecta at constant velocity. We use the velocity at the time the material
first exceeds the local escape velocity as computed by MESA, and it is
typically a few thousand ${\mathrm{km\ s^{-1}}}$. We return to the ejecta
velocity in Section 6. Assuming a constant velocity for the propagation of the
ejecta corresponds to neglecting radiative cooling, internal collisions of the
ejecta, and multi-dimensional effects (e.g., Chen & Woosley, 2019), and we
discuss it here only for illustration purposes. Our output files666Publicly
available at https://doi.org/10.5281/zenodo.3406356. contain the amount of
mass ejected, its center-of-mass velocity, chemical composition and thermal
state (averaged by mass over all the mesh points that exceed the local escape
velocity in the current timestep), which could be used as input for more
detailed simulations to predict the CSM structure around PPI+CC models. This
ejecta output neglects the possibility of fallback, which however could be
implemented when using these files as input for hydrodynamical simulations of
the CSM.

The CSM structure shown in Figure 5 for our example model shows H-poor/He-rich
CSM starting from $\sim$ $10^{13}\,\mathrm{cm}$ and extending out to $\sim$
$10^{14}\,\mathrm{cm}$. The CSM densities reach $10^{-5}-10^{-4}\,\mathrm{g\
cm^{-3}}$. These value are typical for the models in our grid, and fall in the
range of CSM distances and densities inferred from transient observations
(e.g., Gomez et al., 2019).

However, mass ejections that happen in subsequent timesteps (possibly with no
mass ejected in between) in MESA might not be physically distinct events. To
count the mass-ejection events, we need to group mass ejections in timesteps
separated by less than a dynamical timescale in one individual event.
Moreover, one mass ejection event can last several dynamical timescales, for
example the final mass ejection and full disruption of a PISN is expected to
produce a long transient (e.g., Gal-Yam et al., 2009). We estimate the
dynamical timescale as the free-fall timescale
$\tau_{\mathrm{ff}}=2\pi\sqrt{R_{\mathrm{photo}}^{3}/GM_{\mathrm{bound}}}$
where G is Newton’s constant, and $R_{\mathrm{photo}}$ and
$M_{\mathrm{bound}}$ are the (time-dependent) photospheric radius and mass
gravitationally bound to the star. We define the beginning of a mass loss
event as the timestep during which at least $10^{-6}\,M_{\odot}$ has been
removed from the star since the last mass loss event (either in one timestep,
or cumulatively). We require each mass ejection event to last at least one
free fall timescale (calculated at the beginning of the pulse), and define its
end as soon as the amount of mass to be ejected in the following 100
$\tau_{\mathrm{ff}}$ (now calculated at the end of the pulse) is less than
$10^{-7}\,M_{\odot}$. This last condition allows us to count as a single event
mass ejection episodes that last longer than a dynamical timescale. All
together, these requirements enforce that ejections which numerically happen
in different timesteps separated by less than a free fall timescale are not
counted as separate events.

Adopting this criterion to count the mass ejection events, our
$M_{\mathrm{He,init}}=50\,M_{\odot}$ model only has one mass-ejection episode
(cf. tens of core ignitions, see Section 5.1), which starts roughly
$\sim\,0.015\,\mathrm{years}\simeq 130\,\mathrm{hours}$ before CC.

We can tentatively apply the same threshold to define the beginning of an
explosion for a terminal PISN with full disruption. In this case, the duration
of the PISN events in our grid exceeds months even for the least massive PISN
model with $M_{\mathrm{He,init}}=80.75\,M_{\odot}$ in agreement with previous
studies. Our stopping conditions do not allow models to reach what would
appear as the observational end of a PISN.

##  6 Pulsational pair-instability-generated CSM

Figure 6: Number of mass-ejection events caused by pair instability as a
function of CO core mass. The color shading indicates the approximate range
for each behavior: core collapse without experiencing PPI-driven mass loss
(CC, blue), PPI-driven mass loss (PPI+CC, green), or full disruption in a PISN
(yellow), which we define as one mass loss event. The noisiness is caused by
the occurrence of mass loss event right at the time of the final core
collapse. Figure 7: Velocity profile at the onset of core-collapse for a 50
$M_{\odot}$ He core, which undergoes a PPI mass-loss event while collapsing.
The dashed red line shows the escape velocity profile, the thick blue line
indicates the profile of the bound material, while the thinner line represents
material still on the Lagrangian mass grid, but already beyond the escape
velocity. The gray area indicates the model dependent amount of mass lost to
stellar winds.

We now discus the CSM that can be created by PPI evolution across our grid of
models. In principle, this CSM can be probed by time-domain observations, if
the final core-collapse causes an associated explosion (see Section 7), or
because of collisions between shells ejected at different times (e.g., Woosley
et al., 2007; Woosley, 2017, 2019). Most of the CSM properties do not depend
on which criterion is used to define the beginning or end of the pulses,
except the number of pulses and their duration. For these quantities, we adopt
the definition of Section 5.3 based on the ejection of matter, which is the
most relevant for discussing the CSM structure.

Table C summarizes the time, duration, and amount of mass loss in each event
for all the PPI+CC models in our grid, and Figure 6 shows the number of pulses
as mass-ejection events contributing to the CSM. Models evolving to CC without
any mass ejection have zero pulses. We define full disruption in a PISN as a
one-pulse event, although these would not contribute to the CSM itself. The
color in the background emphasizes the various evolutionary behaviors, but
using the definition from Section 5.3 the mass threshold separating CC (blue)
from PPI+CC (green) evolution is well defined at
$M_{\mathrm{He,init}}=40.5\,M_{\odot}$, correponding to
$M_{\mathrm{CO}}=30\,M_{\odot}$.

The number of pulses is zero at the lower end, and increases up to three
distinct mass-ejection events for the central part of the PPI+CC mass range.
At even higher masses, approaching the PPI+CC/PISN boundary, the number of
pulses decreases again, although the amount of mass ejected increases (see
also Figure 8): this is because pulses become more energetic and consume more
nuclear fuel at once (e.g., Woosley et al., 2007; Chatzopoulos & Wheeler,
2012b; Woosley, 2017, 2019).

The green region shows some noise in the number of pulses: the reason for this
is illustrated in Figure 7, which shows the velocity as a function of
Lagrangian mass coordinate for our $M_{\mathrm{He,init}}=50\,M_{\odot}$ core
at the onset of core collapse. Many models exhibit a similar behavior, with an
outgoing pulse wave at the onset of CC: this means that the PPI mass ejection
is still going on while the Fe core starts collapsing.

Figure 8 summarizes the amount of mass lost to PPI-driven pulses across our
model grid. The bottom panel shows the amount of mass lost per individual
pulse, the pulse number is represented by the color of the filled circles.

The typical amount of mass lost varies from $\lesssim 10^{-3}\,M_{\odot}$ for
the lowest-mass models ejecting some mass, up to $\simeq 20\,M_{\odot}$ (of
He-rich material) at the upper mass end, just below the minimum mass for PISN.
For models producing more than one mass ejection event (i.e., the models for
which also a purple and possibly a red dot are shown), the amount of mass lost
per pulse does not behave monotonically with the pulse number. For
$50\,M_{\odot}\lesssim M_{\mathrm{He,init}}\lesssim 62\,M_{\odot}$,
corresponding to $36\,M_{\odot}\lesssim M_{\mathrm{CO}}\lesssim 43\,M_{\odot}$
the second pulse (purple) ejects more mass than the first (blue), while for
higher masses the first pulse removes more mass than the second.

Figure 8: PPI-driven mass loss as a function of the CO core mass. The top
panel shows the total mass ejected in pulses, the bottom panel shows the mass
lost in individual pulses. The amount of mass lost does not have a monotonic
behavior with pulse number, and spans a wide range of values. The first pulse
is shown as a blue dot, and the second and third, if they occur, are shown as
a purple plus and red cross, respectively. Thin vertical lines connect
multiple pulses for the same $M_{\mathrm{CO}}$. Figure 9: Delay time between
the final core-collapse and the end of each mass loss event in a PPI pulse.
Typical delays are on the order of few months, but they increase steeply with
the initial He core mass of the PPI progenitor, which produce fewer but more
energetic pulses, up to $10^{4}$ years. The top inset magnifies the range of
half a year delays on a days scale. The first pulse is shown as a blue dot,
and the second and third, if occurring, are shown as a purple plus and red
cross, respectively, with a thin blue line connecting the pulses of the same
model.

The top panel of Figure 8 shows the total amount of mass lost to PPI ejecta,
i.e., the sum of the mass ejected in each individual pulse. A trend of more
massive models producing more energetic pulses and driving more mass loss is
evident. We provide a simple fitting formula for the total amount of He-rich
mass lost in PPI-driven events for $33\,M_{\odot}\leq M_{\mathrm{CO}}\leq
56.5\,M_{\odot}$ which produce a CSM mass larger than $\sim$ $0.2\,M_{\odot}$,
shown as a dashed gray line in the top panel of Figure 8:

$\frac{\Delta M_{\mathrm{tot}}}{M_{\odot}}=0.000147\times
10^{0.098M_{\mathrm{CO}}/M_{\odot}}\ \ \mathrm{for}\ \ 33\,M_{\odot}\leq
M_{\mathrm{CO}}\leq 56.5\,M_{\odot}\ \ .$ (2)

The total mass lost to pulses should be added to the amount of mass lost due
to winds to calculate the total mass in the CSM. The density distribution of
the CSM generated by PPI-driven pulses and wind mass loss are likely to be
very different from each other. In cases where stars retain a (loosely bound)
H-rich envelope at the onset of the first pulse, the mass of such an envelope
at the start of the pulses should also be added to the total mass lost in the
first mass loss event.

Figure 9 shows the delay time between the end times mass ejections and the
final CC, as a function of the CO core mass of our models. The timing of the
mass-loss events also spans a large range, from zero (see also Figure 7 for
example) to $\sim$ $10^{4}$ years, corresponding roughly to the Kelvin-
Helmholtz timescale of the most massive PPI+CC. We emphasize that many models
in our grid show submonth delays between the last pulse and the final CC,
which makes them candidates for the detection of CSM interactions in early
observations of the SN explosions (see also inset in Figure 9).

The delay (in between pulses and between each pulse and CC) increases with the
core mass, because more massive models produce more energetic pulses that
drive the star farther from gravo-thermal equilibrium, increasing the amount
of time needed to return to equilibrium after a pulse and resume the final
evolution. While typically for very massive stars the neutrino luminosity
greatly exceeds their photon luminosity $L_{\nu}\gg L$ (they are “neutrino
stars”, Fraley, 1968), this is not always true for the most massive PPI+CC
models. For these, the adiabatic expansion of the core can leads to central
temperature and densities too low for significant neutrino cooling to occur.
Thus, after a pulse begins, these models transition from evolving on a
neutrino-mediated thermal timescale ($\propto GM^{2}/RL_{\nu}$) to a photon-
mediated thermal timescale ($\propto GM^{2}/RL$) in between pulses, which
increases their interpulse time.

Figure 10: Center-of-mass velocity of the layers ejected at each pulse. The
first pulse (blue filled circles) produces typically higher ejection
velocities, and the ejecta move at $\sim$ few thousand ${\mathrm{km\ s^{-1}}}$
which suggests a connection with (some subclasses) of SN Ibn.

Figure 10 shows the center-of-mass velocity of the layers ejected (which is
calculated as the mass-weighted average of the center-of-mass velocity of the
layers ejected at each timestep over the duration of the mass ejection event).
Unlike the other quantities characterizing a pulse, the ejecta velocities we
find do not span orders of magnitude, and are typically a few $\sim$
$1000\,{\mathrm{km\ s^{-1}}}$. This suggests that mass ejection during a PPI
might explain the He-rich circumstellar material required to explain at least
some of the spectra of SN Ibn showing narrow He emission lines (e.g.,
Pastorello et al., 2008), provided that there is a way to excite the ejected
shells. This could be due to collisions between the third and second pulse
within the first ejected shell, or because of a successful explosion at the
final CC (see also Section 7).

The first pulse (blue filled circles) are almost always faster than the later
pulses. Conversely, when they occur, the third pulse (red crosses) is often
faster than the second (purple pluses). Therefore, many models might result in
collisions in between the ejecta which can appear as SN impostors, as noted by
Woosley et al. (2007) and Woosley (2017).

We only report the center-of-mass velocity for the ejected shells, because our
treatment of the ejecta is presently very simplified. Radiation-hydrodynamics
simulations of the ejecta propagation would be desirable to properly quantify
the velocity distribution of the ejecta, and in particular quantify the low-
velocity tail which will appear first in observations of CSM interactions.

##  7 Explodability of CC and PPI+CC models

Does the final CC of a post-PPI star result in a successful explosion? This
question remains relevant because a terminal explosion would potentially allow
us to probe all the layers ejected previously. Whether the final collapse of
PPI+CC evolution is accompanied by a successful, albeit possibly weak,
explosion likely constitutes the biggest uncertainty underlying this study,
and it is not a question that we can settle using only our stellar structure
models. We emphasize that even in the absence of a terminal explosion,
electromagnetic transients from PPI are possible and expected because of
collisions between shells (e.g., Woosley et al., 2007; Woosley, 2017, 2019,
and Sec. 6).

Based on the extrapolation of one-dimensional parametric CC simulations, the
typical expectation is that He cores with masses larger than $\gtrsim
10\,M_{\odot}$ fail to explode (e.g., O’Connor & Ott, 2011; Ugliano et al.,
2012; Ertl et al., 2016), therefore a successful terminal explosion in PPI+CC
models appears unlikely. However, the observations of (BH) X-ray binary
kinematics (e.g., Brandt et al., 1995; Fragos et al., 2009; Atri et al., 2019,
and references therein), and possible spin misalignment in gravitational wave
events (e.g., O’Shaughnessy et al., 2017) might require BH natal kicks. These
can occur only with an explosion and some level of asymmetry in the ejecta
(e.g., Janka, 2013, 2017; Chan et al., 2018, 2020, or possibly in the neutrino
flux). Thus, the current understanding of BH formation cannot yet be
considered final.

As far as we are aware, there are no published multi-dimensional calculations
investigating the final “explodability” of PPI+CC structures. Most of the
available studies (e.g., Ott et al., 2018; Kuroda et al., 2018) target much
lower mass progenitors (up to $\sim$$70\,M_{\odot}$ including the H-rich
envelope) and/or rely on artificially imposed large perturbations or enhanced
neutrino-nucleon scattering to obtain explosions (e.g., Chan et al., 2018,
2020). Therefore, these results may not meaningfully extrapolate to the PPI+CC
regime.

Even in the case of successful but weak explosions, these studies either do
not follow the explosion long enough to study the amount of mass ejected (if
any), or find that only a fraction of the envelope is ejected (e.g., Chan et
al., 2020). While our models do not have any H-rich envelope to begin with,
the He core at precollapse can reach radii similar to that of a H-envelope in
yellow-supergiant stars. Specifically, the radii of bound material
($v<v_{\mathrm{esc}}$) at the onset of CC range between a few to a few
hundreds of solar radii ($\sim 10^{11}-10^{13}\,\mathrm{cm}$). Such extended
He cores, effectively akin to He envelopes in some cases, have significantly
different density profiles compared to the stellar structures explored in the
aforementioned studies. It is unclear whether the usual definition of the
core-envelope boundary, based on abundances, makes sense from a core-collapse
dynamics perspective for our models.

Although also unexplored and uncertain in this regime, even very weak
explosions due to the initial neutrino emission before the formation of an
event horizon (e.g., Nadezhin 1980; Lovegrove & Woosley 2013; Coughlin et al.
2018 and possibly observed by Adams et al. 2017) might potentially be
sufficient to start an electromagnetic transient when the outermost layers hit
the previously ejected mass. Once again, this mechanism is thought to only
eject loosely-bound H-rich outer layers at low velocity, and in red supergiant
progenitors of much lower mass than we explore here. The large radial
extension of some of our precollapse models might possibly result in the
ejection of some equally loosely bound material, but the typical amount of
mass with binding energy sufficiently low is only of order $0.01\,M_{\odot}$.
Should this scenario produce some ejecta, the amount of kinetic energy
available seems unlikely to produce transients detectable at large distances
(Chevalier & Irwin, 2012).

Explosion mechanisms relying on rotational energy (e.g., Mösta et al., 2015)
and jets (e.g., Gilkis et al., 2016; Soker, 2019), which have also been
proposed to explain type Ic SNe showing broad lines (e.g., Barnes et al.,
2018) might also deserve attention in this context. In Marchant et al. (2019)
we showed that PPI-driven mass loss does not dramatically decrease the core
angular momentum. However, we are not aware of any “explodability” criterion
for this kind of core-collapse engine.

One key difference between PPI+CC cores and lower mass stars contributing to
the differences in density profiles is the production of ${}^{56}\mathrm{Ni}$
during pulses, before the final CC. We emphasize that the energy released by
${}^{56}\mathrm{Ni}$ decay is not expected to play a role _during_ core-
collapse, as the half-life of ${}^{56}\mathrm{Ni}$ is much longer than the few
hundred millisecond timescale of CC. Its effect is to modify the initial
conditions for CC during the interpulse evolution of the models.

The typical outermost mass coordinate where we find a mass fraction of
${}^{56}\mathrm{Ni}$ larger than 0.01 is less than $5\,M_{\odot}$, thus we
expect most of the ${}^{56}\mathrm{Ni}$ will eventually be accreted into the
final BH. However, the ${}^{56}\mathrm{Ni}$ is produced weeks, months, or up
to several years before CC. This provides sufficient time for the
${}^{56}\mathrm{Ni}$ to decay and causes the core to expand due to this
heating. This makes the precollapse core structure of PPI+CC models
qualitatively different from lower mass models routinely used for CC
simulations. Detailed simulations of the final CC structure (including large
nuclear networks to capture the precollapse deleptonization, Farmer et al.
2016; Renzo et al. 2017) and explosion are needed to shed light on what the
these qualitative differences mean for the “explodability” of PPI+CC models.
In the most massive PPI+CC models, the decay of the ${}^{56}\mathrm{Ni}$
produced during pulses can ultimately unbind what is left of the core
resulting in a minimum BH mass that can be obtained via PPI+CC (see also
Marchant et al., 2019).

Figure 11 shows the amount of ${}^{56}\mathrm{Ni}$ present in our PPI+CC and
PISN models at the end of the evolution: for PPI+CC models the total mass of
${}^{56}\mathrm{Ni}$ is typically $M_{\mathrm{Ni}}\simeq 0.2-0.4\,M_{\odot}$,
i.e., about one order of magnitude more than what is produced in typical core-
collapse SNe (e.g., Wongwathanarat et al., 2013). This value increases steeply
in the PISN range, reaching about $\sim$ $60\,M_{\odot}$ at the upper end, in
good agreement with Heger & Woosley (2002) and Woosley et al. (2002) results.
However, our calculations are based on a 22-isotope nuclear reaction network
which is known to produce $M_{\mathrm{Ni}}$ deviating by up to a factor of
$\sim$ 1.5$\times$ in either direction from results computed with larger
nuclear reaction networks (regardless of the final fate of the models between
CC, PPI+CC, or PISN).

Figure 11: ${}^{56}\mathrm{Ni}$ mass present in the deep interior of the star at the onset of CC or PISN. The color background is the same as Figure 6 and indicates the approximate evolution (CC for blue, PPI+CC for green, and PISN for yellow). We computed these models using a 22-isotope nuclear reaction network. Table 1: He core mass at the beginning of our simulations ($M_{\mathrm{He,init}}$) and at He depletion ($M_{\mathrm{He,depl}}$), carbon-oxygen core mass ($M_{\mathrm{CO}}$), iron core mass ($M_{\mathrm{Fe}}$), mass location $M_{4}$ where the specific entropy is lower than $4k_{B}N_{A}$, and mass gradient $\mu_{4}\equiv dM/dr|_{s=4k_{B}N_{A}}$ for a representative subset of our models eventually forming a BH at $Z=0.001$. All the models we computed are available at https://doi.org/10.5281/zenodo.3406356.  | $M_{\mathrm{He,init}}$ | $M_{\mathrm{He,depl}}$ | $M_{\mathrm{CO}}$ | $M_{\mathrm{Fe}}$ | $M_{4}$ | $\mu_{4}$
---|---|---|---|---|---|---
| $[M_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$
CC | 35.00 | 30.03 | 25.96 | 1.82 | 2.22 | 0.16
36.00 | 30.78 | 26.64 | 1.92 | 2.30 | 0.18
37.00 | 31.53 | 27.33 | 1.75 | 2.33 | 0.18
38.00 | 32.27 | 28.00 | 1.59 | 2.34 | 0.19
39.00 | 33.02 | 28.71 | 1.65 | 2.39 | 0.21
40.00 | 33.75 | 29.38 | 2.07 | 2.44 | 0.25
41.00 | 34.49 | 30.05 | 1.96 | 2.48 | 0.28
PPI+CC | 42.00 | 35.22 | 30.75 | 1.60 | 2.53 | 0.32
43.00 | 35.95 | 31.41 | 1.67 | 2.57 | 0.36
44.00 | 36.68 | 32.07 | 2.17 | 2.62 | 0.11
45.00 | 37.40 | 32.73 | 2.18 | 2.68 | 0.26
46.00 | 38.12 | 33.42 | 1.76 | 2.72 | 0.19
47.00 | 38.83 | 34.07 | 1.92 | 2.77 | 0.20
48.00 | 39.55 | 34.72 | 1.92 | 2.82 | 0.28
48.25 | 39.73 | 34.90 | 2.17 | 2.82 | 0.23
48.75 | 40.08 | 35.22 | 1.72 | 2.86 | 0.27
49.00 | 40.27 | 35.38 | 1.56 | 2.86 | 0.33
50.00 | 40.98 | 36.04 | 2.33 | 2.90 | 0.28
51.00 | 41.69 | 36.67 | 1.69 | 2.64 | 0.56
52.00 | 42.38 | 37.33 | 2.09 | 2.77 | 0.39
53.00 | 43.09 | 37.98 | 1.96 | 2.57 | 0.21
54.00 | 43.79 | 38.63 | 1.95 | 2.56 | 0.24
54.50 | 44.14 | 38.94 | 1.97 | 2.54 | 0.21
54.75 | 44.31 | 39.09 | 1.84 | 2.51 | 0.29
55.00 | 44.49 | 39.25 | 1.80 | 2.42 | 0.30
56.00 | 45.19 | 39.89 | 2.20 | 2.25 | 0.52
57.00 | 45.88 | 40.53 | 2.18 | 1.90 | 0.45
58.00 | 46.57 | 41.17 | 2.21 | 2.56 | 0.45
59.00 | 47.26 | 41.80 | 2.25 | 2.26 | 0.49
60.25 | 48.13 | 42.59 | 2.18 | 2.31 | 0.48
61.75 | 49.14 | 43.51 | 1.65 | 2.24 | 0.29
62.75 | 49.80 | 44.14 | 1.99 | 2.25 | 0.42
63.75 | 50.48 | 44.76 | 2.03 | 2.28 | 0.32
64.25 | 50.81 | 45.05 | 1.68 | 2.24 | 0.41
65.75 | 51.81 | 45.96 | 2.04 | 2.25 | 0.33
66.25 | 52.15 | 46.28 | 1.63 | 2.24 | 0.39
67.25 | 52.80 | 46.89 | 1.94 | 2.25 | 0.29
68.25 | 53.46 | 47.50 | 2.05 | 2.18 | 0.27
69.25 | 54.12 | 48.10 | 2.04 | 2.21 | 0.31
70.25 | 54.78 | 48.70 | 2.02 | 2.17 | 0.27
71.00 | 55.27 | 49.16 | 2.03 | 2.25 | 0.32
72.00 | 55.91 | 49.72 | 2.31 | 2.16 | 0.23
73.00 | 56.56 | 50.35 | 2.16 | 2.19 | 0.42
74.75 | 57.68 | 51.39 | 1.76 | 2.02 | 0.73
75.00 | 57.85 | 51.52 | 2.02 | 2.12 | 0.25
76.25 | 58.64 | 52.24 | 2.15 | 2.15 | 0.39
77.00 | 59.13 | 52.72 | 2.19 | 2.15 | 0.52
80.00 | 61.03 | 54.45 | 2.09 | 2.31 | 0.53
PISN
CC | 200.00 | 125.26 | 114.23 | 31.50 | – | –
220.00 | 134.22 | 122.64 | 10.78 | – | –
250.00 | 147.03 | 134.60 | 5.77 | – | –

Table 1 lists, for a representative subset of models, quantities commonly used
to determine the “explodability” of a stellar model. These focus on the
innermost layers of the star, but already include the effect of the
precollapse decay of ${}^{56}\mathrm{Ni}$ produced during pulses in this
region of the star. Specifically, we report the initial He core mass
$M_{\mathrm{He,init}}$ and its value at the end of He core burning777Since our
models do not have a H-rich envelope, this is the total mass when the central
He mass fraction reaches $X_{c}(^{4}\mathrm{He})<10^{-5}$. We emphasize that
this value is sensitive to the adopted wind mass loss rates.
$M_{\mathrm{He,depl}}$, the CO core mass $M_{\mathrm{CO}}$, the final iron
core mass $M_{\mathrm{Fe}}$, defined as the outermost location where the mass
fraction of ${}^{28}\mathrm{Si}\leq 0.01$ and the mass fraction of Fe-group
elements, i.e., with more than 46 nucleons, exceeds 0.1; and the two
parameters proposed by Ertl et al. (2016). These are the mass coordinate
$M_{4}$ at which the specific entropy reaches $4k_{B}N_{A}$, where $k_{B}$ is
Boltzmann’s constant and $N_{A}$ is Avogadro’s number, and the mass gradient
at this location $\mu_{4}$ (Eq. 6 in Ertl et al., 2016). Since the nuclear
reaction network we employ does not allow for detailed treatment of the
electron captures and $\beta$ decays, which determine the final electron-to-
baryon ratio, we avoid listing the compactness parameter (e.g., O’Connor &
Ott, 2011), which is sensitive to these modeling assumptions (Farmer et al.,
2016; Renzo et al., 2017). We caution that three-dimensional simulations
(e.g., Ott et al., 2018; Kuroda et al., 2018; Chan et al., 2018, 2020) might
give a different outcome than 1D parametric simulations used to assess the
explodability of grids of models (e.g., O’Connor & Ott, 2011; Ugliano et al.,
2012; Müller, 2019; Couch et al., 2019).

Figure 12: CC and PPI+CC models on the Ertl criterion plane for
”explodability”. The red area indicates the uncertain boundary region between
successful explosions and neutron star formation (white) and collapse with BH
formation (gray). The dotted rectangle in the bottom left indicates the range
originally showed in Ertl et al. (2016), most of our models require
extrapolating outside this range. Blue dots correspond to CC models, while
green dots show PPI+CC models.

The iron core mass is typically below $\sim$ $2.5\,M_{\odot}$ below the PISN
BH mass gap. The values for models above the PISN BH mass gap are sensitive to
the amount of nuclear burning going on before the stopping criterion based on
the infall velocity is reached, and because of their large mass, the entropy
is larger than $4k_{B}N_{A}$ throughout these stars. Because of the as of yet
insufficient understanding of black hole formation it is hard to predict
whether these models would give a successful, albeit possibly weak, explosion.
We would expect that if successful explosion can happen, they would be powered
by fallback accretion.

Figure 12 shows our CC (blue) and PPI+CC (green) on the plane used by Ertl et
al. (2016) to determine the “explodability”. They determine it using 1D
parametric simulations aiming at reproducing SN1987A starting with single star
progenitors of initial mass $15-20\,M_{\odot}$. The red area indicates the
region where some of their engines produce explosions and neutron star
formation, while others produce failed explosions and BHs (see their Tab. 2
and Fig. 8). Above this (gray area), Ertl et al. (2016) predict BH formation
with a failed explosion, while below they predict successful explosions with
neutron star formation. Applying this criterion to our models requires
extrapolating significantly outside the range originally explored by Ertl et
al. (2016), shown by the dotted rectangle in the bottom left corner. We omit
from the plot the models models above the PISN BH mass gap in our grid for
which $M_{4}$ and $\mu_{4}$ are zero. Most of our PPI+CC models fall in the
failed explosion region according to the Ertl et al. (2016) criterion. A few
enter in the region that depends on the assumed calibration for their model,
and some cross marginally into the successful explosion and neutron star
formation region. The more explodable models correspond to lower initial He
core masses, but there is a large scatter in the trend. All of our CC models
fall in the region for BH formation without an explosion.

Figure 13: Surface composition at the onset of core-collapse for the PPI+CC
and CC models. ${}^{4}\mathrm{He}$, ${}^{12}\mathrm{C}$, ${}^{16}\mathrm{O}$
are shown as red filled circles, cyan triangles, and blue crosses,
respectively. The dot-dashed horizontal lines of the same color mark the
initial mass fraction for these elements for our $Z=0.001$ model grid. The
colored background indicates the approximate evolution of the star (cf. Figure
2 and Section 3).

Figure 13 shows the surface mass fractions of ${}^{4}\mathrm{He}$,
${}^{12}\mathrm{C}$, ${}^{16}\mathrm{O}$ at the onset of core collapse for our
CC and PPI+CC models. Should the CSM evolve to be (partially) optically thin
at the time of the final CC and in case this results in a successful
explosion, the surface composition of the star and the mixing processes (e.g.,
Dessart et al., 2012) would determine the spectral type of the SN. Most PPI+CC
models experiencing a significant amount of mass loss show He-poor surfaces,
and enhanced carbon (and to a smaller extent) oxygen mass fractions,
corresponding to type Ic SNe. Because of the radial expansion caused by the
pulses, some of these progenitors might look like extended and cool objects at
the onset of collapse, rather than compact and hot progenitors. Conversely, if
the He-rich CSM is optically thick, it might obscure the progenitor star and
the embedded explosion, and we expect that re-processing of the photons by the
shell would produce He lines (possibly narrow and in emission) corresponding
to a type Ib(n) SN.

We omit from Figure 13 the PISN models because our stopping condition for
these conservatively ensures the full star is unbound, but it might not
correspond to the “beginning” of the explosion. Nevertheless, most PISN models
start exploding while still retaining some He at their surface (corresponding
again to type Ib SNe) with the wind mass-loss rates and metallicity adopted
here.

Figure 13 also shows a trend with the initial He core mass: the larger the
initial $M_{\mathrm{He,init}}$, the more mass is lost to winds and PPI, the
lower the surface He mass fraction and the higher the carbon and oxygen mass
fractions. In the mass range $50\,M_{\odot}\lesssim
M_{\mathrm{He,init}}\lesssim 75\,M_{\odot}$, corresponding roughly to the
region where we find three distinct PPI-driven mass loss events, the predicted
surface abundances appear more noisy.

##  8 Comparison to selected supernovae

Stars experiencing the PPI+CC evolution should intrinsically be rare because
of the large initial mass necessary to build up a sufficiently massive core.
To produce a significant amount of CSM via PPI-driven pulses, our results
suggest that the He core mass needs to initially exceed
$M_{\mathrm{He,init}}\gtrsim 42\,M_{\odot}$, corresponding to
$M_{\mathrm{CO}}\gtrsim 31\,M_{\odot}$. Therefore, the rate of observed
transients that can be interpreted as signatures of PPI evolution should be
small. Possibly for this reason an unambiguous detection of PPI+CC/PISN in
time-domain surveys is not yet available, although the physical mechanism
underlying this phenomenon is well understood. We consider here a few notable
and recent H-less type I SNe that have been proposed as PPI+CC candidates.

#### PTF12dam:

Tolstov et al. (2017) modeled the H-less (type I) superluminous supernova
PTF12dam as powered by the combination of ${}^{56}\mathrm{Ni}$ decay and CSM
interaction. They proposed that combination of energy sources invoking the
following scenario: first the H-rich envelope is removed by stellar winds,
then the PPI pulses produce $\sim$ $20-40\,M_{\odot}$ of CSM before the final
CC synthesizes and ejects $M_{\mathrm{{}^{56}Ni}}\simeq 6\,M_{\odot}$ of
radioactive material. Our results, albeit computed with a small nuclear
reaction network, never produce this combination of CSM mass and
$M_{\mathrm{{}^{56}Ni}}$: PPI ejecta exceeding $20\,M_{\odot}$ are found only
for $M_{\mathrm{He,init}}\gtrsim 75\,M_{\odot}$ or equivalently
$M_{\mathrm{CO}}\gtrsim 51\,M_{\odot}$ (cf. Figure 8), and only about $\sim
0.2\,M_{\odot}$ of ${}^{56}\mathrm{Ni}$ is synthesized for PPI+CC models.
Assuming that the final CC proceeds similarly as for lower mass stars, we
expect it would add $\sim$ 0.03-0.05$\,M_{\odot}$ of ${}^{56}\mathrm{Ni}$
(e.g., Wongwathanarat et al., 2013), which does not help to reach the
$M_{\mathrm{{}^{56}Ni}}$ claimed. An initial He core mass exceeding
$M_{\mathrm{He,init}}\gtrsim 140\,M_{\odot}$, or $M_{\mathrm{CO}}\gtrsim
87\,M_{\odot}$, is required to reach the amount of radioactive material
required by Tolstov et al. (2017), which would put the model in the PISN range
where we do not expect CSM from PPI.

#### iPTF16eh:

Lunnan et al. (2018) detected a time and frequency varying MgII line in the
spectrum of the type I superluminous supernova iPTF16eh. They interpreted it
as a light-echo of the explosion bouncing off a layer of CSM at $r\simeq
3.5\cdot 10^{17}\,\mathrm{cm}$ moving at $\sim$ $3\,300\,{\mathrm{km\
s^{-1}}}$, implying an ejection $\sim$ $30\,\mathrm{years}$ before the final
CC. Based on these CSM properties and the models from Woosley (2017), they
inferred a progenitor with $M_{\mathrm{He,init}}\simeq 50-55\,M_{\odot}$ (or
equivalently an initial total mass $\sim$ 115$\,M_{\odot}$). Our models are in
overall agreement with the results from Woosley (2017) used by Lunnan et al.
(2018) to interpret iPTF16eh, although the delay time and ejecta velocity
would agree better with a slightly more massive progenitor, with
$M_{\mathrm{He,init}}\simeq 60-65\,M_{\odot}$, i.e., $M_{\mathrm{CO}}\simeq
43-46\,M_{\odot}$.

#### SN2016iet:

Gomez et al. (2019) analyzed the double-peaked peculiar type I SN 2016iet.
They explored several scenarios (PISN, CSM interaction, and central engine) to
power its light curve. This event showed an unusually high Ca/O ratio, and
extreme offset from the nearest galaxy of $\sim$ 16 kpc, however H$\alpha$
lines appear in the spectra beyond 400 days, possibly indicating local star
formation activity. They also detected a possible light-echo from a H- and He-
poor shell moving at few thousand ${\mathrm{km\ s^{-1}}}$. Regardless of the
scenario assumed, they inferred a large progenitor mass with a CO core
$55\,M_{\odot}\lesssim M_{\mathrm{CO}}\lesssim 120\,M_{\odot}$. The model they
favor to explain the light curve combines the signal from the shock cooling of
the prompt explosion (first peak) and CSM interactions (second peak), but
requires $\sim$ $35\,M_{\odot}$ of CSM, although this value is considerably
uncertain (Gomez et al. private communication). Both the inferred presence of
a shell of H- and He-poor material and the claimed progenitor and CSM masses
suggest PPI+CC as a viable scenario for the progenitor of SN2016iet. Several
models with initial He core mass $M_{\mathrm{He,init}}\gtrsim 50\,M_{\odot}$
produce PPI-driven pulses with mass, timing, and velocity within a factor of
about two from the values inferred by Gomez et al. (2019). However, reaching
that total amount of CSM would either require the progenitor to be at the very
edge of the PISN regime, for which we find long interpulse delays (cf. Figure
9) and also the last pulse tends to produce little mass loss (cf. Figure 8).
Alternatively, relaxing the requirement to match the total CSM mass budget by
allowing for a contribution of the stellar wind to the CSM (e.g., because of
the wind in between pulses running into a slower-moving previously-ejected
shell), models with $60\,M_{\odot}\lesssim M_{\mathrm{He,init}}\lesssim
70\,M_{\odot}$, i.e., $43\,M_{\odot}\lesssim M_{\mathrm{CO}}\lesssim
49\,M_{\odot}$, produce pulses removing larger amounts of mass in the final
few years of the progenitor’s life. This might produce a better agreement with
the observed features of SN2016iet. If that were the case, this event might be
the birth of one of the most massive BHs predicted below the PISN mass gap,
cf. Figure 2. At https://doi.org/10.5281/zenodo.3406356, we provide models
computed at our fiducial metallicity value and at the metallicity of the
galaxy at $16\,$kpc from SN2016iet ($Z=0.00198\simeq 0.14Z_{\odot}$), although
it is likely that a dimmer galaxy, possibly with different $Z$, is coincident
with and presently outshined by SN itself. These can provide input for more
detailed calculations of the CSM structure needed to compare with SN2016iet.

#### SN2006jc, PS15dpn and other narrow-line SNe:

Two out of the three SNe we considered above are super-luminous, however the
final collapse of a PPI+CC progenitor or PISNe does not need to be
superluminous (Woosley, 2017). The PPI is just one possible mechanism to
create CSM, which can produce extreme luminosities by generating radiation
from the kinetic energy of the ejecta and/or narrow emission lines (even if
the luminosity does not reach extreme values). The detection of narrow H lines
determines the classification of a SN as a type IIn, while the detection of
narrow He emission lines determines the classification as type Ibn. Both kinds
of event are too common to be entirely explained with PPI+CC progenitors, and
it is likely that both classes contain events with a diversity of physical
mechanisms (e.g., Pastorello et al. 2008 but see also Hosseinzadeh et al.
2017). Nevertheless, it is possible that at least some of these events might
correspond to the observational counterpart of the death of PPI+CC
progenitors. In particular, our simulations can produce several solar masses
of H-free CSM moving at a few thousand ${\mathrm{km\ s^{-1}}}$, which
correspond to the width of the He lines detected in some SN Ibn without any
fine-tuning required. Even if the detection of narrow lines alone is not
sufficient to associate a specific SN to a PPI event, combining evidences from
previous coincident transients, large ejecta masses or long lightcurve
durations, large ${}^{56}\mathrm{Ni}$ yields, an extremely young surrounding
stellar population, and/or nucleosynthetic signatures might strengthen the
case for associating specific event with this scenario. Possible examples of
SN Ibn that might correspond to PPI+CC are SN2006jc and PS15dpn. The former
showed relatively narrow He lines possibly hinting to asphericity of the CSM
(Foley et al., 2007) and was spatially coincident with an unexplained outburst
two years earlier (e.g., Pastorello et al., 2007; Foley et al., 2007). For the
latter, Wang & Li (2019) proposed to fit the light curve by combining CSM
interaction and radioactive decay, and inferred CSM and ${}^{56}\mathrm{Ni}$
masses of $\sim$$0.8\,M_{\odot}$ and $\sim$ $0.1\,M_{\odot}$, respectively, in
good agreement with our models.

##  9 Limitations and caveats

The stellar evolution simulations presented here require a large number of
assumptions. Work to assess the robustness of these calculations has been
carried out by Marchant et al. (2019); Farmer et al. (2019); Renzo et al.
(2020) (see also Appendix A), to which we refer the readers for more details.

###  9.1 Ignition location and spherical symmetry

Figure 14: Top panel. Temperature and density profiles for example models
approaching the instability, i.e., the first time
$\langle\Gamma_{1}\rangle-4/3=0.01$. The filled circles mark the central
conditions at this stage. The green, yellow, and blue lines show examples of
PPI+CC, PISN, and CC above the mass gap, respectively. All three examples are
labeled according to their $M_{\mathrm{He,init}}$. The red area marks the
region of the EOS where pair-production results in an (local) instability.
Bottom panel. Values of the adiabatic index in the center when approaching the
instability for the entire grid. The thin vertical lines mark the
$M_{\mathrm{CO}}$ of the examples models of the corresponding color in the top
panel.

One of the key assumptions is that spherical symmetry is maintained during the
evolution. Chen et al. (2014); Chen & Woosley (2019) showed that if a pulse
starts symmetrically, hydrodynamic instabilities only weakly deform the pulse.
However, the first stellar layer to become unstable due to pair production in
a pulsating model is not necessarily at the very center, especially at the
lower mass end of the PPI+CC regime. The top panel of Figure 14 shows the
temperature and density profile of three examples with
$M_{\mathrm{He,init}}=50,\,81,\,250\,M_{\odot}$, representative of PPI+CC,
PISN, and CC above the mass gap, respectively. These $M_{\mathrm{He,init}}$
correspond to $M_{\mathrm{CO}}=36.04,\,55.02,\,134.60\,M_{\odot}$,
respectively. The stellar tracks are plotted at the time when the volumetric
pressure-weighted average $\langle\Gamma_{1}\rangle$ first approaches the
instability value 4/3, i.e., $\langle\Gamma_{1}\rangle-4/3=0.01$. The red
shade emphasizes the instability region (neglecting its weak dependence on the
details of the chemical composition), and the text annotations indicate the
physical ingredients that stabilize the structure outside of this region
(Zeldovich & Novikov, 1999; Kippenhahn et al., 2013).

The middle panel of Figure 14 shows the local value of the adiabatic index in
the center $\Gamma_{1,c}$ across our grid, also plotted when each model first
reaches $\langle\Gamma_{1}\rangle-4/3=0.01$. The bottom panel of Figure 14
shows the innermost unstable mass coordinate $M_{u}$, i.e. the innermost
location where $\Gamma_{1}<4/3$. The colors in the middle and bottom panel
have the same meaning as in Figure 2. Overall, as $M_{\mathrm{CO}}$ increases,
the central value of the adiabatic index $\Gamma_{1,c}$ at the beginning of
the instability decreases, and the location of the instability moves inward.
However, this trend is not completely monotonic in the central part of the
PPI+CC regime.

Two example models that ultimately result in a core collapse are shown in the
top panel of Figure 14. The green line corresponds to our
$M_{\mathrm{He,init}}=50\,M_{\odot}$ example for PPI+CC, which shows more
features compared to the other models, because the chemical stratification is
more important in lower mass models. The central region of the $50\,M_{\odot}$
model has a local value of the adiabatic index $\Gamma_{1,c}-4/3>0.01$ when
$\langle\Gamma_{1}\rangle\simeq 4/3$, i.e., the center is stable when the
averaged $\langle\Gamma_{1}\rangle$ approaches instability. The deepest
interior is too dense to become unstable: the contribution of radiation
pressure to the total pressure decreases and the $e^{\pm}$ pairs fill the
available continuum energy levels, raising the Fermi energy
$E^{e^{\pm}}_{\mathrm{Fermi}}$ and consequently the minimum energy photons
need to produce a pair, preventing layers from undergoing the runaway
instability (e.g., Zeldovich & Novikov, 1999). This is true for all our PPI+CC
models: in the bottom panel of Figure 14, all the green points corresponding
to PPI+CC have central values $\Gamma_{1,c}-4/3>0.01$ when
$\langle\Gamma_{1}\rangle-4/3=0.01$.

Since all our models are strongly radiation pressure dominated, they evolve
with most of the mass along the $\Gamma_{1}\simeq 4/3$ locus on the
$(\rho,T)$-plane, and when the instability starts, a nearly homologous
contraction ensues. However, in the densest part of the core, neutrino cooling
dominates over the energy release of the burning. Therefore, the net energy
release starts off-center in our $50\,M_{\odot}$ example, and because of the
assumption of spherical symmetry, the energy release occupies a spherical mass
shell. However, in nature the ignition might not happen simultaneously across
the entire spherical shell, and this could potentially seed an asymmetric
explosion. If asymmetries can build up rapidly during the pair-instability
driven explosion (possibly aided by rotation), they could also lead to orbital
“kicks” when PPI happens in a binary (Marchant et al., 2019).

The location of the unstable layer at the onset of the instability and of the
location where the instability triggers net energy release both move inward
toward the center as $M_{\mathrm{He,init}}$ and $M_{\mathrm{CO}}$ increase.
The least massive models to go PISN are characterized by having their center
close to the instability, i.e., $\Gamma_{1,c}-4/3\leq 0.01$, when
$\langle\Gamma_{1}\rangle-4/3=0.01$ for the first time, as shown by the yellow
filled circles in the bottom panel and the yellow solid line corresponding to
$M_{\mathrm{He,init}}=81\,M_{\odot}$ in the top panel of Figure 14. Our least
massive PISN model has $M_{\mathrm{He,init}}=80.75\,M_{\odot}$ (corresponding
to $M_{\mathrm{CO}}=54.89\,M_{\odot}$) and a central value
$\Gamma_{1,c}-4/3=0.01$ at this evolutionary stage. This model has
$M_{u}\simeq 0.05\,M_{\odot}$ (shown as a dashed horizontal line in the bottom
panel of Figure 14), i.e. the unstable layer extends almost all the way to the
center.

Models forming a BH above the PISN mass gap (cf. blue line in the top panel of
Figure 14 for a $M_{\mathrm{He,init}}=250\,M_{\odot}$ He core) are also
unstable in their very center when $\langle\Gamma_{1}\rangle$ reaches 4/3, but
the ensuing thermonuclear explosion does not cause either pulses or full
disruption. The different outcome is not caused by lack of energy released in
the explosions, but rather by the inefficient use of that energy (Bond et al.,
1984).

To summarize, the temperature, composition, and density profile of the star
when it approaches the instability (i.e., $\langle\Gamma_{1}\rangle-4/3=0.01$
for the first time) are indicative of its future evolution. In particular, the
local value of the adiabatic index in the center $\Gamma_{1,c}$ at this point
can be used to approximately distinguish PISN evolution with no BH remnant (if
also $\Gamma_{1,c}-4/3\leq 0.01$) or PPI+CC evolution with a BH remnant (if
instead the center is safely stable with $\Gamma_{1,c}-4/3>0.01$ when the star
as a whole is becoming unstable). This provides a criterion to decide the
final fate of a star without having to compute the hydrodynamical phase.

###  9.2 CSM structure and composition

We described in Section 6 the amount of mass ejected, its initial velocity,
and the ejection timing resulting from our simulations. We typically keep the
ejecta on our Lagrangian grid for several timesteps after ejection (until
either the bound layers have recovered hydrostatic equilibrium or the onset of
CC is reached). These ejected layers are moving significantly faster than the
escape velocity and the sound speed, and our PPI+CC models exhibit an overall
velocity gradient increasing outward¸ so the ejecta do not cause any back-
reaction on the inner layers that remain bound. Based on the initial mass,
velocity, and time of the ejections, Section 6 illustrates the main features
we expect in the CSM structure surrounding these stars with a toy model
assuming propagation at constant velocity of the ejecta. This is an
oversimplification, since the low density ejecta are likely to be optically
thin and thus can lose energy radiatively. Moreover, as previously noted by
Woosley (2017), the ejected shells can in many cases collide with each other,
and this could also significantly change the CSM structure at the end.
Multidimensional radiation hydrodynamics calculations using our results as
input for the mass, chemical composition, and thermal state of the ejecta
could be used to predict more robustly the CSM structure around PPI+CC models
for comparison with observed transients, and to address the question of how
many progenitors might reach the final CC embedded in a optically-thick layer
of previously ejected material.

Another assumption in our calculations is that the presence of a H-rich
envelope can be neglected to study the dynamics of PPI (Woosley, 2017, 2019)
and that, even if present, such a envelope would be removed early in the
evolution by winds or binary interactions. Should a star retain some H-rich
material until the onset of the first pulse, we can estimate if it would be
detectable in the CSM surrounding these stars assuming that (_i_) the PPI-
driven mass loss timing is unaffected by the presence of the H-rich envelope
and (ii) the entire envelope is ejected in the first pulse (Woosley, 2017).

To estimate the ejection velocities and radii of the H-rich material, we ran a
H-rich $140\,M_{\odot}$ model with initial He abundance $Y=0.27$ and
metallicity $Z=0.001$ with the same setup as in our grids. In Appendix B, we
compare this model to a He core of $55.25\,M_{\odot}$ which produces a
$M_{\mathrm{CO}}\simeq 39.41\,M_{\odot}$, similar in mass to the CO core
produced by the H-rich $140\,M_{\odot}$ star. The H-rich model reaches the
onset of the PPI ($\langle\Gamma_{1}\rangle-4/3=0.01$) with
$\Gamma_{1,c}-4/3=0.05$, so we expect it to follow the PPI+CC evolutionary
path based on Section 9.1. This expectation is confirmed by our results
presented in Appendix B. At the onset of the instability and with our assumed
wind mass loss, this model has a total mass of
$M_{\mathrm{tot}}=83.2\,M_{\odot}$, a He core of $M_{\mathrm{He}}\simeq
50\,M_{\odot}$, and a CO core of $M_{\mathrm{CO}}\simeq 39\,M_{\odot}$. The
remaining envelope has a mass of $M_{\mathrm{env}}\equiv
M_{\mathrm{tot}}-M_{\mathrm{He}}\simeq 33\,M_{\odot}$, however, the
composition of this envelope is dominated by He, with a H mass fraction of
$X\simeq 0.14$, since the winds have carved out material down to the initial
location of the main sequence core of the star. At this stage, the envelope
spans from the He core edge at $R_{\mathrm{He\,core}}=0.58\,R_{\odot}$ to
$R_{*}=3637\,R_{\odot}$, so it is significantly extended.

If we assume propagation at constant velocity of this envelope, we can
estimate the minimum and maximum radii ($R_{\mathrm{min}}$ and
$R_{\mathrm{max}}$, respectively) of this H-rich material at the time of the
final CC and its average density ($\langle\rho\rangle$) with

$\displaystyle\begin{split}R_{\mathrm{min}}=R_{\mathrm{He\,core}}+v_{\mathrm{esc,in}}\times(t_{\mathrm{CC}}-t_{\mathrm{pulse\
end}})\ \,\\\
R_{\mathrm{max}}=R_{*}+v_{\mathrm{esc,out}}\times(t_{\mathrm{CC}}-t_{\mathrm{pulse\
end}})\ \,\\\
\langle\rho\rangle=\frac{3M_{\mathrm{env}}}{4\pi(|R_{\mathrm{max}}^{3}-R_{\mathrm{min}}^{3}|)}\
\ ,\end{split}$ (3)

where $t_{\mathrm{CC}}-t_{\mathrm{pulse\ end}}\simeq 10\,\mathrm{years}$ based
on the $M_{\mathrm{He}}$ and $M_{\mathrm{CO}}$ of this model and Figure 9. If
we assume both escape velocities to be the stellar surface escape velocity,
i.e.,
$v_{\mathrm{esc,out}}=v_{\mathrm{esc,in}}=\sqrt{2GM_{\mathrm{tot}}/R_{*}}$,
then the CSM layer containing H would be at $R_{\mathrm{min}}\simeq
10^{15.46}\,\mathrm{cm}\lesssim r\lesssim
R_{\mathrm{max}}\simeq\,10^{15.51}\,\mathrm{cm}$ and have an average density
of $\langle\rho\rangle\simeq 2.2\cdot 10^{-12}\,\mathrm{g\ cm^{-3}}$.

Instead, if we assume,
$v_{\mathrm{esc,in}}=\sqrt{2GM_{\mathrm{He}}/R_{\mathrm{He\,core}}}$ is the
velocity from the He core edge, while keeping the same $v_{\mathrm{esc,out}}$,
then we obtain $R_{\mathrm{max}}<R_{\mathrm{min}}$ (meaning there would
necessarily be collisions internal to the ejecta that our simplistic toy model
ignores). Nevertheless, the minimum radius at which we would expect CSM with H
in it form the PPI mass ejection is $R_{\mathrm{max}}\simeq 3.2\cdot
10^{16}\,\mathrm{cm}$ and such layer would have an average density of
$\langle\rho\rangle\simeq 3\cdot 10^{-18}\,\mathrm{g\ cm^{-3}}$. The large
difference between these two estimates is mostly due to the large radial
extent of the envelope.

Figure 15: Summary of the pair-instability driven behavior of models as a
function of their initial He core mass ($M_{\mathrm{He,init}}$) and maximum
carbon-oxygen core mass ($M_{\mathrm{CO}}$). The approximate SN type in the
top row speculates on what the resulting SN would appear like, if the final
collapse produces an explosion.

##  10 Conclusions

The broad theoretical understanding of the predicted pair-instability driven
transients has been well established for several decades, however they remain
somewhat elusive from an observational perspective. Recent developments in
stellar evolution calculations allow for the exploration of synergies between
gravitational waves and time-domain observations to better understand the
formation process of the most massive stellar BHs.

We have computed a grid of naked He star models in the mass range
$35\,M_{\odot}\leq M_{\mathrm{He,init}}\lesssim 250\,M_{\odot}$ to investigate
whether these would experience phases of global dynamical instability and
pulsational mass loss due to the pair-production instability. We have computed
grids at two different metallicities, $Z=0.001$ and
$Z=0.00198=14\%\,Z_{\odot}$, although the main features we discuss are not
significantly dependent on $Z$ (except for the wind mass loss rate, see also
Farmer et al. 2019). All our input files and numerical results are available
at https://doi.org/10.5281/zenodo.3406356.

Figure 15 summarizes our main results across the mass range considered. We
find, in agreement with previous studies, that stars enter into the PPI regime
progressively. The production of $e^{\pm}$ initially causes “oscillations” of
the core temperature and nuclear luminosity at the lowest mass end. The least
massive models experiencing an explosive thermonuclear ignition
($M_{\mathrm{He,init}}\gtrsim 37.5\,M_{\odot}$) do not suffer significant
global consequences (“weak pulses”). Increasing further the initial He core
mass, the pulses become progressively stronger, causing at first large radial
expansions (for initial $41\lesssim M_{\mathrm{He,init}}\lesssim
42\,M_{\odot}$), and finally (for initial $M_{\mathrm{He,init}}\gtrsim
42\,M_{\odot}$) also the ejection of matter. The values quoted here are for
the initial He core mass of our models, which can be interpreted as the core
mass at the end of the main sequence of the star. The mapping of these values
to the final (preinstability) He core mass is mass loss and metallicity
dependent (see also Farmer et al. 2019).

The different effects of a pulsational pair instability event on the star
allow for (at least) three different physically-motivated definitions of a
“pulse”, depending on which observable is considered (Section 5). The number
of pulses and which He and CO core masses produce pulses might vary
significantly depending on the observable of interest.

The first definition of pulse we consider (Section 5.1) is based on the core
thermonuclear ignition, following the historical development of studies of
pair-instability evolution. We find that in the lowest mass models the core
ignition does not produce an observable electromagnetic signal or a
significant impact on the final BH mass: the nuclear energy released in the
burning is redistributed and stored in the star without affecting
significantly the outermost layers. The most promising way to detect directly
these core-ignition events is through the variations in the neutrino
luminosity.

The second definition is based on the radial expansion of the models in
response to the core ignition (Section 5.2): this definition shifts the lower
edge of the pulsating regime upward in mass. For the most massive pulsating
models (which also eject mass), the radial expansion itself might be hidden
behind a pseudo-photosphere in the ejected layers. Because of the rarity of
these stars in the local Universe, the most promising way to detect these
radius variations may be through their enhancement of the rate of binary
interactions (e.g., Marchant et al., 2019), although nearby stars that might
undergo this evolution in the future exist (e.g., Bestenlehner et al., 2011;
Crowther et al., 2016; Renzo et al., 2019).

The third definition is based on the ejection of mass (Section 5.3), which
impacts both the circumstellar material around these stars and the final BH
mass they produce. The ejected matter creates shells of ejecta surrounding the
star. If the final collapse results in a successful (even if weak) explosion,
the final SN ejecta can hit this PPI-produced CSM and convert kinetic energy
into radiation.

In our grid, we find full disruption in a PISN for an initial He core mass of
$M_{\mathrm{He,init}}\simeq 80\,M_{\odot}$ corresponding to a final He core
mass of $M_{\mathrm{He}}\simeq 60\,M_{\odot}$ and
$M_{\mathrm{CO}}=55\,M_{\odot}$ after the wind mass loss. With our assumptions
for the wind mass loss and metallicity, most PISN models would still retain
He-rich material at their surface at the onset of the explosion. We propose a
simplified criterion to distinguish full disruption in a PISN from pulsational
behavior producing a final BH based on the adiabatic index at the center of
the star $\Gamma_{1,c}\lesssim 0.01$ at the onset of the instability, defined
as the first moment when the volumetric pressure-weighted average of the
adiabatic index $\langle\Gamma_{1}\rangle-4/3=0.01$ (Section 9.1).This
threshold value allows for the approximately estimate the fate of a stellar
model without the need to compute the hydrodynamical evolution.

The signature on the final BH masses is potentially detectable with a
population of gravitational wave sources (Fishbach & Holz, 2017; Talbot &
Thrane, 2018; Stevenson et al., 2019; Mangiagli et al., 2019). Only for
$M_{\mathrm{He,init}}\gtrsim 42\,M_{\odot}$ ($M_{\mathrm{CO}}\simeq
30.75\,M_{\odot}$) is the stellar core mass significantly reduced and the
circumstellar material significantly affected, as shown in Figure 15. The
maximum BH mass below the PISN mass gap that we find is $\sim$ $45\,M_{\odot}$
and it is formed by the collapse of an initially $M_{\mathrm{He,init}}\simeq
62\,M_{\odot}$ ($M_{\mathrm{CO}}\simeq 43.5\,M_{\odot}$) He core that went
through pulsational mass loss. More massive He cores also produce BHs, but
because of the stronger mass loss due to winds and pair-instability driven
pulses, the resulting BH masses are smaller.

Pair instability does not result in full disruption of the star for an initial
He core mass of $M_{\mathrm{He,init}}\simeq 200\,M_{\odot}$, which forms a BH
of mass $M_{\mathrm{BH}}=125\,M_{\odot}$ after wind mass loss. Above this He
core mass, the photodisintegration of newly synthesized heavy elements during
the thermonuclear explosion prevents the disruption of the entire star (e.g.,
Bond et al., 1984). The boundaries between PISN and BH formation we find are
in very good agreement with previously published results.

We have characterized the CSM properties around the pulsating models resulting
in mass ejections by assuming unperturbed constant velocity propagation of the
ejecta. Under this simplifying assumption, we find that the CSM mass grows
almost monotonically with the initial $M_{\mathrm{He,init}}$ from $\sim$
$10^{-6}\,M_{\odot}$ (for $M_{\mathrm{He,init}}\simeq 42\,M_{\odot}$) to
$\sim$ $20\,M_{\odot}$ at the edge of the PISN range. For initial
$M_{\mathrm{He,init}}\gtrsim 50\,M_{\odot}$, the combined ejection of matter
and mixing during a pulse propagation make He less abundant than C and O at
the stellar surface at the onset of the final core-collapse. The stellar
surface at the onset of core collapse might be obscured by the previously
ejected layers.

The velocity of the ejecta is a few thousand ${\mathrm{km\ s^{-1}}}$, with the
first mass-loss event often resulting in larger velocities. Nevertheless, with
our assumptions, we find numerous self-collisions of the ejecta with
previously ejected layers in agreement with the predictions of Woosley et al.
(2007) and Woosley (2017). This velocity range is close to the width of narrow
He lines detected in some SN Ibn.

The timing of pair-instability driven mass ejections also spans a large range
of values, with a systematic trend of longer delays between pulses for the
more massive models. This is because more massive models produce more
energetic pulses that require a longer time (up to $\sim$ $10^{4}$ years) for
the star to recover its equilibrium.

With the velocity and timing of the ejecta produced by our models, we expect
the PPI-produced circumstellar material to be at $\sim$
$10^{12}-10^{16}\,\mathrm{cm}$ away from the collapsing star at the end of its
evolution. This range covers the distances inferred in observational
candidates for pulsational pair instability evolution.

Gravitational wave detections of merging binary black holes are rapidly
accumulating during the third LIGO/Virgo observing run, and currently
available constraints on their mass will soon become statistically stringent.
Together with the ongoing and upcoming large time-domain survey which will
reveal a plethora of transient, including rare and exotic ones, this will
provide direct constraints on pair-instability evolution of the most massive
stars. Thus, in the future, gravitational and transient observations will soon
shed light on the pair-instability evolution of the most massive stars and the
BHs these produce.

###### Acknowledgements.

We are grateful for the in depth reading of the referee and the suggested
improvements to this manuscript. We acknowledge helpful discussions with M.
Cantiello, D. Hendricks, E. Laplace, I. Mandel, B. Paxton, F. Timmes, L. van
Son, and A. Vigna-Gómez. MR, SJ, and SdM acknowledge funding by the European
Union’s Horizon 2020 research and innovation programme from the European
Research Council (ERC) (Grant agreement No. 715063), and by the Netherlands
Organisation for Scientific Research (NWO) as part of the Vidi research
program BinWaves with project number 639.042.728. RF is supported by the
Netherlands Organisation for Scientific Research (NWO) through a top module 2
grant with project number 614.001.501 (PI de Mink). EZ acknowledges support
from the the Swiss National Science Foundation Professorship grant (project
number PP00P2 176868) and from the Federal Commission for Scholarships for
Foreign Students for the Swiss Gov- ernment Excellence Scholarship (ESKAS No.
2019.0091) for the academic year 2019-2020. Simulations were carried out on
the Dutch national e-infrastructure (Cartesius, project number 16343) with the
support of the SURF Cooperative.

## References

  * Abbott et al. (2019) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Physical Review X, 9, 031040
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Phys. Rev. Lett., 118, 221101
  * Adams et al. (2017) Adams, S. M., Kochanek, C. S., Gerke, J. R., Stanek, K. Z., & Dai, X. 2017, MNRAS, 468, 4968
  * Aoki et al. (2014) Aoki, W., Tominaga, N., Beers, T. C., Honda, S., & Lee, Y. S. 2014, Science, 345, 912
  * Arcavi et al. (2017) Arcavi, I., Hosseinzadeh, G., Howell, D. A., et al. 2017, Nature, 551, 64
  * Arnett (1969) Arnett, W. D. 1969, Ap&SS, 5, 180
  * Atri et al. (2019) Atri, P., Miller-Jones, J. C. A., Bahramian, A., et al. 2019, MNRAS, 489, 3116
  * Barkat et al. (1967) Barkat, Z., Rakavy, G., & Sack, N. 1967, Physical Review Letters, 18, 379
  * Barnes et al. (2018) Barnes, J., Duffell, P. C., Liu, Y., et al. 2018, ApJ, 860, 38
  * Belczynski et al. (2016) Belczynski, K., Heger, A., Gladysz, W., et al. 2016, A&A, 594, A97
  * Bellm (2014) Bellm, E. 2014, in The Third Hot-wiring the Transient Universe Workshop, ed. P. R. Wozniak, M. J. Graham, A. A. Mahabal, & R. Seaman, 27–33
  * Bestenlehner et al. (2011) Bestenlehner, J. M., Vink, J. S., Gräfener, G., et al. 2011, A&A, 530, L14
  * Bond et al. (1984) Bond, J. R., Arnett, W. D., & Carr, B. J. 1984, ApJ, 280, 825
  * Brandt et al. (1995) Brandt, W. N., Podsiadlowski, P., & Sigurdsson, S. 1995, MNRAS, 277, L35
  * Chan et al. (2020) Chan, C., Mueller, B., & Heger, A. 2020, arXiv e-prints, arXiv:2003.04320
  * Chan et al. (2018) Chan, C., Müller, B., Heger, A., Pakmor, R., & Springel, V. 2018, ApJ, 852, L19
  * Chatzopoulos & Wheeler (2012a) Chatzopoulos, E. & Wheeler, J. C. 2012a, ApJ, 748, 42
  * Chatzopoulos & Wheeler (2012b) Chatzopoulos, E. & Wheeler, J. C. 2012b, ApJ, 760, 154
  * Chatzopoulos et al. (2013) Chatzopoulos, E., Wheeler, J. C., & Couch, S. M. 2013, ApJ, 776, 129
  * Chen et al. (2014) Chen, K.-J., Woosley, S., Heger, A., Almgren, A., & Whalen, D. J. 2014, The Astrophysical Journal, 792, 28
  * Chen & Woosley (2019) Chen, K.-J. & Woosley, S. E. 2019, arXiv e-prints [arXiv:1904.12877]
  * Chevalier & Irwin (2012) Chevalier, R. A. & Irwin, C. M. 2012, ApJ, 747, L17
  * Couch et al. (2019) Couch, S. M., Warren, M. L., & O’Connor, E. P. 2019, arXiv e-prints, arXiv:1902.01340
  * Coughlin et al. (2018) Coughlin, E. R., Quataert, E., Fernández, R., & Kasen, D. 2018, MNRAS, 477, 1225
  * Crowther et al. (2016) Crowther, P. A., Caballero-Nieves, S. M., Bostroem, K. A., et al. 2016, MNRAS, 458, 624
  * de Mink et al. (2009) de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243
  * de Mink et al. (2014) de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7
  * Dessart et al. (2012) Dessart, L., Hillier, D. J., Li, C., & Woosley, S. 2012, MNRAS, 424, 2139
  * Ertl et al. (2016) Ertl, T., Janka, H.-T., Woosley, S. E., Sukhbold, T., & Ugliano, M. 2016, ApJ, 818, 124
  * Farmer et al. (2016) Farmer, R., Fields, C. E., Petermann, I., et al. 2016, ApJS, 227, 22
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53
  * Farr et al. (2011) Farr, W. M., Sravan, N., Cantrell, A., et al. 2011, ApJ, 741, 103
  * Fishbach & Holz (2017) Fishbach, M. & Holz, D. E. 2017, ApJ, 851, L25
  * Foley et al. (2007) Foley, R. J., Smith, N., Ganeshalingam, M., et al. 2007, ApJ, 657, L105
  * Fowler & Hoyle (1964) Fowler, W. A. & Hoyle, F. 1964, ApJS, 9, 201
  * Fragos et al. (2009) Fragos, T., Willems, B., Kalogera, V., et al. 2009, ApJ, 697, 1057
  * Fraley (1968) Fraley, G. S. 1968, Ap&SS, 2, 96
  * Fryer et al. (2001) Fryer, C. L., Woosley, S. E., & Heger, A. 2001, ApJ, 550, 372
  * Gal-Yam et al. (2009) Gal-Yam, A., Mazzali, P., Ofek, E. O., et al. 2009, Nature, 462, 624
  * Georgy et al. (2017) Georgy, C., Meynet, G., Ekström, S., et al. 2017, A&A, 599, L5
  * Gilkis et al. (2016) Gilkis, A., Soker, N., & Papish, O. 2016, ApJ, 826, 178
  * Glatzel et al. (1985) Glatzel, W., Fricke, K. J., & El Eid, M. F. 1985, A&A, 149, 413
  * Gomez et al. (2019) Gomez, S., Berger, E., Nicholl, M., et al. 2019, ApJ, 881, 87
  * Hamann & Koesterke (1998) Hamann, W. R. & Koesterke, L. 1998, A&A, 335, 1003
  * Hamann et al. (1995) Hamann, W.-R., Koesterke, L., & Wessolowski, U. 1995, A&A, 299, 151
  * Heger et al. (2003) Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288
  * Heger & Woosley (2002) Heger, A. & Woosley, S. E. 2002, ApJ, 567, 532
  * Hosseinzadeh et al. (2017) Hosseinzadeh, G., Arcavi, I., Valenti, S., et al. 2017, ApJ, 836, 158
  * Janka (2013) Janka, H.-T. 2013, MNRAS, 434, 1355
  * Janka (2017) Janka, H.-T. 2017, ApJ, 837, 84
  * Kasen et al. (2011) Kasen, D., Woosley, S. E., & Heger, A. 2011, ApJ, 734, 102
  * Kippenhahn & Weigert (1967) Kippenhahn, R. & Weigert, A. 1967, ZAp, 65, 251
  * Kippenhahn et al. (2013) Kippenhahn, R., Weigert, A., & Weiss, A. 2013, Stellar Structure and Evolution (Springer-Verlag)
  * Kozyreva et al. (2018) Kozyreva, A., Kromer, M., Noebauer, U. M., & Hirschi, R. 2018, MNRAS, 479, 3106
  * Kuroda et al. (2018) Kuroda, T., Kotake, K., Takiwaki, T., & Thielemann, F.-K. 2018, MNRAS, 477, L80
  * Langer et al. (2007) Langer, N., Norman, C. A., de Koter, A., et al. 2007, A&A, 475, L19
  * Leung et al. (2019) Leung, S.-C., Nomoto, K., & Blinnikov, S. 2019, ApJ, 887, 72
  * Lovegrove & Woosley (2013) Lovegrove, E. & Woosley, S. E. 2013, ApJ, 769, 109
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints [arXiv:0912.0201]
  * Lunnan et al. (2018) Lunnan, R., Chornock, R., Berger, E., et al. 2018, ApJ, 852, 81
  * Maeder & Meynet (2000) Maeder, A. & Meynet, G. 2000, ARA&A, 38, 143
  * Mandel & de Mink (2016) Mandel, I. & de Mink, S. E. 2016, MNRAS, 458, 2634
  * Mangiagli et al. (2019) Mangiagli, A., Bonetti, M., Sesana, A., & Colpi, M. 2019, arXiv e-prints, arXiv:1907.12562
  * Marchant et al. (2016) Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50
  * Marchant et al. (2019) Marchant, P., Renzo, M., Farmer, R., et al. 2019, ApJ, 882, 36
  * Moriya & Langer (2015) Moriya, T. J. & Langer, N. 2015, A&A, 573, A18
  * Mösta et al. (2015) Mösta, P., Ott, C. D., Radice, D., et al. 2015, Nature, 528, 376
  * Müller (2019) Müller, B. 2019, MNRAS, 1527
  * Nadezhin (1980) Nadezhin, D. K. 1980, Ap&SS, 69, 115
  * Nieuwenhuijzen & de Jager (1990) Nieuwenhuijzen, H. & de Jager, C. 1990, A&A, 231, 134
  * Nugis & Lamers (2000) Nugis, T. & Lamers, H. J. G. L. M. 2000, A&A, 360, 227
  * O’Connor & Ott (2011) O’Connor, E. & Ott, C. D. 2011, ApJ, 730, 70
  * O’Shaughnessy et al. (2017) O’Shaughnessy, R., Gerosa, D., & Wysocki, D. 2017, arXiv:1704.03879 [arXiv:1704.03879]
  * Ott et al. (2018) Ott, C. D., Roberts, L. F., da Silva Schneider, A., et al. 2018, ApJ, 855, L3
  * Pastorello et al. (2008) Pastorello, A., Mattila, S., Zampieri, L., et al. 2008, MNRAS, 389, 113
  * Pastorello et al. (2007) Pastorello, A., Smartt, S. J., Mattila, S., et al. 2007, Nature, 447, 829
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Rakavy & Shaviv (1967) Rakavy, G. & Shaviv, G. 1967, ApJ, 148, 803
  * Regos et al. (2020) Regos, E., Vinko, J., & Ziegler, B. L. 2020, arXiv e-prints, arXiv:2002.07854
  * Renzo et al. (2019) Renzo, M., de Mink, S. E., Lennon, D. J., et al. 2019, MNRAS, 482, L102
  * Renzo et al. (2020) Renzo, M., Farmer, R. J., Justham, S., et al. 2020, MNRAS, 493, 4333
  * Renzo et al. (2017) Renzo, M., Ott, C. D., Shore, S. N., & de Mink, S. E. 2017, A&A, 603, A118
  * Soker (2019) Soker, N. 2019, Research in Astronomy and Astrophysics, 19, 095
  * Spera & Mapelli (2017) Spera, M. & Mapelli, M. 2017, MNRAS, 470, 4739
  * Spera et al. (2019) Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889
  * Stevenson et al. (2019) Stevenson, S., Sampson, M., Powell, J., et al. 2019, ApJ, 882, 121
  * Stothers (1999) Stothers, R. B. 1999, MNRAS, 305, 365
  * Takahashi (2018) Takahashi, K. 2018, ApJ, 863, 153
  * Talbot & Thrane (2018) Talbot, C. & Thrane, E. 2018, ApJ, 856, 173
  * Tolstov et al. (2017) Tolstov, A., Nomoto, K., Blinnikov, S., et al. 2017, ApJ, 835, 266
  * Toro et al. (1994) Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25
  * Tramper et al. (2016) Tramper, F., Sana, H., & de Koter, A. 2016, ApJ, 833, 133
  * Ugliano et al. (2012) Ugliano, M., Janka, H.-T., Marek, A., & Arcones, A. 2012, ApJ, 757, 69
  * van den Heuvel & Portegies Zwart (2013) van den Heuvel, E. P. J. & Portegies Zwart, S. F. 2013, ApJ, 779, 114
  * Vigna-Gómez et al. (2019) Vigna-Gómez, A., Justham, S., Mandel, I., de Mink, S. E., & Podsiadlowski, P. 2019, ApJ, 876, L29
  * Vink & de Koter (2005) Vink, J. S. & de Koter, A. 2005, A&A, 442, 587
  * Vink et al. (2000) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
  * Wang & Li (2019) Wang, S.-Q. & Li, L. 2019, arXiv e-prints, arXiv:1905.12623
  * Whalen et al. (2013) Whalen, D. J., Even, W., Frey, L. H., et al. 2013, ApJ, 777, 110
  * Wongwathanarat et al. (2013) Wongwathanarat, A., Janka, H.-T., & Müller, E. 2013, A&A, 552, A126
  * Woosley (2017) Woosley, S. E. 2017, ApJ, 836, 244
  * Woosley (2018) Woosley, S. E. 2018, ApJ, 863, 105
  * Woosley (2019) Woosley, S. E. 2019, ApJ, 878, 49
  * Woosley et al. (2007) Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390
  * Woosley et al. (2002) Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Mod. Phys., 74, 1015
  * Wyrzykowski et al. (2016) Wyrzykowski, Ł., Kostrzewa-Rutkowska, Z., Skowron, J., et al. 2016, MNRAS, 458, 3012
  * Wyrzykowski & Mandel (2019) Wyrzykowski, Ł. & Mandel, I. 2019, arXiv e-prints, arXiv:1904.07789
  * Yoon et al. (2017) Yoon, S.-C., Dessart, L., & Clocchiatti, A. 2017, ApJ, 840, 10
  * Yoon et al. (2012) Yoon, S.-C., Dierks, A., & Langer, N. 2012, A&A, 542, A113
  * Yoon et al. (2006) Yoon, S.-C., Langer, N., & Norman, C. 2006, A&A, 460, 199
  * Yoshida et al. (2016) Yoshida, T., Umeda, H., Maeda, K., & Ishii, T. 2016, MNRAS, 457, 351
  * Zackay et al. (2019) Zackay, B., Venumadhav, T., Dai, L., Roulet, J., & Zaldarriaga, M. 2019, arXiv e-prints, arXiv:1902.10331
  * Zeldovich & Novikov (1999) Zeldovich, Y. B. & Novikov, I. D. 1999, Stars and relativity

##  Appendix A Resolution study

Figure 16: Mass (blue, right axis) and central temperature (red, left axis)
evolution for our example $50\,M_{\odot}$ example He core. The scale on the
time axis emphasizes the short lived final phases. The differences in the
final mass due to variations in the spatial and temporal mesh arise mostly
from the dynamical phase of evolution and are limited to $\Delta M\lesssim
0.5\,M_{\odot}$.

We present here a study of the impact of the numerical resolution in
Lagrangian mass coordinate and time on our results. We refer the interested
readers to Marchant et al. (2019) for a study on the relaxation procedure and
of the numerical resolution on the pulse mass loss, to Farmer et al. (2019)
for a more comprehensive study of the impact of the numerical resolution and
input physics variations on our PPI models, and Renzo et al. (2020) for a
study of the impact of different treatments of time-dependent convection.

Figure 16 shows the evolution in time of the amount of mass bound to the star
(blue) and its core temperature (red) for two $50.0\,M_{\odot}$ He core models
computed with different resolutions. MESA offers many controls to fine-tune
the resolution (see also the provided inlists and run_star_extras.f), here we
vary only three parameters governing the overall variations of averaged
quantities across adjacent mesh points and across timesteps. Our fiducial
(higher) resolution uses mesh_delta_coeff=0.8 (0.6) and
split_merge_amr_nz_baseline=8 000 (10 000) for the spatial resolution during
the hydrostatic and hydrodynamical phases of evolution, respectively. The time
discretization is controlled through varcontrol_target=5d-5 (1d-5). The
largest differences in the evolution are found after the onset of the PPI
pulses, during the dynamical phase after
$\log_{10}\\{(t_{\mathrm{CC}}-t)/\mathrm{[yr]}\\}\lesssim-2$. These result in
a $\Delta M=0.43\,M_{\odot}$ difference in the mass remaining bound (and on
the amount of mass ejected). We emphasize that even our fiducial value
provides a resolution significantly higher than the MESA defaults, with a
number of mesh points $1289\lesssim N\lesssim 6311$ and 87900 timesteps from
the onset of He core burning to the onset of core-collapse. For comparison,
the higher-resolution model shown in Figure 16 has $1583\lesssim N\lesssim
7908$, however it is able to finish the evolution using a slightly smaller
number of timesteps, 86165. This likely indicates that at the higher spatial
resolution the most stringent condition on the timesteps is not
varcontrol_target, and that the higher spatial resolution provides more
numerical stability of the solution allowing for longer timesteps.

##  Appendix B Comparison to H-rich model

A full characterization of PPI+CC evolution for stars with a H-rich envelope
is beyond the scope of this study (we refer interested readers to Woosley
2017; Leung et al. 2019). Here, we present a brief comparison between a
$140\,M_{\odot}$ H-rich model (computed assuming and initial He abundance of
$Y=0.27$) and a $M_{\mathrm{He,init}}=55.25\,M_{\odot}$. Both these models
produce $M_{\mathrm{CO}}\simeq 40\,M_{\odot}$.

Figure 17: Comparison of the evolution of the nuclear luminosity (top) and
neutrino luminosity (center) integrated throughout the star, and central
temperature (bottom) for a H-rich initially $140\,M_{\odot}$ star (red solid
lines) to a $M_{\mathrm{He,init}}=55.25\,M_{\odot}$ core (blue dashed lines).
Both models form a $M_{\mathrm{CO}}\simeq 40\,M_{\odot}$.

Figure 17 shows the evolution of the nuclear and neutrino luminosities
(integrated throughout the entire star) and the central temperature of these
two models. Even if there are some differences, the rough timing and amplitude
of the oscillations of the core quantities are not significantly different
during the pulses, i.e. during the last year of evolution. This is expected,
since in general for all evolved massive stars the core evolution is not
determined by the envelope. Instead, the evolution is driven by the neutrino
losses from the core itself, rather than the photon luminosity at the surface
(e.g., Fraley 1968).

Figure 18: Comparison of the evolution of the total mass (blue), He core mass
(red) and CO core mass (purple) for a $140\,M_{\odot}$ H-rich model (solid)
and a $M_{\mathrm{He,init}}=55.25\,M_{\odot}$ He-core (dashed).

Figure 18 shows the evolution of the total, He and CO core masses for these
two models. For the H-rich model, we define the He core as the outermost
location where the H mass fraction is below 0.5 and the He mass fraction
exceeds 0.1. If this yields zero, we set the He core mass equal to the total
core mass. To avoid noise in the curve caused by mixing episodes at the core
boundary at the end of the main sequence, we only show the He core mass once
its value has stabilized. The H-rich model expells all its envelope at the
first pulse about one year before CC. At this point the He core and the total
mass become the same.

As shown in Table C, our $55.25\,M_{\odot}$ model has three distinct mass-loss
events. The comparison H-rich model also exhibits three mass-loss pulses.
However, the timing in between the events is slightly different in the two
models. This is partly due to the extreme sensitivity to the core masses (or,
more precisely, the entropy profile) at the onset of the pulses. As well as
the presence of the H-rich envelope, especially if extended, which increases
the dynamical timescale and the delay between the core thermonuclear ignition
and shock breakout and thus the mass ejection.

These two models, characterized by similar $M_{\mathrm{CO}}$ (defined in
Section 2 as the maximum throughout the evolution of the CO core mass) yield
similar BH masses of $\simeq 43.5\,M_{\odot}$ for the He core and $\simeq
45\,M_{\odot}$ for the H-rich model. The difference between the BH mass values
found here are smaller than those introduced by other physical and numerical
uncertainties (Farmer et al. 2019).

##  Appendix C Mass loss per pulse

Table C provides data for each individual mass ejection event in our grid of
models. We list the total mass at the beginning of the pulse
$M_{\mathrm{tot}}^{\mathrm{pre-pulse}}$, the amount of mass lost $\Delta
M_{\mathrm{pulse}}$, the delay time between core-collapse and the end of the
pulse, and the pulse duration $\Delta t_{\mathrm{pulse}}$, respectively (as
defined in Section 5.3), and the center-of-mass velocity of the ejected layer
$\langle v\rangle$. This table is also available at
https://doi.org/10.5281/zenodo.3406356.

$M_{\mathrm{He,init}}$ | pulse | $M_{\mathrm{tot}}^{\mathrm{pre-pulse}}$ | $\Delta M_{\mathrm{pulse}}$ | $\log_{10}(t_{\mathrm{CC}}-t_{\mathrm{pulse\ end}})$ | $\Delta t_{\mathrm{pulse}}$ | $\langle v\rangle$
---|---|---|---|---|---|---
$[M_{\odot}]$ |  | $[M_{\odot}]$ | $[M_{\odot}]$ | [yr] | [hours] | $\mathrm{[10^{3}km\ s^{1}]}$
41.25 | 1 | 34.60 | 0.00 | -6.14 | 0 | 3.83
42.00 | 1 | 35.14 | 0.00 | -4.79 | 0 | 2.84
42.25 | 1 | 35.33 | 0.00 | -4.21 | 1 | 2.45
42.50 | 1 | 35.51 | 0.01 | -3.84 | 1 | 1.97
42.75 | 1 | 35.69 | 0.02 | -3.75 | 2 | 1.87
43.00 | 1 | 35.87 | 0.01 | -3.74 | 2 | 1.85
43.25 | 1 | 36.05 | 0.02 | -3.59 | 2 | 1.77
43.50 | 1 | 36.23 | 1.52 | -1.92 | 105 | 3.87
43.75 | 1 | 36.42 | 0.55 | -2.32 | 33 | 3.52
44.00 | 1 | 36.59 | 0.87 | -2.06 | 73 | 3.63
44.25 | 1 | 36.77 | 0.54 | -2.21 | 53 | 3.27
44.50 | 1 | 36.95 | 0.05 | -3.37 | 4 | 1.72
44.75 | 1 | 37.13 | 0.06 | -3.15 | 6 | 1.64
45.00 | 1 | 37.31 | 0.07 | -3.08 | 7 | 1.64
45.25 | 1 | 37.50 | 0.07 | -3.01 | 8 | 1.61
45.50 | 1 | 37.67 | 0.09 | -3.01 | 9 | 1.57
45.75 | 1 | 37.85 | 0.05 | -2.89 | 11 | 1.46
46.00 | 1 | 38.03 | 0.06 | -2.85 | 12 | 1.55
46.25 | 1 | 38.21 | 0.03 | -2.81 | 11 | 1.37
46.50 | 1 | 38.39 | 0.03 | -2.80 | 10 | 1.40
46.75 | 1 | 38.56 | 0.13 | -2.74 | 16 | 1.59
47.00 | 1 | 38.74 | 0.09 | -2.72 | 17 | 1.54
47.25 | 1 | 38.92 | 0.29 | -2.62 | 21 | 1.66
47.50 | 1 | 39.10 | 0.44 | -2.46 | 30 | 1.64
47.75 | 1 | 39.28 | 0.03 | -2.73 | 15 | 1.36
48.00 | 1 | 39.45 | 0.61 | -2.23 | 51 | 1.78
48.25 | 1 | 39.64 | 0.48 | -2.36 | 38 | 1.56
48.75 | 1 | 39.98 | 0.89 | -2.06 | 76 | 1.61
49.00 | 1 | 40.17 | 1.08 | -1.76 | 152 | 1.63
49.25 | 1 | 40.34 | 1.14 | -1.74 | 160 | 1.62
49.50 | 1 | 40.52 | 0.52 | -2.30 | 43 | 1.78
49.75 | 1 | 40.70 | 1.64 | -1.69 | 180 | 1.66
50.00 | 1 | 40.88 | 1.24 | -1.84 | 127 | 1.67
50.25 | 1 | 41.05 | 1.84 | -1.63 | 208 | 1.70
50.50 | 1 | 41.23 | 1.69 | -1.78 | 145 | 1.78
50.75 | 1 | 41.41 | 1.97 | -1.62 | 213 | 1.75
51.00 | 1 | 41.58 | 0.40 | -1.27 | 30 | 2.21
2 | 41.18 | 1.61 | -2.74 | 16 | 1.45
51.50 | 1 | 41.93 | 1.85 | -1.58 | 231 | 1.88
51.75 | 1 | 42.11 | 0.43 | -1.25 | 35 | 2.11
2 | 41.68 | 1.70 | -2.39 | 36 | 1.64
52.00 | 1 | 42.28 | 0.74 | -1.10 | 65 | 2.22
2 | 41.53 | 1.49 | -2.78 | 15 | 1.51
52.25 | 1 | 42.46 | 0.15 | -1.90 | 110 | 1.92
52.50 | 1 | 42.64 | 0.09 | -1.83 | 131 | 1.60
52.75 | 1 | 42.81 | 0.39 | -1.71 | 169 | 1.86
53.00 | 1 | 42.98 | 0.24 | -1.86 | 120 | 1.51
53.25 | 1 | 43.16 | 0.12 | -2.13 | 65 | 1.59
53.50 | 1 | 43.33 | 1.41 | -1.50 | 280 | 1.99
53.75 | 1 | 43.50 | 1.23 | -1.56 | 240 | 1.95
54.00 | 1 | 43.68 | 1.56 | -1.32 | 417 | 1.81
54.50 | 1 | 44.03 | 1.75 | -1.43 | 328 | 1.95
54.75 | 1 | 44.20 | 1.29 | -1.23 | 203 | 2.12
2 | 42.90 | 0.56 | -2.57 | 24 | 1.51
55.00 | 1 | 44.37 | 1.19 | -1.37 | 372 | 1.97
55.25 | 1 | 44.55 | 0.32 | -1.24 | 3 | 2.19
2 | 44.23 | 0.78 | -1.24 | 282 | 1.85
3 | 43.45 | 0.25 | -2.26 | 48 | 1.49
55.50 | 1 | 44.72 | 1.32 | -1.40 | 351 | 1.88
55.75 | 1 | 44.89 | 0.44 | -1.21 | 4 | 2.14
2 | 44.45 | 1.06 | -1.76 | 152 | 1.68
56.00 | 1 | 45.07 | 0.50 | -1.06 | 3 | 2.16
2 | 44.57 | 1.10 | -1.58 | 231 | 1.76