Approximating double n-capture

Approximating double n-capture#

We want to approximate the sequence A(n,γ)X(n,γ)B as A(nn,γ)B by combining the two successive neutron captures into a single rate, assuming equilibrium of nucleus X.

This type of approximation is commonly done in the “aprox” family of networks, for example, approximating the sequence: 52Fe(n,γ)53Fe(n,γ)54Fe.

import pynucastro as pyna

rl = pyna.ReacLibLibrary()

Unapproximated version#

We’ll create a simple network of iron isotopes linked by neutron capture.

lib = rl.linking_nuclei(["fe52", "fe53", "fe54", "n"])

rc = pyna.PythonNetwork(libraries=[lib])
rc.write_network("ncapture.py")

import ncapture

fig = rc.plot(rotated=True, curved_edges=True)
_images/a64968593abab72f06ff63c6175657851a01892b497468cb9f4d79d829eca266.png

Now we’ll integrate it, starting with a mix of 52Fe and n.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

rho = 1.e9
T = 3.e9

X0 = np.zeros(ncapture.nnuc)
X0[ncapture.jn] = 0.5
X0[ncapture.jfe52] = 0.5
Y0 = X0 / ncapture.A

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

fig = plt.figure()
ax = fig.add_subplot(111)

for i in range(ncapture.nnuc):
    ax.loglog(sol.t, sol.y[i,:] * ncapture.A[i], label=f"X({ncapture.names[i].capitalize()})")

ax.set_xlim(right=1.e-11)
ax.set_ylim(1.e-6, 1)

ax.legend(fontsize="small")
ax.set_xlabel("t (s)")
ax.set_ylabel("X")

Text(0, 0.5, 'X')
_images/ec5d78e1d15172af8e96c7292f5bf6b441c851cba76d45387bf983b10a15916f.png

We see that the long term evolution is that all the iron winds up as 54Fe.

Approximating the double n capture#

We can use the helper function create_double_neutron_capture() to create a pair (forward and reverse) of ApproximateRates that approximate out the intermedia nucleus:

rf, rr = pyna.rates.create_double_neutron_capture(rl, pyna.Nucleus("fe52"), pyna.Nucleus("fe54"))

rf

Fe52 + n + n ⟶ Fe54 + 𝛾

We can see what rates this carries internally:

rf.hidden_rates

[Fe52 + n ⟶ Fe53 + 𝛾, Fe53 + n ⟶ Fe54 + 𝛾, Fe54 ⟶ n + Fe53, Fe53 ⟶ n + Fe52]

Tip

Alternately, we could use the make_nn_g_approx() method on the PythonNetwork once we create it.

Now we’ll build a new network with just these approximate rates:

rc2 = pyna.PythonNetwork(rates=[rf, rr])
rc2.write_network("ncapture_approx.py")

import ncapture_approx

fig = rc2.plot(rotated=True, curved_edges=True)
_images/e11c6b733b410fabc89ebda5c6fb86c4811268a8e822ffb147317bd007dd8022.png 
rho = 1.e9
T = 3.e9

X0 = np.zeros(ncapture_approx.nnuc)
X0[ncapture_approx.jn] = 0.5
X0[ncapture_approx.jfe52] = 0.5
Y0 = X0 / ncapture_approx.A

tmax = 1.e-5
approx_sol = solve_ivp(ncapture_approx.rhs, [0, tmax], Y0, method="BDF", jac=ncapture_approx.jacobian,
                       dense_output=True, args=(rho, T), rtol=1.e-6, atol=1.e-8)

fig = plt.figure()
ax = fig.add_subplot(111)

for i in range(ncapture_approx.nnuc):
    ax.loglog(approx_sol.t, approx_sol.y[i,:] * ncapture_approx.A[i],
              label=f"X({ncapture_approx.names[i].capitalize()})")

ax.set_xlim(right=1.e-11)
ax.set_ylim(1.e-6, 1)

ax.legend(fontsize="small")
ax.set_xlabel("t (s)")
ax.set_ylabel("X")

Text(0, 0.5, 'X')
_images/534036212c126da4af84797ab837ef429f56cb8a0e4b10baa077166caa3f1f39.png

Comparison#

Finally, let’s plot the nuclei in common from both nets together

fig = plt.figure()
ax = fig.add_subplot(111)

for i in range(ncapture_approx.nnuc):
    ax.loglog(approx_sol.t, approx_sol.y[i,:] * ncapture_approx.A[i],
              linestyle=":", color=f"C{i}")

    idx = ncapture.names.index(ncapture_approx.names[i])
    ax.loglog(sol.t, sol.y[idx,:] * ncapture.A[idx],
              label=f"X({ncapture.names[idx].capitalize()})",
              linestyle="-", color=f"C{i}")

ax.legend()
ax.set_xlim(right=1.e-11)
ax.set_ylim(1.e-6, 1)

(1e-06, 1)
_images/166f9ae2fc7a678683278457ebfcae0f2c429ccb5f246df881038f9b74625cb4.png

Here the dotted line is the approximate version. We see that the 52Fe is consumed at the same rate. The approximate net produces 54Fe faster simply because there is no intermediate 53Fe in the network.

Contents