SkyNet: A modular nuclear reaction network library

TL;DR

SkyNet is a flexible, modular, open-source nuclear reaction network library designed for simulating nucleosynthesis in various astrophysical environments, supporting diverse reactions and physics with high accuracy.

Contribution

It introduces a comprehensive, adaptable nuclear reaction network framework with detailed physics implementations and extensive validation against existing codes.

Findings

SkyNet achieves agreement with existing codes within a few percent.

It supports a wide range of isotopes and reactions.

The code is modular and easily extendable.

Abstract

Almost all of the elements heavier than hydrogen that are present in our solar system were produced by nuclear burning processes either in the early universe or at some point in the life cycle of stars. In all of these environments, there are dozens to thousands of nuclear species that interact with each other to produce successively heavier elements. In this paper, we present SkyNet, a new general-purpose nuclear reaction network that evolves the abundances of nuclear species under the influence of nuclear reactions. SkyNet can be used to compute the nucleosynthesis evolution in all astrophysical scenarios where nucleosynthesis occurs. SkyNet is free and open-source and aims to be easy to use and flexible. Any list of isotopes can be evolved and SkyNet supports various different types of nuclear reactions. SkyNet is modular so that new or existing physics, like nuclear reactions or equations of state, can easily be added or modified. Here, we present in detail the physics implemented in SkyNet with a focus on a self-consistent transition to and from nuclear statistical equilibrium (NSE) to non-equilibrium nuclear burning, our implementation of electron screening, and coupling of the network to an equation of state. We also present comprehensive code tests and comparisons with existing nuclear reaction networks. We find that SkyNet agrees with published results and other codes to an accuracy of a few percent. Discrepancies, where they exist, can be traced to differences in the physics implementations.

Tables

Table 1. Table 1 Errors between the final mass-summed abundances between the different networks ( Equation   115 ).

Test case	Screening	Self-heating	S – SnoDB	S – W	SnoDB – W	S – X	SnoDB – X	X – W
Neutron-rich r-process	yes	yes/no∗	03.20	05.80	3.40	03.20	0.06600	3.40
(§6.2.1, Figure 8)	yes	no	33.00	35.00	4.50	32.00	1.80000	3.30
	no	yes/no∗	03.40	05.90	3.40	03.40	0.04200	3.40
	no	no	33.00	36.00	3.30	33.00	0.10000	3.20
Explosive X-ray burst	yes	yes/no∗	00.39	02.30	2.10	01.30	1.10000	1.00
(§6.2.2, Figure 9)	yes	no	13.00	14.00	4.00	11.00	3.50000	3.30
	no	yes/no∗	00.38	02.30	2.00	00.47	0.08600	2.00
	no	no	17.00	18.00	1.50	18.00	1.20000	0.40
Hydrostatic C/O burn	yes	yes	00.64	00.68	0.57	03.90	4.10000	3.90
(§6.2.3, Figure 10)	yes	no	00.10	00.22	0.18	00.16	0.22000	0.13
	no	yes	00.68	01.50	1.90	04.00	4.40000	2.50
	no	no	00.10	00.27	0.17	00.10	0.00036	0.17
Notes. S: SkyNet with detailed balance, SnoDB: SkyNet without detailed balance, W: WinNet, X: XNet. The error measures the average fractional difference between the final abundances. The numbers shown are in percent. The three networks generally agree very well with each other. The error between SkyNet without detailed balance (SnoDB), WinNet (W), and XNet (X) are usually of similar magnitude and on the few percent level. We also see that using detailed balance for the inverse rates has a big impact in the first two test cases, especially when self-heating is turned off. Since WinNet and XNet do not use detailed balance, the error is bigger when they are compared to SkyNet with detailed balance (S).
∗ Self-heating is turned on in SkyNet and WinNet, but not in XNet, because its self-heating method does not apply in these test cases. Instead, XNet uses the temperature computed by the self-heating SkyNet run without detailed balance (SnoDB).

Equations

$$\sum _{\alpha \in {\tilde{\mathcal{R}}}_n}[\alpha ]\rightleftharpoons \sum _{\beta \in {\tilde{\mathcal{P}}}_n}[\beta ],$$

$$\left({\partial }_t+\frac{\partial k_{ϵ}^0}{\partial {\mathbf{k}}_{ϵ}}\cdot \nabla \right)f_{ϵ}(x^{\mu },k_{ϵ}^{\mu })$$

$$=$$

$$+(1\pm f_{ϵ})\sum _n{\mathcal{N}}_n^{\text{inverse}}\left[\prod _{\alpha \in {\tilde{\mathcal{R}}}_n,\alpha \ne ϵ}{\int }_{[\alpha ]}(1\pm f_{\alpha })\right]\left[\prod _{\beta \in {\tilde{\mathcal{P}}}_n}{\int }_{[\beta ]}{f}_{\beta }\right]{\delta }^4\left(\sum _{\alpha \in {\tilde{\mathcal{R}}}_n}{k}_{\alpha }^{\mu }-\sum _{\beta \in {\tilde{\mathcal{P}}}_n}{k}_{\beta }^{\mu }\right)\text{mathcalligra}{r}_n\left(k_{\{\alpha \}}^{\mu },k_{\{\beta \}}^{\mu }\right),$$

$${\int }_{[ϵ]}=g_{ϵ}\int \frac{d^3{k}_{ϵ}}{(2\pi {)}^3}$$

$$\text{mathcalligra}r\left(k_{\{\alpha \}}^{\mu },k_{\{\beta \}}^{\mu }\right)=(2\pi {)}^4\frac{⟨|{\mathcal{T}}^2|⟩}{2{\prod }_{\alpha \in \tilde{\mathcal{R}}}2k_{\alpha }^0{\prod }_{\beta \in \tilde{\mathcal{P}}}2k_{\beta }^0},$$

$$\sum _{j\in {\mathcal{R}}_n}{N}_j^n[j]\rightleftharpoons \sum _{l\in {\mathcal{P}}_n}{N}_l^n[l],$$

$${\mathcal{N}}_n^{\text{forward}}=\prod _{j\in {\mathcal{R}}_n}\frac{1}{N_j^n!}\text{and}{\mathcal{N}}_n^{\text{inverse}}=\prod _{l\in {\mathcal{P}}_n}\frac{1}{N_l^n!}.$$

$${\lambda }_n=n_B^{-1}{\mathcal{N}}_n^{\text{forward}}\left[\prod _{j\in {\mathcal{R}}_n}{\left(\frac{n_B}{n_j}\right)}^{N_j^n}\right]\left[\prod _{\alpha \in {\tilde{\mathcal{R}}}_n}{\int }_{[\alpha ]}{f}_{\alpha }\right]\left[\prod _{\beta \in {\tilde{\mathcal{P}}}_n}{\int }_{[\beta ]}(1\pm f_{\beta })\right]{\delta }^4\left(\sum _{\alpha \in {\tilde{\mathcal{R}}}_n}{k}_{\alpha }^{\mu }-\sum _{\beta \in {\tilde{\mathcal{P}}}_n}{k}_{\beta }^{\mu }\right)\text{mathcalligra}{r}_n\left(k_{\{\alpha \}}^{\mu },k_{\{\beta \}}^{\mu }\right),$$

$$n_m={\int }_{[m]}{f}_m=g_m\int \frac{d^3{k}_m}{(2\pi {)}^3}{f}_m$$

$${\partial }_t{n}_{ϵ}+\nabla \cdot \left(⟨{\mathbf{v}}_{ϵ}⟩n_{ϵ}\right)=\sum _n\left[-{\lambda }_n{n}_B^{1-N_{\mathcal{R}}^n}\prod _{j\in {\mathcal{R}}_n}{n}_j^{N_j^n}+{\lambda }_n^{\prime }{n}_B^{1-N_{\mathcal{P}}^n}\prod _{l\in {\mathcal{P}}_n}{n}_l^{N_l^n}\right].$$

$$N_{\mathcal{R}}^n=\sum _{j\in {\mathcal{R}}_n}{N}_j^n\text{and}{N}_{\mathcal{P}}^n=\sum _{l\in {\mathcal{P}}_n}{N}_l^n.$$

$${\partial }_t{n}_i+\nabla \cdot \left(⟨{\mathbf{v}}_i⟩n_i\right)=\sum _{\alpha }{\lambda }_{\alpha }(-R_i^{\alpha }+P_i^{\alpha })N_i^{\alpha }{n}_B^{1-N_{\mathcal{R}}^{\alpha }}\prod _{m\in {\mathcal{R}}_{\alpha }}{n}_m^{N_m^{\alpha }},$$

$$\displaystyle R^{\alpha}_{i}=\left\{\begin{array}[]{@{}l@{\ \ }l@{}}1&\text{if $i\in\mathcal{R}_{\alpha}$}\\ 0&\text{otherwise}\end{array}\right.\qquad\text{and}\qquad P^{\alpha}_{i}=\left\{\begin{array}[]{@{}l@{\ \ }l@{}}1&\text{if $i\in\mathcal{P}_{\alpha}$}\\ 0&\text{otherwise.}\end{array}\right.$$

$$Y_i\equiv \frac{n_i}{n_B}=\frac{N_i/V}{N_B/V}=\frac{N_i}{N_B},$$

$$\frac{d{Y}_i}{dt}$$

$$=\frac{1}{n_B}\left[{\partial }_t{n}_i-\frac{n_i}{n_B}{\partial }_t{n}_B+\nabla \cdot (\mathbf{v}{n}_i)-n_i\nabla \cdot \mathbf{v}-\frac{n_i}{n_B}\left(\nabla \cdot (\mathbf{v}{n}_B)-n_B\nabla \cdot \mathbf{v}\right)\right]$$

$$=\frac{1}{n_B}\left[{\partial }_t{n}_i+\nabla \cdot (\mathbf{v}{n}_i)-\frac{n_i}{n_B}\left({\partial }_t{n}_B+\nabla \cdot (\mathbf{v}{n}_B)\right)\right]$$

$$=\frac{1}{n_B}\left[{\partial }_t{n}_i+\nabla \cdot (\mathbf{v}{n}_i)\right],$$

$$\frac{d{Y}_i}{dt}=\sum _{\alpha }{\lambda }_{\alpha }(-R_i^{\alpha }+P_i^{\alpha })N_i^{\alpha }\prod _{m\in {\mathcal{R}}_{\alpha }}{Y}_m^{N_m^{\alpha }},$$

$$\frac{d{Y}_{\text{ce}{}^{12}C}}{dt}$$

$$\frac{d{Y}_{\text{ce}{}^{4}He}}{dt}$$

$$\frac{d{Y}_{\text{ce}{}^{16}O}}{dt}$$

$$\frac{d{Y}_{\text{ce}{}^{12}C}}{dt}$$

$$\frac{d{Y}_{\text{ce}{}^{20}Ne}}{dt}$$

$$\frac{d{Y}_{\text{ce}{}^{4}He}}{dt}$$

$${\sigma }_{\alpha }(k_1^{\mu }-k_2^{\mu })=\frac{1}{v_{\text{rel}}}\left[\prod _{l\in {\mathcal{P}}_{\alpha }}{\int }_{[l]}\right]{\delta }^4\left(k_1^{\mu }+k_2^{\mu }-\sum _{l\in {\mathcal{P}}_{\alpha }}{k}_l^{\mu }\right)\text{mathcalligra}{r}_{\alpha }\left(k_1^{\mu },k_2^{\mu },k_{\{l\}}^{\mu }\right),$$

$${\sigma }_{\alpha }$$

$${\lambda }_{\alpha }$$

$$r_{i,j}=n_B^{-1}{\lambda }_{\alpha }{n}_1{n}_2=n_1{n}_2{\mathcal{N}}_{\alpha }⟨{\sigma }_{\alpha }{v}_{\text{rel}}⟩.$$

$$(Z,N)\rightleftharpoons Z[\text{p}]+N[\text{n}],$$

$${\mu }_i=Z_i{\mu }_{\text{p}}+N_i{\mu }_{\text{n}},$$

$${\lambda }_{\alpha }\prod _{j\in {\mathcal{R}}_{\alpha }}{Y}_{j,\text{eq}}^{N_j^{\alpha }}={\lambda }_{{\alpha }^{\prime }}\prod _{l\in {\mathcal{P}}_{\alpha }}{Y}_{l,\text{eq}}^{N_l^{\alpha }},$$

$$Y_i=\frac{G_i(T)}{n_B}{\left(\frac{m_{i}T}{2\pi }\right)}^{3/2}\exp [({\mu }_i-m_i)/T],$$

$$\frac{{\lambda }_{{\alpha }^{\prime }}}{{\lambda }_{\alpha }}=\frac{{\prod }_{j\in {\mathcal{R}}_{\alpha }}[G_j/n_B(m_{j}T/2\pi {)}^{3/2}{]}^{N_j^{\alpha }}}{{\prod }_{l\in {\mathcal{P}}_{\alpha }}[G_l/n_B(m_{l}T/2\pi {)}^{3/2}{]}^{N_l^{\alpha }}}\exp \left[\frac{1}{T}\sum _{j\in {\mathcal{R}}_{\alpha }}{N}_j^{\alpha }({\mu }_{j,\text{eq}}-m_j)-\frac{1}{T}\sum _{l\in {\mathcal{P}}_{\alpha }}{N}_l^{\alpha }({\mu }_{l,\text{eq}}-m_l)\right].$$

Full text

Contents

\twocolumngrid@starttoc

toc \onecolumngrid

1. Introduction

Nuclear and weak reactions play a crucial role in many astrophysical scenarios. Nuclear reactions typically occur at high temperatures and densities, because a large amount of energy is required to overcome the Coulomb repulsion between positively charged nuclei. Inside the cores of main-sequence stars, nuclear fusion converts hydrogen into helium, which releases nuclear binding energy as heat, keeping the star from collapsing (Bethe 1939). When massive stars undergo core collapse at the end of their lives, nuclear and weak reactions serve as important energy sinks and sources (neutrino cooling and heating). In core-collapse and Type Ia supernovae (SNe Ia), explosive nuclear burning mainly creates iron-group elements that are ejected into the interstellar medium (Nomoto et al. 1997; Woosley et al. 2002). According to our current understanding, a weak version of the rapid neutron capture process (r-process) can also occur (Wanajo 2013) in core-collapse supernovae. However, it appears that the full r-process that can synthesize all of the heavy elements predominantly happens during neutron star mergers (Freiburghaus et al. 1999b). Heavy elements up to bismuth can also be created in stars via the slow neutron capture process (s-process; Burbidge et al. 1957). Finally, when hydrogen and helium gas accretes onto a white dwarf, the accumulated material can undergo a thermonuclear explosion that creates a short-lived bright flash of light called a nova. If the accretor is a neutron star instead, the thermonuclear explosion results in an X-ray burst (José 2016, ch. 6) or a superburst (Strohmayer & Brown 2002).

To adequately understand these astrophysical objects and phenomena, one needs to account for the nuclear reactions that drive them. In some cases, reaction networks are mainly used to track the nuclear energy generation (e.g., Weaver et al. 1978; Müller 1986; Timmes et al. 2000). But in most cases, the evolution of the entire composition due to nuclear reactions is of interest. Because of the many ways that nuclei can react with each other to form other nuclides, the number of nuclear species that are relevant for many astrophysical processes ranges from dozens to thousands, and the number of nuclear reactions involved is hundreds to tens of thousands. For this reason, a large number of variables (i.e., nuclide abundances) that are coupled together by non-linear nuclear reaction rates (see §2) need to be evolved. Mathematically and computationally, the ensemble of coupled nuclear reactions is described by a nuclear reaction network.

Large-scale (several dozen species or more) nuclear reaction networks were first developed in the late 1960s and early 1970s (e.g., Truran et al. 1966, 1967; Arnett & Truran 1969; Woosley et al. 1973). These first networks were mainly for explosive nuclear burning in massive stellar evolution and supernovae, although earlier stellar evolution models also took nuclear reactions into account and evolved a handful of nuclear species (e.g., Hayashi et al. 1962; Hofmeister et al. 1964; Iben 1967). Early networks ranged in size from a few dozen to around a hundred species, with up to a few hundred reactions connecting the nuclei. Since then, nuclear reaction networks of different sizes have been used to study various astrophysical scenarios. Big bang nucleosynthesis calculations require the smallest networks with typically fewer than a dozen nuclear species, although some authors utilize much bigger networks of up to several dozen species (e.g., Wagoner 1973; Nollett & Burles 2000; Orlov et al. 2000; Coc et al. 2012; Cyburt et al. 2016). Networks with dozens to hundreds of species are employed in stellar evolution codes (e.g., Arnett 1977; Weaver et al. 1978; Maeder & Meynet 2000; Paxton et al. 2011; Bressan et al. 2012; Jones et al. 2015). Similar-sized or larger (up to hundreds of species) networks are also used to compute explosive nucleosynthesis in SNe Ia (e.g., Thielemann et al. 1986a; Iwamoto et al. 1999; Hillebrandt et al. 2013; Seitenzahl et al. 2013; Leung et al. 2015), core-collapse supernovae (e.g., Thielemann et al. 1996; Limongi & Chieffi 2003; Heger & Woosley 2010; Harris et al. 2014), novae (e.g., Weiss & Truran 1990; José & Hernanz 1998; Iliadis et al. 2002; Starrfield et al. 2016), and X-ray bursts (e.g., Schatz et al. 2001; Woosley et al. 2004; Cyburt et al. 2010; Parikh et al. 2013).

The largest nuclear networks are needed to simulate neutron capture processes. For the s-process in massive stars, it may be sufficient to use a few hundred to about a thousand nuclei (e.g., Prantzos et al. 1990; Käppeler et al. 2011; Nishimura et al. 2017). Larger nuclear reaction networks (typically thousands of isotopes) have been used for r-process nucleosynthesis calculations in neutrino-driven winds from core-collapse supernovae (e.g., Woosley & Hoffman 1992; Hoffman et al. 1997; Freiburghaus et al. 1999a; Arcones et al. 2010; Farouqi et al. 2010; Roberts et al. 2010; Wanajo 2013), in the jets of magnetorotational core-collapse supernovae (e.g., Winteler et al. 2012; Nishimura et al. 2015), in the dynamical ejecta of neutron star mergers (e.g., Freiburghaus et al. 1999b; Goriely et al. 2011; Bauswein et al. 2013; Wanajo et al. 2014; Eichler et al. 2015; Just et al. 2015; Fernández et al. 2017), in accretion disk ejecta following neutron star mergers (e.g., Surman et al. 2008; Perego et al. 2014; Martin et al. 2015; Lippuner et al. 2017), and in broader astrophysical contexts (e.g., Panov et al. 1995, 2001; Blinnikov & Panov 1996; Mumpower et al. 2012).

To evolve a nuclear reaction network, the rates of all reactions in the network are required. Most reaction rates, e.g. interactions between two or more nuclides, depend strongly on the energies of the incoming particles, due to Coulomb barrier penetration, resonances in the compound nuclear system, and other effects (e.g., Clayton 1968, §4). The rates of reactions only involving a single particle in the entrance channel, like $\beta$-decays and spontaneous fission, are constant.111Strictly speaking, $\beta$-decay rates are only constant in vacuum. In the medium, electron phase-space blocking introduces a dependency on the electron chemical potential and temperature (e.g., Arnett 1996, §11.3). Some reaction rates involving nuclides sufficiently close to the valley of stability can be measured experimentally as a function of energy, although in many cases astrophysical reactions occur at energies that are much lower than the experimentally accessible energy ranges (e.g., Rolfs & Rodney 1988, §4). Furthermore, most astrophysical processes involve unstable nuclei that may be very far away from stability and are not experimentally accessible for rate measurements. Therefore, theoretical models are necessary to compute the reaction rates needed by the reaction network (e.g., Cyburt et al. 2010). The Hauser–Feshbach approach, which assumes that the reactants form a single compound nucleus that subsequently decays into the reaction products, has been used extensively to compute nuclear reaction rates for astrophysics applications (e.g., Hauser & Feshbach 1952; Holmes et al. 1976; Woosley et al. 1978; Thielemann et al. 1986b; Rauscher & Thielemann 2000; Goriely et al. 2008).

Nuclear reaction networks also require the properties, such as masses and internal nuclear partition functions (e.g., Arcones & Martínez-Pinedo 2011; Brett et al. 2012; Mendoza-Temis et al. 2015; Mumpower et al. 2016), of all nuclides evolved in the network. These properties are needed to compute nuclear statistical equilibrium (NSE) and inverse reaction rates (see §2.3), as well as $\beta$-decay rates. Some of these nuclear properties, such as the masses, also enter the theoretical reaction rate calculations. Since many of the nuclides of interest are extremely unstable, special radioactive ion beam facilities are needed to produce and measure these exotic nuclei (see, e.g., Lunney et al. 2003; Schatz 2013; Mumpower et al. 2016, and references therein). Current radioactive beam facilities have made great progress in measuring unstable nuclei and new facilities or upgrades to current facilities are being built and planned. These new facilities will extend the reach to more exotic nuclei that are highly relevant to nuclear astrophysics scenarios (e.g., Schatz 2013, 2016; Mumpower et al. 2016). For the foreseeable future, however, it is necessary to use theoretical models to compute nuclear masses and $\beta$-decay properties for a large fraction of the nuclear species present in r-process networks (e.g., Lunney et al. 2003; Möller et al. 2003; Mumpower et al. 2016, and references therein).

Many authors who use nuclear reaction networks do not make the code of these networks publicly available. This makes it hard to reproduce and verify published results and also presents a barrier to new researchers joining the field since they first have to write their own nuclear reaction network. Notable exceptions of nuclear reaction networks that are publicly available are the various networks by Timmes (1999),222http://cococubed.asu.edu/code_pages/burn.shtml XNet by Hix & Thielemann (1999),333http://eagle.phys.utk.edu/xnet/trac and NucNet by Meyer & Adams (2007).444https://sourceforge.net/projects/nucnet-tools In this paper, we present a new nuclear reaction network called SkyNet that is publicly available as an open-source software at https://bitbucket.org/jlippuner/skynet (Lippuner & Roberts 2017b). This paper is based on version 1.0 of SkyNet (Lippuner & Roberts 2017a).

SkyNet was initially designed for evolving large reaction networks for r-process nucleosynthesis calculations, but thanks to its modularity and flexibility, SkyNet can easily be used for nucleosynthesis computations in many other astrophysical situations. Besides correctness, the main design goals behind SkyNet are usability and flexibility, making SkyNet an easy to use and versatile nuclear reaction network that is available for anyone to use. SkyNet can evolve an arbitrary set of nuclear species under various different types of nuclear reactions (§5.2). SkyNet can also compute NSE compositions (§2.3) and switch between evolving NSE and the full network in an automated and self-consistent way (§3.4). SkyNet contains electron screening corrections (§4) and an equation of state (EOS) that takes the entire composition into account (§A.2). For ease of use, SkyNet can be used from within Python (§5), and SkyNet can make movies of the nucleosynthesis evolution (see examples at http://stellarcollapse.org/lippunerroberts2015).

SkyNet has been used for r-process nucleosynthesis calculations in different scenarios by various authors: Lippuner & Roberts (2015), Radice et al. (2016), Roberts et al. (2017), Lippuner et al. (2017), Siegel & Metzger (2017), Vlasov et al. (2017), and Fernández et al. (2017). Here, we discuss the physics used in SkyNet, provide details on how it is implemented, and perform code tests and comparisons with other nuclear reaction networks.

This paper is organized as follows. In §2, we derive the equations that govern nuclear abundance evolution and equilibrium. §3 deals with the numerical implementation of the reaction network. We discuss in detail the electron screening corrections implemented in SkyNet in §4. In §5, we describe code implementation details. The code tests and comparisons are the subject of §6. We summarize in §7. In Appendix A, we briefly present the physics of an ideal Boltzmann gas and the EOS implemented in SkyNet. We show how SkyNet computes NSE in Appendix B, and in Appendix C we discuss neutrino interaction reactions.

Throughout this paper, we set the Boltzmann constant $k_{B}=1$ (i.e., all temperatures are measured in energy), the speed of light $c=1$, and the reduced Planck constant $\hbar=1$.

2. Nuclear reaction network basics

Astrophysical nuclear reaction networks track the composition of a system containing many species of nuclei, electrons, positrons, photons, and sometimes neutrinos. Essentially, they evolve the numbers of different nuclei in a system given a set of reactions and rates for those reactions that transmute nuclei into other nuclei. Although it is straightforward to heuristically write down a system of rate equations (Hix & Meyer 2006), it is useful to start from kinetic theory to tie the rate equations to the microscopic processes driving the nuclear transmutations.

2.1. Kinetic theory

Consider a homogeneous system of different species of particles (including nuclei, electrons, etc.) connected by a set of interactions, a subset of which changes particles of one type into another. A reaction indexed by $n$ converts a set of reactants into a set of products and vice versa. We write reaction $n$ as

[TABLE]

where $[\alpha]$ is a reactant, $[\beta]$ is a product, and $\mathcal{\tilde{R}}_{n}$ and $\mathcal{\tilde{P}}_{n}$ are the sets of all reactants and products, respectively. We emphasize that all particles are individually labeled, even if they are of the same species. For example, for the reaction $\smash{\ce{{}^{12}C}}+\smash{\ce{{}^{4}He}}\rightleftharpoons\smash{\ce{{}^{16}O}}+\gamma$, we have $\mathcal{\tilde{R}}_{n}=\{0,1\}$ and $\mathcal{\tilde{P}}_{n}=\{2,3\}$ with $[0]=\smash{\ce{{}^{12}C}}$, $[1]=\smash{\ce{{}^{4}He}}$, $[2]=\smash{\ce{{}^{16}O}}$, and $[3]=\gamma$. For $\smash{\ce{{}^{12}C}}+\smash{\ce{{}^{12}C}}\rightleftharpoons\smash{\ce{{}^{20}Ne}}+\smash{\ce{{}^{4}He}}+\gamma$, we have $\mathcal{\tilde{R}}_{n}=\{0,1\}$ and $\mathcal{\tilde{P}}_{n}=\{2,3,4\}$ with $[0]=[1]=\smash{\ce{{}^{12}C}}$, $[2]=\smash{\ce{{}^{20}Ne}}$, $[3]=\smash{\ce{{}^{4}He}}$, and $[5]=\gamma$. Note that reaction $n$ includes both the forward process (going from reactants to products) and the inverse process (going from products to reactants). Of course, it is arbitrary which set we call reactants and which we call products. In the following, we use the convention that if we consider particle $[\epsilon]$, then we choose the reactants and products such that $\epsilon\in\mathcal{\tilde{R}}_{n}$.

If the particles are uncorrelated, the system can be described in terms of the individual particle distribution functions $f_{\epsilon}$. Kinetic theory then gives the time evolution of the distribution functions as (e.g., Danielewicz & Bertsch 1991; Buss et al. 2012)

[TABLE]

where the sum over $n$ only includes interactions that have $\epsilon\in\mathcal{\tilde{R}}_{n}$. $\mathcal{N}^{\text{forward}}_{n}$ and $\mathcal{N}^{\text{inverse}}_{n}$ are factors that avoid double counting if the interaction involves multiple particles of the same species. These will be defined later after introducing some additional notation. $k_{\epsilon}^{\mu}=(k^{0}_{\epsilon},\bm{k}_{\epsilon})$ is the four-momentum of particle $[\epsilon]$ and $\delta^{4}$ is the 4-dimensional delta function that enforces conservation of momentum. $\mathcalligra{r}_{n}$ denotes the differential rate of reaction $n$. The upper ($+$) signs are for bosons and the lower ($-$) signs are for fermions. We use the shorthand

[TABLE]

for the phase-space integral of particle $[\epsilon]$, where $g_{\epsilon}$ is the spin degeneracy factor. Note that the differential rate $\mathcalligra{r}_{n}$ depends on the momenta of all particles (reactants and products). The second line in Section 2.1 is due to the forward process of interaction $n$, and the third line is due to the inverse process. The third line is required by the assumption of reversibility of interactions. Note that the differential rate is the same for the forward and inverse processes, and the delta function is also identical since it is an even function.

For simple interactions, e.g., weak interactions between nucleons and neutrinos, $\mathcalligra{r}$ is given by

[TABLE]

where $\langle|\mathcal{T}^{2}|\rangle$ is the spin-averaged reduced matrix element (averaged over the spins of both the initial and final states) of the interaction (see Brown & Sawyer (1997) for a discussion in the non-relativistic context). For more complicated interactions between nuclei, $\mathcalligra{r}_{n}$ could include transition probabilities between multiple internal states.

Generally, a reaction network only evolves some subset of the particles present in the system. For instance, photons are assumed to be in equilibrium at all times, and the electron and positron densities are determined by charge neutrality, so their number evolution does not need to be tracked explicitly. Therefore, it is useful to define the part of a reaction that only includes particles that will be tracked by the network as

[TABLE]

where $[j]$ is a reactant species, $[l]$ is a product species, and $\mathcal{R}_{n}$ and $\mathcal{P}_{n}$ (without tildes) are the sets of distinct reactant species and product species including only species that are present in the network, respectively. The positive integers $N_{j}^{n}$ and $N_{l}^{n}$ are the numbers of particles of reactant species $[j]$ destroyed and the number of particles of product species $[l]$ created, respectively. Note that we use Latin indices to refer to particle species, and we use Greek indices to refer to individual particles. In the earlier example of $\smash{\ce{{}^{12}C}}+\smash{\ce{{}^{4}He}}\rightleftharpoons\smash{\ce{{}^{16}O}}+\gamma$, we now have $\mathcal{R}_{n}=\{0,1\}$ and $\mathcal{P}_{n}=\{2\}$ with $[0]=\smash{\ce{{}^{12}C}}$, $[1]=\smash{\ce{{}^{4}He}}$, $[2]=\smash{\ce{{}^{16}O}}$, and $N_{0}^{n}=N_{1}^{n}=N_{2}^{n}=1$. The photon that is contained in $\mathcal{\tilde{P}}_{n}$ is not in $\mathcal{P}_{n}$, because it is not explicitly tracked in the network. Similarly, for $\smash{\ce{{}^{12}C}}+\smash{\ce{{}^{12}C}}\rightleftharpoons\smash{\ce{{}^{20}Ne}}+\smash{\ce{{}^{4}He}}+\gamma$, we get $\mathcal{R}_{n}=\{0\}$ and $\mathcal{P}_{n}=\{1,2\}$ with $[0]=\smash{\ce{{}^{12}C}}$, $[1]=\smash{\ce{{}^{20}Ne}}$, and $[2]=\smash{\ce{{}^{4}He}}$. But in this case, we have $N_{0}^{n}=2$ and $N_{1}^{n}=N_{2}^{n}=1$.

Since Section 2.1 essentially counts the pairs (or triplets, etc.) of reactants that can interact with each other (or pairs of products for the inverse process), we need to be careful to avoid double counting if the interaction involves multiple particles of the same species. If two distinct particles $[0]$ and $[1]$ interact with each other, then there are $N_{0}N_{1}$ distinct pairs, where $N_{0}$ and $N_{1}$ are the numbers of particles $[0]$ and $[1]$, respectively. But if we have $N$ particles of the same species where two react with each other, then the total number of distinct pairs is $N^{2}/2$ and not $N^{2}$. If it is three identical particles that react with each other, the number of distinct triplets is $N^{3}/6$, since there are $6=3!$ ways to order a set of three items. Thus, we need to divide by the product of factorials of the multiplicities of the interacting particle species. With the notation introduced in Equation 5, we can write this as

[TABLE]

We can now define the reaction rate of the forward process as

[TABLE]

where

[TABLE]

is the number density of species $[m]$ and $n_{B}$ is the total baryon number density. $\lambda_{n}$ is the forward process term of reaction $n$ on the right-hand side of Section 2.1 integrated over the phase space of particle $[\epsilon]$ and normalized by the number densities of the particles in the entrance channel. The reaction rate is just the rate at which a reaction proceeds per particle in the entrance channel. These reaction rates are only non-zero when the particles in the entrance channel differ from those in the exit channel. The other interactions included in Section 2.1 may change the momentum space distribution of the particles in the system, but they cannot change the total number of particles of any species. The reaction rate of the inverse process $\lambda_{n}^{\prime}$ is defined analogously to Equation 7 with the reactant and product sets switched.

Now, the standard reaction network equations follow from integrating over the phase space of particle $[\epsilon]$ in Section 2.1 to find

[TABLE]

Here, $\langle{\bm{v}}_{\epsilon}\rangle=(2\pi)^{-3}n_{\epsilon}^{-1}\int d^{3}k_{\epsilon}f_{\epsilon}\partial_{\bm{k}_{\epsilon}}k^{0}_{\epsilon}$ is the average velocity of particle $[\epsilon]$, and we define

[TABLE]

Note that the left-hand side of Equation 9 is for an individual particle, not a particle species, but we need the evolution equations for the particle species. A reaction $n$ that involves $N_{i}^{n}$ particles of species $[i]$ will contribute the right-hand side in Equation 9 $N_{i}^{n}$ times to the derivative of $n_{i}$ and so we multiply the right-hand side by $N_{i}^{n}$. Furthermore, due to the symmetry between the forward and inverse processes in Equation 9, it makes sense to treat the forward and inverse processes separately. So far, we have indexed reactions with $n$ and each reaction consisted of the forward and inverse directions. Let us now index reactions with $\alpha$, where the forward and inverse processes are counted separately. The set of reactions $\{\alpha\}$ is thus twice as big as the set of reactions $\{n\}$, although some inverse reactions may be ignored since they are extremely unlikely to occur. Equation 9 thus becomes

[TABLE]

where

[TABLE]

For every interaction $n$, there is a forward reaction $\alpha$ that has $\lambda_{\alpha}=\lambda_{n}$, $\mathcal{R}_{\alpha}=\mathcal{R}_{n}$, and $N_{\mathcal{R}}^{\alpha}=N_{\mathcal{R}}^{n}$. For that reaction, Equation 11 contributes the forward part of Equation 9 with a multiplicative factor $(-R_{i}^{\alpha}+P^{\alpha}_{i})N_{i}^{\alpha}=-N_{i}^{\alpha}$, since $N_{i}^{\alpha}$ particles of species $[i]$ are destroyed. Similarly, there is an inverse reaction $\alpha^{\prime}$ for the same interaction $n$ that has $\lambda_{\alpha^{\prime}}=\lambda_{n}^{\prime}$, $\mathcal{R}_{\alpha^{\prime}}=\mathcal{P}_{n}$, and $N_{\mathcal{R}}^{\alpha^{\prime}}=N_{\mathcal{P}}^{n}$. This reaction contributes the inverse part of Equation 9 with a factor of $(-R_{i}^{\alpha}+P^{\alpha}_{i})N_{i}^{\alpha}=N_{i}^{\alpha}$, since $N_{i}^{\alpha}$ particles of species $[i]$ are created in the inverse reaction. Note that $N_{m}^{\alpha}=N_{m}^{\alpha^{\prime}}=N_{m}^{n}$ for all reactants and products $[m]$.

Finally, it is useful to define the abundance $Y_{i}$ as

[TABLE]

where $V$ is the volume of the fluid element, and $N_{i}$ and $N_{B}$ are the total numbers of particles of species $[i]$ and baryons, respectively. Since the number density $n_{i}$ of species $[i]$ changes with both the number of particles of species $[i]$ and the total volume of the system, which is often changing in astrophysical systems undergoing nuclear burning, it is convenient to evolve the abundances $Y_{i}$ instead of the number densities $n_{i}$. Assuming that all of the species move as a single fluid, i.e., $\langle{\bm{v}}_{i}\rangle={\bm{v}}$, and using the Lagrangian time derivative, $d/dt=\partial_{t}+{\bm{v}}\cdot\nabla$, we find

[TABLE]

where we used the identity $\bm{v}\cdot\nabla f=\nabla\cdot(\bm{v}f)-f\nabla\cdot\bm{v}$ and the baryon number continuity equation $\partial_{t}n_{B}+\nabla\cdot(\bm{v}n_{B})=0$. But since $\langle\bm{v}_{i}\rangle=\bm{v}$, the above is the left-hand side of Equation 11 and so we get

[TABLE]

which is the familiar abundance evolution equation (e.g., Hix & Thielemann 1999). Essentially, for a given set of rates $\lambda_{\alpha}$, SkyNet solves this coupled, first-order, non-linear system of equations.

Even though it might look somewhat complicated, Equation 19 is very easy to understand. It says that the total time derivative of a species $[i]$ is the sum over all reactions that involve that species, and each reaction contributes a term consisting of the reaction rate multiplied by the abundances of the reactants and a factor that gives the number of particles of species $[i]$ destroyed or created in the reaction. For example, for the forward and inverse reactions $\smash{\ce{{}^{12}C}}+\smash{\ce{{}^{4}He}}\rightleftharpoons\smash{\ce{{}^{16}O}}$, Equation 19 says

[TABLE]

For the reactions $\smash{\ce{{}^{12}C}}+\smash{\ce{{}^{12}C}}\rightleftharpoons\smash{\ce{{}^{20}Ne}}+\smash{\ce{{}^{4}He}}$, we get

[TABLE]

2.2. Reaction rates and velocity-averaged cross-sections

Specializing to astrophysical systems consisting of a range of nuclear species, scattering reactions mediated by the nuclear and electromagnetic forces bring particles into thermal equilibrium at temperature $T$ on a much shorter timescale than nuclear reactions bring particles into chemical equilibrium. In that case, the distribution functions only depend on temperature and the chemical potentials, i.e., $f_{i}=f_{i}(T,\mu_{i})$. As written, the rates defined in Equation 7 depend on the momentum space distribution functions of the particles involved in the reaction and may be quite complicated. Nevertheless, in thermal equilibrium, the reaction rates depend only on the parameters of the distribution functions, i.e., $\lambda_{\alpha}=\lambda_{\alpha}(T,n_{B},\mu_{\{m\}})$, where the index $m$ ranges over reactant and products.

If all of the particles involved in a reaction obey Boltzmann statistics or are photons with chemical potential zero, then we find that $\lambda_{\alpha}=\lambda_{\alpha}(T,n_{B})$ does not depend on the chemical potentials of the particles. This is because we can set the blocking factors $(1\pm f_{i})$ of Boltzmann particles to 1, since quantum effects for Boltzmann particles are negligible. For Boltzmann particles, we have $n_{i}\propto\exp(\mu_{i}/T)$ and $f_{i}\propto\exp(\mu_{i}/T)$, which means that the dependence on $\mu_{i}$ cancels for Boltzmann particles because every $f_{\alpha}$ in Equation 7 is divided by an $n_{\alpha}$.555Actually, the $f_{\alpha}$’s are divided by $n_{j}^{N_{j}^{n}}$, but this is just an $n_{\alpha}$ for every particle of species $[j]$. In most astrophysical scenarios, we can assume that the distribution functions $f_{i}$ follow a thermal Maxwell–Boltzmann distribution, so the rates of reactions involving only nuclei and photons only depend on the temperature and density.

Nuclei can also undergo weak interactions that may involve leptons with non-zero chemical potentials. The leptons are generally not evolved in the network, but rather the electron (and positron) chemical potential is determined by the requirement of charge neutrality or by the number of electrons per baryon $Y_{e}$. Neutrinos are also not evolved in reaction networks since their distribution functions are often non-thermal in astrophysical scenarios in which they play a role. Therefore, weak interaction rates have a dependence $\lambda_{\alpha,\text{weak}}=\lambda_{\alpha,\text{weak}}(T,n_{B},Y_{e},f_{\{\nu\}})$ where $f_{\{\nu\}}$ are the externally specified neutrino distribution functions of the relevant neutrino species (see Appendix C). However, weak decay rates of nuclei are just constants when final state blocking by leptons can be safely ignored.

For two particle reactions, it is common to define the cross-section as (Peskin & Schroeder 1995, §4)

[TABLE]

where $v_{\text{rel}}$ is the relative velocity between particles [1] and [2]. Adopting the viewpoint that the $[i]$ are stationary targets and the $[j]$ are incoming projectiles impinging on the targets, the cross-section $\sigma_{\alpha}$ can be interpreted as (e.g., Clayton 1968, §4)

[TABLE]

where $R_{i,j}$ is the number of reactions per second, $r_{i,j}=R_{i,j}/V$ is the number of reactions per second per volume, and $n_{i}$, $n_{j}$ are the number densities of $[i]$ and $[j]$. Assuming Boltzmann statistics so that $(1\pm f_{l})\to 1$ for the products, Equation 7 gives

[TABLE]

where $\mathcal{N}_{\alpha}$ is the double counting factor from Equation 6. Since the distribution functions are normalized by the densities of species [1] and [2], this expression shows that $\lambda_{\alpha}$ is proportional to the cross-section averaged over the relative velocities between the two particles (after transforming to the center of mass frame; see e.g., Clayton 1968). Therefore, using Equation 27, one arrives at the standard relation between the reaction rate and the velocity-averaged cross-section (e.g., Clayton 1968, §4; Rolfs & Rodney 1988, §3),

[TABLE]

2.3. Nuclear Statistical Equilibrium (NSE) and inverse reaction rates

Section 2.1 shows that for every reaction, there is an inverse reaction. The relationship between the forward and reverse rates, which only depends on the density, temperature, and the internal properties of the nuclei, is generally called detailed balance. In some cases, for example for $\beta$-decays or fission reactions, the inverse reactions are extremely unlikely to occur and can be ignored. For other reactions, e.g., a neutron capture reaction, the inverse reaction can occur very frequently and sometimes even more often than the forward reaction. Thus, it is important to take inverse reactions into account. At temperatures of about 5 GK and above, inverse strong reactions such as photodissociation of nuclides can be in equilibrium with their forward reactions. In that case, the reaction is said to be in chemical equilibrium. For example, the reactions $\smash{\ce{{}^{196}Au}}+\text{n}\rightleftharpoons\smash{\ce{{}^{197}Au}}+\gamma$ and $\smash{\ce{{}^{20}Ne}}+\gamma\rightleftharpoons\smash{\ce{{}^{16}O}}+\smash{\ce{{}^{4}He}}$ can be in chemical equilibrium at sufficiently high temperatures. The situation of all strong reactions being in equilibrium is called NSE. This situation can also be thought of as an equilibrium between the reaction of forming a nucleus $(Z,N)$ from $Z$ free protons and $N$ free neutrons, and its inverse reaction, namely, completely dissociating a nucleus $(Z,N)$ into $Z$ protons and $N$ neutrons. In other words, if NSE holds, then the forward and inverse reactions,

[TABLE]

are in equilibrium for all nuclides that are part of the composition. Of course, there are no reactions that directly create a nuclide $(Z,N)$ out of $Z$ protons and $N$ neutrons. But there is a chain of strong reactions that connects $(Z,N)$ to free neutrons and protons. So if all strong reactions are in equilibrium, then we effectively have the reactions shown above and they are also in equilibrium. When nucleons are in chemical equilibrium with all other nuclear species, the energetic cost of turning $Z_{i}$ protons and $N_{i}$ neutrons into a single nucleus must be zero, which requires

[TABLE]

where $\mu_{i}$ is the chemical potential of species $[i]$.

When the composition moves into NSE, the forward and inverse strong reactions approach equilibrium. In order to ensure that the equilibrium composition determined by the forward and inverse reaction rates is the same as the NSE composition computed from the equality of the chemical potentials, we need to compute the inverse reaction rates directly from the forward rates and nuclide properties. Consider the reaction $\alpha$ and its inverse reaction $\alpha^{\prime}$. In equilibrium, each set of terms on the right-hand side of Section 2.1 must be zero. Then, by the symmetry of the differential rate $\mathcalligra{r}_{\alpha}=\mathcalligra{r}_{\alpha^{\prime}}$ and casting Section 2.1 into Equation 9, we have

[TABLE]

where $Y_{i,\text{eq}}$ is the abundance of species $[i]$ in chemical equilibrium. For a Boltzmann gas, the abundance is given by (Equation A12)

[TABLE]

where $G_{i}(T)$ is the internal partition function of species $[i]$ (see Equation A4) and $m_{i}$ is its rest mass. Substituting the above into Equation 32 yields

[TABLE]

Since the forward and reverse reactions are in equilibrium, the chemical potentials on both sides are equal, hence

[TABLE]

and the chemical potentials in the exponential of Equation 34 cancel. Then, the inverse reaction rate $\lambda_{\alpha^{\prime}}$ is

[TABLE]

where we define

[TABLE]

Although the above is derived under the assumption that the abundances are such that the two reactions are in equilibrium, it still holds for any abundances because the reaction rates only depend on the temperature and density. In Appendix B, we show how calculating NSE is implemented in SkyNet.

3. Network evolution

In this section, we focus on the specific implementation of the physics described in the previous section in SkyNet. In essence, evolving the reaction network means solving the large coupled system of first-order, non-linear ordinary differential equations (ODEs) given in Equation 19. But there are various other pieces that are needed to make the evolution robust and efficient, which are also discussed this section.

3.1. Implicit integration method

The system of ODEs (Equation 19) we need to solve is extremely stiff because of the enormous range of reaction rates, which can span many orders of magnitude (e.g., Timmes 1999; Hix & Meyer 2006). Thus, an explicit integration method would be constrained to taking extremely small time steps. This is why nuclear reaction networks are typically integrated with an implicit method (e.g., Arnett & Truran 1969; Woosley et al. 1973; Timmes 1999; Winteler 2013; Longland et al. 2014). However, some authors have proposed various explicit methods specifically tuned to integrate stiff reaction networks (e.g., Feger 2011; Guidry 2012; Guidry & Harris 2013; Guidry et al. 2013a, b; Brock et al. 2015). Currently, SkyNet uses the first-order implicit backward Euler method (e.g., Gear 1971; Hix & Thielemann 1999), but it is straightforward to implement higher-order implicit integration methods as well. We plan to do this in the future, since Timmes (1999) recommends using higher-order methods, such as the variable-order Bader–Deuflhard method (Bader & Deuflhard 1983). Let the vector $\bm{Y}(t)=Y_{i}(t)$ denote the composition at time $t$. If we want to take a time step $\Delta t$ using the first-order implicit backward Euler method, we write

[TABLE]

where $\bm{x}=\bm{Y}(t+\Delta t)$ are the unknown abundances at the end of the time step we are trying to find and $\bm{Y}(t)$ are the known abundances at the beginning of the time step. $\dot{\bm{Y}}(\bm{Y},T,\rho)$ is the function defined in Equation 19 that gives the time derivatives of the abundances as a function of a given set of abundances $\bm{Y}$, a temperature $T$, and density $\rho$. Note that we need to know the temperature $T$ and density $\rho$ as a function of time. The function $\bm{F}$ is simply the right-hand side of the above equation. We thus have a multidimensional root-finding problem $\bm{0}=\bm{F}(\bm{x},T,\rho)$, where $T$ and $\rho$ are known functions of time. SkyNet uses the Newton–Raphson (NR) method to find the root. For every time step, the following iteration is performed to find the solution $\bm{x}$ to Equation 41:

[TABLE]

where the iteration starts with $\bm{x}_{0}=\bm{Y}(t)$, the abundances from the previous time step, and

[TABLE]

is the Jacobian matrix and $\delta_{ij}$ is the Kronecker delta. The partial derivatives $\partial\dot{Y}_{i}/\partial Y_{j}$ are computed from Equation 19 as

[TABLE]

The size of the Jacobian matrix $J_{\bm{F}}$ is $N\times N$, where $N$ is the number of nuclear species in the network. For large networks ($N\sim 8000$), inverting this large matrix can be quite costly, because the linear system in Equation 41 has to be solved for every NR iteration. Since the Jacobian matrix depends on the unknown abundances at the end of the time step, we have to recompute the Jacobian matrix after every NR iteration that updates our guess for the abundances at the end of the time step. Fortunately, however, the Jacobian matrix is very sparse (only 0.24% of the $N^{2}$ entries are non-zero for large networks), because most nuclear species are not directly connected by a single reaction. By exploiting the sparseness of the Jacobian, we drastically reduce the memory requirement to store the Jacobian and we can also use matrix solver packages that are specifically designed for sparse systems (see §5). The most expensive parts of the evolution are computing the Jacobian entries from Equation 44 and then solving the sparse linear system given by Equation 41. These two operations have to be performed at every NR iteration. The method for choosing the time step $\Delta t$ and when to terminate the NR iteration will be discussed in §3.3.

3.2. Self-heating evolution

In the previous section, we showed how SkyNet integrates the nuclear abundances forward in time if both the temperature and density are given as functions of time. In most applications, the density history $\rho(t)$ is given (for example, when SkyNet is used to post-process nucleosynthesis for tracer particles from a hydrodynamical simulation), but the temperature is not necessarily known as a function of time. Even if we do have a temperature history available, it most likely would not include heating due to the nuclear reactions that SkyNet evolves. But based on the kinetic theory description of the reaction network equations in §2.1, as well as the discussion of detailed balance in §2.3, it is clear that the reaction network and the thermodynamic state of the fluid are intimately related and need to be treated consistently. Therefore, we want the temperature to be evolved in the network under the influence of nuclear reactions, which release nuclear binding energy as heat. This is referred to as a self-heating network evolution (e.g., Freiburghaus et al. 1999b). We still require knowing the density as a function of time, though.

Recall the first law of thermodynamics,

[TABLE]

where $dU$ is the infinitesimal change in internal energy, $\delta Q$ is the infinitesimal heat added to the system from the surroundings, and $\delta W$ is the infinitesimal mechanical work performed by the system. If the system expands or contracts, the work performed is $\delta W=PdV$, where $P$ is the pressure and $dV$ is the infinitesimal change in volume. We use the entropy $S$, volume $V$, and composition $\{N_{k}\}$ as our independent thermodynamic variables. Note that the index $k$ ranges over all particles in the system, not just nuclides. The total differential of the internal energy is thus

[TABLE]

where $T$ is the temperature and $\mu_{k}$ is the chemical potential of particle $k$. Equating Equation 45 and Equation 46 yields

[TABLE]

If we divide the above by $N_{B}$, the total number of baryons, and replace the infinitesimal changes by the differences of the quantities from one time step to the next (i.e., $\Delta X=X(t+\Delta t)-X(t)$), we find

[TABLE]

since $Y_{k}=N_{k}/N_{B}$ (Equation 17). Note that the sum over $k$ includes all particles, and hence also electrons, which can be created or destroyed in weak nuclear reactions. $\Delta q$ is the change in the heat per baryon due to external heating sources. Let $\dot{\epsilon}_{\text{ext}}$ be an imposed external heating rate (per baryon), then

[TABLE]

where $\dot{\epsilon}_{\nu}$ is the neutrino heating/cooling rate of the system on the environment given by Equation C19. Since we do not include neutrinos in the internal state of the system, the neutrino heating/cooling must be treated as an external source of heat. And because we define $\dot{\epsilon}_{\nu}$ as the heating/cooling on the environment, it has a minus sign in Equation 48. Combining Equation 48 with Equation 49 and solving for $\Delta s$ yields

[TABLE]

where we used Equation A12 and switched to index $i$ that runs only over the nuclides. So, we explicitly include the contribution from the $Z_{i}$ electrons that come with nuclide $i$. To make the rest mass terms in the sum closer to unity, we define the mass excess $\mathcal{M}_{i}$ as

[TABLE]

where $A_{i}$ is the number of neutrons and protons of species $i$, and $m_{u}$ is the atomic mass unit, defined such that the mass excess of \ce^12C is exactly 0 (i.e., $m_{u}=m_{\smash{\ce{{}^{12}C}}}/12$). Since $Y_{i}=N_{i}/N_{B}$ (with $N_{i}$ being the number of particles of species $[i]$; see Equation 17), we find

[TABLE]

because species $[i]$ is made up of $A_{i}$ neutrons and protons, and hence $A_{i}$ baryons. Thus, we have $\sum_{i}Y_{i}(t)A_{i}=1$ for all times $t$, which is another way of saying that the total baryon number $N_{B}$ is conserved. Using this, we find

[TABLE]

Using the above and the fact that $n_{i}=Y_{i}n_{B}$, we can write Equation 50 as

[TABLE]

Note that the external heating is accounted for with a first-order forward Euler method. We plan to improve this in the future when we implement higher-order integration methods for the network itself. With the above, SkyNet can update the entropy after every time step and then obtain the new temperature at the end of the time step from the EOS. Therefore, we only need to know the initial entropy (or temperature, from which the entropy is determined). This changes the evolution (Equation 41) slightly, because we now have to use the entropy at the beginning of the time step to estimate the temperature at the end of the time step. That is, we solve

[TABLE]

where $T^{\ast}$ is given by the EOS as

[TABLE]

Equation 55 is solved with the NR method, as described in the previous section. Note that $\Delta t$ is fixed during the NR iterations, which means that the temperature and density are also fixed. After the NR iterations have converged, we have found the new abundances $Y_{i}(t+\Delta t)$ and then we can compute $\Delta s$ according to Equation 54 and update the entropy as

[TABLE]

Hence, we have a hybrid implicit/explicit scheme where the abundances are evolved implicitly but the entropy is evolved explicitly. One could also evolve the entropy implicitly together with the abundances, which would require computing $\partial\dot{Y}_{i}/\partial s$ and $\partial\dot{s}/\partial Y_{i}$ and adding these terms to the Jacobian. We may extend SkyNet in the future to support such a fully implicit scheme, but for now, we have achieved good results with the hybrid approach.

The energy released as nuclear binding energy due to nuclear reactions is

[TABLE]

where we used Equation 53 and $\Delta Y_{i}/\Delta t$ is an approximation for $\dot{Y}_{i}$ over the time step. Note that the minus sign comes from the fact that some rest mass (or mass excess) is converted into energy, and hence the heating rate is positive if there is a net reduction in the total rest mass. Some authors (e.g., Hix & Thielemann 1999) treat $\dot{\epsilon}_{\text{nuc}}$ as an external heat source. This is necessary if the EOS does not depend on the entire composition but only on $\bar{A}$ and $\bar{Z}$, the average mass and charge numbers, for example. In order to compare the nuclear heating in SkyNet to other codes, SkyNet computes and records a total heating rate $\dot{\epsilon}_{\text{tot}}$, regardless of whether self-heating is enabled. $\dot{\epsilon}_{\text{tot}}$ is computed as

[TABLE]

Note that the above has units of erg s*-1* baryon*-1*. To convert it to the more commonly used units of erg s*-1* g*-1*, we simply multiply $\dot{\epsilon}_{\text{tot}}$ by the Avogadro constant $N_{A}$.

Currently, SkyNet records the total heating rate shown above. In reality, this heating rate is composed of multiple components that are thermalized in the material in different ways (e.g., Barnes et al. 2016). For example, emitted electrons and positrons, as well as the kinetic energy of fission fragments, thermalize with very high efficiency, while only a small fraction of the energy released as neutrinos might thermalize. In a future version of SkyNet, we plan to record the different heating rate components, so that thermalization can be taken into account in kilonova light curve calculations, for example.

3.3. Convergence criteria and time stepping

The time step for the network evolution needs to be adjusted depending on how well the NR iteration (Equation 42) converges. All of the default values and thresholds mentioned in this section are adjustable by the user. To check that the NR iterations have completely converged, a standard criterion is (Press et al. 2007)

[TABLE]

where $x_{i}^{(n+1)}$ is the $i$th component of the vector $\bm{x}_{n+1}$, and $\bm{x}=\bm{Y}(t+\Delta t)$ are the unknown abundances at the end of the current time steps that we want to find. The sum only runs over the indices $i$ for which $x^{(n+1)}_{i}\geq Y_{\text{thr}}$ for some abundance threshold $Y_{\text{thr}}$, which we usually set to $10^{-20}$. The default value is $\varepsilon_{\text{tol},\Delta x}=10^{-6}$. Although this convergence criterion ensures that any subsequent NR iterations would no longer change the solution $\bm{x}_{n+1}$, we found that this criterion is too strict in practice. Instead, SkyNet typically uses mass conservation as a heuristic convergence criterion (which is also used by Hix & Thielemann (1999)), which takes the form

[TABLE]

where we usually use $\varepsilon_{\text{tol,mass}}=10^{-10}$. Note that the sum now runs over all nuclear species, and there is no threshold for $x_{i}^{(n+1)}$. Since $x_{i}^{(n+1)}=Y_{i}(t+\Delta t)$, this convergence criterion is simply the conservation of the total baryon number. The user of SkyNet can choose to use Equation 60, Equation 61, or both as the convergence criterion for Equation 42.

Figure 1 shows an r-process evolution with the two different convergence criteria using $\varepsilon_{\text{tol},\Delta x}=10^{-6}$ and $\varepsilon_{\text{tol,mass}}=\smash{10^{-10}}$. These convergence thresholds result in almost exactly the same time-step sizes, but if we made $\varepsilon_{\text{tol},\Delta x}$ smaller, that would result in much smaller time steps. However, using the $\Delta x$ convergence criterion requires an average of 3.1 NR iterations per time step, while mass conservation only needs 1.2 NR iterations per time step. Since the total number of time steps is almost the same, using mass conservation as the convergence criterion is about 2.1 times faster for this particular case. As Figure 1 shows, however, the nucleosynthesis evolution is identical in the two cases. No differences are visible in the entropy, electron fraction, neutron abundance, or average mass number as a function of time. The maximum absolute difference in the final abundances of the two cases (using the $\Delta x$ or mass conservation convergence criterion) is about $1.6\times 10^{-7}$.

SkyNet adjusts the time-step size $\Delta t$ dynamically. Once the NR iterations have converged according to the chosen criterion, SkyNet checks that the composition did not change too much over the last time step. The temperature and entropy are allowed to change by at most 1%. If either of them changes by more than this threshold, then the time step is considered failed and $\Delta t$ is reduced by a factor of two, and the whole step is attempted again with the reduced time-step size. SkyNet also considers the time step as failed if the NR iterations do not converge after 10 iterations, or if the error measure used for the NR convergence criterion increases compared to the error of the previous iteration, or if the error decreases by less than 10%. In all of those cases, $\Delta t$ is reduced and the time step is attempted again. A simplified schematic of this mechanism is shown in Figure 2.

After a successful time step, SkyNet attempts to increase the step size for the next time step. SkyNet tries to double $\Delta t$ after every successful step, but this new time step can be limited if the abundance of a particular nuclide changed by a large amount in the previous time step. If that is the case, then the new time step is limited to the approximate step size necessary to keep the abundance of the nuclide that changed the most over the last time step from changing by more than 10%. Hence, the new time step is computed as

[TABLE]

where $\Delta t_{\text{max}}$ is the maximum allowed time step and $\Delta t^{(n)}$ is the previous time step. Using this adaptive time-stepping mechanism, we typically get a time-step size that grows exponentially with time in freely expanding trajectories while keeping the error measure used for the convergence criterion below its prescribed tolerance.

Very rarely, it is necessary to renormalize the composition. In that case, every abundance is divided by the total mass, i.e.,

[TABLE]

and then the new composition satisfies $\sum_{i}A_{i}Y_{i}=1$ exactly. Although this artificially injects or removes energy from the system, it is useful as a last resort if the time-step size is kept small because the composition is far away from mass conservation (but still within the error tolerance). After renormalization, the evolution usually proceeds normally with a larger time step than before. We renormalize if the time step falls below a certain limit (usually $10^{-16}$), or if there are more than 25 time steps in a row that tried to increase the step size but subsequently failed and had to keep the step size constant. In such cases, it could be that the time step is small because the mass conservation convergence criterion is preventing the time step from increasing. If this is the case, then renormalizing the abundances usually helps to increase the time step, because after renormalizing, mass conservation in Equation 61 is fulfilled exactly. But in some cases, for example when trying to evolve the network near NSE with reaction rates that are inconsistent with NSE (see §6.1.1), the time step is small because the abundances are changing rapidly and so renormalizing the composition does not help.

3.4. NSE evolution mode

If the abundances approach the NSE composition, the forward and inverse strong rates exactly balance (§2.3). In that case, all the partial derivatives in the Jacobian (Equation 44) would be zero, resulting in a singular Jacobian. The Jacobian is not exactly singular, however, because the weak reactions (that are not in equilibrium with their inverses) contribute non-zero derivatives to the Jacobian. Nevertheless, as the strong reactions move into equilibrium, the network time step becomes very small as the Jacobian becomes close to being numerically singular. To alleviate this problem, SkyNet automatically switches from a full network evolution to an NSE evolution scheme, if the strong nuclear reaction time scale becomes shorter than the time scale over which the density changes, and if the temperature is above some threshold (a user setting with a default value of 7 GK). The full network is turned back on when these conditions are no longer satisfied. A similar approach is used by other groups (e.g. Iwamoto et al. 1999; Brachwitz et al. 2000).

If SkyNet determines that switching to NSE evolution is appropriate, it computes the NSE composition from the current internal energy, density, and electron fraction. If the entropy and temperature of that NSE composition differs by less than 1% (user setting) from the current network entropy and temperature, then the switch to NSE is allowed. Otherwise, the full network evolution will continue and SkyNet will try to switch to the NSE evolution mode again after the next step. A test of the NSE evolution mode that demonstrates its necessity and consistency is presented in §6.3.

In the NSE evolution mode, SkyNet no longer evolves the abundances of all nuclear species. Instead, SkyNet only evolves the entropy $s$ and electron fraction $Y_{e}$ of the composition, which can change due to weak reactions, such as $\beta$-decays or neutrino interactions, which can change the charge of nuclides and heat the material. Recall that the electron fraction is $Y_{e}=\sum_{i}Z_{i}Y_{i}$, and so

[TABLE]

where $\dot{Y}_{i}$ is given by Equation 19 as a function of $T$, $\rho$, and $\bm{Y}$. The temperature is given by the EOS as a function of $s$, $\rho$, and $\bm{Y}$. $\bm{Y}$ is given by NSE as a function of $s$, $\rho$, and $Y_{e}$ (see §B.2). Thus, we have

[TABLE]

The rate of change of the entropy is obtained from dividing Equation 54 by $\Delta t$, namely,

[TABLE]

Since $\dot{\epsilon}_{\nu}(t)$, $\dot{\epsilon}_{\text{ext}}(t)$, and $\rho(t)$ are known, we thus have two coupled ODEs for $Y_{e}$ and $s$, which SkyNet integrates with the Runge–Kutta–Fehlberg 4(5) method (e.g., Burden et al. 2015, §5.5). This is a fourth-order explicit integration method that also computes a fifth-order error estimate that is used to adaptively control the integration time step. The heating rate can be calculated analogously to Equation 59 as

[TABLE]

Note that in the NSE evolution mode, we only evolve two variables and they are changing on similar timescales because they are both influenced by the weak reactions. In this case, however, even though the weak reactions span a large range of timescales, this does not introduce any stiffness because we only deal with the sum of the abundance derivatives in both Equations 64 and 69. Thus, we can safely use an explicit integration method.

4. Electron screening

Nuclear reaction rates strongly depend on the Coulomb interaction between the nuclides in the entrance channel (e.g., Salpeter 1954). In conditions where nuclear burning can occur, the nuclides are almost all fully ionized. Therefore, nuclear reaction rates can be computed assuming that bare nuclei of charge $Z_{i}$ interact with each other. At high temperature, the Coulomb interaction energy of the electrons with the nuclei is small compared to the thermal energies of the particles and so the screening effect is negligible. But if the density of the medium is sufficiently high and the temperature sufficiently low, electron screening is likely important. In these conditions, a nucleus will repel neighboring nuclei and attract nearby electrons, thus creating an electron charge cloud around the nucleus. This charge cloud partially screens or shields the nuclear charge $Z_{i}e$, where $e$ is the elementary charge unit. Thus, the Coulomb repulsion between the two positively charged nuclei is reduced by the screening effect, which can enhance the nuclear reaction rates, which depend strongly on the probability of Coulomb barrier penetration. Obviously, screening corrections are only important for charged particle reactions. Neutron capture reactions are unaffected by the polarization of the electron gas.

In this section, we present how electron screening is implemented in SkyNet. Our focus will be on writing down the equations that SkyNet uses to compute the screening corrections in a useful way with adequate justification. We will not develop the screening theory from first principles but refer the reader to the established literature on this subject (e.g., Salpeter 1954; Dewitt et al. 1973; Graboske et al. 1973; Itoh et al. 1979; Ichimaru & Utsumi 1984; Brown & Sawyer 1997; Bravo & García-Senz 1999; Yakovlev et al. 2006, to name but a few). For a handful of reactions, screening has been investigated experimentally (e.g., Engstler et al. 1988; Rolfs & Somorjai 1995; Chen et al. 2004; Gatu Johnson et al. 2017). We use Gaussian cgs units throughout this section.

The strength of the electron screening effect depends mainly on the ratio of the Coulomb interaction energy between a nucleus and the nearby electrons to the thermal energy. If the thermal energy is large compared to the Coulomb interaction energy, then the electron charge cloud around the nucleus will be large and diffuse, providing less screening to the nuclear charge. We define the ion density

[TABLE]

The average interionic spacing is

[TABLE]

Now, we define the dimensionless screening parameter $\Lambda_{0}$ as

[TABLE]

where $\beta=1/T$. For some average (dimensionless) charge per ion $\zeta$ (defined later on in Equation 96), let

[TABLE]

Thus, $\Lambda$ is a measure of the ratio of the average Coulomb interaction energy $\zeta^{2}e^{2}/a$ to the average thermal energy $T$. If $\Lambda$ is large, we expect the screening effect to be strong, and when $\Lambda$ is small, screening should be weak. The screening corrections in these two regimes, as well as the intermediate regime where $\Lambda\sim 1$, are the subject of the following sections, after we introduce the general screening factor that modifies the nuclear reaction rate.

4.1. General screening factor

For the two-body reaction $[1]+[2]\to[3]$, the total Coulomb potential can be written as

[TABLE]

where $U(r_{12})$ is a potential correction to the bare Coulomb interaction between the two nuclei due to screening, $r_{12}$ is the separation of the two reactants, and $Z_{1}$ and $Z_{2}$ are the charge numbers. Salpeter (1954) showed that the screening correction to the nuclear reaction rate $\lambda_{12}$ is given by

[TABLE]

where $U_{0}=U(r_{12}=0)$, $\lambda_{12,\text{no-sc}}$ is the unscreened reaction rate, and $f_{\text{sc}}=e^{-U_{0}/T}$ is the general screening factor. Note that $f_{\text{sc}}$ is sometimes written as $\exp(H_{12}(0))$ in the literature. An approximation for the screening factor $f_{\text{sc}}$ can be found in terms of $Z_{1}$, $Z_{2}$, density, temperature, and other quantities determined by the composition. However, it is advantageous to (equivalently) write down the screening factor in terms of a chemical potential correction $\mu_{\text{sc}}(Z)$ that depends on the charge of a single nucleus. This way, we can apply the screening corrections in NSE as well, and we are not constrained to only correcting reaction rates with two reactants. In this section, we show how to compute the reaction screening factor if we have the chemical potential correction $\mu_{\text{sc}}(Z)$. In the following sections, we will show how to compute $\mu_{\text{sc}}(Z)$ in different screening regimes.

For the same two-body reaction $[1]+[2]\to[3]$, Dewitt et al. (1973) found that the screening factor can also be written as

[TABLE]

where $\beta=1/T$, $\mu_{\text{sc}}(Z)$ is the correction to the chemical potential for the addition of a charge $Z$ to the system when the electron gas is non-uniform, and $Z_{3}=Z_{1}+Z_{2}$ is the charge of the product. It is straightforward to generalize this result to an arbitrary (i.e. not just two-body) reaction $\alpha$ as

[TABLE]

where $\mathcal{R}_{\alpha}$ is the set of reactant species of the reaction and $N_{i}^{\alpha}$ is the number of species $[i]$ destroyed or produced in the reaction (see §2.1 and Equation 5). We also define

[TABLE]

The above expression for the screening factor $f_{\text{sc}}$ for reactions with an arbitrary number of particles in the entrance channel holds for multistep reactions. This can be shown by breaking up the multibody reaction into two-body reactions, e.g., treating $[1]+[2]+[3]\to X$ as $[1]+[2]\to[12]$ followed by $[12]+[3]\to X$, and then calculating the overall screening factor. To justify the general form of $f_{\text{sc}}$, consider a two-body reaction that produces an arbitrary set of products $\mathcal{P}_{\alpha}$, i.e.,

[TABLE]

where $Z_{1}+Z_{2}=\sum_{j\in\mathcal{P}_{\alpha}}N_{j}^{\alpha}Z_{j}$ due to charge conservation. Let $\lambda_{\alpha,\text{no-sc}}$ and $\lambda_{\alpha^{\prime},\text{no-sc}}$ be the reaction rates of the forward and inverse reactions without the screening corrections. The corrected forward rate is

[TABLE]

Equation A12 gives the chemical potential of nuclear species $[i]$ without the screening correction, which we call $\mu_{i,\text{no-sc}}$. The corrected chemical potential is $\mu_{i}=\mu_{i,\text{no-sc}}+\mu_{\text{sc}}(Z_{i})$. Note that the chemical potential correction $\mu_{\text{sc}}(Z_{i})$ enters on the same level as the term $m_{i}$ in Equation A12 and therefore, $\mu_{\text{sc}}(Z_{i})$ can be absorbed into $Q_{\alpha}$, the rest mass difference between the reactants and products, defined in Equation 37. $Q_{\alpha}$ thus becomes

[TABLE]

Substituting the above corrected expressions for $\lambda_{\alpha}$ and $Q_{\alpha}$ into Equation 36 yields

[TABLE]

where we used charge conservation, i.e., $Z_{1}+Z_{2}=\sum_{j\in\mathcal{P}_{\alpha}}N_{j}^{\alpha}Z_{j}=Z_{\mathcal{P}}^{\alpha}$. Thus, the screening factor for the inverse reaction (whose set of reactants is $\mathcal{P}_{\alpha}$) is indeed exactly the generalized screening factor postulated in Equation 78. Thus, we only need to know how to compute $\mu_{\text{sc}}(Z)$, which is the subject of the following sections. Note that in terms of the screened forward rate, we combine the above with detailed balance (Equation 36) to obtain the screened inverse rate

[TABLE]

where we used $Z_{\mathcal{P}}^{\alpha}=Z_{\mathcal{R}}^{\alpha}$ since strong reactions conserve charge.

It now remains to find the screening chemical potential correction $\mu_{\text{sc}}(Z)$ for a nuclide of charge $Z$. Thanks to SkyNet’s modularity, arbitrary expressions for $\mu_{\text{sc}}(Z)$ can be plugged into the general screening framework. In the following sections, we describe the current implementation of $\mu_{\text{sc}}(Z)$ in SkyNet.

4.2. Weak screening

The weak screening limit is the limiting case where the Coulomb interaction energy is much lower than the thermal energy, hence, where $\Lambda\ll 1$. In this case, the electrostatic Poisson–Boltzmann equation describing the screening can be solved approximately to find (Salpeter 1954)

[TABLE]

where $\lambda_{D}$ is the Debye screening length. $\kappa_{D}=\lambda_{D}^{-1}$ is called the Debye wave number, and it is given by (Brown & Sawyer 1997)

[TABLE]

where the Debye wave number of species $[i]$ is

[TABLE]

where $Z_{i}e$ is the charge, $\mu_{i}$ is the chemical potential (defined in Equation A12), and $f_{i}(p,\mu_{i})$ is the distribution function of species $[i]$.

SkyNet assumes that the ions (nuclides) are non-degenerate and non-relativistic. Therefore, they obey Boltzmann statistics with $E(p)=p^{2}/(2m_{i})+m_{i}$, where $p$ is the momentum and $m_{i}$ is the mass of species $[i]$. Hence,

[TABLE]

and so

[TABLE]

and hence

[TABLE]

where $n_{i}$ is the number density of ion species $[i]$.

The electrons and positrons are allowed to be arbitrarily degenerate and arbitrarily relativistic. They are both fermions and thus follow Fermi-Dirac statistics, so

[TABLE]

where $E(p)$ is the total energy

[TABLE]

We find

[TABLE]

and so

[TABLE]

where the extra factor of 2 in the phase-space integral comes from the fact that fermions have two spin states, and we used $Z_{e^{\pm}}=\pm 1$.

The Debye length is thus

[TABLE]

where

[TABLE]

Note that the factor of $n_{B}/n_{I}$ in $\zeta$ occurs because we define $\zeta$ as the RMS charge per ion rather than per nucleon. The distribution functions $f_{e^{\pm}}$ (Equation 91) depend on the electron and positron chemical potentials $\mu_{e^{\pm}}$, which are given by the EOS (Equations A18 and A19). SkyNet uses the adaptive QAG integration routines provided by the GNU Scientific Library666https://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html to evaluate the above integral numerically. The two-body screening potential in Equation 85 becomes

[TABLE]

and so the screening factor is

[TABLE]

This is consistent with the result from Dewitt et al. (1973), who showed that in the weak screening case, the chemical potential correction due to screening is

[TABLE]

Combining the above with the expression for the general screening factor (Equation 78) yields

[TABLE]

as shown in Equation 98.

4.3. Strong and intermediate screening

In the strong screening limit, the Poisson–Boltzmann equation governing screening has to be solved numerically. Dewitt et al. (1973) found that

[TABLE]

where $c_{0}=9/10$, $c_{1}=0.2843$, $c_{2}=-0.054$, $d_{0}=-9/16$, and $d_{1}=0.4600$. The parameter $\Gamma_{0}$ is

[TABLE]

and the applicable average charge per ion is the arithmetic mean $\bar{Z}$, instead of the RMS $\zeta$, which is given by

[TABLE]

We can write Equation 101 in terms of $\Lambda_{0}$ as

[TABLE]

which is valid when $\Gamma=\bar{Z}^{2}\Gamma_{0}\gg 1$. Alternatively, if $\zeta\sim\bar{Z}$, then $\Gamma\sim\Lambda^{2/3}$ and so the strong screening limit applies if $\Lambda\gg 1$. The strong screening limit is also applicable if $\Lambda\ll 1$ but the charge $Z$ is such that $Z\zeta^{2}\Lambda_{0}\gg 1$. Then strong screening applies for charge $Z$ (Dewitt et al. 1973).

The intermediate screening regime is where $\Lambda=\zeta^{3}\Lambda_{0}\sim 1$, in which case Dewitt et al. (1973) found

[TABLE]

where $b=0.860$ and

[TABLE]

4.4. Combining the different screening regimes

Equations 99, 104, and 105 give the chemical potential corrections in the limits of weak, strong, and intermediate screening, respectively. To smoothly transition between these three regimes, we need a single parameter that determines which regime is applicable and a function that smoothly interpolates between the regimes based on the said parameter.

Weak screening applies if $\Lambda=\zeta^{3}\Lambda_{0}\ll 1$ and strong screening if $\Lambda\gg 1$ or $\Lambda\ll 1$ but $Z\zeta^{2}\Lambda_{0}=\Lambda Z/\zeta\gg 1$. We thus define the dimensionless parameter

[TABLE]

Note that $p(Z)\ll 1$ if and only if $\Lambda\ll 1$ and $\Lambda Z/\zeta\ll 1$, in which case weak screening applies. If and only if $\Lambda\gg 1$ or $\Lambda Z/\zeta\gg 1$ is $p(Z)\gg 1$, in which case we are in the strong screening regime. If $p(Z)\sim 1$, then intermediate screening applies. To ensure a smooth transition of $\beta\mu_{\text{sc}}$ from one regime to another, we will compute the screening correction as a weighted sum of the corrections computed in the different regimes. Hence, we compute

[TABLE]

where each coefficient $f_{j}$ is between 0 and 1 and defined as

[TABLE]

Note that $f_{w}(p)+f_{s}(p)+f_{i}(p)=1$ for all values of $p$. The transition from weak to intermediate screening occurs at $p(Z)=1/\sqrt{25}=0.2$, and the transition from intermediate to strong screening occurs at $p(Z)=\sqrt{25}=5$.

Figure 3 demonstrates the transition between the different screening regimes. We choose a composition consisting of 59% neutrons, 40% protons, and 1% \ceU^238 (by mass) at $T=3$ GK and $\rho=2\times 10^{9}\ \text{g cm}^{-3}$. In this composition, all three screening regimes occur as the charge number of a test particle ranges from $Z=1$ to $Z\sim 100$. In order to show the transitions more clearly, we plot the chemical potential corrections from the different regimes for $Z$ ranging from 0.5 to 500 in Figure 3. Such charge numbers are not expected to occur in reality, but the expressions for the chemical potential corrections are valid nonetheless. The chosen composition results in $\Lambda_{0}=0.05089$, $\zeta=0.9104$, $\bar{Z}=0.4079$, $\eta_{b}=0.6214$, and $\Lambda=0.03840$.

Whereas Figure 3 shows the screening correction for different ion charges in a fixed composition, Figure 4 shows the screening correction for a fixed charge ($Z=26$) in different compositions determined by the temperature $T$. Also shown is the Boltzmann chemical potential (Equation A12 without the rest mass) for \ce^56Fe. The composition is determined from NSE with screening corrections (§4.5) using a temperature ranging from 3 to 13 GK, $\rho=10^{8}\ \text{g cm}^{-3}$, and $Y_{e}=0.4$. At low temperatures ($T\lesssim 7$ GK), strong screening is applicable since the thermal energy is small and Coulomb interactions dominate. Electron screening provides a 10% correction over the unscreened Boltzmann chemical potential. Intermediate screening is applicable between 7 and 9 GK. At 9 GK, \ce^4He nuclei are broken up and free neutrons and protons start to dominate the composition. Thus, $\zeta$ remains roughly constant at 1, and so $\Lambda$ and $p(26)$ also stop changing rapidly because $\Lambda_{0}$ does not depend strongly on the temperature. $T\sim 9$ GK is also where weak screening sets in, because the thermal energy now overcomes the Coulomb interaction energy, which has been reduced by the smaller average charge per ion provided by the free protons. In that regime, the screening effect is at the 4% to 0.5% level.

Figure 4 also shows that our screening corrections are indeed approximate. The weak and intermediate screening corrections never meet, so any scheme to transition between them is necessarily approximate and arbitrary to some degree. But considering the fact that our transition scheme needs to work robustly for a wide range of compositions and ion charges, it seems to do reasonably well in interpolating between the different (somewhat disjoint) screening regimes. Although progress has been made in improving screening calculations in various regimes, a unifying theory for screening across all regimes is still elusive (e.g., Itoh et al. 1977; Shaviv & Shaviv 1996, 2000; Chugunov et al. 2007).

4.5. NSE with screening

As was noted above, in addition to screening nuclear reactions, electronic correlations also change the free energetic cost of adding or removing a charged particle from the medium. Therefore, the chemical potentials of the nuclides have an additional correction due to screening, i.e., $\mu_{i}=\mu_{i,\text{no-sc}}+\mu_{\text{sc}}(Z_{i})$, where $\mu_{i,\text{no-sc}}$ is the unscreened chemical potential given by Equation A12 and $\mu_{\text{sc}}(Z_{i})$ is the screening correction given by Equation 108 with $Z_{i}$ being the charge number of nuclide $i$. Appendix B presents in detail how NSE is computed. Here, we briefly discuss how the method shown in Appendix B is modified to take chemical potential corrections into account. Equation B2 becomes

[TABLE]

which means Equation B8 gives

[TABLE]

where $\beta=1/T$. Note that $\mu_{\text{sc}}(Z_{i})$ depends on all $Y_{i}$ because the screening corrections depend on different types of average ion charges ($\bar{Z}$, $\zeta$, and $\eta_{b}$). Computing NSE involves an NR iteration that requires partial derivatives of $Y_{i}$ (see Appendix B). Thus, the screening corrections in Equation 113 introduce a large number of complicated partial derivatives to the Jacobian, which we will not attempt to write down. We experimented with using numerical derivatives to compute the Jacobian, with limited success. Another complication is that Equation 113 itself depends on all $Y_{i}$ on the right-hand side and thus has to be solved iteratively.

We find that it is much more robust to keep the screening corrections fixed during the NR iterations. This does not introduce any additional derivatives in the Jacobian. NSE is computed exactly as shown in Appendix B, with the only difference that Equation 113 is used to compute $Y_{i}$ from $\eta_{i}$, but the terms $\beta\mu_{\text{sc}}(Z_{i})$ are constant throughout the NR iterations. In this case, to obtain an NSE composition that is self-consistent with the screening corrections it is based on, the NSE computation itself needs to be iterated.

We start by computing the NSE composition without screening (i.e., $\mu_{\text{sc}}(Z_{i})=0$). We denote the resulting composition by $\bm{Y}^{(0)}$. Then, we compute $\mu_{\text{sc}}(Z_{i})$ based on $\bm{Y}^{(0)}$ and use these as the constant screening corrections for the next NSE computation that yields $\bm{Y}^{(1)}$. From this new composition we compute new screening corrections, which are used to compute $\bm{Y}^{(2)}$ and so on. This iteration stops once $\max_{i}(|Y_{i}^{(n+1)}-Y_{i}^{(n)}|)<10^{-12}$ or $n$ reaches 20 (both of these criteria can be changed by the user). This method of iteratively updating the screening corrections and computing NSE with them being fixed is not guaranteed to converge (but neither is the NR method itself). However, in practice, we find that this method works very well and converges quite quickly in a large region of parameter space.

Figure 5 shows the NSE abundance distribution as a function of mass number $A$ for three different temperatures. For all temperatures, $\rho=10^{8}\ \text{g cm}^{-3}$ and $Y_{e}=0.4$, so this is the same composition as the one shown in Figure 4. To show the impact of the screening corrections, the NSE compositions with and without screening are shown in the left panel. Screening is strongest for $T=3$ GK and the effect of screening is to reduce the abundances below $A\sim 80$ and slightly increase them above that mass number. Screening is weaker at $T=8$ GK, but still enhances the high-mass abundances above $A\sim 60$. Finally, for $T=13$ GK, screening has virtually no effect. In the right panel of Figure 5, we show the ratio of the screening chemical potential $\mu_{\text{sc}}$ to the Boltzmann chemical potential $\mu_{\text{BZ}}$ (Equation A12 without the rest mass). $\mu_{\text{sc}}$ depends on the composition and the charge, hence it is the same for all isotopes of a given element. But the Boltzmann chemical potential also depends on the mass, partition function, and the abundance of a given isotope. Thus, for a given charge number $Z$, the ratio $\mu_{\text{sc}}/\mu_{\text{BZ}}$ varies for the isotopes of that charge. In Figure 5, we show the range of the chemical potential ratios as colored bands. The bands collapse to a line at large $Z$, because there the abundances of all isotopes are essentially zero. In the strongest screening case ($T=3$ GK), the screening effect ranges from 0.1% to 70%. For $T=8$ GK, it ranges from about 0.02% to 10%, and for $T=13$ GK, the screening effect is much less than 1% except for very large $Z$.

5. Implementation details

The main design goals of SkyNet are usability and flexibility. Since SkyNet is built in a modular fashion, different physics implementations can easily be switched out or new physics can be added, making SkyNet very flexible (see the next sections for details). To achieve the modularity, SkyNet is written in object-oriented C++ and makes use of some C++11 features. SkyNet contains a small amount of Fortran code to provide a minimal interface for SkyNet to be called from Fortran. CMake (http://www.cmake.org) provides a cross-platform, compiler-independent build system for SkyNet that automatically finds the required external libraries. CMake also provides an automated testing facility and SkyNet comes with a suite of tests that check basic functionality and correctness of SkyNet.

To make it easy to use, SkyNet comes with Python bindings that make it possible to use all parts of SkyNet from an interactive Python shell or a Python script. Therefore, one can use standard Python libraries like NumPy (http://www.numpy.org) to read in and manipulate input data, like the list of nuclides to be evolved, the initial composition, density vs. time history, etc., and these data can be passed to SkyNet using standard Python data structures. This means that one does not have to deal with C++ to run SkyNet. But of course, SkyNet can also be used from a C++ or Fortran application. One can even run multiple copies of SkyNet in parallel within Python, using Python’s multiprocessing module—a facility we use extensively whenever post-processing nucleosynthesis on many tracer particles form hydrodynamical simulations (e.g., Lippuner et al. 2017; Roberts et al. 2017). An example of how to run SkyNet in parallel with Python is included with the SkyNet source code available at https://bitbucket.org/jlippuner/skynet. The SkyNet Python bindings are provided by SWIG (http://www.swig.org) and using Python is the most convenient and most flexible way to run SkyNet.

5.1. Modularity and extendability

SkyNet is a modular library of different C++ classes rather than a monolithic program. Some of the most important classes in SkyNet are the various reaction library classes that contain different types of nuclear reactions (see next section), a nuclide library class that contains all nuclear data, and a reaction network class that implements the actual nuclear reaction network. There are various other types of classes that implement specific functionalities. For example, there are different function interpolation classes, ODE integrators, and general numerical method classes (bisection, line search). On the physics side, there are different classes that are responsible for different pieces of physics. The NSE class computes NSE given an electron fraction and two of the following properties: temperature, density, entropy, or internal energy. There are also separate classes that are responsible for the EOS and screening corrections.

Since SkyNet is built in an object-oriented fashion, different parts of the code are separated from each other and only interact via well-defined interfaces. This makes SkyNet extremely modular because the implementation of a certain class can be changed or extended, without having to modify the rest of the code. For example, one could easily extend the nuclide library class to support reading nuclear data from a different file format. Since all of the nuclear data are handled by this one class, only this class has to be modified to support the new file format. Furthermore, some classes are implemented as abstract base classes, meaning they only specify the interface for a particular physics module without tying it to a specific implementation. Examples of this are the EOS class and the screening corrections class. For both of these, SkyNet currently has one implementation, namely the extended Timmes EOS described in §A.2 and the screening corrections discussed in §4. One can easily add a new EOS class that implements a different EOS but has the same interface as the abstract EOS base class. This new EOS class then plugs into the existing SkyNet framework. In a similar way, one can add additional screening implementations to SkyNet.

The various classes provided in SkyNet can be used individually through the Python bindings. For instance, one can use the NSE class in Python to compute NSE in various conditions or use the nuclide library class to access the nuclear data and partition functions from Python.

5.2. Nuclear reaction libraries

SkyNet supports different types of nuclear reactions. Reactions of the same type or from the same data source are grouped into reaction library classes. The network class contains an arbitrary list of reaction library classes that collectively contain all of the reactions that are evolved in the network. The reaction library classes have a common interface that allows the network to be agnostic as to how the reaction rate is determined. Via this interface, the network can tell the reaction libraries to recompute the reaction rates for a given thermodynamic state (temperature, density, electron fraction, electron degeneracy parameter, etc.) to get the contributions to all $\dot{Y}_{i}$ from the reactions in the network and to get the contributions to $\partial\dot{Y}_{i}/\partial Y_{j}$. This makes SkyNet extremely flexible because many different types of reactions can be evolved at the same time, and furthermore, the data for reactions of the same type can be split across multiple files, allowing the user to quickly switch out certain reactions. Finally, thanks to the abstract interface of reaction library classes, it is easy to add new types of reactions to SkyNet.

The following reaction types of nuclear reactions are currently implemented in SkyNet.

• •

Constant: These reactions have a constant rate that does not depend on any properties of the thermodynamic state.

• •

REACLIB: These are reactions that come from the REACLIB database (Cyburt et al. 2010). The rates of these reactions are given by parametric fitting formulae that depend on temperature and density.

• •

Tabulated: This reaction library contains tabulated $\beta$-decay (both $\beta^{-}$ and $\beta^{+}$) and electron/positron capture rates based on Fuller et al. (1982) and Langanke & Martínez-Pinedo (2001). The rates are tabulated as a function of temperature and $Y_{e}\rho$.

• •

Neutrino interactions: These are neutrino emission and absorption reactions on free neutrons and protons. The rates are calculated according to Equations C9, C10, C11, and C12 given the electron neutrino and electron antineutrino distribution functions.

• •

Arbitrary rate functions: This reaction library contains reactions whose rates are given by arbitrary, user-specified functions. This can be used to quickly test a new or modified reaction rate that can depend on various thermodynamic quantities and also time.

Since the different reaction types and rate sources can be used concurrently in SkyNet, care must be taken to ensure that no reaction rate is contained multiple times in the network, since that would effectively multiply the reaction rate by the number of times it occurs. SkyNet provides a facility to remove all reactions in one reaction library that also occur in another library. However, in some cases, there are supposed to be multiple rates for the same reaction. In this case, the total reaction rate is the sum of the individual rates. REACLIB uses this mechanism to capture the different resonant and non-resonant parts of a reaction rate with its limited fitting formula.

6. Code verification and tests

In order to verify the correctness of SkyNet, we compare its results to the results of other existing reaction network codes, specifically WinNet (Winteler 2013) and XNet (Hix & Thielemann 1999), and also to results published in the literature. The scripts and input files to reproduce these code tests are distributed with SkyNet in the directory examples/code_tests. SkyNet also has a test suite that contains simple code tests, regression tests, and tests that compare very simple networks to analytic solutions. The main purpose of that test suite is to ensure that changes to the code do not break the functionality or correctness of SkyNet.

6.1. Nuclear statistical equilibrium

To verify the NSE solver in SkyNet, we perform a consistency test and comparison to literature results of the NSE abundances computed by SkyNet. We use nuclear masses and partition functions distributed with REACLIB (Cyburt et al. 2010), which contains experimental data where available and finite-range droplet macroscopic model (FRDM; see, e.g., Möller et al. 2016) data otherwise. The temperature-dependent partition functions are from Rauscher & Thielemann (2000).

6.1.1 Consistency test

In this section, we verify that the abundances computed with the NSE solver in SkyNet are consistent with the strong reactions. We perform a test evolution starting with purely free neutrons and protons and let only strong reactions take place. We use the strong reaction rates from REACLIB and the default fission rates distributed with SkyNet. The network contains 7824 nuclear species, ranging from free neutron and protons to \ce^337Cn ($Z=112$). We keep the temperature constant at $T=7$ GK and the density $\rho=10^{9}\ \text{g cm}^{-3}$ is also fixed. We pick these values to ensure that the composition achieves NSE within a reasonable amount of time. We set $Y_{e}=0.4$, so the initial composition is $Y_{\text{n}}=0.6$ and $Y_{\text{p}}=0.4$, and since the network does not contain any weak reactions, the electron fraction remains constant at $Y_{e}=0.4$. Screening corrections are enabled for both the NSE solver and the network evolution. At every step, we compare the network abundances to the NSE abundances and compute the error $\Delta Y$ as

[TABLE]

We perform two different network evolutions: one where the strong inverse rates are computed from detailed balance (§2.3) and another one where the inverse rates from REACLIB are used. We never use inverse fission reactions.

The results are shown in Figure 6. After evolving the network using detailed balance to compute the inverse rates for 850 time steps, corresponding to $t=1$ s, the network reaches the NSE composition. The error between NSE and the network composition is $\Delta Y\sim 10^{-11}$ and the deviation of $Y_{e}$ from 0.4 is on the same level. We note that this is comparable to the mass conservation limit of $10^{-10}$ that SkyNet uses as the NR iteration convergence criterion. The neutron and proton abundances also match the values from the NSE composition with very high precision. In the subsequent 250 time steps, the network reaches $t=10^{10}$ s and $\Delta Y$ decreases another order of magnitude. This demonstrates that the NSE solver in SkyNet and the implementation of detailed balance for the inverse rates are consistent. The NSE compositions computed with SkyNet are indeed the compositions that the network produces if the strong reactions are in equilibrium.

Furthermore, Figure 6 also shows that the inverse rates provided in REACLIB are not completely consistent with the NSE composition that is computed from the nuclear data (masses and partition functions) distributed together with REACLIB. If the REACLIB inverse rates are used (“SkyNet no DB” in Figure 6), the network evolution is extremely slow after about 900 steps. We stopped the network evolution after 130,000 steps at $t\sim 0.053$ s, when it became clear that $\Delta Y$ converged to $2.5\times 10^{-4}$. The evolution without detailed balance becomes very slow because the inverse rates from REACLIB try to push the composition into a certain equilibrium configuration, but chemical potential balance predicts a different equilibrium composition. This makes the evolution very difficult and keeps the time step between $10^{-15}$ and $10^{-5}$ s. Computing the inverse rates from detailed balance so that the inverse rates exactly cancel the forward rates when the chemical potentials balance is therefore necessary for the network evolution to be consistent with NSE (see §2.3).

A detailed investigation of the reverse rates in REACLIB reveals that a significant fraction of them were computed with nuclear masses different from the mass model distributed with REACLIB, which we use in SkyNet to compute the inverse rates from detailed balance and to compute NSE. We observe this for the latest REACLIB version (REACLIB V2.2 from 2016-11-14 and also newer pre-release version that is not publicly available yet) as well as previous REACLIB versions. For 74% of the inverse rates, the $Q$ values used in the REACLIB inverse rates differ by 1% or less compared to the $Q$ values computed from our mass model (Equation 37). About 19% of the $Q$ values agree within 1% to 10%, 5% within 10% to 30%, and almost all the $Q$ values agree within a factor of 2. However, there are some rare cases where the $Q$ value used in REACLIB differs by up to two orders of magnitude from the one given by the mass model used in SkyNet. Fortunately, the impact of the discrepancy between the reverse REACLIB rates and the mass model on the final abundances of a network evolution is very small, as we shall see in §6.2. Nevertheless, the test presented in this section highlights the advantage of using detailed balance with the mass model used in the network to compute inverse rates. This is especially important to keep the time step from becoming very small when the composition moves into NSE. Some of the inconsistencies in the inverse rates will be corrected in a future REACLIB version (Schatz 2017). However, because REACLIB is a collection of rates from different sources that use different mass models to compute the rates and their inverses, it seems unlikely that all inverse rates in REACLIB will ever be consistent with a single mass model. We therefore recommend to always directly compute inverse rates from detailed balance with the masses used in the network when consistency with NSE is desired.

6.1.2 Comparison with published results

We compare the SkyNet NSE solver to the NSE results with and without electron screening by Seitenzahl et al. (2009). The composition includes 443 nuclides ranging from free neutrons and protons to multiple isotopes of krypton (see Figure 1 in Seitenzahl et al. 2009). We compute NSE with and without screening at a fixed density of $\rho=5\times 10^{8}\ \text{g cm}^{-3}$ and fixed $Y_{e}=0.5$ for temperatures ranging from 4 to 10 GK. The results are shown in Figure 7 along with the results from Seitenzahl et al. (2009, Figures 3 to 6, data used with permission).

We find excellent agreement between the SkyNet results and those presented in Seitenzahl et al. (2009). When screening corrections are included, the deviation between the two results is slightly larger, because Seitenzahl et al. (2009) use a different screening implementation. The effect of the different screening implementation is most pronounced for the proton mass fraction at $T\gtrsim 9$ GK, but even then, the difference is less than 5%. Seitenzahl et al. (2009) use a fit for the screening corrections that has a significant correction even in the weak screening regime at $T\gtrsim 8$ GK. On the other hand, the weak screening correction in SkyNet becomes much smaller at those temperatures. Hence, the differences in the NSE mass fractions are due to the increased disagreement between the different screening corrections as the temperature increases. Furthermore, Seitenzahl et al. (2009) use different nuclear masses, so that could account for the small differences between our results and theirs when screening is turned off.

6.2. Network evolution

In this section, we present comparisons of nucleosynthesis evolutions with SkyNet and other nuclear reaction networks. We compare SkyNet to WinNet and XNet. WinNet was originally developed at the University of Basel by Winteler (2013) based on the earlier BasNet by Thielemann et al. (2011). WinNet has been used by various authors for r-process nucleosynthesis calculations in core-collapse supernovae and neutron star mergers, and to investigate the impact of nuclear physics on the r-process (e.g., Korobkin et al. 2012; Winteler et al. 2012; Eichler et al. 2015; Martin et al. 2015, 2016). XNet was developed at Oak Ridge National Laboratories by Hix & Thielemann (1999) and has been used for r-process nucleosynthesis in accretion disk outflows and neutron star mergers, and for explosive nucleosynthesis in type I X-ray bursts and core-collapse supernovae (e.g., Surman et al. 2006; Fisker et al. 2008; Roberts et al. 2011; Harris et al. 2017).

Since nuclear physics data, such as nuclear masses, partition functions, and nuclear reaction rates, have a significant influence on the nucleosynthesis calculations, we take care in ensuring that exactly the same nuclear physics input data are used for all of the different codes that we consider. However, this means that we are restricted to using the greatest common denominator of nuclear physics data sources that can be used by all codes. For the comparisons in this section, we use the strong and weak reaction rates distributed in REACLIB (Cyburt et al. 2010), neutron-induced fission reactions with symmetric fission fragments from Panov et al. (2010), and spontaneous fission rates calculated from the approximation of Frankel & Metropolis (1947) using the spontaneous fission barriers of Mamdouh et al. (2001). In future versions of SkyNet, we plan to add additional fission reactions and fission fragment distributions. The nuclear masses and partition functions are again the ones distributed with REACLIB, as in §6.1.

6.2.1 Neutron-rich r-process

We run an r-process nucleosynthesis calculation in a neutron-rich environment with all three networks. We use 7836 nuclear species and about 93,000 reactions. The density history is a trajectory from the ejecta of a black hole–neutron star merger (Roberts et al. 2017). The initial composition is NSE with $T=6.1$ GK, $\rho=7.4\times 10^{9}\ \text{g cm}^{-3}$, and $Y_{e}=0.07$. We run all combinations of screening and self-heating turned on and off. For each case, we perform two separate SkyNet evolutions: one where the inverse rates are computed from detailed balance, and another where the inverse rates from REACLIB are used. We consider these two cases because SkyNet is usually run with inverse rates computed from detailed balance, but WinNet and XNet use the inverse rates from REACLIB, and so we also run SkyNet with those inverse rates for a more direct comparison.

The self-heating method currently implemented in XNet only applies in the case of constant density (Harris 2017). Hence, for this r-process computation with an evolving density, we cannot use the self-heating capability of XNet. Instead, to compare XNet to SkyNet and WinNet when self-heating is turned on, we use the SkyNet temperature history (from the SkyNet run without detailed balance) in XNet. However, the temperature provided to XNet has to be limited to a lower bound of 0.01 GK. At temperatures lower than that, the reaction rate fits from REACLIB no longer apply and some of the rates blow up. SkyNet internally also uses a lower bound 0.01 GK for the REACLIB reactions, but the network temperature used in the EOS is allowed to drop below that bound, until the lower limit of the EOS is reached at around $4\times 10^{5}$ K (this lower bound is due to the tabulated electron/positron part of the Timmes EOS).

Figure 8 shows the results of running the neutron-rich r-process with the three reaction networks with screening and self-heating turned on. We find excellent agreement between the different networks with the temperature evolution of SkyNet and WinNet being virtually indistinguishable (for XNet, we prescribe the SkyNet temperature). The final abundances (at $t=5\times 10^{8}$ s) also agree very well. To compare the results of the three networks quantitatively, we compute the numeric error between the final mass-summed abundances as

[TABLE]

where the mass-summed abundance $Y(A)$ is given by

[TABLE]

This error measure is the average absolute difference between the abundance results divided by the average abundances. This effectively measures the fractional error in the final abundances, averaged over all mass numbers $A$. However, we compute the quotient of the sums rather than the sum of quotients, because the latter would be dominated by tiny abundances ($\sim 10^{-20}$) that may differ by a factor of several between the two networks. This would result in a large overall error, but abundance differences at the $10^{-20}$ level are not important, even if it is by a factor of several.

Table 1 shows the errors in percent between the different networks. Since this is a neutron-rich environment, we expect that screening plays no important role in the nucleosynthesis evolution. The fact that the errors between the different networks are almost the same regardless whether screening is turned on or off confirms that screening is not important in this case. Furthermore, the error between SkyNet without detailed balance and WinNet or XNet are comparable to the errors between WinNet and XNet. For example, with screening and self-heating turned on (first row in Table 1), the error between SkyNet without detailed balance and WinNet is 3.4%, while the error between WinNet and XNet is 3.4%, and the error between XNet and SkyNet without detailed balance is 0.066%. SkyNet is closer to XNet in this case because XNet uses the same temperature as SkyNet whereas WinNet evolves its own temperature. This demonstrates that SkyNet produces results that are compatible with WinNet and XNet. The errors between SkyNet and WinNet or XNet are larger if detailed balance is used to compute the inverse rates in SkyNet. For example, again for the first row in Table 1, the error between WinNet and SkyNet with detailed balance is 5.8%, but only 3.4% when compared to SkyNet without detailed balance. This is not surprising, because SkyNet is effectively evolving slightly different reaction rates when inverse rates are computed from detailed balance. It does illustrate, however, that using detailed balance, which produces inverse rates that are consistent with the nuclear masses and partition functions, has a measurable impact and might be a reasonable standard practice. In the self-heating runs, SkyNet is much closer to XNet than WinNet, because XNet uses the temperature from SkyNet.

6.2.2 X-ray burst

In Figure 9, we show a comparison between the three networks for a different type of trajectory. This trajectory captures the situation of unstable hydrogen burning on the surface of a neutron star, which produces a type I X-ray burst, presented in Schatz et al. (2001). The density and temperature histories as well as the initial and final abundances were graciously provided by Schatz et al. (2001). The temperature starts at 0.2 GK and peaks at 1.9 GK during the burst. The density starts at $1.1\times 10^{6}\ \text{g cm}^{-3}$. The initial composition is 66.0% hydrogen (by mass), 33.6% helium, and 0.4% heavier elements, mostly oxygen. For this test, we use a small network containing only 686 species going up to \ce^136Xe and 8400 reactions.

We again see good agreement in Figure 9 between the three networks in the proton and helium abundance evolution (left panel of Figure 9). But there are a handful of mass numbers at which SkyNet with detailed balance produces final abundances that deviate from the other networks by up to two orders of magnitude (right panel of Figure 9). This indicates that it is vital to compute the inverse reaction rates correctly with detailed balance in this scenario. However, the average fractional abundance errors shown in Table 1 are still at the few percent level, and the differences between SkyNet and WinNet or XNet are comparable to the differences between WinNet and XNet. We note that the errors between SkyNet without detailed balance and WinNet or XNet are smaller if screening is turned off. This indicates that screening is somewhat important in this case and we expect the discrepancy between the different codes to increase if screening is turned on due to the different screening implementations. The effect is especially noticeable when comparing SkyNet to XNet with self-heating turned on, because in this case XNet uses the temperature computed from SkyNet, but if screening is switched off, the error between the two codes decreases from 1.1% to 0.086%.

SkyNet’s screening implementation is presented in §4. WinNet uses a single fit of the two-body screening factor by Chugunov et al. (2007) across all screening regimes. XNet computes the two-body screening function provided by Graboske et al. (1973) for weak and intermediate screening and the one provided by Dewitt & Slattery (1999) for strong screening. XNet then uses a selection rule to select one of the three screening regimes without interpolating between them.

6.2.3 Hydrostatic carbon-oxygen burning

In order to compare the self-heating methods implemented in the three networks, we perform a hydrostatic burn at constant density. We keep $\rho=10^{7}\ \text{g cm}^{-3}$ fixed and start with $T=3$ GK and the initial composition consists of half \ce^12C and half \ce^16O (by mass). As a baseline comparison, we also perform non-self-heating runs, where we keep the temperature fixed at 3 GK. The increase in the errors between the networks with self-heating enabled compared to with it disabled must then be due to the difference in the self-heating implementations in the codes. We use a mid-size network containing all nuclides from the full network (§6.2.1) with $A\leq 100$. This results in a network with 1530 species and 20,000 reactions.

Figure 10 shows the results of this test case with self-heating and screening turned on. We find good qualitative agreement between the three networks, but quantitatively, the differences between XNet and the other two codes are much larger than in the previous test cases. These discrepancies come from the different self-heating implementations in the three codes. WinNet uses the same self-heating method as SkyNet (described in §3.2), but WinNet does not include the entropy change due to the electron chemical potential. Also, WinNet uses the original Timmes EOS, which computes the entropy with a single representative heavy ion species, while SkyNet computes the entropy by considering all species in the network separately (§A.2). However, these two differences have virtually no impact on the temperature evolution or the final abundances produced by WinNet and SkyNet. XNet, on the other hand, uses a different self-heating method. It evolves the temperature directly using (Harris 2017)

[TABLE]

where $c_{V}$ is the specific heat capacity at constant volume (provided by the Timmes EOS), $\dot{Y}_{i}=dY_{i}/dt$ is the abundance time derivative of species $i$, and $\mathcal{M}_{i}$ is the mass excess. As can be seen in the left panel of Figure 10, this method is comparable to the methods in SkyNet and WinNet. Note that XNet is using exactly the same EOS as WinNet. We see only very small differences in the final abundances of the three networks in the right panel of Figure 10.

The differences in the self-heating implementations are also apparent in Table 1. When self-heating is turned on, XNet differs from SkyNet and WinNet by about 4%, and WinNet differs from SkyNet by about 1%. However, when self-heating is disabled, the three networks agree at the 0.2% level, and if screening is turned off, SkyNet without detailed balance and XNet agree to an astounding precision of 0.0004%.

6.3. NSE evolution test

To ensure that the NSE evolution mode in SkyNet produces the correct results, we perform a test that evolves a trajectory with and without the NSE evolution mode. Of course, the trajectory must experience some heating that forces the composition into NSE at some point during the evolution, otherwise the NSE evolution mode would not be triggered. We first attempt this test with a trajectory from a neutron star merger accretion disk outflow simulation (Lippuner et al. 2017). That trajectory experiences late-time fallback, which causes a spike in the density that results in late-time heating and forces the composition into NSE. Although SkyNet is able to evolve this trajectory without issues using the NSE evolution mode, when the NSE evolution mode is turned off, SkyNet gets stuck with a time step of $\sim 10^{-16}$ s for at least 350,000 steps at the time when the heating occurs. Thus, we cannot use this trajectory for this test, since we cannot evolve it without the NSE evolution mode. However, this trajectory serves as an illustration of the necessity of the NSE evolution mode in order to evolve certain trajectories.

Since it is challenging to evolve a trajectory without the NSE evolution mode that moves into NSE during the evolution, we use an artificial trajectory that has a temperature peak that is less than 8 GK. In practice, we found that a trajectory with a peak of 7.6 GK can be evolved without the NSE evolution mode, but anything hotter becomes problematic. So, we use an imposed temperature history that starts out at 5.5 GK and remains constant at that value, except for a short peak up to 7.6 GK at 10 s. We keep the density fixed at $\rho=10^{8}\ \text{g cm}^{-3}$. We evolve SkyNet without screening or self-heating, but with inverse rates from detailed balance, as these are important for consistency with NSE (§6.1.1). We use the same full network (7836 species and 93,000 reactions) as in §6.2.1. The initial composition is NSE with $Y_{e}=0.1$. Figure 11 shows the results of this test. From the left panel, we see that the free neutrons are captured in the first second to make heavy nuclei, raising the average mass number $\bar{A}$ to about 60. The nuclei synthesized are mainly around the first r-process peak and the composition decays toward stability until the temperature starts to rise at around 8 s. The rising temperature forces the material into NSE, which liberates neutrons and protons from nuclei, quickly reducing $\bar{A}$ to around 10. When the temperature drops back to 5.5 GK, the free neutrons and protons are absorbed into nuclei again and the iron peak is formed. In the left panel of Figure 11, there are no visible differences between the SkyNet evolution with and without NSE evolution mode. They produce exactly the same results. The right panel of Figure 11 contains the time step $\Delta t$ and temperature as a function of time. We see that the time steps with and without the NSE evolution mode are the same until SkyNet turns on the NSE evolution mode when the temperature is sufficiently high. With the NSE evolution mode, SkyNet evolves with $\Delta t\sim 10^{-1}$ s until the temperature drops again and the NSE evolution mode is turned off. Without the NSE evolution mode, however, the time step drops to a $\text{few}\times 10^{-5}$ s. After the temperature has returned to 5.5 GK, the time steps in both cases are virtually identical again. Not shown are the final abundances of the two SkyNet evolutions, but we find no discernible differences down to a level of $10^{-30}$. The error according to Equation 115 between the two SkyNet evolutions is 0.0098%, indicating that the NSE evolution mode in SkyNet produces correct results, i.e., exactly the same results as would be obtained with the full network, but that the NSE evolution mode prevents the time step from getting stuck at a very small value.

7. Summary and future work

We presented the new nuclear reaction network SkyNet and the physics that it currently implements. Details are provided of how the abundance evolution equations implemented in SkyNet are derived from kinetic theory. We discuss how inverse reaction rates are computed with detailed balance and how this is related to NSE. A detailed description is given of the numerical methods used in SkyNet for the network integration and the self-heating evolution that accounts for heating due to nuclear reactions. Further, we show how SkyNet automatically transitions between evolving the full network and evolving only the entropy and electron fraction under the influence of weak reactions when the composition is in NSE and all the strong reactions are in equilibrium. A general treatment of electron screening that computes the screening factors for arbitrary strong reactions from chemical potential corrections that only depend on the nucleus charge and the composition is presented. We then show in detail how SkyNet computes these chemical potential corrections in the weak, strong, and intermediate screening regimes and how it smoothly transitions between the regimes. These screening corrections are also taken into account when computing NSE compositions.

After providing some code implementation details that highlight the modularity and expandability of SkyNet, we present comprehensive code tests and comparisons. The NSE compositions computed by SkyNet are shown to be consistent with the evolved strong reactions, but only if detailed balance is used to compute the inverse rates. Furthermore, we show that SkyNet’s NSE results are compatible with results in the published literature. SkyNet is compared to two other nuclear reaction networks in three different cases: a neutron-rich r-process, a proton-rich explosive X-ray burst, and hydrostatic carbon/oxygen burning. All three tests are conducted with and without electron screening and with and without self-heating. We find that all three networks agree with each other at the few percent level in most cases, although there are some situations where the disagreement is larger due to somewhat different implementations of the physics in the codes. Finally, in the appendices, we discuss the physics of ideal Boltzmann gases and how this is implemented in the SkyNet EOS that accounts for all nuclear species individually. Technical details of how NSE is calculated in SkyNet and how the neutrino interactions are implemented are also presented in the appendices.

We hope that SkyNet will be a useful tool for the nuclear astrophysics community to compute nucleosynthesis in various scenarios. We also hope that the theoretical and experimental nuclear physics communities will find SkyNet useful as a low-barrier entry point to running nucleosynthesis models. SkyNet can be used for testing the impact of newly measured or calculated reaction rates or nuclear properties, or to conduct sensitivity studies in order to determine which nuclides should be the focus of future experiments. SkyNet is available as an open-source software at https://bitbucket.org/jlippuner/skynet. We value any feedback, whether it be bug reports, feature requests, or code contributions.

In the future, we plan to extend the electron screening implementation in SkyNet in a way that screening corrections can be consistently accounted for in the EOS as well. We are also thinking of investigating different screening prescriptions and making them available in SkyNet. Currently, only the first-order backward Euler method is implemented in SkyNet, but we intend to add higher-order integration methods in the future. SkyNet is limited to evolving a one-zone model right now (i.e., there is only one density and temperature). We plan to add support for multiple zones in SkyNet and couple it to existing hydrodynamical simulations. Finally, we plan to investigate offloading some or most of the computations in SkyNet to GPUs, which could significantly speed up the network evolution and might be necessary to efficiently evolve large hydrodynamical simulations that are coupled to SkyNet.

We thank J. Austin Harris and W. Raph Hix for giving us access to XNet, for assisting us with running the code comparison tests with XNet, and for helpful discussion of the test results. We are also very grateful to Moritz Reichert, Dirk Martin, Oleg Korobkin, and Marius Eichler for giving us access to WinNet and for helping us run the code comparisons with it. We thank Hendrik Schatz for providing the X-ray burst trajectories and Ivo Seitenzahl for allowing us to use his NSE data. We are very grateful to Christian D. Ott for carefully reading the manuscript and providing helpful comments and suggestions. We also thank our referee, Friedel Thielemann, for the thoughtful comments and additional reference suggestions. L.R. thanks W. Raph Hix and Stan Woosley for many useful discussions concerning nuclear reaction networks and nucleosynthesis. We are indebted to the authors of the sparse linear solver package PARDISO (Schenk & Gärtner 2004, 2006; Karypis & Kumar 1998) for making their code available for free for academic uses. This work was supported in part by the Sherman Fairchild Foundation and NSF awards NSF CAREER PHY-1151197, TCAN AST-1333520, and AST-1205732. The calculations were performed on the Zwicky compute cluster at Caltech, supported by NSF under MRI award PHY-0960291 and by the Sherman Fairchild Foundation. This work benefited from access to the NSF XSEDE computing network under allocation TG-PHY100033 and to NSF/NCSA Blue Waters under NSF PRAC award ACI-1440083. This work was supported in part by NSF grant PHY-1430152 (JINA Center for the Evolution of the Elements).

Appendix A Equation of state (EOS)

SkyNet requires an EOS in order to relate different thermodynamic quantities, such as temperature, entropy, chemical potential, etc., to each other. Ions behave as non-relativistic, non-degenerate particles in the majority of situations where nuclear burning occurs and so the EOS in SkyNet treats all ions as non-interacting, non-degenerate, non-relativistic ideal Boltzmann gases. Electrons and positrons, on the other hand, can be both degenerate and relativistic. An important exception to the assumption that the ions are non-interacting is electron screening, which is discussed in detail in §4.

In this section, we present a brief summary of the most relevant properties of ideal Boltzmann gases and introduce the notation used in this paper. We also describe the EOS implemented in SkyNet.

A.1. Ideal Boltzmann gas

The grand partition function $\mathcal{Z}$ of a Maxwell–Boltzmann gas is given by (e.g., Reichl 1980, §9.D.3)

[TABLE]

where the sum runs over all single-particle states $k$ with energy $\epsilon_{k}$, $\beta=1/T$, and $\mu$ is the chemical potential. Let the particles be non-relativistic with rest mass $m$ and kinetic energy $E=p^{2}/(2m)$. Also, let the particles have internal states with excitation energies $\Delta_{l}$ (with respect to the ground state $l=0$, so $\Delta_{0}=0$) and spins $J_{l}$, so that the multiplicity factor is $2J_{l}+1$. Thus, a state $k$ is described by the internal state label $l$ and momentum $p$, and the energy is given by

[TABLE]

Recalling that the momentum phase-space volume element is $V/h^{3}d^{3}p$, we find

[TABLE]

where $(2J_{l}+1)$ is the multiplicity of the internal state $l$, we used $\int d^{3}p=4\pi\int_{0}^{\infty}dp\,p^{2}$ and $h=2\pi\hbar=2\pi$ since $\hbar=1$, and we define the internal partition function

[TABLE]

Note that the internal partition function is sometimes given normalized to the ground-state spin factor, i.e.,

[TABLE]

where $J_{0}$ is the ground-state spin of the nuclide and $g(T)$ is a tabulated function (e.g., Rauscher & Thielemann 2000).

The grand potential $\Omega$ is defined as

[TABLE]

and the particle number $N$, pressure $P$, entropy $S$, and internal energy $U$ are given by (e.g., Reichl 1980, §9.B.3)

[TABLE]

We find

[TABLE]

which we can solve for the chemical potential $\mu$ to find

[TABLE]

where $n=N/V$ is the number density. The pressure $P$ becomes

[TABLE]

since $\Omega=-NT$. For the entropy $S$ we obtain

[TABLE]

and since $N=-\Omega/T$, the specific entropy per particle $s=S/N$ is

[TABLE]

where we used Equation A12. Finally, the internal energy per particle is

[TABLE]

A.2. Modified Timmes EOS

In the previous section, we found the most relevant thermodynamic properties of a non-interacting, non-relativistic, non-degenerate Boltzmann gas. In this section, we describe the complete EOS implemented in SkyNet. SkyNet uses a modified Timmes EOS developed in Timmes & Arnett (1999) and Timmes & Swesty (2000). The Timmes EOS consists of three independent parts: a photon gas, an arbitrarily degenerate and relativistic electron/positron gas, and a non-degenerate, non-relativistic Boltzmann gas for the heavy ions (Timmes & Arnett 1999). The electron/positron part is implemented via table interpolation of the Helmholtz free energy (Timmes & Swesty 2000).

For the photon gas and electron/positron gas, the code from the original author of the Timmes EOS, which is available at http://cococubed.asu.edu/codes/eos/helmholtz.tbz, is used. That code also provides the electron/positron chemical potential $\eta_{e^{-},\text{Timmes}}=\mu_{e^{-},\text{Timmes}}/T$, which we need to compute neutrino interactions (Appendix C) and electron screening corrections (§4). Note that the electron/positron chemical potential $\mu_{e^{-},\text{Timmes}}$ in the Timmes EOS is defined with the electron rest mass subtracted out (Timmes & Arnett 1999, §2). The positron chemical potential is

[TABLE]

So, the electron and positron chemical potentials that include the rest masses are

[TABLE]

For the heavy ions, the Timmes EOS implementation uses a single representative ion species with mass $\bar{A}$ and charge $\bar{Z}$, which are the average mass and charge of all nuclides, respectively. Since SkyNet has the complete composition information at all times, we decided to extend the original Timmes EOS to take into account all ion species individually. Furthermore, we use the expressions derived in the previous section for the ion quantities, which take the internal nuclear partition functions into account.

The overall specific entropy of the system is computed in units of $k_{B}\,\text{baryon}^{-1}$ as

[TABLE]

where $s_{e^{\pm},\text{Timmes}}$ and $s_{\gamma,\text{Timmes}}$ are the electron/positron and photon entropies provided by the Timmes EOS, respectively. We divide them by $k_{B}N_{A}$, where $N_{A}\approx 6.022\times 10^{23}\ \text{baryon g}^{-1}$ is the Avogadro constant, because the Timmes EOS returns the entropies in units of $\text{erg g}^{-1}\,\text{K}^{-1}$. The specific entropy of the ions $s_{\text{ions}}$ is calculated by SkyNet itself according to Equation A15 as

[TABLE]

where the sum runs over all nuclear species labeled by $i$, $N_{i}$ is the number of particles of species $i$, $N_{B}$ is the total number of baryons, and $s_{i}$ is the entropy per particle of species $i$ given by Equation A15. Recall that the abundance $Y_{i}$ is (Equation 17)

[TABLE]

where $N_{i}$ and $N_{B}$ are the total number of particles species $i$ and the total number of baryons, respectively, and $V$ is the volume. Thus, the abundance $Y_{i}$ is the fraction of particles of species $i$ compared to the total number of baryons. Note that $N_{i}s_{i}$ is the total entropy contribution of species $i$ and so $N_{i}s_{i}/N_{B}$ is the entropy per baryon contribution of species $i$. Also note that

[TABLE]

where $\rho$ is the mass density.

In §B.2, we require the partial derivative of the entropy with respect to temperature. From Equation A15, we get

[TABLE]

The first two partial derivatives are provided by the Timmes EOS, and from Equation A21 we find

[TABLE]

since $\partial f/\partial\ln T=T\partial f\partial T$. In the current SkyNet implementation, however, we ignore the second derivative of $\ln G_{i}(T)$, because the partition functions we currently have available do not have continuous second derivatives. This will be fixed in a future version of SkyNet.

Similar to Equation A20, the specific internal energy is computed as

[TABLE]

in units of $\text{erg g}^{-1}$. The ion internal energy is computed as

[TABLE]

where $m_{n}$ and $m_{p}$ are the neutron and proton mass, respectively, $Z_{i}$ is the charge number of nuclide $i$, and $\text{BE}_{i}$ is its binding energy. Note that the sum in the above expression gives the internal energy per baryon, so we multiply by $N_{A}$ to convert this to the internal energy per gram. We need to add the electron rest mass because it is not accounted for in $u_{e^{\pm},\text{Timmes}}$ (Timmes & Arnett 1999). The binding energy $\text{BE}_{i}$ is defined as

[TABLE]

with $N_{i}$ being the number of neutrons of nuclide $i$ (not to be confused with the same symbol used above for the total number of particles of species $i$ in the composition) and $m_{i}$ being the rest mass of nuclide $i$. Thus, the term $-\text{BE}_{i}-Z_{i}(m_{n}-m_{p})$ in Equation A27 is

[TABLE]

where $A_{i}=Z_{i}+N_{i}$ is the mass number of nuclide $i$. Thus, Equation A27 differs from the expression for the specific internal energy of a single species (Equation A16) only by the subtraction of $A_{i}m_{n}$ from the particle rest mass $m_{i}$. Thus, the specific internal energy we calculate is relative to the neutron rest mass, which has the advantage that the numerical value of the specific internal energy is not too large but may be comparable to the thermal energy. However, if we ignored the rest mass altogether, the internal energy would not be conserved under nuclear reactions. Nuclear reactions change particles from one species to another that have different binding energies, and so we have to account for the binding energy and mass difference between neutrons and protons, as we do in Equation A27. Other EOSs use the same convention (e.g., Lattimer & Swesty 1991). Note that we have

[TABLE]

where $u_{i}$ is given by Equation A16. But since $\sum_{i}Y_{i}A_{i}=1$ (Equation 52), we find $u_{\text{ions}}=N_{A}\sum_{i}Y_{i}u_{i}-N_{A}m_{n}$, which means our definition of specific internal energy only differs by a constant ($N_{A}m_{n}c^{2}\approx 9.065\times 10^{20}\,\text{erg g}^{-1}$) from the specific internal energy we would calculate by using Equation A16 directly.

Currently, the electron screening corrections implemented in SkyNet (§4) are not yet included in the modified Timmes EOS. Since screening is implemented as a correction to the ion chemical potential (Equation A12), it is not straightforward to propagate those corrections to the other thermodynamic quantities of the ions, let alone the electron/positron gas. We plan to incorporate screening into the EOS in a future version of SkyNet.

Appendix B Calculating NSE

In this section, we show in detail how NSE is computed in SkyNet given a temperature and density (§B.1) and with an unknown temperature (§B.2).

B.1. From temperature and density

The NSE evolution mode requires a robust method for calculating NSE. We have a list of nuclides for which we want to calculate the NSE composition, given a temperature $T$, density $\rho$, and electron fraction $Y_{e}$. Recall that NSE is governed by Equation 31:

[TABLE]

where $\mu_{i}$, $\mu_{p}$, and $\mu_{n}$ are the chemical potentials of the nuclide $i$, protons, and neutrons, respectively. To make the values of the chemical potentials closer to unity, we introduce a renormalized chemical potential $\hat{\mu}$ given by

[TABLE]

where $\mathrm{BE}_{i}$ is the binding energy of the nuclide $i$ defined in Equation A28. Recall that

[TABLE]

where the proton and neutron masses $m_{p}$ and $m_{n}$ are generally chosen such that the binding energies of the neutron and proton are exactly zero. Equation B1 now becomes

[TABLE]

where we used the definition of $\mathrm{BE}_{i}$, Equation B3, and the fact that $\mathrm{BE}_{p}=\mathrm{BE}_{n}=0$.

It is always possible to choose any two chemical potentials as the basis vectors and express all other chemical potentials in terms of those two. In terms of the species $l$ and $m$ that have $Z_{l}$ and $Z_{m}$ protons and $N_{l}$ and $N_{m}$ neutrons, we find

[TABLE]

And so

[TABLE]

where

[TABLE]

Clearly, the species $l$ and $m$ must be chosen to have different proton fractions $Z/N$. However, they do not have to be nuclei in the network. We have had reasonable success with $Z_{l}=0$, $N_{l}=1$ and $Z_{m}=-1$, $N_{m}=1$.

Using Equation A12, the abundance $Y_{i}=n_{i}/n_{B}$ in terms of the chemical potential $\hat{\mu}_{i}$ is given by

[TABLE]

where we define $\eta_{i}=\hat{\mu}_{i}/T$.

This system of equations is subject to baryon number conservation and charge conservation,

[TABLE]

These equations can be zeroed by a two-dimensional NR iteration, where the Jacobian is given by

[TABLE]

so that the chemical potential updates are given by

[TABLE]

From Equations B9 and B10, we find

[TABLE]

and similarly

[TABLE]

The trickiest part about calculating NSE is choosing the basis nuclides $l$ and $m$ and the starting guess for $\bm{\eta}=(\eta_{l},\eta_{m})$. If one of the basis nuclides is $(Z,N)=(-1,1)$, then the corresponding $\eta$ can be set to zero in most cases. To choose the second basis nuclide, compute the chemical potential of the most bound nuclide if it had a mass fraction of 1 and compare this to the chemical potential of the neutron if it had mass fraction 1. The second basis nuclide will be the nuclide corresponding to the larger of those two chemical potentials and the starting guess for that $\eta$ comes from that chemical potential. Although this method may not work in every case, we have found it to work robustly in a large region of parameter space.

B.2. With an unknown temperature

In the previous section, we described how to calculate NSE given a temperature $T$, density $\rho$ (from which we get the baryon number density $n_{B}=N_{A}\rho$), and electron fraction $Y_{e}$. However, there are cases where the temperature is unknown but the entropy $s_{0}$ is given instead. In that case, we have an additional unknown variable $T$, and the additional constraint equation

[TABLE]

where $s_{0}$ is the given target entropy and $s(T,\eta_{l},\eta_{m})$ is the entropy given by the EOS from the current guess for $T$, $\eta_{l}$, and $\eta_{m}$. The Jacobian becomes

[TABLE]

Note that Equation B8 gives

[TABLE]

and so

[TABLE]

Recall that the entropy is calculated as (Equation A20)

[TABLE]

The EOS provides $\partial s/\partial T$ (Equation A24), so

[TABLE]

The electron entropy $s_{e^{\pm}}$ and photon entropy $s_{\gamma}$ only depend on the temperature and electron fraction. Thus, they do not depend on the composition $\bm{Y}$ and thus also not on $\eta_{l}$ and $\eta_{m}$. The ion entropy is calculated as (Equation A21)

[TABLE]

so

[TABLE]

From Equation B8, we get

[TABLE]

Combining the above with Equation B19 yields

[TABLE]

where the first partial derivative is provided by the EOS (Equation A24). We also find

[TABLE]

and similarly

[TABLE]

Unfortunately, the NR iterations with three variables are much less stable than if the temperature is fixed, unless a good initial guess for the temperature is available. For this reason, if NSE is computed from a given entropy and density, SkyNet first uses the bisection method (e.g., Burden et al. 2015, §2.1) to find a good guess for the temperature. The bisection attempts to find a guess temperature such that Equation B17 is close to zero. Then, the NR iterations are performed as described in this section, which may lead to a temperature that satisfies the three constraint equations better than the guess temperature found by the bisection method. However, it can also happen that the bisection already found the best temperature.

If NSE needs to be calculated from a given temperature, entropy, and electron fraction, then the bisection method is used to find the density that produces the desired entropy. Similarly, if the internal energy, density, and electron fraction are given, SkyNet uses the bisection method to find the temperature that produces an NSE distribution with the desired internal energy.

Appendix C Neutrino interaction reactions

The rate for a two-particle charged current weak interaction is given by (see Equation 7)

[TABLE]

where particles [1] and [3] are the incoming and outgoing neutrino and lepton, and particles [2] and [4] are the incoming and outgoing nucleons. Under the astrophysical conditions relevant for reaction networks (see, e.g., Reddy et al. 1998), we have

[TABLE]

where $G_{F}$ is the Fermi coupling constant, $g_{V}$ and $g_{A}$ are the vector and axial vector couplings of the weak current to the nucleons, and the energy-dependent correction factor $(1+h_{wm}E_{3})$ comes from weak magnetism and recoil corrections (Horowitz 2002). Due to the large mass of the nucleons relative to the energy scale of neutrinos emitted from sites undergoing nuclear burning, we can assume there is no momentum transfer from the nucleons to the leptons. Also neglecting final state nucleon blocking, we then have

[TABLE]

where $q_{0}$ is the energy difference between the incoming and outgoing nucleons. The angular integrals in momentum space are trivially integrated. We have $E_{i}dE_{i}=k_{i}dk_{i}$ and the delta function gets rid of the integral over $k_{3}$. This gives

[TABLE]

where $E_{1}=E_{3}-q_{0}$, $k_{1}=\sqrt{E_{1}^{2}-m_{1}^{2}}$, and $\theta(x)=1$ if $x>0$ and 0 otherwise.

SkyNet contains neutrino interactions on free nucleons. Currently, the following reactions are implemented:

[TABLE]

Thus, we compute

[TABLE]

where $Q_{\text{ec}}=-Q_{\text{pc}}=m_{n}-m_{p}=1.29333\,\text{MeV}$, $\omega_{x}=\max(m_{e},Q_{x})$ for $x=\text{ec}$ or $x=\text{pc}$, $f_{e}$ and $f_{p}$ are the electron and positron distribution functions with the electron chemical potential $\mu_{e}$ given by Equation A18, and $f_{\nu_{e}}$ and $f_{\bar{\nu}_{e}}$ are the electron neutrino and antineutrino distributions functions, which are usually Fermi-Dirac with $\mu_{\nu_{e}}=\mu_{\bar{\nu}_{e}}=0$, but could be set to something else, too. We also use the Fermi-Dirac distribution for the electrons and positrons, hence

[TABLE]

where $\beta=1/T$. The factor $C$ is given by (Arcones et al. 2010)

[TABLE]

where the matrix element is $B=g_{V}^{2}+3g_{A}^{2}=5.76$ and $K=6144$ s.

Note that the factor $(E-Q_{x})$ in Equations C9, C10, C11, and C12 is the (anti) neutrino energy. We neglect nucleon recoils so that $E_{\nu_{e}}=E+m_{p}-m_{n}=E-Q_{\text{ec}}$ for electron capture and neutrino absorption, and for positron capture and antineutrino absorption, we get $E_{\bar{\nu}_{e}}=E+m_{n}-m_{p}=E+Q_{\text{ec}}=E-Q_{\text{pc}}$. Therefore, to compute the neutrino heating or cooling rates due to these reactions, we simply multiply the integrand by another factor of $(E-Q_{x})$ and we get

[TABLE]

where the negative sign for the (anti) neutrino absorption reactions comes from the fact that in those reactions, the neutrino energies are absorbed, and they thus provide cooling instead of heating. The total neutrino heating/cooling rate is thus

[TABLE]

The integrals shown in this section are evaluated numerically in SkyNet using the adaptive QAG integration routines provided by the GNU Scientific Library777[https://www.gnu.org/software/gsl/manual/html_node/

Numerical-Integration.html](https://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html). As usual, all of the reaction rates and heating rates need to be multiplied by the product of reactant abundances. So, the electron capture and antineutrino absorption rates are multiplied by $Y_{p}$ and the positron capture and neutrino absorption rates are multiplied by $Y_{n}$.

Alternatively, instead of computing the rates from the integrals provided here, one can also specify the rates $\lambda_{\text{ec}}$, $\lambda_{\text{pc}}$, $\lambda_{\nu_{e}}$, and $\lambda_{\bar{\nu}_{e}}$ directly as a function of time, and SkyNet will use these externally given rates. In that case, the neutrino heating/cooling rate also has to be specified directly. This capability is useful if the (anti) neutrino absorption and emission rates are computed in a hydrodynamical simulation and the same rates should be used for the nucleosynthesis calculations in SkyNet for consistency.