pynucastro 2.1: an update on the development of a python library for nuclear astrophysics

pynucastro1 is an open-source python library that provides visualization and analyze techniques to classify, construct, and evaluate nuclear reaction rates and networks. It provides tools that allow users to determine the importance of each rate in the network, based on a specified list of thermodynamic properties. Additionally, pynucastro can output a network in C++ or python for use in simulation codes, include the AMReX-Astrophysics simulation suite. We describe the changes in pynucastro since the last major release, including new capabilities that allow users to generate reduced networks and thermodynamic tables for conditions in nuclear statistical equilibrium.


Introduction
In stars, nuclear reactions are the primary source of energy, and their composition drives their evolution [1,2].Understanding stars means understanding how nuclear reactions work, and their impact on the whole network under different environments and conditions.Simulations of stars model their steady evolution and explosive deaths, requiring accurate and efficient nuclear reaction networks.In the current exascale era, GPU-enabled reaction networks are also required.pynucastro was developed to meet this need-it provides a framework for exploring reaction rates and networks interactively in python and for exporting a reaction network to python or (GPUenabled) C++ code for use in a simulation.pynucastro complements other reaction libraries in the field, for example, SkyNet [3] and WinNet [4], but its goals differ slightly: the emphasis in pynucastro is on interactive exploration of a network, to allow the user to create a network tailored to a particular science problem, and on the ability to export the network in a format that can be used in a multidimensional simulation code.This also allows us to easily explore the effect of the network size on a simulation, as done, for example, in a recent study of flames in X-ray bursts [5].
Recently pynucastro hit the 2.0 milestone [6], with the addition of weak rates, a nuclear statistical equilibrium solver, the ability to derive reverse rates through detailed balance, support for approximate rates, C++ output for the AMReX-Astrophysics suite of simulation codes, and some basic methods to help understand the importance of different rates in a network.The development has since continued, with a lot of new features leading to the 2.1 release, which we summarize here.

Expanded weak rates
The sensitivity of weak rates in core-collapse supernovae, novae, X-ray bursts, and massive stars astrophysical objects and environments plays a significant role in their nucleosynthesis and energy generation [7,8].In pynucastro 2.0, we included the weak-rate tabular data [9,10] for nuclei A = 17 to 28.This was driven by our desire to simulate convective Urca [11].From these tables, we extracted the β − -decay rate (β − ), electron-capture rate (e − ), gamma-energy rate (Γ e ), and the neutrino-energy rate loss (Γ ν ), reformatted in cgs units.We know that both weak reactions, β − -decay and electron-capture, produce the following changes in a nucleus (Z, A): However, the positron-capture and the β + -decay reactions produce the same effects on the nuclei (Z, A): respectively.These additional channels become important at hotter and denser conditions than seen in the Urca simulations, and some rate sources tabulate them in addition to their electron counterparts.In pynucastro 2.1, we combine the bare β − -decay and the positron-capture rates in a single rate, and similarly, the bare electron-capture and the β + -decay rates.This means that each sequence (Z, A) → (Z + 1, A) and (Z, A) → (Z − 1, A) are represented by a single effective link in a network.In pynucastro 2.1, we have improved our weak-rate tabular data availability by including the A = 45 to 65 rates in [12][13][14] to support simulations of massive stars.We also reformatted our tables to consider the original logged tabular entries instead, in cgs units, and in our effective rate convention: log(β − ), log(e − ), log(Γ e ), and log(Γ ν ) to each tabulated pair (log ρY e , log T ).A comparison between pynucastro 2.0 and pynucastro 2.1 weak rate nuclei availability is depicted in Fig 1 .Our addition of rates is driven by our science goals, and more rates will be added as needed.

A new bilinear interpolation scheme for the tabular rates implementation
In pynucastro 2.0, the tabular rate evaluation in a PythonNetwork simply returned the closest tabular entry for the (ρY e , T ) input pair, while the AmrexAstroCxxNetwork C++ network class performed a bi-linear interpolation over the tabular input pairs.Thus, the PythonNetwork functionality was made consistent with the more accurate approach in AmrexAstroCxxNetwork.
In our new release, we reformatted our tabular entries in terms of log 10 (ρY e ) and log 10 (T ) and implemented a bi-linear interpolation scheme for both network classes.For the PythonNetwork class, we implemented a bi-linear interpolation scheme stored in the TableInterpolator class that can be used both for interactive evaluation of the rates and when the network is written out to a python module for use elsewhere.An example of the new interpolation interface is shown in the following script: Note, that we construct first the TableInterpolator object for the na23 ne23 rate, and then we evaluate the rate.

A Simpler C++ Network
The original C++ network support in pynucastro is through AmrexAstroCxxNetwork, and utilizes the data structures in the AMReX library [15] to provide the ability to offload the reaction network onto GPUs.This is compatible with any simulation code that uses AMReX, like Castro [16] or MAESTROeX [17].pynucastro now supports a pure C++ network, through the SimpleCxxNetwork class, that can output the RHS and Jacobian of the network.This was developed to allow other simulation codes to use pynucastro networks.At the moment, screening is not implemented (but it is easy to add a hook for a user-supplied screening function) and tabulated rates are not supported.For codes written in C, the C++ interfaces can be accessed via an extern "C" {} interface.This network output will be expanded as community demand dictates and we hope this new interface will help spur the adoption of pynucastro-generated networks by other simulation codes.

Support for custom rates
The base Rate class in pynucastro was improved to make it easier to add a custom rate.This can allow a user to explore variations or approximations to the standard rates in the pynucastro library.For example, suppose we want to add a simple powerlaw rate: This can be accomplished by creating a new class that inherits from Rate, stores the parameters r 0 and T 0 , and calls the super-class initialization, e.g.This sets up all of the logic needed to understand how to use our new rate, MyRate, in a network.Then we only need to add a few additional methods.If we only want to use the rate in python, then we just need to write an eval() method that computes the rate given temperature, such as: If we want to export the network to python or C++ code, then we need to add one additional method that writes the function string.This custom rate can then be used in any function in pynucastro that takes a Rate object.

The DRGEP modified reduction method implementation
In pynucastro 2.0, we provided a simple find unimportant rates method to help reduce the size of a network.In the latest release, a more extensive reduction method framework is included.The idea of selecting the most important nodes of a graph is also a problem with important applications in combustion and chemical kinetics [18][19][20].The general idea behind these methods starts with a large-sized directed graph G representing the network and a user-defined selection of fixed target nodes T .For each pair of nodes (A, B) of the graph G, one can define the following matrix coefficients r A,B : Here, P A are all the rates on which A is a product, similarly C A are all the rates on which A is a reactant, ν A,i is the stoichiometric coefficient of A in the reaction i, ω i is the i-th rate evaluation, and δ i,B is the Kronecker delta object.From each target nucleus T , a weight R T B is defined for each connected nuclei B to the target T as follows: where the index 1 ≤ i ≤ n represents all the nuclei between T and B for a given path p, including both.The weights R T B can be calculated using a modified Dijkstra's algorithm [20][21][22].From ( 5), all of the paths from a nucleus B to a given target T will have weights less than R T B .For each target T , we can flag whether a nucleus B is considered unimportant by comparing to a user-specified tolerance choice ϵ T : R T B < ϵ T .
We can then remove nuclei with R T B less than ϵ T for all possible targets.The reduced network will be a sub-graph G ′ ⊂ G that will contain at least T nodes: Following the particular case of the DRGEP modified method described in [20][21][22], we have implemented the DRGEP modified reduction method in pynucastro.A direct application of this method is to identify all of the connections that are important for a particular subset of nuclei from a known reaction network set of nuclei.By setting each nucleus from the known network as a target, we can explore which rates are significant from all the currently supported rate libraries, up to an endpoint nucleus given by the user.For example, let us explore how several nuclei may be connected from the known CNO cycle to all the species up to 56 Ni.The importance of 56 Ni endpoint species is explained, for example in [23].We start by creating a pool of rates with all the currently supported REACLIB rates up to 56 Ni with a fixed solar composition.This two-step process is achieved by the load network() method with a fixed network endpoint in 56 Ni.This will create a PythonNetwork containing all the nuclei with A and Z less than the endpoint and all the rates connecting them.Then, by selecting 25 representative nuclei of the original hydrogen-CNO burning reaction network, we construct the targets list.Finally, each link that is ultimately connected to each target, will be evaluated by a weight-constraint function.The tols list contains a tolerance that will determine whether each link should be severed or not based on the previous evaluation.From the original pool network composed of a total of 6758 rates with 584 nuclei, we can drastically reduce it to 607 rates with 63 nuclei.Thus, from all the rates REACLIB rates with nuclei up to 56 Ni, we were able to extract the most important rates connected to the original network nuclei targets.targets = [ " p " , " d " , " he3 " , " he4 " , " li7 " , " be7 " , " be8 " , " b8 " , " c12 " , " n13 " , " n14 " , " n15 " , " o14 " , " o15 " , " o16 " , " o17 " , " o18 " , " f17 " , " f18 " , " f19 " , " f20 " , " ne18 " , " ne19 " , " ne20 " , " ne21 " ] An exploration of alternative reduction method implementations, energy-based weightconstraints functions, and the analysis of their implications in large nuclear reaction networks should be studied further.These algorithms may prove to be useful in highlighting the most important species and rates from an extensive nuclear reaction network, improving our understanding of the group of reactions that dominate our problem and the efficiency of our calculations.

Interfacing with Julia
The SciPy library [24] provides an ODE integration function solve ivp() that provides many different ODE solvers, including BDF solvers, tailored for stiff systems.This is the main way we suggest how a PythonNetwork network should be integrated.However, there are many ODE methods that are not implemented in SciPy, some of which can be more efficient for our problems.
The Julia library DifferentialEquations.jl[25] expands the possibilities to consider.In pynucastro 2.1, we demonstrated how to integrate a python network using the Julia solvers.Performance relies heavily on Numba acceleration [26].This opens the door to further exploration of different ODE solvers.
This interface makes it straightforward to test the variety of solvers on the network of interest.Below is an example using a simple CNO nuclear reaction network with an initial composition of a standard solar mix.This code integrates the system over the specified time using the FBDF solver in the Julia package.Note that the Julia interface requires that the righthand side and Jacobian functions should be combined into a single interface via ODEFunction.After integration, sol.stats can provide a variety of different statistics including the number of function evaluations iterations, and linear algebra operations it performs.These can be useful in determining the performance of the solver.

NSE table construction
pynucastro 2.0 saw the addition of a robust NSE solver.In version 2.1, the NSE solver has been split off into its own network type, NSENetwork, to allow for the addition of new features that support creating NSE tables for use in simulation codes.To support this capability, a new method Composition.binas() has been added to bin a large Composition, for example, representing the NSE state from 100 nuclei, down to a smaller Composition that will be what is carried on the grid in a simulation code.We've also added the ability to supply an optional rate filter function when evaluating the RHS of a network.For the NSE evolution, this can be used to only compute the evolution due to weak interactions, which cause an evolution in Y e further evolving the NSE state.Likewise, the neutrino loss from weak rates can be returned separately when evaluating the energy production in a network.Together these new methods allow us to easily create an NSE table like that in [27], and an interface for this table creation is under development.

Summary and Future Work
The new release of pynucastro added many new features: an expanded weak rate library with a new interpolation scheme over the tabular rates, a new template to cover bare C++ support, the ability to construct a custom rate, network reduction capabilities, a new Julia integration interface, and the capability of creating an NSE table from the network.All of these new features are driven by the science applications that the developers are working on, including studies of convective burning in massive stars leading up to core-collapse.
The next major developments will likely be the addition of the remaining pieces needed to support a self-heating burn completely in python.This includes the addition of an equation of state (to get the temperature evolution) and neutrino cooling sources.Pull requests with initial implementations of each are under review, with the equation of state based on the Helmholtz equation of state [28,29].Once the previous step is completed, we may explore the addition of other equations of state, for example [30].
Finally, all of the major machine learning libraries have python interfaces and much AI research relies on these libraries.By providing a clean, extensible python interface to nuclear reaction rates and networks, we enable the use of pynucastro for machine-learning applications.In particular, since reaction network integration can be one of the most time-consuming parts of a simulation, using neural networks to augment or replace the ODE integration is a promising area [31].
Figure 1: (a) Left: the original weak rate group of nuclei supported pynucastro 2.0.(b) Right: the original weak rate group of nuclei plus the contribution from the [12-14] rates.

Figure 2 :
Figure 2: (a) Left: the available rate nuclei up to 56 Ni, extracted from REACLIB.(b) Right: the reduced network availability, by selecting a network up to 21 Ne.