Explosive Nucleosynthesis

Here we’ll explore burning up to iron group. To do this, we’ll create a moderately-sized network that includes links from He up to iron-group. And to make the integration easier, we’ll do a few approximations along the way.

import pynucastro as pyna

Assembling our library

We start with a list of nuclei, including all \(\alpha\)-nuclei up to \({}^{56}\mathrm{Ni}\). We also add the intermediate nuclei that would participate in \((\alpha,p)(p,\gamma)\) reactions, as well as a few more nuclei, including those identified by Shen & Bildsten 2009 for bypassing the \({}^{12}\mathrm{C}(\alpha,\gamma){}^{16}\mathrm{O}\) rate.

nuclei = ["p",
          "he4", "c12", "o16", "ne20", "mg24", "si28", "s32",
          "ar36", "ca40", "ti44", "cr48", "fe52", "ni56", "zn60",
          "al27", "p31", "cl35", "k39", "sc43", "v47",
          "mn51", "co55", "cu59",
          "n13", "n14", "f18", "ne21", "na22", "na23"]

For our initial library, we take all of the ReacLib rates that link these.

reaclib_lib = pyna.ReacLibLibrary()
core_lib = reaclib_lib.linking_nuclei(nuclei)

Now we’ll add some more nuclei in the iron group.

iron_peak = ["n", "p", "he4",
             "mn51",
             "fe52", "fe53", "fe54", "fe55", "fe56",
             "co55", "co56", "co57",
             "ni56", "ni57", "ni58",
             "cu59", "zn60"]

We want the rates that connect these from ReacLib as well as tabulated weak rate compilations.

iron_reaclib = reaclib_lib.linking_nuclei(iron_peak)

weak_lib = pyna.TabularLibrary()
iron_weak_lib = weak_lib.linking_nuclei(iron_peak)

warning: Cu59 was not able to be linked in TabularLibrary
warning: He4 was not able to be linked in TabularLibrary
warning: Fe53 was not able to be linked in TabularLibrary
warning: Mn51 was not able to be linked in TabularLibrary
warning: Ni58 was not able to be linked in TabularLibrary
warning: Fe52 was not able to be linked in TabularLibrary
warning: Fe54 was not able to be linked in TabularLibrary
warning: Zn60 was not able to be linked in TabularLibrary

all_lib = core_lib + iron_reaclib + iron_weak_lib

Detailed balance

Finally, we replace the reverse rates from ReacLib by rederiving them via detailed balance, and including the partition functions.

rates_to_derive = all_lib.backward().get_rates()

# now for each of those derived rates, look to see if the pair exists

for r in rates_to_derive:
    fr = all_lib.get_rate_by_nuclei(r.products, r.reactants)
    if fr:
        all_lib.remove_rate(r)
        d = pyna.DerivedRate(rate=fr, compute_Q=True, use_pf=True)
        all_lib.add_rate(d)

Removing duplicates

There will be some duplicate rates now because we pulled rates both from ReacLib and from tabulated sources. Here we keep the tabulated version of any duplicate rates.

all_lib.eliminate_duplicates()

Creating the network and rate approximations

Now that we have our library, we can make a network.

net = pyna.PythonNetwork(libraries=[all_lib])

Next, we will do the \((\alpha,p)(p,\gamma)\) approximation and eliminate the intermediate nuclei

net.make_ap_pg_approx(intermediate_nuclei=["cl35", "k39", "sc43", "v47"])
net.remove_nuclei(["cl35", "k39", "sc43", "v47"])

We will also approximate some of the neutron captures:

net.make_nn_g_approx(intermediate_nuclei=["fe53", "fe55", "ni57"])
net.remove_nuclei(["fe53", "fe55", "ni57"])

and finally, make some of the protons into NSE protons

# make all rates with A >= 48 use NSE protons
net.make_nse_protons(48)

Visualizing the network

Let’s visualize the network

fig = net.plot(rotated=True, hide_xalpha=True,
               size=(1500, 450),
               node_size=550, node_font_size=11)

Burning to iron-group

We’ll write out the network to a python module and do a test integration with it.

net.write_network("he_burn.py")

/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: C12 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: N13 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: N14 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: p_nse partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))

import he_burn

from scipy.integrate import solve_ivp
import numpy as np

We’ll pick thermodynamic conditions that are very dense and hot–we should enter NSE here. And weak rates will also have a big effect, meaning that \(Y_e\) should evolve.

Initially, we have \(Y_e = 1/2\).

rho = 1.e9
T = 5.e9

X0 = np.zeros(he_burn.nnuc)
X0[he_burn.jhe4] = 0.98
X0[he_burn.jc12] = 0.02

Y0 = X0/he_burn.A

from pynucastro.screening import chugunov_2007

We’ll specific the time-points where we want to store the solution.

tmax = 10.0
tmin = 1.e-11
factor = 6
n = factor * int(np.log10(tmax) - np.log10(tmin)) + 1
t_eval = np.logspace(np.log10(tmin), np.log10(tmax), n, endpoint=True)

sol = solve_ivp(he_burn.rhs, [0, tmax], Y0,
                method="BDF", jac=he_burn.jacobian,
                t_eval=t_eval, args=(rho, T, chugunov_2007),
                rtol=1.e-5, atol=1.e-8)

Was the integration successful?

sol.success

True

Plotting the evolution

import matplotlib.pyplot as plt

Here we only plot the most abundant species

fig, ax = plt.subplots()

for i in range(he_burn.nnuc):
    Xmax = sol.y[i, :].max() * he_burn.A[i]
    if Xmax > 1.e-1:
        ax.loglog(sol.t, sol.y[i,:] * he_burn.A[i], linewidth=2,
                  label=f"X({he_burn.names[i].capitalize()})")
    elif Xmax > 1.e-2:
        ax.loglog(sol.t, sol.y[i,:] * he_burn.A[i], linewidth=1,
                  label=f"X({he_burn.names[i].capitalize()})")
    elif Xmax > 1.e-3:
        ax.loglog(sol.t, sol.y[i,:] * he_burn.A[i], linewidth=1, ls="--",
                  label=f"X({he_burn.names[i].capitalize()})")
        
ax.set_ylim(1.e-5, 1.1)
ax.legend(fontsize="small", ncol=4, loc=3)
ax.set_xlabel("t (s)")
ax.set_ylabel("X")
ax.margins(0)
ax.grid(ls=":")
fig.set_size_inches(8, 6)

\(Y_e\) evolution

X_final = sol.y[:, -1] * he_burn.A
comp = pyna.Composition(net.unique_nuclei)
comp.set_array(X_final)

Notice that \(Y_e\) has changed a lot due to electron captures.

comp.ye

0.4750014049084172

Visualizing network flow

Finally, we can visualize the evolution of the flow through the network. We’ll just work with a few of the times that were saved. We’ll also show the net rate between nuclei (forward - reverse rate) to make it clearer.

for idx in range(0, n, 3*factor):
    t = sol.t[idx]
    Y = sol.y[:, idx]
    comp.set_array(Y * he_burn.A)
    fig = net.plot(rho, T, comp,
                   rotated=True, hide_xalpha=True,
                   size=(1500, 600),                
                   ydot_cutoff_value=1.e-15,
                   use_net_rate=True,
                   color_nodes_by_abundance=True,
                   Z_range=(1, 31),
                   node_size=500, node_font_size=10)
    fig.suptitle(f"t = {t} s")

/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: C12 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: N13 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: N14 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: p_nse partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))

/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: C12 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: N13 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: N14 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pynucastro/rates/derived_rate.py:135: UserWarning: p_nse partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))