Z Dong

L Mascotto
email: lorenzo.mascotto@unimib.it

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

Keywords: AMS subject classification: 65N12, 65N30, 65N50 discontinuous Galerkin, optimal convergence, a priori error estimation, p-version

est destinée au dépôt

Introduction

The numerical approximation of solutions to the biharmonic problem is a challenging task tracing back to the 60ties [START_REF] Argyris | The TUBA family of plate elements for the matrix displacement method[END_REF]; see also the later work [START_REF] Douglas | A family of finite elements with optimal approximation properties for various Galerkin methods for 2nd and 4th order problems[END_REF]. 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., [START_REF] Morley | The triangular equilibrium element in the solution of plate bending problems[END_REF][START_REF] Wang | The Morley element for fourth order elliptic equations in any dimensions[END_REF].

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., [START_REF] Beirão Da Veiga | A family of C 0 finite elements for Kirchhoff plates I: Error analysis[END_REF][START_REF] Brenner | An a posteriori error estimator for a quadratic C 0 -interior penalty method for the biharmonic problem[END_REF][START_REF] Brenner | C 0 interior penalty methods for fourth order elliptic boundary value problems on polygonal domains[END_REF][START_REF] Engel | 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[END_REF]. The continuity constraint can we 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 [START_REF] Baker | Finite element methods for elliptic equations using nonconforming elements[END_REF] in the late 70ties; amongst others, we recall the following contributions, including the p-and hp-versions of the scheme: [START_REF] Dong | Discontinuous Galerkin methods for the biharmonic problem on polygonal and polyhedral meshes[END_REF][START_REF] Dong | Residual-based a posteriori error estimates for hp-discontinuous Galerkin discretizations of the biharmonic problem[END_REF][START_REF] Georgoulis | Discontinuous Galerkin methods for the biharmonic problem[END_REF][START_REF] Gudi | Mixed discontinuous Galerkin finite element method for the biharmonic equation[END_REF][START_REF] Mozolevski | A priori error analysis for the hp-version of the discontinuous Galerkin finite element method for the biharmonic equation[END_REF][START_REF] Mozolevski | hp-version a priori error analysis of interior penalty discontinuous Galerkin finite element approximations to the biharmonic equation[END_REF][START_REF] Süli | hp-version interior penalty DGFEMs for the biharmonic equation[END_REF].

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 [START_REF] Perugia | An hp-analysis of the local discontinuous Galerkin method for diffusion problems[END_REF] using the error analysis in the L 2 norm for the local discontinuous Galerkin method and later analyzed in [START_REF] Stamm | hp-optimal discontinuous Galerkin methods for linear elliptic problems[END_REF] 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 [START_REF] Brenner | C 0 interior penalty methods for fourth order elliptic boundary value problems on polygonal domains[END_REF], 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 [START_REF] Suri | The p-version of the finite element method for elliptic equations of order 2l[END_REF] 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 [START_REF] Georgoulis | On the suboptimality of the p-version interior penalty discontinuous Galerkin method[END_REF]. 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 R d , d = 1, 2, 3. We denote the Sobolev space of order s ∈ 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 (•, •) s,D , which induces the norm and seminorm • s,D and |•| s,D . Further, given ∈ R + and k ∈ N, we introduce the anisotropic Sobolev space:

H ,k (D) := v ∈ H (D) ∂ α1 x1 . . . ∂ α d x d u 0,D < ∞ ∀α = (α 1 , . . . , α d ) ∈ N d s. th.
|α| = k, α j ≤ 2 ∀j ∈ {1, . . . , d} \ {J} for one J ∈ {1, . . . , d} .

(

The space H ,k (D) contains H 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: ∇ as the gradient; ∆ = ∇•∇ as the Laplacian; ∆ 2 as the bilaplacian; D 2 = ∇∇ T as the Hessian.

Given p ∈ N, let P p (D) and 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 u and v in R d by u • v. Further, we denote the entrywise product between two tensors u and v in R d×d by u : v.

Given two positive quantities a and b, we use the short-hand notation "a b" if there exists a positive constant c independent of the discretization parameters such that "a ≤ c b". Moreover, we write "a ≈ b" if and only if "a b" and "b a".

The continuous problem.

Given Ω ⊂ R d , d = 2, 3, a bounded polygonal/polyhedral domain with unit outward normal vector n Ω and f ∈ L 2 (Ω), we consider the following biharmonic problem:

find u : Ω → R such that ∆ 2 u = f in Ω, u = 0, n Ω • ∇u = 0 on ∂Ω.
(2)

Introduce V := H 2 0 (Ω), B(u, v) = Ω D 2 u : D 2 v ∀u, v ∈ V.
A weak formulation of problem (2) reads

find u ∈ V such that B(u, v) = (f, v) 0,Ω ∀v ∈ V. (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., [START_REF] Beirão Da Veiga | A posteriori error analysis for the Morley plate element with general boundary conditions[END_REF].

2 Interior penalty discontinuous Galerkin methods for the biharmonic problem

In this section, we introduce the IPDG for problem [START_REF] Argyris | The TUBA family of plate elements for the matrix displacement method[END_REF]. 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.

Meshes

We consider sequences {T n } n of decompositions of Ω 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 T n by F n . We further distinguish interior and boundary faces introducing F I n and F B n , respectively. The set of vertices of T n is V n . With each element K ∈ T n , we associate its diameter h K and collect all the diameters in h : Ω → R + defined as

h := h K if x ∈ K ∈ T n , h F if x ∈ F ∈ F n .
We denote the unit normal vector to any interior face F with respect to the global orientation of F by n. We always fix a outward pointing vector for boundary faces. Throughout, we assume quasi-uniformity of the mesh T n , i.e., there exists c mesh > 0 such that

c -1 mesh h K1 ≤ h K2 ≤ c mesh h K1 ∀K 1 , K 2 ∈ T n .
In what follows, we consider a uniform distribution of the polynomial degree p ∈ N over T n .

Variable polynomial degrees can be used [START_REF] Dong | Residual-based a posteriori error estimates for hp-discontinuous Galerkin discretizations of the biharmonic problem[END_REF][START_REF] Ch | p-and hp-Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics[END_REF] but are not the main focus of this work. Define the scalar broken Sobolev space

H 1 (T n , Ω) := {v ∈ L 2 (Ω) | v |K ∈ H 1 (K)}, H ,k (T n ) := {v ∈ L 2 (Ω) | v |K ∈ H ,k (K)} ∀ ∈ R + , k ∈ N.
We introduce the jump • : L 2 (F ) → R and the average {•} : L 2 (F ) → R operators over a face F ∈ F n as follows: if the face F ∈ F I n is an interior face shared by the elements K + and K -, given v ∈ H 1 (T n , Ω) with restrictions over the elements K + and K -denoted by v + and v -, respectively,

v F (x) = v (x) := v + (x) -v -(x), {v} F (x) = {v}(x) := 1 2 (v + (x) + v -(x)) ∀x ∈ F.
On boundary faces F ∈ F B n , 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 T n by adding a subscript n; for instance, ∇ n is the piecewise gradient over T n .

The method

For each mesh T n and polynomial degree p, we introduce the space

V n := {v n ∈ L 2 (Ω) | v n|K ∈ P p (K) ∀K ∈ T n }.
We endow V n with the broken Sobolev norm and seminorm of order s > 0

• 2 s,Tn := K∈Tn • 2 s,K , |•| 2 s,Tn := K∈Tn |•| 2 s,K ,
and introduce the lifting operator L :

V n ⊕ V → [V n ] d×d , d = 2, 3, as Ω L(u n ) : v n := Fn {n • (∇ • v n )} u n - Fn {(v n )n} • ∇u n ∀v n ∈ [V n ] d×d . (4)
Let σ : F n → R + and τ : F n → R + two piecewise constant stabilization parameters, which we shall fix in what follows; see Section 2.3 below. Given D 2 n the piecewise Hessian over T n , we define the bilinear form B n :

V n × V n → R as B n (u n , v n ) := Ω D 2 n u n : D 2 n v n + Ω L(u n ) : D 2 n v n + L(v n ) : D 2 n u n + Fn (σ u n v n + τ ∇u n • ∇v n ) ∀u n , v n ∈ V n . (5)
We consider the following IPDG for the approximation of solutions to (3); see [START_REF] Dong | Residual-based a posteriori error estimates for hp-discontinuous Galerkin discretizations of the biharmonic problem[END_REF]:

find u n ∈ V n such that B n (u n , v n ) = (f, v n ) ∀v n ∈ V n . (6)
Other dG formulations might be used [START_REF] Brenner | C 0 interior penalty methods for fourth order elliptic boundary value problems on polygonal domains[END_REF][START_REF] Georgoulis | Discontinuous Galerkin methods for the biharmonic problem[END_REF]; yet, we stick to (6) as it is the most similar to that of [START_REF] Stamm | hp-optimal discontinuous Galerkin methods for linear elliptic problems[END_REF]. Eventually, we define the natural norm associated with the bilinear form

B n (•, •): for all v ∈ V n ⊕ V , v 2 dG := D 2 n v 2 0,Ω + F ∈Fn τ 1 2 ∇v 2 0,F + F ∈Fn σ 1 2 v 2 0,F =: D 2 n v 2 0,Ω + τ 1 2 ∇v 2 0,Fn + σ 1 2 v 2 0,Fn , (7)
where

(u, v) 0,Fn := F ∈Fn (u, v) 0,F , u 2 0,Fn := (u, u) 0,Fn ∀u, v ∈ V n ⊕ V.

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 L(•) in (4).

Lemma 2.1. There exists a positive constant c S independent of h and p such that

L(v n ) 0,Ω ≤ c S   p 6 h 3 1 2 v n 0,Fn + p 2 h 1 2 ∇v n 0,Fn   ∀v n ∈ V n .
Consequences of Lemma 2.1 are the coercivity and continuity of the dG bilinear form B n (•, •) under suitable choices of the dG parameters σ and τ , namely, given c σ , c τ > 2c S + 1 2 ,

σ := c σ p 6 h 3 , τ := c τ p 2 h . (8)
See, e.g., [START_REF] Georgoulis | Discontinuous Galerkin methods for the biharmonic problem[END_REF]Lemma 5.2].

Lemma 2.2. The bilinear form B n (•, •) is coercive and continuous:

B n (v, v) ≥ 1 2 v 2 dG , B n (u, v) ≤ 2 u dG v dG ∀u, v ∈ V n ⊕ V. (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 h and p such that

(f, v) 0,Ω ≤ c C f 0,Ω v dG ∀v ∈ V n ⊕ V.
Proof. The assertion follows from the Cauchy-Schwarz inequality and Poincaré-Friedrichs type inequalities for piecewise H 2 functions; see, e.g., [START_REF] Brenner | Poincaré-Friedrichs inequalities for piecewise H 2 functions[END_REF]Corollary 4.2].

The well posedness of method (6) follows combining the coercivity of B n (•, •) and the continuity of (f, •) 0,Ω with respect to the • dG norm.

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

u n dG ≤ 2c C f 0,Ω .

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., [START_REF] Canuto | Approximation results for orthogonal polynomials in Sobolev spaces[END_REF][START_REF] Ch | p-and hp-Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics[END_REF] for details. Given I := (-1, 1), let Π 0 p : L 2 ( I) → P p ( I) and Π 1 p : H 1 ( I) → P p ( I) be the L 2 and H 1 projectors onto P p ( I), respectively: for all q p ∈ P p ( I),

(q p , v -Π 0 p v) 0, I = 0 ∀v ∈ L 2 ( I), (q p , (v -Π 1 p v) ) 0, I = 0 (v -Π 1 p v)(±1) = 0 ∀v ∈ H 1 ( I).
The following approximation estimates are valid; see, e.g., [START_REF] Ch | p-and hp-Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics[END_REF]Chapter 3] and [START_REF] Georgoulis | Discontinuous Galerkin methods on shape-regular and anisotropic meshes[END_REF]Lemma 3.5].

Lemma 3.1. For all p ∈ N, given u ∈ H k ( I) with k ≥ 1 and s := min(k, p + 1), there exists a positive constant c independent of p and k such that

u -Π 0 p u 0, I ≤ cp -s |u| s, I , (u -Π 1 p u) 0, I ≤ cp -s+1 |u| s, I , (u -Π 0 p u) 0, I ≤ cp -s+ 3 2 |u| s, I . (10)
Next, we recall best approximation estimates of the elementwise L 2 projector Π 0 p : L 2 (K) → Q p (K) with respect to the boundary L 2 norm on a d-dimensional tensor product-element K; see, e.g., [START_REF] Houston | Discontinuous hp-finite element methods for advection-diffusion-reaction problems[END_REF]Lemma 3.9]. Lemma 3.2. Let K ⊂ 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 ∈ H k (K), k ≥ 1. Then, for all 1 ≤ s ≤ min(p + 1, k), p ∈ N, there exists a positive constant c > 0 independent of h K , k, and p such that

u -Π 0 p u 0,∂K ≤ c h K p + 1 s-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 [START_REF] Beirão Da Veiga | Some estimates for h-p-k-refinement in isogeometric analysis[END_REF]Corollary 2].

Lemma 3.3. Let I := (-1, 1) and k ≥ 0. Then, for all p ≥ 3, there exists a projection operator P 1 p : H k+2 ( I) → P p ( I) such that the following continuity properties are valid: for all u ∈ H k+2 ( I)

(P 1 p u)(±1) = u(±1), (P 1 p u) (±1) = u (±1). ( 11
)
Moreover, the following upper bounds are valid: given s := min(k, p -1),

(u -P 1 p u) 2 0, I ≤ (p -s -1)! (p + s -1)! u (s+2) 2 0, I , (u -P 1 p u) 2 0, I ≤ (p -s -1)! (p + s -1)! 1 (p -1)p u (s+2) 2 0, I , u -P 1 p u 2 0, I ≤ (p -s -1)! (p + s -1)! c(p) u (s+2) 2 0, I , (12)
where c(p) is independent of k and

c(p) ≈ p -4 .
Proof. We provide the details of the proof in Appendix A for completeness.

An application of Stirling's formula, see, e.g., [START_REF] Ch | p-and hp-Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics[END_REF]Corollary 3.12], yields the following p-optimal approximation result for the projector P 1 p .

Corollary 3.4. Let I := (-1, 1) and p ≥ 3. Then, there exists a positive constant c 1D > 0 independent of p and k such that the projection operator P 1 p : H k+2 ( I) → P p ( I), k ≥ 0, from Lemma 3.3 satisfies, for s := min(k, p -1),

c 1D (p -1) 2 u -P 1 p u 0, I + (p -1) (u -P 1 p u) 0, I + (u -P 1 p u) 0, I ≤ e 2 s p-1 (p -1) -s u (s+2) 0, I ≤ e 2 (p -1) -s u (s+2) 0, I . (13)
Based on this, we derive stability properties for the projection operator constructed in Lemma 3.3.

Corollary 3.5. The projection operator P 1 p • 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 ∈ H 2 ( I), (P 1 p u) 0, I ≤ u 0, I , (P 1 p u) 0, I ≤ u 0, I + p -1 u 0, I , P 1 p u 0, I ≤ u 0, I + cp -2 u 0, I .

Proof. The first inequality follows from the definition of the operator P 1 p 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 [START_REF] Brenner | Poincaré-Friedrichs inequalities for piecewise H 2 functions[END_REF] with s = 0.

H 2 estimates on tensor product elements. Given P x p , P y p , and P z p 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

P p • := P x p P y p •, P p • := P x p P y p P z p • .
By tensor product arguments, we deduce the two and three dimensional counterparts of Lemma 3.3.

Theorem 3.6. Let 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 n Q . Further, let p ≥ 3 and u ∈ H 2,k+2 ( Q), k ≥ 0.

Then, the following continuity properties on ∂ Q are valid:

P p u |ei = P x p (u |ei ) if i is even, P y p (u |ei ) if i is odd, (n Q • ∇P p u) |ei = n Q • ∇P y p (u |ei ) if i is even, n Q • ∇P x p (u |ei ) if i is odd. (15)
In particular, P p u and u share the same value at the vertices of Q. If 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 Q = (-1, 1) 2 , there exists a positive constant c independent of p and k such that 3 , there exists a positive constant c independent of p and k such that

u -P p u 0, Q ≤ cp -s-2 ∂ s+2 x u 0, I + ∂ s+2 y u 0, I + ∂ 2 x ∂ s y u 0, I ; (16)
∇(u -P p u) 0, Q ≤ cp -s-1 ∂ s+2 x u 0, Q + ∂ x ∂ s+1 y u 0, Q + ∂ 2 x ∂ s y u 0, Q (17)
+ ∂ s x ∂ 2 y u 0, Q + ∂ s+1 x ∂ y u 0, Q + ∂ s+2 y u 0, Q ; D 2 (u -P p u) 0, Q ≤ cp -s ∂ s+2 x u 0, Q + ∂ 2 x ∂ s y u 0, Q + ∂ s+2 y u 0, Q (18)
+ ∂ s x ∂ 2 y u 0, Q ∂ s+1 x ∂ y u 0, Q + ∂ x ∂ s+1 y u 0, Q + ∂ 2 x ∂ s y u 0, Q . Instead, if Q = (-1, 1)
u -P p u 0, Q ≤ cp -s-2 ∂ s+2 x u 0, I + ∂ s+2 y u 0, I + ∂ s+2 z u 0, I (19)
+ ∂ 2 x ∂ s z u 0, I + ∂ 2 y ∂ s z u 0, I + ∂ 2 x ∂ 2 y ∂ s-2 z u 0, I ; ∇(u -P p u) 0, Q ≤ cp -s-1 ∂ s+2 x u 0, Q + ∂ s+2 y u 0, Q + ∂ s+2 z u 0, Q + ∂ x ∂ s+1 y u 0, Q (20)
+ ∂ y ∂ s+1 z u 0, Q + ∂ z ∂ s+1 x u 0, Q + ∂ x ∂ s+1 y u 0, Q + ∂ y ∂ s+1 z u 0, Q + ∂ z ∂ s+1 x u 0, Q + ∂ x ∂ 2 y ∂ s-1 z u 0, Q + ∂ y ∂ 2 z ∂ s-1 x u 0, Q + ∂ z ∂ 2 x ∂ s-1 y u 0, Q + ∂ 2 x ∂ s y u 0, Q + ∂ 2 y ∂ s z u 0, Q + ∂ 2 z ∂ s x u 0, Q + ∂ 2 x ∂ s z u 0, Q + ∂ 2 y ∂ s x u 0, Q + ∂ 2 z ∂ s y u 0, Q + ∂ 2 x ∂ 2 y ∂ s-2 z u 0, Q + ∂ 2 y ∂ 2 z ∂ s-2 x u 0, Q + ∂ 2 z ∂ 2 x ∂ s-2 y u 0, Q ; D 2 (u -P p u) 0, Q ≤ cp -s ∂ s+2 x u 0, Q + ∂ s+2 y u 0, Q + ∂ s+2 z u 0, Q + ∂ 2 x ∂ s y u 0, Q (21)
+ ∂ 2 y ∂ s z u 0, Q + ∂ 2 z ∂ s x u 0, Q + ∂ 2 x ∂ s z u 0, Q + ∂ 2 y ∂ s x u 0, Q + ∂ 2 z ∂ s y u 0, Q + ∂ s+1 x ∂ y u 0, Q + ∂ s+1 y ∂ z u 0, Q + ∂ s+1 z ∂ x u 0, Q + ∂ x ∂ s+1 y u 0, Q + ∂ y ∂ s+1 z u 0, Q + ∂ z ∂ s+1 x u 0, Q + ∂ x ∂ y ∂ s z u 0, Q + ∂ y ∂ z ∂ s x u 0, Q + ∂ z ∂ x ∂ s y u 0, Q + ∂ x ∂ 2 y ∂ s-1 z u 0, Q + ∂ y ∂ 2 z ∂ s-1 x u 0, Q + ∂ z ∂ 2 x ∂ s-1 y u 0, Q + ∂ 2 x ∂ y ∂ s-1 z u 0, Q + ∂ 2 y ∂ z ∂ s-1 x u 0, Q + ∂ 2 z ∂ x ∂ s-1 y u 0, Q + ∂ 2 x ∂ 2 y ∂ s-2 z u 0, Q + ∂ 2 y ∂ 2 z ∂ s-2 x u 0, Q + ∂ 2 z ∂ 2 x ∂ s-2 y u 0, Q .
Proof. As for the two dimensional case, we refer to [START_REF] Beirão Da Veiga | Some estimates for h-p-k-refinement in isogeometric analysis[END_REF] and present some details in Appendix B for completeness. We report the proof for the three dimensional estimates in Appendix C.

Remark 2. In Theorem 3.6, the p-explicit approximation bounds are derived for functions in the anisotropic Sobolev space H 2,k+2 ( Q) defined in [START_REF] Ainsworth | H 2 -stable polynomial liftings on triangles[END_REF]. 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 P 1 p ; 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 n K , p ≥ 3, and u ∈ H k+2 (K), k ≥ 0. Then, the following continuity properties on ∂K are valid:

P p u |Fi = P x p (u |Fi ) if i is even, P y p (u |Fi ) if i is odd, (n K • ∇P p u) |Fi = n K • ∇P y p (u |Fi ) if i is even, n K • ∇P x p (u |Fi ) if i is odd.
In particular, 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 independent of h, k, and p such that

p 2 h 2 K u -P p u 0,K + p h K ∇(u -P p u) 0,K + D 2 (u -P p u) 0,K ≤ c h s p s |u| s+2,K . (22)
Global H 2 estimates. Corollary 3.7 implies hp-optimal error estimates by means of global H 2 piecewise polynomials over the given tensor product-type meshes. More precisely, we have the following result.

P p : H 2 (Ω) ∩ H 2,2d (T n ) → H 2 (Ω) ∩ V n and a positive constant c independent of h, k, and p such that, for all u ∈ H 2 (Ω) ∩ H k+2 (T n ), p 2 h -2 (u -P p u) 0,Ω + p h -1 ∇(u -P p u) 0,Ω + D 2 (u -P p u) 0,Ω ≤ cp -s h s D s+2 u 0,Tn .
The operator P p is defined on the space H 2 (Ω)∩H 2,2d (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.

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.

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 2 polynomial projector from Theorem 3.8.

Henceforth, given D ⊂ R d , d = 2, 3, either a parallelogram (in two dimensions) or a parallelepiped (in three dimensions), we introduce Π 0 p as the composition of the one dimensional L 2 projectors in Lemma 3.1 along the coordinate directions. In other words, we let

Π 0 p = Π 0,x p Π 0,y p if d = 2, Π 0 p = Π 0,x p Π 0,y p Π 0,z p if d = 3.
Theorem 4.1. Let 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 ∈ H s+2 (Ω), s := min(k, p -1), p ≥ 2. Then, there exists a positive constant c independent of h, k, and p such that

u -u n dG ≤ c h p s u s+2,Ω .
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 dG ≤ 5 inf vn∈Vn u -v n dG + 2 sup ξ∈Vn |B n (u -u n , ξ)| ξ dG . ( 23
)
The first term represents a best polynomial approximation result; the second measures the inconsistency of the scheme.

We have standard hp-optimal (piecewise) polynomial estimates on the first term:

inf vn∈Vn u -v n dG h p s u s+2,Ω . (24)
To show [START_REF] Houston | Discontinuous hp-finite element methods for advection-diffusion-reaction problems[END_REF], 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 [START_REF] Gudi | Mixed discontinuous Galerkin finite element method for the biharmonic equation[END_REF]. Proceeding as in [START_REF] Georgoulis | Discontinuous Galerkin methods for the biharmonic problem[END_REF]Lemma 5.4], we set η := u -u n and write

B n (u -u n , ξ) = Fn n • (∇ • D 2 n η) -n • (∇ • Π 0 p D 2 n η) ξ - Fn (D 2 n η)n -(Π 0 p D 2 n η)n ∇ξ =: A 1 + A 2 . ( 25
)
Since Π 0 p • preserves polynomials of degree p, we have

D 2 n η -Π 0 p D 2 n η = D 2 n (u -u n ) -Π 0 p (D 2 n u) + Π 0 p (D 2 n u n ) = D 2 n u -Π 0 p D 2 n u. (26)
Using ( 7), [START_REF] Beirão Da Veiga | A posteriori error analysis for the Morley plate element with general boundary conditions[END_REF], and ( 26), we write

A 1 h 3 2 p 3 ∇ • (D 2 n η -Π 0 p D 2 n η) 0,Fn σ 1 2 ξ 0,Fn h 3 2 p 3 ∇ • (D 2 n u -Π 0 p D 2 n u) 0,Fn ξ dG . (27)
Next, we focus on the term A 2 and use similar manipulations:

A 2 h 1 2 p D 2 n η -Π 0 p D 2 n η 0,Fn τ 1 2 ξ 0,Fn h 1 2 p D 2 n u -Π 0 p D 2 n u 0,Fn ξ dG . (28)
We insert ( 27) and ( 28) in [START_REF] Kanschat | Divergence-conforming discontinuous Galerkin methods and C 0 interior penalty methods[END_REF] and get

|B n (u -u n , ξ)| h 3 2 p 3 ∇ • (D 2 n u -Π 0 p D 2 n u) 0,Fn + h 1 2 p D 2 n u -Π 0 p D 2 n u 0,Fn ξ dG .
Thus, we arrive at the following bound on the inconsistency term in [START_REF] Gudi | Mixed discontinuous Galerkin finite element method for the biharmonic equation[END_REF]:

sup ξ∈Vn |B n (u -u n , ξ)| ξ dG h 3 2 p 3 ∇ • (D 2 n u -Π 0 p D 2 n u) 0,Fn + h 1 2 p D 2 n u -Π 0 p D 2 n u 0,Fn =: T 1 +T 2 . ( 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:

∇ • (D 2 n u -Π 0 p D 2 n u) = ∂ x (∂ 2 x u -Π 0 p ∂ 2 x u) + ∂ y (∂ xy u -Π 0 p ∂ xy u) ∂ x (∂ yx u -Π 0 p ∂ yx u) + ∂ y (∂ 2 y u -Π 0 p ∂ 2 y u)
.

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

∂ x (∂ 2 x u -Π 0 p ∂ 2 x u) 0,F ≤ ∂ 3 x u -Π 0 p-1 ∂ 3 x u 0,F + Π 0 p-1 ∂ 3 x u -∂ x Π 0 p ∂ 2 x u 0,F =: Z 1 + Z 2 . ( 30
)
As for the term Z 1 , Lemma 3.2 yields

Z 1 h p s-3 2 ∂ 3 x u s-1,K ≤ h p s-3 2 u s+2,K ∀s ≥ 0. ( 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:

Z 2 p h 1 2 Π 0 p-1 ∂ 3 x u -∂ x Π 0 p ∂ 2 x u 0,K ≤ p h 1 2 ∂ 3 x u -Π 0 p-1 ∂ 3 x u 0,K + ∂ x (∂ 2 x u -Π 0 p ∂ 2 x u) 0,K = Z 2,1 + Z 2,2 . (32)
Using the approximation properties of the L 2 projector yields

Z 2,1 p h 1 2 h p s-1 ∂ 3 x u s-1,K ≤ p h 1 2 h p s-1 u s+2,K .
On the other hand, standard manipulations entail

Z 2,2 p h 1 2 ∂ x (∂ 2 x u -Π 0,y p ∂ 2 x u) 0,K + ∂ x (Π 0,y p ∂ 2 x u -Π 0,y p Π 0,x p ∂ 2 x u) 0,K ≤ p h 1 2 ∂ x (∂ 2 x u -Π 0,y p ∂ 2 x u) 0,K + ∂ x (∂ 2 x u -Π 0,x p ∂ 2 x u) 0,K .
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 p 3 2 h 1 2 h p s-1 ∂ 3 x u s-1,K p 3 2 h 1 2 h p 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 p 3 2 h 1 2 h p s-1 u s+2,K , s ≥ 0. (33)
Inserting next [START_REF] Perugia | An hp-analysis of the local discontinuous Galerkin method for diffusion problems[END_REF] and ( 33) in [START_REF] Mozolevski | hp-version a priori error analysis of interior penalty discontinuous Galerkin finite element approximations to the biharmonic equation[END_REF] gives

∂ x (∂ 2 x u -Π 0 p ∂ 2 x u) 0,F p 3 2 h 1 2 h p s-1 u s+2,K = p 5 2 h 3 2 h p s u s+2,K , s ≥ 0.
In overall, also recalling the factor h

3 2
p 3 in front of T 1 , we deduce

T 1 p -1 2 h p s u s+2,Ω . (34)
As for the term T 2 on the right-hand side of (29), we use Lemma 3.2 and obtain

T 2 p -1 2 h p s u s+2,Ω . (35)
Inserting ( 34) and ( 35) in ( 29) leads to

sup ξ∈Vn |B n (u -u n , ξ)| ξ dG p -1 2 h p s u s+2,Ω ≤ h p s u s+2,Ω . (36)
A combination of ( 25), [START_REF] Houston | Discontinuous hp-finite element methods for advection-diffusion-reaction problems[END_REF], and (36) yields the assertion.

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 [START_REF] Georgoulis | On the suboptimality of the p-version interior penalty discontinuous Galerkin method[END_REF], 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.

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 Ω 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 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 Ω is convex, there exists a positive constant c independent of h and p such that

u -u n 0,Ω ≤ c h p 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 Φ such that

∆ 2 Φ = u -u n in Ω, u = 0, n Ω • ∇u = 0 on ∂Ω.
In weak formulation, this problem reads

find u ∈ V such that B(u, v) = (u -u n , v) 0,Ω ∀v ∈ V. (37)
Due to the convexity of Ω, the solution Φ to [START_REF] Wang | A robust numerical method for Stokes equations based on divergence-free H(div) finite element methods[END_REF] belongs to H 4 (Ω) and the following regularity estimate holds true:

Φ 4,Ω u -u n 0,Ω ; (38)
see, e.g., [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF] and [26, Theorem 6, page 182] for the two and three dimensional cases, respectively. Using [START_REF] Wang | A robust numerical method for Stokes equations based on divergence-free H(div) finite element methods[END_REF], splitting the integrals elementwise, and integrating by parts imply

u -u n 2 0,Ω = (u -u n , ∆ 2 Φ) 0,Ω = K∈Tn (-(∇(u -u n ), ∇∆Φ) 0,K + (u -u n , ∇∆Φ • n K ) 0,∂K ) .
Integrating by parts, and using (4) and the smoothness of Φ give

u -u n 2 0,Ω = K∈Tn (D 2 (u -u n ), D 2 Φ) 0,K + ( u -u n , {∇∆Φ • n}) 0,Fn + ( Φ , {∇∆(u -u n ) • n}) 0,Fn -( ∇(u -u n ) , {D 2 Φ • n}) 0,Fn -( ∇Φ , {D 2 (u -u n ) • n}) 0,Fn + (σ 1 2 u n , v n ) 0,Fn + (τ 1 2 ∇u n , ∇v n ) 0,Fn .
As, e.g., in [21, Theorem 4.5], we write the right-hand side as a combination of B n (u -u n , Φ) and a reminder R. Standard computations reveal that

R = ( u -u n , {∇∆Φ • n}) 0,Fn -( ∇(u -u n ) , {D 2 Φ • n}) 0,Fn -(L(u -u n ), Π 0 p D 2 n Φ) 0,Ω + ( Φ , {∇∆(u -u n ) • n}) 0,Fn -( ∇Φ , {D 2 (u -u n ) • n}) 0,Fn -(LΦ, Π 0 p (D 2 n (u -u n ))) 0,Ω = A + B.
Due to the regularity of Φ, we have B = 0. We deduce

u -u n 2 0,Ω = B n (u -u n , Φ) + R, (39)
where

R = ( u -u n , {n • ∇(D 2 Φ -Π 0 p D 2 Φ)}) 0,Fn -( ∇(u -u n ) , {n • (D 2 Φ -Π 0 p D 2 Φ)}) 0,Fn .
Thanks to the smoothness of P p Φ, where P p is as in Theorem 3.8, we can write

B n (u -u n , Φ) = B n (u -u n , Φ -P p Φ).
In fact, L(P p Φ) = 0. Using the continuity property (9), the smoothness and the homogeneous essential boundary conditions of Φ -P p Φ, the approximation properties of P p in [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF], and the stability property [START_REF] Wang | The Morley element for fourth order elliptic equations in any dimensions[END_REF], we deduce

B n (u -u n , Φ) u -u n dG D 2 n (Φ -P p Φ) 0,Ω h p 2 u -u n dG Φ 4,Ω h p 2 u -u n dG u -u n 0,Ω . (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 [START_REF] Wang | The Morley element for fourth order elliptic equations in any dimensions[END_REF] to deduce

R u -u n dG h p 2 Φ 4,Ω h p 2 u -u n dG u -u n 0,Ω . (41)
The assertion follows combining (40) and ( 41) in (39).

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

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 [START_REF] Suri | The p-version of the finite element method for elliptic equations of order 2l[END_REF], Suri constructed an interpolant I BS p of H 2 functions into the space of global H 2 piecewise polynomials on triangular meshes. Such an interpolant satisfies the following approximation properties: there exists a positive constant c only depending on T n such that

u -I BS p u ,Ω ≤ cp -(k-) |u| k,Ω = 0, 1, 2. (42)
The interpolant I BS p is constructed generalizing to the H 2 case the techniques for the H 1 case developed by Babuška and Suri in [START_REF] Babuška | The hp version of the finite element method with quasiuniform meshes[END_REF]. 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 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 ∈ H s+2 (Ω), s := min(k, p -1), p ≥ 2. Then, there exists a positive constant c independent of h and p such that

u -u n dG ≤ c h p s u s+2,Ω , u -u n 0,Ω ≤ c h p 2 u -u n 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 [START_REF] Stein | Singular integrals and differentiability properties of functions[END_REF]. 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 [START_REF] Lederer | Polynomial robust stability analysis for H(div)-conforming finite elements for the Stokes equations[END_REF], 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 [START_REF] Ainsworth | Preconditioning high order H 2 conforming finite elements on triangles[END_REF].

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.

C 0 -IPDG for the biharmonic problem

In this section, we consider C 0 -IPDG for the biharmonic problem; see [START_REF] Brenner | C 0 interior penalty methods for fourth order elliptic boundary value problems on polygonal domains[END_REF][START_REF] Engel | 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[END_REF]. Given the space V n of piecewise continuous polynomials over a mesh T n as in Section 2.1, we introduce the lifting operator L :

V n ⊕ V → [V n ] d×d as Ω L(u n ) : v n := - Fn ({(v n )n} • n) ( ∇u n • n) ∀v n ∈ [V n ] d×d
and the discrete bilinear form

B n (u n , v n ) := Ω D 2 n u n : D 2 n v n + Ω L(u n ) : D 2 n v n + L(v n ) : D 2 n u n + Fn τ ( ∇u n • n)( ∇v n • n) ∀u n , v n ∈ V n . (43)
With an abuse of notation, we denote the norm induced by the bilinear form in (43) by • dG . The bilinear form B n (•, •) 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 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 ∈ H s+2 (Ω), s := min(k, p -1), p ≥ 2. Then, there exists a positive constant c independent of h, and p such that

u -u n dG ≤ c h p s u s+2,Ω , u -u n 0,Ω ≤ c h p 2 u -u n dG .

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 [START_REF] Kanschat | Divergence-conforming discontinuous Galerkin methods and C 0 interior penalty methods[END_REF][START_REF] Wang | A robust numerical method for Stokes equations based on divergence-free H(div) finite element methods[END_REF]. We start with the standard formulation of the Stokes problem in two dimensions: given a polygonal domain

Ω ⊂ R 2 and f ∈ [L 2 (Ω)] 2 , we seek u ∈ [H 1 0 (Ω)] 2 and ρ ∈ L 2 0 (Ω) that solve weakly      -∆ u + ∇ρ = f in Ω, ∇ • u = 0 in Ω, u = 0 on ∂Ω. ( 44
)
The velocity u is divergence free. So, there exists ψ ∈ H 2 0 (Ω) such that u = ∇ × ψ. Introduce the Stokes kernel

Z := v ∈ [H 1 0 (Ω)] 2 div v = 0 .
Testing the first equation in (44) with v ∈ Z gives

(∇ u, ∇ v) 0,Ω = (∇ u, ∇ v) 0,Ω + (∇ρ, v) 0,Ω = (f , v) 0,Ω .
Since v ∈ Z, there exists φ ∈ H 2 0 (Ω) such that v = ∇ × φ. We deduce the stream formulation of the Stokes problem in (44):

(∇ (∇ × ψ), ∇ (∇ × φ)) 0,Ω = (f , ∇ × φ) 0,Ω ∀φ ∈ H 2 0 (Ω). ( 45
)
This is a weak formulation of the biharmonic problem (2) which is alternative to that in [START_REF] Argyris | The TUBA family of plate elements for the matrix displacement method[END_REF]. We can rewrite the right-hand side of (45) integrating by parts:

(f , ∇ × φ) 0,Ω = (∇ × f , φ) 0,Ω .
The strong formulation corresponding to (45) reads as follow: we seek ψ ∈ H 2 0 (Ω) satisfying

∆ 2 ψ = ∇ × f in Ω.
We derive the C 0 -IPDG formulation of this problem using standard dG techniques; see, e.g., [START_REF] Kanschat | Divergence-conforming discontinuous Galerkin methods and C 0 interior penalty methods[END_REF].

On piecewise polynomial spaces, we define the discrete bilinear form

B n (ψ n , φ n ) := Tn (∇ × ∇ψ n ) : (∇ × ∇φ n ) - Fn {∇ × ∇ψ n } • ∇φ n × n - Fn {∇ × ∇φ n } • ∇ψ n × n + Fn τ ∇ψ n × n • ∇φ n × n ∀ψ n , φ n ∈ V n .
With an abuse of notation, we denote the norm induced by the bilinear form in (43) by • 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 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 ψ and ψ n be the solutions to (45) and the corresponding C 0 -IPDG, respectively. Assume that ψ ∈ H s+2 (Ω), s := min(k, p -1), p ≥ 2. Then, there exists a positive constant c independent of h, and p such that

ψ -ψ n dG ≤ c h p s ψ s+2,Ω , ψ -ψ n 0,Ω ≤ c h p 2 ψ -ψ n dG .

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 [START_REF] Babuška | The hp version of the finite element method with quasiuniform meshes[END_REF]; 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. Below, we consider the square domain Ω := (-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 Ω := (-1, 1) 2 \ [0, 1) × (-1, 0] and a given uniform triangular mesh of 6 elements in Section 6.4.

Point singularity inside one element

We consider the exact solution u 1 (x, y) := (x -0.5) 2 + (y -0.5)

2 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-ε (Ω) for any arbitrarily small ε > 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.

Point singularity at an internal mesh vertex

We consider the exact solution

u 2 (x, y) := x 2 + y 2 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-ε (Ω) for any arbitrarily small ε > 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 [START_REF] Babuška | The hp version of the finite element method with quasiuniform meshes[END_REF][START_REF] Suri | The p-version of the finite element method for elliptic equations of order 2l[END_REF].

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 3 2 .
The function u 3 has a point singularity at (1, 0.5), which lies on the boundary of Ω; see the right panel in Figure 1). In particular, it belongs to H 4-ε (Ω) for any arbitrarily small ε > 0. We list the convergence rates under p-refinement in Table 3 and observe a suboptimal by 1.5 order rate O(p -1 2 ) for the DG-norm error. This suboptimality is the counterpart of what is observed in [START_REF] Georgoulis | On the suboptimality of the p-version interior penalty discontinuous Galerkin method[END_REF] for elliptic partial differential equations of second order.

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, θ) := r 5 3 sin 5 3 θ ,
where (r, θ) are the polar coordinates at the re-entrant corner (0, 0). We consider the triangular mesh in Figure 2 Figure 2: Uniform triangular mesh on the L-shape domain.

The function u 4 has a point singularity at (0, 0) and belongs to H 8 3 -ε (Ω) for any arbitrarily small ε > 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 4 3 ) as similarly observed for the square domain in Section 6.3; see Table 4.

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.

A Proof of Lemma 3.3

For the sake of presentation, throughout the proof, we shall write u p instead of P 1 p 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 (ξ), j ∈ N 0 ,

u (ξ) = +∞ i=0 b j L j (ξ), (46)
we define

u p (ξ) := p-2 i=0 b j L j (ξ). ( 47
)
Standard properties of Legendre polynomials, see, e.g., [32, Appendix C], imply

b j = 2j + 1 2 1 -1 u (ξ)L j (ξ)dξ. (48)
For k ≥ 0, we recall [32, Lemma 3.10] that

+∞ i=k 2 2i + 1 (i + k)! (i -k)! |b i | 2 = I (1 -ξ 2 ) k |u (k+2) (ξ)| 2 dξ ≤ u (k+2) 2 0, I . (49)
Using the orthogonality properties of Legendre polynomials [32, eq. (C.24)], the fact that s ≤ p-2, and (49), we obtain

u -u p 2 0, I = +∞ i=p-1 2 2i + 1 |b i | 2 = +∞ i=p-1 2 2i + 1 (i + s)! (i -s)! |b i | 2 (i -s)! (i + s)! ≤ (p -s -1)! (p + s -1)! +∞ i=s 2 2i + 1 (i + s)! (i -s)! |b i | 2 ≤ (p -s -1)! (p + s -1)! u (s+2) 2 0, I .
The two above equations yield the first bound in [START_REF] Brenner | Poincaré-Friedrichs inequalities for piecewise H 2 functions[END_REF]. Next, we introduce

u p (ξ) = ξ -1 u p (t)dt + u (-1). ( 50
)
We have u p (-1) = u (-1). Moreover, recalling that u p and u have the same average over I, we also have

u p (1) = 1 -1 u p (t)dt + u (-1) = 1 -1 u (t)dt + u (-1) = u (1),
which proves [START_REF] Brenner | C 0 interior penalty methods for fourth order elliptic boundary value problems on polygonal domains[END_REF] for the derivative of u p at the endpoints of I.

At this point, we observe

u (ξ) -u p (ξ) = ξ -1 (u (t) -u p (t))dt = +∞ i=p-1 b i ξ -1 L i (t)dt =: +∞ i=p-1 b i φ i (ξ). ( 51
)
Recall the Legendre differential equation [32, eq. (C.2.3)]

((1 -ξ 2 )L i (ξ)) + i(i + 1)L i (ξ) = 0 ∀i ∈ N 0 .
Integrating the above identity over (-1, ξ), ξ ∈ (-1, 1), yields

φ i (ξ) = - 1 i(i + 1) (1 -ξ 2 )L i (ξ). ( 52
)
Recall the orthogonality property of the derivatives of Legendre polynomials [32, eq. (3.39)]:

1 -1 (1 -ξ 2 )L i (ξ)L j (ξ)dξ = 2δ i,j 2i + 1 (i + 1)! (i -1)! . ( 53
)
Combining ( 52) and ( 53), for all i, j ∈ N 0 , we deduce

1 -1 (1 -ξ 2 ) -1 φ i (ξ)φ j (ξ)dξ = 1 i(i + 1)j(j + 1) 1 -1 (1 -ξ 2 )L i (ξ)L j (ξ)dξ = 1 i 2 (i + 1) 2 2δ i,j 2i + 1 (i + 1)! (i -1)! = 2δ i,j i(i + 1)(2i + 1) . (54)
Using s ≤ p -2, we write

1 -1 |u (ξ) -u p (ξ)| 2 dξ = 1 -1 | +∞ i=p-1 b i φ i (ξ)| 2 dξ ≤ 1 -1 (1 -ξ 2 ) -1 | +∞ i=p-1 b i φ i (ξ)| 2 dξ (54) ≤ +∞ i=p-1 2|b i | 2 i(i + 1)(2i + 1) = +∞ i=p-1 2|b i | 2 2i + 1 (i + s)! (i -s)! (i -s)! (i + s)! 1 i(i + 1) ≤ (p -s -1)! (p + s -1)! 1 (p -1)p +∞ i=s 2|b i | 2 2i + 1 (i + s)! (i -s)! (49) ≤ (p -s -1)! (p + s -1)! 1 (p -1)p u (s+2) 2 0, I
, which is the second bound in [START_REF] Brenner | Poincaré-Friedrichs inequalities for piecewise H 2 functions[END_REF]. Finally, we introduce

u p (ξ) = ξ -1 u p (t)dt + u(-1). (55)
We observe that u p (-1) = u(-1). Since L 1 (t) = 1 = L 0 (t), standard manipulations imply

u p (1) -u p (-1) = 1 -1 u p (t)dt (50) = 2u (-1) + 1 -1 t -1 u p (x)dx dt (46),(51) = 2u (-1) + 1 -1 p-2 i=0 b i φ i (t)dt = 2u (-1) + p-2 i=0 b i 1 -1 φ i (t)dt (52) = 2u (-1) - p-1 i=0 b i 1 -1 1 i(i + 1) (1 -t 2 )L i (t)L 1 (t)dt (53) = 2u (-1) - +∞ i=0 b i 1 -1 1 i(i + 1) (1 -t 2 )L i (t)L 1 (t)dt (52),(46) = 2u (-1) + 1 -1 t -1 u (x)dx dt = 1 -1 u (t)dt = u(1) -u(-1).
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(ξ) -u p (ξ) = ξ -1
(u (t) -u p (t))dt.

We arrive at u(ξ) -u p (ξ)

(51) = +∞ i=p-1 b i ξ -1 φ i (x) dx dt =: +∞ i=p-1 b i ψ i (ξ). ( 56
)
We prove certain orthogonality properties of the ψ i functions. Recall the identity [32, eq. (C.2.5)]:

L i (ξ) = L i+1 (ξ) -L i-1 (ξ) 2i + 1 ∀i ≥ 2, L 0 (ξ) = L 1 (ξ).
Integrating over (-1, t), t ∈ (-1, 1), and using L i+1 (-1) = L i-1 (-1), see [32, eq. (C.2.6)], we deduce

φ i (t) = L i+1 (ξ) -L i-1 (ξ) 2i + 1 .
Upon integrating the above identity over (-1, ξ), ξ ∈ (-1, 1), we arrive at

ψ i (ξ) = ξ -1 L i+1 (t) dt - ξ -1 L i-1 (t) dt 2i + 1 (51) = φ i+1 (ξ) -φ i-1 (t) 2i + 1 .
With the notation φ -1 (ξ) = 0, we get

ψ i (ξ)ψ j (ξ) = 1 (2i + 1)(2j + 1) (φ i+1 (ξ)φ j+1 (ξ) + φ i-1 (ξ)φ j-1 (ξ) -φ i+1 (ξ)φ j-1 (ξ) -φ i-1 (ξ)φ j+1 (ξ)) .
In (54), we proved that the φ i functions are orthogonal with respect to the (1 -ξ 2 )-weighted L 2 inner product. Therefore, testing the above identity by (1 -ξ 2 ) -1 and integrating over (-1, 1) we arrive at

1 -1 (1 -ξ 2 ) -1 ψ i (ξ)ψ j (ξ) dξ = 1 (2i + 1)(2j + 1) (ℵ i+1,j+1 + ℵ i-1,j-1 -ℵ i+1,j-1 -ℵ i-1,j+1 ) , (57)
where, from (54),

ℵ ,k = 0 if ≤ 0 or k ≤ 0, 2δ ,k ( +1)(2 +1) if , k > 0.
We provide explicit values for the expression in (57). If j ∈ {i -2, i, i + 2}, then

1 -1 (1 -ξ 2 ) -1 ψ i (ξ)ψ j (ξ) dξ = 0. If j = i, then 1 -1 (1 -ξ 2 ) -1 ψ i (ξ)ψ i (ξ) dξ = i 1 + i 2 ,
where

i 1 := 1 (2i + 1) 2 2 (i + 1)(i + 2)(2i + 3) ∀i ∈ N; i 2 := 0 if i = 1, 1 (2i+1) 2 2 (i-1)i(2i-1) if i ≥ 2. ( 58
) If j = i + 2, then 1 -1 (1 -ξ 2 ) -1 ψ i (ξ)ψ i+2 (ξ) dξ = - 1 (2i + 1)(2i + 5) • 2 (i + 1)(i + 2)(2i + 3) =: -i 3 . ( 59
) If j = i -2, then 1 -1 (1 -ξ 2 ) -1 ψ i (ξ)ψ i-2 (ξ) dξ = -i 4 := - 1 (2i+1)(2i-3) • 2 (i-1)i(2i-1) if i ≥ 3, 0 if i ≤ 2. ( 60
)
In light of the above orthogonality properties, we can estimate the L 2 approximation error as follows:

u -u p 2 0, I (56)
= +∞ i=p-1 b i ψ i 2 0, I ≤ +∞ i,j=p-1 b i b j 1 -1
(1 -ξ 2 ) -1 ψ i (ξ)ψ j (ξ)dξ (58),( 59),( 60)

≤ +∞ i=p-1 |b i | 2 i 1 + i 2 + i 3 + i 4 .
We cope with the four terms on the right-hand side separately. We begin with i 1 :

+∞ i=p-1 |b i | 2 i 1 (58) = +∞ i=p-1 |b i | 2 1 (2i + 1) 2 • 2 (i + 1)(i + 2)(2i + 3) = +∞ i=p-1 2 2i + 1 |b i | 2 (i + s)! (i -s)! (i -s)! (i + s)! 1 (i + 1)(i + 2)(2i + 1)(2i + 3) (49) ≤ (p -s -1)! (p + s -1)! 1 p(p + 1)(2p -1)(2p + 1) u (s+2) 2 0, I . Next, we focus on i 2 , i ≥ 2: +∞ i=max(2,p-1) |b i | 2 i 2 (58) = +∞ i=max(2,p-1) |b i | 2 1 (2i + 1) 2 • 2 (i -1)i(2i -1) = +∞ i=max(2,p-1) 2 2i + 1 |b i | 2 (i + s)! (i -s)! (i -s)! (i + s)! 1 (i -1)i(2i -1)(2i + 1) (49) ≤ (p -s -1)! (p + s -1)! 1 (p -2)(p -1)(2p -3)(2p -1) u (s+2) 2 0, I .
As for the term i 3 , we proceed as follows:

+∞ i=p-1 |b i | 2 i 3 (59) = +∞ i=p-1 |b i | 2 1 (2i + 1)(2i + 5) • 2 (i + 1)(i + 2)(2i + 3) = +∞ i=p-1 2 2i + 1 |b i | 2 (i + s)! (i -s)! (i -s)! (i + s)! 1 (i + 1)(i + 2)(2i + 3)(2i + 5) (49) ≤ (p -s -1)! (p + s -1)! 1 p(p + 2)(2p + 1)(2p + 3) u (s+2) 2 0, I .
Eventually, we cope with the term i 4 , i ≥ 3:

+∞ i=max(3,p-1) |b i | 2 i 4 (60) = +∞ i=max(3,p-1) |b i | 2 1 (2i + 1)(2i -3) • 2 (i -1)i(2i -1) = +∞ i=p-1 2 2i + 1 |b i | 2 (i + s)! (i -s)! (i -s)! (i + s)! 1 (i -1)i(2i -3)(2i -1) ≤ (p -s -1)! (p + s -1)! 1 (p -2)(p -1)(2p -5)(2p -3) u (s+2) 2 0, I .
Collecting the five bounds above concludes the proof.

B Proof of Theorem 3.6 (2D case)

The continuity properties [START_REF] Dong | Residual-based a posteriori error estimates for hp-discontinuous Galerkin discretizations of the biharmonic problem[END_REF] 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 x p and P y p are the projections in Lemma 3.3 along the x and y directions, respectively. Some identities. The following identities are valid:

∂ y P x p u = P x p (∂ y u), ∂ x P y p u = P y p (∂ x u). (61)
We only show the first one as the second can be proven analogously. For all y ∈ (-1, 1), after writing P x p u(x, y) in integral form with respect to the x variable, we have

P x p u(x, y) = x -1 t -1 ∂ 2 x P x p u(x, y)dx + ∂ x u(-1, y) dt + u(-1, y) (11) = x -1 t -1 ∂ 2 x P x p u(x, y)dx + ∂ x u(-1, y) dt + u(-1, y) (47),(48) = x -1 t -1 ∂ 2 x p-2 i=0 2i + 1 2 1 -1 u(s, y)L i (s)ds L i (x)dx + ∂ x u(-1, y) dt + u(-1, y).
Taking the partial derivative with respect to y on both sides yields

∂ y P x p u(x, y) = x -1 t -1 ∂ 2 x p-2 i=0 2i + 1 2 1 -1 ∂ y u(s, y)L i (s)ds L i (x)dx + ∂ x ∂ y u(-1, y) dt + ∂ y u(-1, y) (47)
,( 50),(55) = P x p ∂ y u(x, y).

L 2 estimates. Using the definition of P p , the triangle inequality, the one dimensional approximation properties [START_REF] Canuto | Approximation results for orthogonal polynomials in Sobolev spaces[END_REF], the third stability property in [START_REF] Dong | Discontinuous Galerkin methods for the biharmonic problem on polygonal and polyhedral meshes[END_REF], and the identities (61) we write

u -P p u 0, Q ≤ u -P x p u 0, Q + P x p (u -P y p u) 0, Q p -s-2 ∂ s+2 x u 0, Q + u -P y p u 0, Q + p -2 ∂ 2 x (u -P y p u) 0, Q = p -s-2 ∂ s+2 x u 0, Q + u -P y p u 0, Q + p -2 ∂ 2 x u -P y p ∂ 2 x u 0, Q .
The assertion follows using again the one dimensional approximation properties [START_REF] Canuto | Approximation results for orthogonal polynomials in Sobolev spaces[END_REF].

H 1 estimates. First, we cope with the bound on the derivative with respect to x. Adding and subtracting P x p u, and using the triangle inequality yield ∂

x (u -P p u) 0, Q ≤ ∂ x (u -P x p u) 0, Q + ∂ x P x p (u -P p u) 0, Q =: T 1 + T 2 . ( 62
)
As for the term T 1 , we use the one dimensional approximation properties (13) and get

T 1 p -s-1 ∂ s+2 x u 0, Q . (63)
As for the term T 2 , we use the second stability estimate in ( 14) and get

T 2 ∂ x (u -P y p u) 0, Q + p -1 ∂ 2 x (u -P y p u) 0, 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 ∂ x u -P y p ∂ x u 0, Q + p -1 ∂ 2 x u -P y p ∂ 2 x u 0, Q p -s-1 ∂ x ∂ s+1 y u 0, Q + ∂ 2 x ∂ s y u 0, Q . ( 64
)
Collecting the estimates (63) and ( 64) in (62), we arrive at

∂ x (u -P p u) 0, Q p -s-1 ∂ s+2 x u 0, Q + ∂ x ∂ s+1 y u 0, Q + ∂ 2 x ∂ s y u 0, Q .
With similar arguments for the y derivative term, we deduce [START_REF] Engel | 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[END_REF].

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 [START_REF] Dong | Discontinuous Galerkin methods for the biharmonic problem on polygonal and polyhedral meshes[END_REF], and the identities (61), we obtain

∂ 2 x (u -P p u) 0, Q ≤ ∂ 2 x (u -P x p u) 0, Q + ∂ 2 x P x p (u -P y p u) 0, Q p -s ∂ s+2 x u 0, Q + ∂ 2 x (u -P y p u) 0, Q = p -s ∂ s+2 x u 0, Q + ∂ 2 x u -P y p ∂ 2 x u 0, Q p -s ∂ s+2 x u 0, Q + ∂ 2 x ∂ s y u 0, Q .
Analogously, we can prove

∂ 2 y (u -P p u) 0, Q p -s ∂ s+2 y u 0, Q + ∂ s x ∂ 2 y u 0, Q .
Eventually, we cope with the mixed derivative term. Using the triangle inequality and the identities (61), we get

∂ x ∂ y (u -P p u) 0, Q ≤ ∂ x ∂ y (u -P x p u) 0, Q + ∂ x ∂ y P x p (u -P y p u) 0, Q = ∂ x (∂ y u -P x p ∂ y u) 0, Q + ∂ x P x p ∂ y (u -P y p u) 0, Q =: S 1 + S 2 . (65)
We estimate the two terms S 1 and S 2 separately. Using the one dimensional approximation properties (13), we can write

S 1 p -s ∂ s+1 x ∂ y u 0, Q . (66)
On the other hand, using the one dimensional approximation properties ( 13) and the stability properties ( 14), we get

S 2 ∂ x ∂ y (u -P y p y) 0, Q + p -1 ∂ 2 x ∂ y (u -P y p u) 0, Q = ∂ y (∂ x u -P y p ∂ x y) 0, Q + p -1 ∂ y (∂ 2 x u -P y p ∂ 2 x u) 0, Q p -s ∂ x ∂ s+1 y u 0, Q + ∂ 2 x ∂ s y u 0, Q . (67)
Collecting the estimates (66) and (67) in (65), we arrive at

∂ x ∂ y (u -P p u) 0, Q p -s ∂ s+1 x ∂ y u 0, Q + ∂ x ∂ s+1 y u 0, Q + ∂ 2 x ∂ s y u 0, Q .
Combining the estimates on all second derivatives terms, we obtain [START_REF] Georgoulis | Discontinuous Galerkin methods on shape-regular and anisotropic meshes[END_REF].

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 P x p , P y p , and P z p 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:

∂ y P x p u = P x p (∂ y u), ∂ x P y p u = P y p (∂ x u), ∂ y P z p u = P z p (∂ y u), ∂ z P y p u = P y p (∂ z u), ∂ x P z p u = P z p (∂ x u), ∂ z P x p u = P x p (∂ z u). (68)
L 2 estimates. The triangle inequality and the one dimensional approximation properties (13) imply u -P p u 0, Q ≤ u -P x p u 0, Q + P x p (u -P y p P z p u) 0, Q p -s-2 ∂ s+2 x u 0, Q + P x p (u -P y p P z p u) 0, Q .

We focus on the second term. To this aim, we use the stability properties ( 14), the triangle inequality, the one dimensional approximation properties [START_REF] Canuto | Approximation results for orthogonal polynomials in Sobolev spaces[END_REF], and the identities (68), and deduce P x p (u -P y p P z p u) 0, Q ≤ u -P y p P z p u 0, Q + p -2 ∂ 2 x u -P y p P z p ∂ 2 x u 0, Q u -P y p u 0, Q + P y p (u -P z p u) 0, Q + p -2 ∂ 2 x u -P y p ∂ 2 x u 0, Q + p -2 P y p (∂ 2 x u -P z p ∂ 2 x u) 0, Q p -s-2 ∂ s+2 y u 0, Q + u -P z p u 0, Q + p -2 ∂ 2 y u -P z p ∂ 2 y u 0,

Q + p -s-4 ∂ 2 x ∂ 2 y ∂ s-2 z u 0, Q + p -2 ∂ 2 x u -P z p ∂ 2 x u 0, Q + p -4 ∂ 2 x ∂ 2 y u -P z p ∂ 2 x ∂ 2 y u 0, Q .
The L 2 estimates eventually follow from the one dimensional approximation properties [START_REF] Canuto | Approximation results for orthogonal polynomials in Sobolev spaces[END_REF].

H 1 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 ∂ x (u -P p u) 0, Q ∂ x (u -P x p u) 0, Q + ∂ x P x p (u -P y p P z p u) 0, Q p -s-1 ∂ s+2

x u 0, Q + ∂ x P x p (u -P y p P z p u) 0, Q .

We focus on the second term on the right-hand side. The stability properties ( 14) and the identities (68) entail ∂ x P x p (u -P y p P z p u) 0, Q ∂ x u -P y p P z p ∂ x u 0, Q + p -2 ∂ 2 x u -P y p P z p ∂ 2 x u 0, 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 ∂ x u-P y p ∂ x u 0, Q + P y p (∂ x u-P z p ∂ x u) 0, Q p -s-1 ∂ x ∂ s+1 y u 0, Q + P y p (∂ x u-P z p ∂ x u) 0, Q .
The second term on the right-hand side can be estimated using the stability properties [START_REF] Dong | Discontinuous Galerkin methods for the biharmonic problem on polygonal and polyhedral meshes[END_REF], the identities (68), and the one dimensional approximation properties (13):

P y p (∂ x u-P z p ∂ x u) 0, Q ∂ x u-P z p ∂ x u 0, Q + p -2 ∂ x ∂ 2 y u-P z p ∂ x ∂ 2 y u) 0, Q p -s-1 ∂ x ∂ s+1 z u 0, Q + ∂ x ∂ 2 y ∂ s-1 x u 0, Q .
With similar arguments based on substituting ∂ x u by ∂ 2 x u, we find an upper bound for T 2 :

T 2 p -s-1 ∂ 2 x ∂ s y u 0, Q + ∂ 2 x ∂ s z u 0, Q + ∂ 2 x ∂ 2 y ∂ s-2 z u 0, Q .
Collecting the above estimates, we arrive at

∂ x (u -P p u) 0, Q p -s-1 ∂ s+2 x u 0, Q + ∂ x ∂ s+1 y u 0, Q + ∂ x ∂ s+1 z u 0, Q + ∂ x ∂ 2 y ∂ s-1 z u 0, Q + ∂ 2 x ∂ s y u 0, Q + ∂ 2 x ∂ s z u 0, Q + ∂ 2 x ∂ 2 y ∂ s-2 z u 0, Q .
By similar arguments on the y and x partial derivatives, we deduce [START_REF] Georgoulis | Discontinuous Galerkin methods for the biharmonic problem[END_REF].

H 2 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 [START_REF] Canuto | Approximation results for orthogonal polynomials in Sobolev spaces[END_REF], and the identities (68) imply ∂ 2

x (u -P p u) 0, Q ≤ ∂ 2

x (u -P x p u) 0, Q + ∂ 2 x P x p (u -P y p P z p u) 0, Q p -s ∂ s+2

x u 0, Q + ∂ 2 x u -P y p P z p ∂ 2 x u 0, Q .

We are left with estimating the second term on the right-hand side. Applying a further triangle inequality, the one dimensional approximation properties [START_REF] Canuto | Approximation results for orthogonal polynomials in Sobolev spaces[END_REF], and the stability properties [START_REF] Dong | Discontinuous Galerkin methods for the biharmonic problem on polygonal and polyhedral meshes[END_REF], we arrive at

∂ 2 x u -P y p P z p ∂ 2 x u 0, Q ≤ ∂ 2 x u -P y p ∂ 2 x u 0, Q + P y p (∂ 2 x u -P z p ∂ 2 x u) 0, Q p -s ∂ 2 x ∂ s y u 0, Q + ∂ 2 x u -P z p ∂ 2 x u 0, Q + p -2 ∂ 2 x ∂ 2 y u -P z p ∂ 2 x ∂ 2 y u 0, Q p -s ∂ 2 x ∂ s y u 0, Q + ∂ 2 x ∂ s z u 0, Q + ∂ 2 x ∂ 2 y ∂ s-2 z u 0, Q .
Collecting the two above estimate above gives

∂ 2 x (u -P p u) 0, Q p -s ∂ s+2 x u 0, Q + ∂ 2 x ∂ s y u 0, Q + ∂ 2 x ∂ s z u 0, Q + ∂ 2 x ∂ 2 y ∂ s-2 z u 0, Q . (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 [START_REF] Canuto | Approximation results for orthogonal polynomials in Sobolev spaces[END_REF], the identities (68), and the stability properties ( 14) leads us to

∂ x ∂ y u -P p ∂ x ∂ y u 0, Q ≤ ∂ x ∂ y u -P x p ∂ x ∂ y u 0, Q + P x p (∂ x ∂ y u -P y p P z p ∂ x ∂ y u) 0, Q p -s ∂ s+1 x ∂ y u 0, Q + ∂ y ∂ x u -∂ y P y p P z p ∂ x u 0, Q + p -1 ∂ y ∂ 2 x u -∂ y P y p P z p ∂ 2 x u 0, Q = p -s ∂ s+1
x ∂ y u 0, Q + T 3 + T 4 .

We estimate the terms A and B separately. First, we focus on A. Using the triangle inequality, the one dimensional approximation properties [START_REF] Canuto | Approximation results for orthogonal polynomials in Sobolev spaces[END_REF], the identities (68), and the stability properties ( 14), entail

T 3 ≤ ∂ y ∂ x u -∂ y P y p ∂ x u 0, Q + ∂ y P y p ∂ x u -∂ y P y p ∂ x u 0, Q p -s ∂ x ∂ s+1 y u 0, Q + ∂ y ∂ x u -P z p ∂ y ∂ x u 0, Q + p -1 ∂ 2 y ∂ x u -P z p ∂ 2 y ∂ x u 0, Q p -s ∂ x ∂ s+1 y u 0, Q + ∂ x ∂ y ∂ s z u 0, Q + ∂ x ∂ 2 y ∂ s-1 z u 0, Q .
Next, we focus on the term B. Using the triangle inequality, the one dimensional approximation properties [START_REF] Canuto | Approximation results for orthogonal polynomials in Sobolev spaces[END_REF], the identities (68), and the stability properties [START_REF] Dong | Discontinuous Galerkin methods for the biharmonic problem on polygonal and polyhedral meshes[END_REF], leads to

T 4 ≤ p -1 ∂ y ∂
p -s ∂ 2 x ∂ s y 0, Q + ∂ 2 x ∂ y ∂ s z 0, Q + ∂ 2 x ∂ 2 y ∂ s-2 z 0, Q .
Collecting the above estimates gives

∂ x ∂ y u -P p ∂ x ∂ y u 0, Q p -s ∂ s+1 x ∂ y u 0, Q + ∂ x ∂ s+1 y u 0, Q + ∂ x ∂ y ∂ s z u 0, Q ∂ x ∂ 2 y ∂ s-1 z u 0, Q + ∂ 2 x ∂ s y u 0, Q + ∂ 2 x ∂ y ∂ s-1 z u 0, Q + ∂ 2 x ∂ 2 y ∂ s-2 z u 0, Q . (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.

Theorem 3 . 8 .

 38 Let T n be a tensor product-type mesh over Ω ⊂ R d , d = 2, 3, of either parallelograms (d = 2) or parallelepipeds (d = 3) as in Section 2.1; k ≥ 0; p ∈ N; s := min(p, k -1). Then there exists

Figure 1 :

 1 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.

Table 1 p

 1 -convergence for the exact solution u 1 .

	p	u -u n dG Convergence rate	p	u -u n dG Convergence rate
	2	1.46E+01		3	1.18E+01
	4	5.07E+00	1.52	5	9.67E-01	4.90
	6	1.42E-01	8.83	7	4.19E-02	9.33
	8	2.38E-02	6.18	9	1.62E-02	3.79
	10	1.11E-02	3.41	11	9.29E-03	2.77
	12	6.76E-03	2.73	13	6.12E-03	2.49
	14	4.67E-03	2.40	15	4.39E-03	2.33
	16	3.48E-03	2.21	17	3.33E-03	2.20
	18	2.73E-03	2.05	19	3.33E-03	2.09
	20	2.23E-03	1.93	21	2.16E-03	1.99
	22	1.84E-03	2.03	23	1.80E-03	2.01
	24	1.55E-03	1.96	25	1.53E-03	1.92

Table 2 p

 2 -convergence for the exact solution u 2 .

	p	u -u n dG Convergence rate	p	u -u n dG Convergence rate
	2	5.16E+00		3	4.39E+00
	4	2.98E+00	0.79	5	9.61E-01	2.97
	6	8.64E-02	8.74	7	1.08E-02	13.33
	8	6.22E-03	9.14	9	3.84E-03	4.14
	10	2.47E-03	4.14	11	1.67E-03	4.13
	12	1.18E-03	4.07	13	8.52E-04	4.04
	14	6.34E-04	3.99	15	4.84E-04	3.96
	16	3.76E-04	3.92	17	2.90E-04	4.09
	18	2.31E-04	4.12	19	1.87E-04	3.92
	20	1.53E-04	3.90	21	1.27E-04	3.86
	22	1.05E-04	3.95	23	8.82E-05	4.01
	24	7.44E-05	3.99	25	6.28E-05	4.06

Table 3 p

 3 -convergence for the exact solution u 3 .

	p	u -u n dG Convergence rate	p	u -u n dG Convergence rate
	2	8.61E+00		3	2.57E+00
	4	6.39E-01	3.75	5	1.01E+00	1.83
	6	4.90E-01	0.65	7	7.48E-01	0.89
	8	4.33E-01	0.44	9	6.37E-01	0.64
	10	3.91E-01	0.46	11	5.71E-01	0.55
	12	3.56E-01	0.52	13	5.24E-01	0.51
	14	3.31E-01	0.46	15	4.88E-01	0.50
	16	3.13E-01	0.42	17	4.59E-01	0.49
	18	2.97E-01	0.47	19	4.36E-01	0.48
	20	2.82E-01	0.48	21	4.15E-01	0.48
	22	2.69E-01	0.49	23	3.97E-01	0.49
	24	2.58E-01	0.49	25	3.81E-01	0.49

Table 4 p

 4 -convergence for the exact solution u 1 .

	p	u -u n dG Convergence rate	p	u -u n dG Convergence rate
	2	1.80E+00		3	9.38E-01
	4	6.93E-01	1.82	5	5.50E-01	1.35
	6	4.53E-01	1.27	7	3.84E-01	1.26
	8	3.32E-01	1.25	9	2.91E-01	1.25
	10	2.58E-01	1.26	11	2.32E-01	1.26
	12	2.09E-01	1.27	13	1.91E-01	1.27
	14	1.75E-01	1.27	15	1.62E-01	1.28
	16	1.50E-01	1.28	17	1.39E-01	1.29
	18	1.30E-01	1.29	19	1.22E-01	1.28
	20	1.15E-01	1.29	21	1.08E-01	1.30
	22	1.02E-01	1.31	23	9.60E-02	1.31
	24	9.22E-02	1.31	25	8.77E-02	1.32

  2 x u -∂ y P y p ∂ 2 x u 0, Q + p -1 ∂ y P y p ∂ 2 x u -∂ y P y p P z p ∂ 2 x u 0, Q p -s ∂ 2 x ∂ s y 0, Q + p -1 ∂ y ∂ 2 x u -P z p ∂ y ∂ 2 x u 0, Q + p -2 ∂ 2

	x ∂ 2 y u -P z p ∂ 2 x ∂ 2 y u 0, Q