*p*-robust equivalence of global continuous constrained
and local discontinuous unconstrained approximation, a
*p*-stable local commuting projector, and optimal
elementwise *hp* approximation estimates in *H*(div)

Leszek F. Demkowicz, Martin Vohralík

► To cite this version:

Leszek F. Demkowicz, Martin Vohralík. *p*-robust equivalence of global continuous constrained and local discontinuous unconstrained approximation, a *p*-stable local commuting projector, and optimal elementwise *hp* approximation estimates in *H*(div). 2024. hal-04503603v2

HAL Id: hal-04503603

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

Preprint submitted on 2 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
p-robust equivalence of global continuous constrained and local
discontinuous unconstrained approximation, a p-stable local
commuting projector, and optimal elementwise hp approximation
estimates in H(div)*

Leszek Demkowicz† Martin Vohralík‡§

July 2, 2025

Abstract

Let an open bounded Lipschitz polygon or polyhedron Ω, a function **v** in the Sobolev space *H*(div, Ω), and a simplicial mesh of Ω be given. We prove the equivalence of two piecewise (Raviart–Thomas) polynomial best approximations of **v** in the *L*²-norm: 1) globally on the whole computational domain Ω, with the normal trace continuity requirement and a divergence constraint; 2) locally on each mesh element, without any interelement continuity requirement and without any constraint on the divergence. The former (global-best continuous constrained piecewise polynomial approximation) arises in numerical methods for partial differential equations related to the *H*(div, Ω) space, whereas the latter (local-best discontinuous unconstrained piecewise polynomial approximation) is a key quantity in approximation theory. Crucially, we establish *p*-robustness in that the equivalence constant only depends on the mesh shape regularity and the spatial dimension. This improves the recent result of [IMA J. Numer. Anal. **42** (2022), 1023–1049], where the equivalence constant was possibly dependent on the underlying polynomial degree. Consequently, we obtain fully *h*- and *p*- (mesh-size- and polynomial-degree-) optimal approximation estimates under the minimal Sobolev regularity only requested separately on each mesh element. These two results immediately follow by our construction of an operator from the infinite-dimensional Sobolev space *H*(div, Ω) to its finite-dimensional Raviart–Thomas subspace that has the following properties: 1) it is defined over the entire *H*(div, Ω) and preserves boundary conditions imposed on a part of the boundary of Ω; 2) it is defined locally in a neighborhood of each mesh element; 3) it is based on elementwise *L*²-orthogonal polynomial projections; 4) it is a projector, i.e., it leaves intact objects that are already in the Raviart–Thomas piecewise polynomial space; 5) it is locally and *p*-robustly stable in the *L*²-norm, up to *hp* data oscillation; 6) its approximation property is locally and *p*-robustly equivalent to that of the discontinuous unconstrained (elementwise *L*²-orthogonal) projection; 7) it satisfies the commuting property with the *L*²-orthogonal projection onto piecewise polynomials.

**Key words:** Sobolev space *H*(div), best approximation, continuous approximation, discontinuous approximation, Raviart–Thomas space, local–global equivalence, constrained–unconstrained equivalence, minimal regularity, elementwise regularity, commuting projector, *hp* finite elements, error bound, polynomial-degree robustness

# 1 Introduction

For the space dimension $d = 2, 3$, let $\Omega \subset \mathbb{R}^d$ be an open, bounded, Lipschitz, polygonal or polyhedral domain. Let $H(\text{div}, \Omega)$ be the Sobolev space of functions square-integrable together with their weak divergences, cf. Girault and Raviart [38], Ern and Guermond [31], or Demkowicz [23]. Let a shape-regular simplicial mesh $\mathcal{T}_h$ of $\Omega$ and a polynomial degree $p \ge 0$ be fixed (details on the setting and notation are given in Section 2 below).

<footnote>*The second author acknowledges the J. Tinsley Oden Faculty Fellowship Research Program for his research stay in Austin, Texas.</footnote>

<footnote>†Oden Institute for Computational and Engineering Sciences, 1 University Station, C0200, The University of Texas at Austin, Texas 78712, U.S.A. (leszek@oden.utexas.edu).</footnote>

<footnote>‡Inria, 48 rue Barrault, 75647 Paris, France (martin.vohralik@inria.fr).</footnote>

<footnote>§CERMICS, Ecole nationale des ponts et chaussées, IP Paris, 77455 Marne la Vallée, France.</footnote>
## 1.1 Commuting projectors under minimal regularity

In analysis of numerical methods related to the $H(\text{div}, \Omega)$ space, a crucial role is played by the design of operators $P_{hp}^{\text{div}}$ and $\Pi_{hp}$ such that

$$ \begin{gathered} \mathbf{H}_{0,N}(\text{div}, \Omega) \xrightarrow{\nabla} \mathbf{L}^2(\Omega) \\ \mathbf{P}_{hp}^{\text{div}} \xrightarrow{\nabla} \mathcal{P}_p(\mathcal{T}_h) \cap \mathcal{L}_*^2(\Omega). \end{gathered} \tag{1.1} $$

Here, particularly, $P_{hp}^{\text{div}}$ needs to be defined over the entire infinite-dimensional space $H(\text{div}, \Omega)$, which excludes the so-called canonical Raviart-Thomas projector, cf. [8, 24, 46]. Moreover, $P_{hp}^{\text{div}}$ needs to be a projector, i.e., leave intact objects that are already in the Raviart-Thomas piecewise polynomial space $RT_p(\mathcal{T}_h) \cap H_{0,N}(\text{div}, \Omega)$, and commute in the form expressed in (1.1), which excludes Clément-type [19] (quasi-)interpolation. Moreover, $P_{hp}^{\text{div}}$ should be defined locally, in a neighborhood of each mesh element, and $\Pi_{hp}$ should be the $L^2$-orthogonal projection onto piecewise $p$-degree polynomials. A seminal contribution in this direction is that Falk and Winther [34], following Christiansen and Winther [18] (locality is not devised), and followed by Ern and Guermond [29, 30] (locality or commuting is not devised), Licht [41] (locality is not devised), Arnold and Guzmán [3] ($\Pi_{hp}$ is not the $L^2$-orthogonal projection onto piecewise $p$-degree polynomials), and Gawlik et al. [37] (commuting is not devised). As stated, (1.1) is achieved in Ern et al. [28, Theorem 3.2].

## 1.2 hp approximation estimates in $H(\text{div}, \Omega)$

In addition the properties discussed above, $P_{hp}^{\text{div}}$ from (1.1) should also have correct approximation properties, both with respect to the mesh size $h$ and the polynomial degree $p$. Here, $h$-approximation is customary but $p$-approximation is much more seldom, and more difficult. Up to logarithmic factors in $p$, the latter was achieved in particular in Demkowicz and Buffa [25] and Demkowicz [22]. These logarithmic factors were removed in Bespalov and Heuer [6] and then in Melenk and Rojik [43] when working with weaker norms/higher regularity. In these references, in any case, $P_{hp}^{\text{div}}$ is not defined over the entire $H(\text{div}, \Omega)$. This is rectified in [28, Theorem 3.6].

## 1.3 Local-best-global-best and constrained-unconstrained equivalences

Following the seminal contribution by Veeser [48], with some predecessor results in Carstensen et al. [12, Theorem 2.1 and inequalities (3.2), (3.5), and (3.6)] and Aurada et al. [4, Proposition 3.1], there holds an equivalence between the best approximation of an $H^1(\Omega)$ function globally on the whole computational domain $\Omega$, with the trace continuity requirement, and locally on each mesh element, without any interelement continuity requirement. This result has been recently extended to the $H(\text{div}, \Omega)$-case in Ern et al. [28, Theorem 3.3] (see also Gawlik et al. [37] in the general context and [14, Theorems 1 and 2] and [17, Theorem 3.8] in the $H(\text{curl}, \Omega)$ context). Let $v \in H_{0,N}(\text{div}, \Omega)$, with, for simplicity for the moment, a piecewise polynomial divergence, $\nabla \cdot v \in P_p(\mathcal{T}_h)$. Then, the equivalence writes

$$ v_{hp} \in \min_{\substack{\nabla \cdot v_{hp} = \nabla \cdot v \\ K \in \mathcal{T}_h}} \|v - v_{hp}\|^2 \approx \sum_{v_p \in \mathcal{R}T_p(K)} \|v - v_p\|_K^2, \tag{1.2} $$

(global continuous constrained – local discontinuous unconstrained equivalence)

It is to be noted that $H_{0,N}(\text{div}, \Omega)$-conformity (normal trace is continuous over mesh faces and vanishes on $\Gamma_N \subset \partial\Omega$) and divergence constraints are requested on the left-hand side of (1.2), which is a global-best approximation over the entire $\Omega$. In contrast, crucially, the local-best approximation on the right-hand side of (1.2) is discontinuous and unconstrained. The generic equivalence constant in (1.2) from [28, Theorem 3.3] depends on the mesh shape-regularity and the space dimension $d$, but, unfortunately also (unfavorably) on the polynomial degree $p$. In the $H^1(\Omega)$ context, similar (algebraic) $p$-dependence is obtained in [48] and has been improved to logarithmic in two space dimensions in Canuto et al. [11, Theorem 4]. The concurrent work [50] establishes the equivalent of (1.2) in the $H^1(\Omega)$-case with a $p$-independent (robust) equivalence constant.
## 1.4 Main results of this manuscript

The main results of this manuscript is a construction of an operator $P_{hp}^{\text{div}}$ as in (1.1) such that: 1) it is defined over the entire $H_{0,N}(\text{div}, \Omega)$ and preserves boundary conditions imposed on the Neumann part $\Gamma_N$ of the boundary of $\Omega$; 2) it is defined locally in a neighborhood of each mesh element $K \in \mathcal{T}_h$; 3) it is based on elementwise $L^2$-orthogonal polynomial projections; 4) it is a projector, i.e., it leaves intact objects that are already in the Raviart-Thomas piecewise polynomial subspace of $H_{0,N}(\text{div}, \Omega)$,

$$ P_{hp}^{\text{div}}(\mathbf{v}) = \mathbf{v} \quad \forall \mathbf{v} \in \mathcal{RT}_p(\mathcal{T}_h) \cap H_{0,N}(\text{div}, \Omega); \qquad (1.3) \\ \text{(projection)} $$

5) it is locally and p-robustly stable in the $L^2$-norm, up to hp data oscillation,

$$ \|P_{hp}^{\text{div}}(\mathbf{v})\|_{K}^{2} \lesssim \sum_{L \in \bar{\mathcal{T}}_{K}} \left\{ \|\mathbf{v}\|_{L}^{2} + \left( \frac{h_{L}}{p+1} \|\nabla \cdot \mathbf{v} - \Pi_{hp}(\nabla \cdot \mathbf{v})\|_{L} \right)^{2} \right\}, \qquad (1.4) \\ \text{(L}^2\text{-stability up to data oscillation)} $$

where $\tilde{\mathcal{T}}_K$ is an extended element patch consisting of two layers of vertex neighbors of $K \in \mathcal{T}_h$; note that the second term on the above right-hand side (called hp data oscillation) vanishes if $\nabla \cdot \mathbf{v} \in \mathcal{P}_p(\mathcal{T}_h)$, yielding full $L^2$-stability in this case; 6) its approximation property is locally and p-robustly equivalent to that of the discontinuous unconstrained (elementwise $L^2$-orthogonal) projection:

$$ \|\mathbf{v} - P_{hp}^{\text{div}}(\mathbf{v})\|_{K}^{2} + \left( \frac{h_K}{p+1} \|\nabla \cdot (\mathbf{v} - P_{hp}^{\text{div}}(\mathbf{v}))\|_{K} \right)^{2} \lesssim \sum_{L \in \bar{\mathcal{T}}_K} \left\{ \min_{\mathbf{v}_p \in \mathcal{RT}_p(L)} \|\mathbf{v} - \mathbf{v}_p\|_{L}^{2} + \left( \frac{h_L}{p+1} \|\nabla \cdot \mathbf{v} - \Pi_{hp}(\nabla \cdot \mathbf{v})\|_{L} \right)^{2} \right\}; \qquad (1.5) \\ \text{(approximation equivalent to elementwise L}^2\text{-orthogonal projector)} $$

7) it satisfies the commuting property (1.1) with the $L^2$-orthogonal projection onto piecewise polynomials $\Pi_{hp}$. Crucially, the constant hidden in $\lesssim$ in inequalities (1.4) and (1.5) above only depends on the local mesh shape-regularity given by $\max_{L \in \bar{\mathcal{T}}_K} \kappa_L$ with $\kappa_L$ given by (2.1) below and on the space dimension $d$, in contrast to all references discussed in Section 1.1. All details are presented in Definition 3.3 and Theorem 3.4 below. In three space dimensions, we need Assumption 3.2 on the existence of a p-stable decomposition.

The properties of $P_{hp}^{\text{div}}$ immediately lead to two important consequences. Let $\mathbf{v} \in H_{0,N}(\text{div}, \Omega)$. First, (1.5) immediately implies (1.2) with the hidden constant independent of the polynomial degree $p$, crucially improving [28, Theorem 3.3]. The full version of this result, considering $\mathbf{v} \in H_{0,N}(\text{div}, \Omega)$ with general non-polynomial divergence, is stated in Theorem 3.5.

Second, we establish

$$ \|\mathbf{v} - P_{hp}^{\text{div}}(\mathbf{v})\|_{K}^{2} + \left( \frac{h_K}{p+1} \|\nabla \cdot (\mathbf{v} - P_{hp}^{\text{div}}(\mathbf{v}))\|_{K} \right)^{2} \lesssim \sum_{L \in \bar{\mathcal{T}}_K} \left\{ \left( \frac{h_L^{\min(s_L, p+1)}}{(p+1)^{s_L}} \|\mathbf{v}\|_{H^{s_L}(L)} \right)^2 + \left( \frac{h_L}{p+1} \frac{h_L^{\min(t_L, p+1)}}{(p+1)^{t_L}} \|\nabla \cdot \mathbf{v}\|_{H^{t_L}(L)} \right)^2 \right\} \qquad (1.6) \\ \text{(optimal elementwise hp approximation estimate)} $$

whenever the function $\mathbf{v}$ and its divergence $\mathbf{v}$ additionally have, separately on each mesh element $K \in \mathcal{T}_h$, the Sobolev regularity

$$ \mathbf{v}|_{K} \in H^{s_K}(K) \quad \text{and} \quad (\nabla \cdot \mathbf{v})|_{K} \in H^{t_K}(K) \qquad (1.7) $$

with Sobolev regularity exponents $s_K, t_K \ge 0$ (down to the minimal regularity $s_K = t_K = 0$). The bound (1.6) holds up to a constant that only depends on the mesh shape-regularity, the space dimension $d$, and the regularity exponents $s_K, t_K$; details form the content of Theorem 3.6. This improves the results discussed in Section 1.2 in several directions: no logarithmic factors in $p$ appear; no minimal regularity such as $\mathbf{v} \in H^s(\Omega)$ with $s > 0$ is imposed; no global regularity over the entire $\Omega$ or over patches appears:
(1.7) only requests additional Sobolev regularity separately on each mesh element $K \in \mathcal{T}_h$; in particular, (1.6) improves [28, Theorem 3.6] where the regularity exponents $s_K$ had to be constant over the entire computational domain $\Omega$ (and where a somewhat less sharper treatment of the divergence has been applied).

## 1.5 Crucial tools: polynomial extension operators and p-stable decompositions

There are two crucial tools used to obtain the above results. First, these are polynomial extension operators in the $H(\operatorname{div}, \Omega)$ context, namely that of Ainsworth and Demkowicz [2] for a normal trace lifting on a triangle, that of Demkowicz et al. [26] for a normal trace lifting on a tetrahedron (cf. also the recent work of Falk and Winther [35] for a $d$-simplex), and finally that of Costabel and McIntosh [21] for a divergence lifting on a $d$-simplex. We more precisely employ their broken extensions on patches of elements, obtained in Ern and Vohralík [33, Theorems 2.3 and 2.5, Corollaries 3.3 and 3.8], following Braess et al. [9]. We then generalize these results further to larger (extended) patches and no trace boundary conditions. Second, these are $p$-robustly stable decompositions, where we will namely use that of Schöberl et al. [47] in two space dimensions.

## 1.6 Organization of this manuscript

We set up the notation in Section 2. We then present our main results in full details in Section 3, also including a quick numerical illustration. The more involved proofs are subsequently collected in Sections 4 and 5. We finally present four independent results in the appendices. We first formulate the $p$-stable decomposition result from [47] in a form suitable for us in Appendix A. Next, we introduce the notion of suitable patch enumeration and show its equivalence with shellability of simplicial complexes in Appendix B. We then generalize the results from [33] to larger (extended) patches and no trace boundary conditions in Appendix C. Finally, in Appendix D, we similarly extend the results of [15, Appendix A] concerning seemingly overconstrained minimizations on patch subdomains.

This contribution only concerns the $H(\operatorname{div}, \Omega)$ case. The $H^1$ case is studied in Vohralík [50], whereas the $H(\operatorname{curl}, \Omega)$ case will be addressed in Vohralík [51]. Both these references study locally varying polynomial degrees. For the sake of readability, we only present here the uniform polynomial degree case; all the present results extend to the varying polynomial degree case as in [50, 51].

# 2 Setting and notation

We set here the context and notation.

## 2.1 Domain $\Omega$, simplicial mesh $\mathcal{T}_h$, and patch subdomains $\omega$

Let the computational domain $\Omega \subset \mathbb{R}^d$, $d = 2, 3$, be an open, bounded, and connected Lipschitz polygon or polyhedron. Let $\mathcal{T}_h$ be a simplicial mesh of $\Omega$, i.e., a collection of nontrivial closed triangles or tetrahedra $K$ covering $\bar{\Omega}$, where the intersection of two different simplices is either empty or their common vertex, edge, or face. The shape-regularity parameters of the element and of the entire mesh $\mathcal{T}_h$ are respectively given by

$$ \kappa_K := \frac{h_K}{\rho_K}, \quad \kappa_h := \max_{K \in \mathcal{T}_h} \kappa_K, \qquad (2.1) $$

where $h_K$ is the diameter of the simplex $K$ and $\rho_K$ that of the largest ball contained in $K$. Uniformly bounded $\kappa_h$ allows for families of strongly graded meshes with local refinements but not for anisotropic elements. Let the piecewise constant mesh-size function $h$ be given by $h_K$ on each $K \in \mathcal{T}_h$. Below, we reserve the notation $\omega \subset \mathbb{R}^d$, possibly with subscripts, for open, bounded, Lipschitz, and polygonal or polyhedral subdomains of $\Omega$ corresponding to a set of mesh elements from $\mathcal{T}_h$ such that $\bar{\omega}$ is contractible (homotopic to a ball); we stress this in Assumption 3.1 below. The diameter of $\omega$ is denoted by $h_\omega$.

## 2.2 Vertices, edges, faces, and patches of mesh elements

For a simplex $K \in \mathcal{T}_h$, denote by $\mathcal{V}_K$ the set of its vertices, and let $\mathcal{V}_h$ collect all mesh vertices. Generic vertices will be denoted by $\mathbf{a}$ and $\mathbf{b}$. We will also work with mesh faces $F$, where, henceforth, “face” means “($d-1$)-dimensional face”, i.e., a face in three space dimensions and an edge in two space dimensions.
<figure><img src="image_2.png" /><figcaption>Figure 1: Vertex patch $\mathcal{T}_a$ for a vertex $a \in V_h$ in the interior of $\Omega$ (left) and on the boundary of $\Omega$ (right), $d = 2$</figcaption></figure>

For a vertex $a \in V_h$, denote by $\mathcal{T}_a$ the patch of the elements of $\mathcal{T}_h$ that share $a$ and $\omega_a$ the corresponding open subdomain. Illustration in two space dimensions is provided in Figure 1. We will also need the extended vertex patch $\tilde{\mathcal{T}}_a$ and the corresponding subdomain $\tilde{\omega}_a$; this includes $\mathcal{T}_a$ and elements of all vertex patches $\mathcal{T}_b$ of vertices $b$ from $\mathcal{T}_a$, see Figure 2 (left). Equivalently, $\tilde{\mathcal{T}}_a$ is formed by those elements $L$ from the mesh $\mathcal{T}_h$ that share at least a vertex with an element $K \in \mathcal{T}_a$. Similarly, for a simplex $K \in \mathcal{T}_h$, let $\tilde{\mathcal{T}}_K$ be the extended element patch given by the union of $\tilde{\mathcal{T}}_a$ over all vertices $a$ of the simplex $K$; this comprises $K$ and all elements $L$ sharing a vertex with $K$ or with its vertex neighbor. The corresponding subdomain is denoted by $\tilde{\omega}_K$; an illustration is provided in Figure 2 (right). There is a variety of scenarios that might occur; for instance, for a vertex/element in the interior of $\Omega$, the (extended) vertex/element patch may touch the boundary $\partial\Omega$/be “cropped” by the boundary $\partial\Omega$. All these cases are covered in our construction. In the sequel, we will need to collect the vertices from respectively $\tilde{\mathcal{T}}_a$ and $\tilde{\mathcal{T}}_K$ in the sets $\tilde{V}_a$ and $\tilde{V}_K$. Diameters of respectively $\omega_a$, $\tilde{\omega}_a$, and $\tilde{\omega}_K$ are denoted by $h_{\omega_a}$, $h_{\tilde{\omega}_a}$, and $h_{\tilde{\omega}_K}$.

## 2.3 Hat functions and the partition of unity

Let $a \in V_h$ be an arbitrary mesh vertex. Then the continuous, piecewise first-order polynomial (affine) “hat” function $\psi^a$ takes value 1 at the vertex $a$ and zero at all the other vertices. We note that $\omega_a$ corresponds to the support of $\psi^a$, and that these functions form the partition of unity in that

$$ \sum_{a \in V_h} \psi^a = 1. \qquad (2.2) $$

## 2.4 Boundary subsets $\Gamma_D$ and $\Gamma_N$

Let $\Gamma_D$ and $\Gamma_N$ be two disjoint, relatively open, and possibly empty subsets of the computational domain boundary $\partial\Omega$ such that $\partial\Omega = \overline{\Gamma_D} \cup \overline{\Gamma_N}$. We also require that $\Gamma_D$ and $\Gamma_N$ have polygonal Lipschitz boundaries and we assume that each boundary face of the mesh $\mathcal{T}_h$ lies entirely either in $\overline{\Gamma_D}$ or in $\overline{\Gamma_N}$.

## 2.5 The spaces $H(\text{div})$ on the entire computational domain and on its subdomains

Let $\omega \subseteq \Omega$ be as in Section 2.1. We let $L^2(\omega)$ be the space of scalar-valued square-integrable functions defined on $\omega$. We denote by $(v, w)_\omega := \int_\omega v(x)w(x)dx$ the $L^2(\omega)$ scalar product and by $\|\cdot\|_\omega$ the corresponding norm; we drop the index when $\omega = \Omega$. We also use the notation $L^2(\omega) := [L^2(\omega)]^d$ for vector-valued functions with each component in $L^2(\omega)$. This is equipped with the scalar product $(v, v)_\omega := \int_\omega v(x) \cdot w(x) dx$ and the corresponding norm. We again drop the index when $\omega = \Omega$. The central space for this study is $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)\}$, see Girault and Raviart [38], Ern and Guermond [31],
<figure><img src="image_3.png" /><figcaption>Figure 2: Extended vertex patch $\tilde{T}_a$ for a vertex $a \in V_h$ sufficiently in the interior of $\Omega$ (generated by the vertex patches $T_b$ of all vertices $b$ from $T_a$, marked by a square or a circle) (left) and extended element patch $\tilde{T}_K$ for an element $K \in T_h$ on the boundary of $\Omega$ (generated by the vertex patches $T_b$ of all vertices $b$ marked by a square or a circle) (right), $d = 2$</figcaption></figure>

or Demkowicz [23]. We will employ the notation $\langle \cdot, \cdot \rangle_S$ for the integral product on boundary (sub)sets $S \subset \partial\omega$ or on mesh faces $F$, as well as for duality pairing when $S = \partial\omega$.

Let $n_\omega$ be the unit normal vector on $\partial\omega$, outward to $\omega$. If $\partial\omega$ does not contain any face from $\bar{\Gamma}_N$, cf. Figure 2 (left), let $H_{0,N}(\text{div}, \omega) := H(\text{div}, \omega)$. In general, cf. Figure 2 (right) for an example of $\partial\omega$ containing faces from $\bar{\Gamma}_N$, we let $H_{0,N}(\text{div}, \omega)$ be the subspace of $H(\text{div}, \omega)$ formed by functions with vanishing normal trace on the faces in $\partial\omega \cap \bar{\Gamma}_N$,

$$ H_{0,N}(\text{div}, \omega) := \{v \in H(\text{div}, \omega); v \cdot n_\omega = 0 \text{ on } (\partial\omega \cap \bar{\Gamma}_N)^\circ\}, \quad (2.3) $$

which is understood by duality,

$$ v \cdot n_{\omega} = 0 \text{ on } (\partial\omega \cap \Gamma_N)^{\circ} \iff (v, \nabla\varphi)_{\omega} + (\nabla \cdot v, \varphi)_{\omega} = 0 \quad \forall \varphi \in H_{0,\partial\omega\setminus\Gamma_N}^{1}(\omega). \quad (2.4) $$

Here $H_{0,\partial\omega\setminus\Gamma_N}^1(\omega)$ stands for all functions $\varphi$ from the first-order Sobolev space $H^1(\omega)$ which vanish on the interior of $\partial\omega \setminus \Gamma_N$ in the sense of traces.

Finally, for a vertex patch subdomain $\omega_a$, cf. Figure 1, we will employ the notation $H_{0,N,\psi^a}(\text{div}, \omega_a)$ for the subspace of $H(\text{div}, \omega_a)$ with zero normal trace on those faces in $\partial\omega_a$ where the hat function $\psi^a$ vanishes (all $\partial\omega_a$ for interior vertices) and which lie in the Neumann boundary $\bar{\Gamma}_N$,

$$ H_{0,N,\psi^a}(\text{div}, \omega_a) := \{v \in H(\text{div}, \omega_a); v \cdot n_{\omega_a} = 0 \text{ on } \partial\omega_a \cap \{\psi^a = 0\} \\ \text{and } (\partial\omega_a \cap \Gamma_N)^{\circ}\}. \quad (2.5) $$

In Figure 1, this respectively corresponds to the double line (for interior patches $T_a$, left) or to the double and zigzag lines (for boundary patches $T_a$, right). Similarly, for an arbitrary patch subdomain $\omega$ and a vertex $a \in \bar{\omega}$, $H_{0,N,\psi^a}(\text{div}, \omega_a \cap \omega)$ stands for the subspace of $H(\text{div}, \omega_a \cap \omega)$ with zero normal trace on those faces in $\partial(\omega_a \cap \omega)$ where the hat function $\psi^a$ vanishes or which lie in the Neumann boundary $\bar{\Gamma}_N$,

$$ H_{0,N,\psi^a}(\text{div}, \omega_a \cap \omega) := \{v \in H(\text{div}, \omega_a \cap \omega); v \cdot n_{\omega_a \cap \omega} = 0 \\ \text{on } \partial(\omega_a \cap \omega) \cap \{\psi^a = 0\} \text{ and } (\partial(\omega_a \cap \omega) \cap \Gamma_N)^{\circ}\}. \quad (2.6) $$

This is as above in (2.5) with the exception of vertices $a$ on the boundary of $\omega$: the normal trace of the functions from $H_{0,N,\psi^a}(\text{div}, \omega_a \cap \omega)$ does not have to vanish on $\partial(\omega_a \cap \omega)$ unless this is a part of $\bar{\Gamma}_N$, see $\omega_{b2}$ highlighted by green north-western lines in Figure 3.

## 2.6 Piecewise polynomial spaces

Let $p \ge 0$ be a fixed polynomial degree. For a single simplex $K \in T_h$, we denote by $P_p(K)$ the space of scalar-valued polynomials on $K$ of total degree at most $p$. The notation $[P_p(K)]^d$ then stands for the
space of vector-valued polynomials on K with each component in $\mathcal{P}_p(K)$. The Raviart-Thomas [8, 24, 46] space of degree p on K is given by

$$ \mathcal{RT}_p(K) := [\mathcal{P}_p(K)]^d + \mathcal{P}_p(K)\mathbf{x} = [\mathcal{P}_p(K)]^d \oplus \tilde{\mathcal{P}}_p(K)\mathbf{x}, \qquad (2.7) $$

where $\tilde{\mathcal{P}}_p(K)$ stands for homogeneous polynomials of degree p on K.

We will below extensively use the broken, piecewise polynomial spaces formed from the above element spaces

$$ \begin{aligned} \mathcal{P}_p(\mathcal{T}_h) &:= \{v_{hp} \in L^2(\Omega); v_{hp}|_K \in \mathcal{P}_p(K) \quad \forall K \in \mathcal{T}_h\}, \\ \mathcal{RT}_p(\mathcal{T}_h) &:= \{v_{hp} \in L^2(\Omega); v_{hp}|_K \in \mathcal{RT}_p(K) \quad \forall K \in \mathcal{T}_h\}. \end{aligned} $$

To form the usual finite element $H(\text{div}, \Omega)$-conforming, normal-trace-continuous spaces, we will write $\mathcal{RT}_p(\mathcal{T}_h) \cap H(\text{div}, \Omega)$ and similarly for the subspaces reflecting the different boundary conditions. The same notation will also be used on the patches $\mathcal{T}_a, \tilde{\mathcal{T}}_a$, and $\tilde{\mathcal{T}}_K$.

## 2.7 $L^2$-orthogonal projector onto piecewise polynomials and the elementwise canonical Raviart-Thomas projector

Let $\Pi_{hp}$ denote the elementwise $L^2(\Omega)$-orthogonal projector onto $\mathcal{P}_p(\mathcal{T}_h)$: for $v \in L^2(\Omega)$, $\Pi_{hp}(v) \in \mathcal{P}_p(\mathcal{T}_h)$ is prescribed locally on each element $K \in \mathcal{T}_h$, $\Pi_{hp}(v)|_K \in \mathcal{P}_p(K)$, by

$$ (\Pi_{hp}(v), v_p)_K = (v, v_p)_K \quad \forall v_p \in \mathcal{P}_p(K). \qquad (2.8) $$

Next, we will use the elementwise canonical p-degree Raviart-Thomas projector $\mathcal{I}_{h,p}^{\mathcal{RT}}$: for $v \in \Pi_{K \in \mathcal{T}_h} [C^1(K)]^d$, a function of the $C^1$ regularity in each component, separately on each mesh element $K \in \mathcal{T}_h$, $\mathcal{I}_{h,p}^{\mathcal{RT}}(v) \in \mathcal{RT}_p(\mathcal{T}_h)$ is following [8, 24, 46] prescribed locally on each element $K \in \mathcal{T}_h$, $\mathcal{I}_{h,p}^{\mathcal{RT}}(v)|_K \in \mathcal{RT}_p(K)$, by

$$ \langle \mathcal{I}_{h,p}^{\mathcal{RT}}(v) \cdot \mathbf{n}_K, r_p \rangle_F = \langle v \cdot \mathbf{n}_K, r_p \rangle_F \quad \forall r_p \in \mathcal{P}_p(F), \quad \text{for all faces } F \text{ of } K, \qquad (2.9a) $$

$$ (\mathcal{I}_{h,p}^{\mathcal{RT}}(v), r_p)_K = (v, r_p)_K \quad \forall r_p \in [\mathcal{P}_{p-1}(K)]^d, \qquad (2.9b) $$

where $\mathbf{n}_K$ is the unit outer normal vector of the element K. This projector, crucially, satisfies the commuting property, locally on each $K \in \mathcal{T}_h$,

$$ \nabla \cdot \mathcal{I}_{h,p}^{\mathcal{RT}}(v) = \Pi_{hp}(\nabla \cdot v) \quad \forall v \in \Pi_{K \in \mathcal{T}_h} [C^1(K)]^d. \qquad (2.10) $$

We will only apply $\mathcal{I}_{h,p}^{\mathcal{RT}}$ to piecewise (discontinuous) polynomials which have the requested elementwise $[C^1(K)]^d$ regularity; recall from [8, 24] that one cannot use $\mathcal{I}_{h,p}^{\mathcal{RT}}$ directly on $H_{0,N}(\text{div}, \Omega)$.

## 2.8 Notation $\lesssim$

We will use the notation $a \lesssim b$ when there holds $a \le Cb$ for a positive constant $C$ and $a \approx b$ when $a \lesssim b$ and $b \lesssim a$ hold simultaneously. All dependencies of the hidden constant $C$ will systematically be given. In any case, all constants in this manuscript are independent of the mesh size $h$ and of the polynomial degree $p$.

# 3 Main results

We present here our main results. We rely on two assumptions:

**Assumption 3.1 (Patch subdomains).** For any vertex $\mathbf{a} \in \mathcal{V}_h$ and element $K \in \mathcal{T}_h$, consider the extended vertex patch $\tilde{\mathcal{T}}_a$ or the extended element patch $\tilde{\mathcal{T}}_K$ as per Section 2.2, denoted by $\mathcal{T}_\omega$, with the associated open subdomain $\omega$. We suppose that $\omega$ is Lipschitz, the closure $\bar{\omega}$ of $\omega$ is contractible, and the boundary of $\omega$ does not coincide with the whole Neumann boundary, $\partial\omega \neq \partial\omega \cap \Gamma_N$.

**Assumption 3.2 (A p-stable $\mathcal{RT}_p \cap H(\text{div})$ decomposition on three-dimensional patches).** Suppose that Theorem A.1 also holds for $d = 3$.
Assumption 3.1 is similar to those in [34, 3]. Please remark that it does not request the whole computational domain $\Omega$ to be contractible but merely the patch subdomains. For example for $\Omega$ with a hole, Assumption 3.1 may not be satisfied for a coarse triangulation, but typically will be satisfied for a finer mesh. The same holds true for the assumption that $\partial\hat{\omega}_K$ does not coincide with the whole Neumann boundary. We refer for further details to the recent discussion in [32, Remark 2.1].

We only make Assumption 3.2 in three space dimensions. In two space dimensions, the p-robustly stable $H(\text{div})$ patch decomposition of Theorem A.1 is a simple consequence of [47]. In three space dimensions, Assumption 3.2 holds with the constant in (A.5) below possibly depending on the polynomial degree $p$, which is shown in [15, Appendix B]. Recently, [36] has established the extension of the result of [47] to any space dimension and any differential form, including the result in $H(\text{div})$. Unfortunately, we crucially need a stable decomposition with vanishing global low-order component, which does not seem to be easily available from [36], so that we at present need to rely on Assumption 3.2.

## 3.1 A p-stable local commuting projector in $H_{0,N}(\text{div}, \Omega)$

We first define our p-stable local commuting projector in $H_{0,N}(\text{div}, \Omega)$ and state its properties.

### 3.1.1 Definition of the projector

Our construction extends and builds on [28, Definition 3.1] for equilibration and on [15, Appendix A] for imposing of an additional orthogonality constraint that enables to employ the p-stable decomposition of [47] in a correction stage. Prior to stating it, let us recall the basic steps from [28] and highlight the differences.

The construction in [28, Definition 3.1] proceeds in three steps: 1) elementwise $L^2$-orthogonal projection (local-best approximation) (with a divergence constraint); 2) patchwise equilibration; this crucially employs the hat functions $\psi^a$ from (2.2) and the canonical projector $I_{h,p}^{RT}$ from (2.9) (which in turn prevents proving a p-robust local-best and global-best equivalence as in Theorem 3.5); and 3) gluing of the patchwise contributions. The present construction is slightly more involved but leads to better approximation properties, namely yielding the p-robust local-best and global-best equivalence of Theorem 3.5 and p-robust approximation property (3.14) below. It proceeds in four stages: 1) elementwise $L^2$-orthogonal projection (local-best approximation) (without the divergence constraint); 2) patchwise equilibration and gluing of the patchwise contributions, like above in steps 2) and 3), but with an additional orthogonality constraint; this stage still employs the hat functions $\psi^a$ from (2.2) as well as the canonical projector $I_{h,p}^{RT}$ from (2.9); its main purpose is to cut off the divergence and to impose an elementwise $L^2$-orthogonality with respect to constant vectors; 3) patchwise equilibration of the remainder (with an additional orthogonality constraint) followed by a p-stable decomposition of Appendix A and gluing of the patchwise contributions into a correction; here, crucially, no hat functions $\psi^a$ from (2.2) and no projector such as $I_{h,p}^{RT}$ from (2.9) are used; and 4) combination of the previous steps.

Recall the notation from Section 2. The definition reads:

**Definition 3.3** (A p-stable local commuting projector in $H_{0,N}(\text{div}, \Omega)$). Let a simplicial mesh $\mathcal{T}_h$ of $\Omega$ and a polynomial degree $p \ge 0$ be given. We define

$$ P_{hp}^{\text{div}} : H_{0,N}(\text{div}, \Omega) \to \mathcal{RT}_p(\mathcal{T}_h) \cap H_{0,N}(\text{div}, \Omega) \quad (3.1) $$
*(defined over the entire $H_{0,N}(\text{div}, \Omega)$)*

as follows. Let a function $v \in H_{0,N}(\text{div}, \Omega)$ be given.

**I.** On each mesh element $K \in \mathcal{T}_h$, consider the $L^2(K)$-orthogonal projection of $v$ onto $\mathcal{RT}_p(K)$ (without any normal trace prescription nor any constraint on the divergence)

$$ \tau_{hp}|_K := \arg \min_{v_p \in \mathcal{RT}_p(K)} \|v - v_p\|_K. \quad (3.2) $$
*(elementwise projection $\tau_{hp}$)*

This gives the broken Raviart-Thomas piecewise polynomial

$$ \tau_{hp} \in \mathcal{RT}_p(\mathcal{T}_h). \quad (3.3) $$

**II. Starting from $\tau_{hp}$:**
(a) On each vertex patch $T_a$, $a \in V_h$, see Figure 1, define the Raviart-Thomas piecewise polynomial $\sigma_p^a \in \mathcal{RT}_p(T_a) \cap H_{0,N,\psi^a}(\text{div}, \omega_a)$ via

$$ \sigma_p^a := \arg \min_{\substack{v_p \in \mathcal{RT}_p(T_a) \cap H_{0,N,\psi^a}(\text{div}, \omega_a) \\ \nabla \cdot v_p = \Pi_{hp}(\psi^a \nabla \cdot v + \nabla \psi^a \cdot v) \\ (v_p, r_h)_K = (I_{h,p}^{\mathcal{RT}}(\psi^a \tau_{hp}), r_h)_K \quad \forall r_h \in [P_0(K)]^d, \forall K \in T_a}} \| I_{h,p}^{\mathcal{RT}}(\psi^a \tau_{hp}) - v_p \|_{\omega_a}; \quad (3.4a) $$

(patchwise “no flux” equilibration, with an additional orthogonality constraint if $p \ge 1$)

recall from (2.5) that $H_{0,N,\psi^a}(\text{div}, \omega_a)$ is the subspace of $H(\text{div}, \omega_a)$ with zero normal trace on those faces in $\partial\omega_a$ where the hat function $\psi^a$ vanishes or which lie in the Neumann boundary $\bar{\Gamma}_N$. See Figure 3 (left) for illustration.

(b) Extending $\sigma_p^a$ by zero outside of the patch subdomain $\omega_a$, define

$$ \sigma_{hp} := \sum_{a \in V_h} \sigma_p^a. \quad (3.4b) $$

(gluing patchwise contributions)

This gives the intermediate Raviart-Thomas piecewise polynomial

$$ \sigma_{hp} \in \mathcal{RT}_p(T_h) \cap H_{0,N}(\text{div}, \Omega) \text{ with } \nabla \cdot \sigma_{hp} = \Pi_{hp}(\nabla \cdot v) \quad (3.5) $$

(seemingly overconstrained if $p \ge 1$) equilibration $\sigma_{hp}$ with divergence and projection properties

and the broken Raviart-Thomas piecewise polynomial

$$ \tau_{hp} - \sigma_{hp} \in \mathcal{RT}_p(T_h) \quad (3.6a) $$

$$ (\tau_{hp} - \sigma_{hp}, r_h)_K = 0 \quad \forall r_h \in [P_0(K)]^d, \forall K \in T_h \quad \text{if } p \ge 1. \quad (3.6b) $$

(remainder $\tau_{hp} - \sigma_{hp}$, with vanishing lowest-order moments if $p \ge 1$)

### III. If $p=0$, set $\zeta_{hp} := 0$. Otherwise, if $p \ge 1$, starting from $\tau_{hp} - \sigma_{hp}$:

(a) On each extended vertex patch $\tilde{T}_a$, $a \in V_h$, see Figure 2 (left), define the Raviart-Thomas piecewise polynomial $\zeta_p^a \in \mathcal{RT}_p(\tilde{T}_a) \cap H_{0,N}(\text{div}, \tilde{\omega}_a)$ via

$$ \zeta_p^a := \arg \min_{\substack{v_p \in \mathcal{RT}_p(\tilde{T}_a) \cap H_{0,N}(\text{div}, \tilde{\omega}_a) \\ \nabla \cdot v_p = 0}} \| \tau_{hp} - \sigma_{hp} - v_p \|_{\tilde{\omega}_a}; \quad (3.7a) $$

(seemingly overconstrained divergence-free remainder equilibration)

recall from (2.3) that $H_{0,N}(\text{div}, \tilde{\omega}_a)$ is the subspace of $H(\text{div}, \tilde{\omega}_a)$ with zero normal trace on those boundary faces in $\partial\tilde{\omega}_a$ which lie in $\bar{\Gamma}_N$. See Figure 3 (right) for illustration.

(b) On each extended vertex patch $\tilde{T}_a$, $a \in V_h$, employ to $\zeta_p^a$ the p-stable decomposition of Theorem A.1 or Assumption 3.2 (with $T_\omega = \tilde{T}_a$ and $V_\omega = \tilde{V}_a$),

$$ \zeta_p^a = \sum_{b \in \tilde{V}_a} \zeta_p^{a,b} \quad \text{with in particular } \zeta_p^{a,a} \in \mathcal{RT}_p(T_a) \cap H_{0,N,\psi^a}(\text{div}, \omega_a), \nabla \cdot \zeta_p^{a,a} = 0. \quad (3.7b) $$

(patchwise p-stable equilibrated remainder decomposition)

See Figure 3 (right) for illustration.

(c) Extending the “interior” component $\zeta_p^{a,a}$ by zero outside of the patch subdomain $\omega_a$, define

$$ \zeta_{hp} := \sum_{a \in V_h} \zeta_p^{a,a}. \quad (3.7c) $$

(gluing patchwise correction contributions)

This gives the intermediate Raviart-Thomas piecewise polynomial

$$ \zeta_{hp} \in \mathcal{RT}_p(T_h) \cap H_{0,N}(\text{div}, \Omega) \text{ with } \nabla \cdot \zeta_{hp} = 0. \quad (3.8) $$

(p-robust correction $\zeta_{hp}$ by treatment of the remainder $\tau_{hp} - \sigma_{hp}$ without $\psi^a$ and $I_{h,p}^{\mathcal{RT}}$)
<figure><img src="image_4.png" /><figcaption>Figure 3: The standard non-p-robust equilibration component $\sigma_p^a$ from (3.4a) (left) and the p-robust correction $\zeta_p^a$ from (3.7a) together with its p-stable decomposition (3.7b); only the “interior” component $\zeta_p^{a,a}$ is used (right); $d = 2$, interior of the domain</figcaption></figure>

## IV. Define

$$ P_{hp}^{\text{div}}(\mathbf{v}) := \sigma_{hp} + \zeta_{hp}. \qquad (3.9) $$
$$ (\text{combining the previous steps}) $$

This gives the final Raviart-Thomas piecewise polynomial

$$ P_{hp}^{\text{div}}(\mathbf{v}) \in \mathcal{RT}_p(\mathcal{T}_h) \cap \mathcal{H}_{0,N}(\text{div}, \Omega) \quad \text{with} \quad \nabla \cdot P_{hp}^{\text{div}}(\mathbf{v}) = \Pi_{hp}(\nabla \cdot \mathbf{v}). \qquad (3.10) $$

We will verify the correctness of Definition 3.3 in Section 4 below.

### 3.1.2 Design principles

Let us discuss in detail the design principles of Definition 3.3.

1. The construction of $\tau_{hp}$ in Step I. sets our local-best discontinuous unconstrained projection “target”. There holds $\tau_{hp} \in \mathcal{RT}_p(\mathcal{T}_h)$ but in general $\tau_{hp} \notin \mathcal{H}_{0,N}(\text{div}, \Omega)$. In the rest of Definition 3.3, we search to stay in $\mathcal{RT}_p(\mathcal{T}_h)$, as close as possible to $\tau_{hp}$, keeping its approximation power, but recovering $\mathcal{H}_{0,N}(\text{div}, \Omega)$-conformity.

2. The construction of $\sigma_{hp}$ in Step II. is similar to [28, Definition 3.1, steps 2–3], with the incorporation of the additional orthogonality constraint from [15, Appendix A] if $p \ge 1$. The proof of its p-robustness is obstructed by the presence of the cut-off by the hat functions $\psi^a$ from (2.2) and by the use of the canonical elementwise projector $I_{h,p}^{\mathcal{RT}}$ from (2.9), which brings the polynomial degree increased by $\psi^a$ to $p+1$ back down to $p$. The purpose here is to design a projector capturing the correct divergence as per (3.5) and to obtain the “remainder” $\tau_{hp} - \sigma_{hp}$ with vanishing lowest-order moments as per (3.6b), if $p \ge 1$. The projection property is actually already established here, as $\sigma_{hp}$ is such that if $\mathbf{v} \in \mathcal{RT}_p(\mathcal{T}_h) \cap \mathcal{H}_{0,N}(\text{div}, \Omega)$, then $\tau_{hp} = \sigma_{hp} = \mathbf{v}$ and the remainder $\tau_{hp} - \sigma_{hp}$ vanishes.

3. The construction of $\zeta_{hp}$ in Step III., only nontrivial if $p \ge 1$, is the salient feature for the theoretical proof of p-robustness. Neither the hat functions $\psi^a$ nor the elementwise projector $I_{h,p}^{\mathcal{RT}}$ are present. At the first stage in (3.7a), we employ an equilibration similar to (3.4a) which however 1) does not employ the hat functions $\psi^a$ from (2.2) or the canonical elementwise projector $I_{h,p}^{\mathcal{RT}}$ from (2.9); 2) is divergence-free; and 3) does not impose zero normal trace on $\partial\tilde{\omega}_a$ (except for $(\partial\tilde{\omega}_a \cap \Gamma_N)^\circ)$. At the second stage (3.7b), a p-stable decomposition is applied (note that this cannot be applied directly to the remainder $\tau_{hp} - \sigma_{hp}$ which “broken”, i.e., lies in $\mathcal{RT}_p(\tilde{\mathcal{T}}_a)$ but not in $\mathcal{H}_{0,N}(\text{div}, \tilde{\omega}_a)$). At this stage, the additional orthogonality constraint in (3.7a) (note that $(\tau_{hp} - \sigma_{hp}, \mathbf{r}_h)_K = 0$ follows from (3.6b)) plays a crucial role since it enables to employ the p-stable decomposition of [47] with vanishing lowest-order moments. Note that we merely access the integral volumetric (lowest-order) moments $(v_p, \mathbf{r}_h)_K$ which are available under the the $\mathcal{H}(\text{div}, K)$ regularity, in contrast to the (lowest-order) normal trace face moments such as $\langle v_p \cdot \mathbf{n}, 1 \rangle_F$ (any use of trace face moments
is only possible at the discrete level and also typically spoils p-robustness). Note that in (3.7c), we merely employ the “interior” or “middle” components which do have zero normal trace on $\partial\omega_a$ (for interior vertices) or on $\partial\omega_a \cap \{\psi^a = 0\}$ and $(\partial\omega_a \cap \Gamma_N)^\circ$ (for boundary vertices) as per the definition of $H_{0,N,\psi^a}(\text{div}, \omega_a)$ in (2.5), see Figure 3 for illustration.

4. In Step IV., $P_{hp}^{\text{div}}(v)$ is defined as $\sigma_{hp}$ corrected by $\zeta_{hp}$.

5. The construction relies on local energy minimization problems (3.2), (3.4a), (3.7a) and the p-stable decomposition (3.7b).

6. In comparison to [50, Definition 3.5], the orthogonality constraints with respect to vector-valued piecewise constants are imposed directly in the local minimization problems (3.4a) and (3.7a) and not in a correction stage after local minimization. This seems compulsory to satisfy the (divergence) constraint, not present in [50].

### 3.1.3 Properties of the projector

The following theorem summarizes the properties of the projector from Definition 3.3, improving the results in [25, 22, 18, 6, 34, 29, 30, 41, 43, 3, 37, 28].

**Theorem 3.4** (Commutativity, projection, approximation, and stability of $P_{hp}^{\text{div}}$). Let a simplicial mesh $\mathcal{T}_h$ of $\Omega$ and a polynomial degree $p \ge 0$ be given. Let Assumptions 3.1 and 3.2 hold. The operator $P_{hp}^{\text{div}}: H_{0,N}(\text{div}, \Omega) \to \mathcal{RT}_p(\mathcal{T}_h) \cap H_{0,N}(\text{div}, \Omega)$ from Definition 3.3 satisfies

$$ \nabla \cdot P_{hp}^{\text{div}}(v) = \Pi_{hp}(\nabla \cdot v) \quad \forall v \in H_{0,N}(\text{div}, \Omega), \qquad (3.11) $$
$$ (\text{commutativity}) $$

$$ P_{hp}^{\text{div}}(v) = v \quad \forall v \in \mathcal{RT}_p(\mathcal{T}_h) \cap H_{0,N}(\text{div}, \Omega). \qquad (3.12) $$
$$ (\text{projection}) $$

Moreover, for any function $v \in H_{0,N}(\text{div}, \Omega)$ and any mesh element $K \in \mathcal{T}_h$, there holds

$$ \|v - P_{hp}^{\text{div}}(v)\|_K^2 + \left(\frac{h_K}{p+1} \|\nabla \cdot (v - P_{hp}^{\text{div}}(v))\|_K\right)^2 \lesssim \sum_{L \in \tilde{\mathcal{T}}_K} \left\{ \min_{v_p \in \mathcal{RT}_p(L)} \|v - v_p\|_L^2 + \left(\frac{h_L}{p+1} \|\nabla \cdot v - \Pi_{hp}(\nabla \cdot v)\|_L\right)^2 \right\}, \qquad (3.13) $$
$$ (\text{approximation equivalent to elementwise } \mathbf{L}^2\text{-orthogonal projector}) $$

$$ \|P_{hp}^{\text{div}}(v)\|_K^2 \lesssim \sum_{L \in \tilde{\mathcal{T}}_K} \left\{ \|v\|_L^2 + \left( \frac{h_L}{p+1} \|\nabla \cdot v - \Pi_{hp}(\nabla \cdot v) \|_L \right)^2 \right\}, \qquad (3.15) $$

(L$^2$-stability up to data oscillation)

$$ \|P_{hp}^{\text{div}}(v)\|_K^2 + h_{\Omega}^2 \|\nabla \cdot P_{hp}^{\text{div}}(v)\|_K^2 \lesssim \sum_{L \in \tilde{\mathcal{T}}_K} \{\|v\|_L^2 + h_{\Omega}^2 \|\nabla \cdot v\|_L^2\}, \qquad (3.16) $$

(H(div)-stability)

where, recall from Section 2.2, $\tilde{\mathcal{T}}_K$ collects the elements $L$ of $\mathcal{T}_h$ sharing a vertex with $K$ or with its vertex neighbor. In (3.16), we employ the dimensionally correct scaling by $h_\Omega$, the diameter of $\Omega$ (element diameters $h_L$ can also be used). The constant hidden in $\lesssim$ only depends on the local mesh shape-regularity given by $\max_{L \in \tilde{\mathcal{T}}_K} \kappa_L$ with $\kappa_L$ given by (2.1) and the space dimension $d$.

## 3.2 p-robust equivalence of global continuous constrained and local discontinuous unconstrained approximation in $H_{0,N}(\text{div}, \Omega)$

The following result improves importantly [28, Theorem 3.3], removing the possible dependence of the equivalence constant on the polynomial degree $p$. It is an immediate consequence of Definition 3.3 and Theorem 3.4.
**Theorem 3.5** (p-robust equivalence of local-best and global-best approximations in $H_{0,N}(\text{div}, \Omega)$). Let $v \in H_{0,N}(\text{div}, \Omega)$, a simplicial mesh $\mathcal{T}_h$ of $\Omega$, and a polynomial degree $p \ge 0$ be given. Let Assumptions 3.1 and 3.2 hold. Then

$$ \begin{aligned} & \min_{\substack{v_{hp} \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_h) \cap H_{0,N}(\text{div}, \Omega) \\ \nabla \cdot v_{hp} = \Pi_{hp}(\nabla \cdot v)}} \|v - v_{hp}\|^2 + \sum_{K \in \mathcal{T}_h} \left( \frac{h_K}{p+1} \|\nabla \cdot v - \Pi_{hp}(\nabla \cdot v)\|_K \right)^2 \\ & \approx \sum_{K \in \mathcal{T}_h} \left\{ \min_{v_p \in \mathcal{R}\mathcal{T}_p(K)} \|v - v_p\|_K^2 + \left( \frac{h_K}{p+1} \|\nabla \cdot v - \Pi_{hp}(\nabla \cdot v)\|_K \right)^2 \right\}, \end{aligned} \quad (3.17) $$

(p-robust global continuous constrained – local discontinuous unconstrained equivalence)

where the hidden constant only depends on the mesh shape-regularity parameter $\kappa_h$ given by (2.1) and the space dimension $d$.

*Proof.* Please first note that the second terms are identical on both sides of (3.17); also recall from [28, Remark 3.4] that they have to be included for the equivalence to hold. Then, since the minimization set on the right-hand side of (3.17) is (seemingly much) bigger than that on the left-hand side, the inequality $\gtrsim$ (actually $\ge$) follows. For the $\lesssim$ inequality, we bound the minimum by employing the projector $\mathbf{P}_{hp}^{\text{div}}(v)$ from Definition 3.3. The commuting property (3.11) and elementwise use of (3.14) from Theorem 3.4 below together with a finite overlap argument following from the mesh shape regularity yield the claim:

$$ \begin{aligned} & \min_{\substack{v_{hp} \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_h) \cap H_{0,N}(\text{div}, \Omega) \\ \nabla \cdot v_{hp} = \Pi_{hp}(\nabla \cdot v)}} \|v - v_{hp}\|^2 + \sum_{K \in \mathcal{T}_h} \left( \frac{h_K}{p+1} \|\nabla \cdot v - \Pi_{hp}(\nabla \cdot v)\|_K \right)^2 \\ & \stackrel{(3.1)}{\le} \sum_{K \in \mathcal{T}_h} \left\{ \|v - \mathbf{P}_{hp}^{\text{div}}(v)\|_K^2 + \left( \frac{h_K}{p+1} \|\nabla \cdot (v - \mathbf{P}_{hp}^{\text{div}}(v))\|_K \right)^2 \right\} \\ & \stackrel{(3.11)}{\lesssim} \sum_{K \in \mathcal{T}_h} \left\{ \sum_{L \in \mathcal{T}_K} \left\{ \min_{v_p \in \mathcal{R}\mathcal{T}_p(L)} \|v - v_p\|_L^2 + \left( \frac{h_L}{p+1} \|\nabla \cdot v - \Pi_{hp}(\nabla \cdot v)\|_L \right)^2 \right\} \right\} \\ & \stackrel{(3.14)}{\lesssim} \sum_{K \in \mathcal{T}_h} \left\{ \min_{v_p \in \mathcal{R}\mathcal{T}_p(K)} \|v - v_p\|_K^2 + \left( \frac{h_K}{p+1} \|\nabla \cdot v - \Pi_{hp}(\nabla \cdot v)\|_K \right)^2 \right\}. \end{aligned} $$

□

### 3.3 Optimal local hp approximation estimates under minimal elementwise Sobolev regularity in $H_{0,N}(\text{div}, \Omega)$

Finally, we show how Definition 3.3 and Theorem 3.4 yield optimal local hp approximation estimates under minimal elementwise Sobolev regularity in $H_{0,N}(\text{div}, \Omega)$.

For any element $K \in \mathcal{T}_h$, let $H^{s_K}(K)$, $s_K \ge 0$, denote the space of vector-valued fields in $L^2(K)$ with each component in $H^{s_K}(K)$. We now focus on functions with additional regularity $H^{s_K}(K)$ requested locally on each mesh element. Moreover, we consider the divergence separately: piecewise polynomial (for simplicity of exposition) first, and then of $H^{t_K}(K)$ regularity, $t_K \ge 0$. Here, the Sobolev regularity exponents $s_K$ and $t_K$ can be different for different mesh elements $K \in \mathcal{T}_h$ and also arbitrarily close, and possibly equal to, 0. The following theorem is a fully h- and p- (mesh-size- and polynomial-degree-) optimal approximation estimate. It improves [28, Theorem 3.6] where the Sobolev regularity exponent $s_K$ can also be arbitrarily close (and possibly equal to) 0 but where it is constant, $s = s_K$ for all mesh elements $K \in \mathcal{T}_h$ and where less attention has been paid to the divergence. Theorem 3.6 can be directly used in a priori error analysis of numerical methods for partial differential equations related to the $H(\text{div}, \Omega)$ space; some examples for (least-squares) mixed finite element methods are given in [28, Section 6]. Locally varying polynomial degree can be addressed as in [50, Theorem 3.4] and [51].

**Theorem 3.6** (hp-optimal approximation estimate in $H_{0,N}(\text{div}, \Omega)$ under minimal elementwise Sobolev regularity). Let $v \in H_{0,N}(\text{div}, \Omega)$, a simplicial mesh $\mathcal{T}_h$ of $\Omega$, and a polynomial degree $p \ge 0$ be given. Let Assumptions 3.1 and 3.2 hold and consider the projector $\mathbf{P}_{hp}^{\text{div}}$ of Definition 3.3. For each mesh element $L \in \mathcal{T}_h$, let

$$ v|_L \in H^{s_L}(L) \quad (3.18) $$
for a Sobolev regularity exponent $s_L \ge 0$. We consider two cases.

Case (i) (piecewise polynomial divergence). Let $\nabla \cdot \mathbf{v} \in \mathcal{P}_p(\mathcal{T}_h)$. Then, for each mesh element $K \in \mathcal{T}_h$,

$$ \|\mathbf{v} - \mathbf{P}_{hp}^{\text{div}}(\mathbf{v})\|_K^2 \lesssim \sum_{L \in \mathcal{T}_K} \left( \frac{h_L^{\min(s_L, p+1)}}{(p+1)^{s_L}} \|\mathbf{v}\|_{H^{s_L}(L)} \right)^2. \quad (3.19) $$

(simplified optimal elementwise hp approximation estimate)

Case (ii) (general case). Let, for each mesh element $L \in \mathcal{T}_h$,

$$ (\nabla \cdot \mathbf{v})|_L \in H^{t_L}(L) \quad (3.20) $$

for a Sobolev regularity exponent $t_L \ge 0$. Then, for each mesh element $K \in \mathcal{T}_h$,

$$ \|\mathbf{v} - \mathbf{P}_{hp}^{\text{div}}(\mathbf{v})\|_K^2 + \left( \frac{h_K}{p+1} \|\nabla \cdot (\mathbf{v} - \mathbf{P}_{hp}^{\text{div}}(\mathbf{v}))\|_K \right)^2 \\ \lesssim \sum_{L \in \mathcal{T}_K} \left\{ \left( \frac{h_L^{\min(s_L, p+1)}}{(p+1)^{s_L}} \|\mathbf{v}\|_{H^{s_L}(L)} \right)^2 + \left( \frac{h_L}{p+1} \frac{h_L^{\min(t_L, p+1)}}{(p+1)^{t_L}} \|\nabla \cdot \mathbf{v}\|_{H^{t_L}(L)} \right)^2 \right\}. \quad (3.21) $$

(optimal elementwise hp approximation estimate)

The constants hidden in $\lesssim$ only depend on the local mesh shape-regularity parameters $\max_{L \in \mathcal{T}_K} \kappa_L$ with $\kappa_L$ given by (2.1), the space dimension $d$, the regularity exponents $s_L$, and, for (3.21), the regularity exponents $t_L$.

*Proof.* We use (3.14), observing from (2.7) that for each mesh element $L \in \mathcal{T}_h$, $[\mathcal{P}_p(L)]^d \subset \mathcal{RT}_p(L)$. Thus, well-known hp-approximation bounds, see e.g. [5, Lemma 4.1], imply that

$$ \min_{\mathbf{v}_p \in \mathcal{RT}_p(L)} \|\mathbf{v} - \mathbf{v}_p\|_L \lesssim \frac{h_L^{\min(s_L, p+1)}}{(p+1)^{s_L}} \|\mathbf{v}\|_{H^{s_L}(L)}, \quad (3.22a) $$

$$ \frac{h_L}{p+1} \|\nabla \cdot \mathbf{v} - \Pi_{hp}(\nabla \cdot \mathbf{v})\|_L \lesssim \frac{h_L}{p+1} \frac{h_L^{\min(t_L, p+1)}}{(p+1)^{t_L}} \|\nabla \cdot \mathbf{v}\|_{H^{t_L}(L)}, \quad (3.22b) $$

with the hidden constants only depending on $\kappa_L$, $d$, $s_L$, and $t_L$. Thus (3.21) follows. As for (3.19), it is a simplification of (3.21) where the divergence terms vanish as $\nabla \cdot \mathbf{v} - \Pi_{hp}(\nabla \cdot \mathbf{v}) = 0$ when $\nabla \cdot \mathbf{v} \in \mathcal{P}_p(\mathcal{T}_h)$. $\square$

## 3.4 Numerical illustration

We provide a quick numerical illustration of the projector $\mathbf{P}_{hp}^{\text{div}}$ from Definition 3.3. We start by the following remark:

**Remark 3.7** (Approximation of $\zeta_p^{a,a}$ by additive Schwarz with line search). In order to make the construction of Step III. easily realizable on a computer, we can replace $\zeta_p^{a,a}$ from (3.7b) by an iterative approximation by additive Schwarz with line search. Let $\zeta_p^a$ be given by (3.7a) and suppose it is nonzero. Let $i$ be an iteration index.

1. Set $i = 0$ and $\zeta_p^{a,i} := 0$.

2. For all vertices $\mathbf{b}$ from the extended vertex patch $\tilde{\omega}_a$, $\mathbf{b} \in \tilde{\mathcal{V}}_a$, consider the (small) vertex patches $\omega_b$ and solve

$$ \begin{array}{l} \delta_p^{a,b,i} := \arg \min_{\substack{\mathbf{v}_p \in \mathcal{RT}_p(\mathcal{T}_b \cap \tilde{\mathcal{T}}_a) \cap \mathcal{H}_{0,N,\psi_b}(\text{div}, \omega_b \cap \tilde{\omega}_a) \\ \nabla \cdot \mathbf{v}_p = 0 \\ (\mathbf{v}_p, \mathbf{r}_h)_K = 0 \quad \forall \mathbf{r}_h \in [\mathcal{P}_0(K)]^d, \forall K \in \mathcal{T}_b}} \|\zeta_p^a - \zeta_p^{a,i} - \mathbf{v}_p\|_{\omega_b}. \end{array} $$

This seemingly overconstrained problem is well posed since it satisfies [15, Assumption A.1].

3. Optimize the descent direction $\delta_p^{a,i} := \sum_{\mathbf{b} \in \tilde{\mathcal{V}}_a} \delta_p^{a,b,i}$ by line search: find

$$ \lambda^i := \arg \min_{\lambda \in \mathbb{R}} \|\zeta_p^a - \zeta_p^{a,i} - \lambda \delta_p^{a,i}\|_{\tilde{\omega}_a}. $$
This gives

$$ \lambda^i = \frac{(\zeta_p^a - \zeta_p^{a,i}, \delta_p^{a,i})_{\tilde{\omega}_a}}{\|\delta_p^{a,i}\|_{\tilde{\omega}_a}^2} $$

and the descent $\lambda^i \delta_p^{a,i} = \lambda^i \sum_{b \in \tilde{\mathcal{V}}_a} \delta_p^{a,b,i}$.

4. If the approximate decomposition is sufficiently precise,

$$ \frac{\|\zeta_p^a - \sum_{b \in \tilde{\mathcal{V}}_a} \sum_{j=0}^i \lambda^j \delta_p^{a,b,j}\|_{\tilde{\omega}_a}}{\|\zeta_p^a\|_{\tilde{\omega}_a}} \le \varepsilon, $$

where $\varepsilon$ is the desired relative tolerance, stop. Otherwise update $\zeta_p^{a,i+1} := \zeta_p^{a,i} + \lambda^i \delta_p^{a,i}$, increase $i := i + 1$, and go back to step (2).

Use $\sum_{j=0}^i \lambda^j \delta_p^{a,a,j}$ for $\zeta_p^{a,a}$.

We consider two space dimensions $d = 2$, $\Omega$ a square $(0, 0.5) \times (0, 0.5)$, $\Gamma_N = \emptyset$, and a fixed triangular mesh composed of 18 right-angled triangles ($\Omega$ is divided into $3 \times 3$ identical sub-squares and each of those into 2 triangles). We consider three divergence-free functions $v = (\partial_y w, -\partial_x w)$, where respectively $w(x, y) = (x + y)^8$, $w(x, y) = \sin(x - 0.5) \cos(y - 0.5)$, and $w(x, y) = \ln(\ln(\sqrt{x^2 + y^2}))$, and let the polynomial degree $p$ vary. For computer implementation of Step III. of Definition 3.3, we proceed following Remark 3.7.

The first two functions $v$ are analytical and enable arbitrary regularity exponents $s_L$ in Theorem 3.6, but the last function, constructed following [8, Section 2.5.1], belongs merely to $H(\text{div}, \Omega)$, with the normal trace $v \cdot n$ only in $H^{-1/2}(\partial\Omega)$. We collect the results in Figure 4. In the left column, we report the approximation errors of the elementwise $L^2(\Omega)$-orthogonal projection, $\|v - \tau_{hp}\| = \{\sum_{K \in \mathcal{T}_h} \min_{v_p \in \mathcal{R}\mathcal{T}_p(K)} \|v - v_p\|_K^2\}^{1/2}$, the approximation error $\|v - \sigma_{hp}\|$ after Step II. of Definition 3.3, and the approximation error $\|v - P_{hp}^{\text{div}}(v)\|$ after the final Step IV. of Definition 3.3. We observe very close results. This is confirmed in the right column of Figure 4, where we plot the ratio $\|v - \sigma_{hp}\|/\|v - \tau_{hp}\|$ and, namely, the ratio $\|v - P_{hp}^{\text{div}}(v)\|/\|v - \tau_{hp}\|$. By the approximation property (3.14), the later is theoretically proven to be independent of the polynomial degree $p$, which is confirmed in all three cases. Numerically, already the intermediate commuting projector $\sigma_{hp}$ obtained in Step II. of Definition 3.3 seems to be $p$-robust.

# 4 Correctness of Definition 3.3 of the projector $P_{hp}^{\text{div}}$

We justify here all steps of Definition 3.3 and summarize the properties of the intermediate objects.

## 4.1 Step I. (construction and properties of the discontinuous projection $\tau_{hp}$)

We start with:

**Lemma 4.1** (Definition (3.2) and property (3.3)). For each mesh element $K \in \mathcal{T}_h$, problem (3.2) for $\tau_{hp}|_K$ is well posed. Moreover, (3.3) holds.

*Proof.* Existence and uniqueness of (3.2) are standard. Note that (3.2) is equivalently stated by the Euler-Lagrange conditions: find $\tau_{hp}|_K \in \mathcal{R}\mathcal{T}_p(K)$ such that

$$ (\tau_{hp} - v, v_p)_K = 0 \quad \forall v_p \in \mathcal{R}\mathcal{T}_p(K). \qquad (4.1) $$

As for (3.3), it follows by definition. $\square$

## 4.2 Step II. (construction and properties of the standard (seemingly over-constrained if $p \ge 1$) equilibration $\sigma_{hp}$)

Let us next address:

**Lemma 4.2** (Definition (3.4) and properties (3.5) and (3.6)). For each mesh vertex $a \in \mathcal{V}_h$, problem (3.4a) for $\sigma_p^a$ is well posed. Moreover, defining $\sigma_{hp}$ by (3.4b), (3.5) and (3.6) hold.
<figure><img src="image_5.png" /><figcaption>Figure 4: Approximation errors $||v - \tau_{hp}||$, $||v - \sigma_{hp}||$, and $||v - P_{hp}^{\text{div}}(v)||$ (left) and equivalence ratios $||v - \sigma_{hp}||/||v - \tau_{hp}||$ and $||v - P_{hp}^{\text{div}}(v)||/||v - \tau_{hp}||$ (right), $d = 2$, $v = (\partial_y w, -\partial_x w)^t$, $\nabla \cdot v = 0$</figcaption></figure>

**Proof.** Problem (3.4a) is in a conventional form from, e.g., [8, 24] for $p = 0$; then, existence and uniqueness of $\sigma_p^a$ follow when the Neumann compatibility condition holds if the normal flux is prescribed all along $\partial\omega_a$, i.e., for interior vertices $a$ and for boundary vertices $a$ such that all faces sharing $a$ lie in $\bar{\Gamma}_N$ ($a \notin \bar{\Gamma}_D$). Taking into account definition (2.5) of $H_{0,N,\psi_a}(\text{div}, \omega_a)$, this is satisfied as $(\Pi_{hp}(\psi^a \nabla \cdot v + \nabla \psi^a \cdot v), 1)_{\omega_a} = (\nabla \cdot (\psi^a v), 1)_{\omega_a} = \langle (\psi^a v) \cdot n, 1 \rangle_{\partial\omega_a} = 0$ when $a \notin \bar{\Gamma}_D$.

When $p \ge 1$, however, (3.4a) features an additional orthogonality constraint. For $d = 3$ (the $d = 2$ case
is easier) it, though, exactly fits the framework of [15, Appendix A] with $q' = q = p$, $g^a = \psi^a \nabla \cdot \mathbf{v} + \nabla \psi^a \cdot \mathbf{v}$, and $\tau_h^a = I_{h,p}^{RT}(\psi^a \tau_{hp})$; actually, in [15, Appendix A], there should be $q' = \max\{q, 1\}$ in place of $q' = \min\{q, 1\}$. Let us check [15, Assumption A.1]. Observe that $g^a \in L^2(\omega_a)$, $\tau_h^a \in \mathcal{RT}_p(\mathcal{T}_a)$, and $(g^a, 1)_{\omega_a} = (\nabla \cdot (\psi^a \mathbf{v}), 1)_{\omega_a} = \langle (\psi^a \mathbf{v}) \cdot \mathbf{n}, 1 \rangle_{\partial \omega_a} = 0$ as above when $\mathbf{a} \notin \bar{\Gamma}_D$. Moreover, let $H_1^*(\omega_a)$ be the subspace of $H^1(\omega_a)$ with mean value zero (when $\mathbf{a} \notin \bar{\Gamma}_D$) or the subspace of $H^1(\omega_a)$ with trace zero on $(\partial \omega_a \cap \bar{\Gamma}_D)^\circ$ (when $\mathbf{a} \in \bar{\Gamma}_D$). Let $q_h \in \mathcal{P}_1(\mathcal{T}_a) \cap H_1^*(\omega_a)$. Then

$$ (\tau_h^a, \nabla q_h)_{\omega_a} + (g^a, q_h)_{\omega_a} = (I_{h,p}^{RT}(\psi^a \tau_{hp}), \underbrace{\nabla q_h}_{|K \in [\mathcal{P}_0(K)]^d \forall K \in \mathcal{T}_a})_{\omega_a} + (\nabla \cdot (\psi^a \mathbf{v}), q_h)_{\omega_a} \\ \stackrel{(2.9b)}{=} (\psi^a \tau_{hp}, \nabla q_h)_{\omega_a} - (\psi^a \mathbf{v}, \nabla q_h)_{\omega_a} \\ = \sum_{K \in \mathcal{T}_a} (\tau_{hp} - \mathbf{v}, \underbrace{\psi^a \nabla q_h}_{|K \in [\mathcal{P}_1(K)]^d})_K \\ \stackrel{(4.1)}{=} 0. $$

We have in particular used the assumption $p \ge 1$, (2.7), (2.9b) with $r_h = (\nabla q_h)|_K \in [\mathcal{P}_0(K)]^d$, and (4.1) with $\mathbf{v}_p = (\psi^a \nabla q_h)|_K \in [\mathcal{P}_1(K)]^d \subset \mathcal{RT}_p(K)$. Existence and uniqueness of $\sigma_p^a$ thus follow from [15, Theorem A.2].

For (3.5), $\sigma_{hp} \in \mathcal{RT}_p(\mathcal{T}_h) \cap \mathcal{H}_{0,N}(\text{div}, \Omega)$ follows by (3.4b) and the definitions in Section 2.5. As for the divergence constraint, as in [9, 28], definition (3.4b), the linearity of the weak divergence, the divergence constraints in (3.4a), the linearity of the elementwise $L^2(\Omega)$-orthogonal projector (2.8), and the partition of unity (2.2) give

$$ \begin{aligned} \nabla \cdot \sigma_{hp} &\stackrel{(3.4b)}{=} \sum_{\mathbf{a} \in \mathcal{V}_h} \nabla \cdot \sigma_p^a &\stackrel{(3.4a)}{=} \sum_{\mathbf{a} \in \mathcal{V}_h} \Pi_{hp}(\psi^a \nabla \cdot \mathbf{v} + \nabla \psi^a \cdot \mathbf{v}) \\ &\stackrel{(2.8)}{=} \Pi_{hp} \left( \sum_{\mathbf{a} \in \mathcal{V}_h} (\psi^a \nabla \cdot \mathbf{v} + \nabla \psi^a \cdot \mathbf{v}) \right) &\stackrel{(2.2)}{=} \Pi_{hp}(\nabla \cdot \mathbf{v}). \end{aligned} $$

Finally, property (3.6a) is immediate from (3.3) and (3.5). As for (3.6b), let $K \in \mathcal{T}_h$ and $r_h \in [\mathcal{P}_0(K)]^d$ be fixed. From the orthogonality constraint in (3.4a) imposed if $p \ge 1$, we have, for any vertex $\mathbf{a} \in \mathcal{V}_K$ of the simplex $K$,

$$ (\sigma_p^a - I_{h,p}^{RT}(\psi^a \tau_{hp}), r_h)_K = 0. $$

Thus, summing over all $\mathbf{a} \in \mathcal{V}_K$ and using the linearity of the canonical projector $I_{h,p}^{RT}$ from (2.9), the partition of unity (2.2), the fact that $I_{h,p}^{RT}(\tau_{hp}) = \tau_{hp}$, and definition (3.4b), we have

$$ 0 = \sum_{\mathbf{a} \in \mathcal{V}_K} (\sigma_p^a - I_{h,p}^{RT}(\psi^a \tau_{hp}), r_h)_K = (\sigma_{hp} - \tau_{hp}, r_h)_K $$

which is the claim (3.6b).

□

## 4.3 Step III. (construction and properties of the p-robust correction $\zeta_{hp}$)

We continue with:

**Lemma 4.3 (Definition (3.7a)).** For each mesh vertex $\mathbf{a} \in \mathcal{V}_h$, problem (3.7a) for $\zeta_p^a$ is well posed.

*Proof.* Problem (3.7a) is again not in a conventional form from, e.g., [8, 24], because of the additional orthogonality constraint. The situation is, though, much easier than for (3.4a) in the proof of Lemma 4.2. Indeed, the minimization (3.7a) is convex and the minimization set not empty, since the zero vector is trivially contained; this comes from the data already satisfying $(\tau_{hp} - \sigma_{hp}, r_h)_K = 0$ for all $r_h \in [\mathcal{P}_0(K)]^d$ and for all $K \in \tilde{\mathcal{T}}_a$. □

**Lemma 4.4 (Decomposition (3.7b)).** For each mesh vertex $\mathbf{a} \in \mathcal{V}_h$, the decomposition (3.7b) is well defined.
**Proof.** By definition from (3.7a), $\zeta_p^a$ lies in $\mathcal{RT}_p(\tilde{T}_a) \cap H_{0,N}(\mathrm{div}, \tilde{\omega}_a)$, is divergence-free, and satisfies $(\zeta_p^a, r_h)_K = 0$ for all $r_h \in [\mathcal{P}_0(K)]^d$ and for all $K \in \tilde{T}_a$. Thus assumption (A.1) below is satisfied with $T_\omega = \tilde{T}_a$ and $\omega = \tilde{\omega}_a$. Then (3.7b) follows immediately from (A.2) (with $V_\omega = \tilde{V}_a$) and (A.4) in Theorem A.1 or Assumption 3.2. Note that we only employ the “interior” component $\zeta_p^{a,a}$; this is from (A.3) supported on the vertex patch subdomain $\omega_a \cap \tilde{\omega}_a$ which is simply $\omega_a$ (no patch truncation happens for the “interior” component, see Figure 3 (right)). $\square$

**Lemma 4.5 (Property (3.8)).** Property (3.8) holds true.

**Proof.** The inclusion $\zeta_{hp} \in \mathcal{RT}_p(\tilde{T}_h) \cap H_{0,N}(\mathrm{div}, \Omega)$ follows immediately by (3.7c) and the definitions in Section 2.5. Note that it is crucial that the components $\zeta_p^{a,a}$ have from (3.7b) zero normal trace on those faces in $\partial\omega_a$ where the hat function $\psi^a$ vanishes or which lie in the Neumann boundary $\bar{\Gamma}_N$. The divergence-free property is evident since all the contributions are divergence-free. $\square$

## 4.4 Step IV. (combining the previous steps)

We finish by:

**Lemma 4.6 (Property (3.10)).** Property (3.10) holds true.

**Proof.** This is an immediate consequence of the definition (3.9) and the property (3.5) together with (3.8) if $p \ge 1$. $\square$

# 5 Proof of Theorem 3.4 on properties of the projector $P_{hp}^{\mathrm{div}}$

Let the assumptions of Theorem 3.4 be satisfied. We prove the claims separately.

## 5.1 Commuting and projection

**Lemma 5.1 (Commuting property (3.11)).** The commuting property (3.11) holds true.

**Proof.** This has been already established in Lemma 4.6. $\square$

**Lemma 5.2 (Projection property (3.12)).** The projection property (3.12) holds true.

**Proof.** Let $v \in \mathcal{RT}_p(\tilde{T}_h) \cap H_{0,N}(\mathrm{div}, \Omega)$. Then clearly $\tau_{hp}$ from (3.2) satisfies $\tau_{hp} = v$. Next, from (3.4a), we see that $\sigma_p^a = I_{h,p}^{\mathcal{RT}}(\psi^a \tau_{hp})$. Indeed, $I_{h,p}^{\mathcal{RT}}(\psi^a \tau_{hp}) \in \mathcal{RT}_p(\tilde{T}_a) \cap H_{0,N,\psi^a}(\mathrm{div}, \omega_a)$ by (2.9a); it is crucial that $\tau_{hp} = v$ is normal-trace continuous here. Moreover,

$$ \nabla \cdot I_{h,p}^{\mathcal{RT}} (\psi^a \tau_{hp}) \stackrel{(2.10)}{=} \Pi_{hp} (\nabla \cdot (\psi^a \tau_{hp})) = \Pi_{hp} (\psi^a \nabla \cdot v + \nabla \psi^a \cdot v) $$

by the commuting property (2.10). Consequently, (3.4b) and the linearity of $I_{h,p}^{\mathcal{RT}}$ as well as its projection property give

$$ \sigma_{hp} = \sum_{a \in \mathcal{V}_h} \sigma_p^a = \sum_{a \in \mathcal{V}_h} I_{h,p}^{\mathcal{RT}} (\psi^a \tau_{hp}) = I_{h,p}^{\mathcal{RT}} \left( \sum_{a \in \mathcal{V}_h} \psi^a \tau_{hp} \right) = I_{h,p}^{\mathcal{RT}} (\tau_{hp}) = \tau_{hp}. $$

Thus, also $\sigma_{hp} = v$. Finally, as Step III. of Definition 3.3 only builds on $\tau_{hp} - \sigma_{hp}$ if $p \ge 1$, all $\zeta_p^a, \zeta_p^{a,a}$, and $\zeta_{hp}$ are zero, whereas $\zeta_{hp} = 0$ by definition if $p = 0$. Then, from (3.9), $P_{hp}^{\mathrm{div}}(v) = \sigma_{hp} + \zeta_{hp} = \sigma_{hp} = v$. $\square$

## 5.2 Approximation

**Lemma 5.3 (Approximation property (3.14)).** The approximation property (3.14) holds true.

**Proof.** The case $p = 0$ is treated as in [28, proof of the approximation property (3.6)]; a $p$-dependent constant is harmless in this lowest-order case. We thus henceforth only consider the case $p \ge 1$. In view of (3.11), the second terms in (3.14) are identical. We thus only have to estimate $\|v - P_{hp}^{\mathrm{div}}(v)\|_K$. Let $K \in \tilde{T}_h$ be fixed and recall the notation from Sections 2.2 and 2.5. We proceed in several steps.
(i) Like in (3.7a), but on the extended element patch $\tilde{T}_K$ in place of the extended vertex patch $\tilde{T}_a$, see Section 2.2 and Figure 2, define

$$ \zeta_p^K := \arg \min_{\substack{v_p \in \mathcal{R}\mathcal{T}_p(\tilde{T}_K) \cap \mathcal{H}_{0,N}(\text{div}, \tilde{\omega}_K) \\ \nabla \cdot v_p = 0 \\ (v_p, r_h)_K = (\tau_{hp} - \sigma_{hp}, r_h)_K = 0 \quad \forall r_h \in [\mathcal{P}_0(K)]^d, \forall K \in \tilde{T}_K}} \|\tau_{hp} - \sigma_{hp} - v_p\|_{\tilde{\omega}_K}. \quad (5.1) $$

This problem is trivially well posed as in Lemma 4.3. Now, as in (3.7b), decompose using Theorem A.1 or Assumption 3.2

$$ \zeta_p^K = \sum_{b \in \tilde{V}_K} \zeta_p^{K,b} \quad \text{with} \quad \zeta_p^{K,b} \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_b \cap \tilde{T}_K) \cap \mathcal{H}_{0,N,\psi^b}(\text{div}, \omega_b \cap \tilde{\omega}_K), \nabla \cdot \zeta_p^{K,b} = 0. \quad (5.2) $$

Note that the assumptions (A.1) are satisfied for the choice $\mathcal{T}_\omega = \tilde{T}_K$ and $\mathcal{V}_\omega = \tilde{V}_K$. Now, crucially, as in (3.7b), the contributions for the vertices $a$ of the element $K$, $a \in \mathcal{V}_K$, actually lie in $\mathcal{R}\mathcal{T}_p(\mathcal{T}_a) \cap \mathcal{H}_{0,N,\psi^a}(\text{div}, \omega_a)$ (as $\mathcal{T}_a$ are included in $\tilde{T}_K$, $\mathcal{T}_a \cap \tilde{T}_K = \mathcal{T}_a$ and no patch truncation happens).

(ii) For each vertex $a \in \mathcal{V}_K$, let us also consider $\zeta_p^K$ from (5.1) restricted to the extended vertex patch $\tilde{\omega}_a$ ($\tilde{\omega}_a$ are included in $\tilde{\omega}_K$ by definition). We again decompose $\zeta_p^K|_{\tilde{\omega}_a}$ using Theorem A.1 or Assumption 3.2

$$ \zeta_p^K|_{\tilde{\omega}_a} = \sum_{b \in \tilde{V}_a} \zeta_p^{K,a,b} \quad \text{with} \quad \zeta_p^{K,a,b} \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_b \cap \tilde{T}_a) \cap \mathcal{H}_{0,N,\psi^b}(\text{div}, \omega_b \cap \tilde{\omega}_a), \\ \nabla \cdot \zeta_p^{K,a,b} = 0. \quad (5.3) $$

Assumptions (A.1) are here satisfied for the choice $\mathcal{T}_\omega = \tilde{T}_a$ and $\mathcal{V}_\omega = \tilde{V}_a$. Crucially, from (A.3), as $\zeta_p^K$ and $\zeta_p^K|_{\tilde{\omega}_a}$ are identical on the extended vertex patches $\tilde{\omega}_a$, the $d+1$ contributions $\zeta_p^{K,a}$ from (5.2) for the vertices $a$ of the element $K$ respectively coincide with the $d+1$ contributions $\zeta_p^{K,a,a}$ from (5.3),

$$ \zeta_p^{K,a} = \zeta_p^{K,a,a} \quad \forall a \in \mathcal{V}_K. \quad (5.4) $$

Indeed, by (A.3), these contributions have the vertex patches $\mathcal{T}_a$ as support and the extended vertex patches $\tilde{T}_a$ as dependency regions and, once again, $\zeta_p^K$ and $\zeta_p^K|_{\tilde{\omega}_a}$ coincide on $\tilde{\omega}_a$. The dependency regions being the extended vertex patches $\tilde{T}_a$ are actually the reason for the remainder equilibration (3.7a) and the decomposition (3.7b) to be performed on the extended vertex patches $\tilde{T}_a$; merely the vertex patches $\mathcal{T}_a$ would not be sufficient. From (5.2)-(5.4), we conclude

$$ \zeta_p^K|_{K} \stackrel{(5.2)}{=} \sum_{a \in \mathcal{V}_K} \zeta_p^{K,a} \stackrel{(5.4)}{=} \sum_{a \in \mathcal{V}_K} \zeta_p^{K,a,a}. \quad (5.5) $$

(iii) Recall the definition of $\tau_{hp}$ from (3.2). We estimate by the triangle inequality and employing the definitions (3.9) and (3.7c) together with the equality (5.5),

$$ \begin{aligned} \|v - P_{hp}^{\text{div}}(v)\|_{K} &\le \|v - \tau_{hp}\|_{K} + \|\tau_{hp} - \sigma_{hp} - \zeta_{hp}\|_{K} \\ &\stackrel{(3.7c)}{=} \|v - \tau_{hp}\|_{K} + \left\| \tau_{hp} - \sigma_{hp} - \zeta_{p}^{K} + \sum_{a \in \mathcal{V}_{K}} (\zeta_{p}^{K,a,a} - \zeta_{p}^{a,a}) \right\|_{K} \quad (5.6) \\ &\le \|v - \tau_{hp}\|_{K} + \|\tau_{hp} - \sigma_{hp} - \zeta_{p}^{K}\|_{K} + \sum_{a \in \mathcal{V}_{K}} \|\zeta_{p}^{K,a,a} - \zeta_{p}^{a,a}\|_{\omega_a}. \end{aligned} $$

From (3.2), the first term above already has the target form. For the last term, we crucially use the linearity of the decomposition (A.3) and its $p$-robust stability (A.5). This gives, for a vertex $a \in \mathcal{V}_K$, recalling (5.3) and (3.7b),

$$ \begin{aligned} \|\zeta_p^{K,a,a} - \zeta_p^{a,a}\|_{\omega_a} &\stackrel{(A.5)}{\lesssim} \|\zeta_p^K - \zeta_p^a\|_{\tilde{\omega}_a} \\ &\le \|\tau_{hp} - \sigma_{hp} - \zeta_p^K\|_{\tilde{\omega}_a} + \|\tau_{hp} - \sigma_{hp} - \zeta_p^a\|_{\tilde{\omega}_a} \quad (5.7) \\ &\le \|\tau_{hp} - \sigma_{hp} - \zeta_p^K\|_{\tilde{\omega}_K} + \|\tau_{hp} - \sigma_{hp} - \zeta_p^a\|_{\tilde{\omega}_a}, \end{aligned} $$
where we have followed by adding and subtracting $\tau_{hp} - \sigma_{hp}$, using the triangle inequality, and extending the integration region. We are thus left estimating $\|\tau_{hp} - \sigma_{hp} - \zeta_p^K\|_{\tilde{\omega}_K}$ for $\zeta_p^K$ from (5.1) and $\|\tau_{hp} - \sigma_{hp} - \zeta_p^a\|_{\tilde{\omega}_a}$ for $\zeta_p^a$ from (3.7a). These take the same form, so that we only show the details for the former.

(iv) Let us thus consider (5.1). Such problems (recall that $\tau_{hp}$ from (3.2) merely belongs to $RT_p(\mathcal{T}_h)$ but not to $H_{0,N}(\text{div}, \Omega)$) have recently been analyzed and p-robust stability has been shown in Braess et al. [9] (for $d=2$) and in [33] (for $d=3$) on: 1) vertex patch subdomains $\omega_a$; 2) with no-flux conditions on $\partial\omega_a$; and 3) without the additional orthogonality constraint. The additional orthogonality constraint has recently been analyzed in [15, Appendix A]. We extend these results to the present setting in Appendices C and D and employ them now here.

(v) Let us first treat the additional orthogonality constraint. Taking $\mathcal{T}_\omega = \tilde{\mathcal{T}}_K$, $r_{hp} = 0$, and $\tau_{hp} = \tau_{hp} - \sigma_{hp}$, we see that (D.6) is trivially satisfied, using in particular (3.6b). Thus, Lemma D.3 yields

$$ \begin{align} & \frac{\|\tau_{hp} - \sigma_{hp} - \zeta_p^K\|_{\tilde{\omega}_K}}{\min_{\substack{v_p \in \mathcal{RT}_p(\tilde{\mathcal{T}}_K) \cap H_{0,N}(\text{div}, \tilde{\omega}_K) \\ \nabla \cdot v_p = 0}} \|\tau_{hp} - \sigma_{hp} - v_p\|_{\tilde{\omega}_K}} \tag{5.8} \\ & \stackrel{(5.1)}{=} \frac{(v_p, r_h)_K = (\tau_{hp} - \sigma_{hp}, r_h)_K = 0}{\forall r_h \in [P_0(K)]^d, \forall K \in \tilde{\mathcal{T}}_K} \nonumber \\ & \stackrel{(D.7)}{\lesssim} \frac{\min_{\substack{v_p \in \mathcal{RT}_p(\tilde{\mathcal{T}}_K) \cap H_{0,N}(\text{div}, \tilde{\omega}_K) \\ \nabla \cdot v_p = 0}} \|\tau_{hp} - \sigma_{hp} - v_p\|_{\tilde{\omega}_K}}{\min_{\substack{v_p \in \mathcal{RT}_p(\tilde{\mathcal{T}}_K) \cap H_{0,N}(\text{div}, \tilde{\omega}_K) \\ \nabla \cdot v_p = \Pi_{hp}(\nabla \cdot v)}} \|\tau_{hp} - v_p\|_{\tilde{\omega}_K}} \tag{5.9} \end{align} $$

Indeed, this follows by the shift by $\sigma_{hp}|_{\tilde{\omega}_K}$ since, by (3.5), it lies in $RT_p(\tilde{\mathcal{T}}_K) \cap H_{0,N}(\text{div}, \tilde{\omega}_K)$; the normal-trace continuity of $\sigma_{hp}$ together with $\sigma_{hp} \cdot n_{\tilde{\omega}_K} = 0$ on $\Gamma_N$ are crucial here. In this important conceptual step, the non p-robust usual equilibration $\sigma_{hp}$ is played out.

(vii) We now finally apply Lemma C.2 with $\mathcal{T}_\omega = \tilde{\mathcal{T}}_K$, $r_{hp} = \Pi_{hp}(\nabla \cdot v)$, and $\tau_{hp} = \tau_{hp}$ to deduce that

$$ \min_{\substack{v_p \in \mathcal{RT}_p(\tilde{\mathcal{T}}_K) \cap H_{0,N}(\text{div}, \tilde{\omega}_K) \\ \nabla \cdot v_p = \Pi_{hp}(\nabla \cdot v)}} \|\tau_{hp} - v_p\|_{\tilde{\omega}_K} \lesssim \min_{\substack{w \in H_{0,N}(\text{div}, \tilde{\omega}_K) \\ \nabla \cdot w = \Pi_{hp}(\nabla \cdot v)}} \|\tau_{hp} - w\|_{\tilde{\omega}_K}. \tag{5.10} $$

This is the crucial p-robust stability bound which makes the power of the infinite-dimensional level of $H_{0,N}(\text{div}, \tilde{\omega}_K)$ appear.

(viii) Let temporarily $v \in H_{0,N}(\text{div}, \Omega)$ from the announcement of Theorem 3.4 have a piecewise polynomial divergence, $\nabla \cdot v \in \mathcal{P}_p(\mathcal{T}_h)$. Then $v$ lies in the minimization set on the right-hand side of (5.10), $v|_{\tilde{\omega}_K} \in H_{0,N}(\text{div}, \tilde{\omega}_K)$ with $\nabla \cdot v = \Pi_{hp}(\nabla \cdot v)$, so that

$$ \min_{\substack{w \in H_{0,N}(\text{div}, \tilde{\omega}_K) \\ \nabla \cdot w = \Pi_{hp}(\nabla \cdot v)}} \|\tau_{hp} - w\|_{\tilde{\omega}_K} \le \|\tau_{hp} - v\|_{\tilde{\omega}_K}. \tag{5.11} $$

In the general case, we need to treat the divergence misfit $\nabla \cdot v - \Pi_{hp}(\nabla \cdot v)$. We proceed as in, e.g., [15, Lemma A.3]. First, we employ the primal-dual equivalence, yielding

$$ \begin{align*} \min_{\substack{w \in H_{0,N}(\text{div}, \tilde{\omega}_K) \\ \nabla \cdot w = \Pi_{hp}(\nabla \cdot v)}} \|\tau_{hp} - w\|_{\tilde{\omega}_K} &= \max_{\substack{v \in H_{0,D}^1(\tilde{\omega}_K) \\ \|\nabla v\|_{\tilde{\omega}_K} = 1}} \{ (\tau_{hp}, \nabla v)_{\tilde{\omega}_K} + (\Pi_{hp}(\nabla \cdot v), v)_{\tilde{\omega}_K} \}, \\ \min_{\substack{w \in H_{0,N}(\text{div}, \tilde{\omega}_K) \\ \nabla \cdot w = \nabla \cdot v}} \|\tau_{hp} - w\|_{\tilde{\omega}_K} &= \max_{\substack{v \in H_{0,D}^1(\tilde{\omega}_K) \\ \|\nabla v\|_{\tilde{\omega}_K} = 1}} \{ (\tau_{hp}, \nabla v)_{\tilde{\omega}_K} + (\nabla \cdot v, v)_{\tilde{\omega}_K} \}. \end{align*} $$

Here, $H_{0,D}^1(\tilde{\omega}_K)$ is the subspace of $H^1(\tilde{\omega}_K)$ with vanishing trace on $(\partial \tilde{\omega}_K \cap \Gamma_D)^\circ$; recall (2.3) and that in Assumption 3.1, we suppose $\partial \tilde{\omega}_K \ne \partial \tilde{\omega}_K \cap \Gamma_N$. Thus, to estimate the right-hand side of (5.10) as in (5.11), we need to bound

$$ \max_{\substack{v \in H_{0,D}^1(\tilde{\omega}_K) \\ \|\nabla v\|_{\tilde{\omega}_K} = 1}} (\nabla \cdot v - \Pi_{hp}(\nabla \cdot v), v)_{\tilde{\omega}_K} = \max_{\substack{v \in H_{0,D}^1(\tilde{\omega}_K) \\ \|\nabla v\|_{\tilde{\omega}_K} = 1}} (\nabla \cdot v - \Pi_{hp}(\nabla \cdot v), v - \Pi_{hp}(v))_{\tilde{\omega}_K}. $$
This is achieved using the hp Poincaré inequality

$$ \|v - \Pi_{hp}(v)\|_{L} \lesssim \frac{h_{L}}{p+1} \|\nabla v\|_{L} $$

for all $L \in \tilde{T}_K$. Altogether, we obtain

$$ \min_{\substack{\mathbf{w} \in \mathcal{H}_{0,N}(\mathrm{div}, \tilde{\omega}_K) \\ \nabla \cdot \mathbf{w} = \Pi_{hp}(\nabla \cdot \mathbf{v})}} \|\boldsymbol{\tau}_{hp} - \mathbf{w}\|_{\tilde{\omega}_K} \lesssim \|\boldsymbol{\tau}_{hp} - \mathbf{v}\|_{\tilde{\omega}_K} + \left\{ \sum_{L \in \tilde{T}_K} \left( \frac{h_L}{p+1} \|\nabla \cdot \mathbf{v} - \Pi_{hp}(\nabla \cdot \mathbf{v})\|_L \right)^2 \right\}^{1/2} . \quad (5.12) $$

Combining the above bounds (5.6)–(5.12) gives the assertion (3.14). □

## 5.3 Stability

**Lemma 5.4** (Stability property (3.15)). The stability property (3.15) holds true.

*Proof.* This follows by the triangle inequality from (3.14). Indeed, let $K \in \mathcal{T}_h$ be fixed. Then

$$ \begin{aligned} \| \mathbf{P}_{hp}^{\mathrm{div}}(\mathbf{v}) \|_{K} & \le \| \mathbf{v} \|_{K} + \| \mathbf{v} - \mathbf{P}_{hp}^{\mathrm{div}}(\mathbf{v}) \|_{K} \\ & \lesssim \sum_{L \in \tilde{T}_K} \left\{ \| \mathbf{v} \|_{L}^{2} + \left( \frac{h_L}{p+1} \| \nabla \cdot \mathbf{v} - \Pi_{hp}(\nabla \cdot \mathbf{v}) \|_{L} \right)^2 \right\}^{1/2}, \end{aligned} $$

where we have also used the trivial $L^2(K)$-orthogonal projection stability

$$ \min_{\mathbf{v}_p \in \mathcal{R}\mathcal{T}_p(L)} \|\mathbf{v} - \mathbf{v}_p\|_L \le \|\mathbf{v}\|_L. $$

**Lemma 5.5** (Stability property (3.16)). The stability property (3.16) holds true.

*Proof.* This is trivial from (3.15), the bound $h_L/p + 1 \le h_\Omega$ (or $h_L/p + 1 \le h_L$), and

$$ \begin{aligned} \| \nabla \cdot \mathbf{P}_{hp}^{\mathrm{div}}(\mathbf{v}) \|_{K}^{(3.11)} &\stackrel{(3.11)}{=} \| \Pi_{hp}(\nabla \cdot \mathbf{v}) \|_{K}^{(2.8)} \le \| \nabla \cdot \mathbf{v} \|_{K}, \\ \| \nabla \cdot \mathbf{v} - \Pi_{hp}(\nabla \cdot \mathbf{v}) \|_{L}^{(2.8)} &\le \| \nabla \cdot \mathbf{v} \|_{L}. \end{aligned} $$

## A A p-stable $\mathcal{R}\mathcal{T}_p \cap H(\mathrm{div})$ decomposition on patch subdomains in two space dimensions

We now state a p-stable decomposition result which follows from Schöberl et al. [47, Section 3]. We consider two-dimensional subdomains $\omega \subset \Omega$ and the corresponding meshes $\mathcal{T}_\omega$; in our applications of Theorem A.1, $\mathcal{T}_\omega$ will be either the extended vertex patch $\tilde{T}_a$ with the corresponding subdomain $\tilde{\omega}_a$, or the extended element patch $\tilde{T}_K$ with $\tilde{\omega}_K$, see Section 2.2. Recall the notation from Section 2.5. There holds:

**Theorem A.1** (A p-stable $\mathcal{R}\mathcal{T}_p \cap H(\mathrm{div})$ decomposition on two-dimensional patches). Let $d = 2$, a simplicial mesh $\mathcal{T}_h$ of $\Omega$, a polynomial degree $p \ge 1$, and $\omega \subset \mathbb{R}^d$ an open and bounded Lipschitz polygonal or polyhedral subdomain of $\Omega$, such that $\bar{\omega}$ is contractible, corresponding to a submesh (patch) of $\mathcal{T}_h$ denoted by $\mathcal{T}_\omega$, with vertex set $\mathcal{V}_\omega$, be given. Let

$$ \delta_p \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_\omega) \cap \mathcal{H}_{0,N}(\mathrm{div}, \omega) \quad \text{with} \quad \nabla \cdot \delta_p = 0, \quad (A.1a) $$

$$ (\delta_p, r_h)_K = 0 \quad \forall r_h \in [\mathcal{P}_0(K)]^d, \forall K \in \mathcal{T}_\omega \quad (A.1b) $$

be a p-degree divergence-free Raviart-Thomas piecewise polynomial on $\omega$ respecting the zero normal trace condition on $\Gamma_N$ if $\partial\omega$ contains faces from $\bar{\Gamma}_N$ and with elementwise vanishing lowest-order moments. Then there exists a decomposition of $\delta_p$ as

$$ \delta_p = \sum_{b \in \mathcal{V}_\omega} \delta_p^b \quad (A.2) $$
where the contributions

$$
\begin{equation}
\begin{aligned}
\delta_p^b \text{ are supported on the vertex patch subdomains } \omega_b \cap \omega, \text{ linearly} \\
\text{depend on } \delta_p \text{ on the extended vertex patch subdomains } \tilde{\omega}_b \cap \omega,
\end{aligned}
\tag{A.3}
\end{equation}
$$

and satisfy

$$
\delta_p^b \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_b \cap \mathcal{T}_\omega) \cap \mathbf{H}_{0,N,\psi^b}(\mathrm{div}, \omega_b \cap \omega) \quad \text{with} \quad \nabla \cdot \delta_p^b = 0, \qquad (A.4)
$$

i.e., recalling (2.6), are divergence-free and such that $\delta_p^b \cdot n_{\omega_b \cap \omega} = 0$ on those faces in $\partial(\omega_b \cap \omega)$ where the hat function $\psi^b$ vanishes or which lie in the Neumann boundary $\bar{\Gamma}_N$. Moreover, the decomposition is $p$-stable in that

$$
\| \delta_p^b \|_{\omega_b \cap \omega} \lesssim \| \delta_p \|_{\tilde{\omega}_b \cap \omega} \quad \forall \mathbf{b} \in \mathcal{V}_{\omega}, \qquad (A.5)
$$

where the constant hidden in ≲ only depends on the local mesh shape-regularity parameter κT_ω given by
κT_ω := max_{K∈T_ω} κ_K.

Proof. (i) Let $\delta_p$ satisfy (A.1a). In two space dimensions, it follows, since $\bar{\omega}$ is contractible, see, e.g. [8, Corollary 2.3.2], that

$$
\delta_p = R_{\frac{\pi}{2}}(\nabla s_p), \qquad (A.6)
$$

where $s_p \in \mathcal{P}_{p+1}(\mathcal{T}_\omega) \cap H_{0, (\partial\omega \cap \Gamma_N)^c}^1(\omega)$ is a $(p+1)$-degree (Lagrange) piecewise polynomial, respecting the zero trace condition on $\Gamma_N$ if $\partial\omega$ contains faces from $\bar{\Gamma}_N$. Here,

$$
R_{\frac{\pi}{2}} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}
$$

is the matrix of rotation by $\frac{\pi}{2}$. Moreover, using (A.1b), we see, for any triangle $K \in \mathcal{T}_\omega$ and any of its vertices, $\mathbf{a} \in \mathcal{V}_K$, that

$$
\begin{align*}
0 &= (R_{\frac{\pi}{2}}(\nabla s_p), \nabla \psi^a)_K = (\nabla s_p, R_{\frac{\pi}{2}}^t (\nabla \psi^a))_K \stackrel{\text{Green}}{=} \langle s_p, R_{\frac{\pi}{2}}^t (\nabla \psi^a) \cdot \mathbf{n}_K \rangle_{\partial K} \\
  &= \langle s_p, \nabla \psi^a \cdot (R_{\frac{\pi}{2}} \mathbf{n}_K) \rangle_{F_1 \cup F_2} = \frac{\langle s_p, 1 \rangle_{F_2}}{|F_2|} - \frac{\langle s_p, 1 \rangle_{F_1}}{|F_1|}
\end{align*}
$$

for the two faces (edges) $F_1, F_2$ that share the vertex $\mathbf{a}$ (numbered in the counterclockwise orientation in the triangle $K$, starting from the vertex $\mathbf{a}$). This means that all mean values of $s_p$ on all faces contained in $\mathcal{T}_\omega$ coincide. Thus, to fix $s_p$ from (A.6) completely when $\partial\omega$ contains no face from $\bar{\Gamma}_N$ (not just its (rotated) gradient), we can set its mean value on any face in $\mathcal{T}_\omega$ to zero, and $s_p$ is independent of which face we have chosen, since then all its mean values on all faces are zero,

$$
\frac{\langle s_p, 1 \rangle_F}{|F|} = 0 \quad \text{for all faces } F \text{ of } \mathcal{T}_\omega. \tag{A.7}
$$

(ii) For the above continuous piecewise polynomial $s_p$, consider the decomposition of Schöberl et al. [47, Section 3]. First, let's choose the “coarse grid contribution” ($u_0$ in [47, equation (2)]) as zero. This is eligible in terms of [47, Lemma 3.1], since

$$
\| \nabla 0 \|_{\omega} \le \| \nabla s_p \|_{\omega}, \qquad (A.8a)
$$

$$
\| \nabla s_p \|_{\omega} = \| \nabla s_p \|_{\omega}, \qquad (A.8b)
$$

$$
\|h^{-1}s_p\|_{\omega}^{2} = \sum_{K \in \mathcal{T}_{\omega}} (h_{K}^{-2} \|s_p\|_{K}^{2}) \leq 6 \sum_{K \in \mathcal{T}_{\omega}} \|\nabla s_p\|_{K}^{2}, \qquad (A.8c)
$$

where we have used the face-mean value Poincaré–Friedrichs inequality, see [49, Lemma 4.1] for the
value 6 of the constant. Consequently, there is no global low order component. The construction of [47,
Section 3.2–3.4] then gives the decomposition, see equation (11) in this reference (after associating the
face and element contributions with the vertex contributions),

$$
s_p = \sum_{b \in \mathcal{V}_\omega} s_p^b, \qquad (A.9a)
$$
<figure><img src="image_6.png" /><figcaption>Figure 5: Examples where Definition B.1, property (ii), point 1 is not satisfied (left), Definition B.1, property (i), point 1 is not satisfied (middle), and Definition B.1, property (ii), point 1 is not satisfied (right). For the marked vertex **a** and the hatched simplex, the already enumerated simplices sharing **a** are dotted.</figcaption></figure>

where

$$s_p^b \in \mathcal{P}_{p+1}(\mathcal{T}_b \cap \mathcal{T}_\omega), \quad s_p = 0 \text{ on faces in } \partial(\omega_b \cap \omega) \text{ where the hat} \\ \text{function } \psi^b \text{ vanishes or which lie in the Neumann boundary } \bar{\Gamma}_N. \tag{A.9b}$$

Moreover, this decomposition is p-robustly stable in that, see [47, Section 3.4],

$$\sum_{b \in \mathcal{V}_\omega} \| \nabla s_p^b \|_{\omega_b \cap \omega}^2 \lesssim \| \nabla s_p \|_{\omega}^2.$$

The inspection of the developments of [47, Section 3.2–3.4] shows that $s_p^b$ are solely constructed from and linearly depend on the values of $s_p$ on the extended patches $\tilde{\omega}_b \cap \omega$ and satisfy more precisely the local stability bounds

$$\|\nabla s_p^b\|_{\omega_b \cap \omega} \lesssim \|\nabla s_p\|_{\tilde{\omega}_b \cap \omega} \quad \forall \mathbf{b} \in \mathcal{V}_\omega. \tag{A.10}$$

Crucially, the constant hidden in $\lesssim$ above only depends on the shape-regularity parameter $\kappa_{\mathcal{T}_\omega}$ of the mesh $\mathcal{T}_\omega$.

(iii) Now take

$$\delta_p^b := R_{\frac{\pi}{2}}(\nabla s_p^b). \tag{A.11}$$

It follows that $\delta_p^b$ satisfies the first line in (A.3) and (A.4). Crucially, also the second line in (A.3) is satisfied. The dependence region is indeed $\tilde{\omega}_b \cap \omega$, since from the face-wise zero mean value property (A.7) (a consequence of assumption (A.1b)), $s_p|_K$ only depends on $\delta_p|_K$ for all $K \in \mathcal{T}_\omega$. Moreover, (A.2) follows immediately from (A.9a) to which we apply the rotated gradient, (A.6), and (A.11). Finally, (A.5) is a direct consequence of (A.10) since

$$\|\nabla s_p^b\|_{\omega_b \cap \omega} = \|R_{\frac{\pi}{2}}(\nabla s_p^b)\|_{\omega_b \cap \omega} \text{ and } \|\nabla s_p\|_{\tilde{\omega}_b \cap \omega} = \|R_{\frac{\pi}{2}}(\nabla s_p)\|_{\tilde{\omega}_b \cap \omega}$$

together with (A.11) and (A.6). $\square$

# B Suitable enumeration/shellability of patches of mesh elements

Recall the notation from Section 2 and also recall that by “face”, we mean “$(d-1)$-dimensional face”. Let $|S|$ denote the cardinality (number of elements) of the set $S$. The following definition will be central:

**Definition B.1** (Suitable patch enumeration). Let $\mathcal{T}_\omega$ be a simplicial mesh with the corresponding open and bounded polygonal or polyhedral domain $\omega \subset \mathbb{R}^d$, $d = 2, 3$, such that $\bar{\omega}$ is contractible. An enumeration $\{K_1, \dots, K_{|\mathcal{T}_\omega|}\}$ of the simplices in $\mathcal{T}_\omega$ is suitable if:
(i) (Only for $d=3$) For all $1 < i \le |\mathcal{T}_\omega|$, if there are at least 2 faces of $K_i$ shared with previously enumerated simplices, intersecting in an edge $e$, then 1) all the simplices sharing the edge $e$ different from $K_i$ come sooner in the enumeration; 2) the edge $e$ lies in the interior of $\omega$.

(ii) For all $1 < i \le |\mathcal{T}_\omega|$, if there are $d$ faces of $K_i$ shared with previously enumerated simplices, intersecting in a vertex $\mathbf{a}$, then 1) all the simplices sharing the vertex $\mathbf{a}$ different from $K_i$ come sooner in the enumeration; 2) the vertex $\mathbf{a}$ lies in the interior of $\omega$.

(iii) For all $1 < i \le |\mathcal{T}_\omega|$, there are between 1 and $d$ face neighbors of $K_i$ which have been already enumerated and correspondingly, there is at least 1 face neighbor which has not been enumerated yet, or $K_i$ has a face on the boundary $\partial\omega$. In particular, there is no enumerated face neighbor only for $K_1$ and all face neighbors are already enumerated for $K_{|\mathcal{T}_\omega|}$, which moreover has a face on the boundary $\partial\omega$.

Illustrations are provided in Figure 5.

In algebraic topology/discrete geometry, there exists a concept of shellability of simplicial complexes. Let $\omega$ be an open and bounded polygon or polyhedron with $\bar{\omega}$ contractible. Following Ziegler [53, Definition 5.1 and Examples 5.2.(iii)], we define a $d$-dimensional simplicial complex $\mathcal{T}_\omega$ as a nonempty finite set composed of closed $d'$-dimensional simplices $0 \le d' \le d$ set in $\mathbb{R}^d$ and covering $\bar{\omega}$ such that (i) for $K \in \mathcal{T}_\omega$, a $d'$-dimensional simplex, all its $d''$-dimensional faces, $0 \le d'' \le d' - 1$, are in $\mathcal{T}_\omega$; (ii) the intersection $K \cap L$ of $K, L \in \mathcal{T}_\omega$ is either empty or a $d''$-dimensional face, $0 \le d'' \le d' - 1$, of both $K$ and $L$. We only consider the so-called pure complexes where every simplex of dimension $d' < d$ is a $d'$-dimensional face of some simplex $K \in \mathcal{T}_\omega$ of dimension exactly $d$. Thus, the present simplicial meshes are the $d$-simplices of a pure $d$-dimensional simplicial complex set in a subdomain $\omega$ of $\mathbb{R}^d$. Following Ziegler [53, Definition 8.1 and Remarks 8.3.(ii)] or Kozlov [40, Definition 12.1], we define:

**Definition B.2 (Shelling of a simplicial complex).** Let $\omega$ be an open and bounded polygon or polyhedron with $\bar{\omega}$ contractible. A shelling of a simplicial complex $\mathcal{T}_\omega$ with domain $\bar{\omega}$ is an enumeration $K_1...K_{|\mathcal{T}_\omega|}$ of the $d$-dimensional simplices of $\mathcal{T}_\omega$ such that for all $1 < i \le |\mathcal{T}_\omega|$, the intersection of $K_i$ with the previously enumerated $d$-dimensional simplices is a nonempty collection of $(d-1)$-dimensional faces of $K_i$.

Definition B.2 means that the intersection of $K_i$ with the previously enumerated $d$-dimensional simplices cannot be and cannot include mere points (if $d=2$) and cannot be and cannot include mere points or edges (if $d=3$). The illustrations of not suitable enumerations of Figure 5 apply here as well: these are not shellings. Indeed, the triangle $K_3$ only shares a point with the previously enumerated triangles $K_1$ and $K_2$ (left), the tetrahedron $K_2$ only shares an edge with the tetrahedron $K_1$ (middle), and similarly in the right figure, where in addition the tetrahedron $K_3$ only shares mere edges with the tetrahedra $K_1$ and $K_2$. Denote by $\mathcal{T}_{\omega,i}$ the simplicial complex formed by the $d$-dimensional simplices enumerated before $K_i$, $K_1, \dots, K_{i-1}$, $1 \le i \le |\mathcal{T}_\omega|$. As a distinctive feature, shellability gives that the closure of the underlying open domain of $\mathcal{T}_{\omega,i}$ is a triangulated manifold with boundary, more precisely a topological $d$-ball (and in particular contractible), for any $1 \le i \le |\mathcal{T}_\omega|$, see Ziegler [52, proof of Proposition 2.4.(iv)] or Chaumont-Frelet et al. [13, Lemma 7.5].

It turns out that the following crucial result holds true:

**Lemma B.3 (Equivalence of suitable patch enumeration with shellability).** *Suitable enumeration of Definition B.1 is equivalent with shelling of Definition B.2.*

*Proof.* (i) Shellability $\implies$ suitable patch enumeration. By definition, a mere point ($d=2,3$) or edge ($d=3$) connection of $K_i$ to the previously enumerated $d$-simplices is forbidden. This implies Definition B.1, properties (i) and (ii), see Figure 5 for illustration. Next, the requirement that the intersection of $K_i$ with $\mathcal{T}_{\omega,i}$ is a nonempty collection of $(d-1)$-dimensional faces of $K_i$ implies that there are between 1 and $d$ face neighbors of $K_i$ which have been already enumerated. Finally, the fact that $K_{|\mathcal{T}_\omega|}$ has a face on the boundary $\partial\omega$ is a consequence of the fact that the last $d$-dimensional simplex in any shelling has to be free, see, e.g. [52, proof of Proposition 2.4.(iv)].

(ii) Suitable patch enumeration $\implies$ shellability. Since there are between 1 and $d$ face neighbors of $K_i$ which have been already enumerated for $1 < i \le |\mathcal{T}_\omega|$, the intersection of $K_i$ with $\mathcal{T}_{\omega,i}$ contains between 1 and $d$ $(d-1)$-dimensional faces of $K_i$. Definition B.1, properties (i) and (ii) then imply that the previously enumerated $d$-simplices did not have a mere point ($d=2,3$) or edge ($d=3$) connection to the $d$-simplices enumerated before. $\square$
In two space dimensions, following Bing [7] or Ziegler [53, Examples 8.4.(i)], cf. also or Chaumont-Frelet et al. [13, Lemma 7.11], the situation is simple:

**Theorem B.4** (Shellability of simplicial complexes for $d = 2$). Let $\omega$ be an open and bounded polygon in $\mathbb{R}^2$ with $\bar{\omega}$ contractible. All simplicial complexes over $\omega$ are shellable.

In three space dimensions, a typical patch $\mathcal{T}_\omega$ from finite element mesh will also be shellable. However, rigorously, the situation is much more complex:

**Remark B.5** (Shellable simplicial complexes for $d = 3$). There holds:

* Any simplicial complex with at most 8 vertices is shellable, see Lutz [42, Corollary 6]. (It is conjectured that any simplicial complex with at most 17 tetrahedra is shellable).

* There exists a nonshellable simplicial complex with 9 vertices and 18 tetrahedra, see Lutz [42]. A nonshellable simplicial complex with 10 vertices and 21 tetrahedra with a graphical visualization is given in Ziegler [52].

* Any Delaunay simplicial mesh (and more generally any regular simplicial mesh) is shellable, see Ziegler [53, Definition 5.3 and Corollary 8.14].

* Simplicial complexes with many tetrahedra with respect to vertices are shellable (more precisely vertex-decomposable), see [27, Theorem 1.1].

* Algorithms are available to decide whether a simplicial complex is shellable, see Moriyama [44], Cook [20], and the web page https://macaulay2.com/doc/Macaulay2/share/doc/Macaulay2/SimplicialDecomposability/html/index.html.

* In general, it is NP-complete to decide whether a given simplicial complex is shellable, see Goaoc et al. [39] and Paták and Tancer [45].

* For every simplicial mesh, there exists a subdivision that is shellable, see Bruggesser and Mani [10, Proposition 1]. Actually, Adiprasito and Benedetti [1, Theorem A] show that the barycentric refinement (inserting a barycenter to each edge, face, and tetrahedron) is sufficient.

# C $p$-stable broken $\mathcal{H}(\text{div})$ polynomial extensions on patch subdomains

We summarize here our results on $p$-stable broken $\mathcal{H}(\text{div})$ polynomial extensions on patch subdomains.

## C.1 Available results

First $p$-stable $\mathcal{H}(\text{div})$ polynomial extensions on a single triangle or tetrahedron have been achieved in Ainsworth and Demkowicz [2], Demkowicz et al. [26], and Costabel and McIntosh [21], see also the references therein. Let $K$ be a triangle or a tetrahedron and let $p \ge 0$. Let $\tau_K \in \mathcal{RT}_p(K)$ be a volume datum, $r_K \in \mathcal{P}_p(K)$ a target divergence, and $r_F \in \mathcal{P}_p(F)$ a target normal trace; the latter is prescribed on $\mathcal{F}_K^N$, a subset of all $(d-1)$-dimensional faces of $K$, possibly empty or containing some or all faces of $K$. The combination of the above-cited normal trace and divergence liftings allows to prove, see [33, Lemma A.3], that

$$ \begin{array}{lll} \min_{\substack{v_p \in \mathcal{RT}_p(K) \\ \nabla \cdot v_p = r_K}} & \|\tau_K - v_p\|_K \lesssim & \min_{\substack{v \in \mathcal{H}(\text{div}, K) \\ \nabla \cdot v = r_K}} \|\tau_K - v\|_K, \quad (C.1) \\ \text{$v_p \cdot n_K = r_F$ on all $F \in \mathcal{F}_K^N$} & & \text{$v \cdot n_K = r_F$ on all $F \in \mathcal{F}_K^N$} \end{array} $$

where the hidden constant only depends on the shape-regularity parameter $\kappa_K$ of the element $K$ and the space dimension $d$ (the form (C.1) follows from [33, Lemma A.3] by a shift by $\tau_K$). On $\mathcal{H}(\text{div}, K)$, the normal trace condition is understood by duality as in (2.4). When $\mathcal{F}_K^N$ is composed of all faces of $K$, the Neumann compatibility condition

$$ \sum_{F \in \mathcal{F}_K^N} \langle r_F, 1 \rangle_F = (r_K, 1)_K $$

needs to be satisfied.

Later, $p$-stable broken polynomial extension achieved similar results as (C.1) but on patches of elements, where, crucially, the datum $\tau_{hp}$ is a piecewise (broken Raviart-Thomas) polynomial. For vertex
patches $\omega_a$ and prescribed normal trace boundary conditions on $\partial\omega_a$, they have been established in Braess et al. [9] in two space dimensions and in [33, Corollaries 3.3 and 3.8] (see also [16, Proposition 3.1 and Corollary 4.1]) in three space dimensions.

## C.2 Larger patches and no boundary conditions

We now extend the above results in two directions: for larger patches $\omega$ and without prescription of normal trace boundary conditions on $\partial\omega$. In our application on step (vii) of the proof of Lemma 5.3, we employ this result for to $\omega = \tilde{\omega}_K$ and $\omega = \tilde{\omega}_a$ when away from the Neumann boundary $\Gamma_N$; treatment of the boundary case is postponed to Section C.3.

**Theorem C.1** (p-stable broken $H(\text{div})$ polynomial extension on larger patches and without boundary conditions). Let $\mathcal{T}_\omega$ be a simplicial mesh with the corresponding open, bounded, and Lipschitz polygon or polyhedron $\omega \subset \mathbb{R}^d$, $d = 2, 3$, with $\bar{\omega}$ contractible. Let $\tau_{hp} \in \mathcal{RT}_p(\mathcal{T}_\omega)$ and $r_{hp} \in \mathcal{P}_p(\mathcal{T}_\omega)$ be respectively a volume datum, a broken Raviart-Thomas vector-valued piecewise polynomial, and a target divergence, a scalar-valued piecewise polynomial, of degree $p \ge 0$. Then

$$ v_p \in \mathcal{RT}_p(\mathcal{T}_\omega) \cap \mathcal{H}(\text{div}, \omega) \quad \min_{\nabla \cdot v_p = r_{hp}} \|\tau_{hp} - v_p\|_{\omega} \lesssim \min_{\nabla \cdot v = r_{hp}} \|\tau_{hp} - v\|_{\omega}, \quad (C.2) $$

where the constant hidden in $\lesssim$ only depends on the mesh shape-regularity parameter $\kappa_{\mathcal{T}_\omega} := \max_{K \in \mathcal{T}_\omega} \kappa_K$, the ratio $h_\omega / \min_{K \in \mathcal{T}_\omega} h_K$, and the space dimension $d$.

*Proof for d = 2.* The case $d = 2$ can be handled as the first case for $d = 3$ below, since in two space dimensions, thanks to Theorem B.4 and Lemma B.3, a suitable enumeration of $\mathcal{T}_\omega$ as per Definition B.1 always exists. $\square$

*Proof for d = 3 when a suitable enumeration as per Definition B.1 exists.* We will follow [33, Section 6], see also [16, Section 6.4]. Let

$$ v^* := \arg \min_{v \in \mathcal{H}(\text{div}, \omega) \atop \nabla \cdot v = r_{hp}} \|\tau_{hp} - v\|_{\omega} \quad (C.3) $$

denote the infinite-dimensional $H(\text{div}, \omega)$ minimizer of the right-hand side of (C.2). We present a constructive proof of (C.2) which proceeds along the enumeration of Definition B.1. On each element $K_i$, $1 \le i \le |\mathcal{T}_\omega|$, we in particular construct a suitable minimizer $\xi_i \in \mathcal{RT}_p(K_i)$ and we gradually set

$$ \xi_{hp}|_{K_i} := \xi_i. \quad (C.4) $$

We then verify that

$$ \xi_{hp} \in \mathcal{RT}_p(\mathcal{T}_\omega) \cap \mathcal{H}(\text{div}, \omega) \quad \text{with} \quad \nabla \cdot \xi_{hp} = r_{hp} \quad (C.5) $$

and that

$$ \|\tau_{hp} - \xi_{hp}\|_{\omega} \lesssim \|\tau_{hp} - v^*\|_{\omega}, \quad (C.6) $$

which establishes (C.2). More precisely, on each step $1 \le i \le |\mathcal{T}_\omega|$, we will verify that

$$ \|\tau_{hp} - \xi_i\|_{K_i} \lesssim \|\tau_{hp} - v^*\|_{\omega}. \quad (C.7) $$

This yields (C.6) up to a constant depending on the shape-regularity parameter $\kappa_{\mathcal{T}_\omega}$ of the mesh $\mathcal{T}_\omega$, the ratio $h_\omega / \min_{K \in \mathcal{T}_\omega} h_K$, and the space dimension $d$. Moreover, as $\nabla \cdot \xi_i = r_{hp}|_{K_i}$ and since $\xi_i$ will have its normal trace prescribed by $\xi_j$ on the previously enumerated $K_j$, $\xi_{hp}$ will have no normal trace jumps and (C.5) follows. We proceed along the enumeration $1 \le i \le |\mathcal{T}_\omega|$ of Definition B.1 and consider different cases.

(i) On the first element $K_1$, let

$$ \xi_1 := \arg \min_{v_p \in \mathcal{RT}_p(K_1) \atop \nabla \cdot v_p = r_{hp}|_{K_1}} \|\tau_{hp} - v_p\|_{K_1}. \quad (C.8) $$

This is a well-posed problem. Crucially, since the data $\tau_{hp}|_{K_1}$ and $r_{hp}|_{K_1}$ in (C.8) are polynomial, we know from (C.1) that we can pass to the infinite-dimensional level,

$$ \|\tau_{hp} - \xi_1\|_{K_1} \lesssim \min_{v \in \mathcal{H}(\text{div}, K_1) \atop \nabla \cdot v = r_{hp}|_{K_1}} \|\tau_{hp} - v\|_{K_1}. \quad (C.9) $$
Finally, since the infinite-dimensional minimizer $v^*$ from (C.3) restricted to the element $K_1$, $v^*|_{K_1}$, belongs to the minimization set on the right-hand side of (C.9) (please note that there are no normal trace conditions in (C.9)), we obtain

$$ \|\tau_{hp} - \xi_1\|_{K_1} \lesssim \|\tau_{hp} - v^*\|_{K_1}, \qquad (C.10) $$

which immediately gives (C.7) for $i=1$.

(ii) On each element $K_i$ with exactly one face shared with some previously enumerated simplex, say $F_{i,j}$ shared with $K_j$, $j < i$, we consider

$$ \begin{array}{ll} \xi_i := \arg \min_{\substack{v_p \in \mathcal{R}\mathcal{T}_p(K_i) \\ \nabla \cdot v_p = r_{hp}|_{K_i} \\ v_p \cdot n_{K_i} = \xi_{hp}|_{K_j} \cdot n_{K_i} \text{ on } F_{i,j}}} \|\tau_{hp} - v_p\|_{K_i}. & (C.11) \\ \text{min} & \|\tau_{hp} - v_p\|_{K_i}. \end{array} $$

Please note that since $j < i$ and by (C.4), $\xi_{hp}|_{K_j}$ is known. Then (C.11) is well-posed; there is in particular no compatibility condition to verify, since the normal trace is only imposed on one face. We now again employ (C.1). This yields

$$ \|\tau_{hp} - \xi_i\|_{K_i} \lesssim \min_{\substack{v \in \mathcal{H}(\text{div}, K_i) \\ \nabla \cdot v = r_{hp}|_{K_i} \\ v \cdot n_{K_i} = \xi_{hp}|_{K_j} \cdot n_{K_i} \text{ on } F_{i,j}}} \|\tau_{hp} - v\|_{K_i}. \qquad (C.12) $$

Unfortunately, now $v^*|_{K_i}$ does not belong to the minimization set on the right-hand side of (C.12) since there is a normal trace condition on the face $F_{i,j}$ imposed. The fix is, for the moment, easy. Consider the face neighbor $K_j$, the function $v^* - \xi_{hp}$ on $K_j$ (note that it is divergence-free), and map it to $K_i$ by the contravariant Piola transformation (see, e.g., [31, Section 9]) preserving the face $F_{i,j}$, say $\psi$, forming

$$ v := v^*|_{K_i} - \psi^{-1}((v^* - \xi_{hp})|_{K_j}), \qquad (C.13) $$

see [33, equation (6.10)] for the details. This removes the normal trace of $v^*$ and brings instead the requested $\xi_{hp}|_{K_j} \cdot n_{K_i}$ (in appropriate weak sense), so that $v$ from (C.13) now crucially belongs to the minimization set on the right-hand side of (C.12). Consequently, we obtain

$$ \|\tau_{hp} - \xi_i\|_{K_i} \lesssim \|\tau_{hp} - v^*|_{K_i} + \psi^{-1}((v^* - \xi_{hp})|_{K_j})\|_{K_i}. \qquad (C.14) $$

Finally, by the triangle inequality and the properties of the Piola transform (recall that we suppose shape regularity of $\mathcal{T}_\omega$)

$$ \begin{align} & \|\tau_{hp} - v^*|_{K_i} + \psi^{-1}((v^* - \xi_{hp})|_{K_j})\|_{K_i} \nonumber \\ & \le \|\tau_{hp} - v^*\|_{K_i} + \|\psi^{-1}((v^* - \xi_{hp})|_{K_j})\|_{K_i} \nonumber \\ & \lesssim \|\tau_{hp} - v^*\|_{K_i} + \|v^* - \xi_{hp}\|_{K_j} \tag{C.15} \\ & \le \|\tau_{hp} - v^*\|_{K_i} + \|\tau_{hp} - v^*\|_{K_j} + \|\tau_{hp} - \xi_{hp}\|_{K_j} \nonumber \\ & \lesssim \|\tau_{hp} - v^*\|_{\omega}, \nonumber \end{align} $$

where, in the last estimate, we have employed (C.7) in $K_j$, which has been established previously since $j < i$. Thus (C.7) is established.

(iii) On each element $K_i$ with exactly two faces shared with some previously enumerated simplices, say $F_{i,j}$ shared with $K_j$, $j < i$, and $F_{i,k}$ shared with $K_k$, $k < i$, we consider

$$ \begin{array}{ll} \xi_i := \arg \min_{\substack{v_p \in \mathcal{R}\mathcal{T}_p(K_i) \\ \nabla \cdot v_p = r_{hp}|_{K_i} \\ v_p \cdot n_{K_i} = \xi_{hp}|_{K_j} \cdot n_{K_i} \text{ on } F_{i,j} \\ v_p \cdot n_{K_i} = \xi_{hp}|_{K_k} \cdot n_{K_i} \text{ on } F_{i,k}}} \|\tau_{hp} - v_p\|_{K_i}. & (C.16) \\ \text{min} & \|\tau_{hp} - v_p\|_{K_i}. \end{array} $$

Again, since $j < i$ and $k < i$ and by (C.4), $\xi_{hp}|_{K_j}$ and $\xi_{hp}|_{K_k}$ are known. Then (C.16) is well-posed; there is again no compatibility condition to verify, since the normal trace is only imposed on two faces. We then again employ (C.1), which now yields

$$ \begin{align} \|\tau_{hp} - \xi_i\|_{K_i} \lesssim \min_{\substack{v \in \mathcal{H}(\text{div}, K_i) \\ \nabla \cdot v = r_{hp}|_{K_i} \\ v \cdot n_{K_i} = \xi_{hp}|_{K_j} \cdot n_{K_i} \text{ on } F_{i,j} \\ v \cdot n_{K_i} = \xi_{hp}|_{K_k} \cdot n_{K_i} \text{ on } F_{i,k}}} \|\tau_{hp} - v\|_{K_i}. & (C.17) \\ \text{min} & \|\tau_{hp} - v\|_{K_i}. \end{align} $$
As above in step (ii), the continuous-level minimizer $\mathbf{v}^*$ from (C.3) restricted to $K_i$ does not belong to the minimization set on the right-hand side of (C.17) since there are two normal trace conditions on the two faces $F_{i,j}$ and $F_{i,k}$ imposed. Crucially, by property (i) of Definition B.1 on the enumeration, all the simplices sharing the edge $e$ common to the two faces $F_{i,j}$ and $F_{i,k}$ come sooner in the enumeration and the edge $e$ lies in the interior of $\omega$. This enables to construct a suitable $\mathbf{v}$ in this sprit of (C.13) but which now involves Piola mappings from all the simplices sharing the edge $e$ except for $K_i$. This is done in a “2-folding” way which replaces $\mathbf{v}^* \cdot \mathbf{n}_{K_i}$ on $F_{i,j}$ and $F_{i,k}$ (in a proper weak sense) by respectively $\xi_{hp}|_{K_j} \cdot \mathbf{n}_{K_i}$ and $\xi_{hp}|_{K_k} \cdot \mathbf{n}_{K_i}$; the precise formula is [33, equation (6.12)]. Existence of a two-color refinement around edges of [33, Lemma B.2] is crucial at this step. Then (C.7) is established similarly to (C.15).

(iv) Finally, on each element $K_i$ with exactly three faces shared with some previously enumerated simplices, say $F_{i,j}$ shared with $K_j$, $j < i$, $F_{i,k}$ shared with $K_k$, $k < i$, and $F_{i,l}$ shared with $K_l$, $l < i$, we consider

$$ \xi_i := \arg \min_{\substack{\mathbf{v}_p \in \mathcal{R}\mathcal{T}_p(K_i) \\ \nabla \cdot \mathbf{v}_p = r_{hp}|_{K_i}}} \|\tau_{hp} - \mathbf{v}_p\|_{K_i}. \quad (C.18) $$
$$ \begin{aligned} \mathbf{v}_p \cdot \mathbf{n}_{K_i} &= \xi_{hp}|_{K_j} \cdot \mathbf{n}_{K_i} \quad \text{on } F_{i,j} \\ \mathbf{v}_p \cdot \mathbf{n}_{K_i} &= \xi_{hp}|_{K_k} \cdot \mathbf{n}_{K_i} \quad \text{on } F_{i,k} \\ \mathbf{v}_p \cdot \mathbf{n}_{K_i} &= \xi_{hp}|_{K_l} \cdot \mathbf{n}_{K_i} \quad \text{on } F_{i,l} \end{aligned} $$

Again, all $\xi_{hp}|_{K_j}$, $\xi_{hp}|_{K_k}$, and $\xi_{hp}|_{K_l}$ are known at this stage. Then (C.18) is well-posed; there is still no compatibility condition to verify, since the normal trace is only imposed on three of the four faces of $K_i$. Employing once more (C.1), we have

$$ \|\tau_{hp} - \xi_i\|_{K_i} \lesssim \min_{\substack{\mathbf{v} \in \mathcal{H}(\text{div}, K_i) \\ \nabla \cdot \mathbf{v} = r_{hp}|_{K_i}}} \|\tau_{hp} - \mathbf{v}\|_{K_i}. \quad (C.19) $$
$$ \begin{aligned} \mathbf{v} \cdot \mathbf{n}_{K_i} &= \xi_{hp}|_{K_j} \cdot \mathbf{n}_{K_i} \quad \text{on } F_{i,j} \\ \mathbf{v} \cdot \mathbf{n}_{K_i} &= \xi_{hp}|_{K_k} \cdot \mathbf{n}_{K_i} \quad \text{on } F_{i,k} \\ \mathbf{v} \cdot \mathbf{n}_{K_i} &= \xi_{hp}|_{K_l} \cdot \mathbf{n}_{K_i} \quad \text{on } F_{i,l} \end{aligned} $$

As above in steps (ii) and (iii), the infinite-dimensional minimizer $\mathbf{v}^*|_{K_i}$ does not belong to the minimization set on the right-hand side of (C.19) since there are three normal trace conditions on the three faces $F_{i,j}$, $F_{i,k}$, and $F_{i,l}$ imposed. Crucially, by property (ii) of Definition B.1 on the enumeration, all the simplices sharing the vertex $\mathbf{a}$ common to the three faces $F_{i,j}$, $F_{i,k}$, and $F_{i,l}$ come sooner in the enumeration and the vertex $\mathbf{a}$ lies in the interior of $\omega$. This enables to construct a suitable $\mathbf{v}$ in this sprit of (C.13) but which now involves Piola mappings from all the simplices sharing the vertex $\mathbf{a}$ except for $K_i$. This is done in a “3-folding” way; the precise formula is the equivalent of [33, equation (5.14)] in the $\mathcal{H}(\text{div})$ case. Existence of a three-color refinement around vertices of [33, Lemma B.3] is crucial at this step. Then (C.7) is established similarly to (C.15). $\square$

*Proof for $d=3$ when a suitable enumeration as per Definition B.1 does not exist.* First recall from Remark B.5 and Lemma B.3 that situations where a suitable enumeration as per Definition B.1 does not exist can arise. For the proof in this case, we will use Adiprasito and Benedetti [1, Theorem A] and the results from [33, 16] on vertex patches (only given by elements sharing a given vertex) to transfer the situation to the above proof with mesh enumerated as per Definition B.1.

Suppose that $\mathcal{T}_\omega$ has no suitable enumeration as per Definition B.1, i.e., invoking Lemma B.3, there exists no shelling as per Definition B.2. Let $\bar{\mathcal{T}}_\omega$ be the barycentric refinement of $\mathcal{T}_\omega$, i.e., the tetrahedral submesh of $\mathcal{T}_\omega$ obtained by inserting a barycenter to each edge, face, and tetrahedron in $\mathcal{T}_\omega$, where the new tetrahedra in each tetrahedron $K \in \mathcal{T}_\omega$ have the barycenter of $K$ as vertex, see Figure 6, left. From Adiprasito and Benedetti [1, Theorem A] we know that $\bar{\mathcal{T}}_\omega$ has a shelling as per Definition B.2. Consequently, from the above proof for shellable meshes, we know that

$$ \min_{\substack{\mathbf{v}_p \in \mathcal{R}\mathcal{T}_p(\bar{\mathcal{T}}_\omega) \cap \mathcal{H}(\text{div}, \omega) \\ \nabla \cdot \mathbf{v}_p = r_{hp}}} \|\tau_{hp} - \mathbf{v}_p\|_{\omega} \lesssim \min_{\substack{\mathbf{v} \in \mathcal{H}(\text{div}, \omega) \\ \nabla \cdot \mathbf{v} = r_{hp}}} \|\tau_{hp} - \mathbf{v}\|_{\omega}. \quad (C.20) $$

Let us denote by $\bar{\mathbf{v}}_p^*$ the (unique) minimizer on the left-hand side of (C.20). From $\bar{\mathbf{v}}_p^*$, piecewise polynomial with respect to the barycentric refined mesh $\bar{\mathcal{T}}_\omega$, we now construct $\xi_{hp} \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_\omega) \cap \mathcal{H}(\text{div}, \omega)$ with $\nabla \cdot \xi_{hp} = r_{hp}$, piecewise polynomial with respect to the original mesh $\mathcal{T}_\omega$, by passing through all vertices $\mathbf{a} \in \mathcal{V}$ of the original mesh $\mathcal{T}_\omega$ and considering the vertex patches $\mathcal{T}_a$, cf. Figure 6, middle. We show that $\xi_{hp}$ has a comparable precision to $\bar{\mathbf{v}}_p^*$ in that

$$ \|\tau_{hp} - \xi_{hp}\|_{\omega} \lesssim \|\tau_{hp} - \bar{\mathbf{v}}_p^*\|_{\omega}. \quad (C.21) $$
<figure><img src="image_7.png" /><figcaption>Figure 6: Barycentric refinement of a tetrahedron $K \in \mathcal{T}_\omega$ with vertices $a_1, a_2, a_3, a_4$ (left), vertex patch $\mathcal{T}_{a_1}$ with vertices of the barycentric refinement indicated in the tetrahedron $K$ (middle), tetrahedron $K$ with 6 subtetrahedra sharing $a_1$ (right)</figcaption></figure>

Thus, from (C.20)–(C.21), we conclude

$$ v_p \in \mathcal{R}\mathcal{T}_p(\bar{T}_{\omega,0}) \cap \mathcal{H}(\mathrm{div},\omega) \qquad \min_{\nabla \cdot v_p = r_{hp}} \|\tau_{hp} - v_p\|_{\omega} \le \|\tau_{hp} - \xi_{hp}\|_{\omega} \lesssim \|\tau_{hp} - \bar{v}_p^*\|_{\omega} \lesssim \min_{\nabla \cdot v = r_{hp}} \|\tau_{hp} - v\|_{\omega}, \quad (C.22) $$

which is the desired result (C.2).

Let us enumerate the vertices from $\mathcal{V}$ as $a_1, \dots, a_{|\mathcal{V}|}$. There is no specific order. We define a sequence of tetrahedral meshes $\bar{T}_{\omega,i}$ of $\omega$, $0 \le i \le |\mathcal{V}|$. We set $\bar{T}_{\omega,0} := \bar{T}_{\omega}$, the barycentric refinement mesh. For all $1 \le i \le |\mathcal{V}|$, $\bar{T}_{\omega,i}$ coincides with $\bar{T}_{\omega,i-1}$ on $\omega \setminus \omega_{a_i}$ and coarsens $\bar{T}_{\omega,i-1}$ inside $\omega_{a_i}$. On the last step, $\bar{T}_{\omega,|\mathcal{V}|} = \bar{T}_{\omega}$, the original mesh. More precisely, let $i = 1$. Then the patch subdomain $\omega_{a_1}$ is formed by tetrahedra from $\mathcal{T}_{a_1}$ sharing the vertex $a_1$, see Figure 6, middle, with barycentric refinement as displayed in Figure 6, left. The mesh $\bar{T}_{\omega,1}$ is created by coarsening of $\bar{T}_{\omega,0}$ inside $\omega_{a_1}$ where each $K \in \mathcal{T}_{a_1}$ is now only refined into 6 tetrahedra having $a_1$ as vertex in place of the barycentric refinement, as indicated in Figure 6, right. Note that $\bar{T}_{\omega,1}$ is indeed a coarsening of $\bar{T}_{\omega,0}$ and keeps the triangles on those faces from $\partial\omega_{a_1}$ that do not share the vertex $a_1$ intact. For $i > 1$, we proceed similarly, always keeping the subtriangulation of those mesh faces on the boundary of $\omega_{a_i}$ that do not share the vertex $a_i$, but coarsening $\bar{T}_{\omega,i-1}$ inside $\omega_{a_i}$. With respect to the original mesh $\mathcal{T}_\omega$, the arising $\bar{T}_{\omega,1}$) adds no barycenter of a tetrahedron from $\mathcal{T}_\omega$ having $a_i$ as vertex; 2) adds no barycenter of a face from $\mathcal{T}_\omega$ having $a_i$ as vertex; 3) adds no barycenter of an edge from $\mathcal{T}_\omega$ having $a_i$ as vertex. Thus, after passing through all vertices $a_i$ from $\mathcal{T}_\omega$, we indeed recover in $\bar{T}_{\omega,|\mathcal{V}|}$ the original mesh $\mathcal{T}_\omega$.

We now construct a sequence of piecewise polynomials $\xi_i \in \mathcal{R}\mathcal{T}_p(\bar{T}_{\omega,i}) \cap \mathcal{H}(\mathrm{div},\omega)$ with $\nabla \cdot \xi_i = r_{hp}$, $0 \le i \le |\mathcal{V}|$. We set $\xi_0 := \bar{v}_p^*$, the minimizer on the barycentric refined mesh from (C.20). On each step $1 \le i \le |\mathcal{V}|$, we set $\xi_i = \xi_{i-1}$ on $\omega \setminus \omega_{a_i}$ (recall that $\bar{T}_{\omega,i}$ coincides with $\bar{T}_{\omega,i-1}$ on $\omega \setminus \omega_{a_i}$) and let

$$ \begin{array}{ll} \xi_i|_{\omega_{a_i}} := \arg \min_{\substack{v_p \in \mathcal{R}\mathcal{T}_p(\bar{T}_{\omega,i}|_{\omega_{a_i}}) \cap \mathcal{H}(\mathrm{div},\omega_{a_i}) \\ \nabla \cdot v_p = r_{hp} \\ v_p \cdot n_{\omega_{a_i}} = \xi_{i-1} \cdot n_{\omega_{a_i}} \text{ on faces from } \partial\omega_{a_i} \text{ not sharing } a_i}} \|\tau_{hp} - v_p\|_{\omega_{a_i}}. & (C.23) \\ & \end{array} $$

Crucially, (C.23) is an energy minimization on a vertex patch $\bar{T}_{\omega,i}|_{\omega_{a_i}}$ with the volume data $\tau_{hp}|_{\omega_{a_i}}$ and $r_{hp}$ from respectively the piecewise polynomial spaces $\mathcal{R}\mathcal{T}_p(\bar{T}_{\omega,i}|_{\omega_{a_i}})$ and $\mathcal{P}_p(\bar{T}_{\omega,i}|_{\omega_{a_i}})$ (since piecewise polynomials on the original mesh $\mathcal{T}_\omega$) and with the Neumann boundary datum $\xi_{i-1} \cdot n_{\omega_{a_i}}$ from the normal trace of $\mathcal{R}\mathcal{T}_p(\bar{T}_{\omega,i}|_{\omega_{a_i}})$ on $\partial\omega_{a_i}$. The Neumann equilibrium condition is also clearly satisfied whenever no face from $\partial\omega_{a_i}$ shares the vertex $a_i$, since it is given by $\xi_{i-1}$ whose divergence is $r_{hp}$. Thus, using [33,
Corollaries 3.3 and 3.8] or [16, Corollary 4.1], we obtain

$$
\begin{align}
\|\tau_{hp} - \xi_i\|_{\omega_{a_i}} &\lesssim \min_{\substack{v \in \mathcal{H}(\text{div}, \omega_{a_i}) \\ \nabla \cdot v = r_{hp}}} \|\tau_{hp} - v_p\|_{\omega_{a_i}} \nonumber \\
&\qquad \text{on faces from } \partial\omega_{a_i} \text{ not sharing } a_i \tag{C.24} \\
&\le \|\tau_{hp} - \xi_{i-1}\|_{\omega_{a_i}}, \nonumber
\end{align}
$$

where the second estimate follows as $\xi_{i-1}|_{\omega_{a_i}}$ lies in $\mathcal{H}(\text{div}, \omega_{a_i})$ and satisfies the divergence and Neumann boundary constraints. We can now conclude by two observations: 1) $\xi_{|V|}$ lies in $\mathcal{RT}_p(\mathcal{T}_\omega) \cap \mathcal{H}(\text{div}, \omega)$ defined over the original mesh $\mathcal{T}_\omega$ and satisfies $\nabla \cdot \xi_{|V|} = r_{hp}$, so that we can take $\xi_{hp} := \xi_{|V|}$; 2) (C.24) gives $\|\tau_{hp} - \xi_i\|_{\omega} \lesssim \|\tau_{hp} - \xi_{i-1}\|_{\omega}$, which yields the requested inequality (C.21) since $|V|$, the number of vertices in $\mathcal{T}_\omega$, is bounded by the shape-regularity parameter $\kappa_{\mathcal{T}_\omega}$ and the ratio $h_\omega/\min_{K \in \mathcal{T}_\omega} h_K$. $\square$

### C.3 Extended vertex and element patches and boundary conditions

We now finally formulate the result precisely in the form needed on step (vii) of the proof of Lemma 5.3.

**Corollary C.2** (p-stable broken $\mathcal{H}(\text{div})$ polynomial extension on extended vertex or element patches). Let $a \in \mathcal{V}_h$ or $K \in \mathcal{T}_h$. Consider the extended vertex patch $\tilde{\mathcal{T}}_a$ or the extended element patch $\tilde{\mathcal{T}}_K$ as per Section 2.2, denoted by $\mathcal{T}_\omega$, with the associated open subdomain $\omega$. Let Assumption 3.1 hold. Let $\tau_{hp} \in \mathcal{RT}_p(\mathcal{T}_\omega)$ and $r_{hp} \in \mathcal{P}_p(\mathcal{T}_\omega)$ be respectively a volume datum, a broken Raviart-Thomas vector-valued piecewise polynomial, and a target divergence, a scalar-valued piecewise polynomial. Then

$$ v_p \in \mathcal{RT}_p(\mathcal{T}_\omega) \cap \mathcal{H}_{0,N}(\text{div}, \omega) \qquad \min_{\nabla \cdot v_p = r_{hp}} \|\tau_{hp} - v_p\|_{\omega} \lesssim \min_{\nabla \cdot v = r_{hp}} \|\tau_{hp} - v\|_{\omega}, \quad (C.25) $$

where the constant hidden in $\lesssim$ only depends on the mesh shape-regularity parameter $\kappa_{\mathcal{T}_\omega} := \max_{K \in \mathcal{T}_\omega} \kappa_K$ and the space dimension $d$.

*Proof.* Let $\omega$ have no face on the Neumann boundary $\Gamma_N$, $|\partial\omega \cap \Gamma_N| = 0$. Then (C.25) follows by Theorem C.1; note that the ratio $h_\omega/\min_{L \in \mathcal{T}_\omega} h_L$ for an extended vertex or element patch $\mathcal{T}_\omega$ only depends on the shape-regularity parameter $\kappa_{\mathcal{T}_\omega}$. In the case $|\partial\omega \cap \Gamma_N| \neq 0$ but when the boundary of $\omega$ does not coincide with the whole Neumann boundary $\Gamma_N$, one can proceed following [33, Section 7], [16, Section 7], and the proof of Theorem C.1, case $d=3$ when $\mathcal{T}_\omega$ has no shelling as per Definition B.2. Here, one designs a sequence of mappings where one can deduce the validity of (C.25) in the case $|\partial\omega \cap \Gamma_N| \neq 0$ from the case $|\partial\omega \cap \Gamma_N| = 0$. $\square$

## D Seemingly overconstrained $\mathcal{RT}_p \cap \mathcal{H}(\text{div})$ minimization on patch subdomains

We summarize here our results on seemingly overconstrained $\mathcal{RT}_p \cap \mathcal{H}(\text{div})$ minimization on patch subdomains.

### D.1 Larger patches and no boundary conditions

We extend here the results of [15, Appendix A] in two directions: for larger patches $\omega$ and without prescription of normal trace boundary conditions on $\partial\omega$. In our application on step (v) of the proof of Lemma 5.3, this corresponds to $\omega = \tilde{\omega}_K$ or $\omega = \tilde{\omega}_a$, the case where $\partial\omega$ does not contain any face from $\Gamma_N$; treatment of the other (boundary) cases is postponed to Section D.2.

**Assumption D.1** (Data for seemingly overconstrained minimization). The volume datum $\tau_{hp}$ and the target divergence $r_{hp}$ satisfy

$$ r_{hp} \in \mathcal{P}_p(\mathcal{T}_\omega), \quad \tau_{hp} \in \mathcal{RT}_p(\mathcal{T}_\omega), \quad (D.1a) $$

$$ (\tau_{hp}, \nabla q_h)_\omega + (r_{hp}, q_h)_\omega = 0 \quad \forall q_h \in \mathcal{P}_1(\mathcal{T}_\omega) \cap H_0^1(\omega), \quad (D.1b) $$

i.e., $r_{hp}$ is a broken (piecewise) p-degree polynomial and $\tau_{hp}$ is a broken Raviart-Thomas piecewise polynomial that are “weakly divergence compatible” for homogeneous continuous piecewise affine polynomials.
**Theorem D.2** (Seemingly overconstrained minimization in the Raviart-Thomas spaces on larger patches and without boundary conditions). Let $\mathcal{T}_\omega$ be a simplicial mesh with the corresponding open, bounded, and Lipschitz polygon or polyhedron $\omega \subset \mathbb{R}^d$, $d = 2, 3$, with $\bar{\omega}$ contractible. Let $p \ge 1$ and let $\tau_{hp}$ and $r_{hp}$ satisfy Assumption D.1. Then

$$ \begin{array}{ll} \min_{\substack{v_p \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_\omega) \cap \mathcal{H}(\mathrm{div},\omega) \\ \nabla \cdot v_p = r_{hp}}} & \|\tau_{hp} - v_p\|_{\omega} \lesssim \min_{\substack{v_p \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_\omega) \cap \mathcal{H}(\mathrm{div},\omega) \\ \nabla \cdot v_p = r_{hp}}} \|\tau_{hp} - v_p\|_{\omega}, \quad (D.2) \\ (v_p, r_h)_K = (\tau_{hp}, r_h)_K & \quad \forall r_h \in [\mathcal{P}_0(K)]^d, \forall K \in \mathcal{T}_\omega \end{array} $$

where both problems have a unique solution and where the constant hidden in $\lesssim$ only depends on the mesh shape-regularity parameter $\kappa_{\mathcal{T}_\omega} := \max_{K \in \mathcal{T}_\omega} \kappa_K$, the ratio $h_\omega / \min_{K \in \mathcal{T}_\omega} h_K$, and the space dimension $d$.

*Proof.* We present (an outline of) the proof for $d = 3$; the two-dimensional case is (much) easier. We follow [15, Appendix A]. Let $\theta_p$ denote the minimizer on the left-hand side of (D.2) and $\bar{\theta}_p$ the minimizer on the right-hand side of (D.2). As the existence and uniqueness of $\bar{\theta}_p$ is standard, cf., e.g., [8, 24], we need to show the existence and uniqueness of $\theta_p$ together with

$$ \|\tau_{hp} - \theta_p\|_{\omega} \lesssim \|\tau_{hp} - \bar{\theta}_p\|_{\omega}. \quad (D.3) $$

(i) Let

$$ \begin{array}{ll} \min_{\substack{v_h \in \mathcal{R}\mathcal{T}_h(\mathcal{T}_\omega) \cap \mathcal{H}(\mathrm{div},\omega) \\ \nabla \cdot v_h = 0}} & \|\tau_{hp} - \bar{\theta}_p - v_h\|_{\omega}. \quad (D.4) \\ (v_h, r_h)_K = (\tau_{hp}, \bar{\theta}_p, r_h)_K & \quad \forall r_h \in [\mathcal{P}_0(K)]^d, \forall K \in \mathcal{T}_\omega \end{array} $$

Note that $\varepsilon_h$ is a low-(first)-order Raviart-Thomas piecewise polynomial. We will show its existence and uniqueness and the stability estimate

$$ \|\varepsilon_h\|_{\omega} \lesssim \|\tau_{hp} - \bar{\theta}_p\|_{\omega} \quad (D.5) $$

below in step (ii). Then, shifting $\bar{\theta}_p$ by $\varepsilon_h$,

$$ \begin{array}{l} \bar{\theta}_p + \varepsilon_h \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_\omega) \cap \mathcal{H}(\mathrm{div}, \omega) \text{ with } \nabla \cdot (\bar{\theta}_p + \varepsilon_h) = r_{hp} \\ (\bar{\theta}_p + \varepsilon_h, r_h)_K = (\tau_{hp}, r_h)_K \quad \forall r_h \in [\mathcal{P}_0(K)]^d, \forall K \in \mathcal{T}_\omega. \end{array} $$

Thus, $\bar{\theta}_p + \varepsilon_h$ belongs to the minimization set on the left-hand side of (D.2). Since this minimization is convex, this establishes the existence and uniqueness of $\theta_p$. Moreover,

$$ \|\tau_{hp} - \theta_p\|_{\omega} \le \|\tau_{hp} - \bar{\theta}_p - \varepsilon_h\|_{\omega} \le \|\tau_{hp} - \bar{\theta}_p\|_{\omega} + \|\varepsilon_h\|_{\omega} \lesssim \|\tau_{hp} - \bar{\theta}_p\|_{\omega}, \quad (D.6) $$

which is the desired result (D.3).

(ii) To establish the existence and uniqueness of $\varepsilon_h$ from (D.4), we need to show that the minimization set in (D.4) is not empty. By (D.1b) and by the Green theorem, recalling that $\bar{\theta}_p \in \mathcal{R}\mathcal{T}_p(\mathcal{T}_\omega) \cap \mathcal{H}(\mathrm{div}, \omega)$ with $\nabla \cdot \bar{\theta}_p = r_{hp}$, we see that the datum $\tau_{hp} - \bar{\theta}_p$ satisfies

$$ (\tau_{hp} - \bar{\theta}_p, \nabla q_h)_\omega = -(r_{hp}, q_h)_\omega + (r_{hp}, q_h)_\omega = 0 \quad \forall q_h \in \mathcal{P}_1(\mathcal{T}_\omega) \cap H_0^1(\omega). $$

This is a set of the form studied in [15, Lemma A.5] on vertex patches $\mathcal{T}_a$ and with zero normal trace boundary conditions on $\partial\omega_a$. This proof generalizes to the current setting just as that of the proof of Theorem C.1. As for the stability estimate (D.5), please note that $\tau_{hp} - \bar{\theta}_p$ is the only datum in problem (D.4), which implies (D.5) up to a generic constant with unknown dependencies. The fact that these dependencies only include the shape-regularity parameter $\kappa_{\mathcal{T}_\omega}$ of the mesh $\mathcal{T}_\omega$, the ratio $h_\omega / \min_{K \in \mathcal{T}_\omega} h_K$, and the space dimension $d$ follows by scaling arguments as in [15, Proof of Lemma A.4]; the fact that $\varepsilon_h$ is merely a first-order Raviart-Thomas piecewise polynomial is decisive for p-robustness. $\square$

## D.2 Extended vertex and element patches and boundary conditions

We now finally formulate the result precisely in the form needed on step (v) of the proof of Lemma 5.3. The proof follows that of Lemma C.2.
**Corollary D.3** (Seemingly overconstrained minimization in the Raviart-Thomas spaces on extended vertex or element patches). Let $\mathbf{a} \in V_h$ or $K \in T_h$. Consider the extended vertex patch $\tilde{T}_\mathbf{a}$ or the extended element patch $\tilde{T}_K$ as per Section 2.2, denoted by $T_\omega$, with the associated open subdomain $\omega$. Let Assumption 3.1 hold. Let, for $p \ge 1$,

$$r_{hp} \in \mathcal{P}_p(T_\omega), \quad \tau_{hp} \in \mathcal{RT}_p(T_\omega), \qquad (\text{D.6a})$$

$$(\tau_{hp}, \nabla q_h)_\omega + (r_{hp}, q_h)_\omega = 0 \quad \forall q_h \in \mathcal{P}_1(T_\omega) \cap H_{0, \partial\omega \setminus \Gamma_N}^1(\omega). \qquad (\text{D.6b})$$

Then

$$\min_{\substack{v_p \in \mathcal{RT}_p(T_\omega) \cap \mathcal{H}_{0,N}(\text{div}, \omega) \\ \nabla \cdot v_p = r_{hp} \\ (v_p, r_h)_K = (\tau_{hp}, r_h)_K \quad \forall r_h \in [P_0(K)]^d, \forall K \in \mathcal{T}_\omega}} \|\tau_{hp} - v_p\|_{\omega} \lesssim \min_{\substack{v_p \in \mathcal{RT}_p(T_\omega) \cap \mathcal{H}_{0,N}(\text{div}, \omega) \\ \nabla \cdot v_p = r_{hp}}} \|\tau_{hp} - v_p\|_{\omega}, \qquad (\text{D.7})$$

where the constant hidden in $\lesssim$ only depends on the mesh shape-regularity parameter $\kappa_{T_\omega} := \max_{K \in T_\omega} \kappa_K$ and the space dimension $d$.

## References

[1] Adiprasito, K. A., and Benedetti, B. Subdivisions, shellability, and collapsibility of products. *Combinatorica* **37** (2017), 1–30. https://doi.org/10.1007/s00493-016-3149-8.

[2] Ainsworth, M., and Demkowicz, L. Explicit polynomial preserving trace liftings on a triangle. *Math. Nachr.* **282** (2009), 640–658. https://doi.org/10.1002/mana.200610762.

[3] Arnold, D., and Guzmán, J. Local $L^2$-bounded commuting projections in FEEC. *ESAIM Math. Model. Numer. Anal.* **55** (2021), 2169–2184. https://doi.org/10.1051/m2an/2021054.

[4] Aurada, M., Feischl, M., Kemetmüller, J., Page, M., and Praetorius, D. Each $H^{1/2}$-stable projection yields convergence and quasi-optimality of adaptive FEM with inhomogeneous Dirichlet data in $\mathbb{R}^d$. *ESAIM Math. Model. Numer. Anal.* **47** (2013), 1207–1235. https://doi.org/10.1051/m2an/2013069.

[5] Babuška, I., and Suri, M. The h-p version of the finite element method with quasi-uniform meshes. *RAIRO Modél. Math. Anal. Numér.* **21** (1987), 199–238. http://dx.doi.org/10.1051/m2an/1987210201991.

[6] Bespalov, A., and Heuer, N. A new $H(\text{div})$-conforming $p$-interpolation operator in two dimensions. *ESAIM Math. Model. Numer. Anal.* **45** (2011), 255–275. https://doi.org/10.1051/m2an/2010039.

[7] Bing, R. H. Some aspects of the topology of 3-manifolds related to the Poincaré conjecture. In *Lectures on Modern Mathematics, Vol. II*. Wiley, New York-London-Sydney, 1964, pp. 93–128.

[8] 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.

[9] 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.

[10] Bruggesser, H., and Mani, P. Shellable decompositions of cells and spheres. *Math. Scand.* **29** (1971), 197–205 (1972). https://doi.org/10.7146/math.scand.a-11045.

[11] Canuto, C., Nochetto, R. H., Stevenson, R. P., and Verani, M. Convergence and optimality of hp-AFEM. *Numer. Math.* **135** (2017), 1073–1119. https://doi.org/10.1007/s00211-016-0826-x.

[12] Carstensen, C., Peterseim, D., and Schedensack, M. Comparison results of finite element methods for the Poisson model problem. *SIAM J. Numer. Anal.* **50** (2012), 2803–2823. https://doi.org/10.1137/110845707.
[13] Chaumont-Frelet, T., Licht, M. W., and Vohralík, M. Computable Poincaré-Friedrichs constants for the $L^p$ de Rham complex over convex domains and domains with shellable triangulations. In preparation, 2025.

[14] Chaumont-Frelet, T., and Vohralík, M. Equivalence of local-best and global-best approximations in $H(\text{curl})$. *Calcolo* **58** (2021), 53. https://doi.org/10.1007/s10092-021-00430-9.

[15] Chaumont-Frelet, T., and Vohralík, M. *p*-robust equilibrated flux reconstruction in $H(\text{curl})$ based on local minimizations. Application to a posteriori analysis of the curl-curl problem. *SIAM J. Numer. Anal.* **61** (2023), 1783-1818. https://doi.org/10.1137/21M141909X.

[16] 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. *Found. Comput. Math.* (2024). DOI 10.1007/s10208-024-09674-7, https://doi.org/10.1007/s10208-024-09674-7.

[17] Chaumont-Frelet, T., and Vohralík, M. A stable local commuting projector and optimal hp approximation estimates in $H(\text{curl})$. *Numer. Math.* **156** (2024), 2293-2342. https://doi.org/10.1007/s00211-024-01431-w.

[18] Christiansen, S. H., and Winther, R. Smoothed projections in finite element exterior calculus. *Math. Comp.* **77** (2008), 813-829. http://dx.doi.org/10.1090/S0025-5718-07-02081-9.

[19] Clément, P. Approximation by finite element functions using local regularization. *RAIRO Anal. Numer.* **9** (1975), 77-84.

[20] Cook, II, D. Simplicial decomposability. *J. Softw. Algebra Geom.* **2** (2010), 20-23. https://doi.org/10.2140/jsag.2010.2.20.

[21] 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.

[22] Demkowicz, L. Polynomial exact sequences and projection-based interpolation with application to Maxwell equations. In *Mixed finite elements, compatibility conditions, and applications*, D. Boffi, F. Brezzi, L. F. Demkowicz, R. G. Durán, R. S. Falk, and M. Fortin, Eds., vol. **1939** of *Lecture Notes in Mathematics*. Springer-Verlag, Berlin, 2008, pp. 101-158. Lectures given at the C.I.M.E. Summer School held in Cetraro, June 26-July 1, 2006, Edited by D. Boffi and L. Gastaldi, https://doi.org/10.1007/978-3-540-78319-0_3.

[23] Demkowicz, L. Lecture notes. *Lecture Notes on Energy Spaces*. The University of Texas at Austin, 2018. https://users.oden.utexas.edu/~leszek/classes/EM394H/book2.pdf.

[24] Demkowicz, L. Lecture notes. *Mathematical Theory of Finite Elements*. The University of Texas at Austin, 2023. https://users.oden.utexas.edu/~leszek/classes/EM394H/book.pdf.

[25] Demkowicz, L., and Buffa, A. $H^1$, $H(\text{curl})$ and $H(\text{div})$-conforming projection-based interpolation in three dimensions. Quasi-optimal *p*-interpolation estimates. *Comput. Methods Appl. Mech. Engrg.* **194** (2005), 267-296. https://doi.org/10.1016/j.cma.2004.07.007.

[26] Demkowicz, L., Gopalakrishnan, J., and Schöberl, J. Polynomial extension operators. Part III. *Math. Comp.* **81** (2012), 1289-1326. http://dx.doi.org/10.1090/S0025-5718-2011-02536-6.

[27] Dochtermann, A., Nair, R., Schweig, J., Van Tuyl, A., and Woodroofe, R. Simplicial complexes with many facets are vertex decomposable. *Electron. J. Combin.* **31** (2024), Paper No. 4.34, 11. https://doi.org/10.37236/12984.

[28] Ern, A., Gudi, T., Smears, I., and Vohralík, M. Equivalence of local- and global-best approximations, a simple stable local commuting projector, and optimal hp approximation estimates in $H(\text{div})$. *IMA J. Numer. Anal.* **42** (2022), 1023-1049. http://dx.doi.org/10.1093/imanum/draa103.

[29] Ern, A., and Guermond, J.-L. Mollification in strongly Lipschitz domains with application to continuous and discrete de Rham complexes. *Comput. Methods Appl. Math.* **16** (2016), 51-75. https://doi.org/10.1515/cmam-2015-0034.
[30] Ern, A., and Guermond, J.-L. Finite element quasi-interpolation and best approximation. *ESAIM Math. Model. Numer. Anal.* **51** (2017), 1367–1385. https://doi.org/10.1051/m2an/2016066.

[31] 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/10.1007/978-3-030-56341-7.

[32] Ern, A., Guzmán, J., Potu, P., and Vohralík, M. Local $L^2$-bounded commuting projections using discrete local problems on Alfeld splits. HAL Preprint 04931497, submitted for publication, https://hal.inria.fr/hal-04931497, 2025.

[33] 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.

[34] Falk, R. S., and Winther, R. Local bounded cochain projections. *Math. Comp.* **83** (2014), 2631–2656. http://dx.doi.org/10.1090/S0025-5718-2014-02827-5.

[35] Falk, R. S., and Winther, R. Construction of polynomial preserving cochain extensions by blending. *Math. Comp.* **92** (2023), 1575–1594. https://doi.org/10.1090/mcom/3819.

[36] Falk, R. S., and Winther, R. Local space-preserving decompositions for the bubble transform. *Found. Comput. Math.* (2025). https://doi.org/10.1007/s10208-025-09700-2.

[37] Gawlik, E., Holst, M. J., and Licht, M. W. Local finite element approximation of Sobolev differential forms. *ESAIM Math. Model. Numer. Anal.* **55** (2021), 2075–2099. https://doi.org/10.1051/m2an/2021034.

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

[39] Goaoc, X., Paták, P., Patáková, Z., Tancer, M., and Wagner, U. Shellability is NP-complete. *J. ACM* **66** (2019), Art. 21, 18. https://doi.org/10.1145/3314024.

[40] Kozlov, D. *Combinatorial algebraic topology*, vol. **21** of *Algorithms and Computation in Mathematics*. Springer, Berlin, 2008. https://doi.org/10.1007/978-3-540-71962-5.

[41] Licht, M. W. Smoothed projections and mixed boundary conditions. *Math. Comp.* **88** (2019), 607–635. https://doi.org/10.1090/mcom/3330.

[42] Lutz, F. H. Combinatorial 3-manifolds with 10 vertices. *Beiträge Algebra Geom.* **49** (2008), 97–106.

[43] Melenk, J. M., and Rojik, C. On commuting $p$-version projection-based interpolation on tetrahedra. *Math. Comp.* **89** (2020), 45–87. https://doi.org/10.1090/mcom/3454.

[44] Moriyama, S. Deciding shellability of simplicial complexes with $h$-assignments. *IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences E94.A* (2011), 1238–1241. https://doi.org/10.1587/transfun.E94.A.1238.

[45] Paták, P., and Tancer, M. Shellability is hard even for balls. In *STOC'23—Proceedings of the 55th Annual ACM Symposium on Theory of Computing* ([2023] ©2023), ACM, New York, pp. 1271–1284. https://doi.org/10.1145/3564246.3585152.

[46] 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.

[47] Schöberl, J., Melenk, J. M., Pechstein, C., and Zaglmayr, S. Additive Schwarz preconditioning for $p$-version triangular and tetrahedral finite elements. *IMA J. Numer. Anal.* **28** (2008), 1–24. http://dx.doi.org/10.1093/imanum/drl046.

[48] Veeser, A. Approximating gradients with continuous piecewise polynomial functions. *Found. Comput. Math.* **16** (2016), 723–750. http://dx.doi.org/10.1007/s10208-015-9262-z.
[49] Vohralík, M. On the discrete Poincaré–Friedrichs inequalities for nonconforming approximations of the Sobolev space $H^1$. *Numer. Funct. Anal. Optim.* **26** (2005), 925–952. http://dx.doi.org/10.1080/01630560500444533.

[50] Vohralík, M. *p*-robust equivalence of global continuous and local discontinuous approximation, a *p*-stable local projector, and optimal elementwise *hp* approximation estimates in $H^1$. HAL Preprint 04436063, submitted for publication, https://hal.inria.fr/hal-04436063, 2024.

[51] Vohralík, M. *p*-robust equivalence of global continuous constrained and local discontinuous unconstrained approximation, a *p*-stable local commuting projector, and optimal elementwise *hp* approximation estimates in *H*(curl). In preparation, 2025.

[52] Ziegler, G. M. Shelling polyhedral 3-balls and 4-polytopes. *Discrete Comput. Geom.* **19** (1998), 159–174. https://doi.org/10.1007/PL00009339.

[53] Ziegler, G. M. Updated seventh printing of the first edition. *Lectures on polytopes*, vol. **152** of *Graduate Texts in Mathematics*. Springer-Verlag, New York, 2007. http://dx.doi.org/10.1007/978-1-4613-8431-1.