A Framework for Exploring Nuclear Physics Sensitivity in Numerical Simulations

We describe the AMReX-Astrophysics framework for exploring the sensitivity of astrophysical simulations to the details of a nuclear reaction network, including the number of nuclei, choice of reaction rates, and approximations used. This is explored by modeling a simple detonation with the Castro simulation code. The entire simulation methodology is open-source and GPU-enabled.


Introduction
Nucleosynthesis is the process of forming new atomic nuclei via nuclear reactions.By studying both stellar and primordial nucleosynthesis, we gain better insights into stellar evolution, the chemical evolution of galaxies, and the element composition of the Universe.Moreover, enormous amounts of energy are released during nuclear fusion reactions in stellar nucleosynthesis, which is the power source leading to thermonuclear explosions such as Type Ia supernovae (SN Ia) [1,2] and X-ray bursts (XRBs) [3].Given the wide range of astrophysical environments in which nucleosynthesis occurs, each characterized by distinct thermodynamic conditions and dynamics, nuclear reaction networks are generally tailored to each case.
Due to limited computational resources, simulations of astrophysical explosions often use small to moderate-sized networks.Generally, the number of isotopes and reaction rates that can be included depends on the dimensionality of a simulation.In 1D, it is possible to incorporate more than a thousand isotopes along with thousands of nuclear reaction rates [4].However, the situation becomes significantly more challenging in 2-and 3D, due to the increased memory and the need to advect each species separately.A contemporary 2D XRB simulation can include 25 isotopes easily [5], while a 3D XRB simulation can be restricted to ∼ 10 isotopes [6,7].
A reaction sensitivity study is essential to assess the reliability of a simulation where reactions provide a significant part of the energy budget.This can explore the significance of each individual reaction rate, the size of the network, the modifications to rates (like screening), and approximations used in the network.For example, through a sensitivity study, we can identify the most critical reaction rates that should be included in the network while eliminating unimportant ones to reduce cost.During this process, we also gain insights into how reaction rates impact the properties of the problem, such as the propagating speed of a thermonuclear flame in an XRB [5] and the light curve of thermonuclear explosions, leading to a better convergence between the numerical models and observations.
Here we present the AMReX-Astrophysics / pynucastro framework for exploring nuclear reaction physics in numerical simulations.The reaction networks are generated via the pynucastro library [8,9], which writes GPU-enabled C++ networks for use in the AMReX-Astrophysics Microphysics library [10].In section 2, we describe the overall procedure and workflow for generating an arbitrary reaction network with its right-hand-side (RHS) and Jacobian using pynucastro.In section 3, we demonstrate a simple sensitivity study using the compressible hydrodynamics code Castro [11,12] to model a 1D 4 He detonation.All of these codes are developed openly and are freely available on GitHub.

Arbitrary Network Generation
The first step in a reaction sensitivity study is the generation of the reaction network.We will use pynucastro for this purpose.pynucastro is a Python library that allows users to interactively explore nuclear reaction rates, construct any arbitrary reaction network, and export the RHS and Jacobian for the network for use in an integrator.The ReacLibLibrary class provides access to the collection of JINA ReacLib [13] rates and the TabularLibrary class provides tabulated weak rates from various sources [14,15].A detailed description of pynucastro can be found in [9], here we only emphasize essential points related to network generation.
To illustrate the process of network generation, we start by explaining the design philosophy and steps involved in building the subch simple network, a network employed in our prior research studies [5,16].Subsequently, we showcase how this methodology can be seamlessly extended to develop more complex networks.
The first step in the network generation process is selecting the nuclei of interest.We will build a network that is appropriate for He burning in a detonation.Both ReacLibLibrary and TabularLibrary provide a method called linking nuclei() that takes a list of nuclei as input and returns a Library object with all possible rates that connect the input nuclei.This is demonstrated in the following code: import pynucastro as pyna # Create a library object that contains all ReacLib rates reaclib_lib = pyna .ReacLibLibrary () # Create a list of nuclei nuc_list = [ " p " , " he4 " , " c12 " , " o16 " , " ne20 " , " mg24 " , " si28 " , " s32 " , " ar36 " , " ca40 " , " ti44 " , " cr48 " , " fe52 " , " ni56 " , " al27 " , " p31 " , " cl35 " , " k39 " , " sc43 " , " v47 " , " mn51 " , " co55 " , " n13 " , " n14 " , " f18 " , " ne21 " , " na22 " , " na23 " ] # Returns a subset of ReacLib Library subch_lib = reaclib_lib .linking_nuclei ( nuc_list ) Next, various fine-tuning procedures can be done on this Library to reduce the overall size of the network.It is often desirable to manually remove rates that have been determined as insignificant for energy generation or nucleosynthesis based on past studies.This can be accomplished using Library built-in method remove rate() by supplying either the Rate object or the short name of the rate in the form of a string: "A(x, y)B".Below we demonstrate this by removing the reverse rates of C and O reactions as well as some links involving heavy nuclei: for r in subch_lib .get_rates (): # Check if the products of the rate are the following if sorted ( r .products ) in [[ pyna .Nucleus ( " c12 " ) , pyna .Nucleus ( " c12 " )] , [ pyna .Nucleus ( " c12 " ) , pyna .Nucleus ( " o16 " )] , [ pyna .Nucleus ( " o16 " ) , pyna .Nucleus ( " o16 " )]]: # Remove the rate subch_lib .remove_rate ( r ) # A list of rates to remove by simply supplying the short -name of the rate rates_to_remove = [ " p31 (p , c12 ) ne20 " , " si28 (a , c12 ) ne20 " , " ne20 ( c12 , p ) p31 " , " ne20 ( c12 , a ) si28 " , " na23 (a , g ) al27 " , " al27 (g , a ) na23 " , " al27 (a , g ) p31 " , " p31 (g , a ) al27 " ] # Remove rates individually from the network for r in rates_to_remove : subch_lib .remove_rate ( r ) We can also directly modify the products of a reaction, to further simplify the network.This is useful, for example, in a neutron-capture sequence like A(B, n)Y(n, γ)C where the neutron capture on Y is significantly faster than the first reaction.In this case, a simple approximation involves modifying the product of the first rate to yield the final product of the two-reaction sequence, C, while retaining the rate of the first reaction and adjusting the overall Q-value accordingly.This is demonstrated in the following code: # Select the reaction rates and end product to modify modify_rates = [( " c12 ( c12 , n ) mg23 " , " mg24 " ) , ( " o16 ( o16 , n ) s31 " , " s32 " ) , ( " o16 ( c12 , n ) si27 " , " si28 " )] for r , mp in modify_rates : # Get the specifc Rate object using the short -name _r = reaclib_lib .get_rate_by_name ( r ) # Modify the product of that rate _r .modify_products ( mp ) # Add the modified rate into the customized Library .subch_lib .add_rate ( _r ) This network describes α-captures up to the iron group (and is the starting point for the subch simple network from [5]).However, for modeling SN Ia, the formation of ironpeak isotopes via electron capture during the explosion is important [17][18][19][20].We can include these by adding tabulated weak rates.The code below obtains a new set of rates from both ReacLibLibrary and TabularLibrary in the iron-peak and combines them with our existing rates to create our custom rate library all lib.Note that while the original set of nuclei did not include neutrons, we do include neutrons in the iron-group reactions.

all_lib = subch_lib + iron_reaclib + iron_weak_lib
There may be duplicated weak rates from the ReacLib library and the tabulated weak rates.To avoid errors when creating the network, these duplicates can be identified and removed via the find duplicated links() method: After we've completed the construction of the customized rate libraries, pynucastro can output functions to compute both the righthand-side (RHS) of the ODE system of evolution equations and the Jacobian matrix.Several formats are possible, each provided by a different class: PythonNetwork for Python, SimpleCxxNetwork for a basic C++ network, and AmrexAstroCxxNetwork for C++ for use in the AMReX-Astrophysics suite.For our simulations, we want an AmrexAstroCxxNetwork, which can be used with Castro, and couples with hydrodynamics through a variety of integration techniques (e.g., our simplified-SDC method [21]).We create a network in this format simply as: # Use previously defined customized library to create our network net = pyna .AmrexAstroCxxNetwork ( libraries =[ all_lib ]) Finally, pynucastro can also approximate rate sequences like A(α, p)X(p,γ)B rates into an effective as A(α, γ)B rate, which is used in the traditional "aprox" [22] networks.This approximation holds valid when the proton-capture rate happens at a much faster timescale compared to the (α, p) rate, causing an overall equilibrium in the proton flow.Similar to modify products(), this approximation reduces the size of the network by removing the intermediate nuclei.This method is demonstrated in the following code: # Perform the (a , p )( p , g ) approximation by specifying the intermediate nuclei net .m ake_ap_pg_approx ( intermediate_nuclei =[ " cl35 " , " k39 " , " sc43 " , " v47 " ])]) # Remove intermediate nuclei since now they ' re approximated out net .remove_nuclei ([ " cl35 " , " k39 " , " sc43 " , " v47 " ])]) Our final network has 36 isotopes and 149 rates including 12 TabularRate rates (compared to the original subch simple with 22 isotopes and 57 rates).Note that subch simple also uses the (α, p)(p, γ) approximation for 35 Cl, 39 K, 43 Sc, 47 V, 51 Mn, and 55 Co.Before writing out the actual RHS and the Jacobian, there are also methods that allow the user to explore the network they created and make adjustments if needed.For example, there is an easy interface to plot each reaction rate with a specified density, temperature, and mass abundances that mimic the simulation conditions to help determine the importance of each rate.An example output of the built-in plotting functions is shown in Figure 1.
Once we're satisfied with the current network, pynucastro can write out all the necessary files, including the RHS and the Jacobian of the network, compatible with Microphysics and Castro with a single function: # Write all files needed by the package , Microphysics , to integrate the network net .write_network ()  The AMReX-Astrophysics Microphysics library contains all of the code needed to support the integration of this network, including an equation of state (we'll use the Helmholtz EOS from [23]), integrators such as VODE [24], as well as routines for neutrino cooling and plasma screening.Microphysics can be used on its own to perform network evolution with self-heating or coupled to a simulation code, like Castro.It provides the necessary interfaces to support several different methods of coupling hydrodynamics and reactions.
Typically, a user creates a dedicated directory in Microphysics under Microphysics/networks (we'll use custom net for the network described above).This directory serves as the repository containing both the network-generating script using pynucastro and all the associated autogenerated files.

Network Sensitivity Testing
We will now show how to do a network sensitivity study similar to the one recently performed for XRBs in [5].Any Castro problem setup can be told to use our new network (in the Microphysics directory hierarchy) at compile time, via: To demonstrate a simple sensitivity testing, we use the Detonation test problem in Castro.Detonation models a 1D propagating detonation wave, offering a way to replicate the thermodynamic conditions of the He detonation wave found in the shell-burning stage of the sub-Chandrasekhar double-detonation model [25], as discussed in previous work [16,21].This sensitivity testing on the Detonation problem helps us choose the appropriate network for a realistic sub-Chandrasekhar double-detonation model, especially when studying nucleosynthesis.
To demonstrate this idea, we will compare subch simple and the classic 19-isotope network, aprox19 [26], with the network we generated in Section 2, which will be referred to as He-C-Fe.Our detonation simulation spans a domain up to x = 4.0 × 10 7 cm, with 256 coarse grid zones.By employing two refinement levels with ratios of 2, we reach 1024 zones at the finest level, providing a resolution of ∼ 0.4 km.A zero-gradient outflow boundary condition is used on both boundaries.The initial thermodynamic conditions are taken from the base of the He layer from the simulations in [16]: a uniform ρ = 1.1 × 10 9 g cm −3 , T pert = 1.1 × 10 9 K from 0 < x pert < 1.2 × 10 7 cm, T = 1.75 × 10 8 K for x > 1.2 × 10 7 cm, and a composition of 99 % 4 He and 1 % 14 N.
Figure 2 shows the time evolution of the temperature and specific nuclear energy generation rate using the He-C-Fe network.The peak temperature and specific energy generation rate reach ∼ 3 × 10 9 K and ∼ 10 21 erg g −1 s −1 , consistent with results of the shell-burning stage in the previous study [16].
We now look at the effect of the network.Figure 3 shows the time evolution of the shock speed for our different networks.It is observed that all models reach a shock speed of ∼ 10 4 km s −1 after ∼ 15 ms, while some noticeable differences for t ≲ 15 ms.Specifically, He-C-Fe and subch simple exhibit similar profiles, with roughly double the shock speed compared to aprox19 for t ≲ 15 ms.This difference in shock speed is due to the missing 12 C(p, γ) 13 N(α, p) 16 O rate in aprox19, resulting in weaker initial carbon-burning [27][28][29].A similar effect was seen with flames in XRBs in the study by [5].
This carbon burning enhancement effect is evident when examining the total mass abundance of 12 C in Figure 4, showing the time evolution of 4 He, 12 C, and iron-peak isotopes available in the network.The enhanced carbon burning path in both He-C-Fe and subch simple leads to a greater 12 C consumption rate, leading to a higher shock speed.After ∼ 15 ms, the total mass of 12 C stabilizes, reaching a steady speed for all three networks.In terms of iron-peak nucleosynthesis, both aprox19 and subch simple indicate that 56 Ni is the primary product, whereas He-C-Fe predominantly yields 58 Ni.This insight is important for determining the correct mass abundance of 56 Ni following the explosion since the beta decay from 56 Ni → 56 Co → 56 Fe primarily shapes the observed SN Ia light curve.
To demonstrate the computational cost associated with adding additional iron-peak isotopes, Figure 5 shows the CPU-hours used with each network.As expected, a larger network size corresponds to a higher CPU hour usage.Specifically, the CPU-hours used by He-C-Fe increased drastically after 15 ms, coinciding with the establishment of the steady propagating detonation wave.In the end, He-C-Fe required nearly 4 times the CPU-hours compared to subch simple, which has ∼ 1/3 of the total number of rates in He-C-Fe.This result suggests the need for further simplification and optimization of the network before employing it in actual simulations such as a 2D sub-Chandrasekhar double detonation model.

Future Development
We've outlined the versatile approach for nuclear reaction sensitivity testing using pynucastro, Microphysics, and Castro.By leveraging pynucastro for network generation, Microphysics for

Figure 1 :
Figure1: A visualization of the custom network that we created using the pynucastro package.The color bar shows the reaction rates with solar composition, ρ = 9×10 7 g cm −3 , and T = 6×10 9 K.The horizontal axis and the vertical axis show the atomic number, Z, and the extra number of neutrons compared to the number of protons, respectively.The gray nodes with dotted gray lines represent the (α, p)(p, γ) approximation, which are not directly in the network.All reactions that have the form, A(X, α)B or A(α, X)B, are hidden for better clarity.

Figure 2 :
Figure 2: Temperature and specific nuclear energy generation evolution for He-C-Fe

Figure 3 :
Figure 3: Time evolution of the shock speed for different networks.

Figure 4 :
Figure 4: Time evolution of the total mass of 4 He,12 C and iron-peak isotopes over time.

Figure 5 :
Figure 5: Total number of CPU hours used for each network.