alternate_rates

Sometimes papers provide new measurements of rates that may be better than the current ReacLib library. These rates are typically provided in one of two forms:

• As a new set of fit parameters in the ReacLib format. An example of this is the deBoer et al. [2017] \({}^{12}\mathrm{C}(\alpha, \gamma){}^{16}\mathrm{O}\) rate. This type of rate can be added by subclassing the ReacLibRate class.

• A tabulation of \(N_A \langle \sigma v\rangle\) in terms of temperature. An example of this is the Iliadis et al. [2022] \({}^{16}\mathrm{O}(p, \gamma){}^{17}\mathrm{F}\) rate. This type of rate can be added by subclassing the TemperatureTabularRate class.

We use the alternate_rates module to hold these rates.

import pynucastro as pyna

Comparison of de Boer et al. \({}^{12}\mathrm{C}(\alpha, \gamma){}^{16}\mathrm{O}\)

Here we compare the ReacLib and deBoer et al. versions of \({}^{12}\mathrm{C}(\alpha, \gamma){}^{16}\mathrm{O}\).

First get the ReacLib version

rl = pyna.ReacLibLibrary()
c12_rl = rl.get_rate_by_name("c12(a,g)o16")

We can see the source of this rate

c12_rl.source

{'Label': 'nac2',
 'Author': 'Xu, Y. et al.',
 'Title': 'NACRE2',
 'Publisher': 'Nucl. Phys. A 918, P. 61',
 'Year': '2013',
 'URL': 'https://reaclib.jinaweb.org/labels.php?action=viewLabel&label=nac2'}

This comes from the NACRE II compilation of rates.

Now the DeBoer version. This is available as DeBoerC12agO16:

from pynucastro.rates.alternate_rates import DeBoerC12agO16
c12_deboer = DeBoerC12agO16()

Note that since DeBoer et al. 2012 give the new rate as a parameterization in the ReacLib form, this new rate just subclasses ReacLibRate

DeBoerC12agO16.__base__

pynucastro.rates.reaclib_rate.ReacLibRate

Comparison with ReacLib

Let’s plot a comparison of the rates

import numpy as np
import matplotlib.pyplot as plt

Ts = np.logspace(7.5, 9.5, 100)

r_rl = np.array([c12_rl.eval(T) for T in Ts])
r_deboer = np.array([c12_deboer.eval(T) for T in Ts])

fig, ax = plt.subplots()
ax.loglog(Ts, r_rl, label="ReacLib")
ax.loglog(Ts, r_deboer, label="deBoer")
ax.set_xlabel("T")
ax.set_ylabel(r"$N_A\langle \sigma v\rangle$")
ax.legend()

<matplotlib.legend.Legend at 0x7fabd4173e00>

Here’s the relative difference

fig, ax = plt.subplots()
ax.loglog(Ts, np.abs((r_rl - r_deboer) / r_rl))
ax.set_xlabel("T")
ax.set_ylabel(r"relative difference")
ax.grid(ls=":")

Using alternate rates in a network

To create a network with this rate, we can swap it in for the ReacLib version when we create a library. Here’s a simple He-burning network.

nuclei = ["p", "he4", "c12", "o16", "ne20", "na23", "mg24"]

comp = pyna.Composition(nuclei)
comp.set_equal()

rho = 1.e6
T = 1.e9

Pure ReacLib

First we’ll create a RateCollection with just the ReacLib rates

lib = rl.linking_nuclei(nuclei, with_reverse=False)
rc = pyna.RateCollection(libraries=[lib])
rates = rc.evaluate_rates(rho, T, comp)
for rr, val in rates.items():
    print(f"{str(rr):25} : {val:12.7e}")

C12 + He4 ⟶ O16 + 𝛾       : 2.7441793e-03
O16 + He4 ⟶ Ne20 + 𝛾      : 1.1385847e+00
Ne20 + He4 ⟶ Mg24 + 𝛾     : 4.1265513e-01
Na23 + p ⟶ Mg24 + 𝛾       : 6.9202360e+05
C12 + C12 ⟶ p + Na23      : 1.1970613e-09
C12 + C12 ⟶ He4 + Ne20    : 1.4808559e-09
O16 + C12 ⟶ He4 + Mg24    : 5.5470507e-15
Na23 + p ⟶ He4 + Ne20     : 8.9374939e+06
3 He4 ⟶ C12 + 𝛾           : 2.5845076e-03

Swapping in DeBoer rate

Now we’ll swap in the deBoer rate, first using remove_rate to remove the old rate and then using add_rate to add the new rate to the library.

lib.remove_rate(c12_rl)
lib.add_rate(c12_deboer)

and recreate the new RateCollection

rc2 = pyna.RateCollection(libraries=[lib])
rates = rc2.evaluate_rates(rho, T, comp)
for rr, val in rates.items():
    print(f"{str(rr):25} : {val:12.7e}")

O16 + He4 ⟶ Ne20 + 𝛾      : 1.1385847e+00
Ne20 + He4 ⟶ Mg24 + 𝛾     : 4.1265513e-01
Na23 + p ⟶ Mg24 + 𝛾       : 6.9202360e+05
C12 + C12 ⟶ p + Na23      : 1.1970613e-09
C12 + C12 ⟶ He4 + Ne20    : 1.4808559e-09
O16 + C12 ⟶ He4 + Mg24    : 5.5470507e-15
Na23 + p ⟶ He4 + Ne20     : 8.9374939e+06
3 He4 ⟶ C12 + 𝛾           : 2.5845076e-03
C12 + He4 ⟶ O16 + 𝛾       : 2.2040004e-03

Comparing to the pure ReacLib version, we see the only difference in the evaluation of the rates is in \({}^{12}\mathrm{C}(\alpha, \gamma){}^{16}\mathrm{O}\)

Comparison of Iliadis et al. 2022 \({}^{16}\mathrm{O}(p,\gamma){}^{17}\mathrm{F}\)

The Iliadis et al. 2022 paper provides a tabulation of the \({}^{16}\mathrm{O}(p,\gamma){}^{17}\mathrm{F}\) rate in terms of \(T_9\). We use the median value of the rate.

First we’ll get the ReacLib version

o16_rl = rl.get_rate_by_name("o16(p,g)f17")
o16_rl.source

{'Label': 'ia08',
 'Author': 'Iliadis, C.',
 'Title': 'New reaction rate for 16O( p,γ)17F and its influence on the oxygen isotopic ratios in massive AGB stars',
 'Publisher': 'PRC77, 045802',
 'Year': '2008',
 'URL': 'https://reaclib.jinaweb.org/labels.php?action=viewLabel&label=ia08'}

We see that the ReacLib version of the rate is from an earlier work, Iliadis et al. 2008.

Now we’ll get the new version of the rate

from pynucastro.rates.alternate_rates import IliadisO16pgF17
o16_iliadis = IliadisO16pgF17()

and we can see the base now is a temperature tabulation:

IliadisO16pgF17.__base__

pynucastro.rates.temperature_tabular_rate.TemperatureTabularRate

The rate data is stored in the form \((\log_{10}(T_9), \log_{10}(\lambda))\)

for lT9, lr in zip(o16_iliadis.log_t9_data, o16_iliadis.log_rate_data):
    print(f"{lT9:11.8f} : {lr:11.8f}")

-2.52287875 : -40.36271045
-2.39794001 : -35.85418229
-2.30103000 : -32.64897715
-2.22184875 : -30.20551195
-2.15490196 : -28.25531594
-2.09691001 : -26.64723881
-2.04575749 : -25.28886793
-2.00000000 : -24.11947222
-1.95860731 : -23.09756194
-1.92081875 : -22.19348087
-1.88605665 : -21.38541913
-1.85387196 : -20.65698550
-1.82390874 : -19.99567863
-1.79588002 : -19.39072583
-1.74472749 : -18.32148162
-1.69897000 : -17.40142834
-1.60205999 : -15.56082526
-1.52287875 : -14.15995667
-1.39794001 : -12.12418661
-1.30103000 : -10.68151927
-1.22184875 : -9.58485965
-1.15490196 : -8.71152720
-1.09691001 : -7.99310629
-1.04575749 : -7.38742805
-1.00000000 : -6.86678054
-0.95860731 : -6.41285050
-0.92081875 : -6.01184253
-0.88605665 : -5.65384270
-0.85387196 : -5.33170729
-0.82390874 : -5.03948166
-0.79588002 : -4.77262756
-0.74472749 : -4.30181243
-0.69897000 : -3.89756629
-0.60205999 : -3.09119340
-0.52287875 : -2.48030323
-0.45593196 : -1.99524884
-0.39794001 : -1.59722293
-0.34678749 : -1.26248931
-0.30103000 : -0.97551433
-0.22184875 : -0.50501103
-0.15490196 : -0.13188479
-0.09691001 :  0.17405981
-0.04575749 :  0.43120288
 0.00000000 :  0.65147185
 0.09691001 :  1.08955188
 0.17609126 :  1.42028588
 0.24303805 :  1.68178377
 0.30103000 :  1.89558828
 0.39794001 :  2.22685757
 0.47712125 :  2.47290265
 0.54406804 :  2.66190729

and we can plot the temperature dependence

fig = o16_iliadis.plot()

Comparison with ReacLib

Ts = np.logspace(np.log10(o16_iliadis.table_Tmin), np.log10(o16_iliadis.table_Tmax), 100)
r_rl = np.array([o16_rl.eval(T) for T in Ts])
r_iliadis = np.array([o16_iliadis.eval(T) for T in Ts])

err = np.abs(r_rl - r_iliadis) / r_rl

fig, ax = plt.subplots()
ax.loglog(Ts, err)
ax.set_xlabel("T")
ax.set_ylabel(r"relative difference")
ax.grid(ls=":")

We see that the maximum difference is about 5%