# Frequency-explicit a posteriori error estimates for discontinuous Galerkin discretizations of Maxwell's equations

Théophile Chaumont-Frelet, Patrick Vega

► To cite this version:

Théophile Chaumont-Frelet, Patrick Vega. Frequency-explicit a posteriori error estimates for discontinuous Galerkin discretizations of Maxwell's equations. SIAM Journal on Numerical Analysis, 2024, 62 (1), pp.400-421. 10.1137/22M1516348. hal-03744230

HAL Id: hal-03744230

https://inria.hal.science/hal-03744230v1

Submitted on 2 Aug 2022

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
FREQUENCY-EXPLICIT A POSTERIORI ERROR ESTIMATES FOR
DISCONTINUOUS GALERKIN DISCRETIZATIONS OF MAXWELL'S
EQUATIONS

T. CHAUMONT-FRELET*,† AND P. VEGA‡

ABSTRACT. We propose a new residual-based a posteriori error estimator for discontinuous Galerkin discretizations of time-harmonic Maxwell's equations in first-order form. We establish that the estimator is reliable and efficient, and the dependency of the reliability and efficiency constants on the frequency is analyzed and discussed. The proposed estimates generalize similar results previously obtained for the Helmholtz equation and conforming finite element discretization of Maxwell's equations. In addition, for the discontinuous Galerkin scheme considered here, we also show that the proposed estimator is asymptotically constant-free for smooth solutions. We also present two-dimensional numerical examples that highlight our key theoretical findings and suggest that the proposed estimator is suited to drive h- and hp-adaptive iterative refinements.

KEY WORDS. a posteriori error estimates, hp-adaptivity, discontinuous Galerkin methods, high-frequency problems, Maxwell's equations

1. INTRODUCTION

Time-harmonic Maxwell's equations are a central model in a variety of applications involving electromagnetic fields [22, 27, 48]. Maxwell's equations cannot be analytically solved in complex settings frequently encountered in applications. Instead, numerical simulation tools based on discretization methods are commonly employed in practice, including boundary element [43], finite difference [46], finite element [15, 36], and discontinuous Galerkin [25, 28, 30, 48] methods. Irrespective of the discretization technique, there is a need to make the simulation as efficient as possible, either to treat problems faster or to tackle larger problems.

In this work, we focus on discontinuous Galerkin (DG) discretizations of Maxwell's equa-
tions in the first-order form, where we look for $E, H : \Omega \to \mathbb{C}^3$ such that

$$
(1.1) \qquad \left\{
\begin{array}{lll}
i\omega\varepsilon E - \nabla \times H = \mathbf{J} & \text{in } \Omega, \\
i\omega\mu\mathbf{H} + \nabla \times E = \mathbf{o} & \text{in } \Omega, \\
E \times \mathbf{n} = \mathbf{o} & \text{on } \partial\Omega, \\
\mu\mathbf{H} \cdot \mathbf{n} = 0 & \text{on } \partial\Omega,
\end{array}
\right.
$$

where $\varepsilon, \mu : \Omega \to \mathbb{C}^{3 \times 3}$ are given coefficients representing the electric permittivity and magnetic permeability of the materials contained in $\Omega$, $\mathbf{J} : \Omega \to \mathbb{C}^3$ is the (known) current density,

<footnote>*Inria, 2004 Route des Lucioles, 06902 Valbonne, France</footnote>

<footnote>†Laboratoire J.A. Dieudonné, Parc Valrose, 28 Avenue Valrose, 06108 Nice Cedex 02, 06000 Nice, France</footnote>

<footnote>‡Instituto de Matemáticas, Pontificia Universidad Católica de Valparaíso, Blanco Viel 596, Valparaíso, Chile</footnote>
and $\omega > 0$ is the frequency. As we elaborate below, this setting is rather general and, in particular, unbounded propagation media fit the framework of (1.1) if the coefficients $\varepsilon, \mu$ are suitably modified using perfectly matched layers [4, 5, 36].

An attractive solution to reduce computational costs is to use a posteriori error estimators coupled with adaptive strategies [7, 9, 47]. In this framework, after a discrete solution $(E_h, H_h)$ has been computed on a given mesh $\mathcal{T}_h$, an error estimator $\eta_K$ is associated with each element $K \in \mathcal{T}_h$. These estimators can then be used to decide where to locally refine the mesh, leading to an adaptive loop procedure. Besides, the cumulated estimator $\eta^2 := \sum_{K \in \mathcal{T}_h} \eta_K^2$ can be used to stop the adaptive loop and, more generally, to assess the error level reliably.

In this context, DG methods are especially interesting, as they allow for an easy implementation of hanging nodes and $p$-adaptivity [9]. Another appealing aspect of these methods that we highlight in Section 2.3 is their ability to impose prescribed jumps in the solution easily. For wave propagation problems, this is especially useful to inject incident fields [48, §4.2.2]. Besides, DG methods are sometimes more stable than conforming finite elements on coarse meshes (compare, e.g. [25, Theorem 3.5] and [6, Lemma 3.7]), which may be relevant in adaptive processes starting with coarse discretizations.

The key contribution of this work is the design of a new residual-based a posteriori error estimator for DG discretizations of (1.1). Specifically, our proposed estimator is reliable

$$ (1.2a) \qquad \|(E - E_h, H - H_h)\|_{\Omega} \le C \left( 1 + \max_{K \in \mathcal{T}_h} \frac{\omega h_K}{p_K \vartheta_K} \right) (1 + \gamma_{\text{ba}}) \eta $$

and locally efficient

$$ (1.2b) \qquad \eta_K \le C p_K^{3/2} \left( 1 + \frac{\omega h_K}{p_K \vartheta_K} \right) \|(E - E_h, H - H_h)\|_{\tilde{K}} + \text{osc}_K, $$

where $\|(\cdot, \cdot)\|$ is a suited “energy norm”, $h_K$ and $p_K$ are the size and the polynomial degree associated with the element $K \in \mathcal{T}_h$, and $\vartheta_K$ is a measure of the wavespeed around $K$. The real number $\gamma_{\text{ba}} > 0$ in (1.2a) is the so-called “approximation factor” [12, 15, 24]. It generally grows with the frequency but tends to zero as $\max_{K \in \mathcal{T}_h} h_K/p_K \to 0$, see Section 2.12 below as well as [15, 16] for more details. The term $\text{osc}_K$ in (1.2b) is a “data oscillation term” customary in efficiency estimates. We refer the reader to Theorems 4.7 and 4.11 below, where the estimates in (1.2) are established. Interestingly, these results generalize similar findings for the Helmholtz equation [24, 42] and conforming Nédélec approximations of Maxwell’s equations [15].

On the analysis side, although the estimates we obtain for DG discretization in (1.2) are rather similar to the ones established in [15] for conforming Nédélec elements, the arguments employed are fairly different. In particular, following [15, 24], our analysis relies on duality arguments that turn out to be substantially complicated in the first-order setting we consider here in (1.1), as can be seen in Lemma 4.5 below. DG schemes also have some unique properties as compared to Nédélec discretizations. For instance, we also show that for smooth solutions (or locally refined meshes), our estimator is asymptotically constant-free, meaning that

$$ (1.3) \qquad \|(E - E_h, H - H_h)\|_{\Omega} \le \left( 1 + \theta \left( \max_{K \in \mathcal{T}_h} h_K / p_K \right) \right) \eta $$

with $\lim_{t \to 0} \theta(t) = 0$. This is detailed in Remark 4.8. Another interesting fact is that for divergence-free right-hand sides $(\nabla \cdot J = 0)$, the estimator is “oscillation free”.
We also present a series of detailed two-dimensional numerical experiments. On the one hand, we showcase a few academic benchmarks with known analytical solutions that we use to highlight the main features of our results. In particular, the interplay between the frequency, the mesh size and polynomial degree, and the constants appearing in the estimates (1.2a), (1.2b), and (1.3) is largely illustrated and thoroughly discussed. On the other hand, we also consider more realistic benchmarks where the analytic solution is unavailable. In these cases, we couple the estimator with Dörfler's marking [23] and newest vertex bisection [8] to drive adaptive mesh refinements. We also consider a simple hp-adaptive procedure based on [35]. For these examples, we observe optimal convergence rates in all cases, which indicates that the estimator is suited to drive this kind of adaptive refinements.

The remainder of our work is organized as follows. In Section 2, we precise the setting and main assumptions, and recall useful standard tools. Section 3 presents the range of DG schemes for which our analysis applies and gives the construction of the estimator. Section 3 contains the proof of our main results. Finally, we report on numerical examples in Section 5 before providing some concluding remarks in Section 6.

## 2. PRELIMINARIES

We start by providing key notations and preliminary results.

### 2.1. Domain and coefficients.
We consider Maxwell's equations (1.1) in a Lipschitz polyhedral (not necessarily simply connected) domain $\Omega \subset \mathbb{R}^3$. We denote by $l_{\Omega} := \sup_{x,y \in \Omega} |x - y|$ the diameter of $\Omega$.

The coefficients $\varepsilon, \mu : \Omega \to S(\mathbb{C}^3)$ are two symmetric (but not necessarily self-adjoint) tensor-valued functions describing the electromagnetic properties of the materials contained inside $\Omega$. For the sake of simplicity, we assume that $\Omega$ can be partitioned into a set $\mathcal{P}$ of non-overlapping polyhedral subdomains $P$ such that $\varepsilon|_P$ and $\mu|_P$ are constant for all $P \in \mathcal{P}$. The short-hand notations $\zeta := \varepsilon^{-1}$ and $\chi := \mu^{-1}$ will also be useful.

If $\phi$ is any of the tensor fields mentioned above, we introduce the notations

$$ \phi_{\min}(\mathbf{x}) := \min_{\substack{\mathbf{u} \in \mathbb{C}^3 \\ |\mathbf{u}|=1}} \mathrm{Re} \, \phi(\mathbf{x})\mathbf{u} \cdot \bar{\mathbf{u}}, \quad \phi_{\max}(\mathbf{x}) := \max_{\substack{\mathbf{u} \in \mathbb{C}^3 \\ |\mathbf{u}|=1}} \max_{\substack{\mathbf{v} \in \mathbb{C}^3 \\ |\mathbf{v}|=1}} \mathrm{Re} \, \phi(\mathbf{x})\mathbf{u} \cdot \bar{\mathbf{v}}, $$

as well as

$$ \phi_{D,\min} := \underset{\mathbf{x} \in D}{\operatorname{ess\,inf}} \phi_{\min}(\mathbf{x}), \quad \phi_{D,\max} := \underset{\mathbf{x} \in D}{\operatorname{ess\,sup}} \phi_{\max}(\mathbf{x}) $$

for any open set $D \subset \Omega$, and we assume that $\phi_{\Omega,\min} > 0$. Finally, for $D \subset \Omega$, the notation

$$ c_{D,\min} := (\varepsilon_{D,\max}\mu_{D,\max})^{-1/2}, \quad c_{D,\max} := (\varepsilon_{D,\min}\mu_{D,\min})^{-1/2} $$

are employed for the smallest and highest wavespeeds in $D$.

As described in [15, Remark 2.1], these assumptions cover most interesting application scenarios, including dissipative materials and perfectly matched layers [4, 5, 36].

### 2.2. Functional spaces.
In the following, if $D \subset \Omega$, $\mathcal{L}^2(D)$ denotes the space of square-integrable complex-valued function over $D$, see e.g. [1], and $\mathcal{L}^2(D) := (\mathcal{L}^2(D))^3$. We equip $\mathcal{L}^2(D)$ with the following (equivalent) norms

$$ \|\mathbf{w}\|_D^2 := \int_D |\mathbf{w}|^2, \quad \|\mathbf{w}\|_{\phi,D}^2 := \mathrm{Re} \int_D \phi \mathbf{w} \cdot \bar{\mathbf{w}}, \quad \mathbf{w} \in \mathcal{L}^2(D) $$
where $\phi \in \{\varepsilon, \mu, \zeta, \chi\}$. We denote by $(\cdot, \cdot)_D$ the inner-product of $\mathcal{L}^2(D)$, and we drop the subscript when $D = \Omega$. If $F \subset \bar{\Omega}$ is a two-dimensional measurable planar subset, $\|\cdot\|_F$ and $\langle \cdot, \cdot \rangle_F$ respectively denote the natural norm and inner-product of both $\mathcal{L}^2(F)$ and $\mathcal{L}^2(F) := (\mathcal{L}(F))^3$.

Classically [1], we employ the notation $\mathcal{H}^1(D)$ for the usual Sobolev space of functions $w \in \mathcal{L}^2(D)$ such that $\nabla w \in \mathcal{L}^2(D)$. We also set $\mathcal{H}^1(D) := (\mathcal{H}^1(D))^3$ and introduce the semi-norms

$$ \|\nabla w\|_D^2 := \sum_{j,k=1}^3 \int_D \left| \frac{\partial w_j}{\partial x_k} \right|^2, \quad \|\nabla w\|_{\phi,D}^2 := \sum_{j,k=1}^3 \int_D \phi_{\max} \left| \frac{\partial w_j}{\partial x_k} \right|^2, $$

for $w \in \mathcal{H}^1(\Omega)$ and $\phi \in \{\varepsilon, \mu, \zeta, \chi\}$.

We shall also need Sobolev spaces of vector-valued functions with well-defined divergence and rotation [26]. Specifically, we denote by $\mathcal{H}(\text{curl}, D)$ the space of functions $w \in \mathcal{L}^2(D)$ with $\nabla \times w \in \mathcal{L}^2(D)$, that we equip with the norms

$$ \|v\|_{\text{curl}, \omega, \phi, \psi, D}^2 := \omega^2 \|v\|_{\phi, D}^2 + \|\nabla \times v\|_{\psi, D}^2, \quad v \in \mathcal{H}(\text{curl}, D), $$

for $\phi, \psi \in \{\varepsilon, \mu, \zeta, \chi\}$. We also introduce the following “energy” norm

$$ \|(e, h)\|_{\text{curl}, \omega, D}^2 := \|e\|_{\text{curl}, \omega, \varepsilon, \chi, D}^2 + \|h\|_{\text{curl}, \omega, \mu, \zeta, D}^2, \quad (e, h) \in \mathcal{H}_0(\text{curl}, D) \times \mathcal{H}(\text{curl}, D), $$

for electromagnetic fields on the product space.

In addition, if $\xi: D \to \mathbb{C}^3$ is a measurable tensor-valued function, we will use the notation $\mathcal{H}(\text{div}, \xi, D)$ for the set of functions $w \in \mathcal{L}^2(D)$ with $\nabla \cdot (\xi w) \in L^2(D)$, and we will write $\mathcal{H}(\text{div}^0, \xi, D)$ for the set of fields $w \in \mathcal{H}(\text{div}, \xi, D)$ such that $\nabla \cdot (\xi w) = 0$ in $D$. When $\xi = \mathbf{I}$, the identity tensor, we simply write $\mathcal{H}(\text{div}, D)$ and $\mathcal{H}(\text{div}^0, D)$.

For any of the aforementioned spaces $\mathcal{V}$, the notation $\mathcal{V}_0$ denotes the closure of smooth, compactly supported functions into $\mathcal{L}^2(D)$ (or $\mathcal{L}^2(D)$) with respect to the norm of $\mathcal{V}$. These spaces also correspond to the kernel of the naturally associated trace operators [1, 26].

Finally, if $\mathcal{D}$ is a collection of disjoint sets $D \subset \Omega$ and $\mathcal{V}(D)$ is any of the aforementioned spaces, $\mathcal{V}(\mathcal{D})$ stands for the “broken” space of functions in $v \in \mathcal{L}^2(\Omega)$ (or $\mathcal{L}^2(\Omega)$) such that $v|_D \in \mathcal{V}(D)$ for all $D \in \mathcal{D}$. We employ the same notation for the inner-products, norms, and semi-norms of $\mathcal{V}(\mathcal{D})$ and $\mathcal{V}(D)$, with the subscript $\mathcal{D}$ instead of $D$.

2.3. **TF-SF formulation.** We consider a TF-SF interface $\Gamma_i$ that is either empty (there are no incident fields) or the boundary of a Lipschitz polyhedral subdomain $\Omega_{tf} \subset \Omega$.

We then consider incident fields $E^{\text{inc}}, H^{\text{inc}} \in \mathcal{H}(\text{curl}, \Omega_{tf}) \cap \mathcal{H}(\text{div}, \Omega_{tf}) \cap \mathcal{L}^2(\Gamma_i)$ such that

$$ (2.1) \qquad \begin{cases} i\omega\varepsilon_0 E^{\text{inc}} - \nabla \times H^{\text{inc}} = \mathbf{o} \\ i\omega\mu_0 H^{\text{inc}} + \nabla \times E^{\text{inc}} = \mathbf{o} \end{cases} $$

in $\Omega_{tf}$, where $\varepsilon_0, \mu_0 > 0$ are arbitrary values that typically correspond to the vacuum electric permittivity and magnetic permeability in applications. In addition, for the sake of simplicity, we assume that $\varepsilon = \varepsilon_0 \mathbf{I}$ and $\mu = \mu_0 \mathbf{I}$ in a neighborhood of $\Gamma_i$.

The TF-SF interface $\Gamma_i$ is then employed to inject the incident fields $E^{\text{inc}}$ and $H^{\text{inc}}$ into the computational domain via jump conditions, see e.g. [48, §4.2.2].
2.4. **Variational formulations.** In the remainder of this work, we assume that $J \in \mathcal{H}(\text{div}, \Omega)$.
Then, we may recast (1.1) into a weak formulation, which consists in finding a pair $(E, H) \in
\mathcal{L}^2(\Omega) \times \mathcal{L}^2(\Omega)$ satisfying

$$
(2.2) \quad b((E, H), (v, w)) = i\omega(J, v) + i\omega\ell(v, w) \quad \forall (v, w) \in \mathcal{H}_0(\text{curl}, \Omega) \times \mathcal{H}(\text{curl}, \Omega)
$$

where

$$
b((e, h), (v, w)) := -\omega^2 (\mu h, w) + i\omega (e, \nabla \times w) - i\omega (h, \nabla \times v) - \omega^2 (\epsilon e, v)
$$

and

$$
\ell(v, w) := \langle \mathbf{H}^{\text{inc}}, v \times n_i \rangle_{\Gamma_i} - \langle \mathbf{E}^{\text{inc}}, w \times n_i \rangle_{\Gamma_i},
$$

for all $e, v \in \mathcal{L}^2(\Omega)$ and $(v, w) \in \mathcal{H}_0(\text{curl}, \Omega) \times \mathcal{H}(\text{curl}, \Omega)$. Notice that the duality pairing in the definition of $\ell$ is well-defined since $\mathbf{E}^{\text{inc}}, \mathbf{H}^{\text{inc}} \in \mathcal{H}(\text{curl}, \Omega_{\text{tf}})$ if we understand $\langle \cdot, \cdot \rangle_{\Gamma_i}$ as the duality pairing introduced in [10].

The variational formulations associated with second-order forms of Maxwell's equations
will also be useful. As a result, we introduce

$$
b_{\mathrm{E}}(e, v) := -\omega^2(\varepsilon e, v) + (\chi \nabla \times e, \nabla \times v) \quad \forall e, v \in \mathcal{H}_0(\mathrm{curl}, \Omega)
$$

and

$$
b_{\mathrm{H}}(h, w) := -\omega^2(\mu h, w) + (\zeta \nabla \times h, \nabla \times w) \quad \forall h, w \in \mathcal{H}(\mathrm{curl}, \Omega),
$$

and observe that the Gårding inequalities

$$
(2.3\text{a}) \quad \|e\|_{\text{curl},\omega,\varepsilon,\chi,\Omega}^2 = \operatorname{Re} b_{\text{E}}(e, e) + 2\omega^2 \|e\|_{\varepsilon,\Omega}^2 \quad \forall e \in \mathcal{H}_0(\text{curl}, \Omega)
$$

and

$$
(2.3\text{b}) \quad \|h\|_{\text{curl},\omega,\mu,\zeta,\Omega}^2 = \operatorname{Re} b_{\text{H}}(h,h) + 2\omega^2 \|h\|_{\mu,\Omega}^2 \quad \forall h \in \mathcal{H}(\text{curl},\Omega)
$$

are satisfied.

2.5. **Computational mesh.** The computational mesh $\mathcal{T}_h$ is a partition of $\Omega$ into non-overlapping (closed) simplicial elements $K$. We denote by $\mathcal{F}_h^e$ the set of exterior faces lying on the boundary $\partial\Omega$, by $\mathcal{F}_h^i$ the remaining (interior) faces, and set $\mathcal{F}_h := \mathcal{F}_h^e \cup \mathcal{F}_h^i$. We associate with each face $F \in \mathcal{F}_h^i$ a unit normal vector $n_F$ whose orientation is arbitrary but fixed. For $F \in \mathcal{F}_h^e$, $n_F := n$ is the outward unit vector normal to $\partial\Omega$. For $K \in \mathcal{T}_h$, $\mathcal{F}_K \subset \mathcal{F}_h$ denotes the faces of $K$. The notations

$$
h_K := \sup_{\mathbf{x},\mathbf{y} \in K} |\mathbf{x} - \mathbf{y}|, \quad \rho_K := \sup \{r > 0 \mid \exists \mathbf{x} \in K : B(\mathbf{x}, r) \subset K\},
$$

stand for the diameter of $K$ and the radius of the largest ball contained in $\bar{K}$, and $\beta_K := h_K/\rho_K$ is its shape-regularity parameter. The (global) mesh size, and shape-regularity pa- rameters are respectively defined as $h := \max_{K \in \mathcal{T}_h} h_K$ and $\beta := \max_{K \in \mathcal{T}_h} \beta_K$.

We assume that $\mathcal{T}_h$ is conforming in the sense of [18], that is, the intersection $\overline{K_+} \cap \overline{K_-}$ of two distinct elements $K_\pm \in \mathcal{T}_h$ is either empty, or a single vertex, edge, or face of both $K_-$ and $K_+$. We further require that the mesh $\mathcal{T}_h$ is conforming with the physical partition $\mathcal{P}$. Namely, we assume that for each $K \in \mathcal{T}_h$, there exists $P \in \mathcal{P}$ such that $K \subset P$, which ensures that coefficients $\varepsilon$ and $\mu$ are constant in each element.
**Remark 2.1** (Hanging nodes). *Discontinuous Galerkin methods allow hanging nodes that violate the above assumption and can be especially beneficial in mesh adaptivity techniques [9]. We believe that the present analysis could extend to meshes featuring hanging nodes, but at the price of increased technicalities in the definition of the quasi-interpolation operators described in Section 2.7 below.*

We follow the standard convention for jumps and averages of functions. Namely, if $v \in \mathcal{H}^1(\mathcal{T}_h)$, we define

$$[v]|_F := v_+(n_+ \cdot n_F) + v_-(n_- \cdot n_F), \quad \{v\}|_F := \frac{1}{2}(v_+ + v_-),$$

for all interior faces $F := \partial K_- \cap \partial K_+ \in \mathcal{F}_h^i$ with $v_\pm$ the trace of $v$ on $F$ from the interior of $K_\pm$ and $n_\pm$ the unit normal pointing outward of $K_\pm$. For exterior faces $F \in \mathcal{F}_h^e$, the definitions $[v]|_F := \{v\}|_F := v|_F$ are convenient.

If $K \in \mathcal{T}_h$ and $F \in \mathcal{F}_h$, the associated mesh patches are defined by

$$\mathcal{T}_{K,h} := \{K' \in \mathcal{T}_h \mid \overline{K} \cap \overline{K'} \neq \emptyset\}, \quad \mathcal{T}_{F,h} := \{K' \in \mathcal{T}_h \mid F \subset \partial K'\},$$

and we respectively use the notations $\tilde{K}$ and $\tilde{F}$ for the open domain covered by the elements of $\mathcal{T}_{K,h}$ and $\mathcal{T}_{F,h}$.

For collections of elements and faces $\mathcal{T} \subset \mathcal{T}_h$ and $\mathcal{F} \subset \mathcal{F}_h$, the following broken inner-product will be useful:

$$(\cdot, \cdot)_{\mathcal{T}} := \sum_{K \in \mathcal{T}} (\cdot, \cdot)_K, \quad \langle \cdot, \cdot \rangle_{\mathcal{T}} := \sum_{K \in \mathcal{T}} \langle \cdot, \cdot \rangle_{\partial K}, \quad \langle \cdot, \cdot \rangle_{\mathcal{F}} := \sum_{F \in \mathcal{F}} \langle \cdot, \cdot \rangle_F.$$

**2.6. Polynomial spaces.** In the following, for all $K \in \mathcal{T}_h$ and $q \ge 0$, $\mathcal{P}_q(K)$ stands for the space of (complex-valued) polynomials defined over $K$ and $\mathcal{P}_q(K) := (\mathcal{P}_q(K))^3$. If $\mathcal{T} \subset \mathcal{T}_h$ and $\mathbf{q} := \{q_K\}_{K \in \mathcal{T}}$, then $\mathcal{P}_q(\mathcal{T})$ and $\mathcal{P}_q(\mathcal{T})$ respectively stand for the space of functions that are piecewise in $\mathcal{P}_{q_K}(K)$ and $\mathcal{P}_{q_K}(K)$ for all $K \in \mathcal{T}$.

In the remaining, we associated with each element $K \in \mathcal{T}_h$ a polynomial degree $p_K \ge 1$, and we set $\mathbf{p} := \{p_K\}_{K \in \mathcal{T}_h}$. For the sake of simplicity, we will assume that there exists a constant $\mathcal{C}_{\text{pol}} > 0$ such that

$$(2.4) \qquad \mathcal{C}_{\text{pol}}^{-1} p_K \le p_{K'} \le \mathcal{C}_{\text{pol}} p_K$$

for all neighboring elements $K$ and $K'$ in $\mathcal{T}_h$. We also set $p_F := \max(p_K, p_{K'})$ for $F := \partial K_- \cap \partial K_+ \in \mathcal{F}_h^i$ and $p_F := p_K$ for exterior faces $F \in \mathcal{F}_h^e \cap \mathcal{F}_K$. Then,

$$V_h := \mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega), \quad \tilde{V}_h := \mathcal{P}_p(\mathcal{T}_h) \cap H^1(\Omega)$$

are the usual Lagrange finite element spaces (with and without essential boundary conditions) and

$$\mathbf{W}_h := \mathcal{P}_p(\mathcal{T}_h) \cap \mathcal{H}_0(\text{curl}, \Omega), \quad \tilde{\mathbf{W}}_h := \mathcal{P}_p(\mathcal{T}_h) \cap \mathcal{H}(\text{curl}, \Omega)$$

are the usual second-family of Nédélec spaces [37].

**2.7. Quasi-interpolation.** There exists two operators $\mathcal{Q}_h : H_0^1(\Omega) \to V_h$ and $\mathcal{R}_h : \mathcal{H}_0(\text{curl}, \Omega) \to \mathbf{W}_h$ and a constant $\mathcal{C}_i$ that only depends on $\beta$ such that

$$(2.5) \qquad \frac{p_K}{h_K} \|w - \mathcal{Q}_h w\|_K + \sqrt{\frac{p_K}{h_K}} \|w - \mathcal{Q}_h w\|_{\partial K} \le \mathcal{C}_i \|\nabla w\|_{\tilde{K}}$$
for all $w \in H_0^1(\Omega)$ and

$$ (2.6) \qquad \frac{p_K}{h_K} \|w - \mathcal{R}_h w\|_K + \sqrt{\frac{p_K}{h_K}} \|(w - \mathcal{R}_h w) \times n\|_{\partial K} \le \mathcal{C}_i \|\nabla w\|_{\tilde{K}} $$

for all $w \in \mathcal{H}^1(\mathcal{T}_h) \cap \mathcal{H}_0(\text{curl}, \Omega)$. We will also use quasi-interpolation operators that operate on spaces without essential boundary conditions, namely, $\tilde{Q}_h : H^1(\Omega) \to \tilde{V}_h$ and $\tilde{\mathcal{R}}_h : \mathcal{H}(\text{curl}, \Omega) \to \tilde{\mathcal{W}}_h$. These operators also satisfy (2.5) and (2.6) for all $w \in H^1(\Omega)$ and $w \in \mathcal{H}^1(\mathcal{T}_h) \cap \mathcal{H}(\text{curl}, \Omega)$. We refer the reader to, e.g, [29, 31] for the construction of $Q_h$ and $\tilde{Q}_h$. $\mathcal{R}_h$ and $\tilde{\mathcal{R}}_h$ are then respectively defined by using $Q_h$ and $\tilde{Q}_h$ componentwise.

**2.8. Bubble functions and inverse inequalities.** Bubble functions constitute a standard tool that we will use to prove efficiency estimates [24, 35]. For all elements $K \in \mathcal{T}_h$ and faces $F \in \mathcal{F}_h$, there exists “bubble” functions $b_K \in C^0(\bar{K})$ and $b_F \in C^0(\bar{F})$ such that the following holds.<fnref n="1" /> We have [35, Lemma 2.5]

$$ (2.7a) \qquad \|w\|_K \le \mathcal{C}_b p_K \|b_K^{1/2} w\|_K, \qquad \|v\|_F \le \mathcal{C}_b p_F \|b_F^{1/2} v\|_F $$

for all $w \in \mathcal{P}_p(K)$ and $v \in \mathcal{P}_p(F)$. Here, $\mathcal{C}_b > 0$ is a constant depending on the shape regularity parameter $\beta$. Besides, [35, Lemma 2.5] shows that

$$ (2.7b) \qquad \|\nabla(b_K w)\|_K \le \mathcal{C}_b \frac{p_K}{h_K} \|b_K^{1/2} w\|_K \quad \forall w \in \mathcal{P}_p(K). $$

Finally, [35, Lemma 2.6] guarantees the existence of an extension operator $\mathcal{E} : \mathcal{P}_p(F) \to H_0^1(\bar{F})$ such that $\mathcal{E}(v)|_F = b_F v$ and

$$ (2.7c) \qquad p_K h_K^{-1/2} \|\mathcal{E}(v)\|_{\mathcal{T}_{F,h}} + p_K^{-1} h_K^{1/2} \|\nabla \mathcal{E}(v)\|_{\mathcal{T}_{F,h}} \le \mathcal{C}_b \|b_F^{1/2} v\|_F \quad \forall v \in \mathcal{P}_p(F). $$

Identical results hold for vector-valued functions, applying the above estimates componentwise.

**2.9. Data oscillation.** Our estimates include a “data oscillation” term which is standard in a posteriori error estimation. Perhaps surprisingly, this oscillation term only involves the divergence of the right-hand side $J$, and not its actual values. In particular, there are no oscillation terms in the common case where $\nabla \cdot J = 0$ and $E^{\text{inc}} = H^{\text{inc}} = o$. We thus set

$$ \varrho_h := \arg \min_{q_h \in \mathcal{P}_p(\mathcal{T}_h)} \|\nabla \cdot J - q_h\|_{\Omega}, $$

$$ E_h^{\text{inc}} := \arg \min_{v_h \in \mathcal{P}_p(\mathcal{F}_h^i)} \|E^{\text{inc}} - v_h\|_{\Gamma_i}, $$

$$ H_h^{\text{inc}} := \arg \min_{w_h \in \mathcal{P}_p(\mathcal{F}_h^i)} \|H^{\text{inc}} - w_h\|_{\Gamma_i}, $$

and

$$ \mathrm{osc}_K := \mathrm{osc}_K(J) + \mathrm{osc}_K(E^{\mathrm{inc}}) + \mathrm{osc}_K(H^{\mathrm{inc}}), \qquad \mathrm{osc}_{\mathcal{T}}^2 := \sum_{K \in \mathcal{T}} \mathrm{osc}_K $$

<footnote n="1">¹The results in [35] are rigorously stated for the two-dimensional case. However, as observed in [24, Theorem 4.12], these results naturally extend to the three-dimensional case.</footnote>
for all $K \in \mathcal{T}_h$ and $\mathcal{T} \subset \mathcal{T}_h$, where

$$ \mathrm{osc}_K(\mathbf{J}) := p_K^{3/2} \frac{1}{\sqrt{\varepsilon_{\tilde{K},\min}}} \frac{h_K}{p_K} \|\nabla \cdot \mathbf{J} - \varrho_h\|_K, $$

$$ \mathrm{osc}_K(E^{\mathrm{inc}}) := p_K^2 \sqrt{\varepsilon_{\tilde{K},\max}} \omega \sqrt{\frac{h_K}{p_K}} \|E^{\mathrm{inc}} - E_h^{\mathrm{inc}}\|_{\partial K \cap \Gamma_i}, $$

$$ \mathrm{osc}_K(H^{\mathrm{inc}}) := p_K^2 \sqrt{\mu_{\tilde{K},\max}} \omega \sqrt{\frac{h_K}{p_K}} \|H^{\mathrm{inc}} - H_h^{\mathrm{inc}}\|_{\partial K \cap \Gamma_i}. $$

**2.10. Regular decomposition.** Regular decompositions play an essential role in the derivation of reliability estimates for a posteriori error estimators in the context of $\mathcal{H}_0(\mathrm{curl}, \Omega)$ problems [3, 17, 20, 39, 44], and we refer the reader to [19, 26, 29] for a thorough discussion of this topic.

The results we need follow from [29, Theorem 2.1] and read as follow: for all $\theta \in \mathcal{H}_0(\mathrm{curl}, \Omega)$, there exist $\phi \in \mathcal{H}_0^1(\Omega)$ and $r \in \mathcal{H}_0^1(\Omega)$ such that $\theta = \phi + \nabla r$ with

$$ (2.8a) \qquad \|\nabla \phi\|_{\chi, \Omega} \le \mathcal{C}_r \|\theta\|_{\mathrm{curl}, \omega, \varepsilon, \chi, \Omega}, \qquad \|\nabla r\|_{\varepsilon, \Omega} \le \mathcal{C}_r \|\theta\|_{\varepsilon, \Omega}, $$

and similarly, for all $\vartheta \in \mathcal{H}(\mathrm{curl}, \Omega)$ there exist $\psi \in \mathcal{H}^1(\Omega)$ and $s \in \mathcal{H}^1(\Omega)$ such that $\vartheta = \psi + \nabla s$ with

$$ (2.8b) \qquad \|\nabla \psi\|_{\zeta, \Omega} \le \mathcal{C}_r \|\vartheta\|_{\mathrm{curl}, \omega, \mu, \zeta, \Omega}, \qquad \|\nabla s\|_{\mu, \Omega} \le \mathcal{C}_r \|\vartheta\|_{\mu, \Omega}. $$

In (2.8), $\mathcal{C}_r$ is a constant possibly depending on the geometry of $\Omega$, and, since the result in [29] is enunciated in non-weighted norms, the material contrasts $\varepsilon_{\Omega,\max}/\varepsilon_{\Omega,\min}$ and $\mu_{\Omega,\max}/\mu_{\Omega,\min}$. In addition, it may also depend on $(\omega\ell_{\Omega}/c_{\Omega,\max})^{-1}$ if the domain is not simply-connected, but the constant can only blow up in the low-frequency regime and remains well-behaved in the high-frequency regime on which we focus here (see the discussion in [15]).

**2.11. Well-posedness.** We will work under the assumption that the (adjoint) problem under consideration is well-posed for the fixed frequency $\omega$.

**Assumption 2.2 (Well-posedness).** For all $j \in \mathcal{L}^2(\Omega)$, there exists a unique $e^*(j) \in \mathcal{H}_0(\mathrm{curl}, \Omega)$ such that

$$ (2.9) \qquad b_E(w, e^*(j)) = \omega(w, \varepsilon j) \quad \forall w \in \mathcal{H}_0(\mathrm{curl}, \Omega). $$

Classically, we can infer from Assumption 2.2 that the first-order variational formulation is well-posed (see, e.g. [11, Lemma 3.1], or the proof of Lemma 4.6 below). Similarly, for all $l \in \mathcal{L}^2(\Omega)$, there exists a unique element $h^*(l) \in \mathcal{H}(\mathrm{curl}, \Omega)$ such that

$$ (2.10) \qquad b_H(w, h^*(l)) = \omega(w, \mu l) \quad \forall w \in \mathcal{H}(\mathrm{curl}, \Omega). $$

Thus, we can introduce the notations

$$ (2.11) \qquad \gamma_{\mathrm{st},\mathrm{E}} := \sup_{\substack{j \in \mathcal{H}(\mathrm{div}^0, \varepsilon, \Omega) \\ \|j\|_{\varepsilon, \Omega}=1}} \|e^*(j)\|_{\mathrm{curl}, \omega, \varepsilon, \chi, \Omega}, \qquad \gamma_{\mathrm{st},\mathrm{H}} := \sup_{\substack{l \in \mathcal{H}_0(\mathrm{div}^0, \mu, \Omega) \\ \|l\|_{\mu, \Omega}=1}} \|h^*(l)\|_{\mathrm{curl}, \omega, \mu, \zeta, \Omega}, $$

and set $\gamma_{\mathrm{st}} := \gamma_{\mathrm{st},\mathrm{E}} + \gamma_{\mathrm{st},\mathrm{H}}$.
2.12. **Approximation factors.** As is now standard for high-frequency wave propagation problems [12, 15, 24, 42], our analysis will rely on some “approximation factor” that is employed in the context of a duality argument. Because we study the problem in first-order form, the definition is slightly different than the one proposed earlier for problems in second-order form.

We start by introducing two approximation factors that respectively describe the ability of the discrete spaces $W_h$ and $\tilde{W}_h$ to approximate solutions to (2.9) and (2.10) and are defined by

$$ (2.12a) \qquad \gamma_{ba,E} := \sup_{\substack{j \in \mathcal{H}(\operatorname{div}^0, \varepsilon, \Omega) \\ \|j\|_{\varepsilon, \Omega}=1}} \inf_{e_h \in W_h} \|e^*(j) - e_h\|_{\mathrm{curl}, \omega, \varepsilon, \chi, \Omega} $$

and

$$ (2.12b) \qquad \gamma_{ba,H} := \sup_{\substack{l \in \mathcal{H}_0(\operatorname{div}^0, \mu, \Omega) \\ \|l\|_{\mu, \Omega}=1}} \inf_{h_h \in \tilde{W}_h} \|h^*(l) - h_h\|_{\mathrm{curl}, \omega, \mu, \zeta, \Omega}, $$

and we set $\gamma_{ba} := \gamma_{ba,E} + \gamma_{ba,H}$.

Notice that $\gamma_{ba} \le \gamma_{st}$, showing that the approximation factor is controlled independently of the mesh size $h$ and the approximation orders $p$. It does, however, in general, depend on the wavenumber $\omega l_\Omega/c_{\Omega,\min}$, the geometry of $\Omega$, and the coefficients $\varepsilon$ and $\mu$. Besides, since our definition of $\gamma_{ba}$ immediately hinges on the “standard” approximation factors $\gamma_{ba,E}$ and $\gamma_{ba,H}$, it automatically recovers their key properties. In particular

$$ (2.13) \qquad \gamma_{ba} \le C(\omega, \Omega, \varepsilon, \mu, s) \left(\frac{h}{p}\right)^s, $$

for some $s \in (0, 1]$, where $p := \min_{K \in \mathcal{T}_h} p_K$. Sharper estimates are available, and we refer the reader to [14, 16, 32, 33, 34, 40], for an in-depth discussion on the dependence of $\gamma_{ba}$ on $\omega$, $h$, and $p$.

2.13. **Coefficient contrasts.** For $K \in \mathcal{T}_h$ and $\phi \in \{\varepsilon, \mu, \zeta, \chi\}$, we employ the notation

$$ \mathcal{C}_{c, \phi, K} := \frac{\max_{\tilde{K}} \phi_{\max}}{\min_{\tilde{K}} \phi_{\min}} $$

for the “contrast” of the coefficient $\phi$ in the patch $\mathcal{T}_{K,h}$. We also set $\mathcal{C}_{c,\phi} := \max_{K \in \mathcal{T}_h} \mathcal{C}_{c,\phi,K}$ and $\mathcal{C}_c := \max_{\phi \in \{\varepsilon, \mu, \zeta, \chi\}} \mathcal{C}_{c,\phi}$. Notice that these quantities are actually independent of the mesh $\mathcal{T}_h$, as long as it fits the partition $\mathcal{P}$ and is only affected by the definition of the coefficients $\varepsilon, \mu$.

2.14. **Notation for generic constants.** In the remaining of this document, if $A, B \ge 0$ are two positive real values, we employ the notation $A \lesssim B$ if there exists a constant $C$ that only depends on $\mathcal{C}_i, \mathcal{C}_{\mathrm{pol}}, \mathcal{C}_b, \mathcal{C}_r$, and $\mathcal{C}_c$ such that $A \le CB$. Importantly, $C$ is independent of $\omega, h$ and $p$. However, $C$ may depend on $\Omega$ through $\mathcal{C}_r$, and it may also depend on $\beta$ through $\mathcal{C}_i$ and $\mathcal{C}_b$. We also employ the notation $A \gtrsim B$ if $B \lesssim A$ and $A \sim B$ if $A \lesssim B$ and $A \gtrsim B$.
### 3. DG DISCRETIZATION AND A POSTERIORI ERROR ESTIMATOR

**3.1. Numerical scheme.** Following [2, 30, 41], the discrete problem consists in finding $(E_h, H_h) \in \mathcal{P}_{k+1}(\mathcal{T}_h) \times \mathcal{P}_{k+1}(\mathcal{T}_h)$ such that

$$ (3.1) \qquad b_h((E_h, H_h), (v_h, w_h)) = i\omega(J, v_h) + \ell_h(v_h, w_h) $$

for all $v_h, w_h \in \mathcal{P}_{k+1}(\mathcal{T}_h)$, where

$$ b_h((E_h, H_h), (v_h, w_h)) := b((E_h, H_h), (v_h, w_h)) + \beta_h((E_h, H_h), (v_h, w_h)), $$

$\beta_h$ is a sesquilinear form over $\mathcal{P}_{k+1}(\mathcal{T}_h) \times \mathcal{P}_{k+1}(\mathcal{T}_h)$ that we call the “flux” form and $\ell_h$ is an antilinear form over $\mathcal{P}_{k+1}(\mathcal{T}_h) \times \mathcal{P}_{k+1}(\mathcal{T}_h)$ designed to impose the jump conditions. We assume that if $e_h, h_h, v_h^\dagger, w_h^\dagger \in \mathcal{P}_{k+1}(\mathcal{T}_h)$, then

$$ (3.2) \qquad \beta_h((e_h, h_h), (v_h^\dagger, w_h^\dagger)) = 0, \qquad \ell_h(v_h^\dagger, w_h^\dagger) = \ell(v_h^\dagger, w_h^\dagger), $$

whenever $v_h^\dagger \in \mathcal{H}_0(\text{curl}, \Omega)$ and $w_h^\dagger \in \mathcal{H}(\text{curl}, \Omega)$. Essentially, we ask for the flux form to vanish for conforming test functions.

In practice, the sesquilinear form $\beta_h(\cdot, \cdot)$ is only employed with discrete arguments to assemble the linear system associated with (3.1). However, in the context of abstract mathematical analysis, it is very convenient to extend the domain of definition of $\beta_h(\cdot, \cdot)$ and to apply it to “continuous” arguments as well. To simplify the discussion, we employ the notation $\mathbf{U} := \mathcal{H}_0(\text{curl}, \Omega) \times \mathcal{H}(\text{curl}, \Omega)$ for the “energy” space of “continuous” functions, $\mathbf{D}_h := \mathcal{P}_{k+1}(\mathcal{T}_h) \times \mathcal{P}_{k+1}(\mathcal{T}_h)$ for the “discrete” space of piecewise polynomial functions, and $\mathbf{S}_h := \mathbf{U} + \mathbf{D}_h$. Because of assumption (3.2), we can consistently extend the domain of definition of $\beta_h(\cdot, \cdot)$ to $\mathbf{D}_h \times \mathbf{D}_h + \mathbf{S}_h \times \mathbf{U}$ by simply setting

$$ (3.3) \qquad \beta_h((e, h), (v^\dagger, w^\dagger)) = 0, \qquad \ell_h(v^\dagger, w^\dagger) = \ell(v^\dagger, w^\dagger), $$

for all $(e, h) \in \mathbf{S}_h$ and $(v^\dagger, w^\dagger) \in \mathbf{U}$.

Then, an immediate consequence of (3.3) is that the discrete form is consistent in the sense that

$$ (3.4) \qquad b_h((e, h), (v, w)) = b((e, h), (v, w)) $$

for all $(e, h) \in \mathbf{S}_h$, and $(v, w) \in \mathbf{U}$. In particular, observing that $\mathbf{W}_h \times \tilde{\mathbf{W}}_h \subset \mathbf{U}$, the Galerkin orthogonality property

$$ (3.5) \qquad b_h((E - E_h, H - H_h), (v_h^\dagger, w_h^\dagger)) = 0 $$

holds true for all $v_h^\dagger \in \mathbf{W}_h$ and $w_h^\dagger \in \tilde{\mathbf{W}}_h$, assuming $(E_h, H_h)$ solves (3.1).

**3.2. Examples of flux form.** After formally multiplying (1.1) by tests function $i\omega\bar{v}_h$, $i\omega\bar{w}_h \in \mathcal{P}_{k+1}(\mathcal{T}_h)$ and integrating by parts locally in each element $K$, one obtains the formulation (3.6)

$$ b((E_h, H_h), (v_h, w_h)) - i\omega \sum_{K \in \mathcal{T}_h} \int_{\partial K} \mathbf{H}_h \cdot \bar{v}_h \times \mathbf{n}_K + i\omega \sum_{K \in \mathcal{T}_h} \int_{\partial K} \mathbf{E}_h \cdot \bar{w}_h \times \mathbf{n}_K = i\omega(J, v_h). $$

Obviously, (3.6) is not a satisfactory discrete formulation since no communication between separate mesh elements occurs, all the considered functions being discontinuous. Following
[2], the solution consists in replacing the traces of $E_h$ and $H_h$ by numerical fluxes $E_h^*$ and $H_h^*$, computed from $E_h$ and $H_h$, leading to (3.7)

$$b((E_h, H_h), (v_h, w_h)) - i\omega \sum_{K \in \mathcal{T}_h} \int_{\partial K} H_h^* \cdot \bar{v}_h \times n_K + i\omega \sum_{K \in \mathcal{T}_h} \int_{\partial K} E_h^* \cdot \bar{w}_h \times n_K = i\omega(J, v_h).$$

If the fluxes are single-valued on every face of the mesh, and if $H_h^* = 0$ on $\partial \Omega$, we may rewrite (3.7) with face-by-face integrals as

$$b((E_h, H_h), (v_h, w_h)) + \beta_h((E_h, H_h), (v_h, w_h)) = i\omega(J, v_h)$$

with

$$(3.8) \quad \beta_h((E_h, H_h), (v_h, w_h)) := -i\omega \sum_{F \in \mathcal{F}_h^i} \int_F H_h^* \cdot \llbracket \bar{v}_h \rrbracket \times n_F + i\omega \sum_{F \in \mathcal{F}_h^i} \int_F E_h^* \cdot \llbracket \bar{w}_h \rrbracket \times n_F.$$

One readily sees that for any single-valued flux, the stabilization form $\beta_h(\cdot, \cdot)$ of (3.8) satisfies (3.2). As a result, the forthcoming analysis applies to a variety of DG schemes.

In particular, a rather general family of fluxes we cover reads

$$(3.9) \quad E_h^* := \frac{1}{\{\llbracket Y \rrbracket\}} \left( \{\llbracket Y E_h \rrbracket\} + \frac{\alpha}{2} \llbracket H_h \rrbracket \times n \right), \quad H_h^* := \frac{1}{\{\llbracket Z \rrbracket\}} \left( \{\llbracket Z H_h \rrbracket\} - \frac{\alpha}{2} \llbracket E_h \rrbracket \times n \right),$$

where $Y, Z$ are piecewise constants weights, and $0 \le \alpha \le 1$, see [48, §3.1.3]. We also refer the reader to [28, 30, 41]. These numerical fluxes are called centered fluxes for $\alpha = 0$ and upwind fluxes for $\alpha = 1$.

### 3.3. Hybridization.
One asset of the scheme associated with any of the fluxes defined by (3.9) is that it is “hybridizable”, which means that it can be equivalently rewritten with a Lagrange multiplier living on the faces of the mesh [25, 30, 38]. Such hybridized form is usually called hybrid discontinuous Galerkin (HDG) and exhibits fewer degrees of freedom than the “naive” formulation (3.1). It is thus well suited to speed up the solve of the associated linear system. Here, for the sake of simplicity, we focus on the formulation (3.1), particularly for symmetry reasons with respect to the analysis of second-order schemes. Notice, however, that since the hybridized system is an equivalent reformulation of (3.1), the proposed estimators apply equally well to HDG discretizations.

### 3.4. Numerical solution.
It is an open question whether the discrete problem (3.1) is well-posed for general meshes. However, the following analysis applies to any pair $(E_h, H_h) \in \mathcal{P}_{k+1}(\mathcal{T}_h) \times \mathcal{P}_{k+1}(\mathcal{T}_h)$ satisfying (3.1), and unique solvability is not required.

We nevertheless mention [25] where the authors analyze (the hybridized version of) the method with upwind fluxes ((3.9) with $\alpha = 1$). They focus on a homogeneous medium enclosed by impedance boundary conditions. In this setting, a key feature of the scheme is that it is well-posed without any assumption on the mesh size. While we work under slightly different assumptions here, we believe that this stability result indicates that the method is interesting for adaptivity techniques since a coarse mesh may be used to start the algorithm.

### 3.5. Error estimators.
For each $K \in \mathcal{T}_h$, the estimator is split into four parts

$$\eta_K^2 := \eta_{\text{div},\varepsilon,K}^2 + \eta_{\text{div},\mu,K}^2 + \eta_{\text{curl},\varepsilon,K}^2 + \eta_{\text{curl},\mu,K}^2,$$
where

$$ \eta_{\text{div},\varepsilon,K} := \frac{1}{\sqrt{\varepsilon_{\tilde{K},\min}}} \left( \frac{h_K}{p_K} \| \nabla \cdot (\mathbf{J} - i\omega\varepsilon\mathbf{E}_h) \|_K + \omega \sqrt{\frac{h_K}{p_K}} \| \llbracket \varepsilon(\mathbf{E}_h - \mathbf{E}^{\text{inc}}) \rrbracket \cdot \mathbf{n} \|_{\partial K \setminus \partial \Omega} \right), $$

$$ \eta_{\text{div},\mu,K} := \frac{1}{\sqrt{\mu_{\tilde{K},\min}}} \left( \omega \frac{h_K}{p_K} \| \nabla \cdot (\mu \mathbf{H}_h) \|_K + \omega \sqrt{\frac{h_K}{p_K}} \| \llbracket \mu(\mathbf{H}_h - \mathbf{H}^{\text{inc}}) \rrbracket \cdot \mathbf{n} \|_{\partial K} \right), $$

$$ \eta_{\text{curl},\varepsilon,K} := \frac{1}{\sqrt{\varepsilon_{\tilde{K},\max}}} \| \mathbf{J} - i\omega\varepsilon\mathbf{E}_h + \nabla \times \mathbf{H}_h \|_K + \sqrt{\mu_{\tilde{K},\max}} \omega \sqrt{\frac{h_K}{p_K}} \| \llbracket (\mathbf{H}_h - \mathbf{H}^{\text{inc}}) \rrbracket \times \mathbf{n} \|_{\partial K \setminus \partial \Omega}, $$

$$ \eta_{\text{curl},\mu,K} := \frac{1}{\sqrt{\mu_{\tilde{K},\max}}} \| i\omega\mu\mathbf{H}_h + \nabla \times \mathbf{E}_h \|_K + \sqrt{\varepsilon_{\tilde{K},\max}} \omega \sqrt{\frac{h_K}{p_K}} \| \llbracket (\mathbf{E}_h - \mathbf{E}^{\text{inc}}) \rrbracket \times \mathbf{n} \|_{\partial K}. $$

We also set $\eta^2 := \eta_{\text{div}}^2 + \eta_{\text{curl}}^2$ with

$$ \eta_{\text{div}}^2 := \sum_{K \in \mathcal{T}_h} (\eta_{\text{div},\varepsilon,K}^2 + \eta_{\text{div},\mu,K}^2), \quad \eta_{\text{curl}}^2 := \sum_{K \in \mathcal{T}_h} (\eta_{\text{curl},\varepsilon,K}^2 + \eta_{\text{curl},\mu,K}^2). $$

## 4. MAIN RESULTS

This section introduces our theoretical results.

### 4.1. Preliminary result for the TF-SF formulation.
We start with a preliminary result concerning the jump term in the right-hand side of the DG formulation.

**Lemma 4.1 (Jump form with gradient arguments).** *The identity*

$$ (4.1) \qquad \ell(\nabla q, \nabla r) = i\varepsilon_0\omega\langle \mathbf{E}^{\text{inc}} \cdot \mathbf{n}_i, q \rangle_{\Gamma_i} + i\mu_0\omega\langle \mathbf{H}^{\text{inc}} \cdot \mathbf{n}_i, r \rangle_{\Gamma_i} $$

holds true for all $(q, r) \in H_0^1(\Omega) \times H^1(\Omega)$.

*Proof.* Let $q \in H_0^1(\Omega)$. We will focus on the $\langle \mathbf{H}^{\text{inc}}, \nabla q \times \mathbf{n}_i \rangle_{\Gamma_i}$ term. The other term in the definition of $\ell$ is treated similarly. We have

$$ \langle \mathbf{H}^{\text{inc}}, \nabla q \times \mathbf{n}_i \rangle_{\Gamma_i} = (\nabla \times \mathbf{H}^{\text{inc}}, \nabla q)_{\Omega_{\text{tf}}} = i\varepsilon_0\omega(\mathbf{E}^{\text{inc}}, \nabla q)_{\Omega_{\text{tf}}} = i\varepsilon_0\omega\langle \mathbf{E}^{\text{inc}} \cdot \mathbf{n}_i, q \rangle_{\Gamma_i}, $$

and (4.1) follows: the last duality pairing may be simplified into a usual $\mathcal{L}^2(\Gamma_i)$ inner-product, as $\mathbf{E}^{\text{inc}} \in \mathcal{L}^2(\Gamma_i)$ by assumption. $\square$

**Remark 4.2 (General jump conditions).** If general jump conditions are employed (i.e., $\nabla \times \mathbf{H}^{\text{inc}} \neq i\omega\mathbf{E}^{\text{inc}}$), then a result similar to (4.1) can still be derived if $\mathbf{E}^{\text{inc}} \cdot \mathbf{n}_i$ and $\mathbf{H}^{\text{inc}} \cdot \mathbf{n}_i$ are replaced by $\text{div}_{\Gamma_i} \mathbf{H}^{\text{inc}}$ and $\text{div}_{\Gamma_i} \mathbf{E}^{\text{inc}}$, where $\text{div}_{\Gamma_i}$ denotes the surface divergence along $\Gamma_i$. A similar modification has to be performed in the estimator as well. Although we lose some generality, we prefer to emphasize the presentation of incident fields solution to free-space Maxwell's equations. Indeed, they are usually employed in practice, and the associated results are easier to understand.
**4.2. Reliability.** We start by two lemmas where we show that the residual is controlled by the estimator.

**Lemma 4.3 (Control of the residual).** *The estimates*

$$ (4.2) \qquad |b_h((E - E_h, H - H_h), (\nabla p, \nabla q))| \lesssim \omega \eta_{\mathrm{div}} (\|\nabla p\|_{\varepsilon, \Omega} + \|\nabla q\|_{\mu, \Omega}) $$

and

$$ (4.3) \quad |b_h((E - E_h, H - H_h), (\phi, \psi))| \lesssim \left( 1 + \max_{K \in \mathcal{T}_h} \frac{\omega h_K}{p_K c_{\tilde{K}, \min}} \right) \eta_{\mathrm{curl}} (\|\nabla \phi\|_{\chi, \Omega} + \|\nabla \psi\|_{\zeta, \Omega}) $$

hold true for all $p \in \mathcal{H}_0^1(\Omega)$, $q \in \mathcal{H}^1(\Omega)$, $\phi \in \mathcal{H}_0^1(\Omega)$ and $\psi \in \mathcal{H}^1(\Omega)$.

**Proof.** We first establish (4.2). We observe that for any $v \in \mathcal{H}_0^1(\Omega)$, $w \in \mathcal{H}^1(\Omega)$, we have

$$ b_h((E - E_h, H - H_h), (\nabla v, \nabla w)) = i\omega(J - i\omega\varepsilon E_h, \nabla v) - i\omega(i\omega\mu H_h, \nabla w) + i\omega l_h(\nabla v, \nabla w) $$

so that, thanks to Lemma 4.1, we get

$$ \begin{aligned} & \frac{1}{i\omega} b_h((E - E_h, H - H_h), (\nabla v, \nabla w)) \\ &= -i\omega \langle \varepsilon E_h \cdot n, v \rangle_{\partial \mathcal{T}_h} - (\nabla \cdot (J - i\omega\varepsilon E_h), v)_{\mathcal{T}_h} - i\omega \langle \mu H_h \cdot n, w \rangle_{\partial \mathcal{T}_h} + i\omega \langle \nabla \cdot (\mu H_h), w \rangle_{\mathcal{T}_h} \\ & \quad + i\varepsilon_0 \omega \langle E^{\mathrm{inc}} \cdot n_i, v \rangle_{\Gamma_i} + i\mu_0 \omega \langle H^{\mathrm{inc}} \cdot n_i, w \rangle_{\Gamma_i} \\ &= -i\omega ( [\varepsilon(E_h - E^{\mathrm{inc}})] \cdot n_i, v )_{\mathcal{F}_h^i} - (\nabla \cdot (J - i\omega\varepsilon E_h), v)_{\mathcal{T}_h} \\ & \quad - i\omega \langle [\mu(H_h - H^{\mathrm{inc}})] \cdot n_i, w \rangle_{\mathcal{F}_h^i} + i\omega (\nabla \cdot (\mu H_h), w)_{\mathcal{T}_h}, \end{aligned} $$

and therefore

$$ \begin{aligned} & (4.4) \quad \frac{1}{\omega} |b_h((E - E_h, H - H_h), (\nabla v, \nabla w))| \\ & \le \omega \sum_{F \in \mathcal{F}_h^i} \| [\varepsilon(E_h - E^{\mathrm{inc}})] \cdot n \|_F \|v\|_F + \omega \sum_{F \in \mathcal{F}_h} \| [\mu(H_h - H^{\mathrm{inc}})] \cdot n \|_F \|w\|_F \\ & \quad + \sum_{K \in \mathcal{T}_h} (\|\nabla \cdot (J - i\omega\varepsilon E_h)\|_K \|v\|_K + \omega \|\nabla \cdot (\mu H_h)\|_K \|w\|_K) \\ & \lesssim \sum_{K \in \mathcal{T}_h} \left( \sqrt{\varepsilon_{\tilde{K}, \min}} \eta_{\mathrm{div}, \varepsilon, K} \left( \frac{p_K}{h_K} \|v\|_K + \sqrt{\frac{p_K}{h_K}} \|v\|_{\partial K \setminus \partial \Omega} \right) \right. \\ & \qquad \left. + \sqrt{\mu_{\tilde{K}, \min}} \eta_{\mathrm{div}, \mu, K} \left( \frac{p_K}{h_K} \|w\|_K + \sqrt{\frac{p_K}{h_K}} \|w\|_{\partial K} \right) \right). \end{aligned} $$
Now, let $p \in \mathcal{H}_0^1(\Omega)$ and $q \in \mathcal{H}^1(\Omega)$. Since $\nabla(\mathcal{Q}_h p) \in \mathbf{W}_h$ and $\nabla(\tilde{\mathcal{Q}}_h q) \in \tilde{\mathbf{W}}_h$, by Galerkin orthogonality (3.5), we can apply (4.4) with $v = p - \mathcal{Q}_h p$ and $w = q - \tilde{\mathcal{Q}}_h q$ to show that

$$
\begin{align*}
& \frac{1}{\omega} |b_h((\mathbf{E} - \mathbf{E}_h, \mathbf{H} - \mathbf{H}_h), (\nabla p, \nabla q))| \\
&= |b_h((\mathbf{E} - \mathbf{E}_h, \mathbf{H} - \mathbf{H}_h), (\nabla(p - \mathcal{Q}_h p), \nabla(q - \tilde{\mathcal{Q}}_h q)))| \\
&\lesssim \sum_{K \in \mathcal{T}_h} \left( \sqrt{\varepsilon_{\tilde{K}, \min}} \eta_{\text{div}, \varepsilon, K} \left( \frac{p_K}{h_K} \|p - \mathcal{Q}_h p\|_K + \sqrt{\frac{p_K}{h_K}} \|p - \mathcal{Q}_h p\|_{\partial K \setminus \partial \Omega} \right) \right. \\
&\qquad \left. + \sqrt{\mu_{\tilde{K}, \min}} \eta_{\text{div}, \mu, K} \left( \frac{p_K}{h_K} \|q - \tilde{\mathcal{Q}}_h q\|_K + \sqrt{\frac{p_K}{h_K}} \|q - \tilde{\mathcal{Q}}_h q\|_{\partial K} \right) \right) \\
&\lesssim \sum_{K \in \mathcal{T}_h} \left( \sqrt{\varepsilon_{\tilde{K}, \min}} \eta_{\text{div}, \varepsilon, K} \|\nabla p\|_{\tilde{K}} + \sqrt{\mu_{\tilde{K}, \min}} \eta_{\text{div}, \mu, K} \|\nabla q\|_{\tilde{K}} \right) \\
&\lesssim \eta_{\text{div}} (\|\nabla p\|_{\varepsilon, \Omega} + \|\nabla q\|_{\mu, \Omega}),
\end{align*}
$$

where we additionally employed (2.5). This shows (4.2).

We now focus on (4.3). Similarly, we start with arbitrary elements $v \in \mathcal{H}^1(\mathcal{T}_h) \cap \mathcal{H}_0(\text{curl}, \Omega)$ and $w \in \mathcal{H}^1(\mathcal{T}_h) \cap \mathcal{H}(\text{curl}, \Omega)$. We have

$$
\begin{align*}
b_h((\mathbf{E} - \mathbf{E}_h, \mathbf{H} - \mathbf{H}_h), (\mathbf{v}, \mathbf{w}))
&= (i\omega \mathbf{J}, \mathbf{v}) + (\omega^2 \varepsilon \mathbf{E}_h, \mathbf{v}) + (i\omega \mathbf{H}_h, \nabla \times \mathbf{v}) + (\omega^2 \mu \mathbf{H}_h, \mathbf{w}) - (i\omega \mathbf{E}_h, \nabla \times \mathbf{w}) + i\omega \ell_h(\mathbf{v}, \mathbf{w}) \\
&= (i\omega \mathbf{J} + \omega^2 \varepsilon \mathbf{E}_h + i\omega \nabla \times \mathbf{H}_h, \mathbf{v})_{\mathcal{T}_h} - i\omega \langle \mathbf{H}_h \times \mathbf{n}, \mathbf{v} \rangle_{\partial \mathcal{T}_h} + \langle \mathbf{H}^{\text{inc}}, \mathbf{v} \times \mathbf{n}_i \rangle_{\Gamma_i} \\
&\quad - i\omega(i\omega \mu \mathbf{H}_h + \nabla \times \mathbf{E}_h, \mathbf{w}) + i\omega \langle \mathbf{E}_h \times \mathbf{n}, \mathbf{w} \rangle_{\partial \mathcal{T}_h} - \langle \mathbf{E}^{\text{inc}}, \mathbf{w} \times \mathbf{n}_i \rangle_{\Gamma_i} \\
&= i\omega(\mathbf{J} - i\omega \varepsilon \mathbf{E}_h + \nabla \times \mathbf{H}_h, \mathbf{v})_{\mathcal{T}_h} - i\omega \langle [\mathbf{H}_h - \mathbf{H}^{\text{inc}}] \times \mathbf{n}, \mathbf{v} \rangle_{\mathcal{F}_h^i} \\
&\quad - i\omega(i\omega \mu \mathbf{H}_h + \nabla \times \mathbf{E}_h, \mathbf{w}) + i\omega \langle [\mathbf{E}_h - \mathbf{E}^{\text{inc}}] \times \mathbf{n}, \mathbf{w} \rangle_{\mathcal{F}_h^i}
\end{align*}
$$

and therefore

$$
\begin{align*}
(4.5) \quad |b_h((\mathbf{E} - \mathbf{E}_h, \mathbf{H} - \mathbf{H}_h), (\mathbf{v}, \mathbf{w}))|
&\le \sum_{K \in \mathcal{T}_h} (\omega \|\mathbf{J} - i\omega \varepsilon \mathbf{E}_h + \nabla \times \mathbf{H}_h\|_K \|\mathbf{v}\|_K + \omega \|i\omega \mu \mathbf{H}_h + \nabla \times \mathbf{E}_h\|_K \|\mathbf{w}\|_K) \\
&\quad + \omega \sum_{F \in \mathcal{F}_h^i} \|\llbracket \mathbf{H}_h - \mathbf{H}^{\text{inc}} \rrbracket \times \mathbf{n}\|_F \|\mathbf{v} \times \mathbf{n}\|_F + \omega \sum_{F \in \mathcal{F}_h^i} \|\llbracket \mathbf{E}_h - \mathbf{E}^{\text{inc}} \rrbracket \times \mathbf{n}\|_F \|\mathbf{w} \times \mathbf{n}\|_F \\
&\lesssim \sum_{K \in \mathcal{T}_h} \left( 1 + \frac{\omega h_K}{p_K c_{\tilde{K}, \min}} \right) \left( \frac{1}{\sqrt{\mu_{\tilde{K}, \max}}} \eta_{\text{curl}, \varepsilon, K} \left( \frac{p_K}{h_K} \|\mathbf{v}\|_K + \sqrt{\frac{p_K}{h_K}} \|\mathbf{v} \times \mathbf{n}\|_{\partial K \setminus \partial \Omega} \right) \right. \\
&\qquad \left. + \frac{1}{\sqrt{\varepsilon_{\tilde{K}, \max}}} \eta_{\text{curl}, \mu, K} \left( \frac{p_K}{h_K} \|\mathbf{w}\|_K + \sqrt{\frac{p_K}{h_K}} \|\mathbf{w} \times \mathbf{n}\|_{\partial K} \right) \right).
\end{align*}
$$

Let now $\phi \in \mathcal{H}_0^1(\Omega)$ and $\psi \in \mathcal{H}^1(\Omega)$. Since $\mathcal{R}_h \phi \in \mathbf{W}_h$ and $\tilde{\mathcal{R}}_h \psi \in \tilde{\mathbf{W}}_h$, by Galerkin orthogonality (3.5), we may employ (4.5) with $v = \phi - \mathcal{R}_h \phi$, $w = \psi - \tilde{\mathcal{R}}_h \psi$ and (2.6), showing that

$$
|b_h((\mathbf{E} - \mathbf{E}_h, \mathbf{H} - \mathbf{H}_h), (\phi, \psi))| = |b_h((\mathbf{E} - \mathbf{E}_h, \mathbf{H} - \mathbf{H}_h), (\phi - \mathcal{R}_h \phi, \psi - \tilde{\mathcal{R}}_h \psi))|
$$
$$
\begin{align*}
&\lesssim \sum_{K \in \mathcal{T}_h} \left( 1 + \frac{\omega h_K}{p_K c_{\tilde{K}, \min}} \right) \left( \frac{1}{\sqrt{\mu_{\tilde{K}, \max}}} \eta_{\mathrm{curl}, \varepsilon, K} \left( \frac{p_K}{h_K} \|\phi - \mathcal{R}_h \phi\|_K + \sqrt{\frac{p_K}{h_K}} \|(\phi - \mathcal{R}_h \phi) \times \mathbf{n}\|_{\partial K \setminus \partial \Omega} \right) \right. \\
&\qquad \left. + \frac{1}{\sqrt{\varepsilon_{\tilde{K}, \max}}} \eta_{\mathrm{curl}, \mu, K} \left( \frac{p_K}{h_K} \|\psi - \tilde{\mathcal{R}}_h \psi\|_K + \sqrt{\frac{p_K}{h_K}} \|(\psi - \tilde{\mathcal{R}}_h \psi) \times \mathbf{n}\|_{\partial K} \right) \right) \\
&\lesssim \sum_{K \in \mathcal{T}_h} \left( 1 + \frac{\omega h_K}{p_K c_{\tilde{K}, \min}} \right) \left( \eta_{\mathrm{curl}, \varepsilon, K} \|\nabla \phi\|_{\chi, \tilde{K}} + \eta_{\mathrm{curl}, \mu, K} \|\nabla \psi\|_{\zeta, \tilde{K}} \right) \\
&\lesssim \left( 1 + \max_{K \in \mathcal{T}_h} \frac{\omega h_K}{p_K c_{\tilde{K}, \min}} \right) \eta_{\mathrm{curl}} (\|\nabla \phi\|_{\chi, \Omega} + \|\nabla \psi\|_{\zeta, \Omega}).
\end{align*}
\hfill $\square$

**Lemma 4.4** (General control of the residual). We have

$$
(4.6) \quad |b_h((E - E_h, H - H_h), (\boldsymbol{\theta}, \boldsymbol{\vartheta}))| \lesssim \left( 1 + \max_{K \in \mathcal{T}_h} \frac{\omega h_K}{p_K c_{\tilde{K}, \min}} \right) \eta \|\|(\boldsymbol{\theta}, \boldsymbol{\vartheta})\|_{\mathrm{curl}, \omega, \Omega}
$$
for all $\boldsymbol{\theta} \in \mathcal{H}_0(\mathrm{curl}, \Omega)$ and $\boldsymbol{\vartheta} \in \mathcal{H}(\mathrm{curl}, \Omega)$.

Proof. Using the results from Section 2.10, we can decompose, any $\boldsymbol{\theta} \in \mathcal{H}_0(\mathrm{curl}, \Omega)$ and $\boldsymbol{\vartheta} \in \mathcal{H}(\mathrm{curl}, \Omega)$ as $\boldsymbol{\theta} = \boldsymbol{\phi} + \nabla r$ and $\boldsymbol{\vartheta} = \boldsymbol{\psi} + \nabla s$ with $(\boldsymbol{\phi}, \boldsymbol{\psi}) \in \mathcal{H}_0^1(\Omega) \times \mathcal{H}^1(\Omega)$ and $(r, s) \in \mathcal{H}_0^1(\Omega) \times \mathcal{H}^1(\Omega)$. By employing (4.2) and (4.3) from Lemma 4.3, we have

$$
\begin{align*}
|b_h((E - E_h, H - H_h), (\boldsymbol{\theta}, \boldsymbol{\vartheta}))| &\le |b_h((E - E_h, H - H_h), (\boldsymbol{\phi}, \boldsymbol{\psi}))| \\
&\quad + |b_h((E - E_h, H - H_h), (\nabla r, \nabla s))| \\
&\lesssim \eta_{\mathrm{curl}} \left( 1 + \max_{K \in \mathcal{T}_h} \frac{\omega h_K}{p_K c_{\tilde{K}, \min}} \right) (\|\nabla \boldsymbol{\phi}\|_{\chi, \Omega} + \|\nabla \boldsymbol{\psi}\|_{\zeta, \Omega}) \\
&\quad + \eta_{\mathrm{div}} \omega (\|\nabla r\|_{\varepsilon, \Omega} + \|\nabla s\|_{\mu, \Omega}),
\end{align*}
$$

and we conclude recalling that $\eta^2 := \eta_{\mathrm{curl}}^2 + \eta_{\mathrm{div}}^2$ since (2.8) imply that

$$
\|\nabla \phi\|_{\chi, \Omega} + \|\nabla \psi\|_{\zeta, \Omega} \lesssim \|(\theta, \vartheta)\|_{\mathrm{curl}, \omega, \Omega}, \quad \omega (\|\nabla p\|_{\varepsilon, \Omega} + \|\nabla s\|_{\mu, \Omega}) \lesssim \|(\theta, \vartheta)\|_{\mathrm{curl}, \omega, \Omega}.
$$

$\square$

The next step is an Aubin-Nitsche type result that controls the $\mathcal{L}^2(\Omega)$-norm of the error
to make up for the lack of coercivity of the sesquilinear form $b$. To this end, we first state a
result concerning the approximation factor for first-order schemes.

**Lemma 4.5 (Approximation factor).** For all $\mathbf{j}, \mathbf{l} \in \mathcal{L}^2(\Omega)$, there exists a unique pair $(\mathbf{e}^*, \mathbf{h}^*)(\mathbf{j}, \mathbf{l}) \in \mathcal{H}_0(\mathrm{curl}, \Omega) \times \mathcal{H}(\mathrm{curl}, \Omega)$ such that

$$
b((\mathbf{v}, \mathbf{w}), (\mathbf{e}^*, \mathbf{h}^*)(\mathbf{j}, \mathbf{l})) = \omega(\mathbf{v}, \varepsilon\mathbf{j}) + \omega(\mathbf{w}, \mu\mathbf{l}).
$$

In addition, if $\mathbf{j} \in \mathcal{H}(\mathrm{div}^0, \varepsilon, \Omega)$ and $\mathbf{l} \in \mathcal{H}_0(\mathrm{div}^0, \mu, \Omega)$, we have

$$
(4.7) \quad \inf_{\substack{e_h \in \tilde{W}_h \\ h_h \in \tilde{W}_h}} \|(e^*, h^*)(\mathbf{j}, \mathbf{l}) - (e_h, h_h)\|_{\mathrm{curl}, \omega, \Omega} \lesssim (1 + \gamma_{\mathrm{ba}}) (\|\mathbf{j}\|_{\varepsilon, \Omega} + \|\mathbf{l}\|_{\mu, \Omega}).
$$
Proof. Let $j, l \in \mathcal{L}^2(\Omega)$ and set $(e, h) := (e^*, h^*)(j, l)$. We first observe that

$$-i\omega\nabla \times h = \omega\epsilon j + \omega^2\bar{\epsilon}e, \quad i\omega\nabla \times e = \omega\mu l + \omega^2\bar{\mu}h.$$

Then, considering $\phi \in \mathcal{H}(\text{curl}, \Omega)$ and $\psi \in \mathcal{H}_0(\text{curl}, \Omega)$, selecting the test functions $v := \zeta\nabla \times \phi$ and $w := \chi\nabla \times \psi$, and integrating by parts, we show that

$$b_E(\psi, e) = \omega(\psi, \epsilon j) + i(\bar{\mu}\chi\nabla \times \psi, l), \quad b_H(\phi, h) = \omega(\phi, \mu l) - i(\bar{\epsilon}\zeta\nabla \times \phi, j),$$

for all $\psi \in \mathcal{H}_0(\text{curl}, \Omega)$ and $\phi \in \mathcal{H}(\text{curl}, \Omega)$. At this point, it is tempting to use the approximation factors $\gamma_{\text{ba,E}}$ and $\gamma_{\text{ba,H}}$. However, recalling their definition in (2.12), it is not possible yet, since the right-hand sides are not in $\mathcal{L}^2(\Omega)$. The key idea then consists in “lifting” the last term in the above identities. To do so, we introduce $e_0$ and $h_0$ as the unique elements of $\mathcal{H}_0(\text{curl}, \Omega)$ and $\mathcal{H}(\text{curl}, \Omega)$ such that

$$2\omega^2(\epsilon\psi, e_0) + b_E(\psi, e_0) = i(\bar{\mu}\chi\nabla \times \psi, l), \quad 2\omega^2(\mu\phi, h_0) + b_H(\phi, h_0) = -i(\bar{\epsilon}\zeta\nabla \times \phi, j),$$

for all $\psi \in \mathcal{H}_0(\text{curl}, \Omega)$ and $\phi \in \mathcal{H}(\text{curl}, \Omega)$. As can be seen from (2.3), the left-hand sides correspond to coercive sesquilinear forms, and we have

$$\|e_0\|_{\text{curl},\omega,\epsilon,\chi,\Omega}^2 = \text{Re} \left( 2\omega^2(\epsilon e_0, e_0) + b_E(e_0, e_0) \right) = \text{Re} i(\bar{\mu}\chi\nabla \times e_0, l) \\ \lesssim \|\nabla \times e_0\|_{\chi,\Omega} \|l\|_{\mu,\Omega} \le \|l\|_{\mu,\Omega} \|e_0\|_{\text{curl},\omega,\epsilon,\chi,\Omega}.$$

As a result, we have $\|e_0\|_{\text{curl},\omega,\epsilon,\chi,\Omega} \lesssim \|l\|_{\mu,\Omega}$. Similar arguments show that $\|h_0\|_{\text{curl},\omega,\mu,\zeta,\Omega} \lesssim \|j\|_{\epsilon,\Omega}$, and therefore

$$\|(e_0, h_0)\|_{\text{curl},\omega,\Omega} \lesssim \|j\|_{\epsilon,\Omega} + \|l\|_{\mu,\Omega}.$$

On the other hand, we see that

$$b_E(\psi, e_0) = i(\bar{\mu}\chi\nabla \times \psi, l) - 2\omega^2(\epsilon\psi, e_0), \quad b_H(\phi, h_0) = -i(\bar{\epsilon}\zeta\nabla \times \phi, j) - 2\omega^2(\mu\phi, h_0),$$

and therefore, letting $(\tilde{e}, \tilde{h}) := (e, h) - (e_0, h_0)$, we have

$$b_E(\psi, \tilde{e}) = \omega(\psi, \epsilon j) + 2\omega^2(\epsilon\psi, e_0) = \omega(\psi, \epsilon \tilde{j}),$$

$$b_H(\phi, \tilde{h}) = \omega(\phi, \mu l) + 2\omega^2(\mu\phi, h_0) = \omega(\phi, \mu \tilde{l}),$$

with $(\tilde{j}, \tilde{l}) := (j, l) + 2\omega(\epsilon^{-1}\bar{\epsilon}e_0, \mu^{-1}\bar{\mu}h_0)$.

Now, we observe that picking a gradient as a test function in the definition of $e_0$ and $h_0$ reveals that $e_0 \in \mathcal{H}(\text{div}^0, \bar{\epsilon}, \Omega)$ and $h_0 \in \mathcal{H}_0(\text{div}^0, \bar{\mu}, \Omega)$. Hence $\epsilon^{-1}\bar{\epsilon}e_0 \in \mathcal{H}(\text{div}^0, \epsilon, \Omega)$ and $\mu^{-1}\bar{\mu}h_0 \in \mathcal{H}(\text{div}^0, \mu, \Omega)$. As a result, we have

$$\inf_{e_h \in \mathcal{W}_h} \|\tilde{e} - e_h\|_{\text{curl},\omega,\epsilon,\chi,\Omega} \le \gamma_{\text{ba,E}} \|\tilde{j}\|_{\epsilon,\Omega} \lesssim \gamma_{\text{ba,E}} (\|j\|_{\epsilon,\Omega} + 2\|e_0\|_{\text{curl},\omega,\epsilon,\chi,\Omega}) \\ \lesssim \gamma_{\text{ba,E}} (\|j\|_{\epsilon,\Omega} + \|l\|_{\mu,\Omega}),$$

and similarly

$$\inf_{h_h \in \mathcal{W}_h} \|\tilde{h} - h_h\|_{\text{curl},\omega,\mu,\zeta,\Omega} \lesssim \gamma_{\text{ba,H}} (\|j\|_{\epsilon,\Omega} + \|l\|_{\mu,\Omega}).$$
Now, (4.7) follows since

$$
\begin{align*}
\inf_{\substack{e_h \in \mathcal{W}_h \\ h_h \in \tilde{\mathcal{W}}_h}} \|(e, h) - (e_h, h_h)\|_{\mathrm{curl}, \omega, \Omega} &= \inf_{\substack{e_h \in \mathcal{W}_h \\ h_h \in \tilde{\mathcal{W}}_h}} \|(e_0, h_0) + (\tilde{e}, \tilde{h}) - (e_h, h_h)\|_{\mathrm{curl}, \omega, \Omega} \\
&\le \|(e_0, h_0)\|_{\mathrm{curl}, \omega, \Omega} + \inf_{\substack{e_h \in \mathcal{W}_h \\ h_h \in \tilde{\mathcal{W}}_h}} \|(\tilde{e}, \tilde{h}) - (e_h, h_h)\|_{\mathrm{curl}, \omega, \Omega} \\
&= \|(e_0, h_0)\|_{\mathrm{curl}, \omega, \Omega} + \inf_{e_h \in \mathcal{W}_h} \|\tilde{e} - e_h\|_{\mathrm{curl}, \omega, \varepsilon, \chi, \Omega} + \inf_{h_h \in \tilde{\mathcal{W}}_h} \|\tilde{h} - h_h\|_{\mathrm{curl}, \omega, \mu, \zeta, \Omega} \\
&\lesssim (1 + \gamma_{\mathrm{ba}, \mathrm{E}} + \gamma_{\mathrm{ba}, \mathrm{H}}) (\|\tilde{j}\|_{\varepsilon, \Omega} + \|l\|_{\mu, \Omega}) = (1 + \gamma_{\mathrm{ba}}) (\|j\|_{\varepsilon, \Omega} + \|l\|_{\mu, \Omega}).
\end{align*}
$$

□

**Lemma 4.6** (Aubin-Nitsche). We have

$$
\omega \| \mathbf{E} - \mathbf{E}_h \|_{\varepsilon, \Omega} + \omega \| \mathbf{H} - \mathbf{H}_h \|_{\mu, \Omega} \lesssim \left( 1 + \max_{K \in \mathcal{T}_h} \frac{\omega h_K}{p_K c_{\tilde{K}, \min}} \right) (1 + \gamma_{\text{ba}}) \eta.
$$

*Proof.* The proof relies on the Helmholtz decomposition of the error. We thus define $p \in \mathcal{H}_0^1(\Omega)$ and $q \in \mathcal{H}^1(\Omega)$ such that

$$
(\varepsilon \nabla p, \nabla v) = (\varepsilon(E - E_h), \nabla v), \quad (\mu \nabla q, \nabla w) = (\mu(H - H_h), \nabla w),
$$

for all $v \in \mathcal{H}_0^1(\Omega)$ and $w \in \mathcal{H}^1(\Omega)$. Notice that $p$ is uniquely defined and that $q$ is defined up to a constant, that does not contribute to its gradient. Then, we have

$$
E - E_h = \nabla p + \theta, \quad H - H_h = \nabla q + \vartheta,
$$

with $p \in \mathcal{H}_0^1(\Omega)$, $q \in \mathcal{H}^1(\Omega)$, $\theta \in \mathcal{H}(\mathrm{div}^0, \varepsilon, \Omega)$ and $\vartheta \in \mathcal{H}_0(\mathrm{div}^0, \mu, \Omega)$. For the gradient terms, we have

$$
\begin{align*}
\omega^2 \|\nabla p\|_{\varepsilon, \Omega}^2 + \omega^2 \|\nabla q\|_{\mu, \Omega}^2 &= \omega^2 \mathrm{Re}(\varepsilon(E - E_h), \nabla p) + \omega^2 \mathrm{Re}(\mu(H - H_h), \nabla q) \\
&= - \mathrm{Re} b_h((E - E_h, H - H_h), (\nabla p, \nabla q)) \lesssim \omega \eta_{\mathrm{div}} (\|\nabla p\|_{\varepsilon, \Omega} + \|\nabla q\|_{\mu, \Omega}),
\end{align*}
$$

so that

$$
(4.8) \qquad \omega \| \nabla p \|_{\varepsilon, \Omega} + \omega \| \nabla q \|_{\mu, \Omega} \lesssim \eta_{\mathrm{div}}.
$$

For the remaining terms, we observe that

$$
\begin{align*}
\omega \|\boldsymbol{\theta}\|_{\varepsilon, \Omega}^2 + \omega \|\boldsymbol{\vartheta}\|_{\mu, \Omega}^2 &= \mathrm{Re}(\omega(\boldsymbol{\varepsilon}\boldsymbol{\theta}, \boldsymbol{\theta}) + \omega(\boldsymbol{\mu}\boldsymbol{\vartheta}, \boldsymbol{\vartheta})) = \mathrm{Re}(\omega(\boldsymbol{\theta}, \boldsymbol{\varepsilon}\boldsymbol{\theta}) + \omega(\boldsymbol{\vartheta}, \boldsymbol{\mu}\boldsymbol{\vartheta})) \\
&= \mathrm{Re}(\omega(E - E_h, \boldsymbol{\varepsilon}\boldsymbol{\theta}) + \omega(H - H_h, \boldsymbol{\mu}\boldsymbol{\vartheta})).
\end{align*}
$$

Then, by Lemma 4.5, we may define $(\xi, \zeta)$ as the unique element of $\mathcal{H}_0(\mathrm{curl}, \Omega) \times \mathcal{H}(\mathrm{curl}, \Omega)$ such that

$$
b((w, v), (\xi, \zeta)) = \omega(w, \varepsilon\theta) + \omega(v, \mu\vartheta)
$$

for all $w, v \in \mathcal{L}^2(\Omega)$. Using consistency property (3.4), Galerkin orthogonality (3.5) and (4.6), we have

$$
\begin{align*}
\mathrm{Re}(\omega((E - E_h), \varepsilon\theta) + \omega((H - H_h), \mu\vartheta)) &= \mathrm{Re} b((E - E_h, H - H_h), (\xi, \zeta)) \\
&= \mathrm{Re} b_h((E - E_h, H - H_h), (\xi, \zeta)) = \mathrm{Re} b_h((E - E_h, H - H_h), (\xi - \xi_h, \zeta - \zeta_h)) \\
&\lesssim \left( 1 + \max_{K \in \mathcal{T}_h} \frac{\omega h_K}{p_K c_{\tilde{K}, \min}} \right) \eta \|(\xi - \xi_h, \zeta - \zeta_h)\|_{\mathrm{curl}, \omega, \Omega}
\end{align*}
$$
for all $\xi_h \in W_h$ and $\zeta_h \in \tilde{W}_h$. Then, recalling (4.7), we deduce that

$$ \mathrm{Re}(\omega((E - E_h), \varepsilon\theta) + \omega((H - H_h), \mu\vartheta)) \\ \lesssim \left(1 + \max_{K \in \mathcal{T}_h} \frac{\omega h_K}{p_K c_{\tilde{K}, \min}}\right) (1 + \gamma_{\mathrm{ba}}) \eta (\|\theta\|_{\varepsilon, \Omega} + \|\vartheta\|_{\mu, \Omega}). $$

Hence,

$$ \omega\|\theta\|_{\varepsilon, \Omega}^2 + \omega\|\vartheta\|_{\mu, \Omega}^2 \lesssim \left(1 + \max_{K \in \mathcal{T}_h} \frac{\omega h_K}{p_K c_{\tilde{K}, \min}}\right) (1 + \gamma_{\mathrm{ba}}) \eta (\|\theta\|_{\varepsilon, \Omega} + \|\vartheta\|_{\mu, \Omega}), $$

and the result follows from (4.8). □

We are now ready to establish the main result of this section. Notice that in contrast to second-order schemes [15], the estimate does not stem from a continuous-level Gårding inequality. The “electric-magnetic mismatch” part of the estimator is employed instead.

**Theorem 4.7 (Reliability).** *The estimate*

$$ (4.9) \qquad \|(E - E_h, H - H_h)\|_{\mathrm{curl}, \omega, \mathcal{T}_h} \lesssim \left(1 + \max_{K \in \mathcal{T}_h} \frac{\omega h_K}{p_K c_{\tilde{K}, \min}}\right) (1 + \gamma_{\mathrm{ba}}) \eta $$

holds true.

*Proof.* We start by the observation that

$$ \begin{aligned} \nabla \times (E - E_h) &= -i\omega\mu(H - H_h) - (i\omega\mu H_h + \nabla \times E_h), \\ \nabla \times (H - H_h) &= i\omega\varepsilon(E - E_h) - (J - i\omega\varepsilon E_h + \nabla \times H_h), \end{aligned} $$

which immediately yields the estimates

$$ \begin{aligned} \|\nabla \times (E - E_h)\|_{\chi, \Omega} &\lesssim \omega \|H - H_h\|_{\mu, \Omega} + \sum_{K \in \mathcal{T}_h} \frac{1}{\sqrt{\mu_{\tilde{K}, \max}}} \|i\omega\mu H_h + \nabla \times E_h\|_K, \\ \|\nabla \times (H - H_h)\|_{\zeta, \Omega} &\lesssim \omega \|E - E_h\|_{\varepsilon, \Omega} + \sum_{K \in \mathcal{T}_h} \frac{1}{\sqrt{\varepsilon_{\tilde{K}, \max}}} \|J - i\omega\varepsilon E_h + \nabla \times H_h\|_K. \end{aligned} $$

Adding the last two inequalities, we have

$$ \|\nabla \times (E - E_h)\|_{\chi, \Omega} + \|\nabla \times (H - H_h)\|_{\zeta, \Omega} \lesssim \omega \|E - E_h\|_{\varepsilon, \Omega} + \omega \|H - H_h\|_{\mu, \Omega} + \eta_{\mathrm{curl}}, $$

and (4.9) follows since we already estimated the $\mathcal{L}^2(\Omega)$ terms in Lemma 4.6. □

**Remark 4.8 (Asymptotic estimate for smooth solutions).** If $\varepsilon$ and $\mu$ are real scalars, we can actually rewrite that last line of the above proof as

$$ \|\nabla \times (E - E_h)\|_{\chi, \Omega} + \|\nabla \times (H - H_h)\|_{\zeta, \Omega} \le \omega \|E - E_h\|_{\varepsilon, \Omega} + \omega \|H - H_h\|_{\mu, \Omega} + \eta_{\mathrm{curl}}, $$

without any hidden constant. As a result, assuming that the solution is sufficiently smooth (or that the mesh is locally refined) so that the convergence in $\mathcal{L}^2$ norm happens faster than in the energy norm, we asymptotically have

$$ \|(E - E_h, H - H_h)\|_{\mathrm{curl}, \omega, \Omega} \le \eta_{\mathrm{curl}} = \eta. $$

This behavior is observed several times in the numerical examples reported hereafter.
4.3. **Efficiency.** We now show that the estimator proposed for DG discretizations is efficient. Classically, the proofs of this section hinge on the “bubble” functions introduced in Section 2.8.

We start by showing an upper bound for the “divergence” parts of the estimator, namely $\eta_{\text{div},\varepsilon,K}$ and $\eta_{\text{div},\mu,K}$.

**Lemma 4.9.** We have

$$ \eta_{\text{div},\varepsilon,K} \lesssim p_K^{3/2} \omega \|\mathbf{E} - \mathbf{E}_h\|_{\varepsilon,\tilde{K}} + \text{osc}_{\mathcal{T}_{K,h}}, \quad \eta_{\text{div},\mu,K} \lesssim p_K^{3/2} \omega \|\mathbf{H} - \mathbf{H}_h\|_{\mu,\tilde{K}} + \text{osc}_{\mathcal{T}_{K,h}} $$
for all $K \in \mathcal{T}_h$.

*Proof.* For the sake of readability, we make a slight abuse of notation in the proof, and set $[\mathbf{E}_h^{\text{inc}}]_F := \mathbf{E}_h^{\text{inc}}$ if $F \subset \Gamma_i$ and $[\mathbf{E}_h^{\text{inc}}]_F = 0$ for the remaining faces.

Let $K \in \mathcal{T}_h$ and $v_K := b_K \nabla \cdot (\mathbf{J}_h - i\omega\varepsilon\mathbf{E}_h)$. After integration by parts, we have

$$ \begin{aligned} \|b_K^{1/2} \nabla \cdot (\mathbf{J}_h - i\omega\varepsilon\mathbf{E}_h)\|_K^2 &= (\nabla \cdot (\mathbf{J}_h - i\omega\varepsilon\mathbf{E}_h), v_K)_K \\ &= -(i\omega\varepsilon(\mathbf{E} - \mathbf{E}_h), \nabla v_K)_K - (\nabla \cdot (\mathbf{J} - \mathbf{J}_h), v_K)_K \\ &\le \omega \|\varepsilon(\mathbf{E} - \mathbf{E}_h)\|_K \|\nabla v_K\|_K + \|b_K^{1/2} \nabla \cdot (\mathbf{J} - \mathbf{J}_h)\|_K \|b_K^{1/2} \nabla \cdot (\mathbf{J}_h - i\omega\varepsilon\mathbf{E}_h)\|_K. \end{aligned} $$

Using (2.7a) and (2.7b), it follows that

$$ (4.10) \qquad \frac{1}{\sqrt{\varepsilon_{\tilde{K},\min}}} \frac{h_K}{p_K} \|\nabla \cdot (\mathbf{J} - i\omega\varepsilon\mathbf{E}_h)\|_K \lesssim p_K \omega \|\mathbf{E} - \mathbf{E}_h\|_{\varepsilon,K} + \frac{p_K}{\sqrt{\varepsilon_{\tilde{K},\min}}} \frac{h_K}{p_K} \|\nabla \cdot (\mathbf{J} - \mathbf{J}_h)\|_K. $$

On the other hand, for $F \in \mathcal{F}_h^i \cap \mathcal{F}_K$, we set $v_F := \mathcal{E}(\omega [[\varepsilon(\mathbf{E}_h - \mathbf{E}_h^{\text{inc}})] \cdot \mathbf{n}_F)$. Since $i\omega\nabla \cdot (\varepsilon\mathbf{E}) = \nabla \cdot \mathbf{J}$ and $b_F = 0$ on $\partial\tilde{F}$, we have

$$ \omega^2 \|b_F^{1/2} [[\varepsilon(\mathbf{E}_h - \mathbf{E}_h^{\text{inc}})] \cdot \mathbf{n}_F\|_F^2 = \omega( [[\varepsilon(\mathbf{E}_h - \mathbf{E})] \cdot \mathbf{n}_F, v_F)_F + \omega( [[\varepsilon(\mathbf{E}^{\text{inc}} - \mathbf{E}_h^{\text{inc}})] \cdot \mathbf{n}_F, v_F)_F, \text{ and} $$

$$ \begin{aligned} \omega( [[\varepsilon(\mathbf{E}_h - \mathbf{E})] \cdot \mathbf{n}_F, v_F)_F &= |(\nabla \cdot (i\omega\varepsilon\mathbf{E}_h), v_F)_{\mathcal{T}_{F,h}} + (i\omega\varepsilon\mathbf{E}_h, \nabla v_F)_{\mathcal{T}_{F,h}}| \\ &= | - (\nabla \cdot (\mathbf{J} - i\omega\varepsilon\mathbf{E}_h), v_F)_{\mathcal{T}_{F,h}} - (i\omega\varepsilon(\mathbf{E} - \mathbf{E}_h), \nabla v_F)_{\mathcal{T}_{F,h}}| \\ &\le \|\nabla \cdot (\mathbf{J} - i\omega\varepsilon\mathbf{E}_h)\|_{\mathcal{T}_{F,h}} \|v_F\|_{\tilde{F}} + \omega \|\varepsilon(\mathbf{E} - \mathbf{E}_h)\|_{\tilde{F}} \|\nabla v_F\|_{\tilde{F}}. \end{aligned} $$

Then, thanks to (2.7c), we bound the terms depending on $v_F$, and using both (2.7a) and (2.4), we get

$$ \begin{aligned} (4.11) \qquad & \frac{\omega}{\sqrt{\varepsilon_{\tilde{K},\min}}} \sqrt{\frac{h_K}{p_K}} \| [[\varepsilon(\mathbf{E}_h - \mathbf{E}_h^{\text{inc}})] \cdot \mathbf{n}_F \|_F \\ & \lesssim \frac{\omega}{\sqrt{\varepsilon_{\tilde{K},\min}}} \sqrt{\frac{h_K}{p_K}} p_F \| b_F^{1/2} [[\varepsilon(\mathbf{E}_h - \mathbf{E}_h^{\text{inc}})] \cdot \mathbf{n}_F \|_F \\ & \lesssim \frac{p_K^{1/2}}{\sqrt{\varepsilon_{\tilde{K},\min}}} \frac{h_K}{p_K} \| \nabla \cdot (\mathbf{J} - i\omega\varepsilon\mathbf{E}_h) \|_{\mathcal{T}_{F,h}} + p_K^{1/2} \omega \| \mathbf{E} - \mathbf{E}_h \|_{\varepsilon,\tilde{F}} + \text{osc}_{\mathcal{T}_{F,h}} (\mathbf{E}^{\text{inc}}). \end{aligned} $$

The bound associated with $\eta_{\text{div},\varepsilon,K}$ follows from (4.10) and (4.11). For the sake of shortness, we do not write the proofs for $\eta_{\text{div},\mu,K}$, as it follows from the same arguments. $\square$

We now turn to the “rotation” parts of the estimator, which require increased attention.
**Lemma 4.10.** We have

(4.12a)

$$ \eta_{\text{curl},\varepsilon,K} \lesssim \omega \|\mathbf{E} - \mathbf{E}_h\|_{\varepsilon,K} + p_K \left(1 + \frac{\omega h_K}{p_K c_{\tilde{K},\min}}\right) \|\mathbf{H} - \mathbf{H}_h\|_{\text{curl},\omega,\mu,\zeta,\mathcal{T}_{K,h}} + \text{osc}_{\mathcal{T}_{F,h}}(\mathbf{H}^{\text{inc}}) $$

and

(4.12b)

$$ \eta_{\text{curl},\mu,K} \lesssim p_K \left(1 + \frac{\omega h_K}{p_K c_{\tilde{K},\min}}\right) \|\mathbf{E} - \mathbf{E}_h\|_{\text{curl},\omega,\varepsilon,\chi,\mathcal{T}_{K,h}} + \omega \|\mathbf{H} - \mathbf{H}_h\|_{\mu,K} + \text{osc}_{\mathcal{T}_{F,h}}(\mathbf{E}^{\text{inc}}) $$

for all $K \in \mathcal{T}_h$.

*Proof.* We employ the same notation for $\llbracket \mathbf{E}_h^{\text{inc}} \rrbracket$ as in the previous proof, and we only detail the proof of (4.12b) since (4.12a) is established similarly, given the “symmetry” of the formulation with respect to the electric and magnetic fields. We have

$$ i\omega\mu\mathbf{H}_h + \nabla \times \mathbf{E}_h = -i\omega\mu(\mathbf{H} - \mathbf{H}_h) - \mu\chi\nabla \times (\mathbf{E} - \mathbf{E}_h). $$

As a result, it holds that

$$ \|i\omega\mu\mathbf{H}_h + \nabla \times \mathbf{E}_h\|_K \le \sqrt{\mu_{K,\max}}\omega \|\mathbf{H} - \mathbf{H}_h\|_{\mu,K} + \frac{\mu_{K,\max}}{\sqrt{\mu_{K,\min}}} \|\nabla \times (\mathbf{E} - \mathbf{E}_h)\|_{\chi,K} $$

and

$$ (4.13) \qquad \frac{1}{\sqrt{\mu_{\tilde{K},\max}}} \|i\omega\mu\mathbf{H}_h + \nabla \times \mathbf{E}_h\|_K \lesssim \omega \|\mathbf{H} - \mathbf{H}_h\|_{\mu,K} + \|\nabla \times (\mathbf{E} - \mathbf{E}_h)\|_{\chi,K}. $$

On the other hand, for a face $F \in \mathcal{F}_K$, we set $\mathbf{w}_F := \mathcal{E}(\omega \llbracket \mathbf{E}_h - \mathbf{E}_h^{\text{inc}} \rrbracket \times \mathbf{n}_F)$. With this notation, since $b_F = 0$ on $\partial \tilde{F}$ and $\mathbf{E} \in \mathcal{H}_0(\text{curl}, \Omega)$, we have

$$ \begin{aligned} & \omega^2 \|b_F^{1/2} \llbracket \mathbf{E}_h - \mathbf{E}_h^{\text{inc}} \rrbracket \times \mathbf{n}_F \|_F^2 \\ &= |\omega(\llbracket \mathbf{E}_h - \mathbf{E}^{\text{inc}} \rrbracket \times \mathbf{n}_F, \mathbf{w}_F)_F + \omega(\llbracket \mathbf{E}^{\text{inc}} - \mathbf{E}_h^{\text{inc}} \rrbracket \times \mathbf{n}_F, \mathbf{w}_F)_F| \\ &= |\omega(\mathbf{E} - \mathbf{E}_h, \nabla \times \mathbf{w}_F)_{\mathcal{T}_{F,h}} + (\nabla \times (\mathbf{E} - \mathbf{E}_h), \mathbf{w}_F)_{\mathcal{T}_{F,h}} + \omega(\llbracket \mathbf{E}^{\text{inc}} - \mathbf{E}_h^{\text{inc}} \rrbracket \times \mathbf{n}_F, \mathbf{w}_F)_F| \\ &\le \omega |(\mathbf{E} - \mathbf{E}_h, \nabla \times \mathbf{w}_F)_{\mathcal{T}_{F,h}} + (\nabla \times (\mathbf{E} - \mathbf{E}_h), \mathbf{w}_F)_{\mathcal{T}_{F,h}}| + \text{osc}_{\mathcal{T}_{F,h}}(\mathbf{E}^{\text{inc}}). \end{aligned} $$

Then, it follows from (2.7a), (2.7c) and (2.4) that

$$ \begin{aligned} & \sqrt{\varepsilon_{\tilde{K},\max}} \omega \sqrt{\frac{h_K}{p_K}} \|\llbracket \mathbf{E}_h - \mathbf{E}_h^{\text{inc}} \rrbracket \times \mathbf{n}_F\|_F \lesssim \sqrt{\varepsilon_{\tilde{K},\max}} \omega \sqrt{\frac{h_K}{p_K}} p_F \|b_F^{1/2} \llbracket \mathbf{E}_h - \mathbf{E}_h^{\text{inc}} \rrbracket \times \mathbf{n}_F\|_F \\ & (4.14) \qquad \qquad \qquad \qquad \qquad \lesssim p_K \left(1 + \frac{\omega h_K}{p_K c_{\tilde{K},\min}}\right) \|\mathbf{E} - \mathbf{E}_h\|_{\text{curl},\omega,\varepsilon,\chi,\mathcal{T}_{F,h}} + \text{osc}_{\mathcal{T}_{F,h}}(\mathbf{E}^{\text{inc}}). \end{aligned} $$

Finally, (4.12b) follows from the definition of $\eta_{\text{curl},\mu,K}$, (4.13) and (4.14). $\square$

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

Our key efficiency estimate is a direct consequence of Lemmas 4.9 and 4.10.

**Theorem 4.11 (Efficiency).** The estimate

$$ \eta_K \lesssim p_K^{3/2} \left(1 + \frac{\omega h_K}{p_K c_{\tilde{K},\min}}\right) \|\mathbf{E} - \mathbf{E}_h, \mathbf{H} - \mathbf{H}_h\|_{\text{curl},\omega,\mathcal{T}_{K,h}} + \text{osc}_{\mathcal{T}_{K,h}}(\mathbf{E}^{\text{inc}}). $$

holds true for all $K \in \mathcal{T}_h$.
## 5. NUMERICAL EXPERIMENTS

5.1. **Settings.** For the sake of simplicity, we introduce the frequency $\nu := \omega/2\pi$, and we will assume that $\varepsilon_0 = \mu_0 = 1$ in all the examples below. $\tilde{\xi}$ and $\tilde{\eta}$ respectively denote the relative error and estimators that have been scaled by the norm of the reference solution, i.e.

$$ \tilde{\xi} := \frac{\|(E - E_h, \mathbf{H} - \mathbf{H}_h)\|_{\mathrm{curl},\omega,\Omega}}{\|(E, \mathbf{H})\|_{\mathrm{curl},\omega,\Omega}}, \quad \tilde{\eta} := \frac{\eta}{\|(E, \mathbf{H})\|_{\mathrm{curl},\omega,\Omega}}. $$

5.1.1. *Structured meshes.* We often employ “structured” meshes in the following numerical examples. By this, we mean that the domain is (up to a translation) a square $\Omega := (-\ell, \ell)^2$ and that it is first partitioned into $N \times N$ identical squares. Each of these squares is then subdivided into four triangles by joining the barycenter with each face. The resulting mesh counts $4N^2$ triangular elements, with mesh size $h = 2\ell/N$.

5.1.2. *Unstructured meshes.* We also use “unstructured” meshes generated with a software package. Specifically, we employ mmg2D [21] with the options `-optim`, `-ar=0` and `-hmax=h` to generate a mesh of size $h$.

5.1.3. *Perfectly matched layers.* Most of our experiments employ perfectly matched layers (PML) to mimic an infinite propagation medium. This approach is standard [4, 36], and we proceed as follows. The computational domain will always be (up to a translation) a square $(-\ell, \ell)^2$, and the coefficients $\varepsilon$ and $\mu$ will take the values 1 and $\mathbf{I}$ in a neighborhood of the square's boundary. This original square is extended into a larger square $(-\ell + \ell_{PML}), \ell + \ell_{PML})^2$ for a fixed $\ell_{PML} > 0$, and the coefficients $\varepsilon$ and $\mu$ are artificially modified outside $(-\ell, \ell)^2$ as

$$ \varepsilon := 1/\nu_1\nu_2, \quad \mu := \begin{pmatrix} \nu_2/\nu_1 & 0 \\ 0 & \nu_1/\nu_2 \end{pmatrix}, $$

where $\nu_j(\mathbf{x}) = 1 + i\chi_{|\mathbf{x}_j|>\ell}$ and $\chi_{|\mathbf{x}_j|>\ell}$ is the characteristic function of the set $\{\mathbf{x} \in \mathbb{R}^2; |\mathbf{x}_j| > \ell\}$, for $j = 1$ or 2.

5.1.4. *h-adaptivity.* In several experiments, we consider h-adaptive iterative refinements, and our strategy combines Dörfler's marking [23] with the newest-vertex bisection [8]. Specifically, once the estimator $\{\eta_K\}_{K \in \mathcal{T}_h}$ has been computed, we order the elements $\{K_j\}_{j=1}^{|\mathcal{T}_h|}$ in such a way that $\eta_{K_{j+1}} \le \eta_{K_j}$, and we then select the smallest number $n$ of elements such that

$$ \sum_{j=1}^{n} \eta_{K_j}^2 \ge \theta \eta^2, $$

where $\theta := 0.05$. The elements $\{K_j\}_{j=1}^n$ are then refined using the newest-vertex bisection, starting from $K_1$ and finishing with $K_n$.

5.1.5. *hp-adaptivity.* We will also consider hp-adaptive refinements. In this case, the elements are still marked using Dörfler marking as above, and we employ an algorithm based on [35] to decide between $h$ and $p$ refinements. Specifically, if an element $K$ with diameter $h_{old}$, order $p_{old}$ and estimator $\eta_{old}$ has been refined into new elements $\{\kappa_j\}_{j=1}^n$ with diameters $\{h_{\kappa_j}\}_{j=1}^n$ order $\{p_{\kappa_j}\}_{j=1}^n$ and estimators $\{\eta_{\kappa_j}\}_{j=1}^n$, we define the “ideal” error reduction to be

$$ \mathrm{red} = \left( \frac{h_{\mathrm{new}}}{p_{\mathrm{new}}} \right)^{p_{\mathrm{new}}} \left( \frac{h_{\mathrm{old}}}{p_{\mathrm{old}}} \right)^{p_{\mathrm{old}}}, $$
where $h_{new} = \max_j h_{\kappa_j}$ and $p_{new} = \min_j p_{\kappa_j}$. Then, letting

$$ \eta_{new}^2 := \sum_{j=1}^{n} \eta_{\kappa_j}^2, $$

we perform a p-refinement if $\eta_{new} \le \text{red}\eta_{old}$ and an h-refinement otherwise. Once the new p-distribution has been obtained, it is “smoothed” to ensure that the polynomial degree of two neighboring only varies by one through an iterative increase in the degree of neighboring elements when required. For the first iteration (where we have no history to compute the ideal reduction), we only employ p-refinements.

### 5.2. Planewave in free space.
For our first example, we consider the computational domain $\Omega_0 := (-1, 1)^2$, that we surround with a PML layer of thickness $\ell_{PML} := 0.25$. The entire domain is thus $\Omega := (-1 + \ell_{PML}), 1 + \ell_{PML})^2$. We also set the TF region $\Omega_{TF} := (-\ell_{TF}, \ell_{TF})^2$ with $\ell_{TF} := 0.75$. $\varepsilon = 1$ and $\mu = I$ in $\Omega_0$ and are modified as explained in Section 5.1.3 in $\Omega \setminus \Omega_0$. We then set $J := 0$, and

$$ E_{\text{inc}} := e^{-i\omega d \cdot x}, \quad H_{\text{inc}} := (i\omega)^{-1} \text{curl} E_{\text{inc}} $$

where $d := (\cos \phi, \sin \phi)$ with $\phi := \pi/3$. The analytic solution is then simply given by $E := E_{\text{inc}}\chi_{\Omega_{\text{TF}}}$ and $H := H_{\text{inc}}\chi_{\Omega_{\text{TF}}}$, where $\chi_{\Omega_{\text{TF}}}$ is the characteristic function of $\Omega_{\text{TF}}$. This experiment employs structured meshes as defined in Section 5.1.1. Our goal is to illustrate the behavior of the estimator as the frequency increases.

Figures 1 and 2 respectively report the errors and effectivity indices for different frequencies $\omega$ and polynomial degrees $p$. The observed results are exactly in line with our theoretical prediction: the error is underestimated for coarse meshes, but this effect disappears asymptotically. The asymptotic regime is achieved faster for higher-order methods. Also, the underestimation is more pronounced for higher frequencies. We finally note that the effectivity indices approach one as the meshes are refined, which is expected since the solution is regular here.

### 5.3. Nearly resonant frequencies.
In this example, we consider the unit square $\Omega := (0, 1)^2$ surrounded by PEC boundary conditions (i.e., we do not use PMLs). The source term is $J := 1$, and no incident fields are injected. The (semi-)analytical solution is obtained via Fourier series as

$$ E(x) := (16i\omega) \sum_{n,m \text{ odd}} \frac{1}{nm\pi^2} \frac{1}{(n^2 + m^2)\pi^2 - \omega^2} \sin(n\pi x_1) \sin(m\pi x_2) $$

which we cut at $n, m \le 500$. The magnetic field is obtained by (analytically) differentiating $E$. Notice that the solution belongs to $H^3(\Omega)$, but not to $H^4(\Omega)$.

The first goal of this example is to highlight the behavior of the estimator when approaching a resonance frequency. We focus on the resonance frequency $\nu_r := \sqrt{2}/2$, and we consider a sequence of frequencies

$$ \nu_{\delta} := (1 + \delta) \frac{\sqrt{2}}{2} = \nu_{r} + \frac{\sqrt{2}}{2} \delta, $$

for decreasing values of $\delta \in \{4^{-r}\}_{r=1}^5$. In contrast to the previous example, we employ unstructured meshes here. Figure 3 presents the behavior of the error and the estimator for different values of $\delta$ and $p$ as the mesh is refined, whereas effectivity indices are given in Figure 4. The behavior is similar to the one observed in Figures 1 and 2 and conforms to our
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Planewave example: error and estimator</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Planewave example: efficiencies</figcaption></figure>
theoretical prediction. Indeed, the error is underestimated pre-asymptotically, and this effect is amplified when nearing $v_r$. The asymptotic range is achieved faster for higher polynomial degrees. In Figure 3, we observe the optimal convergence rates for uniform meshes, namely $N_{\text{dofs}}^{\min(p,2)/2}$, given the finite regularity of the solution. Notice also that the “suboptimal” convergence rates for $p = 4$ are only seen “late” in the convergence curves, which is in agreement with the regularity splitting results of [13]. We also observed in Figure 4 that the effectivity indices asymptotically approach one for $p = 1$ and $p = 2$, which is coherent with Remark 4.8 since $E \in H^3(\Omega)$ but $E \notin H^4(\Omega)$.

The next purpose of the present test is to analyze the ability of the estimator to drive h-adaptive refinements. Results for $p = 3$ and 4 are presented for the frequencies $\nu := 1.1$ and 5.1 (respectively close to the resonance frequencies $\nu_r := \sqrt{1^2 + 2^2/2} \sim 1.118$ and $\sqrt{6^2 + 8^2/2} = 5$). We start with a structured mesh for $p = 3$ and an unstructured mesh for $p = 4$. Figure 5 shows the behavior of the error and estimators where our previous comments on the pre-asymptotic underestimation apply. We further see that in contrast with the uniform meshes, the optimal convergence rates are observed despite the finite regularity of the solution, which is in particular due to the refinements close to the domain's corners that can be seen in Figure 7 when $\nu := 1.1$. Notice that for the higher frequency $\nu := 5.1$, we do not observe local refinements in Figure 8, which is again in agreement with [13] and the comment we made on Figure 3. Finally, Figure 6 presents the effectivity indices. The usual pre-asymptotic underestimation is observed for the higher frequency. Besides, the estimator seems asymptotic exact, which is in agreement with Remark 4.8 since local refinements are performed here. We finally note that Figures 7 and 8 show an excellent agreement between the elementwise error and estimator.

5.4. **Reflection of a Gaussian beam.** In this example, we consider the reflection of a Gaussian beam by two heterogeneous prisms modeled by heterogeneous coefficients. The domain of interest $\Omega_0 := (-\ell, \ell)^2$, where $\ell := 1$, is surrounded by a PML of thickness $\ell_{\text{PML}} := 0.1$, so that $\Omega := (-\ell + \ell_{\text{PML}}, \ell + \ell_{\text{PML}})^2$. We set $J := 0$, and we define the TF region $\Omega_{\text{TF}} := (-\ell_{\text{TF}}, \ell_{\text{TF}})^2$ with $\ell_{\text{TF}} := 0.9$. The incident field is the Gaussian beam

$$E_{\text{inc}}(\boldsymbol{x}) := \frac{1}{\sqrt{4\alpha\beta(y_1)}} \exp\left(-\frac{y_2^2}{4\beta(y_1)}\right) \exp(-i\omega y_1), \quad \beta(y_1) := \frac{1}{4\alpha} + i\frac{y_1}{2\omega},$$

with $\boldsymbol{y} = \boldsymbol{x} - \boldsymbol{x}_0$, which is the solution to paraxial Helmholtz equation

$$-\omega^2 u - i\omega \frac{\partial u}{\partial y_1} - \frac{\partial^2 u}{\partial y_2^2} = 0, \quad u(0, y_2) = g(y_2)$$

with initial data $g(y_2) := \exp(-\alpha y_2)$, for $\boldsymbol{x}_0 := (-1.0, 0.3)$ and $\alpha := w^{-2}$, $w = 0.1$, and we set $H_{\text{inc}} := (i\omega)^{-1}\text{curl}E_{\text{inc}}$. Notice that, strictly speaking, these incident fields do not enter our assumption in (2.1) because $E_{\text{inc}}$ solves the paraxial Helmholtz equation, which is only an approximation of the “true” Helmholtz equation. Nevertheless, the paraxial approximation is fairly good for Gaussian beams, so the resulting discrepancy is not seen numerically for the accuracy level we target. The prisms are the two triangles

$$P_1 := \{(0.3, 0.0), (0.6, 0.6), (0.0, 0.6)\}, \quad P_2 := \{(-0.3, 0.0), (0.0, -0.6), (-0.6, -0.6)\},$$

and the electric permittivity is defined by $\varepsilon := 1/10$ in $P_1 \cup P_2$ and $\varepsilon := 1$ outside, and we apply the usual modifications in the PMLs. We also select the frequency $\nu := 50$. The whole setting is sketched on Figure 9, where the initial mesh is also represented.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Nearly resonant example: error and estimator</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Nearly resonant example: efficiencies with upwind fluxes</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Nearly resonant experiment: errors in *h*-adaptive refinements</figcaption></figure>

We start the adaptive algorithm with the mesh represented in Figure 9, and we employ the fixed polynomial degree $p := 6$. We run the adaptive loop for 30 iterations, and the final discrete solution is represented in Figure 10. In Figure 11, we plot the value of the estimator against the number of degrees of freedom throughout the adaptive loop (we do not represent the true error as it is not available here). We observe an initial stagnation, which is coherent with other experiments of adaptivity for high-frequency waves [12]. This initial stagnation is actually expected since we start with a largely unresolved mesh. We then see that the optimal convergence rate is asymptotically reached.

Figure 12 depicts the solution, the mesh size, and the estimator at various iterations of the adaptive loop. Interestingly, we see that the refinements essentially follow the wavefront until the beam goes through the whole domain. More uniform refinements then occur to capture the diffracted rays. We also observe that the mesh is refined on some corners of the prisms, which agrees with the expected presence of singularities at these points.
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Nearly resonant experiment: effectivity indices in *h*-adaptive refinements</figcaption></figure>

**5.5. Scattering by an aircraft.** Our final example is the scattering of planewave by an aircraft. The incident field is given by $E_{\text{inc}} := e^{i\omega d \cdot x}$ with $d = (\cos \phi, \sin \phi)$ and $\phi := 13\pi/12$ whereas $H_{\text{inc}} := (i\omega)^{-1}\text{curl}E_{\text{inc}}$. The domain of interest is the square $\Omega_0 := (-\ell, \ell)^2 \setminus \bar{A}$, where $A \subset\subset \Omega_0$ is an open Lipschitz polygon with 76 vertices representing an aircraft (see Figure 13) and $\ell := 40$, the TF region is $\Omega_{\text{TF}} := (-\ell_{\text{TF}}, \ell_{\text{TF}})^2 \setminus A$ with $\ell_{\text{TF}}$ and we surround $\Omega_0$ with PMLs of length $\ell_{\text{PML}} := 5$. We impose the PEC condition on the boundary of $\Omega := -(-\ell + \ell_{\text{PML}}), \ell + \ell_{\text{PML}})^2 \setminus \bar{A}$. Figure 13 depicts the whole setting.

We employ the *hp*-adaptive algorithm presented in Section 5.1.5 for 15 iterations, starting with the mesh shown in Figure 13 with a uniform polynomial degree distribution $p = 1$. We consider two frequencies: $\nu := 0.1$ and $\nu := 0.5$. The solutions are represented in Figure 14, whereas Figure 15 shows the convergence history of the *hp*-adaptive loop in both cases. We observe an exponential convergence rate for the two frequencies (following [35] and [45, Theorem 4.63], we plot the estimator $\tilde{\eta}$ in log-scale against $N_{\text{dofs}}^{1/3}$ in linear-scale), indicating that the proposed estimator is suited to drive *hp*-adaptive algorithms.
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Nearly resonant experiment: actual (left) and estimated (error) at iteration # 10 in the *h*-adaptive example with $\nu = 1.1$</figcaption></figure>

6. CONCLUSION

We propose a new residual-based a posteriori error estimator for discontinuous Galerkin discretizations (DG) of time-harmonic Maxwell's equations in first-order form. Our estimator covers a range of numerical DG fluxes, including the so-called “central” and “upwind” fluxes. We rigorously establish that the estimator is reliable and efficient, with constants that are independent of the frequency for sufficiently refined meshes. Besides, we show that our estimator is asymptotically constant-free for a smooth solution.

We also present a set of numerical examples that highlights our key findings. We find an excellent agreement between these examples and the behavior predicted by the theory for our estimator. In addition, we employ the estimator to drive *h*- and *hp*-adaptive refinement processes. These examples suggest that the proposed estimator is capable of driving such refinement processes since, in all cases, we observe optimal convergence rates.
<figure><img src="image_10.png" /><figcaption>FIGURE 8. Nearly resonant experiment: actual (left) and estimated (error) at iteration # 10 in the *h*-adaptive example with $\nu = 5.1$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 9. Setting and initial mesh of the Gaussian beam experiment</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIGURE 10. Real part of the solution in the Gaussian beam experiment</figcaption></figure>

## FUNDING

PV was supported by the Chilean National Research and Development Agency (ANID) through the grant ANID FONDECYT No. 3220858.

## REFERENCES

1. R. Adams and J. Fournier, Sobolev spaces, Academic Press, 2003.

2. D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini, *Unified analysis of discontinuous Galerkin methods for elliptic problems*, SIAM J. Numer. Anal. **39** (2002), 1749–1779.

3. R. Beck, R. Hiptmair, R. H. W. Hoppe, and B. Wohlmuth, *Residual based a posteriori error estimators for eddy current computation*, ESAIM Math. Model. Numer. Anal. **34** (2000), 159–182.

4. J. P. Bérenger, *A perfectly matched layer for the absorption of electromagnetics waves*, J. Comput. Phys. **114** (1994), 185–200.
<figure><img src="image_13.png" /><figcaption>FIGURE 11. Convergence in the Gaussian beam experiment</figcaption></figure>

5. __________, *Three-dimensional perfectly matched layer for the absorption of electromagnetic waves*, J. Comput. Phys. **127** (1996), 363–379.

6. M. Bernkopf, S. Sauter, C. Torres, and A. Veit, *Solvability of discrete Helmholtz equations*, arXiv:2105.02273v2, 2022.

7. A. Bespalov, A. Haberl, and D. Praetorius, *Adaptive FEM with coarse initial mesh guarantees optimal convergence rates for compactly perturbed elliptic problems*, Comput. Methods Appl. Mech. Engrg. **317** (2017), 318–340.

8. P. Binev, W. Dahmen, and R. De Vore, *Adaptive finite element methods with convergence rates*, Numer. Math. **97** (2004), 219–268.

9. A. Bonito and R. H. Nochetto, *Quasi-optimal convergence rate of an adaptive discontinuous Galerkin method*, SIAM J. Numer. Anal. **48** (2010), 734–771.

10. A. Buffa and P. Ciarlet Jr., *On traces for functional spaces related to Maxwell's equations part I: an integration by parts formula in Lipschitz polyhedra*, Math. Meth. Appl. Sci. **24** (2001), 9–30.

11. T. Chaumont-Frelet, *Mixed finite element discretization of acoustic Helmholtz problems with high wavenumbers*, Calcolo **56** (2019).

12. T. Chaumont-Frelet, A. Ern, and M. Vohralík, *On the derivation of guaranteed and p-robust a posteriori error estimates for the Helmholtz equation*, Numer. Math. **148** (2021), 525–573.
<figure><img src="image_14.png" /><figcaption>FIGURE 12. Solution mesh size and estimator in the Gaussian beam experiment</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIGURE 13. Setting and initial mesh of the aircraft experiment</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 14. Real part of the solutions in the aircraft experiment</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 15. Convergence in the aircraft experiment</figcaption></figure>
13. T. Chaumont-Frelet and S. Nicaise, *High-frequency behaviour of corner singularities in Helmholtz problems*, ESAIM Math. Model. Numer. Anal. **5** (2018), 1803–1845.

14. —, *Wavenumber explicit convergence analysis for finite element discretizations of general wave propagation problems*, IMA J. Numer. Anal. **40** (2020), 1503–1543.

15. T. Chaumont-Frelet and P. Vega, *Frequency-explicit a posteriori error estimates for finite element discretizations of Maxwell's equations*, SIAM J. Numer. Anal. **60** (2022), 1774–1798.

16. —, *Frequency-explicit approximability estimates for finite element discretizations of Maxwell's equations*, Calcolo **59** (2022).

17. J. Chen, Y. Xu, and J. Zou, *An adaptive edge element method and its convergence for a saddle-point problem from magnetostatics*, Numer. Methods Partial Differ. Equ. **28** (2011), 1643–1666.

18. P. G. Ciarlet, *The finite element method for elliptic problems*, SIAM, 2002.

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

20. A. Demlow and A. N. Hirani, *A posteriori error estimates for finite element exterior calculus: the de Rham complex*, Found. Comput. Math. **14** (2014), 1337–1371.

21. C. Dobrzynski, *MMG3D: User guide*, Tech. Report 422, Inria, 2012.

22. R. C. Dorf, *Electronics, power electronics, optoelectronics, microwaves, electromagnetics and radar*, Taylor & Francis, 2006.

23. W. Dörfler, *A convergent adaptive algorithm for Poisson's equation*, SIAM J. Numer. Anal. **33** (1996), 1106–1124.

24. W. Dörfler and S. Sauter, *A posteriori error estimation for highly indefinite Helmholtz problems*, Comput. Meth. Appl. Math. **13** (2013), 333–347.

25. X. Feng, P. Lu, and X. Xu, *A hybridizable discontinuous Galerkin method for the time-harmonic Maxwell equations with high wave number*, Comput. Methods Appl. Math. **16** (2016), 429–445.

26. V. Girault and P. A. Raviart, *Finite element methods for Navier-Stokes equations: theory and algorithms*, Springer-Verlag, 1986.

27. D. J. Griffiths, *Introduction to Electrodynamics*, Prentice Hall, 1999.

28. J. S. Hesthaven and T. Warburton, *Nodal high-order methods on unstructured grids. Part I. Time-domain solution of Maxwell's equations*, J. Comput. Phys. **181** (2002), 1266–1288.

29. R. Hiptmair and C. Pechstein, *Discrete regular decompositions of tetrahedral discrete 1-forms*, in *Maxwells Equations*, De Gruyter, Berlin, 2019, pp. 199–258.

30. L. Li, S. Lanteri, and R. Perrussel, *A hybridizable discontinuous Galerkin method combined to a Schwarz algorithm for the solution of 3d time-harmonic Maxwell's equations*, J. Comput. Phys. **256** (2014), 563–581.

31. J. M. Melenk, *hp-interpolation of nonsmooth functions and an application to hp-a posteriori error estimation*, SIAM J. Numer. Anal. **43** (2005), 127–155.

32. J. M. Melenk and S. Sauter, *Wavenumber explicit convergence analysis for Galerkin discretizations of the Helmholtz equation*, SIAM J. Numer. Anal. **49** (2011), 1210–1243.

33. J. M. Melenk and S. A. Sauter, *Wavenumber-explicit hp-FEM analysis for Maxwell's equations with transparent boundary conditions*, Found. Comput. Math. **21** (2021), 125–241.

34. —, *Wavenumber-explicit hp-FEM analysis for Maxwell's equations with impedance boundary conditions*, preprint arXiv:2201.02602, 2022.

35. J. M. Melenk and B.I. Wohlmuth, *On residual-based a posteriori estimation in hp-FEM*, Adv. Comput. Math. **15** (2001), 311–331.

36. P. Monk, *Finite element methods for Maxwell's equations*, Oxford University Press, New York, 2003.

37. J. C. Nédélec, *A new family of mixed finite elements in $\mathbb{R}^3$*, Numer. Math. **50** (1986), 57–81.

38. N. C. Nguyen, J. Peraire, and B. Cockburn, *Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations*, J. Comput. Phy. **230** (2011), 7151–7175.

39. S. Nicaise and E. Creusé, *A posteriori error estimation for the heterogeneous Maxwell equations on isotropic and anisotropic meshes*, Calcolo **40** (2003), 249–271.

40. S. Nicaise and J. Tomezyk, *Convergence analysis of a hp-finite element approximation of the time-harmonic Maxwell equations with impedance boundary conditions in domains with an analytic boundary*, Numer. Methods Partial Differ. Equ. **36** (2020), 1868–1903.

41. I. Perugia and D. Schötzau, *The hp-local discontinuous Galerkin method for low-frequency time-harmonic Maxwell equations*, Math. Comp. **72** (2003), 1179–1214.
42. S. Sauter and J. Zech, *A posteriori error estimation of hp-dg finite element methods for highly indefinite Helmholtz problems*, SIAM J. Numer. Anal. **53** (2015), 2414–2440.

43. S.A. Sauter and C. Schwab, *Boundary element methods*, Springer, 2010.

44. J. Schöberl, *A posteriori error estimates for Maxwell equations*, Math. Comp. **77** (2008), 633–649.

45. C. Schwab, *p- and hp-finite element methods*, Oxford Univ. Press, 1998.

46. A. Taflove and S.C. Hagness, *Computational electrodynamics the finite-difference time-domain method*, Artech house, 2005.

47. R. Verfürth, *A posteriori error estimation and adaptive mesh-refinement techniques*, J. Comput. Appl. Math. **50** (1994), 67–83.

48. J. Viquerat, *Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method*, Ph.D. thesis, Université Nice Sophia-Antipolis and Inria project-team Nachos, 2015.