Entropy-Conserving Scheme for Modeling Nonthermal Energies in Fluid Dynamics Simulations

We compare the performance of energy-based and entropy-conserving schemes for modeling nonthermal energy components, such as unresolved turbulence and cosmic rays, using idealized fluid dynamics tests and isolated galaxy simulations. While both methods are aimed to model advection and adiabatic compression or expansion of different energy components, the energy-based scheme numerically solves the nonconservative equation for the energy density evolution, while the entropy-conserving scheme uses a conservative equation for modified entropy. Using the standard shock tube and Zel'dovich pancake tests, we show that the energy-based scheme results in a spurious generation of nonthermal energy on shocks, while the entropy-conserving method evolves the energy adiabatically to machine precision. We also show that, in simulations of an isolated $L_\star$ galaxy, switching between the schemes results in $\approx 20-30\%$ changes of the total star formation rate and a significant difference in morphology, particularly near the galaxy center. We also outline and test a simple method that can be used in conjunction with the entropy-conserving scheme to model the injection of nonthermal energies on shocks. Finally, we discuss how the entropy-conserving scheme can be used to capture the kinetic energy dissipated by numerical viscosity into the subgrid turbulent energy implicitly, without explicit source terms that require calibration and can be rather uncertain. Our results indicate that the entropy-conserving scheme is the preferred choice for modeling nonthermal energy components, a conclusion that is equally relevant for Eulerian and moving-mesh fluid dynamics codes.


INTRODUCTION
Stellar feedback plays a critical role in shaping properties of simulated galaxies on all relevant scales, from the overall stellar mass and star formation rate to metallicity, galaxy morphology, and chemical abundances in the circumgalactic (CGM) and intergalactic medium (IGM; e.g., Governato et al. 2010;Brook et al. 2012;Agertz et al. 2013;Hopkins et al. 2013Hopkins et al. , 2014Stinson et al. 2013;Agertz & Kravtsov 2015, 2016, see also Naab &Ostriker 2017 andVogelsberger et al. 2020 for recent reviews). Modeling feedback processes in numerical simulations of galaxy formation has been the focus of extensive theoretical efforts.
Over the past decades a number of different methods have been developed to model stellar feedback, ranging from the direct injection of energy and momentum to the interstellar medium (ISM) when the physical scales of this injection can be resolved (e.g., Hopkins et al. 2011Hopkins et al. , 2012Hopkins et al. , 2018Agertz et al. 2013;Marinacci et al. 2019;Jeffreson et al. 2021;Smith et al. 2021), to effective models of the ISM and outflow driving when the multiphase ISM structure cannot be sufficiently resolved (e.g., Yepes et al. 1997;Springel & Hernquist 2003;Braun & Schmidt 2012;Vogelsberger et al. 2013;Springel et al. 2018).
One important class of feedback implementations is based on explicit modeling of a nonthermal energy component sourced by young stars in addition to gas thermal energy (e.g., Springel 2000;Agertz et al. 2013;Teyssier et al. 2013;Benincasa et al. 2016). Common examples of such nonthermal components are cosmic rays (e.g., Jubelgas et al. 2008;Uhlig et al. 2012;Booth et al. 2013;Pakmor et al. 2016;Ruszkowski et al. 2017;Wiener et al. 2017;Chan et al. 2019;Buck et al. 2020) and small-scale ISM turbulence (e.g., Braun et al. 2014;Schmidt et al. 2014;Semenov et al. 2016Semenov et al. , 2021Kretschmer & Teyssier 2020;Kretschmer et al. 2022). Such models differ by their sink, source, and extra transport terms (such as cosmic ray transport or turbulent diffusion). Without these terms, the behavior of the nonthermal energy is analogous to the thermal energy and consists of advection with gas density and the "P dV " work done during gas compression or expansion.
Implementation-wise, nonthermal energies can be modeled using the methods developed for modeling the thermal energy of the gas. Such methods were introduced to follow gas temperature in highly supersonic flows typical for galaxy formation simulations 1 where the modeling of thermal energy requires special attention. Indeed, when the flow is subsonic or only moderately supersonic, the thermal energy can be derived as the difference between the total and kinetic energy, e th = e tot − e kin , which are readily available in conservative methods that follow the gas mass, momentum, and total energy. However, when the flow is highly supersonic, the difference e tot − e kin becomes small and comparable to the truncation errors of e tot and e kin , requiring a more accurate method for modeling e th .
Two such methods for modeling thermal energy have been introduced in the early days of hydrodynamic galaxy formation simulations. In the "dual-energy" formalism proposed by Bryan et al. (1995), thermal energy is modeled as a separate variable in addition to total energy, with the P dV source term being computed explicitly. The second method, proposed by Ryu et al. (1993), models the gas entropy as a separate variable, which then can be used to compute temperature or thermal energy. The key advantage of the latter method is that gas entropy obeys a conservative equation, at least outside of shocks and in regions of the flow where dissipative processes are negligible. The entropy equation can thus be solved accurately without explicit source terms that can introduce significant errors.
Both these methods involve several parameters including the criteria that are used to decide whether e th 1 For example, the motion of the cold ISM due to galactic rotation and translational velocity (cs ∼ 1-3 km s −1 , v ∼ 100-300 km s −1 ), accretion of cold gaseous streams at high redshifts, and other motions of T ∼ 10 4 K gas in the IGM (cs ∼ 10 km s −1 , v ∼ few 100 km s −1 ) are flows with extremely high Mach numbers, v/cs > 30.
should be taken from the explicitly modeled variable or synchronized with the total energy as e th = e tot − e kin (Ryu et al. 1993;Bryan et al. 1995;Springel 2010;Teyssier 2015). In most idealized problems, both methods and different synchronization criteria work comparably well, and their differences tend do be small even in more realistic problems (e.g., Costa et al. 2020). However, in some cases, these choices can have significant effects. For example, Villasenor et al. (2021) showed that the formation of shock-heated gaseous halos in the IGM can be strongly delayed by the choice of the criteria used to synchronize thermal and total energies. The choice of the method for modeling of nonthermal components may have an even more significant impact. Indeed, the special treatment of thermal energy is required only in the extreme regime of highly supersonic flows, which are typically limited to a small fraction of the problem volume and where thermal pressure is by definition negligible. In contrast, nonthermal energies must be followed separately from the total energy in the entire simulation domain. Therefore, the choice of the method can impact results via the dynamical effect of nonthermal pressure and via its effect on other physical processes-such as star formation, radiative cooling, chemistry, etc.-which are sometimes coupled to the nonthermal energy.
In this paper, we explore the impact of different choices in modeling nonthermal energy on the evolution of gas using a series of idealized tests and more realistic simulations of an isolated L galaxy. In particular, we compare the two methods in shock tube and hydrodynamic Zel'dovich pancake tests and show that the energy-based method results in the spurious generation of nonthermal energy on shocks while the entropyconserving method can enforce the intended adiabatic behavior. A similar behavior was found by Kudoh & Hanawa (2016) and Gupta et al. (2021) in the context of numerical modeling cosmic rays as a nonthermal energy component. Here we further show that the entropy errors increase significantly for nonthermal components with larger adiabatic indices (such as subgrid turbulence with γ = 5/3 compared to cosmic rays with γ = 4/3), and for strongly compressive shocks. The latter becomes particularly problematic for modeling nonthermal components in galaxy formation simulations where highly compressive shocks are ubiquitous due to efficient radiative cooling in the dense ISM.
The paper is organized as follows. In Section 2, we review the methods of Bryan et al. (1995) and Ryu et al. (1993) and provide additional details about our implementation of these methods in the ART galaxy formation code. We also outline a simple method that can be used to implement the injection of nonthermal energy on shocks (Section 2.3). In Section 3, we present the results of shock tube and hydrodynamic Zel'dovich pancake tests. Then, in Section 4, we further compare the methods in simulations of an isolated L galaxy and show that the choice of the method can significantly impact galaxy properties. In Section 5, we discuss our results and propose a novel method for modeling unresolved turbulent energy by capturing numerical dissipation with the entropy-conserving scheme. We summarize our conclusions in Section 6.

IMPLEMENTATION OF THERMAL AND NONTHERMAL ENERGIES
To explore the impact of specific choices made in modeling nonthermal energy components, we use idealized tests and more realistic simulations of an isolated galaxy performed with the adaptive mesh refinement (AMR) hydrodynamics and N -body code ART (Kravtsov 1999;Kravtsov et al. 2002;Rudd et al. 2008;Gnedin & Kravtsov 2011). The hydrodynamic fluxes in the ART code are handled by a second-order Godunovtype method (Colella & Glaz 1985) with a piecewise linear reconstruction of states at the cell interfaces (van Leer 1979) and a monotonized central slope limiter based on Colella (1985).
Here we summarize our implementation of different methods for solving the set of advection equations with P dV work source terms for the thermal and nonthermal energy components. The key assumption is that the adiabatic evolution of these components can be modeled separately from any nonadiabatic processes, which can be added via source terms. Thus, during a hydrodynamics step and before the source terms are applied, the internal energies of the thermal and nonthermal components evolve adiabatically: Here and throughout the paper, e th and e nt denote energies per unit volume, and ∇ ≡ (∂/∂x, ∂/∂y, ∂/∂z) is the gradient operator. For clarity, we consider a single nonthermal energy component, e nt , but all of the results and conclusions can be trivially generalized for an arbitrary number of such energy components. 2 In our calculations we assume an ideal gas equation of state for each of the components: P i = (γ i − 1)e i with i ∈ [th, nt]. The above set of equations is coupled with gas dynamics via the total pressure P = P th + P nt . In addition, as thermal and nonthermal energies contribute to the gas total energy, e tot , a suitable method should be adopted to ensure that e tot = e kin +e th +e nt , where e kin = ρu u u 2 /2, when such synchronization is appropriate.
The choice of the synchronization method is dictated by the expected behavior of nonthermal energies across shocks. Indeed, in shocked regions, the difference e tot − e kin contains the adiabatic change in e th and e nt and the kinetic energy converted into thermal energy by the shock. In the absence of a subgrid model for nonthermal energy generation within shocks, the conservative choice is to assume that all nonadiabatic energy increase across shocks is thermalized, while the energies of the nonthermal components change adiabatically.
In this work, we consider the assumption that nonthermal components change adiabatically as the most "basic," conservative choice as this is the behavior implied by the underlying equations ( (1) and (2)) for smoothed flows. Any other, nonadiabatic behavior can be taken into account in the partitioning of e tot − e kin between e th and e nt (Section 2.3). In particular, nonthermal energies can be generated within shocks, e.g., via cosmic ray acceleration and turbulence driving. It is also worth noting that the shock structure and jump conditions are generally expected to be modified by kinetic non-MHD effects and propagation of cosmic rays across the shock (e.g., Voelk et al. 1984;Bret 2020;Haggerty & Caprioli 2020;Tsung et al. 2021). These modifications are uncertain and depend on the details and physical parameters of shocks, which generally are not possible to model self-consistently in cosmological simulations. Thus, the effects of such modifications need to be modeled phenomenologically as a part of the CR shock injection model. This was referred to as the closure problem of the shock solution in the presence of nonthermal energy components in a recent study by Gupta et al. (2021).
In the remainder of this section, we outline two methods implemented in the ART code for solving Equations (1)-(2) and our choices of synchronization criteria that determine when e th is reset using Equation (3). These methods were originally developed for modeling thermal energy in highly supersonic flows by Bryan et al. (1995) and Ryu et al. (1993), and in this paper we explore their behavior in the context of modeling nonthermal energies. In principle, different energy components can be treated using different methods, however, we opt for a consistent treatment of all components using the same method.

Energy-based scheme
In the "dual-energy" formulation proposed by Bryan et al. (1995), Equations (1)-(2) are solved in their original form: the energies e th and e nt are advected using conservative mass fluxes, and the P dV work is added as a source term, making the scheme nonconservative.
In the ART code, we advect e th and e nt as passive scalars: we interpolate specific energies, e i /ρ with i ∈ [th, nt], at the upwind side of the interface using the same reconstruction scheme as for all other hydrodynamic variables. The advection flux is then obtained by multiplying these reconstructed values by the average mass flux across that cell interface computed by the Riemann solver. This method ensures that the advection fluxes of e th and e nt are consistent with the mass flux, which enables the code to advect contact discontinuities to machine precision (see Appendix A and, in particular, the middle panel of Figure 10).
In Godunov-type methods, the P dV source terms can be computed relatively accurately and efficiently because the time-averaged velocity of gas at a given interface is computed during the solution of the Riemann problem across cell interfaces and can be used for an accurate estimate of ∇ ·u u u. In the ART code, ∇ ·u u u is computed by accumulating the fluxes of average velocities at cell interfaces and applying the Gauss-Ostrogradsky divergence theorem for each cell. Using this estimate of ∇ · u u u, the source term is applied to the energy components after the advection step.
This method of advancing e th and e nt is used in regions where the flow is highly supersonic and the values of e th from e tot can become smaller than the truncation error of the scheme, meaning highly inaccurate. In regions of modest Mach numbers, e th is synchronized with e tot using Equation (3). Specifically, the energies are synchronized in cells with with the value of η 1 = 10 −3 suggested by Bryan et al. (1995). This condition is equivalent to the Mach number Bryan et al. (1995) apply a second criterion based on the Mach numbers of adjacent cells: e tot,j − e kin,j max(e tot,j−1 , e tot,j , e tot,j+1 ) > η 2 , where j denotes the cell to which the criterion is applied and j ± 1 are its immediate neighbors. We use the value of η 2 = 0.1 as in the original scheme of Bryan et al. (1995), which corresponds to M 4 and makes this a significantly more conservative criterion than Equation (4).
Overall, these two criteria are used to select the method by which e th is evolved-i.e., either by updating its value from e tot (Equation (3)), or by solving Equation (1) independently of e tot . We choose to always follow e nt using Equation (2). For smaller η 1 and η 2 , e th is evolved synchronously with e tot in a larger fraction of the simulated volume, which in the limiting case of η 1 = η 2 = 0 corresponds to not using the dual-energy formulation at all.
Note that our implementation slightly differs from the original scheme of Bryan et al. (1995), in which the first criterion is only used to identify the cells in which thermal pressure is computed using explicitly advected e th , while only the second criterion is used to decide whether to synchronize e th and e tot . This is done to limit the dynamical effect of explicitly advected e th on the solution. However, given that the first criterion selects only the gas with e th e tot , the modeling of e th in such gas should have a negligible effect on the dynamics of the gas. We thus opt to use the first criterion to synchronize e th with e tot and always compute pressures from explicitly advected e th and e nt .
Finally, we also find that modeling of contact discontinuities in presence of multiple energy components with different adiabatic indices can be significantly improved if the total energy flux is computed consistently with the advection fluxes of e th and e nt described above. The total flux of thermal+nonthermal energy produced by the Riemann solver at a given interface is uP/(γ − 1), where γ = P/e + 1, with P = P th + P nt and e = e th + e nt , is the adiabatic index of gas which is explicitly followed by our Riemann solver (Colella & Glaz 1985). Away from shocks, this flux corresponds to the total advection flux of e th + e nt . However, the numerical value of this flux does not necessarily equal to the total flux computed by the advection scheme. We find that, even though this inconsistency is small, it still can produce significant artifacts in idealized tests like the shock tube test or advection of a pressure balance mode.
To make the total energy flux consistent with the advection scheme, we subtract uP/(γ − 1) from the Rie-mann solution for the flux of e tot and add the total flux of e th + e nt computed by the advection scheme. This correction should be done only outside shocked regions where the Riemann solution for pressure does not contain kinetic energy dissipated by shocks. For this reason, we apply this correction only on the interfaces where these two estimates of the total advection flux are within 3% of each other. By design, this correction results only in a small adjustment of the total energy flux, and the scheme remains explicitly energy conserving. We also find that this correction is necessary only in presence of multiple energy components with different adiabatic indices.
The necessity for the above correction is motivated by the design of the Colella & Glaz (1985) Riemann solver as it operates with the total pressure and effective adiabatic indices on the left and right sides of cell interfaces. Because of this, the total flux of thermal and nonthermal energy produced by the Riemann solver will generally be different from the same flux computed by the advection scheme. An alternative solution is to develop a Riemann solver that would operate on individually reconstructed pressure components and follow their evolution across the interface. To this end, one can use analytic solutions for the evolution of the nonthermal component with or without injection on the shock (see Appendix C in Pfrommer et al. 2017) or even try including kinetic non-MHD effects and propagation of cosmic rays (e.g., Voelk et al. 1984;Haggerty & Caprioli 2020;Bret 2020;Tsung et al. 2021). Such a Riemann solver will provide the solution for fluxes of different energy components without requiring a separate advection scheme. Note, however, that to advance thermal and nonthermal energies one still will be required to compute the source terms, which, as we will show below, can generate significant errors in shocked regions. Ryu et al. (1993) proposed an alternative method, where a modified entropy, ρS i ≡ P i /ρ γi−1 , is followed as a separate variable instead of energy. The equations for ρS i can be derived by combining Equations (1)-(2) with the continuity equation,
These equations are in the conservative form and thus can be solved with a numerical scheme that conserves ρS th and ρS nt to machine precision. 4 Given that the gas density can also be conserved to machine precision, this means that we can get accurate estimates of the entropy of the thermal and nonthermal energy components, S th and S nt , using ρS and ρ and computing the entropy fluxes, Sρu u u, consistently with the mass fluxes, ρu u u.
Implementation-wise, we use the same advection scheme for entropy as for the energy-based scheme described above, with specific entropy being advected as a passive scalar. Specifically, to compute the advection flux of entropy across a given cell interface, we interpolate the value of S i = P i /ρ γi (which can be thought of as a proxy for entropy per unit mass) at the upwind side of the interface using the same reconstruction scheme as for all other variables and multiply it by the average mass flux from the solution of the Riemann problem (Springel 2010). 5 The criteria for synchronizing S i with the total energy of the gas, e tot , also differ qualitatively from those in the energy-based method. Indeed, thermal entropy can be generated in shocks, and Equation (7) therefore becomes invalid in regions with strong shocks, meaning that e tot must be used instead. Thus, apart from the threshold on the Mach number of the flow, additional criteria must be used to identify shocked regions. For example, Ryu et al. (1993) check if the flow is locally converging (∇ · u u u < 0) and the pressure jump across the cell is larger than a chosen threshold value.
An alternative method to identify shocked regions was proposed by Springel (2010, Section 3.5). This method is based on the idea that in Godunov-type methods the production of entropy is taken care of by the shock(s) present in the Riemann solution for the hydrodynamic fluxes across cell interfaces. For weak shocks, the dissipated energy is strongly subdominant to the adiabatic increase in the internal energy, and therefore one can identify the regions with strong production of entropy by defining a threshold in the Mach number of such "Riemann shocks." Springel (2010) suggests the value M R,crit = 1.1, so that e tot is used in the regions where M R > M R,crit on at least one of the cell's interfaces. Al-4 In fact, as was pointed out by Ryu et al. (1993), one can construct an arbitrary number of conservative quantities of the form P α ρ 1−γα with α ∈ (−∞, ∞). This can be easily checked by expanding the conservation equation for such a quantity and using Equations (1), (2), and (6). The entropy conservation form corresponds to α = 1. Another special case of α = 1/γ, leading to a conserved variable P 1/γ , was explored by Kudoh & Hanawa (2016). 5 We also tried using the "raw" value of S i from the upwind cell without reconstruction but found that using reconstructed values of S i better preserves sharp contact discontinuities in different pressure components.
though the relation between the shocks in the Riemann solutions and the physical shocks in the simulated flow is not direct, the total entropy generated by the real shock accumulates from the increments produced by the "Riemann shocks" on the interfaces resolving the real shock, which warrants the above criterion. We determine the value of M R at a given interface from the total pressure jumps present in the solution of the Riemann problem at this interface. Such a solution consists of the left (l) and right (r) states sandwiching the so-called "star" region ( ) and separated by either shock(s) and/or expansion fans (e.g., Toro 2009). The relevant pressure jump, j P , is defined as the maximum of P /P l and P /P r . The corresponding Mach number can then be calculated from the Rankine-Hugoniot jump conditions, with pre-and post-shock regions corresponding to the cell on the side with the largest pressure jump (denoted with a subscript "j") and the "star" region, respectively (see also Pfrommer et al. 2017): where γ eff ≡ (γ th P th,j + γ nt P nt,j )/(P th,j + P nt,j ) is the effective adiabatic index in the cell on the side with the largest pressure jump, and γ i ≡ P i /e i + 1 with i ∈ [j, ], and P and e being the total pressure and ther-mal+nonthermal energy in the corresponding regions. In our implementation, γ j is computed using the values of P and e from the corresponding cell, while the value of γ is taken directly from the solution of the Colella & Glaz (1985) Riemann solver, which solves for the variation in γ explicitly. In the case when γ j = γ ≡ γ, Equation (9) can be simplified as We use this simplified form when 2|γ j − γ |/(γ j + γ ) < 0.01, following Pfrommer et al. (2017). Having an estimate of the Mach number of shocks present in the Riemann solution at all cell interfaces, we define M R for a cell as the maximum M R value on that cell's interfaces.
In our experiments, we find that the M R,crit threshold alone cannot prevent spurious heating in regions with strong velocity gradients. Such gradients can be interpreted by the Riemann solver as a discontinuity producing a sufficiently strong shock with a Mach number that can satisfy the above criterion. To filter out such cases, we require that the velocity of the shock(s) in the Riemann solution must be sufficiently large compared to the local flow velocity. In addition, to avoid injection of entropy in poorly resolved highly supersonic flows, we also require that the right-hand side of Equation (3) is larger than a certain fraction of e tot , which is analogous to the criteria used by Ryu et al. (1993) and Bryan et al. (1995). Note that the latter criterion may not be required in moving-mesh codes (e.g., Springel 2010) and when hydrodynamic equations are solved in a locally comoving frame (e.g., Trac & Pen 2004), where a significant fraction of kinetic energy is absorbed by the bulk flow of gas. For static-mesh Eulerian codes, however, this additional criterion is required to avoid artifacts in strongly supersonic flows.
Based on extensive experiments we find that the following set of criteria works well to identify shock regions: where max(M R ) is the maximal Mach number of the "Riemann shocks" at the interfaces of a given cell, |s R,max | is the corresponding shock velocity, and |v cell | is the gas velocity in the cell. When all of the above criteria are satisfied, the thermal energy is synchronized with total energy according to Equation (3), meaning that thermal entropy is injected into the cell in the amount corresponding to the energy dissipated by the shock. Otherwise, the evolution of thermal energy is followed using the entropy conservation Equation (7). To enforce the adiabatic behavior of nonthermal energy components, their evolution is always followed by Equation (8). As was pointed out by Springel (2010), one significant disadvantage of such a scheme is that it forfeits strict conservation of energy. In structure and galaxy formation simulations, however, the total energy is usually not conserved anyway owing to radiative cooling, radiative or feedback heating, star formation, etc. All of these processes are uncertain and modeled approximately. At the same time, as we show below, in the idealized tests, the entropy-conserving scheme performs either comparably or better than the energy-based scheme.
If the strict conservation of energy is nevertheless desirable, one can disregard the solution for the thermal component and compute e th from e tot while also following entropies of nonthermal components (Kudoh & Hanawa 2016;Gupta et al. 2021). Note that solving the equation for thermal entropy simultaneously with e tot is still desirable, so that gas temperature and pressure can be evaluated in regions with highly supersonic flows. In practice, such a scheme can be achieved by relaxing criteria (11) and (12) and using only Equation (13), so that explicitly modeled thermal entropy is used in highly supersonic regions. Such a scheme will ensure the conservation of both total energy and nonthermal entropy.

Generation of nonthermal energies in shocks
Synchronizing e th with e tot according to Equation (3) implies that all kinetic energy dissipated by shocks is thermalized. However, the synchronization procedure can easily be modified to account for nonthermal energy generation on shocks, such as cosmic ray acceleration or driving of small-scale turbulence when these components are modeled as e nt . Indeed, if the adiabatic evolution of e th and e nt can be enforced during a hydrodynamic step, then at the end of the step the difference e diss ≡ e tot − e kin − e th − e nt in each cell will correspond to the total energy dissipated by shocks during the step. Therefore, to convert a fraction ζ of this energy into a nonthermal component, one only needs to add ζe diss to e nt and the remaining (1 − ζ)e diss to e th . For ζ = 0, this scheme is equivalent to using Equation (3).
As we demonstrate in Section 3.2, this scheme works remarkably well when the adiabatic index of the injected energy is the same as that of the thermal energy, and it produces reasonable results when the adiabatic indices are different (e.g., for cosmic rays with γ nt = 4/3). At the same time, this scheme is trivial to implement within the entropy-conserving scheme as it requires only a minor modification of the energy synchronization algorithm, while a number of previous implementations of energy injection on shocks required on-the-fly shockfinding algorithms. Note that if ζ depends on the properties of the shock, as is generally the case for cosmic rays, the shock-finding algorithm is still needed to measure these properties for physical shocks. For example, one may try to estimate the local properties of the shock by considering several adjacent cells or using more complex algorithms (e.g., Ryu et al. 2003;Pfrommer et al. 2006;Skillman et al. 2008;Schaal & Springel 2015). However, the energy partitioning scheme can still be used to inject the corresponding nonthermal energy.
This algorithm is qualitatively similar to the method outlined by Gupta et al. (2021), where the total pressure is partitioned between the thermal and nonthermal components to maintain prescribed fractions in the regions identified as shocks. The main difference of the method described above is that only the kinetic energy locally dissipated by shocks is partitioned, while adiabatic compression is modeled separately by the entropyconserving scheme. This can be advantageous for modeling weak shocks, in which adiabatic compression of the pre-shock thermal and nonthermal components contributes significantly to the post-shock pressure. In principle, adiabatic compression can also be accounted for in the partitioning of the total pressure, however, the fractions of different pressure components depend on the shock Mach number and therefore it requires a shockfinding algorithm, while with the entropy-conserving method, the adiabatic behavior is built into the scheme.
3. IDEALIZED TESTS 3.1. Shock tube test with adiabatic nonthermal energy Figure 1 shows the results of a shock tube test with an additional nonthermal fluid representing cosmic rays (CRs) with γ nt = 4/3 (Pfrommer et al. 2006. The setup is analogous to the classic Sod (1978) test and has the same structure of the solution except that pressure consists of two components with different adiabatic indices, meaning that the effective γ can vary from region to region. As the figure shows, both energy-based and entropy-conserving methods reproduce the analytic solution for the density, velocity, and total pressure, which demonstrates that the fluxes of the mass, momentum, and total energy are computed correctly in both schemes.
The difference between the methods appears in the solution for the CR entropy and pressure as highlighted in the inset panels in the bottom row. In the energybased scheme, CR entropy is spuriously produced at the shock, leading to a ∼20% excess of the post-shock CR pressure. The post-shock thermal pressure becomes underestimated by the same (absolute) amount because in the post-shock region e th is computed from e tot and e cr using Equation (3). In contrast, the entropy-conserving scheme does not suffer from such errors: CR entropy is exactly conserved across the shock, resulting in correct CR and thermal pressures in the post-shock region.
This type of errors was also reported by Kudoh & Hanawa (2016) and Gupta et al. (2021), who also showed that switching to a conservative equation for nonthermal components (cosmic rays in their studies) helps to avoid these errors. All these results indicate that the entropy errors are related to the treatment of the source term in the nonconservative equation of the energybased scheme.
Indeed, the spurious production of nonthermal entropy in the energy-based scheme stems from the fact that inviscid Equations (1) and (2) are formally invalid at shocks because of divergent gradients. In a numerical solution, however, shocks are smeared over several cells via numerical diffusivity and all gradients are finite. As the gas flows through such shocked regions, the P dV source term computed from the local pressure and velocity gradients mixes the adiabatic compression with a fraction of the energy dissipated by the shock, leading to a nonadiabatic behavior. In other words, solving Equa- Dependence of the nonthermal entropy error produced by the energy-based scheme in the post-shock region on the compression ratio of the shock, rs = ρ1/ρ0, for different values of the adiabatic index of the nonthermal component: γnt = 4/3 (blue squares) and γnt = 5/3 (green circles). The entropy error increases monotonically with rs until ∆Snt/Snt becomes ∼ 1, at which point the error diverges. For a larger γnt, this divergence occurs at a smaller rs. For this test, the nonthermal energy is dynamically decoupled from the hydrodynamics by making its initial value arbitrary small and removing its contribution from the total pressure. The color intensity corresponds to different ways of varying rs. Lighter colors correspond to variations in γ th , keeping the shock Mach number high so that rs ≈ (γ th + 1)/(γ th − 1), with γ th spanning a range from γ th = 5 (leftmost points) to γ th = 1.16 and 1.36 (rightmost points for γnt = 4/3 and γnt = 5/3, respectively). Darker colors show the variation in the shock Mach number at fixed γ th = 5/3 (from left to right, M ≈ 1.3, 1.5, 1.9, 3.0, 9.0), with the vertical gray line showing the high-Mach-number limit of rs = 4.
tions (1) and (2) with explicit source terms results in injection of both thermal and nonthermal energy in the shocked regions in addition to the adiabatic increase due to gas compression. In the absence of nonthermal energies, this effect does not introduce severe errors because, in the post-shock region, the solution for the advected e th is usually disregarded and e th is reset from the difference e tot − e kin . In contrast, if nonthermal energies are present, only the sum e th + e nt can be set to e tot − e kin while the individual values of e th and e nt become affected by the entropy error.
The ∼ 20% error in the above test may seem small and unlikely to be important in real applications, e.g., galaxy formation simulations, due to much more uncertain source terms such as those involved in star formation and feedback modeling. However, we find that the error strongly increases for more compressible shocks and larger adiabatic indices of the nonthermal component. Indeed, as the entropy error originates from the  Results of the entropy-conserving scheme in the regime where the energy-based scheme is highly unstable. The initial conditions are the same as in Figure 1, but the adiabatic indices of the thermal and nonthermal components are set to γ th = 1.1 and γnt = 5/3 to make the shock highly compressible, with a compression ratio of rs ≈ 13.2. Such rs is significantly higher than the value at which the entropy error of the energy-based scheme diverges (rs ∼ 6.5 for γnt = 5/3, see Figure 2). The entropy-conserving scheme still performs well. P dV source term (P nt ∇ · u u u = P nt ∂u/∂x for a onedimensional problem), its magnitude relative to e nt increases with γ nt because P nt = (γ nt − 1)e nt . The magnitude of the P dV term also increases with the compression ratio of the shock (r s = ρ 1 /ρ 0 , where "0" and "1" denote pre-and post-shock regions, respectively) because the magnitude of the velocity gradient increases with r s , and the errors are amplified by stronger compression. Figure 2 shows the nonthermal entropy error in the post-shock region as a function of the compression ratio r s for γ nt = 4/3 (blue points) and γ nt = 5/3 (green points). To exclude the dynamical effect of the error on the solution we set the initial value of e nt to an arbitrary small value (< 10 −11 of e th in the post-shock region) and rerun the shock tube test with varying shock Mach number (i.e., varying the initial pressure jump) and adiabatic indices γ th and γ nt . Dark colors show the variation in the Mach number at fiducial γ th = 5/3, which leads to variation in the compression ratio with the limiting value of r s = (γ th + 1)/(γ th − 1) = 4 (shown with the vertical gray line). To explore variations of r s beyond this limiting value, we also rerun the tests with different γ th . Such a variation is relevant for practical applications, e.g., in galaxy formation simulations, because r s > 4 can be achieved due to strong radiative cooling, even if γ th is formally fixed at 5/3.
The figure shows that, while the entropy errors of the γ nt = 4/3 component on strong shocks with r s ≈ 4 are ∼ 20%, the errors of the γ nt = 5/3 component increase by a factor of 4. Thus, the spurious production of nonthermal entropy becomes more problematic for modeling components such as small-scale ISM turbulence with γ = 5/3 (e.g., Robertson & Goldreich 2012;Schmidt et al. 2014). Moreover, the error quickly increases as the shock becomes more compressive, and this increase saturates only when the post-shock e nt becomes comparable to the total energy dissipated by the shock, independent of the pre-shock value of e nt . As the figure shows, for larger γ nt , this divergence of errors occurs at smaller r s : r s ∼ 14 and 6.5 for γ nt = 4/3 and 5/3, respectively. As we will show in Section 4, such errors can significantly change the behavior of nonthermal components in galaxy formation simulations.
In the entropy-conserving scheme, the nonthermal entropy conservation equation is followed in the entire simulation domain. Figure 3 demonstrates that this scheme can produce an adequate solution even for strongly compressive shocks for which the energy-based scheme is highly unstable.

Shock tube test with sourcing of nonthermal energy by the shock
The shock tube test in the previous section assumes that all energy dissipated by the shock is thermalized, so that nonthermal energy behaves adiabatically. In this section, we switch to the case where nonthermal energy is also sourced by the shock as described in Section 2.3. This method requires an accurate modeling of the adiabatic compression of both e th and e nt . The shock-generated entropy errors discussed in the previ- . Shock tube test with 20% of dissipated energy injected in the nonthermal component using the entropyconserving scheme and the energy partitioning method outlined in Section 2.3. The initial conditions are the same as in Figure 1. The left and right columns show the results for different adiabatic indices of the nonthermal component: γnt = 4/3 and 5/3, respectively. The nonthermal energy injected at the shock is tracked separately and shown with green color in the bottom two panels. The analytic solution for all quantities is computed following Pfrommer et al. (2017, Appendix C2) and is shown with the thick lines. When γnt differs from γ th = 5/3, the injected energy in the post-shock region is underestimated by ∼ 20%, while for γnt = γ th = 5/3 the result is in perfect agreement with the analytical solution.
ous section prevent the application of this method with the energy-based scheme; therefore, we only present the results for the entropy-conserving scheme. Figure 4 shows the results of the test where 20% of the energy dissipated by the shock is converted into the nonthermal component, which is tracked separately and is shown with green lines. Different columns show the results for different values of the adiabatic index of the nonthermal component, namely, γ nt = 4/3 and 5/3. As this figure shows, for γ nt = 4/3 the method underproduces nonthermal pressure and entropy in the postshock region by a small amount, while for γ nt = 5/3 the method works remarkably well.
Qualitatively, this difference can be attributed to the effect that the injection of e nt changes the compressibility of the shock and to the fact that a shock is resolved by several cells, each dissipating shock energy. When γ nt is smaller (larger) than γ th = 5/3, converting a fraction of dissipated energy to e nt results in a lower (higher) effective adiabatic index of gas in the post-shock region, making the shock more (less) compressible. The total change in the injected e nt accumulates across the region over which the shock is numerically smeared and thus the effective adiabatic index also changes cell-bycell gradually changing gas compressibility. As the left column in Figure 4 shows, the total change in e nt does not add up correctly and results in smaller than expected injected e nt for e nt = 4/3 < γ th . We also checked that γ nt > γ th leads to an overestimate of the injected e nt .
In contrast, when γ nt = γ th , the effective adiabatic index stays constant and therefore the injection of γ nt does not affect the compressibility of the gas. As the right column in the figure shows, our injection scheme works well in this case.

Zel'dovich pancake test
The formation of a Zel'dovich pancake is a stringent test for an implementation of thermal and nonthermal energies because it involves both a purely adiabatic stage with extremely supersonic gas motion and the formation of strong shocks after the crossing time.
The solution for the density, ρ(x), and velocity, u(x), of a sine perturbation evolving from an initial redshift z i to some later redshift z can be expressed in a parametric form (Zel'dovich 1970): where q is the Lagrangian coordinate, ρ 0 is the average density, k = 2π/λ is the wavenumber of the initial perturbation, and the amplitude of the wave is described by the redshift of wave crossing, z c . . Hydrodynamic Zel'dovich pancake test with 70% of internal energy in the nonthermal form with γnt = γ th = 5/3. The results are shown before and after the wave crossing (left and right sets of panels, respectively) and for the energy-based and entropy-conserving methods (left and right columns in each set). The connected points show simulation results, while the thick lines show the analytic solution: orange lines show gas density, velocity, and total temperature, while red and blue lines show thermal and nonthermal components. Note that the analytic solution is valid only outside the shocked regions. Thin red lines in the top panels show the grid refinement level, with the 0-th level corresponding to 32 cells per box size. The bottom panels in the left set show the relative entropy error, while in the right set, they show the absolute value of entropy normalized to the initial total entropy. The energy-based scheme conserves nonthermal entropy relatively well during the adiabatic stage (with a relative error of < 1%). After the wave crossing, however, the nonthermal entropy suffers from strong numerical heating at the shocks. The entropy-conserving scheme, in contrast, conserves nonthermal entropy throughout the duration of the run, with the nonthermal "temperature" increasing only due to adiabatic compression.
For ease of comparisons, we choose parameters similar to Bryan et al. (1995), Trac & Pen (2004), and Springel (2010): λ = 64h −1 Mpc and z c = 1, and initialize the test at z i = 100. We add a nonthermal energy component with γ nt = 5/3 and initialize e th and e nt so that their initial entropies are constant with values corresponding to the average temperatures of T i = (µm p /k)S i ρ γ−1 0 = 30 K and 70 K for thermal and nonthermal components, respectively (where µm p is the average particle mass in units of the proton mass).
Both entropy components are expected to be conserved until the formation of the shocks at the wave crossing, and the solution for the temperatures is therefore (17) Figure 5 shows the results of the Zel'dovich test before (z = 2) and after (z = 0) the formation of shocks. Both energy-based and entropy-conserving schemes produce good results for the gas density and velocity: during the adiabatic stage, e th and e nt are negligible and thus their treatment does not affect the solution, while the agreement after the wave crossing indicates that both methods capture the formation of physical shocks at z = z c = 1.
In contrast, the solutions for the thermal and nonthermal components are quite different. As the lower left panel shows, before the crossing redshift, the energybased scheme conserves nonthermal entropy to subpercent level, while the thermal entropy suffers from strong numerical heating at the very center of the wave. The entropy-conserving scheme, in contrast, ensures entropy conservation for both components until the shock forms. After shock formation, both schemes produce similar results for the thermal entropy, while the solution for the nonthermal entropy in the energy-based scheme suffers from spurious heating at the shocks, resulting in an increase in the nonthermal entropy by a factor of ∼2 up to ∼50 at the very center of the caustic. The error has the same origin as the errors discussed in Section 3.1 and is completely absent in the entropy-conserving scheme.
The differences in thermal energy evolution before crossing illustrate the effects of the criteria selected to synchronize e th and e nt with e tot − e kin . Large value of thermal energy in the center in the energy-based scheme occurs because thermal energy is reset there from e tot , while e th is still smaller than e tot by orders of magnitude. The truncation error of the scheme (typically ∼ 10 −3 -10 −2 e tot ) thus propagates into e th as it is reset via e th = e tot − e kin − e nt , leading to orders of magnitude increase in the temperature at the wave center. This resetting in the center occurs because the criterion (e tot − e kin )/e tot > η can be satisfied near the center of the wave, where u ∼ 0 even though the flow is highly supersonic in the adjacent cells. Accounting for neighboring cells in the criterion (Equation (5)) or setting the threshold η to a larger value does not prevent this issue, but only delays its onset because the error in e tot − e kin accumulates with time.
The entropy-conserving scheme adopts additional criteria based on the Mach number and velocity of the shocks in the Riemann solution (Equations (11) and (12)) that filter out such cases, ensuring entropy conservation until the wave crossing. In addition, given that during the adiabatic stage the conservation of entropy is enforced instead of e tot , the errors in e tot − e kin do not accumulate.

ISOLATED GALAXY SIMULATIONS
To test different methods for modeling nonthermal energies in a more realistic environment, we use simulations of an isolated ∼L galaxy. Specifically, for the tests presented below, we use a snapshot from our fiducial simulation explored in Semenov et al. (2017Semenov et al. ( , 2018Semenov et al. ( , 2019 as the initial conditions. Below, we briefly de-scribe the aspects of this simulation that are most relevant to the current study and refer the reader to our previous papers for more details. One of the key ingredients of these simulations is the dynamic model for unresolved turbulence, which is implemented as a nonthermal energy. Our implementation (Semenov et al. 2016) is based on the "shear-improved" model of Schmidt et al. (2014). In this model, the turbulent energy on unresolved scales, e turb , is modeled as a nonthermal energy component with an adiabatic index of γ = 5/3, with source and sink terms that describe the turbulent cascade from the resolved velocity fluctuations and dissipation of turbulence on the local cell-crossing time. Apart from these source and sink terms, the subgrid turbulent energy is equivalent to e nt as introduced and tested in Sections 2 and 3.
The initial conditions for our isolated galaxy simulations are taken from the AGORA code comparison project (Kim et al. 2014). The AMR grid cells are adaptively refined when the gas mass in a cell exceeds ∼ 8300 M until the minimal cell size of ∆ = 40 pc is reached. Radiative cooling and heating of gas are modeled using the Gnedin & Hollon (2012) model, assuming constant metallicity at the solar value and a constant radiation background with the H 2 photodissociation rate in the Lyman-Werner bands of 10 −10 s −1 (Stecher & Williams 1967). We also adopt a model for dense gas shielding, which is calibrated against radiative transfer simulations of the ISM (the "L1a" model from Safranek-Shrader et al. 2017).
Local star formation is parametrized via the star formation efficiency per freefall time,ρ = ff ρ/t ff , with a fixed value of ff = 1%. Star-forming gas is defined by using a threshold either in density (Section 4.1) or in the local virial parameter (Section 4.2), where the latter is computed from the local turbulent velocity, σ turb = 2e turb /ρ, by using the definition from Bertoldi & McKee (1992): ≈ 9.35 (σ tot /10 km s −1 ) 2 (n/100 cm −3 )(∆/40 pc) 2 . (18) Stellar feedback is implemented by injecting radial momentum and thermal energy in the amounts calibrated against high-resolution simulations of SN remnants evolution in a nonuniform medium (Martizzi et al. 2015), with the momentum boosted by a factor of 5 to account for the effects of SN clustering (e.g., Gentry et al. 2017Gentry et al. , 2019 and cosmic rays (Diesing & Caprioli 2018). To approximate the effects of pre-SN feedback, we start the injection of momentum from the moment of star particle formation and continue the injection at a constant rate for 40 Myr. The injection rate is set by 5 kpc 10 −2 10 −1 10 0 10 1 10 2 n (cm −3 ) 10 −1 10 0 10 1 σ turb (km s −1 ) −10 −5 0 5 10 10 6 × (∆S turb /S turb,0 ) Figure 6.
Gas density (left), velocity dispersion of subgrid turbulence (center) and the error of the subgrid turbulence entropy (right) in an isolated disk galaxy simulation after ∼ 300 Myr of evolution. The subgrid turbulence is modeled using the entropy-conserved scheme. In this test, subgrid turbulence is initialized such that its entropy is initially constant, and then the galaxy simulation is evolved with all the source and sink terms of the subgrid turbulence model turned off, except for the P dV term. The right panel illustrates the long-term conservation of the entropy of subgrid turbulence, which remains constant to machine precision (∼ 10 −6 for a single-precision floating point variable) after ∼ 300 Myr of evolution. the total number of SNe occurring for a given star particle, which we compute assuming the Chabrier (2003) IMF. We also account for stellar mass loss by using the Leitner & Kravtsov (2011) model.

Adiabatic subgrid turbulence energy
We start by demonstrating that, in a realistic ISM with a wide variation in density and temperature, our implementation of the entropy-conserving scheme ensures a long-term conservation of nonthermal entropy. To this end, we initialize the subgrid turbulent energy such that its entropy is constant and turn off all its source and sink terms. As such a model for turbulence is not realistic, we use a star formation threshold not in α vir but in density, defining all gas with n > n sf = 100 cm −3 as star-forming. We select this threshold value because, for our simulated galaxy, it results in roughly the same amount of star-forming gas as the fiducial α vir threshold (Semenov et al. 2018). For this test, we also modified the original star formation and refinement algorithms such that the nonthermal entropy of gas in a cell is conserved when a star particle is formed or the cell is refined or de-refined. Figure 6 shows the maps of gas density, subgrid turbulent velocity, and entropy error after 300 Myr of evolution. In the absence of source and sink terms, the nonthermal entropy maintains its initial value to machine precision, ∼ 10 −6 for the single-precision floating point used in this test.
We do not show the results for the energy-based scheme because there the entropy errors discussed in Section 3.1 quickly accumulate resulting in a steady in-crease in the turbulent pressure that eventually stabilizes the disk. In a simulation with a more realistic processes, such runaway heating does not occur because the cooling term in the equation for the nonthermal energy leads to efficient dissipation of erroneously generated turbulent energy. Therefore, in the next subsection, we compare the schemes with all the sink and source terms activated.

Full subgrid turbulence model
As a more realistic test, we rerun our galaxy simulation including all source and sink terms in the subgrid turbulence model. We also use a star formation threshold based on the subgrid virial parameter, α vir < α vir,sf = 10, which was our fiducial choice in Semenov et al. (2017Semenov et al. ( , 2018Semenov et al. ( , 2019. Figure 7 compares the maps of gas density, subgrid turbulent velocity, σ turb , and SFR surface density in simulations with the energy-based and entropyconserving schemes. Although the overall distributions of n and σ turb are quite similar, the figure reveals interesting differences. With the entropy-conserving scheme, the overall density structure becomes more flocculent, especially near the disk center. The difference in theΣ map is even more striking: while with the energy-based scheme theΣ distribution is rather clumpy and has a strong deficiency in the central ∼ 1 kpc, in the simulation with the entropy-conserving scheme,Σ is smoother and is not suppressed near the center. To make the comparison more quantitative and explain these differences, Figure 8 shows joint distributions of gas in the n-σ tot plane colored according to the galac-energy-based entropy-conserving 10 −2 10 −1 10 0 10 1 10 2 n (cm −3 ) 10 0 10 1 10 2 σ turb (km s −1 ) 10 −2 10 −1 10 0 Σ (M yr −1 kpc −2 ) Figure 7. Results of the isolated galaxy simulation with the full subgrid turbulence model and star formation prescription coupled with the predicted turbulent velocities. The columns show the mid-plane slices of density, n, and subgrid turbulent velocity, σ turb , and the surface density of SFR,Σ , computed using particles younger than 30 Myr. Different rows show simulations where thermal energy and subgrid turbulence are modeled by using either the energy-based or entropy-conserving scheme. To guide the eye, black contours in the n and σ turb maps show n = 10 cm −3 density isocontours. TheΣ maps show the most apparent difference between the two schemes: the SFR in the run with the energy-based scheme is more clumpy and strongly suppressed near the galaxy center.
tocentric radius of the cells and with the thick dotted line showing the star formation threshold.
In the warm and diffuse part of the ISM with n < 10 cm −3 , the overall shapes of these distributions are similar. In most of such gas, σ tot is dominated by the sound speed, c s ∼ 4-20 km s −1 , which sets the lower envelope of the distribution at n < 10 cm −3 , reminiscent of the shape of the n-T diagram at such densities. We checked that the distributions of σ turb alone are also quite similar at n < 10 cm −3 in both schemes.
In contrast, there are two interesting differences in the cold and dense ISM, n > 10 cm −3 , where subgrid turbulence dominates over thermal energy. First, the turbulent velocities in the main part of the distribution decrease by a factor of ∼1.5 when the entropyconserving scheme is used. Second, with the energybased scheme, the gas near the disk center manifests itself as a prominent feature with n > 100 cm −3 and σ turb ∼ 30-50 km s −1 , while with the entropyconserving scheme, this feature disappears, and σ turb in the central region becomes a continuation of the rest of the distribution.
The source of these differences is the production of e turb in radiative shocks via the mechanism described in Section 3.1. The average distribution of cold gas in the n-σ tot plane is set by the competition between turbulence production and its decay on the local cell-crossing time. While both simulations include physical production terms-adiabatic compression and sourcing by fluctuating velocities-the numerical heating discussed in Section 3.1 also contributes to the production of e turb when the energy-based scheme is used. The fact that in most of the cold ISM the increase in σ turb is only moderate indicates that this heating mechanism does not dominate over the physical production terms. The only exception is the disk center, where switching to the entropy-conserving scheme leads to a decrease in σ turb by a factor of ∼3, which corresponds to a decrease in e turb by an order of magnitude. . Joint distribution of gas density and total subgrid velocity dispersion, which includes both subgrid turbulence and thermal sound speed, σtot = σ 2 turb + c 2 s . The contours show 68%, 95%, and 99% of the gas mass and the color indicates the average galactocentric radius of cells within a pixel. The thick dotted line shows the star formation threshold, α vir,sf = 10, while the thin lines show the values of αvir = 1 and 100 for reference. To reduce noise, the distributions are averaged over 21 snapshots spaced 10 Myr apart. Switching to the entropy-conserving scheme reduces σtot ≈ σ turb in dense gas (n > 10 cm −3 ) by a factor of ∼1.5 and completely removes the horizontal feature at σ turb > 30 km s −1 and n > 100 cm −3 that corresponds to the central region of the disk.
Even though in the context of our subgrid turbulence model the excess of e turb in the disk center is purely numerical, it may reflect a physical enhancement of turbulence production in that region due to strong shear. Such an enhancement may not be captured by the adopted turbulence production term in the model, but it does appear as an enhancement of numerical dissipation of kinetic energy that is partially converted to e turb (see Section 5.2 for further discussion). Therefore, here we solely present the difference between the two schemes of modeling nonthermal energies and leave the detailed investigation of the behavior in the disk center to a future study.
The decrease in σ turb in the cold and dense gas leads to significant changes in the SFR distribution. For example, as σ turb drops near the disk center, gas can cross the star formation threshold more easily in the run with the entropy-conserving scheme in striking contrast with the strongly suppressed SFR in the central ∼ 1 kpc in the run with the energy-based scheme (recall the right panels in Figure 7). At larger galactocentric radii, distribution of star-forming regions also becomes different. As Figure 8 shows, when the entropy-conserving scheme is used, a larger fraction of gas from the main part of the distribution (i.e., excluding the feature corresponding to the disk center) reaches low virial parameters, α vir < α vir,sf = 10, and the typical densities of such star-forming gas become lower. Figure 9 quantifies the effects of the turbulent energy scheme on the SFR, total mass of star-forming gas, M sf , and its average density,n sf , at R > 1 kpc. These three quantities are related viȧ where ff = 1% and the average density of star-forming gas is defined from 1/t ff −1 sf ≡ 3π/(32Gµm pnsf ), assuming µ = 1 and ... sf denoting mass-weighted averages over star-forming regions.
As the figure shows, the effect of the scheme on these quantities is quite significant, at the level of ∼ 10%-30%. In the entropy-conserving scheme, n sf is smaller and M sf is larger, which is consistent with the change in the PDF in Figure 8. These changes partially cancel in the total SFR, as the larger mass of star-forming gas is partly offset by the lower density and longer freefall time of star-forming regions. The total SFR in the run with the entropy-conserving scheme thus increases only moderately.
The direction of these trends is consistent with the model for the origin of global gas depletion times, τ dep ≡ M g /Ṁ , and SFR from continuous and rapid gas cycling between actively star-forming and diffuse, nonstar-forming states in the ISM (Semenov et al. 2017). In this model, the global depletion time is the sum of the times that gas spends in each of these states, which can be written as τ dep = N c t nsf + τ ff / ff , where N c is the total number of cycles required for gas depletion, t nsf is the average residence time of gas in the non-star-forming state, and τ ff / ff is the average depletion time of the . Effect of thermal and nonthermal energy modeling on the star formation in the average disk, i.e., at R > 1 kpc. The panels show the total star formation rate, mass of star-forming gas, and its average density. When the entropy-conserving scheme is used, a lower σ turb in dense gas (see Figure 8) results in a larger amount of star-forming gas but also in a lower average density of such gas. These trends partially cancel out, leading to only a moderate increase in the SFR.
star-forming gas. Switching to the entropy-conserving scheme reduces the production of e turb , making it easier for gas to cross the α vir ∝ σ 2 turb /n threshold, which leads to shorter t nsf and longer τ ff . These trends counteract in the expression for τ dep , but they add up in the expression for the star-forming gas mass fraction: f sf = (τ ff / ff )/τ dep = ( ff N c t nsf /τ ff + 1) −1 . As a result, the increase inṀ = M g /τ dep is smaller than the increase in M sf = f sf M g as M g stays roughly constant over the considered time interval.

Implications for galaxy formation simulations
The results presented in the previous sections show that the choice of the method to model thermal and nonthermal energy components matters. In particular, the energy-based scheme can lead to a spurious generation of both thermal and nonthermal energy components. This was apparent in the Zel'dovich pancake collapse test, which showed significant spurious temperature increase at the center of the pancake right before crossing. Similarly, Villasenor et al. (2021) showed that the mean IGM temperature in their simulations is sensitive to how thermal energy of the gas is modeled and to the parameters of the dual-energy scheme (see their Section 3.4 and Figures 5 and 6).
This difference arose because shock heating of gas in collapsed regions of virialized halos could be reduced or suppressed by the condition for synchronizing thermal energy followed independently with the value computed from the total energy. The results of Villasenor et al. (2021) indicate that these issues can be corrected by the proper choice of thermal energy synchronization parameters, meaning that both energy-and entropy-based schemes can accurately and robustly model thermal energy. However, our results show that the choice of such scheme for any nonthermal component(s), such as subgrid turbulence or cosmic rays, does make a difference. This is because nonthermal components are necessarily evolved separately from the total energy in the entire simulation domain, while such treatment of the thermal energy alone is only needed in a limited volume corresponding to highly supersonic flows.
Many modern fluid dynamics codes employ different schemes to separately model the coherent bulk motion of the fluid and small-scale motions (Trac & Pen 2004;Springel 2010;Duffell & MacFadyen 2011Hopkins 2015;Duffell 2016). Such "moving-mesh" approaches alleviate the inaccuracy of the thermal energy calculation from the total energy by absorbing a significant fraction of kinetic energy into the motion of the mesh. However, our results are equally relevant for these approaches because they indicate that the entropy-conserving scheme is the preferred choice for modeling the energies of nonthermal components. Indeed, for the energy-based formulation, the errors in the nonthermal entropy originate from the nonadiabatic part of the explicit P dV source term inside numerically smeared shocks which occur both in Eulerian and moving-mesh fluid dynamics schemes.
How much the choice of the scheme matters depends on the problem at hand. In the relatively dense parts of the ISM, where cooling is efficient, the spurious generation of energy near shocks may be inconsequential as it is dissipated efficiently. On the other hand, as we showed in the previous section, this choice may result in nonnegligible changes in the star formation rate of a galaxy, the structure of the ISM, and in particular qualitatively different star formation in the central regions of the simulated galaxy.
The origin of the entropy error is closely related to the closure problem of the shock solution in the presence of nonthermal energy components (Gupta et al. 2021). The classical Rankine-Hugoniot jump conditions become insufficient and require additional relations to define the behavior of each of the extra nonthermal components. In the methods explored in our paper, these closure relations are effectively set by the choice of the method to evolve different components and the algorithm that synchronizes these components with the total energy (e.g., Equation (3) or Section 2.3).
The results of Gupta et al. (2021) and our findings indicate that the energy-based implementation implicitly defines a closure relation that deviates from an adiabatic behavior across shocks despite the fact that the underlying equation describes adiabatic evolution. More importantly, the behavior of nonthermal components becomes dependent on the properties of the shock and on the details of the numerical implementation. This is because entropy errors originate from the nonadiabatic excess of the P dV source term in the region over which the shock is numerically smeared. The structure of the smeared shock depends on both physical properties of the shock and details of the numerical method, and so will the entropy error in the post-shock region. This error thus can be avoided by either getting rid of the source term altogether, as we do in this study (see also Kudoh & Hanawa 2016;Gupta et al. 2021), or by imposing an explicit closure relation at shocks, as suggested by Gupta et al. (2021).
Although the entropy-conserving scheme is a wellposed solution to the above problem, Kudoh & Hanawa (2016) and Gupta et al. (2021) reported issues arising when advecting a "pressure-balance" mode, where gas density, velocity, and total pressure are all constant, while the fractions of thermal and nonthermal pressure components experience a jump. We have repeated such a test and found that our implementation of the entropyconserving scheme performs extremely well in this test: density, velocity, and total pressure are all preserved to machine precision (see Appendix A). Thus, the performance in the "pressure-balance" test likely depends on the details of scheme implementation.
Regardless of numerical issues, a strong argument for the entropy-based scheme for nonthermal energy components can be made in models that include the injection or dissipation of energy in shocks. For example, both turbulence and cosmic ray energy can be generated at shock fronts and, as we discussed in Section 2.3, such generation can be handled with ease in the entropyconserving scheme.
This approach can be especially advantageous for sourcing nonthermal energies on radiative shocks. The existing methods based on on-the-fly shock-finding algorithms estimate dissipated energy by measuring the pressure jump across the shock. This can be inaccurate when the pressure in the post-shock region is subject to significant radiative losses. This issue can be alleviated in the proposed method because the energy dissipated by shocks at all cell interfaces can be converted to nonthermal components before radiative losses occur. We leave a more detailed investigation of this issue and a comparison of different methods to future study.
Finally, the entropy-conserving scheme can provide another significant advantage for modeling subgrid turbulence because the resolved turbulent energy dissipated by numerical viscosity can be captured into the subgrid turbulent energy implicitly, as we outline in the next section. This allows one to avoid explicit modeling via source terms that require calibration and can be rather uncertain.

Sourcing subgrid turbulence via capturing kinetic energy dissipated by numerical viscosity
Explicit modeling of unresolved turbulence in galaxy formation simulations can significantly improve the modeling of processes that depend on the small-scale turbulent structure of the ISM, such as star formation, mixing and transport of metals, and the effect of gas clumping on cooling and chemical reaction rates. Such subgrid models of turbulence require a prescription for the energy transfer from the resolved flow to unresolved scales.
For example, in the Large Eddy Simulation (LES) approach adopted in Section 4.2 (see Garnier et al. 2009 andSagaut 2006 for reviews), such source terms are modeled using explicit closure relations, which can be calibrated against direct high-resolution simulations of turbulence (e.g., Schmidt & Federrath 2011). These closures are the major source of uncertainty in modeling subgrid turbulence, especially when the onset of the turbulent cascade is not sufficiently resolved or the local flow is significantly different from the direct simulations used to calibrate the closures.
On the other hand, the strict enforcement of entropy conservation opens up a new, more generic way for modeling the energy cascade to unresolved scales via tracking the local dissipation of the resolved kinetic energy. Indeed, as our results show, this scheme can be used to accurately follow adiabatic evolution of thermal and nonthermal energies during each hydrodynamics step. In regions where thermal energy can be computed reliably via the difference of the total and kinetic energies, e tot − e kin , the change in this energy during a single step can be compared to the adiabatically evolving sum e th + e turb modeled separately using the entropyconserving scheme. 6 The difference between the two estimates will then correspond to the total kinetic energy dissipated by numerical viscosity in each cell during the step.
In most galaxy simulations, this dissipated energy is converted directly into heat and then quickly radiates away. In the real ISM, however, the energy is transferred to turbulent motions on small scales and decays on a longer timescale, close to the local eddy-turnover time. The above method can be used to capture this energy and to convert it into heat on a physically motivated timescale, e.g., turbulent cell-crossing time with an order-of-unity pre-factor that can be calibrated in direct simulations of turbulence as in LES simulations.
Even though the dissipation of kinetic energy in this method is purely numerical and happens on the scale of a few computational cells, empirical evidence suggests that the rate of dissipation is correct. Indeed, turbulent box simulations carried out with different codes generally find turbulent spectra consistent with theoretical expectations, suggesting that the rate of the energy transfer on resolved scales is modeled correctly (e.g., Kritsuk et al. 2011). Given that in developed turbulence, this energy transfer rate is equal to the dissipation rate on the viscous scale, the rate of numerical dissipation must also be equal to the physical dissipation rate.
One additional issue that needs to be considered in such an approach is that not all kinetic energy dissipated by numerical viscosity should source the subgrid turbulence. For example, numerical viscosity can dissipate energy in laminar shearing flows, but no turbulence should be generated in this case. In the context of the LES approach, this motivates the usage of the "shear-improved" versions of the schemes, in which the bulk flow of the gas is factored out from the source terms of subgrid turbulence (Lévêque et al. 2007;Schmidt et al. 2014). Con-versely, a fraction of the energy dissipated by resolved shocks can drive turbulent motions in the post-shock region instead of being fully thermalized. These examples show that an approach described above should include a model for how the dissipated energy is apportioned between thermal and subgrid turbulence energy based on the local flow properties, such as the shear-improved approach and the scheme outlined in Section 2.3.
The idea of unresolved turbulence sourcing by numerical dissipation is a part of a more general assumption upon which the concept of the Implicit Large Eddy Simulations (ILES) is based (see, e.g., Grinstein et al. 2007 for a review): as long as a hydrodynamics solver satisfies the basic physical properties of the equations of hydrodynamics, such as conservation of mass, momentum, and energy, as well as causality and positivity, the effect of truncation errors at the resolution scale is equivalent to the net effect of unresolved scales. In this context, the approach outlined above captures the energy transfer to unresolved scales for the explicitly modeled subgrid turbulence without explicit and uncertain terms and closure relations. Instead, the energy transfer is implicitly handled by the hydrodynamic solver. The subgrid turbulent energy can then be used to model star formation, metal diffusion, and gas clumping factor in exactly the same way as in the LES simulations.

SUMMARY AND CONCLUSIONS
We have investigated the accuracy of energy-based and entropy-conservative schemes in modeling nonthermal energy components. These schemes were initially introduced for modeling thermal energies of highly supersonic flows where e th is strongly subdominant to the total and kinetic energies and therefore cannot be accurately computed as the difference between the two. The energy-based scheme numerically solves the explicit equation for the energy density evolution, which cannot be written in conservative form. The entropy-conserving scheme, on the other hand, uses a conservative equation for modified entropy and can thus conserve the entropy of an energy component to machine precision.
We have examined the performance of the schemes in following the subgrid turbulence and cosmic ray energy in the standard shock tube and Zel'dovich pancake tests and simulations of a realistic isolated L galaxy. Our results can be summarized as follows.
1. The energy-based scheme results in a spurious generation of nonthermal energy on shocks in the shock tube and Zel'dovich pancake tests, while the entropy-conserving method evolves the energy adiabatically to machine precision (see Figures 1, 3, and 5).
2. The magnitude of the nonthermal entropy error in the post-shock region increases with the adiabatic index of the nonthermal component and with the compression ratio of the shock (see Figure 2). The latter is particularly relevant for the cold ISM where shocks can be strongly compressible due to efficient radiative cooling.
3. In simulations of an isolated L galaxy with a turbulence-based star formation prescription, switching between the energy-based and entropyconserved schemes results in a qualitative change in morphology, in particular, a large qualitative difference in star formation in the center of the galaxy (Figure 7), and ≈ 20-30% change in the star formation rate away from the center (Figure 9).
4. We outline and test a simple, physical method for the injection of nonthermal energy on shocks. This method can be used in conjunction with the entropy-conserving scheme and can simplify the implementation of, e.g., injection of cosmic rays and driving of subgrid turbulence by shocks (see Section 2.3 and Figure 4).
Finally, we discuss how the entropy-conserving scheme can provide a straightforward way to capture the kinetic energy dissipated by numerical viscosity into the subgrid turbulent energy implicitly, without explicit source terms that require calibration and can be rather uncertain. It will be interesting to investigate the performance of such an approach using idealized simulations of turbulence and realistic galaxy formation simulations.
Although both methods have been shown to perform well in modeling the evolution of thermal energy, our results indicate that the entropy-conserving scheme is a preferred choice for modeling nonthermal energy components. This conclusion is equally relevant for Eulerian and moving-mesh fluid dynamics codes.

ACKNOWLEDGMENTS
We would like to thank Nick Gnedin, Peng Oh, Tsun Hin Navin Tsung, and Siddhartha Gupta for their insightful comments. We are also grateful to the anonymous referee whose comments helped to improve this paper. Support for V.S. was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51445.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. A.K. was supported by the National Science Foundation grants AST-1714658 and AST-1911111 and NASA ATP grant 80NSSC20K0512. This work was completed in part with resources provided by the University of Chicago Research Computing Center and by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. The analyses presented in this paper were greatly aided by the following free software packages: NumPy (van der Walt et al. 2011), SciPy (Jones et al. 2001(Jones et al. -2016, Matplotlib (Hunter 2007), and yt (Turk et al. 2011). We have also used the Astrophysics Data Service (ADS) and arXiv preprint repository extensively during this project and the writing of the paper.

A. ADVECTION OF A PRESSURE-BALANCE MODE
Kudoh & Hanawa (2016) and Gupta et al. (2021) compared the performance of nonconservative and conservative schemes for modeling cosmic rays and reported issues of the latter scheme in the "pressure-balance" test. In this test, density, velocity, and total (thermal plus nonthermal) pressure are constant throughout, but the fractions of pressure contributed by the thermal and nonthermal components experience a jump at x = 0. Given that pressure and velocity are constant, no force should be generated and the jump in the pressure fractions should simply be advected, while other quantities should remain constant. However, Kudoh & Hanawa (2016) and Gupta et al. (2021) found that their imple-mentation of the conservative scheme produces numerical artifacts in this test, while Gupta et al. (2021) also found that their nonconservative energy-based method performs better.
The results of this test with our implementation of the energy-based and entropy-conserving schemes are shown in Figure 10. Both our schemes perform extremely well in this test: gas density, velocity, and total pressure are all conserved to high precision, even though the jump in the pressure components becomes smeared over several cells. In particular, our entropy-conserving scheme can model advection of the pressure-balance mode to machine precision, in contrast to the results reported by Kudoh & Hanawa (2016) and Gupta et al. (2021).
Why does our entropy scheme perform in this test better than similar models from previous literature? Not  Figure 10. Results of the "pressure-balance" test at t = 0.25. Different columns correspond to different schemes, from left to right: the energy-based scheme described in Section 2, the energy-based scheme that does not synchronize with the solution for etot, and the entropy-based scheme. The initial conditions are the same as in Gupta et al. (2021): (ρ, u, P th , Pcr) = (1, 1, 0.1, 0.9) and (1, 1, 0.9, 0.1) for the left and right initial states, respectively, and the initial location of the discontinuity is marked by the vertical green dotted line. The adiabatic indices of thermal and cosmic ray components are 5/3 and 4/3, respectively. To avoid the overlap in the top panels, the values of ρ and u are shifted by 0.4 and 0.2, respectively. The shaded gray area in the bottom panels shows the range where the scale of the y-axis is linear; the horizontal gray lines indicate the relative error of 1%.
surprisingly, the results of this test strongly depend on the modeling of energy or entropy advection. Indeed, the P dV source term should stay exactly 0 as long as gas velocity stays constant, and therefore any spurious pressure gradients can originate only from the inconsistencies in the advection fluxes of different energy components. These inconsistencies can arise when different components are modeled using different schemes. For example, in the entropy-based methods discussed by Kudoh & Hanawa (2016) and Gupta et al. (2021), only the nonthermal component is modeled using an entropyconserving scheme, while the thermal component is computed from the solution for the total energy. In contrast, in our implementation of the entropy scheme, both entropy components are modeled consistently, using the same advection scheme, which leads to a significantly better performance in this test.
In the energy-based scheme, significant artifacts may also appear when the advection fluxes of e th and e nt are modeled inconsistently with the Riemann fluxes of mass and total energy. For example, we find that not applying the e tot flux correction described at the end of Section 2.1 produces artifacts in the gas density, velocity, and total pressure at a >10% level. In contrast, this correction reduces these errors to 10 −4 (see the left panel in Figure 10). This error is somewhat larger than in the entropy-conserving scheme because e th is reset from the difference e tot − e kin − e nt . To demonstrate this, the middle panel shows the results of the energybased scheme without this synchronization, which only uses advected e th . 7 In this case, the error has the same small magnitude as in the entropy-conserving scheme.
All in all, these results show that the performance of the schemes in the "pressure-balance" test depends on the implementation of the advection for different energy components. In contrast to the findings of previous studies, our implementation of the entropy-conserving scheme performs well in this test.