Appendix B: Equation of state and nuclear burning

In this section we give basic information about the equation of state that is currently used in our simulations of stellar core collapse and supernovae.

B.1 Numerical handling of the equation of state

We employ the equation of state (EoS) of Lattimer & Swesty (1991) with a value of $K_{\rm s}=180~{{\rm MeV}}$ for the incompressibility modulus of bulk nuclear matter. The other parameters can be found in the paper of Lattimer & Swesty (1991). Since this EoS assumes matter to be in nuclear statistical equilibrium (NSE) it is applicable only for sufficiently high temperatures and densities. A more general equation of state, which consistently extends to lower temperatures and densities not being available, we have supplemented the EoS of Lattimer & Swesty (1991) with an EoS that describes a mixture of ideal gases of nucleons and nuclei with given abundances, plus ideal gases of electrons and positrons of arbitrary degeneracy and relativity, plus photons (Janka 1999). The latter EoS includes Coulomb lattice corrections for the pressure, energy density, entropy, and adiabatic index.

In general, the regime of NSE is bounded from below by a complicated hyperplane in $(\rho,T,{Y_{\rm e}})$-space. For our purposes, however, it is sufficient to assume that the transition between the two regimes (and thus the two equations of state) takes place at a fixed value of the density $\rho _0$: If the density of a fluid element fulfils the condition $\rho > \rho_0 = 6\times 10^7$  ${{\rm g}~{\rm cm}^{-3}}$, the temperature in a supernova core is usually large enough ( $T\gtrsim 0.5$ MeV) to ensure that NSE holds. Therefore the EoS of Lattimer & Swesty (1991) can be used at such conditions. If the density drops below the value $\rho _0$, the non-NSE equation of state is invoked. Since both EoSs employ a similar description of the physical properties of stellar matter at densities in the range between 10⁷  ${{\rm g}~{\rm cm}^{-3}}$ and 10⁸  ${{\rm g}~{\rm cm}^{-3}}$, the merging of them across the $(T,{Y_{\rm e}})$-plane at $\rho _0$ is sufficiently smooth as far as, e.g., the pressure, internal energy density and chemical potentials as functions of density are concerned.

A particular complication, however, arises with the chemical composition of the stellar plasma. Lattimer & Swesty (1991) describe the baryonic part of the EoS as a mixture of nucleons, $\alpha$-particles and a heavy nucleus with in general non-integer mass and charge numbers (A,Z). The latter is considered to be representative of a mixture of heavy nuclei that coexist at given $(\rho,T,{Y_{\rm e}})$. Assuming NSE, the corresponding number densities of the nuclear constituents and the mass and charge numbers of the representative heavy nucleus are given as a function of the local values of $\rho$, T, and ${Y_{\rm e}}$. At the transition from NSE to non-NSE, the nuclear composition freezes out and one would ideally want to retain the baryonic abundances as given by the EoS of Lattimer & Swesty (1991) at $\rho _0$. While this is no problem in case of a Lagrangian hydrodynamics code, where the grid cells follow the evolution of individual fluid elements with specific information about the chemical composition attributed, the use of an Eulerian or moving grid requires a more complicated numerical procedure. Here, additional advection equations (similar to Eq. (4)) for the different nuclear components must be integrated to trace the temporal evolution of the composition on the whole numerical grid. With a finite, fixed number of such conservation equations to be solved in the code, one cannot allow for an arbitrarily large ensemble of different nuclear species with non-integer mass and charge numbers, as provided at freeze-out from NSE at $\rho _0$ by the EoS of Lattimer & Swesty (1991). Instead, we predefine a discrete set of nuclei $\{(A_k,Z_k)\}_{k=1,\dots,N_{\rm nuc}}$ at the beginning of the simulation, with k=1 for neutrons, k=2 for protons, k=3 for $\alpha$-particles and $k=4,\dots,N_{\rm nuc}$ for a number of suitably chosen heavy nuclei. When the conditions in a computational cell change from NSE to non-NSE, we must map the NSE composition $n_{\rm n}(\rho_0), n_{\rm p}(\rho_0), n_{\alpha}(\rho_0), n_{(A(\rho_0),Z(\rho_0))}(\rho_0)$ as given by the EoS of Lattimer & Swesty (1991) to our discrete sample of species $\{(A_k,Z_k)\}_{k=1,\dots,N_{\rm nuc}}$ with the constraints of charge neutrality and baryon number conservation.

The following procedure turned out to yield very satisfactory results. It is applied at each time step and in cells of the computational grid, which are close to a possible breakdown of NSE (i.e. which are close to $\rho _0$). Applying this procedure also in grid cells where NSE still holds and the EoS of Lattimer & Swesty (1991) is used, makes sure that the compositional information is available for a transition to the non-NSE regime, if the freeze-out condition ( $\rho<\rho_0$) is reached during a hydrodynamic time step. First we identify the densities of neutrons and $\alpha$-particles with their NSE values: $n_{\rm n}=n_{\rm n}(\rho)$, $n_{\alpha}=n_{\alpha}(\rho)$. For given density $\rho$ and electron fraction ${Y_{\rm e}}$, the number densities of protons and nuclei, $n_{\rm p}$and n_(Aj,Zj), respectively, can then be determined from the requirements of charge neutrality and baryon number conservation. The index (A_j,Z_j) points to a particular nucleus $(A_j,Z_j)\in \{(A_k,Z_k)\}_{k=4...N_{\rm nuc}}$, which is chosen from the set of non-NSE nuclei. It is associated with the representative heavy nucleus $(A(\rho),Z(\rho))$ according to the conditions:

 

$$\vert A_j-A\vert \le \vert A_k-A\vert \quad \forall k\in \{4,\dots,N_{\rm nuc}\}~,$$

$$A_j-Z_j \ge A-Z ~,$$ (B.1)

$$Z_j \le Z ~.$$

The first condition ensures that (A_j,Z_j) is selected as the nucleus of the ensemble $(A_j,Z_j)\in \{(A_k,Z_k)\}_{k=4...N_{\rm nuc}}$which is closest to $(A(\rho),Z(\rho))$ concerning its mass number. A unique solution requires that no isobars exist in the ensemble. The second and third conditions guarantee that $n_{\rm p}\ge 0$ and $n_{(A_j,Z_j)}\ge 0$.

Both the EoS of Lattimer & Swesty (1991) and the non-baryonic part of our low-density equation of state are stored in tabular form for our calculations. Our table of the EoS of Lattimer & Swesty (1991) has 180 and 120 logarithmically spaced entries within the density range $5\times 10^7$  ${{\rm g}~{\rm cm}^{-3}}$ $~\le\rho\le 3\times 10^{15}$  ${{\rm g}~{\rm cm}^{-3}}$and the temperature range $0.1~{{\rm MeV}}\le T \le 80~{{\rm MeV}}$, respectively, and 50 equally spaced entries for $0.001\le {Y_{\rm e}}\le 0.6$. Similarly, the table for the lepton plus photon EoS Janka (1999) has 441 entries for 10^-10  ${{\rm g}~{\rm cm}^{-3}}$ $~\le \rho {Y_{\rm e}}\le 10^{12}$  ${{\rm g}~{\rm cm}^{-3}}$and 141 entries for $8.6\times 10^{-6}~{{\rm MeV}}\le T \le 86~{{\rm MeV}}$. Intermediate values are obtained by trilinear interpolation in $\log \rho$, $\log T$, ${Y_{\rm e}}$ and bilinear interpolation in $\log (\rho {Y_{\rm e}})$, $\log T$, respectively.

Using EoSs in a discretized form we have performed the test calculation suggested by Swesty (1996, Sect. 4) for exploring the accuracy of the EoS evaluation. From our tests we can exclude any serious "unphysical entropy production or loss in otherwise adiabatic flows'' (Swesty 1996). For several combinations of initial values for ${Y_{\rm e}}$ and s that are typical of conditions encountered in core-collapse simulations, we found the deviations from adiabaticity to be negligibly small. The relative change of the entropy per baryon was less than 10^-3 for a density increase by a factor of 10⁴.

B.2 Nuclear burning, dissociation and recombination

Nuclear burning, dissociation of nuclei and the recombination of free nucleons and $\alpha$-particles to heavy nuclei can change the baryonic composition in the regime below the transition density ( $\rho_0 = 6\times 10^7$  ${{\rm g}~{\rm cm}^{-3}}$) from the EoS of Lattimer & Swesty (1991) to the low-density EoS. We take into account complete burning of carbon, nitrogen, oxygen, magnesium and silicon as well as nuclear dissociation and recombination in approximate ways.

Burning:

If the temperature in a grid cell exceeds a threshold value characteristic of a certain nuclear reaction, the burning element is assumed to be instantaneously converted to the end product of the reaction. This is a reasonably good approximation because the burning time scale, which is an extremely steep function of the temperature, is small compared with the time scale for hydrodynamic changes. Moreover, in dynamical supernova models one is mainly interested in the energy release associated with nuclear reactions, but not in detailed information about the nucleosynthetic yields.

When the condition for silicon burning applies, which is taken to be $T \ge 4.5\times 10^9$ K (Hix & Thielemann 1999; Mezzacappa et al. 2001), we simply add up the local mass fractions of $\element[][28]{Si}$ and $\element[][56]{Ni}$: $X^{n+1}_{\element[][56]{Ni}}=X^{n}_{\element[][56]{Ni}}+ X^{n}_{\element[][28]{Si}}$, and set $X^{n+1}_{\element[][28]{Si}}=0$. Analogously, we treat the burning of $\element[][12]{C}$ to $\element[][24]{Mg}$ with a threshold temperature of $2.5\times 10^9$ K and of $\element[][16]{O}$, $\element[][20]{Ne}$, $\element[][24]{Mg}$ to $\element[][28]{Si}$ at $T=3.5\times 10^9$ K (Woosley et al. 2002).

Nuclear dissociation and recombination:

When the density drops into the range $\rho<\rho_0$ and the temperature fulfills $T>3\times 10^9$ K at the same time the baryonic composition is adjusted to account for the photo-disintegration of nuclei and the recombination of free nucleons and $\alpha$-particles to $\element[][56]{Ni}$. We distinguish between three different regimes (I-III), separated by curves $\rho_{1}(T), \rho_{2}(T)$ in the $\rho$-T-plane (see Fig. B.1):

  

$$\log(\rho_1(T)):=11.62 + 1.5~\log(T_9)-39.17/T_9~,$$ (B.2)

$$\log(\rho_2(T)):=10.60 + 1.5~\log(T_9)-47.54/T_9~,$$ (B.3)

where the logarithms are to base 10 and T₉:=T/10⁹ K. The function $\rho _{1}(T)$ is derived from the equations of Saha equilibrium assuming a mixture of half (by mass) $\element[][56]{Fe}$ and half $\alpha$-particles and free neutrons.

Figure B.1: Area in the $\rho$-T-plane in the vicinity of the transition between the EoS of Lattimer & Swesty (1991) and the low-density EoS. The transition density $\rho _0$ is indicated by the vertical dash-dotted line. The curves $\rho _{1}(T)$ and $\rho _{2}(T)$ bound three separate regions (I-III) in the $\rho$-T-plane, where different prescriptions are used for determining the nuclear composition in the low-density regime. Region I contains heavy nuclei and possibly free nucleons, in region II the composition is determined by $\alpha$-particles and free nucleons, and in region III all nuclei and $\alpha$-particles are dissolved. Below the temperature marked by the solid horizontal line the recombination of $\alpha$-particles and free nucleons is suppressed. The threshold temperatures of the different burning reactions are given by the horizontal dotted lines.

Analogously, $\rho _{2}(T)$ is defined by half the mass being in $\alpha$-particles and half in free nucleons (see Shapiro & Teukolsky 1983, Chap. 18). In detail the nuclear composition in the three regions is obtained from the following prescriptions:

I)

If $\rho > \rho_1(T)$, all free nucleons and a fraction $0 \le f_{\rm I} \le 1$ of possibly existing $\alpha$-particles in a corresponding grid cell are assumed to form $\element[][56]{Ni}$. Pairing a suitable number of free neutrons and protons we therefore set for the new composition:

 

$$X^{n+1}_{\rm n} \max\{X^{n}_{\rm n}-X^{n}_{\rm p},0\} ~,$$ =

$$X^{n+1}_{\rm p} \max\{X^{n}_{\rm p}-X^{n}_{\rm n},0\} ~,$$ =

$$X^{n+1}_{\element[][4]{He}} (1-f_{\rm I})\cdot X^{n}_{\element[][4]{He}} ~,$$ =

$$X^{n+1}_{\element[][56]{Ni}} X^{n}_{\element[][56]{Ni}} +f_{\rm I}\cdot X^{n}_{\element[][4]{He}} +2\cdot\min\{X^{n}_{\rm n},X^{n}_{\rm p}\} ~.$$ = (B.4)

II)

If $\rho_2(T) \le \rho \le \rho_1(T)$, $\alpha$-particles are formed by the disintegration of heavy nuclei and the recombination of free nucleons. The new composition is given by:

$$X^{n+1}_{\rm n} \max\{X^{*}_{\rm n}-X^{*}_{\rm p},0\} ~,$$ =

$$X^{n+1}_{\rm p} \max\{X^{*}_{\rm p}-X^{*}_{\rm n},0\} ~,$$ =

$$X^{n+1}_{\element[][4]{He}} X^{*}_{\element[][4]{He}}+2\cdot\min\{X^{*}_{\rm n},X^{*}_{\rm p}\}~,$$ =

$$(1-f_{\rm II})\cdot X^{n}_{k}\quad \forall k\in \{4,\dots,N_{\rm nuc}\} ~,$$ Xⁿ⁺¹_k =

with

 

$$X^{*}_{\rm n} X^{n}_{\rm n}+f_{\rm II}\cdot\sum_{k=4}^{N_{\rm nuc}} \frac{\max\{A_k-2 Z_k,0\}}{A_k}\cdot X^{n}_{k} ~,$$ =

$$X^{*}_{\rm p} X^{n}_{\rm p}+f_{\rm II}\cdot\sum_{k=4}^{N_{\rm nuc}} \frac{\max\{2 Z_k-A_k,0\}}{A_k}\cdot X^{n}_{k} ~,$$ =

$$X^{*}_{\element[][4]{He}} X^{n}_{\element[][4]{He}}+f_{\rm II}\cdot 2\sum_{k=4}^{N_{\rm nuc}}\frac{\min\{A_k-Z_k,Z_k\}}{A_k}\cdot X^{n}_{k} ~.$$ = (B.5)

The factor $0 \le f_{\rm II} \le 1$ is the degree of dissociation of heavy nuclei which is determined as outlined below.

III)

If $\rho < \rho_2(T)$, all heavy nuclei and a fraction $0 \le f_{\rm III} \le 1$ of the $\alpha$-particles is disintegrated into free nucleons:

 

$$X^{n+1}_{\rm n} X^{n}_{\rm n} +f_{\rm III}\cdot\frac{1}{2}X^{n}_{\element[][4]{He}} +\sum_{k=4}^{N_{\rm nuc}} \frac{A_k-Z_k}{A_k}\cdot X^{n}_{k} ~,$$ =

$$X^{n+1}_{\rm p} X^{n}_{\rm p} +f_{\rm III}\cdot\frac{1}{2}X^{n}_{\element[][4]{He}} +\sum_{k=4}^{N_{\rm nuc}} \frac{Z_k}{A_k}\cdot X^{n}_{k} ~,$$ =

$$X^{n+1}_{\element[][4]{He}} (1-f_{\rm III})\cdot X^{n}_{\element[][4]{He}}~,$$ =

$$0 \quad \forall k\in \{4,\dots,N_{\rm nuc}\} ~.$$ Xⁿ⁺¹_k = (B.6)

Nuclear transmutations release (or consume) nuclear binding energy. Consequences of this effect should be taken into account in the hydrodynamical simulation. Since our EoS defines the internal energy density such that contributions from nuclear rest masses are included, the conversion between nuclear binding energy and thermal energy will automatically lead to corresponding changes of the temperature, pressure, entropy etc. of the stellar gas.

The factors $f_{\rm I}$, $f_{\rm II}$ and $f_{\rm III}$ are equal to unity in the corresponding regions I, II or III, respectively of the $\rho$-T-plane but vary between 0 and 1 at the boundary curves (Eqs. (B.2), (B.3)). This means that the composition changes take place only on the curves $\rho _{1}(T)$ and $\rho _{2}(T)$. The degree of dissociation or recombination (expressed by the values $0 \le f_{\rm I}, f_{\rm II},f_{\rm III}\le 1$) is limited by the available internal energy. Starting out with given internal energy density (at a fixed value of $\rho$) and taking into account temperature variations due to changes of the nuclear composition, the numerical algorithm maximizes $f_{\rm I}$, $f_{\rm II}$ or $f_{\rm III}$, while keeping the temperature on the boundary curves until the dissociation or recombination is complete.