# hp-optimal interior penalty discontinuous Galerkin methods for the biharmonic problem

Zhaonan Dong, Lorenzo Mascotto

► To cite this version:

Zhaonan Dong, Lorenzo Mascotto. hp-optimal interior penalty discontinuous Galerkin methods for the biharmonic problem. Journal of Scientific Computing, In press, 96 (30), pp.30. 10.1007/s10915-023-02253-y. hal-03889072v2

HAL Id: hal-03889072

https://inria.hal.science/hal-03889072v2

Submitted on 3 May 2023

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 interior penalty discontinuous Galerkin methods
for the biharmonic problem

Z. Dong*, L. Mascotto†‡§

Abstract

We prove *hp*-optimal error estimates for interior penalty discontinuous Galerkin methods (IPDG) for the biharmonic problem with homogeneous essential boundary conditions. We consider tensor product-type meshes in two and three dimensions, and triangular meshes in two dimensions. An essential ingredient in the analysis is the construction of a global *H*² piecewise polynomial approximants with *hp*-optimal approximation properties over the given meshes. The *hp*-optimality is also discussed for *C*⁰-IPDG in two and three dimensions, and the stream formulation of the Stokes problem in two dimensions. Numerical experiments validate the theoretical predictions and reveal that *p*-suboptimality occurs in presence of singular essential boundary conditions.

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

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

# 1 Introduction

The numerical approximation of solutions to the biharmonic problem is a challenging task tracing back to the 60ties [3]; see also the later work [16]. The main difficulty in designing standard Galerkin methods for this problem resides in the construction of global $H^2$ conforming polynomial spaces over the given meshes. A possible way to circumvent this issue is based on using nonconforming finite element spaces. The Morley element is the first and one of the most popular among them; see, e.g., [28, 38].

Another option to avoid $H^2$ conformity is to use $H^1$ conforming piecewise polynomial spaces. This approach goes under the name of $C^0$-interior penalty discontinuous Galerkin method ($C^0$-IPDG) and has been a very active area of research over the last 20 years; see, e.g., [7, 9, 11, 17]. The continuity constraint can be further weakened leading to the use of fully discontinuous test and trial functions. This approach goes under the name of interior penalty discontinuous Galerkin method (IPDG) and was first investigated for the biharmonic problem by Baker [5] in the late 70ties; amongst others, we recall the following contributions, including the *p*- and *hp*-versions of the scheme: [14, 15, 20, 23, 29, 30, 35].

Standard results in the theory of the IPDG for the biharmonic problem state that the convergence is optimal with respect to the mesh size, and one order and a half suboptimal with respect to the polynomial degree. The reason for this is the use of polynomial inverse estimates while handling the penalization of the bilinear form. A milder suboptimality of half an order occurs also in the approximation of elliptic partial differential equations of second order but it can be removed under two assumptions: (i) the essential boundary conditions are homogeneous; (ii) the underlying mesh allows for the existence of a global $H^1$ piecewise polynomial with *hp*-optimal approximation properties. This was first observed in [31] using the error analysis in the $L^2$ norm for the local

<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</footnote>

<footnote>‡Fakultät für Mathematik, Universität Wien, 1090 Vienna, Austria, lorenzo.mascotto@univie.ac.at</footnote>

<footnote>§IMATI-CNR, 27100, Pavia, Italy</footnote>
discontinuous Galerkin method and later analyzed in [33] using the error analysis in the energy norm for the IPDG.

The aim of this paper is to prove that an analogous result holds true using the IPDG for the biharmonic problem with homogeneous essential boundary conditions. We shall mainly focus on the Hessian-Hessian formulation in [11], and show that similar hp-optimal results are valid also for the C$^0$-IPDG for the biharmonic problem in two and three dimensions, and the stream formulation of the Stokes problem in two dimensions. More precisely, the main result of this contribution is that the optimality with respect to the polynomial degree can be recovered using conforming tensor product-type meshes in two and three dimensions, and triangular meshes in two dimensions.

The analysis hinges upon showing the existence of global $H^2$ piecewise polynomial hp-optimal approximants over the given meshes. On tensor product-type meshes, we shall exhibit explicit hp-optimal bounds in the anisotropic Sobolev scale that are fully explicit with respect to the Sobolev regularity index. We are not aware of existing results in the literature covering the three dimensional case of this; so, this topic is interesting per se. On triangular meshes, such an $H^2$ approximant can be constructed following [36] but the error estimates depend on the Sobolev regularity index implicitly; we are not aware of a counterpart of this approach on tetrahedral meshes.

If essential boundary conditions are nonhomogeneous, for elliptic partial differential equations of second order, it is known that the method is p-suboptimal if essential boundary conditions are singular; see [19]. A similar result might be proven for the biharmonic problem. Showing this goes beyond the scopes of this work, but we shall provide numerical evidence that indeed p-suboptimality occurs.

**Structure of the paper.** In the remainder of this section, we present the notation we shall employ throughout the paper and recall the Hessian-Hessian formulation of the biharmonic problem. We design the corresponding IPDG scheme in Section 2. Section 3 is devoted to recalling and proving several hp-polynomial approximation and stability properties on tensor product-type meshes, which are instrumental to derive hp-optimal convergence rate of the IPDG in Section 4; here, we also discuss hp-optimality on triangular meshes. Variants of the above scheme are discussed in Section 5, namely we shall comment on C$^0$-IPDG in two and three dimensions, and the stream formulation of the Stokes problem in two dimensions. Numerical experiments are presented in Section 6 and conclusions are drawn in Section 7.

**Notation.** Let $D$ be a domain in $\mathbb{R}^d$, $d = 1, 2, 3$. We denote the Sobolev space of order $s \in \mathbb{N}$ by $H^s(D)$. Fractional order and negative order Sobolev spaces are defined using interpolation theory and duality, respectively. We denote the Lebesgue space of square integrable functions over $D$ by $L^2(D)$. We endow the above spaces with the inner product $(\cdot, \cdot)_{s,D}$, which induces the norm and seminorm $\|\cdot\|_{s,D}$ and $|\cdot|_{s,D}$. Further, given $\ell \in \mathbb{R}^+$ and $k \in \mathbb{N}$, we introduce the anisotropic Sobolev space:

$$ H^{\ell,k}(D) := \left\{ \begin{array}{ll} \displaystyle \| v \|_{H^\ell(D)} = \| \partial_{x_1}^{\alpha_1} \dots \partial_{x_d}^{\alpha_d} u \|_{0,D} < \infty & \forall \alpha = (\alpha_1, \dots, \alpha_d) \in \mathbb{N}^d \text{ s. th.} \\ \\ |\alpha| = k, \alpha_j \le 2 \ \forall j \in \{1, \dots, d\} \setminus \{J\} \text{ for one } J \in \{1, \dots, d\} \end{array} \right\}. \quad (1) $$

The space $H^{\ell,k}(D)$ contains $H^\ell$ functions with square integrable mixed derivatives of order $k$, where such mixed derivatives have fewer than three partial derivatives along all directions but one.

We introduce several weak differential operators: $\nabla$ as the gradient; $\Delta = \nabla \cdot \nabla$ as the Laplacian; $\Delta^2$ as the bilaplacian; $D^2 = \nabla^2$ as the Hessian.

Given $p \in \mathbb{N}$, let $\mathbb{P}_p(D)$ and $\mathbb{Q}_p(D)$ represent the spaces of polynomials of degree at most $p$ over a domain $D$ and degree at most $p$ in each variable over a $d$-dimensional cube $D$, respectively.

We denote the scalar product between two vectors $\mathbf{u}$ and $\mathbf{v}$ in $\mathbb{R}^d$ by $\mathbf{u} \cdot \mathbf{v}$. Further, we denote the entrywise product between two tensors $\mathbf{u}$ and $\mathbf{v}$ in $\mathbb{R}^{d \times d}$ by $\mathbf{u} : \mathbf{v}$.

Given two positive quantities $a$ and $b$, we use the short-hand notation “$a \lesssim b$” if there exists a positive constant $c$ independent of the discretization parameters such that “$a \le c b$”. Moreover, we write “$a \approx b$” if and only if “$a \lesssim b$” and “$b \lesssim a$”.
**The continuous problem.** Given $\Omega \subset \mathbb{R}^d$, $d = 2, 3$, a bounded polygonal/polyhedral domain with unit outward normal vector $\mathbf{n}_\Omega$ and $f \in L^2(\Omega)$, we consider the following biharmonic problem: find $u : \Omega \to \mathbb{R}$ such that

$$ \left\{ \begin{array}{ll} \Delta^2 u = f & \text{in } \Omega, \\ u = 0, & \mathbf{n}_\Omega \cdot \nabla u = 0 \text{ on } \partial\Omega. \end{array} \right. \qquad (2) $$

Introduce

$$ V := H_0^2(\Omega), \quad B(u, v) = \int_{\Omega} D^2 u : D^2 v \quad \forall u, v \in V. $$

A weak formulation of problem (2) reads

$$ \left\{ \begin{array}{ll} \text{find } u \in V \text{ such that} \\ B(u, v) = (f, v)_{0,\Omega} & \forall v \in V. \end{array} \right. \qquad (3) $$

Problem (3) is well posed; see, e.g., [10, Section 5.9].

*Remark 1.* We only consider homogeneous boundary conditions as they allow for fewer cumbersome technicalities when deriving hp-optimal error estimates with the IPDG. However, the hp-optimal results that we shall discuss in this paper can be extended to the case of piecewise smooth boundary conditions. Moreover, handling natural boundary conditions does not hide additional technicalities nor impact the optimal convergence rate in $p$; see, e.g., [8]. ■

# 2 Interior penalty discontinuous Galerkin methods for the biharmonic problem

In this section, we introduce the IPDG for problem (3). More precisely, we introduce regular tensor product-type and triangular meshes in Section 2.1; Section 2.2 is devoted to design the method; in Section 2.3, we discuss some properties of the bilinear forms involved in the method.

## 2.1 Meshes

We consider sequences ${\mathcal{T}_n}_n$ of decompositions of $\Omega$ into parallelograms/parallelepipeds in two and three dimensions, and triangular meshes in two dimensions. Henceforth, we shall refer to the case of parallelograms/parallelepipeds meshes as tensor product-type meshes. We denote the set of $(d-1)$-dimensional faces of each $\mathcal{T}_n$ by $\mathcal{F}_n$. We further distinguish interior and boundary faces introducing $\mathcal{F}_n^I$ and $\mathcal{F}_n^B$, respectively. The set of vertices of $\mathcal{T}_n$ is $\mathcal{V}_n$.

With each element $K \in \mathcal{T}_n$, we associate its diameter $h_K$ and collect all the diameters in $\mathbf{h}$: $\Omega \to \mathbb{R}^+$ defined as

$$ \mathbf{h} := \begin{cases} h_K & \text{if } \mathbf{x} \in K \in \mathcal{T}_n, \\ h_F & \text{if } \mathbf{x} \in F \in \mathcal{F}_n. \end{cases} $$

We denote the unit normal vector to any interior face $F$ with respect to the global orientation of $F$ by $\mathbf{n}$. We always fix an outward pointing vector for boundary faces.

Throughout, we assume quasi-uniformity of the mesh $\mathcal{T}_n$, i.e., there exists $c_{\text{mesh}} > 0$ such that

$$ c_{\text{mesh}}^{-1} h_{K_1} \le h_{K_2} \le c_{\text{mesh}} h_{K_1} \quad \forall K_1, K_2 \in \mathcal{T}_n. $$

In what follows, we consider a uniform distribution of the polynomial degree $p \in \mathbb{N}$ over $\mathcal{T}_n$. Variable polynomial degrees can be used [15, 32] but are not the main focus of this work.

Define the scalar broken Sobolev space

$$ \begin{aligned} H^1(\mathcal{T}_n, \Omega) &:= \{v \in L^2(\Omega) \mid v|_K \in H^1(K)\}, \\ H^{\ell,k}(\mathcal{T}_n) &:= \{v \in L^2(\Omega) \mid v|_K \in H^{\ell,k}(K)\} \quad \forall \ell \in \mathbb{R}^+, k \in \mathbb{N}. \end{aligned} $$

We introduce the jump $[\cdot] : L^2(F) \to \mathbb{R}$ and the average $\{\cdot\} : L^2(F) \to \mathbb{R}$ operators over a face $F \in \mathcal{F}_n$ as follows: if the face $F \in \mathcal{F}_n^I$ is an interior face shared by the elements $K^+$ and $K^-$,
given $v \in H^1(\mathcal{T}_n, \Omega)$ with restrictions over the elements $K^+$ and $K^-$ denoted by $v^+$ and $v^-$, respectively,

$$ \|v\|_F(\mathbf{x}) = \|v\|(\mathbf{x}) := v^+(\mathbf{x}) - v^-(\mathbf{x}), \quad \{v\}_F(\mathbf{x}) = \{v\}(\mathbf{x}) := \frac{1}{2}(v^+(\mathbf{x}) + v^-(\mathbf{x})) \quad \forall \mathbf{x} \in F. $$

On boundary faces $F \in \mathcal{F}_n^B$, we set $\|v\| = \{v\} = v$.

The definition of the two operators above can be extended to vector and tensor functions. We avoid using different notation, for the presentation's sake.

For any differential operator introduced in Section 1, we denote its elementwise version associated with $\mathcal{T}_n$ by adding a subscript $n$; for instance, $\nabla_n$ is the piecewise gradient over $\mathcal{T}_n$.

## 2.2 The method

For each mesh $\mathcal{T}_n$ and polynomial degree $p$, we introduce the space

$$ V_n := \{v_n \in L^2(\Omega) \mid v_{n|K} \in \mathbb{P}_p(K) \forall K \in \mathcal{T}_n\}. $$

We endow $V_n$ with the broken Sobolev norm and seminorm of order $s > 0$

$$ \|\cdot\|_{s, \mathcal{T}_n}^2 := \sum_{K \in \mathcal{T}_n} \|\cdot\|_{s, K}^2, \quad |\cdot|_{s, \mathcal{T}_n}^2 := \sum_{K \in \mathcal{T}_n} |\cdot|_{s, K}^2, $$

and introduce the lifting operator $\mathcal{L}: V_n \oplus V \to [V_n]^{d \times d}$, $d = 2, 3$, as

$$ \int_{\Omega} \mathcal{L}(u_n) : \mathbf{v}_n := \int_{\mathcal{F}_n} \{\mathbf{n} \cdot (\nabla \cdot \mathbf{v}_n)\} \|u_n\| - \int_{\mathcal{F}_n} \{(\mathbf{v}_n)\mathbf{n}\} \cdot \|\nabla u_n\| \quad \forall \mathbf{v}_n \in [V_n]^{d \times d}. \quad (4) $$

Let $\boldsymbol{\sigma}: \mathcal{F}_n \to \mathbb{R}^+$ and $\boldsymbol{\tau}: \mathcal{F}_n \to \mathbb{R}^+$ two piecewise constant stabilization parameters, which we shall fix in what follows; see Section 2.3 below. Given $D_n^2$ the piecewise Hessian over $\mathcal{T}_n$, we define the bilinear form $B_n: V_n \times V_n \to \mathbb{R}$ as

$$ B_n(u_n, v_n) := \int_{\Omega} D_n^2 u_n : D_n^2 v_n + \int_{\Omega} (\mathcal{L}(u_n) : D_n^2 v_n + \mathcal{L}(v_n) : D_n^2 u_n) \\ + \int_{\mathcal{F}_n} (\boldsymbol{\sigma}[u_n][v_n] + \boldsymbol{\tau}[\nabla u_n] \cdot [\nabla v_n]) \quad \forall u_n, v_n \in V_n. \quad (5) $$

We consider the following IPDG for the approximation of solutions to (3); see [15]:

$$ \left\{ \begin{array}{ll} \text{find } u_n \in V_n \text{ such that} & \\ B_n(u_n, v_n) = (f, v_n) & \forall v_n \in V_n. \end{array} \right. \quad (6) $$

Other dG formulations might be used [11, 20]; yet, we stick to (6) as it is the most similar to that of [33].

Eventually, we define the natural norm associated with the bilinear form $B_n(\cdot, \cdot)$: for all $v \in V_n \oplus V$,

$$ \begin{aligned} \|v\|_{\text{dG}}^2 &:= \|D_n^2 v\|_{0,\Omega}^2 + \sum_{F \in \mathcal{F}_n} \|\boldsymbol{\tau}^{\frac{1}{2}}[\nabla v]\|_{0,F}^2 + \sum_{F \in \mathcal{F}_n} \|\boldsymbol{\sigma}^{\frac{1}{2}}[v]\|_{0,F}^2 \\ &=: \|D_n^2 v\|_{0,\Omega}^2 + \|\boldsymbol{\tau}^{\frac{1}{2}}[\nabla v]\|_{0,\mathcal{F}_n}^2 + \|\boldsymbol{\sigma}^{\frac{1}{2}}[v]\|_{0,\mathcal{F}_n}^2, \end{aligned} \quad (7) $$

where

$$ (u, v)_{0, \mathcal{F}_n} := \sum_{F \in \mathcal{F}_n} (u, v)_{0, F}, \quad \|u\|_{0, \mathcal{F}_n}^2 := (u, u)_{0, \mathcal{F}_n}, \quad \forall u, v \in V_n \oplus V. $$
## 2.3 Properties of the lifting operator, bilinear form, and right-hand side

Following [20, Lemma 5.1], we have the following continuity result for the lifting operator $\mathcal{L}(\cdot)$ in (4).

**Lemma 2.1.** There exists a positive constant $c_S$ independent of $\mathbf{h}$ and $p$ such that

$$ \|\mathcal{L}(v_n)\|_{0,\Omega} \le c_S \left( \left\| \left( \frac{p^6}{\mathbf{h}^3} \right)^{\frac{1}{2}} \left[ v_n \right] \right\|_{0,\mathcal{F}_n} + \left\| \left( \frac{p^2}{\mathbf{h}} \right)^{\frac{1}{2}} \left[ \nabla v_n \right] \right\|_{0,\mathcal{F}_n} \right) \quad \forall v_n \in V_n. $$

Consequences of Lemma 2.1 are the coercivity and continuity of the dG bilinear form $B_n(\cdot, \cdot)$ under suitable choices of the dG parameters $\sigma$ and $\tau$, namely, given $c_\sigma$, $c_\tau > 2c_S + \frac{1}{2}$,

$$ \sigma := c_{\sigma} \frac{p^6}{\mathbf{h}^3}, \quad \tau := c_{\tau} \frac{p^2}{\mathbf{h}}. \qquad (8) $$

See, e.g., [20, Lemma 5.2].

**Lemma 2.2.** The bilinear form $B_n(\cdot, \cdot)$ is coercive and continuous:

$$ B_n(v, v) \ge \frac{1}{2} \|v\|_{dG}^2, \quad B_n(u, v) \le 2 \|u\|_{dG} \|v\|_{dG} \quad \forall u, v \in V_n \oplus V. \qquad (9) $$

We also have the continuity of the right-hand side in (6).

**Lemma 2.3.** There exists a positive constant $c_C$ independent of $\mathbf{h}$ and $p$ such that

$$ (f, v)_{0,\Omega} \le c_C \|f\|_{0,\Omega} \|v\|_{dG} \quad \forall v \in V_n \oplus V. $$

*Proof.* The assertion follows from the Cauchy-Schwarz inequality and Poincaré-Friedrichs type inequalities for piecewise $H^2$ functions; see, e.g., [12, Corollary 4.2]. □

The well posedness of method (6) follows combining the coercivity of $B_n(\cdot, \cdot)$ and the continuity of $(f, \cdot)_{0,\Omega}$ with respect to the $\|\cdot\|_{dG}$ norm.

**Theorem 2.4.** Method (6) is well posed with a priori estimates

$$ \|u_n\|_{dG} \le 2c_C \|f\|_{0,\Omega}. $$

# 3 Polynomial approximation results on tensor product-type meshes

In this section, we present several hp-optimal polynomial approximation results. They will be useful in Section 4 below when proving the convergence rate of the IPDG on tensor product-type meshes.

**Preliminary technical results.** We recall a result concerned with approximation properties of the $L^2$ and $H^1$ projectors in one dimension; see, e.g., [13, 32] for details. Given $\hat{I} := (-1, 1)$, let $\Pi_p^0 : L^2(\hat{I}) \to \mathbb{P}_p(\hat{I})$ and $\Pi_p^1 : H^1(\hat{I}) \to \mathbb{P}_p(\hat{I})$ be the $L^2$ and $H^1$ projectors onto $\mathbb{P}_p(\hat{I})$, respectively: for all $q_p \in \mathbb{P}_p(\hat{I})$,

$$ (q_p, v - \Pi_p^0 v)_{0, \hat{I}} = 0 \quad \forall v \in L^2(\hat{I}), \quad \begin{cases} (q'_p, (v - \Pi_p^1 v)')_{0, \hat{I}} = 0 \\ (v - \Pi_p^1 v)(\pm 1) = 0 \end{cases} \quad \forall v \in H^1(\hat{I}). $$

The following approximation estimates are valid; see, e.g., [32, Chapter 3] and [18, Lemma 3.5].
**Lemma 3.1.** For all $p \in \mathbb{N}$, given $u \in H^k(\hat{I})$ with $k \ge 1$ and $s := \min(k, p+1)$, there exists a positive constant $c$ independent of $p$ and $k$ such that

$$ \begin{aligned} \|u - \Pi_p^0 u\|_{0, \hat{I}} &\le cp^{-s}|u|_{s, \hat{I}}, \\ \|(u - \Pi_p^1 u)'\|_{0, \hat{I}} &\le cp^{-s+1}|u|_{s, \hat{I}}, \\ \|(u - \Pi_p^0 u)'\|_{0, \hat{I}} &\le cp^{-s+\frac{3}{2}}|u|_{s, \hat{I}}. \end{aligned} \quad (10) $$

Next, we recall best approximation estimates of the elementwise $L^2$ projector $\Pi_p^0: L^2(K) \to Q_p(K)$ with respect to the boundary $L^2$ norm on a $d$-dimensional tensor product-element $K$; see, e.g., [24, Lemma 3.9].

**Lemma 3.2.** Let $K \subset \mathbb{R}^d$, $d=2,3$ be a tensor product-type element, i.e., a parallelogram if $d=2$ and a parallelepiped if $d=3$, $h_K$ its diameter, and $u \in H^k(K)$, $k \ge 1$. Then, for all $1 \le s \le \min(p+1,k)$, $p \in \mathbb{N}$, there exists a positive constant $c > 0$ independent of $h_K$, $k$, and $p$ such that

$$ \|u - \Pi_p^0 u\|_{0, \partial K} \le c \left( \frac{h_K}{p+1} \right)^{s-\frac{1}{2}} \|u\|_{s, K}. $$

Now, we discuss hp-optimal approximation properties by means of global $H^2$ piecewise polynomials over $T_n$. We proceed in three steps.

**$H^2$ estimates in one dimension.** We recall the existence of a polynomial matching the values of a given function and its derivative at the endpoints of an interval with hp-optimal approximation properties; see [6, Corollary 2].

**Lemma 3.3.** Let $\hat{I} := (-1, 1)$ and $k \ge 0$. Then, for all $p \ge 3$, there exists a projection operator $\mathcal{P}_p^1: H^{k+2}(\hat{I}) \to \mathbb{P}_p(\hat{I})$ such that the following continuity properties are valid: for all $u \in H^{k+2}(\hat{I})$

$$ (\mathcal{P}_p^1 u)(\pm 1) = u(\pm 1), \quad (\mathcal{P}_p^1 u)'(\pm 1) = u'(\pm 1). \quad (11) $$

Moreover, the following upper bounds are valid: given $s := \min(k, p-1)$,

$$ \begin{aligned} \|(u - \mathcal{P}_p^1 u)''\|_{0, \hat{I}}^2 &\le \frac{(p-s-1)!}{(p+s-1)!} \|u^{(s+2)}\|_{0, \hat{I}}^2, \\ \|(u - \mathcal{P}_p^1 u)'\|_{0, \hat{I}}^2 &\le \frac{(p-s-1)!}{(p+s-1)!} \frac{1}{(p-1)p} \|u^{(s+2)}\|_{0, \hat{I}}^2, \\ \|u - \mathcal{P}_p^1 u\|_{0, \hat{I}}^2 &\le \frac{(p-s-1)!}{(p+s-1)!} c(p) \|u^{(s+2)}\|_{0, \hat{I}}^2, \end{aligned} \quad (12) $$

where $c(p)$ is independent of $k$ and

$$ c(p) \approx p^{-4}. $$

*Proof.* We provide the details of the proof in Appendix A for completeness. □

An application of Stirling's formula, see, e.g., [32, Corollary 3.12], yields the following p-optimal approximation result for the projector $\mathcal{P}_p^1$.

**Corollary 3.4.** Let $\hat{I} := (-1, 1)$ and $p \ge 3$. Then, there exists a positive constant $c_{1D} > 0$ independent of $p$ and $k$ such that the projection operator $\mathcal{P}_p^1: H^{k+2}(\hat{I}) \to \mathbb{P}_p(\hat{I})$, $k \ge 0$, from Lemma 3.3 satisfies, for $s := \min(k, p-1)$,

$$ \begin{aligned} &c_{1D}(p-1)^2 \|u - \mathcal{P}_p^1 u\|_{0, \hat{I}} + (p-1) \|(u - \mathcal{P}_p^1 u)'\|_{0, \hat{I}} + \|(u - \mathcal{P}_p^1 u)''\|_{0, \hat{I}} \\ &\le \left(\frac{e}{2}\right)^{\frac{s}{p-1}} (p-1)^{-s} \|u^{(s+2)}\|_{0, \hat{I}} \le \frac{e}{2} (p-1)^{-s} \|u^{(s+2)}\|_{0, \hat{I}}. \end{aligned} \quad (13) $$

Based on this, we derive stability properties for the projection operator constructed in Lemma 3.3.
**Corollary 3.5.** The projection operator $\mathcal{P}_p^1$. constructed in Lemma 3.3 satisfies the following stability properties: there exists a positive constant c independent of p and k such that, for all $u \in H^2(\hat{I})$,

$$ \begin{aligned} \|(\mathcal{P}_p^1 u)''\|_{0,\hat{I}} &\le \|u''\|_{0,\hat{I}}, \\ \|(\mathcal{P}_p^1 u)'\|_{0,\hat{I}} &\le \|u'\|_{0,\hat{I}} + p^{-1}\|u''\|_{0,\hat{I}}, \\ \|\mathcal{P}_p^1 u\|_{0,\hat{I}} &\le \|u\|_{0,\hat{I}} + cp^{-2}\|u''\|_{0,\hat{I}}. \end{aligned} \quad (14) $$

*Proof.* The first inequality follows from the definition of the operator $\mathcal{P}_p^1$ in (47) and the orthogonality properties of Legendre polynomials. As for the other two inequalities, it suffices to add and subtract $u'$ and $u$, respectively, use the triangle inequality, and apply the estimates in (12) with $s = 0$. □

**$H^2$ estimates on tensor product elements.** Given $\mathcal{P}_p^x$, $\mathcal{P}_p^y$, and $\mathcal{P}_p^z$ the one dimensional projectors from Lemma 3.3 applied to the restrictions of functions along the directions $x$, $y$, and $z$, respectively, we define the two and three dimensional operators

$$ \mathcal{P}_p := \mathcal{P}_p^x \mathcal{P}_p^y, \qquad \mathcal{P}_p := \mathcal{P}_p^x \mathcal{P}_p^y \mathcal{P}_p^z. $$

By tensor product arguments, we deduce the two and three dimensional counterparts of Lemma 3.3.

**Theorem 3.6.** Let $\hat{Q} := (-1, 1)^2$ be the reference square, with vertical edges $e_1$ and $e_3$, horizontal edges $e_2$ and $e_4$, and outward unit normal vector $\mathbf{n}_{\hat{Q}}$. Further, let $p \ge 3$ and $u \in H^{2,k+2}(\hat{Q})$, $k \ge 0$. Then, the following continuity properties on $\partial \hat{Q}$ are valid:

$$ \mathcal{P}_p u|_{e_i} = \begin{cases} \mathcal{P}_p^x(u|_{e_i}) & \text{if } i \text{ is even,} \\ \mathcal{P}_p^y(u|_{e_i}) & \text{if } i \text{ is odd,} \end{cases} \quad (\mathbf{n}_{\hat{Q}} \cdot \nabla \mathcal{P}_p u)|_{e_i} = \begin{cases} \mathbf{n}_{\hat{Q}} \cdot \nabla \mathcal{P}_p^y(u|_{e_i}) & \text{if } i \text{ is even,} \\ \mathbf{n}_{\hat{Q}} \cdot \nabla \mathcal{P}_p^x(u|_{e_i}) & \text{if } i \text{ is odd.} \end{cases} \quad (15) $$

In particular, $\mathcal{P}_p u$ and $u$ share the same value at the vertices of $\hat{Q}$. If $\hat{Q} = (-1, 1)^3$, analogous continuity properties are valid.

Let $s = \min(k, p-1)$. Then, the following optimal p-error estimates are valid. If $\hat{Q} = (-1, 1)^2$, there exists a positive constant c independent of p and k such that

$$ \|u - \mathcal{P}_p u\|_{0,\hat{Q}} \le cp^{-s-2} \left( \|\partial_x^{s+2} u\|_{0,\hat{I}} + \|\partial_y^{s+2} u\|_{0,\hat{I}} + \|\partial_x^2 \partial_y^s u\|_{0,\hat{I}} \right); \quad (16) $$

$$ \begin{aligned} \|\nabla(u - \mathcal{P}_p u)\|_{0,\hat{Q}} &\le cp^{-s-1} \left( \|\partial_x^{s+2} u\|_{0,\hat{Q}} + \|\partial_x \partial_y^{s+1} u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0,\hat{Q}} \right. \\ &\qquad \left. + \|\partial_x^s \partial_y^2 u\|_{0,\hat{Q}} + \|\partial_x^{s+1} \partial_y u\|_{0,\hat{Q}} + \|\partial_y^{s+2} u\|_{0,\hat{Q}} \right); \end{aligned} \quad (17) $$

$$ \begin{aligned} \|D^2(u - \mathcal{P}_p u)\|_{0,\hat{Q}} &\le cp^{-s} \left( \|\partial_x^{s+2} u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0,\hat{Q}} + \|\partial_y^{s+2} u\|_{0,\hat{Q}} \right. \\ &\qquad \left. + \|\partial_x^s \partial_y^2 u\|_{0,\hat{Q}} \|\partial_x^{s+1} \partial_y u\|_{0,\hat{Q}} + \|\partial_x \partial_y^{s+1} u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0,\hat{Q}} \right). \end{aligned} \quad (18) $$
Instead, if $\hat{Q} = (-1, 1)^3$, there exists a positive constant $c$ independent of $p$ and $k$ such that

$$ \|u - \mathcal{P}_p u\|_{0,\hat{Q}} \le cp^{-s-2} \left( \|\partial_x^{s+2} u\|_{0,\hat{I}} + \|\partial_y^{s+2} u\|_{0,\hat{I}} + \|\partial_z^{s+2} u\|_{0,\hat{I}} \right. \\ \left. \qquad + \|\partial_x^2 \partial_z^s u\|_{0,\hat{I}} + \|\partial_y^2 \partial_z^s u\|_{0,\hat{I}} + \|\partial_x^2 \partial_y^2 \partial_z^{s-2} u\|_{0,\hat{I}} \right); \quad (19) $$

$$ \|\nabla(u - \mathcal{P}_p u)\|_{0,\hat{Q}} \le cp^{-s-1} \left( \|\partial_x^{s+2} u\|_{0,\hat{Q}} + \|\partial_y^{s+2} u\|_{0,\hat{Q}} + \|\partial_z^{s+2} u\|_{0,\hat{Q}} + \|\partial_x \partial_y^{s+1} u\|_{0,\hat{Q}} \right. \\ \left. \qquad + \|\partial_y \partial_z^{s+1} u\|_{0,\hat{Q}} + \|\partial_z \partial_x^{s+1} u\|_{0,\hat{Q}} + \|\partial_x \partial_y^{s+1} u\|_{0,\hat{Q}} + \|\partial_y \partial_z^{s+1} u\|_{0,\hat{Q}} + \|\partial_z \partial_x^{s+1} u\|_{0,\hat{Q}} \\ \left. \qquad + \|\partial_x \partial_y^2 \partial_z^{s-1} u\|_{0,\hat{Q}} + \|\partial_y \partial_z^2 \partial_x^{s-1} u\|_{0,\hat{Q}} + \|\partial_z \partial_x^2 \partial_y^{s-1} u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0,\hat{Q}} + \|\partial_y^2 \partial_z^s u\|_{0,\hat{Q}} \right. \\ \left. \qquad + \|\partial_z^2 \partial_x^s u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_z^s u\|_{0,\hat{Q}} + \|\partial_y^2 \partial_x^s u\|_{0,\hat{Q}} + \|\partial_z^2 \partial_y^s u\|_{0,\hat{Q}} \right. \\ \left. \qquad + \|\partial_x^2 \partial_y^2 \partial_z^{s-2} u\|_{0,\hat{Q}} + \|\partial_y^2 \partial_z^2 \partial_x^{s-2} u\|_{0,\hat{Q}} + \|\partial_z^2 \partial_x^2 \partial_y^{s-2} u\|_{0,\hat{Q}} \right); \quad (20) $$

$$ \|D^2(u - \mathcal{P}_p u)\|_{0,\hat{Q}} \le cp^{-s} \left( \|\partial_x^{s+2} u\|_{0,\hat{Q}} + \|\partial_y^{s+2} u\|_{0,\hat{Q}} + \|\partial_z^{s+2} u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0,\hat{Q}} \right. \\ \left. \qquad + \|\partial_y^2 \partial_z^s u\|_{0,\hat{Q}} + \|\partial_z^2 \partial_x^s u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_z^s u\|_{0,\hat{Q}} + \|\partial_y^2 \partial_x^s u\|_{0,\hat{Q}} + \|\partial_x^{s+1} \partial_y u\|_{0,\hat{Q}} \right. \\ \left. \qquad + \|\partial_x^{s+1} \partial_z u\|_{0,\hat{Q}} + \|\partial_z^{s+1} \partial_x u\|_{0,\hat{Q}} + \|\partial_x \partial_y^{s+1} u\|_{0,\hat{Q}} + \|\partial_y \partial_z^{s+1} u\|_{0,\hat{Q}} + \|\partial_z \partial_x^{s+1} u\|_{0,\hat{Q}} \right. \\ \left. \qquad + \|\partial_x \partial_y \partial_z^s u\|_{0,\hat{Q}} + \|\partial_y \partial_z \partial_x^s u\|_{0,\hat{Q}} + \|\partial_z \partial_x \partial_y^s u\|_{0,\hat{Q}} + \|\partial_x \partial_y^2 \partial_z^{s-1} u\|_{0,\hat{Q}} + \|\partial_y \partial_z^2 \partial_x^{s-1} u\|_{0,\hat{Q}} \right. \\ \left. \qquad + \|\partial_z \partial_x^2 \partial_y^{s-1} u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y \partial_z^{s-1} u\|_{0,\hat{Q}} + \|\partial_y^2 \partial_z \partial_x^{s-1} u\|_{0,\hat{Q}} + \|\partial_z^2 \partial_x \partial_y^{s-1} u\|_{0,\hat{Q}} \right. \\ \left. \qquad + \|\partial_x^2 \partial_y^2 \partial_z^{s-2} u\|_{0,\hat{Q}} + \|\partial_y^2 \partial_z^2 \partial_x^{s-2} u\|_{0,\hat{Q}} + \|\partial_z^2 \partial_x^2 \partial_y^{s-2} u\|_{0,\hat{Q}} \right). \quad (21) $$

*Proof.* As for the two dimensional case, we refer to [6] and present some details in Appendix B for completeness. We report the proof for the three dimensional estimates in Appendix C. $\square$

**Remark 2.** In Theorem 3.6, the p-explicit approximation bounds are derived for functions in the anisotropic Sobolev space $H^{2,k+2}(\hat{Q})$ defined in (1). This is due to the use of tensor product arguments based on one dimensional approximation results. The upside of considering such a regularity is that the error bounds are fully explicit with respect to the polynomial degree $p$ and the Sobolev regularity index $k$; this is crucial to derive exponential convergence rate of the hp-version Galerkin methods. Moreover, estimates as in Theorem 3.6 naturally appear in the approximation of solutions to partial differential equations with anisotropic Sobolev regularity. We cannot extend the techniques derived so far to the case of triangular meshes. Notably, in that case, the estimates also depend on the Sobolev regularity index $k$; see Section 4.3 below for more details. ■

**Remark 3.** In the error estimates of Theorem 3.6, it is possible to switch the order of the mixed derivatives appearing in the norms on the right-hand sides. This is of paramount importance to notice in case one wishes to approximate functions in anisotropic Sobolev spaces. ■

So far, we proved error estimates in two and three dimensions for tensor product-type meshes. These results can be extended to higher dimensions. The corresponding proof hinges upon the same arguments used in Theorem 3.6: error estimates in one dimension; stability of the $H^2$ projector $\mathcal{P}_p^1$; identities as in (61) in general dimension.

By a scaling argument, we deduce hp-optimal error estimates for functions with isotropic Sobolev regularity on tensor product-type elements.

**Corollary 3.7.** Let $K$ be any parallelogram with vertical edges $e_1$ and $e_3$, horizontal edges $e_2$ and $e_4$, and outward unit normal vector $\mathbf{n}_K$, $p \ge 3$, and $u \in H^{k+2}(K)$, $k \ge 0$. Then, the following continuity properties on $\partial K$ are valid:

$$ \mathcal{P}_p u|_{F_i} = \begin{cases} \mathcal{P}_p^x(u|_{F_i}) & \text{if } i \text{ is even,} \\ \mathcal{P}_p^y(u|_{F_i}) & \text{if } i \text{ is odd,} \end{cases} \quad (\mathbf{n}_K \cdot \nabla \mathcal{P}_p u)|_{F_i} = \begin{cases} \mathbf{n}_K \cdot \nabla \mathcal{P}_p^y(u|_{F_i}) & \text{if } i \text{ is even,} \\ \mathbf{n}_K \cdot \nabla \mathcal{P}_p^x(u|_{F_i}) & \text{if } i \text{ is odd.} \end{cases} $$

In particular, $\mathcal{P}_p u$ and $u$ share the same value at the vertices of $K$. If $K$ is a parallelepipeds, analogous continuity properties hold true.
Besides, the following optimal p-error estimates are valid in two and three dimensions: given s = min(k, p - 1) and either a parallelogram or parallelepiped K, there exists a positive constant c in- dependent of h, k, and p such that

$$
\frac{p^2}{h_K^2} \|u - P_p u\|_{0,K} + \frac{p}{h_K} \|\nabla(u - P_p u)\|_{0,K} + \|D^2(u - P_p u)\|_{0,K} \le c \frac{h^s}{p^s} |u|_{s+2,K}. \quad (22)
$$

**Global H** 2 **estimates.** Corollary 3.7 implies hp-optimal error estimates by means of global H piecewise polynomials over the given tensor product-type meshes. More precisely, we have the following result.

**Theorem 3.8.** Let $\mathcal{T}_n$ be a tensor product-type mesh over $\Omega \subset \mathbb{R}^d$, $d = 2,3$, of either parallelograms ($d=2$) or parallelepipeds ($d=3$) as in Section 2.1; $k \ge 0$; $p \in \mathbb{N}$; $s := \min(p, k-1)$. Then there exists $\mathcal{P}_p: H^2(\Omega) \cap H^{2,2d}(\mathcal{T}_n) \to H^2(\Omega) \cap V_n$ and a positive constant $c$ independent of $\mathbf{h}, k$, and $p$ such that, for all $u \in H^2(\Omega) \cap H^{k+2}(\mathcal{T}_n)$,

$$
p^2 \|\mathbf{h}^{-2}(u - \mathcal{P}_p u)\|_{0,\Omega} + p \|\mathbf{h}^{-1}\nabla(u - \mathcal{P}_p u)\|_{0,\Omega} + \|D^2(u - \mathcal{P}_p u)\|_{0,\Omega} \le cp^{-s} \|\mathbf{h}^s D^{s+2}u\|_{0,\mathcal{T}_n}.
$$

The operator $\mathcal{P}_p$ is defined on the space $H^2(\Omega) \cap H^{2,2d}(\mathcal{T}_n)$. The reason for this will be apparent from the construction of the operator itself in Appendices A, B, and C. In particular, we shall need that mixed derivatives of order 2 in each direction are square integrable.

4 Optimal convergence

In this section, we show the main result of this paper: method (6) is hp-optimal on different types
of meshes if the biharmonic problem (3) is endowed with homogeneous boundary conditions. We
focus on tensor product-type meshes in 2D and 3D in Section 4.1; estimates in weaker norms are
shown in Section 4.2; a guideline on how to show hp-optimal error estimates on triangular meshes
is provided in Section 4.3.

4.1 Optimal convergence on tensor product-type meshes

In this section, we prove hp-optimal error estimates for tensor product-type meshes in two and three dimensions. The analysis bases on the existence of the global H² polynomial projector from Theorem 3.8.

Henceforth, given $D \subset \mathbb{R}^d$, $d = 2,3$, either a parallelogram (in two dimensions) or a paral-
lelepiped (in three dimensions), we introduce $\Pi_p^0$ as the composition of the one dimensional $L^2$
projectors in Lemma 3.1 along the coordinate directions. In other words, we let

$$
\Pi_p^0 = \Pi_p^{0,x} \Pi_p^{0,y} \quad \text{if} \quad d=2, \qquad \Pi_p^0 = \Pi_p^{0,x} \Pi_p^{0,y} \Pi_p^{0,z} \quad \text{if} \quad d=3.
$$

**Theorem 4.1.** Let $\mathcal{T}_n$ be a quasi-uniform tensor product-type mesh of either parallelograms (2D) or parallelepipeds (3D) as in Section 2.1. Let $u$ and $u_n$ be the solutions to (3) and (6), respectively. Assume that $u \in H^{s+2}(\Omega)$, $s := \min(k, p-1)$, $p \ge 2$. Then, there exists a positive constant $c$ independent of $h$, $k$, and $p$ such that

$$
\|u - u_n\|_{\mathrm{dG}} \le c \left(\frac{h}{p}\right)^s \|u\|_{s+2, \Omega}.
$$

In other words, for homogeneous essential boundary conditions, method (6) is hp-optimal.

*Proof.* As in [20, Theorem 5.3], we have the following error splitting:

$$
\|u - u_n\|_{\mathrm{dG}} \le 5 \inf_{v_n \in V_n} \|u - v_n\|_{\mathrm{dG}} + 2 \sup_{\xi \in V_n} \frac{|B_n(u - u_n, \xi)|}{\|\xi\|_{\mathrm{dG}}}. \quad (23)
$$

The first term represents a best polynomial approximation result; the second measures the incon-
sistency of the scheme.
We have standard hp-optimal (piecewise) polynomial estimates on the first term:

$$ \inf_{v_n \in V_n} \|u - v_n\|_{dG} \lesssim \left(\frac{h}{p}\right)^s \|u\|_{s+2, \Omega}. \quad (24) $$

To show (24), it suffices to pick $v_n = P_p u$, where $P_p$ is as in Theorem 3.8. The jump terms on interior edges vanish due to the global $H^2$ regularity of $u$ and $P_p u$; the jump terms on boundary edges vanish since $u$ has homogeneous essential boundary conditions and $P_p u$ preserves polynomial essential boundary conditions.

Thus, we only need to derive estimates for the inconsistency term, i.e., the second term on the right-hand side of (23). Proceeding as in [20, Lemma 5.4], we set $\eta := u - u_n$ and write

$$ \begin{aligned} B_n(u - u_n, \xi) = & \int_{\mathcal{F}_n} \{\mathbf{n} \cdot (\nabla \cdot D_n^2 \eta) - \mathbf{n} \cdot (\nabla \cdot \mathbf{\Pi}_p^0 D_n^2 \eta)\} \|\xi\| \\ & - \int_{\mathcal{F}_n} \{(D_n^2 \eta)\mathbf{n} - (\mathbf{\Pi}_p^0 D_n^2 \eta)\mathbf{n}\} \|\nabla \xi\| =: A_1 + A_2. \end{aligned} \quad (25) $$

Since $\mathbf{\Pi}_p^0$ preserves polynomials of degree $p$, we have

$$ D_n^2 \eta - \mathbf{\Pi}_p^0 D_n^2 \eta = D_n^2(u - u_n) - \mathbf{\Pi}_p^0(D_n^2 u) + \mathbf{\Pi}_p^0(D_n^2 u_n) = D_n^2 u - \mathbf{\Pi}_p^0 D_n^2 u. \quad (26) $$

Using (7), (8), and (26), we write

$$ A_1 \lesssim \frac{h^{\frac{3}{2}}}{p^3} \|\nabla \cdot (D_n^2 \eta - \mathbf{\Pi}_p^0 D_n^2 \eta)\|_{0, \mathcal{F}_n} \|\boldsymbol{\sigma}^{\frac{1}{2}} \|\xi\|_{0, \mathcal{F}_n} \lesssim \frac{h^{\frac{3}{2}}}{p^3} \|\nabla \cdot (D_n^2 u - \mathbf{\Pi}_p^0 D_n^2 u)\|_{0, \mathcal{F}_n} \|\xi\|_{dG}. \quad (27) $$

Next, we focus on the term $A_2$ and use similar manipulations:

$$ A_2 \lesssim \frac{h^{\frac{1}{2}}}{p} \|D_n^2 \eta - \mathbf{\Pi}_p^0 D_n^2 \eta\|_{0, \mathcal{F}_n} \|\boldsymbol{\tau}^{\frac{1}{2}} \|\xi\|_{0, \mathcal{F}_n} \lesssim \frac{h^{\frac{1}{2}}}{p} \|D_n^2 u - \mathbf{\Pi}_p^0 D_n^2 u\|_{0, \mathcal{F}_n} \|\xi\|_{dG}. \quad (28) $$

We insert (27) and (28) in (25) and get

$$ |B_n(u - u_n, \xi)| \lesssim \left[ \frac{h^{\frac{3}{2}}}{p^3} \|\nabla \cdot (D_n^2 u - \mathbf{\Pi}_p^0 D_n^2 u)\|_{0, \mathcal{F}_n} + \frac{h^{\frac{1}{2}}}{p} \|D_n^2 u - \mathbf{\Pi}_p^0 D_n^2 u\|_{0, \mathcal{F}_n} \right] \|\xi\|_{dG}. $$

Thus, we arrive at the following bound on the inconsistency term in (23):

$$ \sup_{\xi \in V_n} \frac{|B_n(u - u_n, \xi)|}{\|\xi\|_{dG}} \lesssim \frac{h^{\frac{3}{2}}}{p^3} \|\nabla \cdot (D_n^2 u - \mathbf{\Pi}_p^0 D_n^2 u)\|_{0, \mathcal{F}_n} + \frac{h^{\frac{1}{2}}}{p} \|D_n^2 u - \mathbf{\Pi}_p^0 D_n^2 u\|_{0, \mathcal{F}_n} =: T_1 + T_2. \quad (29) $$

To conclude, we are left to estimate the terms $T_1$ and $T_2$. First, we consider the term $T_1$ and estimate the $L^2$ norm on each face $F$ of the following vector:

$$ \nabla \cdot (D_n^2 u - \mathbf{\Pi}_p^0 D_n^2 u) = \left[ \frac{\partial_x (\partial_x^2 u - \Pi_p^0 \partial_x^2 u) + \partial_y (\partial_{xy} u - \Pi_p^0 \partial_{xy} u)}{\partial_x (\partial_{yx} u - \Pi_p^0 \partial_{yx} u) + \partial_y (\partial_y^2 u - \Pi_p^0 \partial_y^2 u)} \right]. $$

Using the triangle inequality, we have to estimate the four terms appearing above. Since each of them can be handled similarly, we only show the estimates for one of them. Moreover, we show local estimates, say, on each face $F$ belonging (at least) to the element $K$, as the global counterparts follow immediately. The triangle inequality implies

$$ \|\partial_x (\partial_x^2 u - \Pi_p^0 \partial_x^2 u)\|_{0,F} \le \|\partial_x^3 u - \Pi_{p-1}^0 \partial_x^3 u\|_{0,F} + \|\Pi_{p-1}^0 \partial_x^3 u - \partial_x \Pi_p^0 \partial_x^2 u\|_{0,F} =: Z_1 + Z_2. \quad (30) $$

As for the term $Z_1$, Lemma 3.2 yields

$$ Z_1 \lesssim \left(\frac{h}{p}\right)^{s-\frac{3}{2}} \|\partial_x^3 u\|_{s-1, K} \le \left(\frac{h}{p}\right)^{s-\frac{3}{2}} \|u\|_{s+2, K} \quad \forall s \ge 0. \quad (31) $$
As for the term $Z_2$, we use the polynomial inverse trace inequality [32, eq. (4.6.4)] (which gives an order suboptimality in $p$) and the triangle inequality:

$$ 
\begin{aligned}
Z_2 &\lesssim \frac{p}{h^{\frac{1}{2}}} \|\Pi_{p-1}^0 \partial_x^3 u - \partial_x \Pi_p^0 \partial_x^2 u\|_{0,K} \\
&\le \frac{p}{h^{\frac{1}{2}}} \left[ \|\partial_x^3 u - \Pi_{p-1}^0 \partial_x^3 u\|_{0,K} + \|\partial_x(\partial_x^2 u - \Pi_p^0 \partial_x^2 u)\|_{0,K} \right] = Z_{2,1} + Z_{2,2}.
\end{aligned} 
\quad (32) $$

Using the approximation properties of the $L^2$ projector yields

$$ Z_{2,1} \lesssim \frac{p}{h^{\frac{1}{2}}} \left(\frac{h}{p}\right)^{s-1} \|\partial_x^3 u\|_{s-1,K} \le \frac{p}{h^{\frac{1}{2}}} \left(\frac{h}{p}\right)^{s-1} \|u\|_{s+2,K}. $$

On the other hand, standard manipulations entail

$$ 
\begin{aligned}
Z_{2,2} &\lesssim \frac{p}{h^{\frac{1}{2}}} \left[ \|\partial_x(\partial_x^2 u - \Pi_p^{0,y} \partial_x^2 u)\|_{0,K} + \|\partial_x(\Pi_p^{0,y} \partial_x^2 u - \Pi_p^{0,y} \Pi_p^{0,x} \partial_x^2 u)\|_{0,K} \right] \\
&\le \frac{p}{h^{\frac{1}{2}}} \left[ \|\partial_x(\partial_x^2 u - \Pi_p^{0,y} \partial_x^2 u)\|_{0,K} + \|\partial_x(\partial_x^2 u - \Pi_p^{0,x} \partial_x^2 u)\|_{0,K} \right].
\end{aligned} 
$$

Using the approximation properties of the $L^2$ projection on the first term, see Lemma 3.1, and (10) (which gives half an order suboptimality in $p$) on the second term imply

$$ Z_{2,2} \lesssim \frac{p^{\frac{3}{2}}}{h^{\frac{1}{2}}} \left(\frac{h}{p}\right)^{s-1} \|\partial_x^3 u\|_{s-1,K} \lesssim \frac{p^{\frac{3}{2}}}{h^{\frac{1}{2}}} \left(\frac{h}{p}\right)^{s-1} \|u\|_{s+2,K}. $$

We plug the estimates on $Z_{2,1}$ and $Z_{2,2}$ in (32) and arrive at

$$ Z_2 \lesssim \frac{p^{\frac{3}{2}}}{h^{\frac{1}{2}}} \left(\frac{h}{p}\right)^{s-1} \|u\|_{s+2,K}, \quad s \ge 0. \quad (33) $$

Inserting next (31) and (33) in (30) gives

$$ \|\partial_x(\partial_x^2 u - \Pi_p^0 \partial_x^2 u)\|_{0,F} \lesssim \frac{p^{\frac{3}{2}}}{h^{\frac{1}{2}}} \left(\frac{h}{p}\right)^{s-1} \|u\|_{s+2,K} = \frac{p^{\frac{5}{2}}}{h^{\frac{3}{2}}} \left(\frac{h}{p}\right)^s \|u\|_{s+2,K}, \quad s \ge 0. $$

In overall, also recalling the factor $\frac{h^{\frac{3}{2}}}{p^3}$ in front of $T_1$, we deduce

$$ T_1 \lesssim p^{-\frac{1}{2}} \left(\frac{h}{p}\right)^s \|u\|_{s+2,\Omega}. \quad (34) $$

As for the term $T_2$ on the right-hand side of (29), we use Lemma 3.2 and obtain

$$ T_2 \lesssim p^{-\frac{1}{2}} \left(\frac{h}{p}\right)^s \|u\|_{s+2,\Omega}. \quad (35) $$

Inserting (34) and (35) in (29) leads to

$$ \sup_{\xi \in V_n} \frac{|B_n(u - u_n, \xi)|}{\|\xi\|_{dG}} \lesssim p^{-\frac{1}{2}} \left(\frac{h}{p}\right)^s \|u\|_{s+2,\Omega} \le \left(\frac{h}{p}\right)^s \|u\|_{s+2,\Omega}. \quad (36) $$

A combination of (25), (24), and (36) yields the assertion. $\square$

**Remark 4.** Theorem 4.1 states that the DG scheme (6) delivers hp-optimal error estimates for the biharmonic problem (3) with homogeneous boundary conditions. Following [19], it is possible to show that hp-optimal convergence is achieved with polynomial/smooth boundary conditions as well. On the other hand, singularities arising on the boundary of the domain yield a suboptimal p-convergence, which we shall highlight in the numerical experiments; see Section 6.3 below. ■
## 4.2 Optimal convergence on tensor product-type meshes in weaker norms

For tensor product-type meshes, we discuss hp-optimal error estimates in weaker norms under the assumption that $\Omega$ is convex. Notably, we are interested in the $L^2$ and broken-$H^1$ error estimates. Since analogous techniques are employed, we focus on $L^2$ error estimates only.

**Theorem 4.2.** Let $\mathcal{T}_n$ be a quasi-uniform tensor product-type mesh of either parallelograms (2D) or parallelepipeds (3D) as in Section 2.1. Let $u$ and $u_n$ be the solutions to (3) and (6), respectively. Then, if $\Omega$ is convex, there exists a positive constant $c$ independent of $h$ and $p$ such that

$$ \|u - u_n\|_{0,\Omega} \le c \left(\frac{h}{p}\right)^2 \|u - u_n\|_{dG}. $$

Together with Theorem 4.1, this implies hp-optimal convergence in the $L^2$ norm if additional regularity on $u$ is granted.

*Proof.* Consider the auxiliary problem: find $\Phi$ such that

$$ \begin{cases} \Delta^2 \Phi = u - u_n & \text{in } \Omega, \\ u = 0, \quad \mathbf{n}_{\Omega} \cdot \nabla u = 0 & \text{on } \partial \Omega. \end{cases} $$

In weak formulation, this problem reads

$$ \left\{ \begin{array}{ll} \text{find } u \in V \text{ such that} & \\ B(u, v) = (u - u_n, v)_{0,\Omega} & \forall v \in V. \end{array} \right. \qquad (37) $$

Due to the convexity of $\Omega$, the solution $\Phi$ to (37) belongs to $H^4(\Omega)$ and the following regularity estimate holds true:

$$ \|\Phi\|_{4,\Omega} \lesssim \|u - u_n\|_{0,\Omega}; \qquad (38) $$

see, e.g., [22] and [26, Theorem 6, page 182] for the two and three dimensional cases, respectively. Using (37), splitting the integrals elementwise, and integrating by parts imply

$$ \|u - u_n\|_{0,\Omega}^2 = (u - u_n, \Delta^2 \Phi)_{0,\Omega} = \sum_{K \in \mathcal{T}_n} \left( -(\nabla(u - u_n), \nabla \Delta \Phi)_{0,K} + (u - u_n, \nabla \Delta \Phi \cdot \mathbf{n}_K)_{0,\partial K} \right). $$

Integrating by parts, and using (4) and the smoothness of $\Phi$ give

$$ \begin{aligned} \|u - u_n\|_{0,\Omega}^2 = & \sum_{K \in \mathcal{T}_n} (D^2(u - u_n), D^2\Phi)_{0,K} + (\llbracket u - u_n \rrbracket, \{\nabla \Delta \Phi \cdot \mathbf{n}\})_{0,\mathcal{F}_n} \\ & + (\llbracket \Phi \rrbracket, \{\nabla \Delta(u - u_n) \cdot \mathbf{n}\})_{0,\mathcal{F}_n} - (\llbracket \nabla(u - u_n) \rrbracket, \{D^2\Phi \cdot \mathbf{n}\})_{0,\mathcal{F}_n} \\ & - (\llbracket \nabla \Phi \rrbracket, \{D^2(u - u_n) \cdot \mathbf{n}\})_{0,\mathcal{F}_n} + (\boldsymbol{\sigma}^{\frac{1}{2}}[u_n], \llbracket v_n \rrbracket)_{0,\mathcal{F}_n} + (\boldsymbol{\tau}^{\frac{1}{2}}[\nabla u_n], \llbracket \nabla v_n \rrbracket)_{0,\mathcal{F}_n}. \end{aligned} $$

As, e.g., in [21, Theorem 4.5], we write the right-hand side as a combination of $B_n(u - u_n, \Phi)$ and a reminder $R$. Standard computations reveal that

$$ \begin{aligned} R = & [(\llbracket u - u_n \rrbracket, \{\nabla \Delta \Phi \cdot \mathbf{n}\})_{0, \mathcal{F}_n} - (\llbracket \nabla (u - u_n) \rrbracket, \{D^2 \Phi \cdot \mathbf{n}\})_{0, \mathcal{F}_n} - (\mathcal{L}(u - u_n), \Pi_p^0 D_n^2 \Phi)_{0, \Omega}] \\ & + [(\llbracket \Phi \rrbracket, \{\nabla \Delta (u - u_n) \cdot \mathbf{n}\})_{0, \mathcal{F}_n} - (\llbracket \nabla \Phi \rrbracket, \{D^2 (u - u_n) \cdot \mathbf{n}\})_{0, \mathcal{F}_n} - (\mathcal{L}\Phi, \Pi_p^0 (D_n^2 (u - u_n)))_{0, \Omega}] \\ = & A + B. \end{aligned} $$

Due to the regularity of $\Phi$, we have $B = 0$. We deduce

$$ \|u - u_n\|_{0,\Omega}^2 = B_n(u - u_n, \Phi) + R, \qquad (39) $$

where

$$ R = (\llbracket u - u_n \rrbracket, \{\mathbf{n} \cdot \nabla (D^2 \Phi - \Pi_p^0 D^2 \Phi)\})_{0, \mathcal{F}_n} - (\llbracket \nabla (u - u_n) \rrbracket, \{\mathbf{n} \cdot (D^2 \Phi - \Pi_p^0 D^2 \Phi)\})_{0, \mathcal{F}_n}. $$
Thanks to the smoothness of $\mathcal{P}_p\Phi$, where $\mathcal{P}_p$ is as in Theorem 3.8, we can write

$$B_n(u - u_n, \Phi) = B_n(u - u_n, \Phi - \mathcal{P}_p\Phi).$$

In fact, $\mathcal{L}(\mathcal{P}_p\Phi) = 0$. Using the continuity property (9), the smoothness and the homogeneous essential boundary conditions of $\Phi - \mathcal{P}_p\Phi$, the approximation properties of $\mathcal{P}_p$ in (22), and the stability property (38), we deduce

$$ \begin{aligned} B_n(u - u_n, \Phi) &\lesssim \|u - u_n\|_{\mathrm{dG}} \|D_n^2(\Phi - \mathcal{P}_p\Phi)\|_{0,\Omega} \lesssim \left(\frac{h}{p}\right)^2 \|u - u_n\|_{\mathrm{dG}} \|\Phi\|_{4,\Omega} \\ &\lesssim \left(\frac{h}{p}\right)^2 \|u - u_n\|_{\mathrm{dG}} \|u - u_n\|_{0,\Omega}. \end{aligned} \quad (40) $$

We are left with estimating the remainder $R$ on the right-hand side of (39). To this aim, we proceed as in the proof of Theorem 4.1. More precisely, we apply the estimates (34) and (35) with $s = 2$ and use the stability bound (38) to deduce

$$ R \lesssim \|u - u_n\|_{\mathrm{dG}} \left(\frac{h}{p}\right)^2 \|\Phi\|_{4,\Omega} \lesssim \left(\frac{h}{p}\right)^2 \|u - u_n\|_{\mathrm{dG}} \|u - u_n\|_{0,\Omega}. \quad (41) $$

The assertion follows combining (40) and (41) in (39). $\square$

*Remark 5.* As customary, the convexity assumption in Theorem 4.2 can be relaxed by demanding that the domain $\Omega$ implies $H^4$ regularity on the solution of the auxiliary problem (37). ■

## 4.3 Optimal convergence on triangular meshes

In this section, we provide the guidelines to extend Theorems 4.1 and 4.2 to the case of triangular meshes. As in the analysis above, *hp*-optimal convergence estimates follow from the existence of an $H^2$ conforming polynomial interpolant over a given triangular mesh with *hp*-optimal approximation properties. We undertake a different approach compared to the case of tensor product-type meshes as the arguments therein employed do not apply on simplicial meshes.

In [36], Suri constructed an interpolant $\mathcal{I}_p^{\mathrm{BS}}$ of $H^2$ functions into the space of global $H^2$ piece-wise polynomials on triangular meshes. Such an interpolant satisfies the following approximation properties: there exists a positive constant $c$ only depending on $\mathcal{T}_n$ such that

$$ |u - \mathcal{I}_p^{\mathrm{BS}} u|_{\ell, \Omega} \le cp^{-(k-\ell)} |u|_{k, \Omega} \quad \ell = 0, 1, 2. \quad (42) $$

The interpolant $\mathcal{I}_p^{\mathrm{BS}}$ is constructed generalizing to the $H^2$ case the techniques for the $H^1$ case developed by Babuška and Suri in [4]. First, a polynomial with *p*-optimal approximation properties is constructed on each triangle. Then, global $H^2$ regularity is recovered using an $H^2$ stable lifting operator, which preserves essential polynomial conditions on the boundary of a triangle. So, *p*-optimal error bounds of the IPDG on triangular meshes for the biharmonic problem follow using (42) in the analysis of Section 4.1.

We report the convergence results for triangular meshes.

**Theorem 4.3.** Let $\mathcal{T}_n$ be a quasi-uniform triangular mesh as in Section 2.1. Let $u$ and $u_n$ be the solutions to (3) and (6), respectively. Assume that $u \in H^{s+2}(\Omega)$, $s := \min(k, p-1)$, $p \ge 2$. Then, there exists a positive constant $c$ independent of $h$ and $p$ such that

$$ \|u - u_n\|_{\mathrm{dG}} \le c \left(\frac{h}{p}\right)^s \|u\|_{s+2,\Omega}, \quad \|u - u_n\|_{0,\Omega} \le c \left(\frac{h}{p}\right)^2 \|u - u_n\|_{\mathrm{dG}}. $$

In other words, for homogeneous essential boundary conditions, method (6) is *hp*-optimal in the energy and $L^2$ norms also on triangular meshes.
*Remark 6.* The *p*-optimal interpolant in (42) is derived based on the trigonometric polynomials approximation theory. The error bound is *p*-optimal with respect to the Sobolev index *k* but contains unknown constants depending on *k*; this is due to the use of the Stein's extension operator for Lipschitz domains [34]. The corresponding bound cannot be used to derive exponential convergence for the *hp*-version of Galerkin methods but is enough to study *p*-optimal algebraic convergence. So, one of the reasons why we developed the theory of Section 3 is that the *hp*-optimal error estimates are explicit in the Sobolev index *k* on tensor product-type meshes.

■

*Remark 7.* To the best of our knowledge, so far, the Suri-type *p*-optimal interpolation operator in (42) can be constructed in two dimensions only. There are two main reasons for this fact: (i) an explicit construction of general order $H^2$ conforming piecewise polynomial spaces on tetrahedral meshes is still missing; (ii) $H^2$ stable lifting operators for tetrahedra are also unknown. Notwithstanding, it is possible to use $H^1$-conforming *hp*-optimal polynomial interpolation operator on tetrahedral meshes. The resulting error bound for the IPDG is then suboptimal only by half an order of *p*, thus improving the expected one and a half an order suboptimality.

■

*Remark 8.* There exist at least other two different constructions for $H^2$ stable lifting operators. In [27], an $H^2$ stable lifting operator has been derived to prove the *p*-robustness of the inf-sup constant in H(div)-FEM for the Stokes problem in two dimensions. In [1, Theorem 2.2], the existence of an $H^2$ stable lifting operator was shown, which preserves essential polynomial conditions on the boundary of a triangle to design preconditioners for high-order $H^2$ FEM on triangular meshes [2].

■

# 5 Convergence for other IPDG formulations of the biharmonic problem

In this section, we discuss other IPDG formulations that can be proven to be *hp*-optimal using similar techniques to those of Section 4.

## 5.1 $C^0$-IPDG for the biharmonic problem

In this section, we consider $C^0$-IPDG for the biharmonic problem; see [11, 17]. Given the space $V_n$ of piecewise continuous polynomials over a mesh $\mathcal{T}_n$ as in Section 2.1, we introduce the lifting operator $\mathcal{L}: V_n \oplus V \to [V_n]^{d \times d}$ as

$$ \int_{\Omega} \mathcal{L}(u_n) : \mathbf{v}_n := - \int_{\mathcal{F}_n} (\{(\mathbf{v}_n)\mathbf{n}\} \cdot \mathbf{n}) (\llbracket \nabla u_n \rrbracket \cdot \mathbf{n}) \quad \forall \mathbf{v}_n \in [V_n]^{d \times d} $$

and the discrete bilinear form

$$ B_n(u_n, v_n) := \int_{\Omega} D_n^2 u_n : D_n^2 v_n + \int_{\Omega} (\mathcal{L}(u_n) : D_n^2 v_n + \mathcal{L}(v_n) : D_n^2 u_n) + \int_{\mathcal{F}_n} \tau(\llbracket \nabla u_n \rrbracket \cdot \mathbf{n})(\llbracket \nabla v_n \rrbracket \cdot \mathbf{n}) \quad \forall u_n, v_n \in V_n. \quad (43) $$

With an abuse of notation, we denote the norm induced by the bilinear form in (43) by $||\cdot||_{dG}$.

The bilinear form $B_n(\cdot, \cdot)$ is equal to that in (5) without the scalar jump terms. Moreover, the jump of the gradients appear only along the normal direction to the faces. Thus, *p*-optimal error estimates follow as (and are in fact easier to prove than) in Section 4.

For completeness, we report the convergence result.

**Theorem 5.1.** Let $\mathcal{T}_n$ be either a quasi-uniform triangular mesh, or a quasi-uniform tensor product-type mesh of either parallelograms (2D) or parallelepipeds (3D) as in Section 2.1. Let $u$ and $u_n$ be the solutions to (3) and the corresponding $C^0$-IPDG, respectively. Assume that $u \in H^{s+2}(\Omega)$, $s := \min(k, p-1)$, $p \ge 2$. Then, there exists a positive constant $c$ independent of $h$, and $p$ such that

$$ \|u - u_n\|_{dG} \le c \left(\frac{h}{p}\right)^s \|u\|_{s+2, \Omega}, \qquad \|u - u_n\|_{0, \Omega} \le c \left(\frac{h}{p}\right)^2 \|u - u_n\|_{dG}. $$
## 5.2 Divergence-conforming DG methods for the stream formulation of the Stokes problem in two dimensions

In this section, we discuss how to extend the optimal results of Section 4 to divergence-conforming DG methods for the stream formulation of the Stokes problem in two dimensions; see [25, 37].

We start with the standard formulation of the Stokes problem in two dimensions: given a polygonal domain $\Omega \subset \mathbb{R}^2$ and $f \in [L^2(\Omega)]^2$, we seek $\tilde{u} \in [H_0^1(\Omega)]^2$ and $\rho \in L_0^2(\Omega)$ that solve weakly

$$ \left\{ \begin{array}{ll} -\Delta \tilde{u} + \nabla \rho = \mathbf{f} & \text{in } \Omega, \\ \nabla \cdot \tilde{u} = 0 & \text{in } \Omega, \\ \tilde{u} = \mathbf{0} & \text{on } \partial\Omega. \end{array} \right. \qquad (44) $$

The velocity $\tilde{u}$ is divergence free. So, there exists $\psi \in H_0^2(\Omega)$ such that $\tilde{u} = \nabla \times \psi$. Introduce the Stokes kernel

$$ \mathbf{Z} := \{ \mathbf{v} \in [H_0^1(\Omega)]^2 \mid \operatorname{div} \mathbf{v} = 0 \}. $$

Testing the first equation in (44) with $\tilde{\mathbf{v}} \in \mathbf{Z}$ gives

$$ (\nabla \tilde{\mathbf{u}}, \nabla \tilde{\mathbf{v}})_{0,\Omega} = (\nabla \tilde{\mathbf{u}}, \nabla \tilde{\mathbf{v}})_{0,\Omega} + (\nabla \rho, \tilde{\mathbf{v}})_{0,\Omega} = (\mathbf{f}, \tilde{\mathbf{v}})_{0,\Omega}. $$

Since $\tilde{\mathbf{v}} \in \mathbf{Z}$, there exists $\phi \in H_0^2(\Omega)$ such that $\tilde{\mathbf{v}} = \nabla \times \phi$. We deduce the stream formulation of the Stokes problem in (44):

$$ (\nabla (\nabla \times \psi), \nabla (\nabla \times \phi))_{0,\Omega} = (\mathbf{f}, \nabla \times \phi)_{0,\Omega} \quad \forall \phi \in H_0^2(\Omega). \qquad (45) $$

This is a weak formulation of the biharmonic problem (2) which is alternative to that in (3). We can rewrite the right-hand side of (45) integrating by parts:

$$ (\mathbf{f}, \nabla \times \phi)_{0,\Omega} = (\nabla \times \mathbf{f}, \phi)_{0,\Omega}. $$

The strong formulation corresponding to (45) reads as follows: we seek $\psi \in H_0^2(\Omega)$ satisfying

$$ \Delta^2 \psi = \nabla \times \mathbf{f} \text{ in } \Omega. $$

We derive the $C^0$-IPDG formulation of this problem using standard dG techniques; see, e.g., [25]. On piecewise polynomial spaces, we define the discrete bilinear form

$$ B_n(\psi_n, \phi_n) := \int_{\mathcal{T}_n} (\nabla \times \nabla \psi_n) : (\nabla \times \nabla \phi_n) - \int_{\mathcal{F}_n} \{\nabla \times \nabla \psi_n\} \cdot \llbracket \nabla \phi_n \times \mathbf{n} \rrbracket \\ - \int_{\mathcal{F}_n} \{\nabla \times \nabla \phi_n\} \cdot \llbracket \nabla \psi_n \times \mathbf{n} \rrbracket + \int_{\mathcal{F}_n} \boldsymbol{\tau} \llbracket \nabla \psi_n \times \mathbf{n} \rrbracket \cdot \llbracket \nabla \phi_n \times \mathbf{n} \rrbracket \quad \forall \psi_n, \phi_n \in V_n. $$

With an abuse of notation, we denote the norm induced by the bilinear form in (43) by $\|\cdot\|_{dG}$.

Following the same steps as in Section 4, hp-optimal error estimates can be proved. In fact, the problem is the bilinear form appearing in $C^0$-IPDG (43) with rotated components.

For completeness, we report the convergence result.

**Theorem 5.2.** Let $\mathcal{T}_n$ be either a quasi-uniform triangular mesh, or a quasi-uniform tensor product-type mesh of either parallelograms (2D) or parallelepipeds (3D) as in Section 2.1. Let $\psi$ and $\psi_n$ be the solutions to (45) and the corresponding $C^0$-IPDG, respectively. Assume that $\psi \in H^{s+2}(\Omega)$, $s := \min(k, p-1)$, $p \ge 2$. Then, there exists a positive constant $c$ independent of $h$, and $p$ such that

$$ \|\psi - \psi_n\|_{dG} \le c \left(\frac{h}{p}\right)^s \|\psi\|_{s+2,\Omega}, \qquad \|\psi - \psi_n\|_{0,\Omega} \le c \left(\frac{h}{p}\right)^2 \|\psi - \psi_n\|_{dG}. $$
# 6 Numerical experiments

In this section, we validate the theoretical results of Section 4. More precisely, in Section 6.1, we show *p*-optimal convergence for a test case with internal point singularity lying inside one element; in Section 6.2, we test the convergence for a test case with internal point singularity lying on the skeleton and observe a *p*-optimal doubling of the convergence rate [4]; in Section 6.3, we consider a test case that is singular on the boundary of the domain, i.e., a case not covered in the foregoing analysis, and show a *p*-suboptimal rate of convergence.

<figure><img src="image_2.png" /><figcaption>**Figure 1:** The three different point singularities we shall consider. Left panel: point singularity inside one element. Central panel: point singularity at an internal mesh vertex. Right panel: point singularity on the boundary.</figcaption></figure>

Below, we consider the square domain $\Omega := (-1, 1)^2$ and a given uniform Cartesian mesh of four elements in Sections 6.1, 6.2, and 6.3, and the L-shape domain $\Omega := (-1, 1)^2 \setminus [0, 1) \times (-1, 0]$ and a given uniform triangular mesh of 6 elements in Section 6.4.

## 6.1 Point singularity inside one element

We consider the exact solution

$$u_1(x, y) := ((x - 0.5)^2 + (y - 0.5)^2)^{\frac{3}{2}} (1 - x^2)(1 - y^2).$$

The function $u_1$ has a point singularity at $(0.5, 0.5)$ that lies in the interior of an element of the mesh; see the left panel in Figure 1. It belongs to $H^{4-\epsilon}(\Omega)$ for any arbitrarily small $\epsilon > 0$. We list the convergence rates under *p*-refinement in Table 1. The results confirm the optimal error decay $O(p^{-2})$ for the DG-norm error.

<table><caption>**Table 1** *p*-convergence for the exact solution $u_1$.</caption><thead><tr><th>p</th><th>||u - u<sub>n</sub>||<sub>dG</sub></th><th>Convergence rate</th><th>p</th><th>||u - u<sub>n</sub>||<sub>dG</sub></th><th>Convergence rate</th></tr></thead><tbody><tr><td>2</td><td>1.46E+01</td><td></td><td>3</td><td>1.18E+01</td><td></td></tr><tr><td>4</td><td>5.07E+00</td><td>1.52</td><td>5</td><td>9.67E-01</td><td>4.90</td></tr><tr><td>6</td><td>1.42E-01</td><td>8.83</td><td>7</td><td>4.19E-02</td><td>9.33</td></tr><tr><td>8</td><td>2.38E-02</td><td>6.18</td><td>9</td><td>1.62E-02</td><td>3.79</td></tr><tr><td>10</td><td>1.11E-02</td><td>3.41</td><td>11</td><td>9.29E-03</td><td>2.77</td></tr><tr><td>12</td><td>6.76E-03</td><td>2.73</td><td>13</td><td>6.12E-03</td><td>2.49</td></tr><tr><td>14</td><td>4.67E-03</td><td>2.40</td><td>15</td><td>4.39E-03</td><td>2.33</td></tr><tr><td>16</td><td>3.48E-03</td><td>2.21</td><td>17</td><td>3.33E-03</td><td>2.20</td></tr><tr><td>18</td><td>2.73E-03</td><td>2.05</td><td>19</td><td>3.33E-03</td><td>2.09</td></tr><tr><td>20</td><td>2.23E-03</td><td>1.93</td><td>21</td><td>2.16E-03</td><td>1.99</td></tr><tr><td>22</td><td>1.84E-03</td><td>2.03</td><td>23</td><td>1.80E-03</td><td>2.01</td></tr><tr><td>24</td><td>1.55E-03</td><td>1.96</td><td>25</td><td>1.53E-03</td><td>1.92</td></tr></tbody></table>

## 6.2 Point singularity at an internal mesh vertex

We consider the exact solution

$$u_2(x, y) := (x^2 + y^2)^{\frac{3}{2}} (1 - x^2)(1 - y^2).$$
The function $u_2$ has a point singularity at $(0,0)$, which lies at an internal vertex of the mesh; see the central panel in Figure 1. It belongs to $H^{4-\epsilon}(\Omega)$ for any arbitrarily small $\epsilon > 0$. We list the convergence rates under $p$-refinement in Table 2 and observe a doubling rate of the error decay $O(p^{-4})$ for the DG-norm error. This phenomenon can be proven based on the optimal results in Section 4 and the conforming finite element arguments in [4, 36].

<table><caption>Table 2 $p$-convergence for the exact solution $u_2$.</caption><thead><tr><th>p</th><th>||u - u<sub>n</sub>||<sub>dG</sub></th><th>Convergence rate</th><th>p</th><th>||u - u<sub>n</sub>||<sub>dG</sub></th><th>Convergence rate</th></tr></thead><tbody><tr><td>2</td><td>5.16E+00</td><td></td><td>3</td><td>4.39E+00</td><td></td></tr><tr><td>4</td><td>2.98E+00</td><td>0.79</td><td>5</td><td>9.61E-01</td><td>2.97</td></tr><tr><td>6</td><td>8.64E-02</td><td>8.74</td><td>7</td><td>1.08E-02</td><td>13.33</td></tr><tr><td>8</td><td>6.22E-03</td><td>9.14</td><td>9</td><td>3.84E-03</td><td>4.14</td></tr><tr><td>10</td><td>2.47E-03</td><td>4.14</td><td>11</td><td>1.67E-03</td><td>4.13</td></tr><tr><td>12</td><td>1.18E-03</td><td>4.07</td><td>13</td><td>8.52E-04</td><td>4.04</td></tr><tr><td>14</td><td>6.34E-04</td><td>3.99</td><td>15</td><td>4.84E-04</td><td>3.96</td></tr><tr><td>16</td><td>3.76E-04</td><td>3.92</td><td>17</td><td>2.90E-04</td><td>4.09</td></tr><tr><td>18</td><td>2.31E-04</td><td>4.12</td><td>19</td><td>1.87E-04</td><td>3.92</td></tr><tr><td>20</td><td>1.53E-04</td><td>3.90</td><td>21</td><td>1.27E-04</td><td>3.86</td></tr><tr><td>22</td><td>1.05E-04</td><td>3.95</td><td>23</td><td>8.82E-05</td><td>4.01</td></tr><tr><td>24</td><td>7.44E-05</td><td>3.99</td><td>25</td><td>6.28E-05</td><td>4.06</td></tr></tbody></table>

## 6.3 Point singularity on the boundary for the square domain

We consider the exact solution

$$u_3(x, y) := ((x - 1)^2 + (y - 0.5)^2)^{\frac{3}{2}}.$$

The function $u_3$ has a point singularity at $(1,0.5)$, which lies on the boundary of $\Omega$; see the right panel in Figure 1). In particular, it belongs to $H^{4-\epsilon}(\Omega)$ for any arbitrarily small $\epsilon > 0$. We list the convergence rates under $p$-refinement in Table 3 and observe a suboptimal by 1.5 order rate $O(p^{-2/3})$ for the DG-norm error. This suboptimality is the counterpart of what is observed in [19] for elliptic partial differential equations of second order.

<table><caption>Table 3 $p$-convergence for the exact solution $u_3$.</caption><thead><tr><th>p</th><th>||u - u<sub>n</sub>||<sub>dG</sub></th><th>Convergence rate</th><th>p</th><th>||u - u<sub>n</sub>||<sub>dG</sub></th><th>Convergence rate</th></tr></thead><tbody><tr><td>2</td><td>8.61E+00</td><td></td><td>3</td><td>2.57E+00</td><td></td></tr><tr><td>4</td><td>6.39E-01</td><td>3.75</td><td>5</td><td>1.01E+00</td><td>1.83</td></tr><tr><td>6</td><td>4.90E-01</td><td>0.65</td><td>7</td><td>7.48E-01</td><td>0.89</td></tr><tr><td>8</td><td>4.33E-01</td><td>0.44</td><td>9</td><td>6.37E-01</td><td>0.64</td></tr><tr><td>10</td><td>3.91E-01</td><td>0.46</td><td>11</td><td>5.71E-01</td><td>0.55</td></tr><tr><td>12</td><td>3.56E-01</td><td>0.52</td><td>13</td><td>5.24E-01</td><td>0.51</td></tr><tr><td>14</td><td>3.31E-01</td><td>0.46</td><td>15</td><td>4.88E-01</td><td>0.50</td></tr><tr><td>16</td><td>3.13E-01</td><td>0.42</td><td>17</td><td>4.59E-01</td><td>0.49</td></tr><tr><td>18</td><td>2.97E-01</td><td>0.47</td><td>19</td><td>4.36E-01</td><td>0.48</td></tr><tr><td>20</td><td>2.82E-01</td><td>0.48</td><td>21</td><td>4.15E-01</td><td>0.48</td></tr><tr><td>22</td><td>2.69E-01</td><td>0.49</td><td>23</td><td>3.97E-01</td><td>0.49</td></tr><tr><td>24</td><td>2.58E-01</td><td>0.49</td><td>25</td><td>3.81E-01</td><td>0.49</td></tr></tbody></table>

## 6.4 Point singularity on the boundary for the L-shape domain

On the L-shape domain, we consider the exact solution

$$u_4(x, y) = u_4(r, \theta) := r^{\frac{5}{3}} \sin\left(\frac{5}{3}\theta\right),$$
where $(r, \theta)$ are the polar coordinates at the re-entrant corner (0,0). We consider the triangular mesh in Figure 2

<figure><img src="image_3.png" /><figcaption>Figure 2: Uniform triangular mesh on the L-shape domain.</figcaption></figure>

The function $u_4$ has a point singularity at (0,0) and belongs to $H^{\frac{5}{3}-\varepsilon}(\Omega)$ for any arbitrarily small $\varepsilon > 0$; for this reason, optimal convergence estimates are not guaranteed by the foregoing results, as extra regularity on the solution is required. Despite this, we exhibit the usual doubling rate of convergence, namely the decay rate for the DG-norm error is $O(p^{\frac{4}{3}})$ as similarly observed for the square domain in Section 6.3; see Table 4.

<table><caption>Table 4 p-convergence for the exact solution $u_1$.</caption><thead><tr><th>p</th><th>||u - u<sub>n</sub>||<sub>dG</sub></th><th>Convergence rate</th><th>p</th><th>||u - u<sub>n</sub>||<sub>dG</sub></th><th>Convergence rate</th></tr></thead><tbody><tr><td>2</td><td>1.80E+00</td><td></td><td>3</td><td>9.38E-01</td><td></td></tr><tr><td>4</td><td>6.93E-01</td><td>1.82</td><td>5</td><td>5.50E-01</td><td>1.35</td></tr><tr><td>6</td><td>4.53E-01</td><td>1.27</td><td>7</td><td>3.84E-01</td><td>1.26</td></tr><tr><td>8</td><td>3.32E-01</td><td>1.25</td><td>9</td><td>2.91E-01</td><td>1.25</td></tr><tr><td>10</td><td>2.58E-01</td><td>1.26</td><td>11</td><td>2.32E-01</td><td>1.26</td></tr><tr><td>12</td><td>2.09E-01</td><td>1.27</td><td>13</td><td>1.91E-01</td><td>1.27</td></tr><tr><td>14</td><td>1.75E-01</td><td>1.27</td><td>15</td><td>1.62E-01</td><td>1.28</td></tr><tr><td>16</td><td>1.50E-01</td><td>1.28</td><td>17</td><td>1.39E-01</td><td>1.29</td></tr><tr><td>18</td><td>1.30E-01</td><td>1.29</td><td>19</td><td>1.22E-01</td><td>1.28</td></tr><tr><td>20</td><td>1.15E-01</td><td>1.29</td><td>21</td><td>1.08E-01</td><td>1.30</td></tr><tr><td>22</td><td>1.02E-01</td><td>1.31</td><td>23</td><td>9.60E-02</td><td>1.31</td></tr><tr><td>24</td><td>9.22E-02</td><td>1.31</td><td>25</td><td>8.77E-02</td><td>1.32</td></tr></tbody></table>

# 7 Conclusions

We proved that the IPDG for the biharmonic problem is *hp*-optimal if essential boundary conditions are homogeneous, and triangular and tensor product-type meshes (in two and three dimensions) are employed. The analysis hinges upon proving the existence of a global $H^2$ piecewise polynomial over the given meshes with *hp*-optimal approximation properties. The results discussed for the biharmonic problem were extended to other methods, e.g., $C^0$-IPDG, and problems, e.g., the stream formulation of the Stokes problem. The theoretical results were validated on several test cases, which also showed *p*-suboptimality for solutions with singular essential boundary conditions.

**Competing Interests** The authors have no relevant financial or non-financial interests to disclose.

**Data Availability.** The datasets generated during and/or analysed during the current study are available on request.
References

[1] M. Ainsworth and C. Parker. $H^2$-stable polynomial liftings on triangles. *SIAM J. Numer. Anal.*, 58(3):1867–1892, 2020.

[2] M. Ainsworth and C. Parker. Preconditioning high order $H^2$ conforming finite elements on triangles. *Numer. Math.*, 148(2):223–254, 2021.

[3] J. H. Argyris, I. Fried, and D. W. Scharpf. The TUBA family of plate elements for the matrix displacement method. *Aeronaut. J.*, 72(692):701–709, 1968.

[4] 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.

[5] G. A. Baker. Finite element methods for elliptic equations using nonconforming elements. *Math. Comp.*, 31(137):45–59, 1977.

[6] L. Beirão da Veiga, A. Buffa, J. Rivas, and G. Sangalli. Some estimates for $h-p-k$-refinement in isogeometric analysis. *Numer. Math.*, 118(2):271–305, 2011.

[7] L. Beirão da Veiga, J. Niiranen, and R. Stenberg. A family of $C^0$ finite elements for Kirchhoff plates I: Error analysis. *SIAM J. Numer. Anal.*, 45(5):2047–2071, 2007.

[8] L. Beirão da Veiga, J. Niiranen, and R. Stenberg. A posteriori error analysis for the Morley plate element with general boundary conditions. *Internat. J. Numer. Methods Engrg.*, 83(1):1–26, 2010.

[9] S. C. Brenner, T. Gudi, and L.-Y. Sung. An a posteriori error estimator for a quadratic $C^0$-interior penalty method for the biharmonic problem. *IMA J. Numer. Anal.*, 30(3):777–798, 2010.

[10] S. C. Brenner and L. R. Scott. *The mathematical theory of Finite Element Methods*, volume 15. Texts in Applied Mathematics, Springer-Verlag, New York, third edition, 2008.

[11] S. C. Brenner and L.-Y. Sung. $C^0$ interior penalty methods for fourth order elliptic boundary value problems on polygonal domains. *J. Sci. Comp.*, 22(1):83–118, 2005.

[12] S. C. Brenner, K. Wang, and J. Zhao. Poincaré-Friedrichs inequalities for piecewise $H^2$ functions. *Numer. Funct. Anal. Optim.*, 25(5-6):463–478, 2004.

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

[14] Z. Dong. Discontinuous Galerkin methods for the biharmonic problem on polygonal and polyhedral meshes. *Int. J. Numer. Anal. Model.*, 16(5):825–846, 2019.

[15] Z. Dong, L. Mascotto, and O. J. Sutton. Residual-based a posteriori error estimates for hp-discontinuous Galerkin discretizations of the biharmonic problem. *SIAM J. Numer. Anal.*, 59(3):1273–1298, 2021.

[16] J. Douglas Jr., T. Dupont, P. Percell, and R. Scott. A family of finite elements with optimal approximation properties for various Galerkin methods for 2nd and 4th order problems. *RAIRO. Anal. numér.*, 13(3):227–255, 1979.

[17] G. Engel, K. Garikipati, T. J. R. Hughes, M. G. Larson, L. Mazzei, and R. L. Taylor. Continuous/discontinuous finite element approximations of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity. *Comput. Methods Appl. Mech. Engrg.*, 191(34):3669–3750, 2002.

[18] E. H. Georgoulis. *Discontinuous Galerkin methods on shape-regular and anisotropic meshes*. PhD thesis, University of Oxford, 2003.

[19] E. H. Georgoulis, E. Hall, and J. M. Melenk. On the suboptimality of the $p$-version interior penalty discontinuous Galerkin method. *J. Sci. Comput.*, 42(1):54–67, 2010.

[20] E. H. Georgoulis and P. Houston. Discontinuous Galerkin methods for the biharmonic problem. *IMA J. Numer. Anal.*, 29(3):573–594, 2009.

[21] E. H. Georgoulis and A. Lasis. A note on the design of hp-version interior penalty discontinuous Galerkin finite element methods for degenerate problems. *IMA J. Numer. Anal.*, 26(2):381–390, 2006.

[22] P. Grisvard. *Elliptic problems in nonsmooth domains*. SIAM, 2011.

[23] T. Gudi, N. Nataraj, and A. K. Pani. Mixed discontinuous Galerkin finite element method for the biharmonic equation. *J. Sci. Comp.*, 37(2):139–161, 2008.

[24] 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.

[25] G. Kanschat and N. Sharma. Divergence-conforming discontinuous Galerkin methods and $C^0$ interior penalty methods. *SIAM J. Numer. Anal.*, 52(4):1822–1842, 2014.

[26] V. A. Kozlov and V. G. Mazya. Singularities in solutions to mathematical physics problems in non-smooth domains. In *Partial Differential Equations and Functional Analysis*, pages 174–206. Springer, 1996.

[27] P. L. Lederer and J. Schöberl. Polynomial robust stability analysis for $H(\mathrm{div})$-conforming finite elements for the Stokes equations. *IMA J. Numer. Anal.*, 38(4):1832–1860, 2018.

[28] L. S. D. Morley. The triangular equilibrium element in the solution of plate bending problems. *Aero. Quart.*, 19:149–169, 1968.
[29] I. Mozolevski and E. Süli. A priori error analysis for the hp-version of the discontinuous Galerkin finite element method for the biharmonic equation. *Comput. Methods Appl. Math.*, 3(4):596–607, 2003.

[30] I. Mozolevski, E. Süli, and P. R. Bösing. hp-version a priori error analysis of interior penalty discontinuous Galerkin finite element approximations to the biharmonic equation. *J. Sci. Comp.*, 30(3):465–491, 2007.

[31] I. Perugia and D. Schötzau. An hp-analysis of the local discontinuous Galerkin method for diffusion problems. *J. Sci. Comp.*, 17(1):561–571, 2002.

[32] Ch. Schwab. *p- and hp- Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics*. Clarendon Press Oxford, 1998.

[33] B. Stamm and T. Wihler. hp-optimal discontinuous Galerkin methods for linear elliptic problems. *Math. Comp.*, 79(272):2117–2133, 2010.

[34] E. M. Stein. *Singular integrals and differentiability properties of functions*, volume 2. Princeton University Press, 1970.

[35] E. Süli and I. Mozolevski. hp-version interior penalty DGFEMs for the biharmonic equation. *Comput. Methods Appl. Mech. Engrg.*, 196(13-16):1851–1863, 2007.

[36] M. Suri. The p-version of the finite element method for elliptic equations of order 2l. *RAIRO Modél. Math. Anal. Numér.*, 24(2):265–304, 1990.

[37] J. Wang, Y. Wang, and X. Ye. A robust numerical method for Stokes equations based on divergence-free H(div) finite element methods. *SIAM J. Sci. Comput.*, 31(4):2784–2802, 2009.

[38] M. Wang and J. Xu. The Morley element for fourth order elliptic equations in any dimensions. *Numer. Math.*, 103(1):155–169, 2006.

## A Proof of Lemma 3.3

For the sake of presentation, throughout the proof, we shall write $u_p$ instead of $\mathcal{P}_p^1 u$.

We set $u_p''$ as the Legendre expansion of $u''$ truncated at order $p-2$. In other words, expanding $u''$ as a series of Legendre polynomials $L_j(\xi)$, $j \in \mathbb{N}_0$,

$$u''(\xi) = \sum_{i=0}^{+\infty} b_j L_j(\xi), \qquad (46)$$

we define

$$u_p''(\xi) := \sum_{i=0}^{p-2} b_j L_j(\xi). \qquad (47)$$

Standard properties of Legendre polynomials, see, e.g., [32, Appendix C], imply

$$b_j = \frac{2j+1}{2} \int_{-1}^{1} u''(\xi) L_j(\xi) d\xi. \qquad (48)$$

For $k \ge 0$, we recall [32, Lemma 3.10] that

$$\sum_{i=k}^{+\infty} \frac{2}{2i+1} \frac{(i+k)!}{(i-k)!} |b_i|^2 = \int_{\hat{I}} (1-\xi^2)^k |u^{(k+2)}(\xi)|^2 d\xi \le \|u^{(k+2)}\|_{0,\hat{I}}^2. \qquad (49)$$

Using the orthogonality properties of Legendre polynomials [32, eq. (C.24)], the fact that $s \le p-2$, and (49), we obtain

$$\begin{aligned} \|u'' - u_p''\|_{0,\hat{I}}^2 &= \sum_{i=p-1}^{+\infty} \frac{2}{2i+1} |b_i|^2 = \sum_{i=p-1}^{+\infty} \frac{2}{2i+1} \frac{(i+s)!}{(i-s)!} |b_i|^2 \frac{(i-s)!}{(i+s)!} \\ &\le \frac{(p-s-1)!}{(p+s-1)!} \sum_{i=s}^{+\infty} \frac{2}{2i+1} \frac{(i+s)!}{(i-s)!} |b_i|^2 \le \frac{(p-s-1)!}{(p+s-1)!} \|u^{(s+2)}\|_{0,\hat{I}}^2 \end{aligned}$$

The two above equations yield the first bound in (12).

Next, we introduce

$$u_p'(\xi) = \int_{-1}^{\xi} u_p''(t) dt + u'(-1). \qquad (50)$$
We have $u'_p(-1) = u'(-1)$. Moreover, recalling that $u''_p$ and $u''$ have the same average over $\hat{I}$, we also have

$$u'_{p}(1) = \int_{-1}^{1} u''_{p}(t)dt + u'(-1) = \int_{-1}^{1} u''(t)dt + u'(-1) = u'(1),$$

which proves (11) for the derivative of $u_p$ at the endpoints of $\hat{I}$.

At this point, we observe

$$u'(\xi) - u'_{p}(\xi) = \int_{-1}^{\xi} (u''(t) - u''_{p}(t))dt = \sum_{i=p-1}^{+\infty} b_i \int_{-1}^{\xi} L_i(t)dt =: \sum_{i=p-1}^{+\infty} b_i \phi_i(\xi). \quad (51)$$

Recall the Legendre differential equation [32, eq. (C.2.3)]

$$((1 - \xi^2)L'_i(\xi))' + i(i+1)L_i(\xi) = 0 \quad \forall i \in \mathbb{N}_0.$$

Integrating the above identity over $(-1, \xi)$, $\xi \in (-1, 1)$, yields

$$\phi_i(\xi) = -\frac{1}{i(i+1)}(1 - \xi^2)L'_i(\xi). \quad (52)$$

Recall the orthogonality property of the derivatives of Legendre polynomials [32, eq. (3.39)]:

$$\int_{-1}^{1} (1 - \xi^2)L'_{i}(\xi)L'_{j}(\xi)d\xi = \frac{2\delta_{i,j}}{2i+1} \frac{(i+1)!}{(i-1)!}. \quad (53)$$

Combining (52) and (53), for all $i, j \in \mathbb{N}_0$, we deduce

$$\begin{aligned} \int_{-1}^{1} (1 - \xi^2)^{-1} \phi_i(\xi) \phi_j(\xi) d\xi &= \frac{1}{i(i+1)j(j+1)} \int_{-1}^{1} (1 - \xi^2) L'_i(\xi) L'_j(\xi) d\xi \\ &= \frac{1}{i^2(i+1)^2} \frac{2\delta_{i,j}}{2i+1} \frac{(i+1)!}{(i-1)!} = \frac{2\delta_{i,j}}{i(i+1)(2i+1)}. \end{aligned} \quad (54)$$

Using $s \le p-2$, we write

$$\begin{aligned} \int_{-1}^{1} |u'(\xi) - u'_{p}(\xi)|^2 d\xi &= \int_{-1}^{1} \left| \sum_{i=p-1}^{+\infty} b_i \phi_i(\xi) \right|^2 d\xi \le \int_{-1}^{1} (1 - \xi^2)^{-1} \left| \sum_{i=p-1}^{+\infty} b_i \phi_i(\xi) \right|^2 d\xi \\ &\stackrel{(54)}{\le} \sum_{i=p-1}^{+\infty} \frac{2|b_i|^2}{i(i+1)(2i+1)} = \sum_{i=p-1}^{+\infty} \left( \frac{2|b_i|^2}{2i+1} \frac{(i+s)!}{(i-s)!} \right) \frac{(i-s)!}{(i+s)!} \frac{1}{i(i+1)} \\ &\stackrel{(49)}{\le} \frac{(p-s-1)!}{(p+s-1)!} \frac{1}{(p-1)p} \left( \sum_{i=s}^{+\infty} \frac{2|b_i|^2}{2i+1} \frac{(i+s)!}{(i-s)!} \right) \stackrel{(52)}{\le} \frac{(p-s-1)!}{(p+s-1)!} \frac{1}{(p-1)p} \|u^{(s+2)}\|_{0, \hat{I}}^2 \end{aligned}$$

which is the second bound in (12).

Finally, we introduce

$$u_p(\xi) = \int_{-1}^{\xi} u'_{p}(t)dt + u(-1). \quad (55)$$

We observe that $u_p(-1) = u(-1)$. Since $L'_1(t) = 1 = L_0(t)$, standard manipulations imply

$$\begin{aligned} u_p(1) - u_p(-1) &= \int_{-1}^{1} u'_{p}(t)dt \stackrel{(50)}{=} 2u'(-1) + \int_{-1}^{1} \left( \int_{-1}^{t} u''_{p}(x)dx \right) dt \\ &\stackrel{(46),(51)}{=} 2u'(-1) + \int_{-1}^{1} \sum_{i=0}^{p-2} b_i \phi_i(t)dt = 2u'(-1) + \sum_{i=0}^{p-2} b_i \int_{-1}^{1} \phi_i(t)dt \\ &\stackrel{(52)}{=} 2u'(-1) - \sum_{i=0}^{p-1} b_i \int_{-1}^{1} \frac{1}{i(i+1)}(1-t^2)L'_i(t)L'_1(t)dt \\ &\stackrel{(53)}{=} 2u'(-1) - \sum_{i=0}^{+\infty} b_i \int_{-1}^{1} \frac{1}{i(i+1)}(1-t^2)L'_i(t)L'_1(t)dt \\ &\stackrel{(52),(46)}{=} 2u'(-1) + \int_{-1}^{1} \left( \int_{-1}^{t} u''(x)dx \right) dt = \int_{-1}^{1} u'(t)dt = u(1) - u(-1). \end{aligned}$$
Using that $u_p(-1) = u(-1)$, we deduce $u_p(1) = u(1)$.

We are left with proving error estimates in the $L^2$ norm. To this aim, observe that

$$u(\xi) - u_p(\xi) = \int_{-1}^{\xi} (u'(t) - u'_p(t))dt.$$

We arrive at

$$u(\xi) - u_p(\xi) \stackrel{(51)}{=} \sum_{i=p-1}^{+\infty} b_i \int_{-1}^{\xi} \phi_i(x) \, dx \, dt := \sum_{i=p-1}^{+\infty} b_i \psi_i(\xi). \quad (56)$$

We prove certain orthogonality properties of the $\psi_i$ functions. Recall the identity [32, eq. (C.2.5)]:

$$L_i(\xi) = \frac{L'_{i+1}(\xi) - L'_{i-1}(\xi)}{2i+1} \quad \forall i \ge 2, \qquad L_0(\xi) = L'_1(\xi).$$

Integrating over $(-1, t)$, $t \in (-1, 1)$, and using $L_{i+1}(-1) = L_{i-1}(-1)$, see [32, eq. (C.2.6)], we deduce

$$\phi_i(t) = \frac{L_{i+1}(\xi) - L_{i-1}(\xi)}{2i+1}.$$

Upon integrating the above identity over $(-1, \xi)$, $\xi \in (-1, 1)$, we arrive at

$$\psi_i(\xi) = \frac{\int_{-1}^{\xi} L_{i+1}(t) \, dt - \int_{-1}^{\xi} L_{i-1}(t) \, dt}{2i+1} \stackrel{(51)}{=} \frac{\phi_{i+1}(\xi) - \phi_{i-1}(\xi)}{2i+1}.$$

With the notation $\phi_{-1}(\xi) = 0$, we get

$$\psi_i(\xi)\psi_j(\xi) = \frac{1}{(2i+1)(2j+1)} (\phi_{i+1}(\xi)\phi_{j+1}(\xi) + \phi_{i-1}(\xi)\phi_{j-1}(\xi) - \phi_{i+1}(\xi)\phi_{j-1}(\xi) - \phi_{i-1}(\xi)\phi_{j+1}(\xi)).$$

In (54), we proved that the $\phi_i$ functions are orthogonal with respect to the $(1 - \xi^2)$-weighted $L^2$ inner product. Therefore, testing the above identity by $(1 - \xi^2)^{-1}$ and integrating over $(-1, 1)$ we arrive at

$$\int_{-1}^{1} (1 - \xi^2)^{-1} \psi_i(\xi) \psi_j(\xi) \, d\xi = \frac{1}{(2i+1)(2j+1)} (\aleph_{i+1,j+1} + \aleph_{i-1,j-1} - \aleph_{i+1,j-1} - \aleph_{i-1,j+1}), \quad (57)$$

where, from (54),

$$\aleph_{\ell,k} = \begin{cases} 0 & \text{if } \ell \le 0 \text{ or } k \le 0, \\ \frac{2\delta_{\ell,k}}{\ell(\ell+1)(2\ell+1)} & \text{if } \ell, k > 0. \end{cases}$$

We provide explicit values for the expression in (57). If $j \in \{i-2, i, i+2\}$, then

$$\int_{-1}^{1} (1 - \xi^2)^{-1} \psi_i(\xi) \psi_j(\xi) \, d\xi = 0.$$

If $j = i$, then

$$\int_{-1}^{1} (1 - \xi^2)^{-1} \psi_i(\xi) \psi_i(\xi) \, d\xi = \mathfrak{Z}_1^i + \mathfrak{Z}_2^i,$$

where

$$\mathfrak{Z}_1^i := \frac{1}{(2i+1)^2} \frac{2}{(i+1)(i+2)(2i+3)} \quad \forall i \in \mathbb{N}; \qquad \mathfrak{Z}_2^i := \begin{cases} 0 & \text{if } i=1, \\ \frac{1}{(2i+1)^2} \frac{2}{(i-1)i(2i-1)} & \text{if } i \ge 2. \end{cases} \quad (58)$$

If $j = i + 2$, then

$$\int_{-1}^{1} (1 - \xi^2)^{-1} \psi_i(\xi) \psi_{i+2}(\xi) \, d\xi = - \frac{1}{(2i+1)(2i+5)} \cdot \frac{2}{(i+1)(i+2)(2i+3)} = - \mathfrak{Z}_3^i. \quad (59)$$
If $j = i - 2$, then

$$ \int_{-1}^{1} (1 - \xi^2)^{-1} \psi_i(\xi) \psi_{i-2}(\xi) \, d\xi = -\mathbb{Z}_4^i := \begin{cases} -\frac{1}{(2i+1)(2i-3)} \cdot \frac{2}{(i-1)i(2i-1)} & \text{if } i \ge 3, \\ 0 & \text{if } i \le 2. \end{cases} \quad (60) $$

In light of the above orthogonality properties, we can estimate the $L^2$ approximation error as follows:

$$ \begin{aligned} \|u - u_p\|_{0, \hat{I}}^2 &\stackrel{(56)}{=} \left\| \sum_{i=p-1}^{+\infty} b_i \psi_i \right\|_{0, \hat{I}}^2 \le \sum_{i,j=p-1}^{+\infty} \left| b_i \, b_j \int_{-1}^{1} (1 - \xi^2)^{-1} \psi_i(\xi) \psi_j(\xi) \, d\xi \right| \\ &\stackrel{(58),(59),(60)}{\le} \sum_{i=p-1}^{+\infty} |b_i|^2 (\mathbb{Z}_1^i + \mathbb{Z}_2^i + \mathbb{Z}_3^i + \mathbb{Z}_4^i). \end{aligned} $$

We cope with the four terms on the right-hand side separately. We begin with $\mathbb{Z}_1^i$:

$$ \begin{aligned} \sum_{i=p-1}^{+\infty} |b_i|^2 \mathbb{Z}_1^i &\stackrel{(58)}{=} \sum_{i=p-1}^{+\infty} |b_i|^2 \frac{1}{(2i+1)^2} \cdot \frac{2}{(i+1)(i+2)(2i+3)} \\ &= \sum_{i=p-1}^{+\infty} \left( \frac{2}{2i+1} |b_i|^2 \frac{(i+s)!}{(i-s)!} \right) \left( \frac{(i-s)!}{(i+s)!} \frac{1}{(i+1)(i+2)(2i+1)(2i+3)} \right) \\ &\stackrel{(49)}{\le} \frac{(p-s-1)!}{(p+s-1)!} \frac{1}{p(p+1)(2p-1)(2p+1)} \|u^{(s+2)}\|_{0, \hat{I}}^2. \end{aligned} $$

Next, we focus on $\mathbb{Z}_2^i$, $i \ge 2$:

$$ \begin{aligned} \sum_{i=\max(2,p-1)}^{+\infty} |b_i|^2 \mathbb{Z}_2^i &\stackrel{(58)}{=} \sum_{i=\max(2,p-1)}^{+\infty} |b_i|^2 \frac{1}{(2i+1)^2} \cdot \frac{2}{(i-1)i(2i-1)} \\ &= \sum_{i=\max(2,p-1)}^{+\infty} \left( \frac{2}{2i+1} |b_i|^2 \frac{(i+s)!}{(i-s)!} \right) \left( \frac{(i-s)!}{(i+s)!} \frac{1}{(i-1)i(2i-1)(2i+1)} \right) \\ &\stackrel{(49)}{\le} \frac{(p-s-1)!}{(p+s-1)!} \frac{1}{(p-2)(p-1)(2p-3)(2p-1)} \|u^{(s+2)}\|_{0, \hat{I}}^2. \end{aligned} $$

As for the term $\mathbb{Z}_3^i$, we proceed as follows:

$$ \begin{aligned} \sum_{i=p-1}^{+\infty} |b_i|^2 \mathbb{Z}_3^i &\stackrel{(59)}{=} \sum_{i=p-1}^{+\infty} |b_i|^2 \frac{1}{(2i+1)(2i+5)} \cdot \frac{2}{(i+1)(i+2)(2i+3)} \\ &= \sum_{i=p-1}^{+\infty} \left( \frac{2}{2i+1} |b_i|^2 \frac{(i+s)!}{(i-s)!} \right) \left( \frac{(i-s)!}{(i+s)!} \frac{1}{(i+1)(i+2)(2i+3)(2i+5)} \right) \\ &\stackrel{(49)}{\le} \frac{(p-s-1)!}{(p+s-1)!} \frac{1}{p(p+2)(2p+1)(2p+3)} \|u^{(s+2)}\|_{0, \hat{I}}^2. \end{aligned} $$

Eventually, we cope with the term $\mathbb{Z}_4^i$, $i \ge 3$:

$$ \begin{aligned} \sum_{i=\max(3,p-1)}^{+\infty} |b_i|^2 \mathbb{Z}_4^i &\stackrel{(60)}{=} \sum_{i=\max(3,p-1)}^{+\infty} |b_i|^2 \frac{1}{(2i+1)(2i-3)} \cdot \frac{2}{(i-1)i(2i-1)} \\ &= \sum_{i=p-1}^{+\infty} \left( \frac{2}{2i+1} |b_i|^2 \frac{(i+s)!}{(i-s)!} \right) \left( \frac{(i-s)!}{(i+s)!} \frac{1}{(i-1)i(2i-3)(2i-1)} \right) \\ &\stackrel{(49)}{\le} \frac{(p-s-1)!}{(p+s-1)!} \frac{1}{(p-2)(p-1)(2p-5)(2p-3)} \|u^{(s+2)}\|_{0, \hat{I}}^2. \end{aligned} $$

Collecting the five bounds above concludes the proof.
## B Proof of Theorem 3.6 (2D case)

The continuity properties (15) follow from the definition of the operator $P_p$. Therefore, we only focus on the error estimates. We split the proof into several steps. Recall that $P_p^x$ and $P_p^y$ are the projections in Lemma 3.3 along the $x$ and $y$ directions, respectively.

**Some identities.** The following identities are valid:

$$ \partial_y P_p^x u = P_p^x(\partial_y u), \quad \partial_x P_p^y u = P_p^y(\partial_x u). \qquad (61) $$

We only show the first one as the second can be proven analogously. For all $y \in (-1, 1)$, after writing $P_p^x u(x, y)$ in integral form with respect to the $x$ variable, we have

$$ \begin{aligned} P_p^x u(x, y) &= \int_{-1}^{x} \left[ \left( \int_{-1}^{t} \partial_x^2 P_p^x u(x, y) dx \right) + \partial_x u(-1, y) \right] dt + u(-1, y) \\ &\stackrel{(11)}{=} \int_{-1}^{x} \left[ \left( \int_{-1}^{t} \partial_x^2 P_p^x u(x, y) dx \right) + \partial_x u(-1, y) \right] dt + u(-1, y) \\ &\stackrel{(47),(48)}{=} \int_{-1}^{x} \left[ \left( \int_{-1}^{t} \partial_x^2 \left( \sum_{i=0}^{p-2} \frac{2i+1}{2} \int_{-1}^{1} u(s, y) L_i(s) ds \right) L_i(x) dx \right) + \partial_x u(-1, y) \right] dt + u(-1, y). \end{aligned} $$

Taking the partial derivative with respect to $y$ on both sides yields

$$ \begin{aligned} & \partial_y P_p^x u(x, y) \\ &= \int_{-1}^{x} \left[ \left( \int_{-1}^{t} \partial_x^2 \left( \sum_{i=0}^{p-2} \frac{2i+1}{2} \int_{-1}^{1} \partial_y u(s, y) L_i(s) ds \right) L_i(x) dx \right) + \partial_x \partial_y u(-1, y) \right] dt + \partial_y u(-1, y) \\ &\stackrel{(47),(50),(55)}{=} P_p^x \partial_y u(x, y). \end{aligned} $$

**$L^2$ estimates.** Using the definition of $P_p$, the triangle inequality, the one dimensional approximation properties (13), the third stability property in (14), and the identities (61) we write

$$ \begin{aligned} \|u - P_p u\|_{0, \hat{Q}} &\le \|u - P_p^x u\|_{0, \hat{Q}} + \|P_p^x (u - P_p^y u)\|_{0, \hat{Q}} \\ &\lesssim p^{-s-2} \|\partial_x^{s+2} u\|_{0, \hat{Q}} + \|u - P_p^y u\|_{0, \hat{Q}} + p^{-2} \|\partial_x^2 (u - P_p^y u)\|_{0, \hat{Q}} \\ &= p^{-s-2} \|\partial_x^{s+2} u\|_{0, \hat{Q}} + \|u - P_p^y u\|_{0, \hat{Q}} + p^{-2} \|\partial_x^2 u - P_p^y \partial_x^2 u\|_{0, \hat{Q}}. \end{aligned} $$

The assertion follows using again the one dimensional approximation properties (13).

**$H^1$ estimates.** First, we cope with the bound on the derivative with respect to $x$. Adding and subtracting $P_p^x u$, and using the triangle inequality yield

$$ \|\partial_x(u - P_p u)\|_{0, \hat{Q}} \le \|\partial_x(u - P_p^x u)\|_{0, \hat{Q}} + \|\partial_x P_p^x(u - P_p u)\|_{0, \hat{Q}} =: T_1 + T_2. \qquad (62) $$

As for the term $T_1$, we use the one dimensional approximation properties (13) and get

$$ T_1 \lesssim p^{-s-1} \|\partial_x^{s+2} u\|_{0, \hat{Q}}. \qquad (63) $$

As for the term $T_2$, we use the second stability estimate in (14) and get

$$ T_2 \lesssim \|\partial_x(u - P_p^y u)\|_{0, \hat{Q}} + p^{-1} \|\partial_x^2(u - P_p^y u)\|_{0, \hat{Q}}. $$

Thanks to the identities (61) and the one dimensional approximation properties (13), we can estimate the term $T_2$ from above as follows:

$$ T_2 \lesssim \|\partial_x u - P_p^y \partial_x u\|_{0, \hat{Q}} + p^{-1} \|\partial_x^2 u - P_p^y \partial_x^2 u\|_{0, \hat{Q}} \lesssim p^{-s-1} \left( \|\partial_x \partial_y^{s+1} u\|_{0, \hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0, \hat{Q}} \right). \qquad (64) $$

Collecting the estimates (63) and (64) in (62), we arrive at

$$ \|\partial_x(u - P_p u)\|_{0, \hat{Q}} \lesssim p^{-s-1} \left( \|\partial_x^{s+2} u\|_{0, \hat{Q}} + \|\partial_x \partial_y^{s+1} u\|_{0, \hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0, \hat{Q}} \right). $$

With similar arguments for the $y$ derivative term, we deduce (17).
**H**2 **estimates.** We begin by showing an upper bound on the second derivative with respect to x. Using the triangle inequality, the one dimensional approximation properties (13), the stability properties (14), and the identities (61), we obtain

$$
\begin{align*}
\|\partial_x^2(u - \mathcal{P}_p u)\|_{0,\hat{Q}} &\le \|\partial_x^2(u - \mathcal{P}_p^x u)\|_{0,\hat{Q}} + \|\partial_x^2 \mathcal{P}_p^x (u - \mathcal{P}_p^y u)\|_{0,\hat{Q}} \\
&\lesssim p^{-s} \|\partial_x^{s+2} u\|_{0,\hat{Q}} + \|\partial_x^2(u - \mathcal{P}_p^y u)\|_{0,\hat{Q}} = p^{-s} \|\partial_x^{s+2} u\|_{0,\hat{Q}} + \|\partial_x^2 u - \mathcal{P}_p^y \partial_x^2 u\|_{0,\hat{Q}} \\
&\lesssim p^{-s} \left( \|\partial_x^{s+2} u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0,\hat{Q}} \right).
\end{align*}
$$

Analogously, we can prove

$$
\| \partial_y^2 (u - \mathcal{P}_p u) \|_{0, \hat{Q}} \lesssim p^{-s} \left( \| \partial_y^{s+2} u \|_{0, \hat{Q}} + \| \partial_x^s \partial_y^2 u \|_{0, \hat{Q}} \right).
$$

Eventually, we cope with the mixed derivative term. Using the triangle inequality and the identi-
ties (61), we get

$$
\begin{align}
\|\partial_x \partial_y (u - \mathcal{P}_p u)\|_{0, \hat{Q}} &\le \|\partial_x \partial_y (u - \mathcal{P}_p^x u)\|_{0, \hat{Q}} + \|\partial_x \partial_y \mathcal{P}_p^x (u - \mathcal{P}_p^y u)\|_{0, \hat{Q}} \nonumber \\
&= \|\partial_x (\partial_y u - \mathcal{P}_p^x \partial_y u)\|_{0, \hat{Q}} + \|\partial_x \mathcal{P}_p^x \partial_y (u - \mathcal{P}_p^y u)\|_{0, \hat{Q}} =: S_1 + S_2. \tag{65}
\end{align}
$$

We estimate the two terms $S_1$ and $S_2$ separately. Using the one dimensional approximation prop-
erties (13), we can write

$$
S_1 \lesssim p^{-s} \| \partial_x^{s+1} \partial_y u \|_{0, \hat{Q}}. \quad (66)
$$

On the other hand, using the one dimensional approximation properties (13) and the stability
properties (14), we get

$$
\begin{align}
S_2 &\lesssim \|\partial_x \partial_y (u - \mathcal{P}_p^y u)\|_{0, \hat{Q}} + p^{-1} \|\partial_x^2 \partial_y (u - \mathcal{P}_p^y u)\|_{0, \hat{Q}} \nonumber \\
&= \|\partial_y (\partial_x u - \mathcal{P}_p^y \partial_x y)\|_{0, \hat{Q}} + p^{-1} \|\partial_y (\partial_x^2 u - \mathcal{P}_p^y \partial_x^2 u)\|_{0, \hat{Q}} \tag{67} \\
&\lesssim p^{-s} \left( \|\partial_x \partial_y^{s+1} u\|_{0, \hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0, \hat{Q}} \right). \nonumber
\end{align}
$$

Collecting the estimates (66) and (67) in (65), we arrive at

$$
\| \partial_x \partial_y (u - \mathcal{P}_p u) \|_{0, \hat{Q}} \lesssim p^{-s} \left( \| \partial_x^{s+1} \partial_y u \|_{0, \hat{Q}} + \| \partial_x \partial_y^{s+1} u \|_{0, \hat{Q}} + \| \partial_x^2 \partial_y^s u \|_{0, \hat{Q}} \right).
$$

Combining the estimates on all second derivatives terms, we obtain (18).

C Proof of Theorem 3.6 (3D case)

The continuity properties follow exactly as in the two dimensional case. Thus, we only show the details for the approximation properties. We split the proof in several steps. Recall that $\mathcal{P}_p^x, \mathcal{P}_p^y$, and $\mathcal{P}_p^z$ are the projections in Lemma 3.3 along the $x$, $y$, and $z$ directions, respectively.

**Some identities.** Analogous to their two dimensional counterparts in (61), we have the following identities:

$$
\begin{align}
\partial_y \mathcal{P}_p^x u &= \mathcal{P}_p^x (\partial_y u), & \partial_x \mathcal{P}_p^y u &= \mathcal{P}_p^y (\partial_x u), & \partial_y \mathcal{P}_p^z u &= \mathcal{P}_p^z (\partial_y u), \tag{68} \\
\partial_z \mathcal{P}_p^y u &= \mathcal{P}_p^y (\partial_z u), & \partial_x \mathcal{P}_p^z u &= \mathcal{P}_p^z (\partial_x u), & \partial_z \mathcal{P}_p^x u &= \mathcal{P}_p^x (\partial_z u). \notag
\end{align}
$$

**L**² **estimates.** The triangle inequality and the one dimensional approximation properties (13) imply

$$
\|u - \mathcal{P}_p u\|_{0, \hat{Q}} \le \|u - \mathcal{P}_p^x u\|_{0, \hat{Q}} + \|\mathcal{P}_p^x (u - \mathcal{P}_p^y \mathcal{P}_p^z u)\|_{0, \hat{Q}} \lesssim p^{-s-2} \|\partial_x^{s+2} u\|_{0, \hat{Q}} + \|\mathcal{P}_p^x (u - \mathcal{P}_p^y \mathcal{P}_p^z u)\|_{0, \hat{Q}}.
$$
We focus on the second term. To this aim, we use the stability properties (14), the triangle inequality, the one dimensional approximation properties (13), and the identities (68), and deduce

$$
\begin{align*}
\|\mathcal{P}_p^x(u - \mathcal{P}_p^y \mathcal{P}_p^z u)\|_{0, \hat{Q}} &\le \|u - \mathcal{P}_p^y \mathcal{P}_p^z u\|_{0, \hat{Q}} + p^{-2} \|\partial_x^2 u - \mathcal{P}_p^y \mathcal{P}_p^z \partial_x^2 u\|_{0, \hat{Q}} \\
&\lesssim \|u - \mathcal{P}_p^y u\|_{0, \hat{Q}} + \|\mathcal{P}_p^y (u - \mathcal{P}_p^z u)\|_{0, \hat{Q}} + p^{-2} \|\partial_x^2 u - \mathcal{P}_p^y \partial_x^2 u\|_{0, \hat{Q}} + p^{-2} \|\mathcal{P}_p^y (\partial_x^2 u - \mathcal{P}_p^z \partial_x^2 u)\|_{0, \hat{Q}} \\
&\lesssim p^{-s-2} \|\partial_y^{s+2} u\|_{0, \hat{Q}} + \|u - \mathcal{P}_p^z u\|_{0, \hat{Q}} + p^{-2} \|\partial_y^2 u - \mathcal{P}_p^z \partial_y^2 u\|_{0, \hat{Q}} + p^{-s-4} \|\partial_x^2 \partial_y^2 \partial_z^{s-2} u\|_{0, \hat{Q}} \\
&\qquad + p^{-2} \|\partial_x^2 u - \mathcal{P}_p^z \partial_x^2 u\|_{0, \hat{Q}} + p^{-4} \|\partial_x^2 \partial_y^2 u - \mathcal{P}_p^z \partial_x^2 \partial_y^2 u\|_{0, \hat{Q}}.
\end{align*}
$$

The $L^2$ estimates eventually follow from the one dimensional approximation properties (13).

H¹ estimates. We show the details for the x derivative, as the other two cases can be dealt with analogously. The triangle inequality and the one dimensional approximation properties (13) imply

$$
\begin{align*}
\|\partial_x(u - \mathcal{P}_p u)\|_{0, \hat{Q}} &\lesssim \|\partial_x(u - \mathcal{P}_p^x u)\|_{0, \hat{Q}} + \|\partial_x \mathcal{P}_p^x(u - \mathcal{P}_p^y \mathcal{P}_p^z u)\|_{0, \hat{Q}} \\
&\lesssim p^{-s-1} \|\partial_x^{s+2} u\|_{0, \hat{Q}} + \|\partial_x \mathcal{P}_p^x(u - \mathcal{P}_p^y \mathcal{P}_p^z u)\|_{0, \hat{Q}}.
\end{align*}
$$

We focus on the second term on the right-hand side. The stability properties (14) and the identities (68) entail

$$
\|\partial_x \mathcal{P}_p^x(u - \mathcal{P}_p^y \mathcal{P}_p^z u)\|_{0, \hat{Q}} \lesssim \|\partial_x u - \mathcal{P}_p^y \mathcal{P}_p^z \partial_x u\|_{0, \hat{Q}} + p^{-2} \|\partial_x^2 u - \mathcal{P}_p^y \mathcal{P}_p^z \partial_x^2 u\|_{0, \hat{Q}} =: T_1 + T_2.
$$

As for the term $T_1$, the triangle inequality, the identities (68), and the one dimensional approximation properties (13) imply

$$
T_1 \lesssim \|\partial_x u - \mathcal{P}_p^y \partial_x u\|_{0, \hat{Q}} + \|\mathcal{P}_p^y (\partial_x u - \mathcal{P}_p^z \partial_x u)\|_{0, \hat{Q}} \lesssim p^{-s-1} \|\partial_x \partial_y^{s+1} u\|_{0, \hat{Q}} + \|\mathcal{P}_p^y (\partial_x u - \mathcal{P}_p^z \partial_x u)\|_{0, \hat{Q}}.
$$

The second term on the right-hand side can be estimated using the stability properties (14), the identities (68), and the one dimensional approximation properties (13):

$$
\begin{align*}
\|\mathcal{P}_p^y (\partial_x u - \mathcal{P}_p^z \partial_x u)\|_{0, \hat{Q}} &\lesssim \|\partial_x u - \mathcal{P}_p^z \partial_x u\|_{0, \hat{Q}} + p^{-2} \|\partial_x \partial_y^2 u - \mathcal{P}_p^z \partial_x \partial_y^2 u\|_{0, \hat{Q}} \\
&\lesssim p^{-s-1} \left( \|\partial_x \partial_z^{s+1} u\|_{0, \hat{Q}} + \|\partial_x \partial_y^2 \partial_x^{s-1} u\|_{0, \hat{Q}} \right).
\end{align*}
$$

With similar arguments based on substituting $\partial_x u$ by $\partial_x^2 u$, we find an upper bound for $T_2$:

$$
T_2 \lesssim p^{-s-1} \left( \|\partial_x^2 \partial_y^s u\|_{0, \hat{Q}} + \|\partial_x^2 \partial_z^s u\|_{0, \hat{Q}} + \|\partial_x^2 \partial_y^2 \partial_z^{s-2} u\|_{0, \hat{Q}} \right).
$$

Collecting the above estimates, we arrive at

$$
\begin{align*}
\|\partial_x(u - \mathcal{P}_p u)\|_{0, \hat{Q}} &\lesssim p^{-s-1} \left( \|\partial_x^{s+2} u\|_{0, \hat{Q}} + \|\partial_x \partial_y^{s+1} u\|_{0, \hat{Q}} + \|\partial_x \partial_z^{s+1} u\|_{0, \hat{Q}} \right. \\
&\qquad \left. + \|\partial_x \partial_y^2 \partial_z^{s-1} u\|_{0, \hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0, \hat{Q}} + \|\partial_x^2 \partial_z^2 \partial_y^{s-2} u\|_{0, \hat{Q}} \right).
\end{align*}
$$

By similar arguments on the $y$ and $x$ partial derivatives, we deduce (20).

H² estimates. First, we show the details for the second $x$ derivative, since the second $y$ and $z$ derivatives can be dealt with analogously. The triangle inequality, the one dimensional approximation properties (13), and the identities (68) imply

$$
\begin{align*}
\|\partial_x^2 (u - \mathcal{P}_p u)\|_{0, \hat{Q}} &\le \|\partial_x^2 (u - \mathcal{P}_p^x u)\|_{0, \hat{Q}} + \|\partial_x^2 \mathcal{P}_p^x (u - \mathcal{P}_p^y \mathcal{P}_p^z u)\|_{0, \hat{Q}} \\
&\lesssim p^{-s} \|\partial_x^{s+2} u\|_{0, \hat{Q}} + \|\partial_x^2 u - \mathcal{P}_p^y \mathcal{P}_p^z \partial_x^2 u\|_{0, \hat{Q}}.
\end{align*}
$$

We are left with estimating the second term on the right-hand side. Applying a further triangle inequality, the one dimensional approximation properties (13), and the stability properties (14),
we arrive at

$$
\begin{aligned}
\|\partial_x^2 u - \mathcal{P}_p^y \mathcal{P}_p^z \partial_x^2 u\|_{0,\hat{Q}} &\le \|\partial_x^2 u - \mathcal{P}_p^y \partial_x^2 u\|_{0,\hat{Q}} + \|\mathcal{P}_p^y (\partial_x^2 u - \mathcal{P}_p^z \partial_x^2 u)\|_{0,\hat{Q}} \\
&\lesssim p^{-s} \|\partial_x^2 \partial_y^s u\|_{0,\hat{Q}} + \|\partial_x^2 u - \mathcal{P}_p^z \partial_x^2 u\|_{0,\hat{Q}} + p^{-2} \|\partial_x^2 \partial_y^2 u - \mathcal{P}_p^z \partial_x^2 \partial_y^2 u\|_{0,\hat{Q}} \\
&\lesssim p^{-s} \left( \|\partial_x^2 \partial_y^s u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_z^s u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^2 \partial_z^{s-2} u\|_{0,\hat{Q}} \right).
\end{aligned}
$$

Collecting the two above estimate above gives

$$
\|\partial_x^2(u - \mathcal{P}_p u)\|_{0,\hat{Q}} \lesssim p^{-s} \left( \|\partial_x^{s+2}u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^s u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_z^s u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^2 \partial_z^{s-2} u\|_{0,\hat{Q}} \right). \quad (69)
$$

Next, we focus on the mixed *xy* derivative and observe that the *yz* and *xz* counterparts are dealt with analogously. Using the triangle inequality, the one dimensional approximation properties (13), the identities (68), and the stability properties (14) leads us to

$$
\begin{aligned}
\|\partial_x \partial_y u - \mathcal{P}_p \partial_x \partial_y u\|_{0,\hat{Q}} &\le \|\partial_x \partial_y u - \mathcal{P}_p^x \partial_x \partial_y u\|_{0,\hat{Q}} + \|\mathcal{P}_p^x (\partial_x \partial_y u - \mathcal{P}_p^y \mathcal{P}_p^z \partial_x \partial_y u)\|_{0,\hat{Q}} \\
&\lesssim p^{-s} \|\partial_x^{s+1} \partial_y u\|_{0,\hat{Q}} + \|\partial_y \partial_x u - \partial_y \mathcal{P}_p^y \mathcal{P}_p^z \partial_x u\|_{0,\hat{Q}} + p^{-1} \|\partial_y \partial_x^2 u - \partial_y \mathcal{P}_p^y \mathcal{P}_p^z \partial_x^2 u\|_{0,\hat{Q}} \\
&= p^{-s} \|\partial_x^{s+1} \partial_y u\|_{0,\hat{Q}} + T_3 + T_4.
\end{aligned}
$$

We estimate the terms *A* and *B* separately. First, we focus on *A*. Using the triangle inequality, the one dimensional approximation properties (13), the identities (68), and the stability properties (14), entail

$$
\begin{aligned}
T_3 &\le \|\partial_y \partial_x u - \partial_y \mathcal{P}_p^y \partial_x u\|_{0,\hat{Q}} + \|\partial_y \mathcal{P}_p^y \partial_x u - \partial_y \mathcal{P}_p^y \partial_x u\|_{0,\hat{Q}} \\
&\lesssim p^{-s} \|\partial_x \partial_y^{s+1} u\|_{0,\hat{Q}} + \|\partial_y \partial_x u - \mathcal{P}_p^z \partial_y \partial_x u\|_{0,\hat{Q}} + p^{-1} \|\partial_y^2 \partial_x u - \mathcal{P}_p^z \partial_y^2 \partial_x u\|_{0,\hat{Q}} \\
&\lesssim p^{-s} \left( \|\partial_x \partial_y^{s+1} u\|_{0,\hat{Q}} + \|\partial_x \partial_y \partial_z^s u\|_{0,\hat{Q}} + \|\partial_x \partial_y^2 \partial_z^{s-1} u\|_{0,\hat{Q}} \right).
\end{aligned}
$$

Next, we focus on the term *B*. Using the triangle inequality, the one dimensional approximation properties (13), the identities (68), and the stability properties (14), leads to

$$
\begin{aligned}
T_4 &\le p^{-1} \|\partial_y \partial_x^2 u - \partial_y \mathcal{P}_p^y \partial_x^2 u\|_{0,\hat{Q}} + p^{-1} \|\partial_y \mathcal{P}_p^y \partial_x^2 u - \partial_y \mathcal{P}_p^y \mathcal{P}_p^z \partial_x^2 u\|_{0,\hat{Q}} \\
&\lesssim p^{-s} \|\partial_x^2 \partial_y^s\|_{0,\hat{Q}} + p^{-1} \|\partial_y \partial_x^2 u - \mathcal{P}_p^z \partial_y \partial_x^2 u\|_{0,\hat{Q}} + p^{-2} \|\partial_x^2 \partial_y^2 u - \mathcal{P}_p^z \partial_x^2 \partial_y^2 u\|_{0,\hat{Q}} \\
&\lesssim p^{-s} \left( \|\partial_x^2 \partial_y^s\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y \partial_z^s\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^2 \partial_z^{s-2}\|_{0,\hat{Q}} \right).
\end{aligned}
$$

Collecting the above estimates gives

$$
\begin{aligned}
\|\partial_x \partial_y u - \mathcal{P}_p \partial_x \partial_y u\|_{0,\hat{Q}} &\lesssim p^{-s} \left( \|\partial_x^{s+1} \partial_y u\|_{0,\hat{Q}} + \|\partial_x \partial_y^{s+1} u\|_{0,\hat{Q}} + \|\partial_x \partial_y \partial_z^s u\|_{0,\hat{Q}} \|\partial_x \partial_y^2 \partial_z^{s-1} u\|_{0,\hat{Q}} \right. \\
&\qquad \left. + \|\partial_x^2 \partial_y^s u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y \partial_z^{s-1} u\|_{0,\hat{Q}} + \|\partial_x^2 \partial_y^2 \partial_z^{s-2} u\|_{0,\hat{Q}} \right).
\end{aligned}
\tag{70}
$$

Bound (21) follows combining the estimates on the second derivative term (69), the mixed derivative term (70), and the corresponding estimates for similar derivative terms.