© The Authors 2025

1. Introduction

Massive stars with masses of greater than 10 M_⊙ are rare compared to their low-mass counterparts, but their contribution to the evolution and chemical enrichment of galaxies is disproportionately important (Hopkins 2014). Their strong winds and supernova (SN) explosions create mechanical and energetic feedback and enrich their surroundings with heavier elements, determining the properties and evolution of galaxies and of the next generations of stars (Geen et al. 2023). However, not every massive star explodes; a fraction of them are expected to collapse directly and form black holes (BHs). Core-collapse SNe are not always successful and failed explosions are another common path for BH formation (see, e.g., the recent review by Heger et al. 2023). Even successful explosions can be associated with BH formation through fallback accretion (Burrows et al. 2023). Understanding and predicting the fate of stars remains one of the main unsolved problems in astrophysics.

The fate of stars is determined by the explosion mechanism, fallback dynamics, and their final core structure. In this work, we focus on the latter. Progress in stellar and SN physics in recent years has shown that there is no simple initial mass threshold for a star to form a BH. Moreover, almost 30 years ago, Timmes et al. (1996) pointed out that the final structure of a star, characterized for example by the final iron core mass, is not monotonic with mass. This picture was confirmed by multiple studies (e.g., Brown et al. 2001; O’Connor & Ott 2011; Sukhbold & Woosley 2014; Pejcha & Thompson 2015; Ertl et al. 2016; Sukhbold et al. 2018; Limongi & Chieffi 2018; Schneider et al. 2021; Temaj et al. 2024). Moreover, it is now established that most massive stars live in close binary or multiple systems (e.g., Sana et al. 2012), further complicating this picture. Binary interactions have been shown to affect the precollapse core structures of stars, both based on studies of pure He star models approximating stripped stars in binaries (Brown et al. 2001; Woosley 2019; Aguilera-Dena et al. 2022, 2023) and binary evolution models of stripped stars, accretors, and, mergers (Laplace et al. 2021; Schneider et al. 2021, 2023, 2024). In turn, this affects their explodability (Müller et al. 2019; Vartanyan et al. 2021; Woosley et al. 2020; Antoniadis et al. 2022), nucleosynthesis (Farmer et al. 2021, 2023), significantly reduces the parameter space for the formation of compact-object mergers observable with gravitational-wave (GW) observations (Schneider et al. 2021, 2023), and leads to features in the chirp-mass distribution of binary BH mergers (Schneider et al. 2023). However, to study the properties of BH populations and the formation of GW sources, current state-of-the-art studies necessarily have to make simplifying assumptions regarding the formation of BHs, often using analytical prescriptions solely based on the core mass of single-star progenitors (e.g., Fryer et al. 2012). These can lead to substantially different outcomes compared to models that take the structure of stars into account (Patton et al. 2022).

Observationally, there exist only a few direct hints of the link between the pre-SN structure of massive stellar progenitors and the formation of BHs. The red supergiant (RSG), N6946-BH1, which was observed to suddenly vanish in the optical after a short outburst (Gerke et al. 2015; Adams et al. 2017; Sukhbold & Adams 2020; Basinger et al. 2021), could be the first direct progenitor of a BH ever observed. This event is compatible with model predictions for a single star with high compactness undergoing a failed SN explosion and eventually forming a BH (Lovegrove & Woosley 2013; Sukhbold & Adams 2020; Temaj et al. 2024). Very recent observations with JWST potentially challenge this interpretation by identifying an infrared source at the location of this object, which may correspond to a surviving star enshrouded by dust (Beasor et al. 2024) or to the emission from a failed SN (Kochanek et al. 2024). Future observations are needed to better understand whether or not this event was indeed a RSG forming a BH after a failed explosion.

An additional, indirect observational clue as to the link between RSGs and BH formation comes from observations of hydrogen-rich (type II) core-collapse SNe. Archival data searches have unambiguously identified several RSG progenitors at the location of these SNe. These observed RSG SN progenitors tend to have low luminosities of log L/L_⊙ ≤ 5.1 (Smartt 2009). This is in tension with the observed maximum luminosity of about log L/L_⊙ ≊ 5.5 found for RSGs in the galaxy and in the Magellanic Clouds (Davies & Beasor 2020), and is known as the missing RSG problem. A possible interpretation for these “missing” luminous RSG SN progenitors is that these are high-mass stars that “quietly” form BHs instead of exploding. The exact value of this maximum luminosity remains to be determined because of systematic uncertainties associated with the photometric data available for these objects (Davies et al. 2018; Davies & Beasor 2020). Nonetheless, these observations offer important insights into the structures of stars that may be the progenitors of BHs. In principle, the bolometric pre-SN luminosity of a RSG can be directly linked to its final core mass, independently of the uncertainties in convective boundary mixing or rotation (Temaj et al. 2024)¹. However, the initial mass of these progenitors is very uncertain. This is because, even assuming single-star evolution, variations in internal mixing, mass-loss history, and rotation all affect the relation between the initial and final core mass of stars (Farrell et al. 2020). Taking binary evolution into account further complicates this relation (Zapartas et al. 2019, 2021). Generally, the connection between the final core mass and final observable properties of stars is much better constrained (Temaj et al. 2024). Finally, the lack of stars with high core masses exploding as SNe is also supported by studies of their late-time SN spectra (e.g., Jerkstrand et al. 2012) and by age-dating of supernova remnant (SNR) environments (e.g, Jennings et al. 2014).

Based on detailed stellar models, several summarizing quantities have been defined to evaluate the fate of massive stars. The compactness parameter ξ_m (O’Connor & Ott 2011) is commonly used in recent literature and is defined as

$$\begin{array}{c}\hfill {\xi }_m=\frac{m/{\mathrm{M}}_{\odot }}{R(m)/1000\mathrm{km}},\end{array}$$(1)

where m is the mass coordinate at which the compactness is evaluated, typically at a chosen value of 2.5 M_⊙, and R is the radius at this mass coordinate. Essentially, it is a measure of the density (or mass–radius relation Chieffi & Limongi 2020) outside the iron-rich core. This quantity, though arguably arbitrary in its definition, is known to correlate with other key properties, such as the iron-core mass, and the binding energy above the iron-rich core (e.g., Sukhbold & Woosley 2014; Schneider et al. 2021; Temaj et al. 2024, see also Fig. 1). Stars with large iron core masses tend to have a high binding energy outside this core, and are thus difficult to explode using any explosion mechanism and tend to ultimately form BHs (Brown et al. 2001; Sukhbold & Woosley 2014; Heger et al. 2023; Temaj et al. 2024). In recent literature, the compactness parameter has been used as a predictor for the final remnant expected after core collapse, with high values indicating BH- and low values NS formation (O’Connor & Ott 2011; Ugliano et al. 2012; Sukhbold & Woosley 2014; Limongi & Chieffi 2018; Schneider et al. 2021, 2023, 2024; Heger et al. 2023). Multiple studies have pointed out that more sophisticated metrics are needed to accurately capture the explosion physics and understand the conditions for shock revival (Pejcha & Thompson 2015; Ertl et al. 2016; Müller et al. 2016; Sukhbold et al. 2016; Vartanyan et al. 2021; Burrows & Vartanyan 2021). The explodability of stars is a subject of active discussion in the community and several explodability criteria have been proposed and explored in recent years. These include the two-parameter criterion of Ertl et al. (2016), the presence of a steep density profile with a density discontinuity around the Si/O interface (Vartanyan et al. 2021; Tsang et al. 2022; Wang et al. 2022; Boccioli et al. 2023), or a forced explosion condition (Murphy & Dolence 2017; Gogilashvili et al. 2023).

Independently of the discussion about explodability criteria, several unrelated groups making different assumptions regarding the microphysics and using different methods have reported remarkably similar patterns in the final core structure of stars (often summarized by the compactness parameter) as a function of their core or initial mass (O’Connor & Ott 2011; Sukhbold & Woosley 2014; Sukhbold et al. 2018; Limongi & Chieffi 2018; Chieffi & Limongi 2020; Chieffi et al. 2021; Schneider et al. 2021; Patton & Sukhbold 2020; Takahashi et al. 2023; Temaj et al. 2024). Typically, these patterns consist of two prominent peaks in the compactness parameter, separated by ≈15 M_⊙ in initial mass and ≈7 M_⊙ in CO core mass. The robustness of this pattern points to a common underlying physical process determining the final core structure of stars, and with it, their fate.

In a pioneering study, Sukhbold & Woosley (2014) analyzed the pattern in the compactness parameter for the first time in great detail and linked it to the different nuclear burning conditions in the last evolutionary stages of massive stars post-helium burning. Patton & Sukhbold (2020) demonstrated that variations in the final compactness are linked to the initial conditions for core carbon burning, that is, the mass of the CO core and the initial central carbon abundance at core helium exhaustion. Sukhbold et al. (2018) found that the final compactness of a star is influenced by small variations in physical assumptions and resolution, and interpreted this as a sign of intrinsic randomness in the core structure of stars. However, Chieffi & Limongi (2020) argued that these apparently random variations can be traced back to their assumptions regarding the core helium-burning evolution, where semi-convection in particular can result in late ingestion of helium in the core, generating “breathing pulses”, which change the central carbon abundance and lead to different initial conditions for core carbon burning and ultimately to a different final core structure. In our work, which includes convective boundary mixing above the helium-burning core, we do not find these signs of intrinsic stochasticity either (Schneider et al. 2023; Temaj et al. 2024).

Independent studies connected variations in the final core structure of stars to the number and size of carbon-burning shells and to the transition from convective to radiative carbon burning (Brown et al. 2001; Sukhbold & Woosley 2014; Sukhbold & Adams 2020; Chieffi & Limongi 2020). However, this explanation appears incomplete. Even after the transition from convective to radiative core carbon burning, models of stars that undergo radiative core carbon burning can result in low compactness (see, e.g., Fig. 2 of Sukhbold & Adams 2020). In addition, though the cause of the prominent drop in compactness after the first peak remains unclear, Sukhbold & Woosley (2014) identified a link between the base of the carbon shell exceeding the effective Chandrasekhar mass and a smaller oxygen-burning core.

Schneider et al. (2021) identified a connection between models with high compactness and the mass range for which carbon burning and neon burning become neutrino-dominated. Following these findings, in the present study, we investigate the origin of the observed patterns in the final structures of massive stars. We compute detailed simulations of the common progenitors of core-collapse events, namely massive (17–50 M_⊙), single, nonrotating stars at solar metallicity, and focus on the evolution of the innermost 6 M_⊙.

We present our computational setup in Sect. 2 and the overall properties of our models in Sect. 3. In Sect. 4, we conduct a simplified experiment to identify and understand the general physical mechanisms responsible for the observed trends in final core structure. These insights are then applied to our fiducial set of stellar models in Sect. 5. We summarize the physical mechanisms identified as being responsible for determining the main pattern in the final core structure of stars in Sect. 6. We discuss the implications and uncertainties of our results in Sect. 7 and present our conclusions in Sect. 8.

2. Methods

We compute the interior structure of massive single stars with initial masses between 17 and 50 M_⊙ with the MESA stellar evolution code (version 10398, Paxton et al. 2011, 2013, 2015, 2018). Our models build upon Schneider et al. (2021, 2023), with similar assumptions. Specifically, our models are computed at solar metallicity Z = 0.0142 (Asplund et al. 2009) and are nonrotating. We compute convective mixing with an mixing-length theory (Böhm-Vitense 1958) parameter of α_MLT = 1.8 and assume step overshooting of 0.2 pressure-scale height, which is only applied over the H and He-burning convective cores. We adopt the Ledoux criterion for convection and assume a semi-convection efficiency of α_sc = 1.0 (Schootemeijer et al. 2019). We adopt the rescaled “Dutch” wind mass loss rates of Schneider et al. (2021) and enable the MLT++ method of MESA that boosts the local energy transport for outer layers of our massive star models that locally exceed the Eddington limit. The models are computed with the MESA approx21_cr60_plus_co56.net nuclear network until the onset of core collapse, which is defined as the moment when the infall velocity of the iron-rich core exceeds 900 km s⁻¹. This network effectively sets Y_e in the entire iron-rich core by making the approximation that deleptonizations in this core only occur through electron-captures onto ⁵⁶Fe. However, variations in Y_e between different stellar models are generally small (Woosley et al. 2002) and do not affect the formation of the main qualitative patterns in the final core structure that are the subject of this work, which already appear at the end of core Ne burning (see also Sect. 3.3). Using this nuclear reaction network is sufficient for our purpose, but we caution against employing these models as input for predictions of three-dimensional SN simulations or nucleosynthesis yields, which require a larger nuclear network (Farmer et al. 2016; Renzo et al. 2024). To further verify that the final structure pattern is reproduced in models with larger networks, we compute three additional models with initial masses of 21 M_⊙, 22 M_⊙, and 23 M_⊙, for which we employ a nuclear reaction network of 128 isotopes (as recommended by Farmer et al. 2016). The 21 M_⊙ and 22 M_⊙ models encounter numerical difficulties after the iron-core infall velocity exceeds 250 km s⁻¹. By this point the iron core mass and the central entropy change only slightly (by less than 2%), so we consider these to be comparable to our default collapse models. The reaction rates in our models are based on the JINA REACLIB database version 2.2 (Cyburt et al. 2010).

We ensure a high spatial resolution in our models, in particular in zones of high temperature and density (Farmer et al. 2016). More specifically, we adopt a minimum of 2000 grid points and a grid spacing option mesh_delta_coeff = 0.6, which results in an average of 5000–6000 grid points for each model. We also ensure a high temporal resolution throughout the evolution with a maximum time step of 10⁻⁴ years that is further limited based on changes in composition.

To better disentangle the effects responsible for determining the final structure of massive stars, we perform a controlled experiment (see Sect. 4). These models have the same core masses but a different central carbon mass fraction at the moment of core helium depletion (when the central helium mass fraction is lower than 10⁻⁴). This is achieved by computing additional models using the MESA relax_initial_composition method in which we artificially modify the central ¹²C abundance and the central ¹⁶O abundance while keeping the total mass fractions constant. Our base model for this experiment is our fiducial model with an initial mass of 22 M_⊙, which corresponds to the first compactness peak, at the moment of core helium depletion.

We examine the effect of changing the ¹²C(α, γ)¹⁶O reaction rate and discuss these in Appendix A. This notoriously uncertain nuclear reaction plays a crucial role for the evolution of stars, including their fate (Weaver & Woosley 1993; Austin et al. 2014; Sukhbold & Adams 2020; Farmer et al. 2020), as it determines the final central abundances at the end core helium depletion and the mass of the CO core (see also Sect. 3.3). Our fiducial set of models adopts the default MESA rate from Xu et al. (2013). The other two sets we compute adopt the approximately 15% lower rate from Kunz et al. (2002) which is often used in the literature, and a rate that is 10% higher than that of our default model. We compute additional sets of models for which we vary the semi-convection efficiency, discussed in detail in Appendix B. Finally, we perform a resolution test in Appendix C to explore the effect of numerical uncertainties on the occurrence of shell mergers.

To evaluate the final fate of our massive star models, we compute the expected explosion outcome of our models using the semi-analytical parametric neutrino-driven SN explosion model of Müller et al. (2016), with the same assumptions as in Schneider et al. (2021, 2023, 2024) and Temaj et al. (2024). For simplicity, we do not consider BHs formed by SN fallback in this work.

For comparison, we also employ the two-parameter explodability criterion by Ertl et al. (2016). This criterion depends on the mass M₄, which is mass coordinate m where the entropy reaches a value of s/(N_Ak_B) = 4 and on ${\mu }_4=\frac{\mathrm{\Delta }m/{\mathrm{M}}_{\odot }}{\mathrm{\Delta }r/1000\mathrm{km}}{|}_{s=4}$, the radial mass gradient at M₄. Physically, this location M₄ typically corresponds to the Si/O interface (i.e., the transition point between the Si-rich and O-rich layers), which has been found to be a good predictor for successful multi-dimensional neutrino-driven SN explosions (Ertl et al. 2016; Müller et al. 2016; Ertl et al. 2020). We use the s19.8 calibration of Ertl et al. (2016) to models by Woosley et al. (2002) to distinguish successful and failed explosions.

3. Final stellar properties

3.1. Final stellar structure

We characterize the final stellar structure of our models by several quantities, including the compactness parameter, final specific central entropy, iron-core mass, and final binding energy (Fig. 1). Here, the binding energy E_B of a star of mass M above M₄ is defined as

Fig. 1. Final (a) compactness, (b) specific central entropy, (c) iron core mass, and (d) binding energy above M4 at the onset of core collapse as a function of the initial mass. The top axis shows the CO core mass at core helium exhaustion. Circles and crosses represent BH formation and explosions, respectively, according to the Ertl et al. (2016) criterion, while black and white colors indicate BH formation and explosions based on the Müller et al. (2016) SN model.

$$\begin{array}{c}\hfill E_B=-{\displaystyle {\int }_{M_4}^M}\frac{\mathrm{G}m}{r}\mathrm{d}m.\end{array}$$(2)

We identify specific mass ranges, labeled A, B, C, and D in Fig. 1, during which the compactness parameter, final specific central entropy, iron-core mass, and final binding energy follow specific trends (significant increase or decrease), whose physical origin we investigate in Sects. 5.1, 5.2, 5.3, and 5.4, respectively. These mass ranges are indicated as a function of initial mass (bottom axis of Fig. 1) and of the CO core mass, based on the linear relation we derive between the two for our default assumptions (see also Fig. B.1):

$$\begin{array}{c}\hfill M_{\mathrm{CO}}=0.532M_i-5.31.\end{array}$$(3)

Figure 1 demonstrates that all quantities summarizing the final stellar structure follow a very similar, nonlinear trend as a function of mass, confirming earlier findings (e.g., Timmes et al. 1996; Brown et al. 2001; Sukhbold & Woosley 2014; Chieffi & Limongi 2020; Schneider et al. 2021, 2023, 2024; Takahashi et al. 2023; Temaj et al. 2024). For low masses, the values are approximately constant, with only small variations. A characteristic increase (region A) begins at initial (CO core) masses of about 19 (5.8) M_⊙ and reaches a maximum at CO core (initial) masses of 22 (6.5) M_⊙, commonly referred to as the “compactness peak”. It is followed by a decrease (region B) until a minimum is reached at about 25 (8) M_⊙. After a mass range of about 3 M_⊙ with small variations, all quantities then experience a significant second increase from 30 (11) M_⊙ (region C) until a second peak is reached and the values generally decrease from masses of 40 (16) M_⊙ (region D). As shown in the lowest panel, at this point all models have a high binding energy above M₄. Even if a SN shock were successfully propagating after core collapse in these models, potentially leading to an observable SN, the high binding energy means that the formation of a BH is likely for these mass ranges, independently of the explosion mechanism (Sukhbold & Woosley 2014; Heger et al. 2023).

Our models show a lower intrinsic variability in the final structure than models by Sukhbold & Woosley (2014) and Sukhbold et al. (2016, 2018), just like the studies by Chieffi & Limongi (2020), Chieffi et al. (2021), and Takahashi et al. (2023). This is likely due to our choices of convective boundary mixing that prevent the occurrence of “breathing pulses” during core helium burning (Chieffi & Limongi 2020) that can change the relation between the CO core mass and the core carbon mass fraction X_C at core helium exhaustion. Secondary peaks in the quantities shown in Fig. 1 are not observed for our default assumptions and mass sampling but can appear for different choices of physics. We argue that they are probably real (see Appendices A and B for more details).

All quantities shown in Fig. 1 follow similar patterns because they are intrinsically linked (Fryer 2014; Sukhbold & Woosley 2014; Takahashi et al. 2023; Schneider et al. 2021, 2024; Temaj et al. 2024). For completeness, we repeat these arguments below. The connection between the final central entropy (Fig. 1b) and the final iron-core mass (Fig. 1c) can be understood through the effective Chandrasekhar mass (e.g., Timmes et al. 1996; Woosley et al. 2002; Sukhbold & Woosley 2014; Schneider et al. 2024). In massive stars, core collapse is triggered once the iron-rich core has reached a critical mass. However, in contrast to the cores of low-mass stars, the degenerate iron-rich cores of massive stars are hot, that is, they have a finite temperature (entropy) that needs to be taken into account. This thermal structure leads to additional corrections compared to the classical Chandrasekhar mass $M_{\mathrm{Ch},0}=1.457{\left(\frac{Y_e}{0.5}\right)}^2$ (e.g., Timmes et al. 1996, and references therein), and lead to the approximate effective Chandrasekhar mass (ignoring special and relativistic corrections, and Coulomb corrections)

$$\begin{array}{c}\hfill M_{\mathrm{Ch}}\approx M_{\mathrm{Ch},0}[{\left(\frac{s_e}{\pi Y_e{\mathrm{k}}_{\mathrm{B}}{\mathrm{N}}_{\mathrm{A}}}\right)}^2+1],\end{array}$$(4)

where s_e is the average electronic entropy, which is roughly a third of the average central entropy (Baron & Cooperstein 1990; Timmes et al. 1996). This critical mass is mainly sensitive to changes in entropy because the electron fraction Y_e tends to be very similar between different progenitors² (Sukhbold & Woosley 2014). Thus, the final iron-rich core mass that can form in a massive star before it exceeds M_Ch and its core begins to dynamically collapse is directly linked to the central entropy.

As shown in Fig. 1, the compactness parameter also correlates well with the final entropy and iron-core mass, which may appear surprising at first. This connection can be understood by considering that the compactness is essentially a measure of the mass-radius relation outside the iron core (Chieffi & Limongi 2020). By taking into account that the final iron-rich core is adiabatic, this mass radius relation is directly connected to entropy through polytrope relations (Schneider et al. 2021). Thus the final compactness can be understood as being equivalent to central entropy (see also Fryer 2014), itself linked to the final iron-core mass (and M_Ch) as described above. The mass-radius relation and entropy both enter into the definition of binding energy above M₄ (see Eq. (2)), which explains the observed correlation (Fig. 1d).

3.2. Link to BH formation

To estimate the range of stars for which BH formation is expected, beyond considering the binding energy, we apply the Ertl et al. (2016) explodability criterion to our models. This is shown by the different markers (crosses for explosions and circles for BHs) in Fig. 1. Based on this criterion, stars with initial (core) masses below 30 (12) M_⊙ are expected to preferentially explode. Models in the first compactness peak and beyond the second compactness increase are predicted to form BHs. In addition, we compute the expected explosion outcome of our models using the semi-analytical neutrino-driven SN explosion model of Müller et al. (2016). BH formation is indicated with black and successful explosions with white markers in Fig. 1. We find that the Müller et al. (2016) model also predicts BH formation for stars with a high compactness. The mass range for which successful explosions are expected is more extended compared to the outcome of the Ertl et al. (2020) criterion, reaching up to models with masses of 34 (12.5) M_⊙.

Overall, we find that based on these explosion criteria, BH formation is consistently expected for models with high compactness and high central entropy, including the compactness-peak models at a mass of 22–23 (6.5–7) M_⊙ and models beyond a mass of 34 (14) M_⊙. We confirm that using a compactness threshold for models in which BH formation is expected (e.g., ${\xi }_{2.5}^f>0.4$, see dashed line in Fig. 1a) is too simplified to fully reproduce the predictions by more sophisticated measures of explodability. Considering a larger set of models (Schneider et al. 2021; Temaj et al. 2024; Schneider et al. 2024) we note that there is no clear compactness threshold separating successful and unsuccessful explosions (Maltsev et al. in prep.). For ${\xi }_{2.5}^f\lesssim 0.3$ successful explosions are found, and for ${\xi }_{2.5}^f\gtrsim 0.45$ explosions are unsuccessful, similarly to previous works (Müller et al. 2016; Takahashi et al. 2023; Zha et al. 2023). A compactness threshold can thus only be regarded as a rough first approximation of the core-collapse outcome. Approaches based on the entire final stellar structure, such as the Müller et al. (2016) model, can be applied more generally.

3.3. Emergence of the final structure pattern

To understand at which point the observed pattern in the final stellar structure emerges, we trace the value of the central entropy at key evolutionary stages, shown in Fig. 2. At every stage, there is a general underlying trend of an increasing entropy as function of initial (or CO core) mass, which reflects the increasing core mass of stars as a function of their mass. As a function of time, the central entropy generally decreases. This can be understood as the effect of a changing mass-radius relation of the core, as it becomes increasingly compact due to core contraction (Chieffi & Limongi 2020; Schneider et al. 2021). Central nuclear burning episodes temporarily increase the central entropy, while thermal neutrino losses, which become important after core carbon burning, decrease it. At the end of core C burning (blue diamonds in Fig. 2), a first peak feature in the central entropy appears around 20 (5) M_⊙ initial (CO core) masses. It can be understood as a difference in heat content between models before and after 21.5 (6) M_⊙. Models below this mass experience convective core carbon burning while models above this mass burn carbon radiatively as the burning becomes more neutrino-dominated (see the Kippenhahn diagrams in the online data). From the end of core neon burning on (orange crosses in Fig. 2), the characteristic feature corresponding to the compactness peak emerges around 22 (6.4) M_⊙ and remains until core collapse (black squares in Fig. 2). As we discuss later, this is because central entropy reflects the degree of contraction the core experiences in these models. Between the end of core O burning and core collapse, an additional feature emerges in the central entropy landscape for models above 40 (15) M_⊙. The central entropy drops to lower values, though the general trend of increasing entropy as a function of mass remains. The key processes that determine the final entropy or compactness landscape thus already occur close to the time of core neon burning. For the models above 40 (15) M_⊙ additional processes play a role around the time of core silicon burning.

Fig. 2. Evolution of central specific entropy as a function of the initial mass for different evolutionary stages from core hydrogen depletion until core collapse. Key differences in the final entropy start to appear at the end of core carbon burning. Red markers indicate three models computed with a nuclear network comprising 128 isotopes. Taking into account a greater number of nuclear reactions generally leads to a slightly lower final central entropy, though the general qualitative pattern of the compactness peak is unchanged. The characteristic final landscape is already mostly determined at the end of core Ne burning and remains until core collapse.

The emergence of the final structure landscape at the moment of core neon burning demonstrates that using a small nuclear reaction network is sufficient for characterizing the final structure landscape of stars. However, such models are not well-suited as input for multi-dimensional simulations of core collapse, which require larger reaction networks (Farmer et al. 2016; Renzo et al. 2024). We compare our default models to three additional models with initial masses of 21 M_⊙, 22 M_⊙, and 23 M_⊙, for which we employ a nuclear reaction network of 128 isotopes (Farmer et al. 2016), shown with red markers in Fig. 2. We find the same pattern in the formation of a compactness peak as in models computed with a smaller nuclear reaction network, though the final central entropy is 5% smaller. In addition, the patterns in final central entropy shown here are very similar to the ones found by Takahashi et al. (2023), who used a network comprising 300 isotopes and a finer grid (see their Fig. 11). The exact quantitative details of the final structure landscape, especially at the higher mass end, are thus still affected by changes in Y_e induced by taking into account more nuclear reactions, which ultimately change M_Ch, eff (see Eq. (4)). This can slightly shift the mass range for which a high final central entropy is reached, but does not significantly affect the main qualitative trends described here.

4. A controlled experiment for understanding the origin of the final structure

4.1. Post-core-helium-burning evolution

Beyond core helium burning, the evolution of massive stars proceeds differently compared to earlier stages, as thermal neutrino losses become important. Because of their minuscule interaction cross sections, the vast majority of neutrinos escape the star, causing a tremendous energy loss and greatly accelerating the evolution (e.g., Woosley et al. 2002). At the end of core helium burning, the core contracts. For temperatures around log(T/K)≈8.9 and densities around 10⁵ g cm⁻³, carbon burning ignites in the center (see Fig. 3). The nuclear energy generation rate ϵ_nuc can be expressed as ϵ_nuc ∼ ρT²³X_C during carbon burning, where ρ and T are the density and temperature conditions close to the core, and X_C the initial carbon abundance (Woosley et al. 2002). For high enough temperatures, log(T/K)≳9, pair production, followed by electron-positron pair annihilation, can occur. In rare cases, a neutrino-antineutrino pair is produced, which escapes the star. This process is significant enough to be the dominant energy sink in our models with an energy loss rate ϵ_ν ∼ ρ⁻¹T¹² (Sukhbold & Adams 2020). Photo-neutrinos from electron scattering processes also play a role, but the energy loss they cause is an order of magnitude lower than from pair-annihilation neutrinos. The balance between ϵ_ν and ϵ_nuc determines the final evolution of the core. As shown by the expressions for ϵ_nuc and ϵ_ν, it mainly depends on ρ and T, which are set by the core mass, and the initial central carbon abundance X_C.

Fig. 3. Central temperature as a function of the central density for models with the same core mass and different carbon mass fractions XC after core helium burning (indicated in the color bar). The previous evolution of the core (not shown here for clarity) is indistinguishable for all models, while the post-core-helium-burning evolution is greatly affected by the differences in central carbon abundance. Symbols indicate key evolutionary phases. To the right of the gray dashed line, electron degeneracy dominates the pressure. The inset axis indicates the final compactness of the models as a function of decreasing XC.

To disentangle the effects of the core mass and of X_C on the final structure (Patton & Sukhbold 2020) and understand the formation of the compactness peak identified between regions A and B in Fig. 1, we perform a controlled experiment. For a fixed CO core mass (6.62 M_⊙, corresponding to an initial mass of M_i = 22 M_⊙ for our assumption), we compute models with a modified central carbon abundance. Our fiducial model (X_C = 0.23) corresponds to the compactness peak model, while the models with carbon abundances of 0.28 and 0.17 have final compactness values before and after the peak, respectively. As shown in Fig. 3, small differences in their central abundances cause large deviations in the post core helium evolution. Overall, the X_C = 0.23 model (red line in Fig. 3) reaches higher central temperatures and lower densities after core Ne burning compared to the other models at similar evolutionary stages.

4.2. Carbon burning and beyond

In Fig. 4, we summarize the final evolution of the stellar structure for our three models with the same core mass and a different initial core carbon mass fraction X_C. As indicated in the inset figure in Fig. 4a, these models vary in final compactness and represent models before (X_C = 0.27), in (X_C = 0.23), and after (X_C = 0.18) the compactness peak, where the highest final central entropy and iron core mass are reached. Because they have the same core mass, their central temperature and density conditions at the end of core helium burning are the same (see also Fig. 3) and therefore the neutrino loss rate ϵ_ν is the same. However, at lower X_C, neutrinos increasingly dominate the energetic balance during core carbon burning because the nuclear energy generation rate ϵ_nuc decreases.

Fig. 4. Evolution of key elements of the inner 6 M⊙ core structure as a function of time from carbon ignition to core collapse for models with the same core mass and a different central carbon abundances XC. (a) Kippenhahn diagram with selected elements. Full and dashed lines indicate the carbon and neon-burning fronts, respectively. Shaded (hatched) regions show convective carbon (neon) burning regions in the core. The iron-core masses are indicated with dash-dotted lines. The gray horizontal dashed line indicates the classical Chandrasekhar limit. The inset figure shows the final compactness of these models as a function of decreasing XC. (b) Total mass of 12C in the convective carbon-burning shells for each of the models as a function of time.

In Fig. 4a, we show the progression of the C- and Ne-burning fronts in these models in a Kippenhahn-like diagram as a function of the time from core carbon ignition to core collapse. The location of these nuclear burning fronts are defined as the mass coordinate at which the maximum energy generation rate for a particular burning process (e.g., C or Ne), is reached. In other words, a burning front traces the location at which the maximum burning occurs. These burning fronts are important for the final core structure because their location determines the maximum potential growth of the underlying core mass. A model with carbon burning front that reaches a higher mass coordinate forms a larger carbon-free core, and can eventually form a larger iron-rich core (dashed-dotted lines in Fig. 4).

After core carbon burning ignites in the stellar cores, the C-burning front moves outward in mass for all models as it burns through the available carbon fuel. In the X_C = 0.28 model (dark blue lines in Fig. 4), the high core carbon abundance leads to a nuclear energy generation rate that is large enough to overcome the neutrino energy loss rate, which triggers the development of a relatively small (< 0.2 M_⊙) convective zone (Sukhbold & Adams 2020). It brings in 0.1 M_⊙ of fresh ¹²C (see Fig. 4b) and keeps the carbon-burning front in the center until all the carbon has been burned. Subsequently, the core contracts until the temperature in the carbon-rich layers above the former convective zone reach high enough temperatures for carbon burning to take place and the carbon-burning front moves outward in mass. The location of the burning front thus reflects the amount of core contraction that occurred below it. In this model (X_C = 0.28), this first convective episode is followed by two more successive convective regions that become more extended in mass as the temperature at the burning front increases. At the end of core carbon burning (blue cross in Fig. 4b), the core contracts even further, and core neon burning ignites in the center. The C-burning front moves further out in mass (though its progression is slowed by core neon burning, see also Fig. 5a) until the conditions for convection to occur are reached once again. At this point the carbon-burning front has a high temperature and reaches a region with unburned carbon, which means it reaches a high energy generation rate that highly exceeds the neutrino losses, triggering the formation of a large convective zone with over 0.4 M_⊙ of carbon (see Fig. 4b). Because of the large amount of fuel remaining, the C-burning front keeps producing energy at a high rate and stays at this mass coordinate of 3 M_⊙ until the end of the evolution. As a result, the carbon-free core below can only reach a relatively low mass. Hence, after silicon burning, a relatively low-mass iron core of 1.7 M_⊙ forms.

Fig. 5. Time evolution of the ratio between the specific nuclear energy generation rate ϵnuc and the neutrino loss rate ϵν at the location of the C, Ne, and O-burning fronts for models with the same core mass and different central carbon abundances. The dashed horizontal line indicates where ϵnuc = ϵν. When the energy ratio significantly exceeds this line, convection occurs (Sukhbold & Adams 2020). With a decreasing core carbon abundance, the core becomes more neutrino dominated and neon burning, followed by oxygen burning, occur earlier. As highlighted in the red boxes, this temporarily suppresses carbon burning.

For the other models, which have lower central carbon abundances (red and yellow lines in Fig. 4), central carbon-burning proceeds radiatively. Even though these models have the same mass, and therefore initially the same temperature and density conditions as the X_C = 0.28 model, the lower core carbon abundance means that less energy is generated from carbon burning. Neutrino losses dominate and during the initial central carbon burning, energy is transported solely through radiation (see also Fig. 5). The available carbon in the center is depleted quicker, as no fresh carbon is brought in, and thus the core contracts and the C-burning front moves further out in mass.

The contraction also increases the temperature at the burning front. Once the energy generation rate at the carbon-burning front reaches high enough values to exceed the neutrino losses (see also Fig. 5) two effects occur. First, a convective region forms and the carbon-burning front stays at a constant mass coordinate. Second, this convective carbon-burning shell acts like a mirror between the layers above and the core below. As the core below the shell contracts, the layers above greatly expand (shown in the summary of the gravothermal energy in the stellar structure shown in the online data). The convective region increases in mass, bringing in more fuel for the carbon-burning shell (see Fig. 4b) and prolonging the duration of this burning episode.

For X_C = 0.23 and X_C = 0.17 (red and yellow lines in Fig. 4), the carbon-burning front reaches further out in mass during core carbon burning (time coordinate t ≈ −2) than the X_C = 0.28 model, signifying a larger core contraction. At this point, the mass coordinate of the C-burning front of the X_C = 0.23 model, which corresponds to the compactness peak, is m_C = 1.35 M_⊙, close to the classical Chandrasekhar mass M_Ch, 0 = 1.455 M_⊙, while that of the X_C = 0.17 model exceeds it (m_C = 1.83 M_⊙). As we discuss in Sect. 4.4, this is an important difference because electron degeneracy pressure plays a role in these models, slowing down the core contraction during core carbon burning in the X_C = 0.23 and instead speeding it up in the X_C = 0.17 model.

The end of core carbon burning is the critical moment when the behavior of the compactness peak model (X_C = 0.23) and the model beyond the compactness peak start to differ significantly. In the X_C = 0.17 model, core neon burning already ignites while the C-burning front is still experiencing the first convective episode. In contrast, in the compactness-peak model (X_C = 0.23) by the time neon ignites, convective carbon burning has mostly ended (τ = −2.6). Because the core contraction is slower the model can burn a large fraction of the carbon in the convective region before core Ne burning ignites (see also Fig. 4b). Later on (τ ≈ −3), it quickly burns through the remaining carbon in this region, growing a large C-free core. The C-burning front in the X_C = 0.17 model grows more slowly and less far (m/M_⊙ = 3.18 at τ = −5.15) than in the compactness-peak model (m/M_⊙ = 3.33 at τ = −4). As we show in Sect. 4.3, this is a direct consequence of the earlier core neon and oxygen ignition in the X_C = 0.17 model, which suppresses carbon burning and slows down the progression of the C-burning front.

Moreover, towards the end of the evolution (τ = −5.2) the C-burning front of the X_C = 0.17 model suddenly drops in mass and reaches the same level as the Ne-burning front. As we discuss in more detail in Sect. 4.5, this is caused by a shell merger between the C- and Ne-burning front. Eventually, a lower-mass iron core forms in this X_C = 0.17 model compared to the compactness-peak model.

From this experiment, we have demonstrated that the final core structure of stars in this mass range is already determined at core helium depletion, where the initial conditions for core carbon burning are set. The same effect of a change in the final structure pattern due to a change in X_C can be observed for models with different initial or core masses. This is shown more generally in Appendix A, where we investigate changes induced by the ¹²C(α, γ)¹⁶O rate, which is largely responsible for determining X_C at the end of core helium burning.

We have seen that the transition from a final stellar structure with low compactness and iron core mass to high compactness and iron core mass can be understood as the effect of increasingly neutrino-dominated burning. It results in an overall larger core contraction and the formation of a larger carbon-free core, and eventually, to the formation of a larger iron-rich core. At the compactness peak, the neutrino-dominated carbon-burning conditions are such that most of the available carbon is burned before the end of core neon burning. Afterward, the carbon-burning front quickly moves out in mass and grows a particularly large core C-free core, which eventually forms a particularly large iron-rich core. For the transition between the compactness peak model and the model beyond, we identify three mechanisms that are responsible for the observed drop in compactness: (1) the effect of an earlier central neon and oxygen ignition with decreasing X_C, (2) an increased electron degeneracy, and (3) shell mergers.

4.3. Effect of an early core neon ignition

We have observed that the timing of core neon ignition varies significantly between models in our experiment with the same core mass and a changing initial core carbon abundance X_C (Fig. 4). The consequence of this early neon ignition is shown in Fig. 5, where the evolution of the ratio between the nuclear energy generation ϵ_nuc and the neutrino loss rate ϵ_ν at the location of the C, Ne, and O-burning fronts is plotted as a function of time until core collapse. Phases when ϵ_nuc significantly exceeds ϵ_ν correspond to convective episodes (Barkat & Marom 1990; Barkat 1994; Sukhbold & Adams 2020). In the X_C = 0.28 model (Fig. 5a), core neon burning ignites immediately after core carbon depletion, when a core contraction occurs. Neon burning temporarily slows down the progression of the carbon-burning front and dominates the energy generation rate. However, it does not take long before the slowed down carbon-burning front reaches regions of unburned carbon, quickly becoming the main energy source again and forming a large convective region that remains at the same mass coordinate until the end of the evolution (see also Fig. 4b).

In the compactness-peak model (X_C = 0.23, Fig. 5b), core neon and oxygen burning ignite in rapid succession after the convective carbon-burning region has exhausted its fuel. They dominate the energy generation, slowing down the progression of the carbon burning front until the end of core oxygen burning is reached. The ensuing core contraction helps the carbon-burning front reach layers where the nuclear energy generation rate dominates compared to the neutrino losses, triggering the formation of a large convective zone and marking the final location of this burning front.

In the X_C = 0.17, model (Fig. 5c), while the convective carbon-burning episode is still ongoing, the conditions for neon burning are already reached in the center due to the stronger preceding contraction of the core aided by exceeding M_Ch, 0 (see Sect. 4.4). While core neon burning is more neutrino-dominated in this model and produces less energy, this simultaneous burning reduces the luminosity of the carbon-burning front, impacting the extend of the convective region above. Once oxygen burning ignites shortly after, it counters the contraction and helps drastically slow down the progression of the now very neutrino-dominated carbon-burning front. In the subsequent phases, neon and oxygen burning dominate the energy generation, followed by core silicon burning. After silicon depletion, the neon-burning front, which moves out rapidly in mass, merges with the carbon-burning front (see Sect. 4.5). The slowed C-burning implies a smaller carbon-free core and ultimately, a smaller iron core mass and compactness.

4.4. Origin of the changing timing of core neon and oxygen ignition

Here, we investigate the origin of the systematically earlier ignition of neon and oxygen for lower X_C. We find that it is mainly due to a systematic shortening of the duration of core carbon burning for a decreasing amount of fuel (X_C). This is shown for a large set of models with the same core mass and a varying core carbon abundance in Fig. 6a. For reference, we show the final compactness of these models in Fig. 6b. The contraction phase that follows core carbon depletion (when X_C < 10⁻⁴) before core neon burning ignites is significantly shorter than the duration of core carbon burning, on a neutrino-losses accelerated thermal timescale on the order of years (see Fig. 6c). It therefore plays a smaller role than the duration of core carbon burning in determining the timing of core neon ignition.

Fig. 6. Properties of a series of models with the same core mass and a varying initial core carbon abundance XC, indicated with colors. (a) Duration of core carbon burning. (b) Final compactness as a function of XC. Models that undergo convective core carbon burning are highlighted with black outlines. (c) Duration of the phase between core carbon depletion and core neon ignition as a function of mC, the maximum mass coordinate reached by the carbon-burning front before core neon burning. The classical (MCh, 0) and effective (MCh) Chandrasekhar mass at core C depletion are shown as vertical dashed and full lines, respectively.

Besides the duration of core carbon burning, a secondary effect related to electron degeneracy also appears to play a role in accelerating the evolution after core carbon burning, as first pointed out by Sukhbold & Woosley (2014). This is shown in Fig. 7, where the evolution of the electron degeneracy parameter η = μ/k_BT is displayed as colored contours in a Kippenhahn-like plot for each of our example models. η ≈ 0 indicates partial degeneracy, while η ≫ 0 signifies that electrons are strongly degenerate. We also show the evolution of the classical (M_Ch, 0) and effective Chandrasekhar masses (M_Ch) (see Eq. (4)). For a decreasing X_C, the central region where electrons are degenerate during core carbon burning is more extended in mass. This can be understood as a symptom of the overall stronger core contraction experienced by models that are more neutrino-dominated. For the model with X_C = 0.28 (Fig. 7a), electron degeneracy plays a small role during central carbon burning. The C-burning front (at the base of the convective carbon burning shell) remains significantly below the values of the classical and effective Chandrasekhar mass. For the compactness peak model (X_C = 0.23, Fig. 7b), degeneracy pressure helps support the core at the end of core carbon burning, reaching 15% of the total pressure in the center (η_c = 0.93) before it is lifted during core oxygen burning. The carbon burning front remains close to, but below M_Ch, 0 and M_Ch during core carbon burning. For the X_C = 0.17 model (Fig. 7c) the region where electrons are at least partially degenerate extends up M_Ch, 0. The carbon-burning front exceeds M_Ch, 0 and even reaches the value of M_Ch. This appears to help accelerate the core contraction, which is somewhat surprising since the Chandrasekhar mass is known to mainly play a role at full degeneracy, and deserves further study.

Fig. 7. Kippenhahn diagram of the inner 6 M⊙ core structure of models with the same core mass and different core carbon abundance at core C ignition. Colors indicate the dimensionless electron degeneracy parameter η = μ/kBT (electrons are partially degenerate for η ≈ 0 and very degenerate for η ≫ 0). Convective zones are highlighted by the hatched regions and the dotted vertical lines indicate, from left to right, the moments when 12C, 20Ne, 16O, and 28Si are depleted at the center (central abundance lower than 10−4). The gray full and dashed lines indicate the classical (MCh, 0) and effective (MCh) Chandrasekhar mass, respectively. Electron degeneracy during core carbon burning increases for models with lower central carbon abundances.

To understand if this effect applies more generally, we investigate the behavior of a larger set of models with varying X_C. In Fig. 6c, we show the duration of the phase between core carbon depletion and neon ignition as a function of m_C, which is defined as the maximum mass coordinate of the carbon-burning front m_C during core carbon burning (i.e., the base of the furthest convective carbon-burning shell)³. The changes in this timescale are very close to the compactness pattern shown in Fig. 6b as a function of decreasing core carbon abundance. A similar trend in this timescale was also noted by Chieffi et al. (2021). At the end of core carbon burning, all models reach the same value of the classical Chandrasekhar M_Ch, 0 = 1.455 M_⊙, while the effective Chandrasekhar mass varies slightly between models due to their different central entropy, with values of M_Ch ∼ 1.82 − 1.84 M_⊙. Only in the model with the lowest core carbon abundance (X_C = 0.17, shown in more detail in Fig. 7c), m_C exceeds M_Ch.

From the variation of the timescale between core carbon depletion and core neon ignition, we derive the following interpretation. Models with X_C < 0.25 (m_C > 1.2 M_⊙) are strongly neutrino-dominated and burn carbon radiatively in the core, experiencing a stronger contraction than those with larger X_C that still burn carbon convectively in their cores. However, degeneracy pressure helps to counter the core contraction towards the end of core carbon burning and to delay the ignition of core neon burning. For models where m_C > 1.5 M_⊙, the entire region below M_Ch, 0 is at least partially degenerate and degeneracy pressure no longer helps to counter the core contraction, leading to an earlier core neon ignition. Sukhbold & Woosley (2014), Sukhbold et al. (2018), and Sukhbold & Adams (2020) noted a smaller size of the oxygen-burning shell once the base of the carbon-burning shell m_C exceeds M_Ch and linked it to the role of degeneracy. It is noteworthy that we observe the same phenomenology here (see Fig. 7) despite the different assumptions and methods we use. As we discuss in Sect. 5.4, the same mechanism helps explain the second drop in compactness (region D in Fig.1) though in this case it involves the oxygen-burning front and an earlier core Si ignition.

4.5. Shell mergers

In our experiment, the model with the lowest initial carbon abundance X_C = 0.17 experiences a merger of the C- and Ne-burning layers shortly after core Si burning. This is highlighted within the red rectangle in the Kippenhahn diagram (Fig. 8a). The location of the C-burning and Ne-burning shells are indicated with markers. The origin of this shell merger can be understood by inspecting Fig. 8b, which shows the entropy profile in the highlighted region. The entropy in the layers corresponding to the location of the Ne-burning shell (2.5 to 2.9 M_⊙) rapidly increases. This is because after Si-burning ignites, a rapid contraction of the core below the Ne-burning front leads to enhanced Ne burning, which generates entropy. Eventually, the entropy of these layers exceeds that of the carbon-burning layers above (see the entropy profile at a time of log₁₀[(t_cc − t)/yr]= − 3.33). This entropy contrast means that these layers are unstable against convection (Schwarzschild criterion). No convective boundary mixing is assumed in these models after core helium burning, which means that the mixing is solely caused by this entropy change. In Appendix C, we perform a resolution test and find that this shell merger is barely affected by numerical uncertainties. The mixing leads to the merging of these shells and to a new configuration of the stellar structure.

Fig. 8. Final evolution of the inner stellar structure of the model with XC = 0.17 which experiences a shell merger between the C and Ne-burning shells (highlighted in the red box in panel a). (a) Kippenhahn diagram in which colors indicate the regions dominated by nuclear burning or neutrino losses and gray hatched regions indicate convective areas. Brown lines show contours of constant (logarithmic) density. Colored markers trace the burning front for C, Ne, and O burning. Small colored vertical lines indicate the times at which the entropy profiles are shown in the bottom panel. (b) Entropy profiles in the mass and time range highlighted in the red box in the top panel. The entropy in the neon-burning region (2.5 − 2.9 M⊙) gradually increases until it exceeds the entropy of the C-burning layers above, triggering a merger between these regions.

After the merger, the Ne/C burning shell suddenly drops to a lower mass coordinate and a large convective zone develops above. These layers suddenly expand, leading to a dramatic density drop in these regions (see the density contours in Fig. 8a). This particular model additionally experiences a merger of the Ne-burning and O-burning shells (visible at a mass coordinate of 2 M_⊙ at log₁₀[(t_cc − t)/yr]= − 5.2). These mergers lead to a low final mass coordinate of the C/Ne/and O-burning fronts, eventually limiting the development of the Si-burning front, and thus setting a limit to the maximum mass of the iron-rich core.

5. Final core structures of actual stellar models

In the previous section, we have identified key mechanisms that lead to a change in the final core structure through an experiment in which we varied the core carbon abundance for a fixed core mass. We now apply the insights gained to more realistic simulations of massive stars.

The initial conditions for central carbon burning are set by the core mass, which determines the central density and temperature, and by the initial carbon abundance. Hydrostatic equilibrium implies a characteristic trend in the central density and temperature of stars as a function of mass (e.g., Kippenhahn et al. 2013). While their central temperatures systematically increase with mass, their central densities decrease. The central carbon mass fraction at core helium depletion decreases systematically as a function of mass. This is a signature of the increased importance of the ¹²C(α, γ)¹⁶O reaction rate with respect to the triple-alpha process at the end of core helium burning for cores of higher mass and lower density, which effectively leads to an increased destruction of carbon (e.g., Arnett 1972). For increasing masses, the higher central temperatures, lower central densities, and lower central carbon abundance imply stronger neutrino losses (see Sect. 4.2), with important consequences for the final structure.

5.1. Origin of the increase in final compactness (region A)

We first focus on stars at the lower-mass end before the compactness peak (region A in Fig. 1). The change in interior structure as a function of mass is summarized in Fig. 9, in which we show the progression of the C-burning front, the C-free core mass (which follows the progression of neon burning), and the iron core mass as a function of time.

Fig. 9. Evolution of the C-burning front, the C-free core (which follows the progression of neon burning), and the iron-rich core as a function of time from carbon burning to core collapse for models before and up to the compactness peak (region A of Fig. 1, see inset). As carbon burning becomes increasingly neutrino-dominated, the C-burning front moves further out in mass, ultimately leading to the formation of a larger C-free and iron-rich core.

In the lowest-mass model (20 M_⊙) core carbon burning is convective initially. The C-burning front remains at the center until the carbon fuel is exhausted (at τ ≈ −0.75) and the C-burning front moves further out in mass until the conditions for convective core carbon burning are reached once again. Because of the additional fuel from convective core carbon burning, the progression of the C-burning front proceeds somewhat differently in this 20 M_⊙ model compared to the higher-mass models. The carbon-burning phase lasts longer and as a result, neon burning ignites later while the C-burning front is still moving out in mass (as shown by the growth of the C-free core at τ ≈ −0.75, dotted lines in Fig. 9). Eventually, the C-burning front reaches a higher value (2.3 M_⊙) than in the 21 M_⊙ model, which proceeds radiatively. At the end of the evolution, the model forms a relatively low-mass iron core of 1.59 M_⊙.

The transition from convective to radiative core carbon burning has previously been invoked as the origin of patterns in compactness or iron-core mass (Brown et al. 1999; Sukhbold & Woosley 2014; Sukhbold & Adams 2020; Chieffi & Limongi 2020; Takahashi et al. 2023). However, we can clearly see in Fig. 9 that the 21 M_⊙ model, which burns carbon radiatively, does not immediately develop a larger iron-core mass (1.56 M_⊙) than the 20 M_⊙ model (1.59 M_⊙). This transition in the energy transport mechanism from convection to radiation can be understood as a symptom for neutrino losses becoming increasingly dominant, but it is not the cause of the increase in compactness and iron core mass.

For increasing masses, the models experience larger neutrino losses during core carbon burning because they have higher central temperatures, lower densities, and a lower central carbon abundance (see Sect. 4.2). The neutrino losses imply that the layers above the burning front can contract more significantly than when neutrino losses are not dominant. This causes an overall more significant and earlier core contraction in these stars, which is reflected in the location of the burning fronts. The smaller amount of carbon fuel, together with the increased nuclear energy generation rate in higher-mass stars accelerate the outward progression of the C-burning front. For example, as shown in Fig. 9, at the end of the evolution, the C-burning front reaches 2.3 M_⊙ in the 21 M_⊙ model, while it reaches 3.35 M_⊙ in the 22 M_⊙ model (which corresponds to the compactness peak). This progression of the C-burning front halts when the energy generated by carbon burning greatly exceeds the neutrino losses and convection sets in. For the compactness-peak model (22 M_⊙) the C-burning front reaches a particularly high final mass coordinate. This is because most of the carbon in the large convective shell that forms during core carbon burning is already exhausted by the time of core neon ignition (τ ≈ −2.1), and the carbon-burning front quickly burns through the remaining fuel. A further progression of the C-burning front is a symptom of a larger core contraction which results in the growth of a systematically more massive C-free core as a function of initial mass (dotted lines in Fig. 9). In all models, the C-free core mass grows almost up to the final value of the C-burning front. Consequently, the O- and Si-burning fronts also move further out in mass and eventually, a more massive iron-core grows in the higher-mass models (dash-dotted line in Fig. 9).

5.2. Origin of the decrease in final compactness (region B)

For models in region B, the early behavior of the carbon-burning front is similar to the neutrino-dominated models in region A (see Fig. 10). For higher-mass models, the C-burning front moves systematically further out in mass (e.g., reaching m ≈ 1.35 M_⊙ for the 22 M_⊙ model, compared to m ≈ 1.6 M_⊙ for the 24 M_⊙ model). However, around the time of core carbon depletion (τ ≈ −2 Fig. 10), these models experience a clear reversal of the trends observed in region A. The convective carbon-burning episode at τ ≈ −2, where the C-burning front remains at a fixed mass coordinate, is systematically shorter. Afterward, the C-burning front moves out systematically slower in mass for increasing masses and reaches lower final values (e.g., m ≈ 3.35 M_⊙ for the compactness-peak model with an initial mass of 22 M_⊙ compared to m ≈ 2.28 M_⊙ for the 25 M_⊙ model). Consequently, the final iron-core mass is smaller in the higher-mass models compared to the lower-mass ones.

Fig. 10. Same as Fig. 9 but for models beyond the compactness peak (region B in Fig. 1, see inset). As a function of initial mass, neon burning ignites systematically earlier (as traced by the C-free core mass) and slows down the progression of the C-burning front. Shell mergers between the C/Ne/O-burning fronts are responsible for the sudden drops in the C-burning front at τ ≈ −5. Ultimately, for increasing masses, these models form systematically smaller iron-rich cores.

The observed change in final compactness and iron core mass can be traced back to the same mechanisms as outlined in Sect. 4. Higher neutrino losses in higher-mass models and the ensuing contraction of the core, together with the low initial core carbon abundance and shorter duration of core carbon burning, mean that the conditions for core neon burning to occur are met systematically earlier than in lower-mass models. This can be observed by the timing of the growth of the C-free core (dotted lines in Fig. 10), which occurs systematically earlier in higher-mass models. Additionally, as shown in Fig. 11, the extent of the degenerate region within the core increases and helps accelerate the core contraction in models where the C-burning front (at the base of the C-burning shell) exceeds the classical Chandrasekhar mass. This earlier central neon ignition increases the core luminosity which slows down the core contraction above the neon-burning front. This causes the observed slower increase in the C-burning front which implies a slower growth of the C-free core.

Fig. 11. Same as Fig. 7 but for three stellar models of our grid after the first compactness peak.

In addition, central neon burning is more neutrino-dominated for higher mass-models, and leads to a smaller increase in entropy during these burning phases. After core silicon ignition (as traced by the growth of the iron-rich core, dash-dotted lines in Fig. 10), we observe a phase during which neon and oxygen shell burning become more energetic for higher masses, leading to a higher entropy, which ultimately causes shell mergers just like in our previous experiment (see Sect. 4.5). These shell mergers are found for all models in region B beyond the compact peak and can be identified by a sudden drop in the C-burning front shortly after core silicon burning sets in and the iron-core mass starts to increase (e.g., at τ ≈ 5 in the 24 M_⊙ model). In most cases, this involves a merger of the neon-burning and carbon-burning shells. In the highest mass models (M_i > 23 M_⊙), the oxygen and Ne-burning shell merge additionally. The shell merger implies a reconfiguration of the stellar structure and a large expansion of the layers above the C/O burning shells. The lower final location of the C and O-burning fronts after the shell mergers (e.g., m ≈ 2.28 M_⊙ for the 25 M_⊙ model) limit the progression of the Si-burning front. Consequently, the iron-core mass of these higher-mass models remains small.

5.3. Origin of the second compactness increase (region C)

For models between 25 and 28 M_⊙, neutrino losses become increasingly important as the central temperature increases and the initial central carbon abundance decreases further. The core contracts even further during core carbon burning and the carbon-burning front moves further out in mass. Core neon burning ignites even earlier and becomes more neutrino-dominated, which means neon burning proceeds in smaller convective zones. These models are characterized by a large convective oxygen-burning shell that temporarily stops the progression of the carbon-burning front after it reaches its initial location. The final compactness and iron-core mass of these models remains low (see Fig. 1).

After an initial mass of about 27 M_⊙, the final compactness increases again (see Fig. 1). This is due to core neon burning becoming fully neutrino dominated. This is best seen in the interior structure evolution of the models with initial masses from 27 M_⊙ to 32 M_⊙ (available online). For higher initial masses, neon burning changes from a succession of multiple small convection zones (four in the 27 M_⊙ model), to one or two tiny (< 0.2 M_⊙) convective zones (27.5–30 M_⊙) followed by radiative burning, before eventually transitioning to fully radiative, neutrino-dominated neon burning (from models of 31 M_⊙, where convective core neon burning regions completely disappear). Once again, this transition from convective to radiative burning does not immediately lead to a change in the final compactness. Rather, the compactness gradually increases for fully neutrino-dominated models as the core contracts further to compensate the energy loss due to neutrinos.

In Fig. 12 we summarize the changes in the stellar structure for models in region C. All models show a very similar trend as in the neutrino-dominated models in region A. As the neutrino losses increase and the core contracts further, the C-burning front moves further out in mass until it ignites a convective-burning shell. In these models, neon burning ignites even earlier than in region B (as shown by the progression of the C-free core at τ ≈ −2, dotted line in Fig. 12). However, for these masses, it becomes neutrino-dominated, which means that the contraction of the layers above the Ne-burning front cannot be slowed significantly by this early ignition of neon. Furthermore, as the initial carbon abundance decreases, so does the neon abundance after core carbon burning, which shortens the duration and impact of core neon burning. Thus core neon ignition barely affects the progression of the C-burning front and neon shell burning no longer leads to a high enough increase in entropy to cause shell mergers with the C-burning layers. The Ne-burning front moves systematically further out in mass for higher-mass models as the core contracts and the C-free core almost reaches the C-burning front (e.g., 2.09 M_⊙ in the 28.0 M_⊙ model, and 3.9 M_⊙ in the 40.0 M_⊙ model, see dotted lines in Fig. 12). Eventually, for increasing masses, the Si-burning front can also move further out in mass and a systematically larger iron-core mass forms (e.g., 1.7 M_⊙ in the 28.0 M_⊙ model compared to 2.35 M_⊙ in the 40 M_⊙ model, see dash-dotted lines in Fig. 12).

Fig. 12. Same as Fig. 9 but for models in the mass range of the second compactness increase (region C in Fig. 1, see inset). In these models, both carbon and neon burning are fully neutrino-dominated.

5.4. Origin of the second compactness decrease (region D)

For models beyond 40 M_⊙, we once again observe a drop in the final compactness (see Fig. 1). In these models, both carbon and neon are neutrino-dominated and the respective nuclear burning fronts move out so far in mass from the central regions that they no longer play a large role in setting the final iron core mass (see Fig. 13). Instead, the core oxygen and silicon-burning phase increase in importance, as shown in Fig. 14. Oxygen burning produces a large amount of energy per gram and occurs convectively (Woosley et al. 2002). In analogy to the effects observed for the C- and Ne-burning front in lower-mass models, after core oxygen burning the oxygen-burning front initially moves further out in mass the more neutrinos dominate in the center (as indicated in Fig. 13). The timescale between core oxygen and core silicon ignition shortens for increasing masses (as shown by the timing of the growth of the iron-rich core, dash-dotted lines in Fig. 13). The (at least partially) degenerate core region becomes more extended, as it generally reaches up to the base of the oxygen-burning front (see Fig. 14). The core contraction accelerates for models in which the oxygen-burning front reaches a mass coordinate higher than the (classical) Chandrasekhar mass. For example, in the 42 M_⊙ model, the oxygen-burning front (at the base of the oxygen-burning shell) even reaches the effective Chandrasekhar mass (indicated by a dotted line in Fig. 14c) with m ∼ 1.85 M_⊙ at log(t/yr) = − 2.5. When silicon burning occurs during the oxygen shell-burning phase, it slows down the contraction of the layers above, and with it the progression of the oxygen-burning front. Eventually, this limits the maximum progression of the Si-burning front, which forms a lower iron-core mass than in models in which core silicon burning occurs later. For example, as shown in Fig. 14, in the 43.0 M_⊙ model the oxygen-burning front (dashed lines) remains at a mass coordinate of 2.1 M_⊙ after core silicon burning ignites (as traced by the growth of the iron-rich core, dash-dotted lines). It only moves further out in mass shortly before core collapse (τ = −7), which does not leave enough time for the iron-core mass to grow further in mass than 1.85 M_⊙ in this model.

Fig. 13. Same as Fig. 9 but for models in the mass range of the second decrease in compactness (region D in Fig. 1, see inset). The dashed line traces the O-burning fronts. For increasing masses, the progression of the O-burning front is slowed by a systematically earlier ignition of core silicon burning (as traced by the iron core mass).

Fig. 14. Same as Fig. 7 but for three stellar models of our grid after the second compactness increase. Here, neon burning is fully neutrino-dominated and radiative, and not visible in the figure. After convective core oxygen burning (hatched region in the center) the degenerate core region reaches up to the oxygen burning front (at the base of the convective oxygen-burning shell) for all models. Panel c: When this front exceeds the (effective) Chandrasekhar mass, the degenerate core contracts and core silicon burning ignites earlier, slowing the contraction of the layers above and affecting the extent of the convective region.

Beyond region D, for models with higher masses the final compactness remains lower than during the second peak (see Fig. 1a) but above a value of 0.4. The 44 M_⊙ model has an unusually high compactness (0.7) compared to the other models with initial masses larger than 40 M_⊙. We link this back to the oxygen-burning front remaining at a rather low mass coordinate (close to 1.45 M_⊙, not exceeding the effective Chandrasekhar mass) after core oxygen burning, triggering a small (0.1 M_⊙) convective episode. The core contracts less during core oxygen shell burning than for the adjacent models, which means silicon burning ignites later and does not slow down the progression of the oxygen-burning front as effectively, leading to a larger final compactness and iron core mass. Small changes in the energy generation rate or mixing could influence the existence of this small oxygen-burning episode and we therefore consider this model not to be representative of the more general patterns described here. We note that the emergence and disappearance of such small secondary peaks in the final structure landscape can be observed when physical assumptions are varied, such as the efficiency of semiconvection (see Appendix B) and deserve further study. Despite the observed variations, it should be noted that for all these models, the binding energy above the inner core is so high (see Fig. 1 d) that the formation of a BH is expected even if a successful explosion were triggered in the core.

6. Global picture

Based on our findings, we can derive a global picture of the physical processes that determine the final core structures of massive stars. After core helium burning, thermal neutrino emission becomes important and complicates the final evolution by taking away energy from the interior, leading to an acceleration of the evolution. The transition from convective to of radiative carbon and neon burning is a symptom of neutrino losses becoming dominant, but not the cause of the observed changes in the final structure patterns. For higher masses (lower densities, higher temperatures, and lower central carbon abundance), neutrino losses become increasingly important and affect the efficiency and timing of nuclear burning episodes and with it the associated core growth.

We identified physical mechanisms explaining the emergence and decline of the two prominent “compactness peaks” (see Fig. 1) found in several independent studies, that is, mass ranges for which stars develop a high final compactness, iron core mass, central entropy, and binding energy. These are summarized schematically in Fig. 15, where we show how neutrino-dominated nuclear burning in the cores of massive stars post core helium depletion leads to different stellar structures. These outcomes are “written” in the stellar cores at core helium depletion, when the initial burning conditions of the subsequent burning phases are set.

Fig. 15. Schematic representation of the mechanisms leading to the emergence and decline of peaks in final iron core mass and compactness. Horizontal gray dotted guidelines are included to aid the comparison. Star I – Before peak: Schematic stellar structures of a star with a low final iron core mass at two points in the evolution represented by two adjacent wedges. The initial mass fraction of the fuel is still high (see red area in the rectangle) but the main burning stage is neutrino-dominated (arrows). Because neutrinos take away energy, the core contracts and the burning front (red line) moves out in mass as the burning progresses. The large amount of fuel prevents it from moving far out before a convective zone forms above the burning front once the energy generated is high (shaded red area). Ultimately (right wedge), it grows a small fuel-free core and the silicon-burning front (yellow line) cannot move far out in mass, forming a small iron core (black area). Star II – Peak: With a higher mass and lower initial fuel mass fraction, the burning is even more neutrino dominated and a strong contraction occurs. Hence the burning front moves further out in mass than in I, but stays just below the effective Chandrasekhar mass. Aided by degeneracy pressure support, the star burns through almost all available fuel before the next burning stage ignites, after which it quickly moves out in mass, growing a large fuel-free core and eventually a large iron core. Star III – Beyond peak: For an even lower initial fuel abundance and higher mass, the burning is even more neutrino dominated and the core contraction accelerates. The burning front moves further out and exceeds the effective Chandrasekhar mass. This triggers a fast contraction of the partially degenerate core and the next burning stage (blue area) ignites early. This next stage suppresses nuclear burning at the front above. As a result, the burning front moves out in mass slowly and eventually (right wedge) the star grows a low-mass iron core.

6.1. Increase in final iron core mass (Region A and C):

For increasing stellar masses, the central nuclear burning source (C or Ne) becomes increasingly neutrino-dominated. As observed in multiple studies (e.g. Brown et al. 2001; Sukhbold & Woosley 2014; Sukhbold et al. 2018; Sukhbold & Adams 2020; Chieffi & Limongi 2020; Takahashi et al. 2023), energy transport at the center transitions from convection in a large central region followed by shell burning, to ever smaller and more numerous successive convective burning episodes, until it becomes fully radiative for stars in which nuclear burning is fully neutrino-dominated. In Fig. 15I, the structure of a star with a neutrino-dominated, radiative core is shown schematically mid burning (left wedge) and at the onset of core collapse (right wedge). The decreasing amount of fuel together with the neutrino dominance lead to an acceleration of the core contraction and to a faster and further outward progression of the burning front (the point where the maximum nuclear energy generation rate is reached, red lines in Fig. 15I) and associated growth of the C/Ne-free core below. When the burning front reaches layers with a high C/Ne abundance, the generated nuclear energy becomes large enough to significantly exceed the neutrino losses, and convection sets in above the front (shaded red area in Fig. 15I). The star with the structure I still contains a large amount of fuel, which means (a) that a small core contraction is sufficient for convection to set in and (b) that the newly formed convective region contains enough fuel to remain until core collapse without the need of a large core contraction. This limits the growth of the fuel-free core below as the burning front moves only slowly out in mass. Ultimately (right wedge in Fig. 15I), this structure grows a relatively small fuel-free core, which in turn limits the growth of the Si-rich core, which then forms a relatively low-mass iron-rich core.

The final location of the burning front reflects the contraction of the inner core and sets the maximum core growth for the subsequent nuclear-burning episodes. Through this process, for increasing masses, stars in which the central burning source becomes neutrino-dominated grow increasingly massive C/Ne-free cores, which in turn form massive oxygen-free cores and eventually, more massive iron-rich cores. For higher masses, as the cores contract further and their density increases, electron degeneracy also starts to play a role. Degeneracy pressure contributes to supporting the core during the main burning phase (carbon or oxygen, depending on the mass range), delaying the ignition of the next burning phase as long as the burning front does not exceed the effective Chandrasekhar mass.

At the compactness peak (structure II in Fig. 15), nuclear burning at the center is strongly neutrino-dominated due to a higher mass (lower initial density) and low initial fuel mass fraction. The strong core contraction leads to an outward progression of the burning front until it nearly reaches the effective Chandrasekhar mass. Support from degeneracy pressure helps slow the contraction during the main burning phase such that nearly all the fuel in the convective region is burned before the next burning stage ignites in the center. Afterward, the burning front quickly burns through the former convective region as the core contracts further. It grows a particularly massive fuel-free core and eventually, a massive iron core.

6.2. Decrease in final iron core mass (Region B and D)

For even higher masses, stars develop structures that are even more neutrino-dominated (star III in Fig. 15). Two mechanisms, (a) a lower initial fuel abundance decreasing the duration of the main fuel burning phase (core C-burning or Ne-burning), and (b) electron degeneracy no longer helping support the stellar structure once the main burning front and the fuel-free core it produces exceed the (effective) Chandrasekhar mass, both lead to an acceleration of the core contraction. What follows is an early ignition of the next burning fuel while the main nuclear burning phase is still in progress (e.g., Ne ignition during C burning or Si ignition during Ne/O burning, indicated by the blue central region in Fig. 15III). The luminosity of this new nuclear energy source contributes to slowing down the contraction of the layers above and prevents the further growth of the C/O-free core. Consequently, the Si-burning front cannot move far out in mass and grows a lower-mass iron core than in star II.

Additionally, in particular in region B, mergers of the Ne- and C-burning shells, and, in some cases, of the Ne and O-burning shells, can occur due to the high energy generation rate and entropy at these fronts. These lead to a reconfiguration of the stellar structure in which the burning fronts all merge together at the lowest mass coordinate. If these mergers occur before the end of core silicon burning, they limit the progression of the silicon-burning front, and, with it, the growth of the iron-rich core. Eventually, stars that experience such shell mergers form iron cores of even lower mass. Shell mergers before the end of core Si burning are thus an additional mechanism which produces final stellar structures with low iron core masses.

Because of these physical mechanisms, we expect a nonmonotonic landscape of final compactness, iron core mass, central entropy, and binding energy as a function of the initial and core masses of massive stars. Based on a neutrino-driven SN model we have shown that stars with high compactness, iron core masses and binding energy are predicted to form BHs. Hence, the existence of these robust final structure patterns of massive stars observed by multiple independent studies imply the formation of BHs at characteristic mass ranges.

7. Discussion

7.1. Uncertainties in stellar physics

Our stellar evolution models are subject to uncertainties concerning the physics of massive stars that can affect the lanscape of the final structure of stars, and with it, their fate.

7.1.1. Nuclear reaction rates:

As discussed previously, the final core structures of stars are sensitive to the conditions under which core helium and carbon burning take place. In particular, the ¹²C(α, γ)¹⁶O is a key helium-burning process which is still very uncertain and has a tremendous impact on the final fate of stars (e.g. Austin et al. 2014; Farmer et al. 2020). This nuclear reaction rate sets the amount of carbon left after core helium burning (X_C) by fusing carbon into oxygen. As shown in Appendix A, varying this reaction rate changes the core mass range at which stellar structures reach a high final compactness and iron-core mass and are expected to form BHs. A higher nuclear reaction rate resulting in a lower core carbon abundance means that neutrino-dominated burning occurs at lower CO core masses, systematically shifting the compactness landscape (Sukhbold & Adams 2020; Schneider et al. 2023). For example, a 10% higher carbon alpha-capture rate leads to a compactness peak occurring at 0.4 M_⊙ lower CO core masses (see Appendix A). In turn, this increases the number of stars expected to form BHs, as stars with lower (CO core) masses are more common. It will also affect the final masses of BHs, which is expected to be proportional to the final mass of a star and is generally lower for progenitors with lower masses (see also Schneider et al. 2023).

Similarly, the value of the ¹²C+¹²C cross section also systematically shifts the final core structures of stars (Chieffi et al. 2021). The newly determined, higher rate found by Tumino et al. (2018) due to low-lying resonances shifts the compactness landscape to slightly lower core and initial masses compared to the classical rate we use in the present work (Caughlan & Fowler 1988) by allowing core carbon burning to take place at lower temperatures and densities. Carbon thus ignites under conditions where neutrino losses are less severe, which means that the overall core contraction is slightly less significant. This results in a systematic shift in the final compactness to slightly lower values (0.05 in ξ_2.5, see Fig. 9 in Chieffi et al. 2021). Finally, the triple-alpha reaction rate also plays an important role, though it is less uncertain (Austin et al. 2014).

7.1.2. Mass loss:

Wind mass loss (and equivalently, metallicity) has a large impact on the properties of stars and is still very uncertain. Latest observational constraints suggest that wind mass loss at solar metallicity may be less strong than we assume here. When mass loss is strong enough to remove a large fraction of the hydrogen-rich envelope, it can have a large impact on the final structure of stars as it changes the evolution of the helium core, and with it the initial conditions for core carbon burning. In particular, strong wind mass loss can lead to a less efficient hydrogen-burning shell and to a recession of the convective helium-burning core, which results in a large core carbon abundance X_C at core helium depletion. This in turn shifts the parameter space for neutrino-dominated carbon burning to higher CO core masses. Stars that experience strong wind mass loss are thus generally more explodable (just like binary-stripped stars, see Sect. 7.3). In our models with high enough initial masses (≳35 M_⊙), stars may also enter a Luminous Blue Variable phase in which eruptive mass loss has been observed to occur (Smith et al. 2004; Davies et al. 2018). The mechanism behind these eruptions has been suggested to be related to helium opacity in the outer layers (Jiang et al. 2018; Grassitelli et al. 2021). Unlike line-driven winds, this mechanism would also be important at low metallicity. This highly uncertain mass loss is not taken into account in our models but can lead to the removal of the entire hydrogen-rich envelope, which would significantly affect the final core properties, and lead to a final core structure that explodes more easily.

7.1.3. Mixing:

How mixing proceeds in stellar interiors remains an importance uncertainty. Convective core boundary mixing, computed in our models through step overshooting, influences the size of convective cores. During core hydrogen burning, it helps determine the initial mass to helium core mass relation. Similarly, core overshooting during core helium burning influences the resulting CO core mass and thus has a particularly large influence of their final fate. Temaj et al. (2024) showed that core overshooting systematically shifts the compactness landscape of stars as a function of their initial masses, as a larger overshooting allows a star to behave as if it had the core of an initially higher-mass star. This also implies a larger number of BHs formed, as more stars can reach the mass range for BH formation (see also Fig 6. of Schneider et al. 2023). As a function of the CO core mass, the compactness landscape is similar to that of single stars, though the compactness peak shifts towards higher CO core masses (1.3 M_⊙ by varying the step overshooting scaling parameter α_OV from 0.05 to 0.5, Temaj et al. 2024). This shift is because, as convective overshooting affects the core mass of stars, it also changes their luminosity and therefore mass loss rate. Especially for high initial masses, when the total mass loss is large, X_C increases and thus the compactness landscape shifts to higher CO core masses (Temaj et al. 2024). Overshooting during core helium burning also influences the ratio of the total mass to the CO core mass, and thus the final BH mass. In stripped helium stars, large overshooting implies a lower BH mass (see also Fig. 6 of Schneider et al. 2023).

As we explore in more detail Appendix B, semiconvection is another mixing process that has a large effect on the final core structure of stars. It changes the compactness landscape as a function of the initial mass but only slightly affects the compactness landscape as a function of the CO core mass (see Fig. B.3). Semiconvection affects the stellar structure after core hydrogen exhaustion through the development of intermediate convection zones (Sibony et al. 2023). Furthermore, it can lead to late ingestion of helium during core helium burning depending on the chemical gradient at the edge of the helium core (Langer et al. 1985; Langer 1991) and thus not only shifts the CO core mass but also X_C (see Fig. B.2). This can explain the “noise” or “randomness” reported by Sukhbold et al. (2018) in the compactness landscape (see also Chieffi & Limongi 2020). Schootemeijer et al. (2019) found through a comparison with observations that semi-convection may be an efficient process, that is, favoring high values of α_sc ≥ 1. However, it remains a major uncertainty in our models.

7.2. Shell mergers

In our models beyond the compactness peak, we identify several models that experience mergers of their C-, Ne- and in some cases also O-burning shells. As we show in Sect. 4.5 these are the result of an increasing entropy at a nuclear burning front (e.g., Ne or O) that exceeds the entropy of the layers above, implying that these layers are unstable against convection and must mix. Several studies have reported the occurrence of such shell mergers in their one-dimensional models for a similar parameter space (i.e models before and just beyond the compactness peak, typically in the CO core mass range from 3 to 7 M_⊙, Sukhbold & Woosley 2014; Sukhbold et al. 2018; Sukhbold & Adams 2020; Collins et al. 2018; Laplace et al. 2021). We argue that these shell mergers are not artifacts from our treatment of convective boundary mixing since we do not include overshooting above shells post helium burning. In Appendix C, we verify that these shell mergers are barely affected by our choices of spatial and temporal resolution. The exact predictions of the outcomes of shell mergers will still depend on a more accurate treatment of convective mixing than the mixing-length approximation used here, that is, informed by multi-dimensional simulations (Rizzuti et al. 2023). Several three-dimensional calculations have also observed the occurrence of shell mergers, strengthening the finding from one-dimensional simulations (e.g., Couch & Ott 2013; Collins et al. 2018; Yoshida et al. 2019; Andrassy et al. 2020; Yadav et al. 2020; McNeill & Müller 2020). Collins et al. (2018) shows that shell mergers help develop stellar structures that have extended oxygen shells with high convective Mach numbers. In such stellar structure, perturbations can be induced by this convective oxygen shell burning at the time of core collapse. Three-dimensional SN simulations suggest that such perturbations are crucial for increasing the efficiency of neutrino heating in neutrino-driven SN explosions and for a successful shock revival (e.g., Couch & Ott 2013). Thus, stellar structures experiencing shell mergers may be particularly prone to successful explosions. In addition, due to the mixing process and ensuing nuclear burning episodes, stars experiencing shell mergers are expected to produce peculiar abundance patterns. Typically, these are enhanced alpha-capture products in the oxygen-shell region, which in turn affect the nucleosynthesis, in particular the Ca/O ratio of core-collapse SNe (e.g., Ritter et al. 2018; Dessart & John Hillier 2020; Laplace et al. 2021). In our work, we find that shell mergers are expected in stellar structures with strong neon shell burning. We observe that this systematically occurs in models just beyond the compactness peak (23 < M_i/ M_⊙ < 27, 6.3 < M_CO/ M_⊙ < 8.5), in which carbon burning is neutrino-dominated and core neon burning is becoming neutrino dominated.

7.3. Final structures of stars that experience binary interactions

In this work, we have focused on the core structures of single massive stars. However, it is well established that the majority of massive stars live in binary systems and are expected to interact during their lifetime (Sana et al. 2012). As a result of binary interactions, the core structure of stars can be greatly modified, and lead to a different fate even for stars with the same core mass (Laplace et al. 2021; Schneider et al. 2021, 2023, 2024). In addition, the vast majority of BHs are found in binaries observationally. In these systems their properties, in particular their masses, can be measured most accurately (Casares & Jonker 2014). It is therefore important to discuss the effect of binary interactions on the final structures of stars.

The mechanisms we identified here that set the final structure of stars (i.e., the impact of neutrino-dominated burning) apply generally for any stellar structure, as shown by our experiment in Sect. 4 and are mainly influenced by the initial conditions for core carbon burning (i.e., the CO core mass and central carbon abundance, see also Patton & Sukhbold 2020). These conditions can be influenced by binary interactions or strong wind mass loss, in particular if the structure of a star is significantly affected after core H-burning, which will ultimately change the initial conditions for core carbon burning (Brown et al. 1999, 2001; Woosley 2019; Schneider et al. 2021; Laplace et al. 2021). This implies a shift in the compactness landscape.

When it comes to BH formation, binary-stripped stars are particularly important binary products. In isolated massive binary systems, they are the natural progenitors of both BHs that eventually form binary BH mergers that are now commonly detected as GW sources (Schneider et al. 2023). After mass transfer or strong wind mass loss removes the outer envelope of these stars, the conditions under which core helium burning occur can change drastically compared to single stars. A weaker or even absent hydrogen-burning shell, and strong wind mass loss leads to a shrinking of the convective helium-burning core. This results in a higher final core carbon abundance compared to single stars for the same CO core mass (Brown et al. 2001; Schneider et al. 2021; Laplace et al. 2021). The higher core carbon abundance means that the compactness peak, and the corresponding expected parameter space for BH formation is shifted to higher core masses in these stars (M_CO ≈ 7.5 M_⊙, see Fig. 2 in Schneider et al. 2023). Therefore, binary-stripped stars are expected to be common progenitors of SNe and to only form BHs for high initial masses (of about 30 M_⊙), which significantly reduces the parameter space for BH formation (Schneider et al. 2021). Moreover, the core structures of binary-stripped stars have a weak metallicity dependence, implying the formation of universal BH masses at the compactness peaks, with strong observational signatures (Schneider et al. 2023).

The final core structures of accretors and stellar mergers can also be affected by binary interactions. In our recent work (Schneider et al. 2024) we have shown that mass accretion, in particular onto stars that have completed core hydrogen burning, can lead to systematically different final core structures. In particular, the compactness peak in case B and C accretors shifts systematic to lower core masses (Schneider et al. 2024). Contrary to binary-stripped stars, this means that accretors form BHs in stars with lower initial masses than in single stars.

7.4. Compactness peaks and BH formation

In this work, we estimate the outcome of core-collapse for our models based on the Ertl et al. (2016) criterion and on the 1D semi-analytical neutrino-driven model of Müller et al. (2016). We find a correspondence between the formation of BHs and models with high final compactness, which all have a high final binding energy, iron-core mass, and central entropy (see Fig. 1). To first order, compactness (or equivalently, any of these quantities) captures the information that stars with high binding energy are difficult to explode. Thus, independently of the uncertainties regarding the explosion mechanism of core-collapse supernovae, a signature of the characteristic final structure of stars described here is expected in their final remnants and in the properties of their explosions.

However, while compactness captures the general final structure of the star, it is not a good predictor of the detailed explosion outcome (see Sect. 3.1). Explodability criteria and the remnants expected after the core collapse are the subject of active discussion in the SN modeling community, including within groups studying the neutrino-driven explosion mechanism. For example, some 3D neutrino-driven simulations (e.g., Chan et al. 2018; Kuroda et al. 2018; Ott et al. 2018; Burrows et al. 2023) have recently found successful shock revival in models with high compactness. However, the revival of the SN shock does not necessarily imply a full SN explosion with complete ejection of the stellar layers, or that these will all form neutron stars. Due to the high accretion rate typical for models with high compactness, these are still likely to form BHs through fallback accretion, especially in stars with high core masses (Heger et al. 2023). For example, Chan et al. (2018) found that a 40 M_⊙ progenitor star with high compactness forms a BH early on through fallback accretion. In their simulation, the SN shock succeeds in unbinding parts of the envelope despite the explosion energy being lower than the binding energy of the envelope. More recently, Burrows et al. (2023) found that a 40 M_⊙ model with high compactness results in a successful explosion, though the remnant eventually forms a BH due to the high accretion rate onto the proto neutron star. BHs formed through fallback accretion have systematically lower masses than those formed through direct collapse because parts of the star can be expelled. They therefore have distinct properties, which may be constrained through observations (Chan et al. 2020). However, consistently finding explosion of models with high masses and high final compactness would be at odds with observational constraints, which seem to point to a lack of core-collapse SNe from stars with high core masses (e.g., Smartt 2009; Davies & Beasor 2020).

7.5. Observational constraints on the pre-SN structure of stars

The robustness of the patterns found for the final core structure of massive stars, that is, the “compactness peaks” based on the neutrino-dominated burning mechanism identified here, suggests that signatures of this pattern could be identifiable in observations.

Sukhbold & Adams (2020) showed that the luminosity of immediate BH progenitors could be used to test the existence of characteristic core structures corresponding to the compactness peaks observationally. The observation of the disappearing red-supergiant N6946-BH1 could be a piece of evidence (Sukhbold & Adams 2020), though more observations are needed to understand the latest constraints on this object (Beasor et al. 2024; Kochanek et al. 2024). Finding more disappearing stars or additional direct observational link between BHs and SNe, and their progenitors, will provide crucial constraints. For example, Temaj et al. (2024) showed that the luminosity of SN progenitors in hydrostatic equilibrium can be directly linked to their core mass, independently of uncertainties in interior mixing or wind mass loss.

If there is a direct link between massive star structures in the “compactness peak” and BH formation, that is, if it is not fully obscured by SN physics and fallback dynamics, then BH masses should also encode information about the final structure of stars in the compactness peak. We expect that such a link should be readily identified in binary-stripped stars, which have characteristic final structures with a weak metallicity dependence (see also Sec. 7.3). We predicted that this would lead to the formation of BHs of universal masses (Schneider et al. 2023). Alternatively, if SN and fallback physics “wash out” these signals, then identifying broad features in the BH mass distribution and significant numbers of low-mass BHs may provide further constraints into these processes and their relative importance. The complete absence of features in the BH mass distribution of binary-stripped stars (e.g., a flat BH mass distribution) would also give particularly strong constraints on stellar, binary, and SN physics, as this would require narrow ranges of physical ingredients.

Using GW observations to find features in the chirp-mass distribution of binary BH mergers could provide evidence for the existence of the characteristic final stellar structures of binary-stripped stars (Schneider et al. 2023; Disberg & Nelemans 2023). Considering a population of merging binary BHs resulting from isolated binary evolution, we predicted peaks in the chirp-mass distribution of binary BH mergers at ∼8 M_⊙ and ∼14 M_⊙, and a dearth between ∼10 − 12 M_⊙ (Schneider et al. 2023). Based on the published data from LIGO-Virgo-KAGRA (LVK) observations, a peak at a chirp mass of ∼8 M_⊙ is found to be statistically significant, while the existence of a feature at ∼14 M_⊙ is less robust (Talbot & Thrane 2018; Tiwari & Fairhurst 2021; Tiwari 2022; Abbott et al. 2023; Edelman et al. 2023; Farah et al. 2023). Tentative evidence for the predicted dearth in the distribution between 10 and 12 M_⊙ has recently been reported, though more data, even beyond the upcoming O4a data release by the LVK collaboration, is required to further test its existence (Adamcewicz et al. 2024; Galaudage & Lamberts 2025).

Signatures of the characteristic core structures of massive stars that form BHs may also be identified in the mass distribution of BHs from electromagnetic observations (Nambena et al. in prep.). Studies of BHs in X-ray binaries point to the existence of a characteristic peak in BH masses of about ∼8 M_⊙ (e.g., Kreidberg et al. 2012; Casares et al. 2017), which is similar to the mass we predict for BHs in binary-stripped stars originating from the compactness peak (Schneider et al. 2023). The recent discovery of new X-ray quiet BHs in wide binaries (Shenar et al. 2022; El-Badry et al. 2023a,b; Gaia Collaboration 2024) provides a promising avenue to discover more BHs and better understand their properties, and ultimately, their formation and progenitors.

8. Conclusions

Understanding the fate of massive stars, that is, whether stars form BHs at the end of their lives or successfully explode in a SN, is a major goal of astrophysics. In this work, we explored how stellar physics determines the pre-SN core structure and the role it plays in the fate of massive stars by studying the final evolution of single stars of 17 to 50 M_⊙ at solar metallicity. Our findings can be summarized as follows:

• We confirm that the final structure landscape of massive stars found in multiple independent one-dimensional stellar evolution studies (Timmes et al. 1996; Brown et al. 1999; Sukhbold & Woosley 2014; Sukhbold et al. 2018; Müller et al. 2016; Chieffi & Limongi 2020; Chieffi et al. 2021; Schneider et al. 2021, 2023, 2024; Takahashi et al. 2023; Temaj et al. 2024) is robust. The prominent features in the final compactness and iron core mass can be traced back to the conditions under which late nuclear burning occurs in the cores of stars in this mass range.

• Summarizing quantities such as the final compactness, central entropy, iron core mass, and binding energy are all equivalent and intrinsically linked. Stars with high final compactness all have a high binding energy, and we verify that the formation of a BH from these stars is likely based on a neutrino-driven SN model. Independently from uncertainties in the SN mechanism and explodability of stars, stellar physics thus plays a major role in determining the final remnants of stars, in particular the formation of BHs. We confirm that the final structure of massive stars is already mostly “written” in their cores at the end of core helium burning.

• We find that the final compactness increases when carbon burning becomes neutrino-dominated (region A in Fig. 1). The transition from convective to radiative core carbon burning and the change in the extent and number of carbon-burning shells identified in previous studies is a symptom of this effect, but not its cause. Strong neutrino losses force a larger core contraction and the development of a larger carbon-free core. Eventually, this enables the formation of a larger iron core. The same mechanism leads to the formation of a second compactness increase once neon burning becomes neutrino dominated (region C in Fig. 1).

• We trace the drop in compactness after the peak in compactness (regions B and D in Fig. 1) to the effect of an earlier ignition of the next nuclear fuel (neon or silicon) in these stars, which slows down the core contraction and eventually leads to the formation of smaller iron core masses. We find that the next burning phase ignites earlier in these stars because of two main mechanisms: (1) a shorter duration of the main burning stage due to a decrease in available fuel and increase in temperature, and (2) the role of electron degeneracy, which further accelerates the contraction of the core when the fuel-free core exceeds the effective Chandrasekhar mass.

• Shell mergers between the C- and Ne-burning shells, and in some cases also between the Ne- and O-burning shells, occur in our models after the first compactness peak, as observed in other studies (Sukhbold & Woosley 2014; Sukhbold et al. 2018; Collins et al. 2018). These contribute to the drop in compactness after the compactness peak, lead to smaller iron-core masses, and make these stars more explodable. We find that shell mergers are not numerical artifacts but take place because of energetic neon and oxygen burning, which generates enough entropy to mix with the layers above, independently of convective boundary mixing. Shell mergers not only affect the pre-SN structure but also change the pre-SN composition and may have observable signatures, which calls for further investigation.

• The exact final core structure landscape of stars is subject to uncertainties in stellar evolution, such as semi-convective mixing, mass loss, and uncertain nuclear reactions, in particular the ¹²C(α, γ)¹⁶O reaction rate. The origin of small variations and secondary peaks in compactness and central entropy observed when nuclear burning conditions change remains to be explored further. However, the robustness of the main features identified in the final structures of massive stars as a function of their CO core masses provides a clear theoretical prediction and an opportunity to constrain stellar physics from the observation of BHs and their immediate progenitors.

The final structure landscape discussed here applies for any massive star, including those that interact in binary systems. However, binary interactions shift the compactness landscape to different core masses, which affects the parameter space for expected BH formation and successful SNe explosions (Schneider et al. 2021, 2023, 2024). For accurate predictions of SN and BH formation, in particular for understanding the population of gravitational-wave sources, it is crucial to take the different final core structures of single and binary stars into account.

Data availability

Our final core collapse models, additional information, and figures showing the final structure evolution of all our models are available online at https://zenodo.org/doi/10.5281/zenodo.13645155.

Acknowledgments

We thank the referee for constructive comments that helped improve the manuscript. This work was funded by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant agreement No. 945806) and supported by the Klaus Tschira Stiftung. It was also supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1-390900948 (the Heidelberg STRUCTURES Excellence Cluster).

References

Appendix A: The effect of the ¹²C(α, γ)¹⁶O reaction rate

Among all nuclear reactions taking place in the interiors of stars, the ¹²C(α, γ)¹⁶O reaction rate is particularly important for determining the fate of stars (e.g., Weaver & Woosley 1993). However, it remains notoriously uncertain (Farmer et al. 2020). Towards the end of core helium burning, this reaction becomes the main helium-burning process. It determines the fraction of remaining carbon and oxygen in the core after core helium burning, as shown in Fig. A.1. For a higher reaction rate, the core carbon mass fraction X_C reaches systematically lower values. As this sets the initial conditions for the subsequent core carbon-burning stage, this reaction thus plays a crucial role in determining the final structure and fate of stars. In Fig. A.2, we show the effect of varying this reaction rate, from a ∼15% lower (Kunz et al. 2002) to a 10% higher rate than our fiducial assumptions on the final central entropy and compactness. Both as a function of initial and CO core mass, varying this reaction rate systematically shifts the location of the compactness peak to lower masses for a higher reaction rate. This is because a higher rate leads to a lower core carbon abundance and thus to carbon burning becoming neutrino-dominated at lower masses. This is thus a similar effect as described in our experiment in which we vary the core carbon abundances for models with the same core mass (Sect. 4). A difference comes from the role of the ¹²C(α, γ)¹⁶O reaction rate in setting the mass of the CO core (see the small shift in CO core mass in Fig. A.1). Secondary compactness and central entropy peaks appear at CO core masses between 5 and 6 M_⊙ and 7 to 8 M_⊙ for a varying ¹²C(α, γ)¹⁶O rate. These are discussed in more detail in Appendix B. Finally, shifts in the maximum and minimum compactness values can be observed. These can be understood as a varying degree of maximum or minimum contraction in the stellar structure, induced by the varying effect of neutrino losses for different nuclear burning conditions during core carbon burning.

Fig. A.1. Central carbon mass fraction at core helium depletion as a function of the CO core mass for varying 12C(α, γ)16O reaction rate.

Appendix B: The impact of semi-convective mixing on the final structure

Semi-convection is a particularly uncertain mixing process that occurs in stellar layers that are unstable against convection because of their superadiabaticity (Schwarzschild criterion for convection) but supported by a stabilizing molecular weight gradient (Ledoux criterion). It occurs for example at the end of core hydrogen burning when a H/He composition gradient is created by the receding convective core. This can form a transient intermediate convection zone (ICZ) and hence change the structure of the core outside the inner convective region (Sibony et al. 2023). Semi-convection remains difficult to constrain observationally, though a recent study suggest that it may be an efficient process, with an efficiency parameter of α_sc ≥ 1 (Schootemeijer et al. 2019). It has a large influence on the final core structure of stars because it affects the development of the convective helium-burning core (Langer et al. 1985; Langer 1991), generally leading to smaller core masses for a higher semi-convection efficiency (see Fig. B.1). In addition, even small semi-convective layers above a convective burning core can drastically change the initial conditions for the next burning phase. For example, late ingestion of helium due to semiconvection during core helium burning can strongly affect the core carbon abundance at core helium depletion (Langer 1991). This is shown in Fig. B.2, where the core carbon abundance X_C at the end of core helium depletion is plotted as a function of the CO core mass for varying degrees of semi-convective mixing. Between the models with no semi-convection mixing and the α_sc = 0.1 and α_sc = 10 models, there appears to be a systematic shift to a higher X_C. This can be attributed to a systematically smaller helium core mass after core hydrogen burning with respect to the initial mass, which creates a systematically smaller CO core and generates more carbon overall (see also Fig. B.1). Generally, for lower masses and higher densities, the triple alpha reaction (which depends on the density cubed) is favored during core helium burning and alpha captures onto carbon (which linearly depends on density) play less of a role (Brown et al. 2001), which means less carbon is destroyed at lower masses, leading to a higher X_C. However, this trends is no longer reproduced for the α_sc = 1 models, where X_C decreases for CO core masses larger than 5 M_⊙ compared to the α_sc = 0.1 models. This can be understood as the effect of a late ingestion of helium due to semi-convection during core helium burning in these models, which boost alpha captures onto carbon, which reduce X_C.

For the α_sc = 10 models, we observe a significant jump of ∼0.03 in X_C for models with CO core masses larger than 4.5 M_⊙ in Fig. B.2. This is because the lowest-mass models (M_i < 19 M_⊙, M_CO < 4.5 M_⊙) have a different structure compared to higher- models, which all burn helium as red supergiants. Due to the efficient semiconvection, these low-mass models develop a large ICZ after core hydrogen burning and experience a blue loop, igniting helium as blue supergiants (BSGs). The characteristic strong hydrogen-burning shell of the BSG grows a more massive helium core compared to their initial mass than is the case in RSG structures (Fig. B.1). For example, the 18.5 M_⊙ model, which develops a BSG structure, has the same He core mass of 6.2 M_⊙ as the 19 M_⊙ model which has a RSG structure. This larger helium core mass and the ensuing different burning conditions during core helium burning favor carbon destruction through alpha captures and result in a systematically lower X_C in these models.

Changes in semi-convection lead to systematic variations in the final stellar structure, as shown in Fig. B.3. The shift of the compactness and entropy peaks to higher initial masses (2 M_⊙ for a shift of α_sc from 0 to 10) reflects the change in the initial mass to helium core mass relation (see Fig. B.1). However, the CO core mass at which the compactness peak occurs is remarkably similar and varies by only 0.2 M_⊙ between models with different semi-convective mixing. This is because the onset of neutrino-dominated burning happens at a similar core mass (∼6.5 M_⊙) and similar X_C ∼ 0.23. Small variations in the compactness landscape, that is, the existence of small secondary peaks next to the main compactness peak are observed for different values of semi-convection in Fig. B.3. For example, the small peak in central entropy and compactness at an initial (CO core mass) of 20 (5.6) M_⊙ has been identified out in multiple other studies (Sukhbold & Woosley 2014; Sukhbold et al. 2018; Chieffi & Limongi 2020; Chieffi et al. 2021). This appears to be a recurrent feature in this narrow CO core mass and X_C range. Models with initial (CO core) masses between 24 and 28 M_⊙ (8 and 10 M_⊙) also show variations in compactness and final entropy for a changing semi-convection efficiency. This is the region at the transition when neon burning becomes neutrino-dominated and develops through several convective shells. We observed that for a varying semi-convection efficiency, the development of this neon-burning front is affected, which leads to small variations in compactness. However, the exact mechanisms behind these secondary peaks remains to be understood.

Uncertainties in semi-convective mixing can thus lead to variations in the final stellar structure by influencing the initial conditions for core helium and carbon burning. This creates in part the“noise” observed by Sukhbold et al. 2018 in their study of the compactness of massive stars. However, it is noteworthy that semi-convection does not affect the existence of the main features observed in the compactness landscape and does not significantly change the CO core mass range at which the main compactness peak occurs.

Appendix C: Resolution test for the occurrence of shell mergers

To test how the occurrence of shell mergers is affected by numerical uncertainties, we perform a resolution test for the model presented in Sect. 4.5, that is, our default 22 M_⊙ model with a modified central carbon abundance X_C = 0.17, starting from core oxygen exhaustion. We double the spatial and temporal resolution, and in pariticular, double the number of grid points in zones containing composition gradients and nuclear burning regions of C, O, and Ne (MESA settings mesh_logX_species and mesh_dlog_burn_c_dlogP_extra and equivalent). In Fig. C.1 we show the comparison of the Kippenhahn diagrams of the default model and the model with double the resolution. We find no difference in the occurrence of the shell merger between the C and Ne shell (highlighted by the red box). The timing and extent is nearly identical, and can be traced back to the same mechanism, that is, the increase in entropy in the Ne-burning shell. We therefore verify that these shell mergers between the C and Ne shell do not originate from poor resolution. The second shell merger between the C/Ne and the O-burning shells (highlighted by the blue box in Fig. C.1) occurs 5 min earlier for the model with higher resolution. Since this shell merger occurs towards the end of Si shell-burning, after most of the iron-rich core has already formed, it barely affects the final core structure. Numerical uncertainties mainly influence the timing of this shell merger event.

Fig. C.1. Final central entropy (top) and compactness (bottom) as a function of the core mass (left) and initial mass (right) for a varying 12C(α, γ)16O rate.

In Fig. C.2 we compare the final density profiles at the onset of core collapse. The density profiles are almost identical, except for small variations of order ρ ∼ 10g cm⁻³ in the He-rich layer. The final iron core mass is a little larger for the higher-resolution model (M_Fe core = 1.80 compared to M_Fe core = 1.89). Similarly, the final compactness values is a little larger (ξ_2.5 = 0.357 for the default model and ξ_2.5 = 0.364 for the models with double the resolution). We conclude that the shell mergers identified in this work are barely affected by numerics. The C and Ne shell mergers are not numerical artifacts, but instead originate from the high entropy contrast between the Ne-burning shell and the C shell in this model (see Sec. 4.5).

Fig. C.2. He (full lines) and CO (dashed lines) core mass at core helium depletion for models with a varying semiconvection efficiency. A higher efficiency generally leads to a smaller core mass.

Fig. C.3. Central carbon mass fraction at core helium depletion as a function of the CO core mass for varying semi-convection efficiency.

Fig. C.4. Final central entropy and compactness as a function of the initial mass for varying semi-convection efficiency.

Fig. C.5. (a) Same as Fig. 8 but from core oxygen depletion to core collapse. (b) The same model in which we double the spatial and temporal resolution. The red and blue boxes highlight the shell merger between the C and Ne-burning shells, and the C/Ne and O-burning shells, respectively.

Fig. C.6. Final density profile of the models in Fig. C.1 with our default assumptions compared to a model with double the resolution at the onset of core collapse. From left to right, vertical lines indicate the Fe, Si, CO, and He core masses in black and orange for the default and higher resolution model, respectively.

All Figures

Fig. 1. Final (a) compactness, (b) specific central entropy, (c) iron core mass, and (d) binding energy above M4 at the onset of core collapse as a function of the initial mass. The top axis shows the CO core mass at core helium exhaustion. Circles and crosses represent BH formation and explosions, respectively, according to the Ertl et al. (2016) criterion, while black and white colors indicate BH formation and explosions based on the Müller et al. (2016) SN model.

In the text

Fig. 2. Evolution of central specific entropy as a function of the initial mass for different evolutionary stages from core hydrogen depletion until core collapse. Key differences in the final entropy start to appear at the end of core carbon burning. Red markers indicate three models computed with a nuclear network comprising 128 isotopes. Taking into account a greater number of nuclear reactions generally leads to a slightly lower final central entropy, though the general qualitative pattern of the compactness peak is unchanged. The characteristic final landscape is already mostly determined at the end of core Ne burning and remains until core collapse.

In the text

Fig. 3. Central temperature as a function of the central density for models with the same core mass and different carbon mass fractions XC after core helium burning (indicated in the color bar). The previous evolution of the core (not shown here for clarity) is indistinguishable for all models, while the post-core-helium-burning evolution is greatly affected by the differences in central carbon abundance. Symbols indicate key evolutionary phases. To the right of the gray dashed line, electron degeneracy dominates the pressure. The inset axis indicates the final compactness of the models as a function of decreasing XC.

In the text

Fig. 4. Evolution of key elements of the inner 6 M⊙ core structure as a function of time from carbon ignition to core collapse for models with the same core mass and a different central carbon abundances XC. (a) Kippenhahn diagram with selected elements. Full and dashed lines indicate the carbon and neon-burning fronts, respectively. Shaded (hatched) regions show convective carbon (neon) burning regions in the core. The iron-core masses are indicated with dash-dotted lines. The gray horizontal dashed line indicates the classical Chandrasekhar limit. The inset figure shows the final compactness of these models as a function of decreasing XC. (b) Total mass of 12C in the convective carbon-burning shells for each of the models as a function of time.

In the text

Fig. 5. Time evolution of the ratio between the specific nuclear energy generation rate ϵnuc and the neutrino loss rate ϵν at the location of the C, Ne, and O-burning fronts for models with the same core mass and different central carbon abundances. The dashed horizontal line indicates where ϵnuc = ϵν. When the energy ratio significantly exceeds this line, convection occurs (Sukhbold & Adams 2020). With a decreasing core carbon abundance, the core becomes more neutrino dominated and neon burning, followed by oxygen burning, occur earlier. As highlighted in the red boxes, this temporarily suppresses carbon burning.

In the text

Fig. 6. Properties of a series of models with the same core mass and a varying initial core carbon abundance XC, indicated with colors. (a) Duration of core carbon burning. (b) Final compactness as a function of XC. Models that undergo convective core carbon burning are highlighted with black outlines. (c) Duration of the phase between core carbon depletion and core neon ignition as a function of mC, the maximum mass coordinate reached by the carbon-burning front before core neon burning. The classical (MCh, 0) and effective (MCh) Chandrasekhar mass at core C depletion are shown as vertical dashed and full lines, respectively.

In the text

Fig. 7. Kippenhahn diagram of the inner 6 M⊙ core structure of models with the same core mass and different core carbon abundance at core C ignition. Colors indicate the dimensionless electron degeneracy parameter η = μ/kBT (electrons are partially degenerate for η ≈ 0 and very degenerate for η ≫ 0). Convective zones are highlighted by the hatched regions and the dotted vertical lines indicate, from left to right, the moments when 12C, 20Ne, 16O, and 28Si are depleted at the center (central abundance lower than 10−4). The gray full and dashed lines indicate the classical (MCh, 0) and effective (MCh) Chandrasekhar mass, respectively. Electron degeneracy during core carbon burning increases for models with lower central carbon abundances.

In the text

Fig. 8. Final evolution of the inner stellar structure of the model with XC = 0.17 which experiences a shell merger between the C and Ne-burning shells (highlighted in the red box in panel a). (a) Kippenhahn diagram in which colors indicate the regions dominated by nuclear burning or neutrino losses and gray hatched regions indicate convective areas. Brown lines show contours of constant (logarithmic) density. Colored markers trace the burning front for C, Ne, and O burning. Small colored vertical lines indicate the times at which the entropy profiles are shown in the bottom panel. (b) Entropy profiles in the mass and time range highlighted in the red box in the top panel. The entropy in the neon-burning region (2.5 − 2.9 M⊙) gradually increases until it exceeds the entropy of the C-burning layers above, triggering a merger between these regions.

In the text

Fig. 9. Evolution of the C-burning front, the C-free core (which follows the progression of neon burning), and the iron-rich core as a function of time from carbon burning to core collapse for models before and up to the compactness peak (region A of Fig. 1, see inset). As carbon burning becomes increasingly neutrino-dominated, the C-burning front moves further out in mass, ultimately leading to the formation of a larger C-free and iron-rich core.

In the text

Fig. 10. Same as Fig. 9 but for models beyond the compactness peak (region B in Fig. 1, see inset). As a function of initial mass, neon burning ignites systematically earlier (as traced by the C-free core mass) and slows down the progression of the C-burning front. Shell mergers between the C/Ne/O-burning fronts are responsible for the sudden drops in the C-burning front at τ ≈ −5. Ultimately, for increasing masses, these models form systematically smaller iron-rich cores.

In the text

	Fig. 11. Same as Fig. 7 but for three stellar models of our grid after the first compactness peak.
In the text

	Fig. 12. Same as Fig. 9 but for models in the mass range of the second compactness increase (region C in Fig. 1, see inset). In these models, both carbon and neon burning are fully neutrino-dominated.
In the text

	Fig. 13. Same as Fig. 9 but for models in the mass range of the second decrease in compactness (region D in Fig. 1, see inset). The dashed line traces the O-burning fronts. For increasing masses, the progression of the O-burning front is slowed by a systematically earlier ignition of core silicon burning (as traced by the iron core mass).
In the text

Fig. 14. Same as Fig. 7 but for three stellar models of our grid after the second compactness increase. Here, neon burning is fully neutrino-dominated and radiative, and not visible in the figure. After convective core oxygen burning (hatched region in the center) the degenerate core region reaches up to the oxygen burning front (at the base of the convective oxygen-burning shell) for all models. Panel c: When this front exceeds the (effective) Chandrasekhar mass, the degenerate core contracts and core silicon burning ignites earlier, slowing the contraction of the layers above and affecting the extent of the convective region.

In the text

Fig. 15. Schematic representation of the mechanisms leading to the emergence and decline of peaks in final iron core mass and compactness. Horizontal gray dotted guidelines are included to aid the comparison. Star I – Before peak: Schematic stellar structures of a star with a low final iron core mass at two points in the evolution represented by two adjacent wedges. The initial mass fraction of the fuel is still high (see red area in the rectangle) but the main burning stage is neutrino-dominated (arrows). Because neutrinos take away energy, the core contracts and the burning front (red line) moves out in mass as the burning progresses. The large amount of fuel prevents it from moving far out before a convective zone forms above the burning front once the energy generated is high (shaded red area). Ultimately (right wedge), it grows a small fuel-free core and the silicon-burning front (yellow line) cannot move far out in mass, forming a small iron core (black area). Star II – Peak: With a higher mass and lower initial fuel mass fraction, the burning is even more neutrino dominated and a strong contraction occurs. Hence the burning front moves further out in mass than in I, but stays just below the effective Chandrasekhar mass. Aided by degeneracy pressure support, the star burns through almost all available fuel before the next burning stage ignites, after which it quickly moves out in mass, growing a large fuel-free core and eventually a large iron core. Star III – Beyond peak: For an even lower initial fuel abundance and higher mass, the burning is even more neutrino dominated and the core contraction accelerates. The burning front moves further out and exceeds the effective Chandrasekhar mass. This triggers a fast contraction of the partially degenerate core and the next burning stage (blue area) ignites early. This next stage suppresses nuclear burning at the front above. As a result, the burning front moves out in mass slowly and eventually (right wedge) the star grows a low-mass iron core.

In the text

	Fig. A.1. Central carbon mass fraction at core helium depletion as a function of the CO core mass for varying 12C(α, γ)16O reaction rate.
In the text

	Fig. C.1. Final central entropy (top) and compactness (bottom) as a function of the core mass (left) and initial mass (right) for a varying 12C(α, γ)16O rate.
In the text

	Fig. C.2. He (full lines) and CO (dashed lines) core mass at core helium depletion for models with a varying semiconvection efficiency. A higher efficiency generally leads to a smaller core mass.
In the text

	Fig. C.3. Central carbon mass fraction at core helium depletion as a function of the CO core mass for varying semi-convection efficiency.
In the text

	Fig. C.4. Final central entropy and compactness as a function of the initial mass for varying semi-convection efficiency.
In the text

	Fig. C.5. (a) Same as Fig. 8 but from core oxygen depletion to core collapse. (b) The same model in which we double the spatial and temporal resolution. The red and blue boxes highlight the shell merger between the C and Ne-burning shells, and the C/Ne and O-burning shells, respectively.
In the text

	Fig. C.6. Final density profile of the models in Fig. C.1 with our default assumptions compared to a model with double the resolution at the onset of core collapse. From left to right, vertical lines indicate the Fe, Si, CO, and He core masses in black and orange for the default and higher resolution model, respectively.
In the text