Practical Effects of Integrating Temperature with Strang Split Reactions

For astrophysical reacting flows, operator splitting is commonly used to couple hydrodynamics and reactions. Each process operates independent of one another, but by staggering the updates in a symmetric fashion (via Strang splitting) second order accuracy in time can be achieved. However, approximations are often made to the reacting system, including the choice of whether or not to integrate temperature with the species. Here we demonstrate through a simple convergence test that integrating an energy equation together with reactions achieves the best convergence when modeling reactive flows with Strang splitting. Additionally, second order convergence cannot be achieved without integrating an energy or temperature equation.


INTRODUCTION
Simulations of stellar flows require solving the equations of hydrodynamics coupled with a nuclear reaction network. The equations of hydrodynamics with reacting sources are: where ρ is the density, U is the velocity vector, X k are the species mass fractions with creation ratesω k , p is the pressure, E is the specific total energy, andṠ is the nuclear energy generation rate. When we are reacting, we can look at internal energy where e is the specific internal energy or alternately, we can evolve the temperature, T , where c v = ∂e/∂T | ρ (this form neglects composition changes, see Almgren et al. 2008).

NUMERICAL METHODOLOGY
We use the freely-available Castro simulation code (Almgren et al. 2010;Almgren et al. 2020) to solve the equations of hydrodynamics, using an unsplit piecewise parabolic method coupled with reactions. Castro uses either Strang splitting or spectral deferred corrections (SDC) to couple the hydrodynamics and reactions (Zingale et al. 2019). Here we focus on the Strang splitting.
In a Strang split evolution (Strang 1968), we update the full hydrodynamics state, U , as: where R ∆t/2 is the reaction update through a timestep ∆t/2 and A ∆t is the advective update through ∆t. We see with this splitting, each process operates on the state left behind by the previous operation, and the staggering of the physics is done to give second order accuracy in time.
During reactions, we neglect the hydrodynamics terms, so the reactive system updates according to: There are several different approaches taken in the literature to this operator-split reacting system, including some approximations that make integrating the reaction system computationally less expensive: • Evolve (X k ) only. This neglects temperature evolution completely in the burn, only evolving Eq. 9. This is the method used in Fryxell et al. (2000).
• Evolve (X k , T ). This is used in Pakmor et al. (2012) and García-Senz et al. (2013). To avoid expensive equation of state calls in getting the specific heat, we can optionally "freeze" the value of c v at the start of integration. This was discussed in Bell et al. (2004) and until recently was the default method in Castro.
• Evolve (X k , e), and get T from e using the equation of state. This was discussed in Fryxell et al. (1989) and is the current default method in Castro. Raskin et al. (2010) also propose a hybrid system where the first approach is used in most cases, switching to the second approach only near NSE. For all of these cases, density is constant during the reaction operation.
Depending on how vigorous the burning is and how much the temperature changes during a hydrodynamic timestep, one or more of these methods may be reasonable. For explosive reactions, we expect that evolving the full system will be required. The goal of this note is to try to quantify the convergence of a reacting hydrodynamics code with these different approximations. In Zingale et al. (2019), we introduced a test problem where we could measure the convergence of a reacting hydrodynamics problem via Richardson extrapolation-this was an acoustic pulse with helium burning via 3-α and 12 C(α, γ) 16 O. Initially, the domain is pure 4 He, but both 12 C and 16 O are created as time evolves. The published tests showed that we can get overall 4th order in space and time convergence with SDC coupling. Here we run the same test with Strang coupling, looking at the various approaches at incorporating a temperature / energy equation in the reactive portion of the update.
For each method, we run the Castro reacting_convergence test problem at 5 resolutions: 64 2 , 128 2 , 256 2 , 512 2 , and 1024 2 , with the timestep kept fixed in proportion to the grid resolution. We then compute four errors between adjacent resolutions by coarsening the finer resolution run, and computing the L 1 norm over all zones. Figure 1 shows the norm of the error vs. the coarse run resolution. The slope of these lines is a measure of the convergence rate (Oberkampf & Roy 2010). We see that all methods converge at least second order for density, but for the thermodynamic quantities, (ρe) and T , the method where only X k is evolved during reactions has larger errors and much worse convergence than the other methods. Looking at the trace nuclei generated in the burning, 12 C and 16 O, we see a large difference between the two methods that evolve some sort of energy and the one method where only X k is evolved-the latter converging essentially first order at high resolution.

SUMMARY
Looking at global convergence of a reacting hydrodynamics problem we see that second order convergence is only realized when temperature or energy is evolved alongside the composition when using a Strang-split approach to reactions. This is just a single, rather simple problem, but this suggests that reactive hydrodynamics simulations should switch to integrating temperature or another energy equation together with reactions to yield better overall convergence and accuracy. This complements the work of Müller (1986) which showed that when evolving near nuclear statistical equilibrium, evolving entropy with the system is needed for stability. We expect that for problems with vigorous burning, such as detonations, directly coupling the composition and thermodynamic evolution will be especially important.  NumPy (Oliphant 2007), VODE (Brown et al. 1989)