Solving for nuclear statistical equilibrium

pynucastro is able to calculate the abundances at nuclear statistical equilibrium of a given reaction network.

The constrained equations are setup following Calder et al. [2007] and Seitenzahl et al. [2009] and solved using scipy.optimize.fsolve().

Here we show how to find the NSE state of a set of nuclei.

import pynucastro as pyna
import numpy as np
import matplotlib.pyplot as plt

We start by creating a Library object that reads all the ReacLib rates and link different nuclei of choice.

Then initialize NSENetwork by using the the rates created by the Library object.

library = pyna.ReacLibLibrary()

subch_full = library.linking_nuclei(["p", "he4",
                                     "c12", "n13", "n14",
                                     "o16", "f18", "ne20",
                                     "ne21", "na22", "na23",
                                     "mg24", "al27", "si28",
                                     "p31", "s32", "cl35",
                                     "ar36", "k39", "ca40",
                                     "sc43", "ti44", "v47",
                                     "cr48", "mn51",
                                     "fe52","co55","ni56"])

nse = pyna.NSENetwork(libraries=subch_full)

get_comp_nse(rho,T,ye) is a method of NSEState that returns composition at NSE as an object Composition by providing density, temperature, and a presribed electron fraction.

There is pre-set initial guess for the proton and neutron chemical potentials, however one should adjust the initial guess, accordingly if no solutions are found or taking a long time as the algorithm is searching for the correct initial guess. To pass in the initial guess, set initial_guess=guess. guess is a user-supplied list of the form [mu_p, mu_n], where mu_p and mu_n are the chemical potential for proton and neutron, respectively.

It has options of incorporating coulomb correction terms by let use_coulomb_corr=True, and return the solution of proton and neutron chemical potentials by return_sol=True.

comp, sol = nse.get_comp_nse(1e7, 6e9, 0.50, use_coulomb_corr=True, return_sol=True)
fig = comp.plot()

We can also explicitly print the composition

for c, v in comp.X.items():
    print(f"{str(c):6} : {v:20.12g}")

p      :     0.00934989821059
He4    :       0.433640616227
C12    :    7.72859959965e-06
N13    :    2.59190269463e-10
N14    :    1.36631324171e-09
O16    :    1.56975241119e-05
F18    :    2.41138161836e-11
Ne20   :    3.01874050001e-07
Ne21   :    3.95638586959e-09
Na22   :    9.68462282544e-10
Na23   :    7.29695210583e-08
Mg24   :    3.47019416232e-05
Al27   :     3.3525821184e-05
Si28   :      0.0104391634336
P31    :     0.00213500726862
S32    :     0.00916060183422
Cl35   :     0.00400094977488
Ar36   :     0.00487944875964
K39    :     0.00555275469044
Ca40   :     0.00477002847089
Sc43   :      0.0003726716421
Ti44   :    0.000229271130459
V47    :     0.00490466077421
Cr48   :      0.0011465823045
Mn51   :      0.0689585743911
Fe52   :     0.00576589690229
Co55   :       0.415686733967
Ni56   :      0.0189151049137

print(f"After solving, the chemical potential for proton and neutron are {sol[0]:10.7f} and {sol[1]:10.7f}")

After solving, the chemical potential for proton and neutron are -7.7386289 and -9.9990481

For \(\rho = 10^7~\mathrm{g~cm^{-3}}\) and electron fraction, \(Y_e = 0.5\), we can explore how the NSE state changes with temperature.

Here we used a guess for chemical potential of proton and neutron of -6.0 and -11.5, respectively.

rho = 1e7
ye = 0.5
temps = np.linspace(2.5, 12.0, 100)
X_s = []

guess = [-6.0, -11.5]
for i, temp in enumerate(temps):
    nse_comp, sol = nse.get_comp_nse(rho, temp*1.0e9, ye, init_guess=guess,
                                     use_coulomb_corr=True, return_sol=True)
    guess = sol
    nse_X_s = [nse_comp.X[nuc] for nuc in nse_comp.X]
    X_s.append(nse_X_s)

X_s = np.array(X_s)
nuc_names = nse.get_nuclei()

low_limit = 1e-2
fig, ax = plt.subplots()

for k in range(len(nuc_names)):
    line, = ax.plot(temps, X_s[:,k])
    if (max(X_s[:,k]) > low_limit):
        line.set_label(str(nuc_names[k]))

ax.legend(loc = "upper right")
ax.set_xlabel(r'Temp $[\times 10^9 K]$')
ax.set_ylabel('Mass Fraction')
ax.set_yscale('log')
ax.set_ylim([low_limit,1.1])
ax.set_title(rf"$Y_e = {ye}$, $\rho = {rho:6.3g}$")

fig.set_size_inches((7, 6))

For \(\rho = 10^7~\mathrm{g~cm^{-3}}\) and a relatively high temperature of \(T = 6 \times 10^9~\mathrm{K}\), we can look at how the NSE state changes with varying \(Y_e\).

Note: there will be a lowest \(Y_e\) that is attainable by the network, depending on what species are represented. As long as protons are in the network, we can reach a \(Y_e\) of \(1\).

ye_low = min(nuc.Z/nuc.A for nuc in nse.unique_nuclei)

rho = 1e7
ye_s = np.linspace(ye_low, 0.65, 100)
temp = 6.0e9
X_s_1 = []

guess = [-6.0, -11.5]
for i, ye in enumerate(ye_s):
    nse_comp_1, sol = nse.get_comp_nse(rho, temp, ye, init_guess=guess,
                                       use_coulomb_corr=True, return_sol=True)
    guess = sol
    nse_X_s_1 = [nse_comp_1.X[nuc] for nuc in nse_comp_1.X]
    X_s_1.append(nse_X_s_1)

X_s_1 = np.array(X_s_1)
nuc_names = nse.get_nuclei()

fig, ax = plt.subplots()

for k in range(len(nuc_names)):
    line, = ax.plot(ye_s, X_s_1[:,k])
    if (max(X_s_1[:,k]) > 0.01):
        line.set_label(str(nuc_names[k]))

ax.legend(loc="best")
ax.set_xlabel('ye')
ax.set_ylabel('Mass Fraction')
ax.set_yscale('log')
ax.set_ylim([0.01,1])
ax.set_title(rf"$\rho = {rho:6.3g}$ and $T={temp:6.3g}$")

fig.set_size_inches((7, 6))