2. Changing Nuclear Reaction Rates

In this section, we will graduate from using our Intro_MESA_model model directory, and instead use the provided work directory. Changing_rates MESA work directory.

Science goal

In this lab, we will focus on the evolution of a 20 M$_{\odot}$ stellar model from the zero age main-sequence through to core-Helium depletion. We will explore which nuclear reaction rates MESA uses, how to change the rates MESA reads in, and finally we explore the impact of changing the $^{12}$C$(\alpha,\gamma)^{16}$O helium burning reaction rate, and its effect on the core properties of stellar model at core-Helium depletion. These properties are extremely important as the Carbon/Oxygen core mass, $M_{CO}$, is often used in population synthesis calculations as an input into expressions for determining if the stellar model will form a neutron star to a black hole, and the corresponding remenant mass after the model undergoes core-collapse. See Fryer et al. 2012,Ilya et al. 2020, Fryer et al. 2022.

In this lab we only evolve our stellar models to core-Helium depletion, hence we are $M_{CO}$ at core-Helium depletion as a proxy for $M_{CO}$ at core-collapse. These two are not identical, and one must evolve their stellar model to core-collapse if they intend to generate a realistic presupernova stellar structures which can be used to determine the model’s final fate. Note that, most population synthesis prescriptions for stellar model fates are based on stellar models run to core collapse with small nuclear reaction networks, which could problematic as the structure of stellar models computed with small approximate networks are not identical those computed with larger, more physically accurate networks.

Evolving to Core-Helium Depletion

To begin, please make sure to download a copy of the desired changing_rates MESA work directory. This work directory is a slightly modified version of the $MESA_DIR/star/test_suite/20M_pre_ms_to_cc test_suite.

Once downloaded, you can decompress the file by

unzip changing_rates.zip

tree ./changing_rates should return the following.

├── clean
├── history_columns.list
├── profile_columns.list
├── inlist
├── inlist_common
├── inlist_mass_Z_wind_rotation
├── inlist_pgstar
├── inlist_project
├── make
│   └── makefile
├── mk
├── rate_tables
│   ├── c12ag_deboer_sigma_0p0_2000_Tgrid.dat
│   ...
├── re
├── re_nomodfiles
├── README.rst
├── rn_nomodfiles
├── run_lab
├── src
│   ├── run_star_extras.f90
│   └── run.f90
└── zams.mod

All relevent files are briefly described in the table below

MESA STAR work directory

Filename	Description
clean	A bash file for cleaning the model directory.
inlist	The header inlist which points to all other inlists to determine which inlists are read and in what order.
inlist_mass_Z_wind_rotation	A supplemental inlist where the mass, metallicity, and mass loss are set.
inlist_project	The main inlist which contains controls for the stellar evolution of the model.
inlist_common	The common inlist which contains most of the stellar model physics.
inlist_pgstar	The inlist which controls the pgstar output for each single star.
make/	A directory containing the makefile.
mk	A bash file for compiling MESA Star in the model directory.
history_columns.list	A log file which determines which history values are saved in data files as a function of model timestep.
profile_columns.list	A log file which determines which profiles values are saved in data files as a function of Mass/radius.
re	A bash file for restarting the star model executable from photos
run_lab	A bash file for running the star model executable.
src/	A directory containing the three files listed below.
run_star_extras.f90	A fortran file which can be modified to agument the stellar evolution routines.

inlist_project, inlist_common, and inlist_mass_Z_wind_rotation are the three main files that contain the microphysics information of our stellar evolution simulation.

Helium Burning Nuclear Reactions

Nuclear Reaction rate uncertainties have a non-negliggable impact on stellar evolution models and their nucleosynthetic products. See Iliadis et al. 2011,Laird, Fields et al. 2018,Farmer et al. 2019 , Farag et al. 2022, Pignatari et al. 2023

For massive star evolution, four of the most uncertain rates are illustrated below (Figure from Fields et al. 2018)

Some of the most uncertain nuclear reaction rates in nuclear astrophysics are the rates that occur during Helium burning Deboer et al. 2017

That is the triple-$\alpha$ nuclear reaction rates converting Helium into Carbon competes with the c12ag nuclear reaction rate to set the final Carbon and Oxygen abundance in the stellar core at the end of Helium burning

The amount of Carbon and Oxygen in the stellar core is important as it fundamentally changes the nature of Carbon and Oxygen burning, and can have far reaching implications for the resulting presupernova stellar structure and the resulting neutron star, black hole initial mass function.

The uncertainties in the $^{12}$C$(\alpha,\gamma)^{16}$O reaction rate are among the largest (visualized below)

Changing the $^{12}$C$(\alpha,\gamma)^{16}$O rate

MESA’s nuclear reaction rates are a combination of NACRE Angulo et al. 1999 and JINA REACLIB Cyburt et al. 2010. Reaction rate screening corrections are from Chugunov et al. 2007, a dynamic screening method which includes a physical parametrization for the intermediate screening regime and reduces to the weak DeWitt et al. 1973, Graboske et al. 1973 and strong Alastuey and Jancovici 1978, Itoh et al. 1979 screening limits at small and large values of the plasma coupling parameter. Weak reaction rates are based, in order of precedence, on Langanke and Martinez-Pinedo 2000, Oda et al. 1994, and Fuller et al. 1985.

For hardcoded nuclear reaction rates, the definitions are set and called inside $MESA_DIR/rates/private/raw_rates.f90, with definitions for each rate contained inside $MESA_DIR/rates/private/ratelib.f90, or drawn from JINA Reaclib / weaklib.

We would like to change our stellar model to adopt one of the $^{12}$C$(\alpha,\gamma)^{16}$O reaction rates provided by Deboer et al. 2017. These high resolution nuclear reaction rates are available inside ./rate_tables or $MESA_DIR/data/rates_data/rates_tables directories. By varying this rate, we can explore the temperature dependant uncertainty in this nuclear reaction.

TASK
Pick a value of $\sigma$ for the $^{12}$C$(\alpha,\gamma)^{16}$O rate shown in the figure above. Use a different value than those sitting next you.
Change the $^{12}$C$(\alpha,\gamma)^{16}$O reaction rate to one of the Deboer et al. 2017 rates shown in the figure above.
Run your model to completion, and report your values in the google spreadsheet document here

HINT
Look inside your local changing_rates/rate_tables directory, which was copied from $MESA_DIR/rates/rates_data/rate_tables.
To change the rate you are reading you’ll have to modify the rate_list.txt file

❓ Questions
Below are some questions to think about using the pgstar movie output from you stellar model.
1. During core-Helium burning, what happens to $^{14}$N leftover from core-H burning?
2. What do you think should happen to $^{14}$N?
3. How does the your chosen c12ag rate effect the final C/O in the core at core-Helium depletion?
4. Are we missing any reactions with our simplified approx21 network?

Answers: Changing the $^{12}$C$(\alpha,\gamma)^{16}$O reaction rate

uncomment the rate file you would like to read into MESA inside change_rates/rate_tables/rate_list.txt

! c12ag rates from debeor et al. 2017 improved with high resolution (mehta et al. 2022) 
! 0 sigma (median c12ag rate)
r_c12_ag_o16   'c12ag_deboer_sigma_0p0_2000_Tgrid.dat'

! positive sigmas (high c12ag)
!r_c12_ag_o16   'c12ag_deboer_sigma_0p5_2000_Tgrid.dat'
!r_c12_ag_o16   'c12ag_deboer_sigma_1p0_2000_Tgrid.dat'
!r_c12_ag_o16   'c12ag_deboer_sigma_1p5_2000_Tgrid.dat'
!r_c12_ag_o16   'c12ag_deboer_sigma_2p0_2000_Tgrid.dat'
!r_c12_ag_o16   'c12ag_deboer_sigma_2p5_2000_Tgrid.dat'
!r_c12_ag_o16   'c12ag_deboer_sigma_3p0_2000_Tgrid.dat'

! negative sigmas  (lo c12ag)
!r_c12_ag_o16   'c12ag_deboer_sigma_m0p5_2000_Tgrid.dat'
!r_c12_ag_o16   'c12ag_deboer_sigma_m1p0_2000_Tgrid.dat'
!r_c12_ag_o16   'c12ag_deboer_sigma_m1p5_2000_Tgrid.dat'
!r_c12_ag_o16   'c12ag_deboer_sigma_m2p0_2000_Tgrid.dat'
!r_c12_ag_o16   'c12ag_deboer_sigma_m2p5_2000_Tgrid.dat'
!r_c12_ag_o16   'c12ag_deboer_sigma_m3p0_2000_Tgrid.dat'

A typical evolution of the stellar model should look something like the following

The $^{14}$N in the stellar core is quickly converted into $^{20}$Ne, however this is not completely correct. A more detailed nuclear reaction network would accurately capture that $^{14}$N should convert to $^{22}$Ne via the reaction sequence $^{14}$N$(\alpha,\gamma)^{18}$F$(e^{+},\nu)^{18}$O$(\alpha,\gamma)^{22}$Ne.

This remaining $^{22}$Ne persists in stellar models that evolve to become C/O white dwarfs, and is potentially detectible in their pulsation periods, see Chidester et al. 2021.

In massive stellar models, this Ne is provides the neutron excess necessary for a variety of s-process reactions, and influences the budget of Neutrons available for weak nuclear reactions during advanced burning.

The results in the exell document should look something like this

Below are profiles of the same stellar model at core-Helium depletion, with differing $^{12}$C$(\alpha,\gamma)^{16}$O rates

It is important to highlight that while the CO and He core masses are not substantially different, the mass fraction of Carbon available for carbon burning is substantially altered, and this could have larger consequences for the manner in which Carbon ignites in these massive stellar cores. It has long been thought that the nature in which Carbon ignites, radiative versus convectively, is coupled intimately with the resulting presupernova structure of stellars models, and ultimately whether they will collapse into a neutron star or black hole. See Timmes et al., Sukhbold & Woosley 2014, Sukhbold et al. 2018, Sukhbold & Adams 2020.

In lower mass models which form C/O white dwarfs, the change in composition resulting from uncertainties in the c12ag nuclear reaction rate can be constrained through asteroseismology of gravity modes in white dwarfs. See Chidester et al. 2022 Chidester et al. 2023.

Connecting to Population synthesis models

Many population-synthesis calculations do not evolve full stellar structure for every stellar model. Instead, they map pre-collapse core properties (typical at core-C depletion) to compact-remnant masses using analytic prescriptions. Here we connect our MESA models to the smooth remnant-mass prescription in Fryer et al. 2022, which is also used in population-synthesis studies such as Olejak et al. 2022.

Definitions:

• $M_{\mathrm{CO}}$: CO-core mass at collapse (from your MESA model),
• $f_{\mathrm{mix}}$: mixing parameter in the Fryer+2022 fit,
• $M_{\mathrm{crit}}$: critical CO-core mass scale (use $M_{\mathrm{crit}} = 5.75\,M_\odot$),
• $M_{\mathrm{collapse}}$: total mass at collapse.

Using Fryer+2022, the baryonic remnant mass is

$ M_{\mathrm{rem}}^{(\mathrm{raw})} = 1.2 + 0.05\,f_{\mathrm{mix}} + 0.01\left(\frac{M_{\mathrm{CO}}}{f_{\mathrm{mix}}}\right)^2 + \exp\left[f_{\mathrm{mix}}(M_{\mathrm{CO}}-M_{\mathrm{crit}})\right] $

and then

$ M_{\mathrm{rem}} = \min\left(M_{\mathrm{rem}}^{(\mathrm{raw})},\,M_{\mathrm{collapse}}\right) $

This is the quantity we plot below. In this lab, changing the $^{12}$C$(\alpha,\gamma)^{16}$O rate shifts the final core structure (especially $M_{\mathrm{CO}}$), which then shifts $M_{\mathrm{rem}}$. That is the direct link from reaction-rate uncertainty to compact-remnant predictions used in population synthesis.

TASK
For your chosen $\sigma_{C12}$ model, extract the final values for co_core_mass, he_core_mass, center_c12, and center_o16 from your MESA output.
Compute $M_{\mathrm{rem}}$ using the equations above for $f_{\mathrm{mix}}=0.7$ and add it to the google spreadsheet.
Look at the plot of $M_{\mathrm{rem}}$ vs. $\sigma_{C12}$ in the google spreadsheet.
Briefly discuss what this implies for NS/BH outcomes in population-synthesis models.
Compare the results from the spreadsheet with the plots in the Answers block below.

Warning!
The actual final value of $M_{CO}$ at core-collapse can be influenced heavily by the shell mergers during Carbon or Oxygen burning, if they occur. See Sukhbold & Woosley 2014 and Laplace et al. 2025. See also the comments for the convective overshooting &controls in inlist_common.
$M_{CO}$ is not a self consistent proxy for the Chandrasekar mass of the collapsing core. Instead one must actually look at the structure, and specifically the electron fraction $Y_{e}$ of the collapsing core Boccioli et al. 2024. The actual relation between $M_{CO}$ and $M_{ch}$ and the final remenant mass $M_{rem}$ is more complicated than the simplified picture presented in this lab (and often assumed in population synthesis codes). Even then, accurately capturing the mapping between $M_{CO}$, $M_{ch}$, and $M_{rem}$ remains an active area of research.

Answers: $M_{\mathrm{rem}}$ for $^{12}$C$(\alpha,\gamma)^{16}$O reaction rate

Bonus TASK
Change the initial mass of your stellar model in the range between (15,25) M_\odot$ and run to core-Helium depletion again.
For your chosen $\sigma_{C12}$ model, extract the final values for co_core_mass, he_core_mass, center_c12, and center_o16 from your MESA output.
Compute $M_{\mathrm{rem}}$ using the equations for $f_{\mathrm{mix}}=0.7$ and add it to the google spreadsheet.

HINT
Since we have been loading in a a 15 M$_\odot$ ZAMS model, you’ll need to instead create a new model by commenting out load_saved_model = .true. or change it to .false. in the &star_job section of inlist_project
you’ll also want to add the following line to your inlist in the &star_job section of inlist_project: create_pre_main_sequence_model = .true.