# Hybrid high-order approximations of div-curl systems on domains with general topology

Jérémy Dalphin, Jean-Pierre Ducreux, Simon Lemaire, Silvano Pitassi

► To cite this version:

Jérémy Dalphin, Jean-Pierre Ducreux, Simon Lemaire, Silvano Pitassi. Hybrid high-order approximations of div-curl systems on domains with general topology. SIAM Journal on Scientific Computing, 2025, 47 (5), pp.A2988-A3016. 10.1137/25M1744356. hal-04663200v3

HAL Id: hal-04663200

https://hal.science/hal-04663200v3

Submitted on 1 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>
# Hybrid high-order approximations of div-curl systems on domains with general topology

Jérémy Dalphin*¹, Jean-Pierre Ducreux†¹, Simon Lemaire‡², and Silvano Pitassi§²,¹

¹ERMES department, EDF R&D, EDF Lab Paris-Saclay, 91120 Palaiseau, France

²Inria, Univ. Lille, CNRS, UMR 8524 - Laboratoire Paul Painlevé, 59000 Lille, France

June 21, 2025

## Abstract

We devise and analyze hybrid polyhedral methods of arbitrary order for the approximation of div-curl systems on three-dimensional domains featuring non-trivial topology. The div-curl systems we are interested in stem from magnetostatics, and can either be first-order (field formulation) or second-order (vector potential formulation). The well-posedness of the resulting discrete problems essentially hinges on recently established, topologically generic, hybrid versions of the (first and second) Weber inequalities. Our error analysis covers the case of regular solutions. Leveraging (co)homology computation techniques from the literature, we perform an in-depth numerical assessment of our approach, covering, in particular, the case of non-simply-connected domains.

**Keywords:** Div-curl systems; Polyhedral meshes; Hybrid methods; de Rham cohomology; Computational topology.

**AMS Subject Classification 2020:** 65N12, 14F40, 35Q60.

## 1 Introduction

Let $\Omega$ be a domain in $\mathbb{R}^3$, i.e. a bounded and connected Lipschitz open set of $\mathbb{R}^3$. In the $L^2$ framework, the standard (primal) de Rham complex reads:

$$ \{0\} \xrightarrow{0} H^1(\Omega) \xrightarrow{\text{grad}} H(\text{curl}; \Omega) \xrightarrow{\text{curl}} H(\text{div}; \Omega) \xrightarrow{\text{div}} L^2(\Omega) \xrightarrow{0} \{0\}, \quad (1) $$

with cohomology spaces

$$ \begin{aligned} \mathfrak{S}^0 &:= \ker(\text{grad})/(\operatorname{Im}(0)), & \mathfrak{S}^1 &:= \ker(\text{curl})/(\operatorname{Im}(\text{grad})), \\ \mathfrak{S}^2 &:= \ker(\text{div})/(\operatorname{Im}(\text{curl})), & \mathfrak{S}^3 &:= (\ker(0)) = L^2(\Omega)/(\operatorname{Im}(\text{div})). \end{aligned} \quad (2) $$

The spaces $\mathfrak{S}^n$ have respective dimensions equal to the Betti numbers $\beta_n$, with $\beta_0 = 1$ (number of connected components of $\Omega$), $\beta_1 \in \mathbb{N}$ (number of tunnels crossing through $\Omega$), $\beta_2 \in \mathbb{N}$ (number of voids encapsulated by $\Omega$), and $\beta_3 = 0$. In turn, the dual de Rham complex reads:

$$ \{0\} \xrightarrow{0} H_0^1(\Omega) \xrightarrow{\text{grad}_0} H_0(\text{curl}; \Omega) \xrightarrow{\text{curl}_0} H_0(\text{div}; \Omega) \xrightarrow{\text{div}_0} L^2(\Omega) \xrightarrow{0} \{0\}, \quad (3) $$

<footnote>*jeremy.dalphin@edf.fr</footnote>

<footnote>†jean-pierre.ducreux@edf.fr</footnote>

<footnote>‡simon.lemaire@inria.fr (corresponding author)</footnote>

<footnote>§silvano.pitassi@inria.fr</footnote>
where the operators **grad**₀, **curl**₀, and div₀ are respectively L²-adjoint to the operators – div, **curl**, and – **grad**. The zero subscripts aim to recall that the domains of the adjoint differential operators embed zero boundary conditions. The corresponding homology spaces are

$$
\begin{align}
\mathfrak{H}_0^0 &:= \ker(\mathbf{grad}_0) / \operatorname{Im}(0), & \mathfrak{H}_0^1 &:= \ker(\mathbf{curl}_0) / \operatorname{Im}(\mathbf{grad}_0), \tag{4} \\
\mathfrak{H}_0^2 &:= \ker(\operatorname{div}_0) / \operatorname{Im}(\mathbf{curl}_0), & \mathfrak{H}_0^3 &:= (\ker(0))^{\mathrm{L}^2(\Omega)} / \operatorname{Im}(\operatorname{div}_0). \nonumber
\end{align}
$$

Importantly, there holds $\mathfrak{H}_0^n = \mathfrak{H}^{3-n}$ (this is the so-called *Poincaré–Lefschetz duality*).

Let $\mu$ be a possibly varying coefficient, assumed to be uniformly bounded by above and by below by some positive constants. In what follows, the symbols $\perp$ and $\perp_\mu$ respectively denote standard and $\mu$-weighted $L^2$-orthogonality. Let the sequence of differential operators $\{\mathrm{Grad}, \mathrm{Curl}, \mathrm{Div}\}$ stand either (i) for $\{\mathrm{grad}, \mathrm{curl}, \mathrm{div}\}$ from the primal de Rham complex (1), or (ii) for $\{\mathrm{grad}_0, \mathrm{curl}_0, \mathrm{div}_0\}$ from the dual de Rham complex (3). In the first case, let $\mathfrak{h}^n := \mathfrak{H}^n$; in the second, let $\mathfrak{h}^n := \mathfrak{H}_0^n = \mathfrak{H}^{3-n}$. We aim to study problems of the following form. Given some datum $j \in \ker(\mathrm{Div})$, find $h \in \mathrm{Dom}(\mathrm{Curl}) \cap \mathrm{Im}(\mathrm{Grad})^{\perp\mu}$ s.t.

$$
\mathrm{Curl}(h) = j. \tag{5}
$$

One can easily verify that **Curl**, seen as a bounded linear mapping from $\mathrm{Dom}(\mathrm{Curl}) \cap \mathrm{Im}(\mathrm{Grad})^{\perp\mu}$ to $\mathrm{Ker}(\mathrm{Div})$, is a Fredholm operator, with null space $\mathrm{Ker}(\mathrm{Curl}) \cap \mathrm{Im}(\mathrm{Grad})^{\perp\mu} \cong \mathfrak{h}^1$, and defect space $\mathrm{Ker}(\mathrm{Div})/\mathrm{Im}(\mathrm{Curl}) = \mathfrak{h}^2$. Thus, the following alternative holds true for Problem (5):

• either $\dim(\mathfrak{h}^1) = \dim(\mathfrak{h}^2) = 0$, then Problem (5) admits a unique solution;

* in the opposite case, for Problem (5) to admit a solution, it is necessary that $j \in \mathrm{Im}(\mathrm{Curl})$, which amounts, for $j \in \mathrm{Ker}(\mathrm{Div})$, to additionally satisfying $j \perp \mathfrak{h}^2$; the uniqueness of the solution is then recovered by imposing $\dim(\mathfrak{h}^1)$ additional constraints to $h$ (which already satisfies $h \perp_\mu \mathrm{Im}(\mathrm{Grad}))$.

Problem (5) is prototypical of the first-order formulation (also known as *field formulation*) of magnetostatics. In practice, $j : \Omega \to \mathbb{R}^3$ is a given electric current density, $\mu : \Omega \to \mathbb{R}_{>0}$ is the magnetic permeability of the medium, and $h : \Omega \to \mathbb{R}^3$ is the sought magnetic field. The model can be endowed with either normal or tangential boundary conditions, depending on which version (respectively, (1) or (3)) of the de Rham complex the sequence of differential operators stems from. We refer the reader to Section 4.1.1 for a somewhat more conventional rewriting of Problem (5) in the case of normal boundary conditions.

Let us now introduce **Curl**\*, L²-adjoint of **Curl**. Recall that $\mathrm{Ker}(\mathrm{Curl})^{\perp} = \mathrm{Im}(\mathrm{Curl}^*)$. If the solution $h \perp_{\mu} \mathrm{Im}(\mathrm{Grad})$ to Problem (5) (assuming it exists) additionally satisfies $h \perp_{\mu} \mathfrak{h}^1$, then it can be written $h = \mu^{-1}\mathrm{Curl}^*(a)$ for some vector potential $a \in \mathrm{Dom}(\mathrm{Curl}^*)$. Obviously, the vector potential $a$ is non-unique; it is only defined up to an element of $\mathrm{Ker}(\mathrm{Curl}^*)$ (gauge). Gauging out to zero, Problem (5) equivalently rewrites: find $a \in \mathrm{Dom}(\mathrm{Curl}^*) \cap \mathrm{Ker}(\mathrm{Curl}^*)^{\perp}$ s.t.

$$
\mathrm{Curl}(\mu^{-1}\mathrm{Curl}^*(a)) = j, \tag{6}
$$

and the (unique) solution to this problem does exist as soon as $j \in \ker(\mathrm{Div})$ additionally satisfies $j \perp \mathfrak{h}^2$. Problem (6) is prototypical of the second-order formulation (also known as *vector potential formulation*) of magnetostatics. In practice, $a : \Omega \to \mathbb{R}^3$ is the magnetic vector potential. When Problem (5) is endowed with normal (resp. tangential) boundary conditions, Problem (6) is endowed with tangential (resp. normal) boundary conditions. In the second scenario, part of the boundary conditions in Problem (6) are actually natural, and encode the tangential prescription on $h$. For a more conventional rewriting of Problem (6) in the case of tangential boundary conditions (first scenario), we refer to Section 4.1.2.
In this work, we aim at devising arbitrary-order *hybrid* polyhedral discretizations of the first- and second-order model problems (5) and (6). By definition, hybrid methods only attach unknowns to the faces and to the cells of the spatial partition at hand. Examples of such approaches are the Hybridizable Discontinuous Galerkin (HDG) [17], the Weak Galerkin (WG) [39], the Hybrid High-Order (HHO) [22, 23], or the non-conforming Virtual Element (ncVE) [35, 5] methods. As might be expected, and as was first documented in [15, 16], all these technologies are very tightly connected. The devising of hybrid polyhedral methods for div-curl systems has already been addressed to some extent in the HDG and WG literatures. For second-order models, these contributions include [37, 13, 12] (cf. also [24]) for HDG, and [36] for WG, all of them restricted to trivial topologies. Regarding now first-order models, the only contributions we are aware of are [34], in which a primal WG method is introduced, and [38, 9, 10], in which alternative (primal-dual) WG approaches are studied. In all these contributions, full face polynomial spaces are employed to discretize the vectorial variable. Also, either the case of non-trivial topologies is not covered (like in [34]), or it is treated but in a non-robust fashion. Recently, the discretization of both first- and second-order models of magnetostatics has been studied in the HHO context [11]. Therein, the vectorial variable is sought into trimmed face polynomial spaces, reusing and extending ideas introduced in [13] (cf. also [32]). The contribution [11], however, solely covers the case of trivial topologies. Our objective in the present work is to fill this gap. We aim at devising, analyzing, and numerically assessing trimmed HHO methods for (first- and second-order) models set in domains with arbitrary topology. With respect to [38, 9, 10], apart from using trimmed face spaces, our approach for first-order models is also based on a different, more suitable variational formulation (inspired from [31]), which allows for both a leaner construction and a robust handling of harmonic fields. In turn, our stability analysis hinges on the systematic use of the (topologically generic) hybrid Weber inequalities recently established in [33] by the last two authors. Finally, our error analysis (cf. Theorems 7 and 12) covers the case of regular solutions. Let us emphasize that [33] is exclusively concerned with general-purpose discrete functional analysis; it does not address the practical devising of hybrid methods, which is the subject of the present work. To the best of our knowledge, in the realm of polyhedral approaches, our work is the first to deal, from the devising to the implementation in 3D, and in a systematic and robust way, with non-trivial topologies. Let us stress that, from an application point of view, non-trivial topologies are ubiquitous (think, e.g., of the toroidal vacuum chamber of a tokamak). Their proper handling is thus of paramount importance. From an algorithmic point of view, arbitrary topologies come with a number of additional difficulties. It is especially true for non-simply-connected domains, for which discrete cutting surfaces might need to be computed from the polyhedral mesh at hand (see Figure 1).

The material is organized as follows. In Section 2, we introduce both the topological and functional frameworks, and we recall the first and second Helmholtz–Hodge decompositions. In Section 3, we introduce the discrete setting, including polyhedral discretizations, polynomial decompositions, and hybrid spaces. In Section 4, we devise and analyze HHO methods for (first- and second-order) model problems of the forms (5) and (6). Finally, in Section 5, we provide a comprehensive set of numerical experiments on non-trivial 3D domains, assessing the relevance of our methodology.

# 2 Preliminaries

## 2.1 Topological framework

We recall that $\Omega$ denotes a domain in $\mathbb{R}^3$, that is a bounded and connected Lipschitz open set of $\mathbb{R}^3$. Let $\Gamma := \partial\Omega$ denote its boundary, and $n: \Gamma \to \mathbb{R}^3$ be the (almost everywhere defined) unit vector field normal to $\Gamma$, pointing outward from $\Omega$. We recall that the Betti numbers $\beta_1$ and $\beta_2$ of $\Omega$ respectively
count the number of tunnels crossing through $\Omega$ ($\beta_1 \in \mathbb{N}$), and the number of voids encapsulated by $\Omega$ ($\beta_2 \in \mathbb{N}$). For $\Omega$ simply-connected, $\beta_1 = 0$. Likewise, when the boundary $\Gamma$ of $\Omega$ is connected, $\beta_2 = 0$. In what follows, when both $\beta_1, \beta_2$ are equal to zero, we say that the topology of $\Omega$ is trivial.

Whenever $\beta_1 > 0$, we make the following classical assumption: there exist $\beta_1$ orientable and connected two-dimensional manifolds with boundary, denoted $\Sigma_1, \dots, \Sigma_{\beta_1}$, called cutting surfaces, satisfying $\Sigma_i \subset \Omega$ and $\partial\Sigma_i \subset \Gamma$ for all $i \in \{1, \dots, \beta_1\}$, such that the open set $\hat{\Omega} := \Omega \setminus \bigcup_{i \in \{1, \dots, \beta_1\}} \Sigma_i$ is simply-connected, i.e. its first Betti number is zero (see Figure 4 for an example). We will assume, in what follows, that $\hat{\Omega}$ is connected (which is generally the case starting from a connected domain $\Omega$), and that the cutting surfaces are sufficiently regular so that the set $\hat{\Omega}$ is pseudo-Lipschitz (cf. [4, Def. 3.2.2]). For any $i \in \{1, \dots, \beta_1\}$, we let $n_{\Sigma_i} : \Sigma_i \to \mathbb{R}^3$ be the (almost everywhere defined) unit vector field normal to $\Sigma_i$, the orientation of which we arbitrarily prescribe (the cutting surfaces $\Sigma_i$ are supposed to be orientable). The orientation of $n_{\Sigma_i}$ being prescribed, we associate the tag “+” (resp. “−”) to the side of $\hat{\Omega}$ (with respect to $\Sigma_i$) for which $n_{\Sigma_i}$ points outward (resp. inward).

**Remark 1 (Cutting surfaces).** In [7], the class of (so-called) weakly-Helmholtz domains is introduced. Its definition does not assume that the cut domain $\hat{\Omega}$ is simply-connected, which allows to cover more exotic topologies (e.g., the domain obtained as the complement in a cube of a trefoil knot). Our analysis could be extended to cover weakly-Helmholtz domains. However, whereas weakly-Helmholtz domains assume non-intersecting cutting surfaces, we do not make such an assumption here. This allows us to consider the non-weakly-Helmholtz, yet physically relevant configuration of a toroidal vacuum chamber (see Figure 4), for which our approach seamlessly applies, as supported by the numerical results in Section 5.3.

Whenever $\beta_2 > 0$, letting $\Gamma_0$ be the (connected) boundary of the only unbounded component of the exterior open set $\mathbb{R}^3 \setminus \bar{\Omega}$, there exist $\beta_2$ (maximally) connected components $\Gamma_1, \dots, \Gamma_{\beta_2}$ of $\Gamma$ such that $\Gamma = \cup_{j \in \{0, \dots, \beta_2\}} \Gamma_j$ (cf. Figure 4). If $\beta_2 = 0$, there holds $\Gamma = \Gamma_0$.

## 2.2 Functional framework

Let $\mu : \Omega \to \mathbb{R}$ be a given function satisfying, for real numbers $0 < \mu_b \le \mu_\# < \infty$,

$$ \mu_b \le \mu(\mathbf{x}) \le \mu_\# \quad \text{for a.e. } \mathbf{x} \in \Omega. \tag{7} $$

For $m \in \{2, 3\}$, and for $X$ an $m$-dimensional, (relatively) open pseudo-Lipschitz subset of $\bar{\Omega}$, we let $L^2(X)$ (respectively, $L^2(X)$) denote the Lebesgue space of square-integrable functions (respectively, $\mathbb{R}^m$-valued vector fields) over $X$. The standard inner products (and norms) in $L^2(X)$ and $L^2(X)$ are irrespectively denoted by $(\cdot, \mathfrak{g})_X := \int_X \mathfrak{f} \cdot \mathfrak{g}$ (and $\|\cdot\|_{0,X} := \sqrt{\langle \cdot, \cdot \rangle_X}$). We also define $L_0^2(X) := \{\nu \in L^2(X) \mid \int_X \nu = 0\}$ and $L_0^2(X) := \{\nu \in L^2(X) \mid \int_X \nu = 0\}$. For $s > 0$, we let $H^s(X)$ (resp. $H^s(X)$) denote the Sobolev space of functions in $L^2(X)$ (resp. $\mathbb{R}^m$-valued vector fields in $L^2(X)$) possessing square-integrable partial weak derivatives up to order $s$ over $X$ (for fractional $s$, we follow the classical Sobolev–Slobodeckij construction). The standard norms (and semi-norms) in $H^s(X)$ and $H^s(X)$ are irrespectively denoted by $\|\cdot\|_{s,X}$ (and $|\cdot|_{s,X}$). Given a disjoint Lipschitz partition $\mathcal{K}$ of $\Omega$ (in the sense that every $K \in \mathcal{K}$ is an open Lipschitz set, $K_1 \cap K_2 = \emptyset$ for all $K_1, K_2 \in \mathcal{K}$ with $K_1 \ne K_2$, and $\cup_{K \in \mathcal{K}} \bar{K} = \bar{\Omega}$), we also define the broken Sobolev space $H^s(\mathcal{K}) := \{\nu \in L^2(\Omega) : \nu|_{K \in \mathcal{K}} \in H^s(K) \forall K \in \mathcal{K}\}$.

Let $Y$ be a three-dimensional, open Lipschitz subset of $\Omega$. Classically, we let

$$ H(\mathrm{curl}; Y) := \{\nu \in L^2(Y) \mid \mathrm{curl} \, \nu \in L^2(Y)\}, $$

$$ H(\mathrm{div}_{\mu}; Y) := \{\nu \in L^2(Y) \mid \mathrm{div}(\mu\nu) \in L^2(Y)\}, $$
with $H(\operatorname{div}; Y) := H(\operatorname{div}_1; Y)$, as well as their two subspaces $H(\operatorname{curl}^0; Y) := \{v \in H(\operatorname{curl}; Y) | \operatorname{curl} v \equiv 0\}$ and $H(\operatorname{div}_\mu^0; Y) := \{v \in H(\operatorname{div}_\mu; Y) | \operatorname{div}(\mu v) \equiv 0\}$. Let $n_{\partial Y} : \partial Y \to \mathbb{R}^3$ denote the (almost everywhere defined) outward unit normal vector to $\partial Y$. For $v \in H(\operatorname{div}_\mu; Y)$, the normal trace of $\mu v$ on $\partial Y$ can be defined as an element of $H^{-\frac{1}{2}}(\partial Y)$ (space of bounded linear forms on $H^{\frac{1}{2}}(\partial Y)$), denoted $(\mu v)|_{\partial Y} \cdot n_{\partial Y}$. Likewise, for $v \in H(\operatorname{curl}; Y)$, one can give a sense to the rotated tangential trace of $v$ on $\partial Y$ as an element of $H^{-\frac{1}{2}}(\partial Y)^3$ (space of bounded linear forms on $H^{\frac{1}{2}}(\partial Y)^3$), denoted $v|_{\partial Y} \times n_{\partial Y}$. Following [33, Rmk. 1], from now on, tangential vector fields are identified to two-dimensional vector fields. The relevant boundary functional spaces then become $L^2(\partial Y) = L^2(\partial Y)^2$, $H^{\frac{1}{2}}(\partial Y) = H^{\frac{1}{2}}(\partial Y)^2$, and $H^{-\frac{1}{2}}(\partial Y) := H^{-\frac{1}{2}}(\partial Y)^2$, and one may abuse the notation and write $v|_{\partial Y} \times n_{\partial Y} \in H^{-\frac{1}{2}}(\partial Y)$. In what follows, the duality pairings between $H^{-\frac{1}{2}}(\partial Y)$ and $H^{\frac{1}{2}}(\partial Y)$ on the one side, and between $H^{-\frac{1}{2}}(\partial Y)$ and $H^{\frac{1}{2}}(\partial Y)$ on the other side, are irrespectively denoted by $\langle \cdot, \cdot \rangle_{\partial Y}$. If $v \in H(\operatorname{curl}; Y) \cap H(\operatorname{div}; Y) \cap H^s(Y)$ for some $s > \frac{1}{2}$, then for a.e. $x \in \partial Y$, $v|_{\partial Y}(x) = (v|_{\partial Y} \cdot n_{\partial Y})(x)n_{\partial Y}(x) + n_{\partial Y}(x) \times (v|_{\partial Y} \times n_{\partial Y})(x)$. In this case, $v|_{\partial Y} \cdot n_{\partial Y} \in L^2(\partial Y)$, and $v|_{\partial Y} \times n_{\partial Y} \in L^2(\partial Y)$.

To account for essential boundary conditions, we shall define the subspaces

$$H_0(\operatorname{curl}; \Omega) := \{v \in H(\operatorname{curl}; \Omega) | v|_{\Gamma} \times n \equiv 0\}$$

(along with $H_0(\operatorname{curl}^0; \Omega) := H_0(\operatorname{curl}; \Omega) \cap H(\operatorname{curl}^0; \Omega)$), and

$$H_0(\operatorname{div}_\mu; \Omega) := \{v \in H(\operatorname{div}_\mu; \Omega) | (\mu v)|_{\Gamma} \cdot n \equiv 0\}$$

(with $H_0(\operatorname{div}_\mu^0; \Omega) := H_0(\operatorname{div}_\mu; \Omega) \cap H(\operatorname{div}_\mu^0; \Omega)$). We also set $H_0^1(\Omega) := \{v \in H^1(\Omega) | v|_{\Gamma} \equiv 0\}$. Last, assume that the first Betti number $\beta_1$ of $\Omega$ is positive. Then, for $v \in L^2(\hat{\Omega})$ (resp. for $v \in L^2(\hat{\Omega})$), we denote by $\check{v}$ (resp. by $\check{v}$) its continuation to $L^2(\Omega)$ (resp. to $L^2(\Omega)$). Also, for any $v : \hat{\Omega} \to \mathbb{R}$, and $i \in \{1, \dots, \beta_1\}$, denoting $v|_{\Sigma_i}$ and $v|_{\bar{\Sigma}_i}$ the traces of $v$ on $\Sigma_i$ defined (if need be, in a weak sense) from both sides of $\hat{\Omega}$ (respectively tagged by “+” and “-”), we define the jump of $v$ across $\Sigma_i$ by

$$[v]_{\Sigma_i} := v|_{\Sigma_i} - v|_{\bar{\Sigma}_i}. \qquad (8)$$

## 2.3 Helmholtz-Hodge decompositions

We collect classical results about Helmholtz-Hodge decompositions; for further details, we refer the reader to [28, 20, 2] (cf. also [29, 4], as well as [33, Sec. 2.3]). Before proceeding, we define $L^2_\mu(\Omega) := (L^2(\Omega), (\mu \cdot, \cdot)|_{\partial \Omega})$.

### 2.3.1 1st Helmholtz-Hodge decomposition

Let us first consider the harmonic space $H_\tau(\Omega) := H_0(\operatorname{curl}^0; \Omega) \cap H(\operatorname{div}_\mu^0; \Omega)$ which, by (4), satisfies $H_\tau(\Omega) \cong \mathfrak{F}_0^1$. In particular, the space $H_\tau(\Omega)$ has dimension $\beta_2$, and it can be proved that vector fields $w \in H_\tau(\Omega)$ are entirely characterized by the data of $(\langle (\mu w)|_{\Gamma_j} \cdot n, 1 \rangle_{\Gamma_j})_{j \in \{1, \dots, \beta_2\}}$, where $\langle \cdot, \cdot \rangle_{\Gamma_j}$ here stands for the duality pairing in $H^{\frac{1}{2}}(\Gamma_j)$. The following $L^2_\mu(\Omega)$-orthogonal Helmholtz-Hodge decomposition holds true:

$$L^2(\Omega) = \operatorname{grad} (H_0^1(\Omega)) \oplus \frac{\mu_\#}{\mu} \operatorname{curl} (H^1(\Omega) \cap L_0^2(\Omega)) \oplus H_\tau(\Omega). \qquad (9)$$

Furthermore, letting

$$H_\Gamma^1(\Omega) := \{v \in H^1(\Omega), v|_{\Gamma_0} \equiv 0 | \exists (\gamma_j) \in \mathbb{R}^{\beta_2}, v|_{\Gamma_j} \equiv \gamma_j \forall j \in \{1, \dots, \beta_2\}\}, \qquad (10)$$
it can be noticed that

$$ \mathrm{grad}(H_{\Gamma}^{1}(\Omega)) = \mathrm{grad}(H_{0}^{1}(\Omega)) \oplus H_{\tau}(\Omega). \quad (11) $$

Remark, also, that $\mathrm{grad}(H_{\Gamma}^{1}(\Omega)) \subset H_{0}(\mathrm{curl}^{0}; \Omega)$.

### 2.3.2 2nd Helmholtz-Hodge decomposition

Let us now consider the harmonic space $H_n(\Omega) := H(\mathrm{curl}^0; \Omega) \cap H_0(\mathrm{div}_\mu^0; \Omega)$ which, by (4), satisfies $H_n(\Omega) \cong \mathfrak{S}_0^2$. In particular, the space $H_n(\Omega)$ has dimension $\beta_1$, and it can be proved that vector fields $w \in H_n(\Omega)$ are entirely characterized by the data of $((\mu w)|_{\Sigma_i} \cdot n_{\Sigma_i}, 1)_{\Sigma_i \in \{1, \dots, \beta_1\}}$, where $\langle \cdot, \cdot \rangle_{\Sigma_i}$ here stands for the duality pairing in $H^{\frac{1}{2}}(\Sigma_i)$. The following $L_\mu^2(\Omega)$-orthogonal Helmholtz-Hodge decomposition holds true:

$$ L^2(\Omega) = \mathrm{grad}(H^1(\Omega) \cap L_0^2(\Omega)) \oplus \frac{\mu_\#}{\mu} \mathrm{curl}(H^1(\Omega) \cap H_0(\mathrm{curl}; \Omega)) \oplus H_n(\Omega). \quad (12) $$

Furthermore, letting

$$ H_{\Sigma}^{1}(\hat{\Omega}) := \{v \in H^{1}(\hat{\Omega}) \cap L_{0}^{2}(\hat{\Omega}) \mid \exists (\sigma_{i}) \in \mathbb{R}^{\beta_{1}}, \|v\|_{\Sigma_{i}} = \sigma_{i} \quad \forall i \in \{1, \dots, \beta_{1}\}\}, \quad (13) $$

it can be noticed that

$$ \mathrm{gr\v{a}d}(H_{\Sigma}^{1}(\hat{\Omega})) = \mathrm{grad}(H^{1}(\Omega) \cap L_{0}^{2}(\Omega)) \oplus H_{n}(\Omega). \quad (14) $$

We remind the reader that, for $v \in H^1(\hat{\Omega})$, $\mathrm{gr\v{a}d} v$ is the continuation to $L^2(\Omega)$ of $\mathrm{grad} v \in L^2(\hat{\Omega})$, and that the jump $\|\cdot\|_{\Sigma_i}$ is defined in (8). Note also that $\mathrm{gr\v{a}d}(H_{\Sigma}^1(\hat{\Omega})) \subset H(\mathrm{curl}^0; \Omega)$.

## 3 Discrete setting

From now on, we assume that the domain $\Omega \subset \mathbb{R}^3$ is a (Lipschitz) polyhedron.

### 3.1 Polyhedral discretizations

We consider discretizations $\mathcal{D} := (\mathcal{T}, \mathcal{F})$ of $\Omega \subset \mathbb{R}^3$ in the sense of [21, Def. 1.4].

The set $\mathcal{T}$ is a finite collection of disjoint open Lipschitz polyhedra $T$ (the mesh cells), which is assumed to form a partition of the domain, that is $\bar{\Omega} = \cup_{T \in \mathcal{T}} \bar{T}$. For all $T \in \mathcal{T}$, we let $h_T := \max_{x,y \in \bar{T}} |x - y|$ denote the diameter of the cell $T$. We also let $h_{\mathcal{T}}$ be s.t. $h_{\mathcal{T}|T} := h_T$ for all $T \in \mathcal{T}$, and we define the mesh size by $h_{\mathcal{D}} := \max_{T \in \mathcal{T}} h_T$. In turn, the set $\mathcal{F}$ is a finite collection of disjoint connected subsets of $\bar{\Omega}$ (the mesh faces) such that, for all $F \in \mathcal{F}$,

(i) $F$ is a relatively open, Lipschitz polygonal subset of an affine hyperplane, and

(ii) either there are two distinct mesh cells $T^+$, $T^-$ in $\mathcal{T}$ s.t. $\bar{F} \subseteq \partial T^+ \cap \partial T^-$ ($F$ is then an interface), or there is a mesh cell $T \in \mathcal{T}$ s.t. $\bar{F} \subseteq \partial T \cap \Gamma$ ($F$ is then a boundary face).

The set of mesh faces is assumed to form a partition of the mesh skeleton, that is to satisfy $\cup_{T \in \mathcal{T}} \partial T = \cup_{F \in \mathcal{F}} \bar{F}$. For all $F \in \mathcal{F}$, we let $h_F := \max_{x,y \in \bar{F}} |x - y|$ denote the diameter of the face $F$. Interfaces are collected in the set $\mathcal{F}^\circ$, whereas boundary faces are collected in the set $\mathcal{F}^\partial$. For all $T \in \mathcal{T}$, we denote by $\mathcal{F}_T$ the subset of $\mathcal{F}$ which collects the mesh faces lying on the boundary of $T$, so that $\partial T = \cup_{F \in \mathcal{F}_T} \bar{F}$. For all $T \in \mathcal{T}$, consistently with our notation so far, we let $n_{\partial T} : \partial T \to \mathbb{R}^3$ denote the (almost everywhere defined) unit vector field normal to $\partial T$, pointing outward from $T$.
For all $F \in \mathcal{F}_T$, we also let $n_{T,F} := n_{\partial T|_F}$ be the (constant) unit vector normal to the hyperplane containing $F$, and pointing outward from $T$. For all $F \in \mathcal{F}$, we define $n_F$ as the (constant) unit vector normal to $F$ such that either $n_F := n_{T^+,F}$ if $F \subset \partial T^+ \cap \partial T^- \in \mathcal{F}^\circ$, or $n_F := n_{T,F}(= n_{|F})$ if $F \subset \partial T \cap \Gamma \in \mathcal{F}^\partial$. For further use, we also let, for all $T \in \mathcal{T}$ and all $F \in \mathcal{F}_T$, $\varepsilon_{T,F} \in \{-1, 1\}$ be such that $\varepsilon_{T,F} := n_{T,F} \cdot n_F$. Finally, for any $X \in \mathcal{F} \cup \mathcal{T}$, we let $x_X \in \mathbb{R}^3$ be some point inside $X$.

Whenever $\beta_1 > 0$, for all $i \in \{1, \dots, \beta_1\}$, we assume that there exists a subset $\mathcal{F}_{\Sigma_i}^\circ$ of $\mathcal{F}^\circ$ such that $\bar{\Sigma}_i = \bigcup_{F \in \mathcal{F}_{\Sigma_i}^\circ} \bar{F}$, and for which $n_F = n_{\Sigma_i|_F}$ for all $F \in \mathcal{F}_{\Sigma_i}^\circ$. This requirement ensures that we do associate the tag “+” to the side of $\hat{\Omega}$ (with respect to $\Sigma_i$) for which $n_{\Sigma_i}$ is outward, consistently with our assumption from Section 2.1. Notice that, since the cutting surfaces are piecewise planar, the set $\hat{\Omega}$ is indeed pseudo-Lipschitz. We also let $\mathcal{F}^\circ := \bigcup_{i=1}^{\beta_1} \mathcal{F}_{\Sigma_i}^\circ$. When $\beta_2 > 0$, for all $j \in \{0, \dots, \beta_2\}$, we let $\mathcal{F}_{\Gamma_j}^\partial$ denote the subset of $\mathcal{F}^\partial$ such that $\Gamma_j = \bigcup_{F \in \mathcal{F}_{\Gamma_j}^\partial} \bar{F}$.

At the discrete level, the parameter $\mu : \Omega \to [\mu_b, \mu_\#]$ introduced in (7) is assumed to be piecewise constant over the partition $\mathcal{T}$ of the domain $\Omega$, and we let

$$ \mu_b \le \mu_T := \mu|_T \le \mu_\# \quad \forall T \in \mathcal{T}. \qquad (15) $$

The role of this assumption will be made clear in Sections 4.2.1 and 4.2.2 (right before Theorems 7 and 12).

When studying asymptotic properties with respect to the mesh size, one has to adopt a measure of regularity for refined sequences of discretizations. Let us thus consider $(\mathcal{D}_n)_{n \in \mathbb{N}}$, a sequence of polyhedral discretizations $\mathcal{D}_n := (\mathcal{T}_n, \mathcal{F}_n)$ such that $h_{\mathcal{D}_n}$ tends to zero as $n$ goes to infinity. We classically follow [21, Def. 1.9], in which regularity for refined mesh sequences is quantified through a uniform-in-$n$ parameter $\varrho \in (0, 1)$ (the mesh regularity parameter). In a nutshell, it is assumed that, for all $n \in \mathbb{N}$, there exists a matching tetrahedral subtessellation of $\mathcal{T}_n$, (i) which is uniformly-in-$n$ shape-regular, and (ii) whose elements have a diameter that is uniformly-in-$n$ comparable to the diameter of the mesh cell in $\mathcal{T}_n$ they belong to. In what follows, we write $a \lesssim b$ (resp. $a \gtrsim b$) in place of $a \le Cb$ (resp. $a \ge Cb$), if $C > 0$ only depends on $\Omega$, on the mesh regularity parameter $\varrho$, and (if need be) on the underlying polynomial degree, but is independent of both $n$ (and thus $h_{\mathcal{D}_n}$) and $\mu$. When $a \lesssim b \lesssim a$, we simply write $a \approx b$. In particular, for regular sequences $(\mathcal{D}_n)_{n \in \mathbb{N}}$ of discretizations, for all $n \in \mathbb{N}$ and $T \in \mathcal{T}_n$, there holds $\mathrm{card}(\mathcal{F}_T) \le 1$, as well as $h_T \lesssim h_F \le h_T$ for all $F \in \mathcal{F}_T$ (cf. [21, Lem. 1.12]). Also, we assume that, for all $n \in \mathbb{N}$ and $X \in \mathcal{F}_n \cup \mathcal{T}_n$, the point $x_X$ is the center of an $X$-inscribed disk/ball of radius $h_X \le r_X \le h_X$. This last assumption can always be satisfied for regular sequences of discretizations.

## 3.2 Polynomial spaces

For $\ell \in \mathbb{N}$ and $m \in \{2, 3\}$, we let $\mathcal{P}_m^\ell$ denote the linear space of $m$-variate polynomials of total degree at most $\ell$, with the convention that $\mathcal{P}_m^0$ is identified to $\mathbb{R}$, and that $\mathcal{P}_m^{-1} := \{0\}$. For any $X \in \mathcal{F} \cup \mathcal{T}$, we let $\mathcal{P}^\ell(X)$ be the linear space spanned by the restrictions to $X$ of the polynomials in $\mathcal{P}_m^\ell$. For $X$ of Hausdorff dimension $m \in \{2, 3\}$, $\mathcal{P}^\ell(X)$ is isomorphic to $\mathcal{P}_m^\ell$ (cf. [21, Prop. 1.23]). We let $\pi_{\mathcal{P},X}^\ell$ denote the $L^2(X)$-orthogonal projector onto $\mathcal{P}^\ell(X)$. For convenience, we also set $\mathcal{P}^\ell(X) := \mathcal{P}^\ell(X)^m$, that is, $\mathcal{P}^\ell(F) = \mathcal{P}^\ell(F)^2$ for all $F \in \mathcal{F}$, and $\mathcal{P}^\ell(T) = \mathcal{P}^\ell(T)^3$ for all $T \in \mathcal{T}$. We finally define $\pi_{\mathcal{P},X}^\ell$ as the $L^2(X)$-orthogonal projector onto $\mathcal{P}^\ell(X)$.

For any $T \in \mathcal{T}$ and $\ell \in \mathbb{N}$, we define $\mathcal{G}^\ell(T) := \mathrm{grad}(\mathcal{P}^{\ell+1}(T))$, and its Koszul complement $\mathcal{G}^{c,\ell}(T) := \mathcal{P}^{\ell-1}(T) \times (x - x_T)$, both subspaces of $\mathcal{P}^\ell(T)$. The (non $L^2(T)$-orthogonal) polynomial decomposition below is a by-product of the homotopy formula (cf. [3, Thm. 7.1]):

$$ \mathcal{P}^\ell(T) = \mathcal{G}^\ell(T) \oplus \mathcal{G}^{c,\ell}(T). \qquad (16) $$
In addition, by exactness of the polynomial de Rham complex, letting $\mathcal{R}^\ell(T) := \text{curl}(\mathcal{P}^{\ell+1}(T))$, the differential mapping $\text{curl} : \mathcal{G}^{c,\ell}(T) \to \mathcal{R}^{\ell-1}(T)$ is an isomorphism.

For any $F \in \mathcal{F}$ now, we let $H_F$ be the affine hyperplane containing $F$, that we orient according to the normal $n_F$. For any function $w : F \to \mathbb{R}$, we let $\text{grad}_F w : F \to \mathbb{R}^2$ denote the (tangential) gradient of $w$. We also let $\text{rot}_F w : F \to \mathbb{R}^2$ be such that $\text{rot}_F w := (\text{grad}_F w)^\perp$, where $z^\perp$ is defined as the rotation of angle $-\pi/2$ of $z$ in the oriented hyperplane $H_F$. For $\ell \in \mathbb{N}$, we define $\mathcal{R}^\ell(F) := \text{rot}_F(\mathcal{P}^{\ell+1}(F))$, and its Koszul complement $\mathcal{R}^{c,\ell}(F) := \mathcal{P}^{\ell-1}(F)(x-x_F)$, both subspaces of $\mathcal{P}^\ell(F)$, where, for $x \in F$, $(x-x_F) \subset H_F$ is identified to its (two-dimensional) tangential counterpart. By the homotopy formula, in that case again, the following (non $L^2(F)$-orthogonal) polynomial decomposition holds true:

$$ \mathcal{P}^{\ell}(F) = \mathcal{R}^{\ell}(F) \oplus \mathcal{R}^{c,\ell}(F). \qquad (17) $$

For all $T \in \mathcal{T}$ and $F \in \mathcal{F}_T$, since for any $v : T \to \mathbb{R}$ we have $n_F \times ((\text{grad } v)|_F \times n_F) = \text{grad}_F(v|_F)$ (identifying $n_F \times ((\text{grad } v)|_F \times n_F)$ to its two-dimensional proxy), there holds

$$ \mathcal{G}^{\ell}(T)|_{F} \times n_{F} = \mathcal{R}^{\ell}(F), \qquad (18) $$

where vectors in $\mathcal{G}^\ell(T)|_F \times n_F$ are also identified to their tangential counterparts. The above identity (18) provides a characterization for the rotated tangential traces of curl-free vector polynomials. This characterization is instrumental in the design of stable and optimally consistent trimmed HHO methods (cf. [33, Rmk. 10] for further insight).

Last, we introduce broken versions of the 3-variate polynomial spaces $\mathcal{P}^\ell$ and $\mathcal{P}^{\ell'}$, namely $\mathcal{P}^\ell(\mathcal{T}) := \{v \in L^2(\Omega) \mid v|_T \in \mathcal{P}^\ell(T) \forall T \in \mathcal{T}\}$, and $\mathcal{P}^{\ell'}(\mathcal{T}) := \{v \in L^2(\Omega) \mid v|_T \in \mathcal{P}^{\ell'}(T) \forall T \in \mathcal{T}\}$. We classically define on $\mathcal{P}^\ell(\mathcal{T})$ the broken gradient operator $\text{grad}_\mathcal{T}$, and on $\mathcal{P}^{\ell'}(\mathcal{T})$ the broken rotational operator $\text{curl}_\mathcal{T}$. We finally let $\pi_{\mathcal{P},\mathcal{T}}^\ell$ (resp. $\pi_{\mathcal{P},\mathcal{T}}^{\ell'}$) denote the $L^2(\Omega)$-orthogonal (resp. $L^2(\Omega)$-orthogonal) projector onto $\mathcal{P}^\ell(\mathcal{T})$ (resp. $\mathcal{P}^{\ell'}(\mathcal{T})$).

## 3.3 Hybrid spaces

We introduce hybrid counterparts of the spaces $H(\text{curl}; \Omega)$ and $H^1(\hat{\Omega})$ (as well as of relevant subspaces thereof). The first discrete spaces will serve for the discretization of the magnetic variable, whereas the second will be related to the pressure-like variable (Lagrange multiplier). Henceforth, let $k \in \mathbb{N}^*$ denote a given polynomial degree.

### 3.3.1 $H(\text{curl})$-like hybrid spaces

**Global spaces:** We define the following discrete counterpart of the space $H(\text{curl}; \Omega)$:

$$ \underline{U}_{\mathcal{D}}^{k} := \left\{ \underline{v}_{\mathcal{D}} := \left( (\underline{v}_{T})_{T \in \mathcal{T}}, (\underline{v}_{F,\tau})_{F \in \mathcal{F}} \right) : \begin{array}{ll} \underline{v}_{T} \in \mathcal{P}^{k}(T) & \forall T \in \mathcal{T} \\ \underline{v}_{F,\tau} \in \mathcal{Q}^{k}(F) & \forall F \in \mathcal{F} \end{array} \right\}, \qquad (19) $$

where the (possibly trimmed) polynomial space $\mathcal{Q}^k(F)$ shall satisfy

$$ \mathcal{R}^k(F) \subseteq \mathcal{Q}^k(F) \subseteq \mathcal{P}^k(F). \qquad (20) $$

The requirement $\mathcal{R}^k(F) \subseteq \mathcal{Q}^k(F)$ is stability-related (cf. [33, Rmk. 10] for further insight). The space $\mathcal{Q}^k(F)$ may also satisfy, whenever explicitly precised, the additional assumption

$$ \mathcal{P}^{k-1}(F) \subset \mathcal{Q}^k(F). \qquad (21) $$
We let $\pi_{Q,F}^k$ denote the $L^2(F)$-orthogonal projector onto $\mathcal{Q}^k(F)$. In (19), $v_{F,\tau}$ stands for the rotated tangential trace of the (magnetic) variable. We also introduce the subspace

$$ \underline{\mathcal{U}}_{\mathcal{D},0}^k := \{ \underline{v}_{\mathcal{D}} \in \underline{\mathcal{U}}_{\mathcal{D}}^k | v_{F,\tau} \equiv \mathbf{0} \ \forall F \in \mathcal{F}^{\partial} \}, \quad (22) $$

discrete counterpart of $H_0(\text{curl}; \Omega)$.

**Local traits:** Given a cell $T \in \mathcal{T}$, we denote by $\underline{\mathcal{U}}_T^k$ the restriction of $\underline{\mathcal{U}}_{\mathcal{D}}^k$ to $T$, and by $\underline{v}_T := (v_T, (v_{F,\tau})_{F \in \mathcal{T}}) \in \underline{\mathcal{U}}_T^k$ the restriction of the generic element $\underline{v}_{\mathcal{D}} \in \underline{\mathcal{U}}_{\mathcal{D}}^k$. For $\underline{v}_{\mathcal{D}} \in \underline{\mathcal{U}}_{\mathcal{D}}^k$, we also let $v_{\mathcal{T}}$ be the broken polynomial vector field in $\mathcal{P}^k(\mathcal{T})$ such that $v_{\mathcal{T}|T} := v_T$ for all $T \in \mathcal{T}$. Let us now define an $\mathcal{H}(\text{curl})$-like hybrid semi-norm on $\underline{\mathcal{U}}_T^k$. We let, for all $\underline{v}_T \in \underline{\mathcal{U}}_T^k$,

$$ |\underline{v}_T|_{\text{curl},T}^2 := \| \text{curl } v_T \|_{0,T}^2 + \sum_{F \in \mathcal{F}_T} h_F^{-1} \| \pi_{Q,F}^k (v_{T|F} \times n_F) - v_{F,\tau} \|_{0,F}^2, \quad (23) $$

where $v_{T|F} \times n_F$ is identified to its tangential proxy. In turn, at the global level, we classically set, for all $\underline{v}_{\mathcal{D}} \in \underline{\mathcal{U}}_{\mathcal{D}}^k$, $|\underline{v}_{\mathcal{D}}|_{\text{curl},\mathcal{D}}^2 := \sum_{T \in \mathcal{T}} |\underline{v}_T|_{\text{curl},T}^2$. Remark that, whenever $|\underline{v}_{\mathcal{D}}|_{\text{curl},\mathcal{D}} = 0$ for some $\underline{v}_{\mathcal{D}} \in \underline{\mathcal{U}}_{\mathcal{D}}^k$, then, by (18) and (20), one has $v_{\mathcal{T}} \in \mathcal{H}(\text{curl}^0; \Omega) \cap \mathcal{G}^k(\mathcal{T})$. We then define the rotational reconstruction operator. For any $T \in \mathcal{T}$, we let $C_T^{k-1} := \underline{\mathcal{U}}_T^k \to \mathcal{P}^{k-1}(T)$ be the operator such that, for all $\underline{v}_T \in \underline{\mathcal{U}}_T^k$, $C_T^{k-1}(\underline{v}_T) \in \mathcal{P}^{k-1}(T)$ is the unique solution to

$$ (C_T^{k-1}(\underline{v}_T), z)_T = (v_T, \text{curl} z)_T - \sum_{F \in \mathcal{F}_T} \epsilon_{T,F} (v_{F,\tau}, n_F \times (z_{|F} \times n_F))_F \quad \forall z \in \mathcal{P}^{k-1}(T), \quad (24) $$

where $n_F \times (z_{|F} \times n_F)$ is identified to its tangential proxy. At the global level, we let $C_T^{k-1} := \underline{\mathcal{U}}_{\mathcal{D}}^k \to \mathcal{P}^{k-1}(\mathcal{T})$ be such that, for all $\underline{v}_{\mathcal{D}} \in \underline{\mathcal{U}}_{\mathcal{D}}^k$, $C_T^{k-1}(\underline{v}_{\mathcal{D}})_T := C_T^{k-1}(\underline{v}_T)$ for all $T \in \mathcal{T}$. It now remains to define an $\mathcal{H}(\text{curl})$-like hybrid stabilizer. For any $T \in \mathcal{T}$, we introduce the following symmetric, positive semi-definite bilinear form: for all $\underline{v}_T, \underline{w}_T \in \underline{\mathcal{U}}_T^k$,

$$ S_{\text{curl},T}(\underline{w}_T, \underline{v}_T) := \sum_{F \in \mathcal{F}_T} h_F^{-1} (\pi_{Q,F}^k (w_{T|F} \times n_F) - w_{F,\tau}, \pi_{Q,F}^k (v_{T|F} \times n_F) - v_{F,\tau})_F. \quad (25) $$

At the global level, we then classically let $S_{\text{curl},\mathcal{D}}(\underline{w}_{\mathcal{D}}, \underline{v}_{\mathcal{D}}) := \sum_{T \in \mathcal{T}} S_{\text{curl},T}(\underline{w}_T, \underline{v}_T)$ for all $\underline{v}_{\mathcal{D}}, \underline{w}_{\mathcal{D}} \in \underline{\mathcal{U}}_{\mathcal{D}}^k$. Given $T \in \mathcal{T}$, under the assumption that $\mathcal{Q}^k(F)$ (already satisfying (20)) additionally satisfies (21), it is an easy matter to prove that, for all $\underline{v}_T \in \underline{\mathcal{U}}_T^k$,

$$ \|C_T^{k-1}(\underline{v}_T)\|_{0,T}^2 + S_{\text{curl},T}(\underline{v}_T, \underline{v}_T) \approx |\underline{v}_T|_{\text{curl},T}^2. \quad (26) $$

**Approximation:** Let $q > 2$. Given $T \in \mathcal{T}$, we introduce the local reduction operator $\underline{\mathcal{I}}_T^k : \mathcal{H}(\text{curl}; T) \cap \mathcal{L}^q(T) \to \underline{\mathcal{U}}_T^k$ such that, for any $v \in \mathcal{H}(\text{curl}; T) \cap \mathcal{L}^q(T)$,

$$ \underline{\mathcal{I}}_T^k(v) := (\pi_{\mathcal{P},T}^k(v), (\pi_{Q,F}^k(v_{|F} \times n_F))_{F \in \mathcal{F}_T}), \quad (27) $$

where, as now standard, $v_{|F} \times n_F$ is identified to its (two-dimensional) tangential counterpart. Following [25, Sec. 5.1], note that the regularity $v \in \mathcal{H}(\text{curl}; T) \cap \mathcal{L}^q(T)$ is sufficient to give a (weak) meaning to the face polynomial projections in (27) (see also [33, Rmk. 11]). At the domain level, now, the global reduction operator $\underline{\mathcal{I}}_{\mathcal{D}}^k : \mathcal{H}(\text{curl}; \Omega) \cap \mathcal{L}^q(\Omega) \to \underline{\mathcal{U}}_{\mathcal{D}}^k$ is defined so that, for any $v \in \mathcal{H}(\text{curl}; \Omega) \cap \mathcal{L}^q(\Omega)$, $\underline{\mathcal{I}}_{\mathcal{D}}^k(v) := ((\pi_{\mathcal{P},T}^k(v_{|T}))_{T \in \mathcal{T}}, (\pi_{Q,F}^k(v_{|F} \times n_F))_{F \in \mathcal{F}})$. Remark that, since $v \in \mathcal{H}(\text{curl}; \Omega)$, the quantity $v_{|F} \times n_F$ is single-valued at interfaces $F \in \mathcal{F}^o$. Notice also that $\underline{\mathcal{I}}_{\mathcal{D}}^k(H_0(\text{curl}; \Omega) \cap \mathcal{L}^q(\Omega)) \subset \underline{\mathcal{U}}_{\mathcal{D},0}^k$. Given $T \in \mathcal{T}$, under the assumption that $\mathcal{Q}^k(F)$ (already
satisfying (20)) additionally satisfies (21), it is an easy matter to prove the following commutation property: for all $v \in \mathbf{H}(\text{curl}; T) \cap \mathbf{L}^q(T)$,

$$ (C_T^{k-1} \circ \underline{\mathbf{I}}_T^k)(v) = \pi_{\mathcal{P},T}^{k-1}(\text{curl } v). \qquad (28) $$

Finally, without this time the need for the additional assumption (21) on $\mathcal{Q}^k(F)$, the following (optimal) polynomial consistency result holds true: for all $\mathbf{p} \in \mathcal{P}^k(T)$,

$$ S_{\text{curl},T}(\underline{\mathbf{I}}_T^k(\mathbf{p}), \underline{v}_T) = 0 \quad \forall \underline{v}_T \in \underline{\mathbf{U}}_T^k. $$

### 3.3.2 $H^1$-like hybrid spaces

**Global spaces:** We define the following discrete counterpart of the space $H^1(\hat{\Omega})$:

$$ \underline{\hat{P}}_{\underline{\mathcal{D}}}^{k} := \left\{ \underline{q}_{\underline{\mathcal{D}}} := ((q_T)_{T \in \mathcal{T}}, (q_F)_{F \in \mathcal{F} \setminus \mathcal{F}_{\underline{\Sigma}}^{\circ}}, (q_F^{+}, q_F^{-})_{F \in \mathcal{F}_{\underline{\Sigma}}^{\circ}}) : \begin{array}{ll} q_T \in \mathcal{P}^{k-1}(T) & \forall T \in \mathcal{T} \\ q_F \in \mathcal{P}^{k}(F) & \forall F \in \mathcal{F} \setminus \mathcal{F}_{\underline{\Sigma}}^{\circ} \\ q_F^{+}, q_F^{-} \in \mathcal{P}^{k}(F) & \forall F \in \mathcal{F}_{\underline{\Sigma}}^{\circ} \end{array} \right\}, \qquad (29) $$

where we recall that $\mathcal{F}_{\underline{\Sigma}}^{\circ} = \bigcup_{i=1}^{\beta_1} \mathcal{F}_{\Sigma_i}^{\circ}$, with $\mathcal{F}_{\Sigma_i}^{\circ}$ subset of $\mathcal{F}^{\circ}$ such that $\overline{\Sigma_i} = \bigcup_{F \in \mathcal{F}_{\Sigma_i}^{\circ}} \overline{F}$. Here, for any $i \in \{1, \dots, \beta_1\}$ and $F \in \mathcal{F}_{\Sigma_i}^{\circ}$, $q_F^{+}$ and $q_F^{-}$ stand for the traces on $F$ of the pressure-like variable defined, respectively, from the “+” and “-” sides of $\hat{\Omega}$ with respect to the cutting surface $\Sigma_i$. Mirroring (8), for $\underline{q}_{\underline{\mathcal{D}}} \in \underline{\hat{P}}_{\underline{\mathcal{D}}}^{k}$, we define its jump through $F \in \mathcal{F}_{\Sigma_i}^{\circ}$ by $\llbracket \underline{q}_{\underline{\mathcal{D}}} \rrbracket_F := q_F^{+} - q_F^{-}$. We then introduce the following subspace of $\underline{\hat{P}}_{\underline{\mathcal{D}}}^{k}$:

$$ \underline{\hat{P}}_{\underline{\mathcal{D}}, \Sigma}^{b,k} := \{ \underline{q}_{\underline{\mathcal{D}}} \in \underline{\hat{P}}_{\underline{\mathcal{D}}}^{k} : \exists (\sigma_i) \in \mathbb{R}^{\beta_1}, \llbracket \underline{q}_{\underline{\mathcal{D}}} \rrbracket_F \equiv \sigma_i \quad \forall F \in \mathcal{F}_{\Sigma_i}^{\circ}, i \in \{1, \dots, \beta_1\} \}, \qquad (30) $$

which is, modulo the zero-mean condition over $\hat{\Omega}$, a hybrid counterpart of the space $H_{\Sigma}^1(\hat{\Omega})$ defined in (13). Next, we define the following discrete counterpart of the space $H^1(\Omega)$:

$$ \underline{P}_{\underline{\mathcal{D}}}^{k} := \left\{ \underline{q}_{\underline{\mathcal{D}}} := ((q_T)_{T \in \mathcal{T}}, (q_F)_{F \in \mathcal{F}}) : \begin{array}{ll} q_T \in \mathcal{P}^{k-1}(T) & \forall T \in \mathcal{T} \\ q_F \in \mathcal{P}^{k}(F) & \forall F \in \mathcal{F} \end{array} \right\}. \qquad (31) $$

Remark that if $\beta_1 = 0$, then $\underline{\hat{P}}_{\underline{\mathcal{D}}, \Sigma}^{b,k} = \underline{\hat{P}}_{\underline{\mathcal{D}}}^{k} = \underline{P}_{\underline{\mathcal{D}}}^{k}$. Whenever $\beta_1 > 0$, the space $\underline{P}_{\underline{\mathcal{D}}}^{k}$ is isomorphic to the subspace of $\underline{\hat{P}}_{\underline{\mathcal{D}}, \Sigma}^{b,k}$ such that $\llbracket \underline{q}_{\underline{\mathcal{D}}} \rrbracket_F \equiv 0$ for all $F \in \mathcal{F}_{\underline{\Sigma}}^{\circ}$. We may now introduce the following subspaces of $\underline{P}_{\underline{\mathcal{D}}}^{k}$:

$$ \underline{P}_{\underline{\mathcal{D}},0}^{k} := \{ \underline{q}_{\underline{\mathcal{D}}} \in \underline{P}_{\underline{\mathcal{D}}}^{k} : q_F \equiv 0 \quad \forall F \in \mathcal{F}_{\Gamma_0}^{\partial} \}, \qquad (32) $$

as well as

$$ \underline{P}_{\underline{\mathcal{D}},\Gamma}^{k} := \{ \underline{q}_{\underline{\mathcal{D}}} \in \underline{P}_{\underline{\mathcal{D}},0}^{k} : \exists (\gamma_j) \in \mathbb{R}^{\beta_2}, q_F \equiv \gamma_j \quad \forall F \in \mathcal{F}_{\Gamma_j}^{\partial}, j \in \{1, \dots, \beta_2\} \}. \qquad (33) $$

The space $\underline{P}_{\underline{\mathcal{D}},\Gamma}^{k}$ is the hybrid counterpart of the space $H_{\Gamma}^1(\Omega)$ in (10). Remark that, whenever $\beta_2 = 0$, then $\underline{P}_{\underline{\mathcal{D}},\Gamma}^{k} = \underline{P}_{\underline{\mathcal{D}},0}^{k}$ with $\underline{P}_{\underline{\mathcal{D}},0}^{k}$ being, in that case, the hybrid counterpart of $H_0^1(\Omega)$.
**Local traits:** Given a cell $T \in \mathcal{T}$, we denote by $\underline{P}_T^k$ the restriction of $\hat{P}_\mathcal{D}^k$ or $\underline{P}_\mathcal{D}^k$ to $T$, and by $\underline{q}_T := (q_T, (q_{F,T})_{F \in \mathcal{F}_T}) \in \underline{P}_T^k$ the restriction of the generic element $\underline{q}_\mathcal{D}$ in $\hat{P}_\mathcal{D}^k$ or $\underline{P}_\mathcal{D}^k$. For $\underline{q}_\mathcal{D} \in \underline{P}_\mathcal{D}^k$, $q_{F,T} := q_F$ for all $F \in \mathcal{F}_T$, whereas for $\underline{q}_\mathcal{D} \in \underline{P}_\mathcal{D}^k$,

* if $F \in \mathcal{F}_T \setminus \mathcal{F}_\Sigma^\circ$, then $q_{F,T} := q_F$;

* if $F \in \mathcal{F}_T \cap \mathcal{F}_\Sigma^\circ$ with $F \subset \partial T^+ \cap \partial T^-$, then $q_{F,T} := q_F^+$ if $T = T^+$, whereas $q_{F,T} := q_F^-$ if $T = T^-$.

For $\underline{q}_\mathcal{D}$ in $\hat{P}_\mathcal{D}^k$ or $\underline{P}_\mathcal{D}^k$, we also let $q_T$ be the broken polynomial function in $\mathcal{P}^{k-1}(\mathcal{T})$ such that $q_{T|T} := q_T$ for all $T \in \mathcal{T}$. Let us now define an $(L^2)$-scaled $H^1$-like hybrid semi-norm on $\underline{P}_T^k$. We let, for all $\underline{q}_T \in \underline{P}_T^k$,

$$|\underline{q}_T|^2_{\text{grad},T} := h_T^2 \| \text{grad } q_T \|_{0,T}^2 + \sum_{F \in \mathcal{F}_T} h_F \| q_T|_F - q_{F,T} \|_{0,F}^2. \quad (34)$$

In turn, at the global level, for all $\underline{q}_\mathcal{D}$ in $\hat{P}_\mathcal{D}^k$ or $\underline{P}_\mathcal{D}^k$, we set $|\underline{q}_\mathcal{D}|_{\text{grad},\mathcal{D}}^2 := \sum_{T \in \mathcal{T}} |\underline{q}_T|^2_{\text{grad},T}$. Remark that, whenever $|\underline{q}_\mathcal{D}|_{\text{grad},\mathcal{D}} = 0$ for some $\underline{q}_\mathcal{D}$ in $\hat{P}_\mathcal{D}^k$ or $\underline{P}_\mathcal{D}^k$, then there exists $c \in \mathbb{R}$ such that $\underline{q}_\mathcal{D} \cong c \underline{1}_\mathcal{D}$, with $\underline{1}_\mathcal{D}$ the vector of $\underline{P}_\mathcal{D}^k$ whose components are all equal to one. We now define the gradient reconstruction operator. For any $T \in \mathcal{T}$, we let $G_T^k : \underline{P}_T^k \to \mathcal{P}^k(T)$ be the operator such that, for all $\underline{q}_T \in \underline{P}_T^k$, $G_T^k(\underline{q}_T) \in \mathcal{P}^k(T)$ is the unique solution to

$$(G_T^k(\underline{q}_T), z)_T = -(q_T, \operatorname{div} z)_T + \sum_{F \in \mathcal{F}_T} (q_{F,T}, z|_F \cdot n_{T,F})_F \quad \forall z \in \mathcal{P}^k(T). \quad (35)$$

At the global level, we let $G_T^k : \{\hat{P}_\mathcal{D}^k, \underline{P}_\mathcal{D}^k\} \to \mathcal{P}^k(\mathcal{T})$ be such that, for all $\underline{q}_\mathcal{D}$ in $\hat{P}_\mathcal{D}^k$ or $\underline{P}_\mathcal{D}^k$, $G_T^k(\underline{q}_\mathcal{D})|_T := G_T^k(\underline{q}_T)$ for all $T \in \mathcal{T}$. It now remains to define an $(L^2)$-scaled $H^1$-like hybrid stabilizer. For any $T \in \mathcal{T}$, we introduce the following symmetric, positive semi-definite bilinear form: for all $\underline{q}_T, \underline{r}_T \in \underline{P}_T^k$,

$$S_{\text{grad},T}(\underline{r}_T, \underline{q}_T) := \sum_{F \in \mathcal{F}_T} h_F (r_T|_F - r_{F,T}, q_T|_F - q_{F,T})_F. \quad (36)$$

At the global level, we then set $S_{\text{grad},\mathcal{D}}(\underline{r}_\mathcal{D}, \underline{q}_\mathcal{D}) := \sum_{T \in \mathcal{T}} S_{\text{grad},T}(\underline{r}_T, \underline{q}_T)$ for all $\underline{q}_\mathcal{D}, \underline{r}_\mathcal{D}$ in $\hat{P}_\mathcal{D}^k$ or $\underline{P}_\mathcal{D}^k$. Given $T \in \mathcal{T}$, it is an easy matter to prove that, for all $\underline{q}_T \in \underline{P}_T^k$,

$$h_T^2 \|G_T^k(\underline{q}_T)\|_{0,T}^2 + S_{\text{grad},T}(\underline{q}_T, \underline{q}_T) \approx |\underline{q}_T|^2_{\text{grad},T}. \quad (37)$$

Besides, under the assumption that $T$ belongs to a tetrahedral mesh sequence, by [11, Lem. 1] (and adapting the arguments from the proof of [14, Lem. 3.2]), the following important equivalence holds true: for all $\underline{q}_T \in \underline{P}_T^k$,

$$h_T^2 \|G_T^k(\underline{q}_T)\|_{0,T}^2 \approx |\underline{q}_T|^2_{\text{grad},T}. \quad (38)$$

**Approximation:** Given $T \in \mathcal{T}$, we introduce the local reduction operator $\underline{I}_T^k : H^1(T) \to \underline{P}_T^k$ such that, for any $q \in H^1(T)$,

$$\underline{I}_T^k(q) := (\pi_{\mathcal{P},T}^{k-1}(q), (\pi_{\mathcal{P},F}^k(q|_F))_{F \in \mathcal{F}_T}). \quad (39)$$

At the domain level, now, the global reduction operator $\underline{I}_\mathcal{D}^k : \{H^1(\hat{\Omega}), H^1(\Omega)\} \to \{\hat{P}_\mathcal{D}^k, \underline{P}_\mathcal{D}^k\}$ is defined in the following way:
* for any $q \in H^1(\hat{\Omega})$,

$$ \underline{I}_{\mathcal{D}}^{k}(q) := \left( (\pi_{\mathcal{P},T}^{k-1}(q|_{T}))_{T \in \mathcal{T}}, (\pi_{\mathcal{P},F}^{k}(q|_{F}))_{F \in \mathcal{F} \setminus \mathcal{F}_{\Sigma}^{\circ}}, (\pi_{\mathcal{P},F}^{k}(q_{|F}^{+}), \pi_{\mathcal{P},F}^{k}(q_{|F}^{-}))_{F \in \mathcal{F}_{\Sigma}^{\circ}} \right), $$

where, for $F \in \mathcal{F}_{\Sigma}^{\circ}$ such that $F \subset \partial T^{+} \cap \partial T^{-}$, we have let $q^{+} := q|_{T^{+}}$ and $q^{-} := q|_{T^{-}}$;

* for any $q \in H^1(\hat{\Omega})$, $\underline{I}_{\mathcal{D}}^{k}(q) := ((\pi_{\mathcal{P},T}^{k-1}(q|_{T}))_{T \in \mathcal{T}}, (\pi_{\mathcal{P},F}^{k}(q|_{F}))_{F \in \mathcal{F}})$.

Remark that, for $F \in \mathcal{F}^{\circ}$ with $F \subset \partial T^{+} \cap \partial T^{-}$, whenever $q \in H^1(\overline{T^{+}} \cup \overline{T^{-}})$, then the quantity $q|_{F}$ is single-valued. Note, in addition, that $\underline{I}_{\mathcal{D}}^{k}(H_{\Sigma}^{1}(\hat{\Omega})) \subset \underline{P}_{\mathcal{D},\Sigma}^{b,k}$ and that $\underline{I}_{\mathcal{D}}^{k}(H_{\Gamma}^{1}(\hat{\Omega})) \subset \underline{P}_{\mathcal{D},\Gamma}^{k}$. Finally, given $T \in \mathcal{T}$, the following commutation property is valid: for all $q \in H^1(T)$,

$$ (\underline{G}_{T}^{k} \circ \underline{I}_{T}^{k})(q) = \pi_{\mathcal{P},T}^{k}(\mathbf{grad} \, q). \qquad (40) $$

Also, the following polynomial consistency result holds true: for all $p \in \mathcal{P}^{k-1}(T)$,

$$ S_{\text{grad},T}(\underline{I}_{T}^{k}(p), \underline{q}_{T}) = 0 \quad \forall \underline{q}_{T} \in \underline{P}_{T}^{k}. $$

# 4 HHO methods

We consider two models of magnetostatics: (i) the (first-order) field formulation, endowed with normal boundary conditions, and (ii) the (second-order) vector potential formulation, endowed with tangential boundary conditions. This last model was already considered in [11, Sec. 3.2], but therein for trivial topology only. These two problems are equivalent. They correspond to the primal instance (that is, based on the primal de Rham complex) of the model problems (5) and (6). The dual instance could be handled similarly.

## 4.1 Magnetostatics models

Recall the definitions (10) and (13) of the spaces $H_{\Gamma}^{1}(\Omega)$ and $H_{\Sigma}^{1}(\hat{\Omega})$. Henceforth, the (variable) coefficient $\mu$, satisfying (15), will denote the magnetic permeability of the medium.

### 4.1.1 Field formulation

Let $j : \Omega \to \mathbb{R}^3$ be a given electric current density satisfying the following compatibility conditions:

$$ \mathrm{div} \, j = 0 \text{ in } \Omega, \quad \langle j|_{\Gamma_j} \cdot n, 1 \rangle_{\Gamma_j} = 0 \text{ for all } j \in \{1, \dots, \beta_2\}. \qquad (41) $$

We seek the magnetic field $h : \Omega \to \mathbb{R}^3$ such that

$$ \mathrm{curl} \, h = j \quad \text{in } \Omega, \qquad (42a) $$

$$ \mathrm{div} \, b = 0 \quad \text{in } \Omega, \qquad (42b) $$

$$ b|_{\Gamma} \cdot n = 0 \quad \text{on } \Gamma, \qquad (42c) $$

$$ \langle b|_{\Sigma_i} \cdot n|_{\Sigma_i}, 1 \rangle_{\Sigma_i} = 0 \quad \text{for all } i \in \{1, \dots, \beta_1\}, \qquad (42d) $$

with constitutive law $b = \mu h$, where $b : \Omega \to \mathbb{R}^3$ represents the magnetic induction.

In what follows, we assume that $j \in L^2(\Omega)$. By (11), remark that (41) is equivalent to

$$ (j, z)_{\Omega} = 0 \quad \forall z \in \mathrm{grad} \, (H_{\Gamma}^{1}(\Omega)). \qquad (43) $$
In the spirit of Kikuchi [31], we consider the following equivalent (cf. Remark 2 below) weak form for Problem (42): Find $(h, p) \in \mathbf{H}(\text{curl}; \Omega) \times H_{\Sigma}^{1}(\hat{\Omega})$ s.t.

$$A(h, v) + B(v, p) = (j, \text{curl } v)_{\Omega} \quad \forall v \in \mathbf{H}(\text{curl}; \Omega), \quad (44a)$$

$$-B(h, q) = 0 \qquad \forall q \in H_{\Sigma}^{1}(\hat{\Omega}), \quad (44b)$$

with bilinear forms $A: \mathbf{H}(\text{curl}; \Omega) \times \mathbf{H}(\text{curl}; \Omega) \to \mathbb{R}$, $B: \mathbf{H}(\text{curl}; \Omega) \times H^{1}(\hat{\Omega}) \to \mathbb{R}$ set to

$$A(w, v) := (\text{curl } w, \text{curl } v)_{\Omega}, \quad B(w, q) := (\mu w, \text{gr}\acute{a}\d{d} q)_{\Omega}, \quad (45)$$

where we recall that, for $q \in H^1(\hat{\Omega})$, $\text{gr}\acute{a}\d{d} q$ is the continuation to $L^2(\Omega)$ of $\text{grad } q \in L^2(\hat{\Omega})$. The pressure-like variable $p$ is the Lagrange multiplier of the constraints (42b)–(42d) on the magnetic induction $b$. Testing (44a) with $v = \text{gr}\acute{a}\d{d} p \in \text{gr}\acute{a}\d{d}(H_{\Sigma}^{1}(\hat{\Omega})) \subset \mathbf{H}(\text{curl}^0; \Omega)$, and recalling that $p \in L_0^2(\hat{\Omega})$, one can actually infer that $p = 0$ in $\hat{\Omega}$ (recall that $\hat{\Omega}$ is assumed to be connected). The well-posedness of Problem (44) is a direct consequence of the second Weber inequality, more precisely of [33, Rmk. 6] with $\eta := \mu$ (combined with (44b) and (14)).

**Remark 2 (Weak-strong equivalence).** Whereas it is clear that any solution to Problem (42) also solves Problem (44), the converse is less straightforward, in particular when it comes to retrieving (42a). To do so, one has to use the first Helmholtz–Hodge decomposition (9) (with $\mu \leftarrow 1$ therein), combine it to (11), then recall the compatibility condition (43), so as to infer that $(\text{curl } h - j, z)_{\Omega} = 0$ for all $z \in L^2(\Omega)$, thereby yielding (42a).

**Remark 3 (Regularity theory).** Recall that $\Omega$ is a Lipschitz polyhedron. When the coefficient $\mu$ is globally smooth (which, under (15), amounts to assuming that $\mu$ is constant in $\Omega$), it is known (cf. [18, Thm. 2] and [2, Prop. 3.7]) that the solution to Problem (42) (and more generally to Problem (5)) satisfies $h \in H^s(\Omega)$ for some $s > \frac{1}{2}$, with $s \ge 1$ if $\Omega$ is convex. When the coefficient $\mu$ is piecewise smooth (here, under (15), piecewise constant), however, the solution $h$ to Problem (42) (and, more generally, to Problem (5)) is only known (cf. [19, 30, 8]) to belong to $H^s(\Omega)$ for some $s > 0$.

### 4.1.2 Vector potential formulation

Let $j: \Omega \to \mathbb{R}^3$ be a given electric current density complying with (41). We seek the magnetic vector potential $a: \Omega \to \mathbb{R}^3$ such that

$$\text{curl}(\mu^{-1} \text{curl } a) = j \quad \text{in } \Omega, \quad (46a)$$

$$\text{div } a = 0 \quad \text{in } \Omega, \quad (46b)$$

$$a|_{\Gamma} \times n = 0 \quad \text{on } \Gamma, \quad (46c)$$

$$\langle a|_{\Gamma_j} \cdot n, 1 \rangle_{\Gamma_j} = 0 \quad \text{for all } j \in \{1, \dots, \beta_2\}. \quad (46d)$$

In what follows, we assume that $j \in L^2(\Omega)$. Recall that (41) is equivalent to (43). We consider the following weak form for Problem (46): Find $(a, p) \in H_0(\text{curl}; \Omega) \times H_{\Gamma}^{1}(\hat{\Omega})$ s.t.

$$A(a, v) + B(v, p) = (j, v)_{\Omega} \quad \forall v \in H_0(\text{curl}; \Omega), \quad (47a)$$

$$-B(a, q) = 0 \qquad \forall q \in H_{\Gamma}^{1}(\hat{\Omega}), \quad (47b)$$

with bilinear forms $A: \mathbf{H}(\text{curl}; \Omega) \times \mathbf{H}(\text{curl}; \Omega) \to \mathbb{R}$, $B: \mathbf{H}(\text{curl}; \Omega) \times H^1(\hat{\Omega}) \to \mathbb{R}$ set to

$$A(w, v) := (\mu^{-1} \text{curl } w, \text{curl } v)_{\Omega}, \quad B(w, q) := (w, \text{grad } q)_{\Omega}. \quad (48)$$
The pressure-like variable $p$ is the Lagrange multiplier of the constraints (46b) and (46d) on the vector potential $a$. Testing (47a) with $v = \text{grad} p \in \text{grad}(H_Γ^1(\Omega)) \subset H_0(\text{curl}^0; \Omega)$, using (43), and recalling the definition (10), one can actually infer that $p = 0$ in $\Omega$. Consequently, Problems (46) and (47) are equivalent (that is, weak-strong equivalence holds true). The well-posedness of Problem (47) is a direct consequence of the first Weber inequality, more precisely of [33, Rmk. 3] with $\eta := 1$ (combined with (47b) and (11)).

**Remark 4 (Equivalence with Problem (42)).** Whereas it is clear that, for any solution $a$ to Problem (46), setting $b := \text{curl } a$, then $h = \mu^{-1}b$ is solution to Problem (42), the converse is less obvious. It is a consequence of the fact that the equations (42b)-(42d) imply, by the second Helmholtz-Hodge decomposition (12), the existence of $a \in H_0(\text{curl}; \Omega)$ such that $h = \mu^{-1} \text{curl } a$. The latter vector potential is non-unique, which is the reason why one has to further impose a gauge condition. In our case, we impose the so-called Coulomb gauge, which amounts to enforcing (46b) and (46d).

**Remark 5 (Regularity theory).** Given that $h = \mu^{-1} \text{curl } a$ is solution to Problem (42), we infer from Remark 3 that $\mu^{-1} \text{curl } a \in H^r(\Omega)$, with $r > \frac{1}{2}$ if $\mu$ is globally smooth, or $r > 0$ if $\mu$ is piecewise smooth. In turn, since $a \in H_0(\text{curl}; \Omega) \cap H(\text{div}^0; \Omega)$, there holds $a \in H^s(\Omega)$ for some $s > \frac{1}{2}$. In particular, if $\mu$ is globally constant, both $\mu^{-1} \text{curl } a$ and $a$ belong to $H^t(\Omega)$ for some $t > \frac{1}{2}$, with $t \ge 1$ if $\Omega$ is convex. In the general case, the solution to Problem (46) (and, more generally, to Problem (6)) satisfies $a \in H^s(\Omega)$ for some $s > \frac{1}{2}$, and $\mu^{-1} \text{curl } a \in H^r(\Omega)$ for some $r > 0$.

## 4.2 HHO schemes

The HHO methods we devise are directly built upon the weak formulations (44) and (47). At the discrete level, we also heavily leverage the fact that the continuous Lagrange multipliers are equal to zero.

### 4.2.1 Field formulation

Recall the definition (19) of the space $\underline{U}_\mathcal{D}^k$, hybrid counterpart of $H(\text{curl}; \Omega)$, and assume that, for all $F \in \mathcal{F}$, the space $Q^k(F)$ is given by

$$Q^k(F) := \mathcal{R}^k(F) \oplus \mathcal{R}^{c,k-1}(F), \quad (49)$$

in such a way that (20) is indeed fulfilled, and that (21) additionally holds true. Recall the definition (30) of the space $\hat{P}_{\mathcal{D},\Sigma}^{b,k}$. We introduce the following hybrid counterpart of $H_\Sigma^1(\hat{\Omega})$:

$$\underline{\hat{P}}_{\mathcal{D},\Sigma}^k := \{ \underline{q}_\mathcal{D} \in \underline{\hat{P}}_{\mathcal{D},\Sigma}^{b,k} : \int_{\hat{\Omega}} q_T = 0 \}. \quad (50)$$

Remark that $I_\mathcal{D}^k(H_\Sigma^1(\hat{\Omega})) \subset \underline{\hat{P}}_{\mathcal{D},\Sigma}^k$, since functions in $H_\Sigma^1(\hat{\Omega})$ also belong to $L_0^2(\hat{\Omega})$.

Given $T \in \mathcal{T}$, we introduce the local bilinear forms $A_T: \underline{U}_T^k \times \underline{U}_T^k \to \mathbb{R}$, $B_T: \underline{U}_T^k \times \underline{P}_T^k \to \mathbb{R}$, and $S_T: \underline{P}_T^k \times \underline{P}_T^k \to \mathbb{R}$ defined by: for all $\underline{v}_T, \underline{w}_T \in \underline{U}_T^k$, and $\underline{q}_T, \underline{r}_T \in \underline{P}_T^k$,

$$A_T(\underline{w}_T, \underline{v}_T) := (\underline{C}_T^{k-1}(\underline{w}_T), \underline{C}_T^{k-1}(\underline{v}_T))_T + S_{\text{curl},T}(\underline{w}_T, \underline{v}_T), \quad (51a)$$

$$B_T(\underline{w}_T, \underline{q}_T) := \mu_T(\underline{w}_T, \underline{G}_T^k(\underline{q}_T))_T, \quad (51b)$$

$$S_T(\underline{r}_T, \underline{q}_T) := \mu_T^2 \underline{h}_T^2 (\text{grad } r_T, \text{grad } q_T)_T + \mu_T^2 S_{\text{grad},T}(\underline{r}_T, \underline{q}_T), \quad (51c)$$

where the local stabilizers $S_{\text{curl},T}$ and $S_{\text{grad},T}$ are respectively defined by (25) and (36). The scaling $\mu_T^2$ (also $\underline{h}_T^2$) in (51c) is chosen so as to restore consistency between the physical units. As far as the
bilinear contribution $S_T$ is concerned, its role will be discussed in Remark 10 below. As now standard, the global bilinear forms $A_D: \underline{U}_D^k \times \underline{U}_D^k \to \mathbb{R}$, $B_D: \underline{U}_D^k \times \hat{P}_D^k \to \mathbb{R}$, and $S_D: \hat{P}_D^k \times \hat{P}_D^k \to \mathbb{R}$ are assembled summing the local contributions. The discrete HHO problem reads: find $(\underline{h}_D, \underline{p}_D) \in \underline{U}_D^k \times \hat{P}_D^k$ s.t.

$$A_D(\underline{h}_D, \underline{v}_D) + B_D(\underline{v}_D, \underline{p}_D) = (\underline{j}, C_T^{k-1}(\underline{v}_D))_{\Omega} \quad \forall \underline{v}_D \in \underline{U}_D^k, \quad (52a)$$

$$-B_D(\underline{h}_D, \underline{q}_D) + S_D(\underline{p}_D, \underline{q}_D) = 0 \quad \forall \underline{q}_D \in \hat{P}_{\underline{D}, \Sigma}^k. \quad (52b)$$

We now establish well-posedness for Problem (52).

**Lemma 6 (Well-posedness).** For all $\underline{v}_D \in \underline{U}_D^k$ and $\underline{q}_D \in \hat{P}_D^k$, there holds:

$$|\underline{v}_D|^2_{\text{curl},D} + \sum_{T \in \mathcal{T}} \mu_T^2 |\underline{q}_T|^2_{\text{grad},T} \lesssim A_D(\underline{v}_D, \underline{v}_D) + S_D(\underline{q}_D, \underline{q}_D). \quad (53)$$

As a consequence, Problem (52) is well-posed.

*Proof.* Let $\underline{v}_T \in \underline{U}_T^k$ and $\underline{q}_T \in \hat{P}_T^k$. Starting from definition (51a), and invoking (26) (which is valid under assumption (49)), there holds

$$A_T(\underline{v}_T, \underline{v}_T) = \|C_T^{k-1}(\underline{v}_T)\|_{0,T}^2 + S_{\text{curl},T}(\underline{v}_T, \underline{v}_T) \approx |\underline{v}_T|^2_{\text{curl},T}. \quad (54)$$

In the same vein, by (51c), $S_T(\underline{q}_T, \underline{q}_T) = \mu_T^2 |\underline{q}_T|^2_{\text{grad},T}$. Summing over $T \in \mathcal{T}$ yields (53).

Now, since the linear system corresponding to Problem (52) is square, proving well-posedness is equivalent to proving the uniqueness of its solution. Let us show that, if $j \equiv 0$, then necessarily $(\underline{h}_D, \underline{p}_D) = (\underline{0}_D, \underline{0}_D)$. First, remark that the map $|\cdot|_{\text{grad},D}$ (cf. (34)) is a norm on $\hat{P}_D^k$ (cf. (50)), as a by-product of the discrete zero-mean condition over $\hat{\Omega}$. Second, let us test Problem (52) (for $j \equiv 0$) with $\underline{v}_D = \underline{h}_D$ and $\underline{q}_D = \underline{p}_D$. Summing (52a) and (52b), and leveraging (53), we directly infer that $|\underline{h}_D|_{\text{curl},D} = 0$ and $|\underline{p}_D|_{\text{grad},D} = 0$. From the second, since $\underline{p}_D \in \hat{P}_D^k$, we deduce that $\underline{p}_D = \underline{0}_D$. Plugging this value into (52b), we get $B_D(\underline{h}_D, \underline{q}_D) = (\mu \underline{h}_T, G_T^k(\underline{q}_D))_{\Omega} = 0$ for all $\underline{q}_D \in \hat{P}_D^k$. Since, for any $q \in H_{\Sigma}^1(\hat{\Omega})$, $\underline{l}_D^k(q) \in \hat{P}_D^k$, letting $\underline{q}_D = \underline{l}_D^k(q)$ in the latter relation, and invoking the local commutation property (40), together with $\mu \underline{h}_T \in \mathcal{P}^k(\mathcal{T})$ (as a by-product of (15)), we infer

$$(\mu \underline{h}_T, z)_{\Omega} = 0 \quad \forall z \in \text{gr}\acute{d}(H_{\Sigma}^1(\hat{\Omega})). \quad (55)$$

The conclusion then follows from the second hybrid Weber inequality of [33, Rmk. 16] (with $\eta := \mu$) which, combined to (14), yields $\|\mu^{\frac{1}{2}} \underline{h}_T\|_{0,\Omega} \lesssim \mu_{\#}^{\frac{1}{2}} |\underline{h}_D|_{\text{curl},D}(\lesssim \mu_{\#}^{\frac{1}{2}} \|\underline{j}\|_{0,\Omega})$. Since $|\underline{h}_D|_{\text{curl},D} = 0$, we infer that $\underline{h}_D = \underline{0}_D$, thereby concluding the proof. $\square$

Remark that the validity of the discrete orthogonality condition (55), on which crucially hinges the stability of the method, requires assumption (15) on the magnetic permeability to be fulfilled. Last, we prove an error estimate in energy-norm for the solution to Problem (52).

**Theorem 7 (Energy-error estimate).** Assume that $\underline{h} \in H(\text{curl}; \Omega)$ solution to Problem (44) further satisfies $\underline{h} \in H^s(\mathcal{T})$ for some $s \in (1, k+1]$, and $\text{curl} \underline{h} \in H^r(\mathcal{T})$ for some $r \in (s-1, k]$. Then, the following estimate holds true:

$$|\underline{h}_D - \underline{I}_D^k(\underline{h})|_{\text{curl},D}^2 + \sum_{T \in \mathcal{T}} \mu_T^2 |\underline{p}_T|^2_{\text{grad},T} \lesssim \sum_{T \in \mathcal{T}} (h_T^{2r} |\text{curl} \underline{h}|_{r,T}^2 + h_T^{2(s-1)} |\underline{h}|_{s,T}^2), \quad (56)$$

where $(\underline{h}_D, \underline{p}_D) \in \underline{U}_D^k \times \hat{P}_D^k$ is the unique solution to Problem (52).
*Proof.* First, by Sobolev embedding, remark that since $h \in H^s(\mathcal{T})$ for some $s > 1$ ($s > 0$ would be enough), then $h \in L^q(\Omega)$ for some $q > 2$, and one can thus give a meaning to $\underline{\mathcal{I}}_{\mathcal{D}}^k(h) \in \underline{\mathcal{U}}_{\mathcal{D}}^k$. Let us set $|(\underline{\mathbf{w}}_{\mathcal{D}}, \underline{\mathbf{r}}_{\mathcal{D}})|_{\mathcal{D}}^2 := A_{\mathcal{D}}(\underline{\mathbf{w}}_{\mathcal{D}}, \underline{\mathbf{w}}_{\mathcal{D}}) + S_{\mathcal{D}}(\underline{\mathbf{r}}_{\mathcal{D}}, \underline{\mathbf{r}}_{\mathcal{D}})$. The starting point of the proof is the following simple inequality: for all $(\underline{\mathbf{w}}_{\mathcal{D}}, \underline{\mathbf{r}}_{\mathcal{D}}) \in \underline{\mathcal{U}}_{\mathcal{D}}^k \times \underline{\hat{\mathcal{P}}}_{\mathcal{D}, \Sigma}^k$,

$$ |(\underline{\mathbf{w}}_{\mathcal{D}}, \underline{\mathbf{r}}_{\mathcal{D}})|_{\mathcal{D}} \leq \max_{\substack{(\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}}) \in \underline{\mathcal{U}}_{\mathcal{D}}^k \times \underline{\hat{\mathcal{P}}}_{\mathcal{D}, \Sigma}^k \\ |(\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}})|_{\mathcal{D}} = 1}} A_{\mathcal{D}}(\underline{\mathbf{w}}_{\mathcal{D}}, \underline{\mathbf{v}}_{\mathcal{D}}) + B_{\mathcal{D}}(\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{r}}_{\mathcal{D}}) - B_{\mathcal{D}}(\underline{\mathbf{w}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}}) + S_{\mathcal{D}}(\underline{\mathbf{r}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}}). $$

Letting then $\underline{\mathbf{w}}_{\mathcal{D}} = (\underline{\mathbf{h}}_{\mathcal{D}} - \underline{\mathcal{I}}_{\mathcal{D}}^k(h)) \in \underline{\mathcal{U}}_{\mathcal{D}}^k$ and $\underline{\mathbf{r}}_{\mathcal{D}} = (\underline{\mathbf{p}}_{\mathcal{D}} - \underline{\mathcal{I}}_{\mathcal{D}}^k(p)) = \underline{\mathbf{p}}_{\mathcal{D}} \in \underline{\hat{\mathcal{P}}}_{\mathcal{D}, \Sigma}^k$ (it is recalled that $p \equiv 0$), developing the expression into the max, and using Problem (52), yields

$$ \left( |\underline{\mathbf{h}}_{\mathcal{D}} - \underline{\mathcal{I}}_{\mathcal{D}}^k(h)|_{\text{curl}, \mathcal{D}}^2 + \sum_{T \in \mathcal{T}} \mu_T^2 |\underline{\mathbf{p}}_T|_{\text{grad}, T}^2 \right)^{1/2} \lesssim \max_{\substack{(\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}}) \in \underline{\mathcal{U}}_{\mathcal{D}}^k \times \underline{\hat{\mathcal{P}}}_{\mathcal{D}, \Sigma}^k \\ |(\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}})|_{\mathcal{D}} = 1}} \mathcal{E}_{\mathcal{D}}((\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}})), $$

where we have also utilized the estimate (53) in the left-hand side, and set $\mathcal{E}_{\mathcal{D}}((\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}})) := (j, C_{\mathcal{T}}^{k-1}(\underline{\mathbf{v}}_{\mathcal{D}}))_{\Omega} - A_{\mathcal{D}}(\underline{\mathcal{I}}_{\mathcal{D}}^k(h), \underline{\mathbf{v}}_{\mathcal{D}}) + B_{\mathcal{D}}(\underline{\mathcal{I}}_{\mathcal{D}}^k(h), \underline{\mathbf{q}}_{\mathcal{D}})$. Let us define

$$ \begin{aligned} \mathfrak{T}_1 &:= (j, C_{\mathcal{T}}^{k-1}(\underline{\mathbf{v}}_{\mathcal{D}}))_{\Omega} - A_{\mathcal{D}}(\underline{\mathcal{I}}_{\mathcal{D}}^k(h), \underline{\mathbf{v}}_{\mathcal{D}}) \\ &= (\text{curl } h - (C_{\mathcal{T}}^{k-1} \circ \underline{\mathcal{I}}_{\mathcal{D}}^k)(h), C_{\mathcal{T}}^{k-1}(\underline{\mathbf{v}}_{\mathcal{D}}))_{\Omega} - S_{\text{curl}, \mathcal{D}}(\underline{\mathcal{I}}_{\mathcal{D}}^k(h), \underline{\mathbf{v}}_{\mathcal{D}}), \end{aligned} $$

and

$$ \mathfrak{T}_2 := B_{\mathcal{D}}(\underline{\mathcal{I}}_{\mathcal{D}}^k(h), \underline{\mathbf{q}}_{\mathcal{D}}) = (\pi_{\mathcal{P}, \mathcal{T}}^k(\mu h), \mathcal{G}_{\mathcal{T}}^k(\underline{\mathbf{q}}_{\mathcal{D}}))_{\Omega}, $$

where we have invoked (42a) to replace $j$ by $\text{curl } h$ in $\mathfrak{T}_1$, and used the fact that $\mu$ is piecewise constant in $\mathfrak{T}_2$. For $\mathfrak{T}_1$, using (i) the local commutation property (28) along with the fact that $\text{curl } h \in H^r(\mathcal{T})$, and (ii) the fact that $h \in H^s(\mathcal{T})$ for some $s > 1$ ($s > \frac{1}{2}$ would be enough), so that the full trace of $h$ can be given a meaning in $L^2(F)^3$ on each interface (on both sides) and boundary face, we infer from the fractional approximation results of [14, Lem. 2.5] that

$$ \mathfrak{T}_1 \lesssim \left( \sum_{T \in \mathcal{T}} (h_T^{2r} |\text{curl } h|_{r,T}^2 + h_T^{2(s-1)} |\mathbf{h}|_{s,T}^2) \right)^{1/2} A_{\mathcal{D}}(\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{v}}_{\mathcal{D}})^{1/2}. $$

For $\mathfrak{T}_2$ now, using the definition (35) of the (local) gradient reconstruction, and integrating by parts the volumetric term therein, we infer that

$$ \mathfrak{T}_2 = (\mu h, \text{grad}_{\mathcal{T}} q_{\mathcal{T}})_{\Omega} + \sum_{T \in \mathcal{T}} \sum_{F \in \mathcal{F}_T} (\pi_{\mathcal{P},T}^k (\mu_T h|_T)_F \cdot n_{T,F}, q_{F,T} - q_T|_F)_F, $$

where we have invoked the fact that $\text{grad}_{\mathcal{T}} q_{\mathcal{T}} \in \mathcal{P}^{k-2}(\mathcal{T})$ to remove the projector in the first term. Integrating by parts (cell by cell) the first term in the right-hand side, and combining $h \in H_0(\text{div}_{\mu}^0; \Omega)$ (owing to (42b) and (42c)) and $\underline{\mathbf{q}}_{\mathcal{D}} \in \underline{\hat{\mathcal{P}}}_{\mathcal{D}, \Sigma}$ (cf. (50) and (30)), yields

$$ \mathfrak{T}_2 = \sum_{T \in \mathcal{T}} \sum_{F \in \mathcal{F}_T} \varepsilon_{T,F} \left( (\pi_{\mathcal{P},T}^k (\mu_T h|_T)_F - (\mu h)_F) \cdot n_F, q_{F,T} - q_T|_F \right)_F + \sum_{i=1}^{\beta_1} \sum_{F \in \mathcal{F}_{\Sigma_i}^c} \sigma_i ((\mu h)_F \cdot n_F, 1)_F, $$
where we have additionally leveraged that $\mu h \in H^s(\mathcal{T})$ for some $s > 1$ ($s > \frac{1}{2}$ would be sufficient) to give a standard meaning to its normal traces. Now, since, for any $i \in \{1, \dots, \beta_1\}$, there holds $\sum_{F \in \mathcal{F}_{\Sigma_i}} \sigma_i((\mu h)|_F \cdot n_F, 1)_F = \sigma_i((\mu h)|_{\Sigma_i} \cdot n_{\Sigma_i}, 1)_{\Sigma_i}$, we infer by (42d) that

$$ \mathfrak{I}_2 = \sum_{T \in \mathcal{T}} \sum_{F \in \mathcal{F}_T} \varepsilon_{T,F} ((\pi_{\mathcal{P},T}^k (\mu_T h|_T)|_F - (\mu h)|_F) \cdot n_F, q_{F,T} - q_{T|F})_F. $$

The approximation results of [14, Lem. 2.5], along with standard arguments, then show that

$$ \mathfrak{I}_2 \lesssim \left( \sum_{T \in \mathcal{T}} h_T^{2(s-1)} |h|_{s,T}^2 \right)^{1/2} \left( \sum_{T \in \mathcal{T}} \mu_T^2 S_{\text{grad},T}(\underline{q}_T, \underline{q}_T) \right)^{1/2}. $$

Since $\sum_{T \in \mathcal{T}} \mu_T^2 S_{\text{grad},T}(\underline{q}_T, \underline{q}_T) \le S_{\mathcal{D}}(\underline{q}_{\mathcal{D}}, \underline{q}_{\mathcal{D}})$, gathering the estimates on $\mathfrak{I}_1$ and $\mathfrak{I}_2$, and using that $|\underline{v}_{\mathcal{D}}, \underline{q}_{\mathcal{D}}|_{\mathcal{D}} = 1$, we conclude the proof. $\square$

From the error estimate (56), we infer that the HHO scheme converges with optimal rate $h_{\mathcal{D}}^{k-1}$ when $h \in H^{k+1}(\mathcal{T})$. Before concluding this section, a few remarks are in order.

**Remark 8 (Low-regularity solutions and error measure).** In order to yield convergence, the energy error estimate from Theorem 7 essentially requires that $h \in H^s(\mathcal{T})$ for some $s > 1$. By Remark 3, this regularity assumption on the solution is not met in general, even for a globally constant coefficient $\mu$ on a convex domain. Recovering convergence for any $s > 0$ and $r \ge 0$ would require to consider a weaker error measure, morally on $h_T C_T^{k-1}$ instead of $C_T^{k-1}$. This is in line with the fact that, for first-order models, the only meaningful error measure is the $L^2$-error on $h$, which formally scales like $h_T C_T^{k-1}$. For $s \in (0, \frac{1}{2}]$, the analysis would also necessitate to give a meaning to the face terms (adapting, for example, the arguments from [25, Sec. 5]). We do not pursue herein further in this direction.

**Remark 9 (Variant of the scheme).** Since the model at hand is first-order, following [11, Sec. 3.1], it is also possible to define the discrete rotational operator as the broken rotational $\text{curl}_T$, instead of $C_T^{k-1}$. In practice, the new HHO scheme is obtained replacing in Problem (52) every occurrence of $C_T^{k-1}$ by $\text{curl}_T$. The main advantage of doing so is that one can then relax the assumption (21), and consider in place of (49) the leaner space $\mathcal{Q}^k(F) := \mathcal{R}^k(F)$. This choice is actually sufficient (and necessary, thus optimal) for ensuring stability (cf. [33, Rmk. 10]). Redefining $A_T(\underline{w}_T, \underline{v}_T)$, there now holds (in place of (54)) that $A_T(\underline{v}_T, \underline{v}_T) = |\underline{v}_T|^2_{\text{curl},T}$, and this identity holds true irrespectively of (21). There is a price to pay for this complexity reduction. The commutation property (28) is no longer valid in that case, and the proxy for $\text{curl} h$ is then given by $\text{curl}_T(\pi_{\mathcal{P},T}^k(h))$ in lieu of $\pi_{\mathcal{P},T}^{k-1}(\text{curl} h)$. As supported by our numerical experiments (cf. Section 5.5), for singular solutions, this alternative proxy might unfortunately retain insufficient approximability properties.

**Remark 10 (About stability).** Assume $\mathcal{T}$ is a member of a tetrahedral mesh sequence. Then, as already pointed out in [11, Rmk. 3.1], the bilinear contribution $S_{\mathcal{D}}(\underline{p}_{\mathcal{D}}, \underline{q}_{\mathcal{D}})$ can be removed from Problem (52) without compromising its stability. In the tetrahedral case, (38) indeed yields $\sum_{T \in \mathcal{T}} \mu_T^2 |q_T|^2_{\text{grad},T} \lesssim \sum_{T \in \mathcal{T}} \mu_T^2 h_T^2 \|G_T^k(q_T)\|_{0,T}^2$ for all $\underline{q}_{\mathcal{D}} \in \hat{\mathcal{P}}_{\mathcal{D}}^k$. Since the discrete saddle-point system enforces a control on the latter right-hand side (cf. [11, Lem. 3.6]), this actually grants stability. When $\mathcal{T}$ belongs to a general mesh sequence, (37) indicates that $\sum_{T \in \mathcal{T}} \mu_T^2 S_{\text{grad},T}(\underline{q}_T, \underline{q}_T)$ needs also be controlled to ensure stability. This motivates the adding of the bilinear term $S_{\mathcal{D}}(\underline{p}_{\mathcal{D}}, \underline{q}_{\mathcal{D}})$, embedding $\sum_{T \in \mathcal{T}} \mu_T^2 S_{\text{grad},T}(\underline{p}_T, \underline{q}_T)$, to Problem (52). Remarkably, the continuous Lagrange multiplier being zero, any choice of $S_{\mathcal{D}}(\underline{p}_{\mathcal{D}}, \underline{q}_{\mathcal{D}})$ is consistent. We here define $S_{\mathcal{D}}$ from (51c), which
also embeds a volumetric contribution, with the idea of improving stability for the discrete problem. This strategy may yield slightly more accurate results in some situations (cf. [11, Rmk. 3.2]). Other choices of the same nature can be made for $S_D$, see e.g. [11, Eq. (3.4c)]. The present remark extends to Problem (58) (vector potential case).

### 4.2.2 Vector potential formulation

Recall the definition (22) of $\underline{U}_{\underline{D},0}^k$, hybrid counterpart of $H_0(\text{curl}; \Omega)$, and assume that the space $Q^k(F)$ is still given by (49) for all $F \in \mathcal{F}$. Recall, in addition, the definition (33) of $\underline{P}_{\underline{D},\Gamma}^k$, hybrid counterpart of $H_1^1(\Omega)$.

Given $T \in \mathcal{T}$, we introduce the local bilinear forms $A_T : \underline{U}_T^k \times \underline{U}_T^k \to \mathbb{R}$, $B_T : \underline{U}_T^k \times \underline{P}_T^k \to \mathbb{R}$, and $S_T : \underline{P}_T^k \times \underline{P}_T^k \to \mathbb{R}$ defined by: for all $\underline{v}_T, \underline{w}_T \in \underline{U}_T^k$, and $\underline{q}_T, \underline{r}_T \in \underline{P}_T^k$,

$$A_T(\underline{w}_T, \underline{v}_T) := \mu_T^{-1}(\underline{C}_T^{k-1}(\underline{w}_T), \underline{C}_T^{k-1}(\underline{v}_T))_T + \mu_T^{-1}S_{\text{curl},T}(\underline{w}_T, \underline{v}_T), \quad (57a)$$

$$B_T(\underline{w}_T, \underline{q}_T) := (\underline{w}_T, \underline{G}_T^k(\underline{q}_T))_T, \quad (57b)$$

$$S_T(\underline{r}_T, \underline{q}_T) := \mu_T h_T^2 (\text{grad } \underline{r}_T, \text{grad } \underline{q}_T)_T + \mu_T S_{\text{grad},T}(\underline{r}_T, \underline{q}_T), \quad (57c)$$

where the local stabilizers $S_{\text{curl},T}$ and $S_{\text{grad},T}$ are respectively defined by (25) and (36). Here as well, the scaling $\mu_T$ (also $h_T^2$) in (57c) is chosen so as to restore consistency between the physical units. In turn, the global bilinear forms $A_D : \underline{U}_D^k \times \underline{U}_D^k \to \mathbb{R}$, $B_D : \underline{U}_D^k \times \underline{P}_D^k \to \mathbb{R}$, and $S_D : \underline{P}_D^k \times \underline{P}_D^k \to \mathbb{R}$ are assembled by summing the local contributions. The discrete HHO problem reads: Find $(\underline{a}_D, \underline{p}_D) \in \underline{U}_{\underline{D},0}^k \times \underline{P}_{\underline{D},\Gamma}^k$ s.t.

$$A_D(\underline{a}_D, \underline{v}_D) + B_D(\underline{v}_D, \underline{p}_D) = (\underline{j}, \underline{v}_T)_\Omega \quad \forall \underline{v}_D \in \underline{U}_{\underline{D},0}^k, \quad (58a)$$

$$-B_D(\underline{a}_D, \underline{q}_D) + S_D(\underline{p}_D, \underline{q}_D) = 0 \quad \forall \underline{q}_D \in \underline{P}_{\underline{D},\Gamma}^k. \quad (58b)$$

Let us establish well-posedness for Problem (58).

**Lemma 11 (Well-posedness).** For all $\underline{v}_D \in \underline{U}_D^k$ and $\underline{q}_D \in \underline{P}_D^k$, there holds:

$$\sum_{T \in \mathcal{T}} \mu_T^{-1} |\underline{v}_T|^2_{\text{curl},T} + \sum_{T \in \mathcal{T}} \mu_T |\underline{q}_T|^2_{\text{grad},T} \lesssim A_D(\underline{v}_D, \underline{v}_D) + S_D(\underline{q}_D, \underline{q}_D). \quad (59)$$

As a consequence, Problem (58) is well-posed.

*Proof.* Let $\underline{v}_T \in \underline{U}_T^k$ and $\underline{q}_T \in \underline{P}_T^k$. Starting from definition (57a), and invoking (26) (which is valid under assumption (49)), there holds

$$A_T(\underline{v}_T, \underline{v}_T) = \mu_T^{-1} \|\underline{C}_T^{k-1}(\underline{v}_T)\|_{0,T}^2 + \mu_T^{-1} S_{\text{curl},T}(\underline{v}_T, \underline{v}_T) \approx \mu_T^{-1} |\underline{v}_T|^2_{\text{curl},T}. \quad (60)$$

Similarly, by (57c), $S_T(\underline{q}_T, \underline{q}_T) = \mu_T |\underline{q}_T|^2_{\text{grad},T}$. Summing over $T \in \mathcal{T}$ yields (59).

Now, since the linear system corresponding to Problem (58) is square, proving well-posedness is equivalent to proving the uniqueness of its solution. Let us show that, if $\underline{j} \equiv \underline{0}$, then necessarily $(\underline{a}_D, \underline{p}_D) = (\underline{0}_D, \underline{0}_D)$. First, remark that the map $|\cdot|_{\text{grad},D}$ (cf. (34)) is a norm on $\underline{P}_{\underline{D},\Gamma}^k$ (cf. (32) and (33)), as a by-product of the discrete zero boundary condition on $\Gamma_0$. Now, let us test Problem (58) (for $\underline{j} \equiv \underline{0}$) with $\underline{v}_D = \underline{a}_D$ and $\underline{q}_D = \underline{p}_D$. Summing (58a) and (58b), and leveraging (59), we directly infer that $|\underline{a}_D|_{\text{curl},D} = 0$ and $|\underline{p}_D|_{\text{grad},D} = 0$. From the second, since $\underline{p}_D \in \underline{P}_{\underline{D},\Gamma}^k$, we deduce that $\underline{p}_D = \underline{0}_D$. Plugging this value into (58b), we get $B_D(\underline{a}_D, \underline{q}_D) = (\underline{a}_T, \underline{G}_T^k(\underline{q}_D))_\Omega = 0$ for all $\underline{q}_D \in \underline{P}_{\underline{D},\Gamma}^k$. Since, for any $q \in H_1^1(\Omega)$, $\underline{I}_D^k(q) \in \underline{P}_{\underline{D},\Gamma}^k$, letting $\underline{q}_D = \underline{I}_D^k(q)$ in the latter
relation, and leveraging the local commutation property (40), along with $\mathbf{a}_{\mathcal{T}} \in \mathcal{P}^k(\mathcal{T})$, we infer $(\mathbf{a}_{\mathcal{T}}, \mathbf{z})_{\Omega} = 0$ for all $\mathbf{z} \in \text{grad}(H_{\Gamma}^1(\Omega))$. The conclusion then follows (recall that $\underline{\mathbf{a}}_{\mathcal{D}} \in \underline{U}_{\mathcal{D},0}^k$) from the first hybrid Weber inequality of [33, Rmk. 14] (with $\eta := 1$) which, combined to (11), yields $\|\mathbf{a}_{\mathcal{T}}\|_{0,\Omega} \lesssim |\underline{\mathbf{a}}_{\mathcal{D}}|_{\text{curl},\mathcal{D}}(\lesssim \mu_{\#}\|\mathbf{j}\|_{0,\Omega})$. Since $|\underline{\mathbf{a}}_{\mathcal{D}}|_{\text{curl},\mathcal{D}} = 0$, we infer $\underline{\mathbf{a}}_{\mathcal{D}} = \underline{\mathbf{0}}_{\mathcal{D}}$. $\square$

Remark that, contrary to what held true for the field formulation, the stability of Problem (58) is not conditional to assumption (15). In order to simplify the analysis below, we nonetheless stick to (15), and refer to Section 5.4 for a numerical assessment in the case of locally variable permeability. We prove an error estimate in energy-norm for the solution to Problem (58).

**Theorem 12 (Energy-error estimate).** Assume that the solution $\mathbf{a} \in H_0(\text{curl}; \Omega)$ to Problem (47) further satisfies $\mathbf{a} \in H^s(\mathcal{T})$ for some $s \in (1, k+1]$, and $\text{curl} \mathbf{a} \in H^r(\mathcal{T})$ for some $r \in (\max(\frac{1}{2}, s-1), k]$. Then, the following estimate holds true:

$$ \sum_{T \in \mathcal{T}} \mu_T^{-1} |\underline{\mathbf{a}}_T - \underline{\mathbf{I}}_T^k (\mathbf{a}_T)|_{\text{curl},T}^2 + \sum_{T \in \mathcal{T}} \mu_T |\underline{p}_T|_{\text{grad},T}^2 \lesssim \sum_{T \in \mathcal{T}} \mu_T^{-1} (h_T^{2r} |\text{curl } \mathbf{a}|_{r,T}^2 + h_T^{2(s-1)} |\mathbf{a}|_{s,T}^2), \quad (61) $$

where $(\underline{\mathbf{a}}_{\mathcal{D}}, \underline{p}_{\mathcal{D}}) \in \underline{U}_{\mathcal{D},0}^k \times \underline{P}_{\mathcal{D},\Gamma}^k$ is the unique solution to Problem (58).

*Proof.* First, by Sobolev embedding, remark that since $\mathbf{a} \in H^s(\mathcal{T})$ for some $s > 1$ ($s > 0$ would be enough), then $\mathbf{a} \in L^q(\Omega)$ for some $q > 2$, and one can thus give a meaning to $\underline{\mathbf{I}}_{\mathcal{D}}^k(\mathbf{a}) \in \underline{U}_{\mathcal{D},0}^k$. Let us set $|(\underline{\mathbf{w}}_{\mathcal{D}}, \underline{r}_{\mathcal{D}})|_{\mathcal{D}}^2 := A_{\mathcal{D}}(\underline{\mathbf{w}}_{\mathcal{D}}, \underline{\mathbf{w}}_{\mathcal{D}}) + S_{\mathcal{D}}(\underline{r}_{\mathcal{D}}, \underline{r}_{\mathcal{D}})$. Proceeding as in the proof of Theorem 7, we infer that

$$ \left( \sum_{T \in \mathcal{T}} \mu_T^{-1} |\underline{\mathbf{a}}_T - \underline{\mathbf{I}}_T^k (\mathbf{a}_T)|_{\text{curl},T}^2 + \sum_{T \in \mathcal{T}} \mu_T |\underline{p}_T|_{\text{grad},T}^2 \right)^{1/2} \lesssim \max_{\substack{(\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}}) \in \underline{U}_{\mathcal{D},0}^k \times \underline{P}_{\mathcal{D},\Gamma}^k \\ |(\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}})|_{\mathcal{D}} = 1}} \mathcal{E}_{\mathcal{D}}((\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}})), $$

where $\mathcal{E}_{\mathcal{D}}((\underline{\mathbf{v}}_{\mathcal{D}}, \underline{\mathbf{q}}_{\mathcal{D}})) := (\underline{j}, \underline{v}_{\mathcal{T}})_{\Omega} - A_{\mathcal{D}}(\underline{\mathbf{I}}_{\mathcal{D}}^k(\mathbf{a}), \underline{\mathbf{v}}_{\mathcal{D}}) + B_{\mathcal{D}}(\underline{\mathbf{I}}_{\mathcal{D}}^k(\mathbf{a}), \underline{\mathbf{q}}_{\mathcal{D}})$. Let us define

$$ \mathfrak{T}_1 := (\underline{j}, \underline{v}_{\mathcal{T}})_{\Omega} - A_{\mathcal{D}}(\underline{\mathbf{I}}_{\mathcal{D}}^k(\mathbf{a}), \underline{\mathbf{v}}_{\mathcal{D}}), \quad \mathfrak{T}_2 := B_{\mathcal{D}}(\underline{\mathbf{I}}_{\mathcal{D}}^k(\mathbf{a}), \underline{\mathbf{q}}_{\mathcal{D}}) = (\underline{\pi}_{\mathcal{P},\mathcal{T}}^k(\mathbf{a}), \underline{\mathbf{G}}_{\mathcal{T}}^k(\underline{\mathbf{q}}_{\mathcal{D}}))_{\Omega}. $$

For $\mathfrak{T}_1$, (i) replacing $\underline{j}$ by $\text{curl}(\mu^{-1} \text{curl } \mathbf{a})$ (leveraging (46a)), then performing a cell-by-cell integration by parts in $(\underline{j}, \underline{v}_{\mathcal{T}})_{\Omega}$, and (ii) using the definition (24) of the (local) rotational reconstruction operator for $A_{\mathcal{D}}(\underline{\mathbf{I}}_{\mathcal{D}}^k(\mathbf{a}), \underline{\mathbf{v}}_{\mathcal{D}})$, we infer that

$$ \mathfrak{T}_1 = (\mu^{-1}(\text{curl } \mathbf{a} - (\mathbf{C}_{\mathcal{T}}^{k-1} \circ \underline{\mathbf{I}}_{\mathcal{D}}^k)(\mathbf{a})), \text{curl}_{\mathcal{T}} \underline{v}_{\mathcal{T}})_{\Omega} - \sum_{T \in \mathcal{T}} \mu_T^{-1} S_{\text{curl},T} (\underline{\mathbf{I}}_T^k (\mathbf{a}_T)|_T, \underline{\mathbf{v}}_T) \\ + \sum_{T \in \mathcal{T}} \sum_{F \in \mathcal{F}_T} \mathcal{E}_{T,F} (\underline{\mathbf{v}}_T|_F \times \mathbf{n}_F - \underline{\mathbf{v}}_{F,\mathcal{T}}, \mathbf{n}_F \times (\mu^{-1} \text{curl } \mathbf{a} - \mu_T^{-1}(\mathbf{C}_T^{k-1} \circ \underline{\mathbf{I}}_T^k)(\mathbf{a}_T))|_F \times \mathbf{n}_F)_F, $$

where we have also utilized the fact that $\mu^{-1} \text{curl } \mathbf{a} \in H(\text{curl}; \Omega)$ with $\mu^{-1} \text{curl } \mathbf{a} \in H^r(\mathcal{T})$ for some $r > \frac{1}{2}$, together with the fact that $\underline{\mathbf{v}}_{\mathcal{D}} \in \underline{U}_{\mathcal{D},0}^k$. Let us first estimate the two first addends in $\mathfrak{T}_1$, whose sum we denote by $\mathfrak{T}_{1,1}$ ($\mathfrak{T}_{1,2}$ will denote the last addend of $\mathfrak{T}_1$). By (i) the local commutation property (28), along with the fact that $\text{curl } \mathbf{a} \in H^r(\mathcal{T})$, and (ii) the fact that $\mathbf{a} \in H^s(\mathcal{T})$ for some $s > 1$ ($s > \frac{1}{2}$ would be enough to give a meaning to the boundary terms), we get from the fractional approximation results of [14, Lem. 2.5] that

$$ \mathfrak{T}_{1,1} \lesssim \left( \sum_{T \in \mathcal{T}} \mu_T^{-1} (h_T^{2r} |\text{curl } \mathbf{a}|_{r,T}^2 + h_T^{2(s-1)} |\mathbf{a}|_{s,T}^2) \right)^{1/2} \left( \sum_{T \in \mathcal{T}} \mu_T^{-1} |\underline{\mathbf{v}}_T|_{\text{curl},T}^2 \right)^{1/2}. $$
In turn, using [33, Lem. 9] (yielding full control of tangential jumps), together with similar arguments, we infer that

$$ \mathfrak{I}_{1,2} \lesssim \left( \sum_{T \in \mathcal{T}} \mu_T^{-1} h_T^{2r} |\mathbf{curl} \mathbf{a}|_{r,T}^2 \right)^{1/2} \left( \sum_{T \in \mathcal{T}} \mu_T^{-1} |\underline{v}_T|^2_{\mathbf{curl},T} \right)^{1/2}. $$

By (60), we finally get $\mathfrak{I}_1 \lesssim (\sum_{T \in \mathcal{T}} \mu_T^{-1} (h_T^{2r} |\mathbf{curl} \mathbf{a}|_{r,T}^2 + h_T^{2(s-1)} |\mathbf{a}|_{s,T}^2))^{1/2} A_{\mathcal{D}}(\underline{v}_{\mathcal{D}}, \underline{v}_{\mathcal{D}})^{1/2}$. Let us now estimate $\mathfrak{I}_2$. Using the definition (35) of the (local) gradient reconstruction (with the fact that $\underline{q}_{\mathcal{D}} \in \underline{P}_{\mathcal{D},\Gamma}^k$), and integrating by parts the volumetric term therein, we infer

$$ \mathfrak{I}_2 = (\mathbf{a}, \mathbf{grad}_{\mathcal{T}} q_{\mathcal{T}})_{\Omega} + \sum_{T \in \mathcal{T}} \sum_{F \in \mathcal{F}_T} (\pi_{\mathcal{P},T}^k (\mathbf{a}|_T)|_F \cdot \mathbf{n}_{T,F}, q_F - q_{T|F})_F. $$

A cell-by-cell integration by parts of the first term in the right-hand side, and the combination of $\mathbf{a} \in \mathbf{H}(\mathrm{div}^0; \Omega)$ (according to (46b)) and $\underline{q}_{\mathcal{D}} \in \underline{P}_{\mathcal{D},\Gamma}^k$ (cf. (32) and (33)), yields

$$ \mathfrak{I}_2 = \sum_{T \in \mathcal{T}} \sum_{F \in \mathcal{F}_T} \varepsilon_{T,F} ((\pi_{\mathcal{P},T}^k (\mathbf{a}|_T)|_F - \mathbf{a}|_F) \cdot \mathbf{n}_F, q_F - q_{T|F})_F + \sum_{j=1}^{\beta_2} \sum_{F \in \mathcal{F}_j^\theta} \gamma_j (\mathbf{a}|_F \cdot \mathbf{n}_F, 1)_F, $$

where we have additionally utilized the fact that $\mathbf{a} \in \mathbf{H}^s(\mathcal{T})$ for some $s > 1$ ($s > \frac{1}{2}$ would be sufficient) to give a standard meaning to its normal traces. Now, since, for any $j \in \{1, \dots, \beta_2\}$, there holds $\sum_{F \in \mathcal{F}_j^\theta} \gamma_j (\mathbf{a}|_F \cdot \mathbf{n}_F, 1)_F = \gamma_j (\mathbf{a}|_{\Gamma_j} \cdot \mathbf{n}_F, 1)_{\Gamma_j}$, we infer by (46d) that

$$ \mathfrak{I}_2 = \sum_{T \in \mathcal{T}} \sum_{F \in \mathcal{F}_T} \varepsilon_{T,F} ((\pi_{\mathcal{P},T}^k (\mathbf{a}|_T)|_F - \mathbf{a}|_F) \cdot \mathbf{n}_F, q_F - q_{T|F})_F. $$

The approximation results of [14, Lem. 2.5], along with standard arguments, then show that

$$ \mathfrak{I}_2 \lesssim \left( \sum_{T \in \mathcal{T}} \mu_T^{-1} h_T^{2(s-1)} |\mathbf{a}|_{s,T}^2 \right)^{1/2} \left( \sum_{T \in \mathcal{T}} \mu_T S_{\mathbf{grad},T} (\underline{q}_T, \underline{q}_T) \right)^{1/2}. $$

Since $\sum_{T \in \mathcal{T}} \mu_T S_{\mathbf{grad},T}(\underline{q}_T, \underline{q}_T) \le S_{\mathcal{D}}(\underline{q}_{\mathcal{D}}, \underline{q}_{\mathcal{D}})$, gathering the estimates on $\mathfrak{I}_1$ and $\mathfrak{I}_2$, and using that $|(\underline{v}_{\mathcal{D}}, \underline{q}_{\mathcal{D}})|_{\mathcal{D}} = 1$, we conclude the proof. $\square$

From the error estimate (61), we infer that the HHO scheme converges with optimal rate $h_{\mathcal{D}}^k$ when $\mathbf{a} \in \mathbf{H}^{k+1}(\mathcal{T})$. A few remarks are necessary before concluding this section.

**Remark 13 (Low-regularity solutions).** In order to yield convergence, the energy-error estimate from Theorem 12 essentially requires that $\mathbf{a} \in \mathbf{H}^s(\mathcal{T})$ for some $s > 1$. By Remark 5, this regularity assumption on the solution is not met in general, even on a convex domain. However, it is always true that $\mathbf{a} \in \mathbf{H}^s(\mathcal{T})$ for some $s > \frac{1}{2}$. Recovering convergence for any $s > \frac{1}{2}$ would then necessitate to consider a weaker error measure for the scheme, formally on $h_T^{1/2} \mathbf{C}_{\mathcal{T}}^{k-1}$. This is precisely the approach pursued in [26, 27] in the DG context. Note that the correct treatment of interface terms also requires that $\mu^{-1} \mathbf{curl} \mathbf{a} \in \mathbf{H}^r(\mathcal{T})$ for some $r > \frac{1}{2}$. By Remark 5, the latter regularity assumption is met for a globally constant coefficient $\mu$, but is not satisfied in general. Adapting the arguments from [25, Sec. 5], it should however be possible to relax this assumption to $r > 0$, which would cover the general case. We do not pursue herein further in this direction.
**Remark 14** (Variant of the scheme). It is possible, instead of (locally) reconstructing the rotational in $\mathcal{P}^{k-1}(T)$ as in (24), to reconstruct it in the smaller space $\mathcal{R}^{k-1}(T) = \text{curl}(\mathcal{P}^k(T))$. This was the strategy advocated in [11, Sec. 3.2]. In this case, the proxy for **curl** **a** is given by $\pi_{\mathcal{R},\mathcal{T}}^{k-1}(\text{curl } \mathbf{a})$ (with obvious notation for the projector) in place of $\pi_{\mathcal{P},\mathcal{T}}^{k-1}(\text{curl } \mathbf{a})$. This proxy has approximability properties no better than $\text{curl}_{\mathcal{T}}(\pi_{\mathcal{P},\mathcal{T}}^k(\mathbf{a}))$ (see Remark 9). Therefore, since this alternative construction does not allow to relax (21), and given that the complexity gain in using $\mathcal{R}^{k-1}(T)$ instead of $\mathcal{P}^{k-1}(T)$ to reconstruct the rotational is marginal (local systems are small, and local computations are embarrassingly parallel), it is definitely a better option to stick to the choice $\mathcal{P}^{k-1}(T)$ as in (24).

**Remark 15** (Robustness to heterogeneity). Let us introduce $\kappa_{\mu} := \mu_{\#}/\mu_{\flat} \ge 1$ the (global) heterogeneity ratio of the magnetic permeability $\mu$. From (61), one can infer that $|\underline{a}_{\mathcal{D}} - \underline{I}_{\mathcal{D}}^{k}(\mathbf{a})|_{\text{curl},\mathcal{D}}$ is bounded by $\kappa_{\mu}^{1/2}$. Since $\mathbf{b} = \text{curl } \mathbf{a}$ (recall that $\mathbf{a}$ is the magnetic vector potential), notice that $|\underline{a}_{\mathcal{D}} - \underline{I}_{\mathcal{D}}^{k}(\mathbf{a})|_{\text{curl},\mathcal{D}}$ is actually an $L^2$-like measure of the scheme error on $\mathbf{b}$. The computed magnetic induction is thus only mildly sensitive to heterogeneity.

# 5 Numerical experiments

We perform in this section a comprehensive numerical assessment of the HHO schemes we devised and analyzed in Section 4.2. One important goal is to validate numerically the theoretical rates of convergence of the methods. Our focus is on relatively simple geometries, and on test-cases for which an analytical solution is available. The computational domains under study include topologically non-trivial geometries, as well as a domain with reentrant edge, for which the corresponding solution is singular.

We consider two families of (regular) mesh sequences, respectively composed of *tetrahedral* and *polyhedral* meshes. The polyhedral meshes are constructed from the tetrahedral ones by agglomerating tetrahedral elements into polyhedra in a random fashion. The agglomeration procedure concerns all tetrahedra surrounding randomly selected mesh vertices, in such a way that one can construct a mesh sequence whose regularity parameter does not decrease upon refinement. We here solve the discrete Problems (52) and (58), respectively referred to as first- and second-order Problems, modulo the following adaptations: (i) building on Remark 9, for the first-order Problem, when the exact solution is regular, we replace the discrete rotational operator $C_{\mathcal{T}}^{k-1}$ by its broken version **curl**$_{\mathcal{T}}$, and the magnetic face space is then set to $\mathcal{Q}^k(F) := \mathcal{R}^k(F)$ for all $F \in \mathcal{F}$; building on Remark 10, (ii) when $\mathcal{T}$ is a tetrahedral mesh, we remove from both the first- and second-order Problems the contribution $S_{\mathcal{D}}(\underline{p}_{\mathcal{D}}, \underline{q}_{\mathcal{D}})$, whereas (iii) when $\mathcal{T}$ is a polyhedral mesh, we remove from both the first- and second-order Problems, in the definition of $S_{\mathcal{D}}(\underline{p}_{\mathcal{D}}, \underline{q}_{\mathcal{D}})$, the contribution stemming from the local bilinear form $h_{\mathcal{T}}^2(\text{grad } p_{\mathcal{T}}, \text{grad } q_{\mathcal{T}})_\mathcal{T}$. Doing so, we achieve, without compromising stability, a somewhat leaner implementation, and we bear consistency with the numerical results from [11], where these variants were actually implemented (in the case of topologically trivial domains only). Importantly, and as already noticed in [11], the numerical tests that we performed (not reported here) show that these variants achieve virtually the same accuracy as the original schemes (and, interestingly, for examples with non-trivial topology, even slightly improve it). As standard with skeletal methods, for each discrete Problem, both the magnetic and pressure-like cell unknowns are locally eliminated (without additional fill-in) from the global linear system using a Schur complement technique. The resulting (condensed) global linear system is solved using the SparseLU direct solver of the Eigen library, on a laptop Dell Precision 5570 equipped with an Intel Core i7-12800H processor clocked at 2.80GHz, and with 64Gb of RAM. The implementation uses the open-source C++ library
ParaSkel++<fnref n="1" /> [6]. For each computed discrete solution $\underline{u}_{\mathcal{D}} \in \{\underline{h}_{\mathcal{D}}, \underline{a}_{\mathcal{D}}\}$, corresponding to the exact solution $\underline{u} \in \{\underline{h}, \underline{a}\}$, we evaluate the relative energy-error and $L^2$-error by

$$ \text{En\_err} := \frac{\left( \sum_{T \in \mathcal{T}} \eta_T^{-1} |\underline{u}_T - \underline{I}_T^k (\underline{u}_T)|_{\text{curl},T}^2 \right)^{1/2}}{\left( \sum_{T \in \mathcal{T}} \eta_T^{-1} |\underline{I}_T^k (\underline{u}_T)|_{\text{curl},T}^2 \right)^{1/2}} \quad \text{and} \quad \text{L2\_err} := \frac{\|\underline{u}_{\mathcal{T}} - \underline{\pi}_{\mathcal{P},\mathcal{T}}^k(\underline{u})\|_{0,\Omega}}{\|\underline{\pi}_{\mathcal{P},\mathcal{T}}^k(\underline{u})\|_{0,\Omega}}, $$

with $\eta := 1$ in the first-order case, and $\eta := \mu$ in the second-order one. For all the numerical examples below, except in Section 5.4, we set $\mu = 1$. Whenever the computational domain is non-simply-connected, and some additional flux constraints need to be enforced through cutting surfaces so as to guarantee uniqueness for the solution, one needs to compute, from the partition at hand, a discrete representative of each cutting surface. To this aim, we implemented a simplified version of the (co)homology computation algorithm proposed in [1]. We refer the reader to Figure 1 for an example, in the case of a torus discretized by a tetrahedral mesh, of a discrete cutting surface computed by the latter algorithm.

<figure><img src="image_2.png" /><figcaption>Figure 1: **Test-case 5.1 (torus).** Example of a discrete cutting surface computed by the algorithm from [1]. The surface is piecewise planar, and the depicted triangles are faces of the tetrahedral mesh discretizing $\Omega$.</figcaption></figure>

## 5.1 Torus

Our first example consists in a toroidal domain $\Omega$, with midline radius $R = 2$ and internal radius $r = 1$. For this geometry, $\beta_1 = 1$, and $\beta_2 = 0$. We consider the first-order Problem (42), with following smooth exact solution, expressed in cylindrical coordinates:

$$ \mathbf{h}(\rho, \varphi, z) := \left( \cos\left(\frac{\pi}{r}\sqrt{(\rho - R)^2 + z^2}\right) - c \right) \hat{\varphi}, \quad (62) $$

where the constant $c \in \mathbb{R}$ is set so that, for some cutting surface $\Sigma$, the condition $\langle \mathbf{h}_{|\Sigma} \cdot \mathbf{n}_{\Sigma}, 1 \rangle_{\Sigma} = 0$ is satisfied. In our numerical experiments, the electric current density $j$ is set according to (42a), whereas the zero normal boundary condition (42c) and zero flux condition (42d) (valid, respectively, on the curved boundary $\Gamma$, and on the $\Gamma$-delineated cutting surface $\Sigma$) are replaced by their non-homogeneous versions computed from (62) on the given piecewise polygonal approximations of $\Gamma$ and $\Sigma$. For each $k \in \{1, 2, 3\}$, we solve Problem (52) (with $\text{curl}_{\mathcal{T}}$), and depict on Figure 2, for both the tetrahedral and (agglomerated) polyhedral mesh families, the relative energy-error and $L^2$-error as functions of the mesh size ($h_{\mathcal{D}}$), and of the number of degrees of freedom (#DoF) after static

<footnote n="1">¹cf. https://gitlab.inria.fr/simlemai/paraskel (under GNU LGPL v3.0).</footnote>
condensation. For the two mesh families, as predicted by Theorem 7, we obtain a convergence rate for the energy-error of order $k$. We also observe a convergence rate of order $k + 1$ for the $L^2$-error on the magnetic field.

<figure><img src="image_3.png" /><figcaption>Figure 2: Test-case 5.1 (torus). Relative energy-error (top row) and $L^2$-error (bottom row) vs. $h_D$ (left column) and \#DoF (right column) on tetrahedral (solid) and polyhedral (dashed) meshes.</figcaption></figure>

## 5.2 Hollow ball

Our second example consists in a ball $\Omega$ of radius 2, encapsulating a concentric hollow spheric cavity of radius 1, so that its boundary $\Gamma$ consists of two connected components $\Gamma_0$ and $\Gamma_1$. For such a geometry, $\beta_1 = 0$, and $\beta_2 = 1$. We consider the second-order Problem (46), with following smooth exact solution, expressed in spherical coordinates:

$$a(r, \theta, \varphi) := \frac{1}{r^2} \cos(\theta) \hat{r}. \qquad (63)$$

Remark that $a$ satisfies $\langle a|_{\Gamma_1} \cdot n, 1 \rangle_{\Gamma_1} = 0$ (as well as $\langle a|_{\Gamma_0} \cdot n, 1 \rangle_{\Gamma_0} = 0$, since $\text{div } a = 0$ in $\Omega$). In our numerical experiments, the electric current density $j$ is set according to (46a), whereas the zero tangential boundary condition (46c) and zero flux condition (46d) (valid on the curved boundaries $\Gamma$ and $\Gamma_1$) are replaced by their non-homogeneous versions computed from (63) on the given piecewise polygonal approximations of $\Gamma$ and $\Gamma_1$. For each $k \in \{1, 2, 3\}$, we solve Problem (58), and plot on Figure 3 the relative energy-error and $L^2$-error as functions of the mesh size ($h_D$) and of the number of degrees of freedom (\#DoF) after static condensation, for both the tetrahedral and polyhedral mesh families. For both mesh families, we obtain, as predicted by Theorem 12, a convergence rate for the energy-error of order $k$. We also observe a convergence rate of order $k + 1$ for the $L^2$-error on the magnetic vector potential.
<figure><img src="image_4.png" /><figcaption>Figure 3: **Test-case 5.2 (hollow ball).** Relative energy-error (top row) and $L^2$-error (bottom row) vs. $h_D$ (left column) and #DoF (right column) on tetrahedral (solid) and polyhedral (dashed) meshes.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: **Test-case 5.3 (toroidal vacuum chamber).** The toroidal domain $\Omega$, with coaxial hollow toric cavity. Depicted are two possible cutting surfaces $\Sigma_1$ and $\Sigma_2$, as well as the two connected components $\Gamma_0$ and $\Gamma_1$ of $\Gamma$.</figcaption></figure>

## 5.3 Toroidal vacuum chamber

This test-case is defined on the same toroidal domain as in Section 5.1, except that $\Omega$ here additionally encapsulates a coaxial (i.e. with midline radius $R$) hollow toric cavity of internal radius $r/2$ (that is, the cavity has half the internal radius of the torus). For this geometry, $\beta_1 = 2$, and $\beta_2 = 1$. We solve here the field formulation, with the same exact solution (62) (also same value for the constant $c$) as in the fully solid torus test-case of Section 5.1. We depict on Figure 4 two possible cutting surfaces $\Sigma_1$ and $\Sigma_2$, as well as the two connected components $\Gamma_0$ and $\Gamma_1$ of the boundary $\Gamma$. The vector field $h$ given by (62) is actually solution to a non-homogeneous version of the first-order Problem (42): it satisfies (42b)–(42c) and $\langle h|_{\Sigma_1} \cdot n|_{\Sigma_1}, 1 \rangle_{\Sigma_1} = 0$, but $\langle h|_{\Sigma_2} \cdot n|_{\Sigma_2}, 1 \rangle_{\Sigma_2} \neq 0$, so that $h$ “embeds” an harmonic field in $H_n(\Omega)$ (cf. (12)). Here again, for $k \in \{1, 2, 3\}$, we solve Problem (52) (with $\mathrm{curl}_T$), and plot on Figure 5 the relative energy-error and $L^2$-error as functions of $h_D$ and #DoF, for both the tetrahedral and polyhedral mesh families. The obtained energy-error convergence rates are again
consistent with our theory, as established in Theorem 7. The observed L²-error convergence rates on the magnetic field are again of order k + 1.

<figure><img src="image_6.png" /><figcaption>Figure 5: Test-case 5.3 (toroidal vacuum chamber). Relative energy-error (top row) and L²-error (bottom row) vs. h_D (left column) and #DoF (right column) on tetrahedral (solid) and polyhedral (dashed) meshes.</figcaption></figure>

## 5.4 Variable permeability

In this example, we let $\Omega := (0, 1)^3$, and we solve the vector potential formulation, with smooth exact solution

$$ a(x, y, z) := \begin{pmatrix} 2\pi r(x) \sin(2\pi y) \cos(2\pi z) \\ -r'(x) \cos(2\pi y) \cos(2\pi z) \\ -2r'(x) \sin(2\pi y) \sin(2\pi z) \end{pmatrix}, \quad r(x) := x^3, \qquad (64) $$

together with smooth, locally variable (inverse) magnetic permeability $\mu^{-1}(x, y, z) := 1 + x^2 y^2$. The vector field $a$ given by (64) is solution to a non-homogeneous version of the second-order Problem (46): it satisfies (46b) (note that (46d) is here trivial), but fulfills in place of (46c) a non-homogeneous tangential boundary condition. In our experiments, the electric current density $j$ is set according to (46a), and the boundary datum in accordance with (64). Remark that $\mu^{-1}$ is not piecewise constant over the partition of the domain. So as to retain optimal convergence rates, we thus had in our implementation to use weighted cell/face inner products (remark that $\mu^{-1}$ is continuous across interfaces). The relative energy-error is computed accordingly. Their accurate computation necessitates to increase the degree of exactness of the quadrature rules. We here solve Problem (58). The convergence results, collected in Figure 6, show essentially the same behavior as those reported in [11, Sec. 3.2.5] for a unit permeability test-case on the same domain.
<figure><img src="image_7.png" /><figcaption>Figure 6: Test-case 5.4 (variable permeability). Relative energy-error (top row) and L²-error (bottom row) versus h_D (left column) and #DoF (right column) on tetrahedral (solid) and polyhedral (dashed) meshes.</figcaption></figure>

## 5.5 Singular solution

In this last test-case, the domain $\Omega$ is a cylinder of unit radius and height, to which a quarter has been removed, so that it presents a reentrant edge (cf. the left panel of Figure 7). We here solve the field formulation, for the following singular exact solution, expressed in cylindrical coordinates:

$$h(\rho, \varphi, z) := \operatorname{grad} \left( \rho^{\frac{2}{3}} \cos \left( \frac{2}{3} \varphi \right) \right). \qquad (65)$$

The vector field $h$ given by (65) is solution to a non-homogeneous version of the first-order Problem (42): it does satisfy (42b) (remark that (42d) is here trivial), but fulfills in place of (42c) a non-homogeneous normal boundary condition on the curved section of $\Gamma$. In our numerical experiments, the electric current density $j$ is set to zero (accordingly to (42a)), and the boundary datum is set in accordance with (65) on the given piecewise polygonal approximation of the curved section of $\Gamma$. It can be verified that $h \in H^{\frac{2}{3}-\varepsilon}(\Omega)$ for all $\varepsilon > 0$, with $h$ being singular along the reentrant edge. We solve Problem (52) (with $C_{\tau}^{k-1}$) for $k = 1$ on the tetrahedral mesh sequence. Note that the regularity assumptions of Theorem 7 are violated, hence our theory does not apply here. We collect on the right panel of Figure 7 the results for the relative L²-error on the magnetic field. We observe that the L²-error attains its maximum possible rate of convergence of 2/3; no convergence is observed in energy-norm (not reported). We verified, for each $k \in \{2, 3\}$, that the L²-error convergence rate does not depend on the polynomial degree. Finally, we performed the exact same test, but with discrete rotational operator given by $\mathrm{curl}_{\tau}$ instead of $C_{\tau}^{k-1}$. We observed no convergence whatsoever, even in L²-norm. This last observation corroborates the discussion in Remark 9.
<figure><img src="image_8.png" /><figcaption>Figure 7: **Test-case 5.5 (singular solution).** The computational domain $\Omega$ (left panel). Relative $L^2$-error and convergence rate on tetrahedral meshes for $k = 1$ (right panel).</figcaption></figure>

<table><thead><tr><td>$h_D$</td><td>$L^2$-error</td><td>Rate</td></tr></thead><tbody><tr><td>4.20E-01</td><td>3.62E-02</td><td>-</td></tr><tr><td>2.37E-01</td><td>2.72E-02</td><td>0.50</td></tr><tr><td>1.44E-01</td><td>1.99E-02</td><td>0.63</td></tr><tr><td>7.00E-02</td><td>1.24E-02</td><td>0.66</td></tr></tbody></table>

## Acknowledgments

The work of the authors is supported by the action "Préservation de l'emploi de R&D" from the "Plan de Relance" (recovery plan) of the French State, as well as by the Agence Nationale de la Recherche (ANR) under the PRCE grant HIPOTHEC (ANR-23-CE46-0013). The last two authors also acknowledge support from the LabEx CEMPI (ANR-11-LABX-0007).

## References

[1] A. Alonso Rodríguez, E. Bertolazzi, R. Ghiloni, and R. Specogna. Efficient construction of 2-chains representing a basis of $H_2(\bar{\Omega}, \partial\Omega; \mathbb{Z})$. *Adv. Comput. Math.*, 44(5):1411–1440, 2018.

[2] C. Amrouche, C. Bernardi, M. Dauge, and V. Girault. Vector potentials in three-dimensional non-smooth domains. *Math. Methods Appl. Sci.*, 21(9):823–864, 1998.

[3] D. N. Arnold. *Finite element exterior calculus*, volume 93 of CBMS-NSF Regional Conference Series in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2018.

[4] F. Assous, P. Ciarlet Jr., and S. Labrunie. *Mathematical foundations of computational electromagnetism*, volume 198 of Applied Mathematical Sciences. Springer, Cham, 2018.

[5] B. Ayuso de Dios, K. Lipnikov, and G. Manzini. The nonconforming virtual element method. *ESAIM Math. Model. Numer. Anal.*, 50(3):879–904, 2016.

[6] L. Beaude and S. Lemaire. ParaSkel++: a C++ platform for the high-performance, arbitrary-order, 2/3D numerical approximation of PDEs on general polytopal meshes using skeletal Galerkin methods, 2021.

[7] R. Benedetti, R. Frigerio, and R. Ghiloni. The topology of Helmholtz domains. *Expo. Math.*, 30(4):319–375, 2012.

[8] A. Bonito, J.-L. Guermond, and F. Luddens. Regularity of the Maxwell equations in heterogeneous media and Lipschitz domains. *J. Math. Anal. Appl.*, 408(2):498–512, 2013.

[9] S. Cao, C. Wang, and J. Wang. A new numerical method for div-curl systems with low regularity assumptions. *Comput. Math. Appl.*, 114:47–59, 2022.

[10] W. Cao, C. Wang, and J. Wang. An $L^p$-primal-dual weak Galerkin method for div-curl systems. *J. Comput. Appl. Math.*, 422:Paper No. 114881, 17, 2023.

[11] F. Chave, D. A. Di Pietro, and S. Lemaire. A discrete Weber inequality on three-dimensional hybrid spaces with application to the HHO approximation of magnetostatics. *Math. Models Methods Appl. Sci.*, 32(1):175–207, 2022.

[12] G. Chen, J. Cui, and L. Xu. Analysis of a Hybridizable Discontinuous Galerkin method for the Maxwell operator. *ESAIM Math. Model. Numer. Anal.*, 53(1):301–324, 2019.
[13] H. Chen, W. Qiu, K. Shi, and M. Solano. A superconvergent HDG method for the Maxwell equations. *J. Sci. Comput.*, 70(3):1010–1029, 2017.

[14] M. Cicuttin, A. Ern, and N. Pignet. *Hybrid high-order methods – A primer with applications to solid mechanics*. SpringerBriefs in Mathematics. Springer, Cham, 2021.

[15] B. Cockburn. Static condensation, hybridization, and the devising of the HDG methods. In *Building bridges: connections and challenges in modern approaches to numerical partial differential equations*, volume 114 of *Lect. Notes Comput. Sci. Eng.*, pages 129–177. Springer, [Cham], 2016.

[16] B. Cockburn, D. A. Di Pietro, and A. Ern. Bridging the Hybrid High-Order and Hybridizable Discontinuous Galerkin methods. *ESAIM Math. Model. Numer. Anal.*, 50(3):635–650, 2016.

[17] B. Cockburn, J. Gopalakrishnan, and R. Lazarov. Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems. *SIAM J. Numer. Anal.*, 47(2):1319–1365, 2009.

[18] M. Costabel. A remark on the regularity of solutions of Maxwell’s equations on Lipschitz domains. *Math. Methods Appl. Sci.*, 12(4):365–368, 1990.

[19] M. Costabel, M. Dauge, and S. Nicaise. Singularities of Maxwell interface problems. *M2AN Math. Model. Numer. Anal.*, 33(3):627–649, 1999.

[20] R. Dautray and J.-L. Lions. *Mathematical analysis and numerical methods for science and technology – Vol. 3: Spectral theory and applications*. Springer-Verlag, Berlin, 1990. With the collaboration of Michel Artola and Michel Cessenat. Translated from the French by John C. Amson.

[21] D. A. Di Pietro and J. Droniou. *The Hybrid High-Order method for polytopal meshes – Design, analysis, and applications*, volume 19 of *Modeling, Simulation and Applications*. Springer, Cham, 2020.

[22] D. A. Di Pietro and A. Ern. A hybrid high-order locking-free method for linear elasticity on general meshes. *Comput. Methods Appl. Mech. Engrg.*, 283:1–21, 2015.

[23] D. A. Di Pietro, A. Ern, and S. Lemaire. An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators. *Comput. Methods Appl. Math.*, 14(4):461–472, 2014.

[24] S. Du and F.-J. Sayas. A unified error analysis of Hybridizable Discontinuous Galerkin methods for the static Maxwell equations. *SIAM J. Numer. Anal.*, 58(2):1367–1391, 2020.

[25] A. Ern and J.-L. Guermond. Quasi-optimal nonconforming approximation of elliptic PDEs with contrasted coefficients and $H^{1+r}$, $r > 0$, regularity. *Found. Comput. Math.*, 22(5):1273–1308, 2022.

[26] A. Ern and J.-L. Guermond. The discontinuous Galerkin approximation of the grad-div and curl-curl operators in first-order form is involution-preserving and spectrally correct. *SIAM J. Numer. Anal.*, 61(6):2940–2966, 2023.

[27] A. Ern and J.-L. Guermond. Spectral correctness of the simplicial discontinuous Galerkin approximation of the first-order form of Maxwell’s equations with discontinuous coefficients. *SIAM J. Numer. Anal.*, 63(2):661–684, 2025.

[28] V. Girault and P.-A. Raviart. *Finite element methods for Navier–Stokes equations – Theory and algorithms*, volume 5 of *Springer Series in Computational Mathematics*. Springer-Verlag, Berlin, 1986.

[29] P. W. Gross and P. R. Kotiuga. *Electromagnetic theory and computation: a topological approach*, volume 48 of *Mathematical Sciences Research Institute Publications*. Cambridge University Press, Cambridge, 2004.

[30] F. Jochmann. Regularity of weak solutions of Maxwell’s equations with mixed boundary conditions. *Math. Methods Appl. Sci.*, 22(14):1255–1274, 1999.

[31] F. Kikuchi. Mixed formulations for finite element analysis of magnetostatic and electrostatic problems. *Japan J. Appl. Math.*, 6(2):209–221, 1989.
[32] C. Lehrenfeld and J. Schöberl. High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows. *Comput. Methods Appl. Mech. Engrg.*, 307:339–361, 2016.

[33] S. Lemaire and S. Pitassi. Discrete Weber inequalities and related Maxwell compactness for hybrid spaces over polyhedral partitions of domains with general topology. *Found. Comput. Math.*, 25(3):725–763, 2025.

[34] J. Li, X. Ye, and S. Zhang. A weak Galerkin least-squares finite element method for div-curl systems. *J. Comput. Phys.*, 363:79–86, 2018.

[35] K. Lipnikov and G. Manzini. A high-order mimetic method on unstructured polyhedral meshes for the diffusion equation. *J. Comput. Phys.*, 272:360–385, 2014.

[36] L. Mu, J. Wang, X. Ye, and S. Zhang. A weak Galerkin finite element method for the Maxwell equations. *J. Sci. Comput.*, 65(1):363–386, 2015.

[37] N. C. Nguyen, J. Peraire, and B. Cockburn. Hybridizable Discontinuous Galerkin methods for the time-harmonic Maxwell's equations. *J. Comput. Phys.*, 230(19):7151–7175, 2011.

[38] C. Wang and J. Wang. Discretization of div-curl systems by weak Galerkin finite element methods on polyhedral partitions. *J. Sci. Comput.*, 68(3):1144–1171, 2016.

[39] J. Wang and X. Ye. A weak Galerkin finite element method for second-order elliptic problems. *J. Comput. Appl. Math.*, 241:103–115, 2013.