Homework 5 solutions

1. approximate rates

We want to look at the approximation of the rate sequence \((\alpha,p)(p,\gamma)\) and \((\alpha,\gamma)\) into a single effective rate.

Our system is:

\[\begin{equation*} A(\alpha,\gamma)C \end{equation*}\]

\[\begin{equation*} A(\alpha,p)B(p,\gamma)C \end{equation*}\]

which evolves according to:

\[\begin{align*} \frac{dY_A}{dt} &= - \rho Y_A Y_\alpha \lambda_{\alpha,\gamma} + Y_C \lambda_{\gamma,\alpha} - \rho Y_A Y_\alpha \lambda_{\alpha,p} + \rho Y_B Y_p \lambda_{p,\alpha} \\ % \frac{dY_C}{dt} &= + \rho Y_A Y_\alpha \lambda_{\alpha,\gamma} - Y_C \lambda_{\gamma,\alpha} + \rho Y_B Y_p \lambda_{p,\alpha} - Y_C\lambda_{\gamma,p} \\ % \frac{dY_\alpha}{dt} &= - \rho Y_A Y_\alpha \lambda_{\alpha,\gamma} + Y_C \lambda_{\gamma,\alpha} - \rho Y_A Y_\alpha \lambda_{\alpha,p} + \rho Y_B Y_p \lambda_{p,\alpha} \\ % \frac{dY_B}{dt} = \frac{dY_p}{dt} &= + \rho Y_A Y_\alpha \lambda_{\alpha,p} - \rho Y_B Y_p \lambda_{p,\gamma} + Y_C \lambda_{\gamma,p} - \rho Y_B Y_p \lambda_{p,\alpha} \end{align*}\]

where the \(Y\) are the molar fractions of the species and \(\lambda = N_A \langle \sigma v \rangle\)

a.

If we assume that the protons are in equilibrium, then \(dY_p/dt = 0\), giving

\[\rho Y_B Y_p = \frac{\rho Y_A Y_\alpha \lambda_{\alpha,p} + Y_B \lambda_{\gamma,p}} {\lambda_{p,\gamma} + \lambda_{p,\alpha}}\]

If we substitute this into the \(dY_A/dt\) equation, we find:

\[\begin{align*} \frac{dY_A}{dt} = &-\rho Y_A Y_\alpha \left [ \lambda_{\alpha,\gamma} + \lambda_{\alpha,p} \left ( 1 - \frac{\lambda_{p,\alpha}}{\lambda_{p,\gamma} + \lambda_{p,\alpha}} \right ) \right ] \\ &+ Y_B \left [ \lambda_{\gamma,\alpha} + \frac{\lambda_{p,\alpha} \lambda_{\gamma,p}}{\lambda_{p,\gamma} + \lambda_{p,\alpha}} \right ] \end{align*}\]

and we can identify the effective forward rate as (after some simplification) the term in the first set of brackets:

\[\lambda^\prime_{\alpha,\gamma} = \lambda_{\alpha,\gamma} + \frac{\lambda_{\alpha,p}\lambda_{p,\gamma}}{\lambda_{p,\gamma} + \lambda_{p,\alpha}}\]

and the effective reverse rate as the term in the second set of brackets:

\[\lambda^\prime_{\gamma,\alpha} = \lambda_{\gamma,\alpha} + \frac{\lambda_{p,\alpha} \lambda_{\gamma,p}}{\lambda_{p,\gamma} + \lambda_{p,\alpha}} \]

b.

In the limit where the forward rates are much faster than the reverse rates, our approximations reduce to:

\[\lambda^\prime_{\alpha,\gamma} = \lambda_{\alpha,\gamma} + \lambda_{\alpha,p}\]

\[\lambda^\prime_{\gamma,\alpha} \approx 0\]

The forward rate is just the sum of the \(\alpha\)-captures, which makes sense, since this approximation assumed that the proton capture is effectively instantaneous.

2. \(\alpha\) network

We want to evolve the system:

\[\begin{align*} \frac{dn_4}{dt} &= -\frac{3n_4^3}{6} \langle 3\alpha \rangle - n_4 n_{12} \langle \alpha 12 \rangle - n_4 n_{16} \langle \alpha 16 \rangle \\ \frac{dn_{12}}{dt} &= \frac{n_4^3}{6} \langle 3\alpha \rangle - n_4 n_{12} \langle \alpha 12 \rangle \\ \frac{dn_{16}}{dt} &= n_4 n_{12} \langle \alpha 12 \rangle - n_4 n_{16} \langle \alpha 16 \rangle \\ \frac{dn_{20}}{dt} &= n_4 n_{16} \langle \alpha 16 \rangle \end{align*}\]

The rates (from CF88) are:

• \(\langle 3\alpha \rangle = \langle \sigma v \rangle\) for the 3-\(\alpha\) reaction. We can find this expressed as \(N_A^2 \langle \sigma v \rangle_{\alpha\alpha\alpha}\):

  \[\begin{equation*} N_A^2 \langle \sigma v \rangle_{\alpha\alpha\alpha} = \frac{2.79\times 10^{-8}}{T_9^3} \exp(-4.4027/T_9) \end{equation*}\]

• \(\langle \alpha 12 \rangle = \langle \sigma v \rangle\) for the \({}^{12}\mathrm{C}(\alpha,\gamma){}^{16}\mathrm{O}\) reaction. We find this expressed as \(N_A \langle \sigma v \rangle_{C\alpha}\):

  \[\begin{align*} N_A \langle \sigma v \rangle_{C\alpha} = & \frac{1.04\times 10^8}{T_9^2 (1.0 + 0.0489/T_9^{2/3})^2} \exp(-32.120/T_9^{1/3} - (T_9/3.496)^2) + \nonumber \\ & \frac{1.76\times 10^8}{T_9^2 (1.0+0.2654/T_9^{2/3})^2} \exp(-32.120/T_9^{1/3}) + \nonumber \\ &\frac{1.25\times 10^3}{T_9^{3/2}} \exp(-27.499/T_9) + 1.43\times 10^{-2} T_9^5 \exp(-15.541/T_9) \end{align*}\]

• \(\langle \alpha 16 \rangle = \langle \sigma v \rangle\) for the \({}^{16}\mathrm{O}(\alpha,\gamma){}^{20}\mathrm{Ne}\) reaction. We find this expressed as \(N_A \langle \sigma v \rangle_{O\alpha}\):

\[\begin{align*} N_A \langle \sigma v \rangle_{O\alpha} = & \frac{9.37\times 10^9}{T_9^{2/3}} \exp(-39.757/T_9^{1/3}-(T_9/1.586)^2) + \nonumber \\ & \frac{6.21\times 10^1}{T_9^{3/2}} \exp(-10.297/T_9) + \nonumber \\ & \frac{5.38\times 10^2}{T_9^{3/2}} \exp(-12.226/T_9) + \nonumber \\ & 1.30\times 10^1 T_9^2 \exp(-20.093/T_9) \end{align*}\]

a. Number density to molar fractions

As written our system is expressed in terms of number density. Recall:

\[n_k = \frac{\rho X_k}{m_u A_k} = \frac{\rho Y_k}{m_u}\]

or using \(N_A = 1/m_u\), we have:

\[n_k = \frac{\rho N_A X_k}{A_k} = \rho N_A Y_k\]

We’ll work in terms of mass fractions, which have the property that

\[\sum_k X_k = 1\]

For example, to rewrite the evolution equation for \({}^{4}\mathrm{He}\), we could have:

\[\begin{align*} \frac{N_A\rho}{A_4} \frac{dX_4}{dt} = &-\frac{1}{2} \left ( \frac{N_A \rho X_4}{A_4} \right )^3 \langle 3\alpha \rangle - \left (\frac{N_A\rho X_4}{A_4} \right ) \left ( \frac{N_A\rho X_{12}}{A_{12}} \right ) \langle \alpha 12 \rangle \\ &- \left ( \frac{N_A\rho X_4}{A_4} \right ) \left (\frac{N_A\rho X_{16}}{A_{16}} \right ) \langle \alpha 16 \rangle \end{align*}\]

or simplifying:

\[\frac{dY_4}{dt} = -\frac{1}{2} \rho^2 Y_4^3 N_A^2 \langle 3\alpha \rangle - \rho Y_4 Y_{12} N_A \langle \alpha 12 \rangle - \rho Y_4 Y_{16} N_A \langle \alpha 16 \rangle\]

Note how \(N_A\) enters with the rate is consistent with the expressions we were given above.

The system in terms of molar fractions becomes:

\[\begin{align*} \frac{dY_4}{dt} &= -\frac{1}{2} \rho^2 Y_4^3 N_A^2 \langle 3\alpha \rangle - \rho Y_4 Y_{12} N_A \langle \alpha 12 \rangle - \rho Y_4 Y_{16} N_A \langle \alpha 16 \rangle \\ \frac{dY_{12}}{dt} &= \frac{1}{6} \rho^2 Y_4^3 N_A^2 \langle 3\alpha \rangle - \rho Y_4 Y_{12} N_A \langle \alpha 12 \rangle \\ \frac{dY_{16}}{dt} &= \rho Y_4 Y_{12} N_A \langle \alpha 12 \rangle - \rho Y_4 Y_{16} N_A \langle \alpha 16 \rangle \\ \frac{dY_{20}}{dt} &= \rho Y_4 Y_{16} N_A \langle \alpha 16 \rangle \end{align*}\]

b. Integrating

Let’s start by writing our RHS function. We’ll rely on the ODE solver to estimate the Jacobian.

import numpy as np
import matplotlib.pyplot as plt

def rhs(t, Y, rho, T, rc12ago16_factor=1.0):

    Y4  = Y[0]
    Y12 = Y[1]
    Y16 = Y[2]
    Y20 = Y[3]

    # temperature factors
    T9 = T/1.e9
    T913 = T9**(1./3.)
    T923 = T913**2
    T932 = T9**1.5

    # N_A**2 <sigma v> for 3-alpha from CF88 for T9 > 0.08
    r3a = 2.79e-8 / T9**3 * np.exp(-4.4027 / T9)

    # N_A <sigma v> for C12(a,g)O16
    rc12ago16 = (1.04e8 / T9**2 / (1.0 + 0.0489 / T923)**2 *
                 np.exp(-32.120 / T913 - (T9/3.496)**2) + 1.76e8 / T9**2 / (1.0 + 0.2654 / T923)**2 *
                 np.exp(-32.120 / T913) + 1.25e3 / T932 * np.exp(-27.499 / T9) +
                 1.43e-2 * T9**5 * np.exp(-15.541 / T9))
    rc12ago16 *= rc12ago16_factor
    
    # N_A <sigma v> for O16(a,g)Ne20
    ro16agne20 = (9.37e9 / T923 * np.exp(-39.757 / T913 - (T9 / 1.586)**2) + 
                  6.21e1 / T932 * np.exp(-10.297 / T9) + 5.38e2 / T932 * np.exp(-12.226 / T9) +
                  1.30e1 * T9**2 * np.exp(-20.093 / T9))

    dY4dt = -(1./2.) * rho**2 * Y4**3 * r3a - rho * Y4 * Y12 * rc12ago16 - rho * Y4 * Y16 * ro16agne20
    dY12dt = (1./6.) * rho**2 * Y4**3 * r3a - rho * Y4 * Y12 * rc12ago16
    dY16dt = rho * Y4 * Y12 * rc12ago16 - rho * Y4 * Y16 * ro16agne20
    dY20dt = rho * Y4 * Y16 * ro16agne20

    return np.array([dY4dt, dY12dt, dY16dt, dY20dt])

We can already investigate timescales for the burning—consider starting off with pure He, the timescale to consume the helium is:

\[\tau_\mathrm{burn} \sim \frac{Y_0({}^{4}\mathrm{He})}{|dY({}^{4}\mathrm{He})/dt|}\]

Let’s plot this for \(\rho = 10^5~\mathrm{g~cm^{-3}}\) as a function of temperature.

Where we note that \(Y_0({}^4\mathrm{He}) = 1/4\), since the atomic weight is 4.

fig, ax = plt.subplots()

rho = 1.e5
Y = np.array([0.25, 0.0, 0.0, 0.0])
T = np.logspace(8., 9, 100)

tau = []
for TT in T:
    dYdt = rhs(0.0, Y, rho, TT)
    tau.append(Y[0] / np.abs(dYdt[0]))

ax.semilogy(T, tau)
ax.set_xlabel("T [K]")
ax.set_ylabel(r"$\tau$")
ax.grid(ls=":")

Note how strong the temperature dependence is. At the low end, \(T_9 = 8\), the timescale to burn is \(> 10^{14}\) s, while at the high end, it is \(\sim 10^2\) s.

Note that this is only an estimate—once we make C and O, more He will be consumed, further shortening this timescale, but those nuclei will have longer timescales to be consumed.

Now we can integrate the system. We’ll use the SciPy solve_ivp function. Note that this network is stiff, so we will use the "BDF" method.

from scipy.integrate import solve_ivp

def integrate(t_end, Y0, rho, T, rc12ago16_factor=1.0):

    sol = solve_ivp(rhs, [0, t_end], Y0, method="BDF",
                    dense_output=True, args=(rho, T, rc12ago16_factor), rtol=1.e-8, atol=1.e-10)

    return sol.t, sol.y

c. Plotting

Let’s write a simple function to plot the species

def plot(t_out, X_out):
    fig, ax = plt.subplots()

    ax.plot(t_out, X_out[0,:], label="He4")
    ax.plot(t_out, X_out[1,:], label="C12")
    ax.plot(t_out, X_out[2,:], label="O16")
    ax.plot(t_out, X_out[3,:], label="Ne20")

    ax.legend(frameon=False, loc="best")
    ax.set_ylim(1.e-10, 2)
    ax.grid()
    ax.set_xscale("log")
    ax.set_yscale("log")
    
    return fig

We want to plot mass fractions, so let’s make an array of atomic weights

A = np.array([4.0, 12.0, 16.0, 20.0])

Y0 = np.array([0.25, 1.e-10, 1.e-10, 1.e-10])

rho = 1.e5
T = 3e8

tmax = 1.e12

t_out, Y_out = integrate(tmax, Y0, rho, T)

# convert to mass fractions
X_out = np.zeros_like(Y_out)
X_out[0, :] = Y_out[0, :] * A[0]
X_out[1, :] = Y_out[1, :] * A[1]
X_out[2, :] = Y_out[2, :] * A[2]
X_out[3, :] = Y_out[3, :] * A[3]

fig = plot(t_out, X_out)

Notice that as He is consumed, initially we build up C, but once there is enough C, we start to make O and then a little Ne.

What happens at higher T?

Y0 = np.array([0.25, 1.e-10, 1.e-10, 1.e-10])

rho = 1.e5
T = 6e8

tmax = 1.e12

t_out, Y_out = integrate(tmax, Y0, rho, T)

# convert to mass fractions
X_out = np.zeros_like(Y_out)
X_out[0, :] = Y_out[0, :] * A[0]
X_out[1, :] = Y_out[1, :] * A[1]
X_out[2, :] = Y_out[2, :] * A[2]
X_out[3, :] = Y_out[3, :] * A[3]

fig = plot(t_out, X_out)

Now we’ve made a lot of \({}^{20}\mathrm{Ne}\)! and the timescale are overall much shorter.

But we note that at these temperatures were are starting to miss some other important rates, especially those involving \({}^{12}\mathrm{C} + {}^{12}\mathrm{C}\)

d. Rate sensitivity

Let’s explore what happens to the \({}^{12}\mathrm{C}(\alpha,\gamma){}^{16}\mathrm{O}\) rate, but multiplying it by 1.7 as suggested in Weaver & Woosley 1993.

Y0 = np.array([0.25, 1.e-10, 1.e-10, 1.e-10])

rho = 1.e5
T = 3e8

tmax = 1.e12

t_out, Y_out = integrate(tmax, Y0, rho, T, rc12ago16_factor=1.7)

# convert to mass fractions
X_out = np.zeros_like(Y_out)
X_out[0, :] = Y_out[0, :] * A[0]
X_out[1, :] = Y_out[1, :] * A[1]
X_out[2, :] = Y_out[2, :] * A[2]
X_out[3, :] = Y_out[3, :] * A[3]

t_orig, Y_orig = integrate(tmax, Y0, rho, T, rc12ago16_factor=1.0)

# convert to mass fractions
X_orig = np.zeros_like(Y_orig)
X_orig[0, :] = Y_orig[0, :] * A[0]
X_orig[1, :] = Y_orig[1, :] * A[1]
X_orig[2, :] = Y_orig[2, :] * A[2]
X_orig[3, :] = Y_orig[3, :] * A[3]

fig, ax = plt.subplots()

ax.plot(t_out, X_out[0, :], label="He", color="C0")
ax.plot(t_orig, X_orig[0, :], linestyle=":", color="C0")

ax.plot(t_out, X_out[1, :], label="C", color="C1")
ax.plot(t_orig, X_orig[1, :], linestyle=":", color="C1")

ax.plot(t_out, X_out[2, :], label="O", color="C2")
ax.plot(t_orig, X_orig[2, :], linestyle=":", color="C2")

ax.plot(t_out, X_out[3, :], label="Ne", color="C3")
ax.plot(t_orig, X_orig[3, :], linestyle=":", color="C3")

ax.set_yscale("log")
ax.set_xscale("log")

ax.set_ylim(1.e-10, 2)

(1e-10, 2)

Here we see the difference when we run with the enhanced rate (solid lines) – we make slightly more oxygen.

e. Energy release

We can get the energy release from the change in molar abundances over the burn.

We’ll do this as:

\[\epsilon = -N_A c^2 \sum_k \frac{\Delta Y_k}{\Delta t} m_k\]

where \(m_k\) is the mass of nucleus \(k\). This is effectively just accumulating the \(\Delta m c^2\) from the different reactions. The end result will be in units of energy / mass / time.

We can get the mass in terms of the mass excess as:

\[m_k = A_k m_u + \frac{\Delta M_k}{c^2}\]

Here are the mass excesses from AME2020

from scipy import constants

keV2erg = (constants.eV * constants.kilo) / constants.erg
m_u = constants.value("unified atomic mass unit") * constants.kilo
c = constants.speed_of_light / constants.centi
N_A = constants.Avogadro

dm = np.array([2424.91587,  # He4
               0.0,  # C12
               -4737.00217,  # O16
               -7041.9322    # Ne20
              ])

# let's get masses in grams
mass = A * m_u + dm * keV2erg / c**2
mass

array([6.64647907e-24, 1.99264688e-23, 2.65601806e-23, 3.31982279e-23])

we’ll redo our original burn

Y0 = np.array([0.25, 1.e-10, 1.e-10, 1.e-10])

rho = 1.e5
T = 3e8

tmax = 1.e12

t_out, Y_out = integrate(tmax, Y0, rho, T)

dY = np.zeros_like(Y0)
for i in range(len(dY)):
    dY[i] = Y_out[i, -1] - Y_out[i, 0]

eps = 0.0
for i in range(len(dY)):
    eps -= N_A * c**2 * mass[i] * dY[i]
    
eps /= tmax
eps

605872.4480471261

This shows that our energy generation rate was about \(6\times 10^5~\mathrm{erg/g/s}\)