The role of the hadron-quark phase transition in core-collapse supernovae

ABSTRACT

The hadron-quark phase transition in quantum chromodynamics has been suggested as an alternative explosion mechanism for core-collapse supernovae. We study the impact of three different hadron-quark equations of state (EoS) with first-order (DD2F_SF, STOS-B145) and second-order (CMF) phase transitions on supernova dynamics by performing 97 simulations for solar- and zero-metallicity progenitors in the range of |$14\tt {-}100\, \text{M}_\odot$|⁠. We find explosions only for two low-compactness models (14 and |$16\, \text{M}_\odot$|⁠) with the DD2F_SF EoS, both with low explosion energies of |${\sim }10^{50}\, \mathrm{erg}$|⁠. These weak explosions are characterized by a neutrino signal with several minibursts in the explosion phase due to complex reverse shock dynamics, in addition to the typical second neutrino burst for phase-transition-driven explosions. The nucleosynthesis shows significant overproduction of nuclei such as ⁹⁰Zr for the |$14\hbox{-} \text{M}_\odot$| zero-metallicity model and ⁹⁴Zr for the |$16\hbox{-}\text{M}_\odot$| solar-metallicity model, but the overproduction factors are not large enough to place constraints on the occurrence of such explosions. Several other low-compactness models using the DD2F_SF EoS and two high-compactness models using the STOS EoS end up as failed explosions and emit a second neutrino burst. For the CMF EoS, the phase transition never leads to a second bounce and explosion. For all three EoS, inverted convection occurs deep in the core of the protocompact star due to anomalous behaviour of thermodynamic derivatives in the mixed phase, which heats the core to entropies up to 4k_B/baryon and may have a distinctive gravitational-wave signature, also for a second-order phase transition.

equation of state, hydrodynamics, transients: supernovae

1 INTRODUCTION

The iron core of massive stars will undergo gravitational collapse once its mass exceeds the effective Chandrasekhar limit. The minimum mass for stars to burn their cores up to iron is estimated to be about |$8\, {\rm M}_\odot$| (Woosley, Heger & Weaver 2002; Ibeling & Heger 2013). The reduction of degeneracy pressure by electron captures and photodisintegration of heavy nuclei eventually trigger runaway collapse on a free-fall time-scale. The core density increases up to the point where repulsive short-range nuclear interactions come into play, leading to a stiffening of the equation of state and a (first) core bounce, leaving a protoneutron star – or more generally a ‘protocompact star’ (PCS) if non-nucleonic particles eventually appear – at the centre of the star. The rebound of the core launches a shock wave that quickly stalls as its initial kinetic energy is drained by photodissociation of heavy nuclei and neutrino losses as it propagates through the outer core. The subsequent evolution of the shock wave and the high-density EoS determine whether the star ends its life as a compact remnant in form of either a black hole (BH), neutron star (NS), quark star or some other compact star containing non-nucleonic matter. Several mechanisms may revive the shock and result in a core-collapse supernova (CCSN) explosion. The two best-explored scenarios are neutrino-driven explosions, which probably account for most explosions of up to |${\sim }10^{51}\, \mathrm{erg}$|⁠, and the magnetorotational mechanism, which has been proposed to explain exceptionally energetic ‘hypernovae’ with explosion energies of up to |${\sim }10^{52}\, \mathrm{erg}$| (for reviews, see Janka 2012; Burrows 2013; Müller 2020). Both of these mechanisms critically rely on multidimensional effects such as rotation, convection, and other hydrodynamic instabilities. A third proposed mechanism is the phase-transition driven (PT-driven) mechanism (Migdal, Chernoutsan & Mishustin 1979), in particular, a transition from hadrons to quarks (Sagert et al. 2009; Fischer et al. 2018).

By and large, supernova matter requires a conformable treatment for numerous intrinsically different thermodynamic regimes. At temperatures below ≈6 GK, the baryon EoS has to account for heavy nuclei and their abundances, which are usually not in equilibrium with each other and need to be calculated by a nuclear reaction network. At higher temperatures and later stages of the evolution, nuclear statistical equilibrium (NSE) can be applied. From the collapse phase onward, nuclear interactions need to be modelled explicitly at high densities instead of assuming non-interacting nucleons and nuclei. While most current CCSN simulations consider purely nucleonic high-density EoSs (with popular choices including Shen et al. 1998b, a; Lattimer & Swesty 1991; Hempel & Schaffner-Bielich 2010; Hempel et al. 2012; Steiner, Hempel & Fischer 2013), the state of matter is highly uncertain at temperatures of several |$10\, \mathrm{GK}$| and densities well exceeding the nuclear saturation density n₀. It is possible that quarks exist in supernova matter (Gentile et al. 1993; Bednarek, Biesiada & Manka 1996; Drago & Tambini 1999); this is not considered in most standard CCSN models.

The relevant degrees of freedom in the high-density and high-temperature phase of the QCD phase diagram are free quarks and gluons (Witten 1984). Perturbative methods in this regime give expressions for the thermodynamical potential Ω in powers of the QCD strong coupling constant α_s. The latter is density- and temperature-dependent, and becomes small when the quark chemical potential is high or the temperature is large compared to the QCD energy scale Λ_QCD, which corresponds to large chemical potentials of several GeV (Baym et al. 2018a). This is called asymptotic freedom; quarks and gluons in this regime barely interact and can move freely (Peskin & Schroeder 1995). This regime is beyond the densities reached in CCSN. Regions of the QCD phase diagram with lower temperatures and chemical potentials prohibit the use of perturbative methods since the running coupling α_s becomes large and the series expansion of the thermodynamical potential Ω does not converge. Lattice QCD (lQCD) calculations provide insight into matter at near-zero densities and finite temperature (Gattringer & Lang 2010a, b). The solution for QCD predicts a smooth crossover transition near vanishing baryon density at a pseudo-critical temperature in the range |$150\tt {-}160\, \mathrm{MeV}$| (Bazavov et al. 2014, 2019; Borsányi et al. 2014). Several decades of high-energy heavy-ion collision analysis further advanced our understanding on QCD matter (Stöcker & Greiner 1986; Gazdzicki & Gorenstein 1999; Csernai, Magas & Wang 2013; Vovchenko et al. 2020). Heavy-ion collision data indicate a transition from bound hadrons to deconfined quarks at non-vanishing baryon density (Adams et al. 2005; Arsene et al. 2005; Adamczyk et al. 2017). The order of the phase transition at finite densities is, however, not yet fully understood (Annala et al. 2020; Cuteri, Philipsen & Sciarra 2021). With the advancement of multimessenger astronomy (Aartsen et al. 2017; Abbott et al. 2017), astrophysical observations have the potential to further constrain the EoS in the QCD regime as the physical conditions in relativistic heavy-ion collisions are closely linked, e.g. to those in NS mergers (Hades Collaboration et al. 2019; Hanauske et al. 2017, 2019; Most et al. 2022). The extreme environment at the centre of PCSs can exceed densities of |$\rho _\text{centre} \gtrsim 10^{15}\, \mathrm{g}\, \mathrm{cm}^{-3}$| with maximum temperatures of the order |$10^{12}\, \mathrm{K}$|⁠. Phenomenological models are the basis on which predictions about the phase transition from hadrons to quarks rely in these intermediate regimes of the QCD phase diagram. An illustration of the different regimes is shown in Fig. 1.

Figure 1.

Illustrative QCD phase diagram. The x-, y-, and z-axes represent density, temperature, and isospin asymmetry density n_n−n_p. The regime of late post-bounce phase of CCSN and of neutron–star mergers are indicated by the orange band and grey shade.

The hadron-quark phase transition usually softens the EoS during the mixed phase, which lowers the maximal supported mass M_max of the PCS. When the PCS exceeds M_max, it collapses and can perform damped oscillations around its new equilibrium position. If rapid BH formation does immediately follow, the core can bounce a second time, which can, depending on the released binding energy, lead to the formation of a shock wave, which then again can lead to expelling outer layers of the star, leaving behind a stable compact star (Sagert et al. 2009; Fischer et al. 2018).

It is assumed that less massive progenitors with zero-age main sequence (ZAMS) mass |${\lesssim } 20$| M_⊙ explode by the standard NDE-mechanism (Nakamura et al. 2015; Müller et al. 2016; Sukhbold et al. 2016; Burrows et al. 2020), leaving behind NSs, while more massive stars usually collapse into BHs, even if shock revival may sometimes occur and be followed by a weak fallback explosion (e.g. Chan et al. 2018; Kuroda et al. 2018; Ott et al. 2018; Chan, Müller & Heger 2020; Powell, Müller & Heger 2021; Rahman et al. 2021). So far, three-dimensional models of neutrino-driven explosions can only account for CCSNe of normal energies, not significantly exceeding |$10^{51}\, \mathrm{erg}$| (Burrows et al. 2020; Powell & Müller 2020; Bollig et al. 2021). A different mechanism for the most energetic observed CCSNe is probably needed. While the magnetorotational mechanism has long been investigated as an explanation for the most powerful explosions (for recent results on explosion energies, see Kuroda et al. 2020; Obergaulinger & Aloy 2021; Jardine, Powell & Müller 2022), recent studies discussed PT-driven explosions of very massive progenitors as a scenario for various unusually energetic events. For a |$50{\hbox{-}\rm M}_\odot$| progenitor, Fischer et al. (2018) obtained an explosion energy of |$E_\text{exp} = 3\times 10^{51}\, \mathrm{erg}$| in spherical symmetry. The possibility of PT-driven explosions is also relevant for nucleosynthesis. Successful PT-driven explosions have been proposed a site for heavy r-process elements (Nishimura et al. 2012; Fischer et al. 2020). Furthermore, a Galactic PT-supernova would be a promising target for multimessenger observations in neutrinos and gravitational waves. The neutrino signal would provide a characteristic fingerprint for the QCD PT in the form of an electron antineutrino burst (Sagert et al. 2009; Fischer et al. 2018; Zha et al. 2020) that would clearly be observable by present and future detectors such as IceCube, Super-Kamiokande, and Hyper-K. A first-order QCD phase transition is also expected to produce a strong and characteristic gravitational-wave signal peaking at several kHz, regardless of whether a successful explosion ensues (Zha et al. 2020; Kuroda et al. 2021; see also Yasutake et al. 2007 for other effects on the gravitational-wave signal).

The robustness of the PT-driven mechanism, however, is far from clear yet. At this stage, it is important to more systematically scan the parameter space for PT-driven explosions using larger sets of progenitors (⁠|${\gg } 20$|⁠) than in the currently available literature. Furthermore, only first-order phase transitions were considered so far. Recent studies motivate further investigation into the robustness of the PT-driven scenario. Zha, O’Connor & da Silva Schneider (2021) found no successful explosions in 1D simulations using the STOS-B145 EoS, although they found instances of second bounces to a more compact and (transiently) stable PCS, with a strong dependence of the dynamics on the compactness parameter. Using the DD2F_SF EoS, Fischer (2021) recently found explosions only at solar metallicity, but not at low metallicity for two |$75\hbox{-}{\rm M}_\odot$| models.

This paper aims to shed more light on the progenitor dependence of the post-bounce evolution of hybrid PCSs containing quark matter. We perform 97 general-relativistic hydrodynamic simulations with neutrino transport in spherical symmetry for up to 40 progenitors in the mass range |$M\in [14,100]\, {\rm M}_\odot$| with solar and zero metallicity using three different hybrid EoS. We especially focus on the thermodynamic features of the mixed phase and crossover regions and how these affect the post-bounce dynamics once the PCS reaches the threshold density for the appearance of quarks. We find that only two models using the DD2F_SF-1.4 EoS explode. For those two cases, we perform a detailed nucleosynthesis analysis. We also study the neutrino signals of exploding and non-exploding models confirming a similar phenomenology as found by Zha et al. (2020).

This work is structured as follows. In Section 2, we discuss the set of equations of state (EoS) and give an overview of numerical methods including progenitor setup and the nucleosynthesis post-processing. In Section 3, we present the results of our simulations. We interpret the progenitor dependence based on the DD2F_SF EOS and discuss detailed hydrodynamic post-bounce dynamics for two exploding models in the DD2F_SF setup. We analyse the effect of phase transitions by comparing hydrodynamic simulation outcomes for different EoS, including the neutrino signals. Lastly, we review nucleosynthesis results in the ejecta for two exploding models. We summarize and our findings and their implications in Section 4.

2 SIMULATION SETUP

2.1 EoS with quark matter

We study the effect of a quark-hadron phase transition in CCSNe using three high-density EoS with different treatments for the hadronic phase, the quark phase, and the phase transition between them.

2.1.1 DD2-RMF EoS

The DD2F_SF EoS is a hadron-quark EoS featuring a first-order phase transition to deconfined quark matter. The EoS belongs to a new class of hybrid EoS, using a relativistic density-functional formalism (Fischer et al. 2018; Bastian 2021), originally adopted from a relativistic mean-field theory with density-dependent meson-nucleon coupling constants (DD2; Typel & Wolter 1999; Hempel & Schaffner-Bielich 2010; Typel et al. 2010) with a string-flip microscopic quark-matter model (Kaltenborn, Bastian & Blaschke 2017). Repulsive higher order quark-quark interactions give rise to additional pressure contributions with increasing densities (Klähn & Fischer 2015). A vector interaction potential in the quark phase supports high maximum masses for NSs (Kaltenborn et al. 2017) and twin stars (Benic et al. 2015). The phase transition between hadronic and quark matter is modelled via the Gibbs construction (global charge neutrality in the mixed phase).

2.1.2 STOS-B145 EoS

As a second EoS, we use the Shen Bag model (in the following we interchangeably use the abbreviation STOS/STOS-B145) from COMPOSE (Shen et al. 1998b, c; Sagert et al. 2010, 2009; Sugahara & Toki 1994). The STOS-B145 EoS uses a relativistic mean-field approach for the hadronic phase. Quark matter is described by the thermodynamic Bag model containing u, d, and s quarks (Farhi & Jaffe 1984; Greiner, Koch & Stocker 1987). The Bag model is extended by the inclusion of first-order corrections to the strong coupling constant α_s (Sagert et al. 2009, 2010). The strong coupling constant α_s, bag parameter B, and strange quark mass together determine the critical density for the mixed phase. The phase transition region is constructed via the Gibbs condition where both phases in the phase co-existing region have globally conserved charge. We use a bag parameterizations |$B=145\, \mathrm{MeV}$| and a coupling constant α_s = 0.7. The squared speed of sound of the Bag model in the pure quark phase is |$1/3\, c^2$|⁠. The maximum gravitational mass for cold matter in β-equilibrium is 2.01 M_⊙. The corresponding M–R curve has two maxima, which are connected by an unstable branch leading to the twin star phenomenon (Alford, Han & Prakash 2014). The first maximum (at larger radii) can reach up to M ∼ 2.5 M_⊙ (Zha et al. 2021).¹

2.1.3 CMF EoS

The Chiral SU(3)-flavour parity-doublet Polyakov-loop quark-hadron mean-field model (CMF) combines a mean-field description of the interaction between the lowest baryon octet (p, n, Λ, Ξ⁻, Ξ⁰, Σ⁻, Σ⁰, Σ⁺), the three light quark flavours u, d, s, and gluons, as well as contributions of the full hadron-resonance list. Both the lowest octet hadrons and quarks interactions are modelled by a chiral Lagrangian (Papazoglou et al. 1999; Steinheimer, Schramm & Stocker 2011; Motornenko et al. 2020) which allows for chiral symmetry restoration in the hadronic sector as well as in the quark sector. In addition, the thermal contributions of all other known hadronic species, including mesons and baryonic resonances, are included to properly describe QCD matter at intermediate energy densities. The CMF model thus provides a most complete description of the interactions in both the hadronic and the deconfined phases of QCD. In the CMF model, the transition between hadronic matter and quark matter is introduced by an excluded-volume formalism. The parameters of the model are chosen such that properties of nuclear matter are reproduced and the model describes lattice-QCD thermodynamics results. The model therefore incorporates a first-order nuclear liquid-vapour phase transition at densities |${\sim } \rho _\mathrm{sat}$|⁠; a second, but weak first-order phase transition occurs, due to chiral symmetry restoration, at about 4ρ_sat with a critical endpoint at |$T_\mathrm{CeP}\approx 15\, \mathrm{MeV}$|⁠. The transition to quark matter at higher densities occurs as a smooth crossover (Motornenko et al. 2020). At asymptotically high densities, the squared speed of sound approaches the Stefan–Boltzmann limit of |${\sim }1/3 \, c^2$|⁠. The CMF model predicts hybrid NSs with gravitational masses up to |${\sim } 2\,\text{M}_\odot$| for cold NS matter in β-equilibrium. The smooth nature of the crossover from hadrons to quarks does not lead to a third family branch of compact stars. All these components allow the CMF model to be applied for modelling of heavy-ion collisions (Steinheimer et al. 2010; Omana Kuttan et al. 2022), analysis of lattice QCD data (Steinheimer & Schramm 2011, 2014; Motornenko et al. 2020, 2021a; Motornenko et al. 2021b), as well as studies of cold NSs and their mergers (Most et al. 2022). Since the CMF EoS does not include heavy nuclei at sub-saturation density, we extend it to low densities using the SFHx EoS (Steiner et al. 2013), which is matched to the CMF table at densities less than |$8\times 10^{13}\, \mathrm{g}\, \mathrm{cm}^{-3}$|⁠. Note that the CMF model with a similar matching has been recently applied to describe the dynamical evolution of binary NS mergers as well as heavy-ion collisions at the SIS18 accelerator (Most et al. 2022).

2.2 Supernova simulations: numerical methods

For our supernova simulations, we employ the finite-volume neutrino hydrodynamics code CoCoNuT-FMT for solving the general-relativistic equations of hydrodynamic in spherical symmetry in Eulerian form (Müller, Janka & Dimmelmeier 2010; Müller & Janka 2015). The hydrodynamics module CoCoNuT uses higher rder piecewise parabolic reconstruction (Colella & Woodward 1984) and the relativistic HLLC Riemann solver (Mignone & Bodo 2005). We treat convection in 1D using mixing-length theory as applied previously in supernova simulations (Wilson & Mayle 1988; Müller 2015; Mirizzi et al. 2016).

For the neutrino transport, we use the fast multigroup transport fmt scheme of Müller & Janka (2015) which solves the energy-dependent neutrino zeroth moment equation for electron neutrinos, electron antineutrinos, and heavy-flavour neutrinos in the stationary approximation using a one-moment closure from a two-stream Boltzmann equation and an analytic closure at low optical depth. Neutrino interaction rates include absorption and scattering on nuclei and nucleons and bremsstrahlung for heavy-flavour neutrinos in a one-particle rate approximation; see Müller & Janka (2015), Müller et al. (2019) for details. The appearance of quarks is expected to decrease neutrino opacities in the core (Pons et al. 2001; Steiner, Prakash & Lattimer 2001; Colvero & Lugones 2014), which will impact the neutrino emission and any wind outflows from the PCS on longer time-scales. However, one can argue (Fischer et al. 2011) that over short time-scales after the second bounce, neutrino trapping is still effective in the deconfined core region. As a pragmatic solution, we therefore use the nucleonic opacities throughout, assuming a nucleonic composition compatible with charge neutrality. A proper neutrino treatment for the mixed phase and pure quark phase should be further explored in future studies.

Different treatments for the EoS and nuclear reactions are applied in various regimes. At low densities, the matter is treated as a mixture of electrons, positrons, photons, and a perfect gas of nucleons and nuclei. At temperatures below |$5\, \mathrm{GK}$|⁠, we employ a flashing treatment for nuclear reactions following (Rampp & Janka 2002); above |$5\, \mathrm{GK}$|⁠, we assume NSE. At high densities, a tabulated nuclear EoS is used. We track the mass fractions of protons, neutrons, α-particles, and 17 nuclear species, and the electron fraction Y_e in all EoS regimes. The transition density between the low- and high-density EoS regime is set to |$5\times 10^8\, \mathrm{g}\, \mathrm{cm}^{-3}$| during the collapse phase and changed to |$10^{11}\, \mathrm{g}\, \mathrm{cm}^{-3}$| after the collapse. Since our neutrino transport presently does not use consistent opacities for the quark phase, we do not add separate advection equations of the mass fractions of u, d, and s quarks. This is possible since the mass fractions merely act as passive scalars in the high-density EoS regime and do not influence the solution of the equations of hydrodynamics. Instead, we map u, d, and s quarks into neutrons and protons such as to ensure charge neutrality and baryon number conservation, i.e.

$$\begin{eqnarray} \tilde{X}_\mathrm{n} &= X_\mathrm{n} + \frac{1}{3} X_\mathrm{u} + \frac{2}{3} (X_\mathrm{d}+X_\mathrm{s}), \\ \tilde{X}_\mathrm{p} &= X_\mathrm{p} + \frac{2}{3} X_\mathrm{u} + \frac{1}{3} (X_\mathrm{d}+X_\mathrm{s}), \end{eqnarray}$$

where X_u, X_d, X_s, X_n, and X_p are the mass fractions of u, d, s quarks, neutrons, and protons in the EoS table, respectively, and |$\tilde{X}_\text{n}$| and |$\tilde{X}_\text{p}$| are the mass fractions of neutrons and protons used in the code, respectively.

2.3 Progenitor models

We use 40 progenitors in the mass range 14–|$100\, \mathrm{M}_\odot$| with two different metallicities Z = 0 and 0.012 (solar metallicity). Models that start with the letter s (e.g. s14) have solar metallicity, models that start with z have zero metallicity. The number in the model label denotes the ZAMS mass in solar masses. The progenitors have been calculated with the stellar evolution code Kepler (Weaver, Zimmerman & Woosley 1978; Heger & Woosley 2010). The solar-metallicity models are a subset of those in Müller et al. (2016).

Our progenitors differ in several respects from the stellar evolution models used in Fischer et al. (2018, 2020) and Fischer (2021), which were taken from Umeda & Nomoto (2008). Their models were based on stellar evolution calculations using a different treatment of mixing processes (Schwarzschild criterion instead of Ledoux). Furthermore, all of the progenitors from Umeda & Nomoto (2008) have very low but non-zero metallicity (Z = Z_⊙/200). In terms of mass-loss, there is no appreciable difference to our zero-metallicity models; mass-loss will be negligible in both cases. The more relevant parameters for comparison are the He, CO, and Fe core masses, which have significant influence on the supernova dynamics and are strongly dependent on the physics treatment during stellar evolution calculations. The aforementioned differences shift the relationship of core masses and the corresponding ZAMS mass. Our progenitors show a trend towards lower He, CO, and Fe core masses for a given ZAMS mass. More specifically, the |$50\, \mathrm{M}_\odot$| progenitor z50 in our study has a He core mass of |$17.78\, \mathrm{M}_\odot$| as opposed to |$21.8\, \mathrm{M}_\odot$| for the same ZAMS mass in Umeda & Nomoto (2008), Fischer et al. (2018), a CO core mass of |$11.34\, \mathrm{M}_\odot$| as opposed to |$19.3\, \mathrm{M}_\odot$|⁠, a Fe core mass of |$1.86\, \mathrm{M}_\odot$| as opposed to |$2.21\, \mathrm{M}_\odot$|⁠, and a core binding energy |$2.10\times 10^{51}\, \mathrm{erg}$| instead of |$3.67\times 10^{51}\, \mathrm{erg}$|⁠. We find that the |$50\, \mathrm{M}_\odot$| model of Umeda & Nomoto (2008) and Fischer et al. (2018) corresponds most closely to our z60 model with a He, CO and Fe mass of 23.90, 20.90, and |$2.01\, \mathrm{M}_\odot$|⁠, respectively, and a binding energy of |$3.42\times 10^{51}\, \mathrm{erg}$|⁠.

2.4 Nucleosynthesis post-processing

To follow the nucleosynthesis, we post-process the recorded temperature, density, and radius trajectories along with the neutrino fluxes and energies for ν_e, |$\bar{\nu }_{\mathrm{e}}$|⁠, and ν_x neutrinos. The ν_x species stands for the sum of the ν_μ, |$\bar{\nu }_{\mu }$|⁠, ν_τ, and |$\bar{\nu }_\tau$| contributions. The local neutrino energy density and mean energy at the current location of each mass shell are taken directly from the hydrodynamics code.

After the end of the hydrodynamic simulations at time t = t_h, we extrapolate the trajectories to time |$t_\mathrm{f}=1\, \mathrm{yr}$| after the explosion assuming adiabatic homologous expansion where we crudely estimate the velocity from the current radial coordinate and the time since the onset of core collapse:

$$\begin{eqnarray} T(t) = T(t_\mathrm{h}) \times \left(\frac{t}{t_\mathrm{h}}\right)^{-1} \,\,,\quad \rho (t) = \rho (t_\mathrm{h}) \times \left(\frac{t}{t_\mathrm{h}}\right)^{-3} \,\,. \end{eqnarray}$$

(1)

We additionally impose a minimum temperature of |$10^6\, \mathrm{K}$|⁠. At that stage, all regular nuclear reactions are frozen out, only radioactive decays still occur, which we do follow.

For the neutrinos, we assume that the neutrino luminosities, L_ν, and energies, E_ν, decay exponentially after time t_h,

$$\begin{eqnarray} L_\nu (t) &=& L_\nu (t_\mathrm{h}) \times \exp \left(\frac{t_\mathrm{h}-t}{\tau _\nu }\right) \,\,,\quad \end{eqnarray}$$

(2)

$$\begin{eqnarray} E_\nu (t) &=& E_\nu (t_\mathrm{h}) \times \exp \left(\frac{t_\mathrm{h}-t}{4\, \tau _\nu }\right) \,\,, \end{eqnarray}$$

(3)

using a characteristic time-scale of |$\tau _{\nu }=3\, \mathrm{s}$|⁠. The neutrino temperatures (in MeV) are approximated from the neutrino energies as T_ν = E_ν/3.15. The neutrino temperatures in the network are explicitly limited to a range of 3–|$8\, \mathrm{MeV}$| for ν_e, to 3–|$12\, \mathrm{MeV}$| for |$\bar{\nu }_\mathrm{e}$|⁠, and to 4–|$12\, \mathrm{MeV}$| for ν_x.

We use a modified standalone version of the first-order implicit adaptive nuclear reaction network from Kepler (Rauscher et al. 2002), which includes all nuclear species and reactions up to astatine except fission. The network also includes neutrino-induced spallation as described in (Heger et al. 2005). Compared to the Kepler version, the standalone version (Burn code) of the network adds new features for higher computational accuracy. Due to the potentially large time-step in the trajectories, we implement a new iterative adaptive network that repeats a time step until there is no more addition of new species rather than only adjusting the network at the end of a time-step in preparation for the next time-step. We also iteratively sub-cycle the network calculation when the abundance changes are too large or when mass conservation is violated by more than one in 10⁻¹⁴, and linearly interpolate thermodynamic and neutrino quantities from the recorded or extrapolated trajectory grid points.

For temperatures exceeding |$10^{10}\, \mathrm{K}$|⁠, we impose the Y_e from the hydrodynamics code assuming a composition with free nucleons only. When the temperature drops below this threshold temperature, we follow the full network using this new Y_e as a starting point. This accounts for the limited temperature range for neutrino interactions in Kepler and takes advantage of the full non-thermal neutrino energy distribution in CoCoNuT.

3 RESULTS

We summarize key outcomes from all simulations in Tables 1–3, including the time of the phase transition, the presence or absence of a second bounce and neutrino burst, the central lapse function α at the end of the simulation, and the occurrence or absence of an explosion. In addition to the ZAMS mass of the models, it is useful to consider the compactness parameter ξ_2.5 at t = 0, where

$$\begin{eqnarray} \xi _M=\frac{M/\text{M}_\odot }{r(M)/1000\, \text{km}}, \end{eqnarray}$$

(4)

for any mass coordinate M. The compactness parameter is one empirical predictor for the likelihood of a star to explode by the neutrino-driven mechanism (O’Connor & Ott 2011), but it will also prove useful as a stellar structure metric in the context of PT-driven explosions. The compactness parameters for all progenitor models are also plotted in Fig. 2. In addition, the tables show the explodability parameters M₄ and μ₄ of Ertl et al. (2016), which are related to the mass of the iron-silicon core (using the mass shell where the entropy s reaches 4k_B/baryon) and the density outside the Si/O shell interface.

Figure 2.

Compactness parameter ξ_2.5 for the set of progenitors used in our simulations. Yellow stars denote successful explosions with the DD2F_SF EoS, black dots with green/pink edging denote failed explosions for DD2F_SF/STOS where the shock wave is not powerful enough to propagate through the infalling material or the maximum compact star mass is exceeded due to accretion during the explosion phase. Blue circles comprise all other non-exploding progenitors.

Table 1.

Progenitor properties and results for the 1D simulations using the DD2-1.4 EoS. M is the ZAMS mass, M_PCS is the mass of the PCS at the time of the second collapse, t_BH is the time of BH formation (models that did not form BH after t ≥ >6 s are marked with ‘–’), t_{1, b} is the time of the first bounce, t_MP is the time when the progenitor reaches the mixed-phase region, t_{2, b} is the time of the second bounce, α is the central lapse function at the end of the simulation, ξ_2.5 is the progenitor compactness, μ₄ and M₄μ₄ are the explodability parameters from Ertl et al. (2016) (evaluated at t = 0). The column |$C^\mathrm{inv}_\mathrm{L}$| refers to the presence or absence of inverted convection after the phase transition.

Metallicity	M	MPCS(tb, 2)	EoS	Explosion	tBH	t1, b	tMP	t2, b	α	ξ2.5	μ4	M4μ4	|$C^\text{inv}_\text{L}$|
	(M⊙)	(M⊙)			(s)	(s)	(s)	(s)
Solar	14	1.79	DD2F_SF	No	2.056	0.223	2.034	2.056	0.397	0.115	0.079	0.124	Yes
Solar	15	1.87	DD2F_SF	No	1.795	0.240	1.782	1.795	0.315	0.165	0.082	0.141	Yes
Solar	16	1.76	DD2F_SF	Yes	–	0.176	1.971	2.007	0.488	0.142	0.073	0.108	Yes
Solar	17	1.75	DD2F_SF	No	2.187	0.218	2.129	2.187	0.423	0.127	0.052	0.084	Yes
Solar	18	1.95	DD2F_SF	No	1.542	0.252	1.541	1.542	0.211	0.248	0.101	0.184	Yes
Solar	19	1.94	DD2F_SF	No	1.599	0.263	1.560	1.599	0.227	0.246	0.095	0.167	Yes
Solar	20	1.97	DD2F_SF	No	1.516	0.264	1.515	1.516	0.258	0.237	0.103	0.195	Yes
Solar	21	2.19	DD2F_SF	No	1.048	0.312	1.048	1.048	0.200	0.432	0.152	0.321	Yes
Solar	22	2.15	DD2F_SF	No	1.126	0.307	1.126	1.126	0.186	0.402	0.144	0.295	Yes
Solar	23	2.04	DD2F_SF	No	1.368	0.290	1.368	1.368	0.222	0.312	0.109	0.208	Yes
Solar	24	1.97	DD2F_SF	No	1.560	0.272	1.368	1.560	0.174	0.265	0.092	0.168	Yes
Solar	25	2.03	DD2F_SF	No	1.388	0.285	1.387	1.388	0.228	0.304	0.110	0.205	Yes
Solar	30	2.21	DD2F_SF	No	1.022	0.322	1.021	1.022	0.264	0.474	0.194	0.405	Yes
Solar	35	2.33	DD2F_SF	No	0.876	0.351	0.876	0.876	0.225	0.597	0.261	0.597	Yes
Primordial	14	1.71	DD2F_SF	Yes	–	0.186	2.220	2.336	0.513	0.047	0.044	0.071	Yes
Primordial	15	1.79	DD2F_SF	Failed	–	0.156	3.320	2.554	0.495	0.072	0.039	0.060	Yes
Primordial	16	1.77	DD2F_SF	No	2.063	0.206	2.045	2.063	0.414	0.136	0.074	0.117	Yes
Primordial	17	1.91	DD2F_SF	No	1.656	0.226	1.654	1.656	0.297	0.213	0.103	0.180	Yes
Primordial	18	1.76	DD2F_SF	Failed	–	0.192	2.227	2.467	0.458	0.114	0.051	0.076	Yes
Primordial	19	1.87	DD2F_SF	Failed	1.727	0.214	1.726	1.727	0.317	0.197	0.090	0.147	Yes
Primordial	20	1.84	DD2F_SF	No	1.730	0.192	1.726	–	0.354	0.170	0.102	0.150	Yes
Primordial	21	1.86	DD2F_SF	No	1.680	0.198	1.677	–	0.347	0.202	0.099	0.153	Yes
Primordial	22	1.83	DD2F_SF	No	1.746	0.182	1.743	–	0.346	0.206	0.094	0.142	Yes
primordial	23	1.86	DD2F_SF	No	1.769	0.224	1.768	1.769	0.323	0.190	0.086	0.140	Yes
primordial	24	2.16	DD2F_SF	No	1.105	0.308	1.105	1.105	0.258	0.413	0.163	0.325	Yes
Primordial	25	2.21	DD2F_SF	No	0.996	0.304	0.996	0.996	0.265	0.446	0.153	0.332	Yes
Primordial	30	1.95	DD2F_SF	No	1.539	0.244	1.538	1.539	0.268	0.237	0.106	0.185	Yes
Primordial	35	2.25	DD2F_SF	No	0.944	0.320	0.944	0.944	0.248	0.487	0.190	0.407	Yes
Primordial	40	2.37	DD2F_SF	No	0.816	0.333	0.816	0.816	0.127	0.645	0.366	0.764	Yes
Primordial	45	2.36	DD2F_SF	No	0.826	0.333	0.826	0.826	0.226	0.635	0.345	0.740	Yes
Primordial	50	2.28	DD2F_SF	No	0.864	0.332	0.826	0.864	0.215	0.593	0.261	0.584	Yes
Primordial	55	2.34	DD2F_SF	No	0.808	0.305	0.807	0.808	0.233	0.635	0.385	0.679	Yes
Primordial	60	2.30	DD2F_SF	No	0.873	0.299	0.872	0.873	0.223	0.585	0.255	0.489	Yes
Primordial	65	2.32	DD2F_SF	No	0.888	0.316	0.888	0.888	0.229	0.587	0.248	0.490	Yes
Primordial	70	2.40	DD2F_SF	No	0.860	0.350	0.860	0.860	0.212	0.639	0.244	0.533	Yes
Primordial	75	2.44	DD2F_SF	No	0.825	0.352	0.824	0.825	0.198	0.678	0.279	0.601	Yes
Primordial	80	2.48	DD2F_SF	No	0.820	0.378	0.819	0.820	0.257	0.711	0.283	0.640	Yes
Primordial	85	2.64	DD2F_SF	No	0.752	0.412	0.752	0.752	0.220	0.872	0.348	0.844	Yes
Primordial	90	2.67	DD2F_SF	No	0.725	0.416	0.706	0.725	0.209	0.922	0.405	0.981	Yes
Primordial	95	2.63	DD2F_SF	No	0.769	0.433	0.765	0.769	0.232	0.875	0.347	0.854	Yes

Table 2.

Properties of the 1D simulations for 21 different progenitors using the STOS B145 EoS. M is the ZAMS mass, M_PCS is the mass of the PCS at the time of the second collapse, t_{1, b} is the time of the first bounce, t_MP is the time when the progenitor reaches the mixed-phase region, t_{2, b} is the time of the second core bounce, α is the central lapse function at the end of the simulation, ξ_2.5 is the progenitor compactness, μ₄ and M₄μ₄ are the explodability parameters from Ertl et al. (2016) (evaluated at t = 0). |$C^\mathrm{inv}_\mathrm{L}$| refers to the presence (+) or absence (-) of inverted convection after the phase transition.

Metallicity	M	MPCS(tb, 2)	EoS	Explosion	tBH	t1, b	tMP	t2, b	α	ξ2.5	μ4	M4μ4	|$C^\text{inv}_\text{L}$|
	(M⊙)	(M⊙)			(s)	(s)	(s)	(s)
Solar	21	2.19	STOS	no	1.037	0.262	0.356	–	0.176	0.432	0.152	0.321	Yes
Solar	22	2.19	STOS	No	1.289	0.258	0.374	–	0.150	0.402	0.144	0.295	Yes
Solar	30	2.20	STOS	No	0.961	0.272	0.390	–	0.173	0.474	0.194	0.405	Yes
Solar	35	2.21	STOS	No	0.731	0.296	0.404	–	0.185	0.597	0.261	0.597	Yes
Primordial	23	2.26	STOS	No	4.486	0.190	0.310	–	0.214	0.190	0.086	0.140	Yes
Primordial	24	2.19	STOS	No	1.220	0.260	0.374	–	0.191	0.413	0.163	0.325	Yes
Primordial	25	2.20	STOS	No	0.939	0.258	0.370	–	0.168	0.446	0.153	0.332	Yes
Primordial	35	2.20	STOS	No	0.842	0.272	0.372	–	0.153	0.487	0.190	0.407	Yes
Primordial	40	2.21	STOS	No	0.670	0.284	0.390	–	0.214	0.645	0.366	0.764	Yes
Primordial	45	2.21	STOS	No	0.680	0.284	0.402	–	0.157	0.635	0.345	0.740	Yes
Primordial	50	2.21	STOS	No	0.735	0.284	0.386	–	0.109	0.593	0.261	0.584	Yes
Primordial	55	2.21	STOS	No	0.694	0.264	0.356	–	0.187	0.635	0.385	0.679	Yes
Primordial	60	2.20	STOS	No	0.755	0.258	0.354	–	0.169	0.585	0.255	0.489	Yes
Primordial	65	2.20	STOS	No	0.746	0.268	0.388	–	0.170	0.587	0.248	0.490	Yes
Primordial	70	2.25	STOS	Failed	0.658	0.304	0.402	0.658	0.177	0.639	0.244	0.533	Yes
Primordial	75	2.28	STOS	Failed	0.650	0.306	0.406	0.650	0.133	0.678	0.279	0.601	Yes
Primordial	80	2.33	STOS	No	0.642	0.330	0.422	–	0.198	0.711	0.283	0.640	Yes
Primordial	85	2.48	STOS	No	0.612	0.362	0.452	–	0.181	0.872	0.348	0.844	Yes
Primordial	90	2.52	STOS	No	0.607	0.370	0.470	–	0.178	0.922	0.405	0.981	Yes
Primordial	95	2.51	STOS	No	0.620	0.378	0.480	–	0.193	0.875	0.347	0.854	Yes
Primordial	100	2.21	STOS	No	0.731	0.296	0.418	–	0.172	0.597	0.261	0.597	Yes

Table 3.

Summary of models (analogous to Table 1) using the CMF EoS. M is the ZAMS mass, M_PCS is the mass of the PCS at the time of the second collapse, t_BH is the final time of the simulation (t_BH is the time of BH formation (models which cooled before reaching the crossover region and did not collapse to a BH within t ≤ 6 s are marked with ‘–’), t_{1, b} is the time of the first bounce, t_{2, b} is the time of the second core bounce, α is the central lapse function at the end of the simulation, ξ_2.5 is the progenitor compactness, μ₄ and M₄μ₄ are the explodability parameters from Ertl et al. (2016) (evaluated at t = 0). |$C^\mathrm{inv}_\mathrm{L}$| refers to the presence (+) or absence (-) of inverted convection after the phase transition. Some progenitors show very slight convection (with a negative entropy gradient).

Metallicity	M	MPCS(tb, 2)	EoS	Explosion	tBH	t1, b	t2, b	α	ξ2.5	μ4	M4μ4	|$C^\text{inv}_\text{L}$|
	(M⊙)	(M⊙)			(s)	(s)	(s)
Solar	16	2.22	CMF	No	–	0.196	–	0.529	0.142	0.073	0.108	Yes
Solar	18	2.49	CMF	No	4.641	0.272	–	0.229	0.248	0.101	0.184	Yes
Solar	19	2.48	CMF	No	4.683	0.282	–	0.230	0.246	0.095	0.167	Yes
Solar	20	2.49	CMF	No	5.345	0.282	–	0.185	0.237	0.103	0.195	Yes
Solar	21	2.41	CMF	No	1.885	0.326	–	0.177	0.432	0.152	0.321	Yes
Solar	22	2.41	CMF	No	2.213	0.322	–	0.221	0.402	0.144	0.295	Yes
Solar	23	2.45	CMF	No	3.472	0.312	–	0.203	0.312	0.109	0.208	Yes
Solar	24	2.48	CMF	No	4.450	0.294	–	0.188	0.265	0.092	0.168	Yes
Solar	25	2.45	CMF	No	3.481	0.306	–	0.185	0.304	0.110	0.205	Yes
Solar	30	2.40	CMF	No	1.868	0.336	–	0.201	0.474	0.194	0.405	Yes
Primordial	14	1.90	CMF	No	–	0.208	–	0.602	0.047	0.044	0.071	Yes
Primordial	15	1.83	CMF	No	–	0.176	–	0.616	0.072	0.039	0.060	Yes
Primordial	16	2.10	CMF	No	–	0.228	–	0.559	0.136	0.074	0.117	Yes
Primordial	17	2.40	CMF	No	–	0.250	–	0.474	0.213	0.103	0.180	Yes
Primordial	18	1.99	CMF	No	–	0.212	–	0.582	0.114	0.051	0.076	Yes
Primordial	19	2.45	CMF	No	–	0.238	–	0.447	0.197	0.090	0.147	Yes
Primordial	20	2.36	CMF	No	–	0.214	–	0.487	0.170	0.102	0.150	Yes
Primordial	21	2.51	CMF	No	5.861	0.220		0.217	0.202	0.099	0.153	Yes
Primordial	22	2.51	CMF	No	5.978	0.202	–	0.200	0.206	0.094	0.142	Yes
Primordial	23	2.39	CMF	No	5.899	0.246	–	0.477	0.190	0.086	0.140	Yes
Primordial	24	2.41	CMF	No	2.256	0.326	–	0.209	0.413	0.163	0.325	Yes
Primordial	25	2.41	CMF	No	1.722	0.322	–	0.183	0.446	0.153	0.332	Yes
Primordial	30	2.49	CMF	No	5.034	0.264	–	0.230	0.237	0.106	0.185	Yes
Primordial	35	2.42	CMF	No	1.545	0.338	–	0.188	0.487	0.190	0.407	Yes
Primordial	40	2.49	CMF	No	0.960	0.348	–	0.125	0.645	0.366	0.764	Yes
Primordial	45	2.48	CMF	No	0.986	0.348	–	0.077	0.635	0.345	0.740	Yes
Primordial	50	2.47	CMF	No	1.082	0.348	–	0.181	0.593	0.261	0.584	Yes
Primordial	55	2.49	CMF	No	0.959	0.328	–	0.148	0.635	0.385	0.679	Yes
Primordial	60	2.46	CMF	No	1.073	0.322	–	0.131	0.585	0.255	0.489	Yes
Primordial	65	2.46	CMF	No	1.091	0.334	–	0.016	0.587	0.248	0.490	Yes
Primordial	70	2.49	CMF	No	1.009	0.366	–	0.192	0.639	0.244	0.533	Yes
Primordial	75	2.51	CMF	No	0.933	0.370	–	0.190	0.678	0.279	0.601	Yes
Primordial	80	2.54	CMF	No	0.900	0.392	–	0.175	0.711	0.283	0.640	Yes
Primordial	85	2.59	CMF	No	0.744	0.424	–	0.180	0.872	0.348	0.844	Yes
Primordial	95	2.60	CMF	No	0.747	0.440	–	0.210	0.875	0.347	0.854	Yes
Primordial	100	2.41	CMF	No	1.668	0.308	–	0.216	0.404	0.176	0.349	Yes

Among all 97 models, we find only two explosions by the PT-driven mechanism, both for the DD2F_SF EoS, namely the zero-metallicity model z14 and the solar-metallicity model s16. We, therefore, consider the DD2F_SF series first, with a particular focus on the two exploding models.

3.1 DD2F_SF Series

3.1.1 Progenitor-dependent outcomes

The two exploding models z14 and s16 point to a remarkable difference from the scenario of hyperenergetic explosions from very massive progenitors in Fischer et al. (2018). Those two models both have low compactness parameters ξ_2.5 < 0.15 (marked as yellow stars in Fig. 2). Furthermore, they only reach low explosion energies. In Fig. 3, we show the evolution of the diagnostic explosion energy (computed in general relativity following Müller, Janka & Marek 2012) and shock trajectories for these two models. We find values of only |$1.3\times ~10^{50}\, \mathrm{erg}$| for z14 and |$3.6\times ~10^{49}\, \mathrm{erg}$| for s16 at the end of the simulations. Model s16, transiently reaches a diagnostics energy of about |$3\times 10^{50}\, \mathrm{erg}$| immediately after the second bounce, but the diagnostic energy then steadily decreases as the shock scoops up bound outer layers of the star. In Model z14, the shock launched by the second bounce is noticeably weaker; the diagnostic energy never exceeds |$10^{50}\, \mathrm{erg}.$|

Figure 3.

Diagnostic explosion energy E_expl (top panel) and shock radius R_shock (bottom panel) as function of time after the second bounce for models z14 (orange) and s16 (blue) using the DD2F_SF EoS. By the end of the simulation, we find diagnostics energies of |$1.25\times 10^{50}\, \mathrm{erg}$| and |$3.64\times 10^{49}\, \mathrm{erg}$| for s16 and z14, respectively. In both cases, much of the initial explosion energy is drained as the shock scoops up bound shells.

Not all of the other models form BHs quietly, however. Three models in the DD2F_SF series turn out as failed explosions. We call an explosion failed if there is a second bounce after the phase transition and the second shock initially propagates dynamically with positive post-shock velocities through the neutrinosphere that leads to a second neutrino burst, but then proves too weak to propagate through outer infalling material and stalls again. Similar to Zha et al. (2021), we note oscillating behaviour for the failed explosion model z18 in which the PCS oscillates for several milliseconds before it collapses into a BH. The trend towards low compactness parameters for failed explosions is similar to exploding models- we only see failed explosions for ξ_2.5 ≲ 0.2 (models z15, z18, and z19).

In order to gain more insight into the influence of the phase transition after the second collapse, we consider the pressure-weighted mean adiabatic index Γ as a measure for the stability of the PCS (Goldreich & Weber 1980),

$$\begin{eqnarray} \overline{\Gamma } = \frac{\int _\text{PCS} \Gamma \, P\, \text{d}V}{\int _\text{PCS} P\, \text{d}V}, \end{eqnarray}$$

(5)

where the integration volume covers the PCS, which we define by a threshold density of |$10^{11}\, \mathrm{g}\, \mathrm{cm}^{-3}$|⁠. The relation between adiabatic index and stellar stability is well known (Chandrasekhar 1964). For spherical stars in non-isotropic matter, instability can be shown to set in when the pressure-weighted adiabatic index decreases below 4/3 (Goldreich & Weber 1980), though corrections apply in general relativity. The adiabatic index usually decreases during phase transitions. It is apparent that this softening causes the initial contraction of the PCS.

We plot |$\bar{\Gamma }$| as function of time after the first bounce in Fig. 4. The progenitors are colour-coded by their compactness parameter ξ_2.5 at t = 0. High-compactness progenitors reach the phase transition earlier due to the shorter free-fall time-scale of the shells outside the iron core. High-compactness progenitors also show a lower mean adiabatic index compared to low-compactness progenitors. This makes high-compactness progenitors more susceptible to bulk collapse to a BH (without a second bounce) as soon as they hit the phase transition. However, for the DD2F_SF EoS, only 3 out of 40 progenitors do not exhibit a second bounce. The two exploding models s16 and z14 (thick dashed and dotted curves, respectively) exhibit some of the highest values of |$\bar{\Gamma }$|⁠, although not the highest pones altogether. As a reference, we also mark the failed CCSNe as thick dashed/dotted lines. The dependence of |$\bar{\Gamma }$| on the compactness parameter reflects the trend towards higher average PCS entropy in high-compactness parameters noted by Da Silva Schneider et al. (2020) and Zha et al. (2021). In addition, the higher chance of a PT-driven explosion for low-compactness progenitors is also consistent with another phenomenon described by Zha et al. (2021), who noted that the phase transition for oscillating progenitors occurs well below the maximum mass for a stable PCS, and at lower average PCS entropy. At the time of the second collapse, the PCS masses of the exploding models are 1.71 M_⊙ for model z14 and 1.76 M_⊙ for model s16.

Figure 4.

Evolution of the mean adiabatic index |$\bar{\Gamma }$| (defined in equation 5) as function of post-bounce time, colour-coded by the compactness parameter, for the DD2F_SF models. The thick (dashed or dotted) lines show progenitors that either exploded (z14, s16), or failed to explode (z15, z18, and z19) after the emission of a second neutrino burst due to the stalling of the second shock.

3.1.2 DD2F_SF Series – detailed dynamics of selected models

We continue our discussion by the outlining dynamic properties of the exploding models z14 and s16 and the failed explosion model z18. Radial profiles of the density ρ, temperature T, electron fraction Y_e, radial velocity v, and entropy per baryon s² at five epochs from |$2\, \mathrm{ms}$| before the star reaches the mixed phase up to 400 ms after the second bounce are shown in Fig. 5.

Figure 5.

Selected radial profiles for z14, s16, and z18 around the time of second collapse using DD2F_SF. The panels show (from the top to bottom panels) density ρ, temperature T in |$10^{11}\, \mathrm{K}$|⁠, Y_e, radial velocity v, and baryon-specific entropy s as functions of radius at five different time-steps with t = 0 corresponding to when the star enters the mixed phase, which marks the onset of the second collapse.

Shortly before the centre of the PCS enters the mixed phase (purple curves), continuous accretion has brought the central density up to several times saturation density, and the density, temperature, and entropy exhibit typical steady-state accretion profiles. Strong exposure to neutrino heating raises Y_e in the accreted matter to values above 0.5 in much of the gain region, before it drops as material settles in the cooling region. This will later become relevant for nucleosynthesis. Once the mixed phase is reached, the evolution becomes extremely fast. On time-scales of less than a millisecond, the PCS goes through rapid, near-homologous contraction due to the lower adiabatic index in the mixed phase, and due to increasing densities, a major part of the PCS is soon composed of pure quark matter. By the time corresponding to the red curves, the collapsing PCS has already undergone a rebound due to increasing quark-vector interactions, halting the collapse at about |$\sim 2\times 10^{15}\, \, \mathrm{g}\, \mathrm{cm}^{-3}$|⁠. The vector interactions in the quark phase go along with violations of lattice QCD data which we briefly touched upon in Section 2. The shock wave from the second bounce has already broken out of the PCS, overtaken the primary shock, and reached a radius of order |$100\, \mathrm{km}$|⁠, heating the post-shock matter to |$40\tt {-}80 \, k_\mathrm{B}/\mathrm{baryon}$|⁠. The velocity profiles around the time of the second bounce (red, dark yellow) still show strong ringdown oscillations of the PCS, which launch further secondary shocks that somewhat boost the entropy in the ejected matter. Several shock waves from the bounce and ringdown oscillations are formed within the first 2 ms: The shock wave of model z14 (red line in Panel j) at about |$200\, \mathrm{km}$| comes from the initial second collapse and bounce. Another rebound occurs at |$3\, \mathrm{km}$| after |${\sim } 0.5\, \mathrm{ms}$| (on shorter time-scales than displayed here). Multiple core oscillations follow within time-scales |${\lesssim } 2\, \mathrm{ms}$| and lead to the ‘jagged’ velocity profile which we observe for successful as well as failed explosions.

The promptly ejected matter previously located in the gain region mostly maintains Y_e > 0.5 as it is run over by the shock from the second bounce, but the profiles show variations in Y_e in the promptly ejected matter; this can be understood as the result of neutrino irradiation by the intense electron antineutrino burst (Sagert et al. 2009) from the breakout of the second shock through the neutrinosphere.

Already right after the second bounce, the failed explosion model z18 exhibits noticeably smaller positive velocities immediately behind the main (outermost) shock, and smaller amplitudes of the additional shock waves launched by ringdown oscillations. In this model, the shock still propagates out to several hundred kilometres and then stalls again (Panel l, green curve). The shock then recedes to a few tens of kilometres, and the PCS continues to accrete for several hundred milliseconds after the phase transition. In models z14 and s16, the shock continues to propagate outwards, but the post-shock velocity decreases considerably as the initial kinetic energy of the shock is used up to unbind the shells around the PCS. Both of these models develop a neutrino-driven wind from the PCS with high entropy |${\gtrsim }100 \, k_\mathrm{B}/\mathrm{baryon}$| (panels m and n, cyan curves) in the wake of the explosion. The outflow velocities in the wind remain modest, but the developing wind can clearly be seen for model s16 (Panel k, green curve), where the wind crashes into the earlier ejecta in a reverse shock at a few hundred kilometres. We will touch upon the dynamics of the wind and the reverse shock later in Section 3.3, as it leaves interesting traces in the neutrino signal.

With only two exploding models and a few more failed explosions, it is difficult to ascertain the reason for variations in the initial energy of the second shock (and hence the outcome of the phase transition). The trajectories of the central density and temperature in the phase diagram provide some tentative hints, though. All except two progenitors in the DD2F_SF setup undergo a second core bounce. The trajectories of central temperature and density of selected models from the DD2F_SF series are shown in Fig. 6 together with the adiabatic index (background colour), adiabats, and isocontors for the quark fraction at a fixed electron fraction Y_e = 0.25, which approximates the conditions at the centre of the PCS.

Figure 6.

Phase diagram with trajectories of the central density and temperature (thick solid and dashed lines) of selected DD2F_SF models. The two exploding progenitors are thick red and green dashed and the failed explosions of z18 and z15 are shown in thick dashed purple and blue, respectively. The background displays the colour-coded adiabatic index at fixed electron fraction Y_e = 0.25 for the DD2-RMF-1.4 EoS. Two isocontours (dotted, yellow) show quark fractions of X_q = 0.01 and 1. The black dashed lines are isentropes for different entropy values (indicated as contour labels). Most of the progenitors used in this setup with the exception of two models (z20, z21) that underwent a second core bounce.

After the first bounce and before the PCS enters the mixed phase, the central density and temperature of the PCS increase due to mass accretion to |$\log \rho /(\mathrm{g}\, \mathrm{cm}^{-3})\approx 14.9$|⁠, log (T/K) ≈ 11.45. At the onset of the mixed phase, lower, and intermediate progenitors show a sudden jump in temperature which is due to inverse convection as we will discuss in Section 3.2.3. When the PCS enters the mixed phase and the combined mass fraction of u and d-quarks reaches |${\sim } 1{{\ \rm per\ cent}}$| (dotted yellow curve³), the trajectories bend abruptly at approximately constant density |$\log \rho /(\mathrm{g}\, \mathrm{cm}^{-3})=14.8$|⁠. The progenitors are not yet collapsing at this point but the adiabatic index in the mixed phase decreases, softening of the EoS leads to a core contraction of the progenitors within a few milliseconds, and collapse follows once enough matter in the PCS is converted into mixed-phase matter. The adiabatic collapse is accompanied by a decrease in temperature. Such a decrease in temperature under adiabatic compression during a phase transition may seem peculiar, but for the QCD phase transition this behaviour can be connected to the higher entropy in the quark phase under isothermically compression (i.e. along horizontal lines in Fig. 6) as we shall discuss later in Section 3.2.3. When the quark fraction approaches unity, the repulsive quark interactions come into play and significantly stiffen the EoS at densities |$15~\le ~\log \rho /(\mathrm{g}\, \mathrm{cm}^{-3}) \le ~15.1$|⁠. The adiabatic index Γ is highest immediately after the conversion to quark matter. The models with a second bounce overshoot the ‘ridge’ of maximum Γ to various degrees before the central density reaches its maximum – corresponding to the bounce – and then oscillate back and forth. The second bounce occurs at densities of |$\log \rho /(\mathrm{g}\, \mathrm{cm}^{-3}) =15.1$| for z14 and z18 and slightly later at |$\log \rho /(\mathrm{g}\, \mathrm{cm}^{-3})=15.22$| for s16. The latter shows the largest increase in density during the collapse, implying that more gravitational binding energy is released followed by a relatively strong core bounce. This causes a larger explosion energy |$E_\text{exp}=1.25\times 10^{50}\, \mathrm{erg}$|⁠, compared to z14 with |$E_\text{exp}=3.64 \times 10^{49}\, \mathrm{erg}$|⁠. The failed explosion models reach a lower maximum density during the second bounce. This provides a possible explanation for why explodability by the PT mechanism is not related to compactness (or any other obvious structure parameter) in a simple, monotonic fashion: For a successful explosion, the second collapse must neither proceed too violently to high densities since this would result in BH formation, but must concurrently reach sufficiently high densities to launch a strong shock wave. In summary, there is no straightforward way to quantitatively predict the maximum density during the second bounce and the energy of the shock wave solely from progenitor or PCS parameters. Furthermore, CCSNe were explained as a transition from a second to a third family of hot compact stars (Hempel et al. 2016). When the gravitational mass of the PCS exceeds the maximum supported (entropy dependent) maximum PCS mass, the star collapses (Da Silva Schneider et al. 2020). In addition, it was discussed that the mass ratio of both maxima in the mass–radius diagram can determine whether the star collapses into a BH immediately or bounces (Zha et al. 2021). This third family topology in the mass–radius diagram, is entropy dependent.⁴ The STOS-B145 EoS evinces this characteristic of a third family of compact stars at higher entropies.

Whereas we are not able to unambiguously pin down the exact mechanism leading to a successful explosion due to the small number of exploding models, we emphasize that it is the particular – although otherwise unremarkable – core structure of the specific models that lead to the explosion. The core structure strongly varies with mass, as pointed out by Müller et al. (2016), Sukhbold, Woosley & Heger (2018), and Sukhbold & Woosley (2014). Metallicity only plays a minor role, largely leading to some shift of the mass range where a core structure of similar explodability may be incurred. Other parameters, such as stellar rotation, or modelling parameters such as uncertainties in mass-loss, nuclear reaction rates, and mixing processes, may have an even stronger effect. From stellar modelling experience, within reasonable limits, we should expect that all of these predominately shift the mass ranges in similar ways to metallicity. Interestingly, both explosions occur in the intermediate CCSN mass range, for quite similar core structure. The larger explosion energy for s16 is not an effect of metallicity per se.

3.2 EoS comparison

In addition to the DD2F_SF models, we simulated the collapse of 36 progenitors using the CMF EoS and 21 progenitors using the STOS EoS. None of these models developed at PT-driven explosion. We find a second core bounce followed by a second neutrino burst for the STOS EoS for two high-compactness models z70 and z75. None of the CMF models exhibited a second bounce. In the following, we analyse the thermodynamic conditions at the centre of the PCS in these models and identify the reason for the disparate outcomes for the various EoS. Important outcomes and parameters for the simulations using the STOS EoS and CMF EoS are given in Tables 2 and 3. A comparison of the transition density ρ_PT for the different EoS is shown in Table 4, along with other characteristic properties of the phase transition, such as the minimum adiabatic index.

Table 4.

Properties of the mixed phases/crossover region of the DD2, STOS, and CMF EoS at fixed electron fraction Y_e = 0.25. The STOS EoS shows the lowest onset density |$\rho ^\text{onset}_\text{PT}$| of the phase transition, followed by DD2F_SF. The CMF model has the smallest ‘width’, i.e. density range |$\Delta \log \rho ^\text{onset}_\text{PT}$| of the crossover region. STOS features a wide mixed-phase band, about 20 times the width of CMF crossover region. The minimum value |$\Gamma _\mathrm{min}^\mathrm{PT}$| of the adiabatic index for DD2F_SF is |${\sim }0.3$|⁠, about four times smaller than for STOS. DD2F_SF also reaches significantly higher maximum values Γ_max(ρ > ρ_PT) in the mixed phase/pure quark phase compared to STOS and CMF.

	Nature of phase transition	|$\log \rho ^\text{onset}_\text{PT}$| (g cm−3)	Δlog ρPT (g cm−3)	Enthalpic/entropic	|$\Gamma ^\mathrm{PT}_\text{min}$|	Γmax(ρ > ρPT)
DD2F_SF	First order (Maxwell constr)	14.6–15.0	∼0.1–0.5	Yes/no	∼0.3	∼5.2
STOS	First order (Gibbs constr)	14.5–14.1	∼1.1	Yes/yes	∼1.2	∼2.8
CMF	Smooth crossover (‘∞’-order)	14.9–15.0	∼0.05	Yes/no	∼0.5	∼1.2

3.2.1 STOS models

Similar to the DD2F_SF phase diagram in Fig. 6, Fig. 7 shows the evolution of the central temperature and density on top of the adiabatic index (colour-code background) and isocontours for the quark fraction at fixed electron fraction Y_e = 0.25. The appearance of quarks is shown in dotted yellow for X_q = ∑_{i ∈ u, d, s}X_i = 0.01. X_q = 1.0 marks the end of the mixed phase. The STOS EoS features an early onset of the mixed phase at densities ρ_PT with |$14 \lesssim \log (\rho _\text{PT} /\mathrm{g}\, \mathrm{cm}^{-3}) \lesssim 14.5$| for 11.8 ≲ log (T/K) ≲ 12.

Figure 7.

Phase diagram with trajectories of the central density and temperature (thick green and black dashed lines) of two selected STOS-Bag145 models. Out of all simulations in the STOS setup, two progenitors (z70, z75) underwent a second core bounce, followed by a second neutrino burst. Those two progenitors are marked as dashed black (z70) and dash–dotted green (z75). The background displays the colour-coded adiabatic index at fixed electron fraction Y_e = 0.25 for the STOS-Bag145 EoS. Two isocontours (dotted, yellow) show quark fractions of X_q = 0.01 and 1. The black dashed lines are isentropes for different entropy values (indicated as contour labels).

The transition density shifts to lower densities for higher temperatures. This trend is seen in CMF and DD2F_SF as well and is a generic feature of the QCD phase diagram.

The adiabatic index for STOS looks qualitatively different from the DD2F_SF EoS. For STOS, the adiabatic index is roughly constant at the onset of the mixed phase with values Γ ≈ 3 until |$\log \rho /(\mathrm{g}\, \mathrm{cm}^{-3})~\approx ~14.75$| and temperatures below log T/K ≲ 11.5. By contrast, the DD2F_SF EoS shows immediate softening at the onset of the mixed phase. The adiabatic index for the STOS EoS softens for higher densities, but does not reach values below Γ = 1 as is the case for DD2. Γ monotonically decreases with increasing density during the mixed phase and converges towards 1.2 ≲ Γ ≲ 2 in the pure quark phase. This is significantly lower than for DD2F_SF where Γ⪆2.

Isentropes show an increase in temperature at the early onset of the mixed phase and then a decrease at higher densities |$\log \rho / (\text{g {}cm}^{-3}) \gtrsim 15.4$| in the mixed phase. In contrast, DD2F_SF featured only the latter, i.e. a decrease in temperature at constant lines of entropy. The STOS EoS exhibits similar behaviour as DD2F_SF only at temperatures above log (T/K) ≳ 11.8 and entropies s ≳ 5k_B/baryon: The transition density to the mixed-phase transitions to lower densities |$14 \lesssim \log \rho _\text{PT}/ (\text{g {}cm}^{-3}) \lesssim 14.5$| and the adiabatic index significantly softens. In this regime of the mixed phase, the EoS stiffens when quarks become the dominant degrees of freedom.

The gentle softening and the lack of an abrupt increase of Γ along the trajectories of central density and temperature are not conducive to a strong rebound after the initiation of the second collapse in the case of the STOS EoS. The soft-stiff transition due to vector repulsion in the DD2F_SF EoS seems to play a crucial role in the dynamics of the second collapse and (where applicable) the launching of an explosion. It would be interesting to consider PCSs with high initial entropy that would cross the region of low Γ at high temperature in the mixed phase to determine whether PT-driven explosions can occur in this regime,⁵ for greater detail see section 83 in Landau & Lifshitz (1980).

The violence of the collapse and rebound in the STOS model may also be reduced by the fact that the thermodynamic conditions at the centre actually do not follow the adiabats. The central entropy actually increases substantially during the phase transition. The lower and intermediate compactness STOS models ξ_2.5 ≲ 0.7 (red/orange/yellow/light green curves in Fig. 7) jump in temperature during the mixed phase at a constant density of about |$\log \rho /(\mathrm{g}\, \mathrm{cm}^{-3}) \approx 14.8$|⁠. Entropy values accordingly increase from s = 0.9k_B/baryon to 3k_B/baryon. In the higher compactness models (blue), which start with higher central entropy, the increase in temperature and entropy is less pronounced, and the post-collapse entropy is almost progenitor-independent. This phenomenon, which may appear puzzling at first glance, is due to convective mixing as we shall discuss in Section 3.2.3.

None of the low-compactness models experiences a second bounce and contract rather slowly during the mixed phase. A fast collapse follows once the progenitors pass the mixed phase region and contain a pure quark PCS core (see right thick dotted yellow line). The time between core-contraction and BH formation is roughly in between |$5\, \mathrm{ms}\lesssim t_\mathrm{collapse} \lesssim 200\, \mathrm{ms}$|⁠.

The low-compactness model z23 (red) forms an exception and instead collapses to a BH during the mixed phase with a mixed PCS core containing hadrons and quarks.

The evolution of intermediate and more massive progenitors (green/blue) with ξ_2.5 ≳ 0.7 shows a smoother central temperature evolution lacking the sudden jump in T. Furthermore, their collapse occurs faster with |$t_\mathrm{c}\sim 1\, \mathrm{ms}$| and at slightly lower densities during the mixed phase. However, similar to the lower compactness models, BH formation in those models occurs as well beyond the mixed phase and with a pure quark core. Two models in this category, z70 and z75 (thick dashed line), experience a second core bounce and a second neutrino burst.

3.2.2 CMF EoS

We did not obtain any explosions using the CMF EoS, and none of the models exhibits a second core bounce. Again the trajectories of central density and temperature in the phase diagram (Fig. 8) provide clues for this behaviour: The adiabatic index in the CMF phase diagram in Fig. 8 shows a smooth crossover to quarks. The onset density of the mixed phase (appearance of quarks) at |$\log \rho /(\mathrm{g}\, \mathrm{cm}^{-3}) \approx 15$| at fixed electron fraction Y_e = 0.25 is higher than for DD2F_SF (⁠|$\log \rho _\mathrm{onset\, MP} \approx 14.95$| log g cm⁻³) and STOS (⁠|$\log \rho _\mathrm{onset\, MP} \approx 14.5$| log g cm⁻³). Similar to DD2F_SF and STOS, the transition density to quarks shifts to lower density at higher temperatures. The crossover region shows a decrease in temperature with increasing density along isentropes and the change in adiabatic index becomes less pronounced at higher temperatures, which is also an effect of the appearance of mesons in the hadronic and gluons in the deconfined phase. Along the actual trajectories, the adiabatic index decreases gently like for the STOS EoS without substantial stiffening after the phase transition. The width of the crossover region at log (T/K) ≲ 11.8 is significantly shorter than the mixed phases in DD2F_SF and STOS. The evolution of the progenitors shows four distinct types: (1) Low-compactness progenitors ξ_2.5 ≲ 0.3 (red) start to cool after |$t\gt 6\, \mathrm{s}$| and before they reach the onset of the crossover region, and simply do not undergo a second collapse. The CMF EoS leads to a PCS structure that generally requires a larger PCS mass to reach the phase transition, especially for low-compactness models (see the third column for M_PCS in Tables 1 and 3). Therefore, the time to the phase transition is generally longer for the CMF EoS than for DD2F_SF and STOS, which allows neutrino cooling of the PCS to become relevant before the phase transition is reached. (2) Progenitors with slightly higher compactness parameters 0.3 ≲ ξ_2.5 ≲ 0.4 (orange) show a small decrease in central temperature along isentropes during the crossover region, i.e. they cool while they move smoothly through the phase transition. (3) The intermediate- and high-compactness models with ξ_2.5 ≳ 0.5 resemble the evolution of progenitors using the STOS EoS. The Type 3) models are heated at the centre during the PT by inverted convection like the STOS models. Type 3a) ξ_2.5 ≳ 0.5 exhibit a jump in temperature while 3b) show a smoother temperature increase. The collapse time of 3a) and 3b) is similar, contrary to the STOS models where high compactness models with the shape of 3b) collapsed significantly faster and at lower densities.

Figure 8.

Phase diagram with trajectories of the central density and temperature (thick solid lines) of selected CMF models. The background displays the colour-coded adiabatic index at fixed electron fraction Y_e = 0.25 for the CMF EoS. The black dashed lines are isentropes for different entropy values (indicated as contour labels).

3.2.3 Inverted convection and character of the phase transition

The increase (or even jump) in central entropy during the mixed phase sets the STOS and CMF models apart from the DD2F_SF models and may act against a precipitous collapse of the PCS core and a strong rebound. This temperature increase in the STOS and CMF models during the mixed phase is a consequence of ‘inverted convection’ inside the PCS during the phase transition, a phenomenon predicted by Yudin et al. (2016). As Yudin et al. (2016) showed, the (non-relativistic) Ledoux criterion for convective instability

$$\begin{eqnarray} C_\mathrm{L} = \frac{\rm{d}{\rho}}{\rm{d}{\it r}} - \frac{1}{{\it c}_{\it s}^2} \frac{\rm{d}{\it P}}{\rm{d}{\it r}} \gt 0, \end{eqnarray}$$

(6)

can be expressed as

$$\begin{eqnarray} C_\text{L} = -\frac{c^2 T \rho }{P\Gamma c_V} \left(\frac{\partial P}{\partial T}\right)_{\rho ,Y_\mathrm{lep}} \frac{\rm{d}{\it s}}{\rm{d}{\it r}}+ \left(\frac{\partial \rho }{\partial Y_\mathrm{lep}} \right)_{P\!,s} \frac{\rm{d}{\it Y}_{\rm lep}}{\rm{d}{\it r}} \gt 0. \end{eqnarray}$$

(7)

In terms of radial derivatives of entropy s and lepton number Y_lep with the positive heat capacity c_V = T(∂S/∂T)_ρ. Since |$(\partial P/\partial T)_{\rho ,Y_\mathrm{lep}}$| is positive under normal conditions, a negative entropy gradient tends to be destabilizing, but during a phase transition |$(\partial P/\partial T)_{\rho ,Y_\mathrm{lep}}$| can switch sign so that a positive entropy gradient (as usually encountered in the PCS core) acts as destabilizing instead. Qualitatively, the destabilization during a phase transition can also be understood from equation (6). A higher compressibility of mixed-phase matter increases the compactness ∼M/R of the PCS, leading to a steeper pressure gradient and an enlarged destabilizing term |$-c_\mathrm{s}^{-2} \partial P/ \partial r$|⁠. In the relativistic case (Yudin et al. 2016), the situation is qualitatively similar, but additional derivatives appear, and the onset of inverted convection is not strictly limited to |$(\partial P/\partial T)_{\rho ,Y_\mathrm{lep}}\lt 0$|⁠. Furthermore, stabilization or destabilization by lepton number also needs to be taken into account so that |$(\partial P/\partial T)_{\rho ,Y_\mathrm{lep}}\lt 0$| is only an approximate condition for the onset of inverted convection, and the relativistic Ledoux criterion (Thorne 1969) needs to be evaluated directly

$$\begin{eqnarray} C_\text{L} = \frac{\rm{d} {\rho \left(1+\epsilon /{\it c}^2 \right)}}{\rm{d}{\it r}} - \frac{1}{{\it c}_{\it s}^2} \frac{\rm{d}{\it P}}{\rm{d}{\it r}} \gt 0 \end{eqnarray}$$

(8)

with the speed of sound |$c_{\rm s}^2 = \Gamma P / (\rho [1 + \epsilon /c^2 + P/(\rho c^2)])$|⁠. Fig. 9 shows C_L for model z55 with the STOS EoS at selected times around the time of the second collapse along with entropy profiles. As the PCS enters the mixed phase, a wide unstable region appears around |$388\, \mathrm{ms}$| after the first bounce, and a steep entropy gradient at the edge of the low-entropy core is largely eliminated. At |$388\, \mathrm{ms}$|⁠, a small positive entropy gradient remains in the region where equation (8) indicates instability, hinting at the underlying phenomenon of inverted convection. Fig. 9 also shows that the mixing of the PCS core still takes a few milliseconds, i.e. requires several dynamical time-scales

Figure 9.

Relativistic Ledoux criterion (equation 8, top panel) and specific entropy (bottom panel) as function of enclosed mass around the time of the second collapse for STOS model z55. At |$386\, \mathrm{ms}$| after the first bounce, the entropy profile still exhibits the typical post-bounce shape with a positive gradient out to |${\sim }0.7\,\text{M}_\odot$| and a flat gradient in the protoneutron star convection zone further outside. The Ledoux criterion shows a large negative value at the low-entropy core. Another |$2\, \mathrm{ms}$| later, the Ledoux criterion has become positive in most of the formerly stable region inside |${\sim }0.7\,\text{M}_\odot$|⁠, i.e. large parts of the PCS have entered the regime of inverted convection. The entropy gradient has been flattened considerably, but is still positive, and the innermost part of the core is not yet fully mixed. At |$390\, \mathrm{ms}$|⁠, the innermost |${\sim }0.9\,\text{M}_\odot$| of the PCS has formed on a big convection zone with a flat entropy profile.

There is in fact a deeper connection between the properties of the phase transition and the onset of inverted convection (see also Yudin et al. 2016; Hempel et al. 2017). Iosilevskiy (2015) classified first-order PTs into two categories based on the sign of the change in enthalpy or second-order partial derivative of the Gibbs free energy in the mixed phase, i.e.

$$\begin{eqnarray} \Delta H \lt 0 &\rightarrow \left(\frac{\partial P}{\partial T}\right)_{{s}}\, \gt 0 \, \, \, \text{(enthalpic)} \end{eqnarray}$$

(9)

$$\begin{eqnarray} \Delta H \gt 0 &\rightarrow \left(\frac{\partial P}{\partial T}\right)_{{s}}\, \lt 0\, \, \, \text{(entropic)}, \end{eqnarray}$$

(10)

Here, ΔH = T(S_Quarks − S_Hadrons) is the enthalpy difference between the two phases before and after the phase transition. The underlying entropic nature of the hadron-quark phase transition thus can be understood as a consequence of higher entropy in the quark phase.⁶ When hadrons are deconfined into their constituents under adiabatic compression, the degrees of freedom become larger thus the kinetic energy per particle becomes smaller. The behaviour of entropic and enthalpic mixed phases is illustrated in Fig. 10.

Figure 10.

Schematic phase diagram with temperature as function of density for an enthalpic (orange) and entropic (green) EoS. All EoS in our setup feature entropic behaviour during the mixed phase from hadrons to quarks. A mixed phase is entropic if the quark phase has higher entropy after isothermal compression (see grey dotted line). This implies a decrease in temperature for adiabatic compressions (a) → (c). For enthalpic, PT the entropy decreases in the higher phase, which implies an increase in temperature for adiabtic compression (a) → (b).

Furthermore, the negative gradient of adiabats in the P−T plane can be related to the negativity of several other thermodynamic cross derivatives by simple Maxwell relations and particularly to the effect of ‘thermal softening’ (see, e.g. Steiner, Prakash & Lattimer 2000; Nakazato, Sumiyoshi & Yamada 2010; Iosilevskiy 2015; Hempel et al. 2013, 2017),

$$\begin{eqnarray} \left(\frac{\partial T}{\partial \rho}\right)_{s} \le 0 \Leftrightarrow \left(\frac{\partial P}{\partial T}\right)_{\rho } \le 0. \end{eqnarray}$$

(11)

The phase diagrams of all three EoS in Figs 6–8 show a decrease in temperature during adiabatic compression in some parts of the mixed phase and crossover region, and the effect is most pronounced for the DD2F_SF EoS. It is therefore somewhat puzzling that the DD2F_SF models do not show a similar amount of central heating by inverted convection as the STOS and CMF models. This unexpected finding can be resolved, however, by noting that the DD2F_SF models still show some drift to higher core entropy (Fig. 6) during the phase transition, and that the impact of inverted convection also depends on the extent of unstable regions and the available time for mixing.

As discussed before, most of the STOS and CMF models do not evolve through the mixed phase rapidly, in contrast to the DD2F_SF models. This is because the strong entropic phase transition in DD2F_SF is also linked to the thermal softening (equation 11) of matter, i.e. the decrease of Γ, which is related to the behaviour of ∂P/∂T and ∂P/∂S. Entropic phase transitions, therefore, leads to a softer EoS region and is linked to the stability of stars (Steiner et al. 2001). The different dynamics of the DD2F_SF models as opposed to the STOS and CMF models could therefore be interpreted as follows: For a strong entropic phase transition that requires considerable latent heat to deconfine the quarks, significant softening leads to a rapid collapse on a dynamical time-scale that leaves no time for significant mixing by convection and potentially results in a strong bounce. This hypothesis is substantiated by our results for high-compactness progenitors in the STOS setup (blue/green in Fig. 7) showing faster collapse times and lesser mixing in their cores (two of those models exhibiting a second core bounce and neutrino burst). On the other hand, for an entropic phase transition with small latent, the PCS can traverse the mixed phase more slowly, leaving a substantial fraction of its core concurrently in the unstable regime of the mixed phase for inverted convection to become effective, and mitigating the violence of the collapse and a (potential) rebound. We note in passing that the high specific entropies reached in the PCS due to inverted mixing are significantly larger than those reached in a recent study on simulating heavy-ion collisions and binary NS mergers using the same underlying CMF EoS we used (Most et al. 2022).

3.3 Neutrino signals

The neutrino signals of the exploding and non-exploding models show a variety of behaviours, including some features that have not yet been observed in models of PT-driven explosions.

We first consider the two exploding models s16 and z14 using the DD2F_SF EoS. Neutrino luminosities and mean energies for these two models are shown in Fig. 11. Immediately after the second collapse, the neutrino emission conforms to the behaviour known from previous works on PT-driven explosions. Once the shock wave reaches the close-by neutrino sphere, a second burst is released, which is dominated by electron antineutrinos (Sagert et al. 2009; Dasgupta et al. 2010; Fischer et al. 2010a, 2012). After the second neutrino burst, heavy flavour neutrinos exhibit higher luminosity than electron flavour neutrinos due to the quenching of accretion luminosity after successful shock revival. The lack of accretion luminosity is also responsible for the drop in the mean energy of ν_e and |$\bar{\nu }_\mathrm{e}$| compared to the phase prior to the second collapse.

Figure 11.

Neutrino luminosity and neutrino mean energy for eight-selected progenitors as function of time after the first bounce. Vertical dashed lines show onset of mixed phase (dashed red) and collapse (dashed grey). The colour indicates the neutrino flavour. The first and second row shows neutrino signals for the two exploding models z14 and s16 using DD2F_SF (indicated by a blue box). z15 and z18 belong to the class of failed explosion where a second neurino burst is emitted without an explosion following. Furthermore, z15 did not form a BH within the simulation time. The third row shows output for two progenitors in the STOS-B145 setup (indicated by a red box). The two models z70 and z75 exhibit a PT-driven second neutrino burst before they collapse to a BH. The luminosities of these two models show a lower peak luminosity compared to the DD2F_SF luminosities above. The small inserts in the bottom right of the panels for models z70 and z75 magnify the last 2 and |$1.2\, \mathrm{ms}$|⁠, respectively, before BH formation. The y-axis in the inserts of Row 3 ranges from 52.5–|$53.2\, \mathrm{erg}{\,\rm s}^{-1}$| (logarithmic scale); in Row 4, the y-axis of the insert ranges from 14 to |$31.9\, \mathrm{MeV}$| (linear scale). The fourth row shows the two progenitors z21 and z22 in the CMF setup (indicated by a green box). Both models collapse into BHs on a short time-scale without experiencing a second core bounce or second neutrino burst. The rise in the neutrino signal is due to contraction of the neutrinosphere during collapse.

Later on, the luminosities ν_e and |$\bar{\nu }_\mathrm{e}$| show intermittent blips of enhanced neutrino emission. The temporal pattern of these blips appears chaotic, and they could be superficially dismissed as glitches in the neutrino transport solver. These are not unphysical, however, but related to the rather tepid nature of the explosions. Due to relatively slow expansion of the forward shock, the (reverse) termination shock of the neutrino-driven wind from the PCS forms at small radii (see, e.g. Panel k in Fig. 5). The reverse shock quickly starts to propagate backwards, and eventually reaches the PCS. The dense shell ahead of the reverse shock is accreted in the process, which gives rise to a transient enhancement of the emission of ν_e and |$\bar{\nu }_\mathrm{e}$| as accretion luminosity. The accretion is too weak to permanently stifle the outflow, and the neutrino-driven wind is reestablished, which again quenches the accretion luminosity. After a while, a reverse shock is formed again, and the same process repeats a few times at irregular intervals. If observed, such an irregular enhancement of the electron flavour luminosity after a second burst would thus help to differentiate the energetics of PT-driven explosions.

For models that undergo a second bounce, but fail to explode, we find a similar variety of behaviours as Zha et al. (2021). Only a fraction of these also emit a second neutrino burst; in most cases, the shock does not make it to the neutrinosphere. For the DD2F_SF EoS, we find a second bounce for the majority of progenitors with the exception of z20, z21, and z22. Out of all the DD2F_SF models with a second bounce, however, only the models z15, z18, and z19 emit a second neutrino burst. In the STOS-setup only, two progenitors with a high compactness parameter undergo a second core bounce. Both of these models emit a second neutrino burst dominated by |$L_{\overline{\nu }_\mathrm{e}}$|⁠. On the other hand, the CMF EoS does not lead to any second core bounce and consequently no second neutrino burst either.

Fig. 11 shows neutrino luminosities and mean energies for some of these non-exploding cases. The failed explosion models z15, z18, and z19 (DD2F_SF EoS) are characterized by weaker second bursts than the exploding models z14 and s16. In z15, the luminosities of ν_e and |$\bar{\nu }_\mathrm{e}$| and the and mean energies of all flavours settle at higher values than before the second collapse right after the burst due to the smaller radii and higher temperatures of the neutrinospheres. The failed explosion model z18 shows a precipitous drop of electron flavour luminosities and mean energies right after the burst due to the transient quenching of accretion and then another burst before the electron flavour luminosities and all the mean energies settle at higher values than before the second collapse. This behaviour was already seen by Zha et al. (2021) in their simulations. The second peak is stronger with |$L_{\overline{\nu }_\mathrm{e}}\approx 5\times 10^{52} \, \mathrm{erg}\, \mathrm{s}^{-1}$| compared to |$L_{\overline{\nu }_\mathrm{e}}\approx 1\times 10^{53}\, \mathrm{erg}\, \mathrm{s}^{-1}$|⁠.

The two STOS models with a second bounce (z70 and z75) both emit a second neutrino burst right before BH formation. The second neutrino burst is dominated by |$\overline{\nu }_\mathrm{e}$|⁠, followed by ν_{τ, μ}. It is not immediately obvious from Fig. 11 whether the rather small peak in |$L_{\nu _\mathrm{e}}$| and |$L_{\overline{\nu }_\mathrm{e}}$| is indeed from the shock breakout or from increasing temperature of the collapsing neutrinosphere during the collapse to BH. In the inset at the bottom right-hand for both progenitors, we therefore show the last 1.2 ms (z70) and 2 ms (z70). We find an increase (decrease) of the |$\overline{\nu }_\mathrm{e}$| (ν_e) luminosity before collapse, i.e. there is indeed a second burst from shock breakout that is not connected to BH formation per se. The slightly delayed, more sudden increase in ν_τ/μ results from the contraction of the neutrinosphere in the course of BH formation. This effect also explain why (different from the DD2F_SF models) the heavy flavour neutrinos reach higher mean energies than the electron flavour neutrinos during the second burst.

The four bottom right-hand plots show the progenitors z21 and z22 in the CMF setup. Those progenitors do not undergo a second core bounce and immediately form BHs once they enter the crossover region. The slight increase in the neutrino luminosity (and mean neutrino energy) is due to the contraction of the neutrino-sphere, which heats the material and increases neutrino emission rates. The hierarchy is similar to DD2.

3.4 Nucleosynthesis in the ejecta

CCSNe are one of the main production sites for heavy elements in the universe. There are, however, still many open questions about the contribution of CCSNe to heavy elements beyond the iron group. The long-standing notion that CCSNe are the dominant source of rapid neutron-capture process (r-process) elements has been revised in recent years both because modern simulations neither showed the requisite high entropy and low electron fraction in the neutrino-driven wind (Hüdepohl et al. 2010; Fischer et al. 2010b), and because compact binary mergers were identified as a robust site for the r-process (e.g. Freiburghaus, Rosswog & Thielemann 1999; Goriely, Bauswein & Janka 2011). It remains conceivable though, that CCSNe contributes to some r-process production, especially in low-metallicity environments (Spite & Spite 1978; Cowan et al. 1995; Hansen, Montes & Arcones 2014). Jet-driven explosions have been considered extensively as an r-process site (Nishimura et al. 2006; Winteler et al. 2012; Mösta et al. 2018; Grimmett et al. 2021). PT-driven explosions are also a possible candidate because the rapid ejection of material from the PCS surface may provide the requisite neutron-rich conditions for an r-process. Based on 1D explosion models, Fischer et al. (2020) found considerable amounts of r-process material, reaching up to and beyond the third peak. The yields from PT-driven explosions are also of interest because they may potentially provide constraints on the rates or existence of these events.

We identify ejecta based on the same criterion as for the calculation of the diagnostic explosion energy. A mass shell is ejected if the local specific binding energy and radial velocity are both positive. Fig. 12 shows the (I) radius R, (II) entropy, (III) electron fraction Y_e, (IV) density ρ, and (V) temperature for selected ejecta trajectories as functions of time after the second bounce for the exploding model s16 using the DD2F_SF EoS. The ejecta depicted in Fig. 12 cover material in the neutrino-driven wind from mass coordinate |$1.7585\, \mathrm{M}_\odot$| out to shells starting from a few thousand kilometres that were shocked early after the second bounce. The figure does not include all the ejecta that are shocked later after the second bounce. Following Nishimura et al. (2012), we classify ejected matter into:

• prompt ejecta at larger distance |$R \gtrsim 100\, \mathrm{km}$| from the PCS (orange);

• ejecta exposed to ν-heating with Y_e ≥ 0.5 which initially stall after the second core bounce and later get ejected due to neutrino heating (‘Early neutrino-processed ejecta’, yellow); and

• intermediate ejecta exposed to ν-heating with Y_e ≤ 0.5 (blue);

• wind ejecta (turquoise) (turquoise).

Figure 12.

Radial trajectories R(t), entropy s, electron fraction Y_e, density ρ, and temperature T of selected ejected mass shells as function of logarithmic time (s) after the second bounce for the exploding model s16. Matter can be distinguished into (I) prompt ejecta (orange), (II) early neutrino-processed ejecta (yellow), (III) intermediate ejecta (blue), and (IV) neutrino-driven wind material (turquoise).

Prompt ejecta (orange) are expelled by the shock wave with no contact with the neutrino-driven wind. At collapse time, they fall in to radii |$R\approx 100\, \mathrm{km}$| (see I) before they get expelled and keep their initial Y_e ≈ 0.5 from the infall phase. In some of these trajectories, Y_e decreases slightly below 0.5 due to antineutrino captures during the second burst. The early neutrino-processed ejecta (yellow band) are initially propagating outwards after the 2nd collapse, followed by a transient phase of negative velocity before they get re-accelerated by neutrino energy disposition and mechanical work by the wind material inside. Before these ejecta are hit by the second shock, they come close enough to the PCS to reach a relatively low electron fraction of 0.50 ≳ Y_e ≳ 0.18 as a consequence of electron captures. When hit by the second shock, ejection is not fast enough to conserve the low Y_e, and the electron antineutrino burst has no substantial influence on the final Y_e because it occurs when the early neutrino-processed ejecta are still at high densities and charged-current processes have not frozen out yet. Instead, the freeze-out value of Y_e is determined by the neutrino emission during the rapid decay phase of the neutrino luminosities after the second burst. At this early time, the relative difference between electron neutrino and antineutrino luminosities and mean energies is still small enough to make the ejecta proton-rich, analogous to the situation in neutrino-driven explosions (Pruet et al. 2005; Fröhlich et al. 2006).The intermediate ejecta (blue), like the early neutrino-processed ejecta, are driven outwards initially by the second core bounce before they briefly fall in once more to about |$20\, \mathrm{km} \gtrsim R \gtrsim 100\, \mathrm{km}$| (see I) and deleptonize to moderately low Y_e before they are ejected. As a consequence of a slightly increasing difference between electron neutrino and antineutrino mean energies, neutrino processing slowly raises Y_e during ejection. On the other hand, ejection is rather slow and neutrino-processing brings this ejecta component to Y_e ≈ 0.5. The Y_e and entropy in these ejecta are, due to their small amount of mass, subject to a small drift later on because of numerical diffusion and trajectory integration inaccuracies. The mass shells in the neutrino-driven wind (turquoise) are exclusively unbound by neutrino heating. The wind is characterized by moderately high entropies of up to |$120 \, k_\mathrm{B}/\mathrm{baryon}$|⁠. The excess of electron antineutrino emission over electron neutrino emission for an extended period sets favourable conditions for low Y_e in the first few seconds after the second collapse. These conditions reverse again later when the relative difference between the electron antineutrino and neutrino luminosity decreases again (see the graph for model s16 in Fig. 11 at |${\sim }3.5\, \mathrm{s}$|⁠). Hence, neutrino heating turns the neutrino-driven wind slightly proton-rich with Y_e ≈ 0.52 at later times |$t_\mathrm{2,c}\gtrsim 2\, \mathrm{s}$|⁠.

Our findings of early neutron-rich neutrino-driven ejecta for model s16 which, at later times, turns slightly proton-rich were also found (in the context of hybrid EoS) in Nishimura et al. (2012).

In Figs 13 and 14, we plot isotopic production factors for models s16 and z14, respectively. To better understand the nucleosynthesis processes that shape the yield pattern, we also consider spatially resolved information. Fig. 15 shows the abundances of decayed nuclei as function of logarithmic mass coordinate (counting from the mass cut) for s16 (upper plot) and z14 (lower plot). These clearly reflect the different nucleosynthesis regimes that shape the salient features of the overall yields.

Figure 13.

Logarithmic isotopic production factors log (X_i/X_{i, ⊙}) for model s16. Edges connect isotopes which belong to the same element. The two isotopes ⁹⁴Zr (production factor 38) and ¹¹B (production factor 40) show the highest overproduction factor. ¹³⁹La and ¹⁸⁰Ta are produced after the first core bounce and before the second core bounce by the neutrino process in the outer shells.

Figure 14.

Logarithmic isotopic production factors log (X_i/X_{i, ⊙}) for model z14. Edges connect isotopes which belong to the same element. The two isotopes ⁹⁰Zr (production factor 200) and ¹¹B (production factor 13) have the highest overproduction factors.

Figure 15.

Mass numbers A of ejected material for the exploding progenitors s16 (upper panel) and z14 (lower panel) as function of log ejected mass, colour coded by the abundance.

3.4.1 Iron group

The iron group elements around A ∼ 56 are produced in regions with Y_e ≈ 0.5 at moderate entropy s ∼ 40k_B/nuclon within the prompt ejecta (orange band in Fig. 12). They originate from mass coordinates |$10^{-4}\, \text{M}_\odot$| to about |$10^{-2}\, \text{M}_\odot$| outside the mass cut, seen as dark blue bands for both progenitors in Fig. 15. Due to the thermodynamic conditions in this region, the yields are determined by normal or α-rich freeze-out from NSE. In fact, two different regimes can clearly be distinguished in the upper plot of Fig. 15 for model s16, with the inner region showing stripe-like patterns of relatively high mass fractions above A = 56 due to α-captures. In the case of model z14 (Fig. 15), conditions around the mass coordinate 10⁻⁴ M_⊙ also permit iron-group nucleosynthesis, but light-particle capture reactions play a greater role and lead to a more significant production of trans-iron elements as we shall discuss below. We tabulate individual decayed stable iron group nuclides in Table 5. The total ejected mass of those iron group nuclides with 45 ≤ Z ≤ 56 is 3.623 × 10⁻²M_⊙ for model s16 and 4.960 × 10⁻³M_⊙ for z14.

Table 5.

Left-hand panel: absolute masses (solar mass) for (stable decayed) LEPP material for the exploding models z14 (second column) and s16 (third column). Right-hand panel: absolute masses (solar mass) for (stable decayed) iron group nuclides.

LEPP	z14	s16	Iron	z14	s16
			group
78Kr	2.145e-07	2.003e-08	55Mn	1.649e-05	2.780e-04
80Kr	5.820e-07	1.109e-07	54Fe	6.961e-05	2.313e-03
82Kr	1.883e-07	5.005e-07	56Fe	3.756e-03	2.727e-02
83Kr	1.942e-06	1.794e-06	57Fe	2.090e-04	1.351e-03
84Kr	1.261e-05	1.299e-05	58Fe	1.926e-06	5.601e-04
86Kr	2.035e-05	6.922e-06	59Co	1.916e-05	3.110e-04
85Rb	2.079e-06	2.133e-06	58Ni	2.424e-04	1.793e-03
87Rb	3.526e-06	8.081e-07	60Ni	2.472e-04	1.174e-03
84Sr	5.079e-08	2.479e-08	61Ni	2.532e-05	1.489e-04
86Sr	1.194e-07	1.325e-07	62Ni	1.434e-04	4.197e-04
87Sr	3.426e-08	8.903e-08	64Ni	4.831e-05	1.068e-04
88Sr	5.524e-05	8.335e-06	63Cu	1.366e-05	1.309e-04
89Y	5.422e-06	1.109e-06	65Cu	9.219e-06	2.692e-05
90Zr	3.086e-05	3.000e-06	64Zn	8.314e-05	1.071e-04
91Zr	1.118e-06	7.637e-07	66Zn	5.407e-05	3.851e-05
92Zr	3.449e-07	1.775e-06	67Zn	4.827e-06	7.280e-06
94Zr	2.487e-07	2.349e-06	68Zn	1.996e-05	3.832e-05
96Zr	1.693e-08	2.415e-07	70Zn	6.634e-07	8.277e-07
93Nb	8.254e-08	4.490e-07

3.4.2 Light element primary process

The process leading to elements beyond iron with A < 130 is usually referred to as light element primary process (LEPP; Travaglio et al. 2004; Arcones & Montes 2011). LEPP isotopes in the upper plot of Fig. 15 are produced in neutron-rich intermediate ejecta (blue band in Fig. 12). The two LEPP-process zirconium isotopes ⁹⁰Zr (model z14) and ⁹⁴Zr (model s16) show significant overproduction factors. Their absolute yields are tabulated along other LEPP process isotopes in Table 5.

3.4.3 Neutrino process in the envelope: ¹¹B, ¹³⁸La, and ¹⁸⁰Ta

Model s16 exhibits significant overproduction of the heavy odd-odd nuclides ¹³⁸La and ¹⁸⁰Ta which are among the rarest Solar system species (Martínez-Pinedo et al. 2017). We find overproduction factors of ∼60 for ¹⁸⁰Ta and ∼20 for ¹³⁸La for the exploding model s16 whose production is however unrelated to the PT-driven explosion mechanism.

These nuclei belong to the neutron-deficient class of nuclei, called p-nuclei (Rayet et al. 1995). One likely origin for these nuclei is the neutrino (ν)-process (Woosley et al. 1990; Heger et al. 2005), which entails charged-current or neutral-current spallation processes induced by supernova neutrinos. This process also operates in model s16 and can be identified as the source of ¹³⁸La and ¹⁸⁰Ta since both isotopes are produced after the neutrino burst from the first core bounce.⁷¹³⁸La and ¹⁸⁰Ta are mainly produced by the charged reactions ¹³⁸Ba(ν_e, e⁻)La and ¹³⁸Hf(ν_e, e⁻)¹⁸⁰Ta (Woosley et al. 1990; Heger et al. 2005). The first neutrino burst, dominated by ν_e, provides favourable conditions for these charged-current reactions.

We do not see ¹³⁸La and ¹³⁸Ta production for model z14 since the parent weak s-process isotopes ¹³⁸Ba and ¹³⁸Hf are not present in the zero metallicity model z14.

It should be noted, however, that the ¹⁸⁰Ta yield for s16 is possibly lower than what our network predicts. The network calculates the sum of isomeric and ground states which have both very distinct half-lives (the lifetime of the ground state is much shorter⁸). This leads to an overestimation of ¹⁸⁰Ta by a factor of 0.3−0.5 (Rauscher et al. 2002). Since ¹⁸⁰Ta and ¹³⁸La are both produced after the first neutrino burst and before the phase transition-induced explosion, the EoS would have little effect on their production. Additional and more detailed discussion about the interesting isotopes ¹³⁸La and ¹⁸⁰Ta can be found here (Arnould & Goriely 2003; Byelikov et al. 2007; Austin, Heger & Tur 2011; Lahkar et al. 2017; Sieverding et al. 2018).

¹¹B is another rare ν-process isotope that shows large production factors for both progenitors in Figs 13 and 14. Contrary to the heavier p-nuclei isotopes ¹³⁸La and ¹⁸⁰Ta it is mainly produced via the neutral current ¹²C(ν, ν′p)¹¹B on the rather abundant ¹²C (Laha, Beacom & Agarwalla 2014). For both progenitor models, it is produced after the second core bounce and subsequent neutrino burst.

3.4.4 Rate constraints

The overproduction factors can, in principle, be used to constrain the fraction of PT-driven explosions (of a similar type as s16 and z14) among all CCSNe by assuming that the most overproduced nucleus solely originates from this particular type of explosion (e.g. Wanajo et al. 2009). Based on Figs 13 and 14, the nuclides with the largest production factor are ⁹⁴Zr for s16 and ⁹⁰Zr for z14. We assume that there is a negligible contribution to galactic ¹⁶O production for each of these two representatives PTD-CCSN and negligible contribution from non-PT driven CCSN to the production of the two zirconium isotopes ⁹⁰Zr and ⁹⁴Zr. An approximate upper rate constraint on the fraction f (Wanajo et al. 2009, 2018; Wanajo, Janka & Mueller 2011) of PTD-CCSNe is then given by

$$\begin{eqnarray} \frac{f_\mathrm{s16/z14}}{1-f_\mathrm{s16/z14}} = \frac{X_\odot (^{A}\text{Zr})/X_\odot \left(^{16}\text{O} \right)}{M(^{A}\text{Zr})/M_\text{CCSN}\left(^{16}\text{O} \right)}, \end{eqnarray}$$

(12)

where f_s16/z14 stands for the ratio of the respective events to all CCSN-events. A denotes the mass number of the Zr isotope. The nuclear reaction network yields an absolute mass for both isotopes of M_z14(⁹⁰Zr) = 3.856 × 10⁻⁴ M_⊙ and M_s16(⁹⁴Zr) = 2.637 × 10⁻⁶ M_⊙. |$M_\text{CCSN}(^{16}\text{O})=1.5\, \text{M}_\odot$| is the average mass production of ¹⁶O by massive CCSN with masses  M ∈ [13−40] M_⊙ averaged over the stellar initial mass function (Nomoto et al. 2006). The other parameters are X_⊙(¹⁶O) = 7.377 × 10⁻³, X_⊙(⁹⁰Zr) = 1.404 × 10⁻⁸ and X_⊙(⁹⁴Zr) = 4.955 × 10⁻⁹ and taken from Lodders (2021). The fraction of these kinds of explosions under the above assumptions is only constrained to be f_s16 ≤ 0.520 for s16 and to f_z14 ≤ 0.007. Thus, because of the large variations in yields between these two PT-driven explosions, no meaningful constraint can be placed based on the frequency or existence of PT-driven explosion yet. Nucleosynthesis from individual PT-driven explosion models is clearly not sufficiently robust for this purpose. Integrating PT-driven explosions into the bigger picture of chemogalactic evolution would clearly require fine representative sampling of this hypothetical explosion channel by many models, which is not realistic at this stage given the current uncertainties on the high-density equation of state.

4 SUMMARY AND OUTLOOK

In this paper, we studied the effects of the QCD phase transition in CCSNe in order to better gauge the robustness and sensitivities of the proposed PT-driven explosion mechanism. We performed spherically symmetric CCSN simulations for up to 40 progenitors in the mass range |$14\tt {-}100\, \text{M}_\odot$| with three different hadron-quark EoS using the general relativistic CoCoNuT-FMT code with an effective 1D mixing-length treatment for convection. We considered two hybrid EoS with a first-order phase transition to quark matter: DD2F_SF contains stiffening repulsive quark interactions in the quark phase (Fischer et al. 2018; Bastian 2021); the STOS-B145 is matched to the Bag model in the quark phase (Sagert et al. 2009, 2010). The third EoS is a chiral mean-field model with a smooth crossover to quarks (CMF, Motornenko et al. 2020).

While PT-driven explosions have been proposed as a scenario for hyperenergetic supernovae from massive progenitors (Fischer et al. 2018), our results present a more nuanced view on the robustness and potential signatures of such explosions. We find only two explosions among our large set of models, namely for an |$14\hbox{-}\text{M}_\odot$| zero-metallicity progenitor and a |$16\hbox{-}\text{M}_\odot$| solar-metallicity progenitors, both using the DD2F_SF EoS. In another three progenitors models (z15, z18, z19), the shock launched by the rebound after the PT-induced collapse propagates out for several |$100\, \mathrm{km}$|⁠, but ultimately fails to explode the star. There are different parameterizations available for the DD2F_SF-EoS, which differ in various parameters for the treatment of hadronic and quark matter, such as the onset density of the mixed phase and the latent heat; these, in turn, influence the maximum cold NS mass and the mass–radius relation (see Kaltenborn et al. 2017; Bastian 2021 for a detailed discussion). The parameterization we chose has a higher onset density and slightly softer hadronic EoS compared to the one utilized in Fischer et al. (2018, 2020), Fischer (2021). For direct comparison, we simulated a small set of three progenitors using the DD2F_SF parameterizations in Fischer et al. (2018), i.e. DD2F_SF RDF 1.2, (which uses an excluded volume version of DD2F in the hadronic phase). We found very similar outcomes to DD2F_SF RDF 1.4, with the explosion of z14 having a low explosion energy |$E_\mathrm{exp} \sim 10^{50}\, \mathrm{erg}$|⁠, BH formation for z85, and a failed explosion for z15. We conclude that the general trend we see towards low explosion energies (for exclusively low compactness models) is robust within small parameter variations in the DD2_SF-EoS set for simulations with our CoCoNuT-FMT code. We leave a more systemic evaluation for the entire DD2F_SF parameter-set open for future studies. We find no successful explosions for the STOS-145 EoS, although two zero-metallicity progenitors of 70 and 75 M_⊙ undergo a rebound after the PT collapse and end up as ‘failed explosions’. None of the models using the CMF EoS explode. Some low-compactness progenitors do not reach the phase transition at all, whereas models with moderate or high compactness quietly cross the phase transition without a second bounce and eventually collapse to BHs.

The higher likelihood of PT-driven second core bounces for low compactness for the DD2F_SF EoS is consistent with trends identified by Zha et al. (2021), (who applied the STOS-B145 EoS). The phase transition is reached significantly below the maximum stable TOV mass after the phase transition for the lower specific core entropy of the PCS in low-compactness models and a direct collapse to a BH can be avoided. For a second bounce to occur, the maximum mass of the twin star branch in the mass–radius relation has to be larger than the entropy-dependent maximum mass in the large-radius branch (Zha et al. 2021).

The collapse and rebound must also be sufficiently violent in order for an explosion to occur. In the lowest compactness parameter models, the central density of the PCS only increases modestly before a rebound occurs, and the resulting shock wave is too weak to trigger an explosion. This leaves only a small window for successful explosions at low, but not excessively low compactness.

Inspecting the phase diagrams of the different EoS suggests that the structure of the phase transition for the DD2F_SF EoS, with a strong softening (as reflected by the drop in the adiabatic index Γ) during the mixed phase and stiffening in the quark phase due to repulsive vector repulsion, provides more favourable conditions for an explosion. By contrast, the STOS EoS shows gradual softening towards higher densities, and the situation for the CMF EoS is similar for the relevant PCS core entropies. In this context, it is important to note that EoS that include any significant repulsive interaction, or stiffening, in the quark phase usually violate constraints from lattice-QCD simulations (Steinheimer & Schramm 2014) as well as the expected Stefan–Boltzmann limit at high densities.

A further dynamical difference between the DD2F_SF models and the STOS and CMF models lies in the occurrence of ‘inverted convection’ (Yudin et al. 2016). During the phase transition, the positive entropy gradient in the PCS core becomes convectively unstable due to the anomalous behaviour of the thermodynamic derivative (∂P/∂T)_ρ, which can be connected to the entropic nature of the phase transition. As a result, the core of PCS is thoroughly mixed, and the centre of the PCS is heated significantly, raising the central entropy to about 4k_B/baryon. The development of inverted convection could have a further damping effect on the collapse and rebound (where applicable) in the STOS and CMF models. For the DD2F_SF EoS, conditions for inverted convection are met in principle, but due to the stronger softening in the mixed phase, the collapse and rebound occur too rapidly to allow significant mixing by convection in the PCS core.

Our neutrino signals largely show similar characteristic signatures of a first-order phase transition as in previous studies (Sagert et al. 2009; Zha et al. 2021). For successful or ‘failed explosions’ (with transient outward propagation of the shock), there is a strong antineutrino burst (Sagert et al. 2009; Dasgupta et al. 2010; Fischer et al. 2010a, 2012). Similar to Zha et al. (2021), some failed explosions are characterized by a short phase of strongly reduced neutrino emission followed by another burst when the shock stalls again and accretion on to the PCS recommences. For our two exploding models, we find a new neutrino signature, namely a series of minibursts during the explosion phase, which are the result of occasional fallback from the weak wind outflow from the PCS.

The dynamics of the two exploding models is markedly different from earlier studies (Fischer et al. 2020) that proposed PT-driven supernovae as a scenario for hyperenergetic explosions of massive progenitors. The two successful explosion models s16 and z14 only reach energies of |${\sim }1.25\times 10^{50}$| and |$3.64\times 10^{49}\, \mathrm{erg}$|⁠, respectively. The nucleosynthesis is also less spectacular than in previous studies that found significant r-process production (Nishimura et al. 2012; Fischer et al. 2020). The models show appreciable overproduction of some LEPP nuclei, in particular for model z14. The largest production factors are found for ⁹⁰Zr for model z14 and for ⁹⁴Zr for model s16. However, the differences between the two models indicate that PT-driven explosions would not be a very robust site for LEPP nucleosynthesis. For this reason, chemogalactic evolution cannot meaningfully constrain the rate of putative PT-driven explosions.

On a more fundamental level, our results cast doubts on the viability of PT-driven explosions. Since conditions for a PT-driven explosion appear most favourable in low-compactness progenitors as for our DD2F_SF models, it is likely that a neutrino-driven explosion would develop before the PCS enters the phase transition. However, the progenitor dependence seen for the DD2F_SF EoS may not be generic as indicated by the STOS models, where a second bounce does occur (albeit without triggering in explosion) for two high-compactness progenitors that likely would not explode by the neutrino-driven mechanism or at best experience an ‘aborted explosion’ with massive fallback and BH formation.

A more serious problem lies in the thermodynamic properties of the EoS and the nuclear physics assumptions that are required for successful explosions in our study. The explosions of the two DD2F_SF models hinge on the combination of significant softening during the mixed phase followed by abrupt and pronounced stiffening in the quark phase. This peculiar behaviour of the DD2F_SF EoS rests on the problematic assumption of vector repulsion in the quark phase. Although repulsive quark vector interactions have been used in the literature and fulfil the two solar mass constraints (Kaltenborn et al. 2017; Bastian 2021), QCD lattice calculations disfavour repulsion among quarks in the other regimes of QCD at high temperatures and vanishing baryon density (Steinheimer & Schramm 2011, 2014).

These problems do not yet spell the end for a dynamical role of the hadron-quark phase transition in CCSNe, however. Further work is still necessary to resolve the tensions between different results on the progenitor dependence of PT-driven explosions. But even if the phase transition does not provide a pathway to a distinct explosion mechanism in its own right, its effects on CCSN dynamics and observables still deserve more nuanced attention. In particular, the phenomenon of inverted convection needs to be investigated in multidimensional simulations. First, two-dimensional simulations of PT-induced collapse with a hybrid STOS-bag EoS (with a bag constant of |$165\, \mathrm{MeV}$|⁠) were already presented by Zha et al. (2020), who found a prominent high-frequency gravitational-wave signal during the phase around the second bounce due to convective activity deep inside the PCS. While Zha et al. (2020) associated the new gravitational-wave signature with the second collapse and bounce in the wake of a first-order phase transition, we speculate that such bursts could be more generic signatures of the hadron-quark phase transition as long as they provide conditions for inverted convection. EoS with a second-order phase transition like CMF might have analogous fingerprints in gravitational waves. While gravitational-wave signatures from the PCS core are beyond the frequency range of current gravitational-wave detectors, concepts for kilohertz detectors that could probe the relevant high-frequency regime are already being explored (Ackley et al. 2020). Future simulations will need to more thoroughly investigate the potential to leverage gravitational-wave and neutrino detection to constrain the physics of the hadron-quark phase transition.

The thermodynamic conditions obtained in extreme astrophysical events like NS mergers and CCSNe can be probed with present beam energies at heavy-ion collision laboratories such as GSI (Most et al. 2022). Due to inverted convection, the conditions for the majority of progenitors in our CCSN simulations using the CMF EoS show high specific entropies (s ∼ 4k_B/baryon), temperatures (⁠|$T \sim 70\, \mathrm{MeV}$|⁠) and densities (n ∼ 7ρ_sat) which exceed the conditions reached in hydrodynamic simulations of heavy-ion collisions at the SIS18 accelerator and NS merger simulations that utilized the same underlying CMF model (Motornenko et al. 2020; Most et al. 2022). This indicates that heavy-ion collisions with higher beam energies than currently available at GSI, e.g. at the SIS100 accelerator of the Facility for Antiproton and Ion Research, FAIR in Europe, which is presently under construction in Darmstadt, would be necessary in order to study conditions applicable to CCSNe.

ACKNOWLEDGEMENTS

PJ thanks R. Mardling for fruitful discussions. BM acknowledges support by ARC Future Fellowship FT160100035. This work is based on simulations performed within computer time allocations from Astronomy Australia Limited’s ASTAC scheme, the National Computational Merit Allocation Scheme (NCMAS), and an Australasian Leadership Computing Grant on the NCI NF supercomputer Gadi. This research was supported by resources provided by the Pawsey Supercomputing Centre, with funding from the Australian Government and the Government of Western Australia. This research was supported, in part, by the Australian Research Council (ARC) Centre of Excellence (CoE) for Gravitational Wave Discovery (OzGrav) through project number CE170100004, by the ARC CoE for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) through project number CE170100013, and by the National Science Foundation under Grant No. PHY-1430152 (JINA Center for the Evolution of the Elements, JINA-CEE). JS thanks the Samson AG for support. AM acknowledges the Stern-Gerlach Postdoctoral fellowship of the Stiftung Polytechnische Gesellschaft. HS acknowledges the Walter Greiner Gesellschaft zur Förderung der physikalischen Grundlagenforschung e.V. through the Judah M. Eisenberg Laureatus Chair at Goethe Universität.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the authors, subject to considerations of intellectual property law.

Footnotes

1

Depending on the specific entropy per Baryon in the core

2

In the mixed phase and quark phase, ‘entropy per baryon’ is to be understood as the entropy of a fluid element divided by its net baryon number, which can still be defined even when quarks are not bound in baryons anymore.

3

Note that the electron fraction in the phase diagram background is fixed to Y_e = 0.25, so the onset of the mixed phase cannot be pinpointed exactly using the isocontours of the quark fraction.

4

As a consequence of unusual thermodynamic properties which we will discuss in more depth in Section 3.2.3

5

Note, however, that the two-EoS description in close vicinity to the phase construction, that is where both EoS intersect, lacks reliability, since either quark or hadronic models lie beyond their regime of validity (Baym et al. 2018b). Furthermore, the two-EoS approach does not admit a continuous phase transformation and therefore no QCD-critical end-point in the QCD phase diagram (Hempel et al. 2013)

6

The Clausius–Clapeyron relation (Landau & Lifshitz 1980) is often referenced for the negative pressure gradient |$\frac {{\rm d}P}{{\rm d}T}$| for Maxwellian phase transitions (local charge neutrality during the mixed phase) where the pressure in the mixed phase is independent of density which is not the case for a Gibbs construction.

7

The p-process (also called γ-process) can be excluded since no shock has passed through the material at this point (Kusakabe, Iwamoto & Nomoto 2010).

8

An isomeric state is an exited meta-stable (⁠|$\tau \le 1\,$|ns than the ground state) nuclear state (Misch et al. 2021).