# $C^0$-hybrid high-order methods for biharmonic problems

Zhaonan Dong, Alexandre Ern

► To cite this version:

Zhaonan Dong, Alexandre Ern. $C^0$-hybrid high-order methods for biharmonic problems. IMA Journal of Numerical Analysis, 2024, 44 (1), pp.24-57. 10.1093/imanum/drad003. hal-03695484v3

HAL Id: hal-03695484

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

Submitted on 11 Feb 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
# $C^0$-hybrid high-order methods for biharmonic problems

Zhaonan Dong*

Alexandre Ern†

February 11, 2023

## Abstract

We devise and analyze $C^0$-conforming hybrid high-order (HHO) methods to approximate biharmonic problems with either clamped or simply supported boundary conditions. $C^0$-conforming HHO methods hinge on cell unknowns which are $C^0$-conforming polynomials of order $(k+2)$ approximating the solution in the mesh cells and on face unknowns which are polynomials of order $k \ge 0$ approximating the normal derivative of the solution on the mesh skeleton. Such methods deliver $O(h^{k+1}) H^2$-error estimates for smooth solutions. An important novelty in the error analysis is to lower the minimal regularity requirement on the exact solution. The technique to achieve this has a broader applicability than just $C^0$-conforming HHO methods, and to illustrate this point, we outline the error analysis for the well-known $C^0$-conforming interior penalty discontinuous Galerkin (IPDG) methods as well. The present technique does not require bubble functions or a $C^1$-smoother to evaluate the right-hand side in case of rough loads. Finally, numerical results including comparisons to various existing methods showcase the efficiency of the proposed $C^0$-conforming HHO methods.

## 1 Introduction

Biharmonic PDEs are used in the modelling of various physical phenomena, such as thin plate elasticity, micro-electromechanical systems, and phase separation, to mention a few examples. In the present work, we consider the following model problem with two types of boundary conditions (BC's):

$$ \Delta^2 u = f \quad \text{in } \Omega, \quad \begin{cases} u = \partial_n u = 0 & \text{type (I)} \\ u = \partial_{nn} u = 0 & \text{type (II)} \end{cases} \quad \text{on } \partial\Omega. \quad (1) $$

Here, $\Omega$ is an open bounded Lipschitz domain in $\mathbb{R}^d$, $d \in \{2, 3\}$, with boundary $\partial\Omega$, $\partial_n$ denotes the normal derivative on $\partial\Omega$ and $\partial_{nn}$ denotes the normal-normal component of the Hessian on $\partial\Omega$. In the context of plate modelling, type (I) BC is referred to as *clamped* BC, and type (II) as *simply supported* BC. Non-homogeneous BC's can be considered, but we focus on homogeneous BC's for simplicity. Instead, dealing with the boundary condition $\partial_{nn}u = 0$ requires further developments. The regularity of the source term $f$ is specified below in Assumption 1.1. We also observe that the present developments hinge on the weak formulation of (1) involving the Hessian.

The goal of the present work is twofold. The first main goal is to devise and analyze a novel $C^0$-conforming approximation method for the above model problem. The proposed method belongs to the class of hybrid high-order (HHO) methods. These methods were

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

<footnote>†CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée, France, and Inria, 2 rue Simone Iff, 75589 Paris, France.
alexandre.ern@enpc.fr.</footnote>
introduced in [11] for linear diffusion and in [10] for locking-free linear elasticity. Since
then, they have undergone a vigorous development, as reflected, e.g., in the two recent
monographs [9, 7]. Moreover, as discussed in [8, 24, 5], HHO methods are closely related
to hybridizable discontinuous Galerkin (HDG) methods, weak Galerkin (WG) methods,
nonconforming virtual element methods (ncVEM), and multiscale hybrid-mixed (MHM)
methods. Interestingly, the present C⁰-conforming HHO method (C⁰-HHO in short) can be
viewed as a simple approach to extend C⁰-finite element methods for biharmonic problems by
simply adding an additional unknown attached to the mesh faces representing the normal
derivative of the solution. C⁰-HHO methods for the biharmonic problem have not yet
been explored in the literature; we refer the reader to [1, 13, 12] for fully discontinuous
HHO methods. The starting point for devising the present C⁰-HHO methods are the HHO
methods from [13].

$C^0$-conforming approximation methods are popular to discretize biharmonic problems since such methods avoid the intricate construction of $C^1$-conforming approximation spaces while also avoiding severe conditioning issues that can arise with fully nonconforming approximation methods. Examples of $C^0$-conforming approximation methods from the literature include the classical Morley [25, 30] and Hsieh-Clough-Tocher (HCT) finite element methods, the $C^0$-interior penalty discontinuous Galerkin (IPDG) method [15, 2], the $C^0$-weak Galerkin (WG) method [26, 6], and the $C^0$-virtual element method (VEM) [31]. The Morley FEM and the HCT FEM for $d = 3$ are lowest-order methods, whereas the HCT FEM for $d = 2$ [14] and the $C^0$-{IPDG,WG,HHO,VEM} for $d \ge 2$ can reach arbitrary approximation order. $C^0$-IPDG attaches discrete unknowns to the mesh cells only, whereas $C^0$-{WG,HHO} attach discrete unknowns to the mesh cells and faces. A numerical comparison between $C^0$-IPDG and $C^0$-HHO is included herein, indicating the computational efficiency of the $C^0$-HHO approach. Concerning $C^0$-WG, we recall that the only relevant difference between HHO and WG lies in the choice of the discrete unknowns and the design of the stabilization operator. In general, WG employs plain least-squares stabilization, leading to suboptimal convergence rates, whereas HHO employs a more elaborate form of stabilization leading to optimal convergence rates. In other words, to achieve the same convergence rate, $C^0$-WG with plain least-squares stabilization requires more discrete unknowns than $C^0$-HHO methods. Furthermore, among the above $C^0$-conforming methods, $C^0$-VEM is currently the only one supporting general (polytopal) meshes. We notice that an interesting perspective to the present work in order to devise a $C^0$-HHO method on general meshes is to use $C^0$-VEM to build the cell basis functions, while keeping the current face basis functions to handle the normal derivative; this perspective is left to future work. Finally, $C^0$-conforming HDG methods have been so far seldom considered in the literature, with the exception of a brief discussion in [20]. Therefore, leveraging on [8], the present work can be viewed as a contribution to the development of $C^0$-HDG methods for the biharmonic problem.

The second main goal of the present work is to improve on the regularity requirement
on the exact solution to lead the error analysis. The technique to achieve this is a non-
trivial extension of ideas from [19] originally developed in the context of second-order elliptic
problems (see also [18, Chap. 40 & 41]). The result derived herein can be applied to the
$C^0$-HHO method when approximating the biharmonic problem, but also, more broadly, to
the other $C^0$-conforming approximation methods discussed above. Recall that the difficulty
comes from the lack of $C^1$-conformity which causes some difficulties in the error analysis
when it comes to bounding the consistency error. A first possibility is to use bubble functions
together with a $C^1$-smoother when evaluating the right-hand side of the discrete problem,
as shown in [29, 28] for the Morley element and the $C^0$-IPDG method (see also [4] for
further results concerning lowest-order methods). However, these techniques so far meet
with difficulties when it comes to devising a $C^1$-smoother of arbitrary order on tetrahedral
meshes. The alternative road, which is the one followed by most of the above works and
also herein, is to require some (mild) additional regularity assumption on the source term
beyond $H^{-2}(\Omega)$ and on the exact solution beyond $H^2(\Omega)$. The present work hinges on the
following original assumption.
**Assumption 1.1 (Regularity)** We assume that there are real numbers $p > 2$ and $q \in (\frac{2d}{d+2}, 2]$ such that $u \in W^{2,p}(\Omega)$, $\Delta u \in W^{1,q}(\Omega)$, and $f \in W^{-1,q}(\Omega)$.

Assumption 1.1 is, to our knowledge, novel in the analysis of biharmonic problems. If one prefers to remain in the Hilbertian setting, a simpler, but less general, assumption is $u \in H^{2+s}(\Omega)$, with $s > 0$, $\Delta u \in H^1(\Omega)$, and $f \in H^{-1}(\Omega)$. Notice that the assumption $u \in H^{2+s}(\Omega)$, $s > 0$, is still tighter than the one usually made in the literature to analyze $C^0$-conforming methods, which is $u \in H^{2+s}(\Omega)$, $s > \frac{1}{2}$; the assumption $s > 0$ is considered in [4]. Finally, we observe that in Assumption 1.1, the Laplacian of the solution is smoother than its Hessian. Since the Laplace operator contains a relatively large kernel, this assumption is reasonable.

The rest of this work is organized as follows. We present the key identities to bound the consistency error in Section 2. The results presented in this section have a wider outreach beyond $C^0$-HHO methods since they can be applied to other $C^0$-conforming methods (but not to fully nonconforming methods). Then, we devise $C^0$-HHO methods for both types of BC's and establish stability and well-posedness of the discrete problems in Section 3. In Section 4, we perform the error analysis under Assumption 1.1. For completeness, we also outline how the $C^0$-IPDG method can be analyzed under this assumption. Finally, we discuss our numerical results in Section 5.

# 2 Key identities to bound the consistency error

In this section, we introduce some basic notation, present the weak formulation for both types of BC's, and derive the key identities to bound the consistency error. The main result of this section, Lemma 2.3, can be applied to many $C^0$-conforming methods.

## 2.1 Basic notation

We use standard notation for the Lebesgue and Sobolev spaces. In particular, for the fractional-order Sobolev spaces, we consider the Sobolev-Slobodeckij seminorm based on the double integral. For an open, bounded, Lipschitz set $S$ in $\mathbb{R}^d$, $d \in \{1, 2, 3\}$, we denote by $(v, w)_S$ the $L^2(S)$-inner product with appropriate Lebesgue measure, and we employ the same notation for vector- or matrix-valued fields (such fields are denoted using boldface notation). We denote by $\nabla w$ the (weak) gradient of $w$ and by $\nabla^2 w$ its (weak) Hessian. Let $n_S$ be the unit outward normal vector on the boundary $\partial S$ of $S$. Assuming that the functions $v$ and $w$ are smooth enough, we have the following integration by parts formula:

$$ (\Delta^2 v, w)_S = (\nabla^2 v, \nabla^2 w)_S + (\nabla \Delta v, n w)_{\partial S} - (\nabla^2 v n, \nabla w)_{\partial S}. \quad (2) $$

Whenever the context is unambiguous, we denote by $\partial_n$ the (scalar-valued) normal derivative on $\partial S$ and by $\partial_t$ the $(\mathbb{R}^{d-1}$-valued) tangential derivative. We also denote by $\partial_{nn} v$ the (scalar-valued) normal-normal second-order derivative, and by $\partial_{nt} v$ the $(\mathbb{R}^{d-1}$-valued) normal-tangential second-order derivative. The integration by parts formula (2) can then be rewritten as

$$ (\Delta^2 v, w)_S = (\nabla^2 v, \nabla^2 w)_S + (\partial_n \Delta v, w)_{\partial S} - (\partial_{nn} v, \partial_n w)_{\partial S} - (\partial_{nt} v, \partial_t w)_{\partial S}. \quad (3) $$

In what follows, the set $S$ is always a polytope so that its boundary can be decomposed into a finite union of planar faces with disjoint interiors. Expressions involving the tangential derivative on $\partial S$ are then implicitly understood to be evaluated as a summation over the faces composing $\partial S$.

Let $\{\mathcal{T}\}_{h>0}$ be a shape-regular family of simplicial meshes such that each mesh covers the domain $\Omega$ exactly. A generic mesh cell is denoted by $K \in \mathcal{T}$, its diameter by $h_K$, and its unit outward normal by $n_K$. We partition the boundary $\partial K$ of any mesh cell $K \in \mathcal{T}$ by means of the two subsets $\partial K^i := \partial K \cap \Omega$ and $\partial K^b := \partial K \cap \partial \Omega$. The mesh faces are collected in the set $\mathcal{F}$, which is split as $\mathcal{F} = \mathcal{F}^i \cup \mathcal{F}^b$, where $\mathcal{F}^i$ is the collection of the interior faces
(shared by two distinct mesh cells) and $\mathcal{F}^b$ the collection of the boundary faces. We orient every mesh interface $F \in \mathcal{F}^i$ by means of the fixed unit normal vector $n_F$ whose direction is arbitrary but fixed once and for all, whereas we orient every mesh boundary face $F \in \mathcal{F}^b$ by means of the vector $n_F := n_\Omega$. For any mesh cell $K \in \mathcal{T}$, the mesh faces composing its boundary $\partial K$ are collected in the set $\mathcal{F}_{\partial K}$, which is partitioned as $\mathcal{F}_{\partial K} = \mathcal{F}_{\partial K}^i \cup \mathcal{F}_{\partial K}^b$ with obvious notation.

For any real number $\theta \ge 0$, we consider the broken Sobolev space

$$H^\theta(\mathcal{T}) := \{v \in L^2(\Omega) \mid v_K := v_{|K} \in H^\theta(K), \forall K \in \mathcal{T}\}. \quad (4)$$

The jump and average of any function $v \in H^\theta(\mathcal{T})$, $\theta > \frac{1}{2}$, across any mesh interface $F = \partial K_1 \cap \partial K_2 \in \mathcal{F}^i$ are defined by setting $[v]_F(x) := v_{|K_1}(x) - v_{|K_2}(x)$ and $\{v\}_F(x) := \frac{1}{2}(v_{|K_1}(x) + v_{|K_2}(x))$ for a.e. $x \in F$, respectively, where $K_1$ is such that its outward unit normal is $n_F$.

## 2.2 Weak formulations

Recall from Assumption 1.1 that $f \in W^{-1,q}(\Omega)$ with $q \in (\frac{2d}{d+2}, 2]$. Let $q' \in [2, \frac{2d}{d-2})$ be such that $\frac{1}{q} + \frac{1}{q'} = 1$. The weak formulation of the biharmonic problem with type (I) BC's is as follows: Find $u^{(1)} \in H_0^2(\Omega)$ such that

$$(\nabla^2 u^{(1)}, \nabla^2 v)_\Omega = \ell(v) := \langle f, v \rangle_{W^{-1,q}, W_0^{1,q'}}, \quad \forall v \in H_0^2(\Omega). \quad (5)$$

Notice that the right-hand side is meaningful since the Sobolev embedding theorem implies that $H^2(\Omega) \hookrightarrow W^{1,q'}(\Omega)$. Moreover, the model problem (5) is well-posed owing to the Lax-Milgram lemma.

The weak formulation of the biharmonic problem with type (II) BC's is as follows: Find $u^{(II)} \in H_0^1(\Omega) \cap H^2(\Omega)$ such that

$$(\nabla^2 u^{(II)}, \nabla^2 v)_\Omega = \ell(v), \quad \forall v \in H_0^1(\Omega) \cap H^2(\Omega). \quad (6)$$

Owing to the Lax-Milgram lemma, this problem is well-posed (recall that the $H^2$-seminorm defines a norm on $H_0^1(\Omega) \cap H^2(\Omega)$). We also observe that the weak formulations (5) and (6) employ the same bilinear form and the same right-hand side; only the trial and test spaces differ.

## 2.3 Key identity under the classical regularity assumption

Before stating our main result based on our new regularity assumption (Assumption 1.1), it is useful to illustrate the main idea under the classical regularity assumption $f \in H^{-2+s}(\Omega)$ and $u \in H^{2+s}(\Omega)$ with $s \in (\frac{1}{2}, 1]$.

**Lemma 2.1 (Key identity)** Let $s \in (\frac{1}{2}, 1]$. The following holds for all $v \in H^{2+s}(\Omega)$ and all $w \in H^2(\mathcal{T}) \cap H_0^{2-s}(\Omega)$:

$$\langle \Delta^2 v, w \rangle_{H^{-2+s}, H_0^{2-s}} = \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 v, \nabla^2 w_K)_K - (\partial_{nn} v, \partial_n w_K)_{\partial K} \right\}. \quad (7)$$

**Proof.** The proof hinges on a density argument. We consider a mollified sequence of functions $v_\delta \in C^\infty(\bar{\Omega})$ such that $\lim_{\delta \to 0} v_\delta = v$ in $H^{2+s}(\Omega)$. We apply (3) to $v_\delta|_K$ and $w_K := w|_K$ for all $K \in \mathcal{T}$, and sum the result cellwise. This gives

$$(\Delta^2 v_\delta, w)_\Omega = \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 v_\delta, \nabla^2 w_K)_K + (\partial_n \Delta v_\delta, w_K)_{\partial K} - (\partial_{nn} v_\delta, \partial_n w_K)_{\partial K} - (\partial_{nt} v_\delta, \partial_t w_K)_{\partial K} \right\}.$$

The left-hand side is equal to $\langle \Delta^2 v_\delta, w \rangle_{H^{-2+s}, H_0^{2-s}}$. Moreover, the second and fourth summations on the right-hand side vanish. For the second summation, this follows from the fact
that $\partial_n \Delta v_\delta$ is single-valued at the mesh interfaces, whereas $w$ is single-valued at the mesh interfaces and vanishes at the mesh boundary faces. Similar arguments are invoked for the fourth summation. Altogether, we obtain

$$ \langle \Delta^2 v_\delta, w \rangle_{H^{-2+s}, H_0^{2-s}} = \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 v_\delta, \nabla^2 w_K)_K - (\partial_{nn} v_\delta, \partial_n w_K)_{\partial K} \right\}. $$

We can now pass to the limit $\delta \to 0$. This is straightforward for the right-hand side, whereas for the left-hand side, we use [22, Theorem 1.4.4.6] which gives $\lim_{\delta \to 0} \Delta^2 v_\delta \to \Delta^2 v$ in $H^{-2+s}(\Omega)$. This proves the identity (7). □

**Remark 2.2 (Literature)** The identity from Lemma 2.1 is classical; see, e.g., [2]. The present proof is, however, different. The advantages are that it avoids the decomposition of the solution into singular and regular parts and that it works seamlessly in any space dimension.

## 2.4 Key identity under the new regularity assumption

Extending the identity from Lemma 2.1 to the more general setting of Assumption 1.1 requires giving a meaning to the trace of second-order derivatives on each mesh face individually. Before doing this, we briefly recall some material from [19] which was originally devised in the context of second-order elliptic problems.

Let $p$ and $q$ be two real numbers as in Assumption 1.1, i.e., $p > 2$ and $q \in (\frac{2d}{d+2}, 2]$. Let $\varrho \in (2, p]$ be such that $q \ge \frac{\varrho d}{\varrho+d}$ (this is indeed possible since the function $x \mapsto \frac{xd}{x+d}$ is increasing on $[2, \infty)$) and let $\varrho' \in [1, 2)$ be such that $\frac{1}{\varrho} + \frac{1}{\varrho'} = 1$. It is shown in [19] that for every mesh cell $K \in \mathcal{T}$, it is possible to give a meaning separately on every face $F \in \mathcal{F}_{\partial K}$ to the normal component of fields in

$$ \mathbf{V}^d(K) := \{\boldsymbol{\sigma} \in \mathbf{L}^p(K) \mid \nabla \cdot \boldsymbol{\sigma} \in \mathbf{L}^q(K)\}. \qquad (8) $$

Specifically, one defines the operator $\gamma_{K,F}^d : \mathbf{V}^d(K) \to (W_{\frac{1}{\varrho}, \varrho'}(F))'$ such that for all $\varphi \in W_{\frac{1}{\varrho}, \varrho'}(F)$,

$$ \langle \gamma_{K,F}^d(\boldsymbol{\sigma}), \varphi \rangle_F := \int_K \left( \boldsymbol{\sigma} \cdot \nabla L_F^K(\varphi) + (\nabla \cdot \boldsymbol{\sigma}) L_F^K(\varphi) \right) dx, \qquad (9) $$

where $L_F^K : W_{\frac{1}{\varrho}, \varrho'}(F) \to W^{1, \varrho'}(K)$ is a face-to-cell lifting operator satisfying $L_F^K(\varphi)|_F = \varphi$ and $L_F^K(\varphi)|_{\mathcal{F}_{\partial K} \setminus F} = 0$ (see [19, Lemma 3.1]). Notice that $\gamma_{K,F}^d(\boldsymbol{\sigma}) = (\boldsymbol{\sigma} \cdot \mathbf{n}_K)|_F$ whenever the field $\boldsymbol{\sigma}$ is smooth. Consider now the functional space

$$ \mathcal{V}^{(2)}(\Omega) := \{v \in W^{1,p}(\Omega) \mid \Delta v \in L^q(\Omega)\}. \qquad (10) $$

(The superscript refers to the context of second-order PDEs.) We have $(\nabla v)|_K \in \mathbf{V}^d(K)$ for all $v \in \mathcal{V}^{(2)}(\Omega)$ and all $K \in \mathcal{T}$. Moreover, since $\rho > 2$, we have $(w_K)|_F \in W_{\frac{1}{\varrho}, \varrho'}(F)$ for all $w \in H^1(\mathcal{T})$, all $K \in \mathcal{T}$, and all $F \in \mathcal{F}_{\partial K}$. Therefore, it is meaningful to define the following bilinear form for all $v \in \mathcal{V}^{(2)}(\Omega)$ and all $w \in H^1(\mathcal{T})$:

$$ n_{\#}^{(2)}(v, w) := \sum_{K \in \mathcal{T}} \left\{ \sum_{F \in \mathcal{F}_{\partial K}^i} \langle \gamma_{K,F}^d(\nabla v), w_{K|F} - \{\{w\}\}_F \rangle_F + \sum_{F \in \mathcal{F}_{\partial K}^b} \langle \gamma_{K,F}^d(\nabla v), w_{K|F} \rangle_F \right\}. \qquad (11) $$

Since $w_{K|F} - \{\{w\}\}_F = \frac{1}{2}(\mathbf{n}_F \cdot \mathbf{n}_K)[\![w_h]\!]_F$ for all $F \in \mathcal{F}_{\partial K}^i$, we observe that the bilinear form defined in (11) coincides with the one defined in [19, Equ. (3.12)] upon setting therein $\theta_{K,F} = \frac{1}{2}$ if $F \in \mathcal{F}_{\partial K}^i$ and $\theta_{K,F} = 1$ otherwise. Hence, invoking [19, Lemma 3.3], we infer that for all $v \in \mathcal{V}^{(2)}(\Omega)$ and all $w \in H^1(\mathcal{T})$, the following important relation holds:

$$ -(\Delta v, w)_{\Omega} = \sum_{K \in \mathcal{T}} (\nabla v, \nabla w_K)_{K} - n_{\#}^{(2)}(v, w). \qquad (12) $$
We are now ready to address the biharmonic problem. We define the functional space

$$ \mathcal{V}^{(4)}(\Omega) := \{v \in W^{2,p}(\Omega) \mid \Delta v \in W^{1,q}(\Omega)\}. \quad (13) $$

(The superscript refers to the context of fourth-order PDEs.) We also define the following bilinear form for all $v \in \mathcal{V}^{(4)}(\Omega)$ and all $w \in H^2(\mathcal{T})$:

$$ n_{\#}^{(4)}(v, w) := \sum_{i \in \{1:d\}} n_{\#}^{(2)}(\partial_i v, \partial_i w). \quad (14) $$

This definition is meaningful since for all $v \in \mathcal{V}^{(4)}(\Omega)$ and all $w \in H^2(\mathcal{T})$, we have $\partial_i v \in \mathcal{V}^{(2)}(\Omega)$ and $\partial_i w \in H^1(\mathcal{T})$, for all $i \in \{1:d\}$.

**Lemma 2.3 (Key identity)** The following holds for all $v \in \mathcal{V}^{(4)}(\Omega)$ and all $w \in H^2(\mathcal{T}) \cap H_0^1(\Omega)$:

$$ \langle \Delta^2 v, w \rangle_{W^{-1,q}, W_0^{1,q'}} = \sum_{K \in \mathcal{T}} (\nabla^2 v, \nabla^2 w_K)_K - n_{\#}^{(4)}(v, w). \quad (15) $$

**Proof.** Since $\Delta v \in W^{1,q}(\Omega)$, we have $\Delta^2 v = \Delta(\Delta v) \in W^{-1,q}(\Omega)$. Moreover, since $w \in H_0^1(\Omega)$, we infer that

$$ \langle \Delta^2 v, w \rangle_{W^{-1,q}, W_0^{1,q'}} = -(\nabla \Delta v, \nabla w)_\Omega = \sum_{i \in \{1:d\}} -(\Delta \partial_i v, \partial_i w)_\Omega. $$

For all $i \in \{1:d\}$, we have $\partial_i v \in \mathcal{V}^{(2)}(\Omega)$ and $\partial_i w \in H^1(\mathcal{T})$. Applying (12) to the right-hand side of the above equation, we obtain

$$ \langle \Delta^2 v, w \rangle_{W^{-1,q}, W_0^{1,q'}} = \sum_{i \in \{1:d\}} \sum_{K \in \mathcal{T}} (\nabla \partial_i v, \nabla \partial_i w_K)_K - \sum_{i \in \{1:d\}} n_{\#}^{(2)}(\partial_i v, \partial_i w). $$

The conclusion is straightforward by definition of the Hessian and of $n_{\#}^{(4)}$.

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

The following reformulation of $n_{\#}^{(4)}$ will be useful in our analysis. Let us set $n_{F,i} := n_F \cdot e_i$ for all $i \in \{1:d\}$, where $(e_i)_{i \in \{1:d\}}$ denotes the canonical Cartesian basis of $\mathbb{R}^d$.

**Lemma 2.4 (Reformulation of $n_{\#}^{(4)}$)** For all $v \in \mathcal{V}^{(4)}(\Omega)$ and all $w \in H^2(\mathcal{T}) \cap H_0^1(\Omega)$, we have

$$ n_{\#}^{(4)}(v, w) = \sum_{K \in \mathcal{T}} \sum_{i \in \{1:d\}} \left\{ \sum_{F \in \mathcal{F}_{\partial K}^i} \langle \gamma_{K,F}^d(\nabla \partial_i v), n_{F,i} n_F \cdot ((\nabla w_K)|_F - \{\{\nabla w\}\}_F) \rangle_F \right. \\ \left. + \sum_{F \in \mathcal{F}_{\partial K}^b} \langle \gamma_{K,F}^d(\nabla \partial_i v), n_{F,i} n_F \cdot (\nabla w_K)|_F \rangle_F \right\}. \quad (16) $$

**Proof.** We first use the definitions of $n_{\#}^{(2)}$ and $n_{\#}^{(4)}$ to write

$$ n_{\#}^{(4)}(v, w) = \sum_{K \in \mathcal{T}} \sum_{i \in \{1:d\}} \left\{ \sum_{F \in \mathcal{F}_{\partial K}^i} \langle \gamma_{K,F}^d(\nabla \partial_i v), (\partial_i w_K)|_F - \{\{\partial_i w\}\}_F \rangle_F + \sum_{F \in \mathcal{F}_{\partial K}^b} \langle \gamma_{K,F}^d(\nabla \partial_i v), (\partial_i w_K)|_F \rangle_F \right\}. $$

We observe that for all $K \in \mathcal{T}$ and all $F \in \mathcal{F}_{\partial K}$,

$$ (\partial_i w_K)|_F = e_i \cdot (\nabla w_K)|_F = n_{F,i} n_F \cdot (\nabla w_K)|_F + e_i \cdot (\nabla_t w_K)|_F, $$

where $(\nabla_t w_K)|_F$ denotes the tangential gradient of $w_K$ on $F$. Since $w \in H_0^1(\Omega)$, this quantity is single-valued on every mesh interface $F \in \mathcal{F}^i$ and vanishes on every mesh boundary face $F \in \mathcal{F}^b$. Moreover, since $\nabla \partial_i v$ has an integrable divergence on $\Omega$ by assumption, we infer that for all $F = \partial K_1 \cap \partial K_2 \in \mathcal{F}^i$,

$$ \sum_{j \in \{1,2\}} \langle \gamma_{K_j,F}^d(\nabla \partial_i v), \phi \rangle_F = 0, \quad \forall \phi \in W_{\theta}^{1,q'}(F). \quad (17) $$
This implies that

$$ \sum_{K \in \mathcal{T}} \sum_{F \in \mathcal{F}_{\partial K}} \langle \gamma_{K,F}^{\mathrm{d}}(\nabla \partial_i v), (\partial_i w_K)|_F \rangle_F = \sum_{K \in \mathcal{T}} \sum_{F \in \mathcal{F}_{\partial K}} \langle \gamma_{K,F}^{\mathrm{d}}(\nabla \partial_i v), n_{F,i} n_F \cdot (\nabla w_K)|_F \rangle_F. $$

A similar reasoning shows that

$$ \sum_{K \in \mathcal{T}} \sum_{F \in \mathcal{F}_{\partial K}} \langle \gamma_{K,F}^{\mathrm{d}}(\nabla \partial_i v), \{\partial_i w\}_F \rangle_F = \sum_{K \in \mathcal{T}} \sum_{F \in \mathcal{F}_{\partial K}} \langle \gamma_{K,F}^{\mathrm{d}}(\nabla \partial_i v), n_{F,i} n_F \cdot \{\nabla w\}_F \rangle_F. \quad (18) $$

Combining the two above identities proves the assertion.

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

**Remark 2.5 (Simplifications)** Both terms in (18) actually vanish owing to (17). This means that for all $v \in \mathcal{V}^{(4)}(\Omega)$ and all $w \in H^2(\mathcal{T})$, we have

$$ n_{\#}^{(4)}(v, w) = \sum_{K \in \mathcal{T}} \sum_{F \in \mathcal{F}_{\partial K}} \sum_{i \in \{1:d\}} \langle \gamma_{K,F}^{\mathrm{d}}(\nabla \partial_i v), (\partial_i w_K)|_F \rangle_F, \quad (19a) $$

and whenever $w \in H^2(\mathcal{T}) \cap H_0^1(\Omega)$, we also have

$$ n_{\#}^{(4)}(v, w) = \sum_{K \in \mathcal{T}} \sum_{F \in \mathcal{F}_{\partial K}} \sum_{i \in \{1:d\}} \langle \gamma_{K,F}^{\mathrm{d}}(\nabla \partial_i v), n_{F,i} n_F \cdot (\nabla w_K)|_F \rangle_F. \quad (19b) $$

However, the expressions (14) and (16) are those that are needed in this work because we will extend the domain of the first argument $v$ to spaces containing functions which are only piecewise smooth and for which (17) no longer holds true. In this situation, it is important to use (14) and (16) and not (19).

# 3 Discrete setting for $C^0$-HHO methods

In this section, we introduce the key ingredients to formulate $C^0$-HHO methods to approximate the biharmonic problem with both types of BC's. The starting point for the present $C^0$-HHO methods are the fully discontinuous HHO methods from [13]. While the latter methods rely on a triple of discrete unknowns, approximating the solution in each mesh cell, its trace on each mesh face, and the trace of its normal derivative (oriented by $n_F$) on each mesh face, the present $C^0$-HHO methods only rely on a pair of discrete unknowns, approximating the solution in each mesh cell and the trace of its normal derivative on each mesh face.

## 3.1 Local reconstruction, stabilization, and stability

Let $k \ge 0$ be the polynomial degree. Recall that we consider a mesh $\mathcal{T}$ from a shape-regular family of simplicial meshes such that $\mathcal{T}$ covers $\Omega$ exactly. For every mesh cell $K \in \mathcal{T}$, the local $C^0$-HHO space is

$$ \hat{V}_K^k := \mathbb{P}_{k+2}(K) \times \mathbb{P}_k(\mathcal{F}_{\partial K}), \quad (20) $$

with the broken polynomial space $\mathbb{P}_k(\mathcal{F}_{\partial K}) := \times_{F \in \mathcal{F}_{\partial K}} \mathbb{P}_k(F)$. A generic element in $\hat{V}_K^k$ is denoted $\hat{v}_K = (v_K, \gamma_{\partial K})$ with $v_K \in \mathbb{P}_{k+2}(K)$ and $\gamma_{\partial K} \in \mathbb{P}_k(\mathcal{F}_{\partial K})$. The first component of $\hat{v}_K$ aims at representing the solution inside the mesh cell and the second the trace of its normal derivative (oriented by $n_K$) on the cell boundary.

Let $K \in \mathcal{T}$. We define the local reconstruction operator $R_K: \hat{V}_K^k \to \mathbb{P}_{k+2}(K)$ such that, for all $\hat{v}_K \in \hat{V}_K^k$ with $\hat{v}_K := (v_K, \gamma_{\partial K})$, the polynomial $R_K(\hat{v}_K) \in \mathbb{P}_{k+2}(K)$ is uniquely defined by solving the following problem with test functions $w \in \mathbb{P}_{k+2}(K)^\perp := \{w \in \mathbb{P}_{k+2}(K) | (w, \xi)_K = 0, \forall \xi \in \mathbb{P}_1(K)\}$:

$$ (\nabla^2 R_K(\hat{v}_K), \nabla^2 w)_K = (\nabla^2 v_K, \nabla^2 w)_K - (\partial_n v_K - \gamma_{\partial K}, \partial_{nn} w)|_{\partial K}, \quad (21) $$
together with the condition $(R_K(\hat{v}_K), \xi)_K = (v_K, \xi)_K$ for all $\xi \in \mathbb{P}_1(K)$. The local stabilization bilinear form is defined such that, for all $(\hat{v}_K, \hat{w}_K) \in \hat{V}_K^k \times \hat{V}_K^k$, with $\hat{v}_K := (v_K, \gamma_{\partial K})$ and $\hat{w}_K := (w_K, \chi_{\partial K})$,

$$S_{\partial K}(\hat{v}_K, \hat{w}_K) := h_K^{-1} (\Pi_{\partial K}^k(\gamma_{\partial K} - \partial_n v_K), \chi_{\partial K} - \partial_n w_K)_{\partial K}, \quad (22)$$

where $\Pi_{\partial K}^k$ denotes the $L^2$-orthogonal projection onto $\mathbb{P}_k(\mathcal{F}_{\partial K})$.

We define the local bilinear form $a_K$ on $\hat{V}_K^k \times \hat{V}_K^k$ such that

$$a_K(\hat{v}_K, \hat{w}_K) := (\nabla^2 R_K(\hat{v}_K), \nabla^2 R_K(\hat{w}_K))_K + S_{\partial K}(\hat{v}_K, \hat{w}_K). \quad (23)$$

We define the local energy seminorm such that, for all $\hat{v}_K := (v_K, \gamma_{\partial K}) \in \hat{V}_K^k$,

$$|\hat{v}_K|_{\hat{V}_K^k}^2 := \|\nabla^2 v_K\|_{K}^2 + h_K^{-1} \|\gamma_{\partial K} - \partial_n v_K\|_{\partial K}^2. \quad (24)$$

**Lemma 3.1 (Local stability and boundedness)** There is a real number $\alpha > 0$, depending only on the mesh shape-regularity and the polynomial degree $k$, such that, for all $h > 0$, all $K \in \mathcal{T}$, and all $\hat{v}_K \in \hat{V}_K^k$,

$$\alpha |\hat{v}_K|_{\hat{V}_K^k}^2 \le \|\nabla^2 R_K(\hat{v}_K)\|_{K}^2 + S_{\partial K}(\hat{v}_K, \hat{v}_K) \le \alpha^{-1} |\hat{v}_K|_{\hat{V}_K^k}^2. \quad (25)$$

**Proof.** The proof proceeds as that of [13, Lemma 4.1]. □

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

**Remark 3.2 (Cell unknowns)** The choice of the polynomial space $\mathbb{P}_{k+2}(K)$ leverages on [13] and has been made to allow for the simple stabilization bilinear form defined in (22). Alternative choices for the cell unknowns are possible, such as taking the polynomial space $\mathbb{P}_{k+1}(K)$ (see [13, Remark 3.2]) and even $\mathbb{P}_k(K)$ for $k \ge 1$ (see [1]). Both choices lead to fewer cell unknowns, which can marginally alleviate the costs of static condensation, but require a more subtle form for stabilization involving the reconstruction operator, and this typically adds some computational costs. Moreover, whatever the choice for the cell unknowns, the size of the linear system after static condensation is the same.

## 3.2 Global discrete spaces and discrete problems

Let $\mathbb{P}_{k+2}^g(\mathcal{T})$ denote the usual $C^0$-conforming finite element space composed of continuous, piecewise polynomials of degree at most $(k+2)$ on the mesh $\mathcal{T}$ (the superscript g refers to the integrability of the gradient of functions in $\mathbb{P}_{k+2}^g(\mathcal{T})$). We set

$$\hat{V}_h^k := \mathbb{P}_{k+2}^g(\mathcal{T}) \times V_{\mathcal{F}}^k, \quad V_{\mathcal{F}}^k := \mathbb{P}_k(\mathcal{F}), \quad (26)$$

with the broken polynomial space $\mathbb{P}_k(\mathcal{F}) := \times_{F \in \mathcal{F}} \mathbb{P}_k(F)$. We also define the following subspaces accounting for (homogeneous) boundary conditions:

$$V_{\mathcal{T},0}^{k+2} := \mathbb{P}_{k+2}^g(\mathcal{T}) \cap H_0^1(\Omega), \quad V_{\mathcal{F},0}^k := \{\gamma_{\mathcal{F}} \in V_{\mathcal{F}}^k \mid \gamma_{\mathcal{F}|F} = 0, \forall F \in \mathcal{F}^b\}. \quad (27)$$

The global $C^0$-HHO spaces used to approximate the biharmonic problems are then defined as follows:

$$\hat{V}_h^{k,(I)} := V_{\mathcal{T},0}^{k+2} \times V_{\mathcal{F},0}^k, \quad \hat{V}_h^{k,(II)} := V_{\mathcal{T},0}^{k+2} \times V_{\mathcal{F}}^k, \quad (28)$$

for type (I) and type (II) BC's, respectively. This means that the boundary condition $u = 0$ is directly enforced on the trace of the cell unknowns at the boundary, the boundary condition $\partial_n u = 0$ for type (I) BC's is directly enforced on the discrete unknowns representing the normal derivative at the mesh boundary faces, and the boundary condition $\partial_{nn} u = 0$ is not enforced directly, but in a natural way (i.e., it results from the discrete problem). Notice also that $\hat{V}_h^{k,(I)} \subset \hat{V}_h^{k,(II)}$.

Any member of $\hat{V}_h^{k,(I)}$ or $\hat{V}_h^{k,(II)}$ is generically denoted by $\hat{v}_h := (v_{\mathcal{T}}, \gamma_{\mathcal{F}})$ with $v_{\mathcal{T}} := (v_K)_{K \in \mathcal{T}}$ and $\gamma_{\mathcal{F}} := (\gamma_F)_{F \in \mathcal{F}}$. For every mesh cell $K \in \mathcal{T}$, the local components of $\hat{v}_h$ are
collected in the pair $\hat{v}_K := (v_K, \gamma_{\partial K}) \in \hat{V}_K^k$ with $\gamma_{\partial K}|_F := (n_F \cdot n_K)\gamma_F$ for all $F \in \mathcal{F}_{\partial K}$.
The global bilinear form is assembled cellwise from the contributions of all the mesh cells,
yielding

$$a_h(\hat{v}_h, \hat{w}_h) := \sum_{K \in \mathcal{T}} a_K(\hat{v}_K, \hat{w}_K). \qquad (29)$$

The discrete problems are as follows:

* For type (I) BC's, one seeks $\hat{u}_h^{(I)} \in \hat{V}_h^{k,(I)}$ such that

$$a_h(\hat{u}_h^{(I)}, \hat{w}_h) = \ell(w_{\mathcal{T}}), \quad \forall w_h \in \hat{V}_h^{k,(I)}. \qquad (30)$$

* For type (II) BC's, one seeks $\hat{u}_h^{(II)} \in \hat{V}_h^{k,(II)}$ such that

$$a_h(\hat{u}_h^{(II)}, \hat{w}_h) = \ell(w_{\mathcal{T}}), \quad \forall w_h \in \hat{V}_h^{k,(II)}. \qquad (31)$$

Notice that the discrete problems (30) and (31) employ the same discrete bilinear form $a_h$ and the same right-hand side; only the discrete trial and test spaces differ. We also observe that the above right-hand sides are meaningful since $V_{\mathcal{T},0}^{k+2} \subset W_0^{1,q'}(\Omega)$.

Both spaces $\hat{V}_h^{k,(I)}$ and $\hat{V}_h^{k,(II)}$ are equipped with the norm

$$\|\hat{v}_h\|_{\hat{V}_h^k}^2 := \sum_{K \in \mathcal{T}} |\hat{v}_K|^2_{\hat{V}_K^k}, \qquad (32)$$

with the local energy seminorm $|\cdot|_{\hat{V}_K^k}$ defined in (24). To verify that this indeed defines a norm on $\hat{V}_h^{k,(II)}$ (and thus also on $\hat{V}_h^{k,(I)}$), we notice that if $\hat{v}_h = (v_{\mathcal{T}}, \gamma_{\mathcal{F}}) \in \hat{V}_h^{k,(II)}$ satisfies $\|\hat{v}_h\|_{\hat{V}_h^k} = 0$, then $v_{\mathcal{T}}$ is a globally affine function in $\Omega$ which vanishes at the boundary $\partial\Omega$, so that $v_{\mathcal{T}} = 0$; moreover, $\gamma_{\mathcal{F}}$ coincides on each mesh face with the trace of the normal derivative of $v_{\mathcal{T}}$, so that $v_{\mathcal{F}} = 0$ as well. A direct consequence of the lower bound in (25) is therefore that the global bilinear form $a_h$ is coercive on $\hat{V}_h^{k,(II)}$ (and thus also on $\hat{V}_h^{k,(I)}$). Hence, both discrete problems (30) and (31) are well-posed owing to the Lax-Milgram lemma.

**Remark 3.3 (Static condensation)** Unlike the fully discontinuous HHO methods from [1, 13, 12] where all the cell unknowns can be eliminated locally by a static condensation procedure, $C^0$-HHO methods are amenable to static condensation only if $(k+2) \ge (d+1)$ on $d$-dimensional simplices, and in this case, only the so-called cell bubble functions vanishing on the cell boundary can be eliminated, as in classical $C^0$-conforming FEM. Notice also that static condensation is more delicate for $C^0$-IPDG methods since only those bubble functions having a zero second-order derivative at the cell boundary can be locally eliminated.

## 3.3 Analysis tools and local HHO interpolation operator

For the reader's convenience, let us briefly restate some classical discrete inverse inequalities and polynomial approximation properties on shape-regular families of simplicial meshes. The results are classical, and we refer the reader, e.g., to [17, Chap. 9-11] for the proofs. In the rest of this paper, we use the symbol $C$ to denote any positive generic constant (its value can change at each occurrence) that is independent of $h > 0$, the considered mesh cell $K \in \mathcal{T}$, and the considered function in the inequality. The value of $C$ can depend on the shape-regularity parameter of the mesh sequence and the polynomial degree (whenever relevant).

**Lemma 3.4 (Discrete inverse inequalities)** Let $l \ge 0$ be a polynomial degree. There is $C$ such that for all $h > 0$, all $K \in \mathcal{T}$, all $p \in [1, \infty]$, and all $v_h \in \mathbb{P}_l(K)$,

$$\|v_h\|_{\partial K} \le Ch_K^{-\frac{1}{2}} \|v_h\|_K, \qquad (33)$$

$$\|\nabla v_h\|_K \le Ch_K^{-1} \|v_h\|_K, \qquad (34)$$

$$\|v_h\|_{L^p(K)} \le Ch_K^{d(\frac{1}{p}-\frac{1}{2})} \|v_h\|_K. \qquad (35)$$
**Lemma 3.5 (Multiplicative trace inequality)** There is C such that for all $h > 0$, all $K \in \mathcal{T}$, and all $v \in H^1(K)$,

$$ \|v\|_{\partial K} \le C(h_K^{-\frac{1}{2}}\|v\|_K + h_K^{\frac{1}{2}}\|\nabla v\|_K). \quad (36) $$

**Lemma 3.6 (Polynomial approximation)** Let $l \ge 1$ be the polynomial degree and let $d \in \{2, 3\}$ (so that $l+1 > \frac{d}{2}$). Let $\mathcal{L}_h^l$ be the $H^1$-conforming Lagrange interpolation operator onto $\mathbb{P}_l^g(\mathcal{T})$ and let $\mathcal{L}_K^l$ be its local version mapping onto $\mathbb{P}_{k+2}(K)$ for all $K \in \mathcal{T}$. There is $C$ such that for all $r > \frac{d}{2}$, all $m \in \{0, \dots, \lfloor r \rfloor\}$, all $p \in [1, \infty]$, all $h > 0$, all $K \in \mathcal{T}$, and all $v \in H^r(K)$,

$$ h_K^m \|v - \mathcal{L}_K^l(v)\|_{H^m(K)} + h_K^{d(\frac{1}{2} - \frac{1}{p})} \|v - \mathcal{L}_K^l(v)\|_{L^p(K)} \le C h_K^t |v|_{H^t(K)}, \quad (37) $$

where $t := \min\{r, l+1\}$.

Let us briefly highlight some useful consequences of Lemma 3.6. First, taking $l = 1$, $r = 2$, and $m \in \{0, 1\}$ in (37) shows that for all $h > 0$, all $K \in \mathcal{T}$, and all $v \in H^2(K)$,

$$ \|v - \mathcal{L}_K^1(v)\|_K + h_K \|\nabla(v - \mathcal{L}_K^1(v))\|_K \le C h_K^2 \|\nabla^2 v\|_K. \quad (38) $$

Then, for all $v \in H^2(K)$, setting $z := v - \mathcal{L}_K^{k+2}(v)$ for any polynomial degree $k \ge 0$, and observing that $\mathcal{L}_K^1(z) = 0$ since $\mathcal{L}_K^1 \circ \mathcal{L}_K^{k+2} = \mathcal{L}_K^1$, we infer from (38) that $\|\nabla z\|_K \le Ch_K \|\nabla^2 z\|_K$. Hence, owing to the multiplicative trace inequality (36), we have $h_K^{-\frac{1}{2}} \|\nabla z\|_{\partial K} \le C \|\nabla^2 z\|_K$, i.e.,

$$ h_K^{-\frac{1}{2}} \|\nabla(v - \mathcal{L}_K^{k+2}(v))\|_{\partial K} \le C \|\nabla^2(v - \mathcal{L}_K^{k+2}(v))\|_K. \quad (39) $$

**Remark 3.7 (H$^1$-conforming quasi-interpolation operator)** The lower bound $r > \frac{d}{2}$ is not really a restriction for the biharmonic problem in space dimension $d \in \{2, 3\}$ since the weak solution always sits in $H^2(\Omega)$. In higher space dimension, the restriction $r > \frac{d}{2}$ can be lifted by invoking a quasi-interpolation operator onto $\mathbb{P}_l^g(\mathcal{T})$ instead of the Lagrange interpolation operator; we refer the reader to [27, 16] for examples of quasi-interpolation operators that can be considered.

For every mesh cell $K \in \mathcal{T}$, we define the local HHO reduction operator $\hat{\mathcal{I}}_K^k : H^2(K) \to \hat{V}_K^k$ such that, for all $v \in H^2(K)$,

$$ \hat{\mathcal{I}}_K^k(v) := (\mathcal{L}_K^{k+2}(v), \Pi_{\partial K}^k(n_K \cdot \nabla v)) \in \hat{V}_K^k. \quad (40) $$

This definition is meaningful since $H^2(K) \hookrightarrow C^0(\bar{K})$ for $d \in \{2, 3\}$. In addition, we define the local HHO interpolation operator

$$ J_K^{\text{HHO}} := R_K \circ \hat{\mathcal{I}}_K^k : H^2(K) \to \mathbb{P}_{k+2}(K). \quad (41) $$

In what follows, we need to measure the HHO interpolation error in some augmented norm which reflects Assumption 1.1. Specifically, we set

$$ \|v\|_{\sharp,K}^2 := \|\nabla^2 v\|_K^2 + h_K^{2d(\frac{1}{2} - \frac{1}{p})} \|\nabla^2 v\|_{L^p(K)}^2 + h_K^{2+2d(\frac{1}{2} - \frac{1}{q})} \|\nabla \Delta v\|_{L^q(K)}^2, \quad (42) $$

for all $v \in W^{2,p}(K)$ with $\Delta v \in W^{1,q}(K)$ (recall that $p > 2$ and $q \in (\frac{2d}{d+2}, 2]$).

**Lemma 3.8 (HHO interpolation error)** In the above setting, the following holds for all $K \in \mathcal{T}$:

$$ \|v - J_K^{\text{HHO}}(v)\|_{\sharp,K}^2 + S_{\partial K}(\hat{\mathcal{I}}_K^k(v), \hat{\mathcal{I}}_K^k(v)) \le C \|v - \mathcal{L}_K^{k+2}(v)\|_{\sharp,K}^2. \quad (43) $$
**Proof.** (1) Let us first bound $\|v - J_K^{\text{HMO}}(v)\|_{\sharp,K}$. The triangle inequality followed by the discrete inverse inequalities from Lemma 3.4 implies that

$$ \begin{aligned} \|v - J_K^{\text{HMO}}(v)\|_{\sharp,K} &\le \|v - \mathcal{L}_K^{k+2}(v)\|_{\sharp,K} + \|J_K^{\text{HMO}}(v) - \mathcal{L}_K^{k+2}(v)\|_{\sharp,K} \\ &\le \|v - \mathcal{L}_K^{k+2}(v)\|_{\sharp,K} + C \|\nabla^2(J_K^{\text{HMO}}(v) - \mathcal{L}_K^{k+2}(v))\|_{K}, \end{aligned} $$

so that we only need to bound the last term on the right-hand side. A straightforward calculation using the definition (21) of $R_K$ shows that, for all $w \in \mathbb{P}_{k+2}(K)^{\perp}$,

$$ \begin{aligned} (\nabla^2 J_K^{\text{HMO}}(v), \nabla^2 w)_K &= (\nabla^2 \mathcal{L}_K^{k+2}(v), \nabla^2 w)_K - (\partial_n \mathcal{L}_K^{k+2}(v) - \Pi_{\partial K}^k(\partial_n v), \partial_{nn} w)_{\partial K} \\ &= (\nabla^2 \mathcal{L}_K^{k+2}(v), \nabla^2 w)_K - (\partial_n(\mathcal{L}_K^{k+2}(v) - v), \partial_{nn} w)_{\partial K}, \end{aligned} $$

together with the relation $(J_K^{\text{HMO}}(v), \xi)_K = (\mathcal{L}_K^{k+2}(v), \xi)_K$, for all $\xi \in \mathbb{P}_1(K)$. Re-arranging the terms, taking $w = J_K^{\text{HMO}}(v) - \mathcal{L}_K^{k+2}(v)$, and using the discrete trace inverse inequality (33) gives

$$ \|\nabla^2(J_K^{\text{HMO}}(v) - \mathcal{L}_K^{k+2}(v))\|_K \le C h_K^{-\frac{1}{2}} \|\partial_n(\mathcal{L}_K^{k+2}(v) - v)\|_{\partial K}. $$

Bounding the right-hand side using (39) leads to

$$ \|\nabla^2(J_K^{\text{HMO}}(v) - \mathcal{L}_K^{k+2}(v))\|_K \le C \|\nabla^2(v - \mathcal{L}_K^{k+2}(v))\|_K. $$

Putting the above bounds together shows that $\|v - J_K^{\text{HMO}}(v)\|_{\sharp,K} \le C \|v - \mathcal{L}_K^{k+2}(v)\|_{\sharp,K}$ since the $\|\cdot\|_{\sharp,K}$-norm controls the $H^2(K)$-seminorm.

(2) Let us now bound $S_{\partial K}(\hat{\mathcal{I}}_K^k(v), \hat{\mathcal{I}}_K^k(v))$. We have

$$ S_{\partial K}(\hat{\mathcal{I}}_K^k(v), \hat{\mathcal{I}}_K^k(v)) = h_K^{-1} \|\Pi_{\partial K}^k(\Pi_{\partial K}^k(\partial_n v) - \partial_n(\mathcal{L}_K^{k+2}(v)))\|_{\partial K}^2 \le h_K^{-1} \|\partial_n(v - \mathcal{L}_K^{k+2}(v))\|_{\partial K}^2, $$

since $\Pi_{\partial K}^k \circ \Pi_{\partial K}^k = \Pi_{\partial K}^k$ and $\Pi_{\partial K}^k$ is $L^2$-stable. The right-hand side is bounded by means of (39), yielding

$$ S_{\partial K}(\hat{\mathcal{I}}_K^k(v), \hat{\mathcal{I}}_K^k(v)) \le C \|\nabla^2(v - \mathcal{L}_K^{k+2}(v))\|_K^2. $$

Since the $\|\cdot\|_{\sharp,K}$-norm controls the $H^2(K)$-seminorm, this concludes the proof. $\square$

It is convenient to define global versions of the above operators and norms. The global HHO reduction operator $\hat{\mathcal{I}}_h^k : H^2(\Omega) \to \hat{V}_h^k$ is defined such that, for all $v \in H^2(\Omega)$,

$$ \hat{\mathcal{I}}_h^k(v) := (\mathcal{L}_h^{k+2}(v), (\Pi_F^k(n_F \cdot (\nabla v)|_F))_{F \in \mathcal{F}}) \in \hat{V}_h^k, \quad (44) $$

so that the local components of $\hat{\mathcal{I}}_h^k(v)$ are $\hat{\mathcal{I}}_K^k(v|_K)$ for all $K \in \mathcal{T}$. We also notice that $\hat{\mathcal{I}}_h^k(v) \in \hat{V}_h^{k,(II)}$ if $v \in H_0^1(\Omega)$ since the Lagrange interpolation operator preserves the homogeneous Dirichlet boundary condition. Moreover, we have $\hat{\mathcal{I}}_h^k(v) \in \hat{V}_h^{k,(I)}$ if $v \in H_0^2(\Omega)$. Furthermore, the global HHO interpolation operator mapping onto $\mathbb{P}_{k+2}(\mathcal{T})$ is defined such that $J_h^{\text{HMO}}(v)|_K := J_K^{\text{HMO}}(v|_K)$ for all $K \in \mathcal{T}$ and all $v \in H^2(\Omega)$. Finally, we define the global augmented norm such that $\|v\|_{\sharp,h}^2 := \sum_{K \in \mathcal{T}} \|v|_K\|_{\sharp,K}^2$ for all $v \in \mathcal{V}^{(4)}(\Omega)$.

# 4 Error analysis

In this section, we perform the error analysis under Assumption 1.1. For completeness, we also outline the main arguments for the error analysis of $C^0$-IPDG methods. In all cases, the main step is to bound the consistency error by making use of the key identity established in Lemma 2.3.
## 4.1 Preliminaries

Before bounding the consistency error, we need to slightly adapt the bilinear form $n_{\#}^{(4)}$ introduced in Section 2.4. Recall the functional space $V^{(4)}(\Omega) := \{v \in W^{2,p}(\Omega) | \Delta v \in W^{1,q}(\Omega)\}$, $p > 2$, $q \in (\frac{2d}{d+2}, 2]$ (see (13)) and that $n_{\#}^{(4)}$ is defined on $V^{(4)}(\Omega) \times H^2(\mathcal{T})$ using (14) or equivalently (16).

The first adaptation is to allow for discrete functions as the first argument of $n_{\#}^{(4)}$. Since functions in $\mathbb{P}_{k+2}^{g}(\mathcal{T})$ are piecewise smooth, $n_{\#}^{(4)}$ can be extended to $(V^{(4)}(\Omega) + \mathbb{P}_{k+2}^{g}(\mathcal{T})) \times H^2(\mathcal{T})$. Moreover, since $\gamma_{K,F}^d(\nabla\partial_i v_h) = (\partial_n\partial_i v_K)|_F$ for all $v_h \in \mathbb{P}_{k+2}^{g}(\mathcal{T})$ (with the notation $v_K := v_{h|_K}$), a straightforward calculation starting from (16) shows that for all $v_h, w_h \in \mathbb{P}_{k+2}^{g}(\mathcal{T})$,

$$n_{\#}^{(4)}(v_h, w_h) = \sum_{F \in \mathcal{F}} (\{\partial_{nn} v_h\}_F, \llbracket \nabla w_h \rrbracket_F \cdot n_F)_F, \quad (45)$$

where the second-normal derivative is understood to act cellwise on $v_h$. This identity is important when analyzing the $C^0$-IPDG method (see Section 4.4).

However, when analyzing the $C^0$-HHO method, a second adaptation is necessary since discrete test functions in the $C^0$-HHO method have two components, one attached to the mesh cells and one to the mesh faces, and not just one as in the $C^0$-IPDG method. Thus, inspired from (16), we now introduce the following bilinear form on $(V^{(4)}(\Omega) + \mathbb{P}_{k+2}^{g}(\mathcal{T})) \times \hat{V}_h^{k,(\mathrm{II})}$ (recall that $\hat{V}_h^{k,(\mathrm{II})} \subset \hat{V}_h^{k,(\mathrm{II})}$):

$$\hat{n}_{\#}^{(4)}(v, \hat{w}_h) := \sum_{K \in \mathcal{T}} \sum_{F \in \mathcal{F}_{\partial K}} \sum_{i \in \{1:d\}} \langle \gamma_{K,F}^{d}(\nabla \partial_i v), n_{F,i}(n_F \cdot (\nabla w_K)|_F - \chi_F) \rangle_F. \quad (46)$$

Recall the notation $\hat{w}_h = (w_{\mathcal{T}}, \chi_{\mathcal{F}})$ with $w_{\mathcal{T}} = (w_K)_{K \in \mathcal{T}}$ and $\chi_{\mathcal{F}} = (\chi_F)_{F \in \mathcal{F}}$ and that the local components of $\hat{w}_h$ associated with the mesh cell $K \in \mathcal{T}$ are $\hat{w}_K = (w_K, \chi_{\partial K} = ((n_K \cdot n_F)\chi_F)_{F \in \mathcal{F}_{\partial K}})$.

**Lemma 4.1 (Identities for $\hat{n}_{\#}^{(4)}$)** The following holds for all $v_h \in \mathbb{P}_{k+2}^{g}(\mathcal{T})$, all $v \in V^{(4)}(\Omega)$, and all $\hat{w}_h \in \hat{V}_h^{k,(\mathrm{II})}$:

$$\hat{n}_{\#}^{(4)}(v_h, \hat{w}_h) = \sum_{K \in \mathcal{T}} (\nabla^2 v_K, \nabla^2 (w_K - R_K(\hat{w}_K)))_K, \quad (47)$$

$$\hat{n}_{\#}^{(4)}(v, \hat{w}_h) = n_{\#}^{(4)}(v, w_{\mathcal{T}}) - \sum_{F \in \mathcal{F}^b} \sum_{i \in \{1:d\}} \langle \gamma_{K,F}^{d}(\nabla \partial_i v), n_{F,i} \chi_F \rangle_F, \quad (48)$$

where in the summation over all $F \in \mathcal{F}^b$, $K \in \mathcal{T}$ denotes the unique mesh cell such that $F = \partial K \cap \partial \Omega$.

**Proof.** (1) Proof of (47). Since $v_h$ is piecewise smooth, we have $\gamma_{K,F}^d(\nabla\partial_i v_K)|_F = (n_K \cdot \nabla\partial_i v_K)|_F$ for all $K \in \mathcal{T}$ and all $F \in \mathcal{F}_{\partial K}$, so that

$$\sum_{i \in \{1:d\}} \langle \gamma_{K,F}^{d}(\nabla \partial_i v_K)|_F, n_{F,i}(n_F \cdot (\nabla w_K)|_F - \chi_F) \rangle_F = \int_F \partial_{nn} v_K (\partial_n w_K - \chi_{\partial K}) \, ds,$$

where all the normal derivatives are understood to be oriented by $n_K$. Therefore, we obtain

$$\hat{n}_{\#}^{(4)}(v_h, \hat{w}_h) = \sum_{K \in \mathcal{T}} (\partial_{nn} v_K, \partial_n w_K - \chi_{\partial K})_{\partial K}.$$

Recalling the definition (21) of the reconstruction operator $R_K$ applied to $\hat{w}_K$ proves (47).

(2) Proof of (48). Using the expression (16) for $n_{\#}^{(4)}$, we infer that

$$\hat{n}_{\#}^{(4)}(v, \hat{w}_h) - n_{\#}^{(4)}(v, w_{\mathcal{T}}) = \sum_{K \in \mathcal{T}} \sum_{i \in \{1:d\}} \left\{ \sum_{F \in \mathcal{F}_{\partial K}} \langle \gamma_{K,F}^{d}(\nabla \partial_i v), n_{F,i}(n_F \cdot \llbracket \nabla w_{\mathcal{T}} \rrbracket_F - \chi_F) \rangle_F - \sum_{F \in \mathcal{F}_{\partial K}^b} \langle \gamma_{K,F}^{d}(\nabla \partial_i v), n_{F,i} \chi_F \rangle_F \right\}.$$
Re-organizing the summation over the mesh interfaces yields

$$
\begin{align*}
\hat{n}_{\#}^{(4)}(v, \hat{w}_h) - n_{\#}^{(4)}(v, w_{\mathcal{T}}) &= \sum_{F \in \mathcal{F}^1} \sum_{i \in \{1:d\}} \sum_{j \in \{1,2\}} \langle \gamma_{K_j, F}^d(\nabla \partial_i v), n_{F,i}(\mathbf{n}_F \cdot \{\nabla w_{\mathcal{T}}\}_F^T - \chi_F) \rangle_F \\
&\quad - \sum_{F \in \mathcal{F}^b} \sum_{i \in \{1:d\}} \langle \gamma_{K, F}^d(\nabla \partial_i v), n_{F,i} \chi_F \rangle_F,
\end{align*}
$$

recalling the notation $F = \partial K_1 \cap \partial K_2$ for all $F \in \mathcal{F}^1$. Since the first term on the right-hand side vanishes owing to (17), we conclude that (48) holds true.
□

If the function $v$ is smooth, we have $\sum_{i \in \{1:d\}} \langle \gamma_{K,F}^d (\nabla \partial_i v), n_{F,i} \chi_F \rangle_F = \int_F (\partial_{nn} v) \chi_F \, \mathrm{d}s$. Thus, we expect that the right-hand side of (48) vanishes when applied with the first argument equal to the weak solution of (6) (type (II) BC's). The same property holds when considering the weak solution of (5) (type (I) BC's) since, in this case, the face component of the test function vanishes at the boundary faces. Les us now formalize these arguments.

**Lemma 4.2 (Identity for weak solution)** (i) Type (I) BC's: Assume the weak solution $u^{(1)}$ to (5) is in $V^{(4)}(\Omega)$. Then, we have

$$
\hat{n}_{\#}^{(4)}(u^{(1)}, \hat{w}_h) = n_{\#}^{(4)}(u^{(1)}, w_{\mathcal{T}}), \quad \forall \hat{w}_h \in \hat{V}_{h}^{k,(1)}. \qquad (49a)
$$

(ii) Type (ii) BC's: Assume the weak solution $u^{(II)}$ to (6) is in $V^{(4)}(\Omega)$. Then, we have

$$
\hat{n}_{\#}^{(4)}(u^{(\mathrm{II})}, \hat{w}_h) = n_{\#}^{(4)}(u^{(\mathrm{II})}, w_{\mathrm{T}}), \quad \forall \hat{w}_h \in \hat{V}_{h}^{k,(\mathrm{II})}. \qquad (49\mathrm{b})
$$

**Proof.** (1) Proof of (49a). The identity is a simple consequence of (48) and the fact that $\chi_F = 0$ for all $F \in \mathcal{F}^b$ whenever $\hat{w}_h = (w_{\mathcal{T}}, \chi_F) \in \hat{V}_h^{k,(1)}$.

(2) Proof of (49b). Let us prove that for all $F \in \mathcal{F}^b$, we have

$$
\sum_{i \in \{1:d\}} \langle \gamma_{K,F}^{d}(\nabla \partial_i u^{(\mathrm{II})}), n_{F,i} \chi_F \rangle_F = 0, \qquad (50)
$$

where $K \in \mathcal{T}$ is the unique mesh cell such that $F = \partial K \cap \partial \Omega$. Invoking [21, p. 17], we infer
that there is $\phi_\chi \in H_0^1(K) \cap H^2(K)$, $(\partial_n \phi_\chi)|_{F'} = 0$ for all $F' \in \mathcal{F}_{\partial K}^1$ and $(\partial_n \phi_\chi)|_F = \chi_F$.
Let $\tilde{\phi}_\chi$ denote the zero-extension of $\phi_\chi$ to $H_0^1(\Omega) \cap H^2(\Omega)$. Using $\tilde{\phi}_\chi$ as test function in the
weak formulation (6) and since $\Delta u \in W^{1,q}(\Omega)$ and $\tilde{\phi}_\chi \in W_0^{1,q'}(\Omega)$, we obtain

$$
\begin{align*}
0 &= (\nabla^2 u^{(\mathrm{II})}, \nabla^2 \tilde{\phi}_\chi)_\Omega + (\nabla \Delta u^{(\mathrm{II})}, \nabla \tilde{\phi}_\chi)_\Omega \\
  &= (\nabla^2 u^{(\mathrm{II})}, \nabla^2 \phi_\chi)_K + (\nabla \Delta u^{(\mathrm{II})}, \nabla \phi_\chi)_K = \sum_{i \in \{1:d\}} \left\{ (\nabla \partial_i u^{(\mathrm{II})}, \nabla (\partial_i \phi_\chi))_K + (\Delta \partial_i u^{(\mathrm{II})}, \partial_i \phi_\chi)_K \right\}.
\end{align*}
$$

Moreover, we have

$$
\begin{align}
\sum_{i \in \{1:d\}} \langle \gamma_{K,F}^d (\nabla \partial_i u^{(\mathrm{II})}), n_{F,i} \chi_F \rangle_F &= \sum_{i \in \{1:d\}} \langle \gamma_{K,F}^d (\nabla \partial_i u^{(\mathrm{II})}), (\partial_i \phi_\chi)|_F \rangle_F \tag{51} \\
&= \sum_{i \in \{1:d\}} \left\{ (\nabla \partial_i u^{(\mathrm{II})}, \nabla L_F^K (\partial_i \phi_\chi))_K + (\Delta \partial_i u^{(\mathrm{II})}, L_F^K (\partial_i \phi_\chi))_K \right\}, \nonumber
\end{align}
$$

where we used that the tangential derivative of $\phi_\chi$ vanishes on $F$ and the definition (9) of
the normal operator $\gamma_{K,F}^d$. Subtracting the above two relations gives

$$
\sum_{i \in \{1:d\}} \langle \gamma_{K,F}^{\mathrm{d}} (\nabla \partial_i u^{(\mathrm{II})}), n_{F,i} \chi_F \rangle_F = \sum_{i \in \{1:d\}} \left\{ (\nabla \partial_i u^{(\mathrm{II})}, \nabla \delta_i)_K + (\Delta \partial_i u^{(\mathrm{II})}, \delta_i)_K \right\}, \quad (52)
$$

with $\delta_i := L_F^K (\partial_i \phi_\xi) - \partial_i \phi_\xi$. We observe that $\delta_i \in W_0^{1,\rho'}(K) \cap L^2(K)$. Considering a sequence $(\delta_{i,\epsilon})_{\epsilon>0}$ in $C_0^\infty(K)$ which converges to $\delta_i$ in $W^{1,\rho'}(K) \cap L^2(K)$ as $\epsilon \to 0$ and since $\sum_{i \in \{1:d\}} \{ (\nabla(\partial_i u^{(\mathrm{II})}), \nabla \delta_{i,\epsilon})_K + (\Delta(\partial_i u^{(\mathrm{II})}), \delta_{i,\epsilon})_K \} = 0$ for all $\epsilon > 0$, we conclude that (50) holds true. Summing over the mesh boundary faces and invoking (48) proves (49b).
We close this section by stating a boundedness estimate on the bilinear form $\hat{n}_{\#}^{(4)}$. We omit the proof since it follows the arguments from [19, Lemma 3.2].

**Lemma 4.3 (Boundedness of $\hat{n}_{\#}^{(4)}$)** The following holds for all $v \in (\mathcal{V}^{(4)}(\Omega) + \mathbb{P}_{k+2}^{\text{g}}(\mathcal{T}))$ and all $\hat{w}_h \in \hat{V}_h^{k,(\text{II})}$:

$$|\hat{n}_{\#}^{(4)}(v, \hat{w}_h)| \le C \left( \sum_{K \in \mathcal{T}} h_K^{2d(\frac{1}{2}-\frac{1}{p})} \|\nabla^2 v\|_{L^p(K)}^2 + h_K^{2+2d(\frac{1}{2}-\frac{1}{q})} \|\nabla \Delta v\|_{L^q(K)}^2 \right)^{\frac{1}{2}} \left( \sum_{K \in \mathcal{T}} h_K^{-1} \|\chi_{\partial K} - \partial_n w_K\|_{\partial K}^2 \right)^{\frac{1}{2}}. \quad (53)$$

## 4.2 Bound on consistency error

We define the consistency errors $\delta_h^{(\text{I})} \in (\hat{V}_h^{k,(\text{I})})'$ and $\delta_h^{(\text{II})} \in (\hat{V}_h^{k,(\text{II})})'$ such that

$$\langle \delta_h^{(\text{I})}, \hat{w}_h \rangle := \ell(w_{\mathcal{T}}) - a_h(\hat{\mathcal{I}}_h^k(u^{(\text{I})}), \hat{w}_h), \quad \forall \hat{w}_h \in \hat{V}_h^{k,(\text{I})}, \quad (54\text{a})$$

$$\langle \delta_h^{(\text{II})}, \hat{w}_h \rangle := \ell(w_{\mathcal{T}}) - a_h(\hat{\mathcal{I}}_h^k(u^{(\text{II})}), \hat{w}_h), \quad \forall \hat{w}_h \in \hat{V}_h^{k,(\text{II})}, \quad (54\text{b})$$

where the brackets refer to the duality pairing between $(\hat{V}_h^{k,(\text{I})})'$ and $\hat{V}_h^{k,(\text{I})}$ or between $(\hat{V}_h^{k,(\text{II})})'$ and $\hat{V}_h^{k,(\text{II})}$ depending on the context. Recall that the spaces $\hat{V}_h^{k,(\text{I})}$ and $\hat{V}_h^{k,(\text{II})}$ are equipped with the norm $\|\cdot\|_{\hat{V}_h^k}$ defined in (32). Recall also that the $\|\cdot\|_{\#,(\text{K})}$-norm is defined in (42).

**Lemma 4.4 (Consistency)** Let $\delta_h$ denote either $\delta_h^{(\text{I})}$ or $\delta_h^{(\text{II})}$, let $u$ denote either $u^{(\text{I})}$ or $u^{(\text{II})}$, and let $\hat{V}_h^k$ denote either $\hat{V}_h^{k,(\text{I})}$ or $\hat{V}_h^{k,(\text{II})}$. Assume that $u \in \mathcal{V}^{(4)}(\Omega)$. The following holds:

$$\langle \delta_h, \hat{w}_h \rangle \le C \left( \sum_{K \in \mathcal{T}} \|u - \mathcal{L}_K^{k+2}(u)\|_{\#,(K)}^2 \right)^{\frac{1}{2}} \|\hat{w}_h\|_{\hat{V}_h^k}, \quad \forall \hat{w}_h \in \hat{V}_h^k. \quad (55)$$

**Proof.** Using the key identity (15) from Lemma 2.3 followed by the identity (49a) or the identity (49b) from Lemma 4.1 depending on the context, we infer that for all $\hat{w}_h \in \hat{V}_h^k$,

$$\ell(w_{\mathcal{T}}) = \langle f, w_{\mathcal{T}} \rangle_{W^{-1,q}, W_0^{1,q'}} = \sum_{K \in \mathcal{T}} (\nabla^2 u, \nabla^2 w_K)_K - n_{\#}^{(4)}(u, w_{\mathcal{T}}) = \sum_{K \in \mathcal{T}} (\nabla^2 u, \nabla^2 w_K)_K - \hat{n}_{\#}^{(4)}(u, \hat{w}_h).$$

Moreover, using the identity (47) from Lemma 4.1 gives

$$a_h(\hat{\mathcal{I}}_h^k(u), \hat{w}_h) = \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 J_K^{\text{HHO}}(u), \nabla^2 w_K)_K + S_{\partial K}(\hat{\mathcal{I}}_K^k(u), \hat{w}_K) \right\} - \hat{n}_{\#}^{(4)}(J_h^{\text{HHO}}(u), \hat{w}_h).$$

Thus, defining the function $\eta := u - J_h^{\text{HHO}}(u)$, i.e., $\eta|_K := u|_K - J_K^{\text{HHO}}(u|_K)$ for all $K \in \mathcal{T}$, we infer that

$$\langle \delta_h, \hat{w}_h \rangle = \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 \eta, \nabla^2 w_K)_K - S_{\partial K}(\hat{\mathcal{I}}_K^k(u), \hat{w}_K) \right\} - \hat{n}_{\#}^{(4)}(\eta, \hat{w}_h). \quad (56)$$

Let us denote by $\mathcal{T}_1$ the first two terms on the right-hand side and by $\mathcal{T}_2$ the third addend. We bound $\mathcal{T}_{1,K}$ by the same arguments as above, yielding

$$|\mathcal{T}_1| \le C \left( \sum_{K \in \mathcal{T}} \|\nabla^2 \eta\|_K^2 + \|\nabla^2 (u - \mathcal{L}_K^{k+2}(u))\|_K^2 \right)^{\frac{1}{2}} \|\hat{w}_h\|_{\hat{V}_h^k}.$$

Moreover, owing to (53), we have

$$|\mathcal{T}_2| \le C \left( \sum_{K \in \mathcal{T}} h_K^{2d(\frac{1}{2}-\frac{1}{p})} \|\nabla^2 \eta\|_{L^p(K)}^2 + h_K^{2+2d(\frac{1}{2}-\frac{1}{q})} \|\nabla \Delta \eta\|_{L^q(K)}^2 \right)^{\frac{1}{2}} \left( \sum_{K \in \mathcal{T}} h_K^{-1} \|\chi_{\partial K} - \partial_n w_K\|_{\partial K}^2 \right)^{\frac{1}{2}}.$$
Altogether, this implies that

$$|\langle \delta_h, \hat{w}_h \rangle| \le C \left( \sum_{K \in \mathcal{T}} \|\eta\|_{\sharp,K}^2 + \|\nabla^2(u - \mathcal{L}_K^{k+2}(u))\|_K^2 \right)^{\frac{1}{2}} \|\hat{w}_h\|_{\hat{V}_h^k}.$$

Invoking Lemma 3.8 completes the proof. □

**Remark 4.5 (Classical regularity assumption)** For completeness, let us briefly sketch how the consistency error is bounded under the classical regularity assumption $u \in H^{2+s}(\Omega)$ and $f \in H^{-2+s}(\Omega)$, $s \in (\frac{1}{2}, 1]$. As above, we let $u$ denote either $u^{(1)}$ or $u^{(\mathrm{II})}$, and $\hat{V}_h^k$ denote either $\hat{V}_h^{k,(\mathrm{I})}$ or $\hat{V}_h^{k,(\mathrm{II})}$. Starting from the identity from Lemma 2.1 and using that $\partial_{nn}u$ is single-valued at every mesh interface $\mathcal{F}_{\partial K}^{\mathrm{I}}$ and that either $\partial_{nn}u$ or $\chi_{\partial K}$ vanish at every mesh boundary face $F \in \mathcal{F}_{\partial K}^{\mathrm{b}}$ depending on the type of BC that is enforced, we infer that for all $\hat{w}_h \in \hat{V}_h^k$,

$$\langle f, w_{\mathcal{T}} \rangle_{H^{-2+s}, H_0^{2-s}} = \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 u, \nabla^2 w_K)_K - (\partial_{nn} u, \partial_n w_K - \chi_{\partial K})_{\partial K} \right\}.$$

Furthermore, using the identity $R_K \circ \hat{\mathcal{I}}_K^k = J_K^{\mathrm{hho}}$ from Lemma 3.8 together with the definition (21) of $R_K(\hat{w}_K)$ leads to

$$a_h(\hat{\mathcal{I}}_h^k(u), \hat{w}_h) = \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 J_K^{\mathrm{hho}}(u), \nabla^2 w_K)_K - (\partial_{nn} J_K^{\mathrm{hho}}(u), \partial_n w_K - \chi_{\partial K})_{\partial K} + S_{\partial K}(\hat{\mathcal{I}}_K^k(u), \hat{w}_K) \right\}.$$

Defining the function $\eta$ cellwise as $\eta|_K := u|_K - J_K^{\mathrm{hho}}(u)$ for all $K \in \mathcal{T}$, we infer that

$$\langle \delta_h, \hat{w}_h \rangle := \langle f, w_{\mathcal{T}} \rangle_{H^{-2+s}, H_0^{2-s}} - a_h(\hat{\mathcal{I}}_h^k(u), \hat{w}_h) \\ = \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 \eta, \nabla^2 w_K)_K - (\partial_{nn} \eta, \partial_n w_K - \chi_{\partial K})_{\partial K} - S_{\partial K}(\hat{\mathcal{I}}_K^k(u), \hat{w}_K) \right\}.$$

Setting $\|v\|_{\sharp,K}^2 := \|\nabla^2 v\|_K^2 + h_K \|\partial_{nn} v\|_{\partial K}^2$, it is then straightforward to establish that

$$\langle \delta_h, \hat{w}_h \rangle \le C \left( \sum_{K \in \mathcal{T}} \|u - \mathcal{L}_K^{k+2}(u)\|_{\sharp,K}^2 \right)^{\frac{1}{2}} \|\hat{w}_h\|_{\hat{V}_h^k}.$$

## 4.3 Energy-error estimate

We are now ready to establish our main convergence result bounding the error in the energy norm. For completeness, improved error estimates in weaker norms are outlined in Section 4.5. Let us define the discrete error $\hat{e}_h := \hat{\mathcal{I}}_h^k(u) - \hat{u}_h \in \hat{V}_h^k$ with local components $\hat{e}_K := \hat{\mathcal{I}}_K^k(u) - \hat{u}_K \in \hat{V}_K^k$ for all $K \in \mathcal{T}$.

**Theorem 4.6 (Error estimate)** Let $u$ denote either the weak solution $u^{(1)}$ of (5) or the weak solution $u^{(\mathrm{II})}$ of (6). Let $\hat{u}_h$ denote either the discrete solution $\hat{u}_h^{(1)}$ of (30) or the discrete solution $\hat{u}_h^{(\mathrm{II})}$ of (31), respectively. Under Assumption 1.1, the following holds:

$$\sum_{K \in \mathcal{T}} \left\{ \|u - R_K(\hat{u}_K)\|_{\sharp,K}^2 + S_{\partial K}(\hat{e}_K, \hat{e}_K) \right\} \le C \sum_{K \in \mathcal{T}} \|u - \mathcal{L}_K^{k+2}(u)\|_{\sharp,K}^2. \quad (57)$$

Moreover, assuming $u \in H^{2+s}(\Omega)$ with $s > 0$, and letting $t := \min\{s, k + 1\}$ as well as $\mathfrak{J}_t := 1$ if $t < 1 + d(\frac{1}{2} - \frac{1}{q})$ and $\mathfrak{J}_t := 0$ otherwise, we have

$$\sum_{K \in \mathcal{T}} \left\{ \|u - R_K(\hat{u}_K)\|_{\sharp,K}^2 + S_{\partial K}(\hat{e}_K, \hat{e}_K) \right\} \le C \sum_{K \in \mathcal{T}} \left( h_K^t |u|_{H^{t+2}(K)} + \mathfrak{J}_t h_K^{1+d(\frac{1}{2}-\frac{1}{q})} \|\nabla \Delta v\|_{L^q(K)}^2 \right)^2. \quad (58)$$
**Proof.** (1) Let $\delta_h$ denote either $\delta_h^{(I)}$ or $\delta_h^{(II)}$, and let $\hat{V}_h^k$ denote either $\hat{V}_h^{k,(I)}$ or $\hat{V}_h^{k,(II)}$ depending on the context. Since $a_h(\hat{e}_h, \hat{e}_h) = -\langle \delta_h, \hat{e}_h \rangle$, the coercivity of the bilinear form $a_h$ implies that

$$ \alpha \|\hat{e}_h\|_{\hat{V}_h^k}^2 \le a_h(\hat{e}_h, \hat{e}_h) = -\langle \delta_h, \hat{e}_h \rangle \le \|\delta_h\|_{(\hat{V}_h^k)'} \|\hat{e}_h\|_{\hat{V}_h^k}, $$

so that $\|\hat{e}_h\|_{\hat{V}_h^k} \le \frac{1}{\alpha} \|\delta_h\|_{(\hat{V}_h^k)'}$. Since $\sum_{K \in \mathcal{T}} \{\|\nabla^2 R_K(\hat{e}_K)\|_{\sharp,K}^2 + S_{\partial K}(\hat{e}_K, \hat{e}_K)\} \le C \sum_{K \in \mathcal{T}} \{\|\nabla^2 R_K(\hat{e}_K)\|_{\sharp,K}^2 + S_{\partial K}(\hat{e}_K, \hat{e}_K)\} \le C \|\hat{e}_h\|_{\hat{V}_h^k}^2$ owing to the inverse inequality (35) and the upper bound in (25), we infer from Lemma 4.4 that

$$ \sum_{K \in \mathcal{T}} \{ \| \nabla^2 R_K(\hat{e}_K) \|_{\sharp,K}^2 + S_{\partial K}(\hat{e}_K, \hat{e}_K) \} \le C \sum_{K \in \mathcal{T}} \| u - \mathcal{L}_K^{k+2}(u) \|_{\sharp,K}^2. $$

Since $u - R_K(\hat{u}_K) = (u - J_K^{\text{HHO}}(u)) + R_K(\hat{e}_K)$, the triangle inequality combined with Lemma 3.8 and the above bound proves (57).

(2) The estimate (58) results from (57) and the approximation properties of $\mathcal{L}_K^{k+2}$ established in Lemma 3.6. $\square$

**Remark 4.7 (Variant)** *Invoking inverse inequalities shows that for all $K \in \mathcal{T}$,*

$$ \|u_K - R_K(\hat{u}_K)\|_{\sharp,K} \le C \|\nabla^2(u_K - R_K(\hat{u}_K))\|_K \le C S_{\partial K}(\hat{u}_K, \hat{u}_K), \quad (59) $$

where the second bound follows from $(\nabla^2(u_K - R_K(\hat{u}_K)), \nabla^2 w)_K = (\partial_n u_K - \gamma_{\partial K}, \partial_{nn} w)_{\partial K}$ for all $w \in \mathbb{P}_{k+2}(K)^{\perp}$, a discrete trace inequality, $\partial_{nn} w \in \mathbb{P}_k(\mathcal{F}_{\partial K})$, and the definition of the stabilization operator. Moreover, the triangle inequality combined with the bound (43) from Lemma 3.8 implies that

$$ \frac{1}{2} S_{\partial K}(\hat{u}_K, \hat{u}_K) \le S_{\partial K}(\hat{\mathcal{I}}_K^k(u), \hat{\mathcal{I}}_K^k(u)) + S_{\partial K}(\hat{e}_K, \hat{e}_K) \le C \|u - \mathcal{L}_K^{k+2}(u)\|_{\sharp,K} + S_{\partial K}(\hat{e}_K, \hat{e}_K). $$

This leads to the following variant of the above energy-error estimate:

$$ \sum_{K \in \mathcal{T}} \{ \| \nabla^2 (u - u_K) \|_K^2 + S_{\partial K}(\hat{u}_K, \hat{u}_K) \} \le C \sum_{K \in \mathcal{T}} \| u - \mathcal{L}_K^{k+2}(u) \|_{\sharp,K}^2, \quad (60) $$

which follows from (57), the triangle inequality, and the above bounds.

## 4.4 Energy-error estimate for $C^0$-IPDG methods

To illustrate that the key identity (15) from Lemma 2.3 has a broader applicability than $C^0$-HHO methods, we briefly outline here how this identity can be used in the error analysis of $C^0$-IPDG methods. Focusing for brevity on type (II) BC's only, the $C^0$-IPDG bilinear form is defined as follows for all $v_h, w_h \in V_{T,0}^{k+2} := \mathbb{P}_{k+2}^g(T) \cap H_0^1(\Omega)$:

$$ a_h^{\text{IPDG}}(v_h, w_h) := \sum_{K \in \mathcal{T}} (\nabla^2 v_h, \nabla^2 w_h)_K + \sum_{F \in \mathcal{F}^i} \varpi_F h_F^{-1} (\llbracket \partial_n v_h \rrbracket_F, \llbracket \partial_n w_h \rrbracket_F)_F \\ - \sum_{F \in \mathcal{F}^i} (\llbracket \partial_{nn} v_h \rrbracket_F, \llbracket \partial_n w_h \rrbracket_F)_F - \sum_{F \in \mathcal{F}^i} (\llbracket \partial_n v_h \rrbracket_F, \llbracket \partial_{nn} w_h \rrbracket_F)_F, \quad (61) $$

where all the normal derivatives are understood to be oriented by $n_F$ and with the user-defined penalty parameter $\varpi_F > 0$. For type (I) BC's, the three summations over the mesh interfaces are realized over the whole set of mesh faces. The discrete problem consists of finding $u_h \in V_{T,0}^{k+2}$ such that

$$ a_h^{\text{IPDG}}(v_h, w_h) = \ell(w_h) := \langle f, w_h \rangle_{W^{-1,q}, W_0^{1,q'}}, \quad \forall w_h \in V_{T,0}^{k+2}. \quad (62) $$
The stability analysis reveals that the bilinear form $a_h^{\text{IPDG}}$ is coercive on $V_{T,0}^{k+2}$ if the penalty parameters $\omega_F$ are large enough; see, e.g., [2]. The coercivity norm is $\|w_h\|_{V_{T,0}^{k+2}}^2 := \sum_{K \in \mathcal{T}} \|\nabla^2 w_h\|_K^2 + \sum_{F \in \mathcal{F}^i} h_F^{-1} \|\llbracket \partial_n w_h \rrbracket_F\|_F^2$.

The novelty here consists in bounding the consistency error under Assumption 1.1. Let the bilinear form $n_{\#}^{\text{IPDG}}$ be defined as $n_{\#}^{(4)}$, except that the contributions of all the mesh boundary faces $F \in \mathcal{F}^b$ are discarded. The arguments in the proof of (49b) show that $n_{\#}^{\text{IPDG}}(u, w_h) = n_{\#}^{(4)}(u, w_h)$ for all $w_h \in V_{T,0}^{k+2}$. Therefore, using the key identity (15) from Lemma 2.3 gives

$$ \langle f, w_h \rangle_{W^{-1,q}, W_0^{1,q'}} = \sum_{K \in \mathcal{T}} (\nabla^2 u, \nabla^2 w_h)_K - n_{\#}^{(4)}(u, w_h) = \sum_{K \in \mathcal{T}} (\nabla^2 u, \nabla^2 w_h)_K - n_{\#}^{\text{IPDG}}(u, w_h), \quad (63) $$

for all $w_h \in V_{T,0}^{k+2}$. Moreover, using the same proof as for (47), we establish the following identity for all $v_h, w_h \in V_{T,0}^{k+2}$:

$$ n_{\#}^{\text{IPDG}}(v_h, w_h) = \sum_{F \in \mathcal{F}^i} (\llbracket \partial_{nn} v_h \rrbracket_F, \llbracket \partial_n w_h \rrbracket_F)_F. \quad (64) $$

Putting everything together and setting $\eta := u - \mathcal{L}_h^{k+2}(u)$, we can express the consistency error as follows:

$$ \begin{aligned} \langle \delta_h, w_h \rangle &:= \ell(w_h) - a_h^{\text{IPDG}}(\mathcal{L}_h^{k+2}(u), w_h) \\ &= \sum_{K \in \mathcal{T}} (\nabla^2 \eta, \nabla^2 w_h)_K + \sum_{F \in \mathcal{F}^i} \omega_F h_F^{-1} (\llbracket \partial_n \eta \rrbracket_F, \llbracket \partial_n w_h \rrbracket_F)_F \\ &\quad - n_{\#}^{\text{IPDG}}(\eta, w_h) - \sum_{F \in \mathcal{F}^i} (\llbracket \partial_n \eta \rrbracket_F, (\llbracket \partial_{nn} w_h \rrbracket_F)_F), \end{aligned} \quad (65) $$

where we used that $\llbracket \partial_n u \rrbracket_F = 0$ for all $F \in \mathcal{F}^i$. Using the same $\|\cdot\|_{\#}$-norm as for the $C^0$-HHO method, and invoking the Cauchy-Schwarz inequality, Lemma 4.3 to bound $n_{\#}^{\text{IPDG}}(\eta, w_h)$, and the estimate (39) to estimate $\partial_n \eta$ at the mesh interfaces, we infer that

$$ \langle \delta_h, w_h \rangle \le C \left( \sum_{K \in \mathcal{T}} \|u - \mathcal{L}_K^{k+2}(u)\|_{\#}^2 \right)^{\frac{1}{2}} \|w_h\|_{V_{T,0}^{k+2}} \quad \forall w_h \in V_{T,0}^{k+2}. \quad (66) $$

Finally, proceeding as above for the $C^0$-HHO method readily leads to the same error estimates.

## 4.5 Improved error estimates in weaker norms

To derive error estimates in weaker norms, we assume that the following regularity pickup holds true: There exists a constant $C_\Omega$ such that for all $g \in H^\sigma(\Omega)$, $\sigma \in \{-1, 0\}$, the adjoint solution such that $\Delta^2 \zeta_g = g$ in $\Omega$ with either type (I) or (II) BC's satisfies the bound

$$ \|\zeta_g\|_{H^{4+\sigma}(\Omega)} \le C_{\Omega} \|g\|_{H^{\sigma}(\Omega)}, \quad \forall \sigma \in \{-1, 0\}. \quad (67) $$

A sufficient condition for (67) to hold with type (I) BC's is that the domain $\Omega$ is convex; see [23, p. 182]. In what follows, we assume $f \in L^2(\Omega)$ so that the above regularity pickup implies that $u \in H^4(\Omega)$. Hence, the setting of Assumption 1.1 is not needed here, and we simply use the $\|\cdot\|_{\#}$-norm already considered in Remark 4.5, that is, $\|v\|_{\#}^2 := \|\nabla^2 v\|_K^2 + h_K \|\partial_{nn} v\|_{\partial K}^2$ for all $K \in \mathcal{T}$. We also set $\|v\|_{\#}^2 := \sum_{K \in \mathcal{T}} \|v|_K\|_{\#}^2$. Let $\tilde{h} := \max_{K \in \mathcal{T}} h_K$ denote the (global) mesh size and let $\tilde{h}|_K$ be the piecewise constant function such that $\tilde{h}|_K := h_K$ for all $K \in \mathcal{T}$.
**Theorem 4.8 (L²-and H¹-error estimate)** Let $u$ denote either the weak solution $u^{(1)}$ of (5) or the weak solution $u^{(II)}$ of (6). Let $\hat{u}_h$ denote either the discrete solution $\hat{u}_h^{(1)}$ of (30) or the discrete solution $\hat{u}_h^{(II)}$ of (31), respectively. Assume $f \in L^2(\Omega)$. Then, if the regularity estimate (67) holds true with $\sigma = 0$, we have

$$ \sum_{K \in \mathcal{T}} \|u - R_K(\hat{u}_K)\|_K^2 \le C h^{2 \min(k+1,2)} \left( \|\hat{e}_h\|_{\hat{V}_h^k}^2 + \|u - \mathcal{L}_h^{k+2}(u)\|_{\sharp,h}^2 + \|\tilde{h}^2(f - \Pi_h^{k-2}(f))\|_{\Omega}^2 \right), \quad (68) $$

and if the regularity estimate (67) holds true with $\sigma = -1$, we have

$$ \sum_{K \in \mathcal{T}} \|\nabla(u - R_K(\hat{u}_K))\|_K^2 \le C h^2 \left( \|\hat{e}_h\|_{\hat{V}_h^k}^2 + \|u - \mathcal{L}_h^{k+2}(u)\|_{\sharp,h}^2 + \|\tilde{h}^2(f - \Pi_h^{k-2}(f))\|_{\Omega}^2 \right). \quad (69) $$

To prove Theorem 4.8, we need to consider a novel HHO interpolation operator. For all $K \in \mathcal{T}$, let $\mathcal{C}_K^{k+2} : H^2(K) \to \mathbb{P}_{k+2}(K)$ be the local interpolation operator associated with the canonical hybrid finite element in $K$; see [17, section 7.6] (this is the finite element considered in the first step of the discrete de Rham diagram). Recall that for $d = 3$ and a function $v \in H^2(K)$, $\mathcal{C}_K^{k+2}(v)$ is uniquely defined by the following properties:

$$ \mathcal{C}_K^{k+2}(v)(z) = v(z), \quad \forall z \in \mathcal{V}_K, \quad (70a) $$

$$ (\mathcal{C}_K^{k+2}(v), \xi_E)_E = (v, \xi_E)_E, \quad \forall \xi_E \in \mathbb{P}_k(E), \forall E \in \mathcal{E}_K, \quad (70b) $$

$$ (\mathcal{C}_K^{k+2}(v), \xi_F)_F = (v, \xi_F)_F, \quad \forall \xi_F \in \mathbb{P}_{k-1}(F), \forall F \in \mathcal{F}_K, k \ge 1, \quad (70c) $$

$$ (\mathcal{C}_K^{k+2}(v), \xi_K)_K = (v, \xi_K)_K, \quad \forall \xi_K \in \mathbb{P}_{k-2}(K), k \ge 2, \quad (70d) $$

where $\mathcal{V}_K$, $\mathcal{E}_K$, and $\mathcal{F}_K$ collect the vertices, edges, and faces of $K$, respectively. For $d = 2$, $\mathcal{C}_K^{k+2}(v)$ is defined similarly by using (70a), (70b), and (70d) with $\xi_K \in \mathbb{P}_{k-1}(K)$, $k \ge 1$. The following holds for all $K \in \mathcal{T}$ and all $v \in H^3(K)$:

$$ \|v - \mathcal{C}_K^{k+2}(v)\|_{\sharp,K} \le C \|v - \mathcal{L}_K^{k+2}(v)\|_{\sharp,K}. \quad (71) $$

This follows from the triangle inequality, inverse inequalities yielding $\|\mathcal{C}_K^{k+2}(v) - \mathcal{L}_K^{k+2}(v)\|_{\sharp,K} \le C \|\nabla^2(\mathcal{C}_K^{k+2}(v) - \mathcal{L}_K^{k+2}(v))\|_K$, and the $H^2$-stability of $\mathcal{C}_K^{k+2}$ together with the fact that $\mathcal{C}_K^{k+2}(v) - \mathcal{L}_K^{k+2}(v) = \mathcal{C}_K^{k+2}(v - \mathcal{L}_K^{k+2}(v))$.

We then define the new local HHO reduction operator $\hat{\mathcal{P}}_K^k : H^2(K) \to \hat{V}_K^k$ such that, for all $v \in H^2(K)$,

$$ \hat{\mathcal{P}}_K^k(v) := (\mathcal{C}_K^{k+2}(v), \Pi_{\partial K}^k(\mathbf{n}_K \cdot \nabla v)) \in \hat{V}_K^k. \quad (72) $$

as well as the new local HHO interpolation operator such that

$$ T_K^{\text{HHO}} := R_K \circ \hat{\mathcal{P}}_K^k : H^2(K) \to \mathbb{P}_{k+2}(K). \quad (73) $$

The main motivation for the above construction is the following result.

**Lemma 4.9 (New HHO interpolation operator)** The following holds for all $K \in \mathcal{T}$ and all $v \in H^2(K)$:

$$ (\nabla^2 T_K^{\text{HHO}}(v), \nabla^2 w)_K = (\nabla^2 v, \nabla^2 w)_K, \quad \forall w \in \mathbb{P}_{k+2}(K)^{\perp}. \quad (74) $$

Moreover, we have for all $K \in \mathcal{T}$ and all $v \in H^3(K)$,

$$ \|v - T_K^{\text{HHO}}(v)\|_{\sharp,K}^2 + S_{\partial K}(\hat{\mathcal{P}}_K^k(v), \hat{\mathcal{P}}_K^k(v)) \le C \|v - \mathcal{L}_K^{k+2}(v)\|_{\sharp,K}^2. \quad (75) $$

**Proof.** (1) Proof of (74). Using the definition of the reconstruction operator $R_K$, the identity $\nabla \cdot \nabla^2 = \nabla \Delta$ and integration by parts gives

$$ \begin{align*} (\nabla^2 T_K^{\text{HHO}}(v), \nabla^2 w)_K & = (\nabla^2 \mathcal{C}_K^{k+2}(v), \nabla^2 w)_K - (\partial_n \mathcal{C}_K^{k+2}(v) - \Pi_{\partial K}^k(\mathbf{n}_K \cdot \nabla v), \partial_{nn} w)_{\partial K} \\ & = -(\nabla \mathcal{C}_K^{k+2}(v), \nabla \Delta w)_K + (\Pi_{\partial K}^k(\mathbf{n}_K \cdot \nabla v), \partial_{nn} w)_{\partial K} + (\partial_t \mathcal{C}_K^{k+2}(v), \partial_{nt} w)_{\partial K} \\ & = (\mathcal{C}_K^{k+2}(v), \Delta^2 w)_K - (\mathcal{C}_K^{k+2}(v), \mathbf{n}_K \cdot \nabla \Delta w)_{\partial K} + (\Pi_{\partial K}^k(\mathbf{n}_K \cdot \nabla v), \partial_{nn} w)_{\partial K} + (\partial_t \mathcal{C}_K^{k+2}(v), \partial_{nt} w)_{\partial K}. \end{align*} $$
Integrating by parts the last term on the right-hand side gives (for $d = 3$)

$$ (\partial_t C_K^{k+2}(v), \partial_{nt} w)_{\partial K} = \sum_{F \in \mathcal{F}_K} \left\{ -(C_K^{k+2}(v), \partial_t(\partial_{nt} w))_F + \sum_{E \subset \partial F} (C_K^{k+2}(v), n_F \cdot \partial_{nt} w)_E \right\}. $$

Owing to (70) and since $w \in \mathbb{P}_{k+2}(K)$, we infer that

$$ (\nabla^2 T_K^{\text{HHO}}(v), \nabla^2 w)_K = (v, \Delta^2 w)_K - (v, n_K \cdot \nabla \Delta w)_{\partial K} + (n_K \cdot \nabla v, \partial_{nn} w)_{\partial K} + (\partial_t v, \partial_{nt} w)_{\partial K}, $$

whence we deduce from integration by parts that (74) holds true.

(2) Proof of (75). Proceeding as in Lemma 3.8 yields

$$ \|v - T_K^{\text{HHO}}(v)\|_{\#K}^2 + S_{\partial K}(\hat{P}_K^k(v), \hat{P}_K^k(v)) \le C \|v - C_K^{k+2}(v)\|_{\#K}^2, $$

and we conclude owing to (71). □

We are now ready to prove Theorem 4.8. The proof proceeds similarly to that of [11, Thm. 10] on the $L^2$-error estimate for second-order elliptic PDEs, but requires some non-trivial adaptations using the above tools. We detail the proof of (68) and only outline the minor changes needed to prove (69).

**Proof of (68)** Let $\hat{P}_h^k : H^2(\Omega) \to \hat{V}_h^k$ be defined such that, for all $v \in H^2(\Omega)$,

$$ \hat{P}_h^k(v) := (C_h^{k+2}(v), (\Pi_F^k(n_F \cdot (\nabla v)|_F))_{F \in \mathcal{F}}) \in \hat{V}_h^k, $$

where $C_h^{k+2} : H^2(\Omega) \to \mathcal{P}_{k+2}(\mathcal{T})$ is the global interpolation operator associated with the canonical hybrid finite element. Notice that the local components of $\hat{P}_h^k(v)$ are $\hat{P}_K^k(v|_K)$ for all $K \in \mathcal{T}$. Let $\hat{u}_h := (u_\mathcal{T}, \gamma_\mathcal{F}) \in \hat{V}_h^k$ be the discrete solution. Let us set

$$ \hat{e}_h := \hat{u}_h - \hat{P}_h^k(u), $$

so that, setting $\hat{e}_h := (e_\mathcal{T}, \theta_\mathcal{F})$, we have $e_K = u_K - C_K^{k+2}(u)$ for all $K \in \mathcal{T}$ and $\theta_F = \gamma_F - \Pi_F^k(n_F \cdot \nabla u)$ for all $F \in \mathcal{F}$.

(1) Let $\zeta_e$ be the adjoint solution such that $\Delta^2\zeta_e = e_\mathcal{T}$ in $\Omega$ with the appropriate type of BC's. Owing to our regularity assumption with $\sigma = 0$, we have $\|\zeta_e\|_{H^4(\Omega)} \le C_\Omega\|e_\mathcal{T}\|_{\Omega}$. Integration by parts gives

$$ \|e_\mathcal{T}\|_{\Omega}^2 = (e_\mathcal{T}, \Delta^2 \zeta_e)_\Omega = -(\nabla e_\mathcal{T}, \nabla \Delta \zeta_e)_\Omega = -(\nabla e_\mathcal{T}, \nabla \cdot \nabla^2 \zeta_e)_\Omega \\ = \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 e_\mathcal{T}, \nabla^2 \zeta_e)_K - (\partial_n e_\mathcal{T}, \partial_{nn} \zeta_e)_{\partial K} - (\partial_t e_\mathcal{T}, \partial_{nt} \zeta_e)_{\partial K} \right\}. $$

The last summation on the right-hand side vanishes for both types of BC's. Therefore, we have

$$ \begin{aligned} \|e_\mathcal{T}\|_{\Omega}^2 &= \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 e_\mathcal{T}, \nabla^2 \zeta_e)_K - (\partial_n e_\mathcal{T}, \partial_{nn} \zeta_e)_{\partial K} \right\} \\ &= \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 e_\mathcal{T}, \nabla^2 \zeta_e)_K - (\partial_n e_\mathcal{T} - \theta_\mathcal{F}, \partial_{nn} \zeta_e)_{\partial K} \right\} \\ &= \sum_{K \in \mathcal{T}} \left\{ (\nabla^2 e_\mathcal{T}, \nabla^2 T_K^{\text{HHO}}(\zeta_e))_K - (\partial_n e_\mathcal{T} - \theta_\mathcal{F}, \partial_{nn} \zeta_e)_{\partial K} \right\}, \end{aligned} $$

where we used that $\sum_{K \in \mathcal{T}} (\theta_\mathcal{F}, \partial_{nn} \zeta_e)_{\partial K} = 0$ for both types of BC's on the second line and the identity (74) on the third line. Since

$$ (\nabla^2 R_K(\hat{e}_K), \nabla^2 T_K^{\text{HHO}}(\zeta_e))_K = (\nabla^2 e_K, \nabla^2 T_K^{\text{HHO}}(\zeta_e))_K - (\partial_n e_K - \theta_{\partial K}, \partial_{nn} T_K^{\text{HHO}}(\zeta_e))_{\partial K}, $$
straightforward algebra using the definition of $\hat{e}_h$, that $\hat{u}_h$ solves the discrete HHO problem, and that $(f, \zeta_e)_{\Omega} = (\nabla^2 u, \nabla^2 \zeta_e)_{\Omega}$ gives $\|e_{\mathcal{T}}\|_{\Omega}^2 = T_1 + T_2 - T_3$ with

$$T_1 := - \sum_{K \in \mathcal{T}} \left\{ S_{\partial K}(\hat{e}_K, \hat{P}_K^k(\zeta_e)) - (\partial_n e_K - \theta_{\partial K}, \partial_{nn}(\zeta_e - T_K^{\text{HHO}}(\zeta_e)))_{\partial K} \right\},$$

$$T_2 := (\nabla^2 u, \nabla^2 \zeta_e)_{\Omega} - a_h(\hat{P}_h^k(u), \hat{P}_h^k(\zeta_e)),$$

$$T_3 := (f, \zeta_e - C_h^{k+2}(\zeta_e))_{\Omega} = (f - \Pi_h^{k-2}(f), \zeta_e - C_h^{k+2}(\zeta_e))_{\Omega},$$

where the last equality for $T_3$ follows from (70d) (with the convention that $\Pi_h^{k-2}(f) = 0$ for $k \in \{0, 1\}$). It remains to bound the above three terms.

(2) Owing to the Cauchy–Schwarz inequality and (75), we have

$$|T_1| \le C \|\hat{e}_h\|_{\hat{V}_h^k} \|\zeta_e - \mathcal{L}_h^{k+2}(\zeta_e)\|_{\sharp,h}.$$

Moreover, since $\|\zeta_e - \mathcal{L}_h^{k+2}(\zeta_e)\|_{\sharp,h} \le Ch^{\min(k+1,2)}|\zeta_e|_{H^{\min(k+3,4)}(\Omega)}$ and $|\zeta_e|_{H^{\min(k+3,4)}(\Omega)} \le C_{\Omega}\|e_{\mathcal{T}}\|_{\Omega}$ by our assumption on regularity pickup, we infer that

$$|T_1| \le Ch^{\min(k+1,2)} \|\hat{e}_h\|_{\hat{V}_h^k} \|e_{\mathcal{T}}\|_{\Omega}.$$

Furthermore, using the definition of $a_h$ and the identity (74) yields

$$T_2 = \sum_{K \in \mathcal{T}} \left\{ (\nabla^2(u - T_K^{\text{HHO}}(u)), \nabla^2(\zeta_e - T_K^{\text{HHO}}(\zeta_e)))_K - S_{\partial K}(\hat{P}_K^k(u), \hat{P}_K^k(\zeta_e)) \right\}.$$

Using the Cauchy–Schwarz inequality and the same arguments as for $T_1$ gives

$$|T_2| \le Ch^{\min(k+1,2)} \|u - \mathcal{L}_h^{k+2}(u)\|_{\sharp,h} \|e_{\mathcal{T}}\|_{\Omega}.$$

Finally, we have

$$|T_3| \le \|\tilde{h}^2(f - \Pi_h^{k-2}(f))\|_{\Omega} \|\tilde{h}^{-2}(\zeta_e - C_h^{k+2}(\zeta_e))\|_{\Omega} \le Ch^{\min(k+1,2)} \|\tilde{h}^2(f - \Pi_h^{k-2}(f))\|_{\Omega} \|e_{\mathcal{T}}\|_{\Omega}.$$

Putting everything together gives

$$\|e_{\mathcal{T}}\|_{\Omega} \le Ch^{\min(k+1,2)} \left( \|\hat{e}_h\|_{\hat{V}_h^k} + \|u - \mathcal{L}_h^{k+2}(u)\|_{\sharp,h} + \|\tilde{h}^2(f - \Pi_h^{k-2}(f))\|_{\Omega} \right). \quad (76)$$

(3) Using the triangle inequality and the definition of $e_{\mathcal{T}}$ gives

$$\sum_{K \in \mathcal{T}} \|u - R_K(\hat{u}_K)\|_K^2 \le C \sum_{K \in \mathcal{T}} \left\{ \|e_K\|_K^2 + \|u - R_K(\hat{u}_K)\|_K^2 + \|u - C_K^{k+2}(u)\|_K^2 \right\}.$$

The first term on the right-hand side is bounded using (76). The second term is bounded by the Poincaré inequality, (59), and the triangle inequality giving

$$\|u_K - R_K(\hat{u}_K)\|_K \le Ch_K^2 \|\nabla^2(u_K - R_K(\hat{u}_K))\|_K \le Ch_K^2 S_{\partial K}(\hat{u}_K, \hat{u}_K) \\ \le Ch_K^2 (S_{\partial K}(\hat{e}_K, \hat{e}_K) + \|u - \mathcal{L}_K^{k+2}(u)\|_{\sharp,K}).$$

Finally, observing that $\mathcal{L}_K^1$ leaves $(u - C_K^{k+2}(u))$ invariant, invoking (38) and (71) gives

$$\|u - C_K^{k+2}(u)\|_K \le Ch_K^2 \|\nabla^2(u - C_K^{k+2}(u))\|_K \le Ch_K^2 \|u - \mathcal{L}_K^{k+2}(u)\|_{\sharp,K}.$$

Putting everything together proves (68) since $h^2 \le Ch^{\min(k+1,2)}$. $\square$

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

**Proof of (69)** The only salient change with respect to the previous proof is the definition of the adjoint solution which now satisfies $\Delta^2 \zeta_e = \Delta e_{\mathcal{T}} \in H^{-1}(\Omega)$ with the appropriate type of BC's (recall that $e_{\mathcal{T}} \in H_0^1(\Omega)$ by definition). Proceeding as above and invoking the regularity pickup assumption (67) with $\sigma = -1$ then gives

$$\|\nabla e_{\mathcal{T}}\|_{\Omega}^2 = \langle \Delta e_{\mathcal{T}}, e_{\mathcal{T}} \rangle_{H^{-1}, H_0^1} \le Ch \left( \|\hat{e}_h\|_{\hat{V}_h^k} + \|u - \mathcal{L}_h^{k+2}(u)\|_{\sharp,h} + \|\tilde{h}^2(f - \Pi_h^{k-2}(f))\|_{\Omega} \right) |\zeta_e|_{H^3(\Omega)} \\ \le Ch \left( \|\hat{e}_h\|_{\hat{V}_h^k} + \|u - \mathcal{L}_h^{k+2}(u)\|_{\sharp,h} + \|\tilde{h}^2(f - \Pi_h^{k-2}(f))\|_{\Omega} \right) \|\Delta e_{\mathcal{T}}\|_{H^{-1}(\Omega)}.$$

The assertion finally follows from $\|\Delta e_{\mathcal{T}}\|_{H^{-1}(\Omega)} = \|\nabla e_{\mathcal{T}}\|_{\Omega}$.
**Remark 4.10 (Variants and decay rates)** Proceeding as in Remark 4.7 shows that the bounds (68) and (69) hold true also with the left-hand side replaced by $\sum_{K \in \mathcal{T}} \|u - u_K\|_K^2$ and $\sum_{K \in \mathcal{T}} \|\nabla(u - u_K)\|_K^2$, respectively. Concerning the right-hand side of (68) and (69), the first term is controlled by the second term owing to Theorem 4.6, so that both terms converge at rate $O(h^{2(k+1)})$. Moreover, the last term converges at rate $O(h^{2\max(2,k+1)})$ (i.e., at the same rate for $k \ge 1$ and at one order faster for $k = 0$). Taking into account the scaling factor in front of the parenthesis and taking the square root, we infer that the $L^2$-error converges at rate $O(h^2)$ for $k = 0$, which is suboptimal by one order in $h$, whereas it converges at the optimal rate $O(h^{k+3})$ for all $k \ge 1$. Finally, the $H^1$-error converges at the optimal rate $O(h^{k+2})$ for all $k \ge 0$.

# 5 Numerical examples

In this section, we present numerical examples illustrating our theoretical results on the convergence of the $C^0$-HHO method. We also compare the numerical performance of the $C^0$-HHO method to other (classical) $C^0$-conforming methods from the literature. All the computations were run with Matlab R2021b on the Cleps platform at INRIA Paris using 12 cores (all offering the same computational performance in terms of RAM and frequency), and all the linear systems after static condensation (if applicable) are solved using the backslash function (invoking Cholesky's factorization).

## 5.1 Convergence rates for the $C^0$-HHO method

We select $f$ and type (I) BC's on $\Omega := (0, 1)^2$ so that the exact solution is

$$u(x, y) := \sin(\pi x)^2 \sin(\pi y)^2 + e^{-(x-0.5)^2-(y-0.5)^2}.$$

We employ the polynomial degrees $k \in \{0, 1, 2, 3, 4\}$ and a sequence of successively refined triangular meshes consisting of $\{32, 128, 512, 2048, 8192, 32768\}$ cells. Despite an hp-error analysis falls beyond the present scope, we weigh the stabilization terms in (22) by replacing $h_K^{-1}$ by $(k+1)^2 h_K^{-1}$ for all $K \in \mathcal{T}$.

Let us first verify the convergence rates obtained with the $C^0$-HHO method with $k \in \{0, 1, 2, 3, 4\}$. We consider a sequence of successively refined triangular meshes. We measure errors in the (broken) $H^2$-seminorm, in the $H^1$-seminorm, in the $L^2$-norm, and in the stabilization seminorm. The first three errors are evaluated using the reconstruction of the HHO solution cellwise. The errors are reported in Figure 1, and the rates are reported in Table 1 as a function of DoFs$^{1/2}$, where DoFs denotes the total number of globally coupled discrete unknowns (that is, the face unknowns in addition to the cell unknowns except the bubble functions). We observe that the $H^2$-error, $H^1$-error, and stabilization error converge at the optimal rates $O(h^{k+1})$, $O(h^{k+2})$, and $O(h^{k+1})$, respectively, as expected. The $L^2$-error converges at the optimal rate $O(h^{k+3})$, except for $k = 0$ where the rate is only $O(h^2)$. All these rates are consistent with the analysis presented in the previous section. We also notice that some errors and decay rates are not reported on the finest meshes for the higher polynomial degrees. The reason is that a stagnation of the error at levels around $10^{-8}$ is observed owing to the poor conditioning of the linear system (see below). Notice that stagnation does not affect the stabilization error.

Table 2 reports the condition number of the full and condensed linear systems arising from the $C^0$-HHO methods with $k \in \{2, 3, 4\}$ on the same meshes as above. For fixed $k$, all the condition numbers scale as $O(h^{-4})$, as expected. In addition, the condition number of the condensed linear system is always smaller than that of the full linear system. The ratio is about 4, 6, 9 for $k = 2, 3, 4$, respectively. Furthermore, the condition number of the condensed linear system is above $10^{10}$ on the finest mesh for $k \in \{2, 3\}$ and on the two finest meshes for $k = 4$, thereby causing an error stagnation at about $10^{-8}$ with our current implementation. Finally, we mention that the above condition numbers correspond
<figure><img src="image_6.png" /><figcaption>Figure 1: Convergence of $C^0$-HHO in $H^2$-, $H^1$-, $L^2$-(semi)norms, and stabilization seminorm on triangular meshes.</figcaption></figure>

to nodal basis functions based on the Fekete points. We also tested modal basis functions, but obtained even larger condition numbers (by a factor of about 2, 4, 8 for $k = 2, 3, 4$, respectively).

## 5.2 Comparison with $C^0$-IPDG, HHO, DG, and Morley and HCT FEM

In this section, we compare the computational performance of $C^0$-HHO with the $C^0$-IPDG, HHO, dG, Morley, and HCT methods on triangular meshes. We point out that we implement the jump of the normal gradient for dG and $C^0$-IPDG methods using the weight $n_\partial(k + 1)^2 h_F^{-1}$ for all $F \in \mathcal{F}$, with $n_\partial = 4$. In addition, we implement the jump of the trace term for dG methods using the weight $n_\partial(k + 1)^6 h_F^{-3}$ for all $F \in \mathcal{F}$.

We first compare the $C^0$-HHO method to the dG, HHO (more precisely, the so-called HHO(A) method from [13]), and $C^0$-IPDG methods. To put the $C^0$-HHO and the other methods on a fair comparison basis, we compare the $C^0$-HHO and HHO methods with face polynomial degree $k \ge 0$ to the dG and $C^0$-IPDG methods with cell polynomial degree $\ell = k + 2$, so that all the methods deliver the same decay rates on the $H^2$-error. Moreover, all the methods are assembled by using affine geometric mappings from a reference triangle. Interestingly, we point out that using affine geometric mappings can accelerate the assembling procedure about 4 to 8 times for DG and HHO methods in contrast with the use of a physical basis in each mesh cell, see [3, Sec. 6] for a more detailed discussion.

A comparison of total DoFs, assembling time (including static condensation if applicable), and condition number of the linear system (again after static condensation if applicable)
<table><caption>Table 1: Convergence rates of $C^0$-HHO in $H^2$-, $H^1$-, $L^2$-(semi)norms, and stabilization semi-norm on triangular meshes.</caption><thead><tr><th>H<sup>2</sup> error</th><th colspan="2">k = 0</th><th colspan="2">k = 1</th><th colspan="2">k = 2</th><th colspan="2">k = 3</th></tr><tr><th>#T</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th></tr></thead><tbody><tr><td>32</td><td>137</td><td>—</td><td>281</td><td>—</td><td>457</td><td>—</td><td>665</td><td>—</td></tr><tr><td>128</td><td>497</td><td>0.99</td><td>1041</td><td>1.62</td><td>1713</td><td>2.69</td><td>2513</td><td>3.66</td></tr><tr><td>512</td><td>1889</td><td>1.00</td><td>4001</td><td>1.91</td><td>6625</td><td>2.87</td><td>9761</td><td>3.86</td></tr><tr><td>2048</td><td>7361</td><td>0.99</td><td>15681</td><td>1.96</td><td>26049</td><td>2.94</td><td>38465</td><td>3.93</td></tr><tr><td>8192</td><td>29057</td><td>0.98</td><td>62081</td><td>1.98</td><td>103197</td><td>2.97</td><td>152705</td><td>3.96</td></tr><tr><td>32768</td><td>115457</td><td>0.99</td><td>247041</td><td>1.99</td><td>411393</td><td>2.99</td><td>608513</td><td>—</td></tr></tbody><thead><tr><th>H<sup>1</sup> error</th><th colspan="2">k = 0</th><th colspan="2">k = 1</th><th colspan="2">k = 2</th><th colspan="2">k = 3</th></tr><tr><th>#T</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th></tr></thead><tbody><tr><td>32</td><td>137</td><td>—</td><td>281</td><td>—</td><td>457</td><td>—</td><td>665</td><td>—</td></tr><tr><td>128</td><td>497</td><td>1.82</td><td>1041</td><td>2.71</td><td>1713</td><td>3.63</td><td>2513</td><td>4.66</td></tr><tr><td>512</td><td>1889</td><td>1.93</td><td>4001</td><td>2.89</td><td>6625</td><td>3.87</td><td>9761</td><td>4.85</td></tr><tr><td>2048</td><td>7361</td><td>1.96</td><td>15681</td><td>2.95</td><td>26049</td><td>3.94</td><td>38465</td><td>4.96</td></tr><tr><td>8192</td><td>29057</td><td>1.98</td><td>62081</td><td>2.97</td><td>103197</td><td>3.98</td><td>152705</td><td>—</td></tr><tr><td>32768</td><td>115457</td><td>1.99</td><td>247041</td><td>2.98</td><td>411393</td><td>—</td><td>608513</td><td>—</td></tr></tbody><thead><tr><th>L<sup>2</sup> error</th><th colspan="2">k = 0</th><th colspan="2">k = 1</th><th colspan="2">k = 2</th><th colspan="2">k = 3</th></tr><tr><th>#T</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th></tr></thead><tbody><tr><td>32</td><td>137</td><td>—</td><td>281</td><td>—</td><td>457</td><td>—</td><td>665</td><td>—</td></tr><tr><td>128</td><td>497</td><td>1.76</td><td>1041</td><td>3.60</td><td>1713</td><td>4.62</td><td>2513</td><td>5.58</td></tr><tr><td>512</td><td>1889</td><td>1.90</td><td>4001</td><td>3.82</td><td>6625</td><td>4.82</td><td>9761</td><td>5.92</td></tr><tr><td>2048</td><td>7361</td><td>1.95</td><td>15681</td><td>3.92</td><td>26049</td><td>4.94</td><td>38465</td><td>5.82</td></tr><tr><td>8192</td><td>29057</td><td>1.97</td><td>62081</td><td>3.95</td><td>103197</td><td>—</td><td>152705</td><td>—</td></tr><tr><td>32768</td><td>115457</td><td>1.99</td><td>247041</td><td>—</td><td>411393</td><td>—</td><td>608513</td><td>—</td></tr></tbody><thead><tr><th>Stab error</th><th colspan="2">k = 0</th><th colspan="2">k = 1</th><th colspan="2">k = 2</th><th colspan="2">k = 3</th></tr><tr><th>#T</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th><th># DoFs</th><th>rate</th></tr></thead><tbody><tr><td>32</td><td>137</td><td>—</td><td>281</td><td>—</td><td>457</td><td>—</td><td>665</td><td>—</td></tr><tr><td>128</td><td>497</td><td>0.95</td><td>1041</td><td>1.72</td><td>1713</td><td>2.70</td><td>2513</td><td>3.66</td></tr><tr><td>512</td><td>1889</td><td>0.98</td><td>4001</td><td>1.90</td><td>6625</td><td>2.88</td><td>9761</td><td>3.85</td></tr><tr><td>2048</td><td>7361</td><td>0.99</td><td>15681</td><td>1.95</td><td>26049</td><td>2.95</td><td>38465</td><td>3.93</td></tr><tr><td>8192</td><td>29057</td><td>0.99</td><td>62081</td><td>1.98</td><td>103197</td><td>2.98</td><td>152705</td><td>3.97</td></tr><tr><td>32768</td><td>115457</td><td>0.99</td><td>247041</td><td>1.99</td><td>411393</td><td>2.99</td><td>608513</td><td>3.99</td></tr></tbody></table>

is presented in Table 3 for the four discretization methods. We consider a triangular mesh with 32768 cells and let the polynomial degree vary from $k = 0$ to $k = 3$. Concerning assembling time (the reported values are meaningful up to 5-10%), we observe that for the lowest-order cases ($k \in \{0,1\}$), $C^0$-HHO methods are the most effective, whereas for the higher-order cases ($k \in \{2,3\}$), $C^0$-IPDG methods take slightly less time. Another interesting observation is that $C^0$-HHO methods spend about 20% to 30% less time for assembling the linear system compared to (fully discontinuous) HHO methods and dG methods. Concerning the condition number, the first observation is that the condition number for $C^0$-HHO and (fully discontinuous) HHO methods is comparable although the latter is about 10% to 30% larger. The second observation is that the condition number of $C^0$-HHO methods is (significantly) smaller than the one produced by $C^0$-IPDG and (fully discontinuous) dG methods. One reason for this difference is the elimination of the higher-order bubble functions from the cell unknowns in $C^0$-HHO methods. We also mention that the condition
<table><caption>Table 2: Condition number of condensed and full linear system for $C^0$-HHO on triangular meshes.</caption><thead><tr><th rowspan="2">Cond No</th><th colspan="2">k = 2</th><th colspan="2">k = 3</th><th colspan="2">k = 4</th></tr><tr><th>condensed</th><th>full</th><th>condensed</th><th>full</th><th>condensed</th><th>full</th></tr></thead><tbody><tr><td>#T</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>32</td><td>2.99e+04</td><td>8.74e+04</td><td>8.56e+04</td><td>3.93e+05</td><td>2.02e+05</td><td>1.53e+06</td></tr><tr><td>128</td><td>2.96e+05</td><td>9.33e+05</td><td>8.82e+05</td><td>4.41e+06</td><td>2.13e+06</td><td>1.77e+07</td></tr><tr><td>512</td><td>3.41e+06</td><td>1.14e+07</td><td>1.05e+07</td><td>5.64e+07</td><td>2.59e+07</td><td>2.34e+08</td></tr><tr><td>2048</td><td>4.96e+07</td><td>1.72e+08</td><td>1.55e+08</td><td>8.65e+08</td><td>3.89e+08</td><td>3.62e+09</td></tr><tr><td>8192</td><td>7.42e+08</td><td>2.63e+09</td><td>2.36e+09</td><td>1.34e+10</td><td>5.95e+09</td><td>5.68e+10</td></tr><tr><td>32768</td><td>1.13e+10</td><td>4.09e+10</td><td>3.63e+10</td><td>2.11e+11</td><td>9.23e+10</td><td>8.97e+11</td></tr></tbody></table>

number of the $C^0$-IPDG and dG methods are about 4 times smaller if one chooses $n_\partial = 1$ instead of $n_\partial = 4$, although the former value is below the minimal coercivity threshold predicted by the theory.

<table><caption>Table 3: Comparison of total DoFs, assembling time, and condition number for the $C^0$-HHO, $C^0$-IPDG, HHO and dG methods. The polynomial degree is chosen so that all the methods deliver the same decay rates on the $H^2$-error.</caption><thead><tr><th colspan="4">C<sup>0</sup>-HHO</th><th colspan="4">C<sup>0</sup>-IPDG</th></tr><tr><th>order</th><th>kDoFs</th><th>assembling (s)</th><th>Cond No</th><th>order</th><th>kDoFs</th><th>assembling (s)</th><th>Ratio</th></tr></thead><tbody><tr><td>k = 0</td><td>114</td><td>75</td><td>2.27e+08</td><td>ℓ = 2</td><td>66</td><td>112</td><td>67.4</td></tr><tr><td>k = 1</td><td>212</td><td>197</td><td>2.50e+09</td><td>ℓ = 3</td><td>148</td><td>238</td><td>168.0</td></tr><tr><td>k = 2</td><td>310</td><td>507</td><td>1.13e+10</td><td>ℓ = 4</td><td>263</td><td>498</td><td>369.0</td></tr><tr><td>k = 3</td><td>407</td><td>1091</td><td>3.63e+10</td><td>ℓ = 5</td><td>411</td><td>1061</td><td>597.8</td></tr></tbody><thead><tr><th colspan="4">HHO</th><th colspan="4">dG</th></tr><tr><th>order</th><th>kDoFs</th><th>assembling (s)</th><th>Ratio</th><th>order</th><th>kDoFs</th><th>assembling (s)</th><th>Ratio</th></tr></thead><tbody><tr><td>k = 0</td><td>147</td><td>118</td><td>1.3</td><td>ℓ = 2</td><td>197</td><td>154</td><td>23.1</td></tr><tr><td>k = 1</td><td>244</td><td>306</td><td>1.2</td><td>ℓ = 3</td><td>328</td><td>330</td><td>70.0</td></tr><tr><td>k = 2</td><td>342</td><td>666</td><td>1.1</td><td>ℓ = 4</td><td>492</td><td>718</td><td>125.7</td></tr><tr><td>k = 3</td><td>440</td><td>1284</td><td>1.1</td><td>ℓ = 5</td><td>688</td><td>1438</td><td>226.2</td></tr></tbody></table>

Finally, we compare in Table 4 the $C^0$-HHO methods to the Morley and HCT methods on a triangulation composed of 32768 cells. The first observation is that the Morley element only takes about 65% of the assembling time of $C^0$-HHO methods with $k = 0$, but the condition number of the $C^0$-HHO method is about half of that obtained with the Morley element. The second observation is that the $C^0$-HHO method with $k = 1$ takes about 60% of the assembling time of the HCT element and, at the same time, delivers a condition number which is 5 times smaller than that obtained with the HCT element. One possible reason for the discrepancy in assembly time is that HCT methods cannot be constructed by using affine geometric maps and that these methods lead to larger stencils compared to $C^0$-HHO methods.

## 5.3 Nonsmooth problem with type (II) BC's

We select $f$ and the boundary conditions on $\Omega := (-1, 1)^2$ so that the exact solution to (1) is

$$ u(x, y) := \begin{cases} \sin(\pi y) \left( \frac{x^2}{4\pi} + \frac{1}{8\pi^3} (\cos(2\pi x) - 1) \right) & (-1, 0] \times (-1, 1) \\ \sin(\pi y) \left( \frac{x}{4\pi^2} - \frac{1}{8\pi^3} \sin(2\pi x) \right) & (0, 1) \times (-1, 1). \end{cases} \quad (77) $$
<table><caption>Table 4: Comparison of total DoFs, assembling time, and condition number for the C⁰-HHO, Morley, and HCT methods. The polynomial degree is chosen so that the two methods in the same column deliver the same decay rates on the H²-error. Triangular mesh composed of 32768 cells, 49408 edges, and 16641 vertices.</caption><thead><tr><th>k = 0</th><th>kDoFs</th><th>assembling (s)</th><th>Cond No</th><th>k = 1</th><th>kDoFs</th><th>assembling (s)</th><th>Cond No</th></tr></thead><tbody><tr><td>Morley</td><td>65</td><td>42</td><td>4.79e+08</td><td>HCT</td><td>97</td><td>354</td><td>1.20e+10</td></tr><tr><td>C<sup>0</sup>-HHO</td><td>114</td><td>65</td><td>2.27e+08</td><td>C<sup>0</sup>-HHO</td><td>212</td><td>197</td><td>2.50e+09</td></tr></tbody></table>

It can be checked that the solution satisfies $u \in H^{3.5-\epsilon}$ with $\epsilon > 0$ arbitrarily close to zero. More precisely, the fourth-order derivative $\partial_{xxxx} u$ contains a Dirac measure supported on the line $\{x = 0\}$. In addition, the forcing function satisfies $f \in H^{-0.5-\epsilon}(\Omega)$. Thus, the regularity of $u$ and $f$ are compatible with Assumption 1.1.

<figure><img src="image_7.png" /><figcaption>Figure 2: Convergence of C⁰-HHO and C⁰-IPDG methods in H²-seminorms on triangular meshes with $k = 0, 1, 2, 3$.</figcaption></figure>

The error decay rates depend on whether one considers meshes aligned with the line $\{x = 0\}$ that supports the singularity. Considering first aligned meshes composed of $\{32, 128, 512, 2048, 8192, 32768\}$ triangular cells, Figure 2 shows that, as expected, the C⁰-HHO and C⁰-IPDG methods deliver optimal convergence rates. Considering now non-aligned meshes composed of $\{50, 162, 578, 2178, 8450, 33282, 132098\}$ triangular cells (see Figure 3 for an illustration), we notice from Figure 4 that both C⁰-HHO (with $k = 0$) and C⁰-IPDG (with $\ell = k + 2 = 2$) deliver the optimal decay rate $O(h)$ for the error, whereas the decay rate levels off, as expected, at $O(h^{1.5})$ in the higher-order case. This result is again in agreement with the above theoretical results.

## References

[1] F. BONALDI, D. A. DI PIETRO, G. GEYMONAT, AND F. KRASUCKI, *A hybrid high-order method for Kirchhoff-Love plate bending problems*, ESAIM Math. Model. Numer. Anal., 52 (2018), pp. 393-421.

[2] S. C. BRENNER AND L.-Y. SUNG, *C⁰ interior penalty methods for fourth order elliptic boundary value problems on polygonal domains*, J. Sci. Comput., 22/23 (2005), pp. 83-118.

[3] A. CANGIANI, Z. DONG, E. H. GEORGOULIS, AND P. HOUSTON, *hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes*, SpringerBriefs in Mathematics, Springer, Cham, 2017.
<figure><img src="image_8.png" /><figcaption>Figure 3: Two examples of triangular (nonaligned) meshes with 50 (left) and 162 (right) triangles.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 4: Convergence of $C^0$-HHO and $C^0$-IPDG methods in $H^2$-seminorms on triangular meshes with $k = 0, 1, 2, 3$.</figcaption></figure>

[4] C. CARSTENSEN AND N. NATARAJ, *Lowest-order equivalent nonstandard finite element methods for biharmonic plates*, ESAIM Math. Model. Numer. Anal., 56 (2022), pp. 41–78.

[5] T. CHAUMONT-FRELET, A. ERN, S. LEMAIRE, AND F. VALENTIN, *Bridging the multiscale hybrid-mixed and multiscale hybrid high-order methods*, ESAIM Math. Model. Numer. Anal., 56 (2022), pp. 261–285.

[6] G. CHEN AND M. FENG, *A $C^0$-weak Galerkin finite element method for fourth-order elliptic problems*, Numer. Methods Partial Differential Equations, 32 (2016), pp. 1090–1104.

[7] M. CICUTTIN, A. ERN, AND N. PIGNET, *Hybrid high-order methods—a primer with applications to solid mechanics*, SpringerBriefs in Mathematics, Springer, Cham, [2021].

[8] B. COCKBURN, D. A. DI PIETRO, AND A. ERN, *Bridging the Hybrid High-Order and hybridizable discontinuous Galerkin methods*, ESAIM Math. Model. Numer. Anal., 50 (2016), pp. 635–650.

[9] D. A. DI PIETRO AND J. DRONIOU, *The hybrid high-order method for polytopal meshes*, vol. 19 of MS&A. Modeling, Simulation and Applications, Springer, Cham, [2020]. Design, analysis, and applications.

[10] D. A. DI PIETRO AND A. ERN, *A Hybrid High-Order locking-free method for linear elasticity on general meshes*, Comput. Meth. Appl. Mech. Engrg., 283 (2015), pp. 1–21.
[11] D. A. DI PIETRO, A. ERN, AND S. LEMAIRE, *An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators*, Comput. Meth. Appl. Math., 14 (2014), pp. 461–472.

[12] Z. DONG AND A. ERN, *Hybrid high-order method for singularly perturbed fourth-order problems on curved domains*, ESAIM Math. Model. Numer. Anal., 55 (2021), pp. 3091–3114.

[13] Z. DONG AND A. ERN, *Hybrid high-order and weak Galerkin methods for the biharmonic problem*, SIAM J. Numer. Anal., 60 (2022), pp. 2626–2656.

[14] J. DOUGLAS, JR., T. DUPONT, P. PERCELL, AND R. SCOTT, *A family of $C^1$ finite elements with optimal approximation properties for various Galerkin methods for 2nd and 4th order problems*, RAIRO Anal. Numér., 13 (1979), pp. 227–255.

[15] 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 (2002), pp. 3669–3750.

[16] A. ERN AND J.-L. GUERMOND, *Finite element quasi-interpolation and best approximation*, ESAIM Math. Model. Numer. Anal. (M2AN), 51 (2017), pp. 1367–1385.

[17] ——, *Finite Elements I: Approximation and Interpolation*, vol. 72 of Texts in Applied Mathematics, Springer Nature, Cham, Switzerland, 2021.

[18] ——, *Finite Elements II: Galerkin Approximation, Elliptic and Mixed PDEs*, vol. 73 of Texts in Applied Mathematics, Springer Nature, Cham, Switzerland, 2021.

[19] ——, *Quasi-optimal nonconforming approximation of elliptic PDEs with contrasted coefficients and $H^{1+r}$, $r > 0$, regularity*, Found. Comput. Math., 22 (2022), pp. 1273–1308.

[20] G. FU, *Uniform auxiliary space preconditioning for HDG methods for elliptic operators with a parameter dependent low order term*, SIAM J. Sci. Comput., 43 (2021), pp. A3912–A3937.

[21] V. GIRAULT AND P.-A. RAVIART, *Finite element methods for Navier-Stokes equations*, vol. 5 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 1986. Theory and algorithms.

[22] P. GRISVARD, *Elliptic problems in nonsmooth domains*, vol. 69 of Classics in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2011.

[23] V. A. KOZLOV AND V. G. MAZ' YA, *Singularities in solutions to mathematical physics problems in non-smooth domains*, in Partial differential equations and functional analysis, vol. 22 of Progr. Nonlinear Differential Equations Appl., Birkhäuser Boston, Boston, MA, 1996, pp. 174–206.

[24] S. LEMAIRE, *Bridging the hybrid high-order and virtual element methods*, IMA J. Numer. Anal., 41 (2021), pp. 549–593.

[25] L. MORLEY, *The triangular equilibrium element in the solution of plate bending problems*, Aero. Quart., 19 (1968), pp. 149–169.

[26] L. MU, J. WANG, X. YE, AND S. ZHANG, *A $C^0$-weak Galerkin finite element method for the biharmonic equation*, J. Sci. Comput., 59 (2014), pp. 473–495.

[27] R. L. SCOTT AND S. ZHANG, *Finite element interpolation of nonsmooth functions satisfying boundary conditions*, Math. Comp., 54 (1990), pp. 483–493.

[28] A. VEEGER AND P. ZANOTTI, *Quasi-optimal nonconforming methods for symmetric elliptic problems. III—Discontinuous Galerkin and other interior penalty methods*, SIAM J. Numer. Anal., 56 (2018), pp. 2871–2894.
[29] ——, *Quasi-optimal nonconforming methods for symmetric elliptic problems. II—Overconsistency and classical nonconforming elements*, SIAM J. Numer. Anal., 57 (2019), pp. 266–292.

[30] M. WANG AND J. XU, *The Morley element for fourth order elliptic equations in any dimensions*, Numer. Math., 103 (2006), pp. 155–169.

[31] J. ZHAO, S. CHEN, AND B. ZHANG, *The nonconforming virtual element method for plate bending problems*, Math. Models Methods Appl. Sci., 26 (2016), pp. 1671–1687.