$hp$-a posteriori error estimates for hybrid high-order methods applied to biharmonic problems

Problem Statement
This paper addresses the a posteriori error estimation for Hybrid High-Order (HHO) methods applied to the biharmonic problem ($\Delta^2 u = f$) with clamped boundary conditions ($u = \partial_n u = 0$ on $\partial\Omega$). The problem is posed on two- and three-dimensional polytopal Lipschitz domains with general topology (potentially containing holes). The goal is to derive residual-based $hp$-error estimators that are robust with respect to the mesh size ($h$), the polynomial degree ($p$ or $k$), and the domain topology.

Methodology
The authors develop a framework to bound the total error by decomposing it into conforming and nonconforming components.

1. Error Decomposition: The error is split into a conforming part (measured in the energy norm of a reconstructed $H^2_0$ function) and a nonconforming part (the distance between the discrete solution and the conforming space).
2. Nonconforming Error Control:
   • To bound the nonconforming error, the authors construct a conforming $H^2_0$ reconstruction of the discrete solution using a $C^1$-partition of unity.
   • This partition is built using $C^1$-conforming finite element spaces (specifically, the construction from [40] based on Alfeld splittings of the mesh).
   • Crucially, the authors employ local Helmholtz decompositions on vertex patches (stars). Unlike previous approaches that relied on global Helmholtz decompositions (which introduce constants depending on the number of holes in the domain), this local approach ensures the stability constant depends only on the mesh shape-regularity.
3. Conforming Error Control:
   • Two distinct upper bounds for the conforming error are derived, each utilizing a specific interpolation operator:
     • Setting 1: Uses the canonical hybrid finite element interpolation operator ($C^{k+2}_T$). This approach relies on the property that the composition of this operator with the HHO reconstruction yields an $H^2$-elliptic projection. The resulting bound involves only the stabilization term and data oscillation, but it hinges on a numerical assumption regarding the interpolation error (Assumption 3.4).
     • Setting 2: Uses the modified Babuška–Suri $hp$-interpolation operator ($I^{k+2}_{mBS}$). This approach does not require additional assumptions but introduces additional error indicator terms, specifically bulk residuals and normal/tangential flux jumps.
4. Residual Estimators: The paper defines local error indicators comprising stabilization terms, bulk residuals, normal flux jumps, and tangential jumps, scaled by appropriate $hp$-factors.

Key Contributions

• First $hp$-estimator for HHO on biharmonic problems: This work provides the first residual-based $hp$-a posteriori error estimator for HHO methods applied to the biharmonic problem on general 2D and 3D domains.
• Local Helmholtz Decomposition with Topological Robustness: The paper establishes a novel stability estimate for local Helmholtz decompositions on vertex patches. It proves that the stability constant depends solely on the mesh shape-regularity, thereby avoiding the degradation of constants associated with the number of holes in the domain, a known issue in global decomposition approaches.
• Dual Upper Bounds: The derivation of two distinct upper bounds for the conforming error offers flexibility. One bound is "clean" (relying only on stabilization and oscillation) but conditional on a numerically verified assumption; the other is unconditional but includes more complex residual terms.
• Extension to General Topology: The methodology is explicitly designed to handle domains with arbitrary topology (multiple holes), a feature often difficult to address in $hp$-estimation for fourth-order problems.

Results

• Theoretical Bounds: The main theorem (Theorem 5.6) establishes that the total error (broken Hessian norm plus stabilization) is bounded by the sum of the proposed error indicators. The bounds are shown to be $h$-optimal and $p$-suboptimal (by at most one order, or potentially $1/2$ order under the verified assumption).
• Numerical Verification:
  • Assumption Verification: Numerical experiments in 2D and 3D confirm the validity of Assumption 3.4 regarding the interpolation error of the canonical hybrid operator.
  • Convergence Rates: Tests on smooth solutions in 3D demonstrate that both the energy error and the estimator converge at the optimal rate $O(\text{DoFs}^{-(k+1)/d})$.
  • Effectivity Indices: The estimators yield moderate effectivity indices (ranging roughly between 1.2 and 2.0) across various polynomial degrees and mesh refinements.
  • Topological Robustness: Experiments on 2D domains with increasing numbers of holes (0, 1, 4, 8) show that the effectivity indices remain stable, confirming the theoretical claim that the estimator is robust with respect to domain topology.
  • Adaptive Performance: $h$-adaptive algorithms driven by the estimator successfully recover optimal convergence rates for singular solutions in both 2D (L-shaped domain) and 3D (point singularity), outperforming uniform refinement.

Significance
The authors position this work as the first to employ a $C^1$-partition of unity combined with local Helmholtz decompositions on vertex stars to bound nonconforming errors in the $hp$-setting for biharmonic problems. The primary significance lies in the ability to control the nonconforming error without incurring stability constants that deteriorate with the domain's topological complexity (number of holes). Furthermore, the authors note that while the approach is developed for HHO methods, the technique for bounding nonconforming errors via local Helmholtz decompositions is potentially applicable to other nonconforming methods, such as Interior Penalty Discontinuous Galerkin (IPDG) methods. The paper does not claim to provide lower error bounds, noting that such bounds would follow from adaptations of existing arguments but would likely exhibit pessimistic $p$-suboptimality due to the limitations of current $hp$-techniques for constructing specific bubble functions and extension operators.