An introduction to
thermonuclear reaction networks


 

Thanassis Psaltis (@psaltistha)
North Carolina State University & Triangle Universities Nuclear Laboratory
 
 

May 16, 2023
2nd Frontiers Summer School


 

Happy to see you! 😁


 

What do I do?

"I use big machines and computers to study how stars that blow up in space create the stuff we are made of."
Brought to you by: https://xkcd.com/simplewriter/

Why listen to me?

I have used many reaction network codes for thousands of hours.

How did I end up here?

🇬🇷→🇨🇦 → 🇩🇪 → 🇺🇸 → 🤷‍♂


 
 

What are my hobbies?

♟ 📚 🏃‍♂ 🎮 ️ 🍪

From a scale 1-5 my experience with thermonuclear reaction networks is...

🫠

Topics we will cover

 

Thermonuclear reaction rates

 

Forming the network

 

Solving the network

 

Sensitivity studies using reaction networks

Lecture notes and slides available at

http://psaltisa.github.io/teaching/reaction-networks

Thermouclear reaction networks


 

are an essential tool for nuclear astrophysics

W. R. Hix and F.-K. Thielemann, J. Comput. Appl. Math 109, 321 (1999)
C. Jiang et al., New J. Phys. 23, 083035 (2021)

A thermonuclear reaction network is a system of first-order coupled differential equations, which describes the abundance evolution and energy generation, and depends only on the nuclear reaction rates and the local thermodynamical conditions.

Thermonuclear reaction networks are
an essential tool for you!

Understand how stars
live and die




Photo credit: NASA/Solar Dynamics Observatory (SDO)

Understand the origin of the elements

C. Kobayashi, A. Karakas and M. Lugaro, Astrophys. J 900, 179 (2020)

Thermonuclear reaction rates

 

C. Iliadis, "Nuclear Physics of Stars", 2nd edition (2017)
R. Longland et al., Nucl. Phys. A 841, 1 (2010)

Why does the Sun shine?

Four protons are converted to one $\alpha$-particle
$\sim 27\ \mathrm{MeV}$ of energy is released

Thermonuclear Reaction Rate

 

Number of reactions per unit volume and time

\begin{aligned} r_{ab} &= \underbrace{n_a n_b}_{\text{number densities}} \langle \sigma v \rangle_{ab} \\ &= \color{#CC0000}{n_a n_b} \int \sigma(v)~P(v)~v~dv \end{aligned}

Reactions are initiated from their thermal motion → thermonuclear

Number density, abundance and mass fraction

 

\[ n_i = \frac{\rho N_A X_i}{A_i} = \rho N_A Y_i \]

$X_i$, mass fraction • $\sum_i X_i =1$

$Y_i$, abundance • $\sum_i A_i Y_i =1$

Is mass conserved during nucleosynthesis? 🤔

Thermonuclear Reaction Rate

 

Number of reactions per unit volume and time

\begin{aligned} r_{ab} &= \underbrace{n_a n_b}_{\text{number densities}} \langle \sigma v \rangle_{ab} \\ &= n_a n_b \int \sigma(v)~\color{#CC0000}{P(v)~v~dv} \end{aligned}

Reactions are initiated from their thermal motion → thermonuclear

Velocity (energy) distribution in a stellar gas

Maxwell-Boltzmann: $ \phi(E) \propto e^{-E/kT} \sqrt{E} $

Thermonuclear reaction rate


 

\[ \langle \sigma v \rangle = \left( \frac{8}{\pi \mu}\right)^{1/2} (kT)^{-3/2} \int_0^\infty \color{#CC0000}{\sigma(E)}~E~e^{-E/kT} dE \]

We only need an expression of $\sigma(E)$ to evaluate it.

Non-resonant reactions • resonant reactions (tomorrow) • R-matrix theory (R.J. deBoer)
Hauser-Feshbach model (G. Perdikakis) • reaction mechanisms (G. Potel)...

Forming the network

 

D. Arnett, "Supernovae and nucleosynthesis" (1996)
W.R. Hix and F. Thielemann, J. Comput. Appl. Math. 109, 321 (1999)
W.R. Hix and B.S. Meyer, Nucl. Phys. A 777, 188 (2006)

How do abundances evolve?


Decays

\[ a \rightarrow b \] \[ \frac{dY_a}{dt} = \color{#CC0000}{\mathbf{-}} \lambda_a Y_a \] \[ \frac{dY_b}{dt} = \color{#CC0000}{\mathbf{+}} \lambda_a Y_a \]

Capture reactions

\[ a + b \rightarrow c+ \gamma \] \[ \frac{dY_{a,b}}{dt} = \color{#CC0000}{\mathbf{-}}\rho N_A Y_a Y_b \langle \sigma v \rangle_{ab} \] \[ \frac{dY_c}{dt} = \color{#CC0000}{\mathbf{+}} \rho N_A Y_a Y_b \langle \sigma v \rangle_{ab} \]

The rate equations: protons


destroyed by $p(p,e^+\nu)d$ and $d(p,\gamma)\mathrm{^{3}He}$

produced by the $\mathrm{^3He(^3He, 2p)^4He}$

$$ \frac{dY_p}{dt} =\underbrace{ - \rho N_A Y_p^2 \langle \sigma v \rangle_{pp} - \rho N_A Y_d Y_p \langle \sigma v \rangle_{dp}}_{\text{destruction}} + \underbrace{ \rho N_A Y_{\mathrm{^3He}}^2 \langle \sigma v \rangle_{\mathrm{^3He}\mathrm{^3He}}}_{\text{production}} $$

Can you write the rate equations for $d$, $\mathrm{^3He}$ and $\alpha$?

Remember that for two-body reactions $\frac{dY_i}{dt} = \pm \rho N_A Y_i Y_j \langle \sigma v \rangle_{ij}$
The three reactions: $p(p,e^+\nu)d$, $d(p,\gamma)\mathrm{^{3}He}$ & $\mathrm{^3He(^3He, 2p)^4He}$

The rate equations: pp-chain

$$ \frac{dY_p}{dt} = - \rho N_A Y_p^2 \langle \sigma v \rangle_{pp} - \rho N_A Y_d Y_p \langle \sigma v \rangle_{dp} + \rho N_A Y_{\mathrm{^3He}}^2 \langle \sigma v \rangle_{\mathrm{^3He}\mathrm{^3He}} $$
$$ \frac{dY_d}{dt} = \rho N_A \frac{Y_p^2}{2} \langle \sigma v \rangle_{pp} - \rho N_A Y_d Y_p \langle \sigma v \rangle_{dp} $$
$$ \frac{dY_{\mathrm{^3He}}}{dt} = \rho N_A Y_d Y_p \langle \sigma v \rangle_{dp} - \rho N_A Y_{\mathrm{^3He}}^2 \langle \sigma v \rangle_{\mathrm{^3He}\mathrm{^3He}} $$
$$ \frac{dY_\alpha}{dt} = \rho N_A \frac{Y_{\mathrm{^3He}}^2}{2} \langle \sigma v \rangle_{\mathrm{^3He ^3He}} $$

The rate equations

$$ \begin{aligned} \frac{dY_i}{dt} &= \underbrace{\sum_j N_i \lambda_j Y_j}_{\text{decays}} \\ &+\underbrace{\sum_{j,k} \frac{N_i}{N_j! N_k!} \rho N_A \langle \sigma v \rangle_{jk} Y_j Y_k}_{\text{2-body reactions}} \\ &+ \underbrace{\sum_{j,k,l} \frac{N_i}{N_j! N_k! N_l!} (\rho N_A)^2 \langle \sigma v \rangle_{jkl} Y_j Y_k Y_l}_{\text{3-body reactions}} \end{aligned} $$

The rate equations

More generally:

\[ \dot{\mathbf{y}} = f(\mathbf{y,t}) \]

$\mathbf{y}$ is the vector that contains all the abundances $Y_i$ at time $t$ - $\mathbf{y}=[Y_1, Y_2,\ldots Y_{i_{max}}]_t$

The network ingredients

Nuclear Physics

Reaction rates $N_A \langle \sigma v \rangle$ (Reaclib, STARLIB,...), reaction Q-values, half-lifes, fission fragment distributions, partition functions, ...

Astrophysics

Time evolution of the temperature, $T(t)$ and density $\rho(t)$, and sometimes the electron mole fraction ($Y_e$). The initial composition of the environment ($Y_i$). "Tracer particles" and parametric profiles.

For our example we can assume constant temperature and density (hydrostatic burning). $\mathrm{T= 15 \times 10^6~K}, \rho = \mathrm{100~g/cm^3}$.

Uncertainties in the nuclear and astrophysics input will affect the abundances!

Solving the network

 

F. Timmes, Astrophys. J., Suppl. Ser. 124, 241 (1999)
R. Longland, D. Martin and J. José, A&A 563, A67 (2014)

Our goal is to solve
$$ \dot{\mathbf{y}} = f(\mathbf{y, t}) $$

Integration methods: implicit vs explicit

\begin{equation} \mathbf{\dot{y}_n} = \frac{\mathbf{y_{n+1}} - \mathbf{y_n}}{h} = (1-\Theta) \mathbf{\dot{y}_{n+1}} + \Theta \mathbf{\dot{y}_n} \end{equation}

$\mathbf{y_n}, \mathbf{y_{n+1}}$ are the abundance vectors in time $t_n$ and $t_{n+1}$
$h$ is the integration time-step.

It all depends on the $\Theta$...

Integration methods: implicit vs explicit


$$ \frac{\mathbf{y_{n+1}} - \mathbf{y_n}}{h} = \begin{cases} \mathbf{\dot{y}_n} & \Theta=1,~\text{explicit forward Euler} \\ \frac{1}{2}[\mathbf{\dot{y}_{n+1}} + \mathbf{\dot{y}_n}] & \Theta= 1/2,~\text{semi-implicit trapezoidal} \\ \mathbf{\dot{y}_{n+1}} & \Theta=0,~\text{implicit backwards} \end{cases} $$

Explicit methods 😭

Implicit/semi-implicit methods 😁

Implicit/semi-implicit integration methods


"Wagoner's method" (implicit)

R.V. Wagoner, Astrophys. J Suppl. Ser. 18, 237 (1969)

Bader-Deuflehard method (semi-implicit)

G. Bader & P. Deuflehard, Numer. Math. 41, 373 (1983)

Gear's method (implicit)

C.W. Gear, Commun. ACM 14, 176 (1971)

R. Longland et al., A&A 563, A67 (2014) compares the performance of these methods for thermonuclear reaction networks.

The solution: enter the Jacobian

$$ \frac{\mathbf{y_{n+1}} - \mathbf{y_n}}{h} = \mathbf{\dot{y}_{n+1}} $$ $$ \mathbf{y_{n+1}} = \mathbf{y_n} + h f(t_{n+1}, \mathbf{y}, \ldots) $$

$$ [\mathbf{I}- h \mathbf{J}]\mathbf{y_{n+1}} = \mathbf{y_n} $$ $\mathbf{I}$ is the unitary matrix
$\mathbf{J} \equiv \partial f/\partial \mathbf{y}$ is the Jacobian matrix

The Jacobian matrix

$$ \mathbf{J}(Y_1, \ldots, Y_{i_{max}})= \begin{vmatrix} \frac{\partial \dot{Y_1}}{\partial Y_1} & \ldots & \frac{\partial \dot{Y_1}}{\partial Y_{i_{max}}} \\ \vdots & \ddots & \vdots \\ \frac{\partial \dot{Y_{i_{max}}}}{\partial Y_1} & \ldots & \frac{\partial \dot{Y_{i_{max}}}}{\partial Y_{i_{max}}}\end{vmatrix} $$

Can you write the Jacobian matrix for the pp-chain?

You can use the following substitution: $R_{ab}=\rho N_A \langle \sigma v \rangle_{ab}$
The four nuclei: $p,d,\mathrm{^3He},\alpha$
The three reactions: $p(p,e^+\nu)d$, $d(p,\gamma)\mathrm{^{3}He}$ & $\mathrm{^3He(^3He, 2p)^4He}$

The Jacobian matrix: pp-chain


$$ \mathbf{J}= \begin{vmatrix} -2Y_p R_{pp} - Y_d R_{pd} & - Y_p R_{dp} & 2 Y_{\mathrm{^3He}} R_{\mathrm{^3He^3He}} & 0 \\ Y_p R_{pp} - Y_d R_{dp} & -Y_p R_{dp} & 0 & 0 \\ Y_d R_{dp} & Y_p R_{dp} & -2Y_{\mathrm{^3He}}R_{\mathrm{^3He^3He}} & 0 \\ 0 & 0 & Y_{\mathrm{^3He}} R_{\mathrm{^3He ^3He}} & 0 \end{vmatrix} $$

The Jacobian matrix

White space indicates zero elements of the Jacobian! (sparse)

Source: F. Timmes' website

The evolution of the pp-chain abundances

Do you notice anything interesting in the ratio of $X_p/X_d$? 🤔

Equilibrium

$$ \frac{dY_d}{dt} = 0 $$
$$ \rho N_A \frac{Y_p^2}{2} \langle \sigma v \rangle_{pp} - \rho N_A Y_d Y_p \langle \sigma v \rangle_{dp} = 0 $$
$$ \frac{Y_d}{Y_p} = \frac{\langle \sigma v \rangle_{pp}}{2\langle \sigma v \rangle_{dp}} $$

For our conditions, the equilibrium is established in around one minute!

Nuclear Statistical Equilibrium (NSE)

Occurs in high temperatures ($\mathrm{T > 3 \times 10^9 K}$), where all reactions are in equilibrium via strong and electromagnetic interactions.

The abundances are determined by the properties of the gas and not by reactions rates!

Energy Generation

\begin{equation} \epsilon= - \sum_i N_A m_i c^2 \frac{dY_i}{dt}-\frac{d\epsilon_\nu}{dt}~\text{(MeV/g s}) \end{equation}

$m_i c^2$ is the rest mass energy of species $i$ in MeV
$d\epsilon_\nu/dt$ are the neutrino losses

More examples using WebNucleo during the hands-on session!

Sensitivity studies using reaction networks

 

C. Iliadis et al., Astrophys. J Suppl. Ser. 142, 105 (2002)
T. Rauscher et al., Mon. Not. Roy. Astron. Soc. 463, 4153 (2016)
and many more!

How do we identify reactions
whose rate uncertainties affect
nucleosynthesis processes?

Sensitivity studies motivate experiments

Same model, but different nuclear input!

S. Wanajo et al., Astrophys. J 729, 46 (2011)

Monte Carlo Studies

Same model, but change all $(\alpha,n)$ rates

A. Psaltis et al., Astrophys. J 935, 27 (2022)

Take-home message

 

Nuclear reaction networks are an indispensable tool for a nuclear astrophysicist. Learn how it works and use it to your advantage!

Questions

Lecture notes and slides available at

http://psaltisa.github.io/teaching/reaction-networks