# Energy conserving and well-balanced discontinuous Galerkin methods for the
Euler–Poisson equations in spherical symmetry⋆

Weijie Zhang1, Yulong Xing2, Eirik Endeve3,4
1School of Mathematical Sciences, University of Science and Technology of
China, Hefei, Anhui 230026, P.R. China
2Department of Mathematics, The Ohio State University, Columbus, OH 43210, USA
3Multiscale Methods and Dynamics Group, Oak Ridge National Laboratory, Oak
Ridge, TN 37831, USA
4Department of Physics and Astronomy, University of Tennessee Knoxville, TN
37996, USA

$\star$This manuscript has been authored in part by UT-Battelle, LLC, under
contract DE-AC05-00OR22725 with the US Department of Energy (DOE).
The US government retains and the publisher, by accepting the article for
publication, acknowledges that the US government retains a nonexclusive,
paid-up, irrevocable, worldwide license to publish or reproduce the published
form of this manuscript, or allow others to do so, for US government
purposes. DOE will provide public access to these results of federally
sponsored research in accordance with the
DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan)
Contact e-mail<EMAIL_ADDRESS>e-mail<EMAIL_ADDRESS>e-mail<EMAIL_ADDRESS>
(Last updated xxx; in original form xxx)

###### Abstract

This paper presents high-order Runge–Kutta (RK) discontinuous Galerkin methods
for the Euler–Poisson equations in spherical symmetry. The scheme can preserve
a general polytropic equilibrium state and achieve total energy conservation
up to machine precision with carefully designed spatial and temporal
discretizations. To achieve the well-balanced property, the numerical
solutions are decomposed into equilibrium and fluctuation components which are
treated differently in the source term approximation. One non-trivial
challenge encountered in the procedure is the complexity of the equilibrium
state, which is governed by the Lane–Emden equation. For total energy
conservation, we present second- and third-order RK time discretization, where
different source term approximations are introduced in each stage of the RK
method to ensure the conservation of total energy. A carefully designed slope
limiter for spherical symmetry is also introduced to eliminate oscillations
near discontinuities while maintaining the well-balanced and total-energy-
conserving properties. Extensive numerical examples — including a toy model of
stellar core-collapse with a phenomenological equation of state that results
in core-bounce and shock formation — are provided to demonstrate the desired
properties of the proposed methods, including the well-balanced property,
high-order accuracy, shock capturing capability, and total energy
conservation.

###### keywords:

methods: numerical, supernovae: general, shock waves, gravitation,
hydrodynamics

††pubyear: 2021††pagerange: Energy conserving and well-balanced discontinuous
Galerkin methods for the Euler–Poisson equations in spherical symmetry⋆–Energy
conserving and well-balanced discontinuous Galerkin methods for the
Euler–Poisson equations in spherical symmetry⋆

##  1 Introduction

In this paper, we present high-order discontinuous Galerkin (DG) methods for
the Euler–Poisson equations in spherical symmetry, which have the well-
balanced property to preserve hydrostatic equilibrium states exactly and total
energy conservation property at the same time.

The Euler equations with gravitation have wide applications in geophysical and
astrophysical flow problems. In the case of a time-dependent gravitational
potential, the model can be coupled with the Poisson equation to represent the
self-gravity, which leads to the Euler–Poisson equations. They play an
important role in many geophysical and astrophysical flows, for example, core-
collapse supernova explosions Müller & Steinmetz (1995); Couch et al. (2013);
Müller (2020), star formation Ostriker et al. (2001); McKee & Ostriker (2007),
planet formation Armitage (2011); Simon et al. (2016), and plasma physics
applications Guo (1998); Suzuki (2011). Self-gravitating astrophysical
dynamics are often physically complex, and numerical methods are usually
employed to simulate such complicated systems.

The Euler equations with gravitation belong to the family of hyperbolic
conservation laws with source terms. One of the most important features of
such systems is that they admit non-trivial time-independent steady state
solutions. Well-balanced schemes are introduced to preserve such steady states
exactly on the discrete level and shown to be efficient and accurate for
capturing small perturbations to such steady states. These perturbations may
be at the level of the truncation error of standard numerical schemes and can
be hard to capture with relatively coarse meshes. The well-balanced methods
have been widely studied in the context of the shallow water equations over a
non-flat bottom topology, see e.g., Bermudez & Vazquez (1994); LeVeque (1998);
Audusse et al. (2004); Xing & Shu (2005); Noelle et al. (2007); Gallardo et
al. (2007); Xing et al. (2010). In recent years, well-balanced methods for the
Euler equations with static gravity have attracted much attention and have
been developed within several different frameworks; see e.g., Xu et al.
(2010); Käppeli & Mishra (2014); Chandrashekar & Klingenberg (2015); Käppeli &
Mishra (2016); Thomann et al. (2019) for first- and second-order schemes, and
Xing & Shu (2013); Li & Xing (2016b); Ghosh & Constantinescu (2016); Li & Xing
(2016a); Chandrashekar & Zenk (2017); Klingenberg et al. (2019); Veiga et al.
(2019); Grosheintz-Laval & Käppeli (2019); Castro & Parés (2020) for high-
order schemes. Some of these works assume that the desired equilibrium is
explicitly known Klingenberg et al. (2019); Wu & Xing (2021), while others
only need a pre-description of the desired equilibrium Li & Xing (2018), and
work for a class of equilibria. Recently, several works are established
without any information of the desired equilibrium state Käppeli & Mishra
(2016); Franck & Mendoza (2016); Berberich et al. (2021). For the
Euler–Poisson equations considered in this paper, the equilibrium states are
more complicated due to the coupling with the Poisson equation.

For the Euler–Poisson equations, another important feature is that they
conserve the total energy, which is defined as the sum of the potential,
internal, and kinetic energies. In the standard formulation of the
Euler–Poisson equations, the effect of gravity is included as source terms,
and the total energy conservation statement is obtained in a non-trivial way.
Thus, conserving the total energy numerically becomes challenging. For some
systems, e.g., in hydrostatic equilibrium, the total energy can be much
smaller than either the potential or internal energies, which means that even
a small truncation error in standard methods for the potential energy can lead
to a large error in the total energy, and eventually the wrong numerical
solution (Jiang & Goodman, 2011). Fully conservative schemes for the
Euler–Poisson equations, which conserve mass, momentum, and total energy, have
been studied under the framework of finite difference methods in the last
fifteen years. One popular technique is to transfer the energy equation to the
equation for total energy and rewrite the governing equations in conservative
form, see e.g., Jiang et al. (2013). Another popular technique does not
involve the reformulation of the unknown variables, but apply integration by
parts and the mass conservation equation to discretize the source term in the
energy equation, see e.g., Mikami et al. (2008); Hanawa (2019); Mullen et al.
(2021). With a careful approximation of the source term in the energy
equation, one can carry out a rigorous proof to show the conservation of total
energy. In this paper, we adopt the second technique and study it in the
framework of high-order finite element DG methods. We note that we solve the
Euler–Poisson equations in spherical symmetry, where we are unable to
formulate the momentum equation in conservative form. For this reason we do
not consider momentum conservation in this paper (cf. Jiang et al., 2013;
Mullen et al., 2021).

The main objective of this paper is to develop high-order DG methods for the
Euler–Poisson equations, which are well-balanced and at the same time have the
total energy conservation property. The well-balanced DG scheme for the Euler
equations with a time-independent gravitational potential was studied in Li &
Xing (2018), where the key component to achieve the well-balanced property is
to decompose the source into equilibrium and fluctuation components and treat
them differently in the source term approximation. Here we consider the
extension of this technique to the Euler–Poisson equations. One non-trivial
difficulty encountered in the procedure is the complexity of the equilibrium
state, which is now governed by the well known Lane–Emden equation. For total
energy conservation, very recent work presented in Mullen et al. (2021), where
a second-order finite difference, fully conservative scheme was proposed and
studied. Here, the extension to the framework of DG methods is studied, which
involves a special integration by parts and novel second- and third-order
Runge–Kutta (RK) time discretization, where different source term
approximations are introduced in each stage of RK method to ensure the
conservation of total energy. A carefully designed slope limiter in spherical
symmetry is also introduced to eliminate oscillations near discontinuities
while still maintaining the well-balanced and total-energy-conserving
properties. To the best of our knowledge, the design of well-balanced methods
for the Euler–Poisson system has not been studied in the context of DG
methods, and there are no existing Runge–Kutta discontinuous Galerkin (RKDG)
schemes which can conserve the total energy for the Euler–Poisson equations.
This is the first paper trying to tackle both challenges simultaneously.

The main motivating astrophysical application for the present work is the
simulation of core-collapse supernovae (CCSNe) in the context of non-
relativistic, self-gravitating hydrodynamics with DG methods (see also Pochik
et al., 2021). After the collapse of the iron core of a massive star, the
inner core settles into an approximate hydrostatic equilibrium, which is not
easily captured by standard numerical methods, unless relatively high spatial
resolution is used (Käppeli & Mishra, 2016). Moreover, conserving the total
energy in CCSN simulations with standard numerical methods and moderate
spatial resolution is challenging (e.g., Müller et al., 2010). The kinetic
energy of the explosion is a key quantity of interest targeted by CCSN
simulation codes, and is typically on the order of $10^{51}$ erg (or less;
e.g., Lentz et al., 2015; Melson et al., 2015; Burrows et al., 2020). Thus,
for reliable estimates of the explosion energy, the total energy should be
conserved to well within this threshold. The use of high-order, well-balanced,
and energy conserving numerical methods, as developed in this paper, may help
provide reliable estimates for quantities of interest from CCSN simulations at
a reduced computational cost.

The rest of the paper is organized as follows. In Section 2, we introduce the
Euler–Poisson equations, their steady-state solutions, and discuss total
energy conservation. In Section 3, we present the structure-preserving
numerical methods for the Euler–Poisson equations. We start by introducing the
conventional DG methods for the Euler–Poisson equations, and then discuss the
well-balanced modifications and total-energy-conserving source term and time
discretization, which leads to our well-balanced and total-energy-conserving
fully discrete RKDG scheme. In Section 4, numerical examples are given to
verify the properties of our proposed methods. Concluding remarks are provided
in Section 5.

##  2 Mathematical model

In this section, we introduce the Euler equations with self-gravity in
spherical symmetry, and discuss the steady-state solutions and total energy
conservation property of the model.

###  2.1 Euler–Poisson equations

The Euler equations in spherical symmetry take the form

$\displaystyle\frac{\partial\rho}{\partial
t}+\frac{1}{r^{2}}\frac{\partial}{\partial r}\Big{(}\,r^{2}\,\rho
u\,\Big{)}=0,$ (2.1) $\displaystyle\frac{\partial\rho u}{\partial
t}+\frac{1}{r^{2}}\frac{\partial}{\partial r}\Big{(}\,r^{2}\,\big{(}\,\rho
u^{2}+p\,\big{)}\,\Big{)}=\frac{2\,p}{r}-\rho\,\frac{\partial\Phi}{\partial
r},$ (2.2) $\displaystyle\frac{\partial E}{\partial
t}+\frac{1}{r^{2}}\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\big{(}\,E+p\,\big{)}\,u\,\Big{)}=-\rho
u\,\frac{\partial\Phi}{\partial r},$ (2.3)

where $r$ is the radial coordinate, $\rho$ is the mass density, $u$ denotes
the fluid velocity, $p$ is the pressure, and $E=\rho
e+\frac{1}{2}\,\rho\,u^{2}$ is the total non-gravitational energy with $e$
being the specific internal energy. An additional thermodynamic equation to
link $p$ with $(\rho,e)$, called the equation of state (EoS), is needed. For
ideal gases, it is given by

$p=(\gamma-1)\,\rho e,$ (2.4)

where $\gamma$ is the (constant) ratio of specific heats. The gravitational
potential $\Phi$ can be obtained from the density $\rho$ via the Poisson
equation

$\frac{1}{r^{2}}\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\frac{\partial\Phi}{\partial r}\,\Big{)}=4\pi\,G\,\rho,$
(2.5)

where $G$ is the gravitational constant. The coupling of these two models
yield the Euler–Poisson equations in spherical symmetry.

###  2.2 Steady states and the Lane–Emden equation

The Euler equations (2.1)-(2.3) admit the following zero-velocity steady
states:

$\rho=\rho(r),\qquad u=0,\qquad\frac{\partial p}{\partial
r}=-\rho\frac{\partial\Phi}{\partial r}.$ (2.6)

Considering the polytropic hydrostatic equilibrium characterized by

$p=\kappa\rho^{\gamma},$ (2.7)

we can combine (2.5), (2.6) and (2.7) to obtain the steady-state equation

$\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(\frac{r^{2}}{\rho}\kappa\gamma\rho^{\gamma-1}\frac{\partial\rho}{\partial
r}\right)=-4\pi\,G\,\rho,$ (2.8)

which is the equation satisfied by $\rho(r)$. By introducing the quantities
$\theta$ and $n$ defined by

$\displaystyle\rho$
$\displaystyle\equiv\lambda\theta^{n},\qquad\gamma\equiv\frac{n+1}{n},$ (2.9)

with $\lambda\equiv\rho_{c}$ being the value of density $\rho$ at the center
$r=0$, the equation (2.8) can be simplified as

$\frac{(n+1)\kappa\lambda^{\frac{1-n}{n}}}{4\pi\,G}\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial\theta}{\partial r}\right)=-\theta^{n}.$ (2.10)

Let us define the scaled radial coordinate $\xi$ as

$\displaystyle\xi\equiv$
$\displaystyle\frac{r}{\alpha},\qquad\alpha\equiv\sqrt{\frac{(n+1)\kappa\lambda^{\frac{1-n}{n}}}{4\pi\,G}},$
(2.11)

and this equation can be non-dimensionalized into the well-known Lane–Emden
equation for the polytropic hydrostatic equilibrium:

$\frac{1}{\xi^{2}}\frac{\partial}{\partial\xi}\left(\xi^{2}\frac{\partial\theta}{\partial\xi}\right)=-\theta^{n}.$
(2.12)

As a second-order ordinary differential equation for $\theta(\xi)$, it
requires two boundary conditions:

  1. 1.

Since $\lambda\equiv\rho_{c}=\left.\rho\right|_{\xi=0}$ and
$\rho=\lambda\theta^{n}$, we have $\left.\theta\right|_{\xi=0}=1$ at the
center $\xi=0$;

  2. 2.

The polytropic equilibrium (2.7) leads to

$\frac{\partial p}{\partial
r}=\kappa\gamma\rho^{\gamma-1}\frac{\partial\rho}{\partial
r}~{}~{}\propto~{}~{}\frac{\partial\theta}{\partial\xi}.$ (2.13)

We have ${\partial p}/{\partial r}=-\rho~{}{\partial\Phi}/{\partial r}=0$ at
$r=0$ (because there is no mass inside zero radius). Therefore, we conclude
that

$\left.\frac{\partial\theta}{\partial\xi}\right|_{\xi=0}=0.$ (2.14)

###### Remark 2.1

The methods presented in this paper are to preserve the steady state (2.7) for
the ideal EoS (2.4) up to round-off errors, but can deal with problems for
general EoS without preserving the steady states up to machine error.

###  2.3 Total energy conservation

The solutions of the Euler–Poisson system (2.1)-(2.5) satisfy the following
conservation law for the total energy:

$\frac{\partial}{\partial
t}\left(E+\frac{1}{2}\,\rho\,\Phi\right)+\frac{1}{r^{2}}\frac{\partial}{\partial
r}\Big{(}\,r^{2}\left(\big{(}\,E+p\,\big{)}\,u+F_{g}\right)\Big{)}=0,$ (2.15)

where

$F_{g}=\frac{1}{8\pi\,G}\left(\Phi\,\frac{\partial^{2}}{\partial r\partial
t}\Phi-\frac{\partial}{\partial t}\Phi\,\frac{\partial}{\partial
r}\Phi\right)+\rho u\,\Phi,$ (2.16)

which leads to the total energy conservation

$\frac{\partial}{\partial
t}\int_{\Omega}\left(E+\frac{1}{2}\,\rho\,\Phi\right)\,r^{2}\,\mathrm{d}r=0,$
(2.17)

if the boundary fluxes are zero. Here $\frac{1}{2}\,\rho\,\Phi$ is the
canonical gravitational energy density of a self-gravitating system.

Below, we sketch the main derivation steps of (2.15), which will be useful in
the derivation of the total-energy-conserving numerical methods. Let us
decompose the time derivative into two terms as

$\displaystyle\frac{\partial}{\partial
t}\left(E+\frac{1}{2}\,\rho\,\Phi\right)r^{2}$
$\displaystyle=\left(\frac{\partial E}{\partial
t}+\frac{1}{2}\frac{\partial\rho}{\partial
t}\Phi+\frac{1}{2}\rho\,\frac{\partial\Phi}{\partial t}\right)r^{2}$
$\displaystyle=\left(\frac{\partial E}{\partial
t}+\frac{\partial\rho}{\partial
t}\,\Phi\right)r^{2}+\frac{1}{2}\left(\rho\,\frac{\partial\Phi}{\partial
t}-\frac{\partial\rho}{\partial t}\,\Phi\right)r^{2}.$ (2.18)

For the first term, we have

$\displaystyle\left(\frac{\partial E}{\partial t}+\frac{\partial\rho}{\partial
t}\,\Phi\right)r^{2}$ $\displaystyle=\left(-\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\big{(}\,E+p\,\big{)}\,u\,\Big{)}-\rho
u\,\frac{\partial\Phi}{\partial r}\,r^{2}-\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\rho u\,\Big{)}\,\Phi\right)$
$\displaystyle=\left(-\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\big{(}\,E+p\,\big{)}\,u\,\Big{)}-\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\rho u\,\Phi\Big{)}\right)$
$\displaystyle=-\frac{\partial}{\partial
r}\Big{(}\,r^{2}\left(\big{(}\,E+p\,\big{)}\,u+\rho u\,\Phi\,\Big{)}\right),$
(2.19)

which follows from Eq. (2.1) and (2.3). For the second term, we have

$\displaystyle\frac{1}{2}\left(\rho\,\frac{\partial\Phi}{\partial
t}-\frac{\partial\rho}{\partial t}\,\Phi\right)r^{2}$
$\displaystyle=\frac{1}{8\pi\,G}\left(\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\frac{\partial\Phi}{\partial
r}\,\Big{)}\,\frac{\partial\Phi}{\partial t}-\frac{\partial}{\partial
t}\frac{\partial}{\partial r}\Big{(}\,r^{2}\,\frac{\partial\Phi}{\partial
r}\,\Big{)}\,\Phi\right)$
$\displaystyle=\frac{1}{8\pi\,G}\left(\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\frac{\partial\Phi}{\partial r}\frac{\partial\Phi}{\partial
t}\Big{)}-r^{2}\,\frac{\partial\Phi}{\partial
r}\frac{\partial^{2}\Phi}{\partial r\partial t}\right.$ $\displaystyle\hskip
34.14322pt\left.-\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\frac{\partial^{2}\Phi}{\partial r\partial
t}\,\Phi\Big{)}+\frac{\partial}{\partial
t}\Big{(}\,r^{2}\,\frac{\partial\Phi}{\partial
r}\,\Big{)}\,\frac{\partial\Phi}{\partial r}\right)$
$\displaystyle=\frac{1}{8\pi\,G}\left(\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\frac{\partial\Phi}{\partial r}\frac{\partial\Phi}{\partial
t}\Big{)}-\frac{\partial}{\partial
r}\Big{(}\,r^{2}\,\frac{\partial^{2}\Phi}{\partial r\partial
t}\,\Phi\Big{)}\right),$ (2.20)

which follows from Eq. (2.5) and integration by parts. The combination of
these leads to the conservative form of the total energy (2.15).

###### Remark 2.2

We note that the form of the energy flux in Eq. (2.16) is not unique (Jiang et
al., 2013; Mullen et al., 2021). The different energy fluxes will not affect
the numerical methods proposed in this paper, which will be derived based on
the original form (2.1)-(2.5). The energy flux in Eq. (2.16) is introduced
only as a tool for the proof of the total energy conservation property.

##  3 Numerical methods

In this section, we present the high-order, total-energy-conserving, and well-
balanced DG scheme for the Euler–Poisson equations (2.1)-(2.5), which
preserves the polytropic equilibrium (2.8), and at the same time has the total
energy conservation property (2.17) on the discrete level.

###  3.1 Notations

Let us divide the computational domain $\Omega=\\{r:r\in[0,R]\\}$ into
computational cells

$K_{j}=\\{r:r\in[r_{j-\frac{1}{2}},r_{j+\frac{1}{2}}]\\}\quad\mbox{and}\quad\Delta
r_{j}=r_{j+\frac{1}{2}}-r_{j-\frac{1}{2}},$ (3.1)

for $j=1,...,N$. We define the finite dimensional function space

$\mathcal{V}_{h}:=\\{v\in L^{2}(\Omega):\,v|_{K_{j}}\in
P^{k}(K_{j}),\,\forall\,1\leq j\leq N\\},$ (3.2)

where $P^{k}$ denotes the polynomial space up to degree $k$, and let

$\bm{\Pi}_{h}:=\\{(\zeta,\psi,\delta)^{T}:\,\zeta,\psi,\delta\in\mathcal{V}_{h}\\}.$
(3.3)

For any unknown variable $u$, we denote its numerical approximation in the DG
method by $u_{h}$, which belongs to the piecewise polynomial space
$\mathcal{V}_{h}$. For $\psi\in\mathcal{V}_{h}$, the limit values at the cell
boundaries $r_{j+\frac{1}{2}}$ from the left and the right are defined by

$\psi_{j+\frac{1}{2}}^{-}:=\lim_{\epsilon\rightarrow
0^{+}}\psi(r_{j+\frac{1}{2}}-\epsilon),\quad\psi_{j+\frac{1}{2}}^{+}:=\lim_{\epsilon\rightarrow
0^{+}}\psi(r_{j+\frac{1}{2}}+\epsilon).$ (3.4)

We introduce the Gauss-Radau projection, to be used later in designing the
well-balanced methods. For a function $u\in L^{2}(\Omega)$ and $k\geq 1$, we
define its projection $Pu$ into the space $\mathcal{V}_{h}$ as

$\int_{K_{j}}Pu\,\psi\,\mathrm{d}r=\int_{K_{j}}u\,\psi\,\mathrm{d}r,\quad\forall\psi|_{K_{j}}\in
P^{k-1}(K_{j}),$ (3.5)

for every cell $K_{j}$ and

$Pu(r_{j-\frac{1}{2}}^{+})=u(r_{j-\frac{1}{2}}^{+}).$ (3.6)

###  3.2 The approximation of the gravitational potential

Compared with the Euler equations with static gravitational field studied in
Li & Xing (2018); Wu & Xing (2021), the Euler–Poisson equations (2.1)-(2.5)
involve the additional Poisson equation (2.5) which governs the relation
between time dependent $\Phi$ and the density $\rho$. There are extensive
numerical methods that could be used to solve the Poisson equation. Here, we
present the following simple approach to compute $\Phi$ numerically.

Note that the source terms in (2.2) and (2.3) involve only the derivative
$\partial\Phi/\partial r$, however, we will compute the numerical
approximation of both $\partial\Phi/\partial r$ and $\Phi$ in this paper,
denoted by $\partial\Phi_{h}/\partial r$ and $\Phi_{h}$ respectively, as the
latter will be used in the design of total-energy-conserving methods.

We can integrate the Poisson equation (2.5) directly and obtain

$\displaystyle\frac{\partial\Phi_{h}}{\partial r}$
$\displaystyle=\frac{4\pi\,G}{r^{2}}\int_{0}^{r}\rho_{h}\tau^{2}\,\mathrm{d}\tau,$
(3.7) $\displaystyle\Phi_{h}$
$\displaystyle=\Phi_{h}(R)-\int_{r}^{R}\frac{\partial\Phi_{h}}{\partial
r}\,\mathrm{d}r,$ (3.8)

with the boundary conditions ${\partial\Phi_{h}(0)}/{\partial r}=0$ and
$\Phi_{h}(R)=\text{constant}$. The equations (3.7) and (3.8) mean that we
calculate $\frac{\partial\Phi_{h}}{\partial r}$ and $\Phi_{h}$ cell by cell
that

$\displaystyle\frac{\partial\Phi_{h}}{\partial r}(r)$
$\displaystyle=\frac{4\pi\,G}{r^{2}}\int_{r_{j-\frac{1}{2}}}^{r}\rho_{h}\tau^{2}\,\mathrm{d}\tau+\frac{r_{j-\frac{1}{2}}^{2}}{r^{2}}\frac{\partial\Phi_{h}}{\partial
r}(r_{j-\frac{1}{2}}),$ (3.9)

for $r\in K_{j}$, $j=1,...,N$ and

$\displaystyle\Phi_{h}(r)$
$\displaystyle=\Phi_{h}(r_{j+\frac{1}{2}})-\int_{r}^{r_{j+\frac{1}{2}}}\frac{\partial\Phi_{h}}{\partial
r}\,\mathrm{d}r,$ (3.10)

for $r\in K_{j}$, $j=N,...,1$. We set $\Phi_{h}(R)=0$ in the numerical tests
of this paper to observe the total energy conservation up to round-off error.
Note that $\rho_{h}$ is a piecewise polynomial of degree $k$, hence the
integrals in (3.9) and (3.10) can be evaluated exactly over each computational
cell $K_{j}$. The detailed procedure is summarized in the following steps.

  1. 1.

Assume $\rho_{h}$ is $P^{k}$ piecewise polynomial taking the form, for $r\in
K_{j}$, $j=1,...,N$,

$\rho_{h}(r)\Big{|}_{K_{j}}=\sum_{i=0}^{k}\rho_{j,i}\,r^{i}.$ (3.11)

  2. 2.

Compute the integration in (3.9) exactly and obtain $\partial\Phi_{h}/\partial
r$ as

$\displaystyle\frac{\partial\Phi_{h}}{\partial r}(r)=$
$\displaystyle\frac{4\pi\,G}{r^{2}}\left.\sum_{i=0}^{k}\frac{\rho_{j,i}\,\tau^{i+3}}{i+3}\right|_{\tau=r_{j-\frac{1}{2}}}^{r}+\frac{r_{j-\frac{1}{2}}^{2}}{r^{2}}\frac{\partial\Phi_{h}}{\partial
r}(r_{j-\frac{1}{2}})$ $\displaystyle:=$
$\displaystyle\sum_{i=1}^{k+1}g_{j,i}\,r^{i}+\frac{g_{j,-2}}{r^{2}},$ (3.12)

for $r\in K_{j}$, $j=1,...,N$.

  3. 3.

Compute the integration in (3.10) exactly and obtain $\Phi_{h}$ as

$\Phi_{h}(r)=\Phi_{h}(r_{j+\frac{1}{2}})-\left.\left(\sum_{i=1}^{k+1}\frac{g_{j,i}\,\tau^{i+1}}{i+1}-\frac{g_{j,-2}}{\tau}\right)\right|_{\tau=r}^{r_{j+\frac{1}{2}}},$
(3.13)

for $r\in K_{j}$, $j=N,...,1$. Here $\rho_{j,i}$ in (3.11) and $g_{j,i}$ in
(3.12) are the polynomial coefficients of degree $i$ ($i\geq 0$) in the $j$-th
cell for $\rho_{h}$ and ${\partial\Phi_{h}}/{\partial r}$ respectively.
$g_{j,-2}$ in (3.12) are the coefficient of the term ${1}/{r^{2}}$ for
${\partial\Phi_{h}}/{\partial r}$.

###  3.3 The standard DG scheme

In this subsection, we will briefly review the standard DG method for the
Euler–Poisson equations (2.1)-(2.5), which will be used in the numerical
section for comparison. For ease of presentation, we denote the equations
(2.1)-(2.3) as:

$\frac{\partial\bm{u}}{\partial t}+\frac{1}{r^{2}}\frac{\partial}{\partial
r}(r^{2}\bm{f}(\bm{u}))=\bm{s}(\bm{u},\Phi),$ (3.14)

where

$\displaystyle\bm{u}=\left(\begin{matrix}\rho\\\ \rho u\\\
E\end{matrix}\right),~{}\bm{f}(\bm{u})=\left(\begin{matrix}\rho u\\\ \rho
u^{2}+p\\\
(E+p)u\end{matrix}\right),~{}\bm{s}(\bm{u},\Phi)=\left(\begin{matrix}0\\\
\frac{2p}{r}-\rho\frac{\partial\Phi}{\partial r}\\\ -\rho
u\frac{\partial\Phi}{\partial r}\end{matrix}\right).$ (3.15)

To derive the semi-discrete DG scheme, we multiply the equations by $r^{2}$
and test functions, apply integration by parts and replace the boundary value
by a monotone numerical flux, which leads to the following DG scheme: find
$\bm{u}_{h}\in\bm{\Pi}_{h}$ such that for any test function
$\bm{v}=(\zeta,\psi,\delta)^{T}\in\bm{\Pi}_{h}$, it holds that

$\displaystyle\partial_{t}\int_{K_{j}}\bm{u}_{h}\cdot\bm{v}\,r^{2}\mathrm{d}r+r_{j+\frac{1}{2}}^{2}\hat{\bm{f}}_{j+\frac{1}{2}}\cdot\bm{v}_{j+\frac{1}{2}}^{-}-r_{j-\frac{1}{2}}^{2}\hat{\bm{f}}_{j-\frac{1}{2}}\cdot\bm{v}_{j-\frac{1}{2}}^{+}$
$\displaystyle\hskip
34.14322pt-\int_{K_{j}}\bm{f}(\bm{u}_{h})\cdot(\partial_{r}\bm{v})r^{2}\mathrm{d}r=\bm{s}_{j},$
(3.16)

where $\bm{s}_{j}$ is the approximation of
$\int_{K_{j}}\bm{s}(\bm{u},\Phi)\cdot\bm{v}r^{2}\mathrm{d}r$ taking the form

$\bm{s}_{j}=\left(\begin{matrix}0\\\ s_{j}^{[2]}\\\
s_{j}^{[3]}\end{matrix}\right)=\left(\begin{matrix}0\\\
\int_{K_{j}}\left(\frac{2p_{h}}{r}-\rho_{h}\frac{\partial\Phi_{h}}{\partial
r}\right)\psi\,r^{2}\mathrm{d}r\\\ \int_{K_{j}}-(\rho
u)_{h}\frac{\partial\Phi_{h}}{\partial
r}\delta\,r^{2}\mathrm{d}r\end{matrix}\right),$ (3.17)

and $\hat{\bm{f}}$ is the monotone numerical flux. In this paper, to have good
performance in capturing shocks and optimal error convergence rate, we
consider the Harten-Lax-van Leer contact (HLLC) flux (Toro, 2013)

$\hat{\bm{f}}=\hat{\bm{f}}(\bm{u}_{h}^{-},\bm{u}_{h}^{+})=\begin{cases}\bm{f}(\bm{u}_{h}^{-})&\text{
if }0\leq S^{-},\\\
\bm{f}(\bm{u}_{h}^{-})+S^{-}\left(\bm{u}_{*}^{-}-\bm{u}_{h}^{-}\right)&\text{
if }S^{-}\leq 0\leq S^{*},\\\
\bm{f}(\bm{u}_{h}^{+})+S^{+}\left(\bm{u}_{*}^{+}-\bm{u}_{h}^{+}\right)&\text{
if }S^{*}\leq 0\leq S^{+},\\\ \bm{f}(\bm{u}_{h}^{+})&\text{ if }0\geq
S^{+},\end{cases}$ (3.18)

where $S^{-}$, $S^{+}$ and $S^{*}$ are the signal speeds

$\displaystyle
S^{-}=\min\\{u_{h}^{-}-c_{h}^{-},~{}u_{h}^{+}-c_{h}^{+}\\},\quad
S^{+}=\max\\{u_{h}^{-}+c_{h}^{-},~{}u_{h}^{+}+c_{h}^{+}\\},$ (3.19)
$\displaystyle
S^{*}=\frac{p_{h}^{+}-p_{h}^{-}+\rho_{h}^{-}u_{h}^{-}\left(S^{-}-u_{h}^{-}\right)-\rho_{h}^{+}u_{h}^{+}\left(S^{+}-u_{h}^{+}\right)}{\rho_{h}^{-}\left(S^{-}-u_{h}^{-}\right)-\rho_{h}^{+}\left(S^{+}-u_{h}^{+}\right)},$
(3.20)

$c_{h}^{-}$, $c_{h}^{+}$ are the sound speeds calculated from
$\bm{u}_{h}^{-}$, $\bm{u}_{h}^{+}$ respectively, and $\bm{u}_{*}^{+}$,
$\bm{u}_{*}^{+}$ denote the intermediate states which can be computed via

$\bm{u}_{*}^{\pm}=\rho_{h}^{\pm}\left(\frac{S^{\pm}-u_{h}^{\pm}}{S^{\pm}-S^{*}}\right)\left(\begin{matrix}1\\\
S^{*}\\\
\frac{E_{h}^{\pm}}{\rho_{h}^{\pm}}+\left(S^{*}-u_{h}^{\pm}\right)\left(S^{*}-\frac{p_{h}^{\pm}}{\rho_{h}^{\pm}\left(S^{\pm}-u_{h}^{\pm}\right)}\right)\end{matrix}\right).$
(3.21)

The initial condition $\bm{u}_{h,0}\in\bm{\Pi}_{h}$ of the numerical method is
given by

$\bm{u}_{h,0}=P\bm{u}_{ex}(r,t=0),$ (3.22)

where $\bm{u}_{ex}(r,t=0)$ is the exact initial data, and $P$ stands for the
Gauss-Radau projection (3.5)-(3.6).

###  3.4 The well-balanced DG scheme

In this subsection, we will introduce the well-balanced DG scheme which
maintains the polytropic equilibrium (2.8), or equivalently the Lane–Emden
equation (2.12). There are some recent works (Xing, 2014; Grosheintz-Laval &
Käppeli, 2020; Parés & Parés-Pulido, 2021) on designing well-balanced methods
for general steady states including non-zero equilibrium, which will be
studied in future work.

####  3.4.1 Solution of Lane–Emden equation

As illustrated in Section 2.2, the polytropic equilibrium state of the
Euler–Poisson equations is based on the solution of the Lane–Emden equation.
The Lane–Emden equation can be analytically solved (Maciel, 2015) only for a
few special integer values of the index $n$, as outlined below:

Analytical solution for n=0 (i.e., $\gamma=\infty$):
$\displaystyle\theta_{0}(\xi)=1-\frac{1}{6}\xi^{2},$ (3.23) Analytical
solution for n=1 (i.e., $\gamma=2$):
$\displaystyle\theta_{1}(\xi)=\frac{\sin(\xi)}{\xi},$ (3.24) Analytical
solution for n=5 (i.e., $\gamma=\frac{6}{5}$):
$\displaystyle\theta_{5}(\xi)=\frac{1}{\sqrt{1+\frac{1}{3}\xi^{2}}}.$ (3.25)

For all other values of $n$, we must resort to numerical solutions. Rewrite
the equation (2.12) as

$\displaystyle\frac{\partial\theta}{\mathrm{d}\xi}$
$\displaystyle=-\frac{\varphi}{\xi^{2}},\qquad\frac{\partial\varphi}{\mathrm{d}\xi}=\theta^{n}\xi^{2},$
(3.26)

coupled with boundary conditions $\theta(0)=1$ and $\varphi(0)=0$. We denote
them in the vector form by

$\frac{\partial}{\partial\xi}\bm{y}=\bm{F}(\xi,\bm{y})\quad\text{with
}\bm{y}=\left(\begin{aligned} \theta\\\
\varphi\end{aligned}\right)\quad\text{and}\quad\bm{F}(\xi,\bm{y})=\left(\begin{aligned}
-\frac{\varphi}{\xi^{2}}\\\ \theta^{n}\xi^{2}\end{aligned}\right).$ (3.27)

Note that when $\xi=0$, we let $\bm{F}(0,\bm{y}(0))=0$ following the given
boundary conditions. The equations (3.27) is a system of ordinary differential
equations, which can solved by various numerical methods. For example, we can
use the fifth-order Runge-Kutta-Fehlberg technique in Norsett & Wanner (1987)

$\bm{y}_{j+1}=\bm{y}_{j}+h\sum_{i=1}^{s}b_{i}\bm{k}_{i},$ (3.28)

where $\bm{y}_{j}$ denotes the numerical solution at the grid $\xi_{j}$,
$h=\xi_{j+1}-\xi_{j}$ and $\bm{k}_{i}$, $i=1,2,\ldots,s$, is given by

$\displaystyle\bm{k}_{i}$
$\displaystyle=\bm{F}(\xi_{j}+c_{i}h,\bm{y}_{j}+h(a_{i1}\bm{k}_{1}+a_{i2}\bm{k}_{2}+\cdots+a_{i,i-1}\bm{k}_{i-1})).$
(3.29)

with the coefficients $a_{ij}$, $b_{i}$ and $c_{i}$ given in the following
Butcher tableau:

$\displaystyle\begin{array}[]{c|cccccc}0\\\ \frac{1}{4}&\frac{1}{4}\\\
\frac{3}{8}&\frac{3}{32}&\frac{9}{32}\\\
\frac{12}{13}&\frac{1932}{2197}&-\frac{7200}{2197}&\frac{7296}{2196}\\\
1&\frac{439}{216}&-8&\frac{3680}{513}&-\frac{845}{4104}\\\
\frac{1}{2}&-\frac{8}{27}&2&-\frac{3544}{2565}&\frac{1859}{4104}&-\frac{11}{40}\\\
\hline\cr&\frac{16}{135}&0&\frac{6656}{12825}&\frac{28561}{56430}&-\frac{9}{50}&\frac{2}{55}\end{array}$
(3.37)

The numerical solution of (2.12) can be solved with enough accuracy by taking
small enough $h$. We note that the solution of the Lane–Emden equation only
depends on $n$ (i.e. $\gamma$). For each computational example, $\gamma$ is
fixed, hence we can pre-calculate and save the numerical solution $\theta_{n}$
at the beginning of the simulation.

####  3.4.2 Decomposition of the numerical solutions

To design the well-balanced method, we follow the approach in Xing (2014)
where well-balanced methods for the moving water equilibrium of the shallow
water equations are designed. The first step is to separate the numerical
solutions into the well-balanced equilibrium component
$\bm{u}_{h}^{e}\in\bm{\Pi}_{h}$ and the fluctuation part
$\bm{u}_{h}^{f}\in\bm{\Pi}_{h}$ at each time step, which will be elaborated
below.

We start by recovering the desired equilibrium state $\bm{u}^{d}$ which
satisfies the polytropic equilibrium (2.8) and usually does not belong to
$\bm{\Pi}_{h}$. For the given $\gamma$ (or $n$), the solution $\theta_{n}$ of
Lane–Emden equation (2.12) can be pre-computed. Then we evaluate the density
and pressure of the numerical solution $\bm{u}_{h}(r,t)$ at the center $r=0$
and denote them by $\rho_{0}$ and $p_{0}$. By setting
$\kappa=p_{0}/\rho_{0}^{\gamma}$ and
$\alpha=\sqrt{\frac{\gamma}{\gamma-1}\kappa\rho_{0}^{\gamma-2}/(4\pi\,G)}$ in
(2.11), we can define the desired equilibrium state $\bm{u}^{d}$ as

$\bm{u}^{d}(r)=\left(\begin{matrix}\rho_{0}\left(\theta_{n}\left(\frac{r}{\alpha}\right)\right)^{\frac{1}{\gamma-1}}\\\
0\\\
\frac{\kappa}{\gamma-1}\rho_{0}^{\gamma}\left(\theta_{n}\left(\frac{r}{\alpha}\right)\right)^{\frac{\gamma}{\gamma-1}}\end{matrix}\right).$
(3.38)

Suppose the initial condition is in the equilibrium state, i.e.,
$\bm{u}_{ex}(r,0)$ satisfies the polytropic equilibrium (2.8). Note that
although $\bm{u}_{h,0}\in\bm{\Pi}_{h}$ defined in (3.22) is not in perfect
equilibrium, the above procedure can recover the exact equilibrium, i.e., we
can compute $\bm{u}^{d}$ from $\bm{u}_{h,0}$ with
$\bm{u}^{d}=\bm{u}_{ex}(r,0)$.

Next we can define $\bm{u}_{h}^{e}\in\bm{\Pi}_{h}$ as the projection of
$\bm{u}^{d}$ into the DG solution space:

$\bm{u}^{e}_{h}=P\bm{u}^{d},$ (3.39)

and also define the fluctuation term $\bm{u}^{f}\in\bm{\Pi}_{h}$ as:

$\bm{u}^{f}_{h}=\bm{u}_{h}-\bm{u}^{e}_{h}.$ (3.40)

For the $\theta_{n}$ explicitly given in (3.23)-(3.25), the integration in the
definition of the projection in Eq. (3.39) can be evaluated exactly.
Otherwise, the integration is computed by using the values at the Gaussian
quadrature points which can be obtained from interpolation.

###### Remark 3.1

When recovering the desired equilibrium state $\bm{u}^{d}$, two practical
issues in the implementation are noted. First, since the density is positive,
$\theta(\xi)$ should also be positive for robustness of the simulation, and
one should pay attention to the range of the solution of $\theta(\xi)$. If the
analytical solution of the Lane-Emden equation is used, there is a constraint
on the range of $\xi$ for $n=0,1$. For example, $\theta_{0}(\xi)>0$ for
$\xi\in[0,\sqrt{6})$ and $\theta_{1}(\xi)>0$ for $\xi\in[0,\pi)$. If the
numerical solution of the Lane-Emden equation is used, $\theta(\xi)$ may
become negative due to numerical integration errors. Therefore, if there is a
range constraint on $\theta(\xi)$ and a cell $K_{j}$ where the value of
$\theta(\xi)$ is outside of this range constraint, we set
$\left.\bm{u}^{d}\right|_{K_{j}}=0$ for robustness of the simulation. Second,
if the solution is too far away from the equilibrium state, for example, for
the cells $K_{j}$ with

$\rho^{d}(r_{j-\frac{1}{2}})>2\rho_{h,j-\frac{1}{2}}^{+}~{}\text{ or
}~{}p^{d}(r_{j-\frac{1}{2}})>2p_{h,j-\frac{1}{2}}^{+},$ (3.41)

we set $\left.\bm{u}^{d}\right|_{K_{j}}=0$ to avoid the accumulation of error
since $\bm{u}^{d}$ is calculated globally.

####  3.4.3 Well-balanced numerical flux and source term approximation

With the decomposition of the numerical solutions into the equilibrium
component $\bm{u}_{h}^{e}$ and the fluctuation part $\bm{u}_{h}^{f}$ at each
time step, we can now present the well-balanced numerical fluxes and the well-
balanced source term approximation.

We can define the modified cell boundary values of $\bm{u}_{h}$ as

$\bm{u}_{h,j+\frac{1}{2}}^{*,-}=\bm{u}^{d}\left(r_{j+\frac{1}{2}}\right)+\bm{u}_{h,j+\frac{1}{2}}^{f,-},\quad\bm{u}_{h,j+\frac{1}{2}}^{*,+}=\bm{u}^{d}\left(r_{j+\frac{1}{2}}\right)+\bm{u}_{h,j+\frac{1}{2}}^{f,+},$
(3.42)

where $\bm{u}^{d}$ is continuous over the whole computational domain and
defined in (3.38), and $\bm{u}^{f}_{h}$ is defined in (3.40). The well-
balanced numerical flux $\hat{\bm{f}}^{*}$ can be evaluated by

$\displaystyle\hat{\bm{f}}^{*}$
$\displaystyle=\hat{\bm{f}}(\bm{u}_{h}^{*,-},\bm{u}_{h}^{*,+}),$ (3.43)

with $\hat{\bm{f}}$ being the HLLC flux defined in (3.18).

For the well-balanced source term approximation, we follow the main idea in
Xing (2014); Li & Xing (2018), but with some modifications introduced below.
As $s_{j}^{[3]}$ in (3.17) equals to zero automatically at the equilibrium
state, we focus only on the term $s_{j}^{[2]}$. Since $\bm{u}^{d}$ is the
equilibrium solution and continuous, we have

$\displaystyle
r_{j+\frac{1}{2}}^{2}\bm{f}\left(\bm{u}^{d}\left(r_{j+\frac{1}{2}}\right)\right)\cdot\bm{v}_{j+\frac{1}{2}}^{-}-r_{j-\frac{1}{2}}^{2}\bm{f}\left(\bm{u}^{d}\left(r_{j-\frac{1}{2}}\right)\right)\cdot\bm{v}_{j-\frac{1}{2}}^{+}$
$\displaystyle\qquad-\int_{K_{j}}\bm{f}\left(\bm{u}^{d}\right)\cdot(\partial_{r}\bm{v})\,r^{2}\mathrm{d}r-\int_{K_{j}}\bm{s}\left(\bm{u}^{d},\Phi^{d}\right)\cdot\bm{v}\,r^{2}\mathrm{d}r=0,$
(3.44)

where $\Phi^{d}$ is solved exactly from $\rho^{d}$ in (2.5). Because
$\bm{u}_{h}^{e}\in\bm{\Pi}_{h}$ is the projection of $\bm{u}^{d}$ with high-
order accuracy, and $\bm{u}^{d}$ is continuous at the cell interfaces, we have

$\displaystyle
r_{j+\frac{1}{2}}^{2}f^{[2]}\left(\bm{u}^{d}\left(r_{j+\frac{1}{2}}\right)\right)\psi_{j+\frac{1}{2}}^{-}-r_{j-\frac{1}{2}}^{2}f^{[2]}\left(\bm{u}^{d}\left(r_{j-\frac{1}{2}}\right)\right)\psi_{j-\frac{1}{2}}^{+}$
$\displaystyle\quad-\int_{K_{j}}f^{[2]}(\bm{u}^{e}_{h})\,(\partial_{r}\psi)\,r^{2}\mathrm{d}r-\int_{K_{j}}\left(\frac{2p^{e}_{h}}{r}-\rho_{h}^{e}\frac{\partial\Phi^{e}_{h}}{\partial
r}\right)\psi\,r^{2}\mathrm{d}r$ $\displaystyle=\mathcal{O}((\Delta
r_{j})^{k+1}),$ (3.45)

where $f^{[2]}$ denotes the second component of $\bm{f}$ and
$\partial\Phi^{e}_{h}/\partial r$ is evaluated as in (3.7):

$\frac{\partial\Phi^{e}_{h}}{\partial
r}=\frac{4\pi\,G}{r^{2}}\int_{0}^{r}\rho^{e}_{h}\tau^{2}\mathrm{d}\tau,$
(3.46)

with $\frac{\partial\Phi^{e}_{h}(0)}{\partial r}=0$. The approximation of the
source term $\bm{s}_{j}^{\text{wb}}$ is then defined as

$\bm{s}_{j}^{\text{wb}}=\left[0,s_{j}^{[2],\text{wb}},s_{j}^{[3]}\right]^{T},\qquad
s_{j}^{[2],\text{wb}}=s_{j}^{[2]}+s_{j}^{[2],\text{cor}},$ (3.47)

where $s_{j}^{[2]}$ and $s_{j}^{[3]}$ are defined in (3.17) and the correction
term $s_{j}^{[2],\text{cor}}$ takes the form

$\displaystyle s_{j}^{[2],\text{cor}}=$
$\displaystyle~{}r_{j+\frac{1}{2}}^{2}p^{d}\left(r_{j+\frac{1}{2}}\right)\psi_{j+\frac{1}{2}}^{-}-r_{j-\frac{1}{2}}^{2}p^{d}\left(r_{j-\frac{1}{2}}\right)\psi_{j-\frac{1}{2}}^{+}$
$\displaystyle-\int_{K_{j}}p^{e}_{h}\,(\partial_{r}\psi)\,r^{2}\mathrm{d}r-\int_{K_{j}}\left(\frac{2p^{e}_{h}}{r}-\rho_{h}^{e}\frac{\partial\Phi^{e}_{h}}{\partial
r}\right)\psi\,r^{2}\mathrm{d}r,$ (3.48)

which will play an important role in the well-balanced proof.

####  3.4.4 Well-balanced semi-discrete DG scheme

The well-balanced semi-discrete DG scheme can be written as: find
$\bm{u}_{h}\in\bm{\Pi}_{h}$ such that for any test function
$\bm{v}=(\zeta,\psi,\delta)^{T}\in\bm{\Pi}_{h}$, it holds that

$\partial_{t}\int_{K_{j}}\bm{u}_{h}\cdot\bm{v}\,r^{2}\mathrm{d}r=\mathcal{L}_{j}(\bm{u}_{h},\bm{v})=\mathcal{F}_{j}(\bm{u}_{h},\bm{v})+\bm{s}_{j}^{\text{wb}}(\bm{u}_{h},\bm{v}),$
(3.49)

with $\bm{s}_{j}^{\text{wb}}$ defined in (3.47) and

$\displaystyle\mathcal{F}_{j}(\bm{u}_{h},\bm{v})=$ $\displaystyle-
r_{j+\frac{1}{2}}^{2}\hat{\bm{f}}^{*}_{j+\frac{1}{2}}\cdot\bm{v}_{j+\frac{1}{2}}^{-}+r_{j-\frac{1}{2}}^{2}\hat{\bm{f}}_{j-\frac{1}{2}}^{*}\cdot\bm{v}_{j-\frac{1}{2}}^{+}$
$\displaystyle+\int_{K_{j}}\bm{f}(\bm{u}_{h})\cdot(\partial_{r}\bm{v})r^{2}\mathrm{d}r,$
(3.50)

with the source term approximation $\bm{s}_{j}^{\text{wb}}$ defined in (3.47),
and the numerical flux $\hat{\bm{f}}^{*}$ defined in (3.43). We have the
following result on its well-balanced property.

###### Proposition 3.2

The semi-discrete DG scheme (3.49), with initial condition defined in (3.22),
maintains the equilibrium state (2.8) exactly.

###### Proof

Suppose the initial condition is at the equilibrium state (2.8). We will
complete the well-balanced proof in three steps. First, we will show that
$\bm{u}_{h}=\bm{u}^{e}_{h}$ and $\bm{f}^{e}_{h}=0$. By the definition of
$\bm{u}^{d}$ in Eq. (3.38), we can conclude that $\bm{u}^{d}=\bm{u}_{ex}$ as
both are the stationary solutions of (2.8) and share the same value at the
center $r=0$. It then follows from (3.39) and (3.40) that
$\bm{u}^{e}_{h}=\bm{u}_{h}$ and $\bm{u}^{f}_{h}=0$. Moreover, we conclude that
$\partial\Phi_{h}/\partial r=\partial\Phi_{h}^{e}/\partial r$, because
$\partial\Phi_{h}/\partial r$ and $\partial\Phi_{h}^{e}/\partial r$ are
calculated from $\rho_{h}$ and $\rho_{h}^{e}$, respectively, using (3.12), and
$\rho_{h}=\rho_{h}^{e}$.

Second, we would like to show that
$\hat{f}_{j+\frac{1}{2}}^{*,[2]}=p^{d}\left(r_{j+\frac{1}{2}}\right)$. Since
$\bm{u}^{f}=0$, we have that $\bm{u}_{h}^{*,-}=\bm{u}_{h}^{*,+}=\bm{u}^{d}$ at
the interface $r_{j+1/2}$, following the definition (3.42). In Eq. (3.43), we
have

$\displaystyle\hat{\bm{f}}_{j+\frac{1}{2}}^{*}$
$\displaystyle=\hat{\bm{f}}(\bm{u}_{h,j+\frac{1}{2}}^{*,-},\bm{u}_{h,j+\frac{1}{2}}^{*,+})=\bm{f}(\bm{u}_{h,j+\frac{1}{2}}^{*,\pm})$
$\displaystyle=\bm{f}\left(\bm{u}^{d}\left(r_{j+\frac{1}{2}}\right)\right)=\left(\begin{matrix}0\\\
p^{d}\left(r_{j+\frac{1}{2}}\right)\\\ 0\end{matrix}\right),$ (3.51)

where the last equality follows from the zero velocity in the vector
$\bm{u}^{d}$.

Lastly, it is easy to observe that the first and third components of
$\mathcal{L}_{j}$ in Eq. (3.49) are zero. With the source term defined in
(3.47)-(3.4.3), the second component of of $\mathcal{L}_{j}$ can be simplified
as

$\displaystyle\mathcal{L}^{[2]}_{j}(\bm{u}_{h},\bm{v})=\int_{K_{j}}f^{[2]}(\bm{u}_{h})(\partial_{r}\psi)r^{2}\mathrm{d}r-r_{j+\frac{1}{2}}^{2}\hat{f}_{j+\frac{1}{2}}^{*,[2]}\psi_{j+\frac{1}{2}}^{-}$
$\displaystyle\qquad+r_{j-\frac{1}{2}}^{2}\hat{f}_{j-\frac{1}{2}}^{*,[2]}\psi_{j-\frac{1}{2}}^{+}+s_{j}^{[2],\text{wb}}$
$\displaystyle\quad=\uline{\int_{K_{j}}f^{[2]}(\bm{u}_{h})(\partial_{r}\psi)r^{2}\mathrm{d}r}-\dotuline{r_{j+\frac{1}{2}}^{2}p^{d}\left(r_{j+\frac{1}{2}}\right)\psi_{j+\frac{1}{2}}^{-}}$
$\displaystyle\qquad+\dashuline{r_{j-\frac{1}{2}}^{2}p^{d}\left(r_{j-\frac{1}{2}}\right)\psi_{j-\frac{1}{2}}^{+}}+\uwave{\int_{K_{j}}\left(\frac{2p_{h}}{r}-\rho_{h}\frac{\partial\Phi_{h}}{\partial
r}\right)\psi\,r^{2}\mathrm{d}r}$
$\displaystyle\qquad-\uwave{\int_{K_{j}}\left(\frac{2p^{e}_{h}}{r}-\rho_{h}^{e}\frac{\partial\Phi^{e}_{h}}{\partial
r}\right)\psi\,r^{2}\mathrm{d}r}-\uline{\int_{K_{j}}p^{e}_{h}\,(\partial_{r}\psi)\,r^{2}\mathrm{d}r}$
$\displaystyle\qquad+\dotuline{r_{j+\frac{1}{2}}^{2}p^{d}\left(r_{j+\frac{1}{2}}\right)\psi_{j+\frac{1}{2}}^{-}}-\dashuline{r_{j-\frac{1}{2}}^{2}p^{d}\left(r_{j-\frac{1}{2}}\right)\psi_{j-\frac{1}{2}}^{+}}$
$\displaystyle\quad=0,$ (3.52)

where different underlines are used in the last equality to highlight the
terms that cancel each other. Therefore, we can conclude that the semi-
discrete scheme (3.49) maintains the equilibrium state (2.8) exactly.

###  3.5 The well-balanced total-energy-conserving RKDG scheme

In this subsection, we present the approach to design a total-energy-
conserving fully discrete DG method to ensure the scheme has the total energy
conservation property (2.17) on the discrete level. This will involve two
components: the approximation $s_{j}^{[3]}$ of the source term in the energy
equation (2.3), and the temporal discretization. To illustrate the idea, we
will start with the semi-discrete method to explain the approximation
$s_{j}^{[3]}$, followed by the forward Euler time discretization, and the
high-order Runge–Kutta method at the end.

####  3.5.1 Semi-discrete total-energy-conserving method

The key idea of designing the total-energy-conserving scheme is on the
approximation of the source term in the energy equation (2.3). Let us apply
integration by parts on the source term approximation $s_{j}^{[3]}$ in (3.17),
which leads to

$\displaystyle s_{j}^{[3]}=$ $\displaystyle\int_{K_{j}}-(\rho
u)_{h}\,\frac{\partial\Phi_{h}}{\partial r}\,\delta\,r^{2}\mathrm{d}r$
$\displaystyle=$ $\displaystyle-\left((\rho
u)_{h}\,\Phi_{h}\,\delta\,r^{2}\right)\Big{|}_{r_{j-\frac{1}{2}}^{+}}^{r_{j+\frac{1}{2}}^{-}}+\int_{K_{j}}\frac{\partial}{\partial
r}\left((\rho u)_{h}\,r^{2}\right)\,\Phi_{h}\,\delta\mathrm{d}r$
$\displaystyle+\int_{K_{j}}(\rho
u)_{h}\,\Phi_{h}\,\frac{\partial\delta}{\partial r}\,r^{2}\mathrm{d}r$
$\displaystyle\approx$
$\displaystyle-\left(\hat{f}^{*,[1]}\,\Phi_{h}\,\delta\,r^{2}\right)\Big{|}_{r_{j-\frac{1}{2}}^{+}}^{r_{j+\frac{1}{2}}^{-}}-\int_{K_{j}}\frac{\partial\rho_{h}}{\partial
t}\,\Phi_{h}\,\delta\,r^{2}\mathrm{d}r$ $\displaystyle+\int_{K_{j}}(\rho
u)_{h}\,\Phi_{h}\,\frac{\partial\delta}{\partial r}\,r^{2}\mathrm{d}r$
$\displaystyle:=$ $\displaystyle
s_{j}^{[3],\text{tec}}\left(\bm{u}_{h},\hat{\bm{f}}^{*},\frac{\partial\rho_{h}}{\partial
t},\Phi_{h},\delta\right),$ (3.53)

where the superscript ‘tec’ stands for total-energy-conserving, $\delta$ is
the test function and $\hat{f}^{*,[1]}$ is the first component of the
numerical flux in (3.43). Equation (2.1) is used to replace
$\frac{\partial}{\partial r}\left((\rho u)_{h}\,r^{2}\right)$ by
$-r^{2}\partial\rho_{h}/\partial t$ (approximately).

With this reformulation of the source term, we can now modify the semi-
discrete well-balanced method (3.49) slightly, and obtain the semi-discrete
well-balanced and total-energy-conserving scheme: find
$\bm{u}_{h}\in\bm{\Pi}_{h}$ such that for any test function
$\bm{v}=(\zeta,\psi,\delta)^{T}\in\bm{\Pi}_{h}$, it holds that

$\displaystyle\partial_{t}\int_{K_{j}}\bm{u}_{h}\cdot\bm{v}\,r^{2}\mathrm{d}r$
$\displaystyle=$
$\displaystyle\mathcal{F}_{j}(\bm{u}_{h},\bm{v})+\mathcal{S}^{[2],\text{wb}}_{j}(\bm{u}_{h},\bm{v})+\mathcal{S}_{j}^{[3],\text{tec}}\left(\bm{u}_{h},\hat{\bm{f}}^{*},\frac{\partial\rho_{h}}{\partial
t},\Phi_{h},\delta\right),$ (3.54)

where

$\mathcal{S}^{[2],\text{wb}}_{j}=\left[0,s_{j}^{[2],\text{wb}},0\right]^{T},\qquad\mathcal{S}^{[3],\text{tec}}_{j}=\left[0,0,s_{j}^{[3],\text{tec}}\right]^{T}.$
(3.55)

###### Proposition 3.3

For the semi-discrete scheme (3.54), we have the following total energy
conservation property

$\displaystyle\frac{\partial}{\partial
t}\int_{K_{j}}\left(E_{h}+\frac{1}{2}\rho_{h}\,\Phi_{h}\right)r^{2}\,\mathrm{d}r+\Bigg{(}\hat{f}^{*,[3]}+\hat{f}^{*,[1]}\Phi_{h}$
$\displaystyle\qquad\quad\left.-\frac{1}{8\pi\,G}\left(\Phi_{h}\frac{\partial}{\partial
t}\left(\frac{\partial\Phi_{h}}{\partial
r}\right)-\frac{\partial\Phi_{h}}{\partial t}\frac{\partial\Phi_{h}}{\partial
r}\right)\Bigg{)}r^{2}\,\right|_{r_{j-\frac{1}{2}}^{+}}^{r_{j+\frac{1}{2}}^{-}}=0,$
(3.56)

which is consistent with the continuous result in Eq. (2.15), and leads to the
conservation of total energy
$\int_{\Omega}(E_{h}+\frac{1}{2}\rho_{h}\,\Phi_{h})r^{2}\,\mathrm{d}r$.

###### Proof

Following the approach used in the proof of (2.15), we decompose the first
term into two parts:

$\frac{\partial}{\partial
t}\int_{K_{j}}\left(E_{h}+\frac{1}{2}\rho_{h}\,\Phi_{h}\right)r^{2}\,\mathrm{d}r=\text{I}+\text{II},$
(3.57)

with

I $\displaystyle=\int_{K_{j}}\left(\frac{\partial E_{h}}{\partial
t}+\frac{\partial\rho_{h}}{\partial t}\Phi_{h}\right)r^{2}\,\mathrm{d}r,$
(3.58) II
$\displaystyle=\int_{K_{j}}\frac{1}{2}\left(\rho_{h}\frac{\partial\Phi_{h}}{\partial
t}-\frac{\partial\rho_{h}}{\partial t}\Phi_{h}\right)r^{2}\,\mathrm{d}r.$
(3.59)

We set the test function $\bm{v}$ as $(0,0,1)^{T}$ in (3.54) to obtain

$\displaystyle\int_{K_{j}}\frac{\partial E_{h}}{\partial
t}\,r^{2}\,\mathrm{d}r$
$\displaystyle=-\left.\left(\hat{\bm{f}}^{*,[3]}+\hat{\bm{f}}^{*,[1]}\Phi_{h}\right)r^{2}\,\right|_{r_{j-\frac{1}{2}}^{+}}^{r_{j+\frac{1}{2}}^{-}}-\int_{K_{j}}\frac{\partial\rho_{h}}{\partial
t}\Phi_{h}\,r^{2}\,\mathrm{d}r,$ (3.60)

which leads to the simplification of part I as

$\text{I}=-\left.\left(\hat{\bm{f}}^{*,[3]}+\hat{\bm{f}}^{*,[1]}\Phi_{h}\right)r^{2}\,\right|_{r_{j-\frac{1}{2}}^{+}}^{r_{j+\frac{1}{2}}^{-}}.$
(3.61)

Next, note that the evaluation of $\Phi_{h}$ in (3.7) and (3.8) are exact,
i.e.,

$\displaystyle 4\pi\,G\,\rho_{h}\,r^{2}=\frac{\partial}{\partial
r}\left(\,r^{2}\,\frac{\partial\Phi_{h}}{\partial r}\,\right),$ (3.62)

therefore, following the exact same step in the proof of (2.15) in Section
2.3, we have

II
$\displaystyle=\frac{1}{8\pi\,G}\left.\left(\frac{\partial\Phi_{h}}{\partial
t}\frac{\partial\Phi_{h}}{\partial r}-\Phi_{h}\frac{\partial}{\partial
t}\left(\frac{\partial\Phi_{h}}{\partial
r}\right)\right)r^{2}\,\right|_{r_{j-\frac{1}{2}}^{+}}^{r_{j+\frac{1}{2}}^{-}}.$
(3.63)

The combination of these two equations leads to the total energy conservation
property, which finishes the proof.

####  3.5.2 Forward Euler time discretization and total energy conservation

The extension of the total energy conservation property in Proposition 3.3 to
fully discrete schemes coupled with high-order RK methods is a non-trivial
task. Let us start with the simpler first-order Euler method, and use it as an
example to illustrate how to obtain the fully discrete second- and third-order
total-energy-conserving schemes.

The straightforward application of the forward Euler method to the semi-
discrete well-balanced and total-energy-conserving scheme (3.54) may not
conserve the total energy automatically. The only term that needs extra care
is the approximation of $\mathcal{S}_{j}^{[3],\text{tec}}$ in (3.53), (3.55),
and the fully discrete scheme with forward Euler discretization is given by

$\displaystyle\int_{K_{j}}\bm{u}_{h}^{n+1}\cdot\bm{v}\,r^{2}\mathrm{d}r$
$\displaystyle\quad=\int_{K_{j}}\bm{u}_{h}^{n}\cdot\bm{v}\,r^{2}\mathrm{d}r+\Delta
t\Bigg{(}\mathcal{F}_{j}(\bm{u}_{h}^{n},\bm{v})+\mathcal{S}_{j}^{[2],\text{wb}}(\bm{u}_{h}^{n},\bm{v})$
$\displaystyle\qquad+\mathcal{S}_{j}^{[3],\text{tec}}\left(\bm{u}_{h}^{n},\hat{\bm{f}}^{*,n},\frac{\rho^{n+1}-\rho^{n}}{\Delta
t},\frac{\Phi_{h}^{n+1}+\Phi_{h}^{n}}{2},\delta\right)\Bigg{)}.$ (3.64)

Note that although the right-hand side of (3.64) contains $\rho_{h}^{n+1}$ and
$\Phi_{h}^{n+1}$, the proposed scheme is still an explicit scheme as outlined
below. First we can use the density equation to explicitly evaluate
$\rho_{h}^{n+1}$, and obtain $\Phi_{h}^{n+1}$ following (3.12)-(3.13). Next
the momentum equation is solved to update $(\rho u)_{h}^{n+1}$. Finally, with
the available $\rho_{h}^{n+1}$ and $\Phi_{h}^{n+1}$, we can solve the energy
equation to compute $E_{h}^{n+1}$ explicitly.

###### Proposition 3.4

The fully discrete forward Euler DG scheme (3.64) conserves total energy:

$\int_{\Omega}\left(E_{h}^{n+1}+\frac{1}{2}\rho_{h}^{n+1}\,\Phi_{h}^{n+1}\right)r^{2}\mathrm{d}r=\int_{\Omega}\left(E_{h}^{n}+\frac{1}{2}\rho_{h}^{n}\,\Phi_{h}^{n}\right)r^{2}\mathrm{d}r,$
(3.65)

with outer boundary conditions $\Phi_{h}^{n}(R)=\Phi_{h}^{n+1}(R)=0$ and
$\hat{\bm{f}}_{N+\frac{1}{2}}^{*,n,[3]}=0$.

###### Proof

The main structure of the proof is similar to that of the semi-discrete method
in Proposition 3.3, with more terms due to the temporal discretization. In
each cell $K_{j}$, we take the difference of the total energy in (3.65) and
separate it into two parts:

$\displaystyle\int_{K_{j}}\left(\frac{1}{2}\rho_{h}^{n+1}\,\Phi_{h}^{n+1}-\frac{1}{2}\rho_{h}^{n}\,\Phi_{h}^{n}\right)r^{2}\mathrm{d}r+\int_{K_{j}}\left(E_{h}^{n+1}-E_{h}^{n}\right)r^{2}\mathrm{d}r$
$\displaystyle\qquad:=\text{I}+\text{II},$ (3.66)

with

$\displaystyle\text{I}=\int_{K_{j}}\frac{1}{2}\left(\rho_{h}^{n+1}-\rho_{h}^{n}\right)\left(\Phi_{h}^{n+1}+\Phi_{h}^{n}\right)r^{2}\mathrm{d}r+\int_{K_{j}}\left(E_{h}^{n+1}-E_{h}^{n}\right)r^{2}\mathrm{d}r,$
(3.67)
$\displaystyle\text{II}=\int_{K_{j}}\frac{1}{2}\left(-\rho_{h}^{n+1}\,\Phi_{h}^{n}+\rho_{h}^{n}\Phi_{h}^{n+1}\right)r^{2}\mathrm{d}r.$
(3.68)

Let us introduce the notation:

$\Phi_{h}^{n+\frac{1}{2}}=\frac{\Phi_{h}^{n+1}+\Phi_{h}^{n}}{2}.$ (3.69)

We note that $\hat{\bm{f}}^{*,n}$, $\frac{\partial\Phi_{h}^{n}}{\partial r}$,
and $\Phi_{h}^{n}$ are single-valued in our schemes. By setting the test
function $\bm{v}=(0,0,1)^{T}$ in (3.64), we can derive

$\displaystyle\int_{K_{j}}E_{h}^{n+1}\,r^{2}\mathrm{d}r$
$\displaystyle=\int_{K_{j}}E_{h}^{n}\,r^{2}\mathrm{d}r-\int_{K_{j}}\left(\rho_{h}^{n+1}-\rho_{h}^{n}\right)\Phi_{h}^{n+\frac{1}{2}}\,r^{2}\mathrm{d}r$
$\displaystyle\quad-\Delta
t\left.\left(r^{2}\left(\hat{f}^{*,n,[3]}+\hat{f}^{*,n,[1]}\Phi_{h}^{n+\frac{1}{2}}\right)\right)\right|_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}},$
(3.70)

where $\hat{f}^{*,[i]}$ is the $i$-th component of the numerical flux
$\hat{\bm{f}}^{*}$. Therefore, we can simplify the term I as

$\displaystyle\text{I}=-\Delta
t\left.\left(r^{2}\left(\hat{\bm{f}}^{*,n,[3]}+\hat{\bm{f}}^{*,n,[1]}\Phi_{h}^{n+\frac{1}{2}}\right)\right)\right|_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}}.$
(3.71)

Following the equality (3.62) in the evaluation of $\Phi_{h}$, we have

$\displaystyle
4\pi\,G\int_{K_{j}}\rho_{h}^{n+1}\,\Phi_{h}^{n}\,r^{2}\mathrm{d}r=\int_{K_{j}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial\Phi_{h}^{n+1}}{\partial
r}\right)\Phi_{h}^{n}\mathrm{d}r$ $\displaystyle\hskip
34.14322pt=\left.\left(r^{2}\frac{\partial\Phi_{h}^{n+1}}{\partial
r}\Phi_{h}^{n}\right)\right|_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}}-\int_{K_{j}}\frac{\partial\Phi_{h}^{n+1}}{\partial
r}\frac{\partial\Phi_{h}^{n}}{\partial r}\,r^{2}\mathrm{d}r,$ (3.72)
$\displaystyle
4\pi\,G\,\int_{K_{j}}\rho_{h}^{n}\,\Phi_{h}^{n+1}\,r^{2}\mathrm{d}r=\int_{K_{j}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial\Phi_{h}^{n}}{\partial
r}\right)\Phi_{h}^{n+1}\mathrm{d}r$ $\displaystyle\hskip
34.14322pt=\left.\left(r^{2}\frac{\partial\Phi_{h}^{n}}{\partial
r}\Phi_{h}^{n+1}\right)\right|_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}}-\int_{K_{j}}\frac{\partial\Phi_{h}^{n}}{\partial
r}\frac{\partial\Phi_{h}^{n+1}}{\partial r}\,r^{2}\mathrm{d}r.$ (3.73)

Therefore, we can simplify term II as

$\text{II}=\frac{1}{8\pi\,G}\left.\left(r^{2}\frac{\partial\Phi_{h}^{n}}{\partial
r}\Phi_{h}^{n+1}-r^{2}\frac{\partial\Phi_{h}^{n+1}}{\partial
r}\Phi_{h}^{n}\right)\right|_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}}.$ (3.74)

We combine Eqs. (3.66)-(3.71) and sum over all the cells $K_{j}$ to obtain

$\displaystyle\int_{\Omega}\left(E_{h}^{n+1}+\frac{1}{2}\rho_{h}^{n+1}\,\Phi_{h}^{n+1}\right)r^{2}\mathrm{d}r-\int_{\Omega}\left(E_{h}^{n}+\frac{1}{2}\rho_{h}^{n}\,\Phi_{h}^{n}\right)r^{2}\mathrm{d}r$
$\displaystyle\qquad=\sum_{j=1}^{N}\frac{1}{8\pi\,G}\left.\left(r^{2}\frac{\partial\Phi_{h}^{n}}{\partial
r}\Phi_{h}^{n+1}-r^{2}\frac{\partial\Phi_{h}^{n+1}}{\partial
r}\Phi_{h}^{n}\right)\right|_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}}$
$\displaystyle\qquad\quad-\Delta
t\left.\left(r^{2}\left(\hat{\bm{f}}^{*,n,[3]}+\hat{\bm{f}}^{*,n,[1]}\Phi_{h}^{n+\frac{1}{2}}\right)\right)\right|_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}}$
$\displaystyle\qquad=\left.\frac{1}{8\pi\,G}\left(r^{2}\frac{\partial\Phi_{h}^{n}}{\partial
r}\Phi_{h}^{n+1}-r^{2}\frac{\partial\Phi_{h}^{n+1}}{\partial
r}\Phi_{h}^{n}\right)\right|_{0}^{R}$ $\displaystyle\qquad\quad\left.-\Delta
t\left(r^{2}\left(\hat{\bm{f}}^{*,n,[3]}+\hat{\bm{f}}^{*,n,[1]}\Phi_{h}^{n+\frac{1}{2}}\right)\right)\right|_{0}^{R}$
$\displaystyle\qquad=~{}0,$ (3.75)

where the last equality is due to the outer boundary condition
$\Phi_{h}^{n}(R)=\Phi_{h}^{n+1}(R)=\Phi_{h}^{n+\frac{1}{2}}(R)=0$ and
$\hat{f}_{N+\frac{1}{2}}^{*,n,[3]}=0$. Therefore, the fully discrete forward
Euler DG scheme (3.64) has the total energy conservation property.

###### Remark 3.5

The assumptions on the outer boundary condition (i.e.,
$\Phi_{h}^{n}(R)=\Phi_{h}^{n+1}(R)=0$ and
$\hat{\bm{f}}_{N+\frac{1}{2}}^{*,n,[3]}=0$) are only used in the last equality
of the proof. We use these assumptions for ease of presentation. The total
energy conservation property of our numerical methods does not depend on these
assumptions. In Section 4.4, we consider a numerical example without the
assumption $\hat{\bm{f}}_{N+\frac{1}{2}}^{*,n,[3]}=0$, and observe
conservation of total energy, after adding correction terms due to the outer
boundary. We can deal with the case without the assumption
$\Phi_{h}^{n}(R)=\Phi_{h}^{n+1}(R)=0$ in a similar way by adding correction
term. We refer to Section 4.4 for the details on these correction terms and
the numerical observation.

###### Remark 3.6

We note that our proposed scheme (3.64) still has the well-balanced property.
The only thing to check is that the source term approximation
$\mathcal{S}_{j}^{[3],\text{\rm tec}}=0$ holds at the steady state. This holds
due to the fact that $\hat{f}^{*,n,[1]}=0$, $u_{h}^{n}=0$, and also
$\rho_{h}^{n}=\rho_{h}^{n+1}$ by updating the density equation with the well-
balanced DG method at the steady state.

####  3.5.3 High-order Runge-Kutta time discretization

In this section, we will extend the well-balanced and total-energy-conserving
method (3.64) coupled with forward Euler discretization to high-order RK
discretization. In Mullen et al. (2021), the fully discrete energy conserving
schemes with second- and third-order RK time discretization are introduced in
the context of finite difference methods. The key idea is to use different
source term approximations for each stage of the Runge–Kutta method, and a
similar idea will be explored here. Comparing with the RK methods in Mullen et
al. (2021) and this paper, the main difference is that we involve additional
terms, such as the approximation of $\frac{\partial\rho}{\partial t}$. This is
because our DG schemes include test functions and the relationship between the
variables $\bm{u}$ is more complicated.

Let us start with the second-order RK method. For the differential equation of
the general form $w_{t}=\mathcal{L}(w)$, a second-order RK method can be
formulated as

$\displaystyle w^{(1)}$ $\displaystyle=w^{n}+\Delta t\,\mathcal{L}(w^{n}),$
$\displaystyle w^{n+1}$ $\displaystyle=w^{n}+\frac{1}{2}\left(w^{(1)}+\Delta
t\,\mathcal{L}(w^{(1)})\right)$ $\displaystyle=w^{n}+\Delta
t\left(\frac{\mathcal{L}(w^{n})+\mathcal{L}(w^{(1)})}{2}\right).$ (3.76)

Starting from the forward Euler method (3.64), the fully discrete total-energy
conserving scheme with second-order RK method (3.76) is given by

$\displaystyle\int_{K_{j}}\bm{u}_{h}^{(1)}\cdot\bm{v}\,r^{2}\mathrm{d}r$
$\displaystyle\qquad=\int_{K_{j}}\bm{u}_{h}^{n}\cdot\bm{v}\,r^{2}\mathrm{d}r+\Delta
t\Bigg{(}\mathcal{F}_{j}(\bm{u}_{h}^{n},\bm{v})+\mathcal{S}_{j}^{[2],\text{wb}}(\bm{u}_{h}^{n},\bm{v})$
$\displaystyle\qquad\quad+\mathcal{S}_{j}^{[3],\text{tec}}\left(\bm{u}_{h}^{n},\hat{\bm{f}}^{*,n},\frac{\rho^{(1)}-\rho^{n}}{\Delta
t},\Phi_{h}^{(0,1)},\delta\right)\Bigg{)},$ (3.77)
$\displaystyle\int_{K_{j}}\bm{u}_{h}^{n+1}\cdot\bm{v}\,r^{2}\mathrm{d}r$
$\displaystyle\qquad=\int_{K_{j}}\bm{u}_{h}^{n}\cdot\bm{v}\,r^{2}\mathrm{d}r+\Delta
t\Bigg{(}\frac{\mathcal{F}_{j}(\bm{u}_{h}^{n},\bm{v})+\mathcal{F}_{j}(\bm{u}_{h}^{(1)},\bm{v})}{2}$
$\displaystyle\qquad\quad+\frac{\mathcal{S}_{j}^{[2],\text{wb}}(\bm{u}_{h}^{n},\bm{v})+\mathcal{S}_{j}^{[2],\text{wb}}(\bm{u}_{h}^{(1)},\bm{v})}{2}$
$\displaystyle\qquad\quad+\mathcal{S}_{j}^{[3],\text{tec}}\left(\bm{u}_{h}^{(0,1)},\hat{\bm{f}}^{*,(0,1)},\frac{\rho^{n+1}-\rho^{n}}{\Delta
t},\Phi_{h}^{(0,2)},\delta\right)\Bigg{)},$ (3.78)

where we introduced the following notations

$\displaystyle\hat{\bm{f}}^{*,(0,1)}=\frac{1}{2}\left(\hat{\bm{f}}^{*,n}+\hat{\bm{f}}^{*,(1)}\right),\qquad\bm{u}^{(0,1)}=\frac{1}{2}\left(\bm{u}_{h}^{n}+\bm{u}_{h}^{(1)}\right),$
$\displaystyle\Phi_{h}^{(0,1)}=\frac{1}{2}\left(\Phi_{h}^{n}+\Phi_{h}^{(1)}\right),\qquad\Phi_{h}^{(0,2)}=\frac{1}{2}\left(\Phi_{h}^{n}+\Phi_{h}^{n+1}\right).$
(3.79)

The third-order strong-stability-preserving RK method for
$w_{t}=\mathcal{L}(w)$ can be formulated as

$\displaystyle w^{(1)}$ $\displaystyle=w^{n}+\Delta t\,\mathcal{L}(w^{n}),$
$\displaystyle w^{(2)}$
$\displaystyle=\frac{3}{4}w^{n}+\frac{1}{4}\left(w^{(1)}+\Delta
t\,\mathcal{L}(w^{(1)})\right)$ $\displaystyle=w^{n}+\frac{\Delta
t}{2}\left(\frac{\mathcal{L}(w^{n})+\mathcal{L}(w^{(1)})}{2}\right),$
$\displaystyle w^{n+1}$
$\displaystyle=\frac{1}{3}w^{n}+\frac{2}{3}\left(w^{(2)}+\Delta
t\,\mathcal{L}(w^{(2)})\right)$ $\displaystyle=w^{n}+\Delta
t\left(\frac{\mathcal{L}(w^{n})+\mathcal{L}(w^{(1)})+4\mathcal{L}(w^{(2)})}{6}\right).$
(3.80)

The fully discrete total-energy conserving scheme with this third-order RK
method is given by

$\displaystyle\int_{K_{j}}\bm{u}_{h}^{(1)}\cdot\bm{v}\,r^{2}\mathrm{d}r=\int_{K_{j}}\bm{u}_{h}^{n}\cdot\bm{v}\,r^{2}\mathrm{d}r$
$\displaystyle\qquad\quad+\Delta
t\Bigg{(}\mathcal{F}_{j}(\bm{u}_{h}^{n},\bm{v})+\mathcal{S}_{j}^{[2],\text{wb}}(\bm{u}_{h}^{n},\bm{v})$
$\displaystyle\qquad\quad+\mathcal{S}_{j}^{[3],\text{tec}}\left(\bm{u}_{h}^{n},\hat{\bm{f}}^{*,n},\frac{\rho^{(1)}-\rho^{n}}{\Delta
t},\Phi_{h}^{n+\frac{1}{2}},\delta\right)\Bigg{)},$ (3.81)
$\displaystyle\int_{K_{j}}\bm{u}_{h}^{(2)}\cdot\bm{v}\,r^{2}\mathrm{d}r=\int_{K_{j}}\bm{u}_{h}^{n}\cdot\bm{v}\,r^{2}\mathrm{d}r$
$\displaystyle\qquad\quad+\frac{\Delta
t}{2}\Bigg{(}\frac{\mathcal{F}_{j}(\bm{u}_{h}^{n},\bm{v})+\mathcal{F}_{j}(\bm{u}_{h}^{(1)},\bm{v})}{2}$
$\displaystyle\qquad\quad+\frac{\mathcal{S}_{j}^{[2],\text{wb}}(\bm{u}_{h}^{n},\bm{v})+\mathcal{S}_{j}^{[2],\text{wb}}(\bm{u}_{h}^{(1)},\bm{v})}{2}$
$\displaystyle\qquad\quad+\mathcal{S}_{j}^{[3],\text{tec}}\left(\bm{u}_{h}^{(0,1)},\hat{\bm{f}}^{*,(0,1)},\frac{\rho^{(2)}-\rho^{n}}{\Delta
t/2},\Phi_{h}^{(0,2)},\delta\right)\Bigg{)},$ (3.82)
$\displaystyle\int_{K_{j}}\bm{u}_{h}^{n+1}\cdot\bm{v}\,r^{2}\mathrm{d}r=\int_{K_{j}}\bm{u}_{h}^{n}\cdot\bm{v}\,r^{2}\mathrm{d}r$
$\displaystyle\qquad\quad+\Delta
t\Bigg{(}\frac{\mathcal{F}_{j}(\bm{u}_{h}^{n},\bm{v})+\mathcal{F}_{j}(\bm{u}_{h}^{(1)},\bm{v})+4\mathcal{F}_{j}(\bm{u}_{h}^{(2)},\bm{v})}{6}$
$\displaystyle\qquad\quad+\frac{\mathcal{S}_{j}^{[2],\text{wb}}(\bm{u}_{h}^{n},\bm{v})+\mathcal{S}_{j}^{[2],\text{wb}}(\bm{u}_{h}^{(1)},\bm{v})+4\mathcal{S}_{j}^{[2],\text{wb}}(\bm{u}_{h}^{(2)},\bm{v})}{6}$
$\displaystyle\qquad\quad+\mathcal{S}_{j}^{[3],\text{tec}}\left(\bm{u}_{h}^{(0,2)},\hat{\bm{f}}^{*,(0,2)},\frac{\rho^{n+1}-\rho^{n}}{\Delta
t},\Phi_{h}^{(0,3)},\delta\right)\Bigg{)},$ (3.83)

with the following notations

$\displaystyle\hat{\bm{f}}^{*,(0,2)}=\frac{1}{6}\left(\hat{\bm{f}}^{*,n}+\hat{\bm{f}}^{*,(1)}+4\hat{\bm{f}}^{*,(2)}\right),$
(3.84)
$\displaystyle\bm{u}^{(0,2)}=\frac{1}{6}\left(\bm{u}_{h}^{n}+\bm{u}_{h}^{(1)}+4\bm{u}_{h}^{(2)}\right),\qquad\Phi_{h}^{(0,3)}=\frac{1}{2}\left(\Phi_{h}^{n}+\Phi_{h}^{n+1}\right).$

Note that different source term approximations of
$\mathcal{S}_{j}^{[3],\text{tec}}$ are employed in the each stage of the RK
method, in order to simultaneously achieve the total energy conservation
property and high-order accuracy. The proofs of the well-balanced property and
total energy conservation of the high-order RKDG methods (3.77)-(3.78) and
(3.81)-(3.83) follow the exact same approach as that of the forward Euler DG
scheme (3.64), and is omitted here to save space.

###  3.6 TVB limiter

For problems containing strong discontinuities, oscillations may develop in
the solutions obtained with DG methods, and in this case nonlinear limiters
are needed after each stage of the Runge–Kutta methods to control these
oscillations. One popular choice is the total variation bounded (TVB) limiter
Cockburn & Shu (1989). Its extension to the system in spherically symmetrical
coordinates has been considered in Pochik et al. (2021), and will be employed
here, provided some modifications to ensure the total-energy-conserving
property.

We start by defining two different cell averages of $\bm{u}_{h}$ in cell
$K_{j}$: the standard and weighted cell averages given by

$\bar{\bm{u}}_{j}=\frac{\int_{K_{j}}\bm{u}_{h}\,\mathrm{d}r}{\int_{K_{j}}1\,\mathrm{d}r},\qquad\qquad\tilde{\bm{u}}_{j}=\frac{\int_{K_{j}}\bm{u}_{h}\,r^{2}\,\mathrm{d}r}{\int_{K_{j}}r^{2}\,\mathrm{d}r},$
(3.85)

respectively. In cell $K_{j}$, the forward and backward slopes are defined as

$\displaystyle\Delta\bm{u}_{j}^{F}=\frac{\bar{\bm{u}}_{j+1}-\bar{\bm{u}}_{j}}{r_{j+1}-r_{j}},$
$\displaystyle\qquad\qquad\Delta\bm{u}_{j}^{B}=\frac{\bar{\bm{u}}_{j}-\bar{\bm{u}}_{j-1}}{r_{j}-r_{j-1}},$
(3.86)

where $r_{j}=(r_{j+\frac{1}{2}}+r_{j-\frac{1}{2}})/2$ denotes the midpoint of
$K_{j}$. Then we apply the minmod function in Cockburn & Shu (1989) to obtain

$\displaystyle\tilde{\Delta}\bm{u}_{j}=\text{minmod}\left(\Delta\bm{u}_{j},~{}\beta\Delta\bm{u}_{j}^{F},~{}\beta\Delta\bm{u}_{j}^{B}\right),$
(3.87)

where

$\Delta\bm{u}_{j}=\frac{\bm{u}_{h,j+\frac{1}{2}}^{-}-\bm{u}_{h,j-\frac{1}{2}}^{+}}{r_{j+\frac{1}{2}}-r_{j-\frac{1}{2}}},$
(3.88)

with $\beta$ being a constant to be specified. In Pochik et al. (2021), it was
shown that $\beta=1.75$ yields good results for a range of problems, and this
value will also be used in this paper. If $\tilde{\Delta}\bm{u}_{j}$ and
$\Delta\bm{u}_{j}$ are the same, this indicates that a limiter is not needed
in this cell. When they are different, we mark this cell $K_{j}$ as a troubled
cell. In such cell, we define a new linear polynomial $\tilde{\bm{u}}_{h,j}$
as

$\tilde{\bm{u}}_{h,j}=\tilde{\bm{u}}_{j}^{0}+\tilde{\Delta}\bm{u}_{j}(r-r_{j}),\qquad\tilde{\bm{u}}_{j}^{0}=\tilde{\bm{u}}_{j}-\tilde{\Delta}\bm{u}_{j}\frac{\int_{K_{j}}(r-r_{j})\,r^{2}\,\mathrm{d}r}{\int_{K_{j}}r^{2}\,\mathrm{d}r},$
(3.89)

which has the updated slope $\tilde{\Delta}\bm{u}_{j}$ while keeping the same
weighted cell average as $\tilde{\bm{u}}_{j}$. In the cells which are not
marked as troubled cells, we simply set $\tilde{\bm{u}}_{h,j}=\bm{u}_{h,j}$.
Finally, we replace the solution $\bm{u}_{h}$ by the updated solution
$\tilde{\bm{u}}_{h}$ and continue the computation with the updated solution.
This finishes the TVB limiter procedure. One can easily verify that the
weighted cell average of $\tilde{\bm{u}}_{h,j}$ are the same as $\bm{u}_{h}$
in each computational cell, which yields the mass conservation property of the
limiter procedure.

Since the total energy depends nonlinearly on the variable $\rho_{h}$, this
TVB limiter may destroy the total energy conservation property, which is
satisfied by the proposed fully discrete method. To ensure the total-energy-
conserving property, we slightly modify the TVB limiter on the variable
$E_{h}$ as outlined below. Since the Euler–Poisson system does not conserve
the non-gravitational energy $E$ in the PDE level, we propose an additional
correction of $\tilde{E}_{h,j}$ as follows

$\tilde{\tilde{E}}_{h,j}=\tilde{E}_{h,j}+\frac{\int_{K_{j}}\frac{1}{2}(\rho_{h}\phi_{h}-\tilde{\rho}_{h}\tilde{\Phi}_{h})\,r^{2}\,\mathrm{d}r}{\int_{K_{j}}r^{2}\,\mathrm{d}r},$
(3.90)

to ensure that the total energy
$\int_{K_{j}}(E_{h}+\frac{1}{2}\rho_{h}\Phi_{h})~{}dr$ is not changed by the
limiting procedure. Here $\tilde{\tilde{E}}_{h,j}$ is the updated numerical
solution of $E$, $\tilde{E}_{h,j}$ is obtained in (3.89), $\rho_{h}$ is the
numerical solution before limiting, $\tilde{\rho}_{h}$ is the numerical
solution after limiting, $\Phi_{h}$ and $\tilde{\Phi}_{h}$ are the
gravitational potential calculated from $\rho_{h}$ and $\tilde{\rho}_{h}$
respectively. Note that $\tilde{\Phi}_{h}$ is evaluated after
$\tilde{\rho}_{h}$ is available in all the cells, hence even though a cell
$K_{j}$ is not marked as troubled cell, the value of $\tilde{\Phi}_{h}$ in
this cell may be different from the original $\Phi_{h}$ due to modified
$\tilde{\rho}_{h}$ in troubled cells in other locations. Therefore, this
correction (3.90) will be applied for every cell regardless of being marked as
troubled cells or not.

The procedure of applying TVB limiter in each stage of Runge-Kutta method is
summarized below, where the forward Euler time discretization is used for ease
of presentation.

  1. 1.

At each time level $t^{n}$ (or every intermediate stage of Runge-Kutta
method), compute $\rho_{h}^{n+1},(\rho u)_{h}^{n+1}$ for all cells $K_{j}$;

  2. 2.

Apply the TVB limiter to obtain $\tilde{\rho}_{h}^{n+1},\widetilde{\rho
u}^{n+1}$;

  3. 3.

Evaluate $\tilde{\Phi}_{h}^{n+1}$ based on the limited
$\tilde{\rho}_{h}^{n+1}$;

  4. 4.

Compute $E_{h}^{n+1}$ (which employs the limited $\tilde{\rho}_{h}^{n+1}$ and
$\tilde{\Phi}_{h}^{n+1}$) and apply TVB limiter with total-energy-conserving
correction to $\tilde{\tilde{E}}_{h}^{n+1}$ (which involves both
$\rho_{h}^{n+1}$, $\Phi_{h}^{n+1}$ and $\tilde{\rho}_{h}^{n+1}$,
$\tilde{\Phi}_{h}^{n+1}$).

###### Remark 3.7

For the purpose of the well-balanced property, we use
$\bm{u}_{h}-\bm{u}_{h}^{e}$ instead of $\bm{u}_{h}$ as an indicator to
identify the troubled cells Xing (2014). If a cell is marked as a troubled
cell, the update procedure is still applied on $\bm{u}_{h}$ as mentioned
above. In the steady state, we have $\bm{u}_{h}-\bm{u}_{h}^{e}=0$, hence the
TVB limiter will not take effect, and the well-balanced property will not be
affected by the limiter.

##  4 Numerical examples

In this section, numerical examples will be provided to verify the properties
of our proposed scheme, including the well-balanced property, total energy
conservation properties and high-order accuracy. We use $P^{2}$ piecewise
polynomial in the DG method and the third-order RK method (3.81)-(3.83) in the
numerical tests, unless otherwise stated. The CFL number is set as 0.16 to
determine the time step size.

###  4.1 Well-balanced and small perturbation tests

In this example, we consider a simple polytropic equilibrium and verify that
our proposed scheme has the well-balanced property to maintain this
equilibrium up to round-off error. We set $G=1/(4\pi)$ in this example, and
choose two cases, $\gamma=2$ and $\gamma=1.2$, along with $\rho_{0}=1$ and
$\kappa=1$. We have the following initial data

$\rho(r,0)=\frac{\sqrt{2}\sin(\frac{r}{\sqrt{2}})}{r},\quad\rho u(r,0)=0,\quad
p(r,0)=\frac{2\sin^{2}(\frac{r}{\sqrt{2}})}{r^{2}},$ (4.1)

if $\gamma=2$, and

$\rho(r,0)=(1+\frac{1}{18}r^{2})^{-2.5},\quad\rho u(r,0)=0,\quad
p(r,0)=(1+\frac{1}{18}r^{2})^{-3},$ (4.2)

if $\gamma=1.2$, on the domain $\Omega=[0,1]$. The reflecting boundary
condition is considered for the inner boundary and we set
$\bm{u}^{+}(1)=\bm{u}^{-}(1)$ at the outer boundary. We set the stopping time
$t=4$ on the mesh with 200 uniform cells, and present the $L^{1}$ errors of
the numerical solutions in Table 1, where both single and double precisions
have been considered in the simulation. We can see that errors stay at the
level of round-off errors for different precision, which verify the desired
well-balanced property.

Table 1: Example 4.1, $L^{1}$ error of the numerical solutions for different precision in the well-balanced test. Case | Precision | $\rho$ | $\rho u$ | $E$
---|---|---|---|---
$\gamma=2$ | double | 3.89E-13 | 2.70E-15 | 6.52E-14
quad | 3.55E-31 | 3.44E-33 | 5.94E-32
$\gamma=1.2$ | double | 6.75E-13 | 8.00E-15 | 6.31E-13
quad | 6.04E-31 | 8.00E-33 | 5.74E-31

Next, we show the advantage of our proposed scheme in capturing a small
perturbation to the equilibrium state. The initial data is given by imposing a
pressure perturbation to the $\gamma=2$ equilibrium

$\displaystyle\rho(r,0)=\frac{\sqrt{2}\sin(\frac{r}{\sqrt{2}})}{r},\quad\rho
u(r,0)=0,$ $\displaystyle
p(r,0)=\frac{2\sin^{2}(\frac{r}{\sqrt{2}})}{r^{2}}+A\exp(-100r^{2}),$ (4.3)

on the domain $\Omega=[0,0.5]$. The pressure is perturbed by a Gaussian bump
of amplitude $A=10^{-6}$ in this test. We compute the solutions until $t=0.2$.
A reference solution is computed with $N=400$ for comparison. We plot the
velocity and pressure perturbation for $N=100$ in Figure 1, compared with the
numerical solution of the non-well-balanced DG scheme from Section 3.3, and
the reference solution. From the figures, we can see that the well-balanced
scheme resolves the perturbation much better on a relatively coarse mesh.
Similar test under the framework of finite difference methods in three
dimensions can also be found in Käppeli & Mishra (2014).

(a) pressure perturbation of wb

(b) pressure perturbation of non-wb

(c) velocity of wb

(d) velocity of non-wb

Figure 1: Example 4.1, numerical results at time $t=0.2$ for the small
perturbation test. “wb” denotes the proposed DG scheme and “non-wb” denotes
the standard DG scheme. The wb result is compared with non-wb result and the
reference solution.

###  4.2 Accuracy test

  1. 1.

The accuracy test near the equilibrium state.

In this example, we test the accuracy of the numerical solution near the
equilibrium state and use the same initial condition in (4.3) with parameter
$A=0.001$. We set the domain $\Omega=[0,0.5]$, polynomial degree $k=2$ and
stopping time $t=0.2$, same as those in Section 4.1. Since the exact solution
is unknown, we use the numerical solution of $N=640$ as a reference solution.
The error table are shown in Table 2. We can observe the optimal convergence
rate for all the variables. In addition, we also list the errors of the
standard DG scheme (3.16) in Table 3 for comparison. We observe that although
both schemes have the optimal convergence order, the errors of our proposed
scheme are much smaller than those of the standard scheme.

Table 2: Example 4.2, accuracy test near the equilibrium state for $k=2$ with our proposed third-order RKDG scheme (3.81)-(3.83). $N$ | $\rho$ | $\rho u$ | $E$
---|---|---|---
10 | 2.62E-07 | - | 1.63E-07 | - | 2.23E-07 | -
20 | 3.09E-08 | 3.08 | 1.71E-08 | 3.25 | 2.41E-08 | 3.21
40 | 3.73E-09 | 3.05 | 2.16E-09 | 2.98 | 3.08E-09 | 2.97
80 | 4.48E-10 | 3.06 | 2.97E-10 | 2.86 | 4.24E-10 | 2.86
Table 3: Example 4.2, accuracy test near the equilibrium state for $k=2$ with the standard DG scheme (3.16) and third-order RKDG time discretization (3.80) $N$ | $\rho$ | $\rho u$ | $E$
---|---|---|---
10 | 1.84E-04 | - | 1.48E-04 | - | 2.19E-04 | -
20 | 2.62E-05 | 2.81 | 2.03E-05 | 2.87 | 2.16E-05 | 3.34
40 | 3.35E-06 | 2.97 | 2.56E-06 | 2.99 | 3.96E-06 | 2.87
80 | 4.34E-07 | 2.95 | 3.33E-07 | 2.94 | 4.25E-07 | 2.80

  1. (ii)

The accuracy test far away from the equilibrium state.

In this example, we provide an accuracy test for solutions far away from the
equilibrium state, to test the high-order convergence rate of the DG methods.
We consider the following “manufactured” exact solutions

$\rho(r,t)=\frac{\exp(t-r)}{r^{2}},\quad u(r,t)=1,\quad
p(r,t)=\frac{1}{r^{2}}.$ (4.4)

As a result, the Euler–Poisson equations (3.14) becomes

$\frac{\partial\bm{u}}{\partial t}+\frac{1}{r^{2}}\frac{\partial}{\partial
r}(r^{2}\bm{f}(\bm{u}))=\bm{s}(\bm{u},\Phi)+\bm{w}(r),$ (4.5)

with an additional source term $\bm{w}(r)$ given by

$\bm{w}(r)=\left(0,-\frac{\exp(2(t-r))+2r}{r^{4}},-\frac{\exp(2(t-r))}{r^{4}}\right)^{T}.$
(4.6)

In this test, we set $\gamma=2$, $G=1/(4\pi)$, the computational domain is
$\Omega=[0.5,1]$, and the stopping time is set to $t=0.1$. The exact solution
is used to provide the boundary condition for the Euler equations, and the
boundary condition for the Poisson equation is set as

$\frac{\partial\Phi_{h}}{\partial r}(0.5)=-4\exp(t-0.5),\quad\Phi_{h}(0.5)=0.$
(4.7)

Since our computational domain does not contain the origin $r=0$, our approach
of recovering the reference equilibrium state $\bm{u}^{d}$ needs an additional
boundary condition instead of (2.14). For simplicity, we skip the steps of
recovering the reference state in Section 3.4.2 and set a global steady state
$\bm{u}^{d}$ explicitly for all cells without using (3.38):

$\rho^{d}(r)=\frac{\sqrt{2}\sin(\frac{r}{\sqrt{2}})}{r},\quad u^{d}(r)=0,\quad
p^{d}(r)=\frac{2\sin^{2}(\frac{r}{\sqrt{2}})}{r^{2}}.$ (4.8)

We have performed the simulations for various mesh size $N$. The results for
$k=1$ with the second-order RKDG scheme (3.77)-(3.78) and $k=2$ with the
third-order RKDG scheme (3.81)-(3.83) are shown in Table 4. We can observe the
optimal convergence rate for all the variables and $k=1,2$, which confirms the
high-order accuracy of the proposed RKDG method. More specifically, the
different source term approximations in each stage of the third-order RK
method (3.81)-(3.83) yields the desired third-order accuracy.

Table 4: Example 4.2, accuracy test far away from the equilibrium state for $k=1,2$ with equations (4.5). Case | $N$ | $\rho$ | $\rho u$ | $E$
---|---|---|---|---
$k=1$ | 25 | 4.12E-04 | - | 5.17E-04 | - | 6.46E-04 | -
50 | 1.04E-04 | 1.98 | 1.31E-04 | 1.98 | 1.63E-04 | 1.99
100 | 2.63E-05 | 1.99 | 3.29E-05 | 1.99 | 4.10E-05 | 1.99
200 | 6.60E-06 | 1.99 | 8.59E-06 | 2.00 | 1.03E-05 | 2.00
$k=2$ | 25 | 1.29E-05 | - | 1.75E-05 | - | 9.69E-06 | -
50 | 1.82E-06 | 2.82 | 2.41E-06 | 2.86 | 1.33E-06 | 2.87
100 | 2.44E-07 | 2.90 | 3.17E-07 | 2.92 | 1.75E-07 | 2.92
200 | 3.16E-08 | 2.95 | 4.08E-08 | 2.96 | 2.25E-08 | 2.96

###  4.3 Explosion

In this example, we validate the shock capturing and total energy conservation
properties of our proposed scheme. The initial data is given by

$\displaystyle\rho(r,0)=\frac{\sin(\sqrt{2\pi/\kappa}r)}{\sqrt{2\pi/\kappa}r},\quad\rho
u(r,0)=0,$ $\displaystyle p(r,0)=\begin{cases}\alpha\kappa\rho(r,0)^{2},&r\leq
r_{1}\cr\kappa\rho(r,0)^{2},&r>r_{1}\end{cases},$ (4.9)

where we set $\kappa=1$, $\gamma=2$, $G=1$ and increase the equilibrium
pressure by a factor $\alpha=10$ for $r\leq r_{1}=0.1$. The computational
domain is set as $\Omega=[0,0.5]$, and discretized with $N=200$ cells. We use
$P^{2}$ piecewise polynomial and the third-order RK method (3.81)-(3.83). We
set the boundary condition of the velocity $u(0.5,t)=0$ at the outer domain
boundary. We perform the simulation up to time $t=0.15$, and the numerical
results are shown in Figure 2. Both the well-balanced scheme and the standard
DG scheme perform similarly in capturing shocks, which means our proposed
scheme does not diminish the robustness of the shock capturing capability.
Moreover, we can observe that our proposed scheme conserves total energy up to
machine precision, while the standard DG scheme produces an error of about
$3.5\times 10^{-6}$ at $t=0.15$.

Figure 2: The solution of well-balanced scheme (blue) and standard DG scheme
(red) by using $N=200$ cells, compared with the reference solution (black)
produced with $N=800$ cells. From left to right: the numerical solutions of
density, velocity, pressure at time $t=0.15$ and the time history of the
changes in total energy. The maximum absolute value of the changes in total
energy is $8.049\times 10^{-15}$ for the proposed scheme.

###  4.4 Yahil-Lattimer collapse

In this section, we consider the Yahil-Lattimer collapse test, which involves
self-gravity and was studied in Endeve et al. (2019), using standard DG
methods. It models the self-similar collapse of a polytropic star, i.e.
$p=\kappa\rho^{\gamma}$. In Yahil (1983), self-similar solutions to the
gravitational collapse problem were constructed for $6/5\leq\gamma<4/3$. With
two dimensional parameters in the model (the gravitational constant $G$ and
the polytropic constant $\kappa$), the dimensionless similarity variable is

$X=\kappa^{-\frac{1}{2}}G^{(\gamma-1)/2}r(-t)^{\gamma-2},$ (4.10)

where the origin of time is the moment of infinite central density. All the
hydrodynamic variables can be expressed as a function of $X$, and the time-
dependent Euler equations can be recast as a system of ODEs (see Yahil, 1983,
for details). Therefore, we use these self-similar solutions solved by the
ODEs given in Yahil (1983) as a reference solution.

We show some numerical results obtained with $\gamma=1.3$. We set the
computational domain to $\Omega=[0,10^{10}]$ cm discretized with $N=256$
cells, and the collapse time to $(-t)=150$ ms. We use a geometrically
increasing cell spacing

$\Delta r_{j}=r_{j+\frac{1}{2}}-r_{j-\frac{1}{2}}=a^{j-1}\Delta r_{1},\qquad
j=1,...,N,$ (4.11)

with the size of the innermost cell set to $\Delta r_{1}=1\times 10^{5}$ cm,
and increasing at a rate $a=1.03203$. The size of the last element is about
$3\times 10^{8}$ cm. The gravitational constant $G$ is set to $6.67430\times
10^{-8}$ $\mbox{cm}^{-3}\ \mbox{g}^{-1}\ \mbox{s}^{-2}$. We use the reference
solution at time $(-t)=150$ ms to compute the initial density and velocity.
The polytropic constant $\kappa=9.54\times 10^{14}$ is used to give the
initial pressure. We use the reflecting boundary condition for the inner
boundary and zeroth-order extrapolation for the outer boundary.

We simulate collapse until $(-t)=0.5$ ms, and the central density increases
from about $10^{9}$ g $\mbox{cm}^{-3}$ to about $10^{14}$ g $\mbox{cm}^{-3}$.
We plot the density $\rho$ and velocity $u$ at different times in Figure 3,
and compare the results with the reference solutions obtained in Yahil (1983).
The figures show that our numerical method performs well during collapse. We
also compare the total energy conservation property between our proposed
scheme and the standard DG scheme. The total energy is defined as
$E_{tot}=\int_{\Omega}\left(E+\frac{1}{2}\,\rho\,\Phi\right)\,r^{2}\,\mathrm{d}r$.
The total energy conservation for RK3 time discretization $\Delta E$ is
defined as follows

$\displaystyle\Delta E(t^{m+1})=$ $\displaystyle
E_{tot}(t^{m+1})-E_{tot}(t^{m})$ $\displaystyle+4\pi\Delta
t\,R^{2}\frac{\hat{\bm{f}}^{n,[3]}_{N+\frac{1}{2}}+\bm{f}^{(1),[3]}_{N+\frac{1}{2}}+4\bm{f}^{(2),[3]}_{N+\frac{1}{2}}}{6},$
$\displaystyle\Delta E=$ $\displaystyle\sum_{m=1}^{M}\Delta E(t^{m+1}),$
(4.12)

where $R$ is the outer boundary, $N$ is the number of cells and $M$ is the
number of time steps. When the time is close to $(-t)=0.5$ ms and the density
grow rapidly to $10^{14}$ g $\mbox{cm}^{-3}$, our proposed scheme maintains
total energy conservation to round-off error while that of the standard scheme
is much larger.

Figure 3: Example 4.4, the figure of numerical solution (blue) of density
$\rho$ (top left) and velocity $u$ (top right) during collapse, compared with
the standard scheme (red) and the reference solution (black). We compared the
solutions at select central densities, approximately
$[10^{10},10^{11},10^{12},10^{13},10^{14}]$ g cm-3, which correspond to
$(-t)=[51.0,15.0,5.0,1.5,0.5]$ ms. Velocity gradually decreases over time. The
comparison of the total energy conservation between our proposed scheme and
standard scheme versus central density shows in the bottom that when the time
is close to $(-t)=0.5$ ms, our proposed scheme has a much smaller total energy
conservation than the standard scheme.

###  4.5 Toy model of stellar core-collapse, bounce, and shock evolution

We consider a toy model of core-collapse supernova as considered in Janka et
al. (1993); Käppeli & Mishra (2016). This test simulates the spherically
symmetric and adiabatic collapse, bounce, shock evolution, and proto-neutron
star formation for a simplified model using a phenomenological EoS. This test
provides a stringent check on the energy conservation properties of our
proposed scheme — especially during core bounce when core-collapse supernova
codes typically exhibit an abrupt change in the total energy (e.g., Skinner et
al., 2019; Bruenn et al., 2020).

The governing equations are given by (2.1)-(2.3) and (2.5) with a non-ideal
EoS. We first set $\gamma=4/3$ and obtain an equilibrium state according to
(2.6) and (2.7) for a central density $\rho_{c}=10^{10}\,\text{g/cm}^{3}$,
polytropic constant $\kappa=4.897\times 10^{14}$ (in cgs units), and
gravitational constant $G=6.67430\times
10^{-8}\,\text{cm}^{-3}\text{g}^{-1}\text{s}^{-2}$. We initialize the collapse
by reducing the adiabatic index from $\gamma=4/3$ to a slightly smaller value
$\gamma_{1}=1.325$. Then the initial internal energy density is set as $\rho
e=\kappa\rho^{\gamma_{1}}/(\gamma_{1}-1)$ where the initial density $\rho$ is
the equilibrium density for $\gamma=\frac{4}{3}$ and the initial momentum is
set to zero.

The EoS in this test consists of two parts, a polytropic part and a thermal
part, taking the form

$\displaystyle p=p_{\rm{p}}+p_{\rm{th}},$ (4.13) $\displaystyle\rho e=(\rho
e)_{\rm{p}}+(\rho e)_{\rm{th}}.$ (4.14)

The polytropic part is given by

$p_{\rm{p}}=p_{\rm{p}}(\rho)=\begin{cases}\kappa_{1}\rho^{\gamma_{1}},&\rho<\rho_{\rm{nuc}},\\\
\kappa_{2}\rho^{\gamma_{2}},&\rho\geq\rho_{\rm{nuc}},\end{cases}$ (4.15)

where $\rho_{\rm{nuc}}=2\times 10^{14}\,\text{g/cm}^{3}$ is the nuclear
density parameter and separates two different regimes with different adiabatic
indexes, $\gamma_{1}=1.325$ and $\gamma_{2}=2.5$ (This mimics the stiffening
observed in more realistic EoSs as the matter composition transitions from
consisting of nucleons and nuclei to bulk nuclear matter.) The polytropic
internal energy density is given by

$(\rho e)_{\rm{p}}=(\rho
e)_{\rm{p}}(\rho)=\begin{cases}E_{1}\rho^{\gamma_{1}},&\rho<\rho_{\rm{nuc}},\\\
E_{2}\rho^{\gamma_{2}}+E_{3}\rho,&\rho\geq\rho_{\rm{nuc}},\end{cases}$ (4.16)

where the parameters $E_{1},E_{2},E_{3},\kappa_{1},\kappa_{2}$ are given by

$\displaystyle
E_{1}=\frac{\kappa}{\gamma_{1}-1},\quad\kappa_{1}=\kappa,\quad\kappa_{2}=(\gamma_{2}-1)E_{2},$
$\displaystyle
E_{2}=\frac{\kappa}{\gamma_{2}-1}\rho_{\rm{nuc}}^{\gamma_{1}-\gamma_{2}},\quad
E_{3}=\frac{\gamma_{2}-\gamma_{1}}{\gamma_{2}-1}E_{1}\rho_{\rm{nuc}}^{\gamma_{1}-1}.$
(4.17)

One can easily check that the polytropic pressure and internal energy density
are both continuous across the density $\rho=\rho_{\rm{nuc}}$. The thermal
part is given by

$\displaystyle p_{\rm{th}}=(\gamma_{\rm{th}}-1)(\rho e)_{\rm{th}},\qquad(\rho
e)_{\rm{th}}=\rho e-(\rho e)_{\rm{p}},$ (4.18)

where $\gamma_{\rm{th}}=1.5$. We note that the initial thermal pressure is
zero in this test. Combining the above expressions, we can write the complete
EoS in this test as

$\displaystyle p=p(\rho,e)=\begin{cases}(\gamma_{\rm{th}}-1)\rho
e+\frac{\gamma_{1}-\gamma_{\rm{th}}}{\gamma_{1}-1}\kappa\rho^{\gamma_{1}},&\rho<\rho_{\rm{nuc}},\\\
(\gamma_{\rm{th}}-1)\rho
e+\frac{\gamma_{2}-\gamma_{\rm{th}}}{\gamma_{2}-1}\kappa\rho_{\rm{nuc}}^{\gamma_{1}-\gamma_{2}}\rho^{\gamma_{2}}\\\
\qquad\quad-\frac{(\gamma_{\rm{th}}-1)(\gamma_{2}-\gamma_{1})}{(\gamma_{2}-1)(\gamma_{1}-1)}\kappa\rho_{\rm{nuc}}^{\gamma_{1}-1}\rho,&\rho\geq\rho_{\rm{nuc}}.\end{cases}$
(4.19)

We note that there may be a different $\gamma$ in different regions of the
computational domain ($\gamma_{1}$ versus $\gamma_{2}$) and we use the
$\gamma$ of the innermost cell to calculate $n$ and the corresponding
numerical solution $\theta_{n}$ in Section 3.4.1.

We set the computational domain as $\Omega=[0,1.5\times 10^{3}]$ km with a
geometrically increasing cell spacing

$\Delta r_{j}=r_{j+\frac{1}{2}}-r_{j-\frac{1}{2}}=a^{j-1}\Delta r_{1},\qquad
j=1,...,N,$ (4.20)

such that the mesh can be defined by specifying the size of the innermost cell
$\Delta r_{1}$ and the increasing rate $a$. Different values of $\Delta r_{1}$
and $a$ have been utilized in the test with values specified in Table 5. We
use the reflective boundary condition for the inner boundary and zeroth-order
extrapolation for the outer boundary. We set $k=2$ and use the third-order RK
method (3.81)-(3.83) in this test. The simulation is performed from $t=0$ to
$t=0.11$ s. According to the description in Janka et al. (1993); Käppeli &
Mishra (2016), the central density will continue to increase until it exceeds
nuclear density $\rho_{\rm{nuc}}$ and the EoS stiffens to form an inner core
that eventually settles to a new equilibrium configuration (the proto-neutron
star). Due to its inertia, the inner core overshoots its equilibrium and
rebounds to form the shock wave. This is the so-called core bounce, and in
this paper the time of bounce is set as the time when the average density
within the innermost 2 km, which is called central density, reaches its
maximum. Due to the absence of energy losses in our model (i.e., from
deleptonization by neutrinos and dissociation of nuclei below the shock), the
shock wave does not stall, but propagates towards the outer boundary of the
domain.

We note that the dynamics before bounce is similar to the case discussed in
Section 4.4. We refer to the top right panel in Figure 3 for the evolution of
the velocity, and the thermal energy ratio $\mathcal{P}_{th}={(\rho
e)_{th}}/{(\rho e)}$ is almost zero across the whole computational domain
before bounce. To illustrate the dynamics after bounce, we refer to Figure 4,
which shows the fluid velocity and thermal energy ratio versus radius for
select time slices. We can see the shock forms at bounce at a radius between
10 and 20 km, and then gradually propagates to the outer boundary. The thermal
energy remains very small in the inner core, below the location where the
shock formed, while it increases sharply across the shock. Behind the initial
shock, several smaller shocks form and propagate radially as a result of
oscillations in the proto-neutron star as it settles to a hydrostatic
equilibrium state.

Figure 4: Example 4.5, fluid velocity and thermal energy ratio versus radius
after bounce. We use $N=256$ cells and select 6 time slices after the bounce.

We test the proposed well-balanced and energy conserving DG method and the
standard DG method with different number of grids and present them in Table 5,
from which we observe that the time of bounce, the central density of the
bounce, and the final central density at $t=110$ ms are very similar for all
the cases $N=128,256,512,1024,2048$. We show the central density as a function
of time in Figure 5. Both the proposed and standard DG schemes simulate this
test well. In the zoom-in figure, the proposed scheme is shown to be slightly
better than the standard scheme for $N=256$ and $t\in[91,94]$. In Figure 6, we
show the density versus radius at $t=0.11$ s for the case $N=128,256$. We can
observe that there are small shocks at the region $r\in[200,1100]$, and our
proposed scheme performs much better than standard scheme in capturing these
shocks (when compare with the high resolution reference simulation),
especially for the case with $N=256$.

At last, we define the energies as follows

$\displaystyle E_{\rm int}=\int_{\Omega}\rho e\,4\pi
r^{2}\,\mathrm{d}r,~{}~{}E_{\rm kin}=\int_{\Omega}\frac{1}{2}\rho u^{2}\,4\pi
r^{2}\,\mathrm{d}r,$ $\displaystyle E_{\rm
grav}=\int_{\Omega}\frac{1}{2}\rho\,\Phi\,4\pi r^{2}\,\mathrm{d}r,$ (4.21)

where $E_{\rm int}$, $E_{\rm kin}$, and $E_{\rm grav}$ denote the internal
energy, kinetic energy, and gravitational energy, respectively. We list these
three energies $E_{\rm int}$, $E_{\rm kin}$, $-E_{\rm grav}$, and the total
energy conservation $\Delta E$ in (4.12) for different number of cells $N$ at
time $t=0.11$ s in Table 6. Our objective is to study how different schemes
and limiters affect the total energy conservation $\Delta E$. Three different
cases are considered in this table: our well-balanced and total-energy-
conserving scheme, the standard RKDG scheme, and the standard scheme with the
new limiter correction (3.90) (results for this latter scheme are also plotted
in the bottom panels in Figure 5). The reason for including the standard
scheme with the correction is motivated by results from Pochik et al. (2021),
which suggest that limiters may negatively impact the energy conservation
properties of the standard DG scheme for the Euler–Poisson system. From Table
6 (rightmost column), we can see that the well-balanced scheme can maintain
the total energy conservation to round-off errors. For the standard scheme,
neither the case with the standard limiter or the case with the correction
term can maintain the round-off errors. However, we note that the standard
scheme with the correction is substantially better than the standard scheme
with the standard limiter. We plot $E_{\rm int}$, $E_{\rm kin}$, $-E_{\rm
grav}$, and total energy conservation $\Delta E$ versus time in Figure 7 for
the simulations with $N=128$ and $N=256$. We can see that the total energy
conservation for the standard scheme increases rapidly near bounce, and
remains relatively constant thereafter, while for our proposed scheme the
change in the total energy remains small and is not affected by core bounce.

Table 5: Example 4.5, the time of bounce, central density at the bounce time, and central density at the final time for different number of cells. The left and right columns below each label represent the result of the proposed scheme and standard scheme, respectively. $N$ |  $\Delta r_{1}$ [km] | $a-1$ |  $t_{b}$ [ms] |  $\rho_{b}$ [$10^{14}\,\text{g/cm}^{3}$] |  $\rho_{f}$ [$10^{14}\,\text{g/cm}^{3}$]
---|---|---|---|---|---
128 | 2 | $2.292\times 10^{-2}$ | 91.10 | 91.09 | 3.65 | 3.66 | 2.87 | 2.81
256 | 1 | $1.136\times 10^{-2}$ | 91.13 | 91.13 | 3.68 | 3.68 | 2.81 | 2.79
512 | 0.5 | $5.659\times 10^{-3}$ | 91.16 | 91.16 | 3.65 | 3.63 | 2.81 | 2.80
1024 | 0.25 | $2.823\times 10^{-3}$ | 91.16 | 91.16 | 3.63 | 3.63 | 2.81 | 2.80
2048 | 0.125 | $1.410\times 10^{-3}$ | 91.17 | 91.17 | 3.62 | 3.62 | 2.81 | 2.80
Table 6: Example 4.5, four energies at time $t=0.11$ s. We compare the results of three schemes in this table for different number of cells $N$: the well-balanced and total-energy-conserving scheme, the standard scheme, the standard scheme with the new limiter correction (3.90). $N$ | Case | $E_{\rm int}\,[10^{51}\,\text{erg}]$ | $E_{\rm kin}\,[10^{51}\,\text{erg}]$ | $-E_{\rm grav}\,[10^{51}\,\text{erg}]$ | $\Delta E\,[10^{51}\,\text{erg}]$
---|---|---|---|---|---
128 | wb | 120.0 | 3.658 | 122.6 | 4.386$\times 10^{-11}$
standard | 117.7 | 4.091 | 119.1 | 1.269
standard with correction | 119.0 | 3.838 | 121.0 | 4.219$\times 10^{-2}$
256 | wb | 117.7 | 3.452 | 120.0 | 2.886$\times 10^{-10}$
standard | 116.8 | 3.681 | 118.8 | 0.425
standard with correction | 117.3 | 3.543 | 119.6 | 5.976$\times 10^{-3}$
512 | wb | 117.2 | 3.509 | 119.7 | 2.395$\times 10^{-10}$
standard | 116.9 | 3.602 | 119.2 | 0.170
standard with correction | 117.1 | 3.546 | 119.5 | 1.448$\times 10^{-3}$
1024 | wb | 117.2 | 3.542 | 119.7 | 5.404$\times 10^{-10}$
standard | 117.1 | 3.584 | 119.5 | 0.112
standard with correction | 117.1 | 3.559 | 119.6 | 3.545$\times 10^{-4}$
2048 | wb | 117.2 | 3.556 | 119.7 | 1.466$\times 10^{-9}$
standard | 117.1 | 3.578 | 119.6 | 0.038
standard with correction | 117.1 | 3.566 | 119.7 | 4.610$\times 10^{-5}$

Figure 5: Example 4.5, central density as a function of time for the proposed
(top two), the standard (mid two), and the standard with correction (3.90)
(bottom two) DG schemes with $N$=128 (blue dashed), 256 (red dash-dotted), 512
(green dotted) and 2048 (black solid). The right figures represent zoomed-in
versions for $t\in[90,110]$.

Figure 6: Example 4.5, the mass density versus radius at $t=0.11$ s of the
proposed (blue dashed) and the standard DG scheme (red dash-dotted) with
$N=128$ (top two), $256$ (bottom two) compared with a reference solution of
$N=2048$ (black solid). The right two figures represent the zoom-in version at
$r\in[200,1100]$ km.

Figure 7: Example 4.5, the time history of the internal energy $E_{\rm int}$
(blue solid), kinetic energy $E_{\rm kin}$ (green dashed), negative
gravitational energy $-E_{\rm grav}$ (red dash-dotted) and change in total
energy $\Delta E$ (black dotted), with $N=128$ (left figures) and $N=256$
(right figures). We compared the solutions of our proposed scheme (in the top
figures) and the standard DG scheme (in the mid figures) and standard DG
scheme with correction term (3.90) (in the bottom figures).

##  5 Summary and Conclusion

We have developed high-order, total-energy-conserving, and well-balanced
discontinuous Galerkin (DG) methods for solving the Euler–Poisson equations in
spherical symmetry. Our proposed scheme can preserve polytropic steady states
and the total energy up to round-off errors. Key to these properties are the
new way of recovering the steady states, the well-balanced numerical flux, the
novel source term approximations (the well-balanced and total energy
conserving parts), the total energy correction term for the limiter, and the
newly defined time discretization. We have compared the performance of our
proposed scheme with the standard scheme in several different situations,
which all demonstrate the benefits of our proposed scheme. In all these
examples, we can observe the round-off errors for the steady state solutions
and total energy conservation, while the standard scheme can not. In our
opinion, the properties of our proposed scheme may be advantageous for
simulating CCSNe in the context of non-relativistic, self-gravitating
hydrodynamics.

There are still challenges that remain to be solved in future works.
Importantly, CCSNe, and related systems where the methods developed here could
be applicable, are inherently multidimensional due to, e.g., rotation,
hydrodynamic instabilities, and magnetic fields (Müller, 2020). The steady
states considered in this work are valid only in spherical symmetry, and it
will likely become much more complicated to generalize the well-balanced
property to multiple spatial dimensions, which is the main reason we did not
consider multidimensional methods in this paper. For extensions to multiple
spatial dimensions, the main difficulty relates to how the desired steady
states are characterized. However, for problems that can be characterized as
being nearly spherically symmetric (i.e., where the gravitational potential is
dominated by the monopole component), such as CCSNe originating from slowly
rotating stars, the methods developed here may potentially still be
beneficial, but this remains to be investigated. The extension of the energy
conservation property to multiple spatial dimensions appears to be more
straightforward, and will be considered in a future study. Another topic to
consider in a future work is the generalization of the well-balanced property
to tabulated nuclear EoSs needed for more physically realistic models.

## Acknowledgements

This work was carried out when W. Zhang was visiting Department of
Mathematics, The Ohio State University under the support of the China
Scholarship Council (CSC NO. 201906340196). The work of Y. Xing was partially
supported by the NSF grant DMS-1753581. E. Endeve acknowledges support from
the NSF Gravitational Physics Theory Program (NSF PHY 1806692 and 2110177) and
the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the US
Department of Energy Office of Science and the National Nuclear Security
Administration.

## Data Availability Statements

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Armitage (2011) Armitage P. J., 2011, Annual Review of Astronomy and Astrophysics, 49, 195
  * Audusse et al. (2004) Audusse E., Bouchut F., Bristeau M.-O., Klein R., Perthame B. t., 2004, SIAM Journal on Scientific Computing, 25, 2050
  * Berberich et al. (2021) Berberich J. P., Käppeli R., Chandrashekar P., Klingenberg C., 2021, Communications in Computational Physics, 30, 666
  * Bermudez & Vazquez (1994) Bermudez A., Vazquez M. E., 1994, Computers & Fluids, 23, 1049
  * Bruenn et al. (2020) Bruenn S. W., et al., 2020, ApJS, 248, 11
  * Burrows et al. (2020) Burrows A., Radice D., Vartanyan D., Nagakura H., Skinner M. A., Dolence J. C., 2020, MNRAS, 491, 2715
  * Castro & Parés (2020) Castro M. J., Parés C., 2020, Journal of Scientific Computing, 82, 1
  * Chandrashekar & Klingenberg (2015) Chandrashekar P., Klingenberg C., 2015, SIAM Journal on Scientific Computing, 37, B382
  * Chandrashekar & Zenk (2017) Chandrashekar P., Zenk M., 2017, Journal of Scientific Computing, 71, 1062
  * Cockburn & Shu (1989) Cockburn B., Shu C.-W., 1989, Mathematics of computation, 52, 411
  * Couch et al. (2013) Couch S. M., Graziani C., Flocke N., 2013, The Astrophysical Journal, 778, 181
  * Endeve et al. (2019) Endeve E., et al., 2019, in Journal of Physics: Conference Series. p. 012014
  * Franck & Mendoza (2016) Franck E., Mendoza L. S., 2016, Journal of Scientific Computing, 69, 314
  * Gallardo et al. (2007) Gallardo J. M., Parés C., Castro M., 2007, Journal of Computational Physics, 227, 574
  * Ghosh & Constantinescu (2016) Ghosh D., Constantinescu E. M., 2016, AIAA Journal, 54, 1370
  * Grosheintz-Laval & Käppeli (2019) Grosheintz-Laval L., Käppeli R., 2019, Journal of Computational Physics, 378, 324
  * Grosheintz-Laval & Käppeli (2020) Grosheintz-Laval L., Käppeli R., 2020, Journal of Computational Physics, 423, 109805
  * Guo (1998) Guo Y., 1998, Communications in mathematical physics, 195, 249
  * Hanawa (2019) Hanawa T., 2019, in Journal of Physics: Conference Series. p. 012015
  * Janka et al. (1993) Janka H.-T., Zwerger T., Moenchmeyer R., 1993, Astronomy and Astrophysics, 268, 360
  * Jiang & Goodman (2011) Jiang Y.-F., Goodman J., 2011, The Astrophysical Journal, 730, 45
  * Jiang et al. (2013) Jiang Y.-F., Belyaev M., Goodman J., Stone J. M., 2013, New Astronomy, 19, 48
  * Käppeli & Mishra (2014) Käppeli R., Mishra S., 2014, Journal of Computational Physics, 259, 199
  * Käppeli & Mishra (2016) Käppeli R., Mishra S., 2016, Astronomy & Astrophysics, 587, A94
  * Klingenberg et al. (2019) Klingenberg C., Puppo G., Semplice M., 2019, SIAM Journal on Scientific Computing, 41, A695
  * LeVeque (1998) LeVeque R. J., 1998, Journal of computational physics, 146, 346
  * Lentz et al. (2015) Lentz E. J., et al., 2015, ApJL, 807, L31
  * Li & Xing (2016a) Li G., Xing Y., 2016a, Journal of Scientific Computing, 67, 493
  * Li & Xing (2016b) Li G., Xing Y., 2016b, Journal of Computational Physics, 316, 145
  * Li & Xing (2018) Li G., Xing Y., 2018, Journal of Computational Physics, 352, 445
  * Maciel (2015) Maciel W. J., 2015, Introduction to stellar structure. Springer
  * McKee & Ostriker (2007) McKee C. F., Ostriker E. C., 2007, Annual Review of Astronomy and Astrophysics, 45, 565
  * Melson et al. (2015) Melson T., Janka H.-T., Bollig R., Hanke F., Marek A., Müller B., 2015, ApJL, 808, L42
  * Mikami et al. (2008) Mikami H., Sato Y., Matsumoto T., Hanawa T., 2008, The Astrophysical Journal, 683, 357
  * Mullen et al. (2021) Mullen P., Hanawa T., Gammie C., 2021, The Astrophysical Journal Supplement Series, 252, 30
  * Müller (2020) Müller B., 2020, Living Rev. Comput. Astrophys., 6, 3
  * Müller & Steinmetz (1995) Müller E., Steinmetz M., 1995, Computer Physics Communications, 89, 45
  * Müller et al. (2010) Müller B., Janka H.-T., Dimmelmeier H., 2010, The Astrophysical Journal Supplement Series, 189, 104
  * Noelle et al. (2007) Noelle S., Xing Y., Shu C.-W., 2007, Journal of Computational Physics, 226, 29
  * Norsett & Wanner (1987) Norsett E. H. S., Wanner G., 1987, Solving Ordinary Differential Equations I: Nonsti Problems
  * Ostriker et al. (2001) Ostriker E. C., Stone J. M., Gammie C. F., 2001, The Astrophysical Journal, 546, 980
  * Parés & Parés-Pulido (2021) Parés C., Parés-Pulido C., 2021, Journal of Computational Physics, 425, 109880
  * Pochik et al. (2021) Pochik D., Barker B. L., Endeve E., Buffaloe J., Dunham S. J., Roberts N., Mezzacappa A., 2021, The Astrophysical Journal Supplement Series, 253, 21
  * Simon et al. (2016) Simon J. B., Armitage P. J., Li R., Youdin A. N., 2016, The Astrophysical Journal, 822, 55
  * Skinner et al. (2019) Skinner M. A., Dolence J. C., Burrows A., Radice D., Vartanyan D., 2019, ApJS, 241, 7
  * Suzuki (2011) Suzuki M., 2011, Kinetic & Related Models, 4, 569
  * Thomann et al. (2019) Thomann A., Zenk M., Klingenberg C., 2019, International Journal for Numerical Methods in Fluids, 89, 465
  * Toro (2013) Toro E. F., 2013, Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science & Business Media
  * Veiga et al. (2019) Veiga M. H., Velasco-Romero D. A., Abgrall R., Teyssier R., 2019, Communications in Computational Physics, 26, 1
  * Wu & Xing (2021) Wu K., Xing Y., 2021, SIAM Journal on Scientific Computing, 43, A472
  * Xing (2014) Xing Y., 2014, Journal of Computational Physics, 257, 536
  * Xing & Shu (2005) Xing Y., Shu C.-W., 2005, Journal of Computational Physics, 208, 206
  * Xing & Shu (2013) Xing Y., Shu C.-W., 2013, Journal of Scientific Computing, 54, 645
  * Xing et al. (2010) Xing Y., Zhang X., Shu C.-W., 2010, Advances in Water Resources, 33, 1476
  * Xu et al. (2010) Xu K., Luo J., Chen S., 2010, Adv. Appl. Math. Mech, 2, 200
  * Yahil (1983) Yahil A., 1983, The Astrophysical Journal, 265, 1047