Hybridizable Staggered Discontinuous Galerkin Methods for Polyharmonic Equations on Polytopes

Abstract: Hybridizable staggered discontinuous Galerkin methods are developed for arbitrary-order polyharmonic equations $(-Δ)^m u=f$ on shape-regular polytopal meshes in $\mathbb R<sup>d$, for any $m\ge1$, $d\ge2$, and polynomial degree $k\ge0$. The method uses the mixed variable $σ=\nabla^m u$ and a staggered primal--dual mesh to impose complementary continuity on scalar and tensor unknowns, without restrictions such as $d\ge m$. Local trace and bubble enrichments stabilize low-order tensor spaces without adding global unknowns. Hybridization localizes the tensor variable and yields an equivalent stabilization-free weak Galerkin formulation. Well-posedness and optimal energy error estimates are proved, and numerical experiments on polygonal and tetrahedral meshes confirm the predicted rates.

• The paper introduces a novel HSDG method leveraging a staggered primal-dual framework to efficiently solve arbitrary-order polyharmonic equations on polytopal meshes.
• It employs local polynomial enrichments and hybridization to achieve stability and optimal convergence without incurring extra global degrees of freedom.
• Numerical experiments in two and three dimensions demonstrate that the method attains order k+1 convergence in both energy and L2 norms for a variety of mesh types.

Hybridizable Staggered Discontinuous Galerkin Methods for Polyharmonic Equations on Polytopes

Introduction and Motivation

This work introduces a unified class of hybridizable staggered discontinuous Galerkin (HSDG) methods for solving arbitrary-order polyharmonic equations, $(-\Delta)^m u = f$, on general shape-regular polytopal meshes in $\mathbb{R}^d$. The formulation is dimension-agnostic ($d \geq 2$), supports arbitrary differential order ($m \geq 1$), and arbitrary polynomial degree ($k \geq 0$). By extending the mixed finite element framework to leverage a staggered primal-dual mesh, the approach circumvents the difficulties in constructing explicit $H(\operatorname{div}^m)$-conforming tensor finite elements for high-order $m$ and polytopal geometries. Local enrichments recover stability for all discretization orders without incurring extra global degrees of freedom, and hybridization yields an efficient, stabilization-free weak Galerkin (WG) scheme.

Staggered Discretization and Mixed Formulation

The fundamental starting point is a mixed system with $\sigma = \nabla^m u$ as a symmetric $m$-th order tensor in $H(\operatorname{div}^m)$: $\mathbb{R}^d$0 posed with homogeneous Dirichlet and normal derivative boundary conditions up to order $\mathbb{R}^d$1. The mixed weak formulation only requires $\mathbb{R}^d$2 and $\mathbb{R}^d$3 with $\mathbb{R}^d$4-regularity. Classical mixed FEM is available for $\mathbb{R}^d$5 (RT, BDM, Nédélec elements), but explicit $\mathbb{R}^d$6-conforming elements become systematically intractable for higher $\mathbb{R}^d$7 and general polytopes.

The SDG approach exploits staggered compatibility: the scalar variable is elementwise $\mathbb{R}^d$8-conforming (on the primal mesh), and the tensor variable achieves $\mathbb{R}^d$9-type conformity on dual elements, inducing complementary local regularity. Both variables are globally discontinuous, yet continuity is imposed locally in a staggered manner, allowing for stability and local conservation without global $d \geq 2$0-conformity.

Construction of Finite Element Spaces

A key innovation lies in the explicit construction and enrichment of local polynomial tensor spaces $d \geq 2$1. The scalar space on each element is $d \geq 2$2. For the tensor variable, a comprehensive geometric and combinatorial analysis leads to the following features:

• Simplicial lattice indexing: Symmetric tensors are indexed by nodes of a lattice, allowing precise tracking of normal and tangential components relative to mesh faces.
• Trace enrichment: For $d \geq 2$3, additional normal trace layers are explicitly constructed to recover essential boundary continuity.
• Bubble enrichment: Sufficient tensor directions (up to order $d \geq 2$4, with $d \geq 2$5 the number of non-parallel face normals of an element $d \geq 2$6) are included to guarantee local inf-sup stability for all $d \geq 2$7.

The approach thus avoids the need for high-degree $d \geq 2$8-conforming polynomial spaces or the virtual element machinery with stabilization terms.

Hybridization and Weak Galerkin Interpretation

Hybridization is performed by relaxing inter-element continuity of tensor traces and introducing scalar trace unknowns on primal (macro) faces. The resulting global system involves only scalar and face unknowns after static condensation—local tensor variables are eliminated element-wise. The final formulation is equivalent to a stabilization-free weak Galerkin method, where the local projection is onto a space richer than standard polynomials (encompassing barycentric splits), removing the need for classical VEM-style stabilization.

This weak Galerkin interpretation generalizes conforming and nonconforming VEMs, sidesteps explicit stabilization, and unifies the staggered DG methodology for arbitrary order.

Well-posedness and Error Analysis

The analysis is performed under standard shape-regularity, with constants uniform in $d \geq 2$9, depending only on the polynomial degree $m \geq 1$0, order $m \geq 1$1, dimension $m \geq 1$2, and mesh regularity. The following properties are established:

• Coercivity of the discrete weak gradient: The discrete weak gradient operator is $m \geq 1$3-norm coercive, guaranteeing the uniqueness and stability of the scalar problem post-elimination.
• Discrete inf-sup condition: The construction provides a uniform inf-sup constant, ensuring mixed method stability.
• Optimal convergence: For sufficiently regular solutions $m \geq 1$4, the energy norm error behaves as $m \geq 1$5 for both $m \geq 1$6 and $m \geq 1$7 in $m \geq 1$8-seminorm, with higher than first order $m \geq 1$9 convergence.

Numerical Results

Extensive numerical experiments are carried out in two and three dimensions, for biharmonic ($k \geq 0$0) and triharmonic ($k \geq 0$1) problems, using both convex and concave polygonal meshes as well as tetrahedral partitions in 3D. The observed convergence rates strictly match the theoretical predictions:

• For the energy-type errors $k \geq 0$2 and $k \geq 0$3, convergence of order $k \geq 0$4 is attained.
• The $k \geq 0$5 error of $k \geq 0$6 consistently exhibits rates exceeding $k \geq 0$7, as expected from regularity theory.

Notably, even for low $k \geq 0$8 or irregular meshes (with minimal regularity/bubble/traces), the method's stability is robust, provided the element enrichment rules are observed.

Implications and Future Directions

The proposed HSDG method constitutes a principled, unified framework for high-order polyharmonic PDEs on general polytopal meshes. Its reliance on staggered continuity, local enrichment, and hybridization can be further generalized:

• Adaptivity and hp-enrichment: The geometric and algebraic flexibility makes it natural to incorporate mesh adaptivity, variable $k \geq 0$9, and $H(\operatorname{div}^m)$0-enrichment strategies.
• Extension to nonlinear or time-dependent problems: The staggered and weak Galerkin approach can be adapted for higher-order nonlinear PDEs or evolutionary settings.
• Application to complex domains: Since the method is polytopal-agnostic, it is well-suited for interface and multi-material problems where mesh generality is essential.

Potential future work might address efficient solver design for the condensed global system and extension to coupled physical systems where high-order regularity arises in a natural fashion.

Conclusion

This paper presents a hybridizable staggered DG approach for polyharmonic equations that is both theoretically rigorous and practically robust across dimensions and mesh types. By foundationally exploiting staggered conformity, explicit local enrichment, and hybridization, the method achieves optimal convergence and stability while circumventing the limitations of classical high-order conforming or nonconforming FEMs on polytopal domains. The resulting framework offers a scalable and highly generalizable route for simulation of high-order elliptic problems in computational mathematics and engineering.