Nuclear Reactions In MESA: Learning By Example


To get a feel for the basics of the nuclear reaction modules in MESA, I've examined the contents of MESA's net module. In particular, the contents of net's test subdirectory mesa/net/test. I try to balance going into detail and outlining the overall picture. It's my hope that these notes can be used to give one an introduction to reactions in MESA such that they'll be capable of using MESA's nuclear reaction modules in their own code.

net's test directory represents a typical "work" directory, similar to what one would use in carrying out a stellar evolution calculation with MESA. It contains make and src subdirectories:

In [22]:
%cd ~/Research/Codebase/mesa/net/test
In [23]:
burn*        inlist_one_zone_burn             rn*
burn_net.rb  make/                            sample*
ck*          mk*                              src/
clean*       mkx*                             tester*
cleanup*     net_support.rb                   test_output
diff.txt     one_zone_burn_T_Rho_history.txt  test_quietly*
export*      plotter*                         tmp.txt
In [26]:
ls src
mod_test_net.f   sample_net.f         test_net.f
one_zone_burn.f  test_burn_const_P.f  test_net_quietly.f
plot_net.f       test_burn.f          test_net_support.f
In [25]:
ls make
makefile                  plot_net.o             test_burn.o
mod_one_zone_support.mod  sample_net.o           test_net.o
mod_test_net.mod          test_burn_const_p.mod  test_net_quietly.o
mod_test_net.o            test_burn_const_P.o    test_net_support.mod
one_zone_burn.o           test_burn.mod          test_net_support.o

To see how the files in these directories relate, what each's purpose is, and how they piece together MESA modules to do a reaction calculation we'll examine the execution of tester which results from making the program test_net.f found in the src subdirectory. Since the tester program only makes use of local files test_net.f, mod_test_net.f, and test_net_support.f we won't be examining the other files in src. tester is the program executed when we do the standard clean, make, run (output from cleaning/making has been suppressed):

In [34]:
hide = !./cleanup
hide = !./mk
  **************** Basic **************** 

                                log temp         7.500000000
                                 log rho         2.000000000

                            log(eps_nuc)         9.360791959

                                 eps_nuc     0.229504898E+10

                            d_lneps_dlnT        12.186236056
                          d_lneps_dlnRho         1.029292229

 energy generation by category

                                category   log rate (erg/g/sec)

                                      pp         9.360507   0.229354E+10
                                     cno         6.178073   0.150686E+07

             log10(-photodisintegration)       -99.000000  -0.000000E+00

  **************** o18_and_ne22 **************** 

                                log temp         7.500000000
                                 log rho         2.000000000

                            log(eps_nuc)         6.177786613

                                 eps_nuc     0.150586699E+07

                            d_lneps_dlnT        14.551318478
                          d_lneps_dlnRho         1.048016090

 energy generation by category

                                category   log rate (erg/g/sec)

                                     cno         6.177744   0.150572E+07
                                      pp         2.164895   0.146182E+03

             log10(-photodisintegration)       -99.000000  -0.000000E+00

  **************** pp_extras **************** 

                                log temp         7.500000000
                                 log rho         2.000000000

                            log(eps_nuc)         9.359238491

                                 eps_nuc     0.228685428E+10

                            d_lneps_dlnT        12.184186914
                          d_lneps_dlnRho         1.029292264

 energy generation by category

                                category   log rate (erg/g/sec)

                                      pp         9.358952   0.228535E+10
                                     cno         6.178073   0.150686E+07

             log10(-photodisintegration)       -99.000000  -0.000000E+00

  **************** cno_extras **************** 

                                log temp         7.500000000
                                 log rho         2.000000000

                            log(eps_nuc)         5.114771971

                                 eps_nuc     0.130248272E+06

                            d_lneps_dlnT         2.873242105
                          d_lneps_dlnRho         0.206951151

 energy generation by category

                                category   log rate (erg/g/sec)

                                     cno         5.114771   0.130248E+06

             log10(-photodisintegration)       -99.000000  -0.000000E+00

This program calculates the rates for several different reaction networks and displays the total energy generation rate $\dot{\epsilon}$ = eps_nuc [erg/g/s], given a certain temperature and density. The details of the various networks can be found in $MESA_DIR/data/net_data/nets. We'll focus on the details of the 'Basic' network to explore how reaction networks are implemented in MESA.

The calculation procedes in two stages: first the network is established, second the network is evaluated to determine the quantities of interest.

Stage 1: Establishing the Network

From inspecting the file $MESA_DIR/data/net_data/nets/ we see that this basic network provides all branches of the PP chains, two of the CNO cycles, and helium/alpha-chain burning up to $^{24}\rm{Mg}$. In addition a few weak reactions and alternative reactions are included. A great resource for details about nuclear reaction networks is the website of one of the developers of MESA, Frank Timmes.

The call tree of the most salient subroutines for this phase is

--> MTN:test()
----> TNS:test_net_setup()
------> NL:net_init()
------> NL:alloc_net_handle()
------> NL:net_start_def()
------> NL:read_net_file()
------> NL:net_finish_def()

------> NL:net_setup_tables()
--------> RL:read_raw_rates_records()
----------> RI:init_raw_rates_records()
----------> RI:set_rate_ids()
--------> NI:alloc_net_general_info()
----------> RL:make_rate_tables()
------------> RS:do_make_rate_tables()
--------------> RS:read_reaction_from_cache()
--------------> RS:get_net_rates_for_tables()
----------------> RR:set_raw_rates()
------------------> RR:set_raw_rate()
--------------------> RR:do1_of_3()
----------------------> RR:eval_which_raw_rate()
------------------------> RR:eval_raw_rate()
--------------------------> RR:eval_table()
--------------------------> RR:rate_fcn(), e.g. R:rate_tripalf_nacre()

------> NL:get_chem_id_table()
------> NL:get_net_iso_table()
------> NL:get_reaction_id_table()
------> NL:get_net_reaction_table()
----> TNS:do_test_net()
------> TNS:set_composition()

In the above tree the capitals are shorthand for the files containing the subroutine. The full file locations are
TN = $MESA_DIR/net/test/src/test_net.f
MTN = $MESA_DIR/net/test/src/mod_test_net.f
TNS = $MESA_DIR/net/test/src/test_net_support.f
NL = $MESA_DIR/net/public/net_lib.f
RL = $MESA_DIR/rates/public/rates_lib.f
RI = $MESA_DIR/rates/private/rates_initialize.f
NI = $MESA_DIR/net/private/net_initialize.f
RS = $MESA_DIR/rates/private/rates_support.f
RR = $MESA_DIR/rates/private/raw_rates.f
R = $MESA_DIR/rates/private/ratelib.f

Let's step through the series of calls establishing the network.

net_init() through net_finish_def() initialize the net module and read in the network file (in this case, which populates a list of isotopes and possible reactions involving those isotopes. This information is associated with a "handle," which is the index in an array of General_Net_Info types maintained in $MESA_DIR/net/public/net_def.f. A given reaction network has an associated handle/General_Net_Info type. When read_net_file() is called the file is used to initialize the handle's associated General_Net_Info type with the isotopes and reactions. Up to this point we've only named the isotopes and reactions. No physical data such as properties of the isotopes or reaction rates have been determined.

net_setup_tables() reads data from
$MESA_DIR/data/rates_data/rate_tables/rate_list.txt which lists the other files in that directory r1.txt, r2.txt, etc... These are dummy examples of rate tables. The basic structure is that rate_list.txt lists all the rate names and their corresponding files. So the rates for reaction r1 is found in r1.txt. Inside that file will be a list of ($T_8$, rate) pairs. I can't find any non-dummy rate tables, maybe tables aren't often used in practice? After that, set_rate_ids() sets a few rate IDs defined in $MESA_DIR/rates/public/rates_def.f. Honestly, both of these steps seem to not accomplish much.

On the other hand, alloc_net_general_info() is perhaps the most important part of establishing the network. This builds the rate tables and associates them with the network's handle. Where does it get the information to build the tables from? Well, often the table data for a given reaction (e.g. the triple alpha reaction, which in MESA is called r_he4_he4_he4_to_c12) will already be cached in binary files found in $MESA_DIR/data/rates_data/cache.

However, we learn more by considering the case that the reaction rate table isn't cached. As can be seen in the call tree, this ultimately leads to a call to set_raw_rate() which consists of a huge select case statement listing almost all of the reactions MESA knows about (1117 in total at time of writing). There's a more complete list in $MESA_DIR/data/rates_data/reactions.list, which contains 1246 rates at time of writing. The reaction rate for each can be calculated using a number of sources: NACRE rates (C. Angulo et al., Nucl. Phys. A656 (1999)3-187), JINA reaclib rates (, or Frank (X.) Timmes' tables/code (, which in turn are based on several sources cited on the website, notably Caughlin, G. R. & Fowler, W. A. 1988, Atom. Data and Nuc. Data Tables, 40, 283). For all of these there exists a function with a suffix of _nacre, _jina OR _reaclib, or _fxt respectively. These functions can be found in $MESA_DIR/rates/private/ratelib.f.

Let's consider the example of the triple alpha reaction rate. In this test code a parameter is set that preferences the NACRE rates over others. The rate is evaluated by calling rate_tripalf_nacre(). The NACRE paper/data provide both numerically computed tables of rates (which will be the most accurate rates) and analytic approximations of the rates. This test, at least, doesn't access any tables for NACRE rates but instead relies on the analytic approximations. I'm not sure if MESA has the NACRE tables somewhere or not (though it certainly has the ability to load in rates from such tables if you make them yourself).

NACRE's Table 3 of analytic approximations serves to provide values for the constants in the general expression

$N_A \langle \sigma v \rangle = a_0 T_9^{2/3} \exp(-a_1 T_9^{1/3} - (a_2 T_9)^2) * (1 + b_0 T_9 + b_1 T_9^2 + b_2 T_9^3 + b_3 T_9^4 + b_4 T_9^5)$
$+ c_0 T_9^{3/2} e^{-c_1/T_9}$
$+ d_0 T_9^{3/2} e^{-d_1/T_9}$
$+ e_0 T_9^{e_1} e^{-e_2/T_9}~[\rm{cm}^{3(b-1)}~\rm{mol}^{-(b-1)}~\rm{s}^{-1}]$,

I'm trying to be general here by including the interger b for a b-body reaction in the units. The calculations for reactions that aren't 2-body typically lead to a rate with these units. The function rate_tripalf_nacre() provides the constants to calculate rates for $^4\rm{He(}\alpha, \gamma\rm{)}^8\rm{Be}$ and $^8\rm{Be(}\alpha, \gamma\rm{)}^{12}\rm{C}$ which are multiplied to get the net triple alpha reaction (in units of $[\rm{cm}^{6}~\rm{mol}^{-2}~\rm{s}^{-1}]$). This rate calculation is carried out for a range of temperatures in do_make_rate_tables(), generating an array of rates for each reaction at discrete temperature steps which is associated with the network handle/General_Net_Info type: g%rattab. A table of interpolants, g%rattab_f1, is also calculated since you'll usually have to interpolate to match the temperature of your fluid element from the discrete temperature steps available in g%rattab.

get_chem_id_table() through get_net_reaction_table() are necessary because MESA's chem, rates, and other modules have their own collections of isotopes and reactions. This series of calls essentially provides a translation from the IDs and such of a particular network to the IDs and such of the larger set of isotopes/reactions MESA knows about.

Finally, set_composition() initializes the composition for the basic net to what I believe is the expected composition in the Sun at the site of active burning (NOT surface composition).

After all this we've allocated room for and initialized some components of a General_Net_Info representation of our basic network as well as for some other shared variables in the test_net_support module. Now for the calculations!

Stage 2: Evaluating the Network

As before, we'll ground ourselves with the call tree:

--> MTN:test()
----> TNS:test_net_setup() (see previous section)
----> TNS:do_test_net()
------> TNS:Do_One_Net()
--------> CL:composition_info()
--------> NL:net_get()
----------> NE:eval_net()

------------> NI:setup_net_info()

------------> NE:get_weaklib_rates()
--------------> WL:eval_weak_reaction_info()
----------------> EW:do_eval_weak_reaction_info()

------------> NE:get_ecapture_rates()
--------------> EL:eval_ecapture_reaction_info()
----------------> EE:do_eval_ecapture_reaction_info()

------------> NE:set_molar_abundances()
------------> NE:get_rates_with_screening()
--------------> RL:eval_using_rate_tables()
----------------> RS:do_get_raw_rates()
------------------> RS:get_rates_from_table()
--------------> NS:screen_net()

------------> ND:get_derivs()
--------------> ND:get1_derivs()
----------------> NDS:do_three_one()
------------------> NDS:do_three_one_neu()

Where: TN = $MESA_DIR/net/test/src/test_net.f
MTN = $MESA_DIR/net/test/src/mod_test_net.f
TNS = $MESA_DIR/net/test/src/test_net_support.f
CL = $MESA_DIR/chem/public/chem_lib.f
NL = $MESA_DIR/net/public/net_lib.f
NE = $MESA_DIR/net/private/net_eval.f
NI = $MESA_DIR/net/private/net_initialize.f
WL = $MESA_DIR/weaklib/public/weak_lib.f
EW = $MESA_DIR/weaklib/private/eval_weak.f
EL = $MESA_DIR/ecapture/public/ecapture_lib.f
EE = $MESA_DIR/ecapture/private/eval_ecapture.f
RL = $MESA_DIR/rates/public/rates_lib.f
RS = $MESA_DIR/rates/private/rates_support.f
NS = $MESA_DIR/net/private/net_screen.f
ND = $MESA_DIR/net/private/net_derivs.f
NDS = $MESA_DIR/net/private/net_derivs_support.f

Be warned that the test_net_support module makes ample use of shared data and function side-effects. So, for example, if you're looking at the code trying to figure out where in the heck the g variable came from in do_test_net(), you have to realize g is a General_Net_Info pointer shared by all subroutines in test_net_support and was initialized in phase 1 (test_net_setup()). Confusingly, even though Do_One_Net() is a part of the test_net_support module and can thus see g, it has g passed to it. Perhaps there are justifiable reasons for this confusing set of programming practices, but I wish shared data and function side-effects weren't relied upon so much and so inconsistently. In any case, when something's confusing it's probably because it's a shared variable and you'll have to figure out how it got initialized and what it got initialized with by combing through the various files in $MESA_DIR/net/test/src.

Warning out of the way, let's step through the process of evaluating the network.

setup_net_info() sets up a Net_Info variable n which serves as working storage during nuclear reaction calculations (whereas General_Net_Info types store information that should be constant for a given network).

I'm not going to detail the calculations of get_weaklib_rates() or get_ecapture_rates(). Instead I'll detail get_rates_with_screening() and hope the others are a similar procedure that instead uses the weaklib and ecapture modules instead of the rates module.

get_rates_with_screening()'s first major step is to make use of the temperature-dependent raw reaction rates ($N_A \langle \sigma v \rangle(T)$) that were stored in tables in the establishment phase of this calculation. This is made into a density-dependent raw rate via

$RR = Y_e^a \rho^{b'} N_A \langle \sigma v \rangle(T) ~[\rm{g}^{b-1}~\rm{mol}^{-(b-1)}~\rm{s}^{-1}]$,

where $Y_e$ is the electron fraction, $\rho$ the density. The coefficients ($a, b'$) are as follows (note that $b' = b-1$, where $b$ is the number of reaction inputs):
(0, 0) for photodisintegration or decay
(0, 1) for standard two-body reactions
(0, 2) for three-body reactions such as triple alpha
(1, 1) for two-body electron captures
(1, 2) for three-body electron captures such as proton-electron-proton (pep)

Though not used in this test, I should note there's also an option to scale the raw reaction rate $RR$ by an arbitrary factor. This is sometimes useful, for example when calculating the $^{12}\rm{C(}\alpha, \gamma\rm{)}^{16}\rm{O}$ rate (see, e.g. Weaver, T. A., & Woosley, S. E. 1993, PhR, 227, 65).

$RR$ is calculated by eval_using_rate_tables() and is stored in the n%rate_raw(num_reactions) array. After this the rates are multiplied by a screening factor as calculated in screen_net() (let's call this screened rate $R$). In this test screening is calculated as detailed in Graboske, H. C., Dewitt, H. E., Grossman, A. S., & Cooper, M. S. 1973, ApJ, 181, 457. These are, appropriately enough, stored in n%rate_screened(num_reactions).

Confusingly, to me at least, the last major part of the calculation of $\dot{\epsilon}$ occurs in get_derivs(). This function calulates various derivatives that are of interest and will later be used if any burning is to occur. But in this test we're just calculating rates, not burning. For this the important aspect of get_derivs() is that it calculates

$\dot{\epsilon}' = (Q - Q_{neu}) y^b~R ~[\rm{MeV}~\rm{mol}~\rm{g}^{-1}~\rm{s}^{-1}]$,

where $Q$ is the energy released by each reaction in MeV, $Q_{neu}$ is the amount of this energy in the form of neutrinos (and will thus mostly be lost, unavailable for heating the local volume), and there's a $y$ factor for each reaction input's molar fraction [mol/g]. This is converted to $[\rm{erg}~\rm{g}^{-1}~\rm{s}^{-1}]$ in eval_net() via

$\dot{\epsilon} = \dot{\epsilon}' \times \frac{10^6 \rm{eV}}{1 \rm{MeV}} \times \frac{1.602176487 \times 10^{-12} \rm{erg}}{1~\rm{eV}} \times N_A~[\rm{erg}~\rm{g}^{-1}~\rm{s}^{-1}]$,

where $N_A$ is Avogadro's number with units $\rm{mol}^{-1}$.

And there we go! We've calculated the nuclear energy generation rate for a given reaction in cgs units. The test calculates these for each reaction in the given reaction network (e.g. and the sum of them all is what is printed out in the test's output as eps_nuc. In this test the rates are calculated for a temperature of about $30~\rm{MK}$, a density of $100~\rm{g/cm^3}$, and what I believe is the composition in the Sun at the site of active hydrogen burning.