Skeletal Kinetics Reduction for Astrophysical Reaction Networks

A novel methodology is developed to extract accurate skeletal reaction models for nuclear combustion. Local sensitivities of isotope mass fractions with respect to reaction rates are modeled based on the forced optimally time-dependent (f-OTD) scheme. These sensitivities are then analyzed temporally to generate skeletal models. The methodology is demonstrated by conducting skeletal reduction of constant density and temperature burning of carbon and oxygen relevant to Type Ia supernovae (SNe Ia). The 495-isotopes Torch model is chosen as the detailed reaction network. A map of maximum production of 56Ni in SNe Ia is produced for different temperatures, densities, and proton-to-neutron ratios. The f-OTD simulations and the sensitivity analyses are then performed with initial conditions from this map. A series of skeletal models are derived and their performances are assessed by comparison against currently existing skeletal models. Previous models have been constructed intuitively by assuming the dominance of α-chain reactions. The comparison of the newly generated skeletal models against previous models is based on the predicted energy release and 44Ti and 56Ni abundances by each model. The consequences of y e ≠ 0.5 in the initial composition are also explored where y e is the electron fraction. The simulated results show that 56Ni production decreases by decreasing y e as expected, and that the 43Sc is a key isotope in proton and neutron channels toward 56Ni production. It is shown that an f-OTD skeletal model with 150 isotopes can accurately predict the 56Ni abundance in SNe Ia for y e ≲ 0.5 initial conditions.


INTRODUCTION
Direct implementation of detailed reaction networks (RNs) containing many isotopes and reactions in hydrodynamic flow solvers is computationally very expensive.In most cases, it is unavoidable to use efficient reaction kinetics models to conduct large-scale hydrodynamic simulations pertaining to astrophysical explosions, such as Type Ia supernovae (SNe Ia).These reaction kinetics models, which are usually extracted from detailed reaction networks, should reasonably estimate the released energy and isotope abundances.Integration of the ordinary differential equations representing the abundance levels of a set of isotopes of reacting nuclei in the continuum limit serves two functions in stellar models.The primary function of hydrodynamics is to provide the magnitude and sign of the nuclear energy generation rate (Weaver et al. 1978;Bravo et al. 2019;Arnould & Goriely 2020).This is usually the largest energy source in regions conducive to nuclear reactions, and its accurate determination is essential for stellar models.These models usually require accurate predictions of the energy generated by nuclear burning over a wide range of temperatures, densities, and compositions (Timmes 1999;Timmes et al. 2000;Röpke & Sim 2018).The other function is to describe the evolution of the composition.In some stellar events, the isotopic abundances themselves are of primary interest for understanding the origin and evolution of the chemical elements (Pagel 2009;Matteucci 2012;Kobayashi et al. 2020).Moreover, matching observational evidence of certain isotopes, e.g. 56Ni in supernova light curves, gives confidence in the underlying computational model (Seitenzahl & Townsley 2017;Bora et al. 2022).
Thousands of isotopes can participate in a reaction network during a stellar phenomenon (Wanajo 2013;Nishimura et al. 2015;Fernández et al. 2017;Lippuner & Roberts 2017;Psaltis et al. 2022).Accurate predictions of the nuclear energy generation rate and the composition changes in such RNs is computationally expensive.The largest block of memory in a stellar model is usually used for storing the isotopic abundances of every computational cell at all time steps.Even with modern methods for solving stiff systems of ordinary differential equations, integration of the evolution equations of the isotopic abundances dominates the total cost of a stellar model when the number of isotopes evolved is O(100) (Arnett 1996;Nouri et al. 2019).To decrease the computational cost, one has to make a choice between having fewer isotopes (order reduction) or less spatial resolution (or mass resolution).The general response to this trade-off has been the order reduction by using simplified RNs within hydrodynamic solvers to calculate an approximate energy generation rate and isotope mass fractions during stellar explosions (Röpke & Sim 2018).As a post-processing step, the detailed nuclear composition of the ejecta is computed using a large RN by employing Lagrangian tracer particles (Thielemann et al. 1986;Townsley et al. 2016;Leung & Nomoto 2018;Bravo et al. 2019;Seitenzahl & Pakmor 2023).These particles represent passive mass elements and can be evolved in situ, within the simulation, with time steps dictated by the hydrodynamic solvers, or off-site, using an additional recontruction step based on the snapshot data (Sieverding et al. 2023).Current large simulations can use at least 10 6 tracer particles to calculate detailed 3D spatial composition (Seitenzahl & Pakmor 2023).
In chemical combustion, recent advances in data-driven techniques and sensitivity analysis have opened the possibility of significantly enhancing the efficiency and flexibility of generating reduced RNs (Lu & Law 2009).Recent skeletal models, i.e. optimized subsets of detailed reaction networks, can be prepared in an optimized and automated manner, with consistent accuracy throughout the evolution of the network, and can be adapted based on the availability of computing resources.Utilizing such capabilities for nuclear combustion can address certain issues regarding the model reduction of nuclear RNs, e.g.scenarios where the experience required to generate quality reduced models might be lacking.There is a history of utilizing sensitivity analysis in nuclear combustion.This includes studies of the Big Bang (Beaudet & Reeves 1983;Delbourgo-Salvador et al. 1985;Krauss & Romanelli 1990;Smith et al. 1993;Nollett & Burles 2000;Cyburt 2004), stellar explosions (Hix et al. 2003;Parikh et al. 2008;Longland et al. 2010;Longland 2012;Bravo & Martinez-Pinedo 2012;Bliss et al. 2020), and the r-process (Mumpower et al. 2012(Mumpower et al. , 2015;;Sprouse et al. 2020;Barnes et al. 2021).Most of these contributions are based on direct Monte Carlo simulations and their focus is on understanding the impact of nuclear reaction rate and/or other nuclear uncertainties on the resulting nucleosynthesis predictions.
The goal of this work is to develop skeletal RNs suitable for situations pertaining to SNe Ia.A skeletal model is a subset of a detailed reaction model which is generated by eliminating unimportant isotopes and reactions (Smooke (1991); Peters & Rogg (1993); Stagni et al. (2016); Li et al. (2020)).The skeletal reduction is usually the first step in developing a model reduction.The next steps in the reduction include time-scale analysis techniques, e.g.quasi steady state approximation (Stiefenhofer (1998); Girimaji & Ibrahim (2014)), partial equilibrium approximation (Rein (1992); Goussis (2012)), and rate controlled constrained equilibrium (Keck (1990); Hadi et al. (2016)) amongst others.To develop skeletal RNs for SNe Ia, first local sensitivities of isotope mass fractions with respect to reaction rates are analyzed during the constant density and temperature (constant-ρT ) burning of carbon and oxygen with different initial conditions.The isotopes are ranked based on their sensitivities (importance), and several sets of skeletal models with different levels of accuracy are generated by selecting different numbers of important isotopes.The sensitivities are computed by the forced optimally time-dependent (f-OTD) methodology (Donello et al. 2022).This is an on-the-fly reduced order modeling (ROM) technique, recently introduced for computing sensitivities in evolutionary dynamical systems.Unlike the traditional ROM techniques, the f-OTD does not require any offline data generation, and all the computations are carried out online.Nouri et al. (2022) and Liu et al. (2024) conducted a similar sensitivity-based skeletal kinetics reduction technique for chemical combustion which automatically eliminates unimportant reactions and species.Time-dependent f-OTD modes are able to capture sudden transitions associated with the largest finitetime Lyapunov exponents (Babaee et al. 2017b).Time-dependent bases have also been used for stochastic reduced order modeling (Sapsis & Lermusiaux 2009;Cheng et al. 2013;Babaee et al. 2017a;Babaee 2019;Patil & Babaee 2020) and on-the-fly reduced order modeling of reacting species transport equation (Ramezanian et al. 2021;Aitzhan et al. 2022).The f-OTD can be formulated as a special case of the dynamical low-rank approximation (Koch & Lubich (2007)).The specific objectives here are i) to introduce the f-OTD technique for computing sensitivities for a nuclear combustion system, and ii) to find skeletal models for thermonuclear burning in SNe Ia.The first set of skeletal models are applicable to both neutron-rich and equal numbers of neutron and proton scenarios.This is facilitated by the f-OTD skeletal reduction technique, without a priori assumptions or expertise, e.g. the assumption of an equal number of protons and neutrons.Section 2 briefly presents the theory behind the f-OTD method for constant-ρT burning in SNe Ia and the automatic process of eliminating unimportant isotope/reaction from a detailed RN.This elimination process is explained in Section 3 with a simple example, starting from a RN with 21 isotopes and reducing it to a skeletal model with 10 isotopes.Section 4 describes the application of the f-OTD skeletal reduction method to the Torch RN1 (Timmes 1999;Timmes & Swesty 2000;Anninos et al. 2019).This RN considers 495 isotopes, up to 91 Tc, and 6012 reactions.Different skeletal models with different levels of accuracy are extracted and their ability to reproduce the energy release and 44 Ti and 56 Ni abundances of the Torch model are analyzed.Section 5 provides the concluding remarks.
2. SKELETAL REDUCTION WITH F-OTD METHOD 2.1.Reduced-order modeling of the sensitivity matrix with f-OTD For the model description, let isotope k have total charge Z k and atomic weight A k .Let the aggregate total of isotope k have a mass density ρ k and a number density n k in a material with the temperature T and the total mass density ρ.The mass fraction of isotope k is defined as x k = ρ k /ρ = n k A k /ρN A where N A = 6.02252 × 10 23 particles/mol is the Avogadro's number.The mean atomic weight is A m = ( x k /A k ) −1 , and is the equivalent of the mixture molar mass from the combustion literature (Williams 1985).The nuclear abundance of isotope k is and the electron abundance, or electron number fraction, is y e = Z m /A m .This is related to the neutron excess, η, by η = 1 − 2y e , so that η = 0 corresponds to y e = 0.5.The total scalar pressure, the total specific internal energy, and the total specific entropy are denoted by p tot , e tot , and s tot , respectively (Nouri et al. 2019).Other quantities such as the specific heats or adiabatic indices can be determined via an equation of state (Timmes & Arnett 1999) once the partial derivatives of the pressure and the specific internal energy with respect to the density and temperature are known.Consider a nuclear system of n s isotopes reacting through n r reactions: where M k is a symbol for isotope k, and ν ′ kj and ν ′′ kj are the stoichiometric coefficients of isotope k in reaction j.Changes of abundances y = [y 1 , y 2 , . . ., y ns ] T in constant-ρT burning within a carbon-oxygen white dwarf (WD) with constant temperature and pressure can be described by the following initial value problem (Timmes et al. 2000): where t ∈ [0, t f ] is time, t f is the final time, and α = [1, 1, . . ., 1] ∈ R nr is the vector of perturbation parameters.In Eq. ( 2), Q j is the progress rate of reaction j and is equal to the following for one and two body reactions: Here, K f,j and K r,j are the forward and reverse rate of reaction j.The quantity kj is the total possible number of elementary reactions per unit volume (obtained by counting the number of possible collisions) and is based on the number density of particles.Because the collision energies are well below the Coulomb barrier, most collisions do not result in nuclear reactions.Thus, the reaction rate is the product of the collision rate and the tunneling probability.
The abundances in Eq. ( 2) are perturbed by infinitesimal variations of α j , by letting α j = 1 + δα j , where δα j ≪ 1 for j = 1, 2, . . ., n r .The perturbation with respect to α j amounts to an infinitesimal perturbation of progress rates Q j .The sensitivity matrix, S(t) = [s 1 (t), s 2 (t), ...s nr (t)] ∈ R ns×nr , contains local sensitivity coefficients, s j = ∂y/∂α j , and can be calculated by solving the sensitivity equation (SE), where L im = ∂fi ∂ym and F ij = ∂fi ∂αj are the Jacobian and the forcing matrices, respectively.The model reduction strategy is based on selecting reactions, whose perturbations grow most intensely in the abundance Eq. ( 2).The selection of important reactions is performed by instantaneous observation of modeled sensitivities.In f-OTD, the sensitivity matrix S(t) is modeled by factorizing it into a time-dependent subspace in the n s -dimensional phase space of abundances represented by a set of f-OTD modes: U (t) = [u 1 (t), u 2 (t), . . ., u r (t)] ∈ R ns×r .These modes are orthonormal u T i (t)u j (t) = δ ij at all t, where δ ij is the Kronecker delta.The rank of S(t) ∈ R ns×nr is d = min{n s , n r } while the f-OTD modes represent a rank-r subspace, where r ≪ d.To this end, the sensitivity matrix is approximated via the f-OTD decomposition (Fig. 1) as This decomposition is a low-rank approximation of the sensitivity matrix S(t).Therefore, U (t)V (t) T closely approximates S(t) and it is not exact.Both U (t) and V (t) are time-dependent, and their explicit time dependency on t is dropped for brevity.Figure 1 shows the schematic of the decomposition of S into f-OTD components U and V .The evolution equation for U and V are obtained by substituting the sensitivity decomposition (S(t) ≈ U (t)V T (t)) into Eq.( 4): where Q = I −U U T is the orthogonal projection onto the space spanned by the complement of U and C = V T V ∈ R r×r is a correlation matrix.Matrix C(t) is, in general, a full matrix implying that the f-OTD coefficients are correlated.
L r = U T LU ∈ R r×r is a reduced linearized operator.Equation ( 5) represents a coupled system of ODEs and constitutes the f-OTD evolution equations.The f-OTD modes align themselves with the most instantaneously sensitive directions of the abundance evolution equation when perturbed by α.It is shown by Babaee et al. (2017b) that when α is the perturbation to the initial condition, the OTD modes converge exponentially to the eigen-directions of the Cauchy-Green tensor associated with the most intense finite-time instabilities.The primary computational advantage of using f-OTD is that the method only evolves two skinny matrices containing (n s + n r ) × r elements as opposed to n s × n r elements in the SE (Eq.( 4)).This reduces the required memory for ODE solvers drastically and facilitates the application of stiff solvers for evolving sensitivities.Moreover, in the f-OTD decomposition, the sensitivities are stored in the compressed form, i.e., matrices U , and V are kept in the memory as opposed to their multiplication U V T , i.e., the decompressed form.Therefore, in comparison to the full SE, f-OTD decomposition results in the memory compression ratio of (n s × n r )/((n s + n r )r).

Identification of important reactions & isotopes
In f-OTD skeletal reduction, modeled sensitivities are computed in a factorized format by solving Eqs. ( 2), (5a), and (5b), and the values of U , V , and y are stored at resolved time steps 2) is initialized with different sets of isotope abundances, temperature and density within their ranges of interest.Each simulation with a different initial condition is identified as a case.At each resolved time step and for each case, the eigen decomposition of S T S ∈ R nr×nr is computed as AΛA T , and the vector w = (Σλ i |a i |)/(Σλ i ) ∈ R nr is basically the average of eigenvectors of S T S matrix weighted based on their associated eigenvalues (λ i ), and prevents dealing with each eigenvector (a i ) separately.Each component of w, i.e. w i , is positive and associated with a certain reaction (ith reaction).The larger the w i value, the more important the reaction i is.The w max,i denotes the highest value of w i through all resolved time steps and cases.The elements of w max vector are sorted in descending order to find the indices of the most important reactions in the detailed model.Isotopes are also sorted based on their first presence in the sorted reactions, i.e. isotopes which first show up in a higher ranked reaction would be more important than an isotope which first participates in a lower ranked reaction.This results in a reaction and isotope ranking based on w max vector.Finally, a set of isotopes are chosen by defining a threshold ϵ on the element of w max vector and eliminating isotopes whose associated w max,i are less than ϵ.This terminates reactions which include the eliminated isotopes from the detailed model.Since the model reduction is reaction based, any non-reactive isotope with a non-zero mass fraction in the initial condition must be manually added to the skeletal model.
In summary, the nuclear combustion system is instantaneously observed, and the reactions are sorted based on their effects on sensitivities to find the most important isotopes.These isotopes and the reactions which connect them together create the skeletal models.In combustion systems, perturbations with respect to "fast" reactions generate very large sensitivities for short time periods which vanish as t → ∞.On the other hand, perturbations with respect to "slow" reactions generate smaller and more sustained sensitivities.The approach here is based on the instantaneous observation of sensitivities, both slow and fast reactions can leave an imprint on the instantaneous normalized reaction vector (w) if their imprints are larger than the threshold value (ϵ).However, if the sensitivities associated with fast and slow reactions from pre-determined times and locations are combined with each other before dimension reduction, as commonly done in principal component analysis (PCA) type schemes, the smaller sensitivities associated with slow reactions would be out-weighted by the large sensitivities associated with fast reactions.

SKELETAL REDUCTION ON THE APPROX21 RN
The process of eliminating unimportant reactions and isotopes from a kinetics model with f-OTD is demonstrated in this section with a simple example.Let us start with the Approx21 model which is the default MESA network (Paxton et al. 2011(Paxton et al. , 2015) ) for alpha chain reactions.Approx21 evolves n s = 21 isotopes: n, p, 1 H, 3 He, 4 He, 12 C, 14 N, 16 O, 20 Ne, 24 Mg, 28 Si, 32 S, 36 Ar, 40 Ca, 44 Ti, 48 Cr, 56 Cr, 52 Fe, 54 Fe, 56 Fe, and 56 Ni through n r = 112 reactions.In this RN, (α,p)(p,γ) and (γ,p)(p,α) links are included in order to obtain reasonably accurate energy generation rates and abundance levels when the temperature exceeds 2.5e9 K.At these elevated temperatures, the flows through the (α,p)(p,γ) sequences are faster than the flows through the (α,γ) channels.An (α,p)(p,γ) sequence is, effectively, an (α,γ) reaction through an intermediate isotope.By assuming steady-state proton flows through intermediate isotopes or intermediate isotopes2 .The skeletal reduction is exemplified here by analyzing only one case for the constant-ρT burning of a mixture of carbon and oxygen in a WD progenitor.The ignition initiates with T 9 = 3, ρ 9 = 1, and V U Figure 2. Model reduction for Approx21: eigenvalues of the S T S matrix.The sensitivity matrix, S, initially evolves exactly via the SE (Eq.( 4)), and then evolves approximately with the f-OTD equations (Eq.( 5)).The f-OTD simulation (t > 10 −9 ) only evolves r modes of the full order model (t < 10 −9 ).Ignition data is gathered from the constant-ρT burning simulation case described in Fig. 3.
x C,0 = x O,0 = 0.5, so y e = 0.5.Here, T 9 ≡ T /(10 9 K) and ρ 9 ≡ ρ/(10 9 g.cm −3 ).The sensitivities evolve in two phases.In the first, the full-dimensional SE (Eq.( 4)) is solved for a very small time period, e.g. until t = 10 −12 s, to generate the initial conditions for the f-OTD simulation.In the second phase, the sensitivity matrix (S) is approximated by evolving the f-OTD equations (Eq.( 5)).The U and Y matrices are initialized by eigenvalue decomposition of the full sensitivity matrix (S) at the end of the first phase.Figure 2 shows the evolution of the eigenvalues of S T S matrix.It is indicated that i) λ 1 is an order of magnitude larger than λ 2 most of the time, and ii) the modeled sensitivities converge by adding more modes.Therefore, f-OTD simulation with r = 1 provides a reasonable estimation of sensitivities, which is enough if the final goal is to determine the importance of reactions/isotopes.

SKELETAL REDUCTION ON THE TORCH RN
In SNe Ia, the carbon and oxygen burn together to produce nuclei from silicon to the iron peak which are being ejected into the interstellar medium (Nomoto 1997;Woosley et al. 2002;Lippuner & Roberts 2017;Johnson 2019).As this type of supernova does not produce a significant amount of free neutrons, it does not synthesize elements beyond the iron peak (Johnson 2019).Nevertheless, some heavier isotopes can contribute to reactions involving isotopes important to light curves observations, such as 56 Ni.The 495-isotope version of the Torch RN (Paxton et al. 2015) is chosen to represent the detailed RN in this section.The Torch RN extends from 1 H to 91 Tc (Timmes 1999) with  n s = 495 and n r = 6012.The weak reactions are turned on, and no screening is performed on the reaction rates (Fuller et al. 1985).The Helmholtz EOS as developed by Timmes & Swesty (2000) is used with Coulomb correction to calculate the internal energy and the pressure.Several skeletal RNs have been previously proposed based on the Torch model for inline calculations (Timmes 1999;Timmes & Swesty 2000;Anninos et al. 2019) and are used in the MESA code (Paxton et al. 2011;Anninos et al. 2019;Paxton et al. 2015).These include a bare minimum model of the α-chain reactions using 13 isotopes, a 19-isotope RN to also accommodate some hydrogen burning (Weaver et al. 1978), and a 21-isotope RN that adds 56 Cr and 56 Fe and respective equilibrium reaction sequences to the 19 isotope network to attain a lower y e value for pre-supernova models (Paxton et al. 2015).Several important isotopes are produced in burning scenarios with y e significantly lower than 0.5.For example, Woosley (1997) suggests that 48 Ca can only be produced in nature in a subset of SNe Ia, with y e in the range 0.41 to 0.42 and high burning density.The Torch RN covers such scenarios.On the other hand, the performances of the existing skeletal models have not been systematically examined to cover both y e = 0.5 and y e < 0.5 scenarios.This is addressed here by considering initial conditions with different y e values and choosing the 21-isotope RN skeletal model (hereafter denoted Approx21) for comparison with the proposed skeletal RNs.
The SNe Ia progenitor population and burning scenarios cover a wide range of temperatures and densities (Hillebrandt & Niemeyer 2000) so that, most likely, a single skeletal model with a limited number of isotopes (n s ) cannot yield accurate predictions over the full range of conditions which would be covered by a detailed RN.Because of the importance of 56 Ni in SNe Ia, the skeletal models in this work are designed to predict the evolution of x56 Ni correctly.For this purpose, a map of maximum production of this isotope is produced during the course of constant-ρT burning in SNe Ia as shown in Fig. 5.The Torch RN is run 10000 times on a 100×100 grid of T 9 and ρ 9 for y e,0 values of 0.4955 and 0.5.It is observed that 56 Ni is only significantly produced within certain ranges of the temperature and the density values, with a noticeable shift of the 56 Ni production near the peak.In particular, the peak occurs around T 9 = 4.0 at lower densities.Moreover, the maximum production of 56 Ni is decreased by decreasing y e,0 (Iliadis 2015).To generate skeletal reactions, simulations are conducted of constant-ρT burning in SNe Ia for 24 cases, with T 9 ∈ {2, 4, 6}, ρ 9 ∈ {0.001, 0.01, 0.1, 1.0}, and y e,0 ∈ {0.4955, 0.5}.To better show how these cases are distributed, a 2D version of Fig. 5 is shown in Fig. 6.Each case (shown as a red circle in Fig. 6) denotes a set of initial condition for the density, the temperature, the energy and isotope mass fractions.These conditions cover the burning stage, and are of interest in multidimensional simulations (Fryxell et al. 2000;Woosley et al. 2007).The final time for each case is when the mass fraction of 56 Ni reaches its maximum.Figure 7 shows the evolution of isotope mass fractions for two different cases.The red lines in Fig. 7 show the 56 Ni mass fraction and blue triangles denote the final time of f-OTD simulations and sensitivity analysis.The initial mass fractions of 12 C, 16 O, and 22 Ne isotopes in cases with y e,0 = 0.5 are x C,0 = x O,0 = 0.5, and in cases with y e,0 = 0.4955 are x C,0 = 0.45, x O,0 = 0.45, and x N e,0 = 0.1.Note that for the conditions considered here, the sharp decline in 56 Ni mass fraction at t ≳ 10s is likely due to electron capture, causing the network composition to shift to more neutron-rich isotopes of nickel and iron, which is not really relevant to carbon/oxygen burning and therefore not used in the f-OTD analysis.The f-OTD method is used to model the sensitivity matrix with r = 1 mode, and the generated skeletal models by sensitivity analysis are denoted by f-OTD-n in which "n" identifies for the number of isotopes.The predictive capabilities of f-OTD models are compared against those obtained via Torch RN and the Approx21.
A ranking is provided of the first 150 important isotopes in Appendix A considering the maximum characteristic value associated with each isotope, i.e. w max,i .Different skeletal models can be generated by applying a threshold ϵ on w max and eliminating isotopes and their associated reactions with w max,i < ϵ from Torch RN.A comprehensive skeletal model capable of reproducing the energy and isotope mass fraction predictions of Torch RN with a certain accuracy can be developed by specifying an acceptable error level, e.g.5%, in energy or 56 Ni mass fraction estimations.Figure 8(a) shows the isotopes in Torch selected for the f-OTD-150.Each square belongs to an isotope, and isotopes with darker colors are ranked higher (selected earlier) than isotopes with lighter colors.The results show at least 114 isotopes are required to accurately produce 56 Ni with a f-OTD skeletal model.This is because of the importance of some isotopes with 21 ≤ N, Z ≤ 24, and especially 43 Sc in bridging between lighter and heavier isotopes and two QSE clusters (Iliadis 2015;Subedi et al. 2020).It is shown in Fig. 8(a) that 43 Sc plays a key role for proton and neutron channels.Figures 8(b) and 8(c) show the isotopes used in the Approx21 model and SK55, a skeletal model used by Townsley et al. (2019).
The Torch RN is user friendly and flexible to work with any subset of its own isotopes.This means that by providing a ranked list of isotopes from Appendix A, e.g.

CONCLUSIONS
A systematic method for skeletal model reduction of nuclear reaction networks is developed for generating models for the carbon-oxygen combustion in SN Ia covering a range of temperatures and electron number fraction, y e = Z m /A m .In this method, the sensitivities of abundances with respect to reaction rates are modeled using the forced optimally time-dependent (f-OTD) method and are analyzed instantaneously.This results in reaction and isotope rankings based on the correlations between their sensitivities.A key feature of this approach is that it factorizes the sensitivity matrix into a multiplication of two low-ranked time-dependent matrices which evolve based on evolution equations derived from the governing equations of the system.The generated skeletal models are comparatively assessed based on their ability to predict the energy and mass fractions.In particular, the skeletal models as derived here are the first to address situations covering both y e = 0.5 and y e < 0.5.To employ any of the skeletal models developed in this work or to create new f-OTD skeletal models with different numbers of isotopes, one only needs to feed a list of more than 114 ranked isotopes from Appendix A into the Torch RN3 .Further reduction in the number of isotopes (e.g., by using equilibrium assumptions) is a potential future follow-up to this work.
The overall costs of generating an f-OTD model depend on solving the mass fraction equations (which is roughly the same as one non-f-OTD Torch simulation, c torch ), the U equation, whose cost scales as r × c torch , and V equation, whose cost scales as r × n r /n s × c torch .For example, the total cost to generate the f-OTD-150 model with r = 1 modes, for which 10 cases with different densities and temperatures were used, was ∼ 120c torch , which would be negligible, for example, compared to a Monte Carlo simulation with ∼ 10, 000 trials.On the other hand, f-OTD-150 should run ∼ 3.3 faster than Torch.The skeletal reduction technique as described here can be readily extended to other situations.For example, lower values of y e and higher densities to examine the production of 48 Ca, 50 Ti, or 54 Cr.It can also be applied to more complex RNs to examine the production of heavier elements in core-collapse supernovae.With respect to the development of the methodology itself, it can be extended by including the sensitivity analysis based also on transport properties, or even the equation of state (Nouri et al. 2019).Most importantly, as shown recently (Donello et al. 2022), the f-OTD methodology can be used for solving PDEs for multi-dimensional combustion problems in a cost-effective manner -by exploiting the correlations between the spatiotemporal sensitivities of different species with respect to  different parameters.This analysis can be especially insightful for problems containing rare events by providing more insights into global phenomena.

Figure 1 .
Figure 1.Modeling sensitivity matrix S(t) as a multiplication of two low-ranked matrices U (t) and Y (t) which evolve according Eq. (5).Reprinted from Nouri et al. (2022) with permission.

Figure 4 .
Figure 4. Model reduction for Approx21: evolution of isotope mass fractions based on Approx21, Approx13, and f-OTD-10.The last two models are generated from Approx21.

Figure 5 .
Figure 5. Maximum mass fraction of 56 Ni during constant-ρT burning of SNe Ia.Red circles show the f-OTD cases.

Figure 6 .
Figure 6.Maximum mass fraction of 56 Ni during constant-ρT burning of SNe Ia.Red circles show the f-OTD cases.
16 O (1 st rank) to 45 Ca (150 th rank), one can use a f-OTD generated network with the associated isotopes.Figures 9 & 10 show the performance of the f-OTD models in reproducing the energy, mass fractions of 12 C, 44 Ti, and 56 Ni, and y e .The radioactive decays of 44 Ti, and 56 Ni have significant observational applications, and the production of these two isotopes is sensitive to the temperature, density, and y e evolution(Magkotsios et al. (2010)).The energy in Fig.10is normalized by its initial value.The predictions via the Approx21 and SK55 are also presented.The second and third row of subfigures correspond to situations exactly similar to the test cases, but the first and last rows portray the estimations for arbitrary situations within the initial condition domain, i.e., T 9 ∈ [2, 6], ρ 9 ∈ [0.001, 1.0], and y e,0 ∈ [0.4955, 0.5].It is apparent that the f-OTD models with n s ≥ 150 exactly predict the energy evolution of Torch RN and 56 Ni mass fraction within their designed ρ 9 -T 9 -y e ranges.The Approx21 and SK55 models usually over-predict the maximum 56 Ni mass fraction, and Approx21 cannot be used when y e,0 ̸ = 0.5.Replacing Torch with f-OTD models using 114-150 isotopes yields compression ratios ranging from 3.3 to 4.3.Figure11compares the evolution of mass fractions as predicted by Torch model, without any approximation, and nuclear statistical equilibrium (NSE) assumption.The NSE results are generated by using the instantaneous values of ρ 9 -T 9 -y e extracted from non-NSE calculations.The Torch and f-OTD-150 are used as the base reaction network for NSE estimations.It is shown in Figs.11(a) & (b) that only at late time for T 9 = 5 the NSE and non-NSE mass fraction predictions of 44 Ti and 56 Ni (but not 12 C) are close to each other.

Figure 7 .
Figure 7. Evolution of isotope mass fractions in Torch with different initial conditions, highlighting different scenarios for 56 Ni evolution.The blue lines indicate the final time of the f-OTD simulations and sensitivity analysis.The red and gray lines represent the 56 Ni mass fraction and the 10 isotopes with the highest final mass fractions observed during the portrayed ignition process.These isotopes are listed in the left-top corner of each sub-figure

Figure 8 .
Figure 8.(a) The first 150 ranked isotopes of Torch RN sufficient for exact calculation of the maximum 56 Ni abundance and energy in SNe Ia.The darker the squares, the more important (higher ranked) isotopes and the empty squares correspond to eliminated isotopes from the Torch RN.The 43 Sc is a key isotope in proton and neutron channels toward 56 Ni production.The f-OTD models with less than 114 isotopes (Appendix A) do not contain 43 Sc and do not produce the correct amount of 56 Ni.(b) The isotopes for Approx21 and (c) SK55.

Figure 10 .
Figure 10.Model reduction on Torch RN: energy and ye estimations via Torch RN, Approx21, SK55, and f-OTD generated models for four different initial conditions of xc12, xo16, T9, and ρ9.

Figure 11 .
Figure 11.Comparison between mass fractions as predicted by Torch model (black dotted lines) without any equilibrium assumptions and NSE mass fractions estimated based on Torch and f-OTD-150 models (green solid and red dashed lines).NSE and non-NSE simulations have same ye at each time.
Rank Isotope wmax Rank Isotope wmax Rank Isotope wmax Rank Isotope wmax