© The Authors 2025

1. Introduction

Massive stars are keystones of the chemical evolution of the Universe, synthesising heavy elements in their interiors and enriching their surroundings through Type II supernovae explosions (e.g. Woosley et al. 2002). The supernova ejecta composition reflects the cumulative effects of successive nuclearburning stages and extensive nucleosynthesis during massive star evolution (e.g. Bennett et al. 2012; Chieffi & Limongi 2013). A detailed understanding of the stellar evolution phases from hydrogen to silicon burning and their role in driving nucleosynthesis remain incomplete, especially for the most advanced phases, i.e. beyond He-exhaustion. The specific density and temperature conditions reached in these phases enable a broad range of nuclear reactions to occur, including s-process nucleosynthesis, contributing to heavy element formation. These reactions also influence the stellar fate by determining the remnant type and mass as well as the composition of the supernova ejecta (e.g. Chieffi & Limongi 2020). Only the most energetic reactions involving the most abundant elements significantly drivestellar evolution and impact the chemical structure. Among these, nuclear reactions involving carbon and oxygen are key to massive star evolution (e.g. Woosley 1986). Ongoing experimental and theoretical efforts aim to determine the rates of crucial reactions at the relatively low energies of astrophysical interest with improved accuracy and reliability (e.g. Broggini et al. 2018).

Four reactions involving carbon and oxygen are of particular interest in the context of advanced evolution. First, the ¹²C(α, γ)¹⁶O reaction, primarily active during the latter part of He-burning phase, governs the ¹²C/¹⁶O ratio in stars and potentially affects the CO core mass at the end of the He-burning phase. In particular, deBoer et al. (2017) have shown that lower reaction rates lead to a higher ¹²C/¹⁶O ratio with consequences on stellar evolution and structure, eventually influencing the explosiveness (see also e.g. Sukhbold & Adams 2020; Xin et al. 2025). The ¹²C/¹⁶O ratio is ultimately propagated throughout the Universe when a massive star explodes, ejecting the carbon and oxygen from its outer layers, which then contribute to the formation of new stars. Second, the fusion reaction ¹²C+¹²​C, the most important reaction during the C-burning phase, impacts both the nucleosynthesis and the final fate of stars (e.g. Bennett et al. 2012; Pignatari et al. 2013; Frischknecht et al. 2016; Limongi & Chieffi 2018; Chieffi et al. 2021; Monpribat et al. 2022; Dumont et al. 2024). Third, the fusion reaction ¹⁶O+¹⁶​O is the main reaction during the O-burning phase. Finally, the fusion reaction ¹²C+¹⁶​O is expected to occur at the transition between carbon and oxygen fusion, mainly during shell burning. This reaction likely contributes to heavy element production from the end of the core C-burning phase to the core O-burning phase.

Stellar evolution models rely on nuclear physics, with both the choice of isotopes and the adopted reaction network strongly influencing stellar evolution and nucleosynthesis. Nuclear reaction rates are obtained from experimental cross-section measurements and theoretical predictions, which are both subject to uncertainties. At the low energies typical of stellar interiors (a few mega electron-volts), nuclear reaction cross-sections are extremely small, making precise and reliable rate measurements particularly challenging. Consequently, variations in these rates propagate into uncertainties in stellar model predictions. While widely used nuclear rate compilations include Caughlan & Fowler (1988, hereafter CF88), Angulo et al. (1999, Nacre), and Xu et al. (2013, NacreII), newer rates are available, offering the potential to improve our modelling of stellar nuclear physics processes. The aim of the this paper is to explore the newly available reaction rates with regard to their potential impact on the advanced stages of stellar evolution, nucleosynthesis, and the ultimate fate of massive stars, with a particular focus on chemical abundances and stellar structure.

In a previous paper, Dumont et al. (2024, hereafter D24) examined the impact of updated nuclear reaction rates for ¹²C+¹²​C during the core C-burning phase. Here, we extend this analysis by incorporating revised nuclear rates for the ¹²C(α, γ)¹⁶O reaction. We also expand our investigation to include the O-burning phase, exploring the consequences of adopting new reaction rates based on the Hindrance scenario or on theoretical nuclear model predictions for the ¹²C+¹⁶​O and ¹⁶O+¹⁶​O reactions¹.

The outline of this paper is as follows. In Sect. 2, we describe the input physics included in the Geneva stellar evolution code (GENEC). In Sect. 3, we describe the nuclear rates that we are testing in our models. In Sect. 4, we present the impact on stellar evolution from carbon burning up until oxygen burning. Sect. 5 presents the resulting heavy-element nucleosynthesis, with particular emphasis on the effects on the s-process. In Sect. 6, we asses the implications of our models for the final fate of massive stars. Finally, we discuss our results and conclude in Sect. 7.

2. The stellar evolution code

We modelled the evolution of massive stars using the stellar evolution code GENEC as described in Eggenberger et al. (2008), Ekström et al. (2012), Griffiths et al. (2025) and D24. We computed both rotating and non-rotating models at solar metallicity (Z = 0.014) for three different initial masses: 15, 17, and 20 M_⊙. Each model was evolved from the zero age main-sequence (ZAMS) to O-exhaustion, defined as when the oxygen core abundance X_¹⁶O drops below 10⁻³, unless otherwise stated.

The latest version of the GENEC code was developed to better capture the advanced phases of massive star evolution (Griffiths et al. 2025). In this version, radiative opacities are generated within the OPAL opacity tables from Iglesias & Rogers (1996) for temperatures lower than 1 Gk. At higher temperatures, the estimates from Poutanen (2017) are used. The equation of state (EoS) combines a low temperature and density regime, which includes partial ionisation effects, with the fully ionised plasma EoS of Timmes & Swesty (2000). These two regimes are matched at conditions log(T) > 7.55 and log(ρ) > 2.8. For a typical 25 M_⊙ star, this corresponds to a transition between the two EOSs in the core during the He-burning phase. We used the grey atmosphere approximation and treated convection according to mixing length theory, with α_MLT = 1.6, calibrated to the Sun with the solar reference abundances from Asplund et al. (2005) and Cunha et al. (2006), as in D24. We used the Schwarzschild criteria for convective stability and applied core overshooting for phases prior to C-ignition with an overshoot parameter d_over = 0.1 Hp, where Hp is the pressure scale height. In the present paper, we use the same formalism for rotation as described in D24 with V_ini/V_crit = 0.4 and follow the shellular rotation hypothesis by Zahn (1992), Maeder & Zahn (1998, 2004). For the sake of comparison with D24, we do not take into account the action of magnetic fields. However, magnetic torques are drivers of angular momentum transport in stellar evolution. The combined impact of the new rates and magnetic angular momentum transport will be the object of a future work.

3. Nuclear reaction network and rates

The nuclear network used in GENEC during the advanced evolutionary phases has been updated from the version employed in D24 to include the primary reactions that drive the energy generation and chemical evolution during C-, Ne-, and O-burning. The current network, GeValNet48, comprises 48 species. It builds on the spine network GeValNet25 introduced in Griffiths et al. (2025), with the addition of secondary alpha chains to enable the tracking of reactions such as ${}^{16}\mathrm{O}{{(}^{16}\mathrm{O},\mathrm{p})}^{31}\mathrm{P}$. The full list of explicitly tracked isotopes is detailed in Appendix A. Furthermore, we have updated the reaction rates of seven key nuclear reactions that significantly affect the most abundant elements and the neutron density following the recommendations of the Netgen database² and as identified by Choplin et al. (2018) (see Sect. 3.3). For reactions not directly involved in our study, rates were taken from the NACRE II database (Xu et al. 2013) or from CF88 when NACRE II data were not available. The main focus of the present study is to explore the impact and sensitivity of the measured or theoretically predicted reaction rates (listed in Table 1) for key nuclear reactions involving ¹²C and ¹⁶O. These include the ¹²C(α, γ)¹⁶O reaction, which plays a critical role during He-burning as well as the ¹²C+¹²​C, ¹²C+¹⁶​O, and ¹⁶O+¹⁶​O fusion reactions that occur during the later stages of stellar evolution.

3.1. ¹²C(α,γ)¹⁶ O

One of the most important nuclear reactions during stellar evolution is ¹²C(α, γ)¹⁶O. This reaction primarily occurs during the He-burning phase, in competition with the 3α process, and drives the ¹²C/¹⁶O ratio at the end of this phase. In the present work, we adopt the updated reaction rates from deBoer et al. (2017) instead of those from Kunz et al. (2002), which were used in D24. The rates provided by deBoer et al. (2017) are more accurate and have smaller uncertainties, as illustrated in their Fig. 29. Figure 1 compares the two sets of rates and shows the predicted He-burning temperatures in both the core and the shell. In the region of interest, the deBoer et al. (2017) rates are systematically lower than those of Kunz et al. (2002). This revision results in a slightly shorter He-burning lifetime (by ∼2%)³, an approximately 15% increase in the C-burning lifetime, shorter Ne- and O-burning lifetimes (according to internal report), and a higher ¹²C/¹⁶O ratio at the end of He-burning, along with a slightly hotter and denser core.

Fig. 1. Reaction rates for the four main nuclear reactions involving 12C and 16O. Top left: Rates from deBoer et al. (2017) normalised to Kunz et al. (2002) reference for the 12C(α, γ)16O reaction. The grey hatched areas indicate the temperature ranges where core and shell He-burning occur in the GENEC rotating model explored mass range. For the rest of the panels, nuclear reaction rates are taken from CF88 (black), HIN(RES (blue), and TDHF (magenta), all normalised to the CF88 rates (NA < σv> CF88). Shaded areas (with corresponding colour codes) indicate the temperature ranges where core and shell C-burning and core O-burning take place in the 15 and 17 M⊙ GENEC rotating models.

According to Patton & Sukhbold (2020), both the CO core mass and the central ¹²C mass fraction at the end of He-burning influence the final fate of massive non-rotating stars. Variations in these parameters between models can thus affect the nature of the compact remnant produced at the end of stellar evolution. This has been recently confirmed, for example, by Xin et al. (2025) (see also references therein), who examined the impact of an arbitrary ±3σ variation in reaction rates, based on the Kunz et al. (2002) data, and showed, using a different methodology than the one employed here, that such variations influence both the mass and type of remnant.

In our models, we find that the CO core mass at the end of He-burning is largely insensitive to changes in the ¹²C(α, γ)¹⁶O reaction rates. For instance, we obtain M_CO ≈ 2.35 M_⊙, 2.92 M_⊙, and 4.10 M_⊙ for our 15, 17, and 20 M_⊙ models, respectively.

However, the ¹²C/¹⁶O ratio is affected, increasing from ¹²C/¹⁶O ≈ 0.55 with the Kunz et al. (2002) rates to ¹²C/¹⁶O ≈ 0.65 with the updated rates from deBoer et al. (2017), a relative increase of about 20%. This is expected, as the lower ¹²C(α, γ)¹⁶O rate leads to fewer α captures by ¹²C, resulting in more ¹²C remaining and less ¹⁶O being synthesised. The small change observed in the CO core mass is also consistent with expectations, as this quantity is primarily determined by the convective core mass during He-burning, which in turn depends on the total energy released by nuclear reactions. Since the ¹²C(α, γ)¹⁶O reaction does not dominate the energy generation in this phase, the overall energy output remains nearly unchanged between the two models. We will further explore the implications for the final fate of these models in Sect. 6.

3.2. Fusion reactions ¹²C+¹²C, ¹²C+¹⁶O and ¹⁶O+¹⁶O

In Monpribat et al. (2022), hereafter M22, and D24, it was shown that recent measurements for the ¹²C+¹²​C reaction exhibit significant deviations from the commonly used CF88 rates. These updated rates account for a suppression trend (or hindrance) of the fusion cross-section at low energies, and also incorporate a resonance at 2.14 MeV, first observed by Spillane et al. (2007) and later supported by STELLA measurements (Fruet et al. 2020). We previously demonstrated that the adoption of these new rates yields a shorter C-burning lifetime, a hotter and denser burning core, and a shift in the transition from a convective C-burning core to a radiative C-burning core to higher initial masses (approximately 29 and 25 M_⊙ for rotating and non-rotating models, respectively, compared to about 27 and 25 M_⊙ with the CF88 rates). Here, we extend the analysis beyond core C-burning to explore the impact of these updated rates on the C-shell burning phase, as well as on core Ne- and core O-burning phases. In these later stages, the ¹⁶O+¹⁶​O fusion reaction becomes increasingly relevant, along with a potential contribution from the ¹²C+¹⁶​O reaction. Given their similarity to the ¹²C+¹²​C reaction, we adopt a consistent modelling approach for all three reactions and compare the resulting outcomes. We focus mainly on three representative stellar models, as described below. In the absence of experimental data within the relevant Gamow windows for these reactions, this study also serves as a preliminary sensitivity study aimed at guiding future experimental efforts.

As the first energetic reaction of interest for massive stars after He-burning, carbon fusion takes place for temperatures of at least 0.6–0.8 GK:

$$\begin{array}{c}\hfill {}^{12}\mathrm{C}{+}^{12}\mathrm{C}{\to }^{24}{\mathrm{Mg}}^{\ast }\left\{\begin{array}{cc}{\to }^{20}\mathrm{Ne}+\alpha \hfill & (\mathrm{Q}=4.62\mathrm{MeV})\hfill \\ {\to }^{23}\mathrm{Na}+\mathrm{p}\hfill & (\mathrm{Q}=2.24\mathrm{MeV})\hfill \\ {\to }^{23}\mathrm{Mg}+\mathrm{n}\hfill & (\mathrm{Q}=-2.62\mathrm{MeV})\hfill \end{array}\right\},\end{array}$$(1)

where Q (or Q-value) refers to the energy released or consumed by the reaction. For this reaction in particular, substantial experimental efforts and advances in measurement techniques have enabled direct measurements to probe increasingly lower energies. However, the available data are still suffering from significant uncertainties at the centre-of-mass at energies below 2 MeV (Spillane et al. 2007; Jiang et al. 2018; Fruet et al. 2020). To overcome the challenges posed by low-energy measurements, alternative indirect methods, such as the so-called Trojan Horse Method (THM) (Mukhamedzhanov et al. 2008), have been developed. This technique provides indirect estimates of the ¹²C+¹²​C reaction cross-section, which are several orders of magnitude higher than those obtained via direct measurements (Tumino et al. 2018). While these data have been employed in stellar evolution studies (Chieffi et al. 2021; Xin et al. 2025), their interpretation is remains contentious within the community (see e.g. Mukhamedzhanov et al. 2019), further work is needed to reconcile the results of these two complementary approaches. In the present paper, we focus exclusively on the impact of the new rates obtained from direct measurements and theoretical nuclear model predictions, within the temperature range relevant for core and shell C-burning phases.

At temperatures above 1 GK, during both hydrostatic and explosive C- and O-burning phases, when the ¹²C+¹²​C fusion reaction becomes subdominant, the following nuclear reaction becomes relevant:

$$\begin{array}{c}\hfill {}^{12}\mathrm{C}{+}^{16}\mathrm{O}{\to }^{28}{\mathrm{Si}}^{\ast }\left\{\begin{array}{cc}{\to }^{24}\mathrm{Mg}+\alpha \hfill & \text{(Q}=6.77\text{MeV)}\hfill \\ {\to }^{27}\mathrm{Al}+\mathrm{p}\hfill & \text{(Q}=5.17\text{MeV)}\hfill \\ {\to }^{27}\mathrm{Si}+\mathrm{n}\hfill & (\mathrm{Q}=-0.42\mathrm{MeV})\hfill \end{array}\right\}.\end{array}$$(2)

Suitable conditions for reaction (2) occur when the ¹⁶O abundance is high enough (typically between core C-exhaustion and the onset of O-burning, as well as in the shell C-burning regions). Although this reaction has a lower probability to occur than reactions (1) and (3), its main impact lies not in the energy production but in its potential to influence nucleosynthesis yields during the C- and O-burning phases. Experimentally, this reaction has recently been revisited (e.g. Fang et al. 2017) at low energies of astrophysical relevance, enriching the dataset that has been available since the 70s. These new measurements provide modest evidence for the existence of resonances in the deep sub-Coulomb regime.

The last major reaction for which we vary the rates is the dominant nuclear process during the O-burning phase, occurring at temperatures of about 1.9 GK (e.g. Maeder 2009):

$$\begin{array}{c}\hfill {}^{16}\mathrm{O}{+}^{16}\mathrm{O}{\to }^{32}{\mathrm{S}}^{\ast }\left\{\begin{array}{cc}{\to }^{28}\mathrm{Si}+\alpha \hfill & \text{(Q}=9.59\text{MeV)}\hfill \\ {\to }^{31}\mathrm{P}+\mathrm{p}\hfill & \text{(Q}=7.68\text{MeV)}\hfill \\ {\to }^{31}\mathrm{S}+\mathrm{n}\hfill & \text{(Q}=1.45\text{MeV)}\hfill \end{array}\right\}.\end{array}$$(3)

Due to the experimental difficulties to measure this reaction accurately, mainly caused by the technical availability of oxygen targets and the large number of exit channels, only scarce and sometimes disagreeing measurements from the 70s and 80s are currently available (Fernandez et al. 1978; Tserruya et al. 1978; Kuronen et al. 1987, and see for instance Fig. 1 by Simenel et al. 2013).

Based on the available experimental data, several models with a small number of parameters have been developed to describe the astrophysical S-factors and derive the corresponding nuclear reaction rates. The standard and most widely used reference is CF88, which was originally implemented in the GENEC code and remains a common choice for stellar evolution modelling. For fusion reactions in particular, a phenomenological model deviating from CF88 at very low energies, the Hindrance (or HIN model), has recently been incorporated in GENEC. Both models are explored in this paper, and their implementation is described in detail in Sect. 3.2.1. Alternatively, modern theoretical approaches can provide insight into the behaviour of the cross-sections at very low energies. These methods help constrain the rate trends or, in the absence of experimental data, offer predictive estimates that are used in stellar evolution codes (such as those available in the REACLIB database⁴).

For the three fusion reactions that are discussed in the present work (reactions 1, 2, and 3), multiple frameworks have been used in the literature: the Time-Dependent Hartree-Fock method (hereafter TDHF, e.g. Bonche et al. 1978; Simenel et al. 2013; Simenel & Umar 2018; Godbey et al. 2019; Scamps & Hashimoto 2019), the Antisymmetrised Molecular Dynamics method (Taniguchi & Kimura 2020, 2021, hereafter AMD), and the microscopic α-cluster model (e.g. Depastas et al. 2023). However, these models are subject to limitations arising from the approximations and assumptions needed to perform the computations. For example, TDHF and microscopic α-cluster models are unable to capture resonances in the cross-section and may omit relevant physical phenomena, especially for complex systems such as ¹²C+¹²​C. A combination of direct measurements of nuclear reactions and model-based constraints on the energy dependence of the reaction rates is therefore necessary to advance towards a more accurate description of the nuclear reaction rates. The approach followed in this work is to compare the empirical (data-driven) models with the TDHF model, which is described in more detail in Sect. 3.2.2.

3.2.1. Data-driven frameworks

We use as reference rates for the three fusion reactions the parametrised analytic formulae from CF88, a standard reference widely used in the literature. In addition, we also employ the rates from the phenomenological hindrance scenario. This scenario accounts for the suppression in the fusion cross-sections at deep sub-barrier energies, a phenomenon observed in medium-mass fusion systems (hereafter simply referred to as systems) (e.g. Jiang et al. 2007). In Q < 0 systems, the hindrance effect is clearly observed as a more rapid decline of the cross-section than expected, suggesting the existence of a maximum in the S-factor. Whether or not a similar suppression occurs in Q > 0 systems is still unclear. Although this trend has been experimentally investigated across a wide range of systems, it remains under debate whether such suppression also affects light fusion systems relevant to our study.

Assuming this interpretation is correct, Jiang et al. (2007) derived an analytical expression for the S-factor, fitted to experimental data across a broad range of systems. However, a complete physical explanation for the hindrance effect remains elusive. Adding a repulsive contribution to the ion-ion potential enables the reproduction of the observed behaviour. This repulsive term has been suggested to arise, at least partially, from a proper treatment of the Pauli repulsion (Simenel et al. 2017; Umar et al. 2021).

Motivated by experimental constraints, we adopt the more realistic Hindrance+Resonance model proposed by M22 for the case of ¹²C+¹²​C. To describe the hindrance scenario for the ¹²C+¹⁶​O and ¹⁶O+¹⁶​O systems, we use the phenomenological formulation of the S-factor maximum provided by Jiang et al. (2007). Starting from these expressions for the total S-factor, we derive the corresponding non-resonant nuclear reaction rates for each channel (proton, alpha and neutron). The branching ratios for each channel are assumed to follow those of CF88, except for ¹²C+¹²​C, where experimentally derived branching ratios of 65–35% are adopted in favour of the α exit channel. We refer to this approach as the HIN(RES) scenario, which differs from the full Hindrance model (as explored in M22 and D24) by including the resonance only for the ¹²C+¹²​C reaction given the lack of clear experimental resonance data within the Gamow window for ¹²C+¹⁶​O and ¹⁶O+¹⁶​O. Figure 1 shows the relative deviations from CF88 of the rates obtained under the HIN(RES) scenario as a function of temperature. The C- and O-burning temperature regimes are given for the 15 and 17 M_⊙ stellar models. For these masses, we expect the hindrance to impact the ¹²C+¹⁶​O reaction while its influence on the ¹⁶O+¹⁶​O reaction is expected to be minor, since the suppression in the cross-section predominantly occurs at temperatures below those relevant for core oxygen burning.

3.2.2. Time-dependent Hartree-Fock framework

As a complementary approach to the previous data-driven models we choose to explore a third scenario: the impact of nuclear reaction rates derived from the fully microscopic TDHF mean-field theory (for a detailed description of this method see e.g. Simenel 2025). TDHF simulations describe the time evolution of the occupied nucleon wave functions in the self-consistent mean field generated by the ensemble of nucleons. The many-body state is constrained to be an antisymmetrised product of single particle states (i.e. a Slater determinant) at all times. The main input is the energy density functional (EDF), which characterises the effective nucleon-nucleon interaction. Although most applications are based on the Skyrme EDF, recent calculations with the finite range Gogny interaction and including pairing correlations have also been performed by Scamps & Hashimoto (2019).

A significant advantage of the TDHF approach is that it describes both nuclear structure and dynamics on the same footing. In particular, EDF parameters are only fitted to nuclear structure data and equation of state properties, without relying on input from nuclear reactions. Modern TDHF calculations can now be performed in full 3 dimensions, without imposing constraints on the systems shape evolution. The weakness of the real-time mean-field method is that it cannot account for the quantum tunnelling expected in the sub-barrier fusion. To overcome this, predictions are obtained by calculating the nucleus-nucleus potential at an energy just above the Coulomb barrier using the so-called density-constrained TDHF method (DC-TDHF: Umar & Oberacker 2006). The resulting potential is then used to compute barrier transmission probabilities at sub-barrier energies through one-barrier penetration methods. Importantly, DC-TDHF potentials account for crucial dynamical effects such as nucleon transfer (Godbey et al. 2017), which may significantly impact sub-barrier fusion cross-sections. DC-TDHF calculations were performed in Godbey et al. (2019) for the case of ¹²C+¹²​C system, which is also considered in the present work. Here, calculations based on the SLy4d (Kim et al. 1997), UNEDF1 (Kortelainen et al. 2012), and QMC (Stone et al. 2016) EDFs are used. The trends obtained are very similar for the three models and show an increase of the S-factor at low energies, in contrast to the Hindrance scenario (see Appendix B), though still several orders of magnitude lower than the predictions from THM discussed above. For our sensitivity study, we choose to adopt the predictions obtained with the Skyrme SLy4d functional. For the two other systems (¹²C+¹⁶​O and ¹⁶O+¹⁶​O), we perform analogous calculations using the same methodology as in Godbey et al. (2019).

The resulting nuclear rates are shown in Fig. 1 and compared to the CF88 and HIN(RES) rates. At the temperatures where core and shell C- and O-burning are predicted to occur, the TDHF rates are slightly higher than those of CF88 by less than 1 dex. In contrast, the Hindrance scenario yields significantly lower rates (up to 2 orders of magnitude at these temperatures) suggesting it may lead to the largest differences in the results.

3.3. Neutron source nuclear reactions

As the main contributor to free neutrons during massive star evolution, the rate of the ²²Ne(α, n)²⁵Mg reaction has also been updated according to Adsley et al. (2021). This reaction is in competition with the ²²Ne(α, γ)²⁶Mg, whose rate has also been updated based on the same reference. These updated rates differ significantly from previous evaluations (Longland et al. 2012; Wiescher et al. 2023), and their reliability is still under debate, with additional experimental measurements expected to clarify the situation (e.g. Santra et al. 2025). For this study, we adopt these new rates to obtain a first estimation of its impact on massive star nucleosynthesis. Rates from Adsley et al. (2021) are approximately twice as high as those from Longland et al. (2012) and are expected to enhance the free neutron density during the C- to O-burning phases.

The second main contributor to free neutrons comes from the ¹³C(α, n)¹⁶O reaction. In D24, we used the rates from the NACRE compilation (Angulo et al. 1999), but we now adopt the updated rates proposed by Ciani et al. (2021). This new reference results in reduced neutron production while preserving a higher ¹³C abundance throughout the He-burning phase. There is broad consensus regarding the reliability of these neutron source reaction rates.

From the original Netgen database and in line with recent references, we have updated four additional linked reactions, as listed in Table 2. Finally, we note that the update of the ¹²C(α, γ)¹⁶O reaction should not alter significantly the s-process nucleosynthesis, as demonstrated by deBoer et al. (2017).

4. Nuclear reaction rates impact: From C-burning to O-burning

4.1. General evolution

In this section, we compare the evolutionary outcomes of the rotating models under each of the three nuclear reaction rate scenarios, CF88, HIN(RES) and TDHF, described in detail in Sect. 3. Figure 2 shows the evolution of the central temperature and density of our 15 M_⊙ models for each of the three nuclear reaction rates references. Models from D24 are also included in the figure for comparison. Compared to the D24 models (computed up to Ne-ignition), the new models exhibit a hotter and denser core. This is a consequence of the updated reaction rates for ¹²C(α, γ)¹⁶​O during the He-burning phase and the revision of the EOS in the advanced phases. As observed in D24, the C-ignition conditions depend on the reaction rates. In the new models, ignition occurs at lower temperatures in the TDHF and at higher temperatures in the HIN(RES), relative to the CF88 model. We note that the overall shift to higher temperatures in the new models, compared to D24 reduces the difference between the CF88 and HIN(RES) models, as the C-burning temperatures of the 15 M_⊙ star are now closer to the peak of the resonance, where the CF88 and HIN(RES) rates converge (see the top right panel in Fig. 1). A similar behaviour can be observed for O-ignition where the two models following the hindrance scenario show a shift towards ignition in a denser and hotter core.

Fig. 2. Evolution of central temperature (TC) vs. central density (ρC) for rotating models of 15 M⊙ at solar metallicity, shown for the three nuclear reaction rate references considered in this work along with the full Hindrance scenario (HIN model) with each model colour-coded accordingly. Nuclear ignition points are marked by crosses. Models from D24 are indicated in dashed-lines for the CF88, HIN, and HIN(RES) models. The dash-dotted grey line indicates the boundary for degenerate conditions inside the core, while the shaded grey area denotes the domain affected by electron-positron pair instability.

Changes in the rates can also be observed from the perspective of the energy output associated with the various reactions during the evolution. In Fig. 3, we show the total energy production profile for a rotating 15 M_⊙ star during the shell C-burning and core O-burning phases, for each of the nuclear rate prescriptions explored in this study. Nuclear energy produced in the core or in shell is compared to the neutrino energy loss, which becomes significant during the core C-burning phase. Moderate differences between models are observed, mainly reflecting the readjustment of the stellar structure which compensates the lower HIN(RES) rates by higher central temperatures and shorter burning lifetimes, resulting in a modest increase in central energy production.

Fig. 3. Energy production rates ϵ [erg.g−1.s−1] inside a 15 M⊙ rotating star for varying reaction rates. Left: Energy during shell C-burning, approximatively 3000 years after core C-ignition. Right: Energy during core O-burning at central X16O = 0.3. The pale areas in orange indicate where the gravitational energy is negative, i.e. corresponding to an expansion.

4.2. Burning lifetimes

The durations of the various nuclear burning phases are especially sensitive to changes in reaction rates. We observe significant differences in the present work compared to the values reported in D24, which we attribute to the updated reference for the ¹²C(α, γ)¹⁶​O reaction. As discussed in Sect. 3.1, the rates from deBoer et al. (2017) increase the ratio of ¹²C/¹⁶O from about 0.45 to about 0.51, in the case of 15 M_⊙ model, as can be seen in Fig. 4. In parallel with this increase, the change in ¹²C(α, γ)¹⁶​O rates also leads to a reduction in the amount of ²⁰Ne produced at He-exhaustion. Overall, the transition from the Kunz et al. (2002) to the deBoer et al. (2017) rate results in shorter He-, Ne-, and O-burning lifetimes, and an extended C-burning phase.

Fig. 4. Mass fraction of central abundances for a rotating 15 M⊙ from He-burning to O-burning, shown for the three nuclear models: CF88, HIN(RES), and TDHF (line-coded). For comparison, a CF88 model using the rates from Kunz et al. (2002) is shown in dotted lines. Shaded areas show the successive nuclear burning phases with colours corresponding to He-, C-, Ne-, and O-burning (grey, red, green, and blue, respectively).

Regarding the impact of fusion reaction rates variations under the three scenarios, Table 3 presents the lifetimes for the C-burning, Ne-burning, and O-burning phases for rotating and non-rotating models across the explored mass-range. To determine these lifetimes we define the ignition point as the moment when the central mass fraction of the main fuelelement decreases by 3 × 10⁻³ from its maximum value X_fuel,max. We define the point of exhaustion as the time when the central mass fraction reaches minimum X_fuel,min set to 10⁻⁵ for C, and 10⁻³ for Ne and O. In other words, the phase lifetime is: $\mathrm{t}({\mathrm{X}}_{\mathrm{fuel},\max }-3\times {10}^{-3})-\mathrm{t}({\mathrm{X}}_{\mathrm{fuel},\min })$. The lower rates from the HIN(RES) result in comparable or slightly shorter C-burning lifetimes (with a numerical uncertainty of ≈100 years). In contrast, the higher reaction rates from the TDHF model lead to a significant increase in C-burning lifetime for all models relative to the CF88 reference. This result may appear counter-intuitive, as one might expect lower rates to yield longer lifetimes due to less efficient burning. However, as observed by M22 and D24, reduced fusion probabilities decreases the energy output in the burning core. To maintain hydrostatic equilibrium, the star then responds by contracting its core, which in turn heats up (see the crosses in Fig. 2), thus boosting the burning and reducing the phase lifetime. Additionally, the lower fusion probability delays ignition to higher densities, further shortening lifetime.

For the Ne-burning phase, differences between models are less significant, and no clear monotonic trend emerges in response to the changes in ¹²C+¹²​C reaction rates, as previously noted by Bennett et al. (2012). Indeed, for the advanced Ne- and O-burning phases, the lifetime is governed by a more complex interplay between the central temperature reached at the end of the C-burning and the residual abundances of Ne and O.

The O-burning lifetimes are notably shorter in the HIN(RES) models, due to both a lower amount of available ¹⁶O and the correspondingly reduced ¹⁶O+¹⁶​O reaction rates. This result reflects the same underlying mechanism discussed for the ¹²C+¹²​C reaction, with ignition conditions shifted to higher temperatures and densities, as shown in Fig. 2. As suggested by Fig. 1, the CF88 and TDHF models yield similar O-burning lifetimes due to their nearly identical rates at O-burning temperatures.

Finally, we note that both rotating and non-rotating models exhibit trends consistent with those reported in D24; we refer the reader there for further details.

4.3. Stellar structure

The chemical structure of the most abundant elements is shown in Fig. 5 at the end of the C-burning and O-burning phases for a 15 M_⊙ star under each of the nuclear reaction rate scenarios considered. At the end of the C-burning phase, all models predict similar CO core sizes of about 2.50 M_⊙. In the centre, chemical abundances are impacted by the change in reaction rates. The most striking difference is observed for ²³Na, whose abundance is three times higher in the TDHF model than in the other two scenarios. This enhancement is attributed to a lower efficiency in the conversion of ²³Na into ²⁰Ne via subsequent (p, α) reactions, due to the lower central temperature. The extent of the carbon shells is also significantly impacted by the choice of reaction rates, as already noted in M22. In Appendix D.1, we show the Kippenhahn diagrams for the 20 M_⊙ model under each scenario. It is clear that the highest rates from TDHF produce the least extended convective C-burning shell while the lowest rates from HIN(RES) result in a more extended shell.

Fig. 5. Profiles of abundances of 4He, 12C, 16O, 20Ne, 23Na, 28Si, 32S, 34S, and 44Ti (colour-coded) for a rotating 15 M⊙, shown for different reaction rates scenarios distinguished by line style. Left panel: End of the C-burning phase. Right panel: End of the O-burning phase.

The right panel of Fig. 5 shows abundance profiles at the end of the O-burning phase. ²⁸Si, ³²S, and ³⁴S are the main products of this phase, formed through secondary reactions involving the primary ²⁸Si and ³¹P products mentioned in Eq. (3), leading notably to a large abundance of ³⁴S, consistent with finding by Farmer et al. (2016). This amount of ³⁴S is lower in the CF88 model compared to the TDHF and HIN(RES) cases. Moreover, the lower HIN(RES) rates result in a reduced central abundance of ²⁸Si relative to CF88 and TDHF. The resulting Si-core mass, M_Si, is smaller for the CF88 and TDHF models (≈1.10 M_⊙) than for HIN(RES) (≈1.20 M_⊙, defined by ${\mathrm{X}}_{\mathrm{Si}}+{\mathrm{X}}_{\mathrm{S}}<0.2$), in agreement with the results of Hirschi et al. (2004) at the end of O-burning. Other elements are also affected, such as ⁴⁴Ti, which is especially less abundant in the TDHF model compared to the other two scenarios, with mass fraction reaching only ∼6 × 10⁻⁹, instead of 4–7 × 10⁻⁸ in the CF88 and HIN(RES) models.

To sum up, ¹²C+¹²​C and ¹⁶O+¹⁶​O rate variations impact the chemical structure, central abundances of several elements, and the Si core mass at the end of the O-burning phase. However, they have a negligible influence on the size of the CO core at this evolutionary phase.

5. Nucleosynthesis

Neutron capture processes are believed to be at the origin of most heavy elements in the Universe (see e.g. the review by Arcones & Thielemann 2023). One of these, the slow neutron-capture process (s-process), develops in massive stars, at neutron densities of the order of ${\mathrm{N}}_{\mathrm{n}}={10}^5-{10}^{10}$ cm⁻³. It occurs mainly during the core He-burning phase (e.g. Langer et al. 1989; Prantzos et al. 1990; Raiteri et al. 1991a) with ²²Ne(α, n)²⁵Mg as the main neutron source. The ¹³C(α, n)¹⁶O reaction produces a short neutron burst at the onset of core He-burning, barely producing heavy elements. These s-process elements are subsequently transported to the outer stellar layers over the relatively long duration of the He-burning phase. During the brief advanced evolutionary stages, the C-burning shell also contributes to the s-process (e.g. Raiteri et al. 1991b; The et al. 2007; Pignatari et al. 2013). Minor contributions are also expected during the core C-burning phase (e.g. Bennett et al. 2012; Pignatari et al. 2013). In this specific phase, rapid transport of nucleosynthetic products from the core to the outer shells should be involved to enable this process to contribute to the chemical galactic enrichment. During the core O-burning phase, trans-iron elements are destroyed via photodisintegration reactions. In addition, rotation is is known to boost the s-process during core He-burning by promoting the synthesis of primary ¹³C and ²²Ne (Pignatari et al. 2008; Frischknecht et al. 2016; Choplin et al. 2018; Limongi & Chieffi 2018; Banerjee et al. 2019). The impact of variations in nuclear reaction rates and rotation parameters manifests mainly through changes in central temperatures throughout stellar evolution and modifications in the burning lifetimes, both of which influence the neutron exposure. Indeed, the two main neutron sources for s-process in massive stars–¹³C and ²²Ne–contribute in different proportions depending on the temperature and evolutionary stage. In this context, we have also updated seven important reactions for s-process nucleosynthesis, as listed in Table 2 and introduced in Sect. 3.

5.1. One-zone code for nucleosynthesis

We explore the implications for nucleosynthesis using the one-zone code described in Choplin et al. (2016), and previously applied in Taggart et al. (2019), Williams et al. (2022), Frost-Schenk et al. (2022), M22 and D24. The nuclear network has been extended to track 1454 isotopes (against 737 in the previous version in D24) in order to capture photodisintegration reactions during the core O-burning phase. In this approach, only a single thermodynamic trajectory in temperature and density is followed (either core or shell), and transport processes, such as convection, are neglected. This simplified method provides reliable predictions for the s-process nucleosynthesis during the core He-burning phase (e.g. Frost-Schenk et al. 2022). During these advanced evolutionary phases, rotational transport is still active but exerts a weaker mixing influence due to the short timescales in late stellar evolution (as shown, e.g. in Fig. 2 by Choplin et al. 2017). Thus, this approach is sufficient sufficient to capture and interpret the relative effects introduced by variations in nuclear reaction rates. Numerous studies have investigated s-process nucleosynthesis in the core and the shell regions during the advanced phases of evolution (Arcoragi et al. 1991; Gasques et al. 2007; The et al. 2007; Nishimura et al. 2017).

5.2. s-process nucleosynthesis

We now explore the isotopic abundance predictions obtained from the one-zone code at the end of core O-burning and shell C-burning in our 15 M_⊙ rotating models, where the slightly lower temperatures reached are expected to enhance the sensitivity to differences among the nuclear reaction rates prescriptions (see Fig. 1 or D24). We first examine the chemical evolution of the core, following the same approach as in M22 and D24, by tracking the central density and temperature throughout the successive core C-, Ne- and O-burning phases. Figure 6 shows theproduction of both direct C-burning products, such as Ne, Na, and Mg, and trans-iron elements during the core C-burning phase. After Ne-burning, the subsequent O-burning phase leads to the destruction of most trans-iron elements, and the synthesis of typical O-burning products, including Si, S, Ar, Ca, and Ti, as also observed in the right panel of Fig. 5. At the end of the core O-burning phase, heavy elements exhibit either low (30 < Z < 40) or almost null (Z > 40) abundances (coloured lines in Fig. 6). This is due to photo-dissociation reactions such as (γ, n), (γ, p), and (γ, α), which become increasingly efficient at T ≳ 1.8 GK (e.g. Arnould & Goriely 2003), typically reached during the O-burning phase. The sawtooth structure reflects odd-even effects in the nuclear binding energies due to pairing, and is clearly visible in the abundance pattern at the end of O-burning. Beyond the low abundance of fluorine (Z = 9), the largest deviations – up to ∼1 dex – relative to the CF88 reference are obtained for the HIN(RES). These feature a modest overproduction of light elements between Ca (Z = 20) and Fe (Z = 26), and an underproduction of s-process elements at Z > 30, such as Ge, Se, or Sr. We mainly observe the consequences of the lower rates associated with the Hindrance scenario for the ¹⁶O+¹⁶​O reaction, which leads to a shorter burning phase and reduced neutron exposure. This mirrors the behaviour seen in D24 for the ¹²C+¹²​C reaction during core C-burning phase. It is worth noting that the cut-off for Z > 40 is unaffected by the choice of reaction rates.

Fig. 6. Mass fractions at the end of the core O-burning from one-zone calculations based on central conditions of a 15 M⊙ star with and without rotation, shown for three nuclear rate references (colour-coded). The dashed grey line shows the abundances at the end of the core C-burning from the one-zone calculation of the non-rotating CF88 model, while the solid grey line indicates the initial composition at the onset of core C-burning. Lower panels display the abundances normalised to those of the non-rotating CF88 case.

In all models, the s-process distribution in the core is further modified during carbon-burning shell episodes. To explore the associated nucleosynthesis, we have extracted the evolution of density and temperature from the first two C-burning shells, which occur between the end of core C-burning and onset of core Ne-ignition (see technical details in Appendix C). Compared to the helium core, the carbon shells overproduce specific isotopes due to the higher neutron densities, reaching up to 10¹⁰ cm⁻³ within the second C-burning shell, as shown in Table 4. For instance, in the Sm, Eu, and Gd region, we observe an overproduction of ^152, 154Gd, corresponding to the opening of branching points at A = 151 and A = 154. This enhancement is accompanied by a reduction in the abundance of intermediate isotopes along the same nucleosynthetic path (e.g. ^151, 153Eu is depleted by an order of magnitude). However, differences between models remain modest. For example, the ¹⁵²Gd isotope is overproduced by a factor of 10 in the HIN(RES) scenario, compared to a factor of 2 in the TDHF scenario, relative to its abundance at the end of the He-burning phase. The reasons for the small changes observed are twofold. First, the corresponding temperature and density trajectories do not differ significantly enough across the three explored cases to drive substantial nucleosynthetic differences (see Appendix C.1). Second, and more generally, for all our 15 M_⊙ rotating models, the s-process is not very efficient, as quantitatively demonstrated by the peak neutron densities and the total neutron exposures given in Table 4, with maximum exposures reaching only ∼0.1 mbarn⁻¹ during the first C-burning shell.

5.3. Variation of the ¹²C+¹⁶O reaction

Finally, this paper is also intended to guide further experimental determinations of reaction cross-sections. To this end, we performed a sensitivity study for the ¹²C+¹⁶​O reaction. Unlike the main fusion reactions (1 and 3), this reaction does not contribute significantly to energy generation. As such, it is only relevant primarily at the end of core C-burning or during shell C-burning, once the required thermodynamic conditions, low ¹²C abundance, high ¹⁶O abundance, and suitable temperature are met (Iliadis et al. 2010). We investigate nucleosynthesis in the first C-burning shell of our 15 M_⊙ model, where the differences among reaction rates are most pronounced. We conduct the sensitivity study by arbitrarily multiplying or dividing the CF88 reaction rates for ¹²C+¹⁶​O by constant factors. We ind that only an extreme rate variation (by a factor ∼10⁴ would significantly affect the nucleosynthesis.

Specifically, increasing the rate by this factor results in an overproduction of Al and Si by factors of 2.40 and 1.25, respectively, while decreasing the rate yields underproduction of Mg, Al, and Si by ∼20. We conclude that the differences in ¹²C+¹⁶​O reaction rates among the CF88, HIN(RES), and TDHF scenarios (see Fig. 1) are too small to significantly affect nucleosynthesis. Given the current measurements of nuclear reaction rates and predictions for this reaction, no strong impact is expected from improved accuracy, unless previously unidentified resonance enhancing the rate by several orders of magnitude is discovered at astrophysical energies.

6. Remnant predictions

A major goal of massive star modelling is to determine which stars explode and which collapse into black holes (e.g. Chieffi & Limongi 2020; Chieffi et al. 2021). It is essential for constraining black hole (BH), neutron star (NS), and supernova (SN) populations synthesis models. Recent works, such as Ugolini et al. (2025) investigate the link between progenitor mass and the remnant mass produced by core collapse supernovae.

Several attempts to estimate the type of remnant (and the associated explosion type) resulting from massive star evolution have been proposed in the literature (e.g. Maltsev et al. 2025). Some approaches require knowledge of quantities near the onset of core collapse (e.g. Ugliano et al. 2012; Ertl et al. 2016; Müller et al. 2016), while others infer the remnant type based on the CO core properties (e.g. Patton & Sukhbold 2020; Maltsev et al. 2025). Since our models are evolved up to the end of the O-burning phase, we estimate their fate based on properties at this state. The work by Patton & Sukhbold (2020) sought to systematically predict whether the remnant would be a BH or NS, and to estimate its mass, using only the CO core mass and central carbon mass fraction at He-exhaustion as diagnostic parameters in non-rotating stars. From these two parameters, they build a remnant classification map of explosion predictions. While such a two-parameter model cannot fully capture the complex dynamics of core collapse, and state-of-the-art multi-dimensional simulations would be required for a detailed outcome (see e.g. Wang et al. 2022; Burrows et al. 2024), it is still possible to use these parameters to explore how variations in the input physics between models could affect the final fate of massive stars.

6.1. Updated ¹²C(α,γ)¹⁶O reaction rates

In Sect. 3 we have shown that updating the reaction rates of ¹²C(α, γ)¹⁶O from those of Kunz et al. (2002) to those of deBoer et al. (2017) changes the amount of ¹²C left in the core at He-exhaustion, while leaving the CO core mass nearly unchanged (see Appendix E for the exact values of M_CO and ${\mathrm{X}}_{{}^{12}\mathrm{C}}$). Using the remnant classification of Patton & Sukhbold (2020) based on the Kunz et al. (2002) rates, we find that all models would lead to NS formation, with the exception of the 20 M_⊙ model, which collapses into a BH. In Fig. 7, we overlay the position $({\mathrm{M}}_{\mathrm{CO}},{\mathrm{X}}_{{}^{12}\mathrm{C}})$ of our models onto the remnant-type map from Patton & Sukhbold (2020). The general trend resulting from adopting the deBoer et al. (2017) rates is clear: they yield a higher central carbon mass fraction ${\mathrm{X}}_{{}^{12}\mathrm{C}}$ and a slightly smaller CO core mass. Additionally, as the initial mass increases, M_CO increases while ${\mathrm{X}}_{{}^{12}\mathrm{C}}$ decreases, owing to He-burning occurring at higher temperatures. From the remnant map of Patton & Sukhbold (2020), we observe that, broadly speaking, a higher M_CO or a lower ${\mathrm{X}}_{{}^{12}\mathrm{C}}$ correlates with an increased likelihood of BH formation. However, the map contains multiple islands of explosiveness as seen in the background of Fig. 7. The variation in the central value of ${\mathrm{X}}_{{}^{12}\mathrm{C}}$ resulting from the updated reaction rates is large enough to change the predicted remnant type according to Patton & Sukhbold (2020). In contrast, using the prescription of Maltsev et al. (2025), all our models would be classified as NS progenitors, since all their CO core masses are below 6.6 M_⊙, the threshold below which they predict NS formation for all models.

Fig. 7. Position of our models in the $({\mathrm{M}}_{\mathrm{CO}},{\mathrm{X}}_{{}^{12}\mathrm{C}})$ parameter space for all initial masses using the rates from Kunz et al. (2002, red) and deBoer et al. (2017, green). The background colours indicate the remnant mass predictions according to Patton & Sukhbold (2020), where dark grey regions correspond to black hole formation.

6.2. Fusion reactions ¹²C+¹²C, ¹²C+¹⁶O, and ¹⁶O+¹⁶O

After the formation of the CO core, the fusion reactions of ¹²C and ¹⁶O will continue to impact its internal structure. The different rates tested for these two fusion reactions lead to variations in the convective core masses during C-, Ne-, and O-burning, as well as in their respective lifetimes (see Table 3). We also note variations in the convective shell burning, both in the shell mass and in the duration of the burning episode. From the Kippenhahn diagrams of the 20 M_⊙ non-rotating models (Appendix D), we observe the resulting differences in shell C-burning, particularly in the extent of the second carbon shell. In the HIN(RES) model, the convective shell extends almost 0.5 M_⊙ further out compared to TDHF, thereby altering the compactness evolution in this layer (see Appendix D).

The contraction of the models during the second carbon shell episode is faster than during the first. The smaller C-burning shell of the THDF framework displays a comparatively shallow contraction relative to, e.g. the HIN(RES) shell (noticeable by the smaller slope of the iso-radii lines in Fig. D.1). The contraction rate differs significantly among models during the Ne-burning phase. A longer Ne-burning phase sustains the contraction of the outer shells longer. The TDHF calculations yield the shortest Ne-burning phase, resulting in a correspondingly brief contraction period of the outer layers.

The results of these differences in core and shell burning shape the compactness profile and the entropy distribution within the CO core at late times. In Fig. 8, we show the distribution of entropy and compactness as a function of Lagrangian mass near the end of core O-burning for the same models as in Appendix D. We note that the peak compactness, reached at about 1.4 M_⊙, is similar across all models. More importantly, the TDHF model shows the lowest compactness between 2 ≲ M_r ≲ 3.5 M_⊙. This corresponds precisely to the region affected by the second C-burning shell, whose extent and duration vary significantly among the TDHF, CF88, and HIN(RES) scenarios. While compactness is commonly used as a proxy for explosiveness (O’Connor & Ott 2011, but see Aloy & Obergaulinger 2021; Wang et al. 2022; Burrows et al. 2024), in this study we have assessed the explosion outcome primarily through the CO core mass and the central ¹²C abundance, finding that all three models considered here result in successful explosions and form NSs. Although these models have not evolved all the way to collapse, we already observe indications that changes in nuclear reaction rates impact the density and compactness structure, with potential consequences for the explosion mechanism. In Fig. 8, we also represent the entropy per baryon S to illustrate the structural differences among the scenarios. We highlight the variation in shape and location of the entropy discontinuities. The first entropy jump in the HIN(RES) model occurs furthest out due to a larger O-burning core. Additionally, the second entropy jump is smoother in HIN(RES) and sharper in CF88. This reflects the fact that the CF88 model retains an active C-burning shell during O-burning, whereas in the HIN(RES) case, the last convective shell vanishes nearly a year earlier, allowing for diffusive smoothing and contraction in that region. The TDHF cross-sections do not induce an active C-burning shell during the O-burning either. Nevertheless, only a few weeks elapse between the shutdown of O-burning and the time at which the compactness and entropy are evaluated in Fig. 8. Therefore, the entropy discontinuity does not have sufficient time to diffuse, unlike in the HIN(RES) case. The location, shape, and magnitude of these entropy jumps are another critical factor in determining the explosion outcome of massive stars (Imasheva et al. 2025), providing further evidence of structural variability in the CO core driven by the changes in fusion reaction rates explored in this study.

Fig. 8. Compactness (full line), defined in Eq. (D.1), and entropy (dashed line) for the 20 M⊙ models using the rates from CF88, TDHF, and HIN(RES). Models are shown at the midway through core O-burning.

7. Summary and outlook

We have computed models of massive stars to explore the sensitivity of four key nuclear reactions from He-burning to O-burning phases. In agreement with previous studies, we confirm that the reactions ¹²C(α, γ)¹⁶​O and ¹²C+¹²​C are among the most influential in shaping the stellar structure during the advanced burning stages. The former reaction, active during He-burning, governs the ¹²C/¹⁶O ratio in the core at He-exhaustion and thus sets the amount of fuel available for the subsequent C-burning phase, which is governed by the latter. C-burning is the longest of advanced burning stages, during which the energy produced in the core is mainly evacuated by neutrinos. A longer C-burning phase implies greater energy and entropy losses, thereby increasing the core’s susceptibility to degeneracy effects. We confirm that lower ¹²C(α, γ)¹⁶​O rates result in more residual carbon at He-exhaustion and to the formation of less massive cores, delaying the onset of core collapse. In our models, the reduced rates from deBoer et al. (2017) result in shorter He-burning phases, longer C-burning phases, high ¹²C/¹⁶O ratios, and slightly smaller CO core masses (Fig. 7). Consequently, the final fate of the star is also impacted by the adopted rates. Using the remnant prediction scheme by Patton & Sukhbold (2020) (developed for non-rotating models), we demonstrate that the differences in nuclear reaction rates can alter the remnant mass and type of non-rotating models. For example, our 20 M_⊙ model is predicted to form a BH when using the rates from Kunz et al. (2002), but an NS when using those from deBoer et al. (2017).

Variations in the ¹²C+¹²​C and ¹⁶O+¹⁶​O fusion rates primarily affect the lifetimes and ignition conditions of C- and O-burning phases, which in turn influence the cores susceptibility to degeneracy effects. Under the HIN(RES) scenario, these rates lead to shorter lifetimes and ignition at higher temperatures and densities, whereas TDHF rates yield longer C-burning lifetimes and ignition at lower temperatures and densities, albeit to a lesser extent. These differences alter both the abundances of the major isotopes and the internal structure. In particular, the extent of the convective C-burning shells responds to the adopted rates: they expand under lower fusion probabilities and contract under higher ones, thereby modifying the mass distribution and the compactness profile. Although the updated ¹²C+¹²​C and ¹⁶O+¹⁶​O rates induce relatively modest differences in the global nucleosynthetic yields, they produce significant changes in the chemical structure during advanced combustion stages. These structural differences are expected to impact the final ejecta composition, motivating future dedicated studies. Regarding the ¹²C+¹⁶​O reaction, we find that only changes in the rates by several orders of magnitude (≈10⁴) would noticeably alter the nucleosynthesis during the shell C-burning phase, where this reaction should operate. The two new rate models explored in this work do not reach such deviations and are therefore not expected to significantly influence the associatednucleosynthesis.

On the one hand, experimental measurements of fusion reactions at astrophysical energies are challenged by the extremely low cross-sections involved. On the other hand, theoretical approaches, such as the TDHF method used in this work, can estimate reaction rates at these low energies, but with the limitation that they do not account for resonances in the compound system or its nuclear level density, particularly relevant for ¹²C+¹²​C fusion. These two approaches are thus complementary and both are essential for improving stellar evolution models. Progress in both experimental methods and theoretical frameworks will be critical to achieve precise and unified reaction rates inputs for stellar modelling. In this context, a recent advancement in TDHF modelling has demonstrated its ability to reproduce certain resonant structures of the ¹²C+¹²​C system, offering promising new insights on them (Close et al. 2025). Finally, it is important to stress that the nuclear reaction rates investigated here represent only one component of a much broader and complex picture. Our results also depend sensitively on the treatment of mixing processes. Uncertainties in angular momentum transport, convective boundary mixing, or mass loss through stellar winds can all influence the outcomes, in some cases more significantly than variations in nuclear rates. Particularly, the prescriptions adopted for stellar rotation and the possible role ofmagnetic fields exert cumulative effects throughout the stellar evolution, ultimately shaping its final fate (e.g. Meynet et al. 2013; Dumont et al. 2021; Nandal et al. 2024). These aspects will be the subject of a forthcoming study (Dumont et al., in prep.).

Acknowledgments

This work was supported by the European Union (ChETEC-INFRA, project no. 101008324). We thank the anonymous referee for helpful comments on the manuscript. T.D. thanks S. Ekström and R. Hirschi for helpful discussions. MAA and AG acknowledge the support through the grants PID2021-127495NB-I00 and PID2021-125485NB-C21 funded by MCIN/AEI/10.13039/501100011033 and by the European Union, and the Astrophysics and High Energy Physics programme of the Generalitat Valenciana ASFAE/2022/026 funded by MCIN and the European Union NextGenerationEU (PRTR-C17.I1) as well as support from the Prometeo excellence programme grant CIPROM/2022/13 funded by the Generalitat Valenciana. A.C. is a Postdoctoral Researcher of the Fonds de la Recherche Scientifique – FNRS. K.G. acknowledges support from DE-SC0023175 (Office of Science, NUCLEI SciDAC-5 collaboration). This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: Nuclear network

The full set of isotopes included in the GenValNet48 network of the GENEC code is listed in Tab. A.1.

Appendix B: Nuclear reaction rates

Appendix B.1: ¹²C+¹²​C theoretical predictions

As discussed in Sect. 3.2.2, several energy density functionals have been used in the TDHF calculations for the ¹²C+¹²​C nuclear reaction. Figure B.1 presents the resulting S-factor predictions obtained using three of these functionals. In the astrophysical energy range of interest, all TDHF predictions lie above those from the other models compared in this work. Since no specific functional is favoured by the authors, we adopt the results from the Skyrme SLy4d functional (shown in red), which gives the most pronounced deviations, for our sensitivity analysis.

Fig. B.1. S-factor model predictions for the 12C+12​C fusion reaction as a function of the relative energy Erel, obtained using different functionals in the TDHF calculations. The SLy4d functional (in red) is the one considered in this work. For comparison, the CF88 reference is indicated by the black dotted line.

Appendix B.2: Neutron source reactions

The two primary nuclear reactions driving the slow neutron-capture process are ¹³C(α, n)¹⁶O and ²²Ne(α, n)²⁵Mg. In rotating stellar models, the ¹³C(α, n)¹⁶O reaction is favoured due to the higher temperatures reached in the core and in the burning shells. In the present work, we have updated the rates for both reactions. Figure B.2 gives the ¹³C(α, n)¹⁶O reaction rates adopted in this study, and based on the evaluation by Ciani et al. (2021), in comparison to the NACRE compilation used in D24. During the He- and C-burning phases, the updated rates are slightly lower than for the NACRE reference.

Fig. B.2. Reaction rates for 13C(α, n)16O from Ciani et al. (2021), normalised to the NACRE reference rates (Angulo et al. 1999). Dark grey shade areas indicate the temperature ranges corresponding to He- and C-burning in a rotating 17 M⊙ GENEC model with an initial rotation rate of Vini/Vcrit = 0.4.

Figure B.3 shows the ²²Ne(α, n)²⁵Mg reaction rates from Adsley et al. (2021) used in the present work, compared to the reference rates from (Iliadis et al. 2010) used in D24. During the He- and C-burning phases, we expect lower and higher reaction rates, respectively. The same figure also displays the reaction rates for the competing ²²Ne(α, γ)²⁶Mg reaction, which acts in opposition to the neutron-producing channel.

Fig. B.3. Reaction rates from Adsley et al. (2021), normalised to the reference rates of (Iliadis et al. 2010), for the 22Ne(α, n)25Mg and 22Ne(α, γ)26Mg nuclear reactions. Dark grey shade areas indicate the temperature ranges corresponding to He- and C-burning in a rotating 17 M⊙ GENEC model with an initial rotation rate of Vini/Vcrit = 0.4.

Appendix C: Nucleosynthesis one-zone code: Shell C-burning trajectories

To select a representative trajectory for the temperature and density conditions during carbon shell burning, we identify, at each step of the evolution, the location where the energy production from the ¹²C+¹²C reaction peaks. This location consistently corresponds to the base of the convective shell burning zone, as illustrated, e.g. in Fig. 3. In most of the models, two successive convective C-burning shells develop. The evolution of temperature and density within these shells, plotted as a function of the ¹²C mass fraction is shown in Fig. C.1 for a rotating 15 M_⊙ star and the three nuclear reaction rates references: CF88, HINRES, and TDHF. The second shell exhibits slightly higher temperatures, while the densities remains comparable across the two shells. It is worth noting that the second shell has a shorter lifetime (about 1/8 shorter), resulting in a correspondingly reduced neutron exposure.

Fig. C.1. Temperature and density conditions within the first (solid lines) and second (dashed lines) carbon burning shells for the 15 M⊙ rotating star, shown as a function of the 12C mass fraction.

Appendix D: Kippenhahn diagrams

Kippenhan diagrams are a useful tool to understand the internal evolution of stellar models under different nuclear reaction rate scenarios. Figure D.1 shows the Kippenhahn diagrams for the 20 M_⊙ non-rotating models explored in this study. Each panel illustrates the convective zones during the core C-burning phase, the successive shell C-burning episodes, core and shell Ne-burning, and core O-burning. A significant effect of the nuclear rates variations is observed in the extent of the C-burning shells. Especially, the second C-burning shell is significantly more compact in the TDHF model compared to the CF88 and HIN(RES) cases.

In the Kippenhahn diagrams, the background is colour-coded to reflect the compactness of each mass shell, defined as

$$\begin{array}{c}\hfill \xi =\frac{M_r/M_{\odot }}{r/1000\text{km}}.\end{array}$$(D.1)

An overall trend of increasing compactness is evident by the darker green shades, especially towards the end of the computed evolution and in the vicinity of the convective O-burning core. To illustrate the contraction of the selected mass shells, we also overlay iso-radius contours in Fig. D.1. The contraction manifests as increasing mass accumulated at each considered radius (hence, iso-radius contours grow secularly).

Fig. D.1. Kippenhahn diagrams of the 20 M⊙ non-rotating stars from the core C-burning to the core O-burning for the three nuclear reaction rate scenarios. Dark blue shade shows the convective zones. Iso-radii contours are shown by dotted red lines. The compactness ξ is indicated by a colour bar. Top: CF88. Middle: HIN(RES). Bottom: TDHF.

Appendix E: Remnant properties

Table E.1 presents the remnant properties predicted by our non-rotating models. The CO core mass is determined at the end of the He-burning phase and is defined as the mass coordinate where the central helium mass fraction drops below X_He < 10⁻¹. In cases where the model predicts black hole formation, the remnant mass corresponds to that of the He core, defined as the mass coordinate where the central hydrogen mass fraction falls below X_H < 10⁻¹. The remnant type is inferred based on the values of ${\mathrm{X}}_{{}^{12}\mathrm{C}}$ and M_CO following the criteria of Patton & Sukhbold (2020). For models resulting in NS formation, the remnant mass M_rem is taken from the tables of Patton & Sukhbold (2020), which define the NS mass as M₄, corresponding to the mass coordinate where the entropy per baryon reaches 4k_B/baryon. The last column displays M₄ for models that have reached the end of O-burning; this value closely approximates the final mass expected at the pre-SN stage. For models that have not reached the end of O-burning, this quantity is not available and is indicated by ”/”.

All Tables

All Figures

Fig. 1. Reaction rates for the four main nuclear reactions involving 12C and 16O. Top left: Rates from deBoer et al. (2017) normalised to Kunz et al. (2002) reference for the 12C(α, γ)16O reaction. The grey hatched areas indicate the temperature ranges where core and shell He-burning occur in the GENEC rotating model explored mass range. For the rest of the panels, nuclear reaction rates are taken from CF88 (black), HIN(RES (blue), and TDHF (magenta), all normalised to the CF88 rates (NA < σv> CF88). Shaded areas (with corresponding colour codes) indicate the temperature ranges where core and shell C-burning and core O-burning take place in the 15 and 17 M⊙ GENEC rotating models.

In the text

Fig. 2. Evolution of central temperature (TC) vs. central density (ρC) for rotating models of 15 M⊙ at solar metallicity, shown for the three nuclear reaction rate references considered in this work along with the full Hindrance scenario (HIN model) with each model colour-coded accordingly. Nuclear ignition points are marked by crosses. Models from D24 are indicated in dashed-lines for the CF88, HIN, and HIN(RES) models. The dash-dotted grey line indicates the boundary for degenerate conditions inside the core, while the shaded grey area denotes the domain affected by electron-positron pair instability.

In the text

	Fig. 3. Energy production rates ϵ [erg.g−1.s−1] inside a 15 M⊙ rotating star for varying reaction rates. Left: Energy during shell C-burning, approximatively 3000 years after core C-ignition. Right: Energy during core O-burning at central X16O = 0.3. The pale areas in orange indicate where the gravitational energy is negative, i.e. corresponding to an expansion.
In the text

Fig. 4. Mass fraction of central abundances for a rotating 15 M⊙ from He-burning to O-burning, shown for the three nuclear models: CF88, HIN(RES), and TDHF (line-coded). For comparison, a CF88 model using the rates from Kunz et al. (2002) is shown in dotted lines. Shaded areas show the successive nuclear burning phases with colours corresponding to He-, C-, Ne-, and O-burning (grey, red, green, and blue, respectively).

In the text

	Fig. 5. Profiles of abundances of 4He, 12C, 16O, 20Ne, 23Na, 28Si, 32S, 34S, and 44Ti (colour-coded) for a rotating 15 M⊙, shown for different reaction rates scenarios distinguished by line style. Left panel: End of the C-burning phase. Right panel: End of the O-burning phase.
In the text

Fig. 6. Mass fractions at the end of the core O-burning from one-zone calculations based on central conditions of a 15 M⊙ star with and without rotation, shown for three nuclear rate references (colour-coded). The dashed grey line shows the abundances at the end of the core C-burning from the one-zone calculation of the non-rotating CF88 model, while the solid grey line indicates the initial composition at the onset of core C-burning. Lower panels display the abundances normalised to those of the non-rotating CF88 case.

In the text

	Fig. 7. Position of our models in the $({\mathrm{M}}_{\mathrm{CO}},{\mathrm{X}}_{{}^{12}\mathrm{C}})$ parameter space for all initial masses using the rates from Kunz et al. (2002, red) and deBoer et al. (2017, green). The background colours indicate the remnant mass predictions according to Patton & Sukhbold (2020), where dark grey regions correspond to black hole formation.
In the text

	Fig. 8. Compactness (full line), defined in Eq. (D.1), and entropy (dashed line) for the 20 M⊙ models using the rates from CF88, TDHF, and HIN(RES). Models are shown at the midway through core O-burning.
In the text

	Fig. B.1. S-factor model predictions for the 12C+12​C fusion reaction as a function of the relative energy Erel, obtained using different functionals in the TDHF calculations. The SLy4d functional (in red) is the one considered in this work. For comparison, the CF88 reference is indicated by the black dotted line.
In the text

	Fig. B.2. Reaction rates for 13C(α, n)16O from Ciani et al. (2021), normalised to the NACRE reference rates (Angulo et al. 1999). Dark grey shade areas indicate the temperature ranges corresponding to He- and C-burning in a rotating 17 M⊙ GENEC model with an initial rotation rate of Vini/Vcrit = 0.4.
In the text

	Fig. B.3. Reaction rates from Adsley et al. (2021), normalised to the reference rates of (Iliadis et al. 2010), for the 22Ne(α, n)25Mg and 22Ne(α, γ)26Mg nuclear reactions. Dark grey shade areas indicate the temperature ranges corresponding to He- and C-burning in a rotating 17 M⊙ GENEC model with an initial rotation rate of Vini/Vcrit = 0.4.
In the text

	Fig. C.1. Temperature and density conditions within the first (solid lines) and second (dashed lines) carbon burning shells for the 15 M⊙ rotating star, shown as a function of the 12C mass fraction.
In the text

	Fig. D.1. Kippenhahn diagrams of the 20 M⊙ non-rotating stars from the core C-burning to the core O-burning for the three nuclear reaction rate scenarios. Dark blue shade shows the convective zones. Iso-radii contours are shown by dotted red lines. The compactness ξ is indicated by a colour bar. Top: CF88. Middle: HIN(RES). Bottom: TDHF.
In the text