# Maximum principle preserving and entropy stable time implicit DGSEM for nonlinear scalar conservation laws

Florent Renac

<figure><img src="image_1.png" /></figure>

## To cite this version:

Florent Renac. Maximum principle preserving and entropy stable time implicit DGSEM for nonlinear scalar conservation laws. ESAIM: Mathematical Modelling and Numerical Analysis, 2025, 59 (5), pp.2583 - 2612. 10.1051/m2an/2025070. hal-05264803

HAL Id: hal-05264803

https://hal.science/hal-05264803v1

Submitted on 17 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MAXIMUM PRINCIPLE PRESERVING AND ENTROPY STABLE TIME IMPLICIT DGSEM FOR NONLINEAR SCALAR CONSERVATION LAWS

FLORENT RENAC*

**Abstract.** This work concerns the analysis of the discontinuous Galerkin spectral element method (DGSEM) with implicit time stepping for the numerical approximation of nonlinear scalar conservation laws in multiple space dimensions. First we consider the DGSEM with a backward Euler time stepping, then a space-time DGSEM discretization to remove the restriction on the time step. We design first-order graph viscosities in space, and in time for the space-time DGSEM, to make the schemes maximum principle preserving and entropy stable for every admissible convex entropy. We also establish well-posedness of the discrete problems by showing existence and uniqueness of the solutions to the nonlinear implicit algebraic relations that need to be solved at each time step. We then use these low-order schemes as building blocks to design a two-step limiter that successfully captures the physical solution, imposes the maximum principle and entropy stability for any convex entropy imposed by the user on the high-order DGSEM scheme, while keeping its accuracy in smooth regions. These properties hold at any approximation order in space and time and without any constraint on the time step. Numerical experiments in one and two space dimensions are presented to illustrate the properties of these schemes.

**Mathematics Subject Classification.** 35L65, 65M70, 65M12.

Received April 15, 2025. Accepted August 4, 2025.

## 1. INTRODUCTION

We are here interested in the accurate and robust approximation of the following Cauchy problem with a hyperbolic scalar nonlinear conservation law in $d \ge 1$ space dimensions:

$$ \partial_t u + \nabla \cdot \mathbf{f}(u) = 0 \quad \text{in } \Omega \times (0, \infty), \qquad (1.1a) $$

$$ u(\cdot, 0) = u_0(\cdot) \quad \text{in } \Omega, \qquad (1.1b) $$

with $\Omega \subset \mathbb{R}^d$ and $u_0$ in $L^\infty(\mathbb{R}^d)$, whose solution may develop discontinuities in finite time and (1.1a) has to be understood in the sense of distributions where we look for weak solutions that are piecewise smooth solutions that satisfy jump relations at a surface of discontinuity. Weak solutions are not necessarily unique and (1.1) must be supplemented with further admissibility conditions to select the physical solution. We here focus on entropy inequalities for any admissible convex entropy $\eta(u)$ and smooth entropy flux $\mathbf{q}(u) = \int^u \eta'(w)\mathbf{f}'(w) dw$:

$$ \partial_t \eta(u) + \nabla \cdot \mathbf{q}(u) \le 0 \quad \text{in } \Omega \times (0, \infty), \qquad (1.2) $$

*Keywords and phrases.* Hyperbolic scalar equations, maximum principle, entropy stability, discontinuous Galerkin method, summation-by-parts, backward Euler, limiter.

DAAA, ONERA, Université Paris Saclay, F-92322 Châtillon, France.

*Corresponding author: florent.renac@onera.fr
in the sense of distributions. This leads to $L^2$ stability for a uniformly convex entropy [11] and, for $d=1$, an inequality for one entropy is enough to select the physical weak solution when the flux $f(u)$ is strictly convex [54]. For a general flux, it is known that there exists a unique weak solution to (1.1) in $L^\infty(\mathbb{R}^+ \times \mathbb{R}^d)$ that satisfies inequalities (1.2) for any admissible entropy pairs [42, 47, 53, 70]. Existence is obtained by the vanishing viscosity method, while uniqueness follows from Kružkov's seminal work [42] (see, e.g., [23], Sect. 3.3) using the so-called family of Kružkov entropy pairs (see Rem. 4.4).

We assume a locally Lipschitz continuous flux $f(\cdot)$, and entropy solutions to (1.1) are known to satisfy a maximum principle:

$$ \underset{\mathbf{x} \in \Omega}{\operatorname{ess\,inf}} u_0(\mathbf{x}) := m \le u(\mathbf{x}, t) \le M := \underset{\mathbf{x} \in \Omega}{\operatorname{ess\,sup}} u_0(\mathbf{x}) \quad \text{a.e. in } \Omega \times (0, \infty). \quad (1.3) $$

Let us introduce the Lipschitz constant $L_f = L_f(m, M)$ of $f(\cdot)$:

$$ L_f := \sup_{m \le u^+ \ne u^- \le M} \sup_{|\mathbf{n}|=1} \frac{|(\mathbf{f}(u^+) - \mathbf{f}(u^-)) \cdot \mathbf{n}|}{|u^+ - u^-|}. \quad (1.4) $$

The design of the numerical schemes we consider in this work conveniently relies on the Riemann problem in the unit direction $\mathbf{n}$ in $S^{d-1}$: given data $u^-$ and $u^+$, solve

$$ \partial_t u + \partial_x f(u) \cdot \mathbf{n} = 0 \quad \text{in } \mathbb{R} \times (0, \infty), \quad (1.5a) $$

$$ u(x, 0) = u_0(x) := \begin{cases} u^{-} & \text{if } x := \mathbf{x} \cdot \mathbf{n} < 0, \\ u^{+} & \text{if } x > 0. \end{cases} \quad (1.5b) $$

The maximum wave speed in the solution to (1.5) is bounded by $L_f(m, M)$ in (1.4) with $m = \min(u^-, u^+)$ and $M = \max(u^-, u^+)$. Integrating (1.5) and the associated entropy inequality over space and time, we have for all $\alpha \ge L_f$ [25, 29]

$$ m \le W(u^{-}, u^{+}, \mathbf{n}) := \frac{u^{-} + u^{+}}{2} - \frac{(\mathbf{f}(u^{+}) - \mathbf{f}(u^{-})) \cdot \mathbf{n}}{2\alpha} \le M, \quad (1.6a) $$

$$ \eta(W(u^{-}, u^{+}, \mathbf{n})) \le \frac{\eta(u^{-}) + \eta(u^{+})}{2} - \frac{(\mathbf{q}(u^{+}) - \mathbf{q}(u^{-})) \cdot \mathbf{n}}{2\alpha}. \quad (1.6b) $$

We are here interested in the approximation of (1.1) with a high-order discretization that satisfies the above properties at the discrete level. We consider the discontinuous Galerkin spectral element method (DGSEM) based on the collocation between interpolation and quadrature points [40] and tensor products of one-dimensional function bases and quadrature rules. Using diagonal norm summation-by-parts (SBP) operators and the entropy conservative (EC) numerical fluxes from Tadmor [65], semi-discrete EC schemes have been derived in [6, 17] for nonlinear conservation laws and one given entropy. Entropy stable (ES) DGSEM on hexahedral [22, 71] and simplex [7, 10] meshes have been proposed for this given entropy by using the same framework. The particular form of the SBP operators allows to take into account the numerical quadratures that approximate integrals in the numerical scheme compared to other techniques that require their exact evaluation to satisfy the entropy inequality [31, 33]. The DGSEM thus provides a general framework for the design of ES schemes for nonlinear systems of conservation laws. Numerical experiments in [7, 9, 22, 60, 71] highlight the benefits on stability and robustness of the computations, though this not guarantees to preserve neither the entropy stability at the fully discrete level, nor positivity of the numerical solution. Designs of fully discrete ES and positivity-preserving DGSEM have been proposed in [1, 5, 12, 13, 35, 55, 58, 59, 61, 68, 69] among others. Likewise, these schemes are ES for one given entropy only, which may be not sufficient to capture the entropy weak solution as will be observed in the numerical experiments of Section 6. Discrete inequalities for some families of entropy pairs are indeed usually required to prove convergence of the discrete solution to the unique entropy weak solution [4, 8, 46, 64, 66].
These schemes are usually analyzed in semi-discrete form for which the time derivative is not discretized, or when coupled with explicit in time discretizations. Time explicit integration may however become prohibitive for long time simulations or when looking for stationary solutions due to the strong CFL restriction on the time step which gets smaller as the approximation order of the scheme increases to ensure either linear stability [2, 21, 41], or positivity of the approximate solution [74, 75]. We here consider and analyze a DGSEM discretization in space associated with a time implicit integration to circumvent the CFL condition. We first consider a first-order backward Euler time integration, which constitutes an efficient time stepping when looking for steady-state solutions. We then consider a DGSEM discretization in time more adapted for time-resolved simulations.

Little is known about the properties of time implicit DGSEM schemes, apart from the entropy stability for one given entropy which holds providing the semi-discrete scheme is ES due to the dissipative character of many time implicit integration techniques. The works in [52, 57] analyze time implicit modal discontinuous Galerkin (DG) and DGSEM for the discretization of a 1D linear scalar hyperbolic equation and show that a lower bound on the time step is required for the cell-averaged solution to satisfy a maximum principle at the discrete level. A linear scaling limiter of the DG solution around its cell-average [74] is then used to obtain a maximum principle preserving (MPP) scheme. Unfortunately, these schemes are not MPP in multiple space dimensions even on Cartesian grids neither for the full polynomial solution, nor for the cell-averaged solution [51, 52]. Related works on time implicit DG schemes concern radiative transfer equations with enriched function space [51], reduced order quadrature rules [72], and limiters [72]; or the formulation of a constrained optimization problem [67].

We aim at designing schemes that satisfy the maximum principle and all entropy inequalities. The MPP property alone may not be sufficient to capture the entropy weak solution, nor may the satisfaction of only one entropy inequality. For instance, the local extrema diminishing method [32] is MPP, but may not converge to the entropy solution ([26], Lem. 3.2). Likewise, the present experiments show that the DGSEM, though ES for one given entropy, may fail to capture the physical solution. We here add a first-order graph viscosity to the DGSEM schemes. Graph viscosity is a common method of artificial dissipation based on the connectivity graph of the degrees of freedom (DOFs) [25, 27, 32, 55]. The graph viscosity is an efficient technique to impose invariant domain properties to the scheme such as positivity, maximum principle, or minimum entropy principle. Some limiters in high-order schemes also use low-order schemes with graph viscosity as building-blocks such as the flux-corrected transport (FCT) limiter [3, 15, 44, 45, 52, 73], or convex limiting [27, 55]. This is the path we follow and propose a two-step limiter that uses the low-order solution obtained with graph viscosity to successfully impose on the high-order scheme the maximum principle and all entropy inequalities that are required by the user, while keeping its accuracy for smooth solutions. The limited scheme also captures the physical solution, when the unlimited high-order scheme fails. Using the DGSEM with graph viscosity as the low-order scheme allows a minimal modification of any existing DGSEM code that is easy to implement and efficient. To the best of our knowledge, the present limited scheme might be the first time implicit DGSEM scheme that is formally high-order, MPP and ES. These properties hold at any approximation order in space and time and without any constraint on the time step.

The graph viscosity is here based on the connectivity graph of the DOFs either within the mesh element when the DGSEM is associated to the backward-Euler time integration, or within the space-time discretization element for the space-time DGSEM. It is local to those elements and do not use connectivity with neighboring elements making its implementation straightforward. We derive simple estimates on the viscosity coefficients to make the schemes MPP and ES for all admissible convex entropies. We also analyze the well-posedness of the discrete problems and prove existence and uniqueness of the solutions to the nonlinear implicit discrete problems that need to be solved at each time step. One key ingredient to prove these properties is the derivation in Lemmas 2.1 and 2.2 of discrete counterparts to (1.6) for the EC two-point fluxes in space. Although those properties are well known for positivity-preserving and entropic three-point schemes and correspond to consistency relations with the integral forms of (1.1a) and (1.2) ([29], Thm. 3.1), they constitute new results for EC two-point fluxes up to our knowledge. Given the importance of discrete counterparts to (1.6) in the analysis of general discretization schemes [29], these new results may also be useful to analyze time implicit and time explicit flux differencing schemes that use EC fluxes [17].
The paper is organized as follows. Section 2 describes the DGSEM discretization in space, recall some properties of the metric terms, and of the numerical fluxes. The properties of the DGSEM with a backward Euler time stepping and graph viscosity in space are analyzed in Section 3, while the space-time DGSEM with space and time graph viscosity is described and analyzed in Section 4. The limiting procedure is introduced in Section 5. The results are assessed by numerical experiments in Section 6 and concluding remarks about this work are given in Section 7.

## 2. DGSEM APPROXIMATION IN SPACE

We here describe the semi-discrete weak formulation of problem (1.1) with the DGSEM. The domain is discretized with a shape-regular mesh $\Omega_h \subset \mathbb{R}^d$ consisting of nonoverlapping and nonempty open elements $\kappa$ and we assume that it forms a partition of $\Omega$. By $h_\Omega := \sup_{\kappa \in \Omega_h} \text{diam} \kappa$ we denote the mesh size. For the sake of clarity, we introduce the DGSEM in two space dimensions $d = 2$, its extension (resp., restriction) to $d = 3$ (resp., $d = 1$) being straightforward.

### 2.1. The DGSEM approximation in space

We look for approximate solutions in the function space of discontinuous polynomials $V_h^p = \{\phi \in L^2(\Omega_h) : \phi|_{\kappa} \circ \mathbf{x}_{\kappa} \in Q^p(I^2) \forall \kappa \in \Omega_h\}$, where $Q^p(I^2)$ denotes the space of functions over the master element $I^2 := \{\xi = (\xi, \eta) : -1 \le \xi, \eta \le 1\}$ formed by tensor products of polynomials of degree $p \ge 1$ in each direction. Each physical element $\kappa$ is the image of $I^2$ through the mapping $\mathbf{x} = \mathbf{x}_{\kappa}(\xi)$. Likewise, each edge $e$ is the image of $I = [-1, 1]$ through the mapping $\mathbf{x} = \mathbf{x}_e(\xi)$. The approximate solution to (1.1) is sought under the form

$$u_h(\mathbf{x}, t) = \sum_{i,j=0}^{p} \phi_{\kappa}^{ij}(\mathbf{x}) U_{\kappa}^{ij}(t) \quad \forall \mathbf{x} \in \kappa, \kappa \in \Omega_h, \forall t \ge 0, \quad (2.1)$$

where $(U_{\kappa}^{ij})_{0 \le i, j \le p}$ are the DOFs in the element $\kappa$. The subset $(\phi_{\kappa}^{ij})_{0 \le i, j \le p}$ constitutes a basis of $V_h^p$ restricted onto the element $\kappa$ and $(p+1)^2$ is its dimension. Let $(\ell_k)_{0 \le k \le p}$ be the Lagrange interpolation polynomials in one space dimension associated to the Gauss-Lobatto quadrature nodes over $I$, $\xi_0 = -1 < \xi_1 < \dots < \xi_p = 1$: $\ell_k(\xi_l) = \delta_{kl}$ for all $0 \le k, l \le p$, with $\delta_{kl}$ the Kronecker delta. We use tensor products of these polynomials and of Gauss-Lobatto nodes (see Fig. 1): $\phi_{\kappa}^{ij}(\mathbf{x}) = \phi_{\kappa}^{ij}(\mathbf{x}_{\kappa}(\xi)) = \ell_i(\xi)\ell_j(\eta)$, which satisfy the following relation at quadrature points $\xi_{i'j'} = (\xi_{i'}, \xi_{j'})$ in $I^2$:

$$\phi_{\kappa}^{ij}(\mathbf{x}_{\kappa}^{i'j'}) = \delta_{ii'}\delta_{jj'}, \quad 0 \le i, j, i', j' \le p, \quad \mathbf{x}_{\kappa}^{i'j'} := \mathbf{x}_{\kappa}(\xi_{i'j'}), \quad (2.2)$$

so the DOFs correspond to point values of the solution: $U_{\kappa}^{ij}(t) = u_h(\mathbf{x}_{\kappa}^{ij}, t)$.

The integrals over elements and faces are approximated by using the Gauss-Lobatto quadrature rules so the quadrature and interpolation nodes are collocated:

$$\int_{-1}^{1} f(\xi) \, d\xi \approx \sum_{k=0}^{p} \omega_k f(\xi_k),$$

with $\omega_k > 0$, satisfying $\sum_{k=0}^{p} \omega_k = \int_{-1}^{1} ds = 2$, the weights and $\xi_k$ the nodes of the quadrature rule over $I$.

The semi-discrete form of the DGSEM in space is now standard and we refer to [71] and references therein for details on its implementation on unstructured grids with high-order mesh elements. It reads

$$\omega_i \omega_j J_{\kappa}^{ij} \frac{dU_{\kappa}^{ij}}{dt} + R_{\kappa}^{ij}(u_h) = 0 \quad \forall \kappa \in \Omega_h, \ 0 \le i, j \le p, t > 0, \quad (2.3)$$
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Inner and outer elements, $\kappa = \kappa^-$ and $\kappa_e^+$, for $d = 2$; definitions of traces $u_h^\pm$ on the interface $e$ and of the unit outward normal vector $\mathbf{n}_e^k = \mathbf{n}_e(\mathbf{x}_e^k)$ with $\mathbf{x}_e^k = \mathbf{x}_e(\xi_k)$; positions of quadrature points in $\kappa^-$, $\mathbf{x}_\kappa^{ij} = \mathbf{x}_\kappa(\xi_{ij})$, and on $e$, $\mathbf{x}_e^k$, for $p = 3$ (bullets ●). The elements $\kappa$ are interpolated on the same grid of quadrature points as the numerical solution, i.e., $\mathbf{x}_\kappa(\xi) = \sum_{i,j=0}^p \ell_i(\xi)\ell_j(\eta)\mathbf{x}_\kappa^{ij}$.</figcaption></figure>

with $J_\kappa^{ij} := \det \mathbf{x}'_\kappa(\xi_i, \eta_j) > 0$ and

$$R_\kappa^{ij}(u_h) = 2 \left( \omega_j \sum_{k=0}^{p} Q_{ik}^p h_{ec} (U_\kappa^{ij}, U_\kappa^{kj}, \mathbf{n}_\kappa^{(ik)j}) + \omega_i \sum_{k=0}^{p} Q_{jk}^p h_{ec} (U_\kappa^{ij}, U_\kappa^{ik}, \mathbf{n}_\kappa^{i(jk)}) \right) \\ + \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \phi_\kappa^{ij}(\mathbf{x}_e^k) \omega_k J_e^k (h(U_\kappa^{ij}, u_h^+(\mathbf{x}_e^k, t), \mathbf{n}_e^k) - \mathbf{f}(U_\kappa^{ij}) \cdot \mathbf{n}_e^k), \quad (2.4)$$

where $\mathbf{n}_e^k = \mathbf{n}_e(\mathbf{x}_e^k)$, $J_e^k = \det \mathbf{x}'_e(\xi_k) > 0$, and by (2.2), $\phi_\kappa^{ij}(\mathbf{x}_e^k) = 1$ if $\mathbf{x}_\kappa^{ij} = \mathbf{x}_e^k$ and $\phi_\kappa^{ij}(\mathbf{x}_e^k) = 0$ otherwise, while the two first arguments in $h(\cdot, \cdot, \cdot)$ correspond to the traces, $u_h^\mp(\mathbf{x}_e^k, t)$, of the numerical solution at $\mathbf{x}_e^k$ on $e$ (see Fig. 1). The quantities $\mathbf{n}_\kappa^{(ik)j}$ and $\mathbf{n}_\kappa^{i(jk)}$ are defined in (2.9) and have been introduced to keep conservation of the scheme [71].

## 2.2. Metric terms

The coefficients $Q_{kl}^p$ in (2.4) are defined from the entries of the discrete derivative matrix, $D_{kl}^p$, [39] as follows:

$$D_{kl}^p = \ell'_l(\xi_k), \quad Q_{kl}^p = \omega_k D_{kl}^p, \quad 0 \le k, l \le p, \qquad (2.5)$$

where the exponent $p$ refers to the polynomial degree. As noticed in [20], the DGSEM satisfies the summation-by-parts (SBP) property [63]:

$$Q_{kl}^p + Q_{lk}^p = \delta_{kp}\delta_{lp} - \delta_{k0}\delta_{l0}, \quad 0 \le k, l \le p, \qquad (2.6)$$

which is the discrete counterpart to integration by parts. We will also use the following relations of partition of unity and integration of the Lagrange polynomial derivatives:

$$\sum_{l=0}^{p} Q_{kl}^{p} = \omega_{k} \sum_{l=0}^{p} D_{kl}^{p} = 0, \quad \sum_{l=0}^{p} Q_{lk}^{p} = \delta_{kp} - \delta_{k0}, \quad 0 \le k \le p; \qquad (2.7)$$

and the metric identities [38], $\sum_{i=1}^d \partial_{\xi_i} (J_\kappa \nabla \xi_i) = 0$, which may be written for $d = 2$ at the discrete level as

$$\sum_{k=0}^{p} D_{ik}^{p} J_{\kappa}^{kj} \nabla \xi(\mathbf{x}_{\kappa}^{kj}) + D_{jk}^{p} J_{\kappa}^{ik} \nabla \eta(\mathbf{x}_{\kappa}^{ik}) = 0 \quad \forall \kappa \in \Omega_h, \ 0 \le i, j \le p.$$
Combining the above relation with (2.7) we get

$$ \sum_{k=0}^{p} D_{ik}^{p} \mathbf{n}_{\kappa}^{(ik)j} + D_{jk}^{p} \mathbf{n}_{\kappa}^{(jk)} = 0 \quad \forall \kappa \in \Omega_h, 0 \le i, j \le p, \qquad (2.8) $$

where

$$ \mathbf{n}_{\kappa}^{(ik)j} = \frac{1}{2} (J_{\kappa}^{ij} \nabla \xi(\xi_{ij}) + J_{\kappa}^{kj} \nabla \xi(\xi_{kj})), \quad \mathbf{n}_{\kappa}^{i(jk)} = \frac{1}{2} (J_{\kappa}^{ij} \nabla \eta(\xi_{ij}) + J_{\kappa}^{ik} \nabla \eta(\xi_{ik})). \qquad (2.9) $$

Further assuming the mesh is watertight ([30], App. B.2), i.e., the geometry is continuous across faces, the volume and face metric terms are related by the following relations when evaluated at faces (with some slight abuse in the notations, but without ambiguity)

$$ \mathbf{n}_{\kappa}^{(pp)j} = J_{e}(\mathbf{x}_{\kappa}^{pj})\mathbf{n}_{e}(\mathbf{x}_{\kappa}^{pj}), \quad \mathbf{n}_{\kappa}^{(00)j} = -J_{e}(\mathbf{x}_{\kappa}^{0j})\mathbf{n}_{e}(\mathbf{x}_{\kappa}^{0j}) \quad \forall 0 \le j \le p, \qquad (2.10a) $$

$$ \mathbf{n}_{\kappa}^{i(pp)} = J_{e}(\mathbf{x}_{\kappa}^{ip})\mathbf{n}_{e}(\mathbf{x}_{\kappa}^{ip}), \quad \mathbf{n}_{\kappa}^{i(00)} = -J_{e}(\mathbf{x}_{\kappa}^{i0})\mathbf{n}_{e}(\mathbf{x}_{\kappa}^{i0}) \quad \forall 0 \le i \le p, \qquad (2.10b) $$

where $\mathbf{n}_e$ denotes the unit normal to $e$ in $\partial\kappa$ pointing outward from $\kappa$ (see Fig. 1).

We assume that the mesh is regular in the following sense: $h_{\Omega} := \sup_{\kappa \in \Omega_h} \text{diam } \kappa < \infty$ and there exists a finite $\beta = \beta(p, h_{\Omega}) > 0$ such that

$$ J_{\kappa}^{ij} \ge \frac{\beta}{2d} h_{\Omega}^{d}, \quad \frac{1}{4d\beta} h_{\Omega}^{d-1} \ge J_{e}^{k} \left| \mathbf{n}_{\kappa}^{(ik)j} \right|, \left| \mathbf{n}_{\kappa}^{(jk)j} \right| > 0 \quad \forall \kappa \in \Omega_h, 0 \le i, j, k \le p, \qquad (2.11) $$

which enforces classical regularity assumptions ([16], Sect. 24):

$$ |\kappa| := \sum_{0 \le i, j \le p} \omega_i \omega_j J_{\kappa}^{ij} \ge \beta h_{\Omega}^{d}, \quad |\partial \kappa| := \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \omega_k J_{e}^{k} \le \frac{1}{\beta} h_{\Omega}^{d-1}. $$

Finally, we define the cell-average operator:

$$ \langle u_h \rangle_{\kappa}(t) := \sum_{0 \le i, j \le p} \omega_i \omega_j \frac{J_{\kappa}^{ij}}{|\kappa|} U_{\kappa}^{ij}(t). \qquad (2.12) $$

## 2.3. Two-point numerical fluxes

### 2.3.1. Two-point fluxes in cells

The numerical flux in the cells is smooth, consistent, $h_{ec}(u, u, \mathbf{n}) = \mathbf{f}(u) \cdot \mathbf{n}$, symmetric, $h_{ec}(u^{-}, u^{+}, \mathbf{n}) = h_{ec}(u^{+}, u^{-}, \mathbf{n})$, linear in $\mathbf{n}$, and is EC [65] for a given pair $(\vartheta(u), \mathbf{g}(u))$ with $\vartheta$ strictly convex:

$$ (v(u^{+}) - v(u^{-}))h_{ec}(u^{-}, u^{+}, \mathbf{n}) = (\psi(u^{+}) - \psi(u^{-})) \cdot \mathbf{n} \quad \forall u^{\pm}, \mathbf{n} \in \mathbb{R}^{d}, \qquad (2.13) $$

with $v(u) = \vartheta'(u)$ the entropy variable and $\psi(u) = v(u)\mathbf{f}(u) - \mathbf{g}(u)$ the entropy flux potential. Lemma 2.1 states discrete counterparts to stability properties (1.6) for the exact Riemann solution (1.5). Together with Lemma 2.2, they are key results to establish the properties of the numerical schemes in Sections 3 and 4.

**Lemma 2.1.** The EC flux (2.13) defined for the entropy pair $(\vartheta(u), \mathbf{g}(u))$ with $\vartheta''(\cdot) > 0$ satisfies the following relations for all $\alpha \ge L_f$ and all convex entropy pairs $(\eta, \mathbf{q})$:

$$ m \le \mathcal{U}(u^{-}, u^{+}, \mathbf{n}) := (1 - \beta)u^{-} + \beta u^{+} - \frac{h_{ec}(u^{-}, u^{+}, \mathbf{n}) - \mathbf{f}(u^{-}) \cdot \mathbf{n}}{2\alpha|\mathbf{n}|} \le M, \qquad (2.14a) $$

$$ \eta(\mathcal{U}(u^{-}, u^{+}, \mathbf{n})) \le \frac{\eta(u^{-}) + \bar{\eta}(u^{-}, u^{+})}{2} - \frac{q_{ec}(u^{-}, u^{+}, \mathbf{n}) - \mathbf{q}(u^{-}) \cdot \mathbf{n}}{2\alpha|\mathbf{n}|}, \qquad (2.14b) $$
with

$$ \beta = \beta(u^{-}, u^{+}) = \frac{1}{2} \int_{0}^{1} \frac{u(\theta v^{+} + (1 - \theta)v^{-}) - u^{-}}{u^{+} - u^{-}} \mathrm{d}\theta \in \left(0, \frac{1}{2}\right], \qquad (2.15\text{a}) $$

$$ \bar{\eta}(u^{-}, u^{+}) = \int_{0}^{1} \eta \circ u(\theta v^{+} + (1 - \theta)v^{-}) \mathrm{d}\theta, \qquad (2.15\text{b}) $$

$$ q_{ec}(u^{-}, u^{+}, \mathbf{n}) = \int_{0}^{1} \mathbf{q} \circ u(\theta v^{+} + (1 - \theta)v^{-}) \cdot \mathbf{n} \mathrm{d}\theta, \qquad (2.15\text{c}) $$

where $v^{\pm} = v(u^{\pm}) = \vartheta'(u^{\pm})$.

**Proof.** Note that Lemma 2.1 holds when $u^{-} = u^{+} = u$ since $\mathcal{U}(u, u, \mathbf{n}) = u$, $\beta(u, u) = \frac{1}{2} \int_{0}^{1} \theta \mathrm{d}\theta = \frac{1}{4}$, and $h_{ec}(u, u, \mathbf{n}) = \mathbf{f}(u) \cdot \mathbf{n}$. We now assume $u^{-} \neq u^{+}$. Using the entropy variable $v(u) = \vartheta'(u)$, we have $\psi \circ u(v) = v\mathbf{f}(u(v)) - \int^{v} w \mathbf{f}'(u(w)) \mathrm{d}w = \int^{v} \mathbf{f} \circ u(w) \mathrm{d}w$ so the EC flux (2.13) has the closed form

$$ h_{ec}(u^{-}, u^{+}, \mathbf{n}) = \int_{0}^{1} \mathbf{f} \circ u(\theta v^{+} + (1 - \theta)v^{-}) \cdot \mathbf{n} \mathrm{d}\theta. $$

By strict convexity of $\vartheta(u)$, $u = u(v)$ is strictly increasing so $u(\theta v^{+} + (1 - \theta)v^{-})$ lies in $[\min(u^{\pm}), \max(u^{\pm})] \subset [m, M]$ and (2.15a) follows directly. Using (2.15a), we rewrite (2.14a) as

$$ \mathcal{U}(u^{-}, u^{+}, \mathbf{n}) = \int_{0}^{1} \frac{u(\theta v^{+} + (1 - \theta)v^{-}) + u^{-}}{2} - \frac{\mathbf{f} \circ u(\theta v^{+} + (1 - \theta)v^{-}) - \mathbf{f}(u^{-})}{2\alpha|\mathbf{n}|} \cdot \mathbf{n} \mathrm{d}\theta $$

which lies in $[m, M]$ from (1.6a). Likewise, using the Jensen's inequality, then (1.6b), we obtain

$$ \begin{aligned} \eta(\mathcal{U}(u^{-}, u^{+}, \mathbf{n})) &\le \int_{0}^{1} \eta \left( \frac{u(\theta v^{+} + (1 - \theta)v^{-}) + u^{-}}{2} - \frac{\mathbf{f} \circ u(\theta v^{+} + (1 - \theta)v^{-}) - \mathbf{f}(u^{-})}{2\alpha|\mathbf{n}|} \cdot \mathbf{n} \right) \mathrm{d}\theta \\ &\le \int_{0}^{1} \frac{\eta \circ u(\theta v^{+} + (1 - \theta)v^{-}) + \eta(u^{-})}{2} - \frac{\mathbf{q} \circ u(\theta v^{+} + (1 - \theta)v^{-}) - \mathbf{q}(u^{-})}{2\alpha|\mathbf{n}|} \cdot \mathbf{n} \mathrm{d}\theta \end{aligned} $$

which completes the proof. □

**Lemma 2.2.** Under the assumptions of Lemma 2.1, $\beta(u^{-}, u^{+})$ in (2.15a) and $\bar{\eta}(u^{-}, u^{+})$ in (2.15b) satisfy

$$ \bar{\eta}(u^{-}, u^{+}) - \eta(u^{-}) \le 2\beta(u^{-}, u^{+})(\eta(u^{+}) - \eta(u^{-})). \qquad (2.16) $$

**Proof.** Set again $v^{\pm} = v(u^{\pm}) = \vartheta'(u^{\pm})$, (2.16) trivially holds for $u^{+} = u^{-}$ and we now assume $u^{+} \neq u^{-}$, thus $v^{+} \neq v^{-}$. Use the change of variables $u(\theta v^{+} + (1 - \theta)v^{-}) = \zeta u^{+} + (1 - \zeta)u^{-}$ in (2.15b), which is one-to-one since $u(v)$ is strictly increasing by strict convexity of $\vartheta(u)$ in Lemma 2.1, to obtain

$$ \begin{aligned} \bar{\eta}(u^{-}, u^{+}) &= \frac{u^{+} - u^{-}}{v^{+} - v^{-}} \int_{0}^{1} \frac{\eta(\zeta u^{+} + (1 - \zeta)u^{-})}{u'(v(\zeta u^{+} + (1 - \zeta)u^{-}))} \mathrm{d}\zeta \\ &\le \frac{u^{+} - u^{-}}{v^{+} - v^{-}} \int_{0}^{1} \frac{\zeta\eta(u^{+}) + (1 - \zeta)\eta(u^{-})}{u'(v(\zeta u^{+} + (1 - \zeta)u^{-}))} \mathrm{d}\zeta \\ &= 2\beta(u^{-}, u^{+})\eta(u^{+}) + (1 - 2\beta(u^{-}, u^{+}))\eta(u^{-}), \end{aligned} $$

where we have used convexity of $\eta(\cdot)$, then $2\beta(u^{-}, u^{+}) = \frac{u^{+} - u^{-}}{v^{+} - v^{-}} \int_{0}^{1} \frac{\zeta \mathrm{d}\zeta}{u'(v(\zeta u^{+} + (1 - \zeta)u^{-}))}$ with the same change of variables. This completes the proof. □

**Remark 2.3 (Square entropy).** For the square entropy, $\vartheta(u) = \frac{u^2}{2}$, the entropy flux reads $g(u) = \int^{u} w \mathbf{f}'(w) \mathrm{d}w$, and we have $v(u) = u$, $h_{ec}(u^{-}, u^{+}, \mathbf{n}) = \int_{0}^{1} \mathbf{f}(\theta u^{+} + (1 - \theta)u^{-}) \cdot \mathbf{n} \mathrm{d}\theta$, $\beta(u^{-}, u^{+}) = \frac{1}{4}$ in Lemma 2.1, while $\bar{\eta}(u^{-}, u^{+}) = \int_{0}^{1} \eta(\theta u^{+} + (1 - \theta)u^{-}) \mathrm{d}\theta \le \frac{\eta(u^{+}) + \eta(u^{-})}{2}$ in agreement with Lemma 2.2.
### 2.3.2. Two-point fluxes at interfaces.

The numerical flux at interfaces is Lipschitz continuous, consistent, $h(u, u, \mathbf{n}) = f(u) \cdot \mathbf{n}$, conservative, $h(u^{-}, u^{+}, \mathbf{n}) = -h(u^{+}, u^{-}, -\mathbf{n})$, and monotone: $h(u^{-}, u^{+}, \mathbf{n})$ is nondecreasing with $u^{-}$ and nonincreasing with $u^{+}$. The Lipschitz constants of $h$ with respect to its two first arguments satisfy (1.4): $L_{1} + L_{2} \geq L_{f}$. The associated three-point scheme is therefore MPP and ES ([23], Chap. 4): for all $m \leq u^{-}, u, u^{+} \leq M$, $\mathbf{n} \in \mathbb{S}^{d-1}$, and $\alpha \geq L_{f}$, we have

$$m \leq u - \frac{h(u, u^{+}, \mathbf{n}) - h(u^{-}, u, \mathbf{n})}{2\alpha} \leq M, \qquad (2.17)$$

and

$$\eta\left(u - \frac{h(u, u^{+}, \mathbf{n}) - h(u^{-}, u, \mathbf{n})}{2\alpha}\right) \leq \eta(u) - \frac{Q(u, u^{+}, \mathbf{n}) - Q(u^{-}, u, \mathbf{n})}{2\alpha}, \qquad (2.18)$$

with the consistent, $Q(u, u, \mathbf{n}) = \mathbf{q}(u) \cdot \mathbf{n}$, and conservative, $Q(u^{-}, u^{+}, \mathbf{n}) = -Q(u^{+}, u^{-}, -\mathbf{n})$, numerical entropy flux. The Godunov [24], Engquist-Osher [14] and the Rusanov [19] schemes satisfy the above properties.

## 3. TIME IMPLICIT DISCRETIZATION AND GRAPH VISCOSITY

We now focus on a time implicit discretization with a backward Euler method in (2.3). The fully discrete scheme in cell $\kappa \in \Omega_h$ at each time step $n \geq 0$ reads

$$\omega_i \omega_j J_{\kappa}^{ij} \frac{U_{\kappa}^{ij, n+1} - U_{\kappa}^{ij, n}}{\Delta t^{(n)}} + R_{\kappa}^{ij}(u_h^{(n+1)}) + V_{\kappa}^{ij}(u_h^{(n+1)}) = 0, \quad 0 \le i, j \le p, \qquad (3.1)$$

where the residuals $R_{\kappa}^{ij}(\cdot)$ are defined by (2.4), $\Delta t^{(n)} = t^{(n+1)} - t^{(n)} > 0$, with $t^{(0)} = 0$, denotes the time step, and $U_{\kappa}^{ij,n} = U_{\kappa}^{ij}(t^{(n)})$. The projection of the initial condition (1.1b) onto the function space reads $U_{\kappa}^{ij,0} = u_0(\mathbf{x}_{\kappa}^{ij})$.

The artificial viscosity term is a graph viscosity [25] of the form

$$V_{\kappa}^{ij}(u_h) := d_{\kappa} \omega_i \omega_j \sum_{k=0}^{p} \frac{\omega_k}{2} |\mathbf{n}_{\kappa}^{(ik)j}| (U_{\kappa}^{ij} - U_{\kappa}^{kj}) + \frac{\omega_k}{2} |\mathbf{n}_{\kappa}^{i(jk)}| (U_{\kappa}^{ij} - U_{\kappa}^{ik}). \qquad (3.2)$$

One recovers the unlimited time implicit high-order in space scheme when $d_{\kappa} = 0$ for all $\kappa \in \Omega_h$. The scheme with $d_{\kappa} > 0$ is formally first-order accurate in space. We will look for conditions on $d_{\kappa}$ for the scheme (3.1) to satisfy the desired properties.

The method is conservative in the sense that by summing (3.1) over $0 \le i, j \le p$ gives for the cell-averaged solution (2.12)

$$\langle u_h^{(n+1)} \rangle_{\kappa} = \langle u_h^{(n)} \rangle_{\kappa} - \frac{\Delta t^{(n)}}{|\kappa|} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \omega_k J_e^k h \left( u_h^{-} (\mathbf{x}_e^k, t^{(n+1)}), u_h^{+} (\mathbf{x}_e^k, t^{(n+1)}), \mathbf{n}_e^k \right), \qquad (3.3)$$

providing that the discrete metric identities (2.8) are satisfied [38].

We now characterize the properties of the DGSEM scheme (3.1) in Lemma 3.1 for the unlimited scheme without artificial viscosity, and in Theorem 3.2 and 3.3 for the scheme with graph viscosity. Lemma 3.1 is not a new result (see, e.g., [71]), but is recalled for the sake of comparison with Theorems 3.2 and 3.3.

**Lemma 3.1.** The scheme (3.1) without artificial viscosity, $d_{\kappa} = 0 \forall \kappa \in \Omega_h$, and with the EC flux (2.13) defined for a given entropy pair $(\vartheta(u), \mathbf{g}(u))$, with $\vartheta''(\cdot) > 0$, satisfies the following discrete entropy inequality for this pair:

$$\langle \vartheta(u_h^{(n+1)}) \rangle_{\kappa} - \langle \vartheta(u_h^{(n)}) \rangle_{\kappa} + \frac{\Delta t^{(n)}}{|\kappa|} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \omega_k J_e^k G_{ec} \left( u_h^{-} (\mathbf{x}_e^k, t^{(n+1)}), u_h^{+} (\mathbf{x}_e^k, t^{(n+1)}), \mathbf{n}_e^k \right) \leq 0$$

$$\text{with } G_{ec}(u^{-}, u^{+}, \mathbf{n}) = \frac{1}{2}h(u^{-}, u^{+}, \mathbf{n})(\vartheta'(u^{-}) + \vartheta'(u^{+})) - \frac{1}{2}(\psi(u^{-}) + \psi(u^{+})) \cdot \mathbf{n}.$$
**Theorem 3.2.** The scheme (3.1) with the EC flux (2.13) defined for a given entropy pair $(\vartheta(u), g(u))$ with $\vartheta''(\cdot) > 0$, and the graph viscosity in (3.2) such that

$$d_{\kappa} \geq 4L_f L_u \max_{k \neq l} \frac{|D_{kl}^p|}{\omega_l} \quad \forall \kappa \in \Omega_h, \quad L_u := \sup_{m \leq u^+ \neq u^- \leq M} \frac{\vartheta''(u^+)(u^+ - u^-)}{v(u^+) - v(u^-)}, \quad (3.4)$$

with $v(u) = \vartheta'(u)$ and $L_f$ defined in (1.4), has a unique solution satisfying the maximum principle: $m \le U_{\kappa}^{ij,n+1} \le M$ for all $\kappa \in \Omega_h$ and $0 \le i, j \le p$.

**Proof. Preliminaries.** We will use short notations for $\mathbf{f}_{\kappa}^{ij} = \mathbf{f}(U_{\kappa}^{ij})$. Let add the following trivial quantity to the LHS of (3.1):

$$-2\omega_i\omega_j J_{\kappa}^{ij} \mathbf{f}_{\kappa}^{ij,n+1} \cdot \left( \sum_{k=0}^{p} D_{ik}^{p} \mathbf{n}_{\kappa}^{(ik)j} + D_{jk}^{p} \mathbf{n}_{\kappa}^{i(jk)} \right) \stackrel{(2.8)}{=} 0$$

with (2.4) and (3.2) and rewrite the result multiplied by $\frac{\Delta t^{(n)}}{\omega_i \omega_j J_{\kappa}^{ij}}$ as

$$U_{\kappa}^{ij,n+1} + \frac{\Delta t^{(n)}}{J_{\kappa}^{ij}} \left( \sum_{k=0}^{p} 2D_{ik}^{p} \left( h_{ec} \left( U_{\kappa}^{ij}, U_{\kappa}^{kj}, \mathbf{n}_{\kappa}^{(ik)j} \right) - \mathbf{f}_{\kappa}^{ij} \cdot \mathbf{n}_{\kappa}^{(ik)j} \right) + \frac{\omega_k d_{\kappa} |\mathbf{n}_{\kappa}^{(ik)j}|}{2} (U_{\kappa}^{ij} - U_{\kappa}^{kj}) \right. \\ + \sum_{k=0}^{p} 2D_{jk}^{p} \left( h_{ec} \left( U_{\kappa}^{ij}, U_{\kappa}^{ik}, \mathbf{n}_{\kappa}^{i(jk)} \right) - \mathbf{f}_{\kappa}^{ij} \cdot \mathbf{n}_{\kappa}^{i(jk)} \right) + \frac{\omega_k d_{\kappa} |\mathbf{n}_{\kappa}^{i(jk)}|}{2} (U_{\kappa}^{ij} - U_{\kappa}^{ik}) \\ \left. + \frac{1}{\omega_i \omega_j} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \phi_{\kappa}^{ij}(\mathbf{x}_e^k) \omega_k J_e^k \left( h(U_{\kappa}^{ij}, u_h^+(\mathbf{x}_e^k, t), \mathbf{n}_e^k) - \mathbf{f}_{\kappa}^{ij} \cdot \mathbf{n}_e^k \right) \right)^{(n+1)} = U_{\kappa}^{ij,n}.$$

**Existence.** In the above scheme, we substitute

$$h(U_{\kappa}^{ij}, u_{h}^{+}(\mathbf{x}_{e}^{k}, t), \mathbf{n}_{e}^{k}) - \mathbf{f}_{\kappa}^{ij} \cdot \mathbf{n}_{e}^{k} = 2L_{f}U_{\kappa}^{ij} - 2L_{f}\left(U_{\kappa}^{ij} - \frac{h(U_{\kappa}^{ij}, u_{h}^{+}(\mathbf{x}_{e}^{k}, t), \mathbf{n}_{e}^{k}) - \mathbf{f}_{\kappa}^{ij} \cdot \mathbf{n}_{e}^{k}}{2L_{f}}\right),$$

and

$$2D_{ik}^{p} \left( h_{ec} \left( U_{\kappa}^{ij}, U_{\kappa}^{kj}, \mathbf{n}_{\kappa}^{(ik)j} \right) - \mathbf{f}_{\kappa}^{ij} \cdot \mathbf{n}_{\kappa}^{(ik)j} \right) + \frac{\omega_k d_{\kappa} |\mathbf{n}_{\kappa}^{(ik)j}|}{2} (U_{\kappa}^{ij} - U_{\kappa}^{kj}) = \frac{d_{\kappa} \omega_k |\mathbf{n}_{\kappa}^{(ik)j}|}{2\beta_{(ik)j}} U_{\kappa}^{ij} \\ - \frac{d_{\kappa} \omega_k |\mathbf{n}_{\kappa}^{(ik)j}|}{2\beta_{(ik)j}} \left( (1 - \beta_{(ik)j}) U_{\kappa}^{ij} + \beta_{(ik)j} U_{\kappa}^{kj} - \frac{h_{ec} \left( U_{\kappa}^{ij}, U_{\kappa}^{kj}, \mathbf{n}_{\kappa}^{(ik)j} \right) - \mathbf{f}_{\kappa}^{ij} \cdot \mathbf{n}_{\kappa}^{(ik)j} \right) \right. \\ \left. = \frac{d_{\kappa} \omega_k |\mathbf{n}_{\kappa}^{(ik)j}|}{2\beta_{(ik)j}} \left( U_{\kappa}^{ij} - U \left( U_{\kappa}^{ij}, U_{\kappa}^{kj}, \frac{\mathbf{n}_{\kappa}^{(ik)j}}{|\mathbf{n}_{\kappa}^{(ik)j}|} \right) \right),$$

with $\beta_{(ik)j} = \beta(U_{\kappa}^{ij}, U_{\kappa}^{kj})$ defined in (2.15a), $U(\cdot, \cdot, \cdot)$ defined in (2.14a) with $2\alpha = \frac{d_{\kappa} \omega_k}{4\beta_{(ik)j} D_{ik}^{p-}}$ which satisfies $\alpha \ge L_f$ from (3.4), $\beta_{(ik)j} \in (0, \frac{1}{2}]$, and $L_u \ge 1$ in (3.4). The latter inequality is due to $\lim_{u^+ \to u^-} v'(u^+) \frac{u^+ - u^-}{v^+ - v^-} = 1$. Now we define $u_h^{(n+1,m=0)} = u_h^{(n)}$ and introduce the following subiterations for $m \ge 0$

$$\left( 1 + \sum_{k=0}^{p} \left( \gamma_{\kappa}^{(ik)j} + \gamma_{\kappa}^{i(jk)} \right) + \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \gamma_{e}^{k} \right)^{(n+1,m)} U_{\kappa}^{ij,(n+1,m+1)} = U_{\kappa}^{ij,n}$$
$$ 
\begin{aligned}
&+ \sum_{k=0}^{p} \left( \gamma_{\kappa}^{(ik)j} \mathcal{U} \left( U_{\kappa}^{ij}, U_{\kappa}^{kj}, \frac{\mathbf{n}_{\kappa}^{(ik)j}}{|\mathbf{n}_{\kappa}^{(ik)j}|} \right) + \gamma_{\kappa}^{i(jk)} \mathcal{U} \left( U_{\kappa}^{ij}, U_{\kappa}^{ik}, \frac{\mathbf{n}_{\kappa}^{(jk)}}{|\mathbf{n}_{\kappa}^{(jk)}|} \right) \right)^{(n+1,m)} \\
&+ \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \gamma_{e}^{k} \left( U_{\kappa}^{ij} - \frac{h(U_{\kappa}^{ij}, u_{h}^{+}(\mathbf{x}_{e}^{k}, t), \mathbf{n}_{e}^{k}) - \mathbf{f}_{\kappa}^{ij} \cdot \mathbf{n}_{e}^{k}}{2L_{f}} \right)^{(n+1,m)}
\end{aligned} 
\quad (3.5) $$

so from (2.14a) and (2.17), $U_{\kappa}^{ij,(n+1,m+1)}$ is a convex combination of quantities in $[m, M]$ with coefficients

$$ \gamma_{\kappa}^{(ik)j} = \frac{\Delta t^{(n)} d_{\kappa} \omega_k |\mathbf{n}_{\kappa}^{(ik)j}|}{2 \beta_{(ik)j}^{n+1,m} J_{\kappa}^{ij}}, \quad \gamma_{\kappa}^{i(jk)} = \frac{\Delta t^{(n)} d_{\kappa} \omega_k |\mathbf{n}_{\kappa}^{i(jk)}|}{2 \beta_{i(jk)}^{n+1,m} J_{\kappa}^{ij}}, \quad \gamma_{e}^{k} = \frac{2 \Delta t^{(n)} L_f \phi_{\kappa}^{ij}(\mathbf{x}_e^k) \omega_k J_e^k}{\omega_i \omega_j J_{\kappa}^{ij}}. \quad (3.6) $$

Since $u_h^{(n+1,m+1)}$ lies in $[m, M]$, (3.5) defines a map $u_h^{(n+1,m)} \mapsto u_h^{(n+1,m+1)}$ of the convex compact $[m, M]^N$, subset of the normed vector space $\mathbb{R}^N$, into itself. This map is continuous by Lipschitz continuity of the interface flux and regularity of the EC flux. By the Schauder fixed point theorem for infinite dimensional spaces [48, 62], this map has at least one fixed point. The fixed point is solution to (3.1).

*Uniqueness.* We here adapt some existing tools ([16], Sect. 26) to our scheme. Let us assume two solutions $u_h^{(n+1)}$ and $w_h^{(n+1)}$ of (3.1) and drop the exponents $(n+1)$ for the sake of clarity, thus we have

$$ \omega_i \omega_j J_{\kappa}^{ij} \frac{U_{\kappa}^{ij} - W_{\kappa}^{ij}}{\Delta t^{(n)}} + R_{\kappa}^{ij}(u_h) + V_{\kappa}^{ij}(u_h) - R_{\kappa}^{ij}(w_h) - V_{\kappa}^{ij}(w_h) = 0 \quad \forall \kappa \in \Omega_h, 0 \le i, j \le p. $$

Using the SBP property (2.6), then symmetry and consistency of $h_{ec}$, we rewrite

$$ 
\begin{aligned}
2\omega_j Q_{ik}^p h_{ec} (U_{\kappa}^{ij}, U_{\kappa}^{kj}, \mathbf{n}_{\kappa}^{(ik)j}) &= \omega_j \left( (Q_{ik}^p - Q_{ki}^p) h_{ec} (U_{\kappa}^{ij}, U_{\kappa}^{kj}, \mathbf{n}_{\kappa}^{(ik)j}) + \delta_{ip}\delta_{kp} \mathbf{f}_{\kappa}^{pj} \cdot \mathbf{n}_{\kappa}^{(pp)j} - \delta_{i0}\delta_{k0} \mathbf{f}_{\kappa}^{0j} \cdot \mathbf{n}_{\kappa}^{(00)j} \right) \\
&= \omega_j \frac{Q_{ik}^p - Q_{ki}^p}{2} h_{ec} (U_{\kappa}^{ij}, U_{\kappa}^{kj}, \mathbf{n}_{\kappa}^{(ik)j}) - \omega_j \frac{Q_{ki}^p - Q_{ik}^p}{2} h_{ec} (U_{\kappa}^{kj}, U_{\kappa}^{ij}, \mathbf{n}_{\kappa}^{(ki)j}) \\
&\quad + \delta_{ip}\delta_{kp}\omega_j J_e(\mathbf{x}_{\kappa}^{pj}) \mathbf{f}_{\kappa}^{pj} \cdot \mathbf{n}_e(\mathbf{x}_{\kappa}^{pj}) + \delta_{i0}\delta_{k0}\omega_j J_e(\mathbf{x}_{\kappa}^{0j}) \mathbf{f}_{\kappa}^{0j} \cdot \mathbf{n}_e(\mathbf{x}_{\kappa}^{0j}),
\end{aligned}
 $$

where we have used (2.10) in the last step. We thus recast the DGSEM residuals (2.4) plus (3.2) as

$$ 
\begin{aligned}
R_{\kappa}^{ij}(u_h) + V_{\kappa}^{ij}(u_h) = \frac{\omega_j}{2} \sum_{k=0}^{p} G_{\kappa}^{(ik)j}(U_{\kappa}^{ij}, U_{\kappa}^{kj}) - G_{\kappa}^{(ik)j}(U_{\kappa}^{kj}, U_{\kappa}^{ij}) &+ \frac{\omega_i}{2} \sum_{k=0}^{p} G_{\kappa}^{i(jk)}(U_{\kappa}^{ij}, U_{\kappa}^{ik}) - G_{\kappa}^{i(jk)}(U_{\kappa}^{ik}, U_{\kappa}^{ij}) \\
&+ \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \phi_{\kappa}^{ij}(\mathbf{x}_e^k) \omega_k J_e^k h(U_{\kappa}^{ij}, u_h^+(\mathbf{x}_e^k, t), \mathbf{n}_e^k),
\end{aligned}
 $$

where

$$ 
\begin{aligned}
G_{\kappa}^{(ik)j}(U_{\kappa}^{ij}, U_{\kappa}^{kj}) &= (Q_{ik}^p - Q_{ki}^p) h_{ec} (U_{\kappa}^{ij}, U_{\kappa}^{kj}, \mathbf{n}_{\kappa}^{(ik)j}) + d_{\kappa} \omega_i \omega_k |\mathbf{n}_{\kappa}^{(ik)j}| (U_{\kappa}^{ij} - U_{\kappa}^{kj}) \\
G_{\kappa}^{(ik)j}(U_{\kappa}^{kj}, U_{\kappa}^{ij}) &= (Q_{ki}^p - Q_{ik}^p) h_{ec} (U_{\kappa}^{kj}, U_{\kappa}^{ij}, \mathbf{n}_{\kappa}^{(ik)j}) + d_{\kappa} \omega_k \omega_j |\mathbf{n}_{\kappa}^{(ik)j}| (U_{\kappa}^{kj} - U_{\kappa}^{ij}) \\
G_{\kappa}^{i(jk)}(U_{\kappa}^{ij}, U_{\kappa}^{ik}) &= (Q_{jk}^p - Q_{kj}^p) h_{ec} (U_{\kappa}^{ij}, U_{\kappa}^{ik}, \mathbf{n}_{\kappa}^{i(jk)}) + d_{\kappa} \omega_j \omega_k |\mathbf{n}_{\kappa}^{i(jk)}| (U_{\kappa}^{ij} - U_{\kappa}^{ik}) \\
G_{\kappa}^{i(jk)}(U_{\kappa}^{ik}, U_{\kappa}^{ij}) &= (Q_{kj}^p - Q_{jk}^p) h_{ec} (U_{\kappa}^{ik}, U_{\kappa}^{ij}, \mathbf{n}_{\kappa}^{i(jk)}) + d_{\kappa} \omega_k \omega_j |\mathbf{n}_{\kappa}^{i(jk)}| (U_{\kappa}^{ik} - U_{\kappa}^{ij})
\end{aligned}
 $$

are skew-symmetric, i.e., $G_{\kappa}^{(ik)j}(U_{\kappa}^{ij}, U_{\kappa}^{kj}) = -G_{\kappa}^{(ik)j}(U_{\kappa}^{kj}, U_{\kappa}^{ij})$, and are nondecreasing functions of their first argument under (3.4). Indeed, let write $G(u) = (Q_{ik}^p - Q_{ki}^p)h_{ec}(u, u^+, \mathbf{n}) + d_{\kappa} \omega_i \omega_k |\mathbf{n}|(u - u^+)$, then $G'(u) =$
$(Q_{ik}^p - Q_{ki}^p) \partial_i h_{ec}(u, u^+, \mathbf{n}) + d_\kappa \omega_i \omega_k |\mathbf{n}|$. Using (2.13) and $\psi'(u) = v'(u)\mathbf{f}(u)^\top$, then $\psi \circ u(v) = \int^v \mathbf{f} \circ u(w) dw$ we obtain

$$ 
\begin{aligned}
|\partial_1 h_{ec}(u, u^+, \mathbf{n})| &= \left| \frac{v'(u)(\psi(u^+) - \psi(u) - (v^+ - v)\mathbf{f}(u)) \cdot \mathbf{n}}{(v^+ - v)^2} \right| \\
&= \frac{v'(u)}{(v^+ - v)^2} \left| (v^+ - v) \left( \int_0^1 \mathbf{f} \circ u(\theta v^+ + (1-\theta)v) \, d\theta - \mathbf{f}(u) \right) \cdot \mathbf{n} \right| \\
&\le \frac{v'(u)}{|v^+ - v|} L_f |u^+ - u| |\mathbf{n}| \le L_f L_u |\mathbf{n}|,
\end{aligned} 
\quad (3.7) $$

from the definition of $L_u$ in (3.4). We rewrite the difference of interface fluxes as

$$ h(U_{\kappa}^{ij}, u_{h}^{+}, \mathbf{n}_{e}^{k}) - h(W_{\kappa}^{ij}, w_{h}^{+}, \mathbf{n}_{e}^{k}) = \frac{1}{2} (h(U_{\kappa}^{ij}, u_{h}^{+}, \mathbf{n}_{e}^{k}) - h(W_{\kappa}^{ij}, u_{h}^{+}, \mathbf{n}_{e}^{k}) + h(W_{\kappa}^{ij}, u_{h}^{+}, \mathbf{n}_{e}^{k}) - h(W_{\kappa}^{ij}, w_{h}^{+}, \mathbf{n}_{e}^{k})) \\ + \frac{1}{2} (h(U_{\kappa}^{ij}, u_{h}^{+}, \mathbf{n}_{e}^{k}) - h(W_{\kappa}^{ij}, w_{h}^{+}, \mathbf{n}_{e}^{k}) + h(U_{\kappa}^{ij}, w_{h}^{+}, \mathbf{n}_{e}^{k}) - h(W_{\kappa}^{ij}, w_{h}^{+}, \mathbf{n}_{e}^{k})), $$

where since $h(\cdot, \cdot, \mathbf{n})$ is nondecreasing with its first argument by monotonicity, the first and last terms in the right-hand side have same sign as $U_{\kappa}^{ij} - W_{\kappa}^{ij}$. Using a triangle inequality, we obtain

$$ 
\begin{aligned}
& \omega_i \omega_j J_{\kappa}^{ij} \frac{|U_{\kappa}^{ij} - W_{\kappa}^{ij}|}{\Delta t^{(n)}} + \frac{\omega_j}{2} \sum_{k=0}^{p} \left| G_{\kappa}^{(ik)j} (U_{\kappa}^{ij}, U_{\kappa}^{kj}) - G_{\kappa}^{(ik)j} (W_{\kappa}^{ij}, U_{\kappa}^{kj}) \right| \\
& \qquad + \frac{\omega_j}{2} \sum_{k=0}^{p} \left| G_{\kappa}^{(ik)j} (U_{\kappa}^{ij}, W_{\kappa}^{kj}) - G_{\kappa}^{(ik)j} (W_{\kappa}^{ij}, W_{\kappa}^{kj}) \right| \\
& \qquad + \frac{\omega_i}{2} \sum_{k=0}^{p} \left| G_{\kappa}^{i(jk)} (U_{\kappa}^{ij}, U_{\kappa}^{ik}) - G_{\kappa}^{i(jk)} (W_{\kappa}^{ij}, U_{\kappa}^{ik}) \right| \\
& \qquad + \frac{\omega_i}{2} \sum_{k=0}^{p} \left| G_{\kappa}^{i(jk)} (U_{\kappa}^{ij}, W_{\kappa}^{ik}) - G_{\kappa}^{i(jk)} (W_{\kappa}^{ij}, W_{\kappa}^{ik}) \right| \\
& \qquad + \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \frac{\phi_{\kappa}^{ij}(\mathbf{x}_e^k) \omega_k J_e^k}{2} \left| h(U_{\kappa}^{ij}, u_h^+, \mathbf{n}_e^k) - h(W_{\kappa}^{ij}, u_h^+, \mathbf{n}_e^k) \right|_{\mathbf{x}_e^k} \\
& \qquad + \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \frac{\phi_{\kappa}^{ij}(\mathbf{x}_e^k) \omega_k J_e^k}{2} \left| h(U_{\kappa}^{ij}, w_h^+, \mathbf{n}_e^k) - h(W_{\kappa}^{ij}, w_h^+, \mathbf{n}_e^k) \right|_{\mathbf{x}_e^k} \\
& \le \frac{\omega_j}{2} \sum_{k=0}^{p} \left| G_{\kappa}^{(ik)j} (W_{\kappa}^{ij}, U_{\kappa}^{kj}) - G_{\kappa}^{(ik)j} (W_{\kappa}^{ij}, W_{\kappa}^{kj}) \right| \\
& \qquad + \frac{\omega_j}{2} \sum_{k=0}^{p} \left| G_{\kappa}^{(ik)j} (U_{\kappa}^{ij}, U_{\kappa}^{kj}) - G_{\kappa}^{(ik)j} (U_{\kappa}^{ij}, W_{\kappa}^{kj}) \right| \\
& \qquad + \frac{\omega_i}{2} \sum_{k=0}^{p} \left| G_{\kappa}^{i(jk)} (W_{\kappa}^{ij}, U_{\kappa}^{ik}) - G_{\kappa}^{i(jk)} (W_{\kappa}^{ij}, W_{\kappa}^{ik}) \right| \\
& \qquad + \frac{\omega_i}{2} \sum_{k=0}^{p} \left| G_{\kappa}^{i(jk)} (U_{\kappa}^{ij}, U_{\kappa}^{ik}) - G_{\kappa}^{i(jk)} (U_{\kappa}^{ij}, W_{\kappa}^{ik}) \right| \\
& \qquad + \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \frac{\phi_{\kappa}^{ij}(\mathbf{x}_e^k) \omega_k J_e^k}{2} \left| h(U_{\kappa}^{ij}, u_h^+, \mathbf{n}_e^k) - h(U_{\kappa}^{ij}, w_h^+, \mathbf{n}_e^k) \right|_{\mathbf{x}_e^k}
\end{aligned}
$$
$$+ \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \frac{\phi_{\kappa}^{ij}(\mathbf{x}_e^k) \omega_k J_e^k}{2} |h(W_{\kappa}^{ij}, u_h^+, \mathbf{n}_e^k) - h(W_{\kappa}^{ij}, w_h^+, \mathbf{n}_e^k)|_{\mathbf{x}_e^k}$$

since all the terms between vertical bars in the LHS are of the same sign. Multiplying the above relation by a quantity $\phi_{\kappa}^{ij} \ge 0$, which will be defined below, and summing over $\kappa \in \Omega_h$ and $0 \le i, j \le p$, then switching indices $i \leftrightarrow k$ and $j \leftrightarrow k$ in the sums over cells $\kappa$, switching the left and right states in the sums over faces $e$ and using conservation to transform

$$\frac{\phi_{\kappa}^{ij}(\mathbf{x}_e^k) \omega_k J_e^k}{2} |h(U_{\kappa}^{ij}, u_h^+, \mathbf{n}_e^k) - h(U_{\kappa}^{ij}, w_h^+, \mathbf{n}_e^k)|_{\mathbf{x}_e^k} \phi_{\kappa}^{ij} = \frac{\omega_k J_e^k}{2} |h(u_h^+, u_h^-, -\mathbf{n}_e^k) - h(w_h^+, w_h^-, -\mathbf{n}_e^k)|_{\mathbf{x}_e^k} \phi_{h}^{-}(\mathbf{x}_e^k),$$

for which the equivalent contribution from the opposite state reads $\frac{\omega_k J_e^k}{2} |h(u_h^-, u_h^+, \mathbf{n}_e^k) - h(u_h^-, w_h^+, \mathbf{n}_e^k)|_{\mathbf{x}_e^k} \phi_{h}^{+}(\mathbf{x}_e^k)$, we finally get

$$
\begin{align}
\sum_{\kappa \in \Omega_h} \sum_{i,j=0}^{p} a_{\kappa}^{ij} |U_{\kappa}^{ij} - W_{\kappa}^{ij}| & \le \sum_{\kappa \in \Omega_h} \sum_{i,j,k=0}^{p} \frac{\omega_j}{2} \left| G_{\kappa}^{(ik)j} (W_{\kappa}^{ij}, U_{\kappa}^{kj}) - G_{\kappa}^{(ik)j} (W_{\kappa}^{ij}, W_{\kappa}^{kj}) \right| |\varphi_{\kappa}^{ij} - \varphi_{\kappa}^{kj}| \nonumber \\
& \quad + \sum_{\kappa \in \Omega_h} \sum_{i,j,k=0}^{p} \frac{\omega_j}{2} \left| G_{\kappa}^{(ik)j} (U_{\kappa}^{ij}, U_{\kappa}^{kj}) - G_{\kappa}^{(ik)j} (U_{\kappa}^{ij}, W_{\kappa}^{kj}) \right| |\varphi_{\kappa}^{ij} - \varphi_{\kappa}^{kj}| \nonumber \\
& \quad + \sum_{\kappa \in \Omega_h} \sum_{i,j,k=0}^{p} \frac{\omega_i}{2} \left| G_{\kappa}^{i(jk)} (W_{\kappa}^{ij}, U_{\kappa}^{ik}) - G_{\kappa}^{i(jk)} (W_{\kappa}^{ij}, W_{\kappa}^{ik}) \right| |\varphi_{\kappa}^{ij} - \varphi_{\kappa}^{ik}| \nonumber \\
& \quad + \sum_{\kappa \in \Omega_h} \sum_{i,j,k=0}^{p} \frac{\omega_i}{2} \left| G_{\kappa}^{i(jk)} (U_{\kappa}^{ij}, U_{\kappa}^{ik}) - G_{\kappa}^{i(jk)} (U_{\kappa}^{ij}, W_{\kappa}^{ik}) \right| |\varphi_{\kappa}^{ij} - \varphi_{\kappa}^{ik}| \nonumber \\
& \quad + \sum_{\kappa \in \Omega_h} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \frac{\omega_k J_e^k}{2} \left( |h(u_h^-, u_h^+, \mathbf{n}_e^k) - h(w_h^-, u_h^+, \mathbf{n}_e^k)| \right. \nonumber \\
& \quad \left. + |h(u_h^-, w_h^+, \mathbf{n}_e^k) - h(w_h^-, w_h^+, \mathbf{n}_e^k)| \right)_{\mathbf{x}_e^k} |\varphi_h^+ - \varphi_h^-|_{\mathbf{x}_e^k} \nonumber \\
& \le \sum_{\kappa \in \Omega_h} \sum_{i,j,k=0}^{p} \omega_j \left| \mathbf{n}_{\kappa}^{(ik)j} \right| \left( |Q_{ik}^p - Q_{ki}^p| L_f L_u + d_{\kappa} \omega_i \omega_k \right) |U_{\kappa}^{kj} - W_{\kappa}^{kj}| |\varphi_{\kappa}^{ij} - \varphi_{\kappa}^{kj}| \nonumber \\
& \quad + \sum_{\kappa \in \Omega_h} \sum_{i,j,k=0}^{p} \omega_i \left| \mathbf{n}_{\kappa}^{i(jk)} \right| \left( |Q_{jk}^p - Q_{kj}^p| L_f L_u + d_{\kappa} \omega_k \omega_j \right) |U_{\kappa}^{ik} - W_{\kappa}^{ik}| |\varphi_{\kappa}^{ij} - \varphi_{\kappa}^{ik}| \nonumber \\
& \quad + \sum_{\kappa \in \Omega_h} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \frac{\omega_k J_e^k}{2} (L_1 + L_2) |u_h^- - w_h^-|_{\mathbf{x}_e^k} |\varphi_h^+ - \varphi_h^-|_{\mathbf{x}_e^k} \nonumber \\
& = \sum_{\kappa \in \Omega_h} \sum_{i,j=0}^{p} b_{\kappa}^{ij} |U_{\kappa}^{ij} - W_{\kappa}^{ij}| \tag{3.8}
\end{align}
$$

where we have used symmetry of $h_{ec}(\cdot, \cdot, \mathbf{n})$ and the mean value theorem to get

$$
\begin{align*}
|h_{ec}(u, u^{+}, \mathbf{n}) - h_{ec}(u, u^{-}, \mathbf{n})| &= |h_{ec}(u^{+}, u, \mathbf{n}) - h_{ec}(u^{-}, u, \mathbf{n})| \\
&= |\partial_1 h_{ec}(u, \theta u^{+} + (1 - \theta)u^{-}, \mathbf{n})| |u^{+} - u^{-}| \stackrel{(3.7)}{\le} L_f L_u |\mathbf{n}| |u^{+} - u^{-}|,
\end{align*}
$$

while by conservation and Lipschitz continuity we have

$$|h(u_h^-, w_h^+, \mathbf{n}_e^k) - h(w_h^-, w_h^+, \mathbf{n}_e^k)| = |h(w_h^+, u_h^-, -\mathbf{n}_e^k) - h(w_h^+, w_h^-, -\mathbf{n}_e^k)| \le L_2 |u_h^- - w_h^-|.$$
The coefficients in (3.8) read $a_{\kappa}^{ij} = \frac{\omega_i \omega_j J_{\kappa}^{ij}}{\Delta t^{(n)}} \varphi_{\kappa}^{ij}$, and from (2.10) we obtain

$$b_{\kappa}^{ij} = \sum_{k=0}^{p} \omega_j |\mathbf{n}_{\kappa}^{(ik)j}| (|Q_{ik}^{p} - Q_{ki}^{p}| L_f L_u + d_{\kappa} \omega_i \omega_k) |\varphi_{\kappa}^{ij} - \varphi_{\kappa}^{kj}| \\ + \sum_{k=0}^{p} \omega_i |\mathbf{n}_{\kappa}^{i(jk)}| (|Q_{jk}^{p} - Q_{kj}^{p}| L_f L_u + d_{\kappa} \omega_k \omega_j) |\varphi_{\kappa}^{ij} - \varphi_{\kappa}^{ik}| \\ + \frac{L_1+L_2}{2} (\omega_j (\delta_{i0} |\mathbf{n}_{\kappa}^{(00)j}| + \delta_{ip} |\mathbf{n}_{\kappa}^{(pp)j}|) + \omega_i (\delta_{j0} |\mathbf{n}_{\kappa}^{i(00)}| + \delta_{jp} |\mathbf{n}_{\kappa}^{i(pp)}|)) |\varphi_{\kappa}^{ij} - \varphi_{h}^{+} (\mathbf{x}_{\kappa}^{ij})|. $$

We now define $\varphi_h$ as the space average of $\varphi(\mathbf{x}) = e^{-\gamma|\mathbf{x}|}$ in subcells: $U_{i,j=0}^p \overline{\kappa_{ij}} = \overline{\kappa}: \varphi_h(\mathbf{x}) = \frac{1}{|\kappa_{ij}|} \int_{\kappa_{ij}} \varphi(\mathbf{x}) \, dV$ $\forall \mathbf{x} \in \kappa_{ij}$. Such subcells $\kappa_{ij}$ of volume $|\kappa_{ij}| = \omega_i \omega_j J_{\kappa}^{ij}$ are associated to the Gauss-Lobatto nodes, $\mathbf{n}_{\kappa}^{ij}$, and we refer to Appendix B.3 of [30] for their precise definition. We thus have $\sum_{\kappa} \langle \varphi_h \rangle_{\kappa} \le \int_{\mathbb{R}^d} e^{-\gamma|\mathbf{x}|} dV$ finite providing $\gamma > 0$ which will be adequately chosen later on to prove that necessarily $u_h \equiv w_h$. The choice for $\varphi_h$ is not unique, but we need the series in (3.8) to converge and the inequalities in (3.9) below to hold.

Using (2.11) and $|Q_{ik}^p - Q_{ki}^p| L_f L_u \le \omega_k \omega_i d_{\kappa}$ from (3.4), we have

$$a_{\kappa}^{ij} \geq \frac{\omega_0^d \beta h_{\Omega}^d}{2d \Delta t^{(n)}} \inf_{\mathbf{x} \in B_{\kappa}(h_{\Omega})} \varphi(\mathbf{x}), \quad b_{\kappa}^{ij} \leq \frac{(8d_{\kappa} + L_1 + L_2)h_{\Omega}^d}{d\beta} \sup_{\mathbf{x} \in B_{\kappa}(2h_{\Omega})} |\nabla \varphi|, \quad (3.9)$$

with $B_{\kappa}(h_{\Omega})$ the ball of radius $h_{\Omega}$ centered on the barycenter of $\kappa$. Therefore choosing $\gamma > 0$ small enough such that $\gamma e^{3\gamma h} < \frac{\omega_0^d \beta^2}{2(8d_{\kappa}+L_1+L_2)\Delta t^{(n)}}$ imposes $a_{\kappa}^{ij} > b_{\kappa}^{ij}$, so inequality (3.8) cannot hold unless $U_{\kappa}^{ij} = W_{\kappa}^{ij}$ for all the DOFs, which completes the proof. $\square$

**Theorem 3.3.** The scheme (3.1) with the EC flux (2.13), defined for a given entropy pair $(\vartheta(u), \mathbf{g}(u))$ with $\vartheta''(\cdot) > 0$, and the artificial viscosity in (3.2) and (3.4) satisfies the following inequality for every entropy pair $(\eta(u), \mathbf{q}(u))$ with $\eta''(\cdot) \ge 0$:

$$\langle \eta(u_h^{(n+1)}) \rangle_{\kappa} - \langle \eta(u_h^{(n)}) \rangle_{\kappa} + \frac{\Delta t^{(n)}}{|\kappa|} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \omega_k J_e^k Q(u_h^{-}(\mathbf{x}_e^k, t^{(n+1)}), u_h^{+}(\mathbf{x}_e^k, t^{(n+1)}), \mathbf{n}_e^k) \le 0,$$

with the numerical entropy flux $Q(\cdot, \cdot, \cdot)$ defined in (2.18).

*Proof.* We again use short notations as $\eta_{\kappa}^{ij} = \eta(U_{\kappa}^{ij})$ and $\mathbf{q}_{\kappa}^{ij} = \mathbf{q}(U_{\kappa}^{ij})$, and drop the time exponent $(n+1)$ for the sake of clarity. We know that $u_h^{(n+1,m+1)} = u_h^{(n+1,m)} = u_h^{(n+1)}$ is solution to (3.5). We apply the entropy $\eta(\cdot)$ to the convex combination in (3.5) and use the Jensen's inequality, then pass to the limit to obtain an inequality for $u_h^{(n+1)} = \lim_{m \to \infty} u_h^{(n+1,m)}$. Using (2.18) and (2.14b), we get

$$\left(1 + \sum_{k=0}^{p} (\gamma_{\kappa}^{(ik)j} + \gamma_{\kappa}^{i(jk)}) + \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \gamma_e^k\right) \eta_{\kappa}^{ij} \le \eta_{\kappa}^{ij,n} \\ + \sum_{k=0}^{p} \gamma_{\kappa}^{(ik)j} \left( \frac{\bar{\eta}(U_{\kappa}^{ij}, U_{\kappa}^{kj}) + \eta_{\kappa}^{ij}}{2} - \frac{q_{ec}(U_{\kappa}^{ij}, U_{\kappa}^{kj}, \mathbf{n}_{\kappa}^{(ik)j}) - \mathbf{q}_{\kappa}^{ij} \cdot \mathbf{n}_{\kappa}^{(ik)j}}{\frac{d_{\kappa} \omega_k |\mathbf{n}_{\kappa}^{(ik)j}|}{4\beta_{(ik)j} D_{ik}^p}} \right) \\ + \sum_{k=0}^{p} \gamma_{\kappa}^{i(jk)} \left( \frac{\bar{\eta}(U_{\kappa}^{ij}, U_{\kappa}^{ik}) + \eta_{\kappa}^{ij}}{2} - \frac{q_{ec}(U_{\kappa}^{ij}, U_{\kappa}^{ik}, \mathbf{n}_{\kappa}^{(jk)}) - \mathbf{q}_{\kappa}^{ij} \cdot \mathbf{n}_{\kappa}^{(jk)j}}{\frac{d_{\kappa} \omega_k |\mathbf{n}_{\kappa}^{(jk)j}|}{4\beta_{(jk)j} D_{jk}^p}} \right) \\ + \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \gamma_e^k \left( \eta_{\kappa}^{ij} - \frac{Q(U_{\kappa}^{ij}, u_h^{+}(\mathbf{x}_e^k, t), \mathbf{n}_e^k) - \mathbf{q}_{\kappa}^{ij} \cdot \mathbf{n}_e^k}{2L_f} \right),$$
with $\bar{\eta}$ defined in (2.15b). Using the metric identities (2.8), the terms in $\mathbf{q}_{\kappa}^{ij}$ in the second and third lines vanish, then using the expressions of the convex coefficients (3.6) and multiplying the result by $\omega_i \omega_j J_{\kappa}^{ij} / \Delta t^{(n)}$ gives

$$
\begin{align*}
& \omega_i \omega_j J_{\kappa}^{ij} \frac{\eta_{\kappa}^{ij} - \eta_{\kappa}^{ij,n}}{\Delta t^{(n)}} + 2 \sum_{k=0}^{p} \omega_j Q_{ik}^{p} q_{ec} (U_{\kappa}^{ij}, U_{\kappa}^{kj}, \mathbf{n}_{\kappa}^{(ik)j}) + \omega_i Q_{jk}^{p} q_{ec} (U_{\kappa}^{ij}, U_{\kappa}^{ik}, \mathbf{n}_{\kappa}^{(jk)j}) \\
& \qquad + \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \phi_{\kappa}^{ij}(\mathbf{x}_e^k) \omega_k J_e^k (Q(U_{\kappa}^{ij}, u_h^+(\mathbf{x}_e^k, t), \mathbf{n}_e^k) - \mathbf{q}_{\kappa}^{ij} \cdot \mathbf{n}_e^k) \\
& \qquad \leq d_{\kappa} \omega_i \omega_j \sum_{k=0}^{p} \frac{\omega_k}{2} \left( \frac{|\mathbf{n}_{\kappa}^{(ik)j}|}{\beta_{(ik)j}} (\bar{\eta}(U_{\kappa}^{ij}, U_{\kappa}^{kj}) - \eta_{\kappa}^{ij}) + \frac{|\mathbf{n}_{\kappa}^{(jk)j}|}{\beta_{(jk)j}} (\bar{\eta}(U_{\kappa}^{ij}, U_{\kappa}^{ik}) - \eta_{\kappa}^{ij}) \right) \\
& \qquad \stackrel{(2.16)}{\leq} d_{\kappa} \omega_i \omega_j \sum_{k=0}^{p} \omega_k \left( |\mathbf{n}_{\kappa}^{(ik)j}| (\eta_{\kappa}^{kj} - \eta_{\kappa}^{ij}) + |\mathbf{n}_{\kappa}^{(jk)j}| (\eta_{\kappa}^{ik} - \eta_{\kappa}^{ij}) \right).
\end{align*}
$$

The LHS corresponds to a time implicit DGSEM discretization of (1.2). By conservation of the DGSEM, summing the above inequality over $0 \le i, j \le p$, the $q_{ec}(\cdot, \cdot, \cdot)$ terms in the first row of the LHS cancel out with the $\mathbf{q}_{\kappa}^{ij}$ terms in the second row from symmetry of $q_{ec}(\cdot, \cdot, \mathbf{n})$ in (2.15c) and the SBP property (2.6). The RHS vanishes since $\sum_{i,j,k=0}^{p} \omega_i \omega_j \omega_k |\mathbf{n}_{\kappa}^{(ik)j}| (\eta_{\kappa}^{kj} - \eta_{\kappa}^{ij}) = 0$ by switching indices $i \leftrightarrow k$ and using symmetry of $\mathbf{n}_{\kappa}^{(ik)j} = \mathbf{n}_{\kappa}^{(ki)j}$ in (2.9). The same holds for the second sum in the RHS which concludes the proof. $\square$

**Remark 3.4.** Properties of Theorems 3.2 and 3.3 also hold without EC fluxes, i.e., with $h_{ec}(u^{-}, u^{+}, \mathbf{n}) = \frac{1}{2}(\mathbf{f}(u^{-}) + \mathbf{f}(u^{+})) \cdot \mathbf{n}$ in (3.1) and (2.4), but we look for graph viscosity (3.2) as a minimal modification of the scheme to make it MPP and ES and want to keep the scheme as close as it is used in practice. The modification (3.2) to the DGSEM scheme is indeed local to the cell, easy to implement and to differentiate. Likewise, $d_{\kappa}$ can be easily computed once at the beginning of the simulation (see Sect. 6).

**Remark 3.5 (Square entropy).** Conditions (3.4) may be relaxed for the square entropy, $\vartheta(u) = \frac{u^2}{2}$ (see Rem. 2.3), since $L_u = 1$, $\beta(u^{-}, u^{+}) = \frac{1}{4}$, and $|h_{ec}(u^{+}, u, \mathbf{n}) - h_{ec}(u^{-}, u, \mathbf{n})| \le L_f \int_{0}^{1} \theta d\theta |u^{+} - u^{-}|$, so $d_{\kappa} \ge 2L_f \max_{k \ne l} \frac{|D_{kt}^{p}|}{\omega_k}$ is enough.

## 4. SPACE-TIME DGSEM DISCRETIZATION

We now look for a numerical solution to (1.1) in the function space $V_h^{p,q}$ of piecewise polynomials of degrees $p \ge 1$ in space (here $d=2$) and $q \ge 1$ in time:

$$ u_h(\mathbf{x}, t) = \sum_{i,j=0}^{p} \sum_{r=0}^{q} \phi_{\kappa}^{ij}(\mathbf{x}) \phi_{n}^{r}(t) U_{\kappa}^{ij}(t_n^r), \quad \forall (\mathbf{x}, t) \in \kappa \times (t^{(n)}, t^{(n+1)}), $$

with $\phi_n^r(t) = \ell_r(2\frac{t-t^{(n)}}{\Delta t^{(n)}} - 1)$, with again the Gauss-Lobatto quadrature rules and associated Lagrange polynomials as basis functions in both space and time. Hence, every DOF $U_{\kappa}^{ij}(t_n^r)$ corresponds to $u_h(\mathbf{x}_{\kappa}^{ij}, t_n^r)$ with $t_n^r = t^{(n)} + \frac{\Delta t^{(n)}}{2}(1 + \xi_r)$, so $t_n^0 = t^{(n)+}$ and $t_n^q = t^{(n+1)-}$. The space-time DGSEM reads: find $u_h \in V_h^{p,q}$ such that for all $\kappa \in \Omega_h$, $0 \le i, j \le p$ and $0 \le r \le q$ we have

$$ \omega_i \omega_j J_{\kappa}^{ij} T_{\kappa}^{r}(u_h(\mathbf{x}_{\kappa}^{ij}, \cdot)) + \frac{\omega_r \Delta t^{(n)}}{2} (R_{\kappa}^{ij}(u_h(\cdot, t_n^r)) + V_{\kappa}^{ij}(u_h(\cdot, t_n^r))) = 0, \quad (4.1) $$

together with the initial condition $U_{\kappa}^{ij}(t_{-1}^q) = u_0(x_{\kappa}^{ij})$. With some slight abuse, we remove the dependence in $p$ and $q$ of the quadrature weights, but indices $r, m$ (resp., $i, j, k, l$) will always refer to indices in the range $[0, q]$.
(resp., $[0, p]$). The space discretization is defined in (2.4) and (3.2), while the time discretization reads

$$ T_{\kappa}^{r}(u_{h}(\mathbf{x}_{\kappa}^{ij}, \cdot)) = \sum_{m=0}^{q} 2Q_{rm}^{q} U_{ec}(U_{\kappa}^{ij}(t_{n}^{r}), U_{\kappa}^{ij}(t_{n}^{m})) + \delta_{r0}(U_{\kappa}^{ij}(t_{n}^{0}) - U_{\kappa}^{ij}(t_{n-1}^{q})) \\ + d_{n}\omega_{r}\sum_{m=0}^{q} \omega_{m}(U_{\kappa}^{ij}(t_{n}^{r}) - U_{\kappa}^{ij}(t_{n}^{m})), \qquad (4.2) $$

where we have used an upwind numerical flux in time to satisfy causality, while the last term is again a graph viscosity term and

$$ U_{ec}(u^{-}, u^{+}) = \int_{0}^{1} u(\theta v(u^{+}) + (1 - \theta)v(u^{-})) d\theta \stackrel{(2.15a)}{=} 2\beta(u^{-}, u^{+})(u^{+} - u^{-}) + u^{-}, \qquad (4.3a) $$

with $v(u) = \vartheta'(u)$, is smooth, consistent, $U_{ec}(u, u) = u$, and EC in the sense [18, 49]

$$ (v(u^{+}) - v(u^{-}))U_{ec}(u^{-}, u^{+}) = \psi_{t}(u^{+}) - \psi_{t}(u^{-}), \quad \psi_{t}(u) = v(u)u - \vartheta(u), \qquad (4.3b) $$

with $\psi_t(u)$ the entropy potential. By (4.3b), we have $\partial_1 U_{ec}(u^{-}, u^{+}) = \frac{v'(u^{-})}{(v^{+} - v^{-})^2} (\vartheta^{-} - \vartheta^{+} - v^{+}(u^{-} - u^{+})) \in [0, L_u]$ by convexity of $\vartheta(u)$, $v^{-}(u^{-} - u^{+}) \ge \vartheta^{-} - \vartheta^{+} \ge v^{+}(u^{-} - u^{+})$, and the definition of $L_u$ in (3.4). By symmetry, $\partial_2 U_{ec}(u^{-}, u^{+}) = \partial_1 U_{ec}(u^{+}, u^{-})$ with $\partial_1$ and $\partial_2$ the derivatives with respect to the first and second arguments, respectively.

Summing (4.1) over $0 \le r \le q$ and $0 \le i, j \le p$, we obtain

$$ \langle u_h(t_n^q) \rangle_\kappa = \langle u_h(t_{n-1}^q) \rangle_\kappa - \frac{\Delta t^{(n)}}{|\kappa|} \sum_{r=0}^{q} \frac{\omega_r}{2} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \omega_k J_e^k h(u_h^-(\mathbf{x}_e^k, t_n^r), u_h^+(\mathbf{x}_e^k, t_n^r), \mathbf{n}_e^k) \qquad (4.4) $$

which may be compared to (3.3). The following lemma has been established in Theorems 1 & 2 of [18] and is given here for the sake of comparison with Theorem 4.3.

**Lemma 4.1.** The space-time DGSEM (4.1) without artificial viscosity, $d_\kappa = d_n = 0 \forall \kappa \in \Omega_h$ and $n \ge 0$, and with the EC fluxes (2.13) and (4.3) defined for a given entropy pair $(\vartheta(u), \mathbf{g}(u))$, with $\vartheta''(\cdot) > 0$, satisfies the following discrete entropy inequality for this pair:

$$ \langle \vartheta(u_h(t_n^q)) \rangle_\kappa - \langle \vartheta(u_h(t_{n-1}^q)) \rangle_\kappa + \frac{\Delta t^{(n)}}{|\kappa|} \sum_{r=0}^{q} \frac{\omega_r}{2} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \omega_k J_e^k G_{ec}(u_h^-(\mathbf{x}_e^k, t_n^r), u_h^+(\mathbf{x}_e^k, t_n^r), \mathbf{n}_e^k) \le 0. $$

$$ \text{with } G_{ec}(u^{-}, u^{+}, \mathbf{n}) = \frac{1}{2}h(u^{-}, u^{+}, \mathbf{n})(\vartheta'(u^{-}) + \vartheta'(u^{+})) - \frac{1}{2}(\psi(u^{-}) + \psi(u^{+})) \cdot \mathbf{n}. $$

**Theorem 4.2.** The space-time DGSEM (4.1) with the EC fluxes (2.13) and (4.3) defined for a given entropy pair $(\vartheta(u), \mathbf{g}(u))$ with $\vartheta''(\cdot) > 0$, and the graph viscosities in (3.2) and (4.2) such that

$$ d_{\kappa} \ge 4L_{f}L_{u} \max_{k \ne l} \frac{|D_{kl}^{p}|}{\omega_{l}} \quad \forall \kappa \in \Omega_{h}, \quad d_{n} \ge 2L_{u} \max_{r \ne m} \frac{|D_{rm}^{q}|}{\omega_{m}} \quad \forall n \ge 0, \qquad (4.5) $$

with $L_f$ defined in (1.4) and $L_u$ defined in (3.4), has a unique solution satisfying the maximum principle: $m \le U_{\kappa}^{ij}(t_n^r) \le M$ for all $\kappa \in \Omega_h$, $0 \le i, j \le p$, and $0 \le r \le q$.

*Proof. Preliminaries.* We here extend to (4.1) techniques used in the proof of Theorem 3.2 for the scheme (3.1). We will again use notations like $U_n^r = U_{\kappa}^{ij}(t_n^r)$ when there is no ambiguity. Let us add the following regularization to $T_{\kappa}^r(u_h)$ in (4.2)

$$ \nu\omega_r\omega_0(U_{\kappa}^{ij}(t_n^r) - U_{\kappa}^{ij}(t_{n-1}^q)), $$
with $\nu > 0$ and we will then pass to the limit $\nu \to 0$ (see last part of the proof). This term extends the stencil of the graph viscosity in $T_{\kappa}^r(u_h)$ and will be used only to prove uniqueness, while existence and stability will be shown to hold whatever $\nu \ge 0$.

**Existence.** Using (4.3a), we rewrite the time discretization terms (4.2) in (4.1) as

$$ T_{\kappa}^r = \sum_{m=0}^{q} \omega_r \left( (d_n\omega_m + 2D_{rm}^q(1-2\beta_{rm}))U_n^r - (d_n\omega_m - 4D_{rm}^q\beta_{rm})U_n^m \right) + \delta_{r0}(U_n^0 - U_{n-1}^q) + \nu\omega_r\omega_0(U_n^r - U_{n-1}^q) \\ \stackrel{(2.7)}{=} \sum_{m=0}^{q} \gamma_n^{rm}(U_n^r - U_n^m) + (\delta_{r0} + \nu\omega_r\omega_0)(U_n^r - U_{n-1}^q), \quad (4.6) $$

with $\gamma_n^{rm} = \omega_r(d_n\omega_m - 4D_{rm}^q\beta_{rm}) \ge 0$ from (4.5), $\beta_{rm} = \beta(U_n^r, U_n^m) \in (0, \frac{1}{2}]$ defined in (2.15a), and $L_u \ge 1$ in (3.4).

We now rewrite the space discretization terms, $R_{\kappa}^{ij} + V_{\kappa}^{ij}$, in (4.1) as in (3.5) and using (4.6) we again introduce the following subiterations for (4.1) divided by $\omega_i\omega_j J_{\kappa}^{ij}$:

$$ \left( \delta_{r0} + \nu\omega_r\omega_0 + \sum_{m=0}^{q} \gamma_n^{rm} + \frac{\omega_r}{2} \left( \sum_{k=0}^{p} \left( \gamma_{\kappa}^{(ik)j} + \gamma_{\kappa}^{i(jk)} \right) + \sum_{e \in \partial\kappa} \sum_{k=0}^{p} \gamma_e^k \right) \right)^{(m)} U_{\kappa}^{ij}(t_n^r)^{(m+1)} = (\delta_{r0} + \nu\omega_r\omega_0) U_{\kappa}^{ij}(t_{n-1}^q) \\ + \sum_{m=0}^{q} \gamma_n^{rm} U_{\kappa}^{ij}(t_n^m)^{(m)} + \frac{\omega_r}{2} \sum_{k=0}^{p} \left( \gamma_{\kappa}^{(ik)j} U \left( U_{\kappa}^{ij}, U_{\kappa}^{kj}, \frac{\mathbf{n}_{\kappa}^{(ik)j}}{|\mathbf{n}_{\kappa}^{(ik)j}|} \right) + \gamma_{\kappa}^{i(jk)} U \left( U_{\kappa}^{ij}, U_{\kappa}^{ik}, \frac{\mathbf{n}_{\kappa}^{i(jk)}}{|\mathbf{n}_{\kappa}^{i(jk)}|} \right) \right)^{(m)} \\ + \frac{\omega_r}{2} \sum_{e \in \partial\kappa} \sum_{k=0}^{p} \gamma_e^k \left( U_{\kappa}^{ij} - \frac{h(U_{\kappa}^{ij}, u_h^+(\mathbf{x}_e^k, t), \mathbf{n}_e^k) - \mathbf{f}_{\kappa}^{ij} \cdot \mathbf{n}_e^k}{2L_f} \right)^{(m)} \quad (4.7) $$

with $\mathbf{m} \ge 0$, the initial data $U_{\kappa}^{ij}(t_n^{0 \le r \le q})^{(m=0)} = U_{\kappa}^{ij}(t_{n-1}^q)$, and the coefficients defined above and in (3.6). The iteration is a convex combination of quantities in $[m, M]$ whatever $\nu \ge 0$ and one may again use the Schauder fixed point theorem to show that the map $u_h(t_n^{0 \le r \le q})^{(m)} \mapsto u_h(t_n^{0 \le r \le q})^{(m+1)}$ defined by (4.7) has at least one fixed point with values in $[m, M]$. This fixed point is solution to (4.1).

**Uniqueness.** We again assume two solutions $u_h$ and $w_h$ of (4.1) with the same initial data at time $t_{n-1}^q$ and show that necessarily $u_h \equiv w_h$. Subtracting both schemes gives

$$ \omega_i \omega_j J_{\kappa}^{ij} (T_{\kappa}^r(u_h) - T_{\kappa}^r(w_h))_{\mathbf{x}_{\kappa}^{ij}} + \frac{\omega_r \Delta t^{(n)}}{2} (R_{\kappa}^{ij}(u_h) + V_{\kappa}^{ij}(u_h) - R_{\kappa}^{ij}(w_h) - V_{\kappa}^{ij}(w_h))_{t_n^r} = 0. $$

Using SBP (2.6), we rewrite the time discretization terms as

$$ T_{\kappa}^{r}(u_h) = \delta_{rq}U_{n}^{q} - \delta_{r0}U_{n-1}^{q} + \nu\omega_{r}\omega_{0}(U_{n}^{r} - U_{n-1}^{q}) + \sum_{m=0}^{q} (Q_{rm}^{q} - Q_{mr}^{q})U_{ec}(U_{n}^{r}, U_{n}^{m}) + d_{n}\omega_{r}\omega_{m}(U_{n}^{r} - U_{n}^{m}), $$

and obtain

$$ \begin{aligned} T_{\kappa}^{r}(u_h) - T_{\kappa}^{r}(w_h) = & (\delta_{rq} + \nu\omega_r\omega_0)(U_n^r - W_n^r) \\ & + \sum_{m=0}^{q} Q_{rm}^q (U_{ec}(U_n^r, U_n^m) - U_{ec}(W_n^r, U_n^m)) + \frac{1}{2} d_n \omega_r \omega_m (U_n^r - W_n^r) \\ & + \sum_{m=0}^{q} Q_{rm}^q (U_{ec}(W_n^r, U_n^m) - U_{ec}(W_n^r, W_n^m)) - \frac{1}{2} d_n \omega_r \omega_m (U_n^m - W_n^m) \\ & - \sum_{m=0}^{q} Q_{mr}^q (U_{ec}(U_n^r, U_n^m) - U_{ec}(U_n^r, W_n^m)) + \frac{1}{2} d_n \omega_r \omega_m (U_n^m - W_n^m) \end{aligned} $$
$$ -\sum_{m=0}^{q} Q_{mr}^{q}(U_{ec}(U_{n}^{r}, W_{n}^{m}) - U_{ec}(W_{n}^{r}, W_{n}^{m})) - \frac{1}{2}d_{n}\omega_{r}\omega_{m}(U_{n}^{r} - W_{n}^{r}). $$

Since $\partial_1 U_{ec}(u^{-}, u^{+}) \in [0, L_u]$, the first, second and last lines in the RHS have the same sign under (4.5), we can thus combine this result together with (3.8) to get

$$ \sum_{\kappa \in \Omega_h} \sum_{i,j=0}^{p} \sum_{m=0}^{q} a_{\kappa}^{ijr} |U_{\kappa}^{ij}(t_n^r) - W_{\kappa}^{ij}(t_n^r)| \le \sum_{\kappa \in \Omega_h} \sum_{i,j=0}^{p} \sum_{m=0}^{q} b_{\kappa}^{ijr} |U_{\kappa}^{ij}(t_n^r) - W_{\kappa}^{ij}(t_n^r)| $$

with

$$ a_{\kappa}^{ijr} = (\delta_{rq} + \nu\omega_{r}\omega_{0})a_{\kappa}^{ij}(t_n^r), \quad b_{\kappa}^{ijr} = \frac{\omega_r}{2}b_{\kappa}^{ij}(t_n^r) + \frac{\omega_i\omega_j J_{\kappa}^{ij}}{\Delta t^{(n)}} \sum_{m=0}^{q} (d_n\omega_r\omega_m + 2|Q_{rm}^q|L_u)|\varphi_{\kappa}^{ijm} - \varphi_{\kappa}^{ijr}|, $$

where the $a_{\kappa}^{ij}$ and $b_{\kappa}^{ij}$ are defined in (3.8). As in (3.8), we define $\varphi_{\kappa}^{ijr} = \varphi_{\kappa}^{ij}(t_n^r) \ge 0$ the average of $\varphi(\mathbf{x}, t_n^r) = e^{-\gamma(|\mathbf{x}|+t_n^r-t_n^0)}$ over the subcell $\kappa_{ij} \subset \kappa$. Using (3.9) and (2.11) together with $\omega_i\omega_j J_{\kappa}^{ij} \le |\kappa| \le h^d$ we obtain

$$ a_{\kappa}^{ijr} \ge \frac{\nu\omega_0^{d+2}\beta h_{\Omega}^d}{2d\Delta t^{(n)}} \inf_{B_{\kappa}(h_{\Omega}) \times [t_n^0, t_n^1]} \varphi(\mathbf{x}, t), \quad b_{\kappa}^{ijr} \le h_{\Omega}^d \left( \frac{8d_{\kappa} + L_1 + L_2}{d\beta} + \frac{8d_n}{\Delta t^{(n)}} \right) \sup_{B_{\kappa}(2h_{\Omega}) \times [t_n^0, t_n^1]} |\nabla_{\mathbf{x},t}\varphi|, $$

and using $|\nabla_{\mathbf{x},t}\varphi| = |(\nabla\varphi^{\top}, \partial_t\varphi)^{\top}| = \sqrt{2}\gamma\varphi$, any choice of $\gamma > 0$ small enough such that $\gamma e^{\gamma(3h+2\Delta t^{(n)})} (\frac{8d_{\kappa}+L_1+L_2}{d\beta} + \frac{8d_n}{\Delta t^{(n)}}) < \frac{\nu\omega_0^{d+2}\beta}{2\sqrt{2}\Delta t^{(n)}}$ imposes $a_{\kappa}^{ijr} > b_{\kappa}^{ijr}$, so $U_{\kappa}^{ij}(t_n^r) = W_{\kappa}^{ij}(t_n^r)$ for all the DOFs.

Uniqueness at $\nu = 0$. Existence and stability hold whatever $\nu \ge 0$, but uniqueness requires $\nu > 0$. We now justify that uniqueness can be extended to $\nu = 0$. The scheme is a Lipschitz continuous function of $u_h$ and $\nu$. Applying the implicit function theorem for non-differentiable functions [36], for every $\nu > 0$ there exists a unique $u_h = u_{h,\nu}$ and neighborhoods of $u_{h,\nu}$ and $\nu$ where $\nu \mapsto u_{h,\nu}$ is one-to-one and continuous. It is even a Cauchy-continuous function over $(0, \infty)$, since for all Cauchy sequence $(\nu_i)_i$ in $(0, \infty)$, $(u_{h,\nu_i})_i$ is a Cauchy sequence in $\ell_{\infty}(\mathbb{N})$. The latter being a Banach space, $u_{h,\nu}$ can be extended to $\nu = 0$ by Cauchy-continuity. Since $u_h$ exists and is bounded at $\nu = 0$, we obtain uniqueness in this limit. This concludes the proof. $\square$

**Theorem 4.3.** Under the assumptions of Theorem 4.2, the space-time DGSEM (4.1) with the EC fluxes (2.13) and (4.3), defined for a given entropy pair $(\vartheta(u), \mathbf{g}(u))$ with $\vartheta''(\cdot) > 0$, and the graph viscosities in (3.2) and (4.2) satisfies the following inequality for every entropy pair $(\eta(u), \mathbf{q}(u))$ with $\eta''(\cdot) \ge 0$:

$$ \langle \eta(u_h(t_n^q)) \rangle_{\kappa} - \langle \eta(u_h(t_{n-1}^q)) \rangle_{\kappa} + \frac{\Delta t^{(n)}}{|\kappa|} \sum_{r=0}^{q} \frac{\omega_r}{2} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \omega_k J_e^k Q(u_h^-(\mathbf{x}_e^k, t_n^r), u_h^+(\mathbf{x}_e^k, t_n^r), \mathbf{n}_e^k) \le 0, \quad (4.8) $$

with the numerical entropy flux $Q(\cdot, \cdot, \cdot)$ defined in (2.18).

*Proof.* The proof follows the lines of the proof of Theorem 3.3 by applying the convex entropy $\eta(\cdot)$ to the limit $\mathbf{m} \to \infty$ of the convex combination (4.7) with $\nu = 0$ and further using (2.18) and (2.14b), so we omit the details. Compared to Theorem 3.3, the only additional terms to the RHS come from the time contributions. Summing them over $0 \le r \le q$ and using the definition of $\gamma_n^{rm} = \omega_r(d_n\omega_m - 4D_{rm}^q\beta_{rm})$, they read

$$ -\eta_n^0 + \eta_{n-1}^q - \sum_{r,m=0}^{q} \omega_r (d_n \omega_m - 4D_{rm}^q \beta_{rm}) (\eta_n^r - \eta_n^m) = -\eta_n^0 + \eta_{n-1}^q + 4 \sum_{r,m=0}^{q} (Q_{rm}^q \beta_{rm} - Q_{mr}^q \beta_{mr}) \eta_n^r \\ \stackrel{(2.7)}{=} \stackrel{(2.15a)}{-\eta_n^0 + \eta_{n-1}^q - 4 \sum_{r,m=0}^{q} (Q_{rm}^q + Q_{mr}^q) \beta_{mr} \eta_n^r} \stackrel{(2.6)}{=} -\eta_n^q + \eta_{n-1}^q, $$

where we have switched indices $r \leftrightarrow m$ and then used $\beta_{rm} = \frac{1}{2} - \beta_{mr}$ and $\beta_{00} = \beta_{qq} = \frac{1}{4}$ from the definition (2.15a) of $\beta_{rm} = \beta(U_n^r, U_n^m)$. The proof is complete. $\square$
**Remark 4.4** (Kružkov's entropies). The DGSEM in Lemmas 3.1 and 4.1 without graph viscosity is not ES for the Kružkov's entropies [42], $\vartheta(u) = |u - K|$, $g(u) = \operatorname{sgn}(u - K)(f(u) - f(K))$, $K \in \mathbb{R}$, since it requires the mapping $v(u) = \vartheta'(u)$ to be one-to-one, but $v'(u) = \delta(u - K)$. Applying (2.13) and (4.3b), we obtain $h_{ec}(u^{-}, u^{+}, \mathbf{n}) = f(K) \cdot \mathbf{n}$ and $U_{ec}(u^{-}, u^{+}) = K$ which are not consistent. Using the graph viscosity and EC fluxes defined for one strictly convex function, e.g., $\vartheta(u) = \frac{u^2}{2}$, however applies to them as it only requires $\eta''(\cdot) \ge 0$ and, in particular, Theorems 3.3 and 4.3 hold for $\eta(u) = |u - K|$. Kružkov's entropies will be used in the experiments to limit the high-order solution (see Sect. 5 and 6).

**Remark 4.5** (Bounds in limiters). Adding graph viscosity to the schemes (3.1) and (4.1) annihilates they high-order accuracy. These low-order schemes may however constitute building-blocks of limiters such as the flux-corrected transport limiter [3, 52, 73], or the convex limiter [27, 55]. The properties of the present schemes indeed allow to apply several ways to limit the solution to impose different bounds to the high-order solution such as the maximum principle in Theorems 3.2 and 4.2, or entropy inequalities based on any convex entropy in Theorems 3.3 and 4.3. In Section 5, we introduce a two-step limiter based on the FCT limiter and a linear scaling limiter.

# 5. GRAPH VISCOSITY AS A BUILDING BLOCK IN A TWO-STEP LIMITER

We illustrate the use of low-order schemes (3.1) and (4.1) as building blocks in a two-step limiter to make the high-order schemes (i.e., the same schemes (3.1) and (4.1), but without graph viscosities, $d_{\kappa} = d_{n} = 0$) MPP and ES for all the entropies imposed by the user. As advocated in Section 4.3 from [52], the two-step limiter applies first the FCT limiter [3, 28, 44, 45, 52, 73] to the cell averages by imposing convex bounds satisfied by the low-order solution. The bounds are then imposed to all DOFs within the elements through the use of the linear scaling limiter [74, 75]. This two-step limiting process avoids undesirable effects on the approximate solution when the FCT limiting is applied directly to all DOFs, such as accuracy deterioration for smooth solutions, or a frequent switching back and forth between the limited and unlimited schemes ([52], Appendix A). Below we describe the limiter for the space-time DGSEM in Section 4, the application of the limiter to the time implicit DGSEM in Section 3 following the same lines.

## 5.1. FCT limiter for the cell averages

By $u_{LO}$ and $u_{HO}$ we denote the solutions to the low-order scheme, (4.1) with $d_{\kappa}$ and $d_{n}$ given by (4.5), and high-order scheme, (4.1) with $d_{\kappa} = d_{n} = 0$, respectively. Both are solutions to the cell-averaged scheme (4.4). Subtracting the cell-averaged scheme for the low-order solution from the one for the high-order solution gives

$$ \langle u_{HO}(t_n^q) \rangle_{\kappa} - \langle u_{LO}(t_n^q) \rangle_{\kappa} = \frac{\Delta t^{(n)}}{|\kappa|} \sum_{r=0}^{q} \frac{\omega_r}{2} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \omega_k J_e^k \left( h(u_{LO}^{-}(\mathbf{x}_e^k, t_n^r), u_{LO}^{+}(\mathbf{x}_e^k, t_n^r), \mathbf{n}_e^k) \right. \\ \left. - h(u_{HO}^{-}(\mathbf{x}_e^k, t_n^r), u_{HO}^{+}(\mathbf{x}_e^k, t_n^r), \mathbf{n}_e^k) \right) =: \frac{\Delta t^{(n)}}{|\kappa|} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \Delta h_{\kappa, e}^{k, n}, \quad (5.1) $$

where $-\Delta h_{\kappa,e}^{k,n}$ correspond to the antidiffusive fluxes [45, 73] integrated over the time step. Note that by conservation of the numerical flux we have $\Delta h_{\kappa,e}^{n} = -\Delta h_{\kappa,e}^{k,n}$, where $\kappa_e^+$ is the neighboring element sharing the face $e \in \partial \kappa$ (see Fig. 1). The limiter is then defined as follows [28, 45, 73]: we first introduce the limiting factors $0 \le l_{\kappa}^{\pm} \le 1$ defined by

$$ l_{\kappa}^{-} = \min \left( 1, \frac{Q_{\kappa}^{-}}{P_{\kappa}^{-}} \right), \quad P_{\kappa}^{-} = \frac{\Delta t^{(n)}}{|\kappa|} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \min(\Delta h_{\kappa,e}^{k,n}, 0) \le 0, \quad Q_{\kappa}^{-} = m_{\kappa} - \langle u_{LO}(t_n^q) \rangle_{\kappa} \le 0, \quad (5.2a) $$

$$ l_{\kappa}^{+} = \min \left( 1, \frac{Q_{\kappa}^{+}}{P_{\kappa}^{+}} \right), \quad P_{\kappa}^{+} = \frac{\Delta t^{(n)}}{|\kappa|} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \max(\Delta h_{\kappa,e}^{k,n}, 0) \ge 0, \quad Q_{\kappa}^{+} = M_{\kappa} - \langle u_{LO}(t_n^q) \rangle_{\kappa} \ge 0, \quad (5.2b) $$
where $m_κ$ and $M_κ$ are local lower and upper bounds we want to impose to the limited solution $\langle \tilde{u}_h(t_n^q) \rangle_κ$.
The solution in cell $κ$ is then modified through flux limiting according to

$$ \frac{\omega_i \omega_j J_{\kappa}^{ij}}{|\kappa|} (\tilde{U}_{\kappa}^{ij}(t_n^q) - U_{\kappa, \text{HO}}^{ij}(t_n^q)) = \frac{\Delta t^{(n)}}{|\kappa|} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} (l_{\kappa, e}^{k,n} - 1) \phi_{\kappa}^{ij}(\mathbf{x}_e^k) \Delta h_{\kappa, e}^{k,n} \quad \forall 0 \le i, j \le p, \quad (5.3a) $$

so that the cell average of the limited solution reads

$$ \langle \tilde{u}_h(t_n^q) \rangle_\kappa = \langle u_{LO}(t_n^q) \rangle_\kappa + \frac{\Delta t^{(n)}}{|\kappa|} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} l_{\kappa, e}^{k,n} \Delta h_{\kappa, e}^{k,n}, \quad l_{\kappa, e}^{k,n} = \begin{cases} \min(l_\kappa^-, l_\kappa^+) & \text{if } \Delta h_{\kappa, e}^{k,n} < 0 \\ \min(l_\kappa^-, l_\kappa^+) & \text{else.} \end{cases} \quad (5.3b) $$

## 5.2. Linear scaling limiter for the DOFs

We then apply a linear scaling limiter [74, 75] that consists in scaling the solution around its cell average to enforce the bounds at all DOFs:

$$ U_{\kappa}^{ij}(t_n^q) = \theta_{\kappa} \tilde{U}_{\kappa}^{ij}(t_n^q) + (1 - \theta_{\kappa}) \langle \tilde{u}_h(t_n^q) \rangle_{\kappa}, \quad 0 \le i, j \le p, \quad \kappa \in \Omega_h, \quad (5.4a) $$

where $0 \le \theta_κ \le 1$ is defined by

$$ \theta_{\kappa} = \min \left( \frac{M_{\kappa} - \langle \tilde{u}_h(t_n^q) \rangle_{\kappa}}{\max_{ij} \tilde{U}_{\kappa}^{ij}(t_n^q) - \langle \tilde{u}_h(t_n^q) \rangle_{\kappa}}, \frac{m_{\kappa} - \langle \tilde{u}_h(t_n^q) \rangle_{\kappa}}{\min_{ij} \tilde{U}_{\kappa}^{ij}(t_n^q) - \langle \tilde{u}_h(t_n^q) \rangle_{\kappa}}, 1 \right). \quad (5.4b) $$

## 5.3. Definition of bounds

The local bounds in (5.2) and (5.4b) are designed to impose the maximum principle and entropy stability at shocks:

$$ m_{\kappa} = \max \left( m, s_{\kappa} \max_{\eta} m_{\kappa}^{\eta} + (1 - s_{\kappa})m \right), \quad M_{\kappa} = \min \left( M, s_{\kappa} \min_{\eta} M_{\kappa}^{\eta} + (1 - s_{\kappa})M \right), \quad (5.5) $$

where $m$ and $M$ are the maximum principle bounds in (1.3), while $m_κ^η$ and $M_κ^η$ are, respectively, lower and upper roots matching the RHS of (4.8):

$$ \eta(m_{\kappa}^{\eta}) = \eta(M_{\kappa}^{\eta}) = \langle \eta(u_h(t_{n-1}^{q})) \rangle_{\kappa} - \frac{\Delta t^{(n)}}{|\kappa|} \sum_{r=0}^{q} \frac{\omega_r}{2} \sum_{e \in \partial \kappa} \sum_{k=0}^{p} \omega_k J_e^k Q(u_{LO}^{-}(\mathbf{x}_e^k, t_n^r), u_{LO}^{+}(\mathbf{x}_e^k, t_n^r), \mathbf{n}_e^k), $$

for any admissible convex entropy $\eta(\cdot)$ required by the user. We use the Kružkov's entropy $\eta(u) = |u|$ (see Rem. 4.4) in the numerical experiments. It is always possible to apply these convex bounds to $\langle u_{HO} \rangle_κ$, even if the high-order scheme is not ES for $\eta(\cdot)$, since $\max(m, m_κ^η) \le \langle u_{LO}(t_n^q) \rangle_κ \le \min(M, M_κ^η)$ by Theorems 4.2 and 4.3 and using Jensen's inequality, $\eta(\langle u_{LO}(t_n^q) \rangle_κ) \le \langle \eta(u_{LO}(t_n^q)) \rangle_κ$. As a consequence, the limited scheme satisfies the maximum principle and entropy inequalities for all entropies imposed by the user in (5.5) in addition to the entropy $\vartheta$ used to define the EC fluxes (2.13) and (4.3).

By $s_κ$ in (5.5) we denote a shock indicator to impose entropy bounds in shock regions only. Strictly imposing entropy stability may indeed result in overlimiting, see, e.g., Section 4.7 from [27] and Section 3.3 from [37], so we now design the indicator so that $s_κ = 1$ in shock regions, and $s_κ \ll 1$ else. We use the indicator from Persson and Peraire [56] evaluated from $u_{HO}(t_n^q)$ expressed in the basis of Legendre polynomials by first computing $S_κ = \log (\|u_{HO} - \Pi_h^{p-1}u_{HO}\|_{L^2(\kappa)}^2 / \|u_{HO}\|_{L^2(\kappa)}^2)$ with $\Pi_h^{p-1}$ the projector onto $V_h^{p-1}$. Using $S_0 = -4 \log(2p)$, we set: $s_κ = 0$ if $S_κ < S_0 - \frac{1}{10}$; $s_κ = 1$ if $S_κ > S_0 + \frac{1}{10}$; $s_κ = \frac{1}{2} + \frac{1}{2} \sin(5\pi(S_κ - S_0))$ else.
## 5.4. Overview of the two-step limiter

The limiter is applied after each time step and Algorithm 1 gives details on the algorithm flowchart. Let us first consider efficiency aspects of the above limiting procedure. The FCT limiter requires the computation of the low-order and high-order solutions, which induces an additional cost. Once the high-order solution has been computed, we can check if the FCT limiter is required, i.e., if some cell averages $\langle u_{HO}(t_n^q) \rangle_\kappa$ violate the bounds in (5.5) (either $m$ and $M$, or $m_\kappa$ and $M_\kappa$ evaluated from $u_{HO}(t_n^q)$). If not, we apply only the linear scaling limiter (5.4) to the DOFs. In the case of smooth solutions, violating cell averages are expected to occur far less frequently than violating DOFs. Indeed, assume for instance that $u_{HO}$ locally violates the lower bound in some cell $\kappa: \langle u_{HO}(t_n^q) \rangle_\kappa \le m_\kappa \le \langle u(\cdot, t_n^q) \rangle_\kappa$ with $u$ the exact solution. Further assuming that the high-order is uniformly accurate, $\|u_{HO} - u\|_{L^\infty(\Omega_h)} = O(h_\Omega^{\min(p,q)+1})$, we deduce $\|u(\cdot, t_n^q) - m_\kappa\|_{L^\infty(\kappa)} = O(h_\Omega^{\min(p,q)+1})$ which corresponds to a flat region in $\kappa$, e.g., at a local extremum, but is rare on coarse meshes as usually employed with high-order schemes. Moreover, the application of the FCT limiter in (5.3a) only involves loops on faces, while the DOFs limiter requires one loop over elements. Both are easy to implement. When the application of the FCT limiter is required, then computing $u_{LO}$ dominates the cost of the two-step limiter and, hence, roughly doubles the global cost of the computation. As a consequence, the cost of the two-step limiter highly depends on the solution to be computed.

<caption>**Algorithm 1.** Algorithm flowchart of the two-step limiter (5.3) and (5.4) applied at each time step $n \ge 0$.</caption>

**Require:** $u_h(t_n^q)$

1: compute $u_{HO}(t_n^q)$ from $u_h(t_{n-1}^q)$ using (3.1) or (4.1) with $d_\kappa = d_n = 0$

2: boundOK ← True

3: **for** $\kappa \in \Omega_h$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\langle u_{HO}(t_n^q) \rangle_\kappa < m$ **or** $\langle u_{HO}(t_n^q) \rangle_\kappa > M$ **then**

5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;boundOK ← False

6: &nbsp;&nbsp;&nbsp;&nbsp;**end if**

7: **end for**

8: **if** boundOK=False **then**

9: &nbsp;&nbsp;&nbsp;&nbsp;compute $u_{LO}(t_n^q)$ using (3.1) or (4.1) with $d_\kappa$ and $d_n$ given by (4.5)

10: compute the local bounds (5.5)

11: **for** $\kappa \in \Omega_h$ **do**

12: &nbsp;&nbsp;&nbsp;&nbsp;apply the cell-average limiter (5.3) to get $\tilde{u}_h(t_n^q)$

13: **end for**

14: **else**

15: &nbsp;&nbsp;&nbsp;&nbsp;$\tilde{u}_h(t_n^q) \leftarrow u_{HO}(t_n^q)$

16: **end if**

17: apply the DOFs limiter (5.4) to get $u_h(t_n^q)$

The two-step limiter (5.3) and (5.4) is conservative since the DOFs limiter (5.4) does not change the cell average of the solution. Then, by symmetry of the limiting coefficients in (5.3), $l_{\kappa,e}^{k,n} = l_{\kappa^+,e}^{k,n}$, and skew-symmetry of the antidiffusive fluxes, $\Delta h_{\kappa^+,e}^n = -\Delta h_{\kappa,e}^{k,n}$, we have

$$ \sum_{\kappa \in \Omega_h} |\kappa| \langle \tilde{u}_h(t_n^q) \rangle_\kappa = \sum_{\kappa \in \Omega_h} |\kappa| \langle u_{HO}(t_n^q) \rangle_\kappa = \sum_{\kappa \in \Omega_h} |\kappa| \langle u_{LO}(t_n^q) \rangle_\kappa = \sum_{\kappa \in \Omega_h} |\kappa| \langle u_h(t_{n-1}^q) \rangle_\kappa. $$

The two-step limiter also imposes the bounds (5.5) to all DOFs in the limited solution: from (5.3b) we deduce

$$ \langle \tilde{u}_h(t_n^q) \rangle_\kappa - \langle u_{LO}(t_n^q) \rangle_\kappa \ge \sum_{e,k} l_{\kappa,e}^{k,n} \min(\Delta h_{\kappa,e}^{k,n}, 0) \ge l_\kappa^- P_\kappa^- \ge Q_\kappa^- = m_\kappa - \langle u_{LO}(t_n^q) \rangle_\kappa, $$

$$ \langle \tilde{u}_h(t_n^q) \rangle_\kappa - \langle u_{LO}(t_n^q) \rangle_\kappa \le \sum_{e,k} l_{\kappa,e}^{k,n} \max(\Delta h_{\kappa,e}^{k,n}, 0) \le l_\kappa^+ P_\kappa^+ \le Q_\kappa^+ = M_\kappa - \langle u_{LO}(t_n^q) \rangle_\kappa, $$
<table><caption>TABLE 1. Definitions of the one-dimensional problems and numerical parameters.</caption><thead><tr><th>#</th><th>f(u)</th><th>u<sub>0</sub>(x)</th><th>u(0, t)</th><th>u(1, t)</th><th>t</th><th>p</th><th>q</th><th>h<sub>Ω</sub></th><th>CFL</th><th>L<sub>f</sub></th><th>L<sub>u</sub></th></tr></thead><tbody><tr><td>0</td><td><sup>u<sup>2</sup></sup>&frasl;<sub>2</sub></td><td>sin(2πx)</td><td>u(1, t)</td><td>—</td><td><sup>1</sup>&frasl;<sub>8</sub></td><td>1–4</td><td>1–4</td><td><sup>1</sup>&frasl;<sub>40</sub></td><td>1</td><td>1</td><td>1</td></tr><tr><td>1</td><td>u</td><td>Example 1 of [34]</td><td>u(1, t)</td><td>—</td><td>1</td><td>3</td><td>3</td><td><sup>1</sup>&frasl;<sub>40</sub></td><td>1</td><td><sup>1</sup>&frasl;<sub>100</sub></td><td>1</td></tr><tr><td>2</td><td><sup>u<sup>2</sup></sup>&frasl;<sub>2</sub></td><td>sin(2πx)</td><td>u(1, t)</td><td>—</td><td><sup>2</sup>&frasl;<sub>5</sub></td><td>3</td><td>3</td><td><sup>1</sup>&frasl;<sub>40</sub></td><td>1</td><td>1</td><td>1</td></tr><tr><td>3</td><td><sup>u<sup>2</sup></sup>&frasl;<sub>2</sub></td><td>1 – 2x</td><td>1</td><td>–1</td><td>—</td><td>3</td><td>0</td><td><sup>1</sup>&frasl;<sub>40</sub></td><td>10<sup>3</sup></td><td>1</td><td>1</td></tr><tr><td>4</td><td><sup>u<sup>2</sup></sup>&frasl;<sub>2</sub></td><td>1 + sin(2πx)</td><td>u(1, t)</td><td>—</td><td><sup>3</sup>&frasl;<sub>4π</sub></td><td>3</td><td>3</td><td><sup>1</sup>&frasl;<sub>40</sub></td><td>1</td><td>2</td><td>1</td></tr><tr><td>5</td><td><sup>u<sup>2</sup></sup>&frasl;<sub>u<sup>2</sup> + (<sup>1–u</sup>&frasl;<sub>2</sub>)<sup>2</sup></sub></td><td>1<sub>x<0.5</sub></td><td>1</td><td>0</td><td><sup>1</sup>&frasl;<sub>5</sub></td><td>3</td><td>3</td><td><sup>1</sup>&frasl;<sub>40</sub></td><td>1</td><td>2.10</td><td>1</td></tr><tr><td>6</td><td><sup>u<sup>2</sup></sup>&frasl;<sub>u<sup>2</sup> + (<sup>1–u</sup>&frasl;<sub>4</sub>)<sup>2</sup></sub></td><td>3(1<sub>x&gt;0.5</sub> – 1<sub>x&lt;0.5</sub>)</td><td>–3</td><td>3</td><td>1</td><td>3</td><td>3</td><td><sup>1</sup>&frasl;<sub>40</sub></td><td>1</td><td>2.35</td><td>1</td></tr></tbody></table>

hence $m_κ \le \langle \tilde{u}_h(t_n^q) \rangle_κ \le M_κ$, so are the whole DOFs after (5.4).

## 6. NUMERICAL EXPERIMENTS

In this section we present numerical experiments on problems in one and two space dimensions. The objective is here to illustrate the properties of the DGSEM schemes (3.1) and (4.1), and of the limiter from Section 5. We consider problems with either a convex, or a non convex flux.

We define the EC space and time fluxes (2.13) and (4.3) with the square entropy $\vartheta(u) = \frac{u^2}{2}$, so $L_u = 1$ in (3.4) and (4.5) (see Rems. 2.3 and 3.5), and use the Godunov solver [24] at interfaces. We use a high-order scheme with $1 \le p = q \le 4$ in (4.1), except one problem that is steady-state and is solved with the backward Euler DGSEM (3.1) and $p = 3$. The values of $\max_{k \ne l} \frac{D_{kl}}{\omega_l}$ in (3.4) and (4.5) are 1, 3, 9.71, 24.8 for $p$ or $q$ values from 1 to 4 ([52], Tab. 2). The values of $L_f$ in (1.4) are given in Tables 1 and 3 and depend on the physical problem (1.1), so the graph viscosities in (3.4) and (4.5) are easily evaluated a priori once the problem is defined.

The time step is computed from $\frac{\Delta t^{(n)}}{h_{\Omega}} L_f \le \text{CFL}$, where the CFL values are given in Tables 1 and 3. The nonlinear algebraic systems (3.1) and (4.1) are solved by using a quasi-Newton iterative algorithm. Let $\mathcal{R}(\mathbf{U}^{(n)}) = 0$ be the full system (3.1) or (4.1) of size $(q+1)(p+1)^d N$ with $N$ the number of mesh elements and $\mathbf{U}^{(n)}$ the global vector of DOFs at times $t_n^{0 \le r \le q}$. we use subiterations of the form $(\mathbf{M} + \partial_{\mathbf{U}} \mathcal{R}(\mathbf{U}^{(n,m)}))(\mathbf{U}^{(n,m+1)} - \mathbf{U}^{(n,m)}) = -\mathcal{R}(\mathbf{U}^{(n,m)})$ with $\mathbf{M}$ the global mass matrix, $\mathbf{U}^{(n,0)}$ defined from $u_h(\cdot, t^{(n)}$); we linearize the residuals exactly at each subiteration for evaluating the Jacobian. We stop the subiterations when $\|\mathcal{R}(\mathbf{U}^{(n,m)})\|_{\ell^2} < 10^{-10}$.

We will evaluate the effect of graph viscosity by comparing three different computations with schemes (4.1) and (3.1):

**HO scheme:** the graph viscosity is removed by setting $d_κ = d_n = 0$;

**LO scheme:** the graph viscosity coefficients are evaluated from their theoretical values (3.4) and (4.5);

**limited scheme:** we first evaluate the high-order solution with the high-order scheme, then applies the two-step limiter in Section 5.

### 6.1. One-dimensional problems

We consider problems of the form

$$ \partial_t u + \partial_x f(u) = 0 \text{ in } [0, 1] \times (0, \infty), \quad u(\cdot, 0) = u_0(\cdot) \text{ in } [0, 1], \qquad (6.1) $$

with appropriate boundary conditions at $x = 0$ and $x = 1$. The different problems are defined in Table 1.
<table><caption>TABLE 2. Burgers' equation (problem 0): $L^2$ and $L^\infty$ norms of the error obtained at $t = \frac{1}{8}$ with the HO scheme, the limited scheme, and the limited scheme while systematically imposing boundOK=False in Algorithm 1.</caption><thead><tr><th rowspan="2">v<br>||<sub>Ω</sub>||</th><th rowspan="2">1/h<sub>Ω</sub></th><th colspan="3">HO scheme</th><th colspan="3">Limited scheme</th><th colspan="3">boundOK=False</th></tr><tr><th>||e<sub>h</sub>||<sub>L<sup>2</sup>(Ω)</sub></th><th>O<sub>2</sub></th><th>||e<sub>h</sub>||<sub>L<sup>∞</sup>(Ω)</sub></th><th>O<sub>∞</sub></th><th>||e<sub>h</sub>||<sub>L<sup>2</sup>(Ω)</sub></th><th>O<sub>2</sub></th><th>||e<sub>h</sub>||<sub>L<sup>∞</sup>(Ω)</sub></th><th>O<sub>∞</sub></th></tr></thead><tbody><tr><td rowspan="6">1</td><td>20</td><td>5.564e-02</td><td>—</td><td>2.257e-01</td><td>—</td><td>5.405e-02</td><td>—</td><td>2.259e-01</td><td>—</td><td>6.815e-02</td><td>—</td><td>2.747e-01</td><td>—</td></tr><tr><td>40</td><td>2.503e-02</td><td>1.15</td><td>1.394e-01</td><td>0.69</td><td>2.497e-02</td><td>1.15</td><td>1.394e-01</td><td>0.69</td><td>3.107e-02</td><td>1.13</td><td>1.605e-01</td><td>0.77</td></tr><tr><td>80</td><td>1.061e-02</td><td>1.23</td><td>7.943e-02</td><td>0.81</td><td>1.061e-02</td><td>1.23</td><td>7.943e-02</td><td>0.81</td><td>1.281e-02</td><td>1.28</td><td>8.842e-02</td><td>0.86</td></tr><tr><td>160</td><td>4.190e-03</td><td>1.34</td><td>4.258e-02</td><td>0.90</td><td>4.191e-03</td><td>1.34</td><td>4.258e-02</td><td>0.90</td><td>4.862e-03</td><td>1.40</td><td>4.650e-02</td><td>0.93</td></tr><tr><td>320</td><td>1.582e-03</td><td>1.41</td><td>2.183e-02</td><td>0.96</td><td>1.582e-03</td><td>1.41</td><td>2.183e-02</td><td>0.96</td><td>1.797e-03</td><td>1.44</td><td>2.370e-02</td><td>0.97</td></tr><tr><td>20</td><td>1.100e-02</td><td>—</td><td>5.936e-02</td><td>—</td><td>1.100e-02</td><td>—</td><td>5.936e-02</td><td>—</td><td>1.116e-02</td><td>—</td><td>5.936e-02</td><td>—</td></tr><tr><td rowspan="6">2</td><td>40</td><td>2.659e-03</td><td>2.05</td><td>1.777e-02</td><td>1.74</td><td>2.659e-03</td><td>2.05</td><td>1.777e-02</td><td>1.74</td><td>2.982e-03</td><td>1.90</td><td>1.777e-02</td><td>1.74</td></tr><tr><td>80</td><td>4.831e-04</td><td>2.46</td><td>5.333e-03</td><td>1.74</td><td>4.831e-04</td><td>2.46</td><td>5.333e-03</td><td>1.74</td><td>5.698e-04</td><td>2.39</td><td>5.333e-03</td><td>1.74</td></tr><tr><td>160</td><td>8.481e-05</td><td>2.51</td><td>1.019e-03</td><td>2.39</td><td>8.481e-05</td><td>2.51</td><td>1.019e-03</td><td>2.39</td><td>1.019e-04</td><td>2.48</td><td>1.019e-03</td><td>2.39</td></tr><tr><td>320</td><td>1.384e-05</td><td>2.61</td><td>2.137e-04</td><td>2.25</td><td>1.384e-05</td><td>2.61</td><td>2.137e-04</td><td>2.25</td><td>1.706e-05</td><td>2.58</td><td>2.137e-04</td><td>2.25</td></tr><tr><td>20</td><td>1.691e-03</td><td>—</td><td>1.016e-02</td><td>—</td><td>1.691e-03</td><td>—</td><td>1.016e-02</td><td>—</td><td>1.691e-03</td><td>—</td><td>1.016e-02</td><td>—</td></tr><tr><td>40</td><td>1.624e-04</td><td>3.38</td><td>1.798e-03</td><td>2.50</td><td>1.624e-04</td><td>3.38</td><td>1.798e-03</td><td>2.50</td><td>1.624e-04</td><td>3.38</td><td>1.798e-03</td><td>2.50</td></tr><tr><td rowspan="6">3</td><td>80</td><td>4.421e-05</td><td>1.88</td><td>6.024e-04</td><td>1.58</td><td>4.421e-05</td><td>1.88</td><td>6.024e-04</td><td>1.58</td><td>4.421e-05</td><td>1.88</td><td>6.024e-04</td><td>1.58</td></tr><tr><td>160</td><td>6.181e-06</td><td>2.84</td><td>1.005e-04</td><td>2.58</td><td>6.181e-06</td><td>2.84</td><td>1.005e-04</td><td>2.58</td><td>6.181e-06</td><td>2.84</td><td>1.005e-04</td><td>2.58</td></tr><tr><td>320</td><td>6.504e-07</td><td>3.25</td><td>1.362e-05</td><td>2.88</td><td>6.504e-07</td><td>3.25</td><td>1.362e-05</td><td>2.88</td><td>6.504e-07</td><td>3.25</td><td>1.362e-05</td><td>2.88</td></tr><tr><td>20</td><td>3.168e-04</td><td>—</td><td>3.678e-03</td><td>—</td><td>3.167e-04</td><td>—</td><td>3.678e-03</td><td>—</td><td>3.168e-04</td><td>—</td><td>3.678e-03</td><td>—</td></tr><tr><td>40</td><td>8.157e-05</td><td>1.96</td><td>9.888e-04</td><td>1.90</td><td>8.157e-05</td><td>1.96</td><td>9.888e-04</td><td>1.90</td><td>8.157e-05</td><td>1.96</td><td>9.888e-04</td><td>1.90</td></tr><tr><td>80</td><td>6.845e-06</td><td>3.58</td><td>8.246e-05</td><td>3.58</td><td>6.845e-06</td><td>3.58</td><td>8.246e-05</td><td>3.58</td><td>6.845e-06</td><td>3.58</td><td>8.246e-05</td><td>3.58</td></tr><tr><td rowspan="3">4</td><td>160</td><td>2.742e-07</td><td>4.64</td><td>5.141e-06</td><td>4.00</td><td>2.742e-07</td><td>4.64</td><td>5.141e-06</td><td>4.00</td><td>2.742e-07</td><td>4.64</td><td>5.141e-06</td><td>4.00</td></tr><tr><td>320</td><td>1.188e-08</td><td>4.53</td><td>2.850e-07</td><td>4.17</td><td>1.188e-08</td><td>4.53</td><td>2.850e-07</td><td>4.17</td><td>1.188e-08</td><td>4.53</td><td>2.850e-07</td><td>4.17</td></tr></tbody></table>

Table 2 displays error levels for problem 0 involving a smooth solution where we compare results obtained with the high-order scheme to results with the limited scheme. We also provide results when systematically applying the FCT limiter (5.3) in the two-step limiter in Algorithm 1 to analyse its effect on the accuracy since the FCT limiter is not applied in practice as the bounds (5.5) are not violated by the cell-averaged high-order solution. We observe that the limited scheme successfully preserves local extrema even if the FCT limiter is applied as expected from Section 5.4. The limited scheme also achieves high-order accuracy on the full solution and displays error levels in close agreement with the high-order solution.

We now consider problems involving discontinuous solutions. Results are displayed in Figure 2 for the linear transport equation, in Figure 3 for the inviscid Burgers' equation and in Figure 4 for the Buckley-Leverett equation. The DGSEM schemes without graph viscosities are ES for the square entropy $\vartheta(u) = \frac{u^2}{2}$ (see Lems. 3.1 and 4.1), so we expect them to capture the entropy weak solution when the flux is convex [54] as highlighted in Figure 3. This is however no longer the case for the Buckley-Leverett equation in Figure 4, where we observe nonphysical weak solutions. Adding graph viscosity allows to capture the entropy weak solution, but results in large dissipation smearing the discontinuities, in particular for problems 4 and 5. This dissipative character may be illustrated from a comparison with the scaling limiter [75]. Consider, e.g., scheme (3.1) that we rewrite in one space dimension as

$$ \left( 1 + \frac{\Delta t^{(n)}}{J_{\kappa}^{i}} d_{\kappa} \right) U_{\kappa}^{i, n+1} = \frac{\Delta t^{(n)}}{J_{\kappa}^{i}} d_{\kappa} \langle u_{h}^{(n+1)} \rangle_{\kappa} + \hat{U}_{\kappa}^{i, n+1} \quad \forall 0 \le i \le p, $$

with $J_{\kappa}^{i} = \frac{1}{2} \text{diam} \kappa$, $\langle u_{h}^{(n+1)} \rangle_{\kappa} = \frac{1}{2} \sum_{l=0}^{p} \omega_{l} U_{\kappa}^{l, n+1}$ the cell-averaged solution, and $\hat{U}_{\kappa}^{i, n+1} = U_{\kappa}^{i, n} - \frac{\Delta t^{(n)}}{J_{\kappa}^{i} \omega_{i}} R_{\kappa}^{i}(u_{h}^{(n+1)})$ which may be viewed as an unlimited solution. We therefore have $U_{\kappa}^{i, n+1} = (1 - \theta_{\kappa}) \langle u_{h}^{(n+1)} \rangle_{\kappa} + \theta_{\kappa} \hat{U}_{\kappa}^{i, n+1}$ with $\theta_{\kappa} = (1 + \frac{\Delta t^{(n)}}{J_{\kappa}^{i}} d_{\kappa})^{-1} \in (0, 1)$, similar to the linear scaling limiter (5.4a) of $u_{h}^{(n+1)}$ around its cell average [75]. The limit for large viscosity, $\lim_{d_{\kappa} \to \infty} \theta_{\kappa} = 0$, indeed corresponds to imposing $U_{\kappa}^{i, n+1} = \langle u_{h}^{(n+1)} \rangle_{\kappa}$ as may be observed in Figures 2–4. Note that the graph viscosity may be relaxed in the one-dimensional case for
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Linear transport with nonsmooth solution: solutions to problem 1 defined in Table 1. The $p + 1 = 4$ DOFs per mesh element are displayed at the final time (bullets) and compared to the exact solution (lines).</figcaption></figure>

<table><caption>TABLE 3. Definitions of the two-dimensional problems and numerical parameters.</caption><thead><tr><th>#</th><th>f(u)</th><th>Ω</th><th>u<sub>0</sub>(x)</th><th>t</th><th>p</th><th>q</th><th>h<sub>Ω</sub></th><th>CFL</th><th>L<sub>f</sub></th><th>L<sub>u</sub></th></tr></thead><tbody><tr><td>7</td><td><sup>1</sup>&frasl;<sub>2</sub>(u<sup>2</sup>, u<sup>2</sup>)<sup>T</sup></td><td>[0, 1]<sup>2</sup></td><td><span style="font-size: 2em; vertical-align: middle;">{</span><table style="display: inline-table; vertical-align: middle; border-collapse: collapse;"><tbody><tr><td style="border: none; padding: 0 0.2em;">1</td><td style="border: none; padding: 0 0.2em;">if ||x - (<sup>5</sup>&frasl;<sub>8</sub>, <sup>5</sup>&frasl;<sub>8</sub>)||<sub>∞</sub> ≤ <sup>1</sup>&frasl;<sub>4</sub>,</td></tr><tr><td style="border: none; padding: 0 0.2em;">-<sup>3</sup>&frasl;<sub>4</sub></td><td style="border: none; padding: 0 0.2em;">else.</td></tr></tbody></table><span style="font-size: 2em; vertical-align: middle;">}</span></td><td><sup>3</sup>&frasl;<sub>8</sub></td><td>3</td><td>3</td><td><sup>1</sup>&frasl;<sub>128</sub></td><td>10</td><td>√2</td><td>1</td></tr><tr><td>8</td><td>(sin u, cos u)<sup>T</sup></td><td>[-2, 2] × [-2.5, 1.5]</td><td><span style="font-size: 2em; vertical-align: middle;">{</span><table style="display: inline-table; vertical-align: middle; border-collapse: collapse;"><tbody><tr><td style="border: none; padding: 0 0.2em;"><sup>7π</sup>&frasl;<sub>2</sub></td><td style="border: none; padding: 0 0.2em;">if |x| ≤ 1,</td></tr><tr><td style="border: none; padding: 0 0.2em;">π&frasl;<sub>4</sub></td><td style="border: none; padding: 0 0.2em;">else.</td></tr></tbody></table><span style="font-size: 2em; vertical-align: middle;">}</span></td><td>1</td><td>3</td><td>3</td><td><sup>1</sup>&frasl;<sub>80</sub></td><td>10</td><td>1</td><td>1</td></tr></tbody></table>

linear equations and is even not needed for a backward-Euler time stepping to ensure the maximum principle for the cell averages ([52], Thm. 3.3). It is however required for high-order time stepping, but with very low values as highlighted in Table 1.

Finally, applying the two-step limiter from Section 5 significantly reduces the dissipation and results in sharp resolution of the solution features, while successfully limiting the solution. In particular the MPP is always ensured at all DOFs and the physical weak solutions is successfully captured. We stress that imposing only the MPP property is not enough to capture the physically admissible solution for problems 5 and 6 and an additional entropy bound in (5.5) is required. Likewise, the cell-average of the high-order solution has been observed to violate the bounds (5.5), so the linear scaling limiter (5.4) alone cannot be applied.

## 6.2. Two-dimensional problems

We now consider problems in two space dimensions, described in Table 3, which we simulate with a fourth-order scheme in space and time and a large CFL value on unstructured meshes with quadratic curved elements (see Fig. 5). We first consider the inviscid Burgers' equation with a discontinuous initial data [26, 50], and then consider the KPP problem [43] defined by a nonconvex flux.

The numerical solutions are depicted in Figures 6 and 7, while Table 4 reports the minimum and maximum values of the DOFs at final time. Once again, the HO scheme results in large oscillations that violate the maximum principle for both the DOFs and the cell averages, thus preventing to apply the linear scaling limiter (5.4). The HO scheme is also unable to capture the physical solution for the KPP problem with a nonconvex flux in Figure 7, though it is ES for the square entropy which is not sufficient. In contrast, the limited scheme successfully captures the physical solution and ensures the maximum principle at the discrete level (see Tab. 4). Note that the iterative solver for the HO scheme was unable to converge for the KPP problem and (5.1) was not satisfied, with the consequence that (5.3a) and (5.3b) were not equivalent. A simple remedy is to limit the DOFs in (5.3a) from the LO solution, instead of the HO solution. This trick was applied for the KPP problem and proved to successfully satisfy the maximum principle as seen in Table 4.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Inviscid Burgers' equation: solutions to problems 2 (top), 3 (middle), and 4 (bottom) defined in Table 1. The $p+1=4$ DOFs per mesh element are displayed at the final time (bullets) and compared to the exact solution (lines).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Buckley-Leverett equation: solutions to problems 5 (top) and 6 (bottom) defined in Table 1. The $p+1=4$ DOFs per mesh element are displayed at the final time (bullets) and compared to the exact solution (lines).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Example of mesh of a square domain with 16 quadratic curved quadrangles used for the numerical tests in two space dimensions.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6. 2D Burgers' equation: solutions to problem 7 in Table 3 at time $t = 0.375$. The solution is post-processed on a fine grid defined by subdividing every cell into four subcells.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 7. KPP problem: solutions to problem 8 in Table 3 at time $t = 1$. The solution is post-processed on a fine grid defined by subdividing every cell into four subcells.</figcaption></figure>

<table><caption>TABLE 4. Minimum, $\min_{\kappa \in \Omega_h} \min_{0 \le ij, \le p} U_{\kappa}^{ij}(t_n^q)$, and maximum, $\max_{\kappa \in \Omega_h} \max_{0 \le ij, \le p} U_{\kappa}^{ij}(t_n^q)$, values of the DOFs at final time for the two-dimensional problems 7 and 8 defined in Table 3. The second and third columns indicate the bounds $m$ and $M$ in (1.3).</caption><thead><tr><th>#</th><th>m</th><th>M</th><th>HO scheme</th><th>LO scheme</th><th>limited scheme</th></tr></thead><tbody><tr><td>7</td><td>-0.75</td><td>1</td><td>-2.386</td><td>2.835</td><td>-0.75</td><td>0.965</td><td>-0.75</td><td>0.993</td></tr><tr><td>8</td><td>0.785</td><td>10.996</td><td>-7.378</td><td>18.028</td><td>0.75</td><td>10.969</td><td>0.785</td><td>10.996</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>FIGURE 8. Evolutions of the entropy inequality (4.8) integrated over the space domain $\Omega_h$ at every iteration for the square, $\eta(u) = \frac{u^2}{2}$, and Kružkov's, $\eta(u) = |u|$, entropies for 2D problems in Table 3: 2D Burgers' equation (**left**) and KPP problem (**right**).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 9. Cell indicator (black cells) where the entropy inequality (4.8) for the Kružkov's entropy, $\eta(u) = |u|$, is not satisfied at the final iteration of the different computations of 2D problems in Table 3: 2D Burgers' equation (**top**) and KPP problem (**bottom**).</figcaption></figure>

Finally, Figures 8 and 9 indicate whether or not the entropy inequalities (4.8) hold globally and locally, respectively. The results are in good agreement with the theoretical results and confirm that the LO scheme is entropy stable for both entropies, while the HO scheme is entropy stable for the square entropy only (see Lem. 4.1) and violates locally and globally the entropy inequality for the Kružkov's entropy. The limited scheme improves entropy dissipation compared to the HO scheme, but may violate entropy stability locally as observed in Figure 9. This may be attributed to the shock indicator used to impose the entropy bound in the FCT limiter (see Sect. 5.3), which may fail to strictly satisfy the entropy inequality. Even though not fully entropy stable, the limited scheme captures the physical solution unlike the HO scheme.
## 7. CONCLUDING REMARKS

We propose and analyze artificial viscosities in DGSEM schemes with implicit time stepping for the discretization of nonlinear scalar conservation laws in multiple space dimensions. We consider both a backward Euler time stepping for steady-state simulations, and a space-time DGSEM for time resolved simulations. The artificial dissipation is a first-order graph viscosity local to the space-time discretization elements. The schemes have no time step restriction, are MPP, satisfy a fully discrete inequality for every admissible convex entropy, and the associated discrete problems are well-posed. The resulting low-order scheme are then used as building blocks to propose a two-step limiter that is designed to impose the maximum principle and entropy stability, for all entropies required by the user, to the high-order solution, while keeping its accuracy for smooth solutions. Numerical experiments are proposed to illustrate the properties of these schemes. Even being ES, the DGSEM without graph viscosity is not MPP and may fail to capture the entropy weak solution. Applying the two-step limiter keeps robustness and stability of the low-order scheme, while it keeps the resolution capabilities of the high-order scheme and captures the physical solution. Future work will focus on the extension of this approach to systems of conservation laws. Though general discrete counterparts (2.14) to (1.6) for the EC fluxes seem difficult to achieve since the EC fluxes are not uniquely defined the case of systems, the analysis of the MPP and ES properties of the low-order scheme can still be carried out without these properties (see Rem. 3.4), which allows to introduce subiterations that preserve convex invariant domains of the PDEs as in (3.5) and (4.7). The design of limiters adapted to systems and avoiding the computation of a low-order solution will constitute another direction of research.

## DATA AVAILABILITY STATEMENT

No new data/codes were created or analyzed in this study.

## REFERENCES

[1] R. Abgrall, P. Öffner and H. Ranocha, Reinterpretation and extension of entropy correction terms for residual distribution and discontinuous Galerkin schemes: application to structure preserving discretization. *J. Comput. Phys.* **453** (2022) 110955.

[2] H.L. Atkins and C.-W. Shu, Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations. *AIAA J.* **36** (1998) 775-782.

[3] J.P. Boris and D.L. Book, Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works. *J. Comput. Phys.* **11** (1973) 38-69.

[4] F. Bouchut and B. Perthame, Kružkov's estimates for scalar conservation laws revisited. *Trans. Amer. Math. Soc.* **350** (1998) 2847-2870.

[5] V. Carlier and F. Renac, Invariant domain preserving high-order spectral discontinuous approximations of hyperbolic systems. *SIAM J. Sci. Comput.* **45** (2023) A1385-A1412.

[6] M.H. Carpenter, T.C. Fisher, E.J. Nielsen and S.H. Frankel, Entropy stable spectral collocation schemes for the Navier-Stokes equations: discontinuous interfaces. *SIAM J. Sci. Comput.* **36** (2014) B835-B867.

[7] T. Chen and C.-W. Shu, Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws. *J. Comput. Phys.* **345** (2017) 427-461.

[8] B. Cockburn, F. Coquel and P.G. LeFloch, Convergence of the finite volume method for multidimensional conservation laws. *SIAM J. Numer. Anal.* **32** (1995) 687-705.

[9] F. Coquel, C. Marmignon, P. Rai and F. Renac, An entropy stable high-order discontinuous Galerkin spectral element method for the Baer-Nunziato two-phase flow model. *J. Comput. Phys.* **431** (2021) 110135.

[10] J. Crean, J.E. Hicken, D.C. Del Rey Fernández, D.W. Zingg and M.H. Carpenter, Entropy-stable summation-by-parts discretization of the Euler equations on general curved elements. *J. Comput. Phys.* **356** (2018) 410-438.

[11] C.M. Dafermos, Hyperbolic conservation laws in continuum physics, in Grundlehren der mathematischen Wissenschaften. Springer Berlin Heidelberg, Berlin, Heidelberg (2016).

[12] B. Després, Entropy inequality for high order discontinuous Galerkin approximation of Euler equations, in Hyperbolic Problems: Theory, Numerics, Applications (Basel), edited by M. Fey and R. Jeltsch. Birkhäuser Basel (1999) 225-231.
[13] B. Després, Discontinuous Galerkin method for the numerical solution of Euler equations in axisymmetric geometry, in Discontinuous Galerkin Methods: Theory, Computation and Applications. Springer (2000) 315–320.

[14] B. Engquist and S. Osher, One-sided difference approximations for nonlinear conservation laws. *Math. Comp.* **36** (1981) 321–351.

[15] A. Ern and J.-L. Guermond, Invariant-domain preserving high-order time stepping: II. IMEX schemes. *SIAM J. Sci. Comput.* **45** (2023) A2511–A2538.

[16] R. Eymard, T. Gallouët and R. Herbin, Finite volume methods, in Solution of Equation in $\mathbb{R}^n$ (Part 3), Techniques of Scientific Computing (Part 3), Handbook of Numerical Analysis. Vol. 7. Elsevier (2000) 713–1018.

[17] T.C. Fisher and M.H. Carpenter, High-order entropy stable finite difference schemes for nonlinear conservation laws: finite domains. *J. Comput. Phys.* **252** (2013) 518–557.

[18] L. Friedrich, G. Schnücke, A.R. Winters, D.C. Del Rey Fernández, G.J. Gassner and M.H. Carpenter, Entropy stable space-time discontinuous Galerkin schemes with summation-by-parts property for hyperbolic conservation laws. *J. Sci. Comput.* **80** (2019) 175–222.

[19] K.O. Friedrichs and P.D. Lax, Systems of conservation equations with a convex extension. *Proc. Nat. Acad. Sci. USA* **68** (1971) 1686–1688.

[20] G.J. Gassner, A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods. *SIAM J. Sci. Comput.* **35** (2013) A1233–A1253.

[21] G. Gassner and D.A. Kopriva, A comparison of the dispersion and dissipation errors of Gauss and Gauss-Lobatto discontinuous Galerkin spectral element methods. *SIAM J. Sci. Comput.* **33** (2011) 2560–2579.

[22] G.J. Gassner, A.R. Winters and D.A. Kopriva, Split form nodal discontinuous Galerkin schemes with summation-by-parts property for the compressible Euler equations. *J. Comput. Phys.* **327** (2016) 39–66.

[23] E. Godlewski and P.-A. Raviart, Hyperbolic Systems of Conservation Laws. Ellipses (1991).

[24] S. Godunov, A difference scheme for numerical computation of discontinuous solutions of equations of fluid dynamics. *Math. USSR Sbornik* **47** (1959) 271–306.

[25] J.-L. Guermond and B. Popov, Invariant domains and first-order continuous finite element approximation for hyperbolic systems. *SIAM J. Numer. Anal.* **54** (2016) 2466–2489.

[26] J.-L. Guermond and B. Popov, Invariant domains and second-order continuous finite element approximation for scalar conservation equations. *SIAM J. Numer. Anal.* **55** (2017) 3120–3146.

[27] J.-L. Guermond, M. Nazarov, B. Popov and I. Tomas, Second-order invariant domain preserving approximation of the Euler equations using convex limiting. *SIAM J. Sci. Comput.* **40** (2018) A3211–A3239.

[28] J.-L. Guermond, M. Maier, B. Popov and I. Tomas, Second-order invariant domain preserving approximation of the compressible Navier–Stokes equations. *Comput. Methods Appl. Mech. Eng.* **375** (2021) 113608.

[29] A. Harten, P.D. Lax and B. van Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. *SIAM Rev.* **25** (1983) 35–61.

[30] S. Hennemann, A.M. Rueda-Ramírez, F.J. Hindenlang and G.J. Gassner, A provably entropy stable subcell shock capturing approach for high order split form DG for the compressible Euler equations. *J. Comput. Phys.* **426** (2021) 109935.

[31] A. Hiltebrand and S. Mishra, Entropy stable shock capturing space-time discontinuous Galerkin schemes for systems of conservation laws. *Numer. Math.* **126** (2014) 103–151.

[32] A. Jameson, Positive schemes and shock modelling for compressible flows. *Int. J. Numer. Meth. Fluids* **20** (1995) 743–776.

[33] G.S. Jiang and C.-W. Shu, On a cell entropy inequality for discontinuous Galerkin methods. *Math. Comp.* **62** (1994) 531–538.

[34] G.-S. Jiang and C.-W. Shu, Efficient implementation of weighted ENO schemes. *J. Comput. Phys.* **126** (1996) 202–228.

[35] Y. Jiang and H. Liu, Invariant-region-preserving DG methods for multi-dimensional hyperbolic conservation law systems, with an application to compressible Euler equations. *J. Comput. Phys.* **373** (2018) 385–409.

[36] K. Jittornerum, An implicit function theorem. *J. Optim. Theory Appl.* **25** (1978) 575–577.

[37] B. Khobalatte and B. Perthame, Maximum principle on the entropy and second-order kinetic schemes. *Math. Comp.* **62** (1994) 119–131.

[38] D.A. Kopriva, Metric identities and the discontinuous spectral element method on curvilinear meshes. *J. Sci. Comput.* **26** (2006) 302–327.

[39] D.A. Kopriva, Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists And Engineers. Springer Dordrecht (2009).
[40] D.A. Kopriva and G. Gassner, On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods. *J. Sci. Comput.* **44** (2010) 136–155.

[41] L. Krivodonova and R. Qin, An analysis of the spectrum of the discontinuous Galerkin method. *Appl. Numer. Math.* **64** (2013) 1–18.

[42] S. Kružkov, First order quasilinear equations in several independent variables. *Math. Ussr Sbornik* **10** (1970) 217–243.

[43] A. Kurganov, G. Petrova and B. Popov, Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws. *SIAM J. Sci. Comput.* **29** (2007) 2381–2401.

[44] D. Kuzmin, Explicit and implicit FEM-FCT algorithms with flux linearization. *J. Comput. Phys.* **228** (2009) 2517–2534.

[45] D. Kuzmin and S. Turek, Flux correction tools for finite elements. *J. Comput. Phys.* **175** (2002) 525–558.

[46] N.N. Kuznetsov, Accuracy of some approximate methods for computing the weak solutions of a first-order quasi-linear equation. *USSR Comp. Math. Math. Phys.* **16** (1976) 105–119.

[47] P.D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical computation. *Comm. Pure Appl. Math.* **7** (1954) 159–193.

[48] H. Le Dret, Fixed Point Theorems and Applications. Springer International Publishing, Cham (2018) 47–68.

[49] P.G. Le Floch, J.-M. Mercier and C. Rohde, Fully discrete, entropy conservative schemes of arbitrary order. *SIAM J. Numer. Anal.* **40** (2002) 1968–1992.

[50] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems. *Cambridge Texts in Applied Mathematics*. Cambridge University Press (2002).

[51] D. Ling, J. Cheng and C.-W. Shu, Conservative high order positivity-preserving discontinuous Galerkin methods for linear hyperbolic and radiative transfer equations. *J. Sci. Comput.* **77** (2018) 1801–1831.

[52] R. Milani, F. Renac and J. Ruel, Maximum principle preserving time implicit DGSEM for linear scalar hyperbolic conservation laws. *J. Comput. Phys.* (2024) 113254.

[53] O.A. Olešnik, Uniqueness and stability of the generalized solution of the Cauchy problem for a quasi-linear equation. *Uspekhi Mat. Nauk* **14** (1959) 165–170.

[54] E.Y. Panov, Uniqueness of the solution of the Cauchy problem for a first order quasilinear equation with one admissible strictly convex entropy. *Math. Notes* **55** (1994) 517–525.

[55] W. Pazner, Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting. *Comput. Methods Appl. Mech. Eng.* **382** (2021) 113876.

[56] P.-O. Persson and J. Peraire, Sub-cell shock capturing for discontinuous Galerkin methods, in 44th AIAA Aerospace Sciences Meeting and Exhibit (2006).

[57] T. Qin and C.-W. Shu, Implicit positivity-preserving high-order discontinuous Galerkin methods for conservation laws. *SIAM J. Sci. Comput.* **40** (2018) A81–A107.

[58] F. Renac, A robust high-order discontinuous Galerkin method with large time steps for the compressible Euler equations. *Commun. Math. Sci.* **15** (2017) 813–837.

[59] F. Renac, A robust high-order Lagrange-projection like scheme with large time steps for the isentropic Euler equations. *Numer. Math.* **135** (2017) 493–519.

[60] F. Renac, Entropy stable DGSEM for nonlinear hyperbolic systems in nonconservative form with application to two-phase flows. *J. Comput. Phys.* **382** (2019) 1–26.

[61] A.M. Rueda-Ramírez, W. Pazner and G.J. Gassner, Subcell limiting strategies for discontinuous Galerkin spectral element methods. *Comput. Fluids* **247** (2022) 105627.

[62] J. Schauder, Der Fixpunktsatz in Funktionalräumen. *Studia Math.* **2** (1930) 171–180.

[63] B. Strand, Summation by parts for finite difference approximations for $d/dx$. *J. Comput. Phys.* **110** (1994) 47–67.

[64] A. Szepessy, Convergence of a streamline diffusion finite element method for scalar conservation laws with boundary conditions. *ESAIM: M2AN* **25** (1991) 749–782.

[65] E. Tadmor, The numerical viscosity of entropy stable schemes for systems of conservation laws. *I. Math. Comp.* **49** (1987) 91–103.

[66] E. Tadmor, Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems. *Acta Numer.* **12** (2003) 451–512.

[67] J.J.W. van der Vegt, Y. Xia and Y. Xu, Positivity preserving limiters for time-implicit higher order accurate discontinuous Galerkin discretizations. *SIAM J. Sci. Comput.* **41** (2019) A2037–A2063.

[68] F. Vilar and R. Abgrall, A posteriori local subcell correction of high-order discontinuous Galerkin scheme for conservation laws on two-dimensional unstructured grids. *SIAM J. Sci. Comput.* **46** (2024) A851–A883.
[69] F. Vilar, C.-W. Shu and P.-H. Maire, Positivity-preserving cell-centered Lagrangian schemes for multi-material compressible flows: from first-order to high-orders. Part II: the two-dimensional case. *J. Comput. Phys.* **312** (2016) 416–442.

[70] A.I. Vol'pert, The spaces BV and quasilinear equations. *Math. USSR Sbornik* **115** (1967) 255–302.

[71] N. Wintermeyer, A.R. Winters, G.J Gassner and D.A. Kopriva, An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry. *J. Comput. Phys.* **340** (2017) 200–242.

[72] Z. Xu and C.-W. Shu, High order conservative positivity-preserving discontinuous Galerkin method for stationary hyperbolic equations. *J. Comput. Phys.* **466** (2022) 111410.

[73] S.T. Zalesak, Fully multidimensional flux-corrected transport algorithms for fluids. *J. Comput. Phys.* **31** (1979) 335–362.

[74] X. Zhang and C.-W. Shu, On maximum-principle-satisfying high order schemes for scalar conservation laws. *J. Comput. Phys.* **229** (2010) 3091–3120.

[75] X. Zhang, Y. Xia and C.-W. Shu, Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes. *J. Sci. Comput.* **50** (2012) 29–62.

<figure><img src="image_12.png" /></figure>

Please help to maintain this journal in open access!

This journal is currently published in open access under the **Subscribe to Open model (S2O)**. We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.