# Unfitted hybrid high-order methods stabilized by polynomial extension for elliptic interface problems

Erik Burman, Alexandre Ern, Romain Mottier

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

## To cite this version:

Erik Burman, Alexandre Ern, Romain Mottier. Unfitted hybrid high-order methods stabilized by polynomial extension for elliptic interface problems. 2025. hal-04990318v2

HAL Id: hal-04990318

https://hal.science/hal-04990318v2

Preprint submitted on 9 Oct 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Unfitted hybrid high-order methods stabilized by polynomial extension for elliptic interface problems

Erik Burman†

Alexandre Ern‡

Romain Mottier§

## Abstract

In this work, we study the design and analysis of a novel hybrid high-order (HHO) method on unfitted meshes. HHO methods rely on a pair of unknowns, combining polynomials attached to the mesh faces and the mesh cells. In the unfitted framework, the interface can cut through the mesh cells in a very general fashion, and the polynomial unknowns are doubled in the cut cells and the cut faces. In order to avoid the ill-conditioning issues caused by the presence of small cut cells, the novel approach introduced herein is to use polynomial extensions in the definition of the gradient reconstruction operator. Stability and consistency results are established, leading to optimally decaying error estimates. The theory is illustrated by numerical experiments.

**Keywords.** Hybrid high-order methods, Unfitted methods, Polynomial extension, Interface problem, Curved boundary.

**Mathematics Subjects Classification.** 65N15, 65N30, 65N85.

## 1 Introduction

The meshing of complex domains or of internal interfaces is an important bottleneck in the large-scale computation of solutions to partial differential equations. It is well-known that, for standard finite element methods, the mesh must match the geometry sufficiently well, or accuracy will be lost. This is particularly important for high-order methods. Recently, many methods have been developed aiming at simplifying the meshing procedure by allowing for a larger set of geometric shapes for the computational cells. Examples of discretization methods supporting fairly general polygonal (or polyhedral) cells are the discontinuous Galerkin (dG) method [22, 17], the virtual element method (VEM) [8], and hybrid nonconforming methods in various flavors, such as the hybridizable dG (HDG) method [19], the weak Galerkin (WG) method [34], the nonconforming VEM (ncVEM) method [2] and the hybrid high-order (HHO) method [23]. The tight connections between HDG, WG, and HHO methods are highlighted in [18]. These three methods are formulated by means of degrees of freedom (dofs) both in the bulk of the elements and on the element faces, whereas dG methods only use bulk dofs. Using face dofs leads to advantageous properties, such as the elimination of bulk dofs by static condensation and the absence of minimal threshold to weigh the stabilization.

Nevertheless, difficulties remain when using hybrid nonconforming methods in the presence of curved boundaries or interfaces, whereas the treatment of cells with curved faces in dG methods is more straightforward. One possibility is to enrich (with non-polynomial functions) the set of unknowns on curved faces, as in [35] for HHO and in [9] for ncVEM. Another possibility is to work with mesh cells having flat faces, but to build the geometric representation directly into

<footnote>†Department of Mathematics, University College London, London WC1E 6BT, UK.</footnote>

<footnote>‡CERMICS, ENPC, Institut Polytechnique de Paris, F-77455 Marne-la-Vallée cedex 2, and Centre Inria de Paris, 48 rue Barrault, CS 61534, F-75647 Paris, France.</footnote>

<footnote>§CEA, DAM, DIF, F-91297 Arpajon, France, and CERMICS, ENPC, Institut Polytechnique de Paris, F-77455 Marne-la-Vallée cedex 2, and Centre Inria de Paris, 48 rue Barrault, CS 61534, F-75647 Paris, France.</footnote>
the discretization, in such a way that solution accuracy is not lost although the mesh does not
respect the geometry. This idea has been developed first in the framework of unfitted finite
element methods. These methods originate from the penalty method [3, 6] and were further
developed for conforming FEM in [5]. Combining the ideas of unfitted FEM with Nitsche's
method for interface problems led to a consistent unfitted finite element method for interface
problems in [25]. These ideas were then generalized to include more complex coupling problems
including PDEs on surfaces and multiphysics coupling, leading to the cutFEM framework [13].
Unfitted methods have also been considered for dG [7, 27, 30], HDG [20, 32], and HHO [16]
methods. A crucial problem for unfitted methods is to handle the instabilities that can occur
for certain intersections of the mesh and the interface. In the cutFEM framework using a
conforming approximation in each subdomain, this is solved by the addition of certain stabilized
terms [11], whereas in dG methods, which are more flexible with respect to cell geometry, this
can be solved by cell agglomeration [33, 27]. This technique can also be used in the context of
conforming finite elements [4] (therein called aggregated FEM).

In the present work, we continue the development of unfitted HHO methods. The combina-
tion of cell agglomeration and a local weak coupling derived from the HHO framework has been
successfully applied to elliptic interface problems [12], incompressible Stokes flows [14], and wave
propagation in heterogeneous media [15]. Nevertheless, if these methods have to be integrated
in an existing HHO code, the cell-agglomeration procedure becomes a disruptive element since
it requires modifying part of the mesh in the interface zone. It is often desirable to keep the
underlying mesh fixed and modify instead the algebraic structure of the bulk unknowns. One
technique that fits this purpose is the idea of polynomial extension from interior elements to
stabilize the method. The idea was first introduced in the context of Lagrange multipliers [26]
and then using isogeometric approximation and Nitsche's method in [10]. The purpose of this
work is to develop and analyze this technique for hybrid nonconforming methods, focusing on
HHO methods. Owing to the close links between HDG, WG, ncVEM and HHO methods, the
present results are relevant to the development of these other methods as well. Moreover, we
focus on elliptic interface problems, but observe that elliptic problems posed on curved domains
can be handled by the same techniques.

Let us briefly outline the main idea. Recall that a central ingredient in HHO methods is a
local gradient reconstruction from the local cell and face unknowns. The stability and consistency
properties of this reconstructed gradient crucially hinge on a discrete trace inequality on the
normal derivatives at the cell faces. In the case of ill-cut cells, the constant in this inequality
degenerates. One remedy, explored in [16, 12], is to use cell agglomeration around ill-cut cells.
This procedure produces an aggregated mesh solely composed of well-cut cells, whose size is
close to that of the original mesh. The remedy explored herein is to use the cell and face
polynomials of ill-cut cells in the gradient reconstruction of a neighboring well-cut cell. With
this approach, the mesh cells remain unmodified, but the stencil associated with the gradient
reconstruction operator is slightly extended. Moreover, the cell polynomials of the ill-cut cells
are stabilized by using the cell polynomials of the associated well-cut cells, in the same spirit as
the so-called direct ghost penalty method devised in [31, 28]. The main question is then whether
the resulting method maintains the above stability and optimal approximation properties. The
proof of these properties is by no means straightforward and is the main contribution of this work.
Interestingly, the main analysis tools (discrete trace inequality and polynomial approximation)
are of independent interest.
This work is organized as follows. In Section 2, we introduce the model problem. In Section 3, we present the unfitted HHO method stabilized using polynomial extensions. The stability and error analysis of the method is covered in Section 4. Numerical results are presented in Section 5. Finally, the proofs of the analysis tools (discrete trace inequality and polynomial approximation) on unfitted meshes is outlined in Section 6.

## 2 Model problem

Let $\Omega$ be a polyhedral domain in $\mathbb{R}^d$, $d \in \{2, 3\}$ (open, bounded, connected, Lipschitz subset of $\mathbb{R}^d$), and consider a partition of $\Omega$ into two disjoint subdomains $\bar{\Omega} = \bar{\Omega}_1 \cup \bar{\Omega}_2$ with interface $\Gamma := \partial\Omega_1 \cap \partial\Omega_2$, as illustrated in Figure 1. For all $i \in \{1, 2\}$, we set $\bar{i} := 3 - i$, so that $\Omega \setminus \Omega^i = \Omega^{\bar{i}} \cup \Gamma$. For simplicity, we assume that the interface $\Gamma$ is of class $C^2$ and that it does not touch the boundary of $\Omega$. The unit normal vector $n_\Gamma$ to $\Gamma$ conventionally points from $\Omega_1$ to $\Omega_2$. For a smooth enough function $v$ defined on $\Omega_1 \cup \Omega_2$, setting $v_i := v|_{\Omega_i}$ for all $i \in \{1, 2\}$, we denote its jump across $\Gamma$ as $[[v]]_\Gamma := v_1|_\Gamma - v_2|_\Gamma$.

<figure><img src="image_3.png" /><figcaption>Fig. 1: Model problem</figcaption></figure>

Our goal is to approximate the exact solution $u^{\text{ex}} \in H^1(\Omega_1 \cup \Omega_2)$ of the following elliptic interface problem:

$$ -\nabla \cdot (\kappa \nabla u^{\text{ex}}) = f \quad \text{in } \Omega_1 \cup \Omega_2, \qquad (1a) $$

$$ [[u^{\text{ex}}]]_{\Gamma} = g_D \quad \text{on } \Gamma, \qquad (1b) $$

$$ [[\kappa \nabla u^{\text{ex}}]]_{\Gamma} \cdot n_{\Gamma} = g_N \quad \text{on } \Gamma, \qquad (1c) $$

$$ u^{\text{ex}} = 0 \quad \text{on } \partial\Omega, \qquad (1d) $$

with $f \in L^2(\Omega)$, $g_D \in H^{\frac{1}{2}}(\Gamma)$ and $g_N \in L^2(\Gamma)$. For simplicity, we consider a homogeneous Dirichlet boundary condition on $\partial\Omega$ (see (1d)). We assume that the diffusion coefficient $\kappa$ is scalar-valued and that $\kappa_i := \kappa|_{\Omega_i}$ is constant for each $i \in \{1, 2\}$. To fix the ideas, we assume that $\kappa_1 \le \kappa_2$. Our analysis covers the strongly contrasted case where $\kappa_1 \ll \kappa_2$. We emphasize that it is the property $\kappa_1 \le \kappa_2$ that fixes the numbering of the subdomains, and not the fact that one of them touches the boundary. Thus, the numbering in Figure 1 is only illustrative. We also notice that the above assumptions on $\Gamma$ and $\kappa$ can be lifted with additional technicalities.

The weak formulation can be written as

$$ \text{Find } u^{\text{ex}} \in V_{g_D}, \quad a(u^{\text{ex}}, w) = \ell(w) \quad \forall w \in V_0, \qquad (2) $$

with $V_{g_D} := \{v \in H^1(\Omega_1 \cup \Omega_2) | [[v]]_\Gamma = g_D \text{ on } \Gamma, v = 0 \text{ on } \partial\Omega\}$. Notice that $V_0 = H_0^1(\Omega)$. The
symmetric bilinear form $a$ and the linear form $\ell$ are defined as

$$ a(v, w) := \sum_{i \in \{1,2\}} \kappa_i (\nabla v_i, \nabla w_i)_{\Omega_i} \quad (3a) $$

$$ \ell(w) := (f, w)_{\Omega} + (g_N, w)_{\Gamma}. \quad (3b) $$

Here and in what follows, for a measurable subset $S \subset \Omega$, we denote by $(\cdot, \cdot)_S$ the $L^2(S)$-inner product with appropriate Lebesgue measure and by $\|\cdot\|_S$ the corresponding norm. We notice that, in the weak formulation (2), the jump condition (1b) and the boundary condition (1d) are enforced explicitly in the definition of the functional space $V_{g_D}$, whereas the jump condition (1c) is enforced weakly (i.e., results from the weak formulation).

# 3 Unfitted HHO methods

In this section, we present the unfitted HHO method with stabilization of ill-cut cells using polynomial extensions.

## 3.1 Unfitted meshes

Let $\mathcal{T}_h$ be a mesh that covers $\Omega$ exactly. A generic mesh cell is denoted $T \in \mathcal{T}_h$, it is considered to be an open set, $h_T$ denotes its diameter, and $n_T$ its unit outward normal. The mesh size is defined as $h := \max_{T \in \mathcal{T}_h} h_T$. For all $T \in \mathcal{T}_h$, its neighboring layers $\Delta_j(T) \subset \mathcal{T}_h$ are defined by induction as $\Delta_0(T) := T$ and $\Delta_j(T) := \{T' \in \mathcal{T}_h \mid \overline{T'} \cap \overline{\Delta_{j-1}(T)} \neq \emptyset\}$ for all $j \in \mathbb{N}$. In what follows, we use $j \in \{1, 2\}$. To alleviate the notation, we set $\Delta(T) := \Delta_1(T)$. We make the mild assumption that, for every mesh cell $T \in \mathcal{T}_h$, its convex hull, $\text{conv}(T)$, satisfies

$$ \text{conv}(T) \subset \Delta(T). \quad (4) $$

For instance, this assumption trivially holds true if all the mesh cells are convex sets, and it remains a reasonable assumption if the mesh cells are nonconvex. One can also assume more generally that $\text{conv}(T) \subset \Delta_{n_0}(T)$ for some fixed integer $n_0$.

The mesh $\mathcal{T}_h$ is typically composed of cells having a simple shape (although this is not a necessary assumption); in our illustrations, we consider a structured mesh composed of squares. Quite importantly, the mesh $\mathcal{T}_h$ is not fitted to the interface $\Gamma$, so that $\Gamma$ can cut arbitrarily across any mesh cell. For all $T \in \mathcal{T}_h$, we set $T^\Gamma := T \cap \Gamma$ and $T^i := T \cap \Omega^i$ for all $i \in \{1, 2\}$ (some of these sets may be empty). The current implementation assumes that both sub-cells $T^i$ are connected, although this assumption is not needed from a theoretical viewpoint.

The first obvious partition among the mesh cells is between uncut and cut cells. Moreover, among cut cells, we distinguish between well-cut and ill-cut cells. Specifically, we fix a parameter $\vartheta \in (0, 1)$ and say that a cut cell $T$ is a well-cut cell if, for all $i \in \{1, 2\}$, $T^i$ contains a ball of radius $\vartheta h_T$, whereas $T$ is an ill-cut cell if this condition fails. It is shown in [16, Lemma 6.2] that, it the mesh size $h$ is small enough (with respect to the curvature of the interface) and the parameter $\vartheta$ is small enough (with respect to the shape-regularity parameter of the mesh), then, for every ill-cut cell, the above condition can fail on one and only one side of $T$. In what follows, we assume that this is always the case. The above considerations lead to the following partitions:

$$ \mathcal{T}_h := \mathcal{T}_h^{\text{uncut}} \cup \mathcal{T}_h^{\text{cut}}, \quad \mathcal{T}_h^{\text{cut}} := \mathcal{T}_h^{\text{OK}} \cup \mathcal{T}_h^{\text{KO}}, \quad (5) $$
where $\mathcal{T}_h^{\text{uncut}}$ collects the uncut cells, $\mathcal{T}_h^{\text{cut}}$ the cut cells, $\mathcal{T}_h^{\text{OK}}$ the well-cut cells, and $\mathcal{T}_h^{\text{KO}}$ the ill-cut cells. Setting $\mathcal{T}_h^i := \{T \in \mathcal{T}_h \mid T \subset \Omega_i\}$ for all $i \in \{1, 2\}$, we have $\mathcal{T}_h^{\text{uncut}} = \mathcal{T}_h^1 \cup \mathcal{T}_h^2$. For every $T \in \mathcal{T}_h^{\text{uncut}}$, we define $\iota(T) \in \{1, 2\}$ as the index such that $T \in \mathcal{T}_h^{\iota(T)}$. Moreover, for all $T \in \mathcal{T}_h^{\text{KO}}$, we define $\iota(T) \in \{1, 2\}$ as the index for which the above ball condition fails. We then consider the partition $\mathcal{T}_h^{\text{KO}} = \mathcal{T}_h^{\text{KO},1} \cup \mathcal{T}_h^{\text{KO},2}$, where $T \in \mathcal{T}_h^{\text{KO},i}$ if $\iota(T) = i$. In conclusion, a mesh partition that is finer than (5) is

$$ \mathcal{T}_h := \underbrace{\mathcal{T}_h^1 \cup \mathcal{T}_h^2}_{\mathcal{T}_h^{\text{uncut}}} \cup \underbrace{\mathcal{T}_h^{\text{OK}} \cup \mathcal{T}_h^{\text{KO},1}}_{=\mathcal{T}_h^{\text{cut}}} \cup \underbrace{\mathcal{T}_h^{\text{KO},2}}_{=\mathcal{T}_h^{\text{cut}}}. \qquad (6) $$

<figure><img src="image_4.png" /><figcaption>**Fig. 2:** Left: Illustration of the different types of cell in the unfitted mesh. Right: Zoom on the local degrees of freedom in a cut cell $T \in \mathcal{T}_h^{\text{KO}}$; here, the approximation is affine in the sub-cells and constant on the sub-faces.</figcaption></figure>

The polynomial extension technique we propose relies on a pairing operator, which maps every ill-cut cell $S \in \mathcal{T}_h^{\text{KO}}$ to an uncut or a well-cut cell that belongs to the neighborhood $\Delta(S)$. We denote the pairing operator as $\mathcal{N}: \mathcal{T}_h^{\text{KO}} \to \mathcal{T}_h$, so that we have

$$ \mathcal{N}: \mathcal{T}_h^{\text{KO}} \ni S \mapsto T \in (\mathcal{T}_h^i \cup \mathcal{T}_h^{\text{OK}} \cup \mathcal{T}_h^{\text{KO},\bar{i}}) \cap \Delta(S), \quad i := \iota(T). \qquad (7) $$

Letting $\mathcal{N}_i := \mathcal{N}|_{\mathcal{T}_h^{\text{KO},i}}$ for all $i \in \{1, 2\}$, we have more specifically

$$ \mathcal{N}_i: \mathcal{T}_h^{\text{KO},i} \ni S \mapsto T \in (\mathcal{T}_h^i \cup \mathcal{T}_h^{\text{OK}} \cup \mathcal{T}_h^{\text{KO},\bar{i}}) \cap \Delta(S). \qquad (8) $$

We notice that the set $\mathcal{T}_h^{\text{KO},i}$ can be partitioned as follows:

$$ \mathcal{T}_h^{\text{KO},i} = \mathcal{N}_i^{-1}(\mathcal{T}_h^i) \cup \mathcal{N}_i^{-1}(\mathcal{T}_h^{\text{OK}}) \cup \mathcal{N}_i^{-1}(\mathcal{T}_h^{\text{KO},\bar{i}}), \quad \forall i \in \{1, 2\}. \qquad (9) $$

The pairing operator is essentially the one constructed in [12, Algorithm 1]. One difference, however, is that after Step 3 of this algorithm, we need to add an additional Step 4 performing the same operations as in Step 3, but for any mesh cell in $\mathcal{T}_h^{\text{KO},2}$. One way of proceeding is to set $\mathcal{N}_2(T) := S$, whenever $T := \mathcal{N}_1(S) \in \mathcal{T}_h^{\text{KO},2}$ for some cell $S \in \mathcal{T}_h^{\text{KO},1}$. Thus, the well-cut part of $S$ in $\Omega_2$ will be used for $T$, and the well-cut part of $T$ in $\Omega_1$ will be used for $S$. The reason why Step 4 is not needed in the context of cell-agglomeration is that by agglomerating $T$ and $S$ in Step 3, $T$ does not require any further stabilization. An illustration of the pairing operators $\mathcal{N}_1$ and $\mathcal{N}_2$ is shown in Figure 3, and two concrete examples are displayed in the central column of Figure 5 below. Notice that it may happen that some cells in $\mathcal{T}_h^{\text{OK}}$ are in the image of both $\mathcal{N}_1$ and $\mathcal{N}_2$. The above choice in Step 4 aims at reducing the possible impact of the pairing
<figure><img src="image_5.png" /><figcaption>Fig. 3: Pairing operators $\mathcal{N}_1$ and $\mathcal{N}_2$ acting on the ill-cut cells.</figcaption></figure>

operator on the resulting stiffness matrix. There may be room for some further improvement, but we do not delve on this aspect here since the impact of the pairing operator on the stencil is already quite modest, as highlighted in Figure 8 below.

The mesh faces are collected in the set $\mathcal{F}_h$, and we set $\mathcal{F}_h^i := \{F \in \mathcal{F}_h | F \subset \Omega_i\}$ for all $i \in \{1, 2\}$. For every mesh face $F \in \mathcal{F}_h^\Gamma := \mathcal{F}_h \setminus (\mathcal{F}_h^1 \cup \mathcal{F}_h^2)$ cut by the interface $\Gamma$, we define its two sub-faces as $F^i := F \cap \Omega_i$ for all $i \in \{1, 2\}$. For all $T \in \mathcal{T}_h$, $\mathcal{F}_{\partial T}$ is the collection of mesh faces composing the boundary $\partial T$. Moreover, for every cut cell $T \in \mathcal{T}_h^{\text{OK}} \cup \mathcal{T}_h^{\text{KO}}$, recalling that $T^\Gamma = T \cap \Gamma$, the boundary $\partial(T^i)$ of the sub-cell $T^i$ is decomposed as

$$ \partial(T^i) := (\partial T)^i \cup T^\Gamma, \quad (\partial T)^i := \partial T \cap (\overline{\Omega_i} \setminus \Gamma), \quad \forall i \in \{1, 2\}. \qquad (10) $$

To unify the notation, for every uncut cell $T \in \mathcal{T}_h^i$ with $i := \iota(T) \in \{1, 2\}$, we set (recall that $\bar{i} = 3 - i$)

$$ T^i := T, \quad T^{\bar{i}} := \emptyset, \quad (\partial T)^i := \partial T, \quad (\partial T)^{\bar{i}} := \emptyset, \quad T^{\Gamma} := \emptyset. \qquad (11) $$

Finally, for all $T \in \mathcal{T}_h$ and all $i \in \{1, 2\}$, the (possibly empty) set $\mathcal{F}_{(\partial T)^i} := \{F^i = F \cap \Omega_i | F \in \mathcal{F}_{\partial T}\}$ is the collection of the (sub-)faces composing $(\partial T)^i$.

## 3.2 Unfitted HHO spaces

HHO methods utilize polynomial unknowns attached to the mesh cells and to the mesh faces. We consider here the mixed-order HHO method, where the degree of the face unknowns is $k \ge 0$ and that of the cell unknowns is $(k+1)$. For the motivation of using a mixed-order setting, we refer the reader to [16, Remark 4.1].

For every uncut cell $T \in \mathcal{T}_h^i$, $i \in \{1, 2\}$, the local discrete HHO unknowns are a pair composed of one cell polynomial and a collection of face polynomials (one for each face of $T$), so that

$$ \hat{u}_T := (u_T, u_{\partial T}) \in \hat{\mathcal{U}}_T := \mathbb{P}^{k+1}(T; \mathbb{R}) \times \mathbb{P}^k(\mathcal{F}_{\partial T}; \mathbb{R}), \quad \mathbb{P}^k(\mathcal{F}_{\partial T}; \mathbb{R}) := \bigtimes_{F \in \mathcal{F}_{\partial T}} \mathbb{P}^k(F; \mathbb{R}). \qquad (12) $$

Here, $\mathbb{P}^{k+1}(T; \mathbb{R})$ (resp., $\mathbb{P}^k(F; \mathbb{R})$) consists of the restriction to $T$ (resp., $F$) of scalar-valued $d$ variate polynomials of degree at most $(k+1)$ (resp., $(d-1)$-variate polynomials of degree at most $k$ composed with any affine geometric mapping from the hyperplane supporting $F$ to $\mathbb{R}^{d-1}$). For every cut cell $T \in \mathcal{T}_h^{\text{OK}} \cup \mathcal{T}_h^{\text{KO}}$, following [16, 12], we double the HHO unknowns so as to have the usual unknowns available in each sub-cell, except on $T^\Gamma$ where there are no
unknowns. Thus, the local HHO unknowns in every cut cell are

$$ \hat{u}_T := (\hat{u}_{T^1}, \hat{u}_{T^2}) := (u_{T^1}, u_{(\partial T)^1}, u_{T^2}, u_{(\partial T)^2}) \in \hat{\mathcal{U}}_T := \hat{\mathcal{U}}_{T^1} \times \hat{\mathcal{U}}_{T^2}, \quad (13) $$

with $\hat{\mathcal{U}}_{T^i} := \mathbb{P}^{k+1}(T^i) \times \mathbb{P}^k(\mathcal{F}_{(\partial T)^i})$ and $\mathbb{P}^k(\mathcal{F}_{(\partial T)^i}) := \times_{F^i \in \mathcal{F}_{(\partial T)^i}} \mathbb{P}^k(F^i)$, for all $i \in \{1, 2\}$. To unify the notation between cut and uncut cells, we write $\hat{u}_T$ in the same format as in (13) for every uncut cell as well. Thus, we write $\hat{u}_T := (u_T, u_{\partial T}, 0, 0)$ for all $T \in \mathcal{T}_h^1$ and $\hat{u}_T := (0, 0, u_T, u_{\partial T})$ for all $T \in \mathcal{T}_h^2$.

The global HHO space on the unfitted mesh is

$$ \hat{\mathcal{U}}_h := \mathcal{U}_{\mathcal{T}_h^1} \times \mathcal{U}_{\mathcal{F}_h^1} \times \mathcal{U}_{\mathcal{T}_h^2} \times \mathcal{U}_{\mathcal{F}_h^2}, \quad (14) $$

where

$$ \mathcal{U}_{\mathcal{T}_h^i} := \times_{T \in \mathcal{T}_h} \mathbb{P}^{k+1}(T^i; \mathbb{R}), \quad \mathcal{U}_{\mathcal{F}_h^i} := \times_{F \in \mathcal{F}_h} \mathbb{P}^k(F^i; \mathbb{R}), \quad \forall i \in \{1, 2\}. \quad (15) $$

To enforce the homogeneous Dirichlet boundary condition on $\partial\Omega$, we consider the subspace $\hat{\mathcal{U}}_{h0} := \mathcal{U}_{\mathcal{T}_h^1} \times \mathcal{U}_{\mathcal{F}_h^1} \times \mathcal{U}_{\mathcal{T}_h^2} \times \mathcal{U}_{\mathcal{F}_h^2}$, where $\mathcal{U}_{\mathcal{T}_h^1}$ and $\mathcal{U}_{\mathcal{T}_h^2}$ are the subspaces of $\mathcal{U}_{\mathcal{F}_h^1}$ and $\mathcal{U}_{\mathcal{F}_h^2}$, respectively, where all the unknowns attached to faces located on $\partial\Omega$ are set to zero. Recall that only one of the two subdomains $\Omega_i$ touches the boundary, so that, in practice, only one of the two face spaces is modified to account for Dirichlet conditions.

## 3.3 Local reconstruction operator

The local reconstruction operator is a central tool in the devising of HHO methods. Recall that, in the fitted case, one possibility (see [1, 21]) is to consider, for all $T \in \mathcal{T}_h$, the operator $\mathbf{G}_T^k : \hat{\mathcal{U}}_T \to \mathbb{P}^k(T; \mathbb{R}^d)$ such that, for all $\hat{u}_T = (u_T, u_{\partial T}) \in \hat{\mathcal{U}}_T$ and all $\mathbf{q} \in \mathbb{P}^k(T; \mathbb{R}^d)$,

$$ (\mathbf{G}_T^k(\hat{u}_T), \mathbf{q})_T := (\nabla u_T, \mathbf{q})_T + (u_{\partial T} - u_T, \mathbf{q} \cdot \mathbf{n}_T)_{\partial T}. \quad (16) $$

In the unfitted case, the idea proposed in [12] is to reconstruct a gradient in each sub-cell $T^i$. Thus, one defines a pair of local operators $\mathbf{G}_{T^i}^k : \hat{\mathcal{U}}_T \to \mathbb{P}^k(T^i; \mathbb{R}^d)$, $i \in \{1, 2\}$, such that, for all $\hat{u}_T = (u_{T^1}, u_{(\partial T)^1}, u_{T^2}, u_{(\partial T)^2}) \in \hat{\mathcal{U}}_T$ and all $\mathbf{q} \in \mathbb{P}^k(T^i; \mathbb{R}^d)$,

$$ (\mathbf{G}_{T^i}^k(\hat{u}_T), \mathbf{q})_{T^i} := (\nabla u_{T^i}, \mathbf{q})_{T^i} + (u_{(\partial T)^i} - u_{T^i}, \mathbf{q} \cdot \mathbf{n}_T)_{(\partial T)^i} - \delta_{i1}(\llbracket u_T \rrbracket_{\Gamma}, \mathbf{q} \cdot \mathbf{n}_{\Gamma})_{T^i}, \quad (17) $$

with $\llbracket u_T \rrbracket_{\Gamma} := u_{T^1}|_{T^\Gamma} - u_{T^2}|_{T^\Gamma}$ and $\delta_{i1}$ denotes the Kronecker delta. Notice that the two gradients are reconstructed independently, but $\mathbf{G}_{T^1}^k$ depends on the cell components on both sides of the interface owing to the last term in (17). This term is only present when $i = 1$ because we assumed $\kappa_1 \le \kappa_2$ (to fix the ideas, otherwise the term is present only when $i = 2$). As discussed in [12, Section 2.5], this is instrumental to achieve robustness of the error estimate in the highly contrasted case $\kappa_1 \ll \kappa_2$.

Owing to the use of polynomial extensions in the present work, the stencil of the local gradient reconstruction operator needs to be extended. Indeed, on a given mesh cell $T \in \mathcal{T}_h$, this operator now involves the local cell and face unknowns in $T$ (as before), but also the cell and face unknowns of the ill-cut cells $S \in \mathcal{N}^{-1}(T)$ (those are the cells $S \in \mathcal{T}_h$ such that there exists $i \in \{1, 2\}$ so that $\mathcal{N}_i(S) = T$). Therefore, for all $T \in \mathcal{T}_h$, we introduce the notation

$$ \hat{u}_T^N := (\hat{u}_T, (\hat{u}_S)_{S \in \mathcal{N}^{-1}(T)}) \in \hat{\mathcal{U}}_T^N := \hat{\mathcal{U}}_T \times \times_{S \in \mathcal{N}^{-1}(T)} \hat{u}_S. \quad (18) $$
We can now define, for all $T \in \mathcal{T}_h$, the local gradient reconstruction operators $G_{T^i}^k : \hat{\mathcal{U}}_T^N \to \mathbb{P}^k(T^i; \mathbb{R}^d)$, for all $i \in \{1, 2\}$. For all $(T, i) \in \mathcal{P}_h$ and every polynomial $q \in \mathbb{P}^k(T^i; \mathbb{R}^d)$, $q^+$ denotes its extension to $\Delta(T)$ (observe that $T^i \cup \bigcup_{S \in \mathcal{N}_i^{-1}(T)} S^i \subset \Delta(T)$). Then, for all $\hat{u}_T^N \in \hat{\mathcal{U}}_T^N$,

i) If $T \in \mathcal{T}_h^{\text{uncut}}$, we set, for all $q \in \mathbb{P}^k(T^i; \mathbb{R}^d)$ with $i := \iota(T)$,

$$ (G_{T^i}^k(\hat{u}_T^N), q)_{T^i} := (\nabla u_{T^i}, q)_{T^i} + (u_{(\partial T)^i} - u_{T^i}, q \cdot n_T)_{(\partial T)^i} \\ + \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ (u_{(\partial S)^i} - u_{S^i}, q^+ \cdot n_S)_{(\partial S)^i} - \delta_{i1}(\llbracket u_S \rrbracket_{\Gamma}, q^+ \cdot n_{\Gamma})_{S^{\Gamma}} \right\}, \quad (19a) $$

$$ G_{T^i}^k(\hat{u}_T^N) := 0, \quad (19b) $$

where the summation over $S \in \mathcal{N}_i^{-1}(T)$ is void if $T$ is not in the image of $\mathcal{N}_i$;

ii) If $T \in \mathcal{T}_h^{\text{OK}}$, we set, for all $i \in \{1, 2\}$ and all $q \in \mathbb{P}^k(T^i; \mathbb{R}^d)$,

$$ (G_{T^i}^k(\hat{u}_T^N), q)_{T^i} := (\nabla u_{T^i}, q)_{T^i} + (u_{(\partial T)^i} - u_{T^i}, q \cdot n_T)_{(\partial T)^i} - \delta_{i1}(\llbracket u_T \rrbracket_{\Gamma}, q \cdot n_{\Gamma})_{T^{\Gamma}} \\ + \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ (u_{(\partial S)^i} - u_{S^i}, q^+ \cdot n_S)_{(\partial S)^i} - \delta_{i1}(\llbracket u_S \rrbracket_{\Gamma}, q^+ \cdot n_{\Gamma})_{S^{\Gamma}} \right\}; \quad (20) $$

iii) If $T \in \mathcal{T}_h^{\text{KO}}$, we set, for all $\bar{q} \in \mathbb{P}^k(T^{\bar{i}}; \mathbb{R}^d)$ with $i := \iota(T)$,

$$ G_{T^i}^k(\hat{u}_T^N) := \nabla u_{T^i}, \quad (21a) $$

$$ (G_{T^i}^k(\hat{u}_T^N), \bar{q})_{T^i} := (\nabla u_{T^i}, \bar{q})_{T^i} + (u_{(\partial T)^i} - u_{T^i}, \bar{q} \cdot n_T)_{(\partial T)^i} - \delta_{i1}(\llbracket u_T \rrbracket_{\Gamma}, \bar{q} \cdot n_{\Gamma})_{T^{\Gamma}} \\ + \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ (u_{(\partial S)^i} - u_{S^i}, \bar{q}^+ \cdot n_S)_{(\partial S)^i} - \delta_{i1}(\llbracket u_S \rrbracket_{\Gamma}, \bar{q}^+ \cdot n_{\Gamma})_{S^{\Gamma}} \right\}. \quad (21b) $$

Since the discretization method is formulated on the mesh sub-cells, we consider the collection thereof, which is conveniently described as a subset of the product set $\mathcal{T}_h \times \{1, 2\}$. Specifically, we set

$$ \mathcal{P}_h^{\text{OK}} := \{(T, \iota(T)) \mid T \in \mathcal{T}_h^{\text{uncut}}\} \cup \{(T, 1), (T, 2) \mid T \in \mathcal{T}_h^{\text{OK}}\} \cup \{(T, \bar{i}), \mid T \in \mathcal{T}_h^{\text{KO}}, i := \iota(T)\}, \quad (22a) $$

$$ \mathcal{P}_h^{\text{KO}} := \{(T, i) \mid T \in \mathcal{T}_h^{\text{KO}}, i := \iota(T)\}. \quad (22b) $$

Notice that $\mathcal{P}_h := \mathcal{P}_h^{\text{OK}} \cup \mathcal{P}_h^{\text{KO}}$ is the collection of all the mesh sub-cells. As a consequence of (9), the subset $\mathcal{P}_h^{\text{KO}}$ can be enumerated as follows:

$$ \mathcal{P}_h^{\text{KO}} = \{(S, i) \in \mathcal{P}_h \mid S \in \mathcal{N}_i^{-1}(T), (T, i) \in \mathcal{P}_h^{\text{OK}}\}. \quad (23) $$

The above notation allows us to rewrite in a synthetic form the formula for the gradient reconstruction while avoiding the proliferation of cases. Indeed, one readily sees that the formulas (19a)-(20)-(21b) are equivalent to setting, for all $(T, i) \in \mathcal{P}_h^{\text{OK}}$, all $\hat{u}_T^N \in \hat{\mathcal{U}}_T^N$ and all $q \in \mathbb{P}^k(T^i; \mathbb{R}^d)$,

$$ (G_{T^i}^k(\hat{u}_T^N), q)_{T^i} := (\nabla u_{T^i}, q)_{T^i} + (u_{(\partial T)^i} - u_{T^i}, q \cdot n_T)_{(\partial T)^i} - \delta_{i1}(\llbracket u_T \rrbracket_{\Gamma}, q \cdot n_{\Gamma})_{T^{\Gamma}} \\ + \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ (u_{(\partial S)^i} - u_{S^i}, q^+ \cdot n_S)_{(\partial S)^i} - \delta_{i1}(\llbracket u_S \rrbracket_{\Gamma}, q^+ \cdot n_{\Gamma})_{S^{\Gamma}} \right\}, \quad (24) $$
<figure><img src="image_6.png" /><figcaption>Fig. 4: Local stencil for gradient reconstruction operator</figcaption></figure>

where the summation over $S \in \mathcal{N}_i^{-1}(T)$ is void if $T$ is not in the image of $\mathcal{N}_i$. Moreover, (21a) means that, for all $(T, i) \in \mathcal{P}_h^{\text{KO}}$, we set

$$ \mathbf{G}_{T^i}^k(\hat{u}_T^N) := \nabla u_{T^i}. \qquad (25) $$

For all $(T, 1) \in \mathcal{P}_h^{\text{OK}}$, it is convenient to define the discrete lifting operator $\mathbf{L}_{T^1}^k : L^2(\Gamma) \to \mathbb{P}^k(T^1; \mathbb{R}^d)$ such that, for all $g \in L^2(\Gamma)$ and all $\mathbf{q} \in \mathbb{P}^k(T^1; \mathbb{R}^d)$,

$$ (\mathbf{L}_{T^1}^k(g), \mathbf{q})_{T^1} := (g, \mathbf{q} \cdot \mathbf{n}_{\Gamma})_{T^1} + \sum_{S \in \mathcal{N}_1^{-1}(T)} (g, \mathbf{q}^+ \cdot \mathbf{n}_{\Gamma})_{S^1}. \qquad (26) $$

It is also convenient to define on $\Gamma$ the function given by the jump of the cell components of a discrete HHO variable. Specifically, we define the function $[u_{\mathcal{T}_h}]_{\Gamma}$ on $\Gamma$ such that, for all $T \in \mathcal{T}_h^{\text{cut}}$, $[u_{\mathcal{T}_h}]_{\Gamma} := u_{T^1}|_{T^{\Gamma}} - u_{T^2}|_{T^{\Gamma}}$. Then, for all $(T, i) \in \mathcal{P}_h^{\text{OK}}$, we can rewrite (24) as follows:

$$ (\mathbf{G}_{T^i}^k(\hat{u}_T^N) + \delta_{i1} \mathbf{L}_{T^1}^k([u_{\mathcal{T}_h}]_{\Gamma}), \mathbf{q})_{T^i} := (\nabla u_{T^i}, \mathbf{q})_{T^i} + (u_{(\partial T)^i} - u_{T^i}, \mathbf{q} \cdot \mathbf{n}_T)_{(\partial T)^i} \\ + \sum_{S \in \mathcal{N}_i^{-1}(T)} (u_{(\partial S)^i} - u_{S^i}, \mathbf{q}^+ \cdot \mathbf{n}_S)_{(\partial S)^i}. \qquad (27) $$

## 3.4 Stabilization

We consider three local stabilization bilinear forms. The first one is the usual HHO stabilization bilinear form in the mixed-order setting (in the spirit of the HDG Lehrenfeld-Schöberl stabilization [29]). Its role is to weakly enforce the matching between cell- and face-based HHO unknowns on the sub-faces of the mesh sub-cells in each sub-domain. For all $\hat{v}_h, \hat{w}_h \in \hat{\mathcal{U}}_{h0}$, we set

$$ s_h^\circ(\hat{v}_h, \hat{w}_h) := \sum_{(T,i) \in \mathcal{P}_h} \kappa_i h_T^{-1} (\Pi_{(\partial T)^i}^k(v_{T^i}) - v_{(\partial T)^i}), \Pi_{(\partial T)^i}^k(w_{T^i}) - w_{(\partial T)^i})_{(\partial T)^i}, \qquad (28) $$

where $\Pi_{(\partial T)^i}^k$ denotes the $L^2$-orthogonal projector onto $\mathbb{P}^k(\mathcal{F}_{(\partial T)^i}; \mathbb{R})$.
The second stabilization bilinear form is the same as the one introduced in [16, 12] to control the jumps across the interface $\Gamma$. Specifically, we set, for all $\hat{v}_h, \hat{w}_h \in \hat{\mathcal{U}}_{h0}$,

$$s_h^{\Gamma}(\hat{v}_h, \hat{w}_h) := \sum_{T \in \mathcal{T}_h^{\text{cut}}} \kappa_1 h_T^{-1}(\llbracket v_T \rrbracket_{\Gamma}, \llbracket w_T \rrbracket_{\Gamma})_{T^{\Gamma}}. \quad (29)$$

The third stabilization bilinear form is specific to the present setting using polynomial extensions. Its aim is to provide some control on the cell components in ill-cut cells and is devised in the same spirit as the so-called direct ghost penalty method from [31, 28]. Specifically, we set, for all $\hat{v}_h, \hat{w}_h \in \hat{\mathcal{U}}_{h0}$,

$$s_h^{\mathcal{N}}(\hat{v}_h, \hat{w}_h) := \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \sum_{S \in \mathcal{N}_i^{-1}(T)} \eta_{\mathcal{N}} \kappa_i h_T^{-2} (v_{S^i} - v_{T^i}^+, w_{S^i} - w_{T^i}^+)_{T^i}, \quad (30)$$

recalling that the superscript $^+$ denotes the extension of a polynomial originally defined on $T^i$ to $\Delta(T)$, and $\eta_{\mathcal{N}} > 0$ is a user-defined parameter.

Putting everything together, we define, for all $\hat{v}_h, \hat{w}_h \in \hat{\mathcal{U}}_{h0}$,

$$s_h(\hat{v}_h, \hat{w}_h) := s_h^{\circ}(\hat{v}_h, \hat{w}_h) + s_h^{\Gamma}(\hat{v}_h, \hat{w}_h) + s_h^{\mathcal{N}}(\hat{v}_h, \hat{w}_h). \quad (31)$$

## 3.5 Discrete problem

The global discrete problem reads as follows: Find $\hat{u}_h \in \hat{\mathcal{U}}_{h0}$ such that,

$$a_h(\hat{u}_h, \hat{w}_h) = \ell_h(\hat{w}_h) \quad \forall \hat{w}_h \in \hat{\mathcal{U}}_{h0}, \quad (32)$$

where the bilinear $a_h$ is defined as follows: For all $\hat{v}_h, \hat{w}_h \in \hat{\mathcal{U}}_{h0}$,

$$a_h(\hat{v}_h, \hat{w}_h) := b_h(\hat{v}_h, \hat{w}_h) + s_h(\hat{v}_h, \hat{w}_h), \quad (33)$$

with

$$b_h(\hat{v}_h, \hat{w}_h) := \sum_{(T,i) \in \mathcal{P}_h} \kappa_i (\mathbf{G}_{T^i}^k(\hat{v}_T^{\mathcal{N}}), \mathbf{G}_{T^i}^k(\hat{w}_T^{\mathcal{N}}))_{T^i}. \quad (34)$$

Moreover, the linear form $\ell_h$ is such that, for all $\hat{w}_h \in \hat{\mathcal{U}}_{h0}$,

$$\ell_h(\hat{w}_h) := \sum_{(T,i) \in \mathcal{P}_h} (f, w_{T^i})_{T^i} + \sum_{T \in \mathcal{T}_h^{\text{cut}}} \left\{ (g_N, w_{T^2})_{T^\Gamma} + \kappa_1 h_T^{-1} (g_D, \llbracket w_T \rrbracket_\Gamma)_{T^\Gamma} \right\} \\ - \sum_{(T,1) \in \mathcal{P}_h^{\text{OK}}} \kappa_1 (\mathbf{L}_{T^1}^k(g_D), \mathbf{G}_{T^1}^k(\hat{w}_T^{\mathcal{N}}))_{T^1}, \quad (35)$$

where the lifting operator $\mathbf{L}_{T^1}^k$ is defined in (26). Notice that, in the discrete problem (32), only the Dirichlet boundary condition is enforced explicitly (on the face unknowns located on the boundary $\partial\Omega$), whereas both jump conditions across the interface, (1b) and (1c), are enforced weakly. We notice that cell unknowns can still be eliminated from (32) by static condensation. Whenever $\mathcal{N}^{-1}(T)$ is nonempty for a mesh cell $T \in \mathcal{T}_h$, all the cell unknowns attached to $T \cup \bigcup_{S \in \mathcal{N}^{-1}(T)} S$ are blocked together.

To prepare for the consistency error analysis, we state a useful result on the discrete bilinear and linear forms. Although this lemma is only used in Section 4, it is given here so as to highlight why the definitions (24)-(25) of the reconstructed gradient are meaningful.
**Lemma 3.1** (Consistency: preparatory identity). Let $u^{\text{ex}}$ be the weak solution to (2). Let $\hat{v}_h \in \hat{\mathcal{U}}_{h0}$ be arbitrary and set

$$ \delta_{T,i} := \mathbf{G}_{T,i}^k(\hat{v}_T^N) + \delta_{i1} \mathbf{L}_{T1}^k(\llbracket u^{\text{ex}} \rrbracket_{\Gamma}) - \nabla u_i^{\text{ex}}|_{T^i}, \quad \forall (T,i) \in \mathcal{P}_h^{\text{OK}}, \quad (36a) $$

$$ \delta_{T,S,i} := \{\mathbf{G}_{T,i}^k(\hat{v}_T^N) + \delta_{i1} \mathbf{L}_{T1}^k(\llbracket u^{\text{ex}} \rrbracket_{\Gamma})\}|_{S^i} - \nabla u_i^{\text{ex}}|_{S^i}, \quad \forall (T,i) \in \mathcal{P}_h^{\text{OK}}, \forall S \in \mathcal{N}_i^{-1}(T), \quad (36b) $$

$$ \delta_{T,i} := \mathbf{G}_{T,i}^k(\hat{v}_T^N) - \nabla u_i^{\text{ex}}|_{T^i} = \nabla v_{T^i} - \nabla u_i^{\text{ex}}|_{T^i}, \quad \forall (T,i) \in \mathcal{P}_h^{\text{KO}}. \quad (36c) $$

Then, the following holds for all $\hat{w}_h \in \hat{\mathcal{U}}_{h0}$,

$$ a_h(\hat{v}_h, \hat{w}_h) - \ell_h(\hat{w}_h) = \sum_{(T,i) \in \mathcal{P}_h} \kappa_i(\delta_{T,i}, \nabla w_{T^i})_{T^i} \\ + \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i \left\{ (\delta_{T,i} \cdot \mathbf{n}_T, w_{(\partial T)^i} - w_{T^i})_{(\partial T)^i} - \delta_{i1}(\delta_{T,1} \cdot \mathbf{n}_\Gamma, \llbracket w_T \rrbracket_\Gamma)_{T^\Gamma} \right. \\ + \sum_{S \in \mathcal{N}_i^{-1}(T)} \left. \{(\delta_{T,S,i} \cdot \mathbf{n}_S, w_{(\partial S)^i} - w_{S^i})_{(\partial S)^i} - \delta_{i1}(\delta_{T,S,1} \cdot \mathbf{n}_\Gamma, \llbracket w_S \rrbracket_\Gamma)_{S^\Gamma}\} \right\} \\ + s_h^\circ(\hat{v}_h, \hat{w}_h) + s_h^N(\hat{v}_h, \hat{w}_h) + \sum_{T \in \mathcal{T}_h^{\text{cut}}} \kappa_1 h_T^{-1}(\llbracket v_T - u^{\text{ex}} \rrbracket_\Gamma, \llbracket w_T \rrbracket_\Gamma)_{T^\Gamma}. \quad (37) $$

*Proof.* (1) Recalling the properties (1a)-(1c) satisfied by the exact solution and the definition (35) of the linear form $\ell_h$, we infer that, for all $\hat{w}_h \in \hat{\mathcal{U}}_{h0}$, $\ell_h(\hat{w}_h) = L_1 + L_2 + L_3$ with

$$ L_1 := \sum_{(T,i) \in \mathcal{P}_h} -(\nabla \cdot (\kappa_i \nabla u_i^{\text{ex}}), w_{T^i})_{T^i} + \sum_{T \in \mathcal{T}_h^{\text{cut}}} ((\kappa_1 \nabla u_1^{\text{ex}} - \kappa_2 \nabla u_2^{\text{ex}}) \cdot \mathbf{n}_\Gamma, w_{T^2})_{T^\Gamma}, $$

$$ L_2 := \sum_{T \in \mathcal{T}_h^{\text{cut}}} \kappa_1 h_T^{-1}(\llbracket u^{\text{ex}} \rrbracket_\Gamma, \llbracket w_T \rrbracket_\Gamma)_{T^\Gamma}, \quad L_3 := - \sum_{(T,1) \in \mathcal{P}_h^{\text{OK}}} \kappa_1 (\mathbf{L}_{T1}^k(\llbracket u^{\text{ex}} \rrbracket_\Gamma), \mathbf{G}_{T1}^k(\hat{v}_T^N))_{T^1}. $$

Integrating by parts and re-organizing the boundary term on $\Gamma$ (recall that $\mathbf{n}_{T^1} = -\mathbf{n}_{T^2} = \mathbf{n}_\Gamma$) gives

$$ \sum_{(T,i) \in \mathcal{P}_h} -(\nabla \cdot (\kappa_i \nabla u_i^{\text{ex}}), w_{T^i})_{T^i} = \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \{(\nabla u_i^{\text{ex}}, \nabla w_{T^i})_{T^i} - (\nabla u_i^{\text{ex}} \cdot \mathbf{n}_T, w_{T^i})_{(\partial T)^i}\} \\ - \sum_{T \in \mathcal{T}_h^{\text{cut}}} \{(\kappa_1 \nabla u_1^{\text{ex}} \cdot \mathbf{n}_\Gamma, w_{T^1})_{T^\Gamma} - (\kappa_2 \nabla u_2^{\text{ex}} \cdot \mathbf{n}_\Gamma, w_{T^2})_{T^\Gamma}\}. $$

Since $\nabla u_i^{\text{ex}}$ is single-valued on $(\partial T)^i \cap \Omega_i$ and $w_{(\partial T)^i}$ vanishes on the mesh boundary faces, we obtain

$$ \sum_{(T,i) \in \mathcal{P}_h} -(\nabla \cdot (\kappa_i \nabla u_i^{\text{ex}}), w_{T^i})_{T^i} = \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \{(\nabla u_i^{\text{ex}}, \nabla w_{T^i})_{T^i} + (\nabla u_i^{\text{ex}} \cdot \mathbf{n}_T, w_{(\partial T)^i} - w_{T^i})_{(\partial T)^i}\} \\ - \sum_{T \in \mathcal{T}_h^{\text{cut}}} \{(\kappa_1 \nabla u_1^{\text{ex}} \cdot \mathbf{n}_\Gamma, w_{T^1})_{T^\Gamma} - (\kappa_2 \nabla u_2^{\text{ex}} \cdot \mathbf{n}_\Gamma, w_{T^2})_{T^\Gamma}\}. $$

Therefore, we have

$$ L_1 = \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \{(\nabla u_i^{\text{ex}}, \nabla w_{T^i})_{T^i} + (\nabla u_i^{\text{ex}} \cdot \mathbf{n}_T, w_{(\partial T)^i} - w_{T^i})_{(\partial T)^i}\} - \sum_{T \in \mathcal{T}_h^{\text{cut}}} \kappa_1 (\nabla u_1^{\text{ex}} \cdot \mathbf{n}_\Gamma, \llbracket w_T \rrbracket_\Gamma)_{T^\Gamma}. $$

A rewriting of the last term on the right-hand side gives

$$ L_1 = \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \{(\nabla u_i^{\text{ex}}, \nabla w_{T^i})_{T^i} + (\nabla u_i^{\text{ex}} \cdot \mathbf{n}_T, w_{(\partial T)^i} - w_{T^i})_{(\partial T)^i} - \delta_{i1}(\nabla u_1^{\text{ex}} \cdot \mathbf{n}_\Gamma, \llbracket w_T \rrbracket_\Gamma)_{T^\Gamma}\}. $$
We now use the enumeration formula (23) for the pairs $(T, i) \in \mathcal{P}_h^{\text{KO}}$. This gives

$$L_1 = \sum_{(T,i) \in \mathcal{P}_h} \kappa_i (\nabla u_i^{\text{ex}}, \nabla w_{T^i})_{T^i} \\ + \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i \left\{ (\nabla u_i^{\text{ex}} \cdot \mathbf{n}_T, w_{(\partial T)^i} - w_{T^i})_{(\partial T)^i} - \delta_{i1} (\nabla u_1^{\text{ex}} \cdot \mathbf{n}_{\Gamma}, \llbracket w_T \rrbracket_{\Gamma})_{T^i} \right. \\ \left. + \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ (\nabla u_i^{\text{ex}} \cdot \mathbf{n}_S, w_{(\partial S)^i} - w_{S^i})_{(\partial S)^i} - \delta_{i1} (\nabla u_1^{\text{ex}} \cdot \mathbf{n}_{\Gamma}, \llbracket w_S \rrbracket_{\Gamma})_{S^i} \right\} \right\}.$$

(2) For all $\hat{v}_h \in \hat{\mathcal{U}}_{h0}$, setting $\gamma_{T,i} := \mathbf{G}_{T^i}^k(\hat{v}_T^N) + \delta_{i1}\mathbf{L}_{T^1}^k(\llbracket u^{\text{ex}} \rrbracket_{\Gamma})$ for all $(T, i) \in \mathcal{P}_h^{\text{OK}}$ and $\gamma_{T,i} := \mathbf{G}_{T^i}^k(\hat{v}_T^N) = \nabla v_{T^i}$ for all $(T, i) \in \mathcal{P}_h^{\text{KO}}$, and using the definition (34) of the bilinear form $b_h$ together with the definition of $L_3$ given above, we infer that

$$b_h(\hat{v}_h, \hat{w}_h) - L_3 = \sum_{(T,i) \in \mathcal{P}_h} \kappa_i (\mathbf{G}_{T^i}^k(\hat{v}_T^N), \mathbf{G}_{T^i}^k(\hat{w}_T^N))_{T^i} + \sum_{(T,1) \in \mathcal{P}_h^{\text{OK}}} \kappa_1 (\mathbf{L}_{T^1}^k(\llbracket u^{\text{ex}} \rrbracket_{\Gamma}), \mathbf{G}_{T^1}^k(\hat{w}_T^N))_{T^1} \\ = \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i (\gamma_{T,i}, \mathbf{G}_{T^i}^k(\hat{w}_T^N))_{T^i} + \sum_{(T,i) \in \mathcal{P}_h^{\text{KO}}} \kappa_i (\gamma_{T,i}, \mathbf{G}_{T^i}^k(\hat{w}_T^N))_{T^i}.$$

Using the definitions (24)-(25) of the gradient reconstruction operator and since $\gamma_{T,i} \in \mathbb{P}^k(T^i; \mathbb{R}^d)$ for all $(T, i) \in \mathcal{P}_h$, we obtain

$$b_h(\hat{v}_h, \hat{w}_h) = \sum_{(T,i) \in \mathcal{P}_h} \kappa_i(\gamma_{T,i}, \nabla w_{T^i})_{T^i} \\ + \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i \left\{ (\gamma_{T,i} \cdot \mathbf{n}_T, w_{(\partial T)^i} - w_{T^i})_{(\partial T)^i} - \delta_{i1} (\gamma_{T,1} \cdot \mathbf{n}_{\Gamma}, \llbracket w_T \rrbracket_{\Gamma})_{T^i} \right. \\ \left. + \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ (\gamma_{T,i}^{+} \cdot \mathbf{n}_S, w_{(\partial S)^i} - w_{S^i})_{(\partial S)^i} - \delta_{i1} (\gamma_{T,1}^{+} \cdot \mathbf{n}_{\Gamma}, \llbracket w_S \rrbracket_{\Gamma})_{S^i} \right\} \right\}.$$

(3) Combining the identities from Steps (1) and (2) and recalling that $a_h = b_h + s_h$ completes the proof. $\square$

# 4 Stability, consistency, and error estimate

In this section, we perform the stability and error analysis for the unfitted HHO method introduced in the previous section. We use the convention $A \lesssim B$ to abbreviate the inequality $A \le CB$ for positive real numbers $A$ and $B$, where the constant $C$ only depends on the polynomial degree, the mesh shape-regularity, the parameter $\vartheta$ used in the definition of the pairing operator, and the space dimension.

## 4.1 Analysis tools

In this section, we present the main tools to perform the error analysis. The proofs are postponed to Section 6. For all $(T, i) \in \mathcal{P}_h^{\text{OK}}$, recall that the superscript $+$ is used here to indicate the extension of a polynomial originally defined on $T^i$ to $\Delta(T)$ (observe that $T \cup \bigcup_{S \in \mathcal{N}_i^{-1}(T)} S \subset \Delta(T)$).
**Lemma 4.1** (Discrete inverse inequalities). The following holds for all $(T, i) \in \mathcal{P}_h^{\text{OK}}$, all $\phi \in \mathbb{P}^{\ell}(T^i; \mathbb{R})$, and all $\ell \ge 0$,

$$ \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \left\{ \|\phi^+\|_S + h_S^{\frac{1}{2}} \|\phi^+\|_{(\partial S)^i \cup S^\Gamma} \right\} \lesssim \|\phi\|_{T^i}, \quad (38a) $$

$$ \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} h_S^{-\frac{1}{2}} \|(I - \Pi_{(\partial S)^i}^k)(\phi^+)\|_{(\partial S)^i} \lesssim \|\nabla \phi\|_{T^i}. \quad (38b) $$

For all $s \ge 0$ and all $i \in \{1, 2\}$, let $E_i^s: H^s(\Omega_i) \to H^s(\mathbb{R}^d)$ be a stable extension operator. For all $v \in H^s(\Omega_1 \cup \Omega_2)$ with $s > \frac{1}{2}$ and all $(T, i) \in \mathcal{P}_h^{\text{OK}}$, we define

$$ I_{T^i}^{k+1}(v_i) := \Pi_T^{k+1}(E_i^s(v_i))|_{T^i} \in \mathbb{P}^{k+1}(T^i; \mathbb{R}), \quad (39) $$

where $\Pi_T^{k+1}$ denotes the $L^2$-orthogonal projection onto $\mathbb{P}^{k+1}(T; \mathbb{R})$. (Notice that the operator $I_{T^i}^{k+1}$ depends on the Sobolev index $s$, but this dependency is not tracked to simplify the notation; notice also that the extension operator is not needed if $T$ is an uncut cell.) Let us set, for all $(T, i) \in \mathcal{P}_h^{\text{OK}}$,

$$ \epsilon_{T,i}(v_i) := \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \left\{ \begin{aligned} & \|v_i - I_{T^i}^{k+1}(v_i)^+\|_{S^i} + h_S^{\frac{1}{2}} \|v_i - I_{T^i}^{k+1}(v_i)^+\|_{(\partial S)^i} \\ & + h_S \|\nabla(v_i - I_{T^i}^{k+1}(v_i)^+)\|_{S^i} + h_S^{\frac{3}{2}} \|\nabla(v_i - I_{T^i}^{k+1}(v_i)^+)\|_{(\partial S)^i} \end{aligned} \right\}, \quad (40a) $$

$$ \epsilon_{T,i}^{\Gamma}(v) := \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} h_S^{\frac{1}{2}} \| [v - I_T^{k+1}(v)]^\Gamma \|_{S^\Gamma}, \quad (40b) $$

where $[v - I_T^{k+1}(v)]^\Gamma|_{S^\Gamma} := (v_1 - I_{T^1}^{k+1}(v_1))|_{S^\Gamma} - (v_2 - I_{T^2}^{k+1}(v_2))|_{S^\Gamma}$.

**Lemma 4.2** (Approximation). For all $v \in H^s(\Omega_1 \cup \Omega_2)$ with $s \in (\frac{3}{2}, k+2]$ and all $(T, i) \in \mathcal{P}_h^{\text{OK}}$, we have

$$ \epsilon_{T,i}(v_i) \lesssim h_T^s |E_i^s(v_i)|_{H^s(\Delta(T))}, \quad (41a) $$

$$ \epsilon_{T,i}^{\Gamma}(v) \lesssim h_T^s \sum_{i \in \{1,2\}} |E_i^s(v_i)|_{H^s(\Delta_2(T))}. \quad (41b) $$

## 4.2 Stability and well-posedness

We introduce the following norm on $\hat{\mathcal{U}}_{h0}$: For all $\hat{v}_h \in \hat{\mathcal{U}}_{h0}$,

$$ \|\hat{v}_h\|_{h0}^2 := \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \|\nabla v_{T^i}\|_{T^i}^2 + |\hat{v}_h|_S^2 + s_h^\Gamma(\hat{v}_h, \hat{v}_h) + s_h^\mathcal{N}(\hat{v}_h, \hat{v}_h), \quad (42) $$

where

$$ |\hat{v}_h|_S^2 := \sum_{(T,i) \in \mathcal{P}_h} \kappa_i h_T^{-1} \|v_{(\partial T)^i} - v_{T^i}\|_{(\partial T)^i}^2, \quad (43) $$

and the stabilization bilinear forms are defined in Section 3.4. We observe that in general $|\hat{v}_h|_S^2 \neq s_h^\circ(\hat{v}_h, \hat{v}_h)$ since the cell component is projected when evaluating the latter. Furthermore,
we notice the following rewriting:

$$
\begin{equation}
\begin{aligned}
\|\hat{v}_h\|_{h0}^2 = & \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \left\{ \|\nabla v_{T^i}\|_{T^i}^2 + h_T^{-1} \|v_{(\partial T)^i} - v_{T^i}\|_{(\partial T)^i}^2 + \delta_{i1} h_T^{-1} \|\llbracket v_T \rrbracket\|_{T^i}^2 \right\} \\
& + \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \sum_{S \in \mathcal{N}_i^{-1}(T)} \eta_N \kappa_i h_T^{-2} \|v_{S^i} - v_{T^i}\|_{T^i}^2.
\end{aligned}
\tag{44}
\end{equation}
$$

It is straightforward to verify that (42) defines a norm on $\hat{\mathcal{U}}_{h0}$. Indeed, if $\hat{v}_h \in \hat{\mathcal{U}}_{h0}$ satisfies $\|\hat{v}_h\|_{h0} = 0$, then all the cell unknowns and all the face unknowns are constant and take the same value inside each sub-domain $\Omega_i$ and globally in the domain $\Omega$. Since the face unknowns on the boundary $\partial\Omega$ vanish, we conclude that all the components of $\hat{v}_h$ are zero.

**Lemma 4.3** (Stability and boundedness). We have, for all $\hat{v}_h \in \hat{\mathcal{U}}_{h0}$,

$$ \|\hat{v}_h\|_{h0}^2 \lesssim a_h(\hat{v}_h, \hat{v}_h) \lesssim \|\hat{v}_h\|_{h0}^2. \tag{45} $$

*Proof.* (1) Bounds on reconstructed gradient. Let $(T, i) \in \mathcal{P}_h^{\text{OK}}$. Then, taking $\boldsymbol{q} := \nabla v_{T^i} \in \mathbb{P}^k(T^i; \mathbb{R}^d)$ in (24), we obtain

$$
\begin{aligned}
\|\nabla v_{T^i}\|_{T^i}^2 = & (\mathbf{G}_{T^i}^k(\hat{v}_T^N), \nabla v_{T^i})_{T^i} - (v_{(\partial T)^i} - v_{T^i}, \nabla v_{T^i} \cdot \mathbf{n}_T)_{(\partial T)^i} + \delta_{i1}(\llbracket v_T \rrbracket_{\Gamma}, \nabla v_{T^i} \cdot \mathbf{n}_{\Gamma})_{T^i} \\
& - \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ (v_{(\partial S)^i} - v_{S^i}, (\nabla v_{T^i})^+ \cdot \mathbf{n}_S)_{(\partial S)^i} - \delta_{i1}(\llbracket v_S \rrbracket_{\Gamma}, (\nabla v_{T^i})^+ \cdot \mathbf{n}_{\Gamma})_{S^i} \right\}.
\end{aligned}
$$

Since $\nabla v_{T^i} \cdot \mathbf{n}_T \in \mathbb{P}^k(\mathcal{F}_{(\partial T)^i}; \mathbb{R})$ and $(\nabla v_{T^i})^+ \cdot \mathbf{n}_S \in \mathbb{P}^k(\mathcal{F}_{(\partial S)^i}; \mathbb{R})$, we infer that

$$
\begin{aligned}
\|\nabla v_{T^i}\|_{T^i}^2 = & (\mathbf{G}_{T^i}^k(\hat{v}_T^N), \nabla v_{T^i})_{T^i} - (v_{(\partial T)^i} - \Pi_{(\partial T)^i}^k(v_{T^i}), \nabla v_{T^i} \cdot \mathbf{n}_T)_{(\partial T)^i} + \delta_{i1}(\llbracket v_T \rrbracket_{\Gamma}, \nabla v_{T^i} \cdot \mathbf{n}_{\Gamma})_{T^i} \\
& - \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ (v_{(\partial S)^i} - \Pi_{(\partial S)^i}^k(v_{S^i}), (\nabla v_{T^i})^+ \cdot \mathbf{n}_S)_{(\partial S)^i} - \delta_{i1}(\llbracket v_S \rrbracket_{\Gamma}, (\nabla v_{T^i})^+ \cdot \mathbf{n}_{\Gamma})_{S^i} \right\}.
\end{aligned}
$$

Invoking the Cauchy–Schwarz inequality and the discrete trace inequality (38a), we obtain

$$
\|\nabla v_{T^i}\|_{T^i}^2 \lesssim \|\mathbf{G}_{T^i}^k(\hat{v}_T^N)\|_{T^i}^2 + \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \left\{ h_S^{-1} \|v_{(\partial S)^i} - \Pi_{(\partial S)^i}^k(v_{S^i})\|_{(\partial S)^i}^2 + \delta_{i1} h_S^{-1} \|\llbracket v_S \rrbracket\|_{S^i}^2 \right\}.
$$

Proceeding similarly proves that

$$
\|\mathbf{G}_{T^i}^k(\hat{v}_T^N)\|_{T^i}^2 \lesssim \|\nabla v_{T^i}\|_{T^i}^2 + \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \left\{ h_S^{-1} \|v_{(\partial S)^i} - \Pi_{(\partial S)^i}^k(v_{S^i})\|_{(\partial S)^i}^2 + \delta_{i1} h_S^{-1} \|\llbracket v_S \rrbracket\|_{S^i}^2 \right\}.
$$

Moreover, for all $(T, i) \in \mathcal{P}_h^{\text{KO}}$, we have $\mathbf{G}_{T^i}^k(\hat{v}_T^N) = \nabla v_{T^i}$ in $T^i$. Therefore, multiplying by $\kappa_i$ and summing over all $(T, i) \in \mathcal{P}_h$, we infer that

$$ \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \|\nabla v_{T^i}\|_{T^i}^2 \lesssim b_h(\hat{v}_h, \hat{v}_h) + s_h^\circ(\hat{v}_h, \hat{v}_h) + s_h^\Gamma(\hat{v}_h, \hat{v}_h), \tag{46a} $$

$$ b_h(\hat{v}_h, \hat{v}_h) \lesssim \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \|\nabla v_{T^i}\|_{T^i}^2 + s_h^\circ(\hat{v}_h, \hat{v}_h) + s_h^\Gamma(\hat{v}_h, \hat{v}_h). \tag{46b} $$

Since $s_h^\circ(\hat{v}_h, \hat{v}_h) \le |\hat{v}_h|_S^2$, the upper bound in (45) readily follows from (46b).

(2) To prove the lower bound in (45), it remains to estimate $|\hat{v}_h|_S^2$. We first observe that

$$ |\hat{v}_h|_S^2 = \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} h_S^{-1} \|v_{(\partial S)^i} - v_{S^i}\|_{(\partial S)^i}^2. $$
Consider first the case $S = T$. The triangle inequality gives

$$ \|v_{(\partial T)^i} - v_{T^i}\|_{(\partial T)^i} \le \|v_{(\partial T)^i} - \Pi_{(\partial T)^i}^k(v_{T^i})\|_{(\partial T)^i} + \|(I - \Pi_{(\partial T)^i}^k)(v_{T^i})\|_{(\partial T)^i}, $$

and owing to (38b), we infer that

$$ \|v_{(\partial T)^i} - v_{T^i}\|_{(\partial T)^i} \lesssim \|v_{(\partial T)^i} - \Pi_{(\partial T)^i}^k(v_{T^i})\|_{(\partial T)^i} + h_T^{\frac{1}{2}} \|\nabla v_{T^i}\|_{T^i}. \quad (47) $$

Let now $S \in \mathcal{N}_i^{-1}(T)$. The triangle inequality implies that

$$ \|v_{(\partial S)^i} - v_{S^i}\|_{(\partial S)^i} \le \|v_{(\partial S)^i} - \Pi_{(\partial S)^i}^k(v_{S^i})\|_{(\partial S)^i} + \|(I - \Pi_{(\partial S)^i}^k)(v_{S^i} - v_{T^i}^+)\|_{(\partial S)^i} \\ + \|(I - \Pi_{(\partial S)^i}^k)(v_{T^i}^+)\|_{(\partial S)^i}. $$

For the second term on the right-hand side, we simply notice that $\|(I - \Pi_{(\partial S)^i}^k)(v_{S^i} - v_{T^i}^+)\|_{(\partial S)^i} \le \|v_{S^i} - v_{T^i}^+\|_{(\partial S)^i}$, whereas we invoke (38b) to bound the third term. Altogether, this gives

$$ \|v_{(\partial S)^i} - v_{S^i}\|_{(\partial S)^i} \lesssim \|v_{(\partial S)^i} - \Pi_{(\partial S)^i}^k(v_{S^i})\|_{(\partial S)^i} + \|v_{S^i} - v_{T^i}^+\|_{(\partial S)^i} + h_S^{\frac{1}{2}} \|\nabla v_{T^i}\|_{T^i}. $$

The first term on the right-hand side is bounded using the $L^2$-stability of $\Pi_{(\partial S)^i}^k$ and (47). The second term is controlled by means of the discrete trace inequality (38a) and the definition (30) of the stabilization bilinear from $s_h^N$. Squaring the resulting inequality, multiplying by $\kappa_i h_S^{-1}$, summing over $S \in \mathcal{N}_i^{-1}(T)$, and finally summing over $(T, i) \in \mathcal{P}_h^{\text{OK}}$, we obtain

$$ |\hat{v}_h|_S^2 \lesssim s_h^\circ(\hat{v}_h, \hat{v}_h) + s_h^N(\hat{v}_h, \hat{v}_h) + \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \|\nabla v_{T^i}\|_{T^i}^2. $$

Combining this bound with (46a) gives

$$ \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \|\nabla v_{T^i}\|_{T^i}^2 + |\hat{v}_h|_S^2 \lesssim b_h(\hat{v}_h, \hat{v}_h) + s_h^\circ(\hat{v}_h, \hat{v}_h) + s_h^\Gamma(\hat{v}_h, \hat{v}_h) + s_h^N(\hat{v}_h, \hat{v}_h) = a_h(\hat{v}_h, \hat{v}_h), $$

whence the lower bound in (45) readily follows. $\square$

## 4.3 Interpolation operator

For all $T \in \mathcal{T}_h$, the local interpolation operator is defined as follows:

$$ \hat{I}_T^k(v) := (\hat{I}_{T^1}^k(v), \hat{I}_{T^2}^k(v)) := (I_{T^1}^{k+1}(v_1), \Pi_{(\partial T)^1}^k(v_1), I_{T^2}^{k+1}(v_2), \Pi_{(\partial T)^2}^k(v_2)) \in \hat{\mathcal{U}}_T, \quad (48) $$

where $I_{T^i}^{k+1}(v_i)$ is defined in (39) for all $(T, i) \in \mathcal{P}_h^{\text{OK}}$, whereas for all $(S, i) \in \mathcal{P}_h^{\text{KO}}$, we set

$$ I_{S^i}^{k+1}(v_i) := I_{T^i}^{k+1}(v_i)^+|_{S^i}, \quad T := \mathcal{N}_i(S), \quad (49) $$

where we recall that the superscript $+$ is used to indicate the extension of a polynomial originally defined on $T^i$ to $\Delta(T)$.

The global interpolation operator is denoted by $\hat{I}_h^k : H^s(\Omega_1 \cup \Omega_2) \to \hat{\mathcal{U}}_h$ and is such that the local components of $\hat{I}_h^k(v)$ on a mesh cell $T \in \mathcal{T}_h$ are those given by $\hat{I}_T^k(v|_T)$. Notice that $\hat{I}_h^k(u) \in \hat{\mathcal{U}}_{h0}$ for the exact solution since $u|_{\partial\Omega} = 0$.
**Lemma 4.4 (Approximation).** For all $v \in H^s(\Omega_1 \cup \Omega_2)$, $s > \frac{3}{2}$, and all $(T, i) \in \mathcal{P}_h^{\text{OK}}$, set

$$ \delta_{T,i}(v) := \mathbf{G}_{T,i}^k (\hat{I}_T^k(v)^\mathcal{N}) + \delta_{i1} \mathbf{L}_{T1}^k (\llbracket v \rrbracket_\Gamma) - \nabla v_i |_{T^i}, \quad (50a) $$

$$ \delta_{T,S,i}(v) := \{ \mathbf{G}_{T,i}^k (\hat{I}_T^k(v)^\mathcal{N}) + \delta_{i1} \mathbf{L}_{T1}^k (\llbracket v \rrbracket_\Gamma) \}^+ |_{S^i} - \nabla v_i |_{S^i}, \quad \forall S \in \mathcal{N}_i^{-1}(T), \quad (50b) $$

and recall the error measures $\epsilon_{T,i}(v_i)$ and $\epsilon_{T,i}^\Gamma(v)$ defined in (40a)-(40b). Then, the following holds:

$$ \|\delta_{T,i}(v)\|_{T^i} + h_T^{\frac{1}{2}} \|\delta_{T,i}(v)\|_{(\partial T)^i \cup T^\Gamma} \lesssim h_T^{-1} \{\epsilon_{T,i}(v_i) + \delta_{i1} \epsilon_{T,i}^\Gamma(v)\}, \quad (51a) $$

$$ \sum_{S \in \mathcal{N}_i^{-1}(T)} \{ \|\delta_{T,S,i}(v)\|_{S^i} + h_S^{\frac{1}{2}} \|\delta_{T,S,i}(v)\|_{(\partial S)^i \cup S^\Gamma} \} \lesssim h_T^{-1} \{\epsilon_{T,i}(v_i) + \delta_{i1} \epsilon_{T,i}^\Gamma(v)\}. \quad (51b) $$

*Proof.* Let $v \in H^s(\Omega_1 \cup \Omega_2)$, $s > \frac{3}{2}$, and let $(T, i) \in \mathcal{P}_h^{\text{OK}}$. Set

$$ \delta'_{T,i}(v) := \mathbf{G}_{T,i}^k (\hat{I}_T^k(v)^\mathcal{N}) + \delta_{i1} \mathbf{L}_{T1}^k (\llbracket v \rrbracket_\Gamma) - \nabla I_{T,i}^{k+1}(v_i) \in \mathbb{P}^{k+1}(T^i; \mathbb{R}). $$

(1) Let $I_{T_h}^{k+1}$ be composed of the two cell components of $\hat{I}_h^k$. We notice that

$$ \begin{aligned} \|\delta'_{T,i}(v)\|_{T^i}^2 &= (\mathbf{G}_{T,i}^k (\hat{I}_T^k(v)^\mathcal{N}) + \delta_{i1} \mathbf{L}_{T1}^k (\llbracket v \rrbracket_\Gamma) - \nabla I_{T,i}^{k+1}(v_i), \delta'_{T,i}(v))_{T^i} \\ &= (\mathbf{G}_{T,i}^k (\hat{I}_T^k(v)^\mathcal{N}) + \delta_{i1} \mathbf{L}_{T1}^k (\llbracket I_{T_h}^{k+1}(v) \rrbracket_\Gamma) - \nabla I_{T,i}^{k+1}(v_i), \delta'_{T,i}(v))_{T^i} \\ &\quad + \delta_{i1} (\mathbf{L}_{T1}^k (\llbracket v - I_{T_h}^{k+1}(v) \rrbracket_\Gamma), \delta'_{T,i}(v))_{T^i}. \end{aligned} $$

Owing to (27) and the definition of the operators $\hat{I}_T^k$ and $\mathbf{L}_{T1}^k$, we infer that

$$ \begin{aligned} \|\delta'_{T,i}(v)\|_{T^i}^2 ={}& (\Pi_{(\partial T)^i}^k(v_i) - I_{T^i}^{k+1}(v_i), \delta'_{T,i}(v) \cdot \mathbf{n}_T)_{(\partial T)^i} + \delta_{i1}(\llbracket v - I_T^{k+1}(v) \rrbracket_\Gamma, \delta'_{T,1}(v) \cdot \mathbf{n}_\Gamma)_{T^\Gamma} \\ &+ \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ \begin{aligned} & (\Pi_{(\partial S)^i}^k(v_i) - I_{T^i}^{k+1}(v_i)^+, \delta'_{T,i}(v)^+ \cdot \mathbf{n}_S)_{(\partial S)^i} \\ &+ \delta_{i1}(\llbracket v - I_T^{k+1}(v)^+ \rrbracket_\Gamma, \delta'_{T,1}(v)^+ \cdot \mathbf{n}_\Gamma)_{S^\Gamma} \end{aligned} \right\}. \end{aligned} $$

Invoking the definition of the $L^2$-orthogonal projections $\Pi_{(\partial T)^i}^k$ and $\Pi_{(\partial S)^i}^k$ gives

$$ \begin{aligned} \|\delta'_{T,i}(v)\|_{T^i}^2 ={}& (v_i - I_{T^i}^{k+1}(v_i), \delta'_{T,i}(v) \cdot \mathbf{n}_T)_{(\partial T)^i} + \delta_{i1}(\llbracket v - I_T^{k+1}(v) \rrbracket_\Gamma, \delta'_{T,1}(v) \cdot \mathbf{n}_\Gamma)_{T^\Gamma} \\ &+ \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ \begin{aligned} & (v_i - I_{T^i}^{k+1}(v_i)^+, \delta'_{T,i}(v)^+ \cdot \mathbf{n}_S)_{(\partial S)^i} \\ &+ \delta_{i1}(\llbracket v - I_T^{k+1}(v)^+ \rrbracket_\Gamma, \delta'_{T,1}(v)^+ \cdot \mathbf{n}_\Gamma)_{S^\Gamma} \end{aligned} \right\}. \end{aligned} $$

Invoking the Cauchy–Schwarz inequality and the discrete trace inequality (38a), we infer that

$$ h_T^{\frac{1}{2}} \|\delta'_{T,i}(v)\|_{T^i} \lesssim \|v_i - I_{T^i}^{k+1}(v_i)\|_{(\partial T)^i} + \delta_{i1} \|\llbracket v - I_T^{k+1}(v) \rrbracket_\Gamma\|_{T^\Gamma} \\ + \sum_{S \in \mathcal{N}_i^{-1}(T)} \left\{ \|v_i - I_{T^i}^{k+1}(v_i)^+\|_{(\partial S)^i} + \delta_{i1} \|\llbracket v - I_T^{k+1}(v)^+ \rrbracket_\Gamma\|_{S^\Gamma} \right\}. $$

Recalling the definition (40a) of $\epsilon_{T,i}(v)$ and the definition (40b) of $\epsilon_{T,i}^\Gamma(v)$ gives

$$ \|\delta'_{T,i}(v)\|_{T^i} \lesssim h_T^{-1} \{\epsilon_{T,i}(v_i) + \delta_{i1} \epsilon_{T,i}^\Gamma(v)\}, \quad (52) $$
and owing to the discrete trace inequality (38a), we obtain

$$ \|\delta'_{T,i}(v)\|_{T^i} + h_T^{\frac{1}{2}} \|\delta'_{T,i}(v)\|_{(\partial T)^i \cup T^\Gamma} \lesssim \|\delta'_{T,i}(v)\|_{T^i} \lesssim h_T^{-1} \{\epsilon_{T,i}(v_i) + \delta_{i1} \epsilon_{T,i}^\Gamma(v)\}. $$

Finally, since $\delta'_{T,i}(v) = \delta'_{T,i}(v) - \nabla(v_i - I_{T^i}^{k+1}(v_i))$, invoking the triangle inequality proves (51a).

(2) Assuming that $N_i^{-1}(T)$ is nonempty, for all $S \in N_i^{-1}(T)$, we have

$$ \delta_{T,S,i}(v) = \delta'_{T,i}(v)^+|_{S^i} + \nabla(I_{T^i}^{k+1}(v_i)^+|_{S^i} - v_i|_{S^i}). $$

The discrete trace inequality (38a) and (52) imply that

$$ \|\delta'_{T,i}(v)^+\|_{S^i} + h_S^{\frac{1}{2}} \|\delta'_{T,i}(v)^+\|_{(\partial S)^i \cup S^\Gamma} \lesssim \|\delta'_{T,i}(v)\|_{T^i} \lesssim h_T^{-1} \{\epsilon_{T,i}(v_i) + \delta_{i1} \epsilon_{T,i}^\Gamma(v)\}. $$

Invoking the triangle inequality and recalling (40a) proves (51b). □

## 4.4 Consistency

**Lemma 4.5 (Consistency).** Let $u^{\text{ex}}$ be the weak solution to (2). Assume that $u^{\text{ex}} \in H^s(\Omega_1 \cup \Omega_2)$, $s > \frac{3}{2}$. Set $\Phi_h(\hat{w}_h) := a_h(\hat{I}_h^k(u^{\text{ex}}), \hat{w}_h) - \ell_h(\hat{w}_h)$ for all $\hat{w}_h \in \hat{U}_{h0}$. The following holds:

$$ |\Phi_h(\hat{w}_h)| \lesssim \left\{ \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i h_T^{-2} \{\epsilon_{T,i}(u_i^{\text{ex}}) + \delta_{i1} \epsilon_{T,i}^{\Gamma}(u^{\text{ex}})\}^2 \right\}^{\frac{1}{2}} \|\hat{w}_h\|_{h0}. \quad (53) $$

*Proof.* Owing to (37) from Lemma 3.1, we have $\Phi_h(\hat{w}_h) = A_1 + A_2 + A_3$ with

$$ A_1 = \sum_{(T,i) \in \mathcal{P}_h} \kappa_i (\delta_{T,i}(u^{\text{ex}}), \nabla w_{T^i})_{T^i}, $$

$$ A_2 = \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i \left\{ (\delta_{T,i}(u^{\text{ex}}) \cdot n_T, w_{(\partial T)^i} - w_{T^i})_{(\partial T)^i} - \delta_{i1} (\delta_{T,1}(u^{\text{ex}}) \cdot n_\Gamma, \llbracket w_T \rrbracket_\Gamma)_{T^\Gamma} \right. \\ \left. + \sum_{S \in \mathcal{N}_i^{-1}(T)} \{(\delta_{T,S,i}(u^{\text{ex}}) \cdot n_S, w_{(\partial S)^i} - w_{S^i})_{(\partial S)^i} - \delta_{i1} (\delta_{T,S,1}(u^{\text{ex}}) \cdot n_\Gamma, \llbracket w_S \rrbracket_\Gamma)_{S^\Gamma}\} \right\}, $$

$$ A_3 = s_h^\circ (\hat{I}_h^k(u^{\text{ex}}), \hat{w}_h) + \sum_{T \in \mathcal{T}_h^{\text{cut}}} \kappa_1 h_T^{-1} (\llbracket I_T^{k+1}(u^{\text{ex}}) - u^{\text{ex}} \rrbracket_\Gamma, \llbracket w_T \rrbracket_\Gamma)_{T^\Gamma} + s_h^{\mathcal{N}} (\hat{I}_h^k(u^{\text{ex}}), \hat{w}_h), $$

where $\delta_{T,i}(u^{\text{ex}})$ and $\delta_{T,S,i}(u^{\text{ex}})$ are defined in (50) for all $(T, i) \in \mathcal{P}_h^{\text{OK}}$, whereas we set $\delta_{T,i}(u^{\text{ex}}) := \nabla(I_{T^i}^{k+1}(u_i^{\text{ex}}) - u_i^{\text{ex}})$ for all $(T, i) \in \mathcal{P}_h^{\text{KO}}$, and where $\llbracket I_T^{k+1}(u^{\text{ex}}) - u^{\text{ex}} \rrbracket_\Gamma$ is defined in Lemma 4.2. Owing to (49), we infer that the term $A_1$ can be rewritten as follows:

$$ A_1 = \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i \left\{ (\delta_{T,i}(u^{\text{ex}}), \nabla w_{T^i})_{T^i} + \sum_{S \in \mathcal{N}_i^{-1}(T)} (\delta_{T,S,i}(u^{\text{ex}}), \nabla w_{S^i})_{S^i} \right\}. $$

Invoking the Cauchy-Schwarz inequality together with the approximation results from Lemma 4.4 gives

$$ \begin{aligned} |A_1| & \lesssim \left\{ \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i h_T^{-2} \{\epsilon_{T,i}(u_i^{\text{ex}}) + \delta_{i1} \epsilon_{T,i}^{\Gamma}(u^{\text{ex}})\}^2 \right\}^{\frac{1}{2}} \left\{ \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \|\nabla w_{T^i}\|_{T^i}^2 \right\}^{\frac{1}{2}} \\ & \le \left\{ \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i h_T^{-2} \{\epsilon_{T,i}(u_i^{\text{ex}}) + \delta_{i1} \epsilon_{T,i}^{\Gamma}(u^{\text{ex}})\}^2 \right\}^{\frac{1}{2}} \|\hat{w}_h\|_{h0}, \end{aligned} $$
where the last bound follows from the definition (42) of the stability norm $||\cdot||_{h0}$. The same arguments prove that

$$|A_2| \lesssim \left\{ \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i h_T^{-2} \{\epsilon_{T,i}(u_i^{\text{ex}}) + \delta_{i1}\epsilon_{T,i}^{\Gamma}(u^{\text{ex}})\}^2 \right\}^{\frac{1}{2}} \|\hat{w}_h\|_{h0}.$$

Finally, let us write $A_3 = A_{31} + A_{32} + A_{33}$ with obvious notation. First, observing that $s_h^\circ(\hat{w}_h, \hat{w}_h) \le |\hat{w}_h|^2$, we obtain

$$|A_{31}| \le \left\{ \sum_{(T,i) \in \mathcal{P}_h} \kappa_i h_T^{-1} \| \Pi_{(\partial T)^i}^k (u_i^{\text{ex}} - I_{T^i}^{k+1} (u_i^{\text{ex}})) \|_{(\partial T)^i}^2 \right\}^{\frac{1}{2}} |\hat{w}_h|_S \\ \le \left\{ \sum_{(T,i) \in \mathcal{P}_h} \kappa_i h_T^{-1} \| u_i^{\text{ex}} - I_{T^i}^{k+1} (u_i^{\text{ex}}) \|_{(\partial T)^i}^2 \right\}^{\frac{1}{2}} |\hat{w}_h|_S,$$

where the second bound follows from the $L^2$-stability of $\Pi_{(\partial T)^i}^k$. Owing to (49), the summation on the right-hand side can be rewritten as

$$|A_{31}| \le \left\{ \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \kappa_i h_S^{-1} \|u_i^{\text{ex}} - I_{T^i}^{k+1} (u_i^{\text{ex}})\|^2 \right\}^{\frac{1}{2}} |\hat{w}_h|_S.$$

Therefore, we conclude that

$$|A_{31}| \le \left\{ \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i h_T^{-2} \epsilon_{T,i} (u_i^{\text{ex}})^2 \right\}^{\frac{1}{2}} |\hat{w}_h|_S.$$

Concerning $A_{32}$, the Cauchy-Schwarz inequality and the same arguments as above give

$$|A_{32}| \le \left\{ \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i h_T^{-2} \{\epsilon_{T,i}(u_i^{\text{ex}}) + \delta_{i1}\epsilon_{T,i}^{\Gamma}(u^{\text{ex}})\}^2 \right\}^{\frac{1}{2}} s_h^{\Gamma}(\hat{w}_h, \hat{w}_h)^{\frac{1}{2}}.$$

Finally, $s_h^N(\hat{I}_h^k(u^{\text{ex}}), \hat{w}_h) = 0$ owing to (49). Putting everything together proves the claim. $\square$

## 4.5 Error estimate

We are now ready to establish the main result of our error analysis.

**Theorem 4.6 (Error estimate).** Let $u^{\text{ex}}$ be the weak solution to (2). Assume that $u^{\text{ex}} \in H^s(\Omega_1 \cup \Omega_2)$ with $s \in (\frac{3}{2}, k+2]$. The following holds:

$$\left\{ \sum_{(T,i) \in \mathcal{P}_h} \kappa_i \|\nabla(u_i^{\text{ex}} - u_{T^i})\|_{T^i}^2 \right\}^{\frac{1}{2}} \lesssim \left\{ \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i h_T^{2(s-1)} |E_i^s(u_i^{\text{ex}})|_{H^s(\Delta_2(T))}^2 \right\}^{\frac{1}{2}} \lesssim h^{s-1} \sum_{i \in \{1,2\}} \kappa_i^{\frac{1}{2}} |u_i^{\text{ex}}|_{H^s(\Omega_i)}.$$
Proof. Define the discrete error as $\hat{e}_h := \hat{I}_h^k(u^{\text{ex}}) - \hat{u}_h \in \hat{\mathcal{U}}_{h0}$. Since $a_h(\hat{e}_h, \hat{e}_h) = \Phi_h(\hat{e}_h)$, invoking the stability result from Lemma 4.3 and the bound on the consistency error from Lemma 4.5 gives

$$ \|\hat{e}_h\|_{h0}^2 \lesssim \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i h_T^{-2} \{\epsilon_{T,i}(u_i^{\text{ex}}) + \delta_{i1} \epsilon_{T,i}^{\Gamma}(u^{\text{ex}})\}^2. $$

Owing to the approximation result from Lemma 4.2 and since $\kappa_1 \le \kappa_2$ and $\Delta(T) \subset \Delta_2(T)$, we infer that

$$ \|\hat{e}_h\|_{h0}^2 \lesssim \sum_{(T,i) \in \mathcal{P}_h^{\text{OK}}} \kappa_i h_T^{2(s-1)} |E_i^s(u_i^{\text{ex}})|_{H^s(\Delta_2(T))}^2. $$

Invoking the triangle inequality then establishes the first bound in (54). The second bound follows from the shape-regularity of the mesh, and the $H^s$-stability of the extension operators $E_i^s$. $\square$

# 5 Numerical results

In this section, we present numerical results to illustrate the convergence rates established in Theorem 4.6. We also compare the present method to the one from [12] stabilized by a cell-agglomeration procedure, and we briefly investigate some aspects related to the implementation (quadrature, local polynomial bases, pairing criterion).

<figure><img src="image_7.png" /><figcaption>Fig. 5: Circular (top row) and flower-like (bottom row) interfaces. Left column: Coarsest mesh. Central column: Outcome of pairing procedure for polynomial extension stabilization, cells in $\mathcal{T}_h^{\text{KO},1}$ are colored in green, and cells in $\mathcal{T}_h^{\text{KO},2}$ are colored in blue, arrows indicate the pairing operator. Right column: Outcome of pairing procedure for cell-agglomeration stabilization, agglomerated cells are colored in dark.</figcaption></figure>

In all cases, we consider the unit square domain $\Omega := (0,1)^2$ discretized with uniform Cartesian meshes of size $h := \sqrt{2} \times 0.1 \times 2^{-\ell}$ with $\ell \in \{0, \dots, 4\}$. The polynomial degree is taken such that $k \in \{0, \dots, 3\}$. The interface is defined using a level-set function $\Phi$, with $\Gamma := \{(x, y) \in \Omega \mid \Phi(x, y) = 0\}$ and $\Omega_i = \{(x, y) \in \Omega \mid (-1)^i \Phi(x, y) > 0\}$ for all $i \in \{1, 2\}$. We
always take $\kappa_1 := 1$, and modify the diffusivity contrast by taking $\kappa_2 = 10^m$, $m \in \{0, \dots, 4\}$, where $m = 0$ corresponds to no contrast and $m = 4$ to a highly contrasted setting. We consider two shapes for the interface: a circular shape and a flower-like shape, defined by the following level-set functions:

$$ \Phi_C(x, y) := (x - a)^2 + (y - b)^2 - R^2, \qquad (55a) $$

$$ \Phi_F(x, y) := (x - a)^2 + (y - b)^2 - R^2 + c \cos(n\theta), \qquad (55b) $$

with $\theta := \arctan(\frac{y-b}{x-a})$ if $x \ge a$, and $\theta := \pi + \arctan(\frac{y-b}{x-a})$ if $x < a$, $a := b := 0.5$, $R = \frac{1}{3}$, $c := 0.03$, and $n = 8$. Figure 5 illustrates on the coarsest mesh ($\ell = 0$) the two interfaces (left column), and the outcome of the pairing procedure for polynomial extension stabilization (central column) and for cell-agglomeration stabilization (right column).

## 5.1 Implementation aspects

Quadratures along the interface and in the cut sub-cells are realized by dividing the interface into $2^r$ segments, $r \in \mathbb{N}$, and creating a sub-triangulation of the two cut sub-cells. The construction of the sub-triangulation is performed as discussed in [12]; see, in particular, Figure 4.1 therein for an illustration. Unless explicitly stated otherwise, we set $r := 8$. We also take $\eta_N := 20$ for the stabilization bilinear form $s_h^N$ in (30). Concerning the flagging of ill-cut cells, we use the parameter $\vartheta := 0.3$. Our numerical experiments below show a marginal impact of the choice of this parameter on stability and accuracy. As this numerical observation is not backed up by theoretical results, we recommend to use a flagging parameter $\vartheta$ as indicated above.

Another important aspect is the choice of the polynomial bases in the sub-cells and the sub-faces. Both bases are realized by considering centered and scaled monomials. In particular, the barycenter of each sub-cell is used for the centering of the corresponding basis. The barycenter is computed using the above sub-triangulation. The scaling of the monomials is isotropic and uses half of the diameter of the original cell. In the ill-cut cells, say $S_i$ with $i \in \{1, 2\}$, the centering and scaling is defined by considering the merged cell $S_i \cup N_i(S_i)$.

## 5.2 Convergence rates

We consider first a test case with no contrast and no jumps ($g_D = g_N = 0$). The exact solution is taken to be (in Cartesian coordinates $(x, y)$)

$$ u^{\text{ex}}(x, y) := \sin(\pi x) \sin(\pi y). \qquad (56) $$

Errors are reported in Figure 6 for the circular interface (left panel) and the flower-like interface (right panel). Here and below, errors are evaluated using the energy norm considered in Theorem 4.6. In both panels, we compare the results obtained using stabilization by polynomial extension (solid lines) and by cell agglomeration (dashed lines). Both approaches yield similar results, with slightly better errors when using polynomial extension, especially on coarse meshes, except for $k = 3$ and the flower-like interface.

We next study the sensitivity of the errors with respect to the value of the pairing parameter $\vartheta$ used for flagging ill-cut cells. We consider again the exact solution (56). Errors are reported in Figure 7 for $k = 3$ and $\vartheta \in \{0, 0.1, 0.2, 0.3\}$ ($\vartheta = 0$ corresponds to no pairing, irrespective of the ill cut). The left (resp., right) panel corresponds to polynomial extension (resp., cell
<figure><img src="image_8.png" /><figcaption>**Fig. 6:** Errors as a function of the mesh-size for the exact solution (56). Comparison between stabilization by polynomial extension (solid lines) and by cell agglomeration (dashed lines). Left: Circular interface. Right: Flower-like interface.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7:** Errors as a function of the mesh size for the exact solution (56) for $k = 3$ and various values of the pairing parameter $\vartheta$ used for flagging ill-cut cells. Left: polynomial extension. Right: cell agglomeration.</figcaption></figure>

agglomeration). Solid lines correspond to polynomial bases centered at the barycenter of each sub-cell (the default choice in our implementation), whereas dashed lines correspond to centering at the barycenter of the original cell for both sides of the cut. We observe that this second choice yields poor results owing to poor conditioning of the local matrices, whereas the first choice is robust to the presence of ill cuts. Somewhat surprisingly, the robustness is such that it allows to flag cells with a very loose criterion ($\vartheta = 0.1$) or even to flag no cells at all ($\vartheta = 0$), and still obtain close-to-optimal errors.

A further comparison between polynomial extension and cell agglomeration is provided in Figure 8 which compares the sparsity profiles of the stiffness matrices obtained with both procedures using the lowest-order polynomial setting. We observe that both approaches lead to a similar sparsity profile, and that the stiffness matrix corresponding to cell agglomeration is
<figure><img src="image_10.png" /><figcaption>**Fig. 8:** Sparsity profiles of the stiffness matrices obtained using polynomial extension (left) and cell agglomeration (right). In the left panel, the red dots indicate the coupling blocks resulting from polynomial extension. In both panels, the global block-decomposition into cell and face unknowns is indicated by black dashed lines. Coarsest mesh ($l = 0$).</figcaption></figure>

slightly smaller (as expected). The sparsity pattern of the matrix corresponding to cell agglomeration directly results from the data structure of the agglomerated mesh, whereas the sparsity pattern of the matrix corresponding to polynomial extension is derived at the assembly stage by means of the pairing operator.

<figure><img src="image_11.png" /><figcaption>**Fig. 9:** Left: Errors for the exact solution (57) as a function of the diffusivity contrast on the finest mesh ($l = 4$) and for various polynomial degrees ($k \in \{0, \dots, 3\}$). Right: Errors as a function of the mesh size for various values of the sub-triangulation parameter, $k = 3$, $r \in \{2, 4, 6, 8\}$, and various diffusivity contrasts, $m \in \{0, \dots, 4\}$.</figcaption></figure>

We now turn our attention to a test case with diffusivity contrast and jumps across the interface. We first consider the exact solution for the circular interface such that (in radial
coordinates $(\rho, \theta))$

$$ u_1^{\text{ex}}(\rho) := \frac{\rho^6}{\kappa_1}, \qquad u_2^{\text{ex}}(\rho) := \frac{\rho^6}{\kappa_2} + R^6 \left( \frac{1}{\kappa_1} - \frac{1}{\kappa_2} \right). \qquad (57) $$

Notice that this solution does not exhibit jumps across the interface ($g_D = g_N = 0$). The errors in the highly contrasted case where $\kappa_2 = 10^4\kappa_1$ exhibit the same optimal convergence rates as those reported in Figure 6 (not shown for brevity). To highlight the robustness of the approach with respect to the diffusivity contrast, we report in the left panel of Figure 9 the errors as a function of $\kappa_2 = 10^m\kappa_1$, $m \in \{0, \dots, 4\}$, obtained on the finest mesh ($l = 4$) and for various polynomial degrees ($k \in \{0, \dots, 3\}$). The right panel of Figure 9 underlines the importance of ensuring sufficient geometric resolution when performing quadratures in the cut sub-cells. Therein, we report the error as a function of the mesh size for the sub-triangulation parameter $r \in \{2, 4, 6, 8\}$ and various diffusivity contrasts. We observe that in all cases, taking $r \ge 6$ is required to avoid that the geometric errors pollute the optimal decay of the discretization errors.

<figure><img src="image_12.png" /><figcaption>Fig. 10: Errors as a function of the mesh size for the exact solution (58a) (left) and the exact solution (58b) (right) for various polynomial degrees $k \in \{0, \dots, 3\}$; diffusivity contrast set to $\kappa_2 = 10^4\kappa_1$.</figcaption></figure>

We now include jumps across the interface while considering the highly contrasted setting where $\kappa_2 = 10^4\kappa_1$. For this purpose, we modify the exact solution in (57) and set

$$ u_1^{\text{ex}}(\rho) := \frac{\rho^6}{\kappa_1}, \qquad u_2^{\text{ex}}(\rho) := \frac{\rho^8 - R^8}{\kappa_2} + \frac{R^6}{\kappa_1}, \qquad (58a) $$

$$ u_1^{\text{ex}}(\rho) := \frac{\rho^6}{\kappa_1}, \qquad u_2^{\text{ex}}(\rho) := \frac{\rho^6}{\kappa_2}. \qquad (58b) $$

Notice that, for the exact solution (58a), we have $g_N = 2R^5(3 - 4R^2)$ and $g_D = 0$, whereas, for the exact solution (58b), we have $g_N = 0$ and $g_D = R^6(\frac{1}{\kappa_1} - \frac{1}{\kappa_2})$. Errors as a function of the mesh size are reported in Figure 10 for polynomial degrees $k \in \{0, \dots, 3\}$. Optimal convergence rates are observed in all cases. Actually, the errors for both exact solutions are very similar.

A more challenging setting is obtained when considering for the exact solution

$$ u_1^{\text{ex}}(x, y) := \cos(y)e^x, \qquad u_2^{\text{ex}}(x, y) := \sin(\pi x) \sin(\pi y), \qquad (59) $$

leading to variable jump data $g_D$ and $g_N$. Here, we do not consider a diffusivity contrast ($\kappa_1 = \kappa_2 = 1$). Errors as a function of the mesh size are reported in the left panel of Figure 11
<figure><img src="image_13.png" /><figcaption>**Fig. 11:** Left: Errors as a function of the mesh size for the exact solution (59) for various polynomial degrees $k \in \{0, \dots, 3\}$ (with sub-triangulation parameter set to $r = 10$). Right: Errors as a function of the mesh size for $k = 3$ and $r \in \{4, 6, 8, 10, 11\}$. No diffusivity contrast.</figcaption></figure>

for various polynomial degrees $k \in \{0, \dots, 3\}$ and sub-triangulation parameter set to $r = 10$. Optimal convergence rates are observed in all cases. To motivate the choice for $r$, we report in the right panel of Figure 11 the errors as a function of the mesh size for $k = 3$ and $r \in \{4, 6, 8, 10, 11\}$. We observe that the geometric error does not pollute the discretization error only for $r = 10$, except for $k = 3$ on the finest mesh where $r = 11$ is actually necessary to suppress the geometric error (notice in passing the very low values attained by the error).

<figure><img src="image_14.png" /><figcaption>**Fig. 12:** Euclidean condition number of the stiffness matrix on the coarsest mesh ($l = 0$) as a function of the radius of the circular interface (left) and as a function of the position parameter $\delta$ for the square interface, using polynomial extension (solid lines) and cell agglomeration (dashed lines).</figcaption></figure>

Finally, we perform a brief study of the conditioning of the stiffness matrix. In the left panel of Figure 12, we consider a circular interface with radius $R := \frac{1}{3} + \frac{i}{32}$, $i \in \{-4, \dots, 4\}$, and we report the Euclidean condition number as a function of the radius, using either polynomial extension (solid lines) or cell agglomeration (dashed lines). We consider the coarsest mesh
(ℓ = 0) and polynomial degrees $k \in \{0, \dots, 3\}$. We observe that, for each polynomial degree, the conditioning remains fairly insensitive to the stabilization method of the ill-cut cells. To exacerbate the effect of ill-cut cells, we now consider a square interface with level-set function $\Phi_s(x, y) := \max(x - a, y - b) - (0.25 + \delta)$ and position parameter $\delta := 0.5 \times 10^{-p}$, $p \in \{2, \dots, 9\}$. The right panel of Figure 12 reports the Euclidean condition number as a function of the position parameter of the square interface. We observe full robustness of the condition number with respect to the cut position, for both polynomial extension and cell agglomeration. Further numerical results (not displayed for brevity) indicate that a variant of the stabilization bilinear form $s_h^N$ where the penalty is enforced on the faces of $T$ with a weight $h_T^{-1}$ does not lead to a robust behavior, but instead to a linear growth in $\delta^{-1}$ for $k \in \{2, 3\}$. We also emphasize that the errors still behave optimally (figure omitted for brevity), in agreement with the above convergence analysis.

# 6 Proofs

This section collects the proofs of the two results stated in Section 4.1.

## 6.1 Proof of Lemma 4.1

*Proof.* (1) Proof of (38a). Let $(T, i) \in \mathcal{P}_h^{\text{OK}}$ and $\phi \in \mathbb{P}^\ell(T^i; \mathbb{R})$. We have

$$ \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \{\|\phi^+\|_S + h_S^{\frac{1}{2}} \|\phi^+\|_{(\partial S)_i \cup S^\Gamma}\} \lesssim \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \|\phi^+\|_S \lesssim \|\phi^+\|_{\Delta(T)}, $$

where the first bound follows from the discrete trace inequality from [16, Lemma 3.4] upon observing that the aggregated cell $T \cup \mathcal{N}_i^{-1}(T)$ satisfies the ball condition invoked therein, and the second bound follows from the construction of the pairing operator. We next observe that there is a ball $\mathcal{B}(\Delta(T))$ of diameter $\gamma h_T$ such that $\Delta(T) \subset \mathcal{B}(\Delta(T))$, where $\gamma$ only depends on the mesh shape-regularity parameter. We also recall that, since $(T, i) \in \mathcal{P}_h^{\text{OK}}$, there is a ball $\mathcal{B}(T, i)$ of diameter $\vartheta h_T$, so that $\mathcal{B}(T, i) \subset T^i$. We then infer that

$$ \|\phi^+\|_{\Delta(T)} \le \|\phi^+\|_{\mathcal{B}(\Delta(T))} \lesssim \|\phi\|_{\mathcal{B}(T,i)} \le \|\phi\|_{T^i}, $$

where the inverse inequality invoked in the second inequality follows by the arguments given in the proof of [12, Lemma 3.4]. This completes the proof of (38a).

(2) Proof of (38b). Let $(T, i) \in \mathcal{P}_h^{\text{OK}}$ and $\phi \in \mathbb{P}^\ell(T^i; \mathbb{R})$. Let $\langle \phi \rangle_{\mathcal{B}(T,i)}$ denote the mean-value of $\phi$ in $\mathcal{B}(T, i)$. Since $(I - \Pi_{(\partial S)_i}^k)(\phi^+) = (I - \Pi_{(\partial S)_i}^k)((\phi - \langle \phi \rangle_{\mathcal{B}(T,i)}))^+$, invoking the $L^2$-stability of $\Pi_{(\partial S)_i}^k$ followed by (38a) gives

$$ \begin{aligned} \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} h_S^{-\frac{1}{2}} \|(I - \Pi_{(\partial S)_i}^k)(\phi^+)\|_{(\partial S)_i} &\le \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} h_S^{-\frac{1}{2}} \|(\phi - \langle \phi \rangle_{\mathcal{B}(T,i)})^+\|_{(\partial S)_i} \\ &\lesssim h_T^{-1} \|\phi - \langle \phi \rangle_{\mathcal{B}(T,i)}\|_{T^i}, \end{aligned} $$

where we also used the mesh shape-regularity which implies that $h_T \lesssim h_S \lesssim h_T$ for all $S \in \Delta(T)$. The bound (38b) now follows from the discrete Poincaré inequality established in [12, Lemma 3.4]. □
## 6.2 Proof of Lemma 4.2

*Proof.* Let $v \in H^s(\Omega_1 \cup \Omega_2)$ with $s \in (\frac{3}{2}, k+2]$ and let $(T, i) \in \mathcal{P}_h^{\text{OK}}$. Set $\tilde{v}_i := E_i^s(v_i)$.

(1) Recall that the ball $\mathcal{B}(T, i)$ of diameter $\vartheta h_T$ is a subset of $T^i$ since $(T, i) \in \mathcal{P}_h^{\text{OK}}$. Let $\langle v_i \rangle_{\mathcal{B}(T, i)}$ denote the mean-value of $v_i$ in $\mathcal{B}(T, i)$. In this step, we prove that, for all $S \in \{T\} \cup \mathcal{N}_i^{-1}(T)$,

$$ \|\tilde{v}_i - \langle v_i \rangle_{\mathcal{B}(T,i)}\|_S \lesssim h_T \|\nabla \tilde{v}_i\|_{\Delta(T)}. \quad (60) $$

We have $\|\tilde{v}_i - \langle v_i \rangle_{\mathcal{B}(T,i)}\|_S \le \|\tilde{v}_i - \langle v_i \rangle_{\mathcal{B}(T,i)}\|_{\Delta(T)}$ since, by construction, $S \subset \Delta(T)$. Let $\langle \tilde{v}_i \rangle_{\Delta(T)}$ denote the mean-value of $\tilde{v}_i$ on $\Delta(T)$. Invoking [24, Lemma 5.7] gives

$$ \|\tilde{v}_i - \langle \tilde{v}_i \rangle_{\Delta(T)}\|_{\Delta(T)} \lesssim h_T \|\nabla \tilde{v}_i\|_{\Delta(T)}. $$

Moreover, letting $\partial \mathcal{B}(T, i)$ denote the boundary of the ball $\mathcal{B}(T, i)$ and since $\mathcal{B}(T, i) \subset \Delta(T)$, we have

$$ \begin{aligned} \|\langle v_i \rangle_{\mathcal{B}(T,i)} - \langle \tilde{v}_i \rangle_{\Delta(T)}\|_{\Delta(T)} &= |\Delta(T)|^{\frac{1}{2}} |\partial \mathcal{B}(T, i)|^{-\frac{1}{2}} \|\langle v_i \rangle_{\mathcal{B}(T,i)} - \langle \tilde{v}_i \rangle_{\Delta(T)}\|_{\partial \mathcal{B}(T,i)} \\ &\le |\Delta(T)|^{\frac{1}{2}} |\partial \mathcal{B}(T, i)|^{-\frac{1}{2}} (\|v_i - \langle v_i \rangle_{\mathcal{B}(T,i)}\|_{\partial \mathcal{B}(T,i)} + \|\tilde{v}_i - \langle \tilde{v}_i \rangle_{\Delta(T)}\|_{\partial \mathcal{B}(T,i)}), \end{aligned} $$

since $\tilde{v}_i|_{\partial \mathcal{B}(T,i)} = v_i|_{\partial \mathcal{B}(T,i)}$. Invoking a multiplicative trace inequality from $\partial \mathcal{B}(T, i)$ to $\mathcal{B}(T, i)$ gives

$$ \|v_i - \langle v_i \rangle_{\mathcal{B}(T,i)}\|_{\partial \mathcal{B}(T,i)} \lesssim h_T^{-\frac{1}{2}} \|v_i - \langle v_i \rangle_{\mathcal{B}(T,i)}\|_{\mathcal{B}(T,i)} + h_T^{\frac{1}{2}} \|\nabla v_i\|_{\mathcal{B}(T,i)}, $$

and the Poincaré inequality in $\mathcal{B}(T, i)$ then gives

$$ \|v_i - \langle v_i \rangle_{\mathcal{B}(T,i)}\|_{\partial \mathcal{B}(T,i)} \lesssim h_T^{\frac{1}{2}} \|\nabla v_i\|_{\mathcal{B}(T,i)}. $$

The same arguments, together with $\mathcal{B}(T, i) \subset \Delta(T)$ and the mesh shape-regularity, lead to

$$ \|\tilde{v}_i - \langle \tilde{v}_i \rangle_{\Delta(T)}\|_{\partial \mathcal{B}(T,i)} \lesssim h_T^{\frac{1}{2}} \|\nabla \tilde{v}_i\|_{\Delta(T)}. $$

Combining the above bounds and since $|\Delta(T)|^{\frac{1}{2}} |\partial \mathcal{B}(T, i)|^{-\frac{1}{2}} \lesssim h_T^{\frac{1}{2}}$ proves (60).

(2) The higher-order version of (60) is established by invoking the Morrey polynomial of $v_i$ based on mean-values of higher-order derivatives of $v_i$ (see, e.g., the proof of [24, Lemma 5.6]). Omitting the details for brevity, this gives, for all $0 \le m < s$, a polynomial $q_m(v_i) \in \mathbb{P}^{k+1}(T^i; \mathbb{R})$ such that, for all $S \in \{T\} \cup \mathcal{N}_i^{-1}(T)$,

$$ |\tilde{v}_i - q_m(v_i)|_{H^m(S)} \lesssim h_T^{s-m} |\tilde{v}_i|_{H^s(\Delta(T))}. \quad (61) $$

(3) Since $I_{T^i}^{k+1}(q_0(v_i)) = q_0(v_i)$, the triangle inequality gives

$$ \|\tilde{v}_i - I_{T^i}^{k+1}(v_i)^+\|_S \le \|\tilde{v}_i - q_0(v_i)\|_S + \|I_{T^i}^{k+1}(v_i - q_0(v_i))^+\|_S. $$

The first term on the right-hand side is estimated by using (61) with $m = 0$. For the second term, the discrete inverse inequality (38a) gives

$$ \|I_{T^i}^{k+1}(v_i - q_0(v_i))^+\|_S \lesssim \|I_{T^i}^{k+1}(v_i - q_0(v_i))\|_{T^i} \le \|I_{T^i}^{k+1}(v_i - q_0(v_i))\|_T \le \|\tilde{v}_i - q_0(v_i)\|_T, $$
where we used the $L^2$-orthogonality of $I_{T^i}^{k+1}$ in $T$. Invoking again (61) with $m=0$ and combining the above two bounds proves that

$$ \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \|\tilde{v}_i - I_{T^i}^{k+1}(v_i)^+\|_S \lesssim h_T^s |\tilde{v}_i|_{H^s(\Delta(T))}. $$

A similar reasoning with $q_1(v_i)$ in (61) also gives

$$ \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} h_S \|\nabla(\tilde{v}_i - I_{T^i}^{k+1}(v_i)^+)\|_S \lesssim h_T^s |\tilde{v}_i|_{H^s(\Delta(T))}. $$

Altogether, this proves that

$$ \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \left\{ \|v_i - I_{T^i}^{k+1}(v_i)^+\|_{S^i} + h_S \|\nabla(v_i - I_{T^i}^{k+1}(v_i)^+)\|_{S^i} \right\} \\ \le \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \left\{ \|\tilde{v}_i - I_{T^i}^{k+1}(v_i)^+\|_S + h_S \|\nabla(\tilde{v}_i - I_{T^i}^{k+1}(v_i)^+)\|_S \right\} \lesssim h_T^s |\tilde{v}_i|_{H^s(\Delta(T))}. $$

(4) Invoking a multiplicative trace inequality on all the faces composing $(\partial S)^i$ for all $S \in \{T\} \cup \mathcal{N}_i^{-1}(T)$ (observe that the aggregated cell $T \cup \mathcal{N}_i^{-1}(T)$ satisfies the required properties), and since $s > \frac{3}{2}$ by assumption, we infer from the above bounds that

$$ \sum_{S \in \{T\} \cup \mathcal{N}_i^{-1}(T)} \left\{ h_S^{\frac{1}{2}} \|v_i - I_{T^i}^{k+1}(v_i)^+\|_{(\partial S)^i} + h_S^{\frac{3}{2}} \|\nabla(v_i - I_{T^i}^{k+1}(v_i)^+)\|_{(\partial S)^i} \right\} \lesssim h_T^s |\tilde{v}_i|_{H^s(\Delta(T))}. $$

This completes the proof of (41a).

(5) To prove (41b), we bound the jump across $S^\Gamma$ by the triangle inequality, so that we need to estimate the traces on $S^\Gamma$ from both sides of $S^\Gamma$, for all $S \in \{T\} \cup \mathcal{N}_i^{-1}(T)$. On each side, we invoke a multiplicative trace inequality on $S^\Gamma$. This inequality is established with a slight adaptation of the arguments in the proof of [16, Lemma 3.3], whereby we make a specific choice for the apex of the cone, say $C(S^\Gamma)$, considered in that proof.

(5a) Assume first that $S = T$. The triangle inequality gives

$$ h_T^{\frac{1}{2}} \| [v - I_T^{k+1}(v)]_\Gamma \|_{T^\Gamma} \le h_T^{\frac{1}{2}} \| v_i - I_{T^i}^{k+1}(v_i) \|_{T^\Gamma} + h_T^{\frac{1}{2}} \| v_i - I_{T^i}^{k+1}(v_i) \|_{T^\Gamma}. $$

For the first term on the right-hand side, the apex of the cone $C(T^\Gamma)$ is taken to be the center of the ball $B(T, i)$. Since $C(T^\Gamma) \subset \text{conv}(T) \subset \Delta(T)$ by the assumption (4), we infer that

$$ h_T^{\frac{1}{2}} \|v_i - I_{T^i}^{k+1}(v_i)\|_{T^\Gamma} \lesssim \|\tilde{v}_i - I_{T^i}^{k+1}(v_i)^+\|_{\Delta(T)} + h_T \|\nabla(\tilde{v}_i - I_{T^i}^{k+1}(v_i)^+)\|_{\Delta(T)}, $$

and adapting the above arguments gives

$$ h_T^{\frac{1}{2}} \|v_i - I_{T^i}^{k+1}(v_i)\|_{T^\Gamma} \lesssim h_T^s |\tilde{v}_i|_{H^s(\Delta(T))}. $$

To bound $h_T^{\frac{1}{2}} \|v_i - I_{T^i}^{k+1}(v_i)\|_{T^\Gamma}$, we can consider the same cone to derive the multiplicative trace inequality, so that

$$ h_T^{\frac{1}{2}} \|v_i - I_{T^i}^{k+1}(v_i)\|_{T^\Gamma} \lesssim h_T^s |\tilde{v}_i|_{H^s(\Delta(T))}, $$
where $\tilde{v}_i := E_i^s(v_i)$. Altogether, this proves that

$$h_T^{\frac{1}{2}} \| [v - I_T^{k+1}(v)]_{\Gamma} \|_{T^{\Gamma}} \lesssim \sum_{i \in \{1,2\}} h_T^s |\tilde{v}_i|_{H^s(\Delta(T))}.$$

(5b) Assume now that $S \in \mathcal{N}_i^{-1}(T)$. Then, the apex of the cone, $C(S^{\Gamma})$, to establish the multiplicative trace inequality is taken to be the center of the ball $B(S, \bar{i})$ (indeed, if $(S, i) \in \mathcal{P}_h^{\text{KO}}$, then $(S, \bar{i}) \in \mathcal{P}_h^{\text{OK}}$). This yields $C(S^{\Gamma}) \subset \text{conv}(S) \subset \Delta(S) \subset \Delta_2(T)$. Invoking the same arguments as above then yields

$$h_S^{\frac{1}{2}} \| [I_T^{k+1}(v)^+ - v]_{\Gamma} \|_{S^{\Gamma}} \lesssim \sum_{i \in \{1,2\}} h_T^s |\tilde{v}_i|_{H^s(\Delta_2(T))}.$$

This completes the proof. □

## Acknowledgments

EB was partially supported by the EPSRC grants EP/T033126/1 and EP/V050400/1. For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising. RM is thankful to Matteo Ciceuttin (Politecnico Torino, Italy) and Guillaume Delay (Sorbonne University, France) for their instrumental help in code development, and to Zhaonan Dong (INRIA Paris) for stimulating discussions.

## References

[1] M. Abbas, A. Ern, and N. Pignet. Hybrid High-Order methods for finite deformations of hyperelastic materials. *Comput. Mech.*, 62(4):909–928, 2018.

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

[3] I. Babuška. The finite element method with penalty. *Math. Comp.*, 27:221–228, 1973.

[4] S. Badia, F. Verdugo, and A. F. Martín. The aggregated unfitted finite element method for elliptic problems. *Comput. Methods Appl. Mech. Engrg.*, 336:533–553, 2018.

[5] J. W. Barrett and C. M. Elliott. A finite-element method for solving elliptic equations with Neumann data on a curved boundary using unfitted meshes. *IMA J. Numer. Anal.*, 4(3):309–325, 1984.

[6] J. W. Barrett and C. M. Elliott. Finite element approximation of the Dirichlet problem using the boundary penalty method. *Numer. Math.*, 49(4):343–366, 1986.

[7] P. Bastian and C. Engwer. An unfitted finite element method using discontinuous Galerkin. *Internat. J. Numer. Methods Engrg.*, 79(12):1557–1576, 2009.

[8] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo. Basic principles of virtual element methods. *Math. Models Methods Appl. Sci.*, 23(1):199–214, 2013.
[9] L. Beirão da Veiga, Y. Liu, L. Mascotto, and A. Russo. The nonconforming virtual element method with curved edges. *J. Sci. Comput.*, 99(1):Paper No. 23, 35, 2024.

[10] A. Buffa, R. Puppi, and R. Vázquez. A minimal stabilization procedure for isogeometric methods on trimmed geometries. *SIAM J. Numer. Anal.*, 58(5):2711–2735, 2020.

[11] E. Burman. Ghost penalty. *C. R. Math. Acad. Sci. Paris*, 348(21-22):1217–1220, 2010.

[12] E. Burman, M. Cicuttin, G. Delay, and A. Ern. An unfitted hybrid high-order method with cell agglomeration for elliptic interface problems. *SIAM J. Sci. Comput.*, 43(2):A859–A882, 2021.

[13] E. Burman, S. Claus, P. Hansbo, M. G. Larson, and A. Massing. CutFEM: discretizing geometry and partial differential equations. *Internat. J. Numer. Methods Engrg.*, 104(7):472–501, 2015.

[14] E. Burman, G. Delay, and A. Ern. An unfitted hybrid high-order method for the Stokes interface problem. *IMA J. Numer. Anal.*, 41(4):2362–2387, 2021.

[15] E. Burman, O. Duran, and A. Ern. Unfitted hybrid high-order methods for the wave equation. *Comput. Methods Appl. Mech. Engrg.*, 389:Paper No. 114366, 23, 2022.

[16] E. Burman and A. Ern. An unfitted hybrid high-order method for elliptic interface problems. *SIAM J. Numer. Anal.*, 56(3):1525–1546, 2018.

[17] A. Cangiani, Z. Dong, E. H. Georgoulis, and P. Houston. hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes. SpringerBriefs in Mathematics. Springer, Cham, 2017.

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

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

[20] B. Cockburn, W. Qiu, and M. Solano. A priori error analysis for HDG methods using extensions from subdomains to achieve boundary conformity. *Math. Comp.*, 83(286):665–699, 2014.

[21] D. A. Di Pietro and J. Droniou. A Hybrid High-Order method for Leray-Lions elliptic equations on general meshes. *Math. Comp.*, 86(307):2159–2191, 2017.

[22] D. A. Di Pietro and A. Ern. *Mathematical aspects of discontinuous Galerkin methods*, volume 69 of *Mathématiques & Applications* [*Mathematics & Applications*]. Springer-Verlag, Berlin, 2012.

[23] D. A. Di Pietro, A. Ern, and S. Lemaire. An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators. *Comput. Methods Appl. Math.*, 14(4):461–472, 2014.
[24] A. Ern and J.-L. Guermond. Finite element quasi-interpolation and best approximation. *ESAIM Math. Model. Numer. Anal.*, 51(4):1367-1385, 2017.

[25] A. Hansbo and P. Hansbo. An unfitted finite element method, based on Nitsche's method, for elliptic interface problems. *Comput. Methods Appl. Mech. Engrg.*, 191(47-48):5537-5552, 2002.

[26] J. Haslinger and Y. Renard. A new fictitious domain approach inspired by the extended finite element method. *SIAM J. Numer. Anal.*, 47(2):1474-1499, 2009.

[27] A. Johansson and M. G. Larson. A high order discontinuous Galerkin Nitsche method for elliptic problems with fictitious boundary. *Numer. Math.*, 123(4):607-628, 2013.

[28] C. Lehrenfeld and M. Olshanskii. An Eulerian finite element method for PDEs in time-dependent domains. *ESAIM Math. Model. Numer. Anal.*, 53(2):585-614, 2019.

[29] C. Lehrenfeld and J. Schöberl. High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows. *Comput. Methods Appl. Mech. Engrg.*, 307:339-361, 2016.

[30] R. Massjung. An unfitted discontinuous Galerkin method applied to elliptic interface problems. *SIAM J. Numer. Anal.*, 50(6):3134-3162, 2012.

[31] J. Preuß. Higher order unfitted isoparametric space-time fem on moving domains. Master's thesis, NAM, University of Göttingen, 2018.

[32] W. Qiu, M. Solano, and P. Vega. A high order HDG method for curved-interface problems via approximations from straight triangulations. *J. Sci. Comput.*, 69(3):1384-1407, 2016.

[33] W. E. H. Sollie, O. Bokhove, and J. J. W. van der Vegt. Space-time discontinuous Galerkin finite element method for two-fluid flows. *J. Comput. Phys.*, 230(3):789-817, 2011.

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

[35] L. Yemm. A new approach to handle curved meshes in the hybrid high-order method. *Found. Comput. Math.*, 24(3):1049-1076, 2024.