Do not forget the electrons: Extending moderately-sized nuclear networks for multidimensional hydrodynamic codes

Abstract

Context. Nuclear networks are widely used coupled with hydrodynamical simulations of explosive scenarios to account for the change of nuclear species and energy generation rate due to nuclear reactions. In this way, there is a feedback mechanism between the hydrodynamical state and the nuclear processes. Unfortunately, the timescale of nuclear reactions is orders of magnitude smaller than the dynamical timescale that drives hydrodynamical simulations. Therefore, these nuclear networks are usually very small, reduced in most cases to a dozen elements, especially when simulations are carried out in more than one dimension.

Aims. We present here an extended nuclear network, with 90 species, designed for being coupled with hydrodynamic simulations, which includes neutrons, protons, electrons, positrons, and the corresponding neutrino and anti-neutrino emission. This network is also coupled with temperature, making it extremely robust and, together with its size, unique of its kind. The inclusion of electron captures on free protons makes the network very appropriate for multidimensional studies of Type Ia supernova explosions, especially when the exploding object is a massive white dwarf.

Methods. We perform several tests that are relevant to simulate explosive scenarios, such as Type Ia supernovae and core-collapse supernovae. We compare the results of the 90 nuclei network with a standard $\alpha$-chain network with 14 elements to evaluate the differences in the energy generation rate. We also evaluate the relevance of including the electrons in the network in terms of generated yields and how it affects the pressure of a degenerate fluid such as that of white dwarfs. The results obtained with the 90-nuclei network have been verified with a much larger 2000-nuclei network built from REACLIB (WinNet), in terms of nuclear energy generation rate, pressure, and produced yields.

Results. The results obtained with the proposed medium-sized network compare fairly well, to a few percent, with those computed with WinNet in scenarios reproducing the gross physical conditions of current Type Ia supernova explosion models. In those cases where the carbon and oxygen fuel ignites at high density, the high-temperature plateau typical of the nuclear statistical equilibrium regime is well defined and stable, allowing large integration time steps. We show that the inclusion of electron captures on free protons substantially improves the estimation of the electron fraction of the mixture. Therefore, the pressure is better determined than in networks where electron captures are excluded, which will ultimately lead to more reliable hydrodynamic models. Explosive combustion of helium at low density, occurring near the surface layer of a white dwarf, is also better described with the proposed network, which gives nuclear energy generation rates much closer to WinNet than typical reduced alpha networks.

Conclusions. A nuclear network with N=90 species, including electrons, aimed at multidimensional calculations of supernova explosions is described and verified. The proposed network is suitable for the study of Type Ia supernova explosions because it provides better values of pressure and electron abundance than other existing networks with smaller or even a similar size but without including electron capture processes.

1 Introduction

Understanding cosmic violent phenomena such as novae or both types of supernova explosions demands not only a good knowledge of nuclear reaction cross sections but a careful implementation of the nuclear networks in the hydrodynamic codes aimed at simulating these explosions. This issue is especially relevant in the case of Type Ia supernova explosions (SNIa), where the disrupting mechanism is of thermonuclear origin and strong feedback between the hydrodynamics and the released nuclear energy is expected (Nomoto et al. 1984; Hillebrandt & Niemeyer 2000).

The question of the optimal size of these nuclear networks has always challenged the community and is, among other things, closely connected to the available computational resources. Nowadays, hydrodynamic codes that simulate SNIa explosions are, for the most part, of multidimensional nature. The number of spatial and temporal steps involved in the simulation of one of these explosions is huge, which poses practical limits to the size of nuclear networks. Pioneering multidimensional simulations of SNIa incorporated around 10 nuclei in an attempt to have, at best, a reasonable depiction of the released nuclear energy. For example, Benz et al. (1989) implemented an $\alpha$ network to keep track of the nuclear evolution during the collision of two white dwarfs. Timmes et al. (2000) analyzed the performance of two networks, with 7 and 14 species, and concluded that these can account for the nuclear energy generation rate within a $20\%$ precision on average. The use of these small networks has been commonplace in multidimensional simulations of supernova explosions around the edge of the last century, as Reinecke et al. (1999); García-Senz et al. (1999); Plewa et al. (2004), to mention a few connected with different SNIa explosion scenarios. Currently, multidimensional simulations incorporate several dozens of nuclei, which can reproduce the released nuclear energy within a narrow deviation, typically a few percent on average when compared to a larger network (Townsley et al. 2019; Gronow et al. 2021).

However, with current computing technology and numerical methods, it is nowadays feasible to extend the number of nuclei in the network to around one hundred (e.g., Harris et al. 2017; Sandoval et al. 2021; Navó et al. 2023). In that case, the matching of the released nuclear energy and other magnitudes compared to the post-processed values would benefit further. Nevertheless, the post-processing step would still be necessary to obtain nucleosynthetic details of the produced yields.

In this work, we propose a network of 90 species, net90, which can be used to model all kinds of SNIa explosions. As we will see, it can be used to assist in the explosion of a massive white dwarf (WD) made of ¹²C$+^{16}$O with central density $\rho_{c}\simeq 2\times 10^{9}$  g$\cdot$cm^-3 in the so-called Chandrasekhar-mass explosion models. It is also adequate to track the detonation of a tiny helium shell at densities $\simeq 5\times 10^{5}$  g$\cdot$cm^-3 located on top of a moderately massive WD, which is appropriate for studying the Sub-Chandrasekhar mass route to SNIa. net90 also improves over existing reduced, $\alpha$-like, networks regarding core collapse supernova (CCSN) scenarios, but its use is restricted to the stages with low or moderate neutronization.

Our proposal has two important features worth mentioning. First, the abundances of the species $Y_{i}=X_{i}/A_{i}$ (where $X_{i}$ is the mass fraction and $A_{i}$ the atomic mass of the species $i$) and the temperature $T$ are found jointly, after implicitly solving the system of nuclear reactions coupled with the energy equation. It is well known that the implicit coupling of $\{Y_{i},T\}$ leads to a very stable behavior around the nuclear statistical equilibrium (NSE) regime; Mueller (1986), Cabezón et al. (2004) (paper I afterward). The implicit coupling between released nuclear energy, nuclear reaction rates, and temperature makes the scheme robust in simulating high-density combustion, allowing it to naturally handle the quasi (QNSE) and complete nuclear statistical equilibrium (NSE) stages smoothly, without switching to specific NSE routines. It also makes unnecessary the implementation of burning limiters for the simulation of detonation waves (Zingale et al. 2021).

The second novelty of this work is that our moderately-sized network includes the capture reactions $e^{-}+p\rightarrow n+\nu_{e}$ and $e^{+}+n\rightarrow p+\bar{\nu}$. As far as we know, this is the first time that reactions of this kind have been taken into account in a reduced network routine belonging to a multi-D hydrodynamic code dealing with SNIa explosions. As we show below, the benefits of including this reaction, namely a more accurate depiction of the electron pressure and neutronized yields, largely compensate for the shortcoming of a tiny increase ($<0.5\%$) in computing time.

In Sect. 2 we summarize the mathematical formalism used to integrate the chemical and temperature equations. The main features of the different nuclear networks used in this work are described in Sect. 3. Section 4 is devoted to verifying net90 and comparing it with a much larger nuclear network in several scenarios. The computational performance of net90 and related issues are discussed in Sect. 5. Finally, we provide a summary and discussion of the results and prospects for the future in Sect. 6.

2 Mathematical formalism

We follow the method described in paper I and in Sanz et al. (2022) to integrate the nuclear equations. We summarize here the main features of the method.

The abundances $Y_{i}$ of the species evolve according to the set of ordinary differential equations,

$$\frac{dY_{i}}{dt}=\sum_{k,l}{r_{kl}Y_{l}Y_{k}}-\sum_{j}{r_{ij}Y_{i}Y_{j}}+\sum_{m}{\lambda_{m}Y_{m}}-\lambda_{i}Y_{i}\,.$$ (1)

The magnitudes $r_{ij}=r_{ij}(\rho,T)=\rho N_{A}\langle\sigma,\nu\rangle_{ij}$ correspond to the nuclear reaction rates, and $\lambda_{i}=\lambda_{i}(T)$ stands for $\gamma$ photodisintegration rates, $\beta^{+}$ decays, and/or $e^{\pm}$ capture processes.

Following paper I, the energy equation can be written as,

$$dQ +\sum_{i}{\left({BE}_{i}-\frac{\partial U}{\partial Y_{i}}\right)dY_{i}}+\left(B_{pn}-\frac{\partial U}{\partial Y_{e}}\right)dY_{e}-\left(\frac{\partial U}{\partial T}-\frac{d\rho}{\rho^{2}}\frac{\partial P}{\partial T}\right)dT$$ (2)

$$=-T\frac{d\rho}{\rho^{2}}\left(\frac{\partial P}{\partial T}\right)\,,$$

where $BE$ is the nuclear binding energy of the nuclei and, in this work, $dQ\equiv-\dot{\nu}\leavevmode\nobreak\ dt$ accounts for the neutrino losses, $U$ and $P$ are the internal energy and pressure of the gas, and $B_{pn}=0.782$ MeV is the proton+electron to neutron mass-difference. This set of equations is usually integrated keeping $\{\rho,T\}$ constant during the current timestep $\Delta t$. Once the new abundances are known, the nuclear energy released in $\Delta t$ is obtained and added to the energy equation. However, this procedure leads to problems during explosive high-temperature combustion, where the rates of direct and inverse reactions are very large, yet they may balance. Close to or at NSE the integration becomes unstable, demanding very small timesteps or even making the calculation non-feasible. However, the joint implicit integration of the abundance and energy equations has been shown to restore the tractability of the problem (Mueller 1986).

We first write the set of chemical equations (Eq. 1) plus the energy equation (Eq. 2) in a suitable form to obtain an implicit solution via a Newton-Raphson (NR) scheme. For example, Eq. 1 is written as,

$$\begin{split}\frac{{Y_{i}}^{n+1}-{Y_{i}}^{n}}{\Delta t}&=\sum_{k,l}{r_{kl}(T^{n+\theta})Y_{l}^{n+\theta}Y_{k}^{n+\theta}}-\sum_{j}{r_{ij}(T^{n+\theta})Y_{i}^{n+\theta}Y_{j}^{n+\theta}}\\ &+\sum_{m}{\lambda_{m}(T^{n+\theta})Y_{m}^{n+\theta}}-\lambda_{i}(T^{n+\theta})Y_{i}^{n+\theta}\,,\end{split}$$ (3)

where the notation $Y_{i}^{n+\theta}=Y_{i}^{n}+\theta\Delta Y_{i}$ and $T^{n+\theta}=T^{n}+\theta\Delta T$ is used to choose among explicit, implicit, or mixed schemes by varying the value of the parameter $\theta\leavevmode\nobreak\ [0,1]$. For convenience, we use the simplified notation $r_{ij}^{n+\theta}\equiv r_{ij}(T^{n+\theta})$ and $\lambda_{i}^{n+\theta}\equiv\lambda_{i}(T^{n+\theta})$, hereafter. Furthermore, we perform a first-order Taylor expansion, $r_{ij}^{n+\theta}\simeq r_{ij}^{n}+\theta\leavevmode\nobreak\ r{{}^{\prime}}^{n}_{ij}\Delta T$, where $r^{\prime}$ is the derivative of the rate with respect to the temperature (and also for $\lambda^{n+\theta}_{i}$). We apply the NR iterative scheme to solve the nonlinear system of Eq. 3 together with the energy equation (Eq. 2). At the integration step $n+1$ and NR $k-$iteration, the corrections $\delta Y_{i}^{n+1}$ to $(Y_{i}^{n+1})^{k-1}$ and $\delta T^{n+1}$ to $(T^{n+1})^{k-1}$ are found from,

Figure 1: Sketch of the architecture of net90. There are shown, in blue, the three additional nuclei with respect net86 plus electrons (in green)

$$\begin{split}&\frac{\delta Y_{i}^{n+1}}{\Delta t}=\theta\sum_{k,l}r_{kl}^{n+\theta}Y_{l}^{n+\theta}\delta Y_{k}^{n+1}+\theta\sum_{k,l}{r_{kl}^{n+\theta}Y_{k}^{n+\theta}\delta Y_{l}^{n+1}}\\ &-\theta\left(\sum_{j}{{r_{ij}^{n+\theta}Y_{j}^{n+\theta}}}\right)\delta Y_{i}^{n+1}-\theta\sum_{j}{{r_{ij}^{n+\theta}Y_{i}^{n+\theta}}\delta Y_{j}^{n+1}}+\theta\sum_{m}\lambda_{m}^{n+\theta}\delta Y_{m}^{n+1}\\ &-\theta\lambda_{i}^{n+\theta}\delta Y_{i}^{n+1}+\theta\leavevmode\nobreak\ \Bigg{(}\sum_{k,l}{r^{\prime}}_{kl}^{n+\theta}Y_{l}^{n+\theta}Y_{k}^{n+\theta}-\sum_{j}{r^{\prime}}_{ij}^{n+\theta}Y_{i}^{n+\theta}Y_{j}^{n+\theta}\\ &+\sum_{m}{{\lambda{{}^{\prime}}}_{m}^{n+\theta}Y_{m}^{n+\theta}}+{\lambda{{}^{\prime}}}_{i}^{n+\theta}Y_{i}^{n+\theta}\Bigg{)}\leavevmode\nobreak\ \delta T^{n+1}+\sum_{k,l}{r_{kl}^{n+\theta}Y_{l}^{n+\theta}Y_{k}^{n+\theta}}\\ &-\sum_{j}{r_{ij}^{n+\theta}Y_{i}^{n+\theta}Y_{j}^{n+\theta}}+\sum_{m}{\lambda_{m}^{n+\theta}Y_{m}^{n+\theta}}-\lambda_{i}^{n+\theta}Y_{i}^{n+\theta}-\frac{Y_{i}^{n+1}-Y_{i}^{n}}{\Delta t}\,,\end{split}$$ (4)

for the chemical equations, and from,

$$\begin{split}&\sum_{k}\left(BE_{k}-\frac{\partial U_{ions}}{\partial Y_{k}}\right)\delta Y_{k}^{n+1}+\left(B_{pn}-\frac{\partial U_{e}}{\partial{Y_{e}}}-\frac{\partial(\dot{\nu}\Delta t)}{\partial{Y_{e}}}\right)\delta Y_{e}^{n+1}\\ &-\left(\frac{\partial U}{\partial T}-\frac{\Delta\rho}{\rho^{2}}\frac{\partial P}{\partial T}+\frac{\partial(\dot{\nu}\Delta t)}{\partial T}\right)\delta T^{n+1}=\left(\frac{\partial U}{\partial T}-\frac{\Delta\rho}{\rho^{2}}\frac{\partial P}{\partial T}\right)(T^{n+1}-T^{n})\\ &-\sum_{k}\left(BE_{k}-\frac{\partial U_{ions}}{\partial Y_{k}}\right)(Y_{k}^{n+1}-Y_{k}^{n})-\left(B_{pn}-\frac{\partial U_{e}}{\partial Y_{e}}\right)(Y_{e}^{n+1}-Y_{e}^{n})\\ &+\dot{\nu}\Delta t-T\frac{\Delta\rho}{\rho^{2}}\left(\frac{\partial P}{\partial T}\right)\,,\end{split}$$ (5)

for the energy equation. Equation 5 includes the neutrino losses, $\dot{\nu}\Delta t$, produced by electron and positron captures, and the remaining symbols have the usual meaning.

Equations 4 and 5 conform to an $(N+1)\times(N+1)$ linear system of equations with $N$ unknown chemical corrections, $\delta Y_{k}$, plus the temperature correction $\delta T$. The system is sparse and is efficiently solved with an improved version of the algorithm described in Prantzos et al. (1987). As an initial guess in the first NR iteration, we take $Y^{n+1}=Y^{n};T^{n+1}=T^{n}$ and the parameter $\theta$ in Eq. 4 was set to $\theta=0.7$ in the tests described in Sect. 4. This choice of $\theta$ was empirically chosen to balance between $\theta=1$ providing optimal stability and $\theta=0.5$ of perfectly centered integration schemes.

3 Nuclear networks

In this section, we introduce a family of nuclear networks, which we named netX, where $X$ indicates the number of species. From simpler to the more complex, net14, net86, and the new net90 presented here. All of these are based on REACLIB. Verification of net90 is performed by comparing its results with those obtained with the state-of-the-art nuclear network WinNet (Reichert et al. 2023), a much larger network that includes many weak interaction processes. The main features of all these nuclear networks are described below.

3.1 net14

This network was introduced in García-Senz & Cabezón Gómez (2003) and it is representative of typical $\alpha$-networks, widely used coupled to hydrodynamic codes. This network follows an alpha-chain of reactions from ⁴He up to ⁶⁰Zn. Many similar networks had either fewer nuclear species such as Timmes et al. (2000) with only seven nuclei or a bit more, as in Weaver et al. (1978) with 19 isotopes. However, the main characteristic of net14 is that it solves the equations for the evolution of nuclear species coupled with the temperature equation. This methodology was developed with net14 and used later in the other two networks, net86 and net90. Small networks such as net14 are suitable to follow the energy generation rate in explosive scenarios such as Type Ia supernovae, having reasonable feedback with the hydrodynamic state of the fluid. Their accuracy in the energy generation rate is expected to be of the order of $\simeq 20\%$ (Timmes et al. 2000) and their nucleosynthetic yields must be post-processed with larger networks in the range of thousands of nuclei. Despite this, their popularity is due to their reduced size, which results in a very small computational footprint, even when they have to be called many times per fluid element and per hydrodynamic timestep.

3.2 net86

This network was born as an evolution of net14, with the objective of including protons and neutrons, and was also described in Cabezón et al. (2004). To do so, many more nuclear species must be included, and if not careful, the size of the network can increase extremely fast. To that extent, we built this network with a set of blocks of neutron and proton captures around the $\alpha$-chain, again from ⁴He up to ⁶⁰Zn. The result is an extended network with 86 nuclei and $\sim 150$ reactions. This is a reasonable compromise between accuracy (for energy generation rate and yields) and computational burden. net86 also solves the abundance equations jointly with the temperature equation, which proved to be extremely robust, something that becomes even more important when we increase the number of species.

Figure 2: Comparison among $e^{-}$ capture rates at temperatures $T_{9}=1$, $3$, $10$, and different values of $\rho Y_{e}$. The rate obtained numerically with the exact formulation of FFN (solid lines, Fuller et al. 1982) is compared to those obtained with the analytical approach by Hansen (1968) (dashed lines) and the zero temperature limit rate by Bludman et al. (1982) (black dotted line).

Figure 3: Comparison between capture rates $\lambda_{e^{-}p}$ and $\lambda_{e^{+}p}$ at temperatures $T_{9}=1$, $3$, $10$, and different values of $\rho Y_{e}$. We note that, out of the low-$\rho$ high-T region, $\lambda_{e^{+}n}<<\lambda_{e^{-}p}$.

Figure 4: Neutrino and anti-neutrino emissivities (in MeV/s) at temperatures $T_{9}=1$, 3, and 10.

Figure 5: Results of the CO test. Evolution of nuclear mass fractions obtained with net14 (top left), net90-noec (bottom left, without $e^{-}$ captures), and net90 (bottom right, with $e^{-}$ captures on protons included). The vertical dotted lines show the approximate limits of the NSE, when the artificial adiabatic expansion starts, and the region where nuclear reactions are quenched. The panel in the upper right shows the evolution of temperature and pressure for net90-noec and net90. The temperature obtained with net14 is also shown as a reference. The evolution of density is the same for all three nuclear networks.

Figure 6: Results of the CO test. Mass fractions in Fig. 5, but grouped in three nuclear families: Intermediate-Mass $\alpha$-Elements (IME, ²⁸Si$+^{32}$S$+^{36}$Ar$+^{40}$Ca$+^{44}$Ti), $\alpha$-Iron nuclei (⁴⁸Cr$+^{52}$Fe$+^{56}$Ni$+^{60}$Zn, labelled as ⁵⁶Ni), and neutronized Iron-Group Elements (IGE, ⁵⁴Fe$+^{58}$Ni$+^{57}$Co$+^{58}$Co. Solid lines are for net90 and dashed lines are for net90-noec. We also show two lines calculated with net14 ($\alpha$-network).

3.3 net90

net90 is an extension of net86 to include electron and positron captures plus three new nuclei and four additional reactions: ¹²C($p,\gamma)^{13}$N, ¹³N($\alpha,p)^{16}$O, ⁵⁴Fe($\alpha,p)^{57}$Co, and ⁵⁷Co($n,\gamma)^{58}$Co. The first two have fast reaction rates and become relevant during the explosive combustion of He and the other two work as sinks for $\alpha$ particles and neutrons. This makes net90 a suitable option for handling moderately asymmetric matter in the context of multidimensional simulations. A sketch of the network is shown in Fig. 1. In the following, we will specify ’net90-noec’ when referring to ’net90 with $e^{\pm}$ captures turned off’, which in turn is very similar to net86.

In scenarios involving objects such as white dwarfs, where the major contribution to the fluid pressure comes from degenerate electrons, including the electrons in the reduced network is extremely important, and it is something that, as far as we know, has never been done before in the framework of multidimensional hydrodynamic simulations. Indeed, having a process that removes electrons from the fluid, such as capture by protons, can have a significant impact on 3D simulations of Type Ia supernovae. Many reactions involve electrons, but we chose to include only capture by protons, as it is expected to be the dominant process in SNIa (Thielemann et al. 1986) and we can benefit from the structure of net86 without major changes. Furthermore, and according to Brachwitz et al. (2000) and Thielemann et al. (2004) the $e^{-}$ captures on nuclei are clearly subdominant at densities $\rho_{9}\simeq 2$, especially since the weak rates compilation of Langanke & Martínez-Pinedo (2001) significantly reduced the $e^{-}$ capture rates on nuclei of previous compilations (Fuller et al. 1982). We showed and confirmed in the present work (see Sect.4.1.1) that the $e^{-}$ captures on free protons alone could amount up to $\simeq 90\%$ of the captures provided that the nuclear network is large enough (f.e. WinNet).

The reaction rate $e^{-}+p\rightarrow n+\nu$, in a degenerate fluid, can be approached analytically (Hansen 1968; Bludman et al. 1982) or numerically (Fuller et al. 1980). The main advantage of an analytical expression is that its evaluation and its derivative are straightforward, whereas the numerical calculation, although more accurate, requires interpolating from a large data table. We have compared both methods and finally decided to rely on the numerical approach. We show in Figure 2 the $e^{-}$ capture rate, $\lambda_{e^{-}p}$, as a function of $(\rho Y_{e})$ for three fiducial temperatures: $T_{9}=1$, 3, and 10. As we can see, below $\rho Y_{e}\simeq 10^{9}$  g$\cdot$cm^-3 the analytical expressions progressively diverge from the exact numerical value. It should be noted that the expression by Hansen (1968) works well at low density and high temperature, (see f.e. the region $\rho Y_{e}\leq 10^{7},T_{9}=10$), the region where $e^{\pm}$ pair production is relevant.

At temperatures $T_{9}\simeq 10$ and sufficiently low densities, the rate $\lambda_{e^{+}n}$ of the reaction $e^{+}+n\rightarrow p+\bar{\nu}$ competes with $\lambda_{e^{-}p}$. We show this in Fig. 3 which was obtained by numerically integrating the equations of Fuller et al. (1980). Although we expect that in a typical SNIa explosion sequence $\lambda_{e^{+}n}<\lambda_{e^{-}p}$, due to the reduction of $T^{max}$ at lower densities (for example, $T_{9}^{max}\simeq 6,10$ at $\rho\simeq 2\times 10^{7}$ and $10^{9}$ g$\cdot$cm^-3, respectively), the contribution of $\lambda_{e^{+}n}$ could be significant around $T^{max}$ and we therefore included it in the network.

We assume that neutrinos and antineutrinos created during reactions $(e^{-},p)$ and $(e^{+},n)$ freely escape the system, without further interactions. We computed the $\{\nu,\bar{\nu}\}$ emissivities with the scheme developed by Fuller et al. (1980) and stored the values in an interpolating table. Figure 4 shows the expected $\{\nu,\bar{\nu}\}$ luminosity in the range $10^{6}$ g$\cdot$cm^-3 $\leq\rho\leq 10^{10}$ g$\cdot$cm^-3.

Table 1: Initial conditions for all tests performed with the nuclear networks presented in this work.

Test name	$X_{0}$	$T_{9}$	$\rho_{9}$	$\eta_{1}$	$\eta_{2}$	Comment
CO	$\left[{}^{12}\textrm{C}\right]=0.5$	1	2	50	2	WD burning of initially symmetric matter
	$\left[{}^{16}\textrm{O}\right]=0.5$
He	$\left[{}^{4}\textrm{He}\right]=1$	1	$5\times 10^{-4}$	2	1	Shell burning on the WD surface
Si	$\left[{}^{28}\textrm{Si}\right]=1$	6	0.01	2	1	Photodisintegration leading to CCSN

Figure 7: Results of the CO test. Top left: Evolution of the nuclear energy rate as a function of the temperature obtained with net14, net90-noec, and WinNet-noec. Top right: Same, but including $e^{-}$ captures. Letters A, B, C, and D refer to normal combustion, relaxation to NSE, full NSE, and adiabatic expansion, respectively. Bottom left: Evolution of the mean molecular weight $\mu_{i}$ of ions as a function of temperature obtained with net14, net90-noec and WinNet-noec. Bottom right: Same, but including $e^{-}$ captures.

Figure 8: Time-evolution of the nuclear generation rate in the CO test, computed without (left) and with (right) electron captures.

3.4 WinNet

The nuclear reaction network WinNet is presented in Reichert et al. (2023). Unlike the networks presented previously in this work, WinNet contains a flexible number of nuclei. For all calculations, we include 2,281 nuclei up to $Z=50$. This includes all relevant nuclei for the test cases considered. We assume nuclear statistical equilibrium for temperatures higher than 10 GK. Reaction rates have been taken from the REACLIB database (Cyburt et al. 2010). We apply electron screening corrections for all reactions as described in Kravchuk & Yakovlev (2014). Electron and positron capture rates, as well as $\beta^{+}$- and $\beta^{-}$-decays that are contained in the REACLIB database, are replaced by those of Fuller et al. (1985); Oda et al. (1994); Langanke & Martínez-Pinedo (2001); Pruet & Fuller (2003); Suzuki et al. (2016). Notably, this also includes the electron- and positron-capture rates on free protons and neutrons that are taken from Langanke & Martínez-Pinedo (2001). In total, WinNet contains $\sim$30,000 reaction rates, of which $\sim$4,000 reactions are weak reactions. Information on the energy of emitted neutrinos was derived from the same tabulation as the theoretical weak rates (Fuller et al. 1985; Oda et al. 1994; Langanke & Martínez-Pinedo 2001; Pruet & Fuller 2003; Suzuki et al. 2016) and from the ENSDF database (Brown et al. 2018), when it was unavailable in previous tables. Also, we consider the energy of thermal neutrinos (Itoh et al. 1996).

4 Tests

In this section, we present a suite of tests in which, given an initial composition, temperature, and density, we iteratively apply our nuclear networks with a self-adaptive timestep until NSE or a stable configuration has been achieved isochorically. This usually happens faster than the hydrodynamic timescale, defined as:

$$t_{hd}=446/\sqrt{\rho_{0}},$$ (6)

where $\rho_{0}$ is the initial density value (corresponding to $\rho_{9}$ in Table 4). Once the physical time is larger than $t=\eta_{1}t_{hd}$, where $\eta_{1}$ is a constant that depends on the particular scenario, we impose an artificial exponential decrease in density, mimicking an adiabatic expansion of the fluid element in which the nuclear reactions are taking place,

$$\rho(t)=\rho_{0}\exp{\left(-\frac{\Delta t}{\eta_{2}t_{hd}}\right)},$$ (7)

where $\Delta t$ is the current timestep and the choice of $\eta_{2}$ depends on the particular calculation. This produces a cooling that naturally leads to the freezeout of the nuclear reactions. In Table 4 we summarize the initial conditions of all the tests presented in this section. We used the EOS of Timmes & Swesty (2000) to monitor the temperature changes and pressure in all tests.

Figure 9: Results of the CO test. Top left: Evolution of the mass fraction $X_{\alpha}$ of alpha particles of net14, net90-noec, and WinNet-noec. Top right: Same, but for neutron and proton mass fractions. Bottom left: Evolution of the mass fraction $X_{\alpha}$ of the alpha particles of net90 and WinNet. Bottom right: Same, but for proton and neutron mass fractions.

4.1 C+O burning

This test is characteristic of a Chandrasekhar mass white dwarf core that undergoes a nuclear runaway that leads to a type Ia supernova explosion. We used an equal amount of ¹²C and ¹⁶O to mimic the gross conditions prior to the explosion. The initial conditions for this test are summarized in Table 4. Under these conditions, NSE is achieved rapidly when the temperature increases and the radiative captures of protons, neutrons, and alpha particles are balanced by photodisintegrations. The chosen value of $\eta_{1}$ in Table 4 is compatible with the duration of the NSE combustion regime in the center of a massive white dwarf that results in explosive ignition of carbon.

Figure 5 presents the results of the calculation. The top-left panel shows the evolution of the alpha-chain elements obtained with net14. NSE is achieved around iteration 380. The adiabatic expansion sets in around iteration 660 until the density and temperature are so low that all nuclear reactions freeze out around iteration 1000, with a final composition dominated by ⁵⁶Ni (followed by ⁵²Fe) as it is the dominant nuclei in Type Ia supernova explosions.

The bottom left panel of Fig. 5 shows the nuclear evolution when we use net90-noec. To ease comparison, the colors of all alpha-chain elements are the same in all plots, except that in the larger networks we also show a few (most abundant) isotopes in colored, dashed lines, along with proton and neutron abundances (in black). For most of the pre-NSE regime, nuclear mass fractions follow an evolution similar to that obtained by net14. Only when the temperature is high enough ($T_{9}\geq 7$), many new production channels open, given the high abundance of fresh protons and neutrons. This leads to a longer quasi-NSE regime (compared to net14) and a completely different distribution of abundances in NSE, dominated by moderately neutronized isotopes such as ⁵⁸Ni and ⁵⁴Fe, as well as protons and neutrons. After the adiabatic expansion sets in, most of ⁵⁸Ni captures a proton to become ⁵⁹Cu, which produces ⁵⁶Ni through an inverse $(\alpha,p)$ reaction. Another path for ⁵⁸Ni is two consecutive inverse $(n,\gamma)$ reactions, but this is subdominant due to the rapid cooling during the expansion and evidenced in the very fast decrease in neutron abundance. At freezeout, the most abundant element is again ⁵⁶Ni, followed by some of the heavier isotopes in the network. The final mass fractions of $\alpha-$particles calculated with net14 and net90-noec are quite similar and low, $X_{\alpha}\simeq 10^{-4}$, slightly larger in the former.

The bottom right panel of Fig. 5 corresponds to the results obtained by net90. That is, including electron captures. As expected, there are no large differences until the temperature rises enough to trigger the electron captures on protons. This has a major impact on the NSE, where the proton abundance drops by more than an order of magnitude, producing a large abundance of neutrons. Reactions with neutrons progressively reduce the mass fraction of $\alpha-$elements favoring the synthesis of moderately neutron-rich elements. Now the NSE is largely dominated by ⁵⁴Fe (shown as ⁵⁴Fe+⁵⁷Co+⁵⁸Co in Fig. 5) followed by ⁵⁸Ni. Unlike nuclei, electron captures never achieve equilibrium at the densities considered in this work, and neither does the temperature. The larger binding energy of the dominant neutron-rich isotopes increases the released nuclear energy so that the temperature actually increases during the NSE stage.

This is clearly shown in the top right panel of Fig. 5. There, we show the evolution of temperature ($T_{9}$), density ($\rho_{9}$), and pressure. The dashed lines correspond to net90-noec, where no electron captures were included. The solid lines correspond to net90, including electron captures, with the exception of the orange solid line, which represents the temperature evolution obtained with net14, as a reference. Comparing lines with the same color shows the influence of including electron captures in the nuclear network. This leads to a relevant increase in temperature (at NSE and a good fraction of the adiabatic expansion), whereas pressure gets lower as electrons are captured by protons. In this test, $Y_{e}$ decreased from 0.5 to 0.472. As discussed in Sect. 4.1.2, this effect may have an appreciable influence on the dynamical evolution of hot and dense systems in multidimensional simulations where nuclear reactions are relevant.

Figure 6 shows the mass fractions of IME, IGE, and ⁵⁶Ni-like elements (see its caption for a definition of these groups). The inclusion of $e^{-}$ captures produces substantial changes in the abundances after iteration 400. In particular, the production of IME and ⁵⁶Ni is severely overestimated when $e^{-}$ captures are not taken into account. The opposite is true for the IGE elements, which are under-produced, especially during the adiabatic expansion. The reduced net14 produces even more IME and ⁵⁶Ni than net90-noec. Without $e^{-}$ captures the abundance at freezing is dominated by ⁵⁶Ni whereas the simple inclusion of $e^{-}$ captures on protons shifts the dominance to neutronized IGE elements.

4.1.1 Verifying net90 with WinNet

We present here a comparison between net90 and the reference nuclear network WinNet in light of the magnitudes with the largest impact on the hydrodynamics, such as the rate of released nuclear energy $\dot{\epsilon}_{n}(t)$, pressure $P(t)$ and values of $Y_{e}(t)$ and mean molecular weight of ions $\mu_{i}(t)$. To facilitate comparisons among net14, net90-noec and the larger network but without $e^{-}$ captures, WinNet-noec, we first calculate the evolution with net90-noec and use the resulting trajectory $\{T(t),\rho(t)\}$ as input to the net14 and WinNet-noec calculations. The same procedure is used to compare net90 and WinNet, now with $e^{-}$ captures. We stress that although net90 provides a general picture of the abundances of many isotopes during the run time of multidimensional hydrodynamic calculations, the post-processing of the data is still necessary to obtain detailed nucleosynthetic yields.

Figure 7 shows the value of $\dot{\epsilon}_{n}$ and the mean molecular weight of ions $\mu_{i}$ as a function of temperature for the different nuclear networks at ignition density $\rho=2\times 10^{9}$ g$\cdot$cm^-3 and initial composition $X(^{12}\textrm{C})=X(^{16}\textrm{O})=0.5$. The evolution of $\dot{\epsilon}_{n}$ calculated without and with $e^{-}$ captures during the initial rapid temperature increase (labeled with letter A in the uppermost panels) is fairly similar. In this region and below $T_{9}=4$ the three networks produce identical results, whereas WinNet-noec releases a bit more energy between $4\leq T_{9}\leq 8$. This is followed by a rapid relaxation to the NSE regime (region B), where endothermic photodisintegrations take over (dashed lines in the figures). During nuclear statistical equilibrium (region C), the temperature remains close to $T_{9}=8.6$ and the nuclear energy generation rate drops several orders of magnitude. In this stage, calculations with $e^{-}$ captures release more nuclear energy than those that do not consider the effect of $e^{-}$ captures (see also Fig. 8). The match between net90 and WinNet is quite good in this region, and net90 is able to reproduce the loop around $\dot{\epsilon}_{n}=10^{18}$ ergs$\cdot$s^-1 in the upper right panel of Fig. 7. The agreement also holds for a large part of the adiabatic expansion (zone D), where both networks give comparable results between $4\leq T_{9}\leq 8.5$. Note that below $T_{9}\simeq 6$ the reduced $\alpha$-network net14 differs considerably from the results of net90-noec and WinNet-noec (top left panel in Fig. 7).

The impact of including $e^{-}$ captures in the network is evident in Fig. 8 which shows the time evolution of the nuclear energy generation rate. Before NSE, at $t\leq 0.01$ s, the rate is practically the same in all cases, but the behavior is different thereafter, with calculations including $e^{-}$ captures providing a considerably larger (and positive) generation rate. Furthermore, net90-noec produces practically no energy during the equilibrium phase when $t>0.01$ s, showing very low and alternate $\pm$ values, as expected. On the other hand, WinNet-noec releases more energy with a positive sign (left panel), while during the NSE WinNet initially releases more energy, but the difference with net90 gradually decreases, becoming negligible when $t>0.1$ s (right panel).

The mean molecular weight of ions $\mu_{i}$, shown in Fig. 7, has a maximum relative error $\simeq 20\%$ between net90 and WinNet during the fast relaxation to NSE at $T_{9}\simeq 9.3$. However, such a relative difference is much lower in other combustion stages. During freeze-out, the relative difference decreases to $\simeq 4\%$ and $\simeq 8\%$ in the calculations without and with electron captures, respectively. During the adiabatic expansion, there is almost a perfect match between net90-noec and WinNet-noec while net14 differs significantly (bottom left panel). According to the bottom right panel of Fig. 7 the value $\mu_{i}$ calculated with net90 remains close but below that of WinNet, which is attributable to the different mass fractions of helium obtained with the two networks, as commented below.

Figure 9 depicts the evolution of light elements, $n,p$, and $\alpha$ particles. The latter (leftmost panels) are the most abundant during NSE in all cases, with net90-noec giving values closer to WinNet-noec than net14. Models implementing $e^{-}$ captures produce more helium during the NSE due to the higher temperature and leave more uncombined $\alpha$ particles at freezing (magenta lines in the leftmost panels in Fig. 9) The evolution of neutrons and protons (rightmost panels) is fairly well reproduced by net90 in all cases.

Figure 10 shows the evolution of the molar fraction of electrons, $Y_{e}$, which has a major impact on pressure, which in turn is a crucial magnitude in hydrodynamic simulations. The evolution of pressure is shown on the right-hand Y-axis of the same figure. As we can see, including $e^{-}$ captures, if only on free protons, is worthwhile, as it accounts for $\simeq 2/3$ of the total amount of captures. Figure  10 also shows that the contribution of $e^{-}$ captures on nuclei is actually quite small $\simeq 6\%$ (distance between the thick and thin green solid lines). This suggests that both, the size and the architecture of the reduced networks, are relevant to depict the evolution of $Y_{e}$ and that the thin green line calculated with WinNet but excluding $e^{-}$ captures on nuclei, represents a limiting case. As shown in the same figure, the strong dependence of the electron pressure on $Y_{e}$, $P\propto Y_{e}^{4/3}$, makes $e^{-}$ captures relevant to hydrodynamic simulations of Chandrasekhar-mass models of SNIa.

The impact of changing the ignition density is explored in Fig. 11. At each density and in the range $1\times 10^{8}\leq\rho\leq 5\times 10^{9}$ g$\cdot$cm^-3 the burning temperature is first estimated with net90 assuming isochoric combustion that lasts the same amount of time as in the reference case with $\rho=2\times 10^{9}$ g$\cdot$cm^-3 studied above. The system is then allowed to expand with the characteristic time of the reference case in Table 4. The resulting $\rho(t),T(t)$ profiles are used as input for WinNet to estimate more accurate values of $Y_{e}$ and pressure. The left panel of Fig. 11 shows the dependence of

$$g=\frac{[Y_{e}(0)-Y_{e}^{\emph{net90}{}}(t_{f})]}{[Y_{e}(0)-Y_{e}^{netWN}(t_{f})]}\,,$$ (8)

as a function of the ignition density, where $t_{f}$ is the freeze-out time. The estimator $g$ gives the relative fraction of $e^{-}$ captures estimated with net90 compared to those with WinNet. Although the value of $g$ decreases with $\rho_{9}$, it always remains above $63\%$. The slight reduction of $g$ at high densities is due to the progressive relevance of electron captures on neutronized nuclei in WinNet. The impact on pressure is shown in the right panel of Fig. 11 where the depressurization is evident (i.e. the separation from the blue line), especially at high densities. Despite its great simplicity, net90 does a good job here because it can account for more than $60\%$ of the pressure reduction predicted by WinNet. We note in passing that the final values of $Y_{e}$ always lie above the theoretical limit obtained by Arcones et al. (2010) assuming complete beta equilibrium.

Figure 10: Results of the CO test. Evolution of the electron molar fraction and pressure calculated with net90-noec, net90, WinNet, and WinNet_e+p with electron captures only on nucleons (i.e. only with $e^{-}+p$ and $e^{+}+n$ reactions and no electron captures on heavier nuclei).

Figure 11: Results of the CO test. Left: profile of $g=[\Delta Y_{e}(\emph{net90}{})/\Delta Y_{e}(\emph{WinNet}{})]$ at the freeze-out time as a function of the ignition density. Right: profile of total pressure at the beginning of NSE, $P_{0}$, and at the end point of NSE, $P_{f}$, when the adiabatic expansion begins.

4.1.2 Impact on Chandrasekhar-mass explosion models of SNIa

Although the explosion of a white dwarf with a mass approaching the Chandrasekhar-mass limit was historically the favored progenitor to explain Type Ia supernova explosions, it is today thought not to be the dominant production channel for these events (Maoz et al. 2014). However, such a channel may still account for a few percent of these explosions, which can increase to $20\%$ provided that the ignition density range is extended to $\simeq 6\times 10^{9}$ g$\cdot$cm^-3 (Bravo et al. 2022), while a more standard ignition value is $1-3\times 10^{9}$ g$\cdot$cm^-3 (Nomoto et al. 1984). In any case, this interval of densities is within the working range of net90.

Any reduction in $Y_{e}$ caused by $e^{-}$ captures has, to a greater or lesser extent, an impact on the dynamics of the explosion. However, the current reduced or medium-sized nuclear networks used to perform multidimensional simulations of SNIa explosions do not include deleptonization reactions (e.g., Röpke & Hillebrandt 2005).

In matter dominated by degenerate electrons and in the ultra-relativistic regime, the electron pressure is dominant and approached by $P_{e}=k(\rho Y_{e})^{4/3}$, thus:

$$\left(\frac{\delta P_{e}}{P_{e}}\right)_{\rho}=\frac{4}{3}\left(\frac{\delta Y_{e}}{Y_{e}}\right)\,,$$ (9)

which for $Y_{e}=0.5;\delta Y_{e}=-0.03$ gives $\delta P_{e}/P_{e}=-0.08$ which although not large, contributes to delay the expansion of the fluid element.

Furthermore, deleptonization affects the propagation of thermonuclear combustion. It is commonly accepted that in Chandrasekhar-mass models the thermonuclear combustion is initially propagated by convective motions, where hot bubbles filled with combustion products and formed at low radius float and disseminate the combustion at higher altitudes (Niemeyer & Woosley 1997). The ascending velocity of these bubbles is governed by the expansion factor $\mu=\rho_{b}/\rho_{u}$ between burnt and unburnt matter. The deleptonization caused by electron captures accelerates the pressure balance between the hot bubble and the cool surroundings and increases the value of $\mu$, thus reducing the rising velocity of the hot blobs.

We have checked the efficiency of net90 in handling these Chandrasekhar-mass models by tracking the thermonuclear combustion of the central layer of a white dwarf that expands with a density profile $\rho(t)$ adopted from the W7 model by Nomoto et al. (1984). To do that, we first calculate the self-consistent NSE state at $\rho(t=0)$ with net90 and follow the evolution of the temperature and nuclear species as the density declines, following W7, until the complete freeze-out of the mixture. The $T(t)$ profile obtained with net90 is used to postprocess the output with WinNet.

We provide a summary of the results in Fig. 12 which shows the evolution of density, temperature, and pressure (top-left panel), released nuclear energy rate $\dot{\epsilon}_{n}$ (top right), mean molecular weight of ions (bottom left), and electron molar fraction (bottom right), respectively. These results are qualitatively similar to those of Fig. 7 and the same comments given in Sect. 4.1.1 regarding $\dot{\epsilon}_{n}$ apply here. The evolution of $\mu_{i}$ matches WinNet fairly well, especially at $t\geq 0.8$ s and net90 is capable of accounting for a good amount of $e^{-}$ captures ($\simeq 2/3$ of those with WinNet).

Table 6 provides quantitative information of several magnitudes at $t=2$ s at which any nuclear activity has ceased. The final abundance of protons and neutrons is negligible in both calculations. There is some mismatch in the mass fraction of $\alpha$ particles, but the final value with net90 is low enough to still not affect the molecular weight of ions and hence the ionic pressure. This tiny amount of remaining alphas is due to the limited number of $\alpha$ reactions with neutron-rich nuclei included in net90. The final value of $\mu_{i}$ is well reproduced with the reduced network, and the most abundant nuclei is ⁵⁷Co with an atomic mass and a mass fraction similar to ⁵⁶Fe obtained with WinNet. The total energy released at $t=2$ s is in agreement within $4\%$. Note that according to the last column in Table 6 the energy released by the freely escaping neutrinos is not negligible and with $E_{\nu}(\emph{net90}{})/E_{\nu}(\emph{WinNet}{})\simeq 67\%$. This allows us to monitor the approximate neutrino display during the run-time calculation of the explosion.

Figure 12: Results of the W7 test. Top left: Evolution of temperature and pressure for the density evolution of the central region of the W7 model, Top right: Evolution of the rate of released nuclear energy obtained with net90 and WinNet as a function of $T_{9}$. Letters A-D refer to normal combustion, relaxation to NSE, full NSE, and adiabatic expansion, respectively. Bottom-left: Evolution of the mean molecular weight $\mu_{i}$ of ions obtained with net90 and WinNet . Bottom right: Evolution of the electron molar fraction $Y_{e}$ $\mu_{i}$ obtained with net90 and WinNet

Table 2: Results of the W7 test.

Network	$X_{p}$	$X_{n}$	$X_{\alpha}$	$X_{max}$	$\mu_{i}$	$Y_{e}$	$\epsilon_{n}$	$E_{\nu}$
							$\times 10^{17}$ ergs$\cdot$g-1	$\times 10^{17}$ ergs$\cdot$g-1
net90	$0$	$0$	$7.2\times 10^{-4}$	0.56 (57Co)	$56.8$	$0.472$	$8.57$	$1.2$
WinNet	$0$	$0$	$10^{-10}$	0.53 (56Fe)	$56.2$	$0.459$	$8.95$	$1.8$

4.1.3 The deflagration of a white dwarf with net90 and calculated in three-dimensions

We use net90 to simulate the deflagration of a WD in a true three-dimensional environment. Our goal is to show how the implicit network can handle the explosion of a massive white dwarf made of $50\%$ carbon and oxygen and to evaluate the impact of the $e^{-}$ captures on the outcome of the deflagration. We simulate the explosion of a $1.376$ M_☉ WD ($\rho_{c}=2.9\times 10^{9}$ g$\cdot$cm^-3) following the artificial incineration of six small bubbles in the central volume of the WD, at random distances $r\leq 80$ km from the center and aligned with the three coordinate axis. We perform three-dimensional low-resolution calculations with the smoothed particle hydrodynamics (SPH) code SPHYNX (Cabezón et al. 2017) adapted to handle nuclear flames and deflagrations, and whose features will be reported elsewhere (see also García-Senz et al. 2016, for a description of the nuclear flame propagation algorithm). We conducted two simulations with 1 million of SPH particles, including and neglecting electron captures in net90.

Figure 13: XY slice showing the distribution of four groups of isotopes during the deflagration of a WD at two elapsed times: $X_{\alpha}$ (upper row), intermediate-mass elements, $X_{IME}$ (second row, see Fig. 6 for a definition of these groups), iron-group elements $X_{IGE}$ (third row), and ⁵⁶Ni (last row). The first two columns were obtained neglecting $e^{-}$ captures whereas the other two include them. Each box has a side of $1300\times 1300$ km.

Table 3: Final yields during the deflagration of a white dwarf calculated in three-dimensions

Network	p	n	$\alpha$	IME	IGE	56Ni
net90	$2.4\times 10^{-6}$	$2.6\times 10^{-13}$	$1.2\times 10^{-4}$	$5.9\times 10^{-5}$	$4.8\times 10^{-2}$	$6.4\times 10^{-2}$
net90-noec	$2.1\times 10^{-4}$	$5.5\times 10^{-14}$	$1.1\times 10^{-4}$	$4.1\times 10^{-4}$	$1.4\times 10^{-5}$	$1.1\times 10^{-1}$

Figure 14: Distribution of $Y_{e}$ in a XY slice (left) at $t=0.73$ s and evolution of neutrino luminosity in the deflagration test.

Figure 15: Evolution of different energies during the deflagration of a WD. Solid lines correspond to net90 and dashed lines to net90-noec.

Figure 13 shows the distribution of several nuclear species in an XY slice of the WD at two elapsed times. Four of the six high-temperature regions or hot bubbles are clearly visible in the slice. These hot bubbles with $T\simeq 8\times 10^{9}$ K rapidly ascend by Archimedian flotation and deform with the overall effect of increasing the contact surface between the hot ashes and the CO fuel, which accelerates the effective nuclear burning rate. After $t\simeq 1.2$ s the expansion of the WD quenches the nuclear combustion. The amount and distribution of the different families of isotopes in the burned regions are dependent on the $e^{-}$ capture rate. Regions where $e^{-}$ captures are relevant, last two columns in Fig. 13, produce more IGE elements and less ⁵⁶Ni and vice-versa which agree with the results found in Sect. 4.1.2. More detailed values of the final yields of the different families of elements are shown in Table 7. Figure 14 depicts the distribution of the electron molar fraction $Y_{e}$ at $t=0.73$ s (left panel) along the rising plumes of burnt material and with minimum values $Y_{e}\simeq 0.47$ at the tip of the bubbles and in the very central region of the WD. The released neutrinos escaping from the WD produce a display, shown in the same figure (right panel) which could be detected if the event occurs at a sufficiently close distance.

The evolution of different energies during the simulations with and without $e^{-}$ captures is shown in Fig. 15. First note that the total amount of released nuclear energy, $E_{nuc}\simeq 1.8\times 10^{50}$ ergs, is not enough to unbind the star. The total amount of burnt material, $M_{b}\simeq 0.12$ M_☉, is well below the limit of $\simeq 0.30$ M_☉ needed to get the WD unbound (Bravo & García-Senz 2009) and within the range of ”failed deflagrations” described in Röpke et al. (2007). In contrast, the WD will undergo a pulsation having a second chance to explode at longer times, during the re-contraction period (Bravo et al. 2009). Interestingly and according to Figs. 15 and 13, the inclusion of $e^{-}$ captures in the hydrodynamic calculation weakens the kinetic energy during the stages of fast combustion at $t\leq 1.2$ s and reduces the vertical development of the rising plumes of burnt material . Hence, we conclude that $e^{-}$ captures is a necessary ingredient to perform high-precision multi-D calculations of the explosion of a WD approaching the Chandrasekhar-mass limit.

The network net90 is therefore capable to handle all nuclear combustion stages: initial explosive induction of the CO fuel, rapid combustion, QNSE and NSE and breaking of the equilibrium during the fast expansion in a multi-D simulation of the explosion of a massive WD. It also provides a reasonable first picture of the produced nucleosynthetic yields. The time-step always remained $\Delta t\simeq 2\times 10^{-4}$ s in the computed models above, an affordable time-step which is not very different from the current hydrodynamic Courant time.

4.2 He burning

Helium combustion might take place in the outer layers of an accreting white dwarf. Material from a companion star might be ripped away and fall onto the surface of the WD, where it piles up instead of burning. When density is high enough, an ignition will take place at one point (or several) of the He shell, triggering a detonation wave that propagates in all directions, eventually igniting the C+O core of the WD. This is the so-called double detonation mechanism (DDT) (Woosley & Weaver 1994).

In this scenario, the densities in the He shell are significantly lower than those in the core of the WD, with $\rho_{6}\sim 0.5-2$, and maximum temperatures of $T_{9}\sim 2.5-5$. With such conditions, electron captures are expected to play no significant role.

4.2.1 Calculations with net90 and net14

Figure 16 shows the result of the He burning test using the initial conditions described in Table 4. The top-left panel shows the outcome using net14, while the lower panels present the results using net90-noec (left) and net90 (right). When comparing all three results, it is evident that net14 provides very different yields, dominated by lighter metals ⁴⁴Ti, ⁴⁸Cr, and ³⁶Ar. The larger networks products are dominated by elements of the iron family, such as ⁵⁶Ni and ⁵²Fe. The reason for this is the architecture of the networks. The only way net14 can reach ⁵⁶Ni is through alpha captures, which are very suppressed at the upper end of the network due to the low density and temperature. On the other hand, net90-noec (left) and net90 (right) have many other paths through proton captures. These results imply that a small network such as net14 is not suitable for explosive scenarios with low densities in terms of the yields produced.

Furthermore, the thermodynamic evolution of net14 and net90 also shows significant differences, as shown in the top-right panel of Fig. 16. Temperature and pressure values during the high-temperature episode are lower than those obtained with net90, especially pressure. However, having a correct representation of pressure is crucial in the development of self-sustained helium detonations that trigger the explosion. Additionally, as expected at this low density, the inclusion of electron captures has no effect on either the temperature or the pressure. Only elements with very low abundance might be affected, such as ⁵⁸Cu. Unlike the precedent tests, the electron mass fraction slightly increases $Y_{e}=0.5+1\times 10^{-8}$ at the freeze-out in net90, due to the positron captures on free neutrons, which are not blocked by degeneracy effects at the densities and temperatures involved in this scenario.

Figure 16: Results of the He test. Evolution of mass fractions obtained with net14 (top left), net90-noec (bottom left), and net90 (bottom right). Vertical dashed lines show the approximate limits of the high-temperature combustion region, when the artificial adiabatic expansion starts, and the region where nuclear reactions freezeout. The panel on the top right shows the evolution of temperature and pressure for net90 and net14. The calculation with net90-noec follows exactly the same trajectory as with net90 and is not shown. The density evolution is the same for all three nuclear networks.

4.2.2 Verifying net90 with WinNet

We postprocessed the calculations of net90 with both net14 and WinNet and analyzed and compared the ensuing results. Figure 17 shows the evolution of the nuclear energy generation rate calculated with the three networks. The agreement between net90 and WinNet is fairly good, both lines stay close but the larger network leads to a more pronounced, albeit narrow, maximum at $t\simeq 0.012$ s. The match of the reduced net14 is clearly worse, which reinforces the conclusion in Sect. 4.2 that reduced $\alpha-$like networks are not adequate to address Sub-Chandrasekhar-mass explosion models of Type Ia supernova.

The mean molecular weight of ions calculated with net90 perfectly fits that of WinNet while net14 leads to somewhat lower values (see Fig. 18). The corresponding value of $\alpha$ particles is shown in the same figure. Free protons evolve similarly until $t\simeq 1.5$ s but the net90 estimate remains below that of WinNet at $t\geq 1.5$, due to the increasing number of p-production channels in the larger network at later times.

Figure 17: Results of the He test. Evolution of the nuclear energy generation rate calculated with net14, net90, and WinNet.

Figure 18: Results of the He test. Evolution of the mean molecular weight of ions, $\mu_{i}$ (leftmost Y-axis), and protons and helium mass fractions (rightmost y-axis) calculated with net14, net90, and WinNet.

4.3 Si burning

The photodesintegration of Si is relevant in scenarios such as a core-collapse supernova, where the core of a massive star implodes due to its own mass and subsequent deleptonization (e.g., Woosley et al. 2002; Janka et al. 2007). The temperature rises and heavy elements are photodissociated. This process happens nearly isothermally because of the large thermal conductivity of the electrons. To mimic such conditions, we artificially increased the specific heat, which is usually provided by the EOS, so that $\Delta T\sim 0$. Table 4 provides the initial conditions and Fig. 19 shows the results of the test. The setting of this test is similar to that described in Timmes et al. (2000) who compared the performance between two reduced nuclear networks in light of the results obtained with a larger network with N=489 isotopes.

4.3.1 Calculations with net90 and net14

The top-left panel of Fig. 19 clearly shows that, with this high temperature, Si is rapidly destroyed, producing a high amount of alpha particles. These are used to produce all the elements of the alpha-chain, reaching NSE very quickly. In particular, ⁵⁶Ni and alpha particles are the most abundant elements. During adiabatic expansion, a large fraction of this reservoir of alpha particles is used to climb the alpha chain using the seeds present at NSE, producing even more ⁵⁶Ni and an overall decrease in all the other elements. Once the temperature has decreased sufficiently ($T_{9}\lesssim 5$), no more ⁵⁶ Ni is produced, leading to a slight increase in all other elements through alpha captures until they slowly quench as reaction channels are closed from heavier to lighter elements. The final result after freezeout is dominated by ⁵⁶Ni and ⁴He, followed by heavy elements and finally lighter species.

The bottom panels of Fig. 19 show the same test performed with net90-noec (left) and net90 (right). First, we notice that there are very small differences among them. This is due to the low initial density of this test ($\rho_{9}=0.01$). At this density electron captures are subdominant, and only slight differences can be noticed in elements with very low mass fraction. The inclusion of proton and neutron captures is what has the greatest impact compared to the results obtained with net14. In net90-noec and net90 there are many more nuclear paths open for the burning of ²⁸Si. This can be seen in a much lower mass fraction of ²⁸Si in the NSE, the dominance of neutronized isotopes, such as ⁵⁸Ni and ⁵⁴Fe, and the high abundance of protons in that regime.

Once the adiabatic expansion sets in the dominant isotopes during the NSE, ⁵⁸Ni and ⁵⁴Fe have their mass fractions severely reduced. The mass fraction of alpha particles also decreases at the beginning of the expansion, facilitating the build up of ⁵⁶Ni. The freezeout is dominated by ⁵⁶Ni and ⁴He, followed with a stratification similar to that of net14 but with more abundance of heavier elements. From the top-right panel we can see that the temperature and pressure evolution is very similar in all cases, and that the inclusion of electrons has no noticeable effect. Figure 20 shows the mass fractions of IME, IGE, and ⁵⁶Ni-like elements. In contrast to the CO test mentioned earlier, the inclusion of $e^{-}$ captures on protons does not significantly affect the results. Thus, net90-noec and net90 virtually produce the same results. However, the differences between these and net14 are significant during the isothermal combustion phase, with the reduced network generating a substantially higher amount of IME and ⁵⁶Ni. Nevertheless, all networks synthesize similar final amounts of nickel and IME.

Figure 19: Results of the Si test. Evolution of nuclear mass fractions obtained with net14 (top-left), net90-noec (bottom-left), and net90 (bottom-right). Vertical dashed lines show the approximate limits of the NSE, when the artificial adiabatic expansion starts, and the region where nuclear reactions freezeout. The panel in the top right shows the evolution of temperature and pressure for net14 and net90. The evolution of density is the same for all three nuclear networks.

Figure 20: Results of the Si test. The mass fractions in Fig. 19 are grouped in three nuclear families (see Fig. 6 for a definition of these groups) for net90-noec (dashed lines), net90 (continuum lines) and the net14 ($\alpha$-network).

4.3.2 Verifying net90 with WinNet

As in the precedent tests, we take advantage of the temperature and density evolution obtained with net90 to ’post-process’ the calculation with net14 and WinNet and obtain the nucleosynthetic yields and release of nuclear energy, which are in turn used to verify the reduced networks. The main results regarding these magnitudes are shown in Figure 21. The upper-left panel depicts the nuclear energy generation rate that proceeds in two separate regimes. For $t\lesssim 0.01$ s nuclear reactions are endoenergetic due to the prevalence of photodisintegration reactions. In this region net14 and net90 match better with each other than with WinNet, although overall net90 evolves closer to WinNet than net14. The evolution in this region agrees with that obtained by Timmes et al. (2000), although their large-network calculation stays closer to their reduced network estimations than in our work. We observe slight differences between the networks, as (n,p) reactions are dominating for the prevailing conditions. Although all networks strive to achieve an NSE composition, WinNet can achieve it fastest as, in contrast to the other networks, it includes (n,p) reactions. At $t\geq 0.01$ s the energy release is positive and the match between net90 and WinNet improves in that region while net14 fails to reproduce $\dot{\epsilon}_{n}$ between $0.01\leq t\leq 0.3$ s. The upper-right panel shows the evolution of the nuclear energy released until the elapsed time $t$, $\epsilon_{n}(t)$ on a linear scale. Although the pure $\alpha-$ chain gives in general poorer results than net90, the error in the released energy at freeze-out of both networks is similar (see also the second to last column in Table 8)

The lower panels of Fig. 21 show the evolution of the mass fraction of light elements and the mean molecular weight of nuclei calculated with the three networks. The calculation with WinNet leads to a larger population of $\alpha$ particles for $t\lesssim 10^{-3}$ s which lowers the value of $\mu_{i}$. However, net90 and WinNet give almost identical results when $t\geq 10^{-3}$ s. Overall, the values of $\alpha$ and $\mu_{i}$ obtained with net14 do not compare so well with WinNet as those computed with net90. The evolution of free neutrons and protons obtained with net90 and WinNet is fairly similar.

Additional quantitative information on the Si test is included in Table 8. The most abundant element when nuclear reactions have stopped is ⁵⁶Ni, which is abundantly produced and with similar values in all three cases (third column). The value of $Y_{e}$ remains very close to $Y_{e}=0.5$ in all calculations (fifth column) and consequently, the energy drained by the freely escaping neutrinos $E_{\nu}$ is only a small fraction of the total released energy at the freeze-out time (last column). However, $E_{\nu}$ is more than ten times larger in the WinNet calculation because $e^{-}$ captures on nuclei proceed from the very beginning while net90 has to wait until there are enough protons in the mixture.

Table 4: Results of the Si test.

Network	$X_{\alpha}$	$X_{max}$	$\mu_{i}$	$Y_{e}$	$\epsilon_{n}$	$E_{\nu}$
					$\times 10^{16}$ ergs$\cdot$g-1	$\times 10^{15}$ ergs$\cdot$g-1
net14	$0.085$	$0.900$ (56Ni)	$26.51$	$0.50000$	$5.80$	$-$
net90	$0.096$	$0.843$ (56Ni)	$24.98$	$0.49994$	$3.98$	$0.13$
WinNet	$0.096$	$0.853$ (56Ni)	$24.88$	$0.49979$	$4.43$	$1.9$

Figure 21: Results of the Si test. The upper panels show the evolution of the nuclear energy generation rate (left) and total released nuclear energy (right), while the lower panels depict the evolution of the mass fraction of the lighter elements $n$, $p$, and $\alpha$ particles (left), and mean molecular weight of ions, $\mu_{i}$ (right)

5 Computational cost

In previous sections we showed that using a medium-sized nuclear network coupled with hydrodynamical simulations has clear advantages in terms of accuracy and reliability of the results compared to small, $\alpha$-like, networks. However, this comes with a cost. In this section, we present the measurements of the computational cost of each nuclear network.

Table  5 shows the average computational performance of net14, net90-noec, and net90 for 100 full runs of the CO test executed with a single thread on an Intel i7-4790 CPU (3.6 GHz), compiled with gfortran and -O3 flag. All runs execute $\sim 1500$ global iterations and $\sim 6000$ NR iterations, which correspond to $\sim 600,000$ calls to each nuclear subroutine. The columns of Table 5 show, from left to right, the name of the network, the number of global iterations, the ratio of NR iterations per global iteration, the average time per global iteration (that is, for the entire cycle of NR iterations), the average time per network call (namely, per individual NR iteration), the average percentage of time spent on the sparse matrix solver, and the scaling factor for each network call with respect to net14.

From these results, the computational cost of the network scales approximately as $N\log N$ with the number of reaction rates (being $N=17$ for net14 and $N=161$ for net90). The overload of including $e^{-}$ and $e^{+}$ captures is very small ($\sim 0.4\%$). Approximately a fourth of the time per network call is spent in the sparse matrix solver. Hence, any improvement in this part of the code would have a noticeable effect overall.

Coupling chemical equations with temperature increases the algebraic complexity of the integration scheme, and one may wonder if the increment in the timestep during QNSE-NSE conditions compensates for the extra algebraic effort. At $\rho_{9}=1,T_{9}=8$ the timestep during the NSE plateau calculated with net90 is $\Delta t\simeq 10^{-4}$ s, not much different than the Courant time at the center of the WD in current multi-D calculations $\tau_{C}=0.4\leavevmode\nobreak\ (\Delta x/c_{s})\simeq 4\times 10^{-4}$ s for a cell size $\Delta x=5$ km and sound speed $c_{s}=5000$ km$\cdot$s^-1. Removing the thermal implicit coupling in net90 drops the time-step to $\Delta t\simeq 10^{-8}-10^{-10}$ s to control the numerical fluctuations in temperature and abundances. Therefore, incorporating temperature into the Jacobian solver of the nuclear network increases the timestep by a factor $10^{4-5}$ when $T_{9}\geq 8$.

A different stabilizing procedure, described in Paxton et al. (2015), consists on raising the floating point limit in arithmetic operations, for example from the commonly used double precision to quadruple precision, which increases the stability of the integration. Both approaches are not exclusive and could be used jointly.

Table 5: Computational cost of each nuclear network.

Network	It.	NR/It.	$\overline{t}$/It.	$\overline{t}$/NR	$\overline{t}_{\mathrm{matrix}}$	scaling
			($\times 10^{-5}$ s)	($\times 10^{-6}$ s)	(%)
net14	1501	3.47	$1.11\pm 0.23$	$3.21\pm 0.67$	$12.35\pm 8.18$	-
net90-noec	1524	3.72	$17.85\pm 0.62$	$47.93\pm 1.67$	$25.02\pm 2.45$	$\times 14.93$
net90	1508	3.85	$18.53\pm 0.66$	$48.10\pm 1.72$	$25.07\pm 2.57$	$\times 14.98$

6 Conclusions

In this work, we propose and verify a medium-sized nuclear network especially addressed to multidimensional simulations of Type Ia supernova explosions. The network net90, with 90 species from ⁴He to ⁶⁰Zn has two distinctive features. First, the integration of nuclear species is carried out using an implicit scheme that also includes temperature in the Jacobian matrix (Mueller 1986; Cabezón et al. 2004). Such thermal coupling stabilizes the integration scheme during the QNSE and NSE stages, allowing to take considerably larger timesteps than schemes where temperature changes are calculated in an explicit form. Furthermore, both the implicit integrator and the matrix solver have been improved with respect to those described in paper I.

The second feature, and a novelty of this work, is the inclusion of weak interaction processes in the reduced network. Electron captures on protons and nuclei have been traditionally incorporated into large nuclear networks and NSE tables supporting spherically symmetric explosion models of SNIa (e.g., Nomoto et al. 1984) and CCSN and in a few multi-D calculations of Type Ia supernova explosions that interpolate the value of $Y_{e}$ from large NSE tables (e.g., Moll & Woosley 2013). Weak processes are also currently incorporated to post-process the output of multi-D calculations of both kinds of supernova explosions. Nevertheless, and as far as we know, these processes have been ignored by all reduced networks devoted to multi-D hydrodynamic simulations so far. In particular, electron and positron captures on protons and neutrons, respectively. Although it may seem that including weak processes only on free particles is too simplistic, we have shown that net90 manages to account for $\simeq 2/3$ of electron captures under the typical conditions prevailing in Type Ia supernova explosions. Therefore, the pressure of electrons is substantially better estimated than if such reactions are ignored, allowing for a more accurate depiction of the dynamics of the explosion.

We verified net90 by direct comparison with the results obtained with a much larger network, WinNet (Reichert et al. 2023), under typical supernova conditions. These include explosive carbon and oxygen mixtures with maximum densities, prior expansion, in the range $10^{8}$ g$\cdot$cm^-3$\leq\rho_{0}\leq 5\times 10^{9}$ g$\cdot$cm^-3 and edge-lit helium ignition at $\rho\simeq 10^{6}$ g$\cdot$cm^-3. These density ranges are representative of many SNIa explosion models currently considered in the literature, either of Chandrasekhar or Sub-Chandrasekhar mass type.

We found that in the high density regime net90 gives a rather good approach to the released nuclear energy, total pressure and mean molecular weight of ions, $\mu_{i}$, and provides a much better depiction of produced nuclear yields than those obtained with small networks such as net14. According to Table 6, the medium-sized network calculation gives a good estimate of $\mu_{i}$ and total released nuclear energy of a toy model that mimics the central shell of the W7 supernova model by Nomoto et al. (1984) at freeze-out time. The relative deviations of both magnitudes with respect to those computed with WinNet are smaller than $4\%$. However, the deviations of the nuclear energy generation rate could occasionally be larger, as shown in Figs. 7 and 8. Table 6 also shows that net90 leaves more uncombined $\alpha$-particles, which is attributable to the progressive loss of $\alpha$-reaction channels as neutronization increases.

Explosive helium ignition at $\rho\simeq 10^{6}$ g$\cdot$cm^-3 is also supported by net90. This corresponds to point-like edge-lit ignitions near the surface of WDs with masses $0.6\leq M_{wd}\leq 1.0$ M_☉, propagating as detonation waves (Gronow et al. 2021). Too small networks, such as net14, are not very adequate here because temperature and pressure jumps are not accurately captured. This is shown in the upper-right panel of Fig. 16. In contrast, the results with net90 are fairly in agreement with those with WinNet in terms of nuclear energy generation rate $\dot{\epsilon}_{n}$, $\mu_{i}$, dominant nuclei, and mass fractions of light particles (Figs. 17 and 18) which ensures pressure compatibility. Electron and positron captures do not play a significant role in this density regime.

Finally, we have applied net90 to the combustion of silicon at densities $\rho_{0}\simeq 10^{7}$ g$\cdot$cm^-3, this being a standard nuclear test (e.g. Timmes et al. 2000) that is more akin to Core Collapse Supernova (CCSN) than to SNIa. In this test, the silicon is rapidly disintegrated by photons, leading to a negative generation of nuclear energy during some time. According to Fig. 21, the late stages of this phase are neither well described by net14 nor by net90, although the latter is slightly closer to the WinNet solution. After a while, the nuclear generation rate turns to positive values and the match between net90 and WinNet improves, while the reduced net14 is still unable to give realistic results. The electron abundance does not change too much during the interval of time considered in the test. However, both the change of $Y_{e}$ and the energy released by neutrinos are much larger when calculated with WinNet. All of this suggests that, in agreement with Navó et al. (2023), caution should be taken when using medium-sized networks for CCSN related scenarios.

In summary, we proposed and verified a nuclear network implicitly coupled with temperature and built around a standard $\alpha$-chain but expanded to host 90 nuclei including free protons, neutrons, and electrons (Fig. 1). This network includes $e^{-}$+p and $e^{+}$+n capture processes, which make it ideally suited for being implemented in multidimensional simulations of Type Ia Supernova. The feasibility of net90 to handle the deflagration of a massive WD in a three-dimensional calculation was confirmed in Sect. 4.1.3.

The price to pay for such a better depiction of the nuclear energy generation rate and the effect on the overall physical state of matter comes in the form of a larger computational cost. This cost scales roughly as $N\log N$, with the number of nuclear reactions implemented. According to our measurements, net90 is a factor $\times 15$ slower than net14. Even if such a factor cannot change, there is still room for future improvements that could be made to reduce the overall wall-clock time per network call, for example in the architecture of the matrix solver which could be of sparse type, such as the widely used PARDISO, or dense matrix solvers such as LAPACK. Admittedly, an increased computational cost is unavoidable when using a larger nuclear network. Yet, the benefits of considerably more accurate simulations outweigh this cost.