*hp*-optimal convergence of the original DG method for
linear hyperbolic problems on special simplicial meshes

Zhaonan Dong, Lorenzo Mascotto

► To cite this version:

Zhaonan Dong, Lorenzo Mascotto. *hp*-optimal convergence of the original DG method for linear hyperbolic problems on special simplicial meshes. IMA Journal of Numerical Analysis, In press, 10.1093/imanum/drae051. hal-04252310v3

HAL Id: hal-04252310

https://inria.hal.science/hal-04252310v3

Submitted on 9 Jul 2024

**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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
*hp*-optimal convergence of the original DG method for linear hyperbolic problems on special simplicial meshes

Z. Dong,\* L. Mascotto†

**Abstract**

We prove *hp*-optimal error estimates for the original discontinuous Galerkin (DG) method when approximating solutions to first-order hyperbolic problems with constant convection fields in the *L*² and DG norms. The main theoretical tools used in the analysis are novel *hp*-optimal approximation properties of the special projector introduced in [Cockburn, Dong, Guzmán, SINUM, 2008]. We assess the theoretical findings on some test cases.

**AMS subject classification:** 65N12; 65N30; 65N50.

**Keywords:** discontinuous Galerkin; optimal convergence; linear hyperbolic problems; a priori error estimation; *p*-version.

# 1 Introduction

The discontinuous Galerkin (DG) method was first introduced in [21] for the approximation of solutions to the neutron transport equation, i.e., a first-order linear hyperbolic problem. This method, which we shall refer to as the original DG method, was later analysed in [18]; advances on the convergence of the scheme were given in [15]. The convergence rates in the *L*² norm of the *h*-version of the method proven therein are suboptimal by half an order; this was also apparent from the numerical experiments in [20]. On special classes of meshes, optimal convergence for the *h*-version was established in [8].

Fewer results are available for the *p*- and *hp*-versions of the method. In [13], *hp* error estimates were derived for a streamline diffusion version of the discontinuous Galerkin method on fairly general quadrilateral meshes. Shortly after, in [14], for the original method and linear polynomial convection fields, *p* convergence was shown suboptimal by half an order for the convection-diffusion case, but optimal for the linear hyperbolic case; for more general vector fields, suboptimality by one order and a half was discussed as well. The main theoretical tools in this reference are *hp* optimal approximation properties of the *L*² projector on the boundary of tensor product elements. More recently [11] and again for tensor product elements, such a suboptimality was reduced to half an order only for a special class of convective fields.

This paper aims at extending even more the knowledge on the convergence analysis of the *hp*-discontinuous Galerkin (*hp*-DG) for first order hyperbolic problems in the following aspects: for constant convection fields,

* we prove *hp*-optimal convergence of the method in the *L*² norm, thus generalizing the results in [8] to the *p*-version of the method;

* we prove *hp*-optimal convergence of the method in the DG norm also on special simplicial meshes as in [8], thus generalizing the results in [13] to the case of simplicial meshes.

To this aim, we analyse the *hp* approximation properties of the Cockburn-Dong-Guzmán (henceforth denoted by CDG) projector introduced in [9]. In particular, we generalise the one dimensional results in [22, Lemmas 3.5 and 3.6] to simplices in two and three dimensions.

<footnote>\*Inria, 2 rue Simone Iff, 75589 Paris, France and CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée 2, France, zhaonan.dong@inria.fr</footnote>

<footnote>†Dipartimento di Matematica e Applicazioni, Università di Milano Bicocca, 20125 Milan, Italy, lorenzo.mascotto@unimib.it; IMATI-CNR, 27100, Pavia, Italy; Fakultät für Mathematik, Universität Wien, 1090 Vienna, Austria, lorenzo.mascotto@univie.ac.at</footnote>
A possible reason why the simplicial mesh case was not contemplated in [13, 14] is that hp-optimal convergence estimates for the trace of the polynomial L²-projection operator on simplices, crucial in the analysis of the CDG projector properties were derived later [7, 19].

**Notation.** Let $D$ be a Lipschitz domain in $\mathbb{R}^d$, $d = 1, 2$, and $3$, with boundary $\partial D$. The space of Lebesgue measurable and square integrable functions over $D$ is $L^2(D)$. The Sobolev space of positive integer order $s$ is $H^s(D)$. We also write $L^2(D) = H^0(D)$. We endow $H^s(D)$ with the inner product, seminorm, and norm

$$ (\cdot, \cdot)_{s,D}, \qquad |\cdot|_{s,D}, \qquad \|\cdot\|_{s,D}. $$

Interpolation theory is used to construct Sobolev spaces of positive noninteger order; duality is used to define negative order Sobolev spaces.

The trace theorem is valid for $H^s(D)$, $1/2 < s < 3/2$. In particular, given $s$ in the above range and $g$ in $H^{s-1/2}(\Gamma_g)$, being $\Gamma_g$ any subset of $\partial D$ with nonzero measure in $\partial D$, we are allowed to define the space

$$ H_g^s(D, \Gamma_g) := \{v \in H^s(D, \Gamma_D) \mid v|_{\Gamma_D} = g\}. $$

The space of polynomials of nonnegative degree $p$ over $D$ is $\mathbb{P}_p(D)$.

**The continuous problem.** Given $\Omega$ a Lipschitz domain in $\mathbb{R}^d$, $d = 1, 2, 3$, with boundary $\Gamma$, consider $\beta$ in $\mathbb{R}^d$ and $c$ in $L^\infty(\Omega)$. We introduce the inflow part $\Gamma_-$ of the boundary of $\Omega$ as follows: given $\mathbf{n}_\Gamma(\mathbf{x})$ the outward normal to $\Gamma$ at $\mathbf{x}$,

$$ \Gamma_- := \{\mathbf{x} \in \partial\Omega \mid \beta \cdot \mathbf{n}_\Gamma(\mathbf{x}) < 0\}. $$

The characteristic part $\Gamma_0$ is analogously defined as

$$ \Gamma_0 := \{\mathbf{x} \in \partial\Omega \mid \beta \cdot \mathbf{n}_\Gamma(\mathbf{x}) = 0\}. $$

The outflow part $\Gamma_+$ of $\Gamma$ is given by $\Gamma \setminus (\Gamma_- \cup \Gamma_0)$. We are interested in the approximation of solutions to convection-reaction problems: given $g$ in $L^2(\Gamma_-)$,

$$ \left\{ \begin{array}{ll} \text{find } u \text{ such that} & \\ \beta \cdot \nabla u + c u = f & \text{in } \Omega \\ u = g & \text{on } \Gamma_- \end{array} \right. \qquad (1) $$

Introduce the graph space

$$ V := \{v \in L^2(\Omega) \mid \beta \cdot \nabla v \in L^2(\Omega)\}, $$

which we endow with the graph norm

$$ \|v\|_V^2 := \|v\|_{0,\Omega}^2 + \|\beta \cdot \nabla v\|_{0,\Omega}^2. \qquad (2) $$

We can define a trace operator from the graph space onto the space

$$ L^2(|\beta \cdot \mathbf{n}_\Gamma|, \Gamma) := \{v \text{ measurable on } \Gamma \mid (|\beta \cdot \mathbf{n}_\Gamma|v, v)_{0,\Gamma} < \infty\}. $$

Let

$$ V_g := \{v \in V \mid g \text{ is the trace of } v \text{ and belongs to } L^2(|\beta \cdot \mathbf{n}_\Gamma|, \Gamma)\}. \qquad (3) $$

We introduce the bilinear form on $V_g \times L^2(\Omega)$ as

$$ b(u, v) = (\beta \cdot \nabla u, v)_{0,\Omega} + (c u, v)_{0,\Omega}. $$

The weak formulation of problem (1) reads

$$ \left\{ \begin{array}{ll} \text{find } u \in V_g \text{ such that} & \\ b(u, v) = (f, v)_{0,\Omega} & \quad \forall v \in L^2(\Omega). \end{array} \right. \qquad (4) $$

Henceforth, we assume that there exists a positive constant $\bar{c}_0$ such that

$$ \bar{c} := c - \frac{1}{2} \nabla \cdot \beta \geq \bar{c}_0. \qquad (5) $$

Problem (4) is well posed with respect to the graph norm in (2); see, e.g., [3].
**Structure of the paper.** In Section 2, we introduce admissible simplicial meshes in the sense of [8] and recall the original DG method from [21]. The main technical result of the paper, i.e., hp approximation properties of the CDG projector are derived in Section 3 in one, two, and three dimensions. Such estimates are used to derive hp-optimal convergence of the method in the L² norm in Section 4; there, we also show hp-optimal convergence of the method in a DG norm. We assess the theoretical results with several numerical experiments in Section 5 and draw some conclusions in Section 6.

# 2 Admissible meshes and the method

We introduce admissible simplicial meshes for the forthcoming analysis, and define associated Sobolev and polynomial broken spaces in Section 2.1, and recall the original DG method for problem (4) in Section 2.2.

## 2.1 Admissible simplicial meshes

We follow [8, Section 1] and introduce admissible simplicial meshes that are instrumental in deriving the main result of the paper; see Section 3 below.

Let $\mathcal{T}_h$ be a simplicial mesh of the domain $\Omega$ and $\mathcal{F}_h$ be its set of $(d-1)$-dimensional facets. We distinguish the facets in $\mathcal{F}_h$ into internal and boundary facets; the former are those $F$ not contained in $\partial\Omega$. The union of the interior facets of $\mathcal{T}_h$ is $\Gamma_h^{\text{in}}$.

With each element $T$ of $\mathcal{T}_h$ we associate its diameter $h_T$ and its outward unit vector $\mathbf{n}_T$, which is defined almost everywhere on $\partial T$. The set of $(d-1)$-dimensional facets of $T$ is $\mathcal{F}_T$. With each facet $F$ in $\mathcal{F}_h$ we associate its diameter $h_F$ and a unit normal vector $\mathbf{n}_F$, which is pointing outward $\Omega$ if $F$ is a boundary facet. We only consider meshes $\mathcal{T}_h$ that are $\sigma$ shape-regular i.e., any element $T$ of $\mathcal{T}_h$ is star-shaped with respect to a ball of radius larger than or equal to $h_T\sigma$. We associate with each mesh $\mathcal{T}_h$ piecewise differential operators by adding a subscript $h$, for instance $\nabla_h$ denotes the piecewise gradient.

We distinguish the facets in $\mathcal{F}_T$ into outflow, inflow, and characteristics (with respect to $\beta$) facets, depending on whether they satisfy either of the two following properties:

$$ \mathbf{n}_{T|F} \cdot \beta > 0, \qquad \mathbf{n}_{T|F} \cdot \beta < 0, \qquad \mathbf{n}_{T|F} \cdot \beta = 0. $$

The union of the outflow and inflow facets is $\Gamma_{\text{out}}^T$ and $\Gamma_{\text{in}}^T$, respectively.

We require that each element $T$ of $\mathcal{T}_h$ satisfy the two following properties:

(A1) $T$ has only one outflow facet, which we denote by $F_T^+$;

(A2) each interior facet $F$ that is an inflow facet for $T$ is included in the outflow facet for another simplex in the mesh.

The second assumption above implies that $\mathcal{T}_h$ can be nonconforming, i.e., hanging facets are allowed. Meshes satisfying the above properties can be always constructed; see, e.g., [8, Appendix].

The focus of the paper is on p-optimal estimates; therefore, we pick quasi-uniform meshes and denote the mesh size of $\mathcal{T}_h$ by $h$.

Given $p$ in $\mathbb{N}$ and $s$ in $\mathbb{R}^+$, we associate with a mesh $\mathcal{T}_h$ as in Section 2.1 the spaces

$$ \mathbb{P}_p(\mathcal{T}_h) := \{q_p \in L^2(\Omega) \mid q_{p|T} \in \mathbb{P}_p(T) \quad \forall T \in \mathcal{T}_h\} $$

and

$$ H^s(\mathcal{T}_h) := \{v \in L^2(\Omega) \mid v|_T \in H^s(T) \quad \forall T \in \mathcal{T}_h\}. $$

Differential operators defined piecewise over $\mathcal{T}_h$ are denoted with the same symbol of the original operator with an extra subscript $h$; for instance, the broken gradient is $\nabla_h$.

Given an interior facet $F$, we denote by $T^+$ and $T^-$ the two elements of $\mathcal{T}_h$ such that $F$ is an outflow and inflow facet for $T^+$ and $T^-$, respectively. Given $v_h$ in $\mathbb{P}_p(\mathcal{T}_h)$ and $\mathbf{x}$ in the face $F$, we write

$$ v_h^{\pm}(\mathbf{x}) = \lim_{\delta \downarrow 0} v_h(\mathbf{x} \pm \delta\beta). $$
We introduce the jump operator $[\cdot] : H^1(\mathcal{T}_h) \to L^2(\Gamma_h^{\text{in}})$ given by $[\cdot]_{|F} := v_h^+ - v_h^-$ for all internal edges $F$.

Given a boundary facet $F$, we denote by $T$ the only element of $\mathcal{T}_h$ such that $F$ belongs to $\mathcal{F}_T$. Since every function $v_h$ in $\mathbb{P}_p(\mathcal{T}_h)$ is single valued on boundary facets, for all $\mathbf{x}$ in the face $F$, we write $v_h(\mathbf{x}) = v_{hT}(\mathbf{x})$, where $T$ is the only element of $\mathcal{T}_h$ such that $F$ belongs to $\mathcal{F}_T$.

Henceforth, given two positive quantities $a$ and $b$, we write $a \lesssim b$ is there exists a positive constant $c$ possibly depending on the shape-regularity parameter $\sigma$ only, such that $a \le c b$. If $a \lesssim b$ and $b \lesssim a$ at once, we write $a \approx b$.

## 2.2 The original DG method

Consider the space $V_h$ given by $\mathbb{P}_p(\mathcal{T}_h)$. The original [21] DG method for (4) reads

$$ \begin{cases} \text{find } u_h \in V_h \text{ such that} \\ \mathcal{B}(u_h, v_h) = (f, v_h)_{0,\Omega} - (g, \mathbf{n}_\Gamma \cdot \boldsymbol{\beta} \, v_h)_{0,\Gamma_-} & \forall v_h \in V_h. \end{cases} \tag{6} $$

where

$$ \mathcal{B}(u_h, v_h) := (\boldsymbol{\beta} \cdot \nabla_h u_h + c \, u_h, v_h)_{0,\Omega} - \sum_{T \in \mathcal{T}_h} (\llbracket u_h \rrbracket, \mathbf{n}_T \cdot \boldsymbol{\beta} \, v_h^+)_{0,\Gamma_{in}^T} - (u_h, \mathbf{n}_\Gamma \cdot \boldsymbol{\beta} \, v_h)_{0,\Gamma_-}. $$

The solution $u_h$ can be computed starting at the inflow boundary $\Gamma_-$; solving local problems derived from (6) on the elements abutting $\Gamma_-$; transmitting the solution to the neighbouring elements in the $\boldsymbol{\beta}$ direction through upwind.

An equivalent alternative formulation used, e.g., in [8] is derived using an integration by parts and the fact that $\boldsymbol{\beta}$ is a constant field. Notably, an integration by parts implies that the bilinear form $\mathcal{B}(\cdot, \cdot)$ can be rewritten as

$$ \mathcal{B}(u_h, v_h) = -(u_h, \boldsymbol{\beta} \cdot \nabla_h v_h)_{0,\Omega} + \sum_{T \in \mathcal{T}_h} (u_h^-, \mathbf{n}_T \cdot \boldsymbol{\beta} \llbracket v_h \rrbracket)_{0,\Gamma_{in}^T \setminus \Gamma_-} + (c \, u_h, v_h)_{0,\Omega} + (u_h, \mathbf{n}_\Gamma \cdot \boldsymbol{\beta} \, v_h)_{0,\Gamma_+}. \tag{7} $$

Method (6) is well posed and coincides with that introduced in [21] and analysed in [18]. A priori estimates with respect to the data are derived in [13, Lemma 2.4 and Section 5]. Previous results requiring stronger assumptions on $c$ are given, e.g., in [15].

**Proposition 2.1.** Let $\bar{c}$ be given in (5). The following bound holds true:

$$ \bar{c} \|u_h\|_{0,\Omega}^2 + \sum_{T \in \mathcal{T}_h} \left\| \mathbf{n}_T \cdot \boldsymbol{\beta}^{\frac{1}{2}} \llbracket u_h \rrbracket \right\|_{0,\Gamma_{in}^T}^2 + \left\| \mathbf{n}_\Gamma \cdot \boldsymbol{\beta}^{\frac{1}{2}} u_h \right\|_{0,\Gamma}^2 \le \bar{c}^{-1} \|f\|_{0,\Omega} + 2 \|g\|_{0,\Gamma_-}^2. $$

All constants appearing in the stability estimates of Proposition 2.1 are explicit, and independent of $h$ and $p$. This is essential to derive the optimal error estimates in Section 4 below. Former stability results display stability estimates that are suboptimal in the polynomial degree; see, e.g., [15, Theorem 2.1].

Method (6) is consistent. The following Galerkin orthogonality property follows: given $u$ and $u_h$ the solutions to (4) and (6),

$$ \mathcal{B}(u - u_h, v_h) = 0 \qquad \forall v_h \in V_h. \tag{8} $$

In particular, the bilinear form $\mathcal{B}(\cdot, \cdot)$ is also well defined on $V_g \times V_h$, where $V_g$ is defined in (3).

*Remark 1.* The case of nonconstant $\boldsymbol{\beta}$ is considered in [10]. Further assumptions on the advection field are needed. The presentation of this manuscript might be generalised to a similar setting at the price of extra technicalities.

# 3 The CDG projector and its approximation properties

We recall some notation and technical results from the theory of orthogonal polynomials in Section 3.1; we define of the Cockburn-Dong-Guzmán (CDG) projector and state the main result in Section 3.2; we prove $p$-optimal approximation properties of the CDG operator in one, two, and three dimensions in Sections 3.3, 3.4, and 3.5, respectively; in Section 3.6 we prove trace-type estimates for the CDG projectors.
## 3.1 Preliminary results

We recall hp approximation properties of the $L^2$ projector on a simplex T in the $L^2$ norm on the boundary of T. This result traces back to [7, Theorem 2.1] and [19, Theorem 1.1]. Let $\Pi_p^0: L^2(T) \to \mathbb{P}_p(T)$ denote the $L^2$ projector defined as

$$ (v - \Pi_p^0 v, q_p)_{0,T} = 0 \qquad \forall q_p \in \mathbb{P}_p(T). \tag{9} $$

**Lemma 3.1.** Under the notation of Section 2.1, there exists a positive C independent of $\sigma$ and d such that

$$ \|u - \Pi_p^0 u\|_{0, F_T^+} \le C \left(\frac{h_T}{p}\right)^{\frac{1}{2}} |u|_{1,T} \qquad \forall u \in H^1(T). $$

As a consequence, for a possibly larger constant C and k positive, we also have

$$ \|u - \Pi_p^0 u\|_{0, F_T^+} \le C \left(\frac{h_T}{p}\right)^{k+\frac{1}{2}} \|u\|_{k+1, T} \qquad \forall u \in H^{k+1}(T). $$

In the one dimensional case, the norm on the right-hand side of the last inequality is in fact a seminorm.

As we are interested in deriving estimates explicit in the polynomial degree and the estimates explicit in the element size are already known from [9], we shall prove approximation properties on reference elements.

Set $\hat{I} := [-1, 1]$ the reference interval. We recall some properties of orthogonal polynomials; see, e.g., [23]. Let $\{L_j\}_{j=0}^{+\infty}$ be the $L^2$ orthogonal set of the Legendre polynomials over $\hat{I}$ satisfying [23, Corollary 3.6, $\alpha = 0, \beta = 0$]

$$ (L_i, L_j)_{0, \hat{I}} = \frac{2}{2j+1} \delta_{i,j} \qquad \forall i, j \in \mathbb{N}_0. \tag{10} $$

Given $\Gamma(\cdot)$ the Gamma function, further let $\{J_j^\ell\}_{j=0}^{+\infty}, \ell > -1$, be the weighted-$L^2$ orthogonal set of Jacobi polynomials over $\hat{I}$ satisfying [23, Corollary 3.6, $\alpha = \ell, \beta = 0$]

$$ ((1-x)^\ell J_i^\ell, J_j^\ell)_{0, \hat{I}} = \frac{2^{\ell+1} \Gamma(j+\ell+1)\Gamma(\ell+1)}{2j+\ell+1 \Gamma(\ell+1)\Gamma(j+\ell+1)} \delta_{j,\ell} = \frac{2^{\ell+1}}{2j+\ell+1} \delta_{j,\ell} \qquad \forall i,j \in \mathbb{N}_0. \tag{11} $$

## 3.2 The CDG projector and the main result

The CDG operator was introduced for the analysis of superconvergent DG methods for second-order elliptic problems in [9]. It can be defined for simplices in any dimensions satisfying assumptions (A1)-(A2). In particular, recall that $F_T^+$ is the only outflow facet of T.

The CDG operator $P_p: H^{\frac{1}{2}+\varepsilon}(T) \to \mathbb{P}_p(T)$ is defined as

$$ (v - P_p v, q_{p-1}^T)_{0,T} = 0 \qquad \forall q_{p-1}^T \in \mathbb{P}_{p-1}(T), \tag{12a} $$

$$ (v - P_p v, q_p^{F_T})_{0, F_T^+} = 0 \qquad \forall q_p^{F_T} \in \mathbb{P}_p(F_T^+). \tag{12b} $$

In [9], the following error estimate was proven for sufficiently smooth functions:

$$ \|v - P_p v\|_{0,T} \le C h^{p+1} |v|_{p+1, \Omega}. $$

The constant C above depends on $\sigma$ and the polynomial degree p. The remainder of the section is devoted to carefully detail the dependence on p. Namely, we shall prove the following result.

**Theorem 3.2.** Let the assumptions (A1) and (A2) be valid, T be in $\mathcal{T}_h$, and $\Pi_p^0$ be the operator in (9). Then, there exist a positive constant C independent of h and p, such that, for any v in $H^{k+1}(T)$, k positive, we have

$$ \|v - P_p v\|_{0,T} \le C \frac{h^{\min(p,k)+1}}{p^{k+1}} \|v\|_{k+1, T}. $$
### 3.3 The 1D case

We prove Theorem 3.2 in the case of one dimensional simplices, i.e., on intervals. This was already done in [22, Lemmas 3.5 and 3.6]. However, we deem that reviewing the main steps from [22] is beneficial for the understanding of the extension to the two and three dimensional cases in Sections 3.4 and 3.5 below, respectively.

Without loss of generality, we pick $T = \hat{I}$ and assume $-1$ to be the outflow facet. In 1D, the operator $\mathcal{P}_p$ is defined imposing orthogonality up to order $p-1$ over $\hat{I}$ and imposing that $\mathcal{P}_p v(-1) = v(-1)$; this last condition replaces the orthogonality on the outflow facet.

The proof of optimal $p$-convergence can be split in the following steps:

1. we write $\mathcal{P}_p v$ as a combination of Legendre polynomials and compute the corresponding coefficients; up to order $p-1$, the coefficients are the same as those of the expansion of $v$; the coefficient of order $p$ is given by a tail of the other Legendre coefficients;

2. we relate such a tail to the coefficients of $v$ minus its $L^2$ projection on the triangle restricted to the outflow facet;

3. we deduce the assertion using trace error estimates as in Lemma 3.1.

**Step 1: writing $\mathcal{P}_p v$ with respect to the Legendre basis and with explicit coefficients.**
Given $v$ in $L^2(\hat{I})$, we consider its expansion with respect to Legendre polynomials

$$v(x) = \sum_{j=0}^{+\infty} \mathbf{v}_j L_j(x).$$

Consider also the truncated expansion up to order $p$ given by

$$\mathcal{P}_p v(x) = \sum_{j=0}^{p} \hat{\mathbf{v}}_j L_j(x).$$

The orthogonality condition (12a) with respect to polynomials of degree $p-1$ gives

$$\mathcal{P}_p v(x) = \sum_{j=0}^{p-1} \mathbf{v}_j L_j(x) + \hat{\mathbf{v}}_p L_p(x).$$

The last coefficient is given imposing the condition $v(-1) = \mathcal{P}_p v(-1)$. More precisely, since $L_j(-1) = (-1)^j$ for all $j$ in $\mathbb{N}$, we write

$$v(-1) = \sum_{j=0}^{+\infty} (-1)^j \mathbf{v}_j, \qquad \mathcal{P}_p v(-1) = \sum_{j=0}^{p-1} (-1)^j \mathbf{v}_j + (-1)^p \hat{\mathbf{v}}_p.$$

We deduce

$$(-1)^p \hat{\mathbf{v}}_p = \sum_{j=p}^{+\infty} \mathbf{v}_j (-1)^j,$$

whence we get

$$\mathcal{P}_p v(x) = \sum_{j=0}^{p-1} \mathbf{v}_j L_j(x) + \left( \sum_{j=p}^{+\infty} (-1)^{j-p} \mathbf{v}_j \right) L_p(x). \qquad (13)$$
Given $\Pi_p^0$ the $L^2$ projection as in (9), we arrive at

$$ 
\begin{aligned}
(v - \mathcal{P}_p v)(x) &= \sum_{j=p}^{+\infty} \mathbf{v}_j L_j(x) - \left(\sum_{j=p}^{+\infty} (-1)^{j-p} \mathbf{v}_j\right) L_p(x) \\
&= \sum_{j=p+1}^{+\infty} \mathbf{v}_j L_j(x) + \mathbf{v}_p L_p(x) - \mathbf{v}_p L_p(x) - \left(\sum_{j=p+1}^{+\infty} (-1)^{j-p} \mathbf{v}_j\right) L_p(x) \tag{14} \\
&= \sum_{j=p+1}^{+\infty} \mathbf{v}_j L_j(x) - \left(\sum_{j=p+1}^{+\infty} (-1)^{j-p} \mathbf{v}_j\right) L_p(x) \\
&= (v - \Pi_p^0 v)(x) - (-1)^p \left(\sum_{j=p+1}^{+\infty} (-1)^j \mathbf{v}_j\right) L_p(x).
\end{aligned}
$$

**Step 2:** relate the tail of the Legendre coefficients to the trace of $v - \Pi_p^0 v$ on the outflow facet. Observe that

$$ (v - \Pi_p^0 v)(x) = \sum_{j=p+1}^{+\infty} \mathbf{v}_j L_j(x), \quad (v - \Pi_p^0 v)(-1) = \sum_{j=p+1}^{+\infty} \mathbf{v}_j L_j(-1) = \sum_{j=p+1}^{+\infty} (-1)^j \mathbf{v}_j. \quad (15) $$

Combining (14) and (15) gives

$$ (v - \mathcal{P}_p v)(x) = (v - \Pi_p^0 v)(x) - (-1)^p [(v - \Pi_p^0 v)(-1)] L_p(x) $$

Recall from (10) that

$$ \|L_p\|_{0,\hat{I}}^2 = \frac{2}{2p+1}. \quad (16) $$

We deduce

$$ \|v - \mathcal{P}_p v\|_{0,\hat{I}}^2 \le \|v - \Pi_p^0 v\|_{0,\hat{I}}^2 + |(v - \Pi_p^0 v)(-1)|^2 \frac{2}{2p+1}. $$

**Step 3:** deduce the desired bound using Lemma 3.1. The 1D version of Lemma 3.1 reads

$$ |(v - \Pi_p^0 v)(-1)|^2 \lesssim (p+1)^{-1} |v|_{1,\hat{I}}^2. $$

Combining the two equations above yields

$$ \|v - \mathcal{P}_p v\|_{0,\hat{I}} \lesssim \|v - \Pi_p^0 v\|_{0,\hat{I}} + (p+1)^{-1} |v|_{1,\hat{I}}. $$

The assertion follows noting that $\mathcal{P}_p$ and $\Pi_p^0$ preserve polynomials of maximum degree $p$ and standard polynomial approximation results.

## 3.4 The 2D case

We prove Theorem 3.2 in two dimensions. Compared to the 1D case in Section 3.3, the role of orthogonal basis functions is now played by tensor product Jacobi polynomials (with different weights) collapsed on the triangle via the Duffy transformation; in other words, we use Koornwinder polynomials following the construction in [4, 7]; see also [12, 17].

First, we introduce a reference 2D simplex $T$:

$$ \hat{T} := \left\{ \left( z \frac{1-y}{2}, y \right) \in \mathbb{R}^2 \mid (z,y) \in [-1,1]^2 \right\}. \quad (17) $$

The reference triangle $\hat{T}$ corresponds to the reference square $\hat{Q} := [-1,1]^2$ collapsed under the Duffy transformation

$$ (z, y) \in \hat{Q} \quad \implies \quad \left( z \frac{1-y}{2}, y \right). $$
This particular choice of $\hat{T}$ is convenient for our purposes, as we shall analyse the approximation properties of $\mathcal{P}_p$ assuming that $F_{\hat{T}}^+$ is the segment $[-1, 1] \times \{-1\}$.

We introduce the orthogonal polynomial basis (known as Koornwinder polynomial basis [17]) over $\hat{T}$ given by

$$ \Phi_{j,\ell}(x,y) := L_j \left( \frac{2x}{1-y} \right) J_{\ell}^{2j+1}(y) \left( \frac{1-y}{2} \right)^j \quad \forall j, \ell \in \mathbb{N}. \quad (18) $$

This is an orthogonal basis over $\hat{T}$. In fact, using the Duffy transformation, we have

$$ \int_{\hat{T}} \Phi_{i,j}(x,y) \Phi_{k,\ell}(x,y) dx dy = \left( \int_{-1}^{1} L_i(z) L_k(z) dz \right) \left( \int_{-1}^{1} J_j^{2i+1}(y) J_{\ell}^{2k+1}(y) \left( \frac{1-y}{2} \right)^{i+k+1} dy \right). $$

If $i \neq k$, the above inner product is zero. Thus, without loss of generality, we can assume $i = k$. Using (11), we write

$$ \begin{aligned} \int_{\hat{T}} \Phi_{i,j}(x,y) \Phi_{k,\ell}(x,y) dx dy &= \left( \int_{-1}^{1} L_i(z) L_k(z) dz \right) \left( \int_{-1}^{1} J_j^{2i+1}(y) J_{\ell}^{2k+1}(y) \left( \frac{1-y}{2} \right)^{2i+1} dy \right) \\ &= \frac{2}{2i+1} 2^{-2j-1} \frac{2^{2j+2}}{2\ell+2j+2} \delta_{j,\ell} = \frac{2}{2i+1} 2^{-2j-1} \frac{2^{2j+1}}{j+\ell+1} \delta_{j,\ell} = \frac{2}{2i+1} \frac{1}{j+\ell+1} \delta_{j,\ell}. \end{aligned} \quad (19) $$

With this at hand, we extend the results in Section 3.3 to 2D simplices. Henceforth, let $v$ denote a sufficiently regular function over $\hat{T}$.

**Step 1: writing $\mathcal{P}_p v$ with respect to the Koornwinder basis and with explicit coefficients.** We expand $v$ with respect to the Koornwinder basis in (18):

$$ v(x, y) = \sum_{j+\ell=0}^{+\infty} \mathbf{v}_{j,\ell} \Phi_{j,\ell}(x, y). \quad (20) $$

In what follows, the indices appearing under the summation symbol are always taken nonnegative; moreover, we set $\mathbf{v}_{j,\ell}$ equal to 0 if $j$ or $\ell$ are negative.

Since $\mathcal{P}_p v$ belongs to $\mathbb{P}_p(\hat{T})$, we expand it with respect to the Koornwinder basis:

$$ \mathcal{P}_p v(x, y) = \sum_{j+\ell=0}^{p} \hat{\mathbf{v}}_{j,\ell} \Phi_{j,\ell}(x, y). $$

By the orthogonality property (12a), we readily get

$$ \mathcal{P}_p v(x, y) = \sum_{j+\ell=0}^{p-1} \mathbf{v}_{j,\ell} \Phi_{j,\ell}(x, y) + \sum_{j+\ell=p}^{p} \hat{\mathbf{v}}_{j,\ell} \Phi_{j,\ell}(x, y). $$

The coefficients corresponding to the Koornwinder's polynomials of order $p$ are found using the facet orthogonality condition (12b) on the facet $F_{\hat{T}}^+$. Notably, we impose the identity

$$ \int_{F_{\hat{T}}^{+}} \left[ \sum_{j+\ell=p+1}^{\infty} \mathbf{v}_{j,\ell} \Phi_{j,\ell}(x, -1) + \sum_{j+\ell=p}^{\infty} (\mathbf{v}_{j,\ell} - \hat{\mathbf{v}}_{j,\ell}) \Phi_{j,\ell}(x, -1) \right] L_k(x) dx = 0 \quad \forall k = 0, \dots, p. $$

Using (18) and the fact that $J_{\ell}^{2i+1}(-1) = (-1)^{\ell}$, this identity can be rewritten as

$$ \begin{aligned} & \sum_{j+\ell=p} \int_{F_{\hat{T}}^{+}} (\mathbf{v}_{j,\ell} - \hat{\mathbf{v}}_{j,\ell})(-1)^{p-j} L_j(x) L_k(x) dx \\ & + \sum_{\tilde{p}=p+1}^{+\infty} \sum_{j+\ell=\tilde{p}} \int_{F_{\hat{T}}^{+}} \mathbf{v}_{j,\ell}(-1)^{\tilde{p}-j} L_j(x) L_k(x) dx = 0 \quad \forall k = 0, \dots, p. \end{aligned} $$
Further using (16) and the orthogonality property (10) of the Legendre polynomials, we end up with

$$(-1)^{p-k} \frac{2}{2k+1} [\mathbf{v}_{k,p-k} - \hat{\mathbf{v}}_{k,p-k}] + \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-k} \frac{2}{2k+1} \mathbf{v}_{k,\tilde{p}-k} = 0 \quad \forall k = 0, \dots, p.$$

In other words, the order p coefficients of $P_p v$ are given by

$$(-1)^{p-k} \hat{\mathbf{v}}_{k,p-k} = \sum_{\tilde{p}=p}^{+\infty} (-1)^{\tilde{p}-k} \mathbf{v}_{k,\tilde{p}-k} \quad \forall k = 0, \dots, p.$$

This is the 2D counterpart of (13). For completeness, we collect the above identities and arrive at

$$\mathcal{P}_p v(x, y) = \sum_{j+\ell=0}^{p-1} \mathbf{v}_{j,\ell} \Phi_{j,\ell}(x,y) + \sum_{j+\ell=p}^{+\infty} \left( \sum_{\tilde{p}=p}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j,\tilde{p}-j} \right) \Phi_{j,\ell}(x,y).$$

As in the 1D case, we deduce

$$\begin{aligned}
(v - \mathcal{P}_p v)(x, y) &= \sum_{j+\ell=p}^{+\infty} \mathbf{v}_{j,p-j} \Phi_{j,\ell}(x,y) - \sum_{j+\ell=p}^{+\infty} \left(\sum_{\tilde{p}=p}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j,\tilde{p}-j}\right) \Phi_{j,\ell}(x,y) \\
&= \sum_{j+\ell=p+1}^{+\infty} \mathbf{v}_{j,p-j} \Phi_{j,\ell}(x,y) - \sum_{j+\ell=p}^{+\infty} \left(\sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j,\tilde{p}-j}\right) \Phi_{j,\ell}(x,y) \quad (21) \\
&= (v - \Pi_p^0 v)(x, y) - \sum_{j+\ell=p}^{+\infty} \left(\sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j,\tilde{p}-j}\right) \Phi_{j,\ell}(x,y).
\end{aligned}$$

**Step 2: relate the tail of the Koornwinder coefficients to the trace of $v - \Pi_p^0 v$ on the outflow facet.** Using that $J_{\ell}^{2j+1}(-1) = (-1)^{\ell}$ for all $\ell$ and $j$ in $\mathbb{N}$, we write explicitly the trace of $v - \Pi_p^0 v$ on $F_{\hat{T}}^+$:

$$(v - \Pi_p^0 v)(x, y) = \sum_{j+\ell=p+1}^{+\infty} \mathbf{v}_{j,\ell} \Phi_{j,\ell}(x,y), \quad (v - \Pi_p^0 v)(x, -1) = \sum_{j+\ell=p+1}^{+\infty} (-1)^{\ell} \mathbf{v}_{j,\ell} L_j(x).$$

The orthogonality property (10) of the Legendre polynomials over $\hat{I} := [-1, 1]$ and (16), we write

$$\begin{aligned}
\|v - \Pi_p^0 v\|_{0, F_{\hat{T}}^+}^2 &= \sum_{j=0}^{+\infty} \|L_j\|_{0, \hat{I}}^2 \left( \sum_{\ell=\max(p+1-j, 0)}^{+\infty} (-1)^{\ell} \mathbf{v}_{j, \ell} \right)^2 \\
&= \sum_{j=0}^{+\infty} \frac{2}{2j+1} \left( \sum_{\ell=\max(p+1-j, 0)}^{+\infty} (-1)^{\ell} \mathbf{v}_{j, \ell} \right)^2.
\end{aligned}$$

Next, we focus on the second term on the right-hand side of (21). Using (19), we can write

$$\begin{aligned}
\left\| \sum_{j+\ell=p} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j,\tilde{p}-j} \right) \Phi_{j,\ell} \right\|_{0, \hat{T}}^2 &= \sum_{j+\ell=p} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j,\tilde{p}-j} \right)^2 \|\Phi_{j,\ell}\|_{0, \hat{T}}^2 \\
&= \sum_{j+\ell=p} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j,\tilde{p}-j} \right)^2 \frac{2}{2j+1} \frac{1}{j+\ell+1} = \frac{1}{p+1} \sum_{j=0}^{p} \frac{2}{2j+1} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j,\tilde{p}-j} \right)^2 \quad (22) \\
&\le \frac{1}{p+1} \sum_{j=0}^{+\infty} \frac{2}{2j+1} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j,\tilde{p}-j} \right)^2 = \frac{1}{p+1} \|v - \Pi_p^0 v\|_{0, F_{\hat{T}}^+}^2.
\end{aligned}$$

The last identity follows from the fact that the quantity in the parenthesis is squared and the only relevant fact is the alternating sign of the $(-1)^{\tilde{p}-p}$ term.
**Step 3:** deduce the desired bound using Lemma 3.1. Combining the above inequality with (21) and using the 2D version of Lemma 3.1, we deduce

$$ \|v - P_p v\|_{0, \hat{T}} \le \|v - \Pi_p^0\|_{0, \hat{T}} + (p+1)^{-\frac{1}{2}} \|v - \Pi_p^0 v\|_{0, F_{\hat{T}}^+} \lesssim \|v - \Pi_p^0\|_{0, \hat{T}} + (p+1)^{-1} |v|_{1, \hat{T}}. $$

The assertion follows noting that $P_p$ and $\Pi_p^0$ preserve polynomials of maximum degree $p$ and standard polynomial approximation results.

## 3.5 The 3D case

We prove Theorem 3.2 in three dimensions. Compared to the 1D and 2D cases of Sections 3.3 and 3.4, the role of orthogonal basis functions is now played by tensor product Jacobi polynomials (of different orders) collapsed on the tetrahedron using the 3D Duffy transformation; in other words, we use 3D Koornwinder polynomials following the construction in [7, 16, 24].

Introduce the reference cube $\hat{Q} = [-1, 1]^3$ and the reference tetrahedron $\hat{T}$ of vertices

$$ (-1, -1, -1); \quad (1, -1, -1); \quad (0, 1, -1); \quad (0, 0, 1). $$

We denote the lower face of $\hat{T}$ by $F_{\hat{T}}^+$; it coincides with the reference triangle in Section 3.4 at $z = 1$.

We consider the 3D Duffy transformation $\mathcal{D}: \hat{Q} \to \hat{T}$ that maps each $(z_1, z_2, z_3) = \mathbf{z}$ in $\hat{Q}$ into $(x_1, x_2, x_3) = \mathbf{x}$ in $\hat{T}$ as follows:

$$ x_1 = z_1 \frac{1-z_2}{2} \frac{1-z_3}{2}, \quad x_2 = z_2 \frac{1-z_3}{2}, \quad x_3 = z_3. $$

The Jacobian of the transformation is

$$ \det \begin{bmatrix} \frac{1-z_2}{2} & \frac{1-z_3}{2} & 0 & 0 \\ -\frac{1}{2}z_1 & \frac{1-z_3}{2} & \frac{1-z_3}{2} & 0 \\ -\frac{1}{2}z_1 & \frac{1-z_2}{2} & -\frac{z_2}{2} & 1 \end{bmatrix} = \frac{1-z_2}{2} \left( \frac{1-z_3}{2} \right)^2. $$

We also compute the formal inverse $\mathcal{D}^{-1}$ of the transformation $\mathcal{D}$:

$$ z_3 = x_3, \quad z_2 = x_2 \frac{2}{1-z_3} = \frac{2x_2}{1-x_3}, \quad z_1 = x_1 \frac{4}{(1-z_2)(1-z_3)} = \frac{4x_1}{1-2x_2-x_3}. $$

We introduce the following $L^2$-orthogonal basis over $\hat{T}$: given a multi-index $\mathbf{j} = (j_1, j_2, j_3)$ in $\mathbb{N}^3$,

$$ \Phi_j(\mathbf{x}) = L_{j_1} \left( \frac{4x_1}{1-2x_2-x_3} \right) J_{j_2}^{2j_1+1} \left( \frac{2x_2}{1-x_3} \right) \left( \frac{1-2x_2-x_3}{2(1-x_3)} \right)^{j_1} J_{j_3}^{2j_1+2j_2+2} (x_3) \left( \frac{1-x_3}{2} \right)^{j_1+j_2}. $$

We only check the norm of the above Koornwinder polynomials. Transforming $\hat{T}$ into $\hat{Q}$ by $\mathcal{D}^{-1}$ and recalling (11), we deduce

$$ \begin{align} & \int_{\hat{T}} \Phi_j(\mathbf{x}) \Phi_j(\mathbf{x}) d\mathbf{x} \notag \\ &= \int_{\hat{Q}} L_{j_1}(z_1)^2 J_{j_2}^{2j_1+1}(z_2)^2 \left(\frac{1-z_2}{2}\right)^{2j_2+1} J_{j_3}^{2j_1+2j_2+2}(z_3)^2 \left(\frac{1-z_3}{2}\right)^{2j_1+2j_2+2} dz \tag{23} \\ &= \int_{\hat{T}} L_{j_1}(z_1)^2 dz_1 \int_{\hat{T}} J_{j_2}^{2j_1+1}(z_2)^2 \left(\frac{1-z_2}{2}\right)^{2j_2+1} dz_2 \int_{\hat{T}} J_{j_3}^{2j_1+2j_2+2}(z_3)^2 \left(\frac{1-z_3}{2}\right)^{2j_1+2j_2+2} dz_3 \notag \\ &= \frac{2}{2j_1+1} \frac{2}{2j_1+2j_2+2} \frac{2}{2j_1+2j_2+2j_3+3} \quad \forall \mathbf{j} = (j_1, j_2, j_3) \in \mathbb{N}^3. \notag \end{align} $$

Next, we prove $p$-optimal approximation properties in the $L^2$ norm of the operator $P_p$ given as in (12) with the role of outflow face played by the facet $F_{\hat{T}}^+$.
**Step 1:** writing $\mathcal{P}_p v$ with respect to the 3D version Koornwinder basis and with explicit coefficients. Consider the following expansion of $v$:

$$v(\mathbf{x}) = \sum_{|j|=0}^{+\infty} \mathbf{v}_j \Phi_j(\mathbf{x}).$$

In what follows, the indices appearing under the summation symbol are always taken nonnegative; $|j|$ stands for $j_1 + j_2 + j_3$; $\mathbf{v}_j$ is set to 0 if any of the indices $j_1, j_2$, and $j_3$ is negative.

Using (12a), we realise that

$$v(\mathbf{x}) = \sum_{|j|=0}^{p-1} \mathbf{v}_j \Phi_j(\mathbf{x}) + \sum_{|j|=p} \hat{\mathbf{v}}_j \Phi_j(\mathbf{x}),$$

where the coefficients $\hat{\mathbf{v}}_j$, $|j| = p$ are to be determined using (12b).

First, we write the trace of any polynomial basis function on $F_{\hat{T}}^{+}$, i.e., we impose the passage through $x_3 = -1$:

$$\Phi_j(\mathbf{x})|_{F_{\hat{T}}^{+}} = L_{j_1} \left(\frac{2x_1}{1-x_2}\right) J_{j_2}^{2j_1+1}(x_2) \left(\frac{1-x_2}{2}\right)^{j_1} (-1)^{j_3}.$$

Let $\Psi_{\ell_1, \ell_2}$ be the 2D Koornwinder basis in (18) on the reference triangle $F_{\hat{T}}^{+}$. We impose (12b):

$$\sum_{|j|=p} \int_{F_{\hat{T}}^{+}} [\mathbf{v}_j - \hat{\mathbf{v}}_j] (-1)^{j_3} L_{j_1} \left(\frac{2x_1}{1-x_2}\right) L_{j_2}(x_2) \left(\frac{1-x_2}{2}\right)^{j_1} \Psi_{\ell_1, \ell_2}(x_1, x_2) dx_1 dx_2 \\ + \sum_{\tilde{p}=p+1}^{+\infty} \sum_{|j|=\tilde{p}} \int_{F_{\hat{T}}^{+}} \mathbf{v}_j (-1)^{j_3} L_{j_1} \left(\frac{2x_1}{1-x_2}\right) L_{j_2}(x_2) \left(\frac{1-x_2}{2}\right)^{j_1} \Psi_{\ell_1, \ell_2}(x_1, x_2) dx_1 dx_2 = 0 \\ \forall \ell_1, \ell_2 \text{ such that } \ell_1 + \ell_2 \le p.$$

We rewrite the above relation as

$$\sum_{|j|=p} \int_{F_{\hat{T}}^{+}} [\mathbf{v}_j - \hat{\mathbf{v}}_j] (-1)^{j_3} \Psi_{j_1, j_2}(x_1, x_2) \Psi_{\ell_1, \ell_2}(x_1, x_2) dx_1 dx_2 \\ + \sum_{\tilde{p}=p+1}^{+\infty} \sum_{|j|=\tilde{p}} \int_{F_{\hat{T}}^{+}} \mathbf{v}_j (-1)^{j_3} \Psi_{j_1, j_2}(x_1, x_2) \Psi_{\ell_1, \ell_2}(x_1, x_2) dx_1 dx_2 = 0 \quad \forall \ell_1, \ell_2 \text{ such that } \ell_1 + \ell_2 \le p.$$

Using the orthogonality property (19) of the Koornwinder polynomials in 2D, we deduce

$$(-1)^{p-\ell_1-\ell_2} \frac{2}{2\ell_1+1} \frac{1}{\ell_1+\ell_2+1} (\mathbf{v}_{\ell_1, \ell_2, p-\ell_1-\ell_2} - \hat{\mathbf{v}}_{\ell_1, \ell_2, p-\ell_1-\ell_2}) \\ + \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-(\ell_1+\ell_2)} \frac{2}{2\tilde{p}+1} \frac{1}{\tilde{p}+\tilde{p}+1} \mathbf{v}_{\ell_1, \ell_2, \tilde{p}-\ell_1-\ell_2} = 0 \quad \forall \ell_1, \ell_2 \text{ such that } \ell_1 + \ell_2 \le p.$$

In other words, we have

$$\hat{\mathbf{v}}_{\ell_1, \ell_2, p-\ell_1-\ell_2} = \sum_{\tilde{p}=p}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{\ell_1, \ell_2, \tilde{p}-\ell_1-\ell_2} \quad \forall \ell_1, \ell_2 \text{ such that } \ell_1 + \ell_2 \le p.$$

This provides us with the representation

$$\mathcal{P}_p v(\mathbf{x}) = \sum_{|j|=0}^{p-1} \mathbf{v}_j \Phi_j(\mathbf{x}) + \sum_{|j|=p} \left( \sum_{\tilde{p}=p}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j_1, j_2, \tilde{p}-j_1-j_2} \right) \Phi_j(\mathbf{x}).$$
We get

$$ 
\begin{aligned}
(\mathbf{v} - \mathcal{P}_p \mathbf{v})(\mathbf{x}) &= \sum_{|j|=p}^{+\infty} \mathbf{v}_{j_1,j_2,p-j_1-j_2} \Phi_j(\mathbf{x}) - \sum_{|j|=p}^{+\infty} \left( \sum_{\tilde{p}=p}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j_1,j_2,\tilde{p}-j_1-j_2} \right) \Phi_j(\mathbf{x}) \\
&= \sum_{|j|=p+1}^{+\infty} \mathbf{v}_{j_1,j_2,p-j_1-j_2} \Phi_j(\mathbf{x}) - \sum_{|j|=p}^{+\infty} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j_1,j_2,\tilde{p}-j_1-j_2} \right) \Phi_j(\mathbf{x}) \quad (24) \\
&= (\mathbf{v} - \Pi_p^0 \mathbf{v})(\mathbf{x}) - \sum_{|j|=p}^{+\infty} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j_1,j_2,\tilde{p}-j_1-j_2} \right) \Phi_j(\mathbf{x}).
\end{aligned}
$$

**Step 2:** relate the tail of the 3D version Koornwinder coefficients to the trace of $v - \Pi_p^0 v$ on the outflow facet. The restriction of $\mathbf{v} - \Pi_p^0 \mathbf{v}$ on the outflow facet $F_{\hat{T}}^+$ reads

$$ 
\begin{aligned}
(\mathbf{v} - \mathcal{P}_p \mathbf{v})(\mathbf{x})|_{F_{\hat{T}}^+} &= \sum_{|j|=p+1}^{+\infty} \mathbf{v}_j L_{j_1} \left( \frac{2x_1}{1-x_2} \right) J_{j_2}^{2j_1+1} (x_2) \left( \frac{1-x_2}{2} \right)^{j_1} (-1)^{j_3} \\
&= \sum_{|j|=p+1}^{+\infty} \mathbf{v}_j \Psi_{j_1,j_2}(\mathbf{x}) (-1)^{j_3}.
\end{aligned}
$$

We take the $L^2$ norm on both sides over $F_{\hat{T}}^+$, use the orthogonality property (19) of Koornwinder polynomials over $F_{\hat{T}}^+$, and get, for $\tilde{j}$ equal to $\max(p+1 - (j_1 + j_2), 0)$,

$$ \|\mathbf{v} - \mathcal{P}_p \mathbf{v}\|_{0, F_{\hat{T}}^+} = \sum_{j_1+j_2=0}^{+\infty} \|\Psi_{j_1,j_2}\|_{0, F_{\hat{T}}^+}^2 (\sum_{j_3=\tilde{j}} (-1)^{j_3} \mathbf{v}_j)^2 = \sum_{j_1+j_2=0}^{+\infty} \frac{2}{2j_2+1} \frac{1}{j_1+j_2+1} (\sum_{j_3=\tilde{j}} (-1)^{j_3} \mathbf{v}_j)^2. $$

We show a bound on the $L^2$ norm over $\hat{T}$ of the second term on the right-hand side of (24) in terms of the $L^2$ norm of $\mathbf{v} - \mathcal{P}_p \mathbf{v}$ on $F_{\hat{T}}^+$. To this aim, we use (23):

$$ 
\begin{aligned}
\left\| \sum_{|j|=p} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j_1,j_2,\tilde{p}-j_1-j_2} \right) \Phi_j(\mathbf{x}) \right\|_{0, \hat{T}}^2 &= \sum_{|j|=p} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j_1,j_2,\tilde{p}-j_1-j_2} \right)^2 \|\Phi_j\|_{0, \hat{T}}^2 \\
&= \sum_{|j|=p} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j_1,j_2,\tilde{p}-j_1-j_2} \right)^2 \frac{2}{2j_1+1} \frac{1}{j_1+j_2+1} \frac{1}{j_1+j_2+j_3+1} \\
&= \frac{1}{p+1} \sum_{j_1+j_2=0}^{p} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j_1,j_2,\tilde{p}-j_1-j_2} \right)^2 \frac{2}{2j_1+1} \frac{1}{j_1+j_2+1} \\
&\le \frac{1}{p+1} \sum_{j_1+j_2=0}^{+\infty} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j_1,j_2,\tilde{p}-j_1-j_2} \right)^2 \frac{2}{2j_1+1} \frac{1}{j_1+j_2+1} = \frac{1}{p+1} \|\mathbf{v} - \Pi_p^0 \mathbf{v}\|_{0, F_{\hat{T}}^+}^2.
\end{aligned}
$$

**Step 3:** deduce the desired bound using optimal trace error estimates. Combining the above inequality with (21) and using the 3D version of Lemma 3.1, we deduce

$$ \|v - \mathcal{P}_p v\|_{0, \hat{T}} \le \|v - \Pi_p^0\|_{0, \hat{T}} + (p+1)^{-\frac{1}{2}} \|v - \Pi_p^0 v\|_{0, F_{\hat{T}}^+} \lesssim \|v - \Pi_p^0\|_{0, \hat{T}} + (p+1)^{-1} |v|_{1, \hat{T}}. $$

The assertion follows noting that $\mathcal{P}_p$ and $\Pi_p^0$ preserve polynomials of maximum degree $p$ and standard polynomial approximation results.

## 3.6 Trace-type estimates for the CDG operator

The CDG operator satisfies hp-optimal approximation properties in the $L^2$ norm; see Theorem 3.2. Here, we investigate hp-optimal approximation properties in the $L^2$ norm on the boundary of a simplex T, i.e., the counterpart of Lemma 3.1 for the projector $\mathcal{P}_p$.
To this aim, we show separate bounds on the facets. First, we consider the case of the (unique) outflow facet $F_T^+$. By (12b) and Lemma 3.1, we get

$$ \|v - \mathcal{P}_p v\|_{0, F_T^+} \le \|v - \Pi_p^0 v\|_{0, F_T^+} \lesssim \left(\frac{h_T}{p}\right)^{\frac{1}{2}} |v|_{1,T} \quad \forall v \in H^1(T). $$

Approximation properties on the inflow facets are more elaborated and are discussed in the next result.

**Proposition 3.3.** Let the assumptions (A1) and (A2) be valid, $T$ be in $T_h$, $F_T^-$ be any of the inflow facets of $T$, and $\Pi_p^0$ be the operator in (9). Then, there exist positive constants $C_1$ and $C_2$ independent of $h$ and $p$ such that, for every $v$ in $H^{k+1}(T)$, $k$ positive, the following estimate holds true:

$$ \|v - \mathcal{P}_p v\|_{0, F_T^-} \le C_1 \|v - \Pi_p^0 v\|_{0, \partial T} \le C_2 \frac{h^{\min(p,k)+\frac{1}{2}}}{p^{k+\frac{1}{2}}} \|v\|_{k+1, T}. \quad (25) $$

*Proof.* The proof of the 1D case is fairly simple. The two and three dimensional cases are more elaborated; we provide details for the 2D case, as the 3D case can be dealt with similarly. Throughout, we employ the notation of Sections 3.2–3.5. We prove the assertion on the reference simplex $\hat{T}$; the general statement follows from a scaling argument.

**Proof of (25) in 1D.** Here $\hat{T} = [-1, 1]$. From (14), (15), and $L_j(1) = 1$ for all $j$ in $\mathbb{N}$, we deduce

$$ |(v - \mathcal{P}_p v)(1)| \le |(v - \Pi_p^0 v)(-1)| + |(v - \Pi_p^0 v)(1)|. $$

The assertion follows from Lemma 3.1.

**Proof of (25) in 2D.** Here, $\hat{T}$ is as in (17). Recall expansion (20) with the caveat that coefficients related to negative indices are set to zero.

Introduce

$$ \mathcal{R} := \sum_{j+\ell=p} \left( \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j, \tilde{p}-j} \right) \Phi_{j, \ell}(x, y) $$

and

$$ \hat{C}_p := \max_{j+\ell=p} \|\Phi_{j,\ell}(x,y)\|_{0, F_T^-}. $$

We have

$$ \|\mathcal{R}\|_{0, F_T^-}^2 \lesssim \sum_{j+\ell=p} \left| \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j, \tilde{p}-j} \right|^2 \|\Phi_{j,\ell}(x,y)\|_{0, F_T^-}^2 \le \hat{C}_p^2 \sum_{j+\ell=p} \left| \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j, \tilde{p}-j} \right|^2, $$

whence we deduce

$$ \begin{aligned} \|\mathcal{R}\|_{0, F_T^-}^2 &\lesssim \hat{C}_p^2 \sum_{j+\ell=p} \left| \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j, \tilde{p}-j} \right|^2 \frac{2}{2j+1} \frac{2j+1}{2} \\ &\le \hat{C}_p^2 \frac{2p+1}{2} \sum_{j+\ell=p} \left| \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j, \tilde{p}-j} \right|^2 \frac{2}{2j+1}. \end{aligned} $$

We proceed as in (22) and arrive at

$$ \|\mathcal{R}\|_{0, F_T^-}^2 \lesssim \hat{C}_p^2 \frac{2p+1}{2} \sum_{j+\ell=p} \left| \sum_{\tilde{p}=p+1}^{+\infty} (-1)^{\tilde{p}-p} \mathbf{v}_{j, \tilde{p}-j} \right|^2 \frac{2}{2j+1} \le \hat{C}_p^2 \frac{2p+1}{2} \|v - \Pi_p^0 v\|_{0, F_T^+}^2. $$

Combining (21) with the inequality above yields

$$ \|v - \mathcal{P}_p v\|_{0, F_T^-} \le \|v - \Pi_p^0 v\|_{0, F_T^-} + \|\mathcal{R}\|_{0, F_T^-} \lesssim \|v - \Pi_p^0 v\|_{0, F_T^-} + \hat{C}_p p^{\frac{1}{2}} \|v - \Pi_p^0 v\|_{0, F_T^+}. \quad (26) $$
If we were able to prove that $\hat{C}_p \lesssim p^{-\frac{1}{2}}$, then the assertion would follow from Lemma 3.1 and (26). In fact, we would write

$$ \|v - \mathcal{P}_p v\|_{0, F_{\hat{T}}^{-}} \le \|v - \Pi_p^0 v\|_{0, F_{\hat{T}}^{-}} + \|\mathcal{R}\|_{0, F_{\hat{T}}^{-}} \lesssim \|v - \Pi_p^0 v\|_{0, F_{\hat{T}}^{-}} + \|v - \Pi_p^0 v\|_{0, F_{\hat{T}}^{+}}. $$

We are left with showing $\hat{C}_p \lesssim p^{-\frac{1}{2}}$. For the sake of presentation, we fix $F_{\hat{T}}^{-}$ to be the inflow facet of $\hat{T}$ with negative $x$ coordinates. If this is the case, then a direct application of definition (18) gives

$$ \Phi_{j, \ell|F_{\hat{T}}^{-}} = (-1)^j J_{\ell}^{2j+1}(y) \left(\frac{1-y}{2}\right)^j. $$

Taking into account the Jacobian of the transformation mapping the facet $F_{\hat{T}}^{-}$ into the “vertical” facet $\{-1\} \times [-1, 1]$ and recalling (11), we deduce

$$ \|\Phi_{j, \ell}\|_{0, F_{\hat{T}}^{-}}^2 = \int_{-1}^{1} \left[ J_{\ell}^{2j+1}(y) \right]^2 \left( \frac{1-y}{2} \right)^{2j+1} = \frac{2^{2j+2}}{2j+1+2\ell+1} \frac{1}{2^{2j+1}} = \frac{1}{1+j+\ell}. $$

Since we are considering couples $(j, \ell)$ such that $j + \ell = p$, we arrive at

$$ \|\Phi_{j, \ell}\|_{0, F_{\hat{T}}^{-}}^2 = \frac{1}{p+1}, $$

i.e., $\hat{C}_p \lesssim p^{-\frac{1}{2}}$.

□

# 4 A priori analysis

We prove hp-optimal a priori error estimates for method (6). The proof is essentially the hp-version of that of [8, Theorem 2.2]. We report here the details for three reasons: for the sake of completeness; since we are interested in hp-optimal a priori bounds; because we have different assumptions, namely we require (5) for well posedness of the method, but we have no restrictions on the mesh size $h$ as for instance required in [15].

**Theorem 4.1.** Let $T_h$ be a shape-regular, quasi uniform mesh satisfying the assumptions (A1)-(A2), and $u$ and $u_h$ be the solutions to (4) and (6). If (5) holds true and $u$ belongs to $H^{k+1}(\Omega)$, $k$ positive, then there exists a positive constant $C$ independent of $h$ and $p$ but possibly depending on $\sigma$ such that

$$ \|u - u_h\|_{0, \Omega} \le C \frac{h^{\min(k,p)+1}}{p^{k+1}} \|u\|_{k+1, \Omega}. $$

*Proof.* Introduce

$$ e_h := u_h - \mathcal{P}_p u. $$

The rewriting of the bilinear form in (7) and the Galerkin orthogonality (8) imply

$$ \mathcal{B}(e_h, v_h) = \mathcal{B}(u - \mathcal{P}_p u, v_h) = \sum_{j=1}^{3} T_j, $$

where

$$ T_1 := -(u - \mathcal{P}_p u, \beta \cdot \nabla_h v_h)_{0, \Omega}, $$

$$ T_2 := \sum_{T \in \mathcal{T}_h} \left( (u - \mathcal{P}_p u)^-, \mathbf{n}_T \cdot \beta \left[ v_h \right] \right)_{0, \Gamma_{\text{in}}^T} + \sum_{T \in \mathcal{T}_h} (u - \mathcal{P}_p u, \mathbf{n}_T \cdot \beta v_h)_{0, \partial T \cap \Gamma_+}, $$

$$ T_3 := (c(u - \mathcal{P}_p u), v_h)_{0, \Omega}. $$

The term $T_1$ vanishes due to (12a).

As for the term $T_2$, some manipulations imply

$$ T_2 = - \sum_{T \in \mathcal{T}_h} (u - \mathcal{P}_p u, \mathbf{n}_T \cdot \beta [v_h])_{0, \Gamma_{\text{out}}^T} + \sum_{T \in \mathcal{T}_h} (u - \mathcal{P}_p u, \mathbf{n}_T \cdot \beta v_h)_{0, \partial T \cap \Gamma_+}. $$
The last identity is a consequence of the definition of outflow facets. Recalling (12b), also $T_2$ vanishes. We end up with

$$ \mathcal{B}(e_h, v_h) = (c(u - \mathcal{P}_p u), v_h)_{0,\Omega}. $$

This means that $e_h$ solves method (6) with data $f$ and $g$ given by $c(u - \mathcal{P}_p u)$ and $0$, respectively. This identity and the stability estimate in Proposition 2.1 imply

$$ \|e_h\|_{0,\Omega} \lesssim \|c(u - \mathcal{P}_p u)\|_{0,\Omega}. $$

The triangle inequality and the above bound entail

$$ \|u - u_h\|_{0,\Omega} \le \|u - \mathcal{P}_p u\|_{0,\Omega} + \|u - u_h\|_{0,\Omega} \lesssim \|u - \mathcal{P}_p u\|_{0,\Omega} + \|c(u - \mathcal{P}_p u)\|_{0,\Omega}. $$

The assertion follows using Theorem 3.2. $\square$

We also prove hp-optimal convergence for method (6) in the full DG norm

$$ \|v_h\|_{\text{DG}}^2 := \sum_{T \in \mathcal{T}_h} \left( \|c v_h\|_{0,T}^2 + \|\mathbf{n}_T \cdot \boldsymbol{\beta}|^{\frac{1}{2}} v_h\|_{0,\Gamma_{\text{in}}^T \cap \Gamma_-}^2 + \|\mathbf{n}_T \cdot \boldsymbol{\beta}|^{\frac{1}{2}} v_h^+\|_{0,\Gamma_{\text{out}}^T \cap \Gamma_+}^2 + \|\mathbf{n}_T \cdot \boldsymbol{\beta}|^{\frac{1}{2}} \llbracket v_h \rrbracket\|_{0,\Gamma_{\text{in}}^T \setminus \Gamma_-}^2 \right). $$

**Theorem 4.2.** Let $\mathcal{T}_h$ be a shape-regular, quasi uniform mesh satisfying the assumptions (A1)-(A2), and $u$ and $u_h$ be the solutions to (4) and (6). If (5) holds true and $u$ belongs to $H^{k+1}(\Omega)$, $k$ positive, then there exists a positive constant $C$ independent of $h$ and $p$ but possibly depending on $\sigma$ such that

$$ \|u - u_h\|_{\text{DG}} \le C \frac{h^{\min(k,p)+\frac{1}{2}}}{p^{k+\frac{1}{2}}} \|u\|_{k+1,\Omega}. $$

*Proof.* This is the simplicial version of [13, Theorem 3.7]. As in the proof of that result, given $\eta := u - \Pi_p^0 u$, it is possible to show that

$$ \begin{aligned} \|u - u_h\|_{\text{DG}} & \lesssim (\sum_{T \in \mathcal{T}_h} \|\eta\|_{0,T}^2)^{\frac{1}{2}} + \left( \sum_{T \in \mathcal{T}_h} \left( \|\mathbf{n}_T \cdot \boldsymbol{\beta}|^{\frac{1}{2}} \eta^+\|_{0,\Gamma_{\text{out}}^T \cap \Gamma_+}^2 + \|\mathbf{n}_T \cdot \boldsymbol{\beta}|^{\frac{1}{2}} \eta^+\|_{0,\Gamma_{\text{in}}^T \cap \Gamma_-}^2 \right. \right. \\ & \qquad \left. + \|\mathbf{n}_T \cdot \boldsymbol{\beta}|^{\frac{1}{2}} \eta^-\|_{0,\Gamma_{\text{out}}^T \cap \Gamma_-}^2 + \|\mathbf{n}_T \cdot \boldsymbol{\beta}|^{\frac{1}{2}} \eta^+\|_{0,\Gamma_{\text{out}}^T \cap \Gamma_-}^2 \right) \right)^{\frac{1}{2}} =: A + B. \end{aligned} $$

The hidden constant above is independent of $h$ and $p$.

The term $A$ is dealt with using standard hp Babuška-Suri type approximation properties of the $L^2$ projection in the $L^2$ norm; see, e.g., [1, Lemma 4.5]. The term $B$ is dealt with using approximation properties of the $L^2$ projection in the boundary $L^2$ norm as discussed in Lemma 3.1. $\square$

Theorems 4.1 and 4.2 improve the current state of the art of the literature along different directions:

* [8, Theorem 2.2] shows h-optimal convergence in the $L^2$ norm; here, we provide full explicit track of the p-convergence;

* [13, Theorem 3.7] and [14, Theorem 5.1] display hp-optimal convergence in the full DG norm on Cartesian-type meshes; here, we show optimal hp-optimal convergence in the full DG norm on special simplicial meshes and also improved convergence (on simplicial meshes) in the $L^2$ norm.

The theoretical limitations of Theorems 4.1 and 4.2 are the use

* of the special meshes in Section 2.1 from [8];

* of assumption (5).

In Section 5 below, we investigate whether the two limitations above can be overcome in practice or are necessary condition for the hp-optimality.
Remark 2. Following the proof of Theorem 4.2, it is possible to derive hp-optimal convergence of DG methods on simplicial meshes for diffusion-advection-reaction problems as in [14]. In words, one needs to care also of an interior penalty discontinuous discretization of a diffusion term; the error of the method is then estimated from above by the distance in a DG type norm between the exact solution and its L² projection. Such a quantity is then estimated optimally using the Babuška-Suri type approximation properties of the L² projection in the L² norm (using [1]), in the H¹ seminorm (using [5]), and in the L² norm on the boundary (using Lemma 3.1).

# 5 Numerical experiments

We assess the convergence rates of method (6) in the L² and full DG norms proven in Theorems 4.1 and 4.2. We focus on the p-version of the method in two dimensions, fix two triangular meshes of 50 and 32 triangles partitioning the square domain $\Omega := (-1, 1)^2$, see Figure 1, and use polynomial degrees p from 1 to 40. Standard composite Gaussian quadrature is employed throughout. We pick an L² orthonormal basis over simplices as that constructed in [16].

We consider two test cases. The first one involves a constant convection field such that assumption (5) holds true and the meshes satisfy the construction in Section 2.1. Instead, a second test case is devoted to check the performance of the scheme even if the convection field and the meshes do not satisfy the assumptions needed in the analysis of Sections 3 and 4.

<figure><img src="image_2.png" /><figcaption>Figure 1: Left-panel: mesh consisting of 50 triangles. Right-panel: mesh consisting of 32 triangles.</figcaption></figure>

**Test case 1.** As in [6, Section 5.2], we consider the convection field $\beta = (1, 1)$, the reaction coefficient $c = 1$, and the exact solution, for given positive $\alpha$,

$$u(x, y) = \begin{cases} \cos(\pi y/2) & \text{in } (-1, 0] \times (-1, 1) \\ \cos(\pi y/2) + x^{\alpha} & \text{in } (0, 1) \times (-1, 1). \end{cases} \quad (27)$$

The right-hand side $f$ is computed accordingly. The two meshes in Figure 1 are such that each element has only one outflow facet for the given convection field $\beta$.

The solution $u$ belongs to $H^{1/2+\alpha-\epsilon}(\Omega)$ for any positive and arbitrarily small $\epsilon$. The singularity is located in the interior of a mesh cell for the mesh in Figure 1 (left-panel) and at the interface of several cell elements for the mesh in Figure 1 (right-panel).

We report the results in Figure 2 and 3 for the 50 and 32 triangles cases, respectively. We pick $\alpha = 0.5, 1.5$, and $2.5$.

The convergence rate of the DG scheme under $p$ refinement in the L²-norm error is of order $O(p^{-(1/2+\alpha)})$, which is optimal. The convergence rate in the DG norm is also of order $O(p^{-(1/2+\alpha)})$, i.e., we observe convergence rates better than those proven in theory.

For the case of 32 triangular elements, the singularity lies at the interface of several mesh cell elements. Therefore, the usual doubling of the convergence takes place, see, e.g., [2], and the convergence rate is $O(p^{-2\alpha})$ in both norms.
<figure><img src="image_3.png" /><figcaption>**Figure 2:** *p*-version for the case of constant convection field $\beta$. We consider the exact solution $u$ in (27) and fix the triangular mesh with 50 elements as in Figure 1 so as the mesh assumptions in Section 2.1 are satisfied. The singularity lies in the interior of a cell element. Left-panel: $L^2$ norm error. Right-panel: DG norm error.</figcaption></figure>

**Test case 2.** As in [14, Example 1], we consider the convection field $\beta = (2 - y^2, 2 - x)$, the reaction coefficient $c = 1+(1+x)(1+y^2)$, and the exact solution given by $u$ in (27). The right-hand side $f$ is computed accordingly. The two meshes in Figure 1 do not satisfy the mesh assumptions in Section 2.1 for the given convection field $\beta$. The regularity of the exact solution $u$ is already discussed in Test case 1.

We report the results in Figure 4 and 5 for the 50 and 32 triangles cases, respectively. We pick $\alpha = 0.5, 1.5$, and $2.5$.

We observe convergence rates as those for the constant convection field case. This is in accordance with the numerical observation in [14, Example 1].

In both test cases, the use of an orthonormal basis implies that the condition number remains moderate, whence the convergence rates do not deteriorate for high $p$.

# 6 Conclusions

We analyzed $p$-optimal approximation properties of the CDG operator in the $L^2$ norm and the $L^2$ norm of the trace on simplices with only one outflow facet in 1D, 2D, and 3D. These results are instrumental in deriving $p$-optimal error estimates for the original DG method on special simplicial meshes. Numerical experiments validate the predicted convergence rates measured in the $L^2$ norm also for convection fields and meshes not satisfying the assumptions used in our proofs; convergence rates are observed for the error measured in the full DG norm that are better than those proven in theory.

**Acknowledgements.** We are very grateful to the anonymous referees who helped us to improve the presentation and readability of the paper. LM is member of the Gruppo Nazionale Calcolo Scientifico-Istituto Nazionale di Alta Matematica (GNCS-INdAM) and was partially supported by the Italian MIUR through the PRIN Grant No. 2022-NAZ-0279.

# References

[1] I. Babuška and M. Suri. The hp version of the finite element method with quasiuniform meshes. *ESAIM Math. Model. Numer. Anal.*, 21(2):199–238, 1987.

[2] I. Babuška and M. Suri. The optimal convergence rate of the $p$-version of the finite element method. *SIAM J. Numer. Anal.*, 24(4):750–776, 1987.
<figure><img src="image_4.png" /><figcaption>**Figure 3:** *p*-version for the case of constant convection field $\beta$. We consider the exact solution $u$ in (27) and fix the triangular mesh with 32 elements as in Figure 1 so as the mesh assumptions in Section 2.1 are satisfied. The singularity lies at the interface of several cell elements. Left-panel: $L^2$ norm error. Right-panel: DG norm error.</figcaption></figure>

[3] C. Bardos. Problèmes aux limites pour les équations aux dérivées partielles du premier ordre à coefficients réels; théorèmes d'approximation; application à l'équation de transport. *Ann. Sci. de l'Ecole Norm. Superieure*, 3(2):185-233, 1970.

[4] S. Beuchler and J. Schöberl. New shape functions for triangular *p*-FEM using integrated Jacobi polynomials. *Numer. Math.*, 103(3):339-366, 2006.

[5] C. Canuto and A. Quarteroni. Approximation results for orthogonal polynomials in Sobolev spaces. *Math. Comp.*, 38(157):67-86, 1982.

[6] P. Castillo, B. Cockburn, I. Perugia, and D. Schötzau. An a priori error analysis of the local discontinuous Galerkin method for elliptic problems. *SIAM J. Numer. Anal.*, 38(5):1676-1706, 2000.

[7] A. Chernov. Optimal convergence estimates for the trace of the polynomial $L^2$-projection operator on a simplex. *Math. Comp.*, 81(278):765-787, 2012.

[8] B. Cockburn, B. Dong, and J. Guzmán. Optimal convergence of the original DG method for the transport-reaction equation on special meshes. *SIAM J. Numer. Anal.*, 46(3):1250-1265, 2008.

[9] B. Cockburn, B. Dong, and J. Guzmán. A superconvergent LDG-hybridizable Galerkin method for second-order elliptic problems. *Math. Comp.*, 77(264):1887-1916, 2008.

[10] B. Cockburn, B. Dong, J. Guzmán, and J. Qian. Optimal convergence of the original DG method on special meshes for variable transport velocity. *SIAM J. Numer. Anal.*, 48(1):133-146, 2010.

[11] Z. Dong and L. Mascotto. On the suboptimality of the *p*-version discontinuous Galerkin methods for first order hyperbolic problems. *Numerical Methods and Algorithms in Science and Engineering*, 700, 2021.

[12] M. Dubiner. Spectral methods on triangles and other domains. *J. Sci. Comput.*, 6:345-390, 1991.

[13] P. Houston, Ch. Schwab, and E. Süli. Stabilized *hp*-finite element methods for first-order hyperbolic problems. *SIAM J. Numer. Anal.*, 37(5):1618-1643, 2000.

[14] P. Houston, Ch. Schwab, and E. Süli. Discontinuous *hp*-finite element methods for advection-diffusion-reaction problems. *SIAM J. Numer. Anal.*, 39(6):2133-2163, 2002.

[15] C. Johnson and J. Pitkäranta. An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. *Math. Comp.*, 46(173):1-26, 1986.

[16] G. E. Karniadakis and S. Sherwin. *Spectral/hp element methods for computational fluid dynamics*. Oxford University Press, 2005.

[17] T. Koornwinder. Two-variable analogues of the classical orthogonal polynomials. In *Theory and application of special functions*, pages 435-495. Elsevier, 1975.

[18] P. Lasaint and P.-A. Raviart. On a finite element method for solving the neutron transport equation. In *Mathematical aspects of finite elements in partial differential equations (Proc. Sympos., Math. Res. Center, Univ. Wisconsin, Madison, Wis., 1974)*, pages 89-123. Academic Press, New York-London, 1974.

[19] J. M. Melenk and T. Wurzer. On the stability of the boundary trace of the polynomial $L^2$-projection on triangles and tetrahedra. *Comp. Math. Appl. Sc.*, 67(4):944-965, 2014.

[20] T. E. Peterson. A note on the convergence of the discontinuous Galerkin method for a scalar hyperbolic equation. *SIAM J. Numer. Anal.*, 28(1):133-140, 1991.
<figure><img src="image_5.png" /><figcaption>**Figure 4:** *p*-version for the case of variable convection field $\beta$. We consider the exact solution $u$ in (27) and fix the triangular mesh with 50 elements as in Figure 1. The mesh assumptions in Section 2.1 are not satisfied. The singularity lies in the interior of a cell element. Left-panel: $L^2$ norm error. Right-panel: DG norm error.</figcaption></figure>

[21] W. H. Reed and T. R. Hill. Triangular mesh methods for the neutron transport equation. Technical report, Los Alamos Scientific Lab., 1973.

[22] D. Schötzau and Ch. Schwab. Time discretization of parabolic problems by the hp-version of the discontinuous Galerkin finite element method. *SIAM J. Numer. Anal.*, 38(3):837–875, 2000.

[23] J. Shen, T. Tang, and L.-L. Wang. *Spectral methods: algorithms, analysis and applications*, volume 41. Springer Science & Business Media, 2011.

[24] T. Warburton and J. S. Hesthaven. On the constants in hp-finite element trace inverse inequalities. *Comput. Methods Appl. Mech. Engrg.*, 192(25):2765–2773, 2003.
<figure><img src="image_6.png" /><figcaption>**Figure 5:** *p*-version for the case of variable convection field $\beta$. We consider the exact solution $u$ in (27) and fix the triangular mesh with 32 elements as in Figure 1. The mesh assumptions in Section 2.1 are not satisfied. The singularity lies at the interface of several cell elements. Left-panel: $L^2$ norm error. Right-panel: DG norm error.</figcaption></figure>