Well-balanced methods for computational astrophysics

We review well-balanced methods for the faithful approximation of solutions of systems of hyperbolic balance laws that are of interest to computational astrophysics. Well-balanced methods are specialized numerical techniques that guarantee the accurate resolution of non-trivial steady-state solutions, that balance laws prominently feature, and perturbations thereof. We discuss versatile frameworks and techniques for generic systems of balance laws for finite volume and finite difference methods. The principal emphasis of the presentation is on the algorithms and their implementation. Subsequently, we specialize in hydrodynamics’ Euler equations to exemplify the techniques and give an overview of the available well-balanced methods in the literature, including the classic hydrostatic equilibrium and steady adiabatic flows. The performance of the schemes is evaluated on a selection of test problems.


Introduction
Numerical methods for the approximate solution of balance laws play a central role in the simulation of many interesting and challenging phenomena in computational astrophysics. Balance laws take the generic form where u is the vector of conserved variables, f the flux tensor and s the vector of source terms, respectively. Examples of balance laws include (ideal) hydrodynamics or Euler equations, (ideal) magnetohydrodynamics, and radiation (magneto)hydrodynamics, as well as their relativistic counterparts. The origin of the source term on the right-hand side may be physical (e.g., chemical reactions, external forces, nonideal effects), geometric (e.g., curvilinear coordinates) or both (e.g., curved spacetime). The faithful modeling of complex astrophysical phenomena with balance laws is generally not feasible with (semi-)analytical methods alone. Hence, solutions can only be sought approximately by numerical means. Numerical methods for (hyperbolic) conservation laws, that is, the homogeneous Eq. (1) with s 0, are in a mature stage of development. We refer, for example, to the recent comprehensive review by Balsara (2017) and references therein. The approximation of balance laws is often not much more involved and can easily be done by supplementing a consistent discretization of the source term s. In this way, highly accurate solution approximations can be obtained efficiently by computational means.
However, there are particular regimes where conventional numerical methods encounter difficulties. 1 Balance laws often possess non-trivial steady-state solutions where the flux divergence exactly balances the source term. Numerical methods do not necessarily satisfy a discrete version of this subtle equilibrium balance. Consequently, steady states are not resolved exactly but are approximated with an error of the order of the method's truncation error. To simulate phenomena near steady states, the numerical resolution needs to be high enough such that the continuous pile-up of these truncation errors does not obscure the processes of interest during the simulation timeframe. Especially in multi-dimensional simulations, the required resolution may entail prohibitively high computational costs. The shortcoming near steady states was realized early on in the development of numerical methods for balance laws. This lead to the suggestion to exploit the steady-state solutions in the discrete representation of the approximate solution (see, e.g., Liu 1979;Glaz and Liu 1984;Glimm et al. 1984;van Leer 1984;Huang and Liu 1986;Roe 1987). The idea is to replace the common (piecewise) polynomial approximate solution representation with a (piecewise) steady one that fulfills the subtle balance Eq. (2) either exactly or approximately. Thereby, only deviations from steady state induce dynamics which is indeed highly desirable. As a matter of fact, this is a generalization of the fundamental property of numerical methods for the homogeneous equations that only deviations from a constant state trigger wave motion. As noted by van Leer (1984), the construction of such steady-state distributions is in general difficult as it requires the local solution of a boundary value problem for Eq. (2). The solvability is challenging and a solution can usually only be obtained numerically. The pioneering work by Mellema et al. (1991), Eulderink and Mellema (1995) suggests constructing an equilibrium subgrid model by locally approximating equivalent initial value problems numerically. For example, such equilibrium subgrid models were successfully constructed for hydrostatic equilibrium by Zingale et al. (2002) and even for general relativity by Kastaun (2006). This led to much improved numerical resolution near equilibrium states.
An additional design principle was introduced by Cargo and LeRoux (1994) to overcome the challenges near steady states. They constructed a scheme for the Euler equations with gravity source terms capable of preserving exactly a discrete form of hydrostatic equilibrium and termed the scheme as well-balanced (or, in French, ''un schéma équilibre''). A well-balanced numerical method satisfies a discrete form of the equilibrium balance Eq. (2) exactly, independent of the resolution. Therefore, these methods can accurately resolve solutions that are small perturbations of equilibrium data. Many such schemes have been developed since, especially for the shallow water equations with bottom topography used in environmental applications. In the context of the shallow water equations, the well-balanced property is also referred to as the exact C-Property put forward in a seminal paper by Bermudez and Vazquez (1994). We refer to the comprehensive reviews by Noelle et al. (2009), Xing and Shu (2014), Kurganov (2018), the textbook by Bouchut (2004) and the references therein for further information. An extensive review of well-balanced and related schemes for many applications can also be found in the textbook by Gosse (2013). Moreover, we refer to Amadori and Gosse (2015) for an extensive theoretical treatment and rigorous numerical analysis of well-balanced schemes on simple balance laws.
A popular framework for the construction of well-balanced numerical methods is rooted in the piecewise steady or subgrid equilibrium representation. The framework combines a piecewise steady reconstruction, consisting of an equilibrium subgrid model and a piecewise polynomial equilibrium-preserving reconstruction, and a well-balanced source term discretization. Many of the aforementioned schemes above have been constructed along these ingredients. In this text, we focus on this framework as it combines conceptual simplicity and versatility in that it applies to a wide range of numerical methods for balance laws ranging from finite volume to discontinuous Galerkin over finite difference methods. For the clarity and conciseness of the presentation, we concentrate on particular flavors of these numerical methods. In particular, we focus the presentation on higher-order Godunov-type finite volume methods and finite difference methods with flux splitting. Moreover, the emphasis of this review is on algorithmic ideas, not necessarily on the underlying theory.
Another versatile framework to construct well-balanced methods is based on the reformulation of Eq. (1) as a homogeneous quasi-linear PDE system of the form This framework is applied in the context of hyperbolic systems with non-conservative products using the path-conservative finite volume methods (see, e.g., Cargo and LeRoux 1994;Greenberg and LeRoux 1996;Greenberg et al. 1997;Gosse 2000Gosse , 2001Parés and Castro 2004;Parés 2006;Castro et al. 2007LeVeque 2010). In this form, a special family of paths in phase space can be constructed such that a well-balanced method results. These paths can be obtained from the explicit knowledge of the solutions of the Riemann problems of Eq. (3), which may be difficult and expensive in general, or through a so-called generalized hydrostatic reconstruction technique by Castro et al. (2007. The latter is closely related to the piecewie steady reconstruction technique (Castro and Parés 2020). Furthermore, the framework is able to deal with singular source terms. However, we shall not pursue further the presentation of this theoretically pleasing and elegant framework in this text and redirect the interested reader to the given references. See also the recent comprehensive review by Castro et al. (2017) on this framework. Before we proceed to the outline, we also mention that well-balanced methods may be considered as part of the family of so-called structure-preserving methods. These methods are designed such that certain properties of the physical model (i.e., the balance law, the partial differential equation, etc.) are fulfilled in some form at the discrete level. Such properties may be in the form of so-called companion balance or conservation laws that are automatically satisfied at the analytical level. For example, the second law of thermodynamics puts admissibility criteria in the form of entropy conditions on flow discontinuities such as shock waves. The preservation of physical states, e.g., positive mass density, pressure or subluminal velocities. The rotational invariance of the equations of (magneto-)hydrodynamics implies the conservation of angular momentum. Faraday's law, together with the fact that magnetic monopoles have never been observed in nature, imply the solenoidal character of the magnetic field in Maxwell's equations and magnetohydrodynamics. Self-gravitating flows conserve total momentum and energy. Although consistent numerical methods may fulfill such structures in the infinite resolution limit, this is often unsatisfactory in practice as the needed resolution may result in unaffordable large computational costs. Also, note that near discontinuities, which solutions of balance laws prominently feature, these errors are inevitably of order one. Hence, structure-preserving methods are often not a luxury choice but a necessity. The design of structure-preserving methods that maintain a discrete form of such structures is a rich and challenging line of research on numerical methods for balance/conservation laws. We refer the interested reader to, e.g., Evans and Hawley (1988), Balsara and Spicer (1999), Tóth (2000), Morton and Roe (2001), Tadmor (2003), Mishra and Tadmor (2011), Jiang et al. (2013), Després and Labourasse (2015), Schaal et al. (2015), Katz et al. (2016), Zanotti and Dumbser (2016), Balsara and Kim (2016), Tang (2017, 2018), Mullen et al. (2021) and references therein.
The text is organized as follows: • Section 2 presents a brief introduction to finite volume methods and motivates the well-balanced methods on the basis of an extremely simple model equation, namely the linear advection-reaction equation. This is followed by a general framework for the construction of well-balanced finite volume methods. The procedure is examplified on the Euler equations in spherical symmetry featuring a geometric source term. The section rounds up with a general discussion of well-balanced methods within finite difference frameworks. • Section 3 focuses on well-balanced methods for the Euler equations. Several flavors of well-balanced methods are presented with differing local steady-state determination strategies. The section completes by a battery of numerical test problems on which the performance of well-balanced methods is commonly assessed.
Before we proceed, let us state that any inadvertent omission or understatement of credit to authors to whom more was due, we humbly offer a sincere apology in advance.
2 Well-balanced discretization 2.1 One-dimensional methods We begin by considering a one-dimensional system of balance laws in the form Here u, f and s are vectors of m components: u ¼ uðx; tÞ is the vector of conserved variables, f ¼ f ðuÞ the vector of flux functions, and s ¼ sðuÞ the vector of source terms. 2 In the following, we will tacitly assume that (i) the system is of hyperbolic nature: the Jacobian of the flux function vector AðuÞ ¼ of ou has real eigenvalues and an associated set of linearly independent eigenvectors for all u of interest, (ii) the source term sðuÞ is not singular: bounded source terms do not change the Rankine-Hugoniot jump conditions of the system.
Next, we outline a standard finite volume discretization of the balance law Eq. (4) in order to introduce our notation and set the stage for the following developments. For further details and precise derivation, we refer to the many excellent textbooks available in the literature, see, e.g., LeVeque (1992LeVeque ( , 2002, Godlewski and Raviart (1996), Laney (1998), Hirsch (2007), Toro (2009).

Finite volume discretization
A standard finite volume method discretizes the spatial domain of interest X ¼ ½0; L into a finite number N of control volumes or cells X i ¼ ½x iÀ1=2 ; x iþ1=2 (i ¼ 1; . . .; N). For the i-th cell, x iAE1=2 denote the left/right cell interfaces and x i ¼ ðx iÀ1=2 þ x iþ1=2 Þ=2 the cell centers. For ease of presentation, we shall assume a uniform discretization with constant cell size Dx ¼ x iþ1=2 À x iÀ1=2 . However, this assumption can easily be relaxed within a finite volume discretization.
Integrating the balance law Eq. (4) over cell X i and dividing by the cell size Dx yields where we introduced the cell averages of conserved variables and source terms. We also introduce the convention that a quantity with an overbar denotes a cell average, while one without a point value. Equation (5) represents an exact evolution equation for the cell-averaged conserved variables. The numerical approximation is introduced by replacing the exact fluxes and source terms by so-called numerical fluxes and source terms: Here, the U i , F iAE1=2 and S i denote approximations of the cell-averaged conserved variables, the fluxes through the cell interfaces and the cell-averaged source terms at time t: In the following, we use the convention that exact solutions are denoted by lower case letters and approximations by upper case letters. Equation (7) is a generic semi-discrete 3 finite volume discretization in one space dimension. Furthermore, in Eq. (7) we denote the so-called spatial discretization operator by LðUÞ i . Next, we briefly describe the individual components of a finite volume method. For ease of notation, we suppress the temporal dependency.

Reconstruction R
The primary unknowns in a finite volume method are the cell averages. To evaluate the numerical fluxes through cell interfaces and compute cell averages of the source terms, within each cell a subcell profile of the conserved variables U i ðxÞ has to be reconstructed from the cell averages fU k g. Because discontinuities may be present in the solution, special care is needed to reconstruct non-oscillatory subcell profiles that avoid spurious Gibbs phenomena.
We denote such a reconstruction procedure R, which recovers an r-th order accurate profile Q i ðxÞ of a quantity q(x) within cell X i from the cell averages fq k g, by with Here S i ¼ f. . .; i À 1; i; i þ 1; . . .g is the so-called stencil of the reconstruction for cell X i , which consists of cell X i and a certain number of neighboring cells. For systems, the reconstruction procedure can be applied component-wise to the cell averages of the conserved variables vector Numerous such reconstruction procedures have been developed in the literature, and a non-exhaustive list includes the Total Variation Diminishing (TVD) and the Monotonic Upwind Scheme for Conservation Laws (MUSCL) methods (see, e.g., van Leer 1979;Harten et al. 1983;Sweby 1984;Laney 1998;LeVeque 2002;Toro 2009), the Piecewise Parabolic Method (PPM) by Colella and Woodward (1984), the Essentially Non-Oscillatory (ENO) (see, e.g., Harten et al. 1987), Weighted ENO (WENO) (see, e.g., Shu 2009 and references therein) and Central WENO (CWENO) methods (see, e.g., Levy et al. 1999Levy et al. , 2000Cravero et al. 2018). For example, a spatially first-order accurate reconstruction consists of a piecewise constant profile A spatially second-order accurate piecewise linear reconstruction à la TVD/MUSCL is given by where DU i are some appropriately limited slopes (to avoid monotonicity violation and ensuing spurious oscillations). A popular example is the so-called generalized minmod slope limiter family where h 2 ½1; 2 is a parameter and the minmod function is defined by minmodða 1 ; a 2 ; :::Þ ¼ Equation (14) has to be understood component-wise.
reproduces the traditional minmod (monotonized centered) limiter (see, e.g., Toro 2009 and references therein for further information). See Fig. 1 for an illustration of a piecewise constant/linear reconstruction. Straightforward component-wise reconstruction for systems of balance laws may sometimes lead to some undesirable oscillations in the results, especially when strong flow discontinuities interact. In that case, it may prove beneficial to perform the reconstruction in local characteristic variables (see, e.g., Harten et al. 1987;Qiu and Shu 2002;Toro 2009) where L i ¼ LðU i Þ and R i ¼ RðU i Þ are the matrices of left and right eigenvectors, respectively. The eigenvectors are typically evaluated at the cell average of the i-th cell X i whose reconstruction is performed, hence the name local.

Numerical fluxes F
The numerical fluxes are obtained by resolving the discontinuities at cell interfaces naturally arising from the per cell reconstruction (see Fig. 1). This is commonly done à la Godunov by solving (approximately) the Riemann problem at cell interfaces where the point values U iþ1=2AE are the cell interface reconstructed conserved variables Notice that the value on the left (L)/right (R) of the interface x iþ1=2 is obtained from the reconstruction in cell X i /X iþ1 . The numerical flux is required to be consistent with the physical flux function f , i.e., F ðU; UÞ ¼ f ðUÞ, and Lipschitz continuous. The latter is required for accuracy reasons (see, e.g., Harten et al. 1987). Moreover, certain numerical fluxes have the ability to exactly recognize isolated discontinuities such as contacts or shocks in (magneto-) hydrodynamics (see, e.g., Toro 2009 for details). A simple and popular choice for the numerical flux is the so-called Rusanov flux (Rusanov 1962;Toro 2009) Fig. 1 Illustration of the reconstruction procedure for the sine function uðxÞ ¼ sinðxÞ (solid black line). Within each control volume or cell is shown the cell average U i (solid blue) and a piecewise linear TVD/ MUSCL reconstruction U i ðxÞ (solid red). Notice the limiter's action near the two extrema, where the slopes are clipped to ensure the monotonicity of the reconstruction. The piecewise constant cell averages correspond to a piecewise constant reconstruction where F L=R ¼ f ðU L=R Þ and S max is an estimate of the largest characteristic speed in the solution of the Riemann problem S max ¼ max m jk m j (k m are the eigenvalues of the flux Jacobian). This numerical flux is sometimes also called a local Lax-Friedrichs (LLF) flux.

Numerical source terms S
For the integration of the source terms, there are essentially two standard methods. The first one is the so-called unsplit method. In this method, the source terms are typically incorporated directly into the spatial discretization operator as tacitly already done in Eq. (7). An accurate approximation of the cell average of the source terms are obtained by numerical integration. Let Q i denote a q-th order accurate quadrature rule over the i-th cell X i . The cell average of the source terms are then computed by where the x i;a 2 X i and x a denote the N q quadrature nodes and weights of Q i , respectively. Assuming that the point values of the source terms can be evaluated with spatial order of accuracy s, then the resulting discretization is spatially minðq; sÞ-th order accurate (provided enough smoothness, of course). 4 A popular example is the second-order accurate midpoint rule Higher-order rules are provided, for example, by the Gauss-Legendre or Gauss-Lobatto quadrature rules (see, e.g., Press et al. 1993). The second family of methods are the so-called splitting or fractional-step methods. In these methods, the original problem Eq. (4) is first split (or fractured) into two subproblems of the form: In this approach, one alternates adroitly between solving the two subproblems A and B. This is indeed a very practical approach: Problem A is a standard (homogeneous) conservation law and Problem B is a simple Ordinary Differential Equation (ODE). For both subproblems, there exist many excellent numerical methods and software libraries which implement them. By extension, this approach allows for straightforward modularization. Next, we catalogue two popular splitting methods. Let S Dt A and S Dt B denote the discrete solution operators that advance a discrete solution U n by a time step Dt for problems A and B where the superscript labels the time step and the ''star'' shall stress the fact that these are only partially evolved states. An obvious splitting method is then given by which is first-order accurate in time (provided that S Dt A and S Dt B are at least of that same temporal order). This splitting is sometimes termed as Godunov splitting. A second-order accurate in time method is given by the so-called Strang splitting where one sandwiches a full step Dt with solution operator A between two half steps Dt=2 with solution operator B. Of course, full second-order accuracy in time requires that the individual solution operators S Dt A and S Dt B possess an equivalent or higher order of accuracy.

Time discretization T
The semi-discrete evolution equations Eq. (7) for the cell averages U i represent a system of ordinary differential equations that has to be approximately integrated in time. For that purpose, the temporal domain of interest T ¼ ½t i ; t f is discretized into time steps Dt n ¼ t nþ1 À t n , where the superscripts label the respective time step.
The simplest time integration method is of course the temporally first-order accurate explicit Euler method For higher-order time integration, there are essentially two large families of methods: Runge-Kutta and predictor-corrector methods. A popular representative of the Runge-Kutta family is the temporally second-order accurate explicit Heun method It is a so-called Strong Stability-Preserving Runge-Kutta (SSP-RK) method, and it is often labeled by SSP-RK2 in the literature (because it is a two-stage SSP-RK method). These methods have certain desirable stability properties when integrating non-linear conservation/balance laws (see, e.g., Gottlieb et al. 2001 and references therein). Another very popular method is the third-order accurate SSP-RK3 method, which is shown below in Eq. (126). A popular representative of the predictor-corrector methods is the temporally second-order accurate MUSCL-Hancock method (van Leer 1984). Possible highorder extensions of this methodology can be achieved by evolution of the solution in the small with help of the Cauchy-Kowalevski procedure (Harten et al. 1987). Another possibility is provided by the so-called Arbitrary DERivative (ADER) methods (see, e.g., Castro and Toro 2008;Toro 2009 and references therein). A distinctive feature of predictor-corrector methods is that they are one-step methods, which makes them extremely attractive in an adaptive mesh refinement context (see, e.g., Balsara 2017).
For time explicit approaches as above, Eqs. (26) and (27), the time step Dt n is in general required to fulfill a so-called CFL condition of the form (Courant et al. 1928) where S n i is the speed of the fastest wave in cell X i at time t n and C CFL is the CFL number. The latter needs to fall within a certain range for linear stability.
Time implicit approaches are also possible and especially adapted for so-called stiff problems involving vastly different timescales (see, e.g., Kwatra et al. 2009;Viallet et al. 2011Viallet et al. , 2013Kifonidis and Müller 2012;Miczek et al. 2015 and references therein). Stiffness may originate in the different wave propagation characteristics (such as advective versus acoustic waves), strong chemical reactions and many more.

Assembling a finite volume scheme
A generic finite volume scheme Eq. (7) for the one-dimensional balance law Eq. (4) is now easily assembled with the previously described components: (1) A spatially r-th order accurate reconstruction R (Eq. (11)).
(3) An unsplit source terms discretization S (Eq. (20)) based on s-th order accurate point value evaluations and a q-th order accurate quadrature rule Q. (4) A s-th order accurate time integrator T . This results in a minðq; r; s; sÞ-th order accurate finite volume scheme (for smooth enough solutions, of course). A similar assemblage can be realized with an appropriate splitting method for the source terms.
However, the approximation of near steady states characterized by the near balance of flux divergence and source term is quite challenging for such a generic finite volume scheme It turns out that the steady states of interest are not exactly representable by polynomials used in the reconstruction procedure in general. Therefore the piecewise polynomial reconstruction will introduce truncation errors at every time step. Another issue is that the flux divergence and source term discretizations are commonly computed independently. This further makes the above discrete near balance unlikely.
In the next section, we present a simple illustrating example followed by a general technique to well-balance such steady states within a finite volume framework.

Example: linear advection-reaction equation
We now illustrate the issues that can arise when numerically approximating balance laws near steady states. Consider the simple linear advection-reaction equation modeling the transport of some radioactive material of concentration u(x, t) with constant advection velocity a [ 0 and decay rate k [ 0. An exact solution is easily derived with the method of characteristics uðx; tÞ ¼ e Àkt u 0 ðx À atÞ; where u 0 ðxÞ is the initial concentration uðx; 0Þ ¼ u 0 ðxÞ. The exact solution reflects the anticipated behavior that the initial concentration is advected to the right with velocity a and decays along the way with rate k. An interesting feature of the above simple model Eq. (29) is that it possesses nontrivial steady-state solutions which are of the form for some constant C. The steady states are characterized by a subtle balance between the advection and decay processes. Their spatial variation is ruled by the ratio between the decay and the advection timescale, commonly known as the Damköhler number Da ¼ k a . Let us solve approximately Eq. (29) over the computational domain X ¼ ½0; 2 discretized by N uniform cells X i (i ¼ 1; . . .; N). For illustration, we choose two first-order accurate finite volume schemes. The first scheme consists of piecewise constant reconstruction Eq. (12), the Rusanov numerical flux Eq. (19), the unsplit source term discretization based on the midpoint rule Eq. (21), and the explicit Euler time integration Eq. (26). Explicitly, this gives the following fully discrete evolution equation for the cell-averaged concentration U n i within cell X i : The second scheme employs Godunov splitting Eq. (24) for the source term discretization and reads: Note that explicit Euler time integration is used in both subproblems. In principle, the exact solution could be used in Eq. (34b), i.e., U nþ1 i ¼ e ÀkDt U Ã i . However, this does not affect the following discussion. Without any surprise, an astute reader will recognize here the classical first-order upwind method for the linear advection equation. Both schemes are first-order accurate in space and time and are linearly stable provided that the time step Dt is chosen such that 0\ Dt Dx a 1 and 0\Dtk\2. We fix a ¼ k ¼ 1 and evolve a slightly perturbed steady state Eq. (32) as shown in Fig. 2a for one time unit. The small perturbation centered around x ¼ 0:5 is advected by one unit to the right and its amplitude decays by a factor e À1 . In the same panel are also shown the approximate results obtained with the unsplit Eq. (33) and split Eq. (34) first-order schemes. Both schemes show qualitatively correct results. More quantitatively, Fig. 2b displays the equilibrium perturbation, that is, the difference between the solution and the background steady state. We observe that the perturbation is indeed advected by the correct distance by both schemes. However, we also observe that both schemes show significant discrepancies with the expected solution away from the perturbation.
To further highlight the issue, we evolve the unperturbed steady state Eq. (32) for one time unit with both schemes. The results are shown in Fig. 3a. By comparison with Fig. 2b, we see clear evidence that the spurious deviations are due to the inability of both schemes to maintain the steady state discretely. To illustrate the origin of the problem, Fig. 3b shows the exact steady state together with the cell averages U 0 i at the initial time for a few cells. The cell averages also correspond to the piecewise constant solution representation within each cell resulting from the first-order reconstruction. It is clear that these piecewise constant subcell profiles are inadequate to represent the steady state within the cells faithfully. More precisely, the piecewise constant approximation of the exponentially varying steady state Eq. (32) inevitably introduces truncation errors of order OðDxÞ.
Likewise, higher-order polynomial reconstruction procedures of order r introduce truncation errors of order OðDx r Þ. Therefore, the schemes will introduce local truncation errors of order OðDx r Þ near non-polynomial steady states. If the goal is to simulate small perturbations on top of a steady state, the numerical resolution needs to be increased to the point that these local truncation errors do not obscure the phenomena of interest. Similarly, if the goal is to simulate phenomena near a steady state for an extended time (compared to a characteristic timescale on which the steady state would react to equilibrium perturbations), the resolution needs to be increased such that the pile-up of these local truncation errors in each time step does not corrupt the phenomena of interest. This increase in resolution may cause prohibitively high computational costs, especially in multiple dimensions.
This inadequacy of standard piecewise reconstruction procedures was realized early on in the development of numerical methods for balance laws. This motivated for example Liu (1979), Glaz andLiu (1984), andvan Leer (1984) to replace the piecewise constant reconstruction within each cell by a piecewise steady reconstruction Note that this equilibrium subcell profile U n eq;i ðxÞ depends on the cell under consideration and may be adapted in each time step. Hence the subscripts ''eq'' and ''i'', and the superscript ''n''. Since the steady states of the considered linear advectionreaction equation are known explicitly Eq. (32), the desired piecewise steady reconstruction is of the form U n eq;i ðxÞ ¼ C i e Àk=aðxÀx i Þ ; x 2 X i ; and the constant C i is simply fixed by matching with the i-th cell's average Eq. (38) Plugging this reconstruction into the unsplit first-order scheme gives Analogously for the split first-order scheme, one obtains Let's evolve the unperturbed steady state with the split and unsplit schemes using the above piecewise steady reconstruction Eq. (36). The resulting equilibrium perturbation at final time t f ¼ 1 is shown in i (blue solid lines) at initial time for a few cells. The latter also corresponds to the piecewise constant reconstruction on which the unsplit/split first-order schemes are based observe that the piecewise steady reconstruction does not improve the situation for the split scheme. Actually, the spurious equilibrium deviations are even slightly worse in this example. In contrast, the unsplit scheme with piecewise steady reconstruction preserves the steady state down to machine precision (% 10 À16 for the double precision floating-point representation used in the computations). Figure 4b displays the results for the slightly perturbed steady state. The split scheme evolves the perturbation faithfully, but is afflicted by the scheme's local truncation errors at the steady state. On the other hand, the unsplit scheme not only advect the bump very well; it additionally relaxes back to the steady state once the perturbation passed through.
A straightforward computation shows that the unsplit first-order scheme with piecewise steady reconstruction is exact for the advection-reaction equation's steady states Eq. (32). Cargo and LeRoux (1994) subsequently coined the term wellbalanced for numerical schemes with the property of preserving a discrete form of certain steady states exactly. In the above derivation, we implicitly fixed some choices within the scheme. When fixing the equilibrium subcell profile U eq;i ðxÞ, we chose exact integration in Eq. (38) to match with the i-th cell average. Instead, a quadrature rule, e.g., the midpoint rule, could be used. Similarly, we chose the exact solution Eq. (39) as the equilibrium subcell profile. As suggested by Roe (1987), Mellema et al. (1991) and Eulderink and Mellema (1995), one could chose an approximate solution Eq. (37) as the equilibrium subcell profile. In the next section, we present a general framework for the construction of well-balanced high-order finite volume schemes. At the root, it is based on a high-order generalization of the piecewise steady reconstruction idea.
We remark that fractional step or splitting methods could also be adapted to improve their performance near steady states. This can be achieved by carefully matching the boundary conditions used in the conservation law evolution Eq. (22a) and the source term integration Eq. (22b). However, we shall not pursue this idea in (a) (b) Fig. 4 The left panel shows the equilibrium perturbation for the numerically evolved steady state with a resolution of N ¼ 64 cells at final time t f ¼ 1. The blue and red pluses are obtained with the unsplit/split first-order schemes using the piecewise steady reconstruction. The zoom in shows that the unsplit scheme is able to preserve the steady state down to machine precision. The right panel shows the equilibrium perturbation for the slightly perturbed steady-state concentration simulated with the unsplit/ split first-order schemes using the piecewise steady reconstruction. The unsplit scheme clearly relaxes back to the steady state away from the concentration perturbation the sequel, and we refer to LeVeque (1986LeVeque ( , 2002 and references therein for a general procedure.

Well-balanced finite volume schemes
From the linear advection-reaction example, we see that the idea of a piecewise steady solution representation can lead to a finite volume scheme capable of preserving a steady state exactly, i.e., a well-balanced finite volume scheme. In the following, we present a simple recipe for constructing arbitrarily high-order wellbalanced finite volume schemes in a systematic manner. The recipe relies on a highorder generalization of the piecewise steady reconstruction and a special discretization of the source terms guaranteeing the discrete preservation of steady states. We stress that the recipe is a humbly distilled version of the methodologies found in the vast literature about well-balanced finite volume schemes given in Sect. 1. The principle of well-balanced finite volume methods based on piecewise steady reconstruction is to decompose the solution into an equilibrium part and a (not necessarily small 5 ) perturbation part where the equilibrium part u eq ðxÞ fulfills the steady-state balance One obvious requirement for the piecewise steady reconstruction is thus the ability to compute such steady states. It turns out that this is difficult in general. However, we will tacitly assume that the differential equation Eq. (44) can be solved (exactly or approximately) for certain steady states U eq ðxÞ ¼ u eq ðxÞ þ OðDx Þ: Here denotes the spatial order of accuracy of the computed equilibrium solution. If Eq. (44) can be solved analytically for certain steady states, we may slightly abuse the notation and set ¼ 1.
Solving for equilibrium is usually the main challenge when designing a wellbalanced scheme. The difficulty depends strongly on the considered balance law and the associated steady states. For the linear advection-reaction equation, there is only one steady state Eq. (39) and it is known analytically. The Euler equations of fluid dynamics feature a myriad of steady states. We will look at one particular class that arises when considering the Euler equations in spherical symmetry in Sect. 2.5. Section 3 will look at the steady states that occur when the considered fluid is subject to gravitational forces. Luckily, one is generally not interested in all possible steady states within one practical simulation. Hence, it is often sufficient to design well-balanced schemes for certain stationary states of practical interest. The solvability of Eq. (44) is then restricted to these cases to which we will refer loosely as the steady states of interest in the following.
Next, we describe the modifications to the standard reconstruction and source term integration procedures of Godunov-type finite volume schemes for the homogeneous equations to build a well-balanced scheme for the steady states of interest. This involves the subtle correction of the reconstruction and source term integration that somehow incorporates the steady states of interest. We repeat the obvious that the following developments evidently hinge on the computability of the steady states of interest.

Piecewise steady reconstruction WR
As in a standard finite volume scheme, within each cell a subcell profile U i ðxÞ has to be reconstructed from the cell averages fU k g. A piecewise steady reconstruction WR is then given by the decomposition where U eq;i ðxÞ and dU i ðxÞ denote the local equilibrium and perturbation reconstruction parts in cell X i , respectively. The stencil of the piecewise steady reconstruction is denoted as previously by S i . We now describe each part in detail. Within each cell X i , the local equilibrium reconstruction U eq;i ðxÞ is determined by fitting an equilibrium solution U eq ðxÞ among the steady states of interest to the cell average U i . Since the cell average U i may be arbitrarily far from a steady state of interest, this is done in two substeps. The first substep consists of projecting U i onto a cell average U eq;i consistent with the steady states of interest. The second substep determines the local equilibrium reconstruction U eq;i ðxÞ in cell X i by matching an equilibrium profile Eq. (45) to the equilibrium projected cell average U eq;i where Q i denotes a q-th order accurate quadrature rule over cell X i . We allow that the matching Eq. (47) is done exactly (using exact integration) and again slightly abuse the notation by setting q ¼ 1. For instance, the matching was done exactly with Eq. (38) in the linear advection-reaction example of Sect. 2.3. This results in a minð; qÞ-th order accurate local equilibrium reconstruction within each cell. However, the difficulty of this equilibrium projection and matching depends strongly on the balance law and steady states of interest. Some concrete examples are provided in Sects. 2.5 and 3. In addition, it is important to realize that not every given cell average must correspond to an equilibrium among the steady states of interest. Indeed, the solution may be far from a steady state. Therefore, the possibility that the local equilibrium reconstruction does not succeed must be taken into account. In that case, the local equilibrium reconstruction is simply set to zero U eq;i ðxÞ 0. The local equilibrium perturbation dU i ðxÞ within each cell X i is obtained by extrapolating the cell's local equilibrium profile U eq;i ðxÞ to neighboring cells, where it is compared with their cell averages. This senses how much the neighboring cells are perturbed with respect to the equilibrium in cell X i . Cell-averages of these equilibrium perturbations can then be fed to any standard r-th order accurate piecewise polynomial reconstruction procedure to recover a local equilibrium perturbation profile as Like for the standard reconstruction procedure Eq. (16), the equilibrium perturbation reconstruction can also be performed in local characteristic variables. The piecewise steady reconstruction WR (Eq. (46)) is illustrated in Fig. 5 for a scalar quantity. The reconstruction is minð; q; rÞ-th order accurate close or far from the steady states of interest (for smooth enough solutions, of course). It is intuitively clear that the local equilibrium perturbation dU i ðxÞ vanishes if cell averages of a steady state of interest are fed to the piecewise steady reconstruction. Hence, a minð; qÞ-th order accurate discrete form of the steady states of interest is exactly reconstructed by the piecewise steady reconstruction. A proof is sketched in Sect. 2.4.3. Also note that if we set U eq;i ðxÞ 0, then WR automatically reduces to the standard piecewise polynomial reconstruction procedure R. This is important in practice when there exists no solution of Eq. (47), i.e., no local equilibrium solution matching with the given cell's average is found.
A possible variation of the piecewise steady reconstruction found in the literature is given by Sketch of the piecewise steady reconstruction of some scalar quantity from the cell averages fU k g.
For the illustration, we assume that the steady-state projection of the cell averages is trivial U eq;i ¼ U i as for the steady states of the linear advection-reaction equation. Left panel: At the beginning of the reconstruction process, we are given the cell averages of the i-th cell and its immediate neighbors (solid black piecewise constant lines). The equilibrium reconstruction U eq;i ðxÞ is built such that it matches the cell average U i by Eq. (47) and extrapolated to the neighboring cells k ¼ i AE 1 (solid blue line). The cell averages of U eq;i ðxÞ are then computed in the neighboring cells k ¼ i AE 1 (solid blue piecewise constant lines) and the cell-averaged equilibrium perturbations as seen from the i-th cell are computed. Right panel: The equilibrium perturbation dU i ðxÞ is reconstructed by a standard reconstruction procedure Eq. (48). Finally, by combining the equilibrium U eq;i ðxÞ and perturbation dU i ðxÞ reconstruction as Eq. (46) one obtains the equilibrium-preserving piecewise steady reconstruction U i ðxÞ in the left panel (solid black line) which separates the solution into local equilibrium and relative perturbation parts (e.g., Chandrashekar and Klingenberg 2015;Berberich et al. 2019). The local equilibrium reconstruction U eq;i ðxÞ is obtained with the same two substeps as above.
The local relative equilibrium perturbation is computed by where the expression is to be understood component-wise. One drawback of this form is that it does not automatically reduce to a standard reconstruction if (some components of) the local equilibrium U eq;i ðxÞ vanishes. In that case, one simply switches to a standard reconstruction (of these components) with some additional implementation logic. If the reconstruction is not sensitive to the shift with a constant, for any constant C and cell averages fQ k g, then Eq. (49) can be rewritten as and the relative equilibrium perturbation Eq. (50) as Most reconstruction methods possess property Eq. (51) because non-oscillating behavior is usually enforced by limiting first and higher derivatives of the reconstruction polynomial and these are not affected by the addition of a constant. Both forms of the relative piecewise steady reconstruction share similar properties to the ''absolute'' one Eq. (46). An example is given in Sect. 3.3.1.

Well-balanced source term discretization WS
A direct numerical integration of the source term as in Eq. (20) will in general not lead to a well-balanced scheme. Instead, one uses the previously introduced piecewise steady reconstruction that decomposes the solution into an equilibrium and a perturbation part to perform the following seemingly frivolous manipulation which simply adds and subtracts the source term evaluated with the local equilibrium reconstruction U eq;i within cell X i . As suggested, e.g., by Huang and Liu (1986), Audusse et al. (2004), and Botta et al. (2004), the equilibrium part fulfills the steady-state balance by construction, and can be trivially integrated by applying the fundamental theorem of calculus. The cell average of the source term Eq. (54) can therefore be approximated by applying exact integration to the equilibrium part and numerical integration to the perturbation part as follows Here Q i denotes a q-th order accurate quadrature rule over cell X i . Note that Q i may be different from the quadrature rule used in the piecewise steady reconstruction. We will refer to it with the same symbol since the same quadrature rule is typically used. Equation (56) results in a minð; q; r; sÞ-th order accurate discretization of the source term. At a steady state of interest (i.e., U i U eq;i ), Eq. (56) reduces to As we will see below, this is crucial for the well-balanced property of the scheme. Moreover, note that if the local equilibrium reconstruction part U eq;i ðxÞ vanishes, the above source term discretization automatically reduces to the standard one in Eq. (20). This is again important in practice when no local equilibrium matching with the given cell average is found (i.e., no solution to Eq. (47) is found). An alternative form of the well-balanced source discretization is based on Richardson extrapolation. The idea is to write the source term as which has to be understood component-wise. The fact that the equilibrium part U eq;i fulfills the steady-state balance by construction is used in the second equality. However, note that rewriting the source term in this way may not be possible for all the components of a particular system of balance laws, because they are trivially fulfilled at the steady states of interest. Hence, they are not relevant for the construction of a well-balanced scheme and can be discretized in a standard way. For the sake of presentation, we ignore this subtlety in the derivation of the wellbalanced source term discretization based on this form. An illustrative example of this alternative well-balanced source term discretization is provided in Sect. 3.3.1. Consider the following second-order approximation of the cell-averaged source term based on the form above where we introduce the symbol T i for this particular quadrature rule (due to its resemblance with the trapezoidal rule). At a steady state of interest U i U eq;i , this clearly reduces to Eq. (57) which is again crucial for well-balancing as we shall see below. Unfortunately, it is still only a second-order source term discretization. To overcome this limitation, Noelle et al. (2006) ingeniously suggest to use Richardson extrapolation. Let us introduce a composite quadrature rule based on Eq. (59). The cell By applying the quadrature rule T i to each subinterval and summing up, we obtain the following composite quadrature rule This again reduces to Eq. (57) at a steady state of interest by telescoping of the sum, but it is still only second-order accurate. However, Noelle et al. (2006) note that the quadrature rule T i is also symmetric and therefore possesses an asymptotic error expansion of the form for any (smooth enough) function f. Richardson extrapolation then combines the T N c i to cancel out increasingly higher error terms in the expansion. For example, fourthand a sixth-order accurate quadrature rules are readily obtained: Thus, arbitrary high-order well-balanced source term discretizations can be obtained from the alternative form Eq. (59). Although it possesses similar properties as the well-balanced source term discretization Eq. (56), one drawback of the alternative form is that it does not automatically reduce to a standard source term discretization in case the local equilibrium part vanishes in the piecewise steady reconstruction (i.e., no solution to Eq. (47) is found). However, this can easily be handled with some additional implementation logic.
This results in a minð; q; r; s; sÞ-th order accurate well-balanced finite volume scheme (for smooth enough solutions). The scheme preserves exactly a minð; qÞ-th order accurate discrete form of the steady states of interest (up to machine precision). Furthermore, such a well-balanced scheme automatically falls back to a standard high-order finite volume scheme if the local equilibrium reconstruction part vanishes. 6 This guarantee is important in practice since one is assured that if the piecewise steady reconstruction fails to determine a local equilibrium profile (because it may not exist), the scheme reduces decently to a standard scheme without any loss of accuracy and robustness.
To round off this section, we demonstrate the well-balanced property of such a scheme, that is, its ability to preserve exactly (up to machine precision) a discrete form U eq ðxÞ of the steady states of interest u eq ðxÞ it was designed for. 7 For simplicity, we assume that both u eq ðxÞ and its approximation U eq ðxÞ are continuous. As we shall see below, this requirement can easily be waived. Let the scheme be given cell averages fU i g of such a steady state of interest. These cell averages are computed from a given steady state u eq ðxÞ approximated discretely by U eq ðxÞ with where Q i denotes the same q-th order quadrature rule as used when matching the local equilibrium profile with the cell averages in the piecewise steady reconstruction Eq. (47). We shall term such initial data as well-prepared initial data. Given such well-prepared initial data, we reciprocally assume that the local equilibrium reconstruction Eq. (47) recovers within every cell X i the restriction of U eq ðxÞ in respective cell, and that its extrapolation over the computational domain X recovers U eq ðxÞ everywhere, i.e., Of course, this assumption needs to be verified for the particular balance law and steady states of interest, and represents the core challenge in the construction of a well-balanced scheme. Taking this assumption for granted, it is obvious that the piecewise steady reconstruction WR (Eq. (46)) recovers the given steady state U eq ðxÞ in every cell as the local equilibrium perturbation vanishes everywhere, i.e., dU i ðxÞ 0, and we have The same holds true for the alternative piecewise steady reconstructions Eq. (49) or Eq. (52)) with slightly adapted arguments. For the numerical flux, we therefore have due to the fitting of the piecewise steady reconstruction WR at every cell interface Similarly for the source term discretization WS Eq. (56), we have The alternative source term discretization Eqs. (60) and (62) likewise reduces to the above expression. Plugging Eqs. (67) and (69) into Eq. (7), one obtains and therefore the scheme is well-balanced as claimed, i.e., it preserves a minð; qÞ-th order accurate discrete form of the steady states of interest exactly (up to machine precision). We assumed that the steady states of interest and their discrete approximation are continuous in Eq. (68). However, it is straightforward to generalize the wellbalanced scheme (and the above demonstration) to steady states with stationary discontinuities located at cell interfaces. For this purpose, some mild additional requirements for the numerical flux function F and the standard piecewise polynomial reconstruction procedure R are necessary: (i) the numerical flux F has to be able to resolve exactly the (stationary) discontinuities allowed by the steady states of interest, and (ii) the reconstruction procedure R has to reduce to a piecewise constant reconstruction near isolated discontinuities at cell interfaces. Both requirements ensure that at the stationary discontinuities the numerical flux agrees with the exact flux like in Eq. (66).

Example: the Euler equations in spherical symmetry
As a simple and practical example of the construction of a well-balanced finite volume scheme, we consider the Euler equations in spherical symmetry expressing the conservation of mass, momentum and energy. Here r is the radial coordinate, q the mass density, v the radial velocity, E ¼ qe þ 1 2 qv 2 the total fluid energy composed of internal and kinetic energy densities, and p the pressure. The latter is related to the density q and specific internal energy e through an equation of state p ¼ pðq; eÞ.
A computationally convenient form of Eq. (70) is given by where the vector of conserved variables, fluxes and source terms are and the area and volume functions are This form is particularly convenient because the fluxes take exactly the same form as in the one-dimensional planar geometry case. Hence, the same numerical fluxes can directly be used. The drawback of this form is the introduction of a geometric source term that becomes singular near the origin. Furthermore, note that the source term may depend non-linearly on the conserved variables through the equation of state.
A particular steady state of Eq. (70) and Eq. (71) is a resting fluid with uniform density and pressure profile. It is of course highly desirable that a numerical scheme faithfully reproduces this seemingly trivial equilibrium. The steady state of interest u eq is therefore simply where q eq ¼ const is the constant density and qe eq ¼ const the constant internal energy density, respectively. It clearly fulfills where p eq ¼ pðq eq ; qe eq Þ ¼ const is the constant equilibrium pressure. A straightforward discretization of Eq. (71) on a spherical domain D ¼ ½R 0 ; R 1 , 0 R 0 \R 1 , with a semi-discrete finite volume method gives for the i-th cell Here U i denotes the approximate cell average of the conserved variables over a (spherical shell) cell X i ¼ ½r iÀ1=2 ; r iþ1=2 with inner/outer radius r iAE1=2 of volume The fluxes through the inner/outer (spherical shell) cell boundary of area A iAE1=2 ¼ Aðr iAE1=2 Þ are approximated by a numerical flux function e.g., the Rusanov flux Eq. (19), and the cell interface extrapolated point values of the conserved variables U iAE1=2À =U iAE1=2þ are obtained from a reconstruction procedure. For simplicity of the example, let's fix spatial accuracy to second order by choosing a piecewise linear reconstruction centered at the (spherical shell) cell where S i ¼ fi À 1; i; i þ 1g is the stencil and the limited slopes can be computed with the generalized minmod slope Eq. (14). Accordingly, we also choose the second-order accurate midpoint quadrature rule 9 for approximating integrals of a function f over a (spherical shell) cell X i : This immediately gives the following (naive) discretization of the geometric source term where p i is the pressure at the cell center. The latter is simply obtained by evaluating the piecewise linearly reconstructed conserved variables Eq. (79) at cell center where the peculiarity that cell averages correspond to point values at cell center up to second-order accuracy is especially apparent. This concludes the description of a plain-vanilla finite volume scheme for the Euler equations in spherical symmetry. We now construct a well-balanced finite volume scheme capable of preserving a resting fluid Eq. (74) exactly following the recipe in Sect. 2.4 based on the just described scheme. First, we need to devise a piecewise steady reconstruction procedure for the resting fluid equilibrium, i.e., our steady state of interest we wish to preserve. We begin with the local equilibrium reconstruction part. The first substep in the local equilibrium reconstruction is the projection of the cell averages onto equilibrium cell averages consistent with the resting fluid equilibrium. This substep is necessary because the averages could be arbitrarily far from the steady state of interest (i.e., non-vanishing radial momentum and kinetic energy densities), and it is simply accomplished by The equilibrium cell average of the density is simply set to the cell-averaged density and the equilibrium cell average of the momentum density is set to zero as is consistent with the steady state of interest. An expression for the cell average of the internal energy density is provided by Eq. (83). Although this is only spatially second-order accurate in general, it becomes exact when the fluid is at rest, thereby establishing consistency with the steady state of interest Eq. (74). The second substep is then to match a local equilibrium reconstruction U eq;i ðrÞ to the cell's X i average equilibrium projected conserved variables U eq;i as in Eq. (47). This is indeed trivial given the constant nature of the considered steady state of interest The local equilibrium perturbation reconstruction Eq. (48) is then where we used that U eq;i ðrÞ is a simple constant, i.e., As result, we obtain the following piecewise steady reconstruction In the last equality, we used again that U eq;i ðrÞ is simply a constant. However, the above piecewise steady reconstruction can be much simplified. The limited slopes DdU i can be reduced with the following observation which means that the equilibrium U eq;i drops out in the computation of the slopes Eq. (14). Hence, we have that the limited slopes of the local equilibrium perturbation reconstruction in Eq. (86) reduce to the slopes used in the standard piecewise linear reconstruction Eq. (79): DdU i ¼ DU i . Now by combining this result with Eq. (86) and plugging it into Eq. (88), we obtain that the piecewise steady reconstruction simplifies to the standard piecewise linear reconstruction: Of course, this is not surprising as we simply subtract a constant from the data to be (piecewise linearly) reconstructed. It is nevertheless a welcome simplification when implementing the present scheme. Finally, we construct the appropriate well-balanced source term discretization with Eq. (56): We substituted the midpoint rule Eq. (80) in the second equality, and we used the fact that the pressure computed from the piecewise steady U i ðrÞ and the local equilibrium reconstruction U eq;i ðrÞ coincide at the cell center r i in the third equality (see Eqs. (82) and (83)).
It is now straightforward to show that the just derived source term discretization Eq. (91) is indeed able to preserve a resting fluid with uniform density and pressure exactly. Hence, we have designed a well-balanced scheme for this particular steady state. This is confirmed by the numerical results displayed in Fig. 6.
We remark that the above spatially second-order accurate source term discretization Eq. (91) is well-known among the practitioners in the field (see, e.g., Mönchmeyer and Müller 1989;Li 2003;Skinner and Ostriker 2010;Wang and Johnsen 2013). The above expression for the geometric source term can also be motivated from the derivation of the momentum equation Eq. (72) which expresses the pressure gradient in Eq. (70b) with the following simple application of the product rule Expressing oA oV in a discrete finite volume sense then immediately gives Eq. (91). However, the here described discretization is in principle extensible to arbitrary spatial orders of accuracy. It would be interesting to combine the above with highorder reconstruction procedures for orthogonal curvilinear coordinates devised by Mignone (2014) and Shadab et al. (2019) together with specifically designed weighted 10 Gauss quadrature rules.

Extension to several space dimensions
We now extend the one-dimensional recipe in Sect. 2.4 to build arbitrarily highorder well-balanced finite volume schemes for multidimensional systems of balance laws where u ¼ uðx; tÞ is the vector of conserved variables, f ¼ fðuÞ the flux tensor and s ¼ sðuÞ the vector of source terms. As in the one-dimensional case, we tacitly assume that (i) the system is of hyperbolic nature (the Jacobian of the flux tensor n Á of ou is diagonalizable with real eigenvalues for any direction n) and that (ii) the source term is not singular. For ease of presentation, we focus on the two-dimensional case in Cartesian coordinates where f ¼ f ðuÞ and g ¼ gðuÞ are the vectors of fluxes in x-and y-direction, i.e., the components of the flux tensor f ¼ ½f; g T in Cartesian coordinates. However, the extension to three dimensions and other coordinate systems is straightforward. In the following subsection, we begin by concisely describing a standard finite volume discretization of the balance law Eq. (93) to introduce our notation. More comprehensive descriptions can be found in the excellent textbooks listed at the end of Sect. 2.1. The extension of the one-dimensional recipe to design well-balanced schemes in several space dimensions is presented in the subsequent subsections.

Finite volume discretization
We consider a rectangular spatial domain X ¼ ½x min ; x max Â ½y min ; y max discretized uniformly (for ease of presentation) by N x and N y cells or finite volumes in x-and y- Fig. 6 The figure shows the radial velocity after one time unit of a constant state with unit density and pressure (q ¼ p ¼ 1). The blue/red line corresponds to the results obtained by the standard naive/wellbalanced second-order schemes from Sect. 2.5 with a resolution of N ¼ 64. In both simulations, solid wall boundary conditions were enforced at the upper boundary. From the plot it is obvious that the standard treatment of the geometric source term results in spurious velocity fluctuations near the origin. In contrast, the well-balanced scheme shows a still standing radial profile, as is expected direction, respectively.
The cells are labeled by X i;j ¼ X i Â X j ¼ ½x iÀ1=2 ; x iþ1=2 Â ½y jÀ1=2 ; y jþ1=2 , the constant cell sizes by Dx ¼ x iþ1=2 À x iÀ1=2 and Dy ¼ y jþ1=2 À y jÀ1=2 , and the cell volumes by X i;j ¼ DxDy. We also introduce a non-directional cell size h ¼ maxðDx; DyÞ for convenience. The A semi-discrete finite volume scheme for the numerical approximation of Eq. (94) then takes the following form where the U i;j denote the approximate cell averages of the conserved variables, the F iAE1=2;j and G i;jAE1=2 are approximate facial averages of the fluxes through the cell boundary in respective direction, gðuðx; y jAE1=2 ; tÞÞ dx; ð97Þ and the S i;j are approximate cell averages of the source term The next paragraphs compactly describe the components of a generic finite volume scheme in several space dimensions. For the sake of presentation, we suppress the temporal dependence of the quantities. Reconstruction The first task is to reconstruct an accurate subcell profile from the cell-averaged conserved variables. We denote such an r-th order accurate piecewise polynomial reconstruction procedure by where S i;j denotes the stencil of the reconstruction for cell X i;j . Many such reconstruction procedures have been developed in the literature and we refer to the references previously mentioned in Sect. 2.2.1. For example, the stencil for a spatially first-order accurate piecewise constant consists only of the cell itself For a spatially second-order accurate piecewise linear reconstruction, the stencil includes the four adjacent cells Numerical fluxes The numerical fluxes through the cell faces are obtained by numerical integration of one-dimensional numerical fluxes. The facially averaged numerical fluxes in x-direction are given by where Q iþ1=2;j denotes a q-th order accurate quadrature rule with N q nodes y j;b 2 X j and weights x b , and F is a one-dimensional numerical flux formula in x-direction (see Sect. 2.2.2). Likewise, the facially averaged numerical fluxes in y-direction are given by where Q i;jþ1=2 denotes a q-th order accurate quadrature rule with N q nodes x i;a 2 X i and weights x a , and G is a one-dimensional numerical flux formula in y-direction.
In general, the numerical flux formulas in respective direction are often obtained by appropriate rotation of a one-dimensional flux formula in x-direction (see, e.g., Toro 2009 for details).
Numerical source terms We shall consider only unsplit methods for the numerical integration of the source terms. The cell-averaged numerical source terms are then given by where Q i;j denotes a q-th order accurate quadrature rule with N q Â N q nodes x i;a 2 X i , y j;b 2 X j and weights x a , x b in x-and y-direction, respectively. If we assume that the point values of the source term can be evaluated with spatial order of accuracy s, then this gives a spatially minðq; sÞ-th order discretization of the source term (provided enough smoothness). In practice, the same quadrature rules are often used in the numerical flux integration along the x-and y-direction. A tensor product quadrature rule is typically used for the numerical source term integration. For first-and second-order accuracy, the midpoint rule is the quadrature rule of choice. Beyond second-order accuracy, N q -point Gauss-Legendre or Gauss-Lobatto quadrature rules are usually used (N q [ 1).

Time discretization
The semi-discrete evolution Eq. (95) can be integrated with s-th order in time as in the one-dimensional case (see Sect. 2.2.4).
This concludes the brief description of a standard minðq; r; s; sÞ-order accurate finite volume scheme in several space dimensions (for smooth enough solutions). We refer again to the excellent textbooks in the literature for precise derivations and generalizations (curvilinear coordinates, unstructured meshes, etc.).

Steady states
Balance laws in several space dimensions too admit non-trivial steady states. As in the one-dimensional case, the numerical approximation of solutions near a steady state characterized by a delicate balance is generally challenging for standard finite volume schemes. The underlying reason is again twofold. First, standard reconstruction procedures are not well suited to represent steady-state solutions and, second, the source term discretization is performed independently from the discrete flux divergence. Both conspire that steady states are only preserved up to truncation errors. Hence, the numerical resolution needs to be high enough over the entire simulation duration such that the physical phenomena of interest are not affected by these truncation errors. The required resolution in several dimensions may then quickly lead to prohibitively high computational costs.
Although truncation errors are at the very essence of numerical approximation, it is again highly desirable to design schemes that preserve exactly a discrete form of the steady states; even more so than in one dimension. The one-dimensional recipe based on piecewise steady reconstruction naturally extends to designing wellbalanced schemes in multiple dimensions. The principle is again to decompose the solution into equilibrium and (not necessarily small) perturbation parts uðx; yÞ ¼ u eq ðx; yÞ þ duðx; yÞ; where the equilibrium part u eq ðx; yÞ fulfills the steady state balance As before, the piecewise steady reconstruction requires the computability of such multi-dimensional steady states. In general, this is an even more difficult undertaking than in one dimension. However, we again tacitly assume that the differential equation Eq. (104) can be solved for certain steady states of interest, either exactly ( ¼ 1) or approximately (\1). We reiterate that this is the main challenge when designing a well-balanced finite volume method for a particular system of balance laws. Some concrete examples are given in Sect. 3 for the Euler equations. Taking the above assumption for granted, the construction of a wellbalanced finite volume scheme follows the same recipe as in one dimension. The following subsections describe the ingredients in detail.
Before we proceed, let us stress once more that the following developments crucially hinge on the (exact or approximate) solvability for the multi-dimensional steady states of interest. We will see some examples for the Euler equations, where this can be achieved for barotropic fluids. In general, however, this is far from obvious. Let us mention here the truly two-dimensional well-balanced schemes developed by Bianchini and Gosse (2018), Gosse and Vauchelet (2020), Gosse (2021) for several balance laws. The latter construct fully two-dimensional steady states profiles by solving elliptic boundary-value problems and may serve as a guidance for the type of equations of interest in computational astrophysics.

Piecewise steady reconstruction WR
The piecewise steady reconstruction WR of a subcell profile U i;j ðxÞ within each cell X i;j from the associated cell averages U i;j is given by the decomposition where U eq;i;j ðx; yÞ and dU i;j ðx; yÞ denote the local equilibrium and perturbation part in cell X i;j , respectively. Next, we present the construction of each part. Within each cell X i;j , the local equilibrium reconstruction U eq;i;j ðx; yÞ is found by fitting an equilibrium solution U eq ðx; yÞ among the steady states of interest to the cell average U i;j . Because the given cell average U i;j may be arbitrarily far from an equilibrium, we first need to project U i;j onto a cell average U eq;i;j that is consistent with the steady states of interest. The local equilibrium reconstruction U eq;i;j ðx; yÞ in cell X i;j is then determined by matching an equilibrium profile Eq. (105) with the cell average of the equilibrium projected conserved variables U eq;i;j where Q i;j denotes a q-th order accurate quadrature rule over cell X i;j . Within every cell, we now have a minð; qÞ-th order accurate local equilibrium profile (if exact integration is chosen, then q ¼ 1). As in the one-dimensional case, the difficulty of this step depends strongly on the considered balance law and associated steady states of interest. Moreover, the possibility that no equilibrium profile is found needs to be considered too. In that case, it is again simply set to zero U eq;i;j ¼ 0. The local equilibrium perturbation dU i;j ðx; yÞ is obtained by extrapolating the local equilibrium profile U eq;i;j from cell X i;j to the neighboring cells, computing the cell average of this extrapolated equilibrium profile, and applying a standard piecewise polynomial reconstruction to the cell-averaged equilibrium perturbations: Intuitively, this senses how far away the states in the neighboring cells are from the equilibrium U eq;i;j in cell X i;j . It is clear that the piecewise steady reconstruction preserves the equilibrium by construction, since the equilibrium perturbation vanishes dU i;j 0, and it is minð; q; rÞ-th order accurate at and away from equilibrium (for sufficiently smooth solutions). Likewise, it is obvious that the piecewise steady reconstruction reduces to a standard reconstruction Eq. (99) if U eq;i;j 0 (e.g., if no local equilibrium is found in Eq. (107)). The alternative formulations of the piecewise steady reconstruction Eqs. (49) and (51) can also be extended to several space dimensions in a straightforward manner. However, we shall leave this to the interested reader.

Well-balanced source term discretization WS
As in the one-dimensional setting, a direct discretization of the source terms as in Eq. (102) will in general not lead to a well-balanced scheme. In lieu thereof, we rewrite the source terms as by trivially adding and subtracting the local equilibrium component and replacing the latter by the equilibrium flux divergence. The well-balanced source term discretization is then obtained by numerical quadrature as follows S i;j ¼ À 1 Dx Q iþ1=2;j ðf ðU eq;i;j ÞÞ À Q iÀ1=2;j ðf ðU eq;i;j ÞÞ À Á À 1 Dy Q i;jþ1=2 ðgðU eq;i;j ÞÞ À Q i;jÀ1=2 ðgðU eq;i;j ÞÞ À Á þ Q i;j ðsðU i;j Þ À sðU eq;i;j ÞÞ; where the divergence theorem is applied to the equilibrium fluxes by numerical integration over the cell boundary. This is clearly a minð; q; r; sÞ-th order accurate discretization of the source term. Moreover, note that the above well-balanced source term discretization seamlessly reduces to a standard discretization if no local equilibrium is found in the piecewise steady reconstruction (i.e., U eq;i;j 0). The alternative well-balanced source term discretization Eq. (60) and its highorder extension via Richardson extrapolation Eq. (62) can also be generalized to multiple space dimensions. We again leave this to the interested reader.

Assembling a well-balanced finite volume scheme
A well-balanced finite volume scheme Eq. (95) for the two-dimensional conservation law Eq. (94) is now readily assembled with the previously described ingredients: (1) A minð; q; rÞ-th order accurate piecewise steady reconstruction WR (Eq. (106)).
The resulting scheme is a min ð; q; r; s; sÞ-th order accurate well-balanced finite volume scheme for balance laws in two dimensions (for smooth enough solutions). The scheme preserves exactly a min ð; qÞ-th order accurate discrete form of the steady states of interest (up to machine precision). The proof follows the same steps as in the one-dimensional case and is not repeated here.
Although we have focussed on the structured two-dimensional Cartesian case, well-balanced schemes can be designed for three-dimensional Cartesian and curvilinear meshes following the same recipe. The same applies to unstructured meshes (see, e.g., Grosheintz-Laval 2021).

Well-balanced finite difference schemes
In this subsection, we focus on the construction of well-balanced finite difference schemes. Before we begin, we concisely summarize a generic finite difference scheme for one-dimensional balance laws. We refer to the excellent available textbooks and review articles for more comprehensive presentations, e.g., Laney (1998), Shu (1998Shu ( , 2009. Equipped with the basic principles and accompanying notation, we then present a (non-exhaustive) selection of frameworks for wellbalancing finite difference schemes. Further frameworks are developed, for instance, by Bermudez and Vazquez (1994), Gascón and Corberán (2001) (2021) and we refer the interested reader to these references.

Conservative finite difference schemes
A semi-discrete conservative finite difference method approximates the differential form of the balance law Eq. (4) directly by where the primal unknowns are the point values U i of the conserved variables at cell centers x i . They approximate the point value of the exact solution uðx; tÞ at cell centers: The F iAE1=2 and S i denote the numerical fluxes through the cell interfaces and the numerical source term at cell centers, respectively. We briefly describe these components in the following paragraphs. Since we consider semi-discrete schemes, we will drop the time dependence for ease of presentation.

Flux reconstruction
The numerical fluxes at cell interfaces F iAE1=2 are assumed to be consistent with the physical flux function f and Lipschitz continuous. The flux difference in Eq. (111) is constructed to be an (r-th order) accurate approximation of the flux divergence in Eq. (4) at cell centers: The numerical fluxes F iAE1=2 approximate the cell interface values hðx iAE1=2 Þ of a certain function hðxÞ implicitly defined by The above definition implies that the cell averages of this function hðxÞ are given by So the approximation of hðxÞ boils down to a reconstruction problem. The same reconstruction procedures R as in finite volume methods can be used. Hence, an r-th order accurate approximation of hðx iAE1=2 Þ is readily computed: One can show that the discrete flux divergence Eq. (113) is then an r-th order accurate approximation of the flux divergence at cell centers. Note that this may seem surprising as Eq. (113) looks like a simple second-order centered finite difference approximation (ensuring the crucial conservative character of the scheme). Unfortunately, beyond second-order accuracy (r [ 2), this accuracy claim is only valid for uniform or smooth grids (see, e.g., Merriman 2003;Shu 2009). For the subsequent developments, we have to slightly refine the so far introduced notation for the reconstruction procedures in Sect. 2.2.1. The reconstruction of some quantity Q within cell X i from cell averages fQ k g over the cell's stencil S i can be written as where the c i;l denote the reconstruction coefficients. In general, these depend on the evaluation location x 2 X i and on certain so-called smoothness measures of the reconstruction data fQ k g k2S i to guarantee non-oscillatory behavior. For smooth data, these coefficients tend towards optimal constants that maximize the accuracy of the reconstruction. Instead of using the data of the quantity to reconstruct, the smoothness measures can be based on the cell averages of some other quantity 11 fW k g 11 To guarantee accurate and non-oscillatory reconstruction properties, the smoothness of quantity Q should imply smoothness of the other quantity W and reciprocally.
where the second argument of the reconstruction procedure is now a set of couples fðQ k ; W k Þg k2S i . The first element is the cell average of the quantity to reconstruct and the second element is the cell average of the quantity on which the smoothness measures are based. For finite difference schemes, a straightforward component-wise flux reconstruction for (systems of) balance laws may also lead to spurious oscillations, especially at high orders of accuracy and when strong flow discontinuities interact. Then a reconstruction in local characteristic variables as described at the end of Sect. 2.2.1 can be advantageous (see, e.g., Jiang and Shu 1996;Balsara and Shu 2000;Qiu and Shu 2002).
Flux splitting For numerical stability, upwinding is introduced by splitting the flux into right(þ) and left(-) propagating contributions The characteristic values of of þ ou are all non-negative and the characteristic values of of À ou are all non-positive, i.e., of þ ou ! 0 and of À ou A popular flux splitting is the (local) Lax-Friedrichs flux splitting The fluxes can then be discretized according to the wave direction as The right going flux contribution F þ iþ1=2À is biased one cell to the left from the cell interface x iþ1=2 by using the reconstruction stencil S i . Accordingly, the left going flux contribution F À iþ1=2þ is biased one cell to the right from the cell interface x iþ1=2 by using the reconstruction stencil S iþ1 . The complete upwinded numerical flux is then simply the sum of both contributions The maximum in Eq. (122) can either be taken globally over the whole grid fU i g or locally, for example, over the flux reconstruction stencils S i [ S iþ1 . Note that a can be interpreted as a parameter introducing artificial viscosity into the scheme, ensuring its numerical stability (see, e.g., Laney 1998).
Numerical source term The numerical source term S i is a point-value evaluated at cell centers such that it is an (s-th order) accurate approximation: Time integration The semi-discrete evolution equations Eq. (111) for the point values U i can be approximately integrated in time by s-th order accurate ODE solvers. Also for high-order finite difference schemes, the SSP Runge-Kutta methods are often used in practice. One popular example is the temporally thirdorder accurate SSP Runge-Kutta (SSP-RK3) method where L denotes the spatial discretization operator in Eq. (111). This completes the short description of a standard minðr; s; sÞ-th order accurate finite difference scheme. The extension to several space dimensions is achieved by discretizing the flux divergence direction by direction. A significant computational advantage of finite difference schemes over finite volume schemes is that the multidimensional reconstruction procedures and quadrature rules are avoided entirely. The price of this advantage is the restriction to uniform or smooth curvilinear grids. We refer to the references given earlier for further details and developments. Next, we have a closer look at what happens near steady states.

Steady states
Standard finite difference schemes alike have troubles in preserving steady-state solutions in general. The fundamental cause of these troubles is analogous to finite volume methods: the source term discretization is independent of the discrete flux divergence. The consequence, in turn, is that steady states are generally preserved only up to truncation errors, requiring the numerical resolution to be high enough that they do not obscure the phenomena of interest.
The objective of well-balanced finite difference schemes is to preserve a discrete form of certain steady-state solutions. To do this, it is tacitly assumed, as for finite volume schemes, that the differential equation Eq. (44) can be solved (exactly or approximately) for the steady states of interest Eq. (45). The following subsections introduce a selection of frameworks for building well-balanced finite difference schemes for one-dimensional balance laws. Multi-dimensional balance laws can be treated by straightforwardly extending the one-dimensional building blocks.

Well-balanced finite difference schemes based on source term decomposition
A framework for the construction of high-order well-balanced finite difference schemes for a class of balance laws was developed by Xing and Shu (2005, 2006a, Li and Xing (2018b). The framework can also be applied to finite volume and discontinuous Galerkin methods (see, e.g., Xing and Shu 2006b;Noelle et al. 2009).
The framework assumes that each source term component can be (analytically) decomposed as follows where s m denotes the m-th component of the source term sðu; xÞ. 12 The functions b m ¼ b m ðxÞ are supposed to depend only on space and the a m ¼ a m m; x ð Þ depend additionally on so-called equilibrium variables denoted by m. The latter are function of the conserved variables and space, i.e., m ¼ mðu; xÞ, and they are required to become constant at the steady states of interest U eq mðU eq ; xÞ ¼ const.; explaining the name equilibrium variables. Similarly, the functions a m are assumed to become spatially constant for U eq : a m ðmðU eq ; xÞ; xÞ ¼ const.
As a result, such a source term decomposition clearly has the following property at a steady state of interest U eq : o ox f m ðU eq Þ ¼ s m ðU eq ; xÞ which obviously implies that f m ðU eq Þ À a m mðU eq ; xÞ; x À Á b m ðxÞ ¼ const.
Although these assumptions may seem rather restrictive, a broad class of balance laws and associated steady states possess such a source term decomposition. We will provide some illustrating examples below together with a general recipe.
The key idea is now to mimic the above property at the discrete level to ensure the exact balance at steady states. We begin by discussing the case without flux splitting. Then, this is achieved by approximating the derivative of the function b m with the same finite difference operator as used for the approximation of the flux divergence Eq. (113) Note that the above reconstruction uses the m-th component of the flux f m as smoothness measures. This ensures the consistency of the difference operators applied to the flux and source term by construction. 13 An r-th order accurate discretization of the m-th source term component is thus given by and the complete source term vector is simply S i ¼ ½Á Á Á ; S m;i ; Á Á Á T . Before we discuss the case with flux splitting, we verify that the scheme Eq. (111) and the above source term discretization Eq. (133) is wellbalanced. Suppose that the conserved variables point values U i are initialized with an equilibrium U eq ðxÞ fulfilling all the necessary assumptions, that is, U i ¼ U eq ðx i Þ. Then Eq. (129) gives immediately at all cell centers x i a m ðmðU i ; x i Þ; x i Þ ¼ a m ¼ const. and one verifies that Eq. (130) holds discretely as follows In the last equality, the consistency of the reconstruction was used (constant functions are reconstructed exactly). Now the semi-discrete update equation Eq. (111) gives for the m-th solution component U m;i Since this is valid for all the solution components, we conclude that the scheme is well-balanced as claimed. Next, we discuss the case with flux splitting Eq. (119) and Eq. (123). Then, the source term is split into a right (þ) and a left (-) contribution as follows Note that the above reconstruction uses the m-th component of the left(þ)/right (-) propagating flux contributions f AE m as smoothness measures. This ensures the consistency of the difference operators applied to the flux and source term by construction.
To guarantee an exact balance, the flux splitting Eq. (121) also needs to be slightly modified. One possibility suggested by Xing and Shu (2006a) Due to Eq. (129), this ensures that at steady state the artificial viscosity tends towards zero. Importantly, this does not interfere with the scheme's design accuracy (see Xing and Shu 2006a). However, the constant e a m should be suitably adjusted to maintain enough artificial viscosity for numerical stability. This depends on the concrete application and may require some fine-tuning. The well-balancing property with flux splitting can now be shown easily as before. This ends the description of the well-balanced finite difference schemes based on source term decomposition. To conclude on this family of schemes, let us give two illustrating examples of balance laws with their corresponding equilibrium variables and source term decomposition.
Example: Linear advection-reaction equation Consider (again) the linear advection-reaction equation which still has the steady-state solutions A particular choice for the equilibrium variable of Eq. (137) is which becomes clearly a constant for the steady states U eq ðxÞ. A source term decomposition of the form Eq. (127) is then given by with a 1 ðm; xÞ ¼ e m and b 1 ðxÞ ¼ ae À k a x : Example: Euler equations As a slightly more complex example, consider the onedimensional Euler equations describing the motion of fluids subject to a gravitational field with potential /. They are given by the conservation of mass, momentum and fluid energy where the vector of conserved variables, fluxes and source terms are Here, q is the mass density, v x the velocity and E ¼ qe þ 1 2 qv 2 x the total fluid energy density. The pressure p is related to the density q and specific internal energy e by an equation of state p ¼ pðq; eÞ. For steady adiabatic flow, the equilibrium variables are given by (see, e.g., Landau and Lifshitz 1987) where h and s denote the fluid's specific enthalpy and entropy, respectively. Let U eq ðxÞ denote the conserved variables corresponding to some given equilibrium variables m eq , i.e., U eq ðxÞ ¼ uðm eq ; xÞ is the inverse transformation of the equilibrium variables to conserved variables (see Grosheintz-Laval and Käppeli 2020 for a way to compute it). Also, let p eq ðxÞ denote the corresponding equilibrium pressure. A source term decomposition of the form Eq. (127) is then given by the following. The first component is trivial where m eq are some fixed equilibrium variables. However, we shall not further discuss the above generic decomposition and refer the interested reader to the original research articles cited at the beginning of this section for the details.

Well-balanced finite difference schemes based on piecewise steady flux reconstruction
Parés and Parés-Pulido (2021) recently proposed a framework for designing highorder well-balanced finite difference schemes for balance laws. It is elegantly based on a piecewise steady flux reconstruction and, as such, is closely related to the finite volume recipe based on piecewise steady reconstruction presented in Sect. 2.4. Unfortunately, the formulation is generally not fully conservative in that if a system of balance laws possesses a conservative subsystem, 14 the method may not reduce to a conservative finite difference method for that subsystem of conservation laws. However, the conservation errors are of the order of the truncation errors and vanish in the infinite resolution limit. Since high-order finite difference schemes are based on flux reconstructions, the key idea is to decompose the flux into equilibrium and (not necessarily small) perturbation parts f ðuðx; tÞÞ ¼ f ðu eq ðxÞÞ þ f ðuðx; tÞÞ À f ðu eq ðxÞÞ À Á where the equilibrium part f ðu eq ðxÞÞ fulfills by definition the steady-state balance Eq. (44). Like in the piecewise steady reconstruction used in finite volume schemes, this requires the ability to compute such steady states of interest Eq. (45). Parés and Parés-Pulido (2021) and suggest the following semi-discrete finite difference approximation Next, we describe the Parés and Parés-Pulido (2021) well-balanced finite difference scheme in detail. This will especially clarify the double index notation for the numerical equilibrium perturbation fluxes dF i;iAE1=2 . One begins with the determination of a local equilibrium profile U eq;i ðxÞ in each cell X i by fitting an equilibrium solution U eq ðxÞ among the steady states of interest to the point value U i . Because U i may be far from a steady state, it is first projected onto a point value U eq;i that is consistent with the steady states of interest. Therewith, the local equilibrium profile U eq;i ðxÞ in cell X i is determined by matching an equilibrium profile U eq ðxÞ to the equilibrium projected point value U eq;i o ox f ðU eq;i ðxÞÞ ¼ sðU eq;i ðxÞÞ ð155Þ such that the profile is anchored at the cell center point value U eq;i Note that this is the main difference with the piecewise steady reconstruction, where the local equilibrium profile is matched with the cell average. The equilibrium perturbation flux divergence df is discretized to r-th order as where the numerical equilibrium perturbation fluxes dF i;iAE1=2 approximate the cell interface values dh i ðx iAE1=2 Þ of a certain function dh i ðxÞ implicitly defined by Following the same principles as in a standard finite difference scheme, the above definition implies that the cell averages of this function dh i ðxÞ are given by Now the double index notation becomes clear. The first index ''i'' indicates the cell where the equilibrium reconstruction U eq;i ðxÞ is anchored. The second index ''k'' indicates the cell center x k where the equilibrium perturbation flux is evaluated. Also notice that the equilibrium reconstruction U eq;i ðxÞ might only be an -th order approximation of a true equilibrium u eq ðxÞ, which means that also the cell averages dh i;k are only accurate to this precision. An r-th order reconstruction procedure is then applied to the cell averages dh i;k to obtain the numerical equilibrium perturbation fluxes at cell interfaces For numerical stability, these numerical fluxes need to be upwinded and this will be explained below. Notice that two reconstructions, dF i;iþ1=2 and dF iþ1;iþ1=2 , have to be computed at every cell interface. In contrast, only one is needed in a standard finite difference scheme. The numerical perturbation source term dS i is discretized as Note that dS i does not vanish in general because the solution could be far away from a steady state. Before we discuss the upwinding of the numerical equilibrium perturbation fluxes, let us verify that the scheme Eq. (154) is indeed well-balanced. Suppose that we initialize the conserved variables point values U i with a steady state of interest U eq ðxÞ that is, we initialize with well-prepared initial data. The equilibrium projection U eq;i ¼ U i is trivial. Assuming that for all cells X i the local equilibrium reconstruction U eq;i ðxÞ matches with U eq ðxÞ at cell centers (i.e., U eq;i ðx i Þ ¼ U eq ðx i Þ), then the equilibrium perturbation flux vanishes at all cell centers (i.e., f ðU k Þ À f ðU eq;i ðx k ÞÞ ¼ 0). Any consistent reconstruction procedure will then produce vanishing numerical equilibrium perturbation fluxes at cell interfaces dF i;iAE1=2 for all cells X i . Along the same lines, the numerical source term perturbation Eq. (161) vanishes identically. Therefore, the semi-discrete finite difference scheme Eq. (154) gives and the scheme is well-balanced as advertised.
For numerical stability, the numerical equilibrium perturbation fluxes have to be properly upwinded. As in a standard finite difference scheme, this is achieved by flux splitting: The complete upwind numerical equilibrium perturbation fluxes are then For example, the (local) Lax-Friedrichs flux splitting Eq. (121) could be used. It is straightforward to show that flux splitting does not interfere with the well-balancedness of the scheme. This concludes the presentation of the high-order well-balanced finite difference schemes of Parés and Parés-Pulido (2021). We want to stress the fact that the numerical equilibrium perturbation fluxes at the cell interfaces are not equal in general: Of course, this is not unexpected for balance laws and their non-conservative character due to their source terms. However, the latter may also be the case for the conservative equations in a system of balance laws. Fortunately, one can show that the conservation errors made by the schemes are of the order of the truncation error and, therefore, tend to vanish with increasing resolution. Moreover, this non-conservation can be avoided for certain balance laws using some fine-tuning of the artificial viscosity terms in the numerical fluxes. We refer to Parés and Parés-Pulido (2021) for the intricate details.

Well-balanced methods for the Euler equations
In this section, we showcase the frameworks introduced in Sect. 2. We consider the Euler equations of compressible hydrodynamics as a prototypical example for a system of balance laws. The section starts by introducing the equations and their steady states. This is followed by a classification of well-balanced schemes for the Euler equations and the presentation of some representative example for each class. The section concludes by presentating several numerical test problems highlithing the performance of the well-balanced schemes.

The Euler equations
The Euler equations describe the motion of ideal fluids subject to (gravitational) forces in which thermal conductivity and viscosity are unimportant. They express the conservation of fluid mass, momentum and energy: Here q is the fluid mass density, v the velocity and p the pressure. The total fluid energy E ¼ qe þ q 2 v 2 is composed of internal and kinetic energy densities. The pressure p is related to the density q and the specific internal energy e through an equation of state (EoS) p ¼ pðq; eÞ. The latter determines all thermodynamic quantities by specifying any two values describing the state. 15 The source terms on the right-hand side of Eqs. (166) model the effect of gravitational forces on the fluid. They are characterized by a gravitational potential, which may either be a given function of the coordinates, / ¼ /ðxÞ, or prescribed by the fluid's self-gravity depending on the concrete application. In the latter case, the potential is determined by the Poisson equation, where G is the gravitational constant. By including the gravitational interaction explicitly into the conservation balance, as opposed to a generic external force, the Euler equations can be written in alternative forms that emphasize their total conservative character. In the case of self-gravitating flows, the momentum source term can be reformulated as the divergence of the gravitational stress tensor, allowing its inclusion into the momentum flux tensor (see, e.g., Shu 1992). Unfortunately, when gravity is described by an external static gravitational potential, the source term in the momentum equations cannot be eliminated. However, the energy Eq. (166c) can be reformulated as in both cases. Here E T ¼ E þ E G is the total (fluid and gravitational) energy density and F G represents a ''gravitational energy flux''. For a static gravitational field, the gravitational energy density is E G ¼ q/. If the gravitational potential is prescribed by the fluid's self-gravity, the gravitational energy density is E G ¼ q/=2, where the factor 1/2 avoids the double-counting of pairs of fluid elements. For time-independent gravitational fields, as is the case at steady states which are our primal concern, this energy flux component becomes F G ¼ qv/. Hence, the energy equation can be written in conservation form Below, we will also use a slightly different form that evolves directly the fluid energy E with the gravity source term expressed as follows For static gravitational fields, the source term can then be discretized such that the total energy E T is preserved in a straightforward manner. We refer to the developments in Springel (2010), Hanawa (2019), Katz et al. (2016), Mullen et al. (2021) for an in-depth discussion and, in particular, the extension to self-gravitating flows. Before we proceed with a discussion of some interesting steady states of the Euler equations, let us rewrite them in canonical balance law form. For simplicity, we restrict the presentation to a one-dimensional setting. The Euler Eqs. (166) then compactly read where the conserved variables, fluxes, and gravitational source term are given by Moreover, the primitive variables are denoted by w ¼ ½q; v x ; p T .

Steady states of the Euler equations
The Euler Eqs. (166) allow for a myriad of non-trivial steady states. Many interesting astrophysical applications occur near or involve such equilibrium flows.
A particular example is hydrostatic equilibrium, which describes the mechanical balance between pressure and gravity forces. In the case of self-gravity, the hydrostatic equation can be combined with the Poisson Eq. (167) to yield (Landau and Lifshitz 1987) It must be emphasized that the hydrostatic equilibrium equation only describes a mechanical equilibrium. A certain thermal stratification has to be supplemented to integrate the equations. Another steady-state example of the Euler Eqs. (166) is provided by steady adiabatic flow. Such steady flows are governed by Bernoulli's equation (Landau and Lifshitz 1987) where h is the specific enthalpy. The above relation holds along each streamline, i.e., lines tangent to the velocity of the flow. In general, the constant may take different values for different streamlines. Flows near such steady states are ubiquitous in nature as in accretion and wind phenomena. However, hydrostatic equilibrium and steady adiabatic flow are just two particular examples. In many astrophysical applications, the flow of interest takes place near more complex equilibrium configurations. One example is rotational hydrostatic equilibrium which describes the balance between pressure, gravitational and centrifugal forces. Under certain conditions, such equilibrium configurations fulfill Eq. (173) with the gravitational potential / replaced by an effective potential U including effects of rotation whereis the cylindrical radius and X (the constant) angular velocity (see, e.g., Tassoul 1978). The Euler equations are an idealized model for the description of flows. They lack any dissipative processes such as thermal conduction and viscous stresses. Including these effects leads to the Navier-Stokes equations (see, e.g., Landau and Lifshitz 1987). In astrophysical flows, a substantial amount of the energy density, momentum density and stress is in the form of radiation (e.g., photons, neutrinos). Such radiating flows are described by the equation of radiation hydrodynamics (see, e.g., Mihalas and Weibel-Mihalas 1984;Castor 2004). Similarly, magnetic fields play a prominent role and a fluid model is provided by the equations of magnetohydrodynamics (MHD). These more sophisticated physical models and resulting equations possess even richer classes of non-trivial steady states. For example, radiative hydrostatic equilibrium or two-dimensional plasma steady states as described by the Grad-Shafranov equation would be of interest, to name a couple.
A (general) relativistic flow description would also be of interest. We hope that the following presentation for the Euler equations can serve as a useful guideline to develop well-balanced schemes for these extended physical models and their more intricate steady states.

Well-balanced methods for the Euler equations
Many well-balanced schemes for the Euler equations have been developed in the literature. Although the schemes follow different design philosophies, they may be broadly classified based on the steady states of interest they preserve:

A priori known steady states
The steady state of interest is assumed to be globally known.

Barotropic steady states
The steady states of interest assume a certain barotropic relation, effectively imposing a thermal stratification of the equilibrium state.

Discrete steady states
The steady states of interest are built from a consistent discretization of the defining PDE.
The order reflects the successive weakening of the made assumptions about the steady states of interest. Below we classify many of the well-balanced schemes for the Euler equations that are available in the literature into one of these categories. A representative is also presented for each category. This classification and representative selection is a more or less arbitrary choice of the author and as such is not perfect. For instance, some schemes may be classified in several categories, e.g., any barotropic or discrete steady states type scheme can be reduced to an a priori type scheme by simply freezing the piecewise steady reconstruction. Indeed, a taxonomy that fits it all may not even exist. Nevertheless, we consider it useful for a first, albeit rough, classification.

A priori known steady states
The first category of well-balanced schemes for the Euler equations that we will discuss assumes that the steady state of interest is globally known. This allows the construction of well-balanced finite volume, finite difference and discontinuous Galerkin schemes that can preserve any known steady state, making these methods extremely versatile. The only caveat is, of course, that the steady state has to be known a priori. However, this is not such a severe restriction as it may seem as long as the phenomena of interest don't deviate too much from the fixed steady state. Several well-balanced schemes of this category have been developed in the literature: • Constantinescu (2015, 2016) developed high-order well-balanced finite difference schemes. The schemes are based on the source term decomposition method of Xing and Shu (2013) which is extended to types of equilibria encountered in atmospheric flow simulations. This includes isentropic atmospheres and stratified atmospheres with specified Brunt-Väisälä frequency. • Li and Xing (2016a) developed high-order well-balanced finite volume schemes for isothermal and polytropic hydrostatic equilibrium on the basis of the source term decomposition method for finite volume schemes of Shu (2006b, 2013). Li and Xing (2018b) simplify the original procedure of Xing and Shu (2013) by saving some costly WENO reconstructions in the discrete source term evaluation. This is possible because the discrete source term can be evaluated once at the beginning of the simulation since the equilibrium state is known initially. Furthermore, they extend the formalism to more general steady states including isothermal and polytropic hydrostatic equilibrium. Robust highorder discontinuous Galerkin schemes are developed by Wu and Xing (2021) following the techniques introduced by Li and Xing (2016b). These schemes are capable of balancing any known hydrostatic equilibrium, have a guaranteed positivity-preserving property for general equation of states, and can handle unstructured meshes. • Touma et al. (2016) elaborated a well-balanced second-order unstaggered central finite volume scheme that is able to preserve isothermal equilibrium. The necessary back-and-forth projections between the unstaggered and staggered cells are made equilibrium-preserving by using a variant of the surface gradient method of Zhou et al. (2001). The method provides discrete solutions on a single grid by using a ''ghost'' staggered grid that is only used in the discrete evolution steps (i.e., it is effectively an unstaggered central scheme). volume schemes on moving nonconforming meshes able to preserve rotational hydrostatic equilibrium. This is accomplished within the well-balanced pathconservative framework of  where the source terms are treated as non-conservative products combined with a well-balanced reconstruction operator. • Veiga et al. (2019) developed discontinuous Galerkin schemes that can exactly balance any known equilibrium. Furthermore, they systematically compare the well-balanced schemes with unbalanced standard schemes regarding the computational cost for increasingly higher-order schemes. They observe that well-balanced schemes pay off, especially in multi-dimensional settings. • Berberich et al. (2019) developed second-order well-balanced finite volume schemes for arbitrary known hydrostatic equilibrium. The schemes can handle curvilinear grids and general equations of state (see also Berberich et al. 2018). They are based on known non-dimensionalized density a and pressure b functions on which a piecewise steady reconstruction is built and combined with a well-balanced source term discretization. Hence, they term their formalism the a-b well-balanced method. Klingenberg et al. (2019) generalize the a-b method to arbitrary orders of accuracy with CWENO reconstruction procedures and Richardson extrapolation for the source term discretization. The latter reconstruction procedures are particularly well suited for the source term discretization as they avoid any negative stencil weights within the cell by construction. Thomann et al. (2020Thomann et al. ( , 2019 combined the a-b method with relaxation Riemann solvers and implicit-explicit (IMEX) time integration tailored for the efficient treatment of low Mach-number flows. • Li and Gao (2021) devised a strategy to build high-order well-balanced finite difference schemes by introducing specialized nonlinear WENO differential operators which fulfill a certain homogenization condition that guarantees the exact balance between the flux and source terms. The latter are discretized with the source term decomposition method. • Berberich et al. (2021a) designed a general framework to construct high-order well-balanced finite volume schemes for any known solution of a given hyperbolic system. Their framework also preserves known time-dependent solutions, which may be interesting for certain applications such as uncertainty quantification. Kanbar et al. (2020) applied the formalism to unstaggered, second-order, central finite volume schemes. • Edelmann et al. (2021) propose an interesting comparison of several wellbalanced finite volume solvers for low Mach-number flows that are relevant for multi-dimensional stellar structure and evolution simulations. Moreover, they developed a multi-dimensional extension of the Cargo and LeRoux (1994) scheme. They emphasize that when combining well-balanced schemes with low Mach number solvers, care must be taken not to introduce spurious numerical artifacts. Low Mach-number solvers rely on special numerical flux functions that reduce the unphysical numerical diffusion in low Mach-number regimes. This can interfere negatively with the vanishing diffusion of well-balanced schemes at steady states.
We next sketch two representative of this category of well-balanced methods. We opted for the a-b well-balanced finite volume method for hydrostatic equilibrium of Berberich et al. (2018Berberich et al. ( , 2019 and collaborators. Subsequently, we also briefly describe another method that relies on a slightly different principle than most of the methods of the present category.

a-b well-balanced method
For simplicity, we consider a one-dimensional setting and limit the spatial accuracy to second-order. The a-b well-balanced method assumes that the hydrostatic equilibrium Eq. (173) to be preserved is explicitly known in terms of two dimensionless scalar functions aðxÞ and bðxÞ, q eq ðxÞ ¼ q 0 aðxÞ and p eq ðxÞ ¼ p 0 bðxÞ; fulfilling 1 q eq dp eq dx The constants q 0 and p 0 anchor the equilibrium density and pressure at some reference coordinate x 0 . The steady state of interest U eq ðxÞ in Eq. (45) is therefore explicitly known U eq ðxÞ ¼ which also highlights the importance of the functions aðxÞ and bðxÞ. We begin by the piecewise steady reconstruction in primitive variables. Berberich et al. (2018Berberich et al. ( , 2019 use the relative form Eq. (52). Since the steady state of interest is explicitly known, the local equilibrium reconstruction within each cell X i is trivial. Indeed, the equilibrium projection and matching steps simplify to a formality W eq;i ¼ q eq;i 0 p eq;i The cell-averaged equilibrium density q eq;i and pressure p eq;i or, equivalently, the a i and b i can be computed over the whole computational domain once at the beginning of the simulation and stored. The local relative equilibrium perturbation reconstruction Eq. (53) gives where the velocity and pressure are computed from the cell-averaged conserved variables The latter choice limits the spatial accuracy to formally second-order, regardless of whether a higher-order reconstruction procedure R is used. Note that the velocity component uses a standard reconstruction since it does not participate in the steady state of interest. If we further assume that the reconstruction procedure fulfills a certain scale invariance property 16 for any cell-averaged quantity Q k and constant C, we obtain the following final expression of the piecewise steady reconstruction It is obvious that if equilibrium cell-averages W eq;i are given to the above piecewise steady reconstruction, it reduces to the exact steady state Eq. (180) by construction. It remains to discuss the well-balanced source term discretization. Berberich et al. (2018Berberich et al. ( , 2019 use the alternative form Eq. (59). Since only the momentum source term is relevant at the steady state of interest, we directly obtain Moreover, Berberich et al. (2018Berberich et al. ( , 2019 use the conservative formulation of the total (fluid and gravitational) energy Eq. (169). The complete source term discretization thus reads It is straightforward to verify that the a-b method is able to preserve any known hydrostatic equilibrium Eq. (179) or Eq. (180). The extension to several space dimension is also uncomplicated (e.g., following the recipe in Sect. 2.6). Extending the method beyond second-order accuracy is slightly more tricky. The issue stems from the pressure reconstruction Eqs. (182) and (185), which estimates the pressure based on the cell-averaged conserved variables Eq. (183). We refer to Klingenberg et al. (2019) for the details.

Equilibrium truncation error annihilation method
Let us mention another approach that exists as ''folklore'' among practitioners. An exact balance of an a priori known steady state is simply obtained by subtracting the discretization error at steady state in each time step. At the analytical level and in the infinite resolution limit, this is tantamount to subtracting a zero from the equation (because the steady state fulfills the balance by definition, of course). In semi-discrete evolution form, this simply reads where L is the spatial discretization operator (see Eqs. (7) and (95) for finite volume) and U eq the explicitly computable cell averages of the known steady state U eq . Similarly, the same can be ported to other spatial discretizations such as finite difference and discontinuous Galerkin methods. For example, this method was used by Dedner et al. (2001) for the simulation of waves in stratified stellar atmospheres. The latter approach is probably easier and computationally cheaper to implement in an existing solver since the equilibrium discretization error LðU eq Þ can be calculated once and for all at the beginning of a simulation. However, this approach may interact in less predictable ways within the reconstruction steps than the previously discussed methods such as the a-b method.

Barotropic steady states
The second category of well-balanced schemes for the Euler equations is designed to preserve barotropic steady states. In barotropic fluids, the density is a function of pressure only. This establishes a thermal equilibrium stratification and that information is directly exploited by the well-balanced schemes in this category. A prominent example is provided by isentropic conditions in which the specific entropy is constant. Consider the fundamental thermodynamic relation where h is the specific enthalpy, T the temperature and s the specific entropy. Hydrostatic equilibrium Eq. (173) under isentropic conditions (ds ¼ 0) then gives which can be trivially integrated to Note that this is an explicit expression for hydrostatic equilibrium under isentropic conditions (assuming the gravitational potential is known). Once the constant fixed at some reference coordinate, the hydrostatic stratification is fully determined. Bernoulli's Eq. (175) is the generalization to steady adiabatic flow. Along the same lines, isothermal hydrostatic equilibrium is derived from the fundamental thermodynamic relation where g is the specific Gibbs free energy, yielding More generally, the expression is integrable for barotropic fluids and hydrostatic equilibrium takes the form Here H ¼ Hð#; pÞ is a thermodynamic potential depending on the natural variables # and p. For example, the general expression encompasses the isentropic (H ¼ h, The above explicit expressions for barotropic hydrostatic equilibrium can then be used for the design of well-balanced schemes. Many schemes of this category have been developed in the literature: • Botta et al. (2004) designed a well-balanced second-order finite volume scheme for isentropic hydrostatic equilibrium as frequently encountered in numerical weather prediction and climate modeling. The scheme is based on a local piecewise steady reconstruction of isentropic hydrostatic states within each grid cell that are adapted to the local thermodynamic conditions at each time step. • Fuchs et al. (2010a) constructed second-order finite volume schemes that preserve isothermal hydrostatic equilibrium. The schemes use a piecewise linear steady reconstruction of isothermal hydrostatic states explicitly exploiting the exponential stratification of density and pressure. Remarkably, the schemes are also well-balanced for certain magnetostatic equilibria such as used in the simulation of waves in stellar atmospheres (Rosenthal et al. 2002). • LeVeque (2010) developed a well-balanced second-order method for polytropic gas dynamics using the f-wave approach of Bale et al. (2002). The particular source term average at cell interfaces to guarantee the exact preservation of steady states is constructed using the theory of path conservative methods. Gundlach and LeVeque (2011) used the approach to study the universality in the run-up of shock waves to stellar surfaces. Ahmad and Lindeman (2007) have applied the f-wave approach to atmospheric flow problems and they report promising results. • Xu et al. (2010), Luo et al. (2011) constructed symplecticity-preserving gaskinetic schemes for compressible Euler and Navier-Stokes equations with gravity. The second-order schemes represent the gravitational potential as a piecewise constant with a potential jump at every cell interface. The schemes are designed such that an isothermal hydrostatic state is preserved during the process of particle transport and collision, which necessitates the use of an exact Maxwellian velocity distribution. See also the recent approach by Chen et al. (2020). • Xing and Shu (2013) designed high-order well-balanced finite difference schemes for isothermal hydrostatic equilibrium. The approach uses the source term decomposition method of Xing and Shu (2006a) (see Sect. 2.7.3) using the special form of isothermal hydrostatic states. Li and Xing (2016b) extend the formalism to discontinuous Galerkin schemes using the source term decomposition of Xing and Shu (2006b). • Käppeli and Mishra (2014), Käppeli (2017) developed second-order finite volume schemes that retain barotropic hydrostatic states exactly (up to machine precision). The schemes are based on the piecewise steady reconstruction of barotropic hydrostatic states using thermodynamic potentials and are capable of dealing with general EoS. Grosheintz-Laval and Käppeli (2019); Grosheintz-Laval (2021) extended the schemes to (spatially) arbitrary order of accuracy and unstructured meshes. Grosheintz-Laval and Käppeli (2020) generalized the second-order schemes to adiabatic steady flows.  constructed well-balanced second-order finite volume schemes for isothermal and polytropic hydrostatic equilibrium. The schemes rewrite the gravitational source terms in a specific form by exploiting the structure of the equilibrium state (similar to Xing and Shu 2013) and a piecewise steady reconstruction that uses equilibrium scaled variables. • Chandrashekar and Zenk (2017) designed high-order nodal discontinuous Galerkin methods for isothermal and polytropic hydrostatic equilibrium. The schemes use a form of the source term decomposition method as Xing and Shu (2013) combined with Gauss-Lobatto-Legendre quadrature rules. The latter choice ensures that at an equilibrium the solution is continuous between the cells and the flux and source discretizations exactly match. • Li and Xing (2018a) developed high-order (modal) discontinuous Galerkin schemes capable of preserving isothermal and polytropic hydrostatic states. The schemes are based on a generalized hydrostatic reconstruction with an equilibrium state recovery technique and a special projection operator guaranteeing the necessary continuity conditions of the numerical fluxes at the cell interfaces and a well-balanced source term discretization following Xing and Shu (2006c). • Gómez-Bueno et al. (2021a, 2021b) developed a general framework for constructing high-order well-balanced finite volume schemes for one-dimensional balance laws. The schemes are based on a piecewise steady reconstruction that constructs local equilibrium by solving the steady-state defining ODEs numerically. (In general, the schemes could therefore also be classified in the discrete steady state category.) In the context of the Euler equations, high-order well-balanced schemes for adiabatic sub-and supersonic steady flows are designed within the framework. By assuming adiabatic flow, they derive a system of ODEs for the steady states that is then solved numerically to derive the piecewise steady reconstruction.
In the following, we briefly outline a representative out of this category of wellbalanced schemes for barotropic hydrostatic equilibrium. We opt for the wellbalanced finite volume schemes of Käppeli and Mishra (2014), Käppeli (2017) and their higher-order extension by Grosheintz-Laval and . This (admittedly not entirely impartial) choice is motivated by the versatility of the approach as it seamlessly adapts to any barotropic relation and relies only on fundamental thermodynamic relations (i.e., it works for any EoS). We now follow the recipe from Sect. 2.4 and prepare the necessary ingredients to design a onedimensional well-balanced finite volume scheme for barotropic hydrostatic equilibrium. The multi-dimensional case is briefly addressed at the end of the description. The starting point is to ensure the computability of the steady states of interest. The barotropic hydrostatic equilibrium Eq. (195) immediately gives the following equilibrium profiles # eq ðxÞ ¼ # 0 and H eq ðxÞ ¼ H 0 þ / 0 À /ðxÞ: Here # 0 , H 0 and / 0 are the equilibrium's thermodynamic natural variable and potential, and the gravitational potential evaluated at some reference coordinate x 0 , respectively. The equilibrium conserved variables are then U eq ðxÞ ¼ where the equilibrium density q eq ðxÞ and internal energy density qe eq ðxÞ are computed through the EoS given the natural variable and thermodynamic potential. Note that we distinguish between the exact u eq ðxÞ and approximate U eq ðxÞ steady states of interest (see Eq. (45)). Although this may seem overly pedantic, it takes into account that the gravitational potential /ðxÞ may only be known approximately (say up to order OðDx Þ). The next ingredient is the piecewise steady reconstruction from the cell-averaged conserved variables fU i g. It consists of reconstructing within each cell X i local equilibrium U eq;i ðxÞ and perturbation dU i ðxÞ parts. The local equilibrium is fixed in two substeps. First, the cell average U i is projected onto a cell average that is consistent with the steady states of interest. This is of course trivial for the density and the momentum components. For the energy component, an estimate of the cellaveraged internal density is needed. A natural choice is provided by Note that this choice is consistent with the steady states of interest (i.e., it is exact at equilibrium when qv x;i 0). With this we then have the local equilibrium projected cell-averaged conserved variables The second substep matches a steady state of interest profile Eq. (197) to the equilibrium projected cell averages U eq;i by Eq. (47). To this end, the local equilibrium profiles are anchored at the cell center x i 18 # eq;i ðxÞ ¼ # 0;i and H eq;i ðxÞ ¼ where the # 0;i , H 0;i and / i ¼ /ðx i Þ are point values of the equilibrium's natural variable, thermodynamic potential and gravitational potential at cell center. It is assumed that the gravitational potential /ðxÞ can be evaluated anywhere needed, either exactly or approximately (up to some order OðDx Þ). The anchor values # 0;i and H 0;i are then set by requiring that x a q eq;i # 0;i ; H eq;i ðx i;a Þ À Á ; x a qe eq;i # 0;i ; H eq;i ðx i;a Þ À Á : In general, Eq. (201) represents a system of two (non-linear) equations for the anchor values # 0;i and H 0;i of the local equilibrium reconstruction profile in cell X i . The system can efficiently be solved by a (hybrid) Newton method (see, e.g., Press et al. 1993;Dennis and Schnabel 1996). We remark that the derivatives needed for the Jacobian matrix computation in Newton's method are of thermodynamic nature and provided by any EoS. Moreover, the initial guesses for # 0;i and H 0;i can directly be computed from the local equilibrium projected cell averages U eq;i through the EoS, i.e., # 0;i ¼ #ðq i ; qe eq;i Þ and H 0;i ¼ Hðq i ; qe eq;i Þ: Note that the latter values are second-order approximation of the values at cell center and are therefore already pretty close to the solution. Furthermore, these initial guesses are already the solution of the system Eq. (201) for a second-order scheme that uses the midpoint rule. Equipped with the local equilibrium reconstruction, the local equilibrium perturbation dU i ðxÞ can be computed following the recipe in Sect. 2.4.1. As a result, we obtain a piecewise steady reconstruction WR for barotropic hydrostatic equilibrium. Likewise, the well-balanced source term discretization is obtained as described in Sect. 2.4.2. The latter can also be constructed such that the scheme is total (fluid and gravitational) energy-conserving. This can be achieved either by evolving the total energy Eq. (169) directly or by discretizing the source term in Eq. (170) appropriately. 19 This concludes the description of the required ingredients to assemble a onedimensional finite volume scheme that preserves a discrete form of barotropic hydrostatic equilibrium. A particularity of that steady state Eq. (195) is its validity in several space dimensions. Therefore, a well-balanced finite volume scheme for multi-dimensional barotropic hydrostatic equilibrium can readily be developed according the recipe in Sect. 2.6.

Discrete steady states
The third category of well-balanced schemes for the Euler equations avoid any assumption on the thermal stratification of the steady states. In effect, they aim at preserving a consistent discretization of the PDE underlying the steady states of interest. For hydrostatic equilibrium, this is rp ¼ Àqr/; which the methods in this category solve numerically (given the density and gravitational acceleration) for the hydrostatic pressure. Therefore, the methods in this category are, in some sense, truly in the spirit of the piecewise steady reconstruction as early advocated by van Leer (1984), Eulderink and Mellema (1995), Mellema et al. (1991). However, this endeavor is difficult in general, especially in several space dimensions. We will come back briefly to this issue at the end of the section.
Many schemes of this type have been developed in the literature: • LeVeque and Bale (1999) developed second-order finite volume schemes that preserve hydrostatic equilibrium and steady adiabatic flow (see also LeVeque et al. 1998). The schemes are based on the quasi-steady wave-propagation algorithm of LeVeque (1998). The method replaces the piecewise constant solution representation within a cell with two constant states separated by a single discontinuity at the middle of the cell. The jump within the cell is chosen such that it exactly cancels out the source term. This leads to modified Riemann problems at cell interfaces that only encode perturbations from the steady state. The steady states are preserved with second-order accuracy by construction, and perturbations are propagated with the same accuracy within the wavepropagation algorithm. • Fuchs et al. (2011) presented well-balanced second-order finite volume schemes for stratified non-isothermal magnetic atmospheres (see also Fuchs et al. 2010b). The schemes are based on an extension beyond the isothermal case developed by Fuchs et al. (2010a) and are capable of preserving certain non-isothermal magnetostatic equilibria. The schemes are applied to the simulation of waves in the outer solar (chromosphere and corona) and other stellar atmospheres. A similar approach is used by Krause (2019). • Vides et al. (2014) developed a second-order finite volume schemes for selfgravitating astrophysical flows. Although the schemes are not strictly wellbalanced 20 (according to the authors), they nevertheless display a substantial improvement over a standard scheme. The schemes are based on a relaxationtype Riemann solver that incorporates the gravity source terms, which effectively couples more strongly the fluid to the gravitational forces. Padioleau et al. (2019) extend the method to low Mach number flow regimes and apply it to compressible convection. • Desveaux et al. (2015) constructed a well-balanced finite volume scheme to capture non-explicit hydrostatic equilibria (see also Desveaux 2013;Desveaux et al. 2014). It is based on a relaxation-type Riemann solver able to resolve hydrostatic equilibrium by directly including the gravity source terms. The resulting scheme preserves a spatially second-order accurate discrete form of hydrostatic equilibrium, but perturbations are only evolved with first order accuracy. • Käppeli and Mishra (2016) designed well-balanced second-order finite volume schemes for hydrostatic sates. They are based on a local discrete hydrostatic reconstruction that directly integrates the equilibrium pressure given the density and gravity forces (see also Käppeli and Mishra 2015). Popov et al. (2019) present some first applications of the schemes to multi-dimensional stellar structure calculations. Moreover, they present an improved discretization of the gravity source term in the energy equation that avoids unphysical energy changes when simulating quasi-stationary convection for extended time periods. Grosheintz-Laval (2021) adapted the schemes to semi-structured icosahedral grids. Berberich et al. (2021b) generalize the second-order schemes to arbitrary orders of accuracy. • Franck and Mendoza (2016) constructed well-balanced asymptotic-preserving finite volume schemes for the Euler equations with gravity and friction source terms. The schemes are based on a Lagrange?remap approach in combination with a relaxation procedure specifically designed to capture the asymptotic limit induced by the friction source term (i.e., the scheme produces consistent and stable approximations for arbitrarily high friction coefficients). Moreover, the schemes are capable of preserving an arbitrary high-order discretization of hydrostatic equilibrium, but perturbations are only evolved with first order accuracy. An extension to two-dimensional unstructured meshes is also presented. • Chertock et al. (2018) developed second-order central-upwind finite volume schemes capable of preserving hydrostatic equilibrium. The schemes are based on a purely conservative reformulation of the equations that avoids the source terms by introducing global fluxes, which essentially corresponds to an integration of the hydrostatic pressure over the whole computational domain. The inherent viscosity of the central-upwind scheme is tweaked by a smooth cutoff function that tends towards zero when the computed solution is locally close to a steady state. See also Gascón and Corberán (2001), Caselles et al. (2009) for a closely related approach. • Varma and Chandrashekar (2019) extended the approach of Chandrashekar and Klingenberg (2015) to arbitrary thermal stratification. To construct the hydrostatic reconstruction and well-balanced source term discretization, they first rewrite the hydrostatic equilibrium in a particular form and discretize the resulting expressions. Effectively, this allows the parametrization of the assumed subcell thermal equilibrium in the piecewise steady reconstruction.
In the following, we briefly sketch in a one-dimensional setting the well-balanced schemes of Käppeli and Mishra (2016) in their arbitrary high-order formulation by Berberich et al. (2021b). Again, this (to be sure not fully impartial) choice is motivated by the flexibility of the approach as it easily adapts to any EoS. The plan is to prepare all the necessary ingredients to apply the recipe outlined in Sect. 2.4. As mentioned earlier, it is currently a challenge to extend the well-balanced schemes of this category to the multidimensional case. Some discussion of this issue is provided at the end of the one-dimensional description.
First of all, we need to be able to compute the steady states of interest. For onedimensional hydrostatic equilibrium, dp eq dx we obtain by direct integration Here, q eq and p eq ðxÞ denote the hydrostatic density and pressure, and p 0 is the pressure at some reference coordinate x 0 . The equilibrium conserved and primitive variables are then where the equilibrium energy density qe eq ðxÞ is computed through the EoS. Note that we are again distinguishing the exact u eq ðxÞ and the approximate U eq ðxÞ hydrostatic equilibrium state (see Eq. (45)). Both, the equilibrium density and the gravitational potential will be approximated by (polynomial) reconstruction and interpolation (up to some desired accuracy OðDx Þ).
The second ingredient is the piecewise steady reconstruction, whose goal is to build accurate equilibrium subcell profiles from the cell-averaged conserved variables fU i g that are consistent with the steady states of interest. This local equilibrium profile is found by fitting an equilibrium Eq. (206) to the cell's average U i . The first substep is to perform an equilibrium projection of U i onto a cell average U eq;i consistent with hydrostatic equilibrium. A convenient option is as in the barotropic equilibrium case provided by where the cell-averaged equilibrium internal energy density is estimated directly from the cell average U i as This choice is exact at hydrostatic equilibrium qv x;i 0 (i.e., it is consistent with the steady states of interest). The second substep matches a hydrostatic equilibrium profile Eq. (205) to the equilibrium projected cell average U eq;i by Eq. (47). For that purpose, the local equilibrium pressure profile p eq;i ðxÞ is anchored at the cell center where p 0;i is a point value of the equilibrium pressure at cell center, q eq;i ðxÞ is the local equilibrium density reconstruction and d/=dx ð Þ i ðxÞ the gravitational 21 Any other anchor point is possible. However, the cell center is again a convenient choice. acceleration in cell X i . Before we can evaluate Eq. (209), we have to choose a subcell representation for these quantities. The gravitational acceleration is computed by (polynomially) interpolating the gravitational potential f/ i g and differentiating (even if the potential is analytically known and the derivative could be evaluated exactly). An obvious choice for the local equilibrium density reconstruction q eq;i ðxÞ is the one provided by the standard reconstruction procedure for the density (Berberich et al. 2021a): However, other choices are possible. See Fig. 7 for a few possibilities. This influences (unsurprisingly) the accuracy OðDx Þ to which the equilibrium profiles are computed and the overall stencil size (Berberich et al. 2021a). Note that since the integrand in Eq. (209) is a simple polynomial, the integral can be easily evaluated analytically. The anchor value p 0;i is then fixed by requiring that the pressure profile Eq. (209) matches with the cell-averaged equilibrium internal energy density x a qe q eq;i ðx i;a Þ; p 0;i À This is a scalar equation for the pressure p 0;i at the cell center and it can be efficiently solved iteratively by, e.g., a (hybrid) Newton method (see, e.g., Press et al. 1993;Dennis and Schnabel 1996). The iteration starts with the pressure estimated from the cell-average mass and internal energy density Eq. (208) through the EoS. The initial guess is a second-order approximation of the pressure point value at the cell center and is, therefore, already quite close to the sought solution. Moreover, it can be solved analytically for second-order accurate schemes or simple EoS like the ideal gas law (arbitrary order). The local equilibrium reconstruction in cell X i is then U eq;i ðxÞ ¼ q eq;i ðxÞ 0 qeðq eq;i ðxÞ; p eq;i ðxÞÞ The local equilibrium perturbation dU i ðxÞ can now be computed as described by the recipe in Sect. 2.4.1. This gives us the complete piecewise steady reconstruction WR for arbitrarily stratified hydrostatic equilibrium. Along the same lines, the wellbalanced source term discretization follows from Sect. 2.4.2: Remarkably, this is the same discretization as provided by a standard (high-order) finite volume scheme. This is a welcome simplification when implemented into an existing solver. As for the barotropic well-balanced schemes, the source term in the energy equation can be avoided by either evolving the total (fluid and gravitational) energy or by discretizing appropriately the source terms in Eq. (170). Such total energy-conserving schemes are advantageous for the long-term simulation of nearequilibrium configurations. Unfortunately, the just described one-dimensional well-balanced scheme does not generalize to the multi-dimensional case in a straightforward manner (e.g., by following the recipe in Sect. 2.6). The culprit is that the discrete (re)construction of the hydrostatic equilibrium profile via approximation of Eq. (205), is typically dependent on the path C (from the reference coordinate x 0 to x). Only under special circumstances, such as constant density or alignment of gravity forces with one coordinate axis, a unique construction of the hydrostatic profile in such a direct way is feasible. As a matter of fact, the numerical integration of the hydrostatic profile via Eq. (214) gives sensible results (i.e., path independent) only if a discrete curl operator applied to the reconstructed integrand (q eq r/) vanishes everywhere (i.e., the integrability condition for hydrostatic equilibrium is fulfilled in a discrete sense). Therefore, the direct (line) integration approach seems hopeless. A more promising approach may rely on trying to solve a certain boundary value problem locally: Fig. 7 Equilibrium pressure profile for a polytrope of index n ¼ 1 (c ¼ 2) over the radial interval [0, 1.2] discretized by 6 cells: Examples for the equilibrium density reconstruction q eq;i ðxÞ (based on the symmetric stencil fq Àk ; . . .; q þk g of size 2k þ 1) and gravitational potential interpolation / i ðxÞ (based on the symmetric stencil f/ ÀkÀ1 ; . . .; / þkþ1 g of the size 2k þ 3) for the construction of the hydrostatic pressure profile p eq;i ðxÞ from Eq. (209). The left panel displays the equilibrium density as reconstructed with piecewise constants (k ¼ 0), the gravitational potential as interpolated with piecewise parabolic (k ¼ 1) and the resulting piecewise linear equilibrium pressure profile. The right panel shows the accuracy of the equilibrium pressure profile for several stencil size choices (k ¼ 0; 1; 2; 3) together with the estimated order of convergence (EOC). We empirically observe that the equilibrium pressure profile has accuracy ¼ 2ðk þ 1Þ, which is increased by one than one would have expected from the order of the equilibrium density reconstruction. We (presumably) attribute this increase in accuracy to a similar phenomena appearing in numerical integration, where even-degree (i.e., symmetric) Newton-Cotes quadrature rules have a by one higher order of accuracy than expected For example, a relaxation method could solve the elliptical PDE within a cell to construct its local hydrostatic pressure distribution (assuming the pressure in the neighboring cells is given). At an equilibrium state, the hope is that the relaxation approach realizes that the cell's pressure is in equilibrium with its surrounding cells and local gravity forces. Away from an equilibrium state, certain criteria have to be devised to select a plausible local equilibrium state. However, this is beyond the scope of the present review. However, let us end the present section on an optimistic note. Although not wellbalanced in a multi-dimensional sense, it turns out that the above well-balanced scheme considerably improves the simulation of nearly hydrostatic configurations even in this suboptimal case.

Numerical examples
In this section, we showcase some typical numerical test problems used to assess the performance of well-balanced schemes in the context of the Euler equations. We here deliberately focus on standard test problems that are easily replicable by an interested reader wishing to test his or her implementation/scheme. Using wellbalanced schemes in concrete astrophysical applications usually brings in a plethora of complications (initial conditions, appropriate boundary conditions, microphysics, etc.) and goes beyond the scope of this review. Before we proceed, we give general comments on how the problems below are initialized, boundary conditions are handled, the timescale the simulations are run, and how the accuracy is measured.
Given a point-wise defined initial condition in conserved variables u 0 ¼ u 0 ðxÞ, the initialization of the simulation depends on the nature of the chosen scheme. For finite volume schemes, the cell-averaged conserved variables at the initial time U 0 i are computed by (sufficiently accurate) numerical integration An appropriate multi-dimensional quadrature rule is used for multi-dimensional problems. Often, it is convenient to formulate the initial conditions in a different set of variables, such as the primitive variables w ¼ ½q; v; p T , and apply a point-wise transformation u 0 ðxÞ ¼ uðw 0 ðxÞÞ in the above initializations. Boundary conditions are usually a delicate and strongly application-dependent matter, even more so when the dynamics of interest occurs close to a steady state. Typically, the local equilibrium profile found in the piecewise steady reconstruction in the last physical cells is extrapolated into an appropriate number of ghost cells. For finite volume schemes, the cell averages in the ghost cells are computed by U 0 k ¼ Q k ðU eq;1 Þ for k\1; Another option is to freeze the data in the ghost cells to the initial equilibrium conditions. In multiple dimensions, the boundary conditions are applied in a direction-by-direction manner.
To characterize a timescale on which a model reacts to perturbations of its equilibrium, we define the sound crossing time where c s denotes the speed of sound and the integral has to be taken over the extent of the steady state of interest C eq . The sound crossing time corresponds to the time it takes for a sound wave to propagate back and forth through the equilibrium configuration. It gives a measure of how quickly a steady configuration reacts to any perturbations of its equilibrium. The accuracy of the schemes is quantified by computing the errors where k:k denotes some norm, usually the L 1 -norm. Here q is a selected quantity of interest (e.g., density, pressure, velocity, ...) and q ref is a reference solution. The reference solution is the stationary state to be maintained discretely or a numerically obtained solution on a very fine grid. Although the comparison with a numerically obtained reference solution does not provide rigorous evidence of convergence, it nevertheless indicates a meaningful measure of the errors.

Hydrostatic atmospheres
The simplest setup one can image is that of a one-dimensional hydrostatic atmosphere subject to a constant gravitational acceleration g: Despite its apparent simplicity, it has applications in situations where gravitational forces change slowly with respect to other quantities of interest, such as numerical weather prediction, climate modeling of (exo-) planets, and simulation of waves in stellar atmospheres. As a matter of fact, Eq. (220) describes only a mechanical equilibrium of the atmosphere. To uniquely integrate Eq. (220), one needs in addition to a density q 0 and pressure p 0 at the base of the atmosphere x 0 also a thermal stratification. This is simply because the pressure depends on one 22 additional thermodynamic quantity besides density such as the temperature T (p ¼ pðq; TÞ) or the specific entropy s (p ¼ pðq; sÞ). We assume a monoatomic ideal gas law EoS where R is the gas constant, c v the specific heat at constant volume and c the ratio of specific heats.
For an isothermal atmosphere T ¼ T 0 ¼ const., one can then immediately integrate Eq. (220) to yield Here H 0 ¼ RT 0 =g is the so-called scale height. Similarly, under isentropic conditions s ¼ s 0 ¼ const., Eq. (220) can be analytically integrated to p eq ðxÞ ¼ e s 0 =c v qðxÞ c and q eq ðxÞ ¼ q cÀ1 with s 0 ¼ c v lnðp 0 =q c 0 Þ. Note that the isentropic atmosphere has a surface at a finite height while the isothermal one extends to infinity.
Following Käppeli and Mishra (2016), we set the computational domain to X ¼ ½0; 2, and the EoS parameters to R ¼ 1, c v ¼ 3R=2 and c ¼ 5=3. For the isentropic atmosphere, we set the density and pressure to q 0 ¼ 1 and p 0 ¼ 1 at the base x 0 ¼ 0. The resulting isentropic atmosphere has a scale height HðxÞ ¼ À dx d ln P decreasing linearly from Hð0Þ ¼ 1 at the bottom to Hð2Þ ¼ 0:2 at the top and a sound crossing time of s sound % 4:3. Similarly for the isothermal atmosphere, we set the density and pressure to q 0 ¼ 1 and p 0 ¼ 1 at the base x 0 ¼ 0. The resulting isothermal atmosphere has a constant scale height of H 0 ¼ 1 and a sound crossing time of s sound % 3:1. The density and pressure profiles of both atmospheres are shown in Fig. 8. Similar setups have been used in most (if not all) publications for designing well-balanced methods for the Euler equations with gravity.
Next, we present several test cases based on the simple setting of isothermal and isentropic hydrostatic atmospheres. The results are obtained with the well-balanced second-order finite volume scheme of Käppeli and Mishra (2016), also outlined in Sect. 3.3.3. Below it will be termed as the ''well-balanced (WB) discrete scheme'' as it belongs to the category of well-balanced schemes with an underlying discrete piecewise steady reconstruction. For the sake of comparison, we also show the results obtained with a standard second-order finite volume scheme obtained by switching off the piecewise steady reconstruction. Below it will be termed the ''standard (STD) scheme.''

Well-balanced property
The first thing to verify is, of course, the well-balanced property of the scheme: Is it able to preserve the steady state it was designed for up to machine precision? Indeed, due to the finite precision of the computer's approximation of real numbers, exact balance can, in general, not be expected.
To this end, one initializes the computation with so-called well-prepared initial data, i.e., the discrete steady state that the well-balanced method is designed to balance. These well-prepared initial conditions are then evolved with the wellbalanced scheme for a certain time characteristic for the considered steady state, e.g., a multiple of the sound crossing time Eq. (218). At the end of the computation, one computes the difference between the initial and final states in some norm (e.g., the L 1 norm).
The results of such one-dimensional computations with the well-balanced and the standard scheme are shown in Fig. 9. From the left panel of the figure, it is clear that the well-balanced scheme maintains the discrete hydrostatic equilibrium down to machine precision. On the other hand, the standard scheme cannot preserve the discrete equilibrium. However, note that the error for the standard scheme gets smaller with increasing resolution. 23 Indeed, in the limit of an extremely high resolution, the unbalanced scheme can also preserve the equilibrium. This is a matter of consistency of the numerical method with the PDE.
The right panel of Fig. 9 shows the evolution of the maximum absolute velocity for a given resolution (N ¼ 512 cells) and up to one hundred sound crossing times. The advantage of well-balanced methods for long-term integration of nearequilibrium flows becomes quite clear: the local truncation errors of the standard scheme piles up in each step and create spurious flow. Conversely, only the (unavoidable) round-off errors accumulate for the well-balanced scheme. This makes such well-balanced schemes well suited to study numerically natural phenomena that occur on a hydrostatic background.

Wave propagation: shake the base
The next test assesses the ability of a scheme to propagate waves on top of a hydrostatic atmosphere. For brevity, only the results for the isentropic atmosphere are shown. Fuchs et al. (2010a) suggest to excite waves at the bottom of the atmosphere by imposing a periodic velocity perturbation in the lower boundary  23 Interestingly, the standard second-order base scheme shows a third-order convergence towards the steady state. Such a superconvergence has been observed before and attributed to the fact that the base scheme may be asymptotically high order, i.e., it has a higher order of accuracy at a steady state than its design order of accuracy (see Fjordholm et al. 2011).
where m\1 is the boundary cell index (i.e., ghost cell index) and A is the amplitude of the perturbation. The simulation is stopped shortly before the excited waves hit the upper boundary at t f ¼ 1:8. Practically, this test is a simplified version of similar setups used in the study of wave propagation in stellar atmospheres (see, e.g., Bogdan et al. 2003;Rosenthal et al. 2002;Fuchs et al. 2010c and references therein). The setup is run for three amplitudes A ¼ 10 À8 ; 10 À6 ; 10 À1 for several resolutions with the well-balanced/standard scheme. The obtained results are compared with a numerical reference solution computed with the well-balanced scheme and an N ¼ 8192 resolution. The left panel of Fig. 10 displays the errors in velocity.
For the small amplitude A ¼ 10 À8 case, we observe the superiority of the wellbalanced versus the standard scheme, i.e., the committed velocity errors are orders of magnitude smaller. The well-balanced scheme shows a rough second-order convergence. Although way off, the standard scheme seems to show third-order (super)convergence already observed in the well-balanced property test. The left panel of Fig. 11 shows the velocity profile for the standard and the well-balanced schemes for N = 512, together with the reference solution. The well-balanced scheme can resolve the wave pattern very accurately. On the other hand, the standard scheme shows spurious deviations because of its inability to resolve the hydrostatic background properly.
The standard and well-balanced schemes do equally well for the large amplitude A ¼ 10 À1 case. Both show an order of convergence close to one, which is expected because the large amplitude waves quickly steepen into saw-tooth waves, propagating up the atmosphere. The velocity profile for both schemes and the reference solution are shown in the right panel of Fig. 11. This case shows that the well-balanced piecewise steady reconstruction does not destroy the shock-capturing properties of the base scheme. The intermediate amplitude A ¼ 10 À6 case is interesting. The well-balanced scheme is clearly superior at low resolutions ( 128). The standard scheme shows superconvergence in this regime with roughly order three. This is the regime where the hydrostatic atmosphere dominates the committed error, while the wave pattern is totally unresolved. At higher resolutions ( [ 128), the wave pattern dominates the committed error, and the expected second-order accuracy is recovered. The wellbalanced scheme shows an order of convergence of roughly two over the entire resolution range.

Wave propagation: pressure bump
The next test was introduced by LeVeque and Bale (1999). A pressure perturbation is added at the center of the atmosphere, pðxÞ ¼ p eq ðxÞ 1 þ Ae À200ðxÀ1Þ 2 ; which excites two acoustic pulses propagating downwards/upwards. The simulations are stopped shortly before the pulses reach the domain boundaries. An advantage of this setup is that boundary conditions do not play a role, which can be delicate, especially for high-order schemes. For this test, we only show results for the isothermal atmosphere. Similar results are obtained for the isentropic atmosphere. The setup is also run until t f ¼ 0:4 for three different amplitudes A ¼ 10 À8 ; 10 À6 ; 10 À1 and several resolution with the well-balanced/standard scheme. The results are compared to a numerically obtained reference solution computed with the well-balanced scheme at an N ¼ 8192 resolution. The right panel of Fig. 10 shows the errors in velocity and Fig. 12 shows the velocity profiles for the small and large amplitudes cases. The figures show clearly that similar conclusions to the preceding test can be drawn. Hence, we do not repeat them here for brevity.

Polytrope
The next series of test problems model considers very simple stellar models known as polytropes. A polytrope is a static configuration of an adiabatic gaseous sphere held together by self-gravitation (see, e.g., Chandrasekhar 1967;Kippenhahn et al. 2012 in spherical symmetry. Here r is the radial coordinate and G is the gravitational constant. The purely mechanical equilibrium constraints are supplemented by a thermal equilibrium in the form of a barotropic relation as A relation of this form is called a polytropic relation with polytropic constant K and polytropic exponent c. Hence the name polytrope. Although a polytropic relation is often only an approximate model, it is an exact relation when the pressure is dominated by a (non-or relativistic) completely degenerate electron gas. For instance, the latter is the case in white dwarfs and iron cores of evolved massive stars. With the polytropic relation Eq. (226), Eqs. (224) and (225) can be combined into a single equation Fig. 11 Plot of the velocity profile for the small (left) and large (right) amplitude waves propagating up the isentropic atmosphere. The solid/dashed red lines are the solutions obtained with the standard/wellbalanced scheme with N ¼ 512, respectively. In solid blue is also shown a reference solution obtained with the well-balanced scheme and N ¼ 8192 which is known as the Lane-Emden equation. For three values c ¼ 6=5; 2; 1, the Lane-Emden equation can be solved analytically. Specifically for c ¼ 2, the density and pressure are given by q eq ðrÞ ¼ q C sinðarÞ ar and p eq ðrÞ ¼ KqðrÞ 2 : Here q C is the central density of the polytrope and The gravitational potential is given by /ðrÞ ¼ À2KqðrÞ: Note that a c ¼ 2 polytrope has a finite surface radius r surface ¼ p=a. Moreover, it fulfills the isentropic equilibrium hðrÞ þ /ðrÞ ¼ const ¼ 0 for any r ! 0. Following the setup of Käppeli and Mishra (2014), we set the computational domain to a three-dimensional cube X ¼ ½À1=2; þ1=2 3 uniformly discretized by N 3 cells. The model constants are set to q C ¼ G ¼ K ¼ 1 and we assume an ideal gas law EoS with c ¼ 2. The resulting radial density, pressure and gravitational potential profiles are shown in the left panel of Fig. 13. The polytrope has a sound crossing time from its center r ¼ 0 to r ¼ 1=2 of s sound % 0:74. Although the static configuration is built from an inherently three-dimensional physical problem, a star, it can also be applied in a Cartesian one-and two-dimensional setting as a test with non-constant gravitational acceleration (see, e.g., Li andXing 2018b, a, 2016a;Qian et al. 2018;Grosheintz-Laval and Käppeli 2019;Berberich et al. 2021b). The Fig. 12 Plot of the velocity profile for the small (left) and large (right) amplitude pressure perturbations on the isothermal atmosphere. The solid/dashed red lines are the solutions obtained with the standard/ well-balanced scheme with N ¼ 512, respectively. In solid blue is also shown a reference solution obtained with the well-balanced scheme and N ¼ 8192 simulation of a polytrope in a general relativistic context with well-balanced schemes was also considered by Gosse (2015).
We present several test cases based on the polytrope. We compare the results of a standard, unbalanced scheme with two well-balanced schemes. The first wellbalanced scheme is the second-order finite volume scheme of Käppeli and Mishra (2014), also outlined in Sect. 3.3.2, based on a barotropic piecewise hydrostatic reconstruction-termed the barotropic well-balanced scheme below. The second well-balanced scheme is the second-order finite volume scheme of Käppeli and Mishra (2016), also outlined in Sect. 3.3.3, based on a discrete piecewise hydrostatic reconstruction-termed the discrete well-balanced scheme below. The standard scheme is obtained by switching off the piecewise steady reconstruction in the discrete well-balanced scheme, resulting in a second-order finite volume scheme based on piecewise linear reconstruction in primitive variables 24 .

Well-balanced property
The first verification is again the well-balanced property of the scheme. The threedimensional simulation is initialized with the hydrostatic polytrope, qðx; y; zÞ ¼ q eq ðrÞ; vðx; y; zÞ ¼ 0; pðx; y; zÞ ¼ p eq ðrÞ; /ðx; y; zÞ ¼ /ðrÞ; with r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x 2 þ y 2 þ z 2 p using the midpoint rule to compute the cell-averaged conserved variables and is evolved for twenty sound crossing times t f ¼ 20 s sound % 14:8 with the three schemes.
We remark that the initial conditions fulfill the isentropic equilibrium h þ / ¼ const. exactly in a point-wise fashion. Therefore, the initial cell-averages Fig. 13 Left panel: Radial density, pressure and gravitational potential profile of the polytrope up to the corners of the three-dimensional domain (r ¼ ffiffi ffi 3 p =2 % 0:87). Right panel: Maximum absolute radial velocity as a function of time in the simulations of the three-dimensional hydrostatic polytrope for the standard unbalanced/barotropic/discrete well-balanced second-order finite volume schemes 24 Similarly, the standard scheme could be obtained by switching off the piecewise steady reconstruction in the barotropic well-balanced scheme, resulting in a second-order finite volume scheme based on piecewise linear reconstruction in conserved variables. However, both standard schemes give similar results and, therefore, only the results obtained with piecewise linear reconstruction in the primitive variables are shown below. correspond exactly to the discrete hydrostatic equilibrium preserved by the barotropic well-balanced scheme.
The errors in density and pressure at final time t f for the standard unbalanced/ barotropic/discrete well-balanced schemes are summarized in Table 1. We first observe that the barotropic well-balanced scheme produces errors on the order of the machine precision (for double-precision) as is expected. On the other hand, the standard unbalanced scheme suffers from large spurious deviations. This corresponds to the pile-up of the truncation errors at each time step. The errors also show the expected behavior with increasing resolution for a second-order scheme. Interestingly, also the discrete well-balanced schemes shows errors on the order of the machine precision. It turns out that for the special ideal EoS with a ratio of specific heats c ¼ 2, the discrete well-balanced scheme preserves isentropic hydrostatic equilibrium (h þ / ¼ const.) exactly, even in multiple dimensions, without any requirement of alignment of grid axes and gravity force. As a result, the truncation error vanishes by design for the well-balanced schemes and only the (unavoidable) round-off errors sum up. This is further highlighted in Fig. 13, where the maximum absolute radial velocity is shown as a function of time. It is clear that the standard scheme produces spurious deviations from the hydrostatic state.

Wave propagation: pressure bump
To test the capability of the schemes to evolve small perturbations of the multidimensional hydrostatic equilibrium, we add a small Gaussian hump in pressure at the center of the model star pðx; y; zÞ ¼ p eq ðrÞ 1 þ Ae Àr 2 =w 2 ; with amplitude A ¼ 10 À3 and width w ¼ 0:1. The problem is stopped at t f ¼ 0:2 just before the excited waves reach the boundary. A reference solution was computed in one-dimensional spherical symmetry using the well-balanced barotropic scheme with resolution N ¼ 8192.
The errors in radial velocity are displayed in the left panel of Fig. 14. We note that all the schemes show the expected second-order accuracy. However, we also observe that both well-balanced schemes show roughly three orders of magnitude smaller errors than the standard unbalanced scheme. Furthermore, we note that the errors for the well-balanced scheme on the coarsest resolution are smaller than the respective errors of the unbalanced scheme on the finest resolution. This fact is further highlighted in the right panel of Fig. 14 showing scatter plots of radial velocity for the standard scheme at the highest resolution (N 3 ¼ 256 3 ), both wellbalanced schemes at the coarsest resolution (N 3 ¼ 32 3 ) and the reference solution. This underlines the superiority and computational efficiency of well-balanced schemes for simulating small disturbances on top of a stationary state, especially in a multi-dimensional setting. The errors in density and pressure show the same trends.

Wave propagation: explosion
To check the schemes' robustness and shock-capturing properties, we add a localized overpressure region at the center of the model star   The triggered outward propagating explosion is evolved until t f ¼ 0:15. A reference solution was computed in one-dimensional spherical symmetry using the wellbalanced barotropic scheme with resolution N ¼ 8192.
The left panel of Fig. 15 displays the errors in radial velocity as a function of resolution. We observe that all schemes show comparable errors and order of convergence close to unity, as expected for discontinuous solutions. Therefore, we conclude that the piecewise steady reconstruction in both well-balanced schemes does not diminish the robustness of the standard, base shock-capturing scheme. The errors in density and pressure show similar tendencies. The right panel of Fig. 15 shows scatter plots of radial velocity as a function of radius for all the schemes, including the one-dimensional reference solution for comparison. The overpressurized central region quickly expands driving a first strong shock wave outward. As the shock wave moves out, gravity starts to pull back some matter behind it, driving a collapse which then eventually leads to a rebound in the center. This rebound then drives another outward running shock wave. At the final time t f ¼ 0:15, this cycle has been repeated twice (hence the two strong shock waves) and is about to happen again as the matter is pulled back (the negative velocities below r\0:1).