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)
../_images/6914e4894c7dca42839ae1c45487c11b47a02871dd835cccd86fdfc7ca81aa64.png

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)
../_images/79f42d208312dce9ca2c2c378345dc56d9c84760358a0cdf9ccdb2587e3f95ff.png

\(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'))
../_images/5a4e691867a1c904294e3c6323ada74b761e9af1bd1989f647a4faa36dd59810.png ../_images/51d24dd1520fd54a9a9206fa25cc6548d34b16d2f07b65202b84355403085c46.png ../_images/0ecacdc37beb673032d12ec9b34351060ceb44e422e865ef52aa21fe4793d206.png ../_images/ff708393f7b9732c42db31a3d8d8b5e40450d28b9c98708647eb874a1be50daa.png ../_images/68cfcfbb011b7b329231a50a7a92cb65b812c6fa69e4e794c2df469ba3d8298c.png