# p-robust equilibrated flux reconstruction in H(curl) based on local minimizations. Application to a posteriori analysis of the curl-curl problem

Théophile Chaumont-Frelet, Martin Vohralík

► To cite this version:

Théophile Chaumont-Frelet, Martin Vohralík. p-robust equilibrated flux reconstruction in H(curl) based on local minimizations. Application to a posteriori analysis of the curl-curl problem. SIAM Journal on Numerical Analysis, 2023, 61 (4), pp.1783–1818. 10.1137/21M141909X. hal-03227570v2

HAL Id: hal-03227570

https://inria.hal.science/hal-03227570v2

Submitted on 20 Aug 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
p-robust equilibrated flux reconstruction in $H(\text{curl})$ based on
local minimizations. Application to a posteriori analysis of the
curl-curl problem*

Théophile Chaumont-Frelet†‡ Martin Vohralík§¶

August 20, 2022

Abstract

We present a local construction of $H(\text{curl})$-conforming piecewise polynomials satisfying a prescribed curl constraint. We start from a piecewise polynomial not contained in the $H(\text{curl})$ space but satisfying a suitable orthogonality property. The procedure employs minimizations in vertex patches and the outcome is, up to a generic constant independent of the underlying polynomial degree, as accurate as the best-approximations over the entire local versions of $H(\text{curl})$. This allows to design guaranteed, fully computable, constant-free, and polynomial-degree-robust a posteriori error estimates of the Prager-Synge type for Nédélec's finite element approximations of the curl-curl problem. A divergence-free decomposition of a divergence-free $H(\text{div})$-conforming piecewise polynomial, relying on over-constrained minimizations in Raviart-Thomas' spaces, is the key ingredient. Numerical results illustrate the theoretical developments.

**Key words:** Sobolev space $H(\text{curl})$, Sobolev space $H(\text{div})$, equilibrated flux reconstruction, p-robustness, a posteriori error estimate, divergence-free decomposition, broken polynomial extension

# 1 Introduction

A posteriori error estimation by equilibrated flux reconstruction has achieved a great attention for elliptic model problems like the Poisson problem. For an $H^1$-conforming discretization whose flux is not in $H(\text{div})$, one has to reconstruct a flux in $H(\text{div})$ satisfying a prescribed divergence constraint. To design high-performance algorithms, the procedure must furthermore be localized and can not involve a solution of any supplementary global problem. Then, a guaranteed, fully computable, and constant-free upper bound on the unknown discretization error follows from the equality of Prager and Synge [35]. There are several techniques of such an equilibrated flux reconstruction. Following Ladevèze and Leguillon [29] and Ainsworth and Oden [2], normal fluxes on mesh faces can first be constructed and then lifted elementwise as in Nicaise et al. [34], dual Voronoï-type grids can be employed for local non-overlapping minimizations in $H(\text{div})$ as in Luce and Wohlmuth [31] or Hannukainen et al. [27], or a localization by the partition of unity via the finite element hat basis functions can be used for an overlapping combination of best-possible vertex-patch fluxes as in Destuynder and Métivet [15] or Braess and Schöberl [8]. This last approach is conceptual and, as established in Braess et al. [7] and Ern and Vohralík [21], it gives estimates robust with respect to the polynomial degree $p$ (henceforth termed p-robust).

In contrast, there is only a handful of results available for the curl-curl problem, where, for an $H(\text{curl})$-conforming discretization whose curl is not in $H(\text{curl})$, one has to locally reconstruct a flux in $H(\text{curl})$ satisfying a prescribed curl constraint. An approach based on patchwise minimizations for the lowest-order case $p = 0$ has been designed in [8]. Its generalization for arbitrary $p \ge 1$, however, turns surprisingly difficult and, to the best of our knowledge, has not been presented yet. Several workarounds appeared in the literature recently, though. A conceptual discussion appears in Licht [30], whereas a

<footnote>*This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 647134 GATIPOR).</footnote>

<footnote>†Inria, 2004 Route des Lucioles, 06902 Valbonne, France (theophile.chaumont@inria.fr).</footnote>

<footnote>‡Laboratoire J.A. Dieudonné, Parc Valrose, 28 Avenue Valrose, 06108 Nice, France</footnote>

<footnote>§Inria, 2 rue Simone Iff, 75589 Paris, France (martin.vohralik@inria.fr).</footnote>

<footnote>¶CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée, France</footnote>
construction following in spirit [29, 2] has been proposed and analyzed in Gedicke et al. [23]. This last approach has been recently modified in Gedicke et al. [24] in order to achieve p-robustness. A broken patchwise equilibration procedure that bypasses the Prager-Synge theorem is proposed and proved p-robust in Chaumont-Frelet et al. [10]; it relies on smaller edge patches, but the arising estimates are not constant-free.

The purpose of this contribution is to design an equilibrated flux reconstruction in $H(\text{curl})$ employing best-possible local fluxes. In doing so, we rely on localization by the partition of unity via the hat functions and overlapping flux combinations, in generalization of the concept of [8] to arbitrary $p \ge 0$. Consequently, we identify the equivalent in $H(\text{curl})$ of the concept of equilibrated flux reconstruction in $H(\text{div})$ from [15, 8, 7, 20, 21]. This is then used for a posteriori error estimation when the Nédélec (edge) finite elements of arbitrary degree $p \ge 0$ are used for approximation of the curl-curl problem. It leads to guaranteed, fully computable, and constant-free a posteriori error estimates that are locally efficient and robust with respect to the polynomial degree $p$; (higher-order) data oscillation terms are rigorously included in our analysis. Our p-robust efficiency proofs are based on the seminal volume and tangential trace p-robust extensions on a single tetrahedron of Costabel and McIntosh [13, Proposition 4.2] and Demkowicz et al. [14, Theorem 7.2]. These results were recently extended into a stable broken polynomial extension for a single tetrahedron in Chaumont-Frelet et al. [9, Theorem 2], an edge patch of tetrahedra in Chaumont-Frelet et al. [10, Theorem 3.1], and for a vertex patch of tetrahedra in Chaumont-Frelet and Vohralík [11, Theorem 3.3, see also Corollary 4.3].

An important step in the construction of our estimators is to decompose the given divergence-free right-hand side into locally supported divergence-free contributions. Starting from the available (lowest-order Galerkin) orthogonality property, we propose a multi-stage procedure relying on two central technical results of independent interest: over-constrained minimization in Raviart-Thomas' spaces leading to suitable elementwise orthogonality properties, and a decomposition of a divergence-free piecewise polynomial with the above elementwise orthogonality properties into local divergence-free contributions. These issues are related to the developments on divergence-free decompositions in Scheichl [37], Alonso Rodríguez et al. [3, 4], and the references therein.

This contribution is organized as follows. Section 2 fixes the notation. Section 3 introduces the curl-curl problem, its Nédélec finite element discretization, and identifies therefrom two abstract assumptions under which our analysis is performed. In Section 4, we motivate our approach at the continuous level. Section 5 then presents our main results: in Section 5.1, we develop a divergence-free decomposition of the given target curl; in Section 5.2, we present the equilibrated flux reconstruction based on local minimization in $H(\text{curl})$, as well as its p-robust stability; finally, these abstract results are applied in Section 5.3 to the Nédélec finite element discretization of the curl-curl problem. Section 6 is dedicated to a numerical illustration, whereas Section 7 collects some technical details and proofs. Finally, in Appendices A and B, we present the two central technical results on over-constrained minimization and divergence-free decomposition.

# 2 Notation

The purpose of this section is to set the necessary notation. Let $\omega, \Omega \subset \mathbb{R}^3$ be open, Lipschitz polyhedra; $\Omega$ will be used to denote the computational domain, while we reserve the notation $\omega \subseteq \Omega$ for its simply connected subsets. Notice that we do not require $\Omega$ to be simply connected.

## 2.1 Sobolev spaces $H^1$, $H(\text{curl})$, and $H(\text{div})$

We let $L^2(\omega)$ be the space of scalar-valued square-integrable functions defined on $\omega$; we use the notation $L^2(\omega) := [L^2(\omega)]^3$ for vector-valued functions with each component in $L^2(\omega)$. We denote by $\|\cdot\|_{\omega}$ the $L^2(\omega)$ or $L^2(\omega)$ norm and by $(\cdot, \cdot)_{\omega}$ the corresponding scalar product; we drop the index when $\omega = \Omega$. We will extensively work with the following three Sobolev spaces: 1) $H^1(\omega)$, the space of scalar-valued $L^2(\omega)$ functions with weak gradients in $L^2(\omega)$, $H^1(\omega) := \{v \in L^2(\omega); \nabla v \in L^2(\omega)\}$; 2) $H(\text{curl}, \omega)$, the space of vector-valued $L^2(\omega)$ functions with weak curls in $L^2(\omega)$, $H(\text{curl}, \omega) := \{v \in L^2(\omega); \nabla \times v \in L^2(\omega)\}$; and 3) $H(\text{div}, \omega)$, the space of vector-valued $L^2(\omega)$ functions with weak divergences in $L^2(\omega)$, $H(\text{div}, \omega) := \{v \in L^2(\omega); \nabla \cdot v \in L^2(\omega)\}$. We refer the reader to Adams [1] and Girault and Raviart [25] for an in-depth description of these spaces. Moreover, component-wise $H^1(\omega)$ functions will be denoted by $H^1(\omega) := \{v \in L^2(\omega); v_i \in H^1(\omega), i = 1, \dots, 3\}$. We will employ the notation $\langle \cdot, \cdot \rangle_S$ for the integral product on boundary (sub)sets $S \subset \partial\omega$.
## 2.2 Sobolev spaces with partially vanishing traces on $\partial\Omega$

Let $\Gamma_D, \Gamma_N$ be two disjoint, relatively open, and possibly empty subsets of the computational domain boundary $\partial\Omega$ such that $\partial\Omega = \bar{\Gamma}_D \cup \bar{\Gamma}_N$. We assume in addition that each boundary face of the mesh $\mathcal{T}_h$ defined below lies entirely either in $\bar{\Gamma}_D$ or in $\bar{\Gamma}_N$. Then $H_{0,D}^1(\Omega)$ is the subspace of $H^1(\Omega)$ formed by functions vanishing on $\Gamma_D$ in the sense of traces, $H_{0,D}^1(\Omega) := \{v \in H^1(\Omega); v = 0 \text{ on } \Gamma_D\}$. Let $n_\Omega$ be the unit normal vector on $\partial\Omega$, outward to $\Omega$. Let $T = D$ or $N$; then $H_{0,T}(\text{curl}, \Omega)$ is the subspace of $H(\text{curl}, \Omega)$ formed by functions with vanishing tangential trace on $\Gamma_T$, $H_{0,T}(\text{curl}, \Omega) := \{v \in H(\text{curl}, \Omega); v \times n_\Omega = 0 \text{ on } \Gamma_T\}$, where $v \times n_\Omega = 0$ on $\Gamma_T$ means that $(\nabla \times v, \varphi) - (v, \nabla \times \varphi) = 0$ for all functions $\varphi \in H^1(\Omega)$ such that $\varphi \times n_\Omega = 0$ on $\partial\Omega \setminus \Gamma_T$. Finally, $H_{0,N}(\text{div}, \Omega)$ is the subspace of $H(\text{div}, \Omega)$ formed by functions with vanishing normal trace on $\Gamma_N$, $H_{0,N}(\text{div}, \Omega) := \{v \in H(\text{div}, \Omega); v \cdot n_\Omega = 0 \text{ on } \Gamma_N\}$, where $v \cdot n_\Omega = 0$ on $\Gamma_N$ means that $(v, \nabla \varphi) + (\nabla \cdot v, \varphi) = 0$ for all functions $\varphi \in H_{0,N}^1(\Omega)$. Fernandes and Gilardi [22] present a thorough characterization of tangential (resp. normal) traces of $H(\text{curl}, \Omega)$ (resp. $H(\text{div}, \Omega)$) on a part of the boundary $\partial\Omega$.

## 2.3 Cohomology space

The space $\mathcal{H}(\Omega, \Gamma_D)$ of functions $v \in H_{0,D}(\text{curl}, \Omega) \cap H_{0,N}(\text{div}, \Omega)$ such that $\nabla \times v = 0$ and $\nabla \cdot v = 0$ is the “cohomology” space associated with the domain $\Omega$ and the partition of its boundary $\partial\Omega = \bar{\Gamma}_D \cup \bar{\Gamma}_N$. When $\Omega$ is simply connected and $\Gamma_D$ is connected, this space is trivial; then the conditions associated with it below can be disregarded. In the general case, $\mathcal{H}(\Omega, \Gamma_D)$ is finite-dimensional, and its dimension depends on the topology of $\Omega$ and $\Gamma_D$, see [22, 26].

## 2.4 Tetrahedral mesh, patches of elements, and the hat functions

Let $\mathcal{T}_h$ be a simplicial mesh of the domain $\Omega$, i.e., $\cup_{K \in \mathcal{T}_h} K = \bar{\Omega}$, where any element $K \in \mathcal{T}_h$ is a closed tetrahedron with nonzero measure, and where the intersection of two different tetrahedra is either empty or their common vertex, edge, or face. The shape-regularity parameter of the mesh $\mathcal{T}_h$ is the positive real number $\kappa_{\mathcal{T}_h} := \max_{K \in \mathcal{T}_h} h_K / \rho_K$, where $h_K$ is the diameter of the tetrahedron $K$ and $\rho_K$ is the diameter of the largest ball contained in $K$. These assumptions are standard and allow for strongly graded meshes with local refinements. We will use the notation $a \lesssim b$ when there exists a positive constant $C$ only depending on $\kappa_{\mathcal{T}_h}$ such that $a \le Cb$.

We denote the set of vertices of the mesh $\mathcal{T}_h$ by $\mathcal{V}_h$; it is composed of interior vertices lying in $\Omega$ and of vertices lying on the boundary $\partial\Omega$. For an element $K \in \mathcal{T}_h$, $\mathcal{F}_K$ denotes the set of its faces and $\mathcal{V}_K$ the set of its vertices. Conversely, for a vertex $a \in \mathcal{V}_h$, $\mathcal{T}_a$ denotes the patch of the elements of $\mathcal{T}_h$ that share $a$, and $\omega_a$ is the corresponding open subdomain with diameter $h_{\omega_a}$. A particular role below will be played by the continuous, piecewise affine “hat” function $\psi^a$ which takes value 1 at the vertex $a$ and zero at the other vertices. We note that $\omega_a$ corresponds to the support of $\psi^a$ and that the functions $\psi^a$ form the partition of unity

$$ \sum_{a \in \mathcal{V}_h} \psi^a = 1. \qquad (2.1) $$

We will also need the patch $\tilde{\mathcal{T}}_a$ extended by one layer of neighbors and the associated subdomain $\tilde{\omega}_a$, corresponding to the supports of the hat functions $\psi^b$ for all vertices $b$ contained in the patch $\mathcal{T}_a$.

## 2.5 Piecewise polynomial spaces

Let $q \ge 0$ be an integer. For a single tetrahedron $K \in \mathcal{T}_h$, denote by $\mathcal{P}_q(K)$ the space of scalar-valued polynomials on $K$ of total degree at most $q$, and by $[\mathcal{P}_q(K)]^3$ the space of vector-valued polynomials on $K$ with each component in $\mathcal{P}_q(K)$. The Nédélec [6, 33] space of degree $q$ on $K$ is then given by

$$ \mathcal{N}_q(K) := [\mathcal{P}_q(K)]^3 + x \times [\mathcal{P}_q(K)]^3. \qquad (2.2) $$

Similarly, the Raviart-Thomas [6, 36] space of degree $q$ on $K$ is given by

$$ \mathcal{RT}_q(K) := [\mathcal{P}_q(K)]^3 + \mathcal{P}_q(K)x. \qquad (2.3) $$

We note that (2.2) and (2.3) are equivalent to the writing with a direct sum and only homogeneous polynomials in the second terms. The second term in (2.2) can also equivalently be given by homogeneous $(q+1)$-degree polynomials $v_h$ such that $x \cdot v_h(x) = 0$ for all $x \in K$.
We will below extensively use the broken, piecewise polynomial spaces formed from the above element spaces

$$
\begin{align*}
\mathcal{P}_q(\mathcal{T}_h) &:= \{v_h \in L^2(\Omega); v_h|_{K} \in \mathcal{P}_q(K) \quad \forall K \in \mathcal{T}_h\}, \\
\mathcal{N}_q(\mathcal{T}_h) &:= \{v_h \in L^2(\Omega); v_h|_{K} \in \mathcal{N}_q(K) \quad \forall K \in \mathcal{T}_h\}, \\
\mathcal{RT}_q(\mathcal{T}_h) &:= \{v_h \in L^2(\Omega); v_h|_{K} \in \mathcal{RT}_q(K) \quad \forall K \in \mathcal{T}_h\}.
\end{align*}
$$

To form the usual finite-dimensional Sobolev subspaces, we will write $\mathcal{P}_q(\mathcal{T}_h) \cap H^1(\Omega)$ (for $q \ge 1$), $\mathcal{N}_q(\mathcal{T}_h) \cap \mathbf{H}(\text{curl}, \Omega)$, $\mathcal{RT}_q(\mathcal{T}_h) \cap \mathbf{H}(\text{div}, \Omega)$ (both for $q \ge 0$), and similarly for the subspaces reflecting the different boundary conditions. The same notation will also be used on the patches of mesh elements $\mathcal{T}_a$.

## 2.6 $L^2$-orthogonal projectors and the Raviart–Thomas interpolator

For $q \ge 0$, let $\Pi_q$ denote the $L^2(K)$-orthogonal projector onto $\mathcal{P}_q(K)$. Since this is an elementwise procedure, we keep the same notation for the $L^2(\Omega)$-orthogonal projector onto $\mathcal{P}_q(\mathcal{T}_h)$ given for $v \in L^2(\Omega)$ as $\Pi_q(v) \in \mathcal{P}_q(\mathcal{T}_h)$ such that $(\Pi_q(v), w_h) = (v, w_h)$ for all $w_h \in \mathcal{P}_q(\mathcal{T}_h)$. Then, $\Pi_q$ is given componentwise by $\Pi_q$.

Let $K \in \mathcal{T}_h$ be a mesh tetrahedron and let $v \in [C^1(K)]^3$ be given. Following [6, 36], the canonical $q$-degree Raviart–Thomas interpolant $I_{K,q}^{\mathcal{RT}}(v) \in \mathcal{RT}_q(K)$, $q \ge 0$, is given by

$$ \langle I_{K,q}^{\mathcal{RT}}(v) \cdot n_K, r_h \rangle_F = \langle v \cdot n_K, r_h \rangle_F \quad \forall r_h \in \mathcal{P}_q(F), \quad \forall F \in \mathcal{F}_K, \quad (2.4a) $$

$$ (I_{K,q}^{\mathcal{RT}}(v), r_h)_K = (v, r_h)_K \quad \forall r_h \in [\mathcal{P}_{q-1}(K)]^3. \quad (2.4b) $$

Less regular functions can be used in (2.4), but $v \in [C^1(K)]^3$ will be sufficient for our purposes; we will actually only employ polynomial $v$ as arguments of $I_{K,q}^{\mathcal{RT}}$. This interpolator crucially satisfies, on the tetrahedron $K$, the commuting property

$$ \nabla \cdot I_{K,q}^{\mathcal{RT}}(v) = \Pi_q(\nabla \cdot v) \quad \forall v \in [C^1(K)]^3. \quad (2.5) $$

## 2.7 Sobolev spaces on the patch subdomains $\omega_a$

Let $a \in \mathcal{V}_h$ be an interior vertex. Then we set 1) $H_0^1(\omega_a) := \{v \in H^1(\omega_a); (v, 1)_{\omega_a} = 0\}$, so that $H_0^1(\omega_a)$ is the subspace of those $H^1(\omega_a)$ functions whose mean value vanishes; 2) $H_0(\text{curl}, \omega_a) := \{v \in \mathbf{H}(\text{curl}, \omega_a); v \times n_{\omega_a} = 0 \text{ on } \partial\omega_a\}$, where the tangential trace is understood by duality as above in Section 2.2; and, similarly, 3) $H_0(\text{div}, \omega_a) := \{v \in \mathbf{H}(\text{div}, \omega_a); v \cdot n_{\omega_a} = 0 \text{ on } \partial\omega_a\}$. We will also need 4) $H^1(\text{curl}, \omega_a) := \mathbf{H}(\text{curl}, \omega_a)$ (the symbol $\dagger$ is used here for notational purposes). The situation is more subtle for boundary vertices. As a first possibility, if $a \in \Gamma_N$ (i.e., $a \in \mathcal{V}_h$ is a boundary vertex such that all the faces sharing the vertex $a$ lie in $\Gamma_N$), then the spaces $H_0^1(\omega_a)$, $H_0(\text{curl}, \omega_a)$, $H_0(\text{div}, \omega_a)$, and $H^1(\text{curl}, \omega_a)$ are defined as above. Secondly, when $a \in \Gamma_D$, then at least one of the faces sharing the vertex $a$ lies in $\Gamma_D$; we denote by $\gamma_D$ the subset of $\Gamma_D$ formed by all mesh faces sharing the vertex $a$ and lying in $\Gamma_D$. In this situation, we let 1) $H_0^1(\omega_a) := \{v \in H^1(\omega_a); v = 0 \text{ on } \gamma_D\}$; 2) $H_0(\text{curl}, \omega_a) := \{v \in \mathbf{H}(\text{curl}, \omega_a); v \times n_{\omega_a} = 0 \text{ on } \partial\omega_a \setminus \gamma_D\}$; 3) $H_0(\text{div}, \omega_a) := \{v \in \mathbf{H}(\text{div}, \omega_a); v \cdot n_{\omega_a} = 0 \text{ on } \partial\omega_a \setminus \gamma_D\}$; and 4) $H^1(\text{curl}, \omega_a) := \{v \in \mathbf{H}(\text{curl}, \omega_a); v \times n_{\omega_a} = 0 \text{ on } \gamma_D\}$.

## 2.8 Functional inequalities

To work with data oscillation terms, we will employ the three following functional inequalities. First, from [12, Theorems 3.4 and 3.5], [28, Theorem 2.1], and the discussion in [10, Section 3.2.1], it follows that there exists a constant $C_L$ such that for all $v \in H_{0,D}(\text{curl}, \Omega)$, there exists $w \in H^1(\Omega) \cap H_{0,D}(\text{curl}, \Omega)$) such that $\nabla \times w = \nabla \times v$, and

$$ \|\nabla w\| \le C_L \|\nabla \times v\|. \quad (2.6) $$

When either $\Gamma_D$ or $\Gamma_N$ has zero measure and if $\Omega$ is convex, one can take $C_L = 1$, see [12] together with [25, Theorem 3.7] for Dirichlet boundary conditions and [25, Theorem 3.9] for Neumann boundary conditions.

Second, for any mesh element $K \in \mathcal{T}_h$ and $v \in H^1(K)$, there holds the Poincaré inequality

$$ \|v - \Pi_0(v)\|_K \le \frac{h_K}{\pi} \|\nabla v\|_K, \quad (2.7) $$
where $\Pi_0(v)$ denotes the componentwise mean value of $v$ on $K$.

Third, the Poincaré–Friedrichs–Weber inequality, see [22, Proposition 7.4] and more precisely [10, Theorem A.1] for the form of the constant, will be useful: for all vertices $a \in V_h$ and all vector-valued functions $v \in H^\dagger(\text{curl}, \omega_a) \cap H_0(\text{div}, \omega_a)$ with $\nabla \cdot v = 0$, we have

$$ \|v\|_{\omega_a} \lesssim h_{\omega_a} \|\nabla \times v\|_{\omega_a}. \qquad (2.8) $$

Strictly speaking, the inequality is established in [10, Theorem A.1] for edge patches, but the proof can be easily extended to vertex patches.

# 3 Setting

The purpose of this section is to introduce the curl-curl problem and its Nédélec finite element approximation. We also identify, in a form of two self-standing assumptions, the kernel properties solely needed for our analysis.

## 3.1 Current density

The following assumption is central for us:

**Assumption 3.1** (Current density $j$). Let $j$ be $H_{0,N}(\text{div}, \Omega)$-conforming, divergence-free, and $L^2(\Omega)$-orthogonal to the cohomology space $\mathcal{H}(\Omega, \Gamma_D)$, i.e.,

$$ j \in H_{0,N}(\text{div}, \Omega), \qquad (3.1a) $$

$$ \nabla \cdot j = 0, \qquad (3.1b) $$

$$ (j, \varphi) = 0 \quad \forall \varphi \in \mathcal{H}(\Omega, \Gamma_D). \qquad (3.1c) $$

Let us recall from Section 2.3 that when $\Omega$ is simply connected and $\Gamma_D$ is connected, then (3.1c) can be disregarded. Sometimes, to illustrate the main ideas, we will additionally suppose that $j$ is a piecewise $p$-degree Raviart–Thomas polynomial, $j \in \mathcal{RT}_p(T_h) \cap H_{0,N}(\text{div}, \Omega)$. Assumption 3.1 equivalently means that $j$ belongs to the range of the curl operator, i.e., there exists $v \in H_{0,N}(\text{curl}, \Omega)$ such that $\nabla \times v = j$.

## 3.2 The curl-curl problem

The curl-curl problem we study here reads: find the magnetic vector potential $A: \Omega \to \mathbb{R}^3$ such that

$$ \nabla \times (\nabla \times A) = j, \quad \nabla \cdot A = 0 \quad \text{in } \Omega, \qquad (3.2a) $$

$$ A \times n_{\Omega} = 0, \quad \text{on } \Gamma_D, \qquad (3.2b) $$

$$ (\nabla \times A) \times n_{\Omega} = 0, \quad A \cdot n_{\Omega} = 0, \quad \text{on } \Gamma_N, \qquad (3.2c) $$

with the additional requirement that $(A, \varphi) = 0$ for all $\varphi$ from the cohomology space $\mathcal{H}(\Omega, \Gamma_D)$ introduced in Section 2.3 to ensure uniqueness. Introducing $K(\Omega) := \{v \in H_{0,D}(\text{curl}, \Omega); \nabla \times v = 0\}$, the weak formulation of problem (3.2), cf., e.g., [6], consists in finding a pair $(A, q) \in H_{0,D}(\text{curl}, \Omega) \times K(\Omega)$ such that

$$ (A, \varphi) = 0 \quad \forall \varphi \in K(\Omega), \qquad (3.3a) $$

$$ (\nabla \times A, \nabla \times v) + (q, v) = (j, v) \quad \forall v \in H_{0,D}(\text{curl}, \Omega). \qquad (3.3b) $$

Picking the test function $v = q$ in (3.3b), we see that $q = 0$. Thus $A$ is such that

$$ A \in H_{0,D}(\text{curl}, \Omega) \qquad (3.4a) $$

$$ (\nabla \times A, \nabla \times v) = (j, v) \quad \forall v \in H_{0,D}(\text{curl}, \Omega). \qquad (3.4b) $$

**Remark 3.2** (Characterization (3.4) of the magnetic vector potential $A$). All the main developments below actually rely solely on (3.4), so that in particular the vector field $A$ can in our setting only be defined up to a curl-free component. Remark that the existence of $A$ satisfying (3.4) is a direct consequence of Assumption 3.1, and that a direct consequence of (3.4) is that $\nabla \times A \in H_{0,N}(\text{curl}, \Omega)$ with $\nabla \times (\nabla \times A) = j$.
## 3.3 Nédélec finite element approximation

For an integer $p \ge 0$ that we consider fixed henceforth, let the Nédélec finite element space be given by $V_h := \mathcal{N}_p(T_h) \cap H_{0,D}(\text{curl}, \Omega)$. The subspace $K_h := \{v_h \in V_h; \nabla \times v_h = 0\}$ is simply $\nabla(\mathcal{P}_{p+1}(T_h) \cap H_{0,D}^1(\Omega))$ when $\Omega$ is simply connected and $\Gamma_D$ is connected, and can be readily identified by introducing “cuts” in the mesh mimicking the construction of the cohomology space $\mathcal{H}(\Omega, \Gamma_D)$, see [26, Chapter 6]. The finite element approximation of (3.3) is a pair $(A_h, q_h) \in V_h \times K_h$ such that

$$ (A_h, \varphi_h) = 0 \quad \forall \varphi_h \in K_h, \qquad (3.5a) $$

$$ (\nabla \times A_h, \nabla \times v_h) + (q_h, v_h) = (j, v_h) \quad \forall v_h \in V_h. \qquad (3.5b) $$

Observing that $K_h \subset K$, this actually leads to $A_h \in V_h$ such that

$$ (\nabla \times A_h, \nabla \times v_h) = (j, v_h) \quad \forall v_h \in V_h. \qquad (3.6) $$

In the developments below, we can actually still weaken (3.6) and rely solely on:

**Assumption 3.3** (Discrete magnetic vector potential $A_h$). Let $A_h$ be a piecewise $p$-degree Nédélec polynomial satisfying a lowest-order Nédélec orthogonality property,

$$ A_h \in \mathcal{N}_p(T_h) \cap H_{0,D}(\text{curl}, \Omega), \qquad (3.7a) $$

$$ (\nabla \times A_h, \nabla \times v_h) = (j, v_h) \quad \forall v_h \in \mathcal{N}_0(T_h) \cap H_{0,D}(\text{curl}, \Omega). \qquad (3.7b) $$

# 4 Motivation

Let $j$ satisfy Assumption 3.1. We motivate here our approach by showing how an equilibrated flux $h$ may be constructed locally from any $A$ satisfying (3.4) at the continuous level. These observations are the basis of the actual flux equilibration procedure involving $A_h$ satisfying Assumption 3.3 at the discrete level that we develop in Sections 5.1 and 5.2 below. We would in particular like to identify a patchwise construction such that

$$ h^a \in H_0(\text{curl}, \omega_a), \qquad (4.1a) $$

$$ h := \sum_{a \in \mathcal{V}_h} h^a \in H_{0,N}(\text{curl}, \Omega), \qquad (4.1b) $$

$$ \nabla \times h = j. \qquad (4.1c) $$

At the continuous level, the solution is trivially

$$ h^a = \psi^a(\nabla \times A). \qquad (4.2) $$

We now rewrite the above definition implicitly. The idea is to introduce

$$ h^a := \arg \min_{\substack{v \in H_0(\text{curl}, \omega_a) \\ \nabla \times v = j^a}} \|v - \psi^a(\nabla \times A)\|_{\omega_a}^2 \quad \forall a \in \mathcal{V}_h \qquad (4.3) $$

with a suitable curl constraint $j^a$. Since

$$ \nabla \times (\psi^a(\nabla \times A)) = \psi^a \underbrace{(\nabla \times (\nabla \times A))}_{j} + \underbrace{\nabla \psi^a \times (\nabla \times A)}_{\theta^a} \qquad (4.4) $$

we have

$$ j^a := \psi^a j + \theta^a, \quad \theta^a := \nabla \psi^a \times (\nabla \times A). \qquad (4.5) $$

Importantly, it holds that

$$ \theta^a \in H_0(\text{div}, \omega_a), \qquad (4.6a) $$

$$ \nabla \cdot \theta^a = \underbrace{\nabla \times \nabla \psi^a}_{0} \cdot (\nabla \times A) - \nabla \psi^a \cdot \underbrace{\nabla \times (\nabla \times A)}_{j} = -\nabla \psi^a \cdot j, \qquad (4.6b) $$

$$ \sum_{a \in \mathcal{V}_h} \theta^a = \sum_{a \in \mathcal{V}_h} \nabla \psi^a \times (\nabla \times A) = 0, \qquad (4.6c) $$
where the last property follows by the partition of unity (2.1). Consequently,

$$j^a = \psi^a j + \theta^a \in H_0(\operatorname{div}, \omega_a), \qquad (4.7a)$$

$$\nabla \cdot j^a = \nabla \psi^a \cdot j + \psi^a \underbrace{\nabla \cdot j}_{0} + \nabla \cdot \theta^a = 0, \qquad (4.7b)$$

$$\sum_{a \in \mathcal{V}_h} j^a = j, \qquad (4.7c)$$

which gives a decomposition of the divergence-free current density j into divergence-free contributions j^a defined over the vertex patch subdomains ω_a. The above auxiliary fields θ^a can also be defined implicitly as the solution to the minimization problems:

$$\theta^a := \arg \min_{\substack{v \in H_0(\operatorname{div}, \omega_a) \\ \nabla \cdot v = -\nabla \psi^a \cdot j}} \|v - \nabla \psi^a \times (\nabla \times A)\|_{\omega_a}^2 \quad \forall a \in \mathcal{V}_h. \qquad (4.8)$$

We shall now mimic (4.3), (4.7), and (4.8) at the discrete level.

# 5 Main results

In this section, we summarize our main results.

## 5.1 Stable divergence-free patchwise decomposition of the given current density j

The central issue for our approach is a stable divergence-free patchwise decomposition of the current density j in the spirit of (4.7). For this purpose, we first design an appropriate discrete variant of (4.8), where we crucially rely on the patchwise orthogonality stemming from Assumption 3.3. We will initially be requested to work with the increased polynomial degree

$$p' := \min\{p, 1\}, \qquad (5.1)$$

recalling that $p \ge 0$ is fixed in Section 3.3. We start with:

**Definition 5.1** (Patchwise contributions j_h^a). Let j and A_h satisfy respectively Assumptions 3.1 and 3.3. Carry out the three following steps:

1. For all vertices a ∈ V_h, consider the p'-degree Raviart-Thomas patchwise (seemingly over-constrained) minimizations

$$\theta_h^a := \arg \min_{\substack{v_h \in \mathcal{R}\mathcal{T}_{p'}(\mathcal{T}_a) \cap \mathcal{H}_0(\operatorname{div}, \omega_a) \\ \nabla \cdot v_h = \Pi_{p'}(-\nabla \psi^a \cdot j) \\ (v_h, r_h)_K = (\nabla \psi^a \times (\nabla \times A_h), r_h)_K \quad \forall r_h \in [P_0(K)]^3, \forall K \in \mathcal{T}_a}} \|v_h - \nabla \psi^a \times (\nabla \times A_h)\|_{\omega_a}^2, \qquad (5.2)$$

where in addition to the usual normal trace and divergence, the constraints additionally also concern elementwise product with piecewise vector-valued constants.

2. Extending θ_h^a by zero outside of the patch subdomains ω_a, set

$$\delta_h := \sum_{a \in \mathcal{V}_h} \theta_h^a. \qquad (5.3)$$

For all tetrahedra K ∈ T_h, consider the (p + 1)-degree Raviart-Thomas elementwise minimizations:

$$\delta_h^a|_K := \arg \min_{\substack{v_h \in \mathcal{R}\mathcal{T}_1(K) \\ \nabla \cdot v_h = 0}} \|v_h - I_{K,1}^{\mathcal{R}\mathcal{T}}((\psi^a \delta_h)|_K)\|_K^2 \quad \forall a \in \mathcal{V}_K \text{ when } p=0, \qquad (5.4a)$$

$$v_h \cdot n_K = I_{K-1}^{\mathcal{R}\mathcal{T}}((\psi^a \delta_h)|_K) \cdot n_K \quad \text{on } \partial K$$

$$\delta_h^a|_K := \arg \min_{\substack{v_h \in \mathcal{R}\mathcal{T}_{p+1}(K) \\ \nabla \cdot v_h = 0 \\ v_h \cdot n_K = \psi^a \delta_h \cdot n_K \text{ on } \partial K}} \|v_h - \psi^a \delta_h\|_K^2 \quad \forall a \in \mathcal{V}_K \text{ when } p \ge 1, \qquad (5.4b)$$

which yields the divergence-free decomposition

$$\delta_h = \sum_{a \in \mathcal{V}_h} \delta_h^a.$$
3. For all vertices $\mathbf{a} \in \mathcal{V}_h$, define

$$j_h^a := \psi^a \mathbf{j} + \theta_h^a - \delta_h^a. \qquad (5.5)$$

For a vertex $\mathbf{a} \in \mathcal{V}_h$ and the extended (second-order) patch $\tilde{\mathcal{T}}_a$, define the data oscillation term

$$\tilde{\eta}_{\text{osc},j}^{a} := \left\{ \sum_{K \in \tilde{\mathcal{T}}_a} \left( \frac{h_K}{\pi} \| \mathbf{j} - \mathbf{\Pi}_{p'}(\mathbf{j}) \|_K \right)^2 \right\}^{1/2}. \qquad (5.6)$$

We crucially have:

**Theorem 5.2** (Stable divergence-free patchwise decomposition of $\mathbf{j}$). Let $\mathbf{j}$ and $\mathbf{A}_h$ satisfy respectively Assumptions 3.1 and 3.3. Let $j_h^a$ be given by Definition 5.1 for all vertices $\mathbf{a} \in \mathcal{V}_h$. Then

$$j_h^a \in \mathbf{H}_0(\text{div}, \omega_a), \qquad (5.7a)$$

$$\nabla \cdot j_h^a = \nabla \psi^a \cdot (\mathbf{j} - \mathbf{\Pi}_{p'}(\mathbf{j})), \qquad (5.7b)$$

$$\sum_{a \in \mathcal{V}_h} j_h^a = \mathbf{j}, \qquad (5.7c)$$

where the extension of $j_h^a$ by zero outside of the patch subdomains $\omega_a$ is understood in the last two properties. Moreover, when $\mathbf{j} \in \mathbf{RT}_p(\mathcal{T}_h) \cap \mathbf{H}_{0,N}(\text{div}, \Omega)$ is piecewise polynomial, then, in strengthening of (5.7a)-(5.7b),

$$j_h^a \in \mathbf{RT}_{p+1}(\mathcal{T}_a) \cap \mathbf{H}_0(\text{div}, \omega_a), \qquad (5.8a)$$

$$\nabla \cdot j_h^a = 0. \qquad (5.8b)$$

Let in addition $\mathbf{A}$ satisfying (3.4) be arbitrary and let, as in (4.5),

$$\mathbf{j}^a := \psi^a \mathbf{j} + \nabla \psi^a \times (\nabla \times \mathbf{A}).$$

Then

$$\| \mathbf{j}^a - j_h^a \|_{\omega_a} \lesssim h_{\omega_a}^{-1} [ \| \nabla \times (\mathbf{A} - \mathbf{A}_h) \|_{\bar{\omega}_a} + \tilde{\eta}_{\text{osc},j}^a ]. \qquad (5.9)$$

## Remarks

Several remarks about Definition 5.1 and Theorem 5.2 are in order:

1. At the discrete level, $\nabla \psi^a \times (\nabla \times \mathbf{A}_h) \notin \mathbf{H}_0(\text{div}, \omega_a)$, in contrast to $\nabla \psi^a \times (\nabla \times \mathbf{A})$, see (4.5)-(4.6). The auxiliary field $\theta_h^a$ from (5.2) is the projection of $\nabla \psi^a \times (\nabla \times \mathbf{A}_h)$ to $\mathbf{RT}_{p'}(\mathcal{T}_a) \cap \mathbf{H}_0(\text{div}, \omega_a)$ satisfying $\nabla \cdot \theta_h^a = \Pi_{p'}(-\nabla \psi^a \cdot \mathbf{j})$. Step 1 of Definition 5.1 thus mimics (4.8) and achieves equivalents to (4.6a) and (4.6b). Unfortunately, $\delta_h = \sum_{a \in \mathcal{V}_h} \theta_h^a$ given by (5.3) typically does not equal 0, which would mimic (4.6c).

2. Step 2 of Definition 5.1 yields the corrected fields $\theta_h^a - \delta_h^a$ which mimic (4.6) entirely in that (see Lemma 7.4 below for details)

$$\theta_h^a - \delta_h^a \in \mathbf{RT}_{p+1}(\mathcal{T}_a) \cap \mathbf{H}_0(\text{div}, \omega_a),$$

$$\nabla \cdot (\theta_h^a - \delta_h^a) = \Pi_{p'}(-\nabla \psi^a \cdot \mathbf{j}),$$

$$\sum_{a \in \mathcal{V}_h} (\theta_h^a - \delta_h^a) = \delta_h - \delta_h = 0.$$

3. Step 3 of Definition 5.1, in view of Theorem 5.2, finally materializes (4.7) at the discrete level.

4. Property (5.9) from Theorem 5.2 shows that the local discrete decomposition (5.5) compares in a $p$-robust way to the continuous-level decomposition (4.7), up to data oscillation given by (5.6).

5. Minimization (5.2) contains an additional constraint on the elementwise product with piecewise vector-valued constants. Existence, uniqueness, and $p$-robust stability theory for such problems is developed in Appendix A. Section 7 shows that this applies to our setting under the orthogonality in Assumption 3.3.

6. The additional constraint in (5.2) also ensures the existence, uniqueness, and $p$-robust stability of the elementwise problems (5.4), where $\delta_h^a$ form a divergence-free local decomposition of $\delta_h$ following Appendix B, see Lemma 7.4 below.
## 5.2 Equilibrated flux reconstruction based on local patchwise minimizations in $H(\text{curl})$

We now identify an appropriate discrete variant of (4.1)–(4.3), giving a locally computable equilibrated flux $h_h$. Let $j_h^a$ be given by Definition 5.1 and set

$$ \bar{j}_h^a := \arg \min_{\substack{v_h \in \mathcal{R}\mathcal{T}_{p+1}(\mathcal{T}_a) \cap \mathcal{H}_0(\text{div}, \omega_a) \\ \nabla \cdot v_h = 0 \\ (v_h, r_h)_K = (j_h^a, r_h)_K \quad \forall r_h \in [\mathcal{P}_0(\mathcal{K})]^3, \forall K \in \mathcal{T}_a}} \|v_h - j_h^a\|_{\omega_a}^2. \quad (5.10) $$

When $j$ is piecewise polynomial, $j \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_h) \cap \mathcal{H}_{0,N}(\text{div}, \Omega)$, then (5.8) implies that $\bar{j}_h^a = j_h^a$, so that there is no need for (5.10). In general, the role of (5.10) is to prepare a piecewise polynomial datum for a discrete version of (4.3): it projects the non-polynomial and non-divergence-free $j_h^a$ to $\bar{j}_h^a \in \mathcal{R}\mathcal{T}_{p+1}(\mathcal{T}_a) \cap \mathcal{H}_0(\text{div}, \omega_a)$ with $\nabla \cdot \bar{j}_h^a = 0$. Problem (5.10) has the same form as problem (5.2) and is well-posed following Appendix A below. With the Raviart–Thomas divergence-free $\bar{j}_h^a$, the following Nédélec local equilibration problem is well posed by standard arguments, see, e.g., [6]:

**Definition 5.3** (Equilibrated flux reconstruction based on local minimization in $H(\text{curl})$). Let $j$ and $A_h$ satisfy respectively Assumptions 3.1 and 3.3 and let, for all vertices $a \in \mathcal{V}_h$, $j_h^a$ be given by Definition 5.1 and $\bar{j}_h^a$ by (5.10). Consider the patchwise minimizations

$$ h_h^a := \arg \min_{\substack{v_h \in \mathcal{N}_{p+1}(\mathcal{T}_a) \cap \mathcal{H}_0(\text{curl}, \omega_a) \\ \nabla \times v_h = \bar{j}_h^a}} \|v_h - \psi^a(\nabla \times A_h)\|_{\omega_a}^2. \quad (5.11a) $$

Extending $h_h^a$ by zero outside of $\omega_a$, define

$$ h_h := \sum_{a \in \mathcal{V}_h} h_h^a. \quad (5.11b) $$

Recall $\tilde{\eta}_{\text{osc},j}^a$ from (5.6) and define

$$ \eta_{\text{osc},j}^a := h_{\omega_a} \|\bar{j}_h^a - j_h^a\|_{\omega_a}. \quad (5.12) $$

Crucially, the construction of Definition 5.3 is a stable equilibration:

**Theorem 5.4** (Equilibrium property and p-robust stability of the flux reconstruction). Let $j$ and $A_h$ satisfy respectively Assumptions 3.1 and 3.3. Then the equilibrated flux reconstruction $h_h$ from Definition 5.3 satisfies

$$ h_h \in \mathcal{N}_{p+1}(\mathcal{T}_h) \cap \mathcal{H}_{0,N}(\text{curl}, \Omega), \quad (5.13a) $$

$$ \nabla \times h_h = j \quad \text{when } j \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_h) \cap \mathcal{H}_{0,N}(\text{div}, \Omega). \quad (5.13b) $$

Let in addition $A$ satisfying (3.4) be arbitrary. Then

$$ \|h_h^a - \psi^a(\nabla \times A_h)\|_{\omega_a} \lesssim \|\nabla \times (A - A_h)\|_{\tilde{\omega}_a} + \eta_{\text{osc},j_h^a}^a + \tilde{\eta}_{\text{osc},j}^a. $$

## 5.3 Guaranteed, fully computable, constant-free, and p-robust a posteriori error estimates for the curl-curl problem

We apply here the results of Sections 5.1 and 5.2 to a posteriori error analysis of the curl-curl problem (3.2).

**Theorem 5.5** (Guaranteed, fully computable, and constant-free upper bound). Let $j$ satisfy Assumption 3.1, let $A$ be the weak solution to the curl-curl problem given by (3.3), and let $A_h$ be its Nédélec finite element approximation given by (3.5). Let $j_h^a$ be given by Definition 5.1 for all vertices $a \in \mathcal{V}_h$, and let $h_h$ be given by Definition 5.3. Then

$$ \|\nabla \times (A - A_h)\| \le \eta_{\text{tot}} := \underbrace{\|h_h - \nabla \times A_h\|}_{\eta} + C_L \underbrace{\left\{ \sum_{K \in \mathcal{T}_h} \underbrace{\frac{h_K^2}{\pi^2} \|j - \nabla \times h_h\|_K^2}_{\eta_{\text{osc},K}^2} \right\}}_{\eta_{\text{osc}}}^{1/2} $$

and

$$ \|h_h - \nabla \times A_h\|_K + \eta_{\text{osc},K} \lesssim \sum_{a \in \mathcal{V}_K} [\|\nabla \times (A - A_h)\|_{\tilde{\omega}_a} + \eta_{\text{osc},j_h^a}^a + \tilde{\eta}_{\text{osc},j}^a]. $$
# Remarks

Several remarks are in order:

1. On the discrete level, $\psi^a(\nabla \times A_h) \notin H_0(\text{curl}, \omega_a)$, in contrast to $\psi^a(\nabla \times A)$ on the continuous level, see (4.1)-(4.2). The equilibrated flux contribution $h_h^a$ from (5.11a) is its constrained projection to $\mathcal{N}_{p+1}(T_a) \cap H_0(\text{curl}, \omega_a)$. It mimics (4.3) at the discrete level.

2. When $j \in \mathcal{RT}_p(T_h) \cap H_{0,N}(\text{div}, \Omega)$, all the data oscillation estimators in Theorem 5.5 vanish. Indeed, (5.13b) implies that $\eta_{\text{osc}} = 0$, whereas $\tilde{j}_h^a = j_h^a$ gives $\eta_{\text{osc},j}^a = 0$, see (5.12). Similarly, $\tilde{\eta}_{\text{osc},j}^a$ from (5.6) vanishes as well (this is actually true up to $j \in \mathcal{RT}_{p'}(T_h) \cap H_{0,N}(\text{div}, \Omega)$, since $\nabla \cdot j = 0$). Moreover, all these terms are higher-order with respect to $\|\nabla \times (A - A_h)\|$ if $j$ is piecewise smooth. We also note that using (5.11b), (5.11a), and (5.7c), the data oscillation term $\eta_{\text{osc},K}$ can equivalently be rewritten with

$$j - \nabla \times h_h = \sum_{a \in \mathcal{V}_h} (j_h^a - \tilde{j}_h^a). \qquad (5.14)$$

3. The equilibration of Definition 5.3 is performed in local Nédélec spaces of order $p+1$. This is in agreement with p-robust flux equilibrations from [7, 20, 21, 11]. Similarly to [7, 19], it is also possible to design a downgrade of the orders of the local problems (5.11a) from $p+1$ to $p$.

Let us first discuss the case $p \ge 1$. The first step is to replace (5.4b) by (5.4a) with $\mathcal{RT}_1$ replaced by $\mathcal{RT}_p$. Then, according to Theorem B.1 with $q' = p$, we obtain $\delta_h^a \in \mathcal{RT}_p(T_a) \cap H_0(\text{div}, \omega_a)$ in place of (7.5a) below. Second, we employ (elementwise) the projector $I_{K,p}^{\mathcal{RT}}(\psi^a j)$ in (5.5) in place of $\psi^a j$. Let $j \in \mathcal{RT}_p(T_h) \cap H_{0,N}(\text{div}, \Omega)$. Then $\sum_{a \in \mathcal{V}_h} j_h^a = j$ and $\nabla \cdot j_h^a = 0$ as in (5.7c), (5.8b), but $j_h^a \in \mathcal{RT}_p(T_a) \cap H_0(\text{div}, \omega_a)$ in place of (5.8a). Consequently, (5.11a) can be brought down to

$$h_h^a := \arg \min_{v_h \in \mathcal{N}_p(T_a) \cap H_0(\text{curl}, \omega_a) \atop \nabla \times v_h = j_h^a} \|v_h - I_p^N(\psi^a(\nabla \times A_h))\|_{\omega_a}^2, \qquad (5.15)$$

where $I_p^N$ is the elementwise canonical p-degree Nédélec interpolate, analogue to (2.4). This leads to a cheaper procedure where the guaranteed estimate of Theorem 5.5 (with $\eta_{\text{osc}} = 0$) still holds true when $j \in \mathcal{RT}_p(T_h) \cap H_{0,N}(\text{div}, \Omega)$; general $j$ can be covered by data oscillation terms. Similarly, the local efficiency of Theorem 5.5 is also preserved, with, however, the p-robustness theoretically lost. In particular, from (B.6b), estimate (7.6) below still holds true, up to a possibly p-dependent constant.

Alternatively, for $p = 0$ in particular, because of $p' = p + 1$ employed in (5.2), we need to replace (5.5) by

$$j_h^a|_K := I_{K,0}^{\mathcal{RT}}((\psi^a j + \theta_h^a - \delta_h^a)|_K) \quad \forall K \in \mathcal{T}_h.$$

Let $j \in \mathcal{RT}_0(T_h) \cap H_{0,N}(\text{div}, \Omega)$. Then clearly $j_h^a \in \mathcal{RT}_0(T_a) \cap H_0(\text{div}, \omega_a)$. Moreover, from (2.5), $\nabla \cdot I_{K,0}^{\mathcal{RT}}(\delta_h^a|_K) = P_0(\nabla \cdot (\delta_h^a)|_K) = 0$, whereas $\nabla \cdot I_{K,0}^{\mathcal{RT}}((\psi^a j)|_K) = P_0(\nabla \cdot (\psi^a j)|_K) = (\nabla \psi^a \cdot j)|_K$, also using that $j \in [P_0(T_a)]^3$ as above in point 2, and similarly $\nabla \cdot I_{K,0}^{\mathcal{RT}}(\theta_h^a|_K) = P_0(\nabla \cdot (\theta_h^a)|_K) = (-\nabla \psi^a \cdot j)|_K$. Thus $\nabla \cdot j_h^a = 0$. Finally, $(\sum_{a \in \mathcal{V}_h} j_h^a)|_K = I_{K,0}^{\mathcal{RT}}(\sum_{a \in \mathcal{V}_h} (\psi^a j + \theta_h^a - \delta_h^a)|_K) = I_{K,0}^{\mathcal{RT}}(j|_K) = j|_K$ by the linearity of the Raviart-Thomas projector $\mathcal{RT}_0$. Then the above discussion for $p \ge 1$ applies.

4. The approach of [23, 24] includes solutions of local, a priori over-determined, problems on vertex patches in a multi-stage procedure. The present (again a priori over-constrained) problems (5.2) and consecutive steps in Definitions 5.1 and 5.3 share this spirit, though the minimizations directly determine the best-possible local energy error estimator contributions.

# 6 Numerical illustration

This section presents some numerical examples illustrating the key features of the estimator of Theorem 5.5. We impose the Dirichlet boundary condition on the whole boundary, i.e., $\Gamma_D := \partial \Omega$. We consider both structured meshes and unstructured meshes. When we speak about a “structured” mesh, we mean a Cartesian partition of $\Omega$ into $N \times N \times N$ cubes where each cube is first subdivided into 6 pyramids (with the basis a face and the apex the barycenter of the cube) and then each pyramid into $\sqrt{3}/(2N)$. On the other hand, the “unstructured” meshes are generated with the software pacakge MMG3D [16], where we simply require a maximum element size. These are typically quasi-uniform, but do not have any particular reapinging structure (every vertex patch is different). For both types of meshes, we consider the Nédélec finite element approximation (3.5) with varying degree $p \ge 1$.
<figure><img src="image_2.png" /><figcaption>Figure 1: [Smooth solution with limited regularity (6.1)] Uniform mesh refinement.</figcaption></figure>

## 6.1 $H^3(\Omega)$ solution with a polynomial right-hand side

We first consider the unit cube $\Omega := (0, 1)^3$ and a polynomial right-hand side $j := (0, 0, 1)$, so that the data oscillation estimator $\eta_{osc}$ vanishes. We can show that the solution is given by $A = (0, 0, A_3)$ with

$$A_3(\boldsymbol{x}) := \frac{16}{\pi^4} \sum_{n,m \ge 1} \frac{1}{nm(n^2 + m^2)} \sin(n\pi\boldsymbol{x}_1) \sin(m\pi\boldsymbol{x}_2). \quad (6.1)$$

This function belongs to $H^3(\Omega)$ but not to $H^4(\Omega)$. In practice, we cut the series at $n = m = 100$, and obtain $\nabla \times A$ by analytically differentiating (6.1).

We first fix the polynomial degree and consider a sequence of meshes. We use $p = 1$ and structured meshes with $N = 1, 2, 4, 8$, and then $p = 2$ and a sequence of unstructured meshes. Figure 1 presents the corresponding errors, estimates, and effectivity indices. We observe the expected convergence rate $h^2$ (recall that $A \in H^3(\Omega)$ merely). The estimator $\eta = \eta_{tot}$ ($\eta_{osc} = 0$ here) closely follows the error $\|\nabla \times (A - A_h)\|$, and the effectivity index given by the ratio $\eta/\|\nabla \times (A - A_h)\|$ is above but close to the optimal value 1; we actually numerically observe asymptotic exactness.

We then fix a mesh and increase the polynomial degree $p$ from 1 to 6. We consider two configurations: a structured mesh where the unit cube is split into 24 tetrahedra as described above and an unstructured mesh consisting of 20 tetrahedra. Figure 2 reports the results. The convergence is not exponential, which is expected because of the solution's finite regularity. Also in this setting, the estimator closely follows the actual error, and the effectivity index always remains above but close to 1. In particular, the effectivity index does not increase with $p$, which illustrates the $p$-robustness of the estimator.

Although this is not reported in the figures, we also numerically check that the reconstructed flux $h_h$ is indeed equilibrated, i.e., $\|j - \nabla \times h_h\| = 0$. Because of finite precision arithmetics, this value is not exactly zero, but ranges between $10^{-15}$ and $10^{-11}$, which is perfectly reasonable compared to the actual error levels.
<figure><img src="image_3.png" /><figcaption>Figure 2: [Smooth solution with limited regularity (6.1)] Uniform polynomial degree refinement.</figcaption></figure>

## 6.2 Analytical solution with a general right-hand side

We consider again the unit cube $\Omega := (0, 1)^3$, this time with a non-polynomial right-hand side $j := 8\pi^2(\sin(2\pi x_2) \sin(2\pi x_3), 0, 0)$. The associated solution is analytic,

$$A := (\sin(2\pi x_2) \sin(2\pi x_3), 0, 0). \qquad (6.2)$$

Figure 3 presents an $h$ convergence experiment with the same settings as above. The optimal convergence rate $h^{p+1}$ is observed for $||\nabla \times (A - A_h)||$. The oscillation-free estimator $\eta$ closely follows the actual error, with a possible slight underestimation, whereas the total estimator including data oscillation $\eta_{tot} = \eta + C_L \eta_{osc}$ from Theorem 5.5 gives a guaranteed upper bound with a slight overestimation; as discussed in Section 2.8, we can take here $C_L = 1$. In agreement with the theory, the influence of $\eta_{osc}$ diminishes with mesh refinement, and we again numerically observe asymptotic exactness. We then consider a $p$ convergence test. In Figure 4, we now observe the expected exponential convergence rate of the error and a perfect behavior of the effectivity indices. More precisely, as the mesh is not refined here, $\eta_{osc}$ does not necessarily go faster to zero than the error; this would be the case if the hp-version of (2.7), with $\Pi_p$ and $Ch_K/(p+1)$ in place of respectively $\Pi_0$ and $h_K/\pi$, was used. We put forward here (2.7), where the is no unknown constant $C$, leading to a fully computable $\eta_{osc}$.

## 6.3 Adaptivity with a singular solution

Our last experiment features a singular solution in a nonconvex domain, following [10, 23]. Specifically, we consider an L-shape example where $\Omega := L \times (0, 1)$, with

$$L := \{x = (r \cos \theta, r \sin \theta); |x_1|, |x_2| \le 1, \ 0 \le \theta \le 3\pi/2\}.$$

The right-hand side $j$ is non-polynomial and chosen such that

$$A(x) = (0, 0, \chi(r)r^{\alpha} \sin(\alpha\theta)), \qquad (6.3)$$
<figure><img src="image_4.png" /></figure>

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

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

<figure><img src="image_7.png" /><figcaption>Figure 3: [Analytical solution (6.2)] Uniform mesh refinement.</figcaption></figure>

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

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

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

<figure><img src="image_11.png" /><figcaption>Figure 4: [Analytical solution (6.2)] Uniform polynomial degree refinement.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 5: [Singular solution (6.3)] Adaptive mesh refinement.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 6: [Singular solution (6.3)]. Estimated (left) and actual (right) error distributions on the initial mesh. Top view (top) and side view (bottom).</figcaption></figure>

where $\alpha := 3/2$, $r^2 := |\mathbf{x}_1|^2 + |\mathbf{x}_2|^2$, $(\mathbf{x}_1, \mathbf{x}_2) = r(\cos \theta, \sin \theta)$, and $\chi : (0, 1) \to \mathbb{R}$ is a smooth cutoff function such that $\chi = 0$ in a neighborhood of 1. One easily checks that $\nabla \cdot \mathbf{A} = 0$. Besides, since $\Delta (r^\alpha \sin(\alpha\theta)) = 0$ near the origin, the right-hand side is non-singular (i.e., $j \in L^2(\Omega)$), and the singularity appearing in the solution is solely due to the re-entrant edge.

We couple our estimator with an adaptive strategy based on Dörfler's marking [17] for $\eta_K := \|\mathbf{h}_h - \nabla \times \mathbf{A}_h\|_K$ and MMG3D [16] to build a series of adaptively refined meshes. We select $p = 2$ and an initial mesh made of 415 elements.

The behaviors of the error and of the estimators $\eta$ and $\eta_{tot}$ with respect to the number of degrees of freedom $N_{dofs}$ are presented in Figure 5. Here we still take $C_L = 1$ in front of $\eta_{osc}$, though we do not anymore have a theoretical support for this. The effectivity indices stay close to one, even on unstructured and locally refined meshes, with $\eta_{tot}/\|\nabla \times (\mathbf{A} - \mathbf{A}_h)\|_K$ always above one. Besides, the optimal convergence rate is observed (it is limited to $-2/3$ when using isotropic elements in the presence of an edge singularity, see [5, Section 4.2.3]). This seems to indicate that the estimator is perfectly suited to drive adaptive mesh refinement, and illustrates the local efficiency of Theorem 5.5.

Finally, Figures 6–7 present the meshes generated by the adaptive algorithm, the estimators $\eta_K = \|\mathbf{h}_h - \nabla \times \mathbf{A}_h\|_K$, and the elementwise errors $\|\nabla \times (\mathbf{A} - \mathbf{A}_h)\|_K$ (the top face and the faces sharing the re-entrant edge). The meshes are refined close to the re-entrant edge, as expected. The estimated error distribution closely matches the actual one, illustrating again the local efficiency of the estimator.
<figure><img src="image_14.png" /><figcaption>Figure 7: [Singular solution (6.3)] Estimated (left) and actual (right) error distributions at adaptive mesh refinement iteration #10. Top view (top) and side view (bottom).</figcaption></figure>

# 7 Technical details and proofs

This section collects some technical details and the proofs of all the claims above.

## 7.1 Equivalent form of Assumption 3.3

Recall from Section 2.4 the piecewise affine “hat” function $\psi^a$ associated with the vertex $a \in \mathcal{V}_h$, as well as the notation $H^1_*(\omega_a)$ from Section 2.7. The following technical result holds true:

**Lemma 7.1** (Equivalence of images by the curl operator). There holds

$$ \nabla \times \left[ \operatorname{span} \left( \psi^a|_{\omega_a} \nabla (\mathcal{P}_1(\mathcal{T}_a) \cap H^1_*(\omega_a)) \right) \right] = \nabla \times [\mathcal{N}_0(\mathcal{T}_h) \cap \mathbf{H}_{0,D}(\operatorname{curl}, \Omega)]. \quad (7.1) $$

*Proof.* Let $a \in \mathcal{V}_h$. For any $q_h \in \mathcal{P}_1(\mathcal{T}_a) \cap H^1_*(\omega_a)$, clearly $\psi^a|_{\omega_a} \nabla q_h$, extended by zero outside of the patch subdomain $\omega_a$, lies in $\mathbf{H}_{0,D}(\operatorname{curl}, \Omega)$ (though in general not in $\mathcal{N}_0(\mathcal{T}_h)$). Moreover, $\nabla \times (\psi^a|_{\omega_a} \nabla q_h) = \nabla \psi^a|_{\omega_a} \nabla q_h$, which is a piecewise constant vector-valued polynomial on the patch $\mathcal{T}_a$ whose extension by zero outside of the patch subdomain $\omega_a$ has a continuous normal trace on interfaces and zero normal trace on $\Gamma_D$. Thus, this extension belongs to the lowest-order divergence-free Raviart-Thomas space, which implies $\nabla \times (\psi^a|_{\omega_a} \nabla q_h) = \nabla \times w_h$ on $\omega_a$ for $w_h$ which belongs to $\mathcal{N}_0(\mathcal{T}_h) \cap \mathbf{H}_{0,D}(\operatorname{curl}, \Omega)$. Thus, in (7.1), there holds the inclusion $\subseteq$.

Conversely, following, e.g., Monk [32, Section 5.5.1] or Ern and Guermond [18, Section 15.1], the space $\mathcal{N}_0(\mathcal{T}_h) \cap \mathbf{H}_{0,D}(\operatorname{curl}, \Omega)$ is spanned by the set of the “edge functions” $\{\psi^e\}_{e \in \mathcal{E}_h^D}$, where $\mathcal{E}_h^D$ denotes the mesh edges not lying in $\bar{\Gamma}_D$. If $e$ is the edge between vertices $a, b \in \mathcal{V}_h$, then $\psi^e = \psi^a \nabla \psi^b - \psi^b \nabla \psi^a$. Moreover, if one of the vertices of $e$ lies in $\bar{\Gamma}_D$, we chose the convention that $a \in \bar{\Gamma}_D$, so that we have $(\psi^b - c_b)|_{\omega_a} \in H^1_*(\omega_a)$ for some constant $c_b$ in all cases. Now, since $\nabla \times \psi^e = 2\nabla \psi^a \times \nabla \psi^b = 2\nabla \times (\psi^a \nabla \psi^b) = 2\nabla \times (\psi^a \nabla (\psi^b - c_b))$, we have found $q_h := (\psi^b - c_b)|_{\omega_a}/2 \in \mathcal{P}_1(\mathcal{T}_a) \cap H^1_*(\omega_a)$ such that, after zero extension, $\nabla \times (\psi^a|_{\omega_a} \nabla q_h) = \nabla \times \psi^e$, and the inclusion $\supseteq$ in (7.1) holds. $\square$

The following alternative formulation of Assumption 3.3 is crucial:

**Lemma 7.2** (Patchwise orthogonality). Let $j$ satisfy Assumption 3.1. Then $A_h$ satisfies Assumption 3.3 if and only if $A_h \in \mathcal{N}_p(\mathcal{T}_h) \cap \mathbf{H}_{0,D}(\operatorname{curl}, \Omega)$ and

$$ (\psi^a j, \nabla q_h)_{\omega_a} + (\nabla \psi^a \times (\nabla \times A_h), \nabla q_h)_{\omega_a} = 0 \quad \forall q_h \in \mathcal{P}_1(\mathcal{T}_a) \cap H^1_*(\omega_a), \forall a \in \mathcal{V}_h. \quad (7.2) $$
Proof. Since $\nabla \psi^a|_{\omega_a} \times \nabla q_h = \nabla \times (\psi^a|_{\omega_a} \nabla q_h)$,

$$(\nabla \psi^a \times (\nabla \times A_h), \nabla q_h)_{\omega_a} = -(\nabla \times A_h, \nabla \psi^a \times \nabla q_h)_{\omega_a} = -(\nabla \times A_h, \nabla \times (\psi^a \nabla q_h))_{\omega_a}.$$

For any $v \in H_{0,N}(\text{curl}, \Omega)$ such that $j = \nabla \times v$, the Green theorem in turn gives

$$(\psi^a j, \nabla q_h)_{\omega_a} = (j, \psi^a \nabla q_h)_{\omega_a} = (v, \nabla \times (\psi^a \nabla q_h))_{\omega_a}.$$

Finally, again by the Green theorem, for any $v_h \in \mathcal{N}_0(\mathcal{T}_h) \cap H_{0,D}(\text{curl}, \Omega)$,

$$(j, v_h) = (\nabla \times v, v_h) = (v, \nabla \times v_h).$$

Applying these identities respectively in (7.2) and (3.7b), the assertion follows from Lemma 7.1. $\square$

## 7.2 Properties of the auxiliary fields $\theta_h^a$, $\delta_h$, and $\delta_h^a$ from Definition 5.1

We collect here some important results on $\theta_h^a$, $\delta_h$, and $\delta_h^a$ from (5.2)–(5.4). We start with the following application of the self-standing result on over-constrained minimization in the Raviart–Thomas spaces that we present in Appendix A below. Let $\eta_{\text{osc},j}^a$ be defined as $\tilde{\eta}_{\text{osc},j}^a$ in (5.6) but on the patch $\mathcal{T}_a$ only.

**Lemma 7.3** (Existence, uniqueness, and stability of $\theta_h^a$ from (5.2)). There exists a unique solution $\theta_h^a$ to problem (5.2) for all $a \in \mathcal{V}_h$. Moreover, it satisfies the stability estimate

$$\|\theta_h^a - \nabla \psi^a \times (\nabla \times A_h)\|_{\omega_a} \lesssim \min_{\substack{v \in H_0(\text{div}, \omega_a) \\ \nabla \cdot v = -\nabla \psi^a \cdot j}} \|v - \nabla \psi^a \times (\nabla \times A_h)\|_{\omega_a} + h_{\omega_a}^{-1} \eta_{\text{osc},j}^a.$$

Proof. We choose $g^a := (-\nabla \psi^a \cdot j)|_{\omega_a}$, $\tau_h^a := (\nabla \psi^a \times (\nabla \times A_h))|_{\omega_a}$, $q := p$ and verify the assumptions of Theorem A.2 in three steps. Note that $\Pi_{p'}(\nabla \psi^a \cdot j) = \nabla \psi^a \cdot \Pi_{p'}(j)$ and that $\|\nabla \psi^a \cdot (j - \Pi_{p'}(j))\|_K \le \|\nabla \psi^a\|_{\infty, \omega_a} \|j - \Pi_{p'}(j)\|_K \lesssim h_{\omega_a}^{-1} \|j - \Pi_{p'}(j)\|_K$, where $\|\nabla \psi^a\|_{\infty, \omega_a} \lesssim h_{\omega_a}^{-1}$ follows from the shape regularity of the mesh, which gives rise to $h_{\omega_a}^{-1} \eta_{\text{osc},j}^a$ from the data oscillation term in Theorem A.2.

**Step 1.** Assumption (A.1a). From (3.1a), $g^a \in L^2(\omega_a)$, so that the first condition in (A.1a) is satisfied. From (3.7a), in turn, on $\omega_a$, it follows that $\nabla \times A_h \in [P_p(\mathcal{T}_a)]^3$, see, e.g., [6, Corollary 2.3.2], so that $\tau_h^a = \nabla \psi^a \times (\nabla \times A_h) \in [P_p(\mathcal{T}_a)]^3 \subset \mathcal{R}T_p(\mathcal{T}_a) \subset \mathcal{R}T_{p'}(\mathcal{T}_a)$. Thus the second (polynomial) condition in (A.1a) is also satisfied.

**Step 2.** Assumption (A.1b). For vertices $a \in \mathcal{V}_h$ such that $a \notin \bar{\Gamma}_D$, the Green theorem and $j \in H_{0,N}(\text{div}, \Omega)$ from (3.1a) together with $\nabla \cdot j = 0$ from (3.1b) imply

$$-(\nabla \psi^a \cdot j, 1)_{\omega_a} = -(\nabla \psi^a, j)_{\omega_a} = (\psi^a, \nabla \cdot j)_{\omega_a} = 0.$$

**Step 3.** Assumption (A.1c). For any $q_h \in P_1(\mathcal{T}_a) \cap H_*^1(\omega_a)$, again the Green theorem yields

$$-(\nabla \psi^a \cdot j, q_h)_{\omega_a} \stackrel{(3.1b)}{=} -(\nabla \cdot (\psi^a j), q_h)_{\omega_a} = (\psi^a j, \nabla q_h)_{\omega_a},$$

so that the patchwise orthogonality property (7.2) implies

$$(\nabla \psi^a \times (\nabla \times A_h), \nabla q_h)_{\omega_a} - (\nabla \psi^a \cdot j, q_h)_{\omega_a} = 0. \qquad (7.3)$$

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

Similarly, an important part of the results of the following lemma are consequences of Appendix B below:

**Lemma 7.4** (Auxiliary correction fields $\delta_h$ and $\delta_h^a$). For $\delta_h$ given by (5.3), there holds

$$\delta_h \in \mathcal{R}T_{p'}(\mathcal{T}_h) \cap H_{0,N}(\text{div}, \Omega) \quad \text{and} \quad \nabla \cdot \delta_h = 0. \qquad (7.4)$$

In addition, there exists a unique solution $\delta_h^a|_K$ to problems (5.4) for all tetrahedra $K \in \mathcal{T}_h$ and all vertices $a \in \mathcal{V}_K$, yielding the local divergence-free decomposition

$$\delta_h^a \in \mathcal{R}T_{p+1}(\mathcal{T}_a) \cap H_0(\text{div}, \omega_a) \quad \text{and} \quad \nabla \cdot \delta_h^a = 0 \quad \forall a \in \mathcal{V}_h, \qquad (7.5a)$$

$$\delta_h = \sum_{a \in \mathcal{V}_h} \delta_h^a. \qquad (7.5b)$$

Moreover, for all tetrahedra $K \in \mathcal{T}_h$ and all vertices $a \in \mathcal{V}_K$, there holds the local stability estimate

$$\|\delta_h^a\|_K \lesssim \|\delta_h\|_K. \qquad (7.6)$$
Proof. The patchwise contributions $\theta_h^a$ extended by zero outside of the patch subdomains $\omega_a$ belong to $RT_{p'}(T_h) \cap H_{0,N}(\text{div}, \Omega)$, so that the first property in (7.4) is immediate. The second property in (7.4) then follows by the divergence constraint in (5.2), the linearity of the projector $\Pi_{p'}$, and the partition of unity (2.1), since

$$ \nabla \cdot \delta_h = \sum_{a \in \mathcal{V}_h} \nabla \cdot \theta_h^a = \sum_{a \in \mathcal{V}_h} \Pi_{p'}(-\nabla \psi^a \cdot j) = \Pi_{p'} \left[ \sum_{a \in \mathcal{V}_h} -\nabla \psi^a \cdot j \right] = \Pi_{p'}(0) = 0. $$

Let $K \in \mathcal{T}_h$ and $r_h \in [\mathcal{P}_0(K)]^3$ be fixed. Then definition (5.3), which gives $\delta_h|_K = \sum_{b \in \mathcal{V}_K} \theta_h^b$, the partition of unity (2.1), which implies $\sum_{b \in \mathcal{V}_K} (\nabla \psi^b \times (\nabla \times A_h))|_K = 0$, and the elementwise orthogonality constraint in (5.2) lead to

$$ (\delta_h, r_h)_K = \sum_{b \in \mathcal{V}_K} (\theta_h^b - \nabla \psi^b \times (\nabla \times A_h), r_h)_K = 0. $$

This is condition (B.2). Thus, Theorem B.1 can be employed, where we choose $q := p'$ together with $q' := p'$ for $p = 0$ and $q' := p' + 1$ for $p \ge 1$. This implies the existence and uniqueness of solutions $\delta_h^a|_K$ to problems (5.4), the properties (7.5a), the decomposition (7.5b), and the stability bound (7.6). Note in particular that we only employ (B.6b) with $q' = q$ in the lowest-order case with $q = 1$, whereas in other cases, we employ (B.6b) with $q' = q+1$, so there is indeed no polynomial degree dependence in (7.6). $\square$

## 7.3 Decomposition of the current density j and its stability from Theorem 5.2

We are now ready to prove Theorem 5.2.

Proof of Theorem 5.2 (decomposition). Property (5.7a) is immediate since $\psi^a j \in H_0(\text{div}, \omega_a)$ in view of assumption (3.1a), from (5.2) which gives $\theta_h^a \in RT_{p'}(T_a) \cap H_0(\text{div}, \omega_a)$, and from the first property in (7.5a). Property (5.7b) follows since $\nabla \cdot (\psi^a j) = \nabla \psi^a \cdot j$ in view of assumption (3.1b) and using $\nabla \cdot \theta_h^a = \Pi_{p'}(-\nabla \psi^a \cdot j) = -\nabla \psi^a \cdot \Pi_{p'}(j)$ from (5.2) and $\nabla \cdot \delta_h^a = 0$, which is the second property in (7.5a). Finally, (5.7c) follows from the partition of unity (2.1) which gives $\sum_{a \in \mathcal{V}_h} \psi^a j = j$ together with (5.3) and (7.5b). When $j \in RT_p(T_h) \cap H_{0,N}(\text{div}, \Omega)$, (5.8) immediately follows from (5.7) and the fact that $\psi^a j \in RT_{p+1}(T_a) \cap H_0(\text{div}, \omega_a)$. $\square$

Proof of Theorem 5.2 (stability). We develop

$$ j^a - j_h^a = \nabla \psi^a \times (\nabla \times A) - \theta_h^a + \delta_h^a = \nabla \psi^a \times (\nabla \times (A - A_h)) - (\theta_h^a - \nabla \psi^a \times (\nabla \times A_h)) + \delta_h^a. $$

For the first term above, we immediately see

$$ \| \nabla \psi^a \times (\nabla \times (A - A_h)) \|_{\omega_a} \le \| \nabla \psi^a \|_{\infty, \omega_a} \| \nabla \times (A - A_h) \|_{\omega_a}. $$

For the second term above, we employ Lemma 7.3 with $v = \nabla \psi^a \times (\nabla \times A)$, which lies in $H_0(\text{div}, \omega_a)$ with divergence equal to $-\nabla \psi^a \cdot j$ by virtue of (4.6), which leads to

$$ \| \theta_h^a - \nabla \psi^a \times (\nabla \times A_h) \|_{\omega_a} \lesssim \| \nabla \psi^a \times (\nabla \times (A - A_h)) \|_{\omega_a} + h_{\omega_a}^{-1} \eta_{\text{osc},j}. \quad (7.7) $$

For the last term, we first recall (7.6), i.e., $\|\delta_h^a\|_K \lesssim \|\delta_h\|_K$ for every $K \in \mathcal{T}_a$. Now definition (5.3), the partition of unity (2.1), and the triangle inequality imply

$$ \| \delta_h \|_K = \left\| \sum_{b \in \mathcal{V}_K} (\theta_h^b - \nabla \psi^b \times (\nabla \times A_h)) \right\|_K \le \sum_{b \in \mathcal{V}_K} \| \theta_h^b - \nabla \psi^b \times (\nabla \times A_h) \|_{\omega_b}, $$

which extends by one layer beyond the patch $\omega_a$ and can be bounded by (7.7). The shape regularity of the mesh ensures that $\|\nabla \psi^a\|_{\infty, \omega_a} \lesssim h_{\omega_a}^{-1}$ and $\|\nabla \psi^b\|_{\infty, \omega_b} \simeq \|\nabla \psi^a\|_{\infty, \omega_a}$ for all vertices $b$ in the patch $\mathcal{T}_a$. Hence, (5.9) follows upon combining the above developments. $\square$
## 7.4 Equilibrated flux reconstruction from Section 5.2 and its stability

To prove Theorem 5.4, we rely on the following crucial result:

**Theorem 7.5** (p-robust $H(\text{curl}, \omega_a)$ stability). For a vertex $\mathbf{a} \in \mathcal{V}_h$, let $\mathbf{A}_h \in \mathcal{N}_p(\mathcal{T}_a) \cap \mathbf{H}(\text{curl}, \omega_a)$ and $\bar{\mathbf{j}}_h^a \in \mathcal{RT}_{p+1}(\mathcal{T}_a) \cap \mathbf{H}_0(\text{div}, \omega_a)$ with $\nabla \cdot \bar{\mathbf{j}}_h^a = 0$ be given. Then

$$ v_h \in \mathcal{N}_{p+1}(\mathcal{T}_a) \cap \mathbf{H}_0(\text{curl}, \omega_a) \qquad \min_{\nabla \times v_h = \bar{\mathbf{j}}_h^a} \|v_h - \psi^a(\nabla \times \mathbf{A}_h)\|_{\omega_a} \lesssim \min_{\nabla \times v = \bar{\mathbf{j}}_h^a} \|v - \psi^a(\nabla \times \mathbf{A}_h)\|_{\omega_a}. \quad (7.8) $$

On a single tetrahedron $K$ in place of the vertex patch $\mathcal{T}_a$, Theorem 7.5 follows by the seminal contributions of Costabel and McIntosh [13, Proposition 4.2] and Demkowicz et al. [14, Theorem 7.2], see [9, Theorem 2]. On an edge patch, such a result has been established in [10, Theorem 3.1]. The further extension to a vertex patch has recently been established in [11, Theorem 3.3, see also Corollary 4.3].

*Proof of Theorem 5.4 (equilibration).* Property (5.13a) follows immediately from $h_h^a \in \mathcal{N}_{p+1}(\mathcal{T}_a) \cap \mathbf{H}_0(\text{curl}, \omega_a)$ of (5.11a) and (5.11b). For piecewise polynomial $j \in \mathcal{RT}_p(\mathcal{T}_h) \cap \mathbf{H}_{0,N}(\text{div}, \Omega)$, $\nabla \times h_h^a = \bar{\mathbf{j}}_h^a = \mathbf{j}_h^a$ from (5.11a) and (5.10). Property (5.13b) is then an direct consequence of (5.7c) and (5.11b). $\square$

*Proof of Theorem 5.4 (stability).* Fix a vertex $\mathbf{a} \in \mathcal{V}_h$ and use $\mathbf{j}^a = \psi^a \mathbf{j} + \nabla \psi^a \times (\nabla \times \mathbf{A}) = \nabla \times (\psi^a(\nabla \times \mathbf{A}))$ as in property (4.4). This implies $(\mathbf{j}^a, v)_{\omega_a} = (\psi^a(\nabla \times \mathbf{A}), \nabla \times v)_{\omega_a}$ for any $v \in \mathbf{H}^\dagger(\text{curl}, \omega_a)$. Then Theorem 7.5 and a primal-dual equivalence as in, e.g., [10, Lemma 5.5] imply

$$ \begin{aligned} \|\mathbf{h}_h^a - \psi^a(\nabla \times \mathbf{A}_h)\|_{\omega_a} &\lesssim \min_{\nabla \times v = \bar{\mathbf{j}}_h^a} \|v - \psi^a(\nabla \times \mathbf{A}_h)\|_{\omega_a} \\ &= \sup_{\substack{v \in \mathcal{H}^\dagger(\text{curl}, \omega_a) \\ \|\nabla \times v\|_{\omega_a} = 1}} \{ (\bar{\mathbf{j}}_h^a, v)_{\omega_a} - (\psi^a(\nabla \times \mathbf{A}_h), \nabla \times v)_{\omega_a} \} \\ &\le \sup_{\substack{v \in \mathcal{H}^\dagger(\text{curl}, \omega_a) \\ \|\nabla \times v\|_{\omega_a} = 1}} (\bar{\mathbf{j}}_h^a - \mathbf{j}^a, v)_{\omega_a} + \|\psi^a(\nabla \times (\mathbf{A} - \mathbf{A}_h))\|_{\omega_a} \\ &\le \sup_{\substack{v \in \mathcal{H}^\dagger(\text{curl}, \omega_a) \\ \|\nabla \times v\|_{\omega_a} = 1}} (\bar{\mathbf{j}}_h^a - \mathbf{j}^a, v)_{\omega_a} + \|\nabla \times (\mathbf{A} - \mathbf{A}_h)\|_{\omega_a}. \end{aligned} $$

We are thus left to treat the first term above.

Fix $v \in \mathbf{H}^\dagger(\text{curl}, \omega_a)$ with $\|\nabla \times v\|_{\omega_a} = 1$. Consider $q \in H_*^1(\omega_a)$ such that $(\nabla q, \nabla w)_{\omega_a} = (v, \nabla w)_{\omega_a} \quad \forall w \in H_*^1(\omega_a)$.

Then $\tilde{v} := v - \nabla q$ lies in both $\mathbf{H}^\dagger(\text{curl}, \omega_a)$ and $\mathbf{H}_0(\text{div}, \omega_a)$ and is divergence-free, $\nabla \cdot \tilde{v} = 0$. Thus, the Poincaré-Friedrichs-Weber inequality (2.8) implies

$$ \|\tilde{v}\|_{\omega_a} \lesssim h_{\omega_a} \|\nabla \times \tilde{v}\|_{\omega_a} = h_{\omega_a} \|\nabla \times v\|_{\omega_a} = h_{\omega_a}. \quad (7.9) $$

Note that $\bar{\mathbf{j}}_h^a - \mathbf{j}^a \in \mathbf{H}_0(\text{div}, \omega_a)$ with $\nabla \cdot (\bar{\mathbf{j}}_h^a - \mathbf{j}^a) = 0$; indeed, this follows from (4.7a)-(4.7b) together with (5.10). Thus, the Green theorem gives

$$ (\bar{\mathbf{j}}_h^a - \mathbf{j}^a, \nabla q)_{\omega_a} = 0. \quad (7.10) $$

Thanks to this crucial property, we can play in $\tilde{v}$ and use (7.9): employing additionally the Cauchy-Schwarz inequality and the triangle inequality, we have

$$ \begin{aligned} (\bar{\mathbf{j}}_h^a - \mathbf{j}^a, v)_{\omega_a} &= (\bar{\mathbf{j}}_h^a - \mathbf{j}^a, \tilde{v})_{\omega_a} \le \|\bar{\mathbf{j}}_h^a - \mathbf{j}^a\|_{\omega_a} \|\tilde{v}\|_{\omega_a} \\ &\lesssim h_{\omega_a} [\|\bar{\mathbf{j}}_h^a - \mathbf{j}_h^a\|_{\omega_a} + \|\mathbf{j}_h^a - \mathbf{j}^a\|_{\omega_a}], \end{aligned} \quad (7.11) $$

and we conclude by (5.9) from Theorem 5.2. $\square$
## 7.5 A posteriori error estimates from Section 5.3

We can finally prove Theorem 5.5.

*Proof of Theorem 5.5 (reliability).* For a piecewise polynomial current density, $j \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_h) \cap \mathbf{H}_{0,N}(\text{div}, \Omega)$, Theorem 5.4 implies $h_h \in \mathbf{H}_{0,N}(\text{curl}, \Omega)$ with $\nabla \times h_h = j$. Thus, in this case the claim follows with $\eta_{\text{osc}} = 0$ by the Prager-Synge theorem [35] in the $\mathbf{H}(\text{curl})$-context, see, e.g., [8, Theorem 10] or [23, Theorem 3.1].

In general, we proceed as follows. Since $\mathbf{A}, \mathbf{A}_h \in \mathbf{H}_{0,D}(\text{curl}, \Omega)$,

$$ \|\nabla \times (\mathbf{A} - \mathbf{A}_h)\| = \max_{\mathbf{v} \in \mathbf{H}_{0,D}(\text{curl}, \Omega) \atop \|\nabla \times \mathbf{v}\| = 1} (\nabla \times (\mathbf{A} - \mathbf{A}_h), \nabla \times \mathbf{v}). $$

Fix $\mathbf{v} \in \mathbf{H}_{0,D}(\text{curl}, \Omega)$ with $\|\nabla \times \mathbf{v}\| = 1$ and consider $\mathbf{w}$ from (2.6). Note that since $h_h \in \mathbf{H}_{0,N}(\text{curl}, \Omega)$ from Theorem 5.4, the Green theorem and $\nabla \times \mathbf{w} = \nabla \times \mathbf{v}$ give

$$ (\nabla \times h_h, \mathbf{w}) = (h_h, \nabla \times \mathbf{w}) = (h_h, \nabla \times \mathbf{v}). $$

Similarly, $\nabla \times \mathbf{w} = \nabla \times \mathbf{v}$ and the weak solution characterization (3.4) lead to

$$ (\nabla \times \mathbf{A}, \nabla \times \mathbf{v}) = (\nabla \times \mathbf{A}, \nabla \times \mathbf{w}) = (j, \mathbf{w}). $$

Thus

$$ (\nabla \times (\mathbf{A} - \mathbf{A}_h), \nabla \times \mathbf{v}) = (j - \nabla \times h_h, \mathbf{w}) + (h_h - \nabla \times \mathbf{A}_h, \nabla \times \mathbf{v}). $$

The second term is trivially bounded by the estimator $\eta$ via the Cauchy-Schwarz inequality, so that we are left with bounding the first one.

Property (5.14) and the additional orthogonality constraint in (5.10) lead to

$$ (j - \nabla \times h_h, \mathbf{w}) = \sum_{K \in \mathcal{T}_h} \left( \sum_{a \in \mathcal{V}_K} (j_h^a - \bar{j}_h^a), \mathbf{w} \right)_K = \sum_{K \in \mathcal{T}_h} (j - \nabla \times h_h, \mathbf{w} - \mathbf{\Pi}_0(\mathbf{w}))_K. $$

Consequently, the Poincaré inequality (2.7), (2.6), and $\|\nabla \times \mathbf{v}\| = 1$ give

$$ (j - \nabla \times h_h, \mathbf{w}) \le \sum_{K \in \mathcal{T}_h} \eta_{\text{osc},K} \|\nabla \mathbf{w}\|_K \le \eta_{\text{osc}} \|\nabla \mathbf{w}\| \le C_L \eta_{\text{osc}} \|\nabla \times \mathbf{v}\| = C_L \eta_{\text{osc}}. $$

□

**Remark 7.6 (Comparison with (7.11)).** Above, we could also write

$$ (j - \nabla \times h_h, \mathbf{w}) = \sum_{a \in \mathcal{V}_h} (j_h^a - \bar{j}_h^a, \mathbf{w})_{\omega_a}, $$

where the terms in the sum are similar to (7.11) from Section 7.4. In contrast to (7.11), it seems that we cannot pass through the Poincaré-Friedrichs-Weber inequality (2.8) in absence of an exactly divergence-free field (recall from (5.7b) that $\nabla \cdot j_h^a = \nabla \psi^a \cdot (j - \mathbf{\Pi}'_p(j))$ only in general), and rather need to resort to the switch from $\mathbf{v} \in \mathbf{H}_{0,D}(\text{curl}, \Omega)$ to $\mathbf{w}$ of (2.6) and to make use of the Poincaré inequality (2.7).

*Proof of Theorem 5.5 (efficiency).* Property (5.14), the triangle inequality, and definition (5.12) immediately lead to $\eta_{\text{osc},K} \le \sum_{a \in \mathcal{V}_K} \eta_{\text{osc},j_h^a}$. Moreover, Definition (5.11b) together with the partition of unity (2.1) imply

$$ \|h_h - \nabla \times \mathbf{A}_h\|_K = \left\| \sum_{a \in \mathcal{V}_K} (h_h^a - \psi^a(\nabla \times \mathbf{A}_h)) \right\|_K \le \sum_{a \in \mathcal{V}_K} \|h_h^a - \psi^a(\nabla \times \mathbf{A}_h)\|_{\omega_a}. $$

Thus employing Theorem 5.4 concludes the proof. □
# A Over-constrained minimization in Raviart–Thomas’ spaces

## A.1 Assumption and statement of the over-constrained minimization

In this appendix, we consider a fixed mesh vertex $\mathbf{a} \in \mathcal{V}_h$. Let an integer $q \ge 0$ be fixed and set $q' := \min\{q, 1\}$. We employ the notation of Section 2 and in particular recall that $\lesssim$ means smaller or equal to up to a constant only depending on the mesh shape-regularity parameter $\kappa_{\mathcal{T}_h}$. We also assume a polynomial form, mean value zero, and patchwise orthogonality conditions on the two data $g^a$ and $\tau_h^a$:

**Assumption A.1 (Data $g^a$ and $\tau_h^a$).** The data $g^a$ and $\tau_h^a$ satisfy

$$g^a \in L^2(\omega_a) \quad \text{and} \quad \tau_h^a \in \mathcal{RT}_{q'}(\mathcal{T}_a), \qquad (A.1a)$$

$$(g^a, 1)_{\omega_a} = 0 \quad \text{when } \mathbf{a} \notin \bar{\Gamma}_D, \qquad (A.1b)$$

$$(\tau_h^a, \nabla q_h)_{\omega_a} + (g^a, q_h)_{\omega_a} = 0 \quad \forall q_h \in \mathcal{P}_1(\mathcal{T}_a) \cap H_*^1(\omega_a). \qquad (A.1c)$$

We consider the following (seemingly over-constrained) minimization problem in the Raviart–Thomas space $\mathcal{RT}_{q'}(\mathcal{T}_a) \cap H_0(\text{div}, \omega_a)$:

$$\theta_h^a := \arg \min_{\substack{v_h \in \mathcal{RT}_{q'}(\mathcal{T}_a) \cap H_0(\text{div}, \omega_a) \\ \nabla \cdot v_h = \Pi_{q'}(g^a) \\ (v_h, r_h)_K = (\tau_h^a, r_h)_K \quad \forall r_h \in [\mathcal{P}_0(K)]^3, \forall K \in \mathcal{T}_a}} \|v_h - \tau_h^a\|_{\omega_a}^2. \qquad (A.2)$$

The following result is of independent interest:

**Theorem A.2 (Over-constrained minimization in the Raviart–Thomas spaces).** Let Assumption A.1 hold. Then there exists a unique solution $\theta_h^a$ to problem (A.2), satisfying the stability estimate

$$\|\theta_h^a - \tau_h^a\|_{\omega_a} \lesssim \min_{\substack{v \in H_0(\text{div}, \omega_a) \\ \nabla \cdot v = g^a}} \|v - \tau_h^a\|_{\omega_a} + \left\{ \sum_{K \in \mathcal{T}_a} \left( \frac{h_K}{\pi} \|\Pi_{q'}(g^a) - g^a\|_K \right)^2 \right\}^{1/2}.$$

## A.2 Auxiliary conventional minimization

In addition to (A.2), it will be useful to also consider

$$\bar{\theta}_h^a := \arg \min_{\substack{v_h \in \mathcal{RT}_{q'}(\mathcal{T}_a) \cap H_0(\text{div}, \omega_a) \\ \nabla \cdot v_h = \Pi_{q'}(g^a)}} \|v_h - \tau_h^a\|_{\omega_a}^2. \qquad (A.3)$$

Minimizations (A.3) are in a conventional format in that the constraints only concern normal trace and divergence. Moreover, they fulfill the following important property:

**Lemma A.3 (Existence, uniqueness, and stability of $\bar{\theta}_h^a$ from (A.3)).** Let Assumption A.1 hold. Then there exists a unique solution $\bar{\theta}_h^a$ to problem (A.3), satisfying the stability estimate

$$\|\bar{\theta}_h^a - \tau_h^a\|_{\omega_a} \lesssim \min_{\substack{v \in H_0(\text{div}, \omega_a) \\ \nabla \cdot v = \Pi_{q'}(g^a)}} \|v - \tau_h^a\|_{\omega_a} + \left\{ \sum_{K \in \mathcal{T}_a} \left( \frac{h_K}{\pi} \|\Pi_{q'}(g^a) - g^a\|_K \right)^2 \right\}^{1/2}. \qquad (A.4)$$

*Proof.* Existence and uniqueness of $\bar{\theta}_h^a$ from (A.3) are classical following, e.g., [6], thanks to the Neumann boundary compatibility condition (A.1b); note that this implies $(\Pi_{q'}(g^a), 1)_{\omega_a} = 0$ when $\mathbf{a} \notin \bar{\Gamma}_D$. Moreover, since $\Pi_{q'}(g^a) \in \mathcal{P}_{q'}(\mathcal{T}_a)$ and $\tau_h^a \in \mathcal{RT}_{q'}(\mathcal{T}_a)$, taking $p = q'$, $\tau_p = -\tau_h^a$, and $r_K = (\Pi_{q'}(g^a) - \nabla \cdot \tau_h^a)|_K$ in [21, Corollaries 3.3 and 3.6] for an interior vertex $\mathbf{a}$ and [21, Corollary 3.8] and [11, Proposition 3.1] for a boundary vertex $\mathbf{a}$ lead to

$$\|\bar{\theta}_h^a - \tau_h^a\|_{\omega_a} \lesssim \min_{\substack{v \in H_0(\text{div}, \omega_a) \\ \nabla \cdot v = \Pi_{q'}(g^a)}} \|v - \tau_h^a\|_{\omega_a} = \|\nabla \tilde{r}^a\|_{\omega_a}.$$

The equality above is a classical primal-dual equivalence, with $\tilde{r}^a \in H_*^1(\omega_a)$ given by

$$(\nabla \tilde{r}^a, \nabla v)_{\omega_a} = (\tau_h^a, \nabla v)_{\omega_a} + (\Pi_{q'}(g^a), v)_{\omega_a} \quad \forall v \in H_*^1(\omega_a).$$
Thus, as in, e.g., [20, Theorem 3.17],

$$
\begin{align*}
\|\nabla \tilde{r}^a\|_{\omega_a} &= \max_{\substack{v \in H_h^1(\omega_a) \\ \|\nabla v\|_{\omega_a} = 1}} \{(\tau_h^a, \nabla v)_{\omega_a} + (\Pi_{q'}(g^a), v)_{\omega_a}\} \\
&= \max_{\substack{v \in H_*^1(\omega_a) \\ \|\nabla v\|_{\omega_a} = 1}} \{(\tau_h^a, \nabla v)_{\omega_a} + (g^a, v)_{\omega_a} + (\Pi_{q'}(g^a) - g^a, v)_{\omega_a}\}.
\end{align*}
$$

The projection orthogonality and the elementwise Poincaré inequality then lead to

$$
\begin{aligned}
| (\Pi_{q'}(g^a) - g^a, v)_{\omega_a} | &= \left| \sum_{K \in \mathcal{T}_a} (\Pi_{q'}(g^a) - g^a, v - \Pi_0 v)_K \right| \\
&\le \left\{ \sum_{K \in \mathcal{T}_a} \left( \frac{h_K}{\pi} \|\Pi_{q'}(g^a) - g^a\|_K \right)^2 \right\}^{1/2} \|\nabla v\|_{\omega_a},
\end{aligned}
$$

and (A.4) follows, since

$$
\max_{\substack{v \in H_h^1(\omega_a) \\ \|\nabla v\|_{\omega_a} = 1}} \{(\tau_h^a, \nabla v)_{\omega_a} + (g^a, v)_{\omega_a}\} = \min_{\substack{v \in H_0(\text{div}, \omega_a) \\ \nabla \cdot v = g^a}} \|v - \tau_h^a\|_{\omega_a}
$$

by the same primal-dual equivalence argument.

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

A.3 Auxiliary first-order over-constrained minimization and proof of Theo-
rem A.2

Let, in addition to (A.2) and (A.3), the first-order Raviart–Thomas piecewise polynomial $\bar{\epsilon}_h^a$ be given by

$$
\begin{equation}
\begin{aligned}
\bar{\epsilon}_h^a &:= \operatorname*{arg\,min}_{\substack{v_h \in \mathcal{R}\mathcal{T}_1(\bar{\tau}_a) \cap \mathcal{H}_0(\operatorname{div}, \omega_a) \\ \nabla \cdot v_h = 0}} \|v_h - \tau_h^a + \bar{\theta}_h^a\|_{\omega_a}^2. \quad (A.5) \\
(v_h, r_h)_K &= (\tau_h^a - \bar{\theta}_h^a, r_h)_K \quad \forall r_h \in [\mathcal{P}_0(K)]^3, \forall K \in \mathcal{T}_a
\end{aligned}
\end{equation}
$$

The field $\bar{\epsilon}_h^a$ can be seen as the correction of $\bar{\theta}_h^a$ from (A.3) necessary to fulfill the constraints on the elementwise product with piecewise vector-valued constants in (A.2). As one might expect, the patchwise orthogonality assumption (A.1c), turns to be the key for the following crucial technical result:

**Lemma A.4** (Existence, uniqueness, and stability of $\bar{\epsilon}_h^a$ from (A.5)). Let Assumption A.1 hold. Then there exists a unique solution $\bar{\epsilon}_h^a$ to problem (A.5), and the following stability estimate holds true:

$$
\| \bar{\epsilon}_h^a \|_{\omega_a} \lesssim \| \tau_h^a - \bar{\theta}_h^a \|_{\omega_a}. \qquad (A.6)
$$

We postpone the proof Lemma A.4 to the sections below; let us now first show that Lemma A.4
implies the results announced in Theorem A.2.

*Proof of Theorem A.2.* It follows straightforwardly from (A.3) and (A.5) that $\bar{\theta}_h^a + \bar{\epsilon}_h^a$ lies in the mini-
mization set of (A.2). Consequently, the existence and uniqueness of (A.2) follows since the minimized
functional in (A.2) is convex. Moreover, the triangle inequality together with Lemma A.4 implies

$$
\| \theta_h^a - \tau_h^a \|_{\omega_a} \le \| \bar{\theta}_h^a + \bar{\epsilon}_h^a - \tau_h^a \|_{\omega_a} \le \| \bar{\epsilon}_h^a \|_{\omega_a} + \| \bar{\theta}_h^a - \tau_h^a \|_{\omega_a} \lesssim \| \bar{\theta}_h^a - \tau_h^a \|_{\omega_a},
$$

and we conclude by Lemma A.3.

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

A.4 Piola mappings

In order to prove the technical results below, we will rely on Piola mappings. An extensive description
may be found in, e.g., [18, Chapters 7.2 and 9.2], and we only list here the essential properties we need.

If $U, V \subset \mathbb{R}^3$ are open sets with Lipschitz boundaries and $\phi: U \to V$ is a bilipschitz mapping, the gradient-, curl-, divergence-preserving and broken Piola mappings are the applications $\phi^g: L^2(U) \to L^2(V)$, $\phi^c, \phi^d: L^2(U) \to L^2(V)$, and $\phi^b: L^2(U) \to L^2(V)$ respectively defined by

$$
\phi^g(v) = v \circ \phi^{-1}, \quad \phi^c(w) = (\mathbb{J}^{-T} w) \circ \phi^{-1}, \quad \phi^d(w) = \left( \frac{\mathbb{J}}{|\mathbb{J}|} w \right) \circ \phi^{-1}, \quad \phi^b(v) = (|\mathbb{J}|v) \circ \phi^{-1}, \quad (A.7)
$$
<figure><img src="image_18.png" /><figcaption>Figure 8: Interior vertex patch (left) and the element $K_1$ with the face $F_1^{\text{ext}}$ (right)</figcaption></figure>

for all $v \in L^2(U)$ and $w \in L^2(V)$, where $\mathbb{J}$ is the Jacobian matrix of $\phi$ and $|\mathbb{J}|$ its determinant. These mappings are invertible. In addition, if $\gamma_U \subset \partial U$ and $\gamma_V := \phi(V)$, with a similar notation as in Section 2.2, $\phi^g$, $\phi^c$, and $\phi^d$

$$H_{0,\gamma_U}^1(U) \to H_{0,\gamma_V}^1(V), \quad \mathbf{H}_{0,\gamma_U}(\text{curl}, U) \to \mathbf{H}_{0,\gamma_V}(\text{curl}, V), \quad \mathbf{H}_{0,\gamma_U}(\text{div}, U) \to \mathbf{H}_{0,\gamma_V}(\text{div}, V) \quad (A.8)$$

are bijections, and more generally, the full, tangential, and normal traces on $\gamma_U$ are respectively transported by $\phi^g$, $\phi^c$, and $\phi^d$ on $\gamma_V$. The commutativity properties

$$\phi^c(\nabla v) = \nabla(\phi^g(v)) \quad \nabla(\phi^d(w)) = \phi^b(\nabla \cdot w) \quad (A.9)$$

for $v \in H^1(U)$ and $w \in \mathbf{H}(\text{div}, U)$ will be useful. We will also need the formula

$$((\phi^d)^{-1}(\mathbf{v}), \mathbf{w})_U = \varepsilon(\mathbf{v}, \phi^c(\mathbf{w}))_V \quad \forall \mathbf{v} \in L^2(V), \forall \mathbf{w} \in L^2(U), \quad (A.10)$$

where $\varepsilon$ is the (constant) sign of the determinant of $\mathbb{J}$. Finally, if $U$ is a polyhedron covered by a tetrahedral mesh $\mathcal{T}_U$ and $\phi|_K$ is affine for each $K \in \mathcal{T}_U$, denoting by $\mathcal{T}_V$ the tetrahedral mesh of $V$ induced by $\phi$, we have the bijections

$$\phi^g : \mathcal{P}_q(\mathcal{T}_U) \to \mathcal{P}_q(\mathcal{T}_V), \quad \phi^c : \mathcal{N}_q(\mathcal{T}_U) \to \mathcal{N}_q(\mathcal{T}_V), \quad \phi^d : \mathcal{R}\mathcal{T}_q(\mathcal{T}_U) \to \mathcal{R}\mathcal{T}_q(\mathcal{T}_V) \quad (A.11)$$

for all integer $q \ge 0$. In addition, in this case, we have

$$\|\phi^d\|^3 \|\phi^d\|^{-1} \le C(\kappa_{\mathcal{T}_U}, \kappa_{\mathcal{T}_V}) \quad (A.12)$$

with $\|\phi^d\|$ denoting the norm operator of $\phi^d : L^2(U) \to L^2(V)$ (and vice-versa for $\|(\phi^d)^{-1}\|$), and $\kappa_{\mathcal{T}_U}$, $\kappa_{\mathcal{T}_V}$ the shape-regularity constants of $\mathcal{T}_U$ and $\mathcal{T}_V$ as in Section 2.4.

## A.5 A preliminary result

Before proving Lemma A.4, we establish the following preliminary result:

**Lemma A.5 (Orthogonalities).** Let $\mu^a \in L^2(\omega_a)$ satisfy $(\mu^a, \nabla q_h)_{\omega_a} = 0$ for all $q_h \in \mathcal{P}_1(\mathcal{T}_a) \cap H_*(\omega_a)$. Then, the following set is non-empty:

$$W_h(\mathcal{T}_a, \mu^a) := \left\{ v_h \in \mathcal{R}\mathcal{T}_1(\mathcal{T}_a) \cap \mathbf{H}_0(\text{div}, \omega_a) \;\middle|\; \begin{array}{l} \nabla \cdot \mathbf{v}_h = 0 \\ (\mathbf{v}_h, \mathbf{r}_h)_K = (\mu^a, \mathbf{r}_h)_K \\ \forall \mathbf{r}_h \in [\mathcal{P}_0(K)]^3, \forall K \in \mathcal{T}_a \end{array} \right\}.$$

*Proof. Step 1: interior patches.* We start with the case where the vertex $a \in \mathcal{V}_h$ does not lie on the boundary $\partial\Omega$, cf. Figure 8, left. We will construct a particular element $w_h \in W_h(\mathcal{T}_a, \mu^a)$ by an explicit run through the patch $\mathcal{T}_a$ of tetrahedra sharing the vertex $a \in \mathcal{V}_h$, similarly as in [7, 21]. Specifically,
following the concept of shelling of a polytopal complex, see [38, Theorem 8.12] and [21, Lemma B.1], there exists an enumeration $K_i$, $1 \le i \le |\mathcal{T}_a|$, of the tetrahedra in the patch $\mathcal{T}_a$ such that, except for the first tetrahedron in the enumeration $K_1$: (i) if there are at least two faces corresponding to the neighbors of $K_i$ which have been already enumerated, then all the tetrahedra of $\mathcal{T}_a$ sharing this edge come sooner in the enumeration; (ii) except for the last element $K_{|\mathcal{T}_a|}$, there are one or two neighbors of $K_i$ which have been already enumerated and correspondingly two or one neighbors of $K_i$ which have not been enumerated yet.

Consider a pass through the patch $\mathcal{T}_a$ in the sense of the above enumeration. For the tetrahedron $K_i$, $1 \le i \le |\mathcal{T}_a|$, let us denote by $\mathcal{F}_i^\sharp$ the faces of $K_i$ corresponding to the neighbors of $K_i$ which have been already passed through and $F^j = \partial K_i \cap \partial K_j \in \mathcal{F}_i^\sharp$ the face corresponding to the neighbor $K_j$. Also, let $F_i^{\text{ext}}$ be the face of $K_i$ lying on the patch boundary $\partial\omega_a$. Consider the problem

$$ w_h^i := \arg \min_{ \substack{ \mathbf{v}_h \in \mathcal{R}\mathcal{T}_1(K_i) \\ \nabla \cdot \mathbf{v}_h = 0 \\ \mathbf{v}_h \cdot \mathbf{n}_{K_i} = 0 \text{ on } F_i^{\text{ext}} \\ \mathbf{v}_h \cdot \mathbf{n}_{K_i} = w_h^j \cdot \mathbf{n}_{K_i} \text{ on all } F^j \in \mathcal{F}_i^\sharp \\ (\mathbf{v}_h, \mathbf{r}_h)_{K_i} = (\boldsymbol{\mu}^a, \mathbf{r}_h)_{K_i} \quad \forall \mathbf{r}_h \in [\mathcal{P}_0(K_i))^3 } } \|\mathbf{v}_h - \boldsymbol{\mu}^a\|_{K_i}^2, \quad (A.13) $$

similar to (A.5) but reduced to the single tetrahedron $K_i$. If a solution to (A.13) exists, $w_h$ defined as $w_h|_{K_i} := w_h^i$ belongs to $W_h(\mathcal{T}_a, \boldsymbol{\mu}^a)$. We are thus left to establish the existence and uniqueness of (A.13).

Step 1a: the first element $K_1$. Let us start with the first element $K_1$, cf. Figure 8, right. Then the set $\mathcal{F}_1^\sharp$ is empty, and the constraints in (A.13) lead us to ask whether in the first-order Raviart-Thomas space $\mathcal{R}\mathcal{T}_1(K_1)$, one can impose simultaneously the divergence, the normal flux through one face, and moments against constant functions. We will reason by the canonical degrees of freedom, see, e.g., [6, Proposition 2.3.4 and Figure 2.14.c] or (2.4), and find a suitable $\mathbf{v}_h \in \mathcal{R}\mathcal{T}_1(K_1)$ satisfying the constraints in (A.13). First, we see that in $\mathcal{R}\mathcal{T}_1(K_1)$, the normal flux $\mathbf{v}_h \cdot \mathbf{n}_{K_1}$ on $F_1^{\text{ext}}$ can be fixed to zero and the moments against constants $(\mathbf{v}_h, \mathbf{r}_h)_{K_1}$ can be fixed as in (A.13). We still have the freedom to choose the normal fluxes $\mathbf{v}_h \cdot \mathbf{n}_{K_1}$ on the faces of $K_1$ different from $F_1^{\text{ext}}$, and the question is whether this can be done so as to fix the divergence of $\mathbf{v}_h$ to zero. By [6, Proposition 2.3.3], there holds

$$ \nabla \cdot \mathbf{v}_h = 0 \Leftrightarrow (\nabla \cdot \mathbf{v}_h, q_h)_{K_1} = 0 \quad \forall q_h \in \mathcal{P}_1(K_1). $$

Employing the Green theorem and the fact that $\mathbf{v}_h \cdot \mathbf{n}_{K_1} = 0$ on $F_1^{\text{ext}}$,

$$ (\nabla \cdot \mathbf{v}_h, q_h)_{K_1} = \langle \mathbf{v}_h \cdot \mathbf{n}_{K_1}, q_h \rangle_{\partial K_1 \setminus F_1^{\text{ext}}} - (\mathbf{v}_h, \nabla q_h)_{K_1}. $$

Now, since $\nabla q_h \in [\mathcal{P}_0(K_1)]^3$, the last term above is fixed from the last constraint in (A.13), so the question becomes: can one choose $\mathbf{v}_h \cdot \mathbf{n}_{K_1}$ on $\partial K_1 \setminus F_1^{\text{ext}}$ such that

$$ \langle \mathbf{v}_h \cdot \mathbf{n}_{K_1}, q_h \rangle_{\partial K_1 \setminus F_1^{\text{ext}}} = (\boldsymbol{\mu}^a, \nabla q_h)_{K_1} \quad \forall q_h \in \mathcal{P}_1(K_1), \quad (A.14) $$

which gives 4 conditions for the 9 remaining degrees of freedom (there are 4 degrees of freedom in $\mathcal{P}_1(K_1)$ and 3 degrees of freedom per face in $\mathcal{R}\mathcal{T}_1(K_1)$ following [6, Proposition 2.3.4]).

We proceed as follows. Out of the three faces of $K_1$ different from $F_1^{\text{ext}}$, choose one and impose $\mathbf{v}_h \cdot \mathbf{n}_{K_1} = 0$ therein. Then we are left to set $\mathbf{v}_h \cdot \mathbf{n}_{K_1}$ on two faces, say $F$ and $\tilde{F}$. For $F$, consider the three hat basis functions $\psi_F^k$, $1 \le k \le 3$, as in Section 2.4, corresponding to its three vertices. Restricted to $\tilde{F}$, which is necessary a face neighboring with $F$, they belong to $\mathcal{P}_1(\tilde{F})$, and one of the restrictions, say $\psi_F^3$, is zero on $\tilde{F}$. Thus, there holds

$$ \langle \mathbf{v}_h \cdot \mathbf{n}_{K_1}, \psi_F^3 \rangle_{\tilde{F}} = 0, $$

and, following [6, Proposition 2.3.4], we can prescribe

$$ \langle \mathbf{v}_h \cdot \mathbf{n}_{K_1}, \psi_F^k \rangle_{\tilde{F}} := 0 \quad 1 \le k \le 2. $$

Moreover, restricted to $F$, $\psi_F^k$ create a basis of $\mathcal{P}_1(F)$, whereas restricted to $K_1$, they belong to $\mathcal{P}_1(K_1)$. Thus, we can also set

$$ \langle \mathbf{v}_h \cdot \mathbf{n}_{K_1}, \psi_F^k \rangle_F := (\boldsymbol{\mu}^a, \nabla \psi_F^k)_{K_1} \quad 1 \le k \le 3. $$

With the choices made so far, we see that (A.14) holds for the three hat functions $\psi_F^k$, $1 \le k \le 3$. Finally, consider $\psi_F^4$, the hat basis function corresponding to the vertex opposite to the face $F$. Restricted to $F$, it is zero, so that

$$ \langle \mathbf{v}_h \cdot \mathbf{n}_{K_1}, \psi_F^4 \rangle_F = 0. $$
Moreover, restricted to $\tilde{F}$, it completes $\psi_F^1$ and $\psi_F^2$ (restricted to $\tilde{F}$) to create a basis of $\mathcal{P}_1(\tilde{F})$, and
restricted to $K_1$, it belongs to $\mathcal{P}_1(K_1)$, so that we can prescribe

$$ \langle v_h \cdot n_{K_1}, \psi_F^4 \rangle_{\tilde{F}} := (\mu^a, \nabla \psi_F^4)_{K_1}. $$

Thus, (A.14) also holds for $\psi_F^4$, and since $\psi_F^k$, $1 \le k \le 4$, restricted to $K_1$ create a basis of $\mathcal{P}_1(K_1)$, (A.14) holds true, and a unique $w_h^1$ from (A.13) exists.

Step 1b: any element $K_i$ with $|\mathcal{F}_i^\sharp| = 1$. We now investigate those consecutive elements $K_i$ which are such that two neighbors of $K_i$ have not been passed through yet. This means that exactly one neighbor of $K_i$, say $K_j$, has already been passed through, so there is one face $F^j$ in the set $\mathcal{F}_i^\sharp$. Since $v_h \cdot n_{K_i} = w_h^j \cdot n_{K_i}$ on $F^j$ is requested in (A.13), (A.14) asks if one can choose $v_h \cdot n_{K_i}$ on $\partial K_i \setminus \{F_i^{\text{ext}}, F^j\}$ such that

$$ \langle v_h \cdot n_{K_i}, q_h \rangle_{\partial K_i \setminus \{F_i^{\text{ext}}, F^j\}} = (\mu^a, \nabla q_h)_{K_i} - \langle w_h^j \cdot n_{K_i}, q_h \rangle_{F^j} \quad (A.15) $$

for all $q_h \in \mathcal{P}_1(K_i)$, which is still undetermined, giving 4 conditions for the 6 remaining degrees of freedom. The reasoning is similar as for $K_1$. Still denoting $F$ and $\tilde{F}$ the two remaining faces and $\psi_F^k$, $1 \le k \le 4$ the hat basis functions, we again have

$$ \langle v_h \cdot n_{K_i}, \psi_F^3 \rangle_{\tilde{F}} = 0, \quad \langle v_h \cdot n_{K_i}, \psi_F^4 \rangle_F = 0. $$

Moreover, imposing

$$ \begin{align*} \langle v_h \cdot n_{K_i}, \psi_F^k \rangle_{\tilde{F}} &:= 0 \quad 1 \le k \le 2, \\ \langle v_h \cdot n_{K_i}, \psi_F^4 \rangle_{\tilde{F}} &:= (\mu^a, \nabla \psi_F^4)_{K_i} - \langle w_h^j \cdot n_{K_i}, \psi_F^4 \rangle_{F^j}, \\ \langle v_h \cdot n_{K_i}, \psi_F^k \rangle_F &:= (\mu^a, \nabla \psi_F^k)_{K_i} - \langle w_h^j \cdot n_{K_i}, \psi_F^k \rangle_{F^j} \quad 1 \le k \le 3 \end{align*} $$

yields (A.15), and $w_h^i$ exists.

Step 1c: any element $K_i$ with $|\mathcal{F}_i^\sharp| = 2$. We now investigate those consecutive elements $K_i$ which are such that only one neighbor of $K_i$ has not been passed through yet, with $K_{j1}$ and $K_{j2}$ already passed through and faces $F^{j1}$, $F^{j2}$ in the set $\mathcal{F}_i^\sharp$. Denote $F$ the only remaining face, so that $F_i^{\text{ext}}$, $F^{j1}$, $F^{j2}$, and $F$ are the four faces of the tetrahedron $K_i$. As in (A.14) and (A.15), we need to ensure that

$$ \langle v_h \cdot n_{K_i}, q_h \rangle_F = (\mu^a, \nabla q_h)_{K_i} - \langle w_h^{j1} \cdot n_{K_i}, q_h \rangle_{F^{j1}} - \langle w_h^{j2} \cdot n_{K_i}, q_h \rangle_{F^{j2}} \quad (A.16) $$

for all $q_h \in \mathcal{P}_1(K_1)$. This time, the system is over-determined in that we request 4 conditions for the 3 remaining degrees of freedom of the normal components $v_h \cdot n_{K_i}$ on the face $F$. As above, we can impose

$$ \langle v_h \cdot n_{K_i}, \psi_F^k \rangle_F := (\mu^a, \nabla \psi_F^k)_{K_i} - \langle w_h^{j1} \cdot n_{K_i}, \psi_F^k \rangle_{F^{j1}} - \langle w_h^{j2} \cdot n_{K_i}, \psi_F^k \rangle_{F^{j2}} \quad 1 \le k \le 3, $$

which fixes $v_h \cdot n_{K_i}$ on the face $F$. Now, noting that $\langle v_h \cdot n_{K_i}, \psi_F^4 \rangle_F = 0$, it follows that to prove (A.16), we need to show that

$$ (\mu^a, \nabla \psi_F^4)_{K_i} - \langle w_h^{j1} \cdot n_{K_i}, \psi_F^4 \rangle_{F^{j1}} - \langle w_h^{j2} \cdot n_{K_i}, \psi_F^4 \rangle_{F^{j2}} = 0. \quad (A.17) $$

To prove (A.17), recall from property (i) of the enumeration (giving that all other elements sharing the edge $e$ common to $F^{j1}$ and $F^{j2}$ have been already passed through) and the previous steps, see (A.14) and (A.15), that

$$ \langle w_h^j \cdot n_{K_j}, \psi_F^4 \rangle_{\partial K_j} = (\mu^a, \nabla \psi_F^4)_{K_j} \quad (A.18) $$

for all the tetrahedra $K_j$ sharing the edge $e$, different from $K_i$. Recalling from assumptions of Lemma A.5, we have

$$ 0 = (\mu^a, \nabla \psi_F^4)_{\omega_a} = (\mu^a, \nabla \psi_F^4)_{\omega_e}, \quad (A.19) $$

where $\omega_e$ is the part of $\omega_a$ corresponding to the elements sharing the edge $e$; the last equality holds since in the vertex patch subdomain $\omega_a$, $\psi_F^4$ is only supported on the edge patch subdomain $\omega_e$. Denote by $\tilde{\omega}_e$ the part of $\omega_e$ without the element $K_i$. Then the normal traces orientation, the zero normal trace boundary conditions $w_h^j \cdot n_{K_j} = 0$ on the faces $F_j^{\text{ext}}$ together with the zero values of $\psi_F^4$ on $\partial \omega_e \setminus \partial \omega_a$, the Green theorem first applied on $\tilde{\omega}_e$ and later individually on $K_j$, and the notation $w_h|_{K_j} = w_h^j$ for the
previous $K_j^o \subset \tilde{\omega}_e$ give

$$
\begin{align}
& - \langle \boldsymbol{w}_h^{j_1} \cdot \boldsymbol{n}_{K_i}, \psi_F^4 \rangle_{F^{j_1}} - \langle \boldsymbol{w}_h^{j_2} \cdot \boldsymbol{n}_{K_i}, \psi_F^4 \rangle_{F^{j_2}} = \langle \boldsymbol{w}_h^{j_1} \cdot \boldsymbol{n}_{\tilde{\omega}_e}, \psi_F^4 \rangle_{F^{j_1}} + \langle \boldsymbol{w}_h^{j_2} \cdot \boldsymbol{n}_{\tilde{\omega}_e}, \psi_F^4 \rangle_{F^{j_2}} \nonumber \\
& = \langle \boldsymbol{w}_h \cdot \boldsymbol{n}_{\tilde{\omega}_e}, \psi_F^4 \rangle_{\partial \tilde{\omega}_e} = (\boldsymbol{w}_h, \nabla \psi_F^4)_{\tilde{\omega}_e} + (\nabla \cdot \boldsymbol{w}_h, \psi_F^4)_{\tilde{\omega}_e} \nonumber \\
& = \sum_{j; K_j^o \subset \tilde{\omega}_e} \{ (\boldsymbol{w}_h^j, \nabla \psi_F^4)_{K_j} + (\nabla \cdot \boldsymbol{w}_h^j, \psi_F^4)_{K_j} \} = \sum_{j; K_j^o \subset \tilde{\omega}_e} \langle \boldsymbol{w}_h^j \cdot \boldsymbol{n}_{K_j}, \psi_F^4 \rangle_{\partial K_j} \tag{A.20} \\
& \stackrel{(A.18)}{=} \sum_{j; K_j^o \subset \tilde{\omega}_e} (\boldsymbol{\mu}^a, \nabla \psi_F^4)_{K_j} = (\boldsymbol{\mu}^a, \nabla \psi_F^4)_{\omega_e} - (\boldsymbol{\mu}^a, \nabla \psi_F^4)_{K_i} \stackrel{(A.19)}{=} -(\boldsymbol{\mu}^a, \nabla \psi_F^4)_{K_i}, \nonumber
\end{align}
$$

which is (A.17). Thus, there exists a unique $w_h^j$ from (A.13) also on this $K_i$.

Step 1d: the last element $K_{|\mathcal{T}_a|}$. According to property (ii) of the enumeration, the last element $K_{|\mathcal{T}_a|}$ is such that $|\mathcal{F}_{|\mathcal{T}_a|}^\sharp| = 3$, so that all the neighbors have been already passed through. Consequently, all the degrees of freedom of $v_h$ are fixed from the last three constraints in (A.13), and we need to show that $\nabla \cdot v_h = 0$, i.e., that

$$ (\nabla \cdot v_h, q_h)_{K_{|\mathcal{T}_a|}} = 0 \quad \forall q_h \in \mathcal{P}_1(K_{|\mathcal{T}_a|}), $$

since $\nabla \cdot v_h \in \mathcal{P}_1(K_{|\mathcal{T}_a|})$. Using the Green theorem and the constraints in (A.13) as above, this is equivalent to verifying that

$$
\begin{align}
0 = (\boldsymbol{\mu}^a, \nabla \psi_F^k)_{K_{|\mathcal{T}_a|}} & - \langle \boldsymbol{w}_h^{j_1} \cdot \boldsymbol{n}_{K_{|\mathcal{T}_a|}}, \psi_F^k \rangle_{F^{j_1}} \tag{A.21} \\
& - \langle \boldsymbol{w}_h^{j_2} \cdot \boldsymbol{n}_{K_{|\mathcal{T}_a|}}, \psi_F^k \rangle_{F^{j_2}} - \langle \boldsymbol{w}_h^{j_3} \cdot \boldsymbol{n}_{K_{|\mathcal{T}_a|}}, \psi_F^k \rangle_{F^{j_3}} \nonumber
\end{align}
$$

for all $1 \le k \le 4$, where $F^{j_1}$, $F^{j_2}$, $F^{j_3}$ are the three faces in $\mathcal{F}_{|\mathcal{T}_a|}^\sharp$ and $\psi_F^k$ are the hat basis functions associated with the four vertices of $K_{|\mathcal{T}_a|}$. As in (A.19), assumptions of Lemma A.5 imply

$$ 0 = (\boldsymbol{\mu}^a, \nabla \psi_F^k)_{\omega_a} \quad 1 \le k \le 4. \tag{A.22} $$

Moreover, as in (A.18), it follows from (A.14), (A.15), and (A.16)

$$ \langle \boldsymbol{w}_h^j \cdot \boldsymbol{n}_{K_j}, \psi_F^k \rangle_{\partial K_j} = (\boldsymbol{\mu}^a, \nabla \psi_F^k)_{K_j} \quad 1 \le k \le 4 \tag{A.23} $$

is satisfied on all elements $K_j$ of the patch $\mathcal{T}_a$ other than $K_{|\mathcal{T}_a|}$. Let $\tilde{\omega}_a$ correspond to the patch subdomain $\omega_a$ without the element $K_{|\mathcal{T}_a|}$. Then, as in (A.20),

$$
\begin{align}
& - \langle \boldsymbol{w}_h^{j_1} \cdot \boldsymbol{n}_{K_{|\mathcal{T}_a|}}, \psi_F^k \rangle_{F^{j_1}} - \langle \boldsymbol{w}_h^{j_2} \cdot \boldsymbol{n}_{K_{|\mathcal{T}_a|}}, \psi_F^k \rangle_{F^{j_2}} - \langle \boldsymbol{w}_h^{j_3} \cdot \boldsymbol{n}_{K_{|\mathcal{T}_a|}}, \psi_F^k \rangle_{F^{j_3}} \nonumber \\
& = \langle \boldsymbol{w}_h^{j_1} \cdot \boldsymbol{n}_{\tilde{\omega}_a}, \psi_F^k \rangle_{F^{j_1}} + \langle \boldsymbol{w}_h^{j_2} \cdot \boldsymbol{n}_{\tilde{\omega}_a}, \psi_F^k \rangle_{F^{j_2}} + \langle \boldsymbol{w}_h^{j_3} \cdot \boldsymbol{n}_{\tilde{\omega}_a}, \psi_F^k \rangle_{F^{j_3}} \nonumber \\
& = \langle \boldsymbol{w}_h \cdot \boldsymbol{n}_{\tilde{\omega}_a}, \psi_F^k \rangle_{\partial \tilde{\omega}_a} = (\boldsymbol{w}_h, \nabla \psi_F^k)_{\tilde{\omega}_a} + (\nabla \cdot \boldsymbol{w}_h, \psi_F^k)_{\tilde{\omega}_a} \nonumber \\
& = \sum_{j; K_j^o \subset \tilde{\omega}_a} \{ (\boldsymbol{w}_h^j, \nabla \psi_F^k)_{K_j} + (\nabla \cdot \boldsymbol{w}_h^j, \psi_F^k)_{K_j} \} = \sum_{j; K_j^o \subset \tilde{\omega}_a} \langle \boldsymbol{w}_h^j \cdot \boldsymbol{n}_{K_j}, \psi_F^k \rangle_{\partial K_j} \tag{A.23} \\
& \stackrel{(A.22)}{=} \sum_{j; K_j^o \subset \tilde{\omega}_a} (\boldsymbol{\mu}^a, \nabla \psi_F^k)_{K_j} = (\boldsymbol{\mu}^a, \nabla \psi_F^k)_{\omega_a} - (\boldsymbol{\mu}^a, \nabla \psi_F^k)_{K_{|\mathcal{T}_a|}} \stackrel{(A.22)}{=} -(\boldsymbol{\mu}^a, \nabla \psi_F^k)_{K_{|\mathcal{T}_a|}} \nonumber
\end{align}
$$

for all $1 \le k \le 4$, i.e., (A.21). Thus, there exists a minimizer $w_h^{|\mathcal{T}_a|}$ of (A.13) on $K_{|\mathcal{T}_a|}$.

Step 2. Boundary patches with flat boundaries. We now investigate the case where the vertex $a \in V_h$ lies on the boundary $\partial \Omega$. We present in this step in detail the case of a boundary patch $\mathcal{T}_a$ for which $\Gamma_a$, the part of $\partial \omega_a$ that contains the faces sharing the vertex $a$, is contained in a plane $H$, which we call a “flat boundary” case. For the sake of simplicity, assume that either $\Gamma_a \subset \Gamma_D$ or $\Gamma_a \subset \Gamma_N$ and, without loss of generality, that $H = \{x \in \mathbb{R}^3; x_3 = 0\}$. The symmetrization operator around the plane $H$, $\phi: x = (x_1, x_2, x_3) \to (x_1, x_2, -x_3)$, as in [21, Section 7] and [11, Section 7], will be instrumental in the proof. Specifically, we introduce the symmetrized patch $\tilde{\mathcal{T}}_a := \mathcal{T}_a \cup \phi(\mathcal{T}_a)$, with the associated domain $\tilde{\omega}_a$, obtained by mapping the elements of $\mathcal{T}_a$ by $\phi$. We will employ the Piola mappings from (A.7) to relate the set $W_h(\mathcal{T}_a, \boldsymbol{\mu}^a)$ from the announcement of Lemma A.5 to a set $W_h(\tilde{\mathcal{T}}_a, \tilde{\boldsymbol{\mu}}^a)$ with an extended datum $\tilde{\boldsymbol{\mu}}^a$. Then, the result will follow by Step 1, since $\tilde{\mathcal{T}}_a$ is an interior patch.

Step 2a. The case $a \in \Gamma_D$. We start by defining the extended datum $\tilde{\boldsymbol{\mu}}^a \in L^2(\tilde{\omega}_a)$ from $\boldsymbol{\mu}^a$: we simply set $\tilde{\boldsymbol{\mu}}^a := \boldsymbol{\mu}^a$ in $\omega_a$ and $\tilde{\boldsymbol{\mu}}^a := \phi^d(\boldsymbol{\mu}^a)$ on $\tilde{\omega}_a \setminus \omega_a$. Let $\tilde{q}_h \in \mathcal{P}_1(\tilde{\mathcal{T}}_a) \cap H_*^1(\tilde{\omega}_a)$. Recalling that
$\phi$ is a symmetrization operator, its (constant) Jacobian matrix has a negative determinant. As a result, using (A.10) and (A.9),

$$ (\tilde{\mu}^a, \nabla \tilde{q}_h)_{\tilde{\omega}_a} = (\mu^a, \nabla \tilde{q}_h)_{\omega_a} + (\phi^d(\mu^a), \nabla \tilde{q}_h)_{\tilde{\omega}_a \setminus \omega_a} = (\mu^a, \nabla \tilde{q}_h)_{\omega_a} - (\mu^a, (\phi^c)^{-1}(\nabla \tilde{q}_h))_{\omega_a} \\ = (\mu^a, \nabla \tilde{q}_h)_{\omega_a} - (\mu^a, \nabla ((\phi^g)^{-1}(\tilde{q}_h)))_{\omega_a} = (\mu^a, \nabla q_h)_{\omega_a} $$

with $q_h := \tilde{q}_h - (\phi^g)^{-1}(\tilde{q}_h)$. Because $(\phi^g)^{-1}$ preserves the trace on $H$, we see that $q_h = 0$ on $H$ (see (A.8)), and since it also maps piecewise polynomials to piecewise polynomials (see (A.11)), $q_h \in P_1(\mathcal{T}_a) \cap H_*^1(\omega_a)$ (recall from Section 2.7 that $H_*^1(\omega_a) = \{v \in H^1(\omega_a); v = 0 \text{ on } \gamma_D = \Gamma_a\}$ here). Hence, $(\tilde{\mu}^a, \nabla \tilde{q}_h)_{\tilde{\omega}_a} = 0$ by our assumption $(\mu^a, \nabla q_h)_{\omega_a} = 0$. Thus $\tilde{\mu}^a$ satisfies the assumption of Lemma A.5 on the interior patch $\tilde{\mathcal{T}}_a$, and therefore Step 1 ensures that the set $W_h(\tilde{\mathcal{T}}_a, \tilde{\mu}^a)$ is non-empty.

We now consider an arbitrary element $\tilde{w}_h \in W_h(\tilde{\mathcal{T}}_a, \tilde{\mu}^a)$ and set $w_h := \tilde{w}_h|_{\omega_a}$. Since $\tilde{w}_h \in \mathcal{RT}_1(\tilde{\mathcal{T}}_a) \cap H_0(\text{div}, \tilde{\omega}_a)$, it is clear that $w_h \in \mathcal{RT}_1(\mathcal{T}_a) \cap H_0(\text{div}, \omega_a)$, namely as no normal trace boundary conditions are required on $\Gamma_a \subset H$. Indeed, in this case, $\Gamma_a = \gamma_D$ in the notation of Section 2.7, so that $H_0(\text{div}, \omega_a) = \{v \in H(\text{div}, \omega_a); v \cdot n_{\omega_a} = 0 \text{ on } \partial\omega_a \setminus \Gamma_a\}$. Moreover, $\nabla \cdot w_h = \nabla \cdot \tilde{w}_h = 0$ on $\omega_a$. Finally, $(w_h, r_h)_K = (\mu^a, r_h)_K$ for all $r_h \in [P_0(K)]^3$ and all $K \in \mathcal{T}_a$ since $\mathcal{T}_a \subset \tilde{\mathcal{T}}_a$ and simply $(w_h, r_h)_K = (\tilde{w}_h, r_h)_K$, $(\mu^a, r_h)_K = (\tilde{\mu}^a, r_h)_K$, and $\tilde{w}_h \in W_h(\tilde{\mathcal{T}}_a, \tilde{\mu}^a)$, so that $(\tilde{w}_h, r_h)_K = (\tilde{\mu}^a, r_h)_K$. This concludes the proof that the set $W_h(\mathcal{T}_a, \mu^a)$ is non-empty in this case.

Step 2b. The case $a \in \Gamma_N$. In this case, we extend the datum $\mu^a$ by setting $\tilde{\mu}^a := \mu^a$ on $\omega_a$ and $\tilde{\mu}^a := 0$ on $\tilde{\omega}_a \setminus \omega_a$. If $\tilde{q}_h \in P_1(\tilde{\mathcal{T}}_a) \cap H_*^1(\tilde{\omega}_a)$, we have

$$ (\tilde{\mu}^a, \nabla \tilde{q}_h)_{\tilde{\omega}_a} = (\mu^a, \nabla \tilde{q}_h)_{\omega_a} = 0 $$

since $\tilde{q}_h|_{\omega_a} \in P_1(\mathcal{T}_a) \cap H^1(\omega_a)$, whose gradients have the same span as those of $P_1(\mathcal{T}_a) \cap H_*^1(\omega_a)$, the zero mean value subspace of $P_1(\mathcal{T}_a) \cap H^1(\omega_a)$ following Section 2.7 in this case. It thus follows from Step 1 that $W_h(\tilde{\mathcal{T}}_a, \tilde{\mu}^a)$ is non-empty.

Consider an element $\tilde{w}_h \in W_h(\tilde{\mathcal{T}}_a, \tilde{\mu}^a)$ and set $w_h := \tilde{w}_h|_{\omega_a} - (\phi^d)^{-1}(\tilde{w}_h|_{\tilde{\omega}_a \setminus \omega_a})$. We need to show that $w_h \in W_h(\mathcal{T}_a, \mu^a)$. Recall that here the functions in $H_0(\text{div}, \omega_a)$ need to satisfy the no-flow boundary condition on the whole patch boundary $\partial\omega_a$ and in particular on $H: H_0(\text{div}, \omega_a) = \{v \in H(\text{div}, \omega_a); v \cdot n_{\omega_a} = 0 \text{ on } \partial\omega_a\}$ from Section 2.7 in this case. Since the Piola mapping $(\phi^d)^{-1}$ maps piecewise Raviart-Thomas polynomials to piecewise Raviart-Thomas polynomials (cf. (A.11)) and preserves the divergence (cf. (A.9)) and the normal trace (cf. (A.8)), it is clear that $w_h \in \mathcal{RT}_1(\mathcal{T}_a) \cap H_0(\text{div}, \omega_a)$ and $\nabla \cdot w_h = 0$. It remains to show that $(w_h, r_h)_K = (\mu^a, r_h)_K$ for all $r_h \in [P_0(K)]^3$ and all $K \in \mathcal{T}_a$. Let $K \in \mathcal{T}_a$ and $r_h \in [P_0(K)]^3$ and let $\tilde{K}$ be the tetrahedron corresponding to $K$ by the symmetry map $\phi$. Then

$$ (w_h, r_h)_K = (\tilde{w}_h, r_h)_K - ((\phi^d)^{-1}(\tilde{w}_h), r_h)_K = (\tilde{w}_h, r_h)_K + (\tilde{w}_h, \phi^c(r_h))_{\tilde{K}} \\ = (\tilde{\mu}^a, r_h)_K + (\tilde{\mu}^a, \phi^c(r_h))_{\tilde{K}} = (\tilde{\mu}^a, r_h)_K = (\mu^a, r_h)_K, $$

where we have used (A.10), that the Piola mapping $\phi^c$ maps piecewise constant vectors onto piecewise constant vectors (this can be seen from the definition (A.7) of $\phi^c$ since its Jacobian matrix is constant here), that $\tilde{w}_h \in W_h(\tilde{\mathcal{T}}_a, \tilde{\mu}^a)$, and finally that $\tilde{\mu}^a$ is the extension of $\mu^a$ by zero to the symmetrized patch. This concludes the proof that $W_h(\mathcal{T}_a, \mu^a)$ is non-empty in this case.

Step 3: general boundary patches. For general boundary patches, the proof follows the lines of Step 2 while employing the extension and restriction operators introduced in [11, Section 7] instead of the (simpler) symmetrization operator $\phi$ of Step 2. We do not give details here. $\square$

## A.6 Proof of Lemma A.4

We can now finally establish a proof of Lemma A.4.

*Proof of Lemma A.4. Step 1. Existence and uniqueness.* The minimization set in (A.5) is the set $W_h(\mathcal{T}_a, \mu^a)$ of Lemma A.5 with $\mu^a := \tau_h^a - \theta_h^a$. Since the minimization functional in (A.5) is convex, it is sufficient to show that $W_h(\mathcal{T}_a, \mu^a)$ is non-empty to ensure the existence and uniqueness of $\epsilon_h^a$ from (A.5). From Lemma A.5, we need to show that $(\mu^a, \nabla q_h)_{\omega_a} = 0$ for all $q_h \in P_1(\mathcal{T}_a) \cap H_*^1(\omega_a)$. This is actually a direct consequence of assumption (A.1c). Indeed, from the divergence constraint in (A.3) and since $q' \ge 1$ and $q_h \in P_1(\mathcal{T}_a) \cap H_*^1(\omega_a)$, we have

$$ (\mu^a, \nabla q_h)_{\omega_a} = (\tau_h^a, \nabla q_h)_{\omega_a} - (\theta_h^a, \nabla q_h)_{\omega_a} = (\tau_h^a, \nabla q_h)_{\omega_a} + (\nabla \cdot \theta_h^a, q_h)_{\omega_a} \\ = (\tau_h^a, \nabla q_h)_{\omega_a} + (\Pi_{q'}(g^a), q_h)_{\omega_a} = (\tau_h^a, \nabla q_h)_{\omega_a} + (g^a, q_h)_{\omega_a} \stackrel{(A.1c)}{=} 0. $$
**Step 2. Stability bound.** We now proceed with the proof of the stability (A.6).

**Step 2a: generic stability bound.** Set again $\mu^a := \tau_h^a - \bar{\theta}_h^a$, and denote by $\mu_h^a$ the $L^2(\omega_a)$-orthogonal projection of $\mu^a$ onto $[P_1(\mathcal{T}_a)]^3$. Considering the Euler(-Lagrange) equations associated with (A.5), it is clear that we can equivalently replace $\mu^a$ by $\mu_h^a$ in the definition (A.5) of $\bar{\epsilon}_h^a$. Furthermore, because (A.5) is a quadratic minimization problem with linear constraints, the operator $T: [P_1(\mathcal{T}_a)]^3 \ni \mu_h^a \to \bar{\epsilon}_h^a \in \mathcal{RT}_1(\mathcal{T}_a) \cap H_0(\text{div}, \omega_a)$ (well-defined from Step 1) is linear. Since both $[P_1(\mathcal{T}_a)]^3$ and $\mathcal{RT}_1(\mathcal{T}_a) \cap H_0(\text{div}, \omega_a)$ are finite-dimensional spaces, the operator $T$ is continuous, and there exists a constant $C(\mathcal{T}_a)$ such that

$$ \|\bar{\epsilon}_h^a\|_{\omega_a} \le C(\mathcal{T}_a) \|\mu_h^a\|_{\omega_a} \le C(\mathcal{T}_a) \|\mu^a\|_{\omega_a}, \quad (A.24) $$

where we used the fact that $\mu_h^a$ is defined from $\mu^a$ by projection in the last inequality. The constant $C(\mathcal{T}_a)$ is independent of the polynomial degree $q$ (recall that (A.5) works with $\mathcal{RT}_1$ elements only) but depends on the patch $\mathcal{T}_a$ in an unspecified way. To make the dependence explicit, we resort in the next step to a reference patch.

**Step 2b: explicit stability bound.** For a fixed shape-regularity parameter $\kappa_{\mathcal{T}_h}$ from Section 2.4, there exists a maximal number of elements $N(\kappa_{\mathcal{T}_h})$ allowed in any patch $\mathcal{T}_a$. In turn, for any $N(\kappa_{\mathcal{T}_h})$, there exists a finite set of reference patches $\{\hat{T}\}$ such that for all vertex patches $\mathcal{T}_a$, there exists a reference patch $\hat{T}$ and a bilipschitz mapping $\phi: \omega_a \to \hat{\omega}$ ($\hat{\omega}$ being the open domain associated with $\hat{T}$) such that $\phi|_K$ is an affine mapping between the tetrahedron $K \in \mathcal{T}_a$ and a tetrahedron $\hat{K} \in \hat{T}$. The associated Piola mapping $\phi^d$ from (A.7) will be useful.

Crucially, we observe that for all $\hat{K} \in \hat{T}$, $v \in L^2(K)$, and $\hat{r}_h \in [P_0(\hat{K})]^3$, there exists $r_h \in [P_0(K)]^3$ such that $(\phi^d(v), \hat{r}_h)_{\hat{K}} = (v, r_h)_K$, since, elementwise, the Piola transform amounts to a multiplication by a constant matrix and a change of coordinates. It follows that $\phi^d$ maps the minimization set of (A.5) on $\mathcal{T}_a$ into the minimization set of the equivalent problem set on $\hat{T}$ with constraints $\phi^d(\mu^a)$.

Now, on the reference patch $\hat{T}$, if $\hat{\epsilon}_h$ is the minimizer of (A.5) with the datum $\phi^d(\mu^a)$, we conclude from Step 2a that

$$ \|\hat{\epsilon}_h\|_{\hat{\omega}} \le C(\kappa_{\mathcal{T}_h}) \|\phi^d(\mu^a)\|_{\hat{\omega}} \le C(\kappa_{\mathcal{T}_h}) \|\phi^d\|_{\omega_a} \|\mu^a\|_{\omega_a}. $$

On the other hand, since $(\phi^d)^{-1}(\hat{\epsilon}_h)$ belongs the minimization set on $\mathcal{T}_a$, we have

$$ \|\bar{\epsilon}_h^a - \mu^a\|_{\omega_a} \le \|(\phi^d)^{-1}(\hat{\epsilon}_h) - \mu^a\|_{\omega_a} \le \|(\phi^d)^{-1}\|\|\bar{\epsilon}_h\|_{\hat{\omega}} + \|\mu^a\|_{\omega_a}, $$

so that

$$ \|\bar{\epsilon}_h^a - \mu^a\|_{\omega_a} \le (1 + C(\kappa_{\mathcal{T}_h}) \|(\phi^d)^{-1}\|\|\phi^d\|) \|\mu^a\|_{\omega_a}. $$

At this point, we conclude the proof since $\|(\phi^d)^{-1}\|\|\phi^d\|$ only depends on $\kappa_{\mathcal{T}_h}$ due to (A.12) and $\|\bar{\epsilon}_h^a\|_{\omega_a} \le \|\mu^a\|_{\omega_a} + \|\bar{\epsilon}_h^a - \mu^a\|_{\omega_a}$. $\square$

## B Decomposition of a divergence-free piecewise polynomial with an elementwise orthogonality into local divergence-free contributions

Let $q \ge 0$ be a fixed integer and recall the notation of Section 2; namely, $I_{K,q}^{\mathcal{RT}}$ is the canonical $q$-degree Raviart-Thomas interpolate on the given mesh element $K \in \mathcal{T}_h$ from (2.4), and $\lesssim$ means smaller or equal to up to a constant only depending on the mesh shape-regularity parameter $\kappa_{\mathcal{T}_h}$. The following result is of independent interest:

**Theorem B.1** (Decomposition of a divergence-free Raviart-Thomas piecewise polynomial with an elementwise orthogonality constraint into local divergence-free contributions). Let

$$ \delta_h \in \mathcal{RT}_q(\mathcal{T}_h) \cap H_{0,N}(\text{div}, \Omega) \quad \text{with} \quad \nabla \cdot \delta_h = 0 \quad (B.1) $$

be a divergence-free $q$-degree Raviart-Thomas piecewise polynomial that is elementwise orthogonal to vector-valued constants,

$$ (\delta_h, r_h)_K = 0 \quad \forall r_h \in [P_0(K)]^3, \forall K \in \mathcal{T}_h. \quad (B.2) $$
Then there exists a unique solution to the $q'$-degree Raviart-Thomas elementwise minimizations, $q' = q$ or $q' = q + 1$,

$$ \delta_h^a := \arg \min_{\substack{v_h \in \mathcal{RT}_{q'}(K) \\ \nabla \cdot v_h = 0 \\ v_h \cdot n_K = \mathcal{IR}_{K,q'}((\psi^a \delta_h)|_K) \cdot n_K \text{ on } \partial K}} \|v_h - \mathcal{IR}_{K,q'}((\psi^a \delta_h)|_K)\|_K^2 \quad (B.3) $$

for all tetrahedra $K \in \mathcal{T}_h$ and all vertices $a \in \mathcal{V}_K$. This yields patchwise divergence-free contributions

$$ \delta_h^a \in \mathcal{RT}_{q'}(\mathcal{T}_a) \cap \mathcal{H}_0(\mathrm{div}, \omega_a) \quad \text{with} \quad \nabla \cdot \delta_h^a = 0 \quad \forall a \in \mathcal{V}_h, \quad (B.4) $$

decomposing $\delta_h$ as

$$ \delta_h = \sum_{a \in \mathcal{V}_h} \delta_h^a. \quad (B.5) $$

Moreover, for all tetrahedra $K \in \mathcal{T}_h$ and all vertices $a \in \mathcal{V}_K$, there hold the local stability estimates

$$ \|\delta_h^a - \mathcal{IR}_{K,q'}((\psi^a \delta_h)|_K)\|_K \lesssim \|\delta_h\|_K, \quad (B.6a) $$

$$ \|\delta_h^a\|_K \lesssim_{q'} \|\delta_h\|_K, \quad (B.6b) $$

where $\lesssim_{q'}$ means $\lesssim$ for $q' = q + 1$ and up to a constant only depending on the mesh shape-regularity parameter $\kappa_{\mathcal{T}_h}$ and the degree $q$ when $q' = q$.

**Remark B.2** (The two settings $q' = q$ or $q' = q + 1$ in Theorem B.1). With the choice $q' = q$, the contributions $\delta_h^a$ in Theorem B.1 stay in the same degree Raviart-Thomas space as the datum $\delta_h$, but, unfortunately, the stability (B.6b) is not necessarily q-robust. For q-robustness, the choice $q' = q + 1$, increasing the degree of $\delta_h^a$ by one, is to be used. Note that in this case, the Raviart-Thomas interpolator $\mathcal{IR}_{K,q'}^{\mathcal{RT}}$ can be disregarded, since then $\mathcal{IR}_{K,q'}^{\mathcal{RT}}((\psi^a \delta_h)|_K) = (\psi^a \delta_h)|_K$.

*Proof.* Let $\delta_h$ satisfy (B.1) and (B.2). We address (B.3)-(B.6) in four steps.

**Step 1.** *Proof of the well-posedness of* (B.3). Fix $K \in \mathcal{T}_h$ and $a \in \mathcal{V}_K$. The existence and uniqueness of $\delta_h^a|_K$ from (B.3) are classical following, e.g., [6], when the Neumann compatibility condition $\langle \mathcal{IR}_{K,q'}^{\mathcal{RT}}((\psi^a \delta_h)|_K) \cdot n_K, 1 \rangle_{\partial K} = 0$ is satisfied. This can be shown via (2.4a), the Green theorem, the assumption $\nabla \cdot \delta_h = 0$ in (B.1), and the elementwise orthogonality assumption (B.2) (note that $(\nabla \psi^a)|_K \in [\mathcal{P}_0(K)]^3$) as

$$ \langle \mathcal{IR}_{K,q'}^{\mathcal{RT}}((\psi^a \delta_h)|_K) \cdot n_K, 1 \rangle_{\partial K} = \langle \psi^a \delta_h \cdot n_K, 1 \rangle_{\partial K} = \langle \delta_h \cdot n_K, \psi^a \rangle_{\partial K} = (\nabla \cdot \delta_h, \psi^a)_K + (\delta_h, \nabla \psi^a)_K = 0. $$

**Step 2.** *Proof of the stability estimates* (B.6). Still for a fixed $K \in \mathcal{T}_h$ and $a \in \mathcal{V}_K$, consider the problem

$$ \hat{\delta}_h^a := \arg \min_{\substack{v_h \in \mathcal{RT}_{q'}(K) \\ \nabla \cdot v_h = (-\nabla \psi^a \cdot \delta_h)|_K \\ v_h \cdot n_K = 0 \text{ on } \partial K}} \|v_h\|_K^2. \quad (B.7) $$

This problem is again well-posed since, from (B.2), $(\nabla \psi^a \cdot \delta_h, 1)_K = (\delta_h, \nabla \psi^a)_K = 0$; moreover, $(\nabla \psi^a \cdot \delta_h)|_K \in \mathcal{P}_q(K) \subset \mathcal{P}_{q'}(K)$, since from $\nabla \cdot \delta_h = 0$, it follows that $\delta_h|_K \in [\mathcal{P}_q(K)]^3$ (see, e.g., [6, Corollary 2.3.1]). It follows that $\hat{\delta}_h^a|_K = \delta_h^a|_K - \mathcal{IR}_{K,q'}^{\mathcal{RT}}((\psi^a \delta_h)|_K)$; indeed, the commuting property (2.5) yields, on the simplex $K$, $\nabla \cdot (\mathcal{IR}_{K,q'}^{\mathcal{RT}}(\psi^a \delta_h)) = \mathcal{P}_{q'}(\nabla \cdot (\psi^a \delta_h)) = \mathcal{P}_{q'}(\nabla \psi^a \cdot \delta_h) = \nabla \psi^a \cdot \delta_h$. Problem (B.7) fits the framework of [21, Lemma A.3] with $r_F = 0$, $r_K = (-\nabla \psi^a \cdot \delta_h)|_K$, and $p = q'$, so that

$$ \|\delta_h^a - \mathcal{IR}_{K,q'}^{\mathcal{RT}}((\psi^a \delta_h)|_K)\|_K = \|\hat{\delta}_h^a\|_K = \min_{\substack{v_h \in \mathcal{RT}_{q'}(K) \\ \nabla \cdot v_h = -\nabla \psi^a \cdot \delta_h \\ v_h \cdot n_K = 0 \text{ on } \partial K}} \|v_h\|_K \lesssim \min_{v \in \mathcal{H}(\mathrm{div}, K)} \|v\|_K = \|\nabla \zeta_K\|_K. $$

Here, by primal-dual equivalence, $\zeta_K \in H_*^1(K)$ is such that

$$ (\nabla \zeta_K, \nabla v)_K = -(\nabla \psi^a \cdot \delta_h, v)_K \quad \forall v \in H_*^1(K) $$

with $H_*^1(K) := \{v \in H^1(K); (v, 1)_K = 0\}$, where the Poincaré inequality gives $\|v\|_K \lesssim h_K \|\nabla v\|_K$. Then the Cauchy-Schwarz inequality and shape regularity yield

$$ \|\nabla \zeta_K\|_K = \max_{\substack{v \in H_*^1(K) \\ \|\nabla v\|_K = 1}} (\nabla \zeta_K, \nabla v)_K = \max_{\substack{v \in H_*^1(K) \\ \|\nabla v\|_K = 1}} -(\nabla \psi^a \cdot \delta_h, v)_K \lesssim \|\nabla \psi^a\|_{\infty,K} \|\delta_h\|_K h_K \lesssim \|\delta_h\|_K. $$
Combining the two above estimates gives the desired stability result (B.6a). The other stability result (B.6b) follows from (B.6a) by the triangle inequality together with the non-q-robust stability bound $||\mathbf{I}_{K,q'}^{RT}((\psi^a \delta_h)|_K)||_K \lesssim q'$, $||\psi^a \delta_h||_K \le ||\delta_h||_K$ when $q' = q$, whereas $||\mathbf{I}_{K,q'}^{RT}((\psi^a \delta_h)|_K)||_K = ||\psi^a \delta_h||_K \le ||\delta_h||_K$ when $q' = q + 1$.

**Step 3.** *Proof of the patchwise properties (B.4).* The first property in (B.4) follows from the prescription of the normal components in (B.3), whereas the second one is the divergence prescription in (B.3).

**Step 4.** *Proof of the decomposition (B.5).* Finally, in order to prove (B.5), set $\tilde{\delta}_h := \sum_{a \in \mathcal{V}_h} \delta_h^a$. Now fix an element $K \in \mathcal{T}_h$ and remark that from the normal trace constraint in (B.3) and the linearity of the interpolator $\mathbf{I}_{K,q'}^{RT}$, on $\partial K$,

$$ \begin{aligned} \tilde{\delta}_h|_K \cdot \mathbf{n}_K &= \sum_{a \in \mathcal{V}_K} \delta_h^a|_K \cdot \mathbf{n}_K = \sum_{a \in \mathcal{V}_K} \mathbf{I}_{K,q'}^{RT}((\psi^a \delta_h)|_K) \cdot \mathbf{n}_K = \mathbf{I}_{K,q'}^{RT} \left[ \sum_{a \in \mathcal{V}_K} (\psi^a \delta_h)|_K \right] \cdot \mathbf{n}_K \\ &= \mathbf{I}_{K,q'}^{RT}(\delta_h)|_K \cdot \mathbf{n}_K = \delta_h|_K \cdot \mathbf{n}_K \end{aligned} $$

also using the partition of unity (2.1). Similarly, by the divergence constraint in (B.3) and $\nabla \cdot \delta_h = 0$ from (B.1), on $K$,

$$ \nabla \cdot \tilde{\delta}_h = \sum_{a \in \mathcal{V}_K} \nabla \cdot \delta_h^a = 0 = \nabla \cdot \delta_h. $$

Consequently, $(\tilde{\delta}_h - \delta_h)|_K \in \mathcal{RT}_{q'}(K)$ has zero normal trace and divergence. Moreover, the Euler conditions of problem (B.3) state

$$ (\delta_h^a - \mathbf{I}_{K,q'}^{RT}((\psi^a \delta_h)|_K), v_h)_K = 0 \quad \forall v_h \in \mathcal{RT}_{q'}(K) \text{ with } \nabla \cdot v_h = 0 \text{ and } v_h \cdot \mathbf{n}_K = 0 $$

on $\partial K$. Summing this over all vertices $a \in \mathcal{V}_K$ and using again the linearity of $\mathbf{I}_{K,q'}^{RT}$,

$$ (\tilde{\delta}_h - \delta_h, v_h)_K = 0 \quad \forall v_h \in \mathcal{RT}_{q'}(K) \text{ with } \nabla \cdot v_h = 0 \text{ and } v_h \cdot \mathbf{n}_K = 0 \text{ on } \partial K, $$

so that indeed $\tilde{\delta}_h = \delta_h$ on any mesh element $K \in \mathcal{T}_h$. $\square$

## References

[1] Adams, R. A. Pure and Applied Mathematics, Vol. 65. Sobolev spaces. Academic Press [A subsidiary of Harcourt Brace Jovanovich, Publishers], New York-London, 1975.

[2] Ainsworth, M., and Oden, J. T. A posteriori error estimation in finite element analysis. Pure and Applied Mathematics (New York). Wiley-Interscience [John Wiley & Sons], New York, 2000.

[3] Alonso Rodríguez, A., Camaño, J., De Los Santos, E., and Rapetti, F. A graph approach for the construction of high order divergence-free Raviart-Thomas finite elements. Calcolo 55 (2018), Paper No. 42, 28. https://doi.org/10.1007/s10092-018-0284-1.

[4] Alonso Rodríguez, A., Camaño, J., Ghiloni, R., and Valli, A. Graphs, spanning trees and divergence-free finite elements in domains of general topology. IMA J. Numer. Anal. 37 (2017), 1986–2003. https://doi.org/10.1093/imanum/drw047.

[5] Apel, T. Anisotropic finite elements: local estimates and applications. Advances in Numerical Mathematics. B. G. Teubner, Stuttgart, 1999.

[6] Boffi, D., Brezzi, F., and Fortin, M. Mixed finite element methods and applications, vol. 44 of Springer Series in Computational Mathematics. Springer, Heidelberg, 2013. https://doi.org/10.1007/978-3-642-36519-5.

[7] Braess, D., Pillwein, V., and Schöberl, J. Equilibrated residual error estimates are p-robust. Comput. Methods Appl. Mech. Engrg. 198 (2009), 1189–1197. http://dx.doi.org/10.1016/j.cma.2008.12.010.

[8] Braess, D., and Schöberl, J. Equilibrated residual error estimator for edge elements. Math. Comp. 77 (2008), 651–672. http://dx.doi.org/10.1090/S0025-5718-07-02080-7.

[9] Chaumont-Frelet, T., Ern, A., and Vohralík, M. Polynomial-degree-robust H(curl)-stability of discrete minimization in a tetrahedron. C. R. Math. Acad. Sci. Paris 358 (2020), 1101–1110. https://doi.org/10.5802/crmath.133.
[10] Chaumont-Frelet, T., Ern, A., and Vohralík, M. Stable broken $H(\text{curl})$ polynomial extensions and $p$-robust a posteriori error estimates by broken patchwise equilibration for the curl-curl problem. *Math. Comp.* **91** (2022), 37–74. https://doi.org/10.1090/mcom/3673.

[11] Chaumont-Frelet, T., and Vohralík, M. Constrained and unconstrained stable discrete minimizations for $p$-robust local reconstructions in vertex patches in the de Rham complex. HAL Preprint 03749682, submitted for publication, https://hal.inria.fr/hal-03749682, 2022.

[12] Costabel, M., Dauge, M., and Nicaise, S. Singularities of Maxwell interface problems. *M2AN Math. Model. Numer. Anal.* **33** (1999), 627–649. https://doi.org/10.1051/m2an:1999155.

[13] Costabel, M., and McIntosh, A. On Bogovskiǐ and regularized Poincaré integral operators for de Rham complexes on Lipschitz domains. *Math. Z.* **265** (2010), 297–320. http://dx.doi.org/10.1007/s00209-009-0517-8.

[14] Demkowicz, L., Gopalakrishnan, J., and Schöberl, J. Polynomial extension operators. Part II. *SIAM J. Numer. Anal.* **47** (2009), 3293–3324. http://dx.doi.org/10.1137/070698798.

[15] Destuynder, P., and Métivet, B. Explicit error bounds in a conforming finite element method. *Math. Comp.* **68** (1999), 1379–1396. http://dx.doi.org/10.1090/S0025-5718-99-01093-5.

[16] Dobrzynski, C. MMG3D: User guide. Tech. Rep. 422, Inria, France, 2012.

[17] Dörfler, W. A convergent adaptive algorithm for Poisson’s equation. *SIAM J. Numer. Anal.* **33** (1996), 1106–1124. http://dx.doi.org/10.1137/0733054.

[18] Ern, A., and Guermond, J.-L. *Finite Elements I. Approximation and Interpolation*, vol. **72** of *Texts in Applied Mathematics*. Springer International Publishing, Springer Nature Switzerland AG, 2021. https://doi-org.ezproxy.is.cuni.cz/10.1007/978-3-030-56341-7.

[19] Ern, A., and Vohralík, M. Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs. *SIAM J. Sci. Comput.* **35** (2013), A1761–A1791. http://dx.doi.org/10.1137/120896918.

[20] Ern, A., and Vohralík, M. Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations. *SIAM J. Numer. Anal.* **53** (2015), 1058–1081. http://dx.doi.org/10.1137/130950100.

[21] Ern, A., and Vohralík, M. Stable broken $H^1$ and $H(\text{div})$ polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions. *Math. Comp.* **89** (2020), 551–594. http://dx.doi.org/10.1090/mcom/3482.

[22] Fernandes, P., and Gilardi, G. Magnetostatic and electrostatic problems in inhomogeneous anisotropic media with irregular boundary and mixed boundary conditions. *Math. Models Methods Appl. Sci.* **7** (1997), 957–991. https://doi.org/10.1142/S0218202597000487.

[23] Gedicke, J., Geevers, S., and Perugia, I. An equilibrated a posteriori error estimator for arbitrary-order Nédélec elements for magnetostatic problems. *J. Sci. Comput.* **83** (2020), Paper No. 58, 23. https://doi.org/10.1007/s10915-020-01224-x.

[24] Gedicke, J., Geevers, S., Perugia, I., and Schöberl, J. A polynomial-degree-robust a posteriori error estimator for Nédélec discretizations of magnetostatic problems. *SIAM J. Numer. Anal.* **59** (2021), 2237–2253. https://doi-org.ezproxy.is.cuni.cz/10.1137/20M1333365.

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

[26] Gross, P. W., and Kotiuga, P. R. *Electromagnetic theory and computation: a topological approach*, vol. **48** of *Mathematical Sciences Research Institute Publications*. Cambridge University Press, Cambridge, 2004. https://doi.org/10.1017/CBO9780511756337.

[27] Hannukainen, A., Stenberg, R., and Vohralík, M. A unified framework for a posteriori error estimation for the Stokes problem. *Numer. Math.* **122** (2012), 725–769. http://dx.doi.org/10.1007/s00211-012-0472-x.
[28] Hiptmair, R., and Pechstein, C. *Discrete regular decompositions of tetrahedral discrete 1-forms*. De Gruyter, 2019, ch. 7, pp. 199–258. https://doi.org/10.1515/9783110543612-007.

[29] Ladevèze, P., and Leguillon, D. Error estimate procedure in the finite element method and applications. *SIAM J. Numer. Anal.* **20** (1983), 485–509.

[30] Licht, M. W. Higher-order finite element de Rham complexes, partially localized flux reconstructions, and applications. Preprint, http://www.math.ucsd.edu/~mlicht/pdf/preprint.hoflux.pdf, 2019.

[31] Luce, R., and Wohlmuth, B. I. A local a posteriori error estimator based on equilibrated fluxes. *SIAM J. Numer. Anal.* **42** (2004), 1394–1414. http://dx.doi.org/10.1137/S0036142903433790.

[32] Monk, P. *Finite element methods for Maxwell's equations*. Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2003. https://doi.org/10.1093/acprof:oso/9780198508885.001.0001.

[33] Nédélec, J.-C. Mixed finite elements in $\mathbb{R}^3$. *Numer. Math.* **35** (1980), 315–341.

[34] Nicaise, S., Witowski, K., and Wohlmuth, B. I. An a posteriori error estimator for the Lamé equation based on equilibrated fluxes. *IMA J. Numer. Anal.* **28** (2008), 331–353. http://dx.doi.org/10.1093/imanum/drm008.

[35] Prager, W., and Synge, J. L. Approximations in elasticity based on the concept of function space. *Quart. Appl. Math.* **5** (1947), 241–269.

[36] Raviart, P.-A., and Thomas, J.-M. A mixed finite element method for 2nd order elliptic problems. In *Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975)*. Springer, Berlin, 1977, pp. 292–315. Lecture Notes in Math., Vol. 606.

[37] Scheichl, R. Decoupling three-dimensional mixed problems using divergence-free finite elements. *SIAM J. Sci. Comput.* **23** (2002), 1752–1776. http://dx.doi.org/10.1137/S1064827500375886.

[38] Ziegler, G. M. *Lectures on polytopes*, vol. **152** of *Graduate Texts in Mathematics*. Springer-Verlag, New York, 1995. http://dx.doi.org/10.1007/978-1-4613-8431-1.