Combining ReacLib Rates with Electron Capture Tables

Here’s an example of using tabulated weak rates from Suzuki et al. (2016) together with rates from the ReacLib database.

We’ll build a network suitable for e-capture supernovae.

import pynucastro as pyna

First create a library from ReacLib using a set of nuclei.

reaclib_library = pyna.ReacLibLibrary()

all_nuclei = ["p", "he4",
              "ne20", "o20", "f20",
              "mg24", "al27", "o16",
              "si28", "s32", "p31"]
ecsn_library = reaclib_library.linking_nuclei(all_nuclei,with_reverse=True)

Here are the rates it chose

print(ecsn_library)

O16 + He4 ⟶ Ne20 + 𝛾           [Q =   4.73 MeV] (O16 + He4 --> Ne20 <co10_reaclib__>)
O16 + O16 ⟶ He4 + Si28         [Q =   9.59 MeV] (O16 + O16 --> He4 + Si28 <cf88_reaclib__>)
O16 + O16 ⟶ p + P31            [Q =   7.68 MeV] (O16 + O16 --> p + P31 <cf88_reaclib__>)
O20 ⟶ F20 + e⁻ + 𝜈             [Q =   3.81 MeV] (O20 --> F20 <wc12_reaclib_weak_>)
F20 ⟶ Ne20 + e⁻ + 𝜈            [Q =   7.02 MeV] (F20 --> Ne20 <wc12_reaclib_weak_>)
Ne20 + He4 ⟶ Mg24 + 𝛾          [Q =   9.32 MeV] (Ne20 + He4 --> Mg24 <il10_reaclib__>)
Mg24 + He4 ⟶ Si28 + 𝛾          [Q =   9.98 MeV] (Mg24 + He4 --> Si28 <st08_reaclib__>)
Al27 + p ⟶ He4 + Mg24          [Q =   1.60 MeV] (Al27 + p --> He4 + Mg24 <il10_reaclib__>)
Al27 + p ⟶ Si28 + 𝛾            [Q =  11.59 MeV] (Al27 + p --> Si28 <il10_reaclib__>)
Al27 + He4 ⟶ P31 + 𝛾           [Q =   9.67 MeV] (Al27 + He4 --> P31 <ths8_reaclib__>)
Si28 + He4 ⟶ S32 + 𝛾           [Q =   6.95 MeV] (Si28 + He4 --> S32 <ths8_reaclib__>)
P31 + p ⟶ He4 + Si28           [Q =   1.92 MeV] (P31 + p --> He4 + Si28 <il10_reaclib__>)
P31 + p ⟶ S32 + 𝛾              [Q =   8.86 MeV] (P31 + p --> S32 <il10_reaclib__>)
Ne20 ⟶ He4 + O16               [Q =  -4.73 MeV] (Ne20 --> He4 + O16 <co10_reaclib__reverse>)
Mg24 + He4 ⟶ p + Al27          [Q =  -1.60 MeV] (Mg24 + He4 --> p + Al27 <il10_reaclib__reverse>)
Mg24 ⟶ He4 + Ne20              [Q =  -9.32 MeV] (Mg24 --> He4 + Ne20 <il10_reaclib__reverse>)
Si28 + He4 ⟶ O16 + O16         [Q =  -9.59 MeV] (Si28 + He4 --> O16 + O16 <cf88_reaclib__reverse>)
Si28 + He4 ⟶ p + P31           [Q =  -1.92 MeV] (Si28 + He4 --> p + P31 <il10_reaclib__reverse>)
Si28 ⟶ He4 + Mg24              [Q =  -9.98 MeV] (Si28 --> He4 + Mg24 <st08_reaclib__reverse>)
Si28 ⟶ p + Al27                [Q = -11.59 MeV] (Si28 --> p + Al27 <il10_reaclib__reverse>)
P31 + p ⟶ O16 + O16            [Q =  -7.68 MeV] (P31 + p --> O16 + O16 <cf88_reaclib__reverse>)
P31 ⟶ He4 + Al27               [Q =  -9.67 MeV] (P31 --> He4 + Al27 <ths8_reaclib__reverse>)
S32 ⟶ He4 + Si28               [Q =  -6.95 MeV] (S32 --> He4 + Si28 <ths8_reaclib__reverse>)
S32 ⟶ p + P31                  [Q =  -8.86 MeV] (S32 --> p + P31 <il10_reaclib__reverse>)

Now let’s specify the weak rates we want from Suzuki et al. – these tables are included with pynucastro:

tabular_lib = pyna.TabularLibrary()

ecsn_tabular = tabular_lib.linking_nuclei(["f20", "o20", "ne20"])

Let’s look at the rates

tr = ecsn_tabular.get_rates()
tr

[F20 ⟶ Ne20 + e⁻ + 𝜈,
 F20 + e⁻ ⟶ O20 + 𝜈,
 Ne20 + e⁻ ⟶ F20 + 𝜈,
 O20 ⟶ F20 + e⁻ + 𝜈]

These tables are in terms of \(T\) and \(\rho Y_e\). We can easily plot the tabulated electron capture rates just like ReacLib rates:

fig = tr[0].plot(rhoYmin=4.e8, rhoYmax=2.e9, Tmin=1.e8, Tmax=2.e9,
                 figsize=(8, 4))

Let’s create a rate collection from just the ReacLib rates and look to see which are actually important

rc = pyna.RateCollection(libraries=[ecsn_library])

Here we’ll pick thermodynamic conditions appropriate to the oxygen burning shell in a white dwarf

comp = pyna.Composition(rc.get_nuclei())
comp.set_nuc("o16", 0.5)
comp.set_nuc("ne20", 0.3)
comp.set_nuc("mg24", 0.1)
comp.set_nuc("o20", 1.e-5)
comp.set_nuc("f20", 1.e-5)
comp.set_nuc("p", 1.e-5)
comp.set_nuc("he4", 1.e-2)
comp.set_nuc("al27", 1.e-2)
comp.set_nuc("si28", 1.e-2)
comp.set_nuc("s32", 1.e-2)
comp.set_nuc("p31", 1.e-2)
comp.normalize()

fig = rc.plot(rho=7.e9, T=1.e9, comp=comp, ydot_cutoff_value=1.e-20)

Note, this rate collection already includes some weak rates from ReacLib – we want to remove those.

We’ll also take the opportunity to remove any rates that are not important, by looking for a \(|\dot{Y}| > 10^{-20}~s^{-1}\).

new_rate_list = []
ydots = rc.evaluate_rates(rho=7.e9, T=1.e9, composition=comp)
for rate in rc.rates:
    if abs(ydots[rate]) >= 1.e-20 and not rate.weak:
        new_rate_list.append(rate)        

Now let’s add the tabular rates to this pruned down rate list

new_rate_list += tr

Finally, we’ll create a new rate collection by combining this new list of rates with the list of tables:

rc_new = pyna.RateCollection(rates=new_rate_list)

fig = rc_new.plot(rho=7.e9, T=1.e9, comp=comp, curved_edges=True, rotated=True)

We can see the values of the rates at our thermodynamic conditions easily:

rc_new.evaluate_rates(rho=7.e9, T=1.e9, composition=comp)

{Ne20 ⟶ He4 + O16: 4.6316695198155e-18,
 O16 + He4 ⟶ Ne20 + 𝛾: 2163.4896055892,
 Ne20 + He4 ⟶ Mg24 + 𝛾: 470.46571662959457,
 Mg24 + He4 ⟶ Si28 + 𝛾: 65.33425925449085,
 Al27 + p ⟶ Si28 + 𝛾: 3786.588954080847,
 Al27 + He4 ⟶ P31 + 𝛾: 0.05763581125914332,
 Si28 + He4 ⟶ S32 + 𝛾: 0.08322738610569612,
 P31 + p ⟶ S32 + 𝛾: 2375.1304964039123,
 O16 + O16 ⟶ p + P31: 7.15493989546395e-17,
 O16 + O16 ⟶ He4 + Si28: 2.5949312230335824e-17,
 Mg24 + He4 ⟶ p + Al27: 0.2543862399303773,
 Al27 + p ⟶ He4 + Mg24: 5920.818449985259,
 Si28 + He4 ⟶ p + P31: 0.00019528201299493384,
 P31 + p ⟶ He4 + Si28: 5491.960442675607,
 F20 ⟶ Ne20 + e⁻ + 𝜈: 1.143049854823845e-14,
 F20 + e⁻ ⟶ O20 + 𝜈: 4.212096660029042e-09,
 Ne20 + e⁻ ⟶ F20 + 𝜈: 4.3100255627801964e-11,
 O20 ⟶ F20 + e⁻ + 𝜈: 2.688308786951095e-28}

Additionally, we can see the specific energy generation rate for these conditions

rc_new.evaluate_energy_generation(rho=7.e9, T=1.e9, composition=comp)

9.667025588667893e+22