# A generalized Hessian-based error estimator for an IPDG formulation of the biharmonic problem in two dimensions

Théophile Chaumont-Frelet, Joscha Gedicke, Lorenzo Mascotto

► To cite this version:

Théophile Chaumont-Frelet, Joscha Gedicke, Lorenzo Mascotto. A generalized Hessian-based error estimator for an IPDG formulation of the biharmonic problem in two dimensions. 2025. hal-05176686

HAL Id: hal-05176686

https://inria.hal.science/hal-05176686v1

Preprint submitted on 22 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
A generalized Hessian-based error estimator for an IPDG formulation of the biharmonic problem in two dimensions

Théophile Chaumont-Frelet,\* Joscha Gedicke† Lorenzo Mascotto‡

**Abstract**

We consider a two dimensional biharmonic problem and its discretization by means of a symmetric interior penalty discontinuous Galerkin method. Based on the “div-div” complex, a novel split of an error measure based on a generalized Hessian into two terms measuring the conformity and nonconformity of the scheme is proven. This splitting is the departing point for the design of a new reliable and efficient error estimator, which does not involve any DG stabilization. Such an error estimator can be bounded from above by the standard DG residual error estimator. Numerical results assess the theoretical predictions, including the efficiency of the proposed estimator.

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

**Keywords:** discontinuous Galerkin method; biharmonic problem; a posteriori error analysis; generalized Hessian; div-div complex

# 1 Introduction

**State-of-the-art.** Residual-type error estimators for discontinuous Galerkin (DG) discretizations of elliptic problems are by now a fairly well-understood topic. In 2003, three works involving second order elliptic problems were published: Rivière and Wheeler [27] investigated different DG schemes and showed a posteriori error bounds in the $L^2$ norm; Becker, Hansbo, and Larson [1] proved a posteriori error estimates in the energy norm and, in the derivation of the upper bound, they used a Helmholtz decomposition to estimate the nonconforming part of the error, cf. [14]; Karakashian and Pascal discussed analogous estimates employing different tools, i.e., averaging operators mapping piecewise polynomials into globally continuous piecewise polynomials. DG residual error estimators typically contain terms that also appear in residual estimators for conforming methods (elemental residuals, normal jumps of the gradients), plus other jump terms that control the nonconformity of the discrete solution.

Such a structure of the error estimator essentially extends to the case of DG discretizations of fourth order problems. Brenner, Gudi, and Sung [3] exhibited a posteriori error estimators for a lowest order $C^0$ IPDG discretization of the biharmonic problem in Hessian-Hessian formulation; a general order IPDG discretization of the same problem in Laplacian-Laplacian formulation was investigated by Georgoulis, Houston, and Virtanen in [20]. In both references, two dimensional problems were considered; essential tools in the analysis were arithmetic averaging operators, mapping continuous piecewise polynomials and piecewise polynomials, respectively, into $C^1$ finite element spaces. The extension to the three dimensional case and polynomial degree explicit estimates was tackled in [16], where the averaging operator was replaced by an elliptic reconstruction operator into $H^2$ functions, which are not necessarily piecewise polynomials. In all cases, the proposed error estimators contained, in addition to terms that one would expect from $H^2$-conforming discretizations (elemental residuals, jumps along certain directions of broken second and third order derivatives), also terms related to the nonconformity of the scheme (jumps and normal jumps of the gradients).

Such stabilization terms are delicate for two main reasons. First, they involve stability constants depending on the degree $p$ of the scheme; this may influence the $p$ “efficiency” of the error

<footnote>\*Inria Univ. Lille and Laboratoire Paul Painlevé, 59655 Villeneuve-d'Ascq, France, theophile.chaumont@inria.fr</footnote>

<footnote>†Institut für Numerische Simulation, Universität Bonn, 53115 Bonn, gedicke@ins.uni-bonn.de</footnote>

<footnote>‡Dipartimento di Matematica e Applicazioni, Università di Milano-Bicocca, 20125 Milan, Italy; IMATI-CNR, 27100, Pavia, Italy; Fakultät für Mathematik, Universität Wien, 1090 Vienna, Austria, lorenzo.mascotto@unimib.it</footnote>
estimators [16]. Second, the choice of the stabilization parameters affects the proof of the quasi-optimality of adaptive schemes driven by that error estimator. The works [24] and [2] discussed that the stabilization parameter for DG discretizations of second order elliptic problems must be even larger than what is needed for the proof of the well-posedness of the scheme, in order to show the convergence and the quasi-optimality of those adaptive algorithms. Effectively, it was proven in [25] for the second order case and [15] for the fourth order case discretized with a $C^0$-IPDG scheme, that this additional restriction is not really necessary for the proof of the convergence of the scheme; yet, it remains open the question whether that condition is necessary or not for the proof of quasi-optimality.

**Contents and contributions.** We consider the biharmonic problem in two dimensions with clamped boundary conditions and its discretization by means of a symmetric interior penalty method, based on piecewise polynomials of general order larger than or equal to 2. We consider a generalized Hessian $\mathbb{H}_h(u_h)$ of the discrete solution $u_h$ to the method, which consists of the broken Hessian and the fourth order DG lifting of $u_h$ in the sense of display (3.4) below. Given $u$ the solution to the continuous problem, for an error measure of the form

$$ \|\mathbb{D}^2 u - \mathbb{H}_h(u_h)\|_{0,\Omega}, $$

which does not depend on the stabilization parameters of the scheme unless through $u_h$, we construct an error estimator computable from $\mathbb{H}_h(u_h)$ and the data of the problem only, which is reliable and efficient for the above error measure. This estimator does not contain explicitly the stabilization terms that typically appear in standard residual a posteriori error bounds.

**Crucial tools in the analysis.** Amongst others, we pinpoint the following results that are needed to derive the a posteriori error bounds:

* trace operators (see Proposition 2.1) and a Green identity (see display (2.3)) for the space of $H(\text{div div})$ symmetric tensors;

* the exactness of the two dimensional div-div complex (see display (2.4));

* $C^1$ quasi-interpolators mapping DG functions into globally $H^2$ piecewise polynomials, which preserves values at the vertices (see display (5.5));

* several bubble-type inverse inequalities [28] (see Section 5.1);

* the existence of a stable right-inverse of the sym curl operator (see Appendix A).

**Comments on the novel error estimator.** The current setting has some differences compared to that of standard residual error estimators as in [3, 16, 20]. First, we consider an error measure (4.4) and design an error estimator (5.2a), which do not involve stabilization terms explicitly. This may lead in the future to different and possibly easier proofs of convergence of related adaptive algorithms. The new error estimator is reliable and efficient, and can be estimated from above by the standard residual one (see Section 5.5). The derivation of the error measure and estimator follows ideas from [8] with significant further elaboration related to the fourth order nature of the problem.

**List of the main results.** For the reader's convenience we spell out a list of the main results of this work:

* a novel split (4.5) of an error measure involving the generalized Hessian in (3.9) into two terms that measure the conforming and nonconforming parts;

* the two essential properties listed in (5.1) and proven in Section 5.4 that a generalized Hessian should satisfy in order to construct a reliable and efficient error estimator;

* the a posteriori error estimates in Theorem 5.3 for the global error estimator in (5.2);

* the comparison with the standard DG residual error estimator in Section 5.5.
**Outline of the paper.** We introduce some Sobolev spaces, the $H(\operatorname{div} \operatorname{div})$ complex, and the continuous problem in Section 2. The interior penalty method is presented in Section 3, including a rewriting based on generalized Hessians. Section 4 is devoted to derive an error measure, which can be split into a combination of a conforming and a nonconforming error measures. An error estimator and a posteriori error bounds for those error measures are detailed in Section 5. Numerical results are presented in Section 6 and some conclusions are drawn in Section 7. Appendix A is concerned with the construction of a stable right-inverse of the sym **curl** operator.

## 2 The continuous setting

**Differential operators.** We denote the gradient, divergence, vector divergence, Laplacian, bi-laplacian, curl, symmetric curl, and Hessian operators by $\nabla$, div, **div**, $\Delta$, $\Delta^2$, **curl**, sym **curl**, and $\mathbb{D}^2$.

**Lebesgue and basic Sobolev spaces.** In what follows, $D$ is a Lipschitz domain with diameter $h_D$ and boundary $\partial D$. The Lebesgue space of measurable functions with finite $\|v\|_{L^2(D)} = \|v\|_{0,D} := (\int_D |v|^2)^{\frac{1}{2}}$ norm is denoted by $L^2(D)$; this norm is induced by the inner product $(u, v)_{0,D} := \int_D u v$.

Given a positive integer $k$, $H^k(D)$ denotes the space of $L^2(D)$ functions with distributional derivatives $D^\alpha$ of order $k$ in $L^2(D)$. We introduce the seminorms and norms

$$|v|_{k,D} := \left( \sum_{|\alpha|=k} \|D^\alpha v\|_{0,D} \right)^{\frac{1}{2}}, \qquad \|v\|_{k,D} := \left( \sum_{\ell=0}^{k} (h_D^{k-\ell} |v|_{H^k(D)})^2 \right)^{\frac{1}{2}}.$$

Noninteger order Sobolev spaces are constructed via interpolation theory [3]. For vector and tensor Lebesgue and Sobolev spaces, we replace the symbols $H$ and $L$ by $\mathbf{H}$ and $\mathbb{H}$, and $\mathbf{L}$ and $\mathbb{L}$. In the tensor case, additional subscripts $\mathcal{S}$ and $\mathbb{T}$ mean that we consider symmetric and trace free tensors. As a matter of notation, we shall employ Latin letters to denote scalars, boldface Greek (on occasions also Latin) letters to denote vector fields, and blackboard bold Latin letters to denote tensors.

A trace theorem holds true for Sobolev spaces [18, Chapter 3]: there exists a continuous operator mapping $H^s(D)$ to $L^2(\partial D)$ for $s$ larger than $1/2$. For $s$ in $(1/2, 3/2)$, we denote the image of $H^s(D)$ through the trace operator by $H^{s-1/2}(\partial D)$. For positive $s$, $H_0^s(D)$ is the closure of $C_0^\infty(D)$ with respect to the $H^s(D)$ norm. The spaces $H_0^1(D)$ and $H_0^2(D)$ coincide with the subspaces of $H^1(D)$ and $H^2(D)$ of functions with zero trace over $\partial D$, and zero trace and trace of the normal component of the gradient over $\partial D$, respectively.

We define the Sobolev space $H^{-1}(\Omega)$ by duality:

$$H^{-1}(D) := [H_0^1(D)]^*, \qquad \|w\|_{-1,D} := \sup_{v \in H_0^1(D)} \frac{\langle w, v \rangle}{\|v\|_{1,D}}.$$

**Some Sobolev boundary spaces.** For a given simply connected, Lipschitz polygonal domain $\Omega$, we define two spaces on its boundary $\partial\Omega$, which will be useful to define trace operators for a certain Sobolev space in (2.2) below. The tangent and (outward) normal vectors of $\partial\Omega$ are $\tau_\Omega$ and $\nu_\Omega$.

For a given subset $F$ of $\partial\Omega$ with nonzero 1-dimensional measure, let

$$H_{0,0}^{\frac{1}{2}}(F) := \{v \in H^{\frac{1}{2}}(F) \mid (\tilde{v} - v)|_F = 0 \text{ where } \tilde{v} \in H^{\frac{1}{2}}(\partial\Omega), \tilde{v}_{\partial\Omega \setminus F} = 0\}.$$

We introduce a boundary space as in [9, Section 2.3]

$$H_{n,0}^{\frac{1}{2}}(\partial\Omega) := \{\partial_{\nu_{\Omega}} v|_{\partial\Omega} \mid v \in H^2(\Omega) \cap H_0^1(\Omega)\} = \{g \in L^2(\partial\Omega) \mid g|_F \in H_{0,0}^{\frac{1}{2}}(F) \quad \forall \text{ facets } F\}$$

and endow it with the norm

$$\|g\|_{H_{n,0}^{\frac{1}{2}}(\partial\Omega)} := \inf_{v \in H^2(\Omega) \cap H_0^1(\Omega), \partial_{\nu_{\Omega}} v|_{\partial\Omega} = g} \|v\|_{2,\Omega}.$$
We also introduce the space

$$H_{e,0}^{\frac{3}{2}}(\partial\Omega) := \{v|_{\partial\Omega} | v \in H^2(\Omega), \partial_{\nu\Omega}v|_{\partial\Omega} = 0, v(\nu) = 0 \quad \forall \nu \text{ vertices of } \Omega\}$$

and endow it with the norm

$$\|g\|_{H_{e,0}^{\frac{3}{2}}(\partial\Omega)} := \inf_{\substack{v \in H^2(\Omega), \partial_{\nu\Omega} v|_{\partial\Omega} = 0, v|_{\partial\Omega} = g \\ v(\nu) = 0 \quad \forall \nu \text{ vertices of } \Omega}} \|v\|_{2,\Omega}.$$

Finally, we consider the spaces

$$H_n^{-\frac{1}{2}}(\partial\Omega) := [H_{n,0}^{\frac{1}{2}}(\partial\Omega)]', \qquad H_e^{-\frac{3}{2}}(\partial\Omega) := [H_{e,0}^{\frac{3}{2}}(\partial\Omega)]', \qquad (2.1)$$

which we can endow with their dual norms.

**The $H(\text{div div}, \Omega)$ space.** Introduce the lowest order Raviart-Thomas space $RT_0(\Omega)$, i.e., the span of $(1, 0)$, $(0, 1)$, and $(x, y)$, and the spaces

$$\mathbb{H}(\text{div div}, \Omega) := \{A \in L^2(\Omega) \mid \text{div div} A \in L^2(\Omega)\}, \\ \mathbb{H}(\text{div div}^0, \Omega) := \{A \in L^2(\Omega) \mid \text{div div} A = 0\}, \qquad (2.2)$$

which we endow with the norm

$$\|\cdot\|_{H(\text{div div}, \Omega)} := (\text{diam}(\Omega)^{-4} \|\cdot\|_{L^2(\Omega)}^2 + \|\text{div div} \cdot\|_{L^2(\Omega)}^2)^{\frac{1}{2}}.$$

We define $\mathbb{H}_S(\text{div div}, \Omega)$ as the subspace of symmetric tensors in $\mathbb{H}(\text{div div}, \Omega)$. In what follows, $\partial_{\tau\Omega}$ and $\partial_{\nu\Omega}$ denote the tangential and normal derivatives along $\partial\Omega$.

A trace theorem for $\mathbb{H}_S(\text{div div})$ spaces is valid; see [9, Lemmas 2.3 and 2.4].

**Proposition 2.1.** For any $\mathbb{B}$ in $\mathbb{H}_S(\text{div div}, \Omega)$, there exists a positive constant $C$ such that

$$\| \nu_{\Omega}^{T} \mathbb{B} \nu_{\Omega} \|_{H_n^{-\frac{1}{2}}(\partial \Omega)} \leq C \| \mathbb{B} \|_{H(\text{div div}, \Omega)}, \qquad \| \partial_{\tau_{\Omega}} (\tau_{\Omega}^{T} \mathbb{B} \tau_{\Omega}) + \nu_{\Omega}^{T} \text{div} \mathbb{B} \|_{H_e^{-\frac{3}{2}}(\partial \Omega)} \leq C \| \mathbb{B} \|_{H(\text{div div}, \Omega)}.$$

The two trace operators above admit continuous right-inverses.

The proof of Proposition 2.1 is strictly related to the validity of the following Green's identity [9, Lemma 2.1]: for all $v$ in $H^2(\Omega)$ and all $\mathbb{B}$ in the space of symmetric tensors in $C^2(\Omega)$,

$$(\text{div div} \mathbb{B}, v)_{0,0} = (\mathbb{B}, \mathbb{D}^2 v)_{0,0} - \sum_{F \in \mathcal{F}^\Omega} \sum_{\nu \in \mathcal{V}^F} \text{sign}_{F,\nu} (\tau_{\Omega}^T \mathbb{B} \nu_{\Omega})(\nu) v(\nu) \\ - \sum_{F \in \mathcal{F}^\Omega} \left[ (\nu_{\Omega}^T \mathbb{B} \nu_{\Omega}, \partial_{\nu_{\Omega}} v)_{0,F} - (\partial_{\tau_{\Omega}} (\tau_{\Omega}^T \mathbb{B} \nu_{\Omega}) + \nu_{\Omega}^T \text{div} \mathbb{B}, v)_{0,F} \right], \qquad (2.3)$$

where

$$\text{sign}_{F,\nu} := \begin{cases} 1 & \text{if } \nu \text{ is the endpoint of } F \\ -1 & \text{if } \nu \text{ is the start endpoint of } F \end{cases} \quad \forall F \in \mathcal{F}^\Omega, \nu \in \mathcal{V}^F.$$

The last two inner products in (2.3) are to be interpreted as suitable duality pairings for tensors $\mathbb{B}$ in $\mathbb{H}_S(\text{div div}, \Omega)$:

$$\sum_{F \in \mathcal{F}^\Omega} (\nu_{\Omega}^T \mathbb{B} \nu_{\Omega}, \partial_{\nu_{\Omega}} v)_{0,F} = -\frac{1}{2}, \partial_{\Omega} \langle \nu_{\Omega}^T \mathbb{B} \nu_{\Omega}, \partial_{\nu_{\Omega}} v \rangle_{\frac{1}{2}, \partial_{\Omega}}, \\ \sum_{F \in \mathcal{F}^\Omega} (\partial_{\tau_{\Omega}} (\tau_{\Omega}^T \mathbb{B} \nu_{\Omega}) + \nu_{\Omega}^T \text{div} \mathbb{B}, v)_{0,F} = -\frac{3}{2}, \partial_{\Omega} \langle \partial_{\tau_{\Omega}} (\tau_{\Omega}^T \mathbb{B} \nu_{\Omega}) + \nu_{\Omega}^T \text{div} \mathbb{B}, v \rangle_{\frac{3}{2}, \partial_{\Omega}}.$$

The bounds in Proposition 2.1 follow from (2.3), picking first $v$ in $H_{n,0}^{\frac{1}{2}}(\partial\Omega)$ and then in $H_{e,0}^{\frac{3}{2}}(\partial\Omega)$.
**Some finite element spaces.** In what follows, $V_h^{H^1}$, $\Sigma_h$, and $V_h$ denote the standard vector La-
grange finite element space of order $p$, the div **div**-conforming tensor finite element space detailed,
e.g., in [9] and [10, Section 5.4] of order $p$, see also Remark 2.2 below, and the standard DG space
of order $p-2$. In Section 5 below, for $p=3,4$, we shall be using the HCT space [17]. None of this
spaces will be employed in the implementation of the scheme (3.8) below, nor in the computation
of the error estimator (5.2); the space $\Sigma_h$ is included for completeness as a pivot space to construct
the FE sequence in (2.5), but is not used in the forthcoming analysis.

**Remark 2.2.** We recall [9, Section 3.2] that the space $\bar{\Sigma}_h$ consists of piecewise tensors polynomials that are $H(\text{div \, div})$-conforming on the whole domain in the sense that they (i) are continuous at the vertices; (ii) admit matching (polynomial) traces in $H_n^{-\frac{1}{2}}$ and $H_e^{-\frac{3}{2}}$ in the sense of (2.1) across all facets.

**Exact sequences.** We are in a position to introduce the two dimensional “div-div” complex,
see [9, 10] and [11, Section 2.4]:

$$
\mathrm{RT}_0(\Omega) \xrightarrow{\subset} \mathrm{H}^1(\Omega) \xrightarrow{\text{sym curl}} \mathbb{H}_S(\operatorname{div \, div}, \Omega) \xrightarrow{\text{div curl}} \mathbb{L}^2(\Omega) \xrightarrow{0} 0. \quad (2.4)
$$

A finite element counterpart of the sequence in (2.4) is

$$
\mathrm{RT}_0(\Omega) \xrightarrow{\subset} \mathbf{V}_h^{H^1} \xrightarrow{\mathrm{sym \, curl}} \bar{\Sigma}_h \xrightarrow{\mathrm{div \, div}} \mathbf{V}_h \xrightarrow{0} 0. \quad (2.5)
$$

The two sequences (2.4) and (2.5) are exact.

**The model problem.** Given $f$ be in $L^2(\Omega)$, we are interested in finding $u : \Omega \to \mathbb{R}$ such that

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

Introduce

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

Noting that $\Delta^2 = \operatorname{div} \mathbf{div}(\mathbb{D}^2)$, a weak formulation of (2.6) reads

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

Problem (2.7) is well-posed; see, e.g., [4, Section 5.9].

3 The method and a generalized Hessian

**Meshes.** We consider a simplicial mesh $\mathcal{T}_h$. Given $h_T$ and $\rho_T$ the diameter and the maximal inradius of each element $T$ in $\mathcal{T}_h$, we introduce the positive shape-regularity parameter $\gamma_T$

$$
\frac{h_T}{\rho_T} \leq \gamma_T. \tag{3.1}
$$

In what follows, the constants appearing in the inequalities are allowed to depend on $\gamma :=$
$\max_{T \in \mathcal{T}_h} \gamma_T$; on occasions, for given positive constants $a$ and $b$, we shall write $a \lesssim b$ if there
exists a positive constant $C$ only depending on $\gamma$.

With each $\mathcal{T}_h$, we associate its set of total $\mathcal{F}_h$, internal $\mathcal{F}_h^I$, and boundary $\mathcal{F}_h^B$ facets. The diameter of a facet $F$ is $h_F$. For a given mesh with given set of vertices $\{v\}$, $h$ denotes a function defined on $\Omega \setminus \{v\}$ such that it equals $h_T$ and $h_F$ in the interior of each element $T$ and in the interior of each facet $F$, respectively. The set of facets of an element $T$ is $\mathcal{F}^T$. Given $T$ and $F$ in $\mathcal{T}_h$ and $\mathcal{F}_h$, we define the patches

$$
\omega^T := \cup \{\tilde{T} \in \mathcal{T}_h \mid \tilde{T} \cap T \neq \emptyset\}, \quad \omega_F := \cup \{\omega^T \mid T \in \mathcal{T}_h, F \in \mathcal{F}^T\}, \quad \tilde{\omega}_F := \cup \{T \mid F \in \mathcal{F}^T\}. \quad (3.2)
$$
For positive s, we introduce the broken Sobolev space

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

For each facet F, we define the jump $\llbracket \cdot \rrbracket_F$ and average $\llbracket \cdot \rrbracket_F$ operators on $L^2(F)$ as follows: for a boundary facet F and v in $L^2(F)$, we set $\llbracket v \rrbracket_F = \llbracket v \rrbracket_F = v$; for an internal facet F with fixed normal vector $\nu_F$, shared by $T^+$ and $T^-$ with outward unit normal vectors $\nu_{T^+}$ and $\nu_{T^-}$, given v in $H^1(\mathcal{T}_h, \Omega)$ with restrictions $v_{T^+}$ and $v_{T^-}$ over $T^+$ and $T^-$, we set

$$\llbracket v \rrbracket_F := (v_{T^+} \nu_F \cdot \nu_{T^+} + v_{T^-} \nu_F \cdot \nu_{T^-})|_F, \qquad \llbracket v \rrbracket_F = \frac{1}{2}(v_{T^+} + v_{T^-})|_F.$$

We employ the same notation for these operators while applied to vector fields and remove the subscript F when it is clear from the context.

We introduce the following skeletal inner product and norm: for any subset $\tilde{\mathcal{F}}_h$ of $\mathcal{F}_h$,

$$(u, v)_{0, \tilde{\mathcal{F}}_h} = \sum_{F \in \tilde{\mathcal{F}}_h} (u, v)_{0, F}, \qquad \|v\|_{0, \tilde{\mathcal{F}}_h}^2 := (v, v)_{0, \tilde{\mathcal{F}}_h} \qquad \forall v \in L^2(F), \forall F \in \tilde{\mathcal{F}}_h.$$

All differential operators defined piecewise over the mesh elements are denoted with the same symbol of the original operator and a subscript h; for instance, the piecewise Hessian over $\mathcal{T}_h$ is denoted by $\mathbb{D}_h^2$.

**The method.** Let $c_\sigma$ and $c_\tau$ be two stabilization parameters to be chosen sufficiently large; see Lemma 3.2 below. We introduce the IPDG bilinear form

$$\begin{aligned} B_h(u_h, v_h) := & (\mathbb{D}_h^2 u_h, \mathbb{D}_h^2 v_h)_{0, \Omega} + (\llbracket u_h \rrbracket, \{\{\nu_F^T \mathbf{div} \mathbb{D}_h^2 u_h\}\})_{0, \mathcal{F}_h} - (\llbracket \nabla u_h \rrbracket, \{\{\nu_F^T \mathbb{D}_h^2 u_h\}\})_{0, \mathcal{F}_h} \\ & + (\{\{\nu_F^T \mathbf{div} \mathbb{D}_h^2 u_h\}\}, \llbracket v_h \rrbracket)_{0, \mathcal{F}_h} - (\{\{\nu_F^T \mathbb{D}_h^2 u_h\}\}, \llbracket \nabla v_h \rrbracket)_{0, \mathcal{F}_h} \\ & + c_\sigma(h^{-3} \llbracket u_h \rrbracket, \llbracket v_h \rrbracket)_{0, \mathcal{F}_h} + c_\tau(h^{-1} \llbracket \nabla u_h \rrbracket, \llbracket \nabla v_h \rrbracket)_{0, \mathcal{F}_h} \quad \forall u_h, v_h \in V_h. \end{aligned} \quad (3.3)$$

We introduce the lifting operator $\mathbb{L}_h: V_h \oplus V \to [V_h]^{2 \times 2}$ as

$$(\mathbb{L}_h(v), \mathbb{B}_h)_{0, \Omega} := (\llbracket v \rrbracket, \{\{\nu_F^T \mathbf{div} \mathbb{B}_h\}\})_{0, \mathcal{F}_h} - (\llbracket \nabla v \rrbracket, \{\{\nu_F^T \mathbb{B}_h\}\})_{0, \mathcal{F}_h} \quad \forall \mathbb{B}_h \in [V_h]^{2 \times 2}, \quad (3.4)$$

which satisfies the following result; see [19, Lemma 5.1].

**Lemma 3.1.** The lifting operator in (3.4) is well-defined. Moreover, there exists a positive constant $c_{\mathbb{L}_h}$ only depending on $c_\sigma$ and $c_\tau$ in (3.3), the shape-regularity parameter $\gamma$, and p such that

$$\|\mathbb{L}_h(v_h)\|_{0, \Omega} \le c_{\mathbb{L}_h} \left( h^{-\frac{3}{2}} \|\llbracket v_h \rrbracket\|_{0, \mathcal{F}_h} + h^{-\frac{1}{2}} \|\llbracket \nabla v_h \rrbracket\|_{0, \mathcal{F}_h} \right) \quad \forall v_h \in V_h. \quad (3.5)$$

If v belongs to V, then $\mathbb{L}_h(v) = [0]^{d \times d}$.

Using (3.4), the bilinear form in (3.3) can be equivalently rewritten as

$$\begin{aligned} B_h(u_h, v_h) := & (\mathbb{D}_h^2 u_h, \mathbb{D}_h^2 v_h)_{0, \Omega} + (\mathbb{L}_h(u_h), \mathbb{D}_h^2 v_h)_{0, \Omega} + (\mathbb{D}_h^2 u_h, \mathbb{L}_h(v_h))_{0, \Omega} \\ & + c_\sigma(\llbracket u_h \rrbracket, \llbracket v_h \rrbracket)_{0, \mathcal{F}_h} + c_\tau(\llbracket \nabla_h u_h \rrbracket, \llbracket \nabla_h v_h \rrbracket)_{0, \mathcal{F}_h} \quad \forall u_h, v_h \in V_h. \end{aligned} \quad (3.6)$$

The bilinear form in (3.6) is well-defined also on $[H_0^2(\Omega)]^2$ whereas that in (3.3) is not. Introduce the DG norm

$$\|v_h\|_{\text{DG}}^2 := \|\mathbb{D}_h^2 v_h\|_{0, \Omega}^2 + \|h^{-\frac{3}{2}} \llbracket v_h \rrbracket\|_{0, \mathcal{F}_h}^2 + \|h^{-\frac{1}{2}} \llbracket \nabla v_h \rrbracket\|_{0, \mathcal{F}_h}^2. \quad (3.7)$$

The following result is valid; see [19, Lemma 5.2].

**Lemma 3.2.** Let the stabilization parameters $\sigma$ and $c_\tau$ in (3.6) be sufficiently large. The following coercivity and continuity properties hold true: there exist positive constant $c_\nu$ and $c_{\#}$ only depending on $c_\sigma$ and $c_\tau$ in (3.3), the shape-regularity parameter $\gamma$, and p such that

$$B_h(v_h, v_h) \ge c_\nu \|v_h\|_{\text{DG}}^2, \qquad B_h(u_h, v_h) \le c_{\#} \|u_h\|_{\text{DG}} \|v_h\|_{\text{DG}}.$$
We are now in a position to introduce an IPDG method for the discretization of problem (2.7):

$$ \left\{ \begin{array}{l} \text{find } u_h \in V_h \text{ such that} \\ B_h(u_h, v_h) = (f, v_h)_{0,\Omega} \quad \forall v_h \in V_h. \end{array} \right. \qquad (3.8) $$

The well-posedness of method (3.8) follows from Lax-Milgram lemma, Lemma 3.2, Cauchy-Schwarz' inequality on the right-hand side of (3.8), and the Poincaré inequality for piecewise $H^2$ functions in [5].

Given $u_h$ the solution to (3.8) and $L_h$ as in (3.4), we introduce the generalized Hessian

$$ \mathbb{H}_h(u_h) := \mathbb{D}_h^2 u_h + \mathbb{L}_h(u_h) \qquad (3.9) $$

and rewrite, for all $u_h$ and $v_h$ in $V_h$, the DG bilinear form in (3.6) as

$$ B_h(u_h, v_h) = (\mathbb{H}(u_h), \mathbb{H}(v_h))_{0,\Omega} - (\mathbb{L}(u_h), \mathbb{L}(v_h))_{0,\Omega} + c_{\sigma}([\mathbb{L}u_h], [\mathbb{L}v_h])_{0,\mathcal{X}_h} + c_{\tau}([\mathbb{V}h u_h], [\mathbb{V}h v_h])_{0,\mathcal{X}_h}. $$

# 4 A splitting of the generalized Hessian error

Given $s$ in $V$ such that

$$ s := \arg \min_{t \in V} \| \mathbb{H}_h(u_h) - \mathbb{D}^2 t \|_{0,\Omega} \iff (\mathbb{D}^2 s - \mathbb{H}_h(u_h), \mathbb{D}^2 v) = 0 \quad \forall v \in V, \qquad (4.1) $$

Pythagoras' theorem implies

$$ \|\mathbb{D}^2 u - \mathbb{H}_h(u_h)\|_{0,\Omega}^2 = \|\mathbb{D}^2(u-s)\|_{0,\Omega}^2 + \|\mathbb{D}^2 s - \mathbb{H}_h(u_h)\|_{0,\Omega}^2 := T_1^2 + T_2^2. \qquad (4.2) $$

As for the term $T_1$, a duality argument reveals that

$$ T_1 = \sup_{v \in V, |v|_{2,\Omega}=1} (\mathbb{D}^2(u-s), \mathbb{D}^2 v)_{0,\Omega} \stackrel{(4.1)}{=} \sup_{v \in V, |v|_{2,\Omega}=1} (\mathbb{D}^2 u - \mathbb{H}_h(u_h), \mathbb{D}^2 v)_{0,\Omega}. \qquad (4.3) $$

We also rewrite the term $T_2$ as a sup over a suitable space. To this aim, we need a technical result.

**Proposition 4.1.** The tensor $\mathbb{D}^2 s - \mathbb{H}_h(u_h)$ belongs to $\mathbb{H}(\text{div } \mathbf{div}^0, \Omega)$ in (2.2).

*Proof.* The proof boils down to using the Green's identity in (2.3) and showing that all the terms on the right-hand side vanish. The latter task is achieved by taking test functions $v$ with a suitable local supports. □

Using Proposition 4.1 and observing that the space $\mathbb{H}(\text{div } \mathbf{div}^0, \Omega)$ is Hilbert with respect to the $L^2(\Omega)$ norm, a duality argument entails the following rewriting of the term $T_2$ in (4.2):

$$ T_2 = \sup_{\substack{\psi \in H(\text{div } \mathbf{div}^0, \Omega), \\ \|\psi\|_{0,\Omega}=1}} (\mathbb{H}_h(u_h) - \mathbb{D}^2 s, \psi)_{0,\Omega}. $$

Using the “divdiv” complex (2.4), the fact that $\mathbb{D}^2 s$ is symmetric, and an integration by parts, we arrive at

$$ \begin{aligned} T_2 &= \sup_{\substack{\psi \in \mathbb{H}^1(\Omega), \\ \|\text{sym curl } \psi\|_{0,\Omega}=1}} (\mathbb{H}_h(u_h) - \mathbb{D}^2 s, \text{sym curl } \psi)_{0,\Omega} \\ &= \sup_{\substack{\psi \in \mathbb{H}^1(\Omega), \\ \|\text{sym curl } \psi\|_{0,\Omega}=1}} (\mathbb{H}_h(u_h), \text{sym curl } \psi)_{0,\Omega}. \end{aligned} \qquad (4.4) $$

Collecting (4.3) and (4.4), we obtain the crucial identity

$$ \begin{aligned} \|\mathbb{D}^2 u - \mathbb{H}_h(u_h)\|_{0,\Omega}^2 &= \left( \sup_{v \in V, |v|_{2,\Omega}=1} (\mathbb{D}^2 u - \mathbb{H}_h(u_h), \mathbb{D}^2 v)_{0,\Omega} \right)^2 \\ &\quad + \left( \sup_{\substack{\psi \in \mathbb{H}^1(\Omega), \\ \|\text{sym curl } \psi\|_{0,\Omega}=1}} (\mathbb{H}_h(u_h), \text{sym curl } \psi)_{0,\Omega} \right)^2. \end{aligned} \qquad (4.5) $$
# 5 A posteriori bounds for a class of generalized Hessians

The forthcoming analysis applies to all piecewise polynomial tensor fields $\mathbb{H}_h$ satisfying the following properties:

$$ (\mathbb{H}_h, \mathbb{D}^2 w_h)_{0,\Omega} = (f, w_h)_{0,\Omega} \qquad \forall w_h \in V_h \cap H_0^2(\Omega), \qquad (5.1a) $$

$$ (\mathbb{H}_h, \mathrm{sym} \, \mathrm{curl} \, \theta_h)_{0,\Omega} = 0 \qquad \forall \theta_h \in \mathbf{V}_h^{\mathrm{H}^1}. \qquad (5.1b) $$

In what follows, we shall refer to such tensor fields as generalized Hessians. In (5.1b), we recall that $\mathbf{V}_h^{\mathrm{H}^1}$ denotes the standard vector $\mathbf{H}^1$-conforming finite element space. In Section 5.4, we shall be proving that $\mathbb{H}_h(u_h)$ in (3.9) satisfies (5.1).

**Remark 5.1.** In (5.1a), the space $V_h \cap H_0^2(\Omega)$ is the finite element space $V_h^{\mathrm{H}^2}$ as in [26] for $p$ larger than or equal to 5. Assuming that a given mesh $\mathcal{T}_h$ is an HCT refinement of a coarser mesh $\tilde{\mathcal{T}}_h$ (i.e., an element $\tilde{T}$ in $\tilde{\mathcal{T}}_h$ is refined into three elements of $\mathcal{T}_h$ obtained connecting the vertices of $\tilde{T}$ to its centroid), the space $V_h \cap H_0^2(\Omega)$ contains the HCT finite element space of degree $p = 3, 4$ on the mesh $\tilde{\mathcal{T}}_h$.

Aim of this section is to derive upper and lower a posteriori error bounds for the error in (4.5) in terms of the error estimator

$$ \eta^2 := \sum_{T \in \mathcal{T}_h} \eta_T^2, \qquad (5.2a) $$

where

$$ \eta_T^2 := \sum_{j=1}^{5} \eta_{j,T}^2 := h_T^4 \|f - \mathrm{div} \, \mathrm{div} \, \mathbb{H}_h\|_{0,T}^2 \\ + \frac{1}{2} \sum_{F \in \mathcal{F}^T \cap \mathcal{F}_h^1} \left( h_F \|\boldsymbol{\nu}_F^T \llbracket \mathbb{H}_h \rrbracket \boldsymbol{\nu}_F \|_{0,F}^2 + h_F^3 \|\partial_{\boldsymbol{\tau}_F}(\boldsymbol{\tau}_F^T \llbracket \mathbb{H}_h \rrbracket \boldsymbol{\nu}_F) + \boldsymbol{\nu}_F^T \llbracket \mathrm{div} \, \mathbb{H}_h \rrbracket \|_{0,F}^2 \right) \\ + h_T^2 \|\mathrm{curl} \, \mathrm{sym} \, \mathbb{H}_h\|_{0,T}^2 + \sum_{F \in \mathcal{F}^T} h_F \|\boldsymbol{\tau}_F^T \llbracket \mathrm{sym} \, \mathbb{H}_h \rrbracket \|_{0,F}^2. \qquad (5.2b) $$

The local error estimators in (5.2b) do not depend explicitly on the stabilization of the DG method. The first two terms and the last two terms on the right-hand side of (5.2b) are residual-type terms related to the conforming and nonconforming parts of the error, respectively, i.e., the two terms on the right-hand side in (4.5).

**Remark 5.2.** All the powers of $h$ in (5.2b) are positive.

The main result of this section is reported here; its proof is a consequence of the computations in Sections 5.2 and 5.3.

**Theorem 5.3.** Given $u$ and $u_h$ the solutions to (2.7) and (3.8), let $\mathbb{H}_h(u_h)$ be the generalized Hessian as in (3.9), and $\eta$ and $\eta_T$ be the global and local error estimators as in (5.2a) and (5.2b). Furthermore, let $f_h$ denote the piecewise $L^2$ projection onto the space of polynomials of degree $p$. Then, the following reliability and efficiency error bounds are valid: there exist positive constants $C_{\mathrm{rel}}$ and $C_{\mathrm{eff}}$ independent of $h$, but depending on $\gamma$ in (3.1), the domain $\Omega$ through $C_{\mathrm{JG}}$ in (5.21), and the polynomial degree $p$, such that

$$ \|\mathbb{D}^2 u - \mathbb{H}_h(u_h)\|_{0,\Omega} \leq C_{\mathrm{rel}}\eta, \qquad (5.3a) $$

$$ \eta_T \leq C_{\mathrm{eff}} \left( \|\mathbb{D}_h^2 u - \mathbb{H}_h(u_h)\|_{0,\omega^T} + \sum_{\tilde{T} \in \omega^T} h_{\tilde{T}}^2 \|f - f_h\|_{0,\tilde{T}} \right) \quad \forall T \in \mathcal{T}_h. \qquad (5.3b) $$

Bound (5.3b) holds true for any polynomial degree $p$ larger than or equal to 2. Bound (5.3a) is valid assuming $p$ larger than or equal to 5 on any mesh, and $p = 3, 4$ on meshes stemming from HCT refinements as discussed in Remark 5.1.
Numerical examples in Section 6 below indicate that bound (5.3a) remains valid for $p = 2, 3, 4$ on any mesh.

The remainder of this section is organized as follows. After introducing technical tools in Section 5.1, including quasi-interpolation results and polynomial inverse estimates, we derive a posteriori error bounds in Sections 5.2 and 5.3, for any generalized Hessian satisfying the properties detailed in (5.1). We prove that the generalized Hessian in (3.9) satisfies (5.1) in Section 5.4. In order to shorten the notation, we shall write $\mathbb{H}_h$. A comparison with standard DG residual error estimators is given in Section 5.5.

## 5.1 Technical tools

This section is devoted to discussing several technical tools, which will be instrumental in the a posteriori error analysis. Notably, we discuss: the approximation properties of $C^0$ and $C^1$ quasi-interpolants, which are essential in deriving the reliability of the error estimator; extension operators, bubble functions, and polynomial inverse estimates, which are needed in deriving the efficiency of the error estimator.

**$C^0$ quasi-interpolants.** Introduce the vector Clément quasi-interpolant $\theta_I$ of $\theta$ in $\mathbf{H}^1(\Omega)$ as in [12]. We have the standard approximation properties

$$ \|\theta - \theta_I\|_{0,T} + h_T |\theta - \theta_I|_{1,T} \lesssim h_T |\theta|_{1,\omega_T}, \qquad \|\theta - \theta_I\|_{0,F} \lesssim h_T^{\frac{1}{2}} |\theta|_{1,\omega_F}. \quad (5.4) $$

**$C^1$ quasi-interpolants preserving nodal values.** For $p$ larger than or equal to 5, consider the quasi-interpolant $v_I$ in $V_h^{H^2}$ of a function $v$ in $V$ as in [22, Theorem 3.2], which is a modification of that in [21, Section 5] and allows for exact pointwise interpolation at the vertices of the mesh; see also [7, Theorem 2]. The following interpolation estimates can be derived from the results in [22, Theorem 3.2] and standard trace inequalities: for all elements $T$ and all $F$ in $\mathcal{F}^T$,

$$ h_T^{-2} \|v - v_I\|_{0,T} + \sum_{F \in \mathcal{F}^T} \left( h_F^{-\frac{3}{2}} \|v - v_I\|_{0,F} + h_F^{-\frac{1}{2}} \|\nabla(v - v_I)\|_{0,F} \right) \lesssim |v|_{2,\omega_T}. \quad (5.5) $$

Compared to other $C^1$ quasi-interpolants, e.g., that in [21], the one we pick has the property of interpolating functions at the vertices of the mesh [22, Theorem 3.2 (a)]:

$$ (v - v_I)(\nu) = 0 \qquad \forall \nu \text{ vertex of } \mathcal{T}_h. \qquad (5.6) $$

With a similar reasoning, for $p = 3, 4$, and meshes as those discussed in Remark 5.1, one can define an HCT quasi-interpolant satisfying (5.5) and (5.6). In fact, the interpolant in [21] is based on nodal averaging employing auxiliary $L^2$ dual basis on patches/elements. If $v$ belongs to $H^2(\Omega)$, we may replace averages with point values at the pointwise vertex degrees of freedom. This newly defined operator preserves polynomials and is continuous in $H^2(\Omega)$, which is enough to derive the interpolation estimates.

**A facet extension operator.** We shall be considering extension operators $\mathcal{E}_F = \mathcal{E} : L^2(F) \to L^2(\tilde{\omega}_F)$ for all $F$ in $\mathcal{F}_h$. Such operators are defined so as to satisfy

$$ \mathcal{E}v(x) = v(x') \text{ in } T, \quad T \subset \omega_F, \text{ if } x' \in F \text{ is such that } \lambda_j(x') = \lambda_j(x), \quad (5.7) $$

where $\lambda_j$ is one of the two barycentric coordinates of $T$. In particular, this operator acts as the identity on the facet $F$.

**$C^0$ bubble functions and inverse estimates.** In the treatment of the lower bound for the nonconforming part of the error, see Section 5.3, we shall be using two different bubbles. The first one is defined on the element $T$ and is given by

$$ \tilde{b}_T := 27 \prod_{j=1}^3 \lambda_j^T, \qquad (5.8) $$
where the $\lambda_j^T$, $j = 1, 2, 3$, are the barycentric coordinates of $T$. We have [28] the standard inverse estimates

$$ \|q_p\|_{0,T} \lesssim \|\tilde{b}_T^{\frac{1}{2}} q_p\|_{0,T}, \quad |\tilde{b}_T q_p|_{1,T} \lesssim h_T^{-1} \|q_p\|_{0,T} \quad \forall q_p \in P_p(T). \quad (5.9) $$

The second one is defined on the patch $\tilde{\omega}_F$, see (3.2), for a given facet $F$ and is given by

$$ \tilde{b}_{F|T} = 4 \prod_{j=1}^{2} \lambda_j^T \quad \forall T \subset \tilde{\omega}_F, \quad (5.10) $$

where $\lambda_j^T$, $j = 1, 2, 3$, are the barycentric coordinates of $T$ subset of $\tilde{\omega}_F$, $j = 3$ being the index related to the facet $F$. Given $\mathcal{E}$ as in (5.7), we have [28] the inverse estimates

$$ \|q_p\|_{0,F} \lesssim \|\tilde{b}_F^{\frac{1}{2}} q_p\|_{0,T}, \quad |\tilde{b}_F \mathcal{E}(q_p)|_{1,T} + h_T^{-1} \|\tilde{b}_F \mathcal{E}(q_p)\|_{0,T} \lesssim h_F^{\frac{1}{2}} \|q_p\|_{0,F} \quad \forall q_p \in P_p(F). \quad (5.11) $$

**C¹ bubble functions and inverse estimates.** In the treatment of the lower bound for the conforming part of the error, see Section 5.2, we shall be using three different bubbles. The first one is defined on the element $T$ and, for $\tilde{b}_T$ as in (5.8), is given by

$$ b_T := \tilde{b}_T^2. \quad (5.12) $$

In particular, this is a $C^1$ bubble vanishing and with vanishing normal component of the gradient over $\partial T$. We have [28] the standard inverse estimates

$$ \|q_p\|_{0,T} \lesssim \|\tilde{b}_T^{\frac{1}{2}} q_p\|_{0,T}, \quad |b_T q_p|_{2,T} \lesssim h_T^{-2} \|q_p\|_{0,T} \quad \forall q_p \in P_p(T). \quad (5.13) $$

The second one is associated with the facet $F$ and, for $\tilde{b}_F$ as in (5.10), is given by

$$ b_F := \tilde{b}_F^2. \quad (5.14) $$

In particular, this is a $C^1$ bubble vanishing and with vanishing normal component of the gradient over $\partial \tilde{\omega}_F$. We have [28] the standard inverse estimates: given $\mathcal{E}$ as in (5.7),

$$ \begin{aligned} \|q_p\|_{0,F} &\lesssim \|\tilde{b}_F^{\frac{1}{2}} q_p\|_{0,F}, \quad h_F^2 |b_F \mathcal{E}(q_p)|_{2,T} + \|b_F \mathcal{E}(q_p)\|_{0,T} \lesssim h_F^{\frac{1}{2}} \|q_p\|_{0,F}, \\ \|\partial_{\nu_F} (b_F \mathcal{E}(q_p))\|_{0,F} &\lesssim h_F^{-1} \|q_p\|_{0,F} \quad \forall q_p \in P_p(F). \end{aligned} \quad (5.15) $$

The third one is less standard, see [3, Lemma 4.3] and [6, Lemma 6], and is defined on the patch $\tilde{\omega}_F$ as in (3.2). Given $\mathcal{E}$ as in (5.7) and

$$ b_F^* \text{ as in [6, eq. (10)] satisfying } b_F^*|_{\partial \tilde{\omega}_F} = 0, \quad \|q_p\|_{0,F} \lesssim \|(b_F^*)^{\frac{1}{2}} q_p\|_{0,F} \quad \forall q_p \in P_p(F), \quad (5.16) $$

for all $q_p$ in $P_p(F)$, there exists a function $w_F = w_F(q_p)$ in $H_0^2(\tilde{\omega}_F) \subset H_0^2(\Omega)$ piecewise in $P_{p+4}(T)$, $T$ in $\tilde{\omega}_F$, depending on $b_F^*$ such that

$$ \partial_{\nu_F} w_F = b_F^* \mathcal{E}(q_p) = b_F^* q_p \text{ on } F, \quad w_F = 0 \text{ on } F, \quad \|w_F\|_{0,\tilde{\omega}_F} + h_F^2 |w_F|_{2,\tilde{\omega}_F} \lesssim h_F^{\frac{3}{2}} \|q_p\|_{0,F}. \quad (5.17) $$

## 5.2 Upper and lower bounds of the conforming part of the error

**The upper bound.** For all $v$ in $V$ with $|v|_{2,\Omega} = 1$, we have

$$ \begin{aligned} & (\mathbb{D}^2(u-s), \mathbb{D}^2 v)_{0,\Omega} \stackrel{(4.1)}{=} (\mathbb{D}^2 u - \mathbb{H}_h, \mathbb{D}^2 v)_{0,\Omega} \stackrel{(5.1a)}{=} (\mathbb{D}^2 u - \mathbb{H}_h, \mathbb{D}^2 (v-v_I))_{0,\Omega} \\ & \stackrel{(2.3)}{=} \sum_{T \in \mathcal{T}_h} (f - \operatorname{div} \operatorname{div} \mathbb{H}_h, v-v_I)_{0,T} - \sum_{F \in \mathcal{F}_h^1} (\nu_F^T [\mathbb{H}_h], \nabla(v-v_I))_{0,F} + \sum_{F \in \mathcal{F}_h^1} (\nu_F^T [\operatorname{div} \mathbb{H}_h], v-v_I)_{0,F} \\ & \stackrel{(5.6)}{=} \sum_{T \in \mathcal{T}_h} (f - \operatorname{div} \operatorname{div} \mathbb{H}_h, v-v_I)_{0,T} - \sum_{F \in \mathcal{F}_h^1} (\nu_F^T [\mathbb{H}_h], \nu_F, \nu_F^T \nabla(v-v_I))_{0,F} \\ & \qquad + \sum_{F \in \mathcal{F}_h^1} (\partial_{\tau_F} (\tau_F^T [\mathbb{H}_h], \nu_F) + \nu_F^T [\operatorname{div} \mathbb{H}_h], v-v_I)_{0,F}. \end{aligned} $$
Further using Cauchy-Schwarz' inequality (first for integrals, then for sequences) and the approximation properties (5.5) of $\mathcal{C}^1$ quasi-interpolants, and recalling that $|v|_{2,\Omega} = 1$ yield

$$ (\mathbb{D}^2(u-s), \mathbb{D}^2 v)_{0,\Omega} \lesssim \left[ \sum_{T \in \mathcal{T}_h} h_T^4 \|f - \operatorname{div} \mathbf{div} \mathbb{H}_h\|_{0,T}^2 + \sum_{F \in \mathcal{F}_h^1} (h_F \|\boldsymbol{\nu}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F\|_{0,F}^2 + h_F^3 \|\partial_{\tau_F}(\boldsymbol{\tau}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F) + \boldsymbol{\nu}_F^T [\mathbf{div} \mathbb{H}_h]\|_{0,F}^2) \right]^{\frac{1}{2}} \stackrel{(5.2b)}{\le} \eta. $$

**The lower bound.** We begin with estimating the elemental terms. Given $b_T$ as in (5.12), we observe that

$$ \begin{aligned} \|f_h - \operatorname{div} \mathbf{div} \mathbb{H}_h\|_{0,T}^2 &\lesssim (f_h - \operatorname{div} \mathbf{div} \mathbb{H}_h, \underbrace{b_T^2(f_h - \operatorname{div} \mathbf{div} \mathbb{H}_h)}_{:=w_T})_{0,T} \\ &\stackrel{(2.7)}{=} (f - f_h, w_T)_{0,T} + (\mathbb{D}^2 u - \mathbb{H}_h, \mathbb{D}^2 w_T)_{0,T} \\ &\le \|f - f_h\|_{0,T} \|w_T\|_{0,T} + \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,T} \|\mathbb{D}^2 w_T\|_{0,T} \\ &\stackrel{(5.13)}{\lesssim} \|f - f_h\|_{0,T} \|w_T\|_{0,T} + \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,T} h_T^{-2} \|w_T\|_{0,T} \\ &\stackrel{(5.8)}{\le} (\|f - f_h\|_{0,T} + h_T^{-2} \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,T}) \|f_h - \operatorname{div} \mathbf{div} \mathbb{H}_h\|_{0,T}. \end{aligned} $$

For all $T$ in $\mathcal{T}_h$, we deduce

$$ h_T^2 \|f_h - \operatorname{div} \mathbf{div} \mathbb{H}_h\|_{0,T} \lesssim h_T^2 \|f - f_h\|_{0,T} + \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,T}, $$

and, by the triangle inequality,

$$ h_T^2 \|f - \operatorname{div} \mathbf{div} \mathbb{H}_h\|_{0,T} \lesssim h_T^2 \|f - f_h\|_{0,T} + \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,T}. \quad (5.18) $$

Next, we focus on the normal-normal component term appearing in (5.2b). Given $\mathcal{E}$ and $b_F^*$ as in (5.7) and (5.16) (here $\boldsymbol{\nu}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F$ plays the role of $q_p$), we write

$$ \begin{aligned} \|\boldsymbol{\nu}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F\|_{0,F}^{2} &\stackrel{(5.16)}{\lesssim} (b_F^* \mathcal{E}(\boldsymbol{\nu}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F), \boldsymbol{\nu}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F)_{0,F} \stackrel{(5.17)}{=} (\partial_{\boldsymbol{\nu}_F} w_F, \boldsymbol{\nu}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F)_{0,F} \\ &\stackrel{(5.16)}{=} \sum_{j=1}^{2} (\partial_{\boldsymbol{\nu}_{T_j}} w_F, \boldsymbol{\nu}_{T_j}^T [\mathbb{H}_h \boldsymbol{\nu}_{T_j})_{0,\partial T_j} \stackrel{(2.3)}{=} -(\operatorname{div} \mathbf{div} \mathbb{H}_h, w_F)_{0,\tilde{\omega}_F} + (\mathbb{H}_h, \mathbb{D}^2 w_F)_{0,\tilde{\omega}_F} \\ &\stackrel{(2.7)}{=} -(\operatorname{div} \mathbf{div} \mathbb{H}_h - f, w_F)_{0,\tilde{\omega}_F} + (\mathbb{H}_h - \mathbb{D}^2 u, \mathbb{D}^2 w_F)_{0,\tilde{\omega}_F} \\ &\stackrel{(5.17)}{\lesssim} (h_F^2 \|f - \operatorname{div} \mathbf{div} \mathbb{H}_h\|_{0,\tilde{\omega}_F} + \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\tilde{\omega}_F}) h_F^{-\frac{1}{2}} \|\boldsymbol{\nu}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F\|_{0,F}. \end{aligned} $$

We simplify the display above into

$$ h^{\frac{1}{2}} \|\boldsymbol{\nu}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F\|_{0,F} \stackrel{(5.18)}{\lesssim} h_T^2 \|f - f_h\|_{0,\tilde{\omega}_F} + \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\tilde{\omega}_F}. \quad (5.19) $$

Finally, we handle the facet terms involving the normal component of the divergence and the normal-tangential component of the generalized Hessian. Let $F$ be shared by the elements $T_1$ and $T_2$. Given $\mathcal{E}$ and $b_F$ as in (5.7) and (5.14), we have

$$ \begin{aligned} & \|\partial_{\tau_F}(\boldsymbol{\tau}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F) + \boldsymbol{\nu}_F^T [\mathbf{div} \mathbb{H}_h]\|_{0,F}^2 \\ &\stackrel{(5.15)}{\lesssim} (\partial_{\tau_F}(\boldsymbol{\tau}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F) + \boldsymbol{\nu}_F^T [\mathbf{div} \mathbb{H}_h]) \underbrace{, b_F \mathcal{E}(\partial_{\tau_F}(\boldsymbol{\tau}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F) + \boldsymbol{\nu}_F^T [\mathbf{div} \mathbb{H}_h]))}_{:=: w_F} \quad (0,F) \\ &\stackrel{(2.3)}{=} \sum_{j=1}^{2} \left[ - (\operatorname{div} \mathbf{div} \mathbb{H}_h, w_F)_{0,T_j} + (\mathbb{H}_h, \mathbb{D}^2 w_F)_{0,T_j} \right. \\ &\qquad \left. - \sum_{\tilde{F} \in \mathcal{F}_h^T} \sum_{\nu \in \nu_{\tilde{F}}} (\boldsymbol{\tau}_{T_j}^T [\mathbb{H}_h \boldsymbol{\nu}_{T_j}) (\nu) \underbrace{w_F(\nu)}_{=0} - \sum_{\tilde{F} \in \mathcal{F}_h^T} (\boldsymbol{\nu}_{T_j}^T [\mathbb{H}_h \boldsymbol{\nu}_{T_j}, \underbrace{\partial_{\nu_{T_j}} w_F}_{=0 \text{ if } \tilde{F} \neq F}) \right] \quad (0,F) \\ &\stackrel{(2.7)}{=} (f - \operatorname{div} \mathbf{div}_{\mathbb{H}_h} [\mathbb{H}_h, w_F])_{0,\omega_F} + (\mathbb{H}_h - \mathbb{D}^2 u, \mathbb{D}^2 w_F)_{0,\omega_F} + (\boldsymbol{\nu}_F^T [\mathbb{H}_h] \boldsymbol{\nu}_F, \partial_{\boldsymbol{\nu}_F} w_F)_{0,F} =: S_1 + S_2 + S_3. \end{aligned} $$
We estimate the three terms on the right-hand side separately. As for the term $S_1$, we write

$$
\begin{align*}
S_1 &\le \|f - \operatorname{div} \mathbf{div}_h \mathbb{H}_h\|_{0,\omega_F} \|w_F\|_{0,\omega_F} \\
&\stackrel{(5.15)}{\lesssim} \|f - \operatorname{div} \mathbf{div}_h \mathbb{H}_h\|_{0,\omega_F} h_F^{\frac{1}{2}} \|\partial_{\tau_F}(\tau_F^T [\mathbb{H}_h] \nu_F) + \nu_F^T [\mathbf{div} \mathbb{H}_h]\|_{0,F} \\
&\stackrel{(5.18)}{\lesssim} h_F^{-\frac{3}{2}} [h_T^2 \|f - f_h\|_{0,\omega_F} + \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\omega_F}] \|\partial_{\tau_F}(\tau_F^T [\mathbb{H}_h] \nu_F) + \nu_F^T [\mathbf{div} \mathbb{H}_h]\|_{0,F}.
\end{align*}
$$

As for the term $S_2$, we have

$$
S_2 \le \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\omega_F} \|\mathbb{D}^2 w_F\|_{0,\omega_F} \stackrel{(5.15)}{\lesssim} h_F^{-\frac{3}{2}} \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\omega_F} \left\| \partial_{\tau_F}(\tau_F^T [\mathbb{H}_h] \nu_F) + \nu_F^T [\mathbf{div} \mathbb{H}_h] \right\|_{0,F}.
$$

We now focus on the term $S_3$:

$$
S_3 \le \| \nu_F^T [\mathbb{H}_h] \nu_F \|_{0,F} \| \partial_{\nu_F} w_F \|_{0,F} \stackrel{(5.15)}{\lesssim} h_F^{\frac{1}{2}} \| \nu_F^T [\mathbb{H}_h] \nu_F \|_{0,F} h_F^{-\frac{3}{2}} \| \partial_{\tau_F} (\tau_F^T [\mathbb{H}_h] \nu_F) + \nu_F^T [\mathbf{div} \mathbb{H}_h] \|_{0,F}.
$$

Collecting the four displays above yields

$$
\begin{align}
& h_F^{\frac{3}{2}} \|\partial_{\tau_F}(\tau_F^T [\mathbb{H}_h] \nu_F) + \nu_F^T [\mathbf{div} \mathbb{H}_h]\|_{0,F} \nonumber \\
&\lesssim h_T^2 \|f - f_h\|_{0,\omega_F} + \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\omega_F} + h_F^{\frac{1}{2}} \|\nu_F^T [\mathbb{H}_h] \nu_F\|_{0,F} \tag{5.20} \\
&\stackrel{(5.19)}{\lesssim} h_T^2 \|f - f_h\|_{0,\omega_F} + \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\omega_F}. \nonumber
\end{align}
$$

The lower bound for the conforming part of the error is then a consequence of (5.18), (5.19), and (5.20).

5.3 Upper and lower bounds of the nonconforming part of the error

A right-inverse of the sym curl operator. There exists a stable right-inverse of the sym curl operator: given $\mathbb{B}$ in $\mathbb{H}(\operatorname{div} \mathbf{div}^0, \Omega)$, it is possible to construct $\theta$ in $\mathbf{H}^1(\Omega)$ satisfying the following: there exists a positive constant $C_{\text{JG}}$ only depending on the shape of $\Omega$ such that

$$
\operatorname{sym} \operatorname{curl}(\theta) = \mathbb{B}, \qquad \|\theta\|_{1,\Omega} \le C_{\text{JG}} \|\mathbb{B}\|_{0,\Omega}. \tag{5.21}
$$

We refer to Appendix A for a proof.

**The upper bound.** Using (5.21) with $\mathbb{B}$ equal to $\operatorname{sym} \operatorname{curl} \psi$, we pick $\theta$ in $\mathbf{H}^1(\Omega)$ as the $\operatorname{sym} \operatorname{curl}$ right-inverse of $\operatorname{sym} \operatorname{curl} \psi$, and given $\theta_I$ the vector Clément quasi-interpolant of $\theta$ in $\mathbf{H}^1(\Omega)$ as in (5.4), we have

$$
\begin{align*}
(\mathbb{H}_h, \operatorname{sym} \operatorname{curl} \psi)_{0,\Omega} &\stackrel{(5.21)}{=} (\mathbb{H}_h, \operatorname{sym} \operatorname{curl} \theta)_{0,\Omega} \\
&\stackrel{(5.1b)}{=} (\mathbb{H}_h, \operatorname{sym} \operatorname{curl}(\theta - \theta_I))_{0,\Omega} = (\operatorname{sym} \mathbb{H}_h, \operatorname{curl}(\theta - \theta_I))_{0,\Omega} \\
&\stackrel{(5.4)}{\lesssim} \left[ \sum_{T \in \mathcal{T}_h} h_T^2 \|\operatorname{curl} \operatorname{sym} \mathbb{H}_h\|_{0,T}^2 + \sum_{F \in \mathcal{F}_h} h_F \|\tau_F^T [\operatorname{sym} \mathbb{H}_h]\|_{0,F}^2 \right]^{\frac{1}{2}} |\theta|_{1,\Omega} \\
&\stackrel{(5.21)}{\lesssim} \left[ \sum_{T \in \mathcal{T}_h} h_T^2 \|\operatorname{curl} \operatorname{sym} \mathbb{H}_h\|_{0,T}^2 + \sum_{F \in \mathcal{F}_h} h_F \|\tau_F^T [\operatorname{sym} \mathbb{H}_h]\|_{0,F}^2 \right]^{\frac{1}{2}} \underbrace{\|\operatorname{sym} \operatorname{curl} \psi\|_{0,\Omega}}_{=1}.
\end{align*}
$$

**The lower bound.** We begin with estimating the elemental terms. For all $T$ in $\mathcal{T}_h$, given $\tilde{b}_T$ as in (5.8), we have

$$
\begin{align*}
\|\operatorname{curl} \operatorname{sym} \mathbb{H}_h\|_{0,T}^2 &\lesssim (\operatorname{curl} \operatorname{sym} \mathbb{H}_h, \underbrace{\operatorname{curl} \operatorname{sym} \mathbb{H}_h}_{=:w_T})_{0,T} \\
&= (\operatorname{curl}(\operatorname{sym} \mathbb{H}_h - \mathbb{D}^2 u), w_T)_{0,T} \\
&= (\operatorname{sym} \mathbb{H}_h - \mathbb{D}^2 u, \operatorname{curl} w_T)_{0,T} \le \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,T} \|\operatorname{curl} w_T\|_{0,T}.
\end{align*}
$$
By noting that

$$ \|\mathrm{curl} \mathbf{w}_T\|_{0,T} \lesssim |\mathbf{w}_T|_{1,T} \stackrel{(5.9)}{\lesssim} h_T^{-1} \|\mathrm{curl} \mathrm{sym} \mathbb{H}_h\|_{0,T}, $$

we deduce, for all elements $T$,

$$ h_T \|\mathrm{curl} \mathrm{sym} \mathbb{H}_h\|_{0,T} \lesssim \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,T}. \qquad (5.22) $$

We now focus on the facet terms. We focus on internal facets $F$ shared by the elements $T_1$ and $T_2$; the analysis for boundary facets is dealt with analogously. Given $\mathcal{E}$ and $\tilde{b}_F$ as in (5.7) and (5.10), we have

$$ \begin{aligned} \|\boldsymbol{\tau}_F^T [\mathrm{sym} \mathbb{H}_h]\|_{0,F}^2 &\stackrel{(5.11)}{\lesssim} (\boldsymbol{\tau}_F^T [\mathrm{sym} \mathbb{H}_h], \tilde{b}_F \mathcal{E}(\boldsymbol{\tau}_F^T [\mathrm{sym} \mathbb{H}_h]))_{0,F} \\ &= \sum_{j=1}^2 [(\mathrm{curl} \mathrm{sym} \mathbb{H}_h, \mathbf{w}_F)_{0,T_j} + (\mathrm{sym} \mathbb{H}_h, \mathrm{curl} \mathbf{w}_F)_{0,T_j}] \\ &= (\mathrm{curl}_h \mathrm{sym} \mathbb{H}_h, \mathbf{w}_F)_{0,\tilde{\omega}_F} + (\mathrm{sym} \mathbb{H}_h, \mathrm{curl} \mathbf{w}_F)_{0,\tilde{\omega}_F} =: S_1 + S_2. \end{aligned} $$

We estimate the two terms on the right-hand side separately. As for the term $S_1$, we note that

$$ \begin{aligned} S_1 &\lesssim h_F \|\mathrm{curl}_h \mathrm{sym} \mathbb{H}_h\|_{0,\omega_F} h_F^{-1} \|\mathbf{w}_F\|_{0,\omega_F} \\ &\stackrel{(5.11),(5.22)}{\lesssim} \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,T} h_F^{-\frac{1}{2}} \|\mathbf{w}_F\|_{0,F} \le \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,T} h_F^{-\frac{1}{2}} \|\boldsymbol{\tau}_F^T [\mathrm{sym} \mathbb{H}_h]\|_{0,F}. \end{aligned} $$

As for the term $S_2$, using the fact that each column of $\mathbb{D}^2 u$ belongs to $\mathbf{H}(\mathrm{curl}^0, \Omega)$ and an integration by parts, we deduce

$$ \begin{aligned} S_2 &= (\mathrm{sym} \mathbb{H}_h - \mathbb{D}^2 u, \mathrm{curl} \mathbf{w}_F)_{0,\omega_F} \le \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\omega_F} \|\mathrm{curl} \mathbf{w}_F\|_{0,\omega_F} \\ &\stackrel{(5.11)}{\lesssim} \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\omega_F} h_F^{-\frac{1}{2}} \|\mathbf{w}_F\|_{0,F} \le \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\omega_F} h_F^{-\frac{1}{2}} \|\boldsymbol{\tau}_F^T [\mathrm{sym} \mathbb{H}_h]\|_{0,F}. \end{aligned} $$

We combine the displays above and get

$$ h_F^{\frac{1}{2}} \|\boldsymbol{\tau}_F^T [\mathrm{sym} \mathbb{H}_h]\|_{0,F} \lesssim \|\mathbb{D}^2 u - \mathbb{H}_h\|_{0,\omega_F}. \qquad (5.23) $$

The lower bound for the nonconforming part of the error follows from (5.22) and (5.23).

## 5.4 Properties of the generalized Hessian

Here, we prove that the generalized Hessian in (3.9) satisfies (5.1a) and (5.1b). Identity (5.1a) is a consequence of the following computations: for all $w_h$ in $V_h \cap H_0^2(\Omega)$,

$$ (\mathbb{H}_h(u_h), \mathbb{D}^2 w_h)_{0,\Omega} \stackrel{(3.9)}{=} (\mathbb{D}_h^2 u_h + \mathbb{L}_h(u_h), \mathbb{D}^2 w_h)_{0,\Omega} \stackrel{w_h \in H_0^2(\Omega), (3.6)}{=} B_h(u_h, w_h) \stackrel{w_h \in V_h, (3.8)}{=} (f, w_h). $$

The proof of identity (5.1b) is more involved and is intrinsically related to the proof of Green's identity (2.3). We shall be equivalently proving that, cf. (2.5),

$$ (\mathbb{D}_h^2 u_h, \mathbb{A}_h)_{0,\Omega} = -(\mathbb{L}_h(u_h), \mathbb{A}_h)_{0,\Omega}, \qquad \mathbb{A}_h = \mathrm{sym} \, \mathrm{curl} \, \theta_h \quad \forall \theta_h \in \mathbf{V}_h^{H_1}. \qquad (5.24) $$

As for the left-hand side of (5.24), we integrate by parts elementwise twice and get

$$ (\mathbb{D}_h^2 u_h, \mathbb{A}_h)_{0,\Omega} = \sum_{T \in \mathcal{T}_h} [(u_h, \mathrm{div} \, \mathrm{div} \, \mathbb{A}_h)_{0,T} + (\nabla u_h, \mathbb{A}_h \, \boldsymbol{\nu}_T)_{0,\partial T} - (u_h, \boldsymbol{\nu}_T^T (\mathrm{div} \, \mathbb{A}_h))_{0,\partial T}]. $$

The bulk terms vanish since $\mathrm{div} \, \mathrm{div} \, \mathbb{A}_h = 0$ over $\Omega$. Consider the decomposition into tangential and normal parts

$$ \mathbb{A}_h \boldsymbol{\nu}_T = (\boldsymbol{\nu}_T^T \mathbb{A}_h \boldsymbol{\nu}_T) \boldsymbol{\nu}_T + (\boldsymbol{\tau}_T^T \mathbb{A}_h \boldsymbol{\nu}_T) \boldsymbol{\tau}_T \qquad (5.25) $$
This allows us to write

$$ (\mathbb{D}_h^2 u_h, \mathbb{A}_h)_{0,\Omega} = \sum_{T \in \mathcal{T}_h} [(\nu_T^T \nabla u_h, \nu_T^T \mathbb{A}_h \nu_T)_{0,\partial T} + (\tau_T^T \nabla u_h, \tau_T^T \mathbb{A}_h \nu_T)_{0,\partial T} - (u_h, \nu_T^T (\mathrm{div} \, \mathbb{A}_h))_{0,\partial T}]. $$

As for the second terms, we observe that

$$ (\tau_T^T \nabla u_h, \tau_T^T \mathbb{A}_h \nu_T)_{0,\partial T} = (\partial_{\tau\Omega} u_h, \tau_T^T \mathbb{A}_h \nu_T)_{0,\partial T}. $$

With the notation as in (2.3), a tangential integration by parts on $\partial T$ gives

$$ (\tau_T^T \nabla u_h, \tau_T^T \mathbb{A}_h \nu_T)_{0,\partial T} = -(u_h, \partial_{\tau_T}(\tau_T^T \mathbb{A}_h \nu_T))_{0,\partial T} + \sum_{F \in \mathcal{F}_h} \sum_{\nu \in \mathcal{V}^F} \mathrm{sign}_{F,\nu} u_h(\nu)(\tau_T^T \mathbb{A}_h \nu_T)(\nu). $$

Overall, we get

$$ (\mathbb{D}_h^2 u_h, \mathbb{A}_h)_{0,\Omega} = \sum_{T \in \mathcal{T}_h} [(\nu_T^T \nabla u_h, \nu_T^T \mathbb{A}_h \nu_T)_{0,\partial T} - (u_h, \nu_T^T \mathrm{div} \, \mathbb{A}_h + \partial_{\tau_T}(\tau_T^T \mathbb{A}_h \nu_T))_{0,\partial T}] \tag{5.26} \\ + \sum_{F \in \mathcal{F}_h} \sum_{\nu \in \mathcal{V}^F} \mathrm{sign}_{F,\nu} u_h(\nu)(\tau_F^T \mathbb{A}_h \nu_F)(\nu). $$

We focus on the term on the right-hand side of (5.24). The definition of the lifting operator in (3.4) and the splitting in (5.25) imply

$$ \begin{aligned} & (\mathbb{L}_h(u_h), \mathbb{A}_h)_{0,\Omega} \\ &= (\llbracket u_h \rrbracket, \{\{\nu_F^T \mathrm{div} \, \mathbb{A}_h\}\}\}_{0,\mathcal{F}_h} - (\llbracket \nabla u_h \rrbracket, \{\{\mathbb{A}_h \nu_F\}\}\}_{0,\mathcal{F}_h} \\ &= (\llbracket u_h \rrbracket, \{\{\nu_F^T \mathrm{div} \, \mathbb{A}_h\}\}\}_{0,\mathcal{F}_h} - (\llbracket \nu_F^T \nabla u_h \rrbracket, \{\{\nu_F^T \mathbb{A}_h \nu_F\}\}\}_{0,\mathcal{F}_h} - (\llbracket \tau_F^T \nabla u_h \rrbracket, \{\{\tau_F^T \mathbb{A}_h \nu_F\}\}\}_{0,\mathcal{F}_h}. \end{aligned} $$

The continuity of $\nu_F^T \mathbb{A}_h \nu_F$ through facets, cf. the $H(\mathrm{div} \, \mathrm{div})$ conformity of $\mathbb{A}_h$ discussed in Remark 2.2, implies

$$ (\mathbb{L}_h(u_h), \mathbb{A}_h)_{0,\Omega} = (\llbracket u_h \rrbracket, \{\{\nu_F^T \mathrm{div} \, \mathbb{A}_h\}\}\}_{0,\mathcal{F}_h} - (\llbracket \nu_F^T \nabla u_h \rrbracket, \nu_F^T \mathbb{A}_h \nu_F)_{0,\mathcal{F}_h} \\ - (\llbracket \tau_F^T \nabla u_h \rrbracket, \{\{\tau_F^T \mathbb{A}_h \nu_F\}\}\)_{0,\mathcal{F}_h}. $$

We rewrite the last term on the right-hand side, also using an integration by parts on the boundaries of the elements and standard DG manipulations, as

$$ \begin{aligned} & - (\llbracket \tau_F^T \nabla u_h \rrbracket, \{\{\tau_F^T \mathbb{A}_h \nu_F\}\}\)_{0,\mathcal{F}_h} = - (\llbracket \partial_{\tau_F} u_h \rrbracket, \{\{\tau_F^T \mathbb{A}_h \nu_F\}\}\)_{0,\mathcal{F}_h} \\ &= (\llbracket u_h \rrbracket, \{\{\partial_{\tau_F}(\tau_F^T \mathbb{A}_h \nu_F)\}\}\)_{0,\mathcal{F}_h} - \sum_{F \in \mathcal{F}_h} \sum_{\nu \in \mathcal{V}^F} \mathrm{sign}_{F,\nu} u_h(\nu)(\tau_F^T \mathbb{A}_h \nu_F)(\nu). \end{aligned} $$

Combining the above displays yields

$$ (\mathbb{L}_h(u_h), \mathbb{A}_h)_{0,\Omega} = (\llbracket u_h \rrbracket, \{\{\nu_F^T \mathrm{div} \, \mathbb{A}_h + \partial_{\tau_F}(\tau_F^T \mathbb{A}_h \nu_F)\}\}_{0,\mathcal{F}_h} \\ - (\llbracket \nu_F^T \nabla u_h \rrbracket, \nu_F^T \mathbb{A}_h \nu_F)_{0,\mathcal{F}_h} - \sum_{F \in \mathcal{F}_h} \sum_{\nu \in \mathcal{V}^F} \mathrm{sign}_{F,\nu} u_h(\nu)(\tau_F^T \mathbb{A}_h \nu_F)(\nu). $$

The continuity of $\nu_F^T \mathrm{div} \, \mathbb{A}_h + \partial_{\tau_F}(\tau_F^T \mathbb{A}_h \nu_F)$ through facets, cf. the $H(\mathrm{div} \, \mathrm{div})$ conformity of $\mathbb{A}_h$ discussed in Remark 2.2, implies

$$ -(\mathbb{L}_h(u_h), \mathbb{A}_h)_{0,\Omega} = -(\llbracket u_h \rrbracket, \nu_F^T \mathrm{div} \, \mathbb{A}_h + \partial_{\tau_F}(\tau_F^T \mathbb{A}_h \nu_F))_{0,\mathcal{F}_h} \\ + (\llbracket \nu_F^T \nabla u_h \rrbracket, \nu_F^T \mathbb{A}_h \nu_F)_{0,\mathcal{F}_h} + \sum_{F \in \mathcal{F}_h} \sum_{\nu \in \mathcal{V}^F} \mathrm{sign}_{F,\nu}(\tau_F^T \mathbb{A}_h \nu_F)(\nu) u_h(\nu). \tag{5.27} $$

Up to rewriting the sum over the facets as sums over the boundary of elements, the right-hand sides in (5.26) and (5.27) coincide. Identity (5.24) and thus also (5.1b) follow.
## 5.5 Comparison with the standard residual error estimator

The aim of this section is to show that the reliable and efficient error estimator in (5.2a) based on $H_h(u_h)$ can be estimated from above by standard DG residual error estimators available in the literature.

To see this, we recall the general order error estimator in [16, Definition 4.1] ([20, eq. (4.1)] is concerned with for biharmonic problem in Laplace-Laplace formulation)<fnref n="1" />:

$$ \mathcal{J}_T^2 := \sum_{\ell=1}^6 \mathcal{J}_{\ell,T}^2 := h_T^4 \|f - \operatorname{div} \operatorname{div} \mathbb{D}^2 u_h\|_{0,T}^2 + \sum_{F \in \mathcal{F}^T} \left( h_F^{-3} \|[[u_h]]\|_{0,F}^2 + h_F^{-1} \|[[\nabla u_h]]\|_{0,F}^2 \right. \\ \left. \qquad + h_F \|[[\mathbb{D}_h^2 u_h]] \nu_F\|_{0,F}^2 + h_F \|[[\mathbb{D}_h^2 u_h]] \tau_F\|_{0,F}^2 + h_F^3 \|[[\operatorname{div}_h \mathbb{D}_h^2 u_h]] \cdot \nu_F\|_{0,F}^2 \right). \quad (5.28a) $$

For future convenience, we also define the corresponding global error estimator

$$ \mathcal{J}^2 := \sum_{T \in \mathcal{T}_h} \mathcal{J}_T^2. \quad (5.28b) $$

We want to prove that $\eta_T \lesssim \mathcal{J}_T$; to this aim, we estimate each $\eta_{j,T}$, $j = 1, \dots, 5$, in (5.2b) in terms of combinations of $\mathcal{J}_{\ell,T}$, $\ell = 1, \dots, 6$, in (5.28a). In all cases, we split the generalized Hessian $H_h(u_h)$ into the standard Hessian $\mathbb{D}_h^2 u_h$ and the lifting $\mathbb{L}_h(u_h)$. As for the latter, using inverse estimates and the stability result (3.4) bounds all the contributions involving $\mathbb{L}_h(u_h)$ in terms of $\mathcal{J}_{2,T}$ and $\mathcal{J}_{3,T}$. Therefore, we only focus on the remainder parts involving the standard broken Hessian.

The resulting term stemming from $\eta_{1,T}$ boils down to $\mathcal{J}_{1,T}$. On the other hand, as for the term resulting from $\eta_{2,T}$, we have

$$ \sum_{F \in \mathcal{F}^T} h_F \| \nu_F^T [[\mathbb{D}_h^2 u_h]] \nu_F \|_{0,F}^2 \le \mathcal{J}_{4,T}^2. $$

As for the term resulting from $\eta_{3,T}$, we write

$$ \begin{aligned} & \sum_{F \in \mathcal{F}^T} h_F^3 \| \partial_{\tau_F} (\tau_F^T [[\mathbb{D}_h^2 u_h]] \nu_F) + \nu_F^T [[\operatorname{div} \mathbb{D}_h^2 u_h]] \|_{0,F}^2 \\ & \lesssim \sum_{F \in \mathcal{F}^T} h_F^3 ( \| \partial_{\tau_F} (\tau_F^T [[\mathbb{D}_h^2 u_h]] \nu_F) \|_{0,F}^2 + \| \nu_F^T [[\operatorname{div} \mathbb{D}_h^2 u_h]] \|_{0,F}^2 ) \\ & \lesssim \sum_{F \in \mathcal{F}^T} h_F \| [[\mathbb{D}_h^2 u_h]] \nu_F \|_{0,F}^2 + \sum_{F \in \mathcal{F}^T} h_F^3 \| \nu_F^T [[\operatorname{div} \mathbb{D}_h^2 u_h]] \|_{0,F}^2 = \mathcal{J}_{4,T}^2 + \mathcal{J}_{6,T}^2. \end{aligned} $$

The term resulting from $\eta_{4,T}$ vanishes as $\operatorname{curl} \mathbb{D}_h^2 u_h$ is piecewise zero. As for the term resulting from $\eta_{5,T}$, we write

$$ \sum_{F \in \mathcal{F}^T} h_F \| \tau_F^T [[\mathbb{D}_h^2 u_h]] \|_{0,F}^2 = \mathcal{J}_{5,T}^2. $$

The considerations above imply that the local error estimator (5.2b) can be estimated from above by the standard DG residual error estimator. In Section 6 below, we shall assess the practical performance of the two error estimators.

# 6 Numerical results

**Test cases.** The first test case, cf. with [23, pp. 107–108], is set on the L-shaped domain $\Omega_1 := (-1, 1)^2 \setminus ([0, 1) \times (-1, 0])$ with exact solution

$$ u_1(r, \theta) = r^{1+z} g(\theta), \quad (6.1) $$

<footnote n="1">¹The first and sixth terms in display (5.28a) are written equivalently in [16] as (up to the scaling factors) $\|f - \Delta^2 u_h\|_{0,T}$ since for smooth functions $\Delta^2 \cdot = \operatorname{div} \operatorname{div} \mathbb{D}^2 \cdot$ and $\|[[\nabla_h \Delta_h^2 u_h]] \cdot \nu_F\|_{0,F}$ since $\nabla_h \Delta_h^2 \cdot = \operatorname{div}_h \mathbb{D}_h^2$.</footnote>
where $(r, \theta)$ are the polar coordinates at the re-entrant corner $(0,0)$, $z = 0.544483736782464$, and, for $\omega = 3\pi/2$,

$$g(\theta) := \left[ \frac{1}{1-z} \sin((z-1)\omega) - \frac{1}{z+1} \sin((z+1)\omega) \right] \cos((z-1)\theta) - \cos((z+1)\theta) - \left[ \frac{1}{1-z} \sin((z-1)\theta) - \frac{1}{z+1} \sin((z+1)\theta) \right] \cos((z-1)\omega) - \cos((z+1)\omega).$$

The function $u_1$ has inhomogeneous clamped boundary conditions on part of $\partial\Omega_1$ (it has effectively homogeneous boundary conditions on the two edges at the re-entrant corner); we refer to the next paragraph for further discussions on this aspect.

The second test case is set on the unit square domain $\Omega_2 = (0, 1)^2$ with exact solution

$$u_2(x, y) = \sin^2(\pi x) \sin^2(\pi y). \qquad (6.2)$$

**Handling inhomogeneous boundary conditions.** The analysis discussed so far considered the biharmonic problem with homogeneous clamped boundary conditions. The inhomogeneous case follows closely the presentation of the previous sections. For the sake of completeness, we sketch here some details on the inhomogeneous boundary conditions

$$u = g_1, \qquad \partial_{\nu\Omega} u = g_2 \qquad \text{on } \partial\Omega.$$

Formal computations reveal that

$$\begin{aligned} \tau_{\Omega}^{T} \mathbb{D}^{2} u_{|\partial\Omega} &= \partial_{\tau_{\Omega}} (\nabla u_{|\partial\Omega}) = \partial_{\tau_{\Omega}} ((\partial_{\tau_{\Omega}} g_{1}) \tau_{\Omega} + g_{2} \nu_{\Omega}) = (\partial_{\tau_{\Omega}^{2}}^{2} g_{1}) \tau_{\Omega} + (\partial_{\tau_{\Omega}} g_{2}) \nu_{\Omega}, \\ \nabla u_{|\partial\Omega} &= (\nu_{\Omega}^{T} \nabla u_{|\partial\Omega}) \nu_{F} + (\tau_{\Omega}^{T} \nabla u_{|\partial\Omega}) \tau_{\Omega} = g_{2} \nu_{\Omega} + (\partial_{\tau_{\Omega}} g_{1}) \tau_{\Omega}. \end{aligned} \qquad (6.3)$$

For the definition of the lifting operators (3.4), the jump operators on boundary facets modify into

$$[\mathbb{u}_h]_F = \mathbb{u}_h|_F - g_1, \qquad [\nabla \mathbb{u}_h]_F \stackrel{(6.3)}{=} \nabla \mathbb{u}_h|_F - [g_2 \nu_{\Omega} + (\partial_{\tau_{\Omega}} g_1) \tau_{\Omega}].$$

Identity (5.1b) holds true only up to some terms involving the boundary data that eventually modifies the local error estimator in (5.2b), and notably the term $\eta_{5,T}$, which on boundary facets $F$ in $\mathcal{E}^T$ now reads

$$\tau_{\Omega}^{T} [\mathrm{sym} \, \mathbb{H}_h(\mathbb{u}_h)]_F = (\tau_{\Omega}^{T} \mathrm{sym} \, \mathbb{H}_h(\mathbb{u}_h) - (\partial_{\tau_{\Omega}^2}^2 g_1) \tau_{\Omega} - (\partial_{\tau_{\Omega}} g_2) \nu_{\Omega})|_F. \qquad (6.4)$$

The reason for this resides in the fact that identity (4.4) is not true any more. In particular, we have

$$\begin{aligned} (\mathbb{D}^2 u, \mathrm{sym} \, \mathrm{curl} \, \theta)_{0,\Omega} &= (\mathbb{D}^2 u, \mathrm{curl} \, \theta)_{0,\Omega} = (\nabla u, \mathrm{curl} \, \theta \cdot \nu_{\Omega})_{0,\partial\Omega} = (\nabla u, \nabla \theta \cdot \tau_{\Omega})_{0,\partial\Omega} \\ &\stackrel{(6.3)}{=} ((\partial_{\tau_{\Omega}} g_1) \tau_{\Omega} + g_2 \nu_{\Omega}, \partial_{\tau_{\Omega}} \theta)_{0,\partial\Omega} = -((\partial_{\tau_{\Omega}^2}^2 g_1) \tau_{\Omega}, \theta)_{0,\partial\Omega} - ((\partial_{\tau_{\Omega}} g_2) \nu_{\Omega}, \theta)_{0,\partial\Omega}, \end{aligned}$$

which is different from zero in general (while it was zero for homogeneous boundary conditions). In fact, the discrete counterpart of the above identities holds, whence (5.1b) is not valid in general, but should involve the boundary data on the right-hand side with an arbitrary test function $\theta_h$. This reasoning eventually combines with the estimates in Section 5.3 and yields (6.4).

**Outline of the remainder of the section.** In Section 6.1, we compare the performance of the error measure on the left-hand side of (4.5) with the standard DG error under uniform refinements and for different orders. We investigate the practical efficiency of the error estimators in (5.2a) and the standard DG-residual error estimator under uniform h- and p-refinements in Section 6.2. Section 6.3 is devoted to introduce an adaptive algorithm driven by the error estimator in (5.2a) and investigate its performance. In (3.6), we shall always consider

$$c_{\sigma} = 3p^{6}, \qquad c_{\tau} = 9p^{2}.$$

Moreover, in the local estimators (5.2) and (5.28), we replace $h^{\lambda}$, $\lambda$ in $\mathbb{Z}$, with $h^{\lambda}/p^{\lambda}$; this poses no modifications in the analysis and gives a more reasonable scaling in terms of the degree of the scheme.
## 6.1 Comparison of two error measures

Given $u$ and $u_h$ the solutions to (2.7) and (3.8), $\mathbb{H}_h(u_h)$ as in (3.9), and $\|\cdot\|_{\text{DG}}$ as in (3.7), we are interested here in comparing the following two error measures:

$$ \|\mathbb{D}^2 u - \mathbb{H}_h(u_h)\|_{0,\Omega}, \qquad \|u - u_h\|_{\text{DG}}. \qquad (6.5) $$

We pick the exact solution $u_2$ in (6.2), $p = 2$ and $p = 5$, sequences of uniformly refined meshes, and report the results in Tables 1 and 2.

<table><caption>**Table 1:** Error measures in (6.5). Exact solution $u_2$ in (6.2); $p = 2$; sequences of uniformly refined meshes.</caption><thead><tr><th>level</th><th>||u - u<sub>h</sub>||<sub>DG</sub></th><th>order</th><th>||D<sup>2</sup> u - H<sub>h</sub>(u<sub>h</sub>)||<sub>0,Ω</sub></th><th>order</th></tr></thead><tbody><tr><td>1</td><td>1.5095 × 10<sup>1</sup></td><td>—</td><td>9.4596 × 10<sup>0</sup></td><td>—</td></tr><tr><td>2</td><td>1.4193 × 10<sup>1</sup></td><td>0.09</td><td>1.1561 × 10<sup>1</sup></td><td>-0.29</td></tr><tr><td>3</td><td>1.1255 × 10<sup>1</sup></td><td>0.33</td><td>1.0065 × 10<sup>1</sup></td><td>0.20</td></tr><tr><td>4</td><td>4.7746 × 10<sup>0</sup></td><td>1.24</td><td>5.3481 × 10<sup>0</sup></td><td>0.91</td></tr><tr><td>5</td><td>2.4031 × 10<sup>0</sup></td><td>0.99</td><td>2.7127 × 10<sup>0</sup></td><td>0.98</td></tr><tr><td>6</td><td>1.2006 × 10<sup>0</sup></td><td>1.00</td><td>1.3597 × 10<sup>0</sup></td><td>1.00</td></tr><tr><td>7</td><td>5.9949 × 10<sup>-1</sup></td><td>1.00</td><td>6.8016 × 10<sup>-1</sup></td><td>1.00</td></tr><tr><td>8</td><td>2.9948 × 10<sup>-1</sup></td><td>1.00</td><td>3.4011 × 10<sup>-1</sup></td><td>1.00</td></tr><tr><td>9</td><td>1.4966 × 10<sup>-1</sup></td><td>1.00</td><td>1.7006 × 10<sup>-1</sup></td><td>1.00</td></tr><tr><td>10</td><td>7.4812 × 10<sup>-2</sup></td><td>1.00</td><td>8.5029 × 10<sup>-2</sup></td><td>1.00</td></tr><tr><td>11</td><td>3.7401 × 10<sup>-2</sup></td><td>1.00</td><td>4.2514 × 10<sup>-2</sup></td><td>1.00</td></tr></tbody></table>

<table><caption>**Table 2:** Error measures in (6.5). Exact solution $u_2$ in (6.2); $p = 5$; sequences of uniformly refined meshes.</caption><thead><tr><th>level</th><th>||u - u<sub>h</sub>||<sub>DG</sub></th><th>order</th><th>||D<sup>2</sup> u - H<sub>h</sub>(u<sub>h</sub>)||<sub>0,Ω</sub></th><th>order</th></tr></thead><tbody><tr><td>1</td><td>9.1635 × 10<sup>0</sup></td><td>—</td><td>8.7192 × 10<sup>0</sup></td><td>—</td></tr><tr><td>2</td><td>1.9076 × 10<sup>0</sup></td><td>2.26</td><td>1.7656 × 10<sup>0</sup></td><td>2.30</td></tr><tr><td>3</td><td>2.8588 × 10<sup>-2</sup></td><td>6.06</td><td>2.5638 × 10<sup>-2</sup></td><td>6.11</td></tr><tr><td>4</td><td>7.0726 × 10<sup>-3</sup></td><td>2.02</td><td>6.6088 × 10<sup>-3</sup></td><td>1.96</td></tr><tr><td>5</td><td>4.4821 × 10<sup>-4</sup></td><td>3.98</td><td>4.1799 × 10<sup>-4</sup></td><td>3.98</td></tr><tr><td>6</td><td>2.8049 × 10<sup>-5</sup></td><td>4.00</td><td>2.6155 × 10<sup>-5</sup></td><td>4.00</td></tr><tr><td>7</td><td>1.7624 × 10<sup>-6</sup></td><td>3.99</td><td>1.6442 × 10<sup>-6</sup></td><td>3.99</td></tr></tbody></table>

The two error measures in (6.5) converge with the expected order and deliver approximately the same accuracy.

## 6.2 Efficiency under uniform refinements

We assess the performance of the effectivity indices

$$ \mathfrak{F}_{\mathbb{H}_h} := \frac{\eta}{\|\mathbb{D}^2 u - \mathbb{H}_h(u_h)\|_{0,\Omega}}, \qquad \mathfrak{F}_{\text{DG}} := \frac{\mathfrak{J}}{\|u - u_h\|_{\text{DG}}}, \qquad (6.6) $$

where $\eta$ and $\mathfrak{J}$ are given in (5.2a) and (5.28b). First, we consider uniform mesh refinements, pick the exact solution $u_1$ in (6.1), $p = 2, 3, 4, 5$, and $p = 6, 7, 8, 9$ and report the results in Figures 1 and 2.

Next, for the test case with exact solution $u_1$ in (6.1) and a fixed structured mesh of 6 triangles, we investigate the performance of the effectivity indices in (6.6) under $p$-refinements in Figure 3 for $p = 2, \dots, 20$.

From Figures 1 and 2, we assess the efficiency of the error estimator in (5.2); for low degree $p$, the standard residual error estimator exhibits lower indices, while for larger values of $p$ the indices tend to be more similar. From Figure 3, we observe that the two error estimators are not $p$-robust (as one could have expected from the use of inverse estimates in the derivation of the lower bound (5.3b)); importantly, for increasing $p$, the effectivity indices converge to one another.
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Effectivity indices in (6.6). Test case with exact solution $u_1$ in (6.1); $p = 2, 3, 4, 5$; uniform mesh refinements. Left-panel: effectivity index $\mathfrak{F}_{\text{th}}$. Right-panel: effectivity index $\mathfrak{F}_{\text{DG}}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2:** Effectivity indices in (6.6). Test case with exact solution $u_1$ in (6.1); $p = 6, 7, 8, 9$; uniform mesh refinements. Left-panel: effectivity index $\mathfrak{F}_{\text{th}}$. Right-panel: effectivity index $\mathfrak{F}_{\text{DG}}$.</figcaption></figure>

## 6.3 Numerical results for an adaptive algorithm

We use here the error estimators in (5.2) and (5.28) to drive an adaptive algorithm of the form

SOLVE → ESTIMATE → MARK → REFINE

where

* in the **SOLVE** step, we employ a direct solver;

* in the **ESTIMATE** step, we employ the error estimators in (5.2) and (5.28)

* in the **MARK** step, we use Dörfler's marking strategy with marking parameter $1/2$;

* in the **REFINE** step, we use newest vertex bisection.

We run the adaptive algorithm on the test case with exact solution $u_1$ in (6.1) for fixed $p = 2$ and 5. In Figure 4, we report the error measures in (6.5), and the error estimators in (5.2) and (5.28) under uniform and adaptive mesh refinements.

Optimal convergence using the adaptive algorithm with both estimators is recovered; the expected suboptimal rates are observed for under uniform mesh refinements.
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Effectivity indices in (6.6). Test case with exact solution $u_1$ in (6.1); fixed structured mesh of 6 triangles; *p*-refinement. Effectivity indices $\mathcal{F}_{H_h}$ and $\mathcal{F}_{DG}$ in (6.6).</figcaption></figure>

# 7 Conclusions

We designed an error estimator for an IPDG discretization of the biharmonic problem, which does not depend explicitly on the stabilization of the method. Upper and lower bounds with respect to a suitable error measure are proven. Crucial tools in the analysis were a rewriting of the error measure based on the two dimensional div-div complex, quasi-interpolants in $C^1$ finite element spaces, and the right-inverse of the sym curl operator. The reliability and efficiency of the error estimator were verified in practice as well. An adaptive algorithm was finally used to construct sequences of geometrically refined meshes to approximate solutions with corner singularities.

Open problems that may be tackled in the future include the extension to the three dimensional case and the convergence analysis of the adaptive algorithm, which may be simpler than that for the standard IPDG error estimator due to the absence of explicit stabilizations in the proposed error estimator; cf. also Remark 5.2.

**Acknowledgments.** The authors are grateful to Johnny Guzmán for his kind advice on the right-inverse of the sym curl operator discussed in Appendix A. LM has been partially funded by MUR (PRIN2022 research grant n. 202292JW3F). LM was also partially supported by the European Union (ERC Synergy, NEMESIS, project number 101115663). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. LM is member of the Gruppo Nazionale Calcolo Scientifico-Istituto Nazionale di Alta Matematica (GNCS-INdAM).

# References

[1] R. Becker, P. Hansbo, and M. G. Larson. Energy norm a posteriori error estimation for discontinuous Galerkin methods. *Comput. Methods Appl. Mech. Engrg.*, 192(5-6):723–733, 2003.

[2] A. Bonito and R. H. Nochetto. Quasi-optimal convergence rate of an adaptive discontinuous Galerkin method. *SIAM J. Numer. Anal.*, 48(2):734–771, 2010.

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

[4] S. C. Brenner and L. R. Scott. *The Mathematical Theory of Finite Element Methods*, volume 15 of *Texts in Applied Mathematics*. Springer, New York, third edition, 2008.

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

[6] C. Carstensen, D. Gallistl, and J. Gedicke. Residual-based a posteriori error analysis for symmetric mixed Arnold-Winther FEM. *Numer. Math.*, 142(2):205–234, 2019.
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Test case with exact solution $u_1$ in (6.1); $p = 2, 5$; uniform and adaptive mesh refinements. Left-panel: $||D^2 u - H_h(u_h)||_{0,Ω}$; the adaptive algorithm is driven by the error estimator in (5.2). Right-panel: $||u - u_h||_{DG}$; the adaptive algorithm is driven by the error estimator in (5.28).</figcaption></figure>

[7] C. Carstensen and J. Hu. Hierarchical Argyris finite element method for adaptive and multigrid algorithms. *Comput. Methods Appl. Math.*, 21(3):529–556, 2021.

[8] T. Chaumont-Frelet. A new family of a posteriori error estimates for non-conforming finite element methods leading to stabilization-free error bounds. http://arxiv.org/abs/2506.23381, 2025.

[9] L. Chen and X. Huang. Finite elements for divdiv-conforming symmetric tensors. http://arxiv.org/abs/2005.01271, 2021.

[10] L. Chen and X. Huang. A finite element elasticity complex in three dimensions. *Math. Comp.*, 91(337):2095–2127, 2022.

[11] L. Chen and X. Huang. Finite element complexes in two dimensions. http://arxiv.org/abs/2206.00851, 2024.

[12] Ph. Clément. Approximation by finite element functions using local regularization. *RAIRO Anal. Numer.*, 9:77–84, 1975.

[13] M. Costabel and A. McIntosh. On Bogovskiǐ and regularized Poincaré integral operators for de Rham complex on Lipschitz domains. *Math. Z.*, 265:297–320, 2010.

[14] E. Dari, R. Durán, C. Padra, and V. Vampa. A posteriori error estimators for nonconforming finite element methods. *RAIRO Anal. Numer.*, 30(4):385–400, 1996.

[15] A. Dominicus, F. D. Gaspoz, and Ch. Kreuzer. Convergence of an adaptive $C^0$-interior penalty Galerkin method for the biharmonic problem. *IMA J. Numer. Anal.*, 2024.

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

[17] J. Douglas, 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. Numer.*, 13(3):227–255, 1979.

[18] A. Ern and J.-L. Guermond. *Finite Elements I—Approximation and Interpolation*, volume 72 of *Texts in Applied Mathematics*. Springer, Cham, 2021.

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

[20] E. H. Georgoulis, P. Houston, and J. Virtanen. An *a posteriori* error indicator for discontinuous Galerkin approximations of fourth-order elliptic problems. *IMA J. Numer. Anal.*, 31(1):281–298, 2011.

[21] V. Girault and L. R. Scott. Hermite interpolation of nonsmooth functions preserving boundary conditions. *Math. Comp.*, 71(239):1043–1074, 2002.

[22] B. Gräßle. Optimal multilevel adaptive FEM for the Argyris element. *Comput. Methods Appl. Mech. Engrg.*, 399:Paper No. 115352, 19, 2022.

[23] P. Grisvard. *Singularities in Boundary Value Problems*, volume 22 of *Recherches en Mathématiques Appliquées [Research in Applied Mathematics]*. Masson, Paris; Springer-Verlag, Berlin, 1992.

[24] O. A. Karakashian and F. Pascal. Convergence of adaptive discontinuous Galerkin approximations of second-order elliptic problems. *SIAM J. Numer. Anal.*, 45(2):641–665, 2007.
[25] Ch. Kreuzer and E. H. Georgoulis. Convergence of adaptive discontinuous Galerkin methods. *Math. Comp.*, 87(314):2611–2640, 2018.

[26] J. Morgan and R. Scott. A nodal basis for $C^1$ piecewise polynomials of degree $n \ge 5$. *Math. Comput.*, 29:736–740, 1975.

[27] B. Rivière and M. F. Wheeler. A posteriori error estimates for a discontinuous Galerkin method applied to elliptic problems. *Comput. Math. Appl.*, 46(1):141–163, 2003.

[28] R. Verfürth. *A Posteriori Error Estimation Techniques for Finite Element Methods*. Numerical Mathematics and Scientific Computation. Oxford University Press, Oxford, 2013.

# A Regularized right-inverse of the sym curl operator

Define

$$ \mathbb{H}(\mathbf{div}^0, \Omega) := \{ \mathbb{B} \in \mathbb{L}^2(\Omega) \mid \mathbf{div} \mathbb{B} = \mathbf{0} \}. $$

Recall the following two estimates from [13, Theorem 4.9] (with $s = -1$, and $s = 0$ row-wise). There exists positive constants $C_1(\Omega)$ and $C_2(\Omega)$ only depending on $\Omega$ such that, for all $\tau$ in $\mathbf{H}^{-1}(\Omega)$ with $\mathrm{div} \tau = 0$ and for all $\mathbb{B}$ in $\mathbb{H}(\mathbf{div}^0, \Omega)$, it is possible to construct $\omega$ in $L^2(\Omega)$ and $\theta$ in $\mathbf{H}^1(\Omega)$ satisfying

$$ \mathrm{curl} \omega = \tau, \qquad \|\omega\|_{0,\Omega} \le C_1(\Omega) \|\tau\|_{-1,\Omega}, \qquad (A.1) $$

and

$$ \mathrm{curl} \theta = \mathbb{B}, \qquad |\theta|_{1,\Omega} \le C_2(\Omega) \|\mathbb{B}\|_{0,\Omega}. \qquad (A.2) $$

We are in a position to prove the following result.

**Theorem A.1** (Regularized right-inverse of the sym curl operator). *The sym curl operator admits a stable right-inverse. More precisely, given $C_1(\Omega)$ and $C_2(\Omega)$ as in (A.1) and (A.2), for all $\mathbb{A}$ in $\mathbb{H}_S(\mathrm{div} \mathbf{div}^0, \Omega)$, it is possible to construct $\theta$ in $\mathbf{H}^1(\Omega)$ satisfying*

$$ \mathrm{sym} \, \mathrm{curl} \, \theta = \mathbb{A}, \qquad |\theta|_{1,\Omega} \le C_2(\Omega)(1+2C_1(\Omega)) \|\mathbb{A}\|_{0,\Omega}. \qquad (A.3) $$

*Proof.* Any $\mathbb{A}$ in $\mathbb{H}_S(\mathrm{div} \mathbf{div}^0, \Omega)$ is such that $\mathrm{div} \mathbb{A}$ is divergence free, and

$$ \mathrm{div} \mathbb{A} \in \mathbf{H}^{-1}(\Omega), \qquad \| \mathrm{div} \mathbb{A} \|_{-1,\Omega} \le \| \mathbb{A} \|_{0,\Omega}. \qquad (A.4) $$

We invoke (A.1), giving the existence of $\omega$ in $L^2(\Omega)$ such that

$$ \mathrm{curl} \omega = \mathrm{div} \mathbb{A}, \qquad \|\omega\|_{0,\Omega} \stackrel{(A.1)}{\le} C_1(\Omega) \| \mathrm{div} \mathbb{A} \|_{-1,\Omega} \stackrel{(A.4)}{\le} C_1(\Omega) \| \mathbb{A} \|_{0,\Omega}. \qquad (A.5) $$

Consider the skew-symmetric matrix

$$ \mathcal{O} := \begin{pmatrix} 0 & \omega \\ -\omega & 0 \end{pmatrix} \in \mathbb{L}^2(\Omega), \qquad \|\mathcal{O}\|_{0,\Omega} \le 2 \|\omega\|_{0,\Omega} \stackrel{(A.5)}{\le} 2C_1(\Omega) \|\mathbb{A}\|_{0,\Omega}, \qquad (A.6) $$

and observe that $\mathrm{div} \mathcal{O} = \mathrm{curl} \omega = \mathrm{div} \mathbb{A}$. We readily deduce that $\mathbb{A} - \mathcal{O}$ belongs to $\mathbb{H}(\mathbf{div}^0, \Omega)$. We invoke (A.2), giving the existence of $\theta$ in $\mathbf{H}^1(\Omega)$ such that

$$ \mathrm{curl} \theta = \mathbb{A} - \mathcal{O}, \qquad |\theta|_{1,\Omega} \stackrel{(A.2)}{\le} C_2(\Omega) \|\mathbb{A} - \mathcal{O}\|_{0,\Omega} \stackrel{(A.6)}{\le} C_2(\Omega)(1+2C_1(\Omega)) \|\mathbb{A}\|_{0,\Omega}. $$

Since $\mathbb{A}$ is symmetric and $\mathcal{O}$ is skew-symmetric, $\mathrm{sym} \, \mathrm{curl} \, \theta = \mathrm{sym}(\mathbb{A} - \mathcal{O}) = \mathbb{A}$. $\square$