# A hybridizable discontinuous Galerkin method with transmission variables for time-harmonic electromagnetic problems

Ari Rappaport, Théophile Chaumont-Frelet, Axel Modave

► To cite this version:

Ari Rappaport, Théophile Chaumont-Frelet, Axel Modave. A hybridizable discontinuous Galerkin method with transmission variables for time-harmonic electromagnetic problems. 2025. hal-05016261v2

HAL Id: hal-05016261

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

Preprint submitted on 9 May 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A HYBRIDIZABLE DISCONTINUOUS GALERKIN METHOD WITH TRANSMISSION VARIABLES FOR TIME-HARMONIC ELECTROMAGNETIC PROBLEMS*

A. E. RAPPAPORT†, T. CHAUMONT-FRELET‡, AND A. MODAVE†

**Abstract.** The CHDG method is a hybridizable discontinuous Galerkin (HDG) finite element method suitable for the iterative solution of time-harmonic wave propagation problems. Hybrid unknowns corresponding to transmission variables are introduced at the element interfaces and the physical unknowns inside the elements are eliminated, resulting in a hybridized system with favorable properties for fast iterative solution. In this paper, we extend the CHDG method, initially studied for the Helmholtz equation, to the time-harmonic Maxwell equations. We prove that the local problems stemming from hybridization are well-posed and that the fixed-point iteration naturally associated to the hybridized system is contractive. We propose a 3D implementation with a discrete scheme based on nodal basis functions. The resulting solver and different iterative strategies are studied with several numerical examples using a high-performance parallel C++ code.

**Key words.** Maxwell equations, Discontinuous finite elements, High-order methods, Hybridization, Iterative solvers

**MSC codes.** 65F08, 65N22, 65N30

**1. Introduction.** The numerical solution of the time-harmonic Maxwell equations is fundamental for simulating electromagnetic wave propagation, with applications ranging from antenna design to photonics and radar systems. In this context, finite element methods are well established as a robust and flexible solution framework, offering advantages such as high-order accuracy and geometric flexibility. However, fine meshes may be required due to pollution effects. For this class of problems, the finite element methods lead to large, sparse, complex-valued linear systems, which are difficult to solve with direct and iterative solution methods. On the one hand, direct approaches are expensive and difficult to parallelize. On the other hand, iterative strategies may require many iterations due to intrinsic properties of the problem.

Discontinuous Galerkin (DG) methods have been developed and analyzed for both time-domain and time-harmonic electromagnetic problems, see e.g. [13, 17, 24] and [27, 5, 26, 2], respectively. From a practical point of view, the cell local nature of the degrees of freedom (DOFs) in DG methods leads to access patterns that are favorable to high performance computing and in particular, to GPUs [17, 13]. However, the number of unknowns in the resulting linear system is larger than for a corresponding continuous Galerkin counterpart, and solving these systems for time-harmonic application problems is still prohibitively expensive.

Several specialized methods have emerged to address these challenges and have been applied specifically to the time-harmonic Maxwell equations, namely hybridizable DG (HDG) methods [23, 18, 3] and Trefftz approaches [15, 14, 9, 25]. The HDG methods consist in introducing hybrid unknowns at the element interfaces and eliminating the physical unknowns. In the Trefftz approach, basis functions tailored for the considered problem are used. In both cases, the number of unknowns are reduced (sometimes drastically). Recent studies have also focused on advanced preconditioning techniques and domain decomposition methods to improve the scalability of DG-based solvers, see, e.g., [7, 6, 8, 19, 11] and references therein. This is still an active area of research and the interplay between the discretization and the preconditioning strategies is very much intertwined.

In the recent work [22], a HDG variant, termed CHDG, is studied in the context of Helmholtz problems. The key idea is to reformulate the local DG flux in terms of the so-called transmission

<footnote>*Distributed under Creative Commons CC BY 4.0 license.</footnote>

<footnote>**Funding:** This work was funded in part by the ANR JCJC project WavesDG (research grant ANR-21CE46-0010) and by the Agence de l'Innovation de Défense [AID] through Centre Interdisciplinaire d'Etudes pour la Défense et la Sécurité [CIEDS] (project 2022 ElectroMath). The authors thank the IDCS research facility of Ecole Polytechnique for providing access to the Cholesky cluster.</footnote>

<footnote>†POEMS, CNRS, Inria, ENSTA, Institut Polytechnique de Paris, 91120 Palaiseau, France (ari.rappaport@ensta.fr, axel.modave@ensta.fr).</footnote>

<footnote>‡Inria Univ. Lille and Laboratoire Paul Painlevé, 59655 Villeneuve-d'Ascq, France (theophile.chaumont@inria.fr).</footnote>
variables (incoming and outgoing) and then to introduce skeleton unknowns based on these variables. Then, an elimination process similar to the one in HDG allows to solve only for the incoming variable. The method naturally results in a system of the form $(I - \Pi S)g = b$ that can be solved with a fixed-point iteration, similar to the ultra weak variational formulation (UWVF) [4]. Numerical results indicate that standard Krylov solvers such as CGNR or GMRES require fewer iterations with CHDG than with other discretization methods.

In this paper, we extend CHDG to the time-harmonic Maxwell equations with constant coefficients, and we propose a 3D implementation based on nodal basis functions. Appropriate transmission variables are defined in this setting, and the good properties obtained in [22] are recovered. In particular, the element-local problems associated with the hybridization are well-posed, and the fixed-point iteration is a contraction. We give implementation details for the resulting scheme with a strategy based on nodal basis functions, which has advantages for the practical parallel implementation. Specific preconditioned CGNR and GMRES iterations of the hybridized system are discussed, which in practice lead to faster convergence.

The rest of this paper is organized as follows. In §2 we introduce the mathematical model, the standard DG method, the transmission variables, and the hybridization technique. In §3 the local and hybridized problems are analyzed. In §4 we present the discrete systems with implementation details, and we discuss linear solver strategies for the hybridized problem. In §5 we present numerical results on a suite of benchmarks that both validate code correctness and compare the solver strategies on reference problems. All experiments are conducted in 3D with a parallel C++ code. Finally, we present conclusion and future perspectives in §6.

## 2. Discontinuous Galerkin method and hybridization.
Let $\Omega \subset \mathbb{R}^3$ be a domain with a polyhedral boundary partitioned into three mutually disjoint polytopal Lipschitz subsets, i.e. $\partial\Omega = \Gamma_E \cup \Gamma_H \cup \Gamma_I$. In this paper, we consider the time-harmonic Maxwell equations completed with standard boundary conditions

$$ (2.1) \qquad \left\{ \begin{array}{ll} \iota\kappa e - \nabla \times h = 0 & \text{in } \Omega, \\ \iota\kappa h + \nabla \times e = 0 & \text{in } \Omega, \\ n \times e = s_E & \text{on } \Gamma_E, \\ n \times h = s_H & \text{on } \Gamma_H, \\ -n \times (n \times e) + n \times h = s_I & \text{on } \Gamma_I, \end{array} \right. $$

where the unknowns $e: \Omega \to \mathbb{C}^3$ and $h: \Omega \to \mathbb{C}^3$ represent the electric and magnetic field respectively, $\iota$ is the imaginary unit, $\kappa > 0$ represents the (constant) wave number, $n$ is the outward unit normal vector with respect to $\partial\Omega$, and $s_E$, $s_E$ and $s_E$ are given surface sources. Volumetric sources are easily handled as well.

### 2.1. Discretization and notation.
We consider a conforming mesh $\mathcal{T}_h$ of the domain $\Omega$ consisting of simplicial elements $K$. The collection of element boundaries is denoted by $\partial\mathcal{T}_h := \{\partial K | K \in \mathcal{T}_h\}$, the collection of faces of the mesh is denoted by $\mathcal{F}_h$, and the collection of faces of an element $K$ is denoted by $\mathcal{F}_K$.

The numerical physical fields, denoted by $e_h$ and $h_h$, belong to the broken space

$$ V_h := \prod_{K \in \mathcal{T}_h} \mathcal{P}_p(K), $$

where $p \ge 0$ is a given polynomial degree and $\mathcal{P}_p(\cdot)$ is the space of complex vector-valued polynomials of degree smaller than or equal to $p$. For a face $F$ of an element $K$, we shall use the tangential trace and tangential component trace, which, for a vector-valued polynomial $v_K \in \mathcal{P}_p(K)$, are defined as

$$ (2.2) \qquad n_{K,F} \times v_K \quad \text{and} \quad \pi_{K,F}^t(v_K) := -n_{K,F} \times (n_{K,F} \times v_K), $$

where $n_{K,F}$ denotes the outward unit normal to $F$. These traces verify the identity $n_{K,F} \times \pi_{K,F}^t(v_K) = n_{K,F} \times v_K$, and they belong to the tangential trace space

$$ (2.3) \qquad \mathcal{P}_p^t(F) := \{v \in \mathcal{P}_p(F) | n_F \cdot v = 0\}. $$
Let us note that the traces (2.2) can be applied to any tangential vector-valued polynomial $v_F^t \in \mathcal{P}_p^t(F)$ without ambiguities. The orientation of the normal $n_F$ depends on the context. We introduce the following sesquilinear forms for complex vector-valued functions

$$ 
\begin{align*}
(v, w)_K &:= \int_K v \cdot \bar{w} \, dx, & (v, w)_{\mathcal{T}_h} &:= \sum_{K \in \mathcal{T}_h} (v, w)_K, \\
\langle v, w \rangle_F &:= \int_F v \cdot \bar{w} \, d\sigma(x), & \langle v, w \rangle_{\mathcal{F}_h} &:= \sum_{F \in \mathcal{F}_h} \langle v, w \rangle_F, \\
\langle v, w \rangle_{\partial K} &:= \sum_{F \in \mathcal{F}_K} \langle v, w \rangle_F, & \langle v, w \rangle_{\partial \mathcal{T}_h} &:= \sum_{K \in \mathcal{T}_h} \langle v, w \rangle_{\partial K}.
\end{align*}
$$

By convention, the quantities used in the surface integral $\langle \cdot, \cdot \rangle_{\partial K}$ are associated to the faces of $K$, e.g. $n_{K,F}$ and $\pi_{K,F}^t(v_K)$ with $F \in \mathcal{F}_K$. The subscripts $K$ and $F$ are dropped when there is no danger of ambiguity from context. For $v, w \in \mathcal{P}_p(K)$, we have the integration by part formula

$$ (\nabla \times w, v)_K - (w, \nabla \times v)_K = \langle n \times w, \pi^t(v) \rangle_{\partial K} $$

and the following identities

$$ 
\begin{align*}
\langle \pi_{K,F}^t(v), n_{K,F} \times w \rangle_F &= -\langle n_{K,F} \times v, \pi_{K,F}^t(w) \rangle_F, \\
\langle n_{K,F} \times v, n_{K,F} \times w \rangle_F &= \langle \pi_{K,F}^t(v), \pi_{K,F}^t(w) \rangle_F.
\end{align*}
$$

**2.2. Standard DG formulation with upwind fluxes.** The standard DG formulation for the time-harmonic Maxwell equations is given by

**PROBLEM 2.1.** Find $(e_h, h_h) \in V_h \times V_h$ such that, for all $(v_h, w_h) \in V_h \times V_h$,

$$ 
\left\{ 
\begin{array}{l}
\iota_K(e_h, v_h)_{\mathcal{T}_h} - (h_h, \nabla \times v_h)_{\mathcal{T}_h} - \langle n \times \hat{h}^t(e_h, h_h), \pi^t(v_h) \rangle_{\partial \mathcal{T}_h} = 0, \\
\iota_K(h_h, w_h)_{\mathcal{T}_h} + (e_h, \nabla \times w_h)_{\mathcal{T}_h} + \langle n \times \hat{e}^t(e_h, h_h), \pi^t(w_h) \rangle_{\partial \mathcal{T}_h} = 0,
\end{array}
\right.
$$

with the numerical fluxes $n \times \hat{h}^t(e_h, h_h)$ and $n \times \hat{e}^t(e_h, h_h)$ defined below.

We choose upwind fluxes, see e.g. [13, 24]. For each interior face $F \not\subset \partial\Omega$ of each element $K$, they are given by

$$ (2.4a) \qquad \left\{ 
\begin{array}{l}
n_{K,F} \times \hat{e}_F^t := n_{K,F} \times (e_K + e_{K'})/2 + \pi_{K,F}^t(h_K - h_{K'})/2, \\
n_{K,F} \times \hat{h}_F^t := n_{K,F} \times (h_K + h_{K'})/2 - \pi_{K,F}^t(e_K - e_{K'})/2,
\end{array}
\right. $$

where $K'$ is the neighboring element of $K$ sharing the face $F$. For a boundary face $F$, they are defined as

$$ (2.4b) \qquad \left\{ 
\begin{array}{l}
n_{K,F} \times \hat{e}_F^t := s_E \\
n_{K,F} \times \hat{h}_F^t := n_{K,F} \times h_K - (\pi_{K,F}^t(e_K) + n_{K,F} \times s_E)
\end{array}
\right. \qquad \text{if } F \subset \Gamma_E, $$

$$ (2.4c) \qquad \left\{ 
\begin{array}{l}
n_{K,F} \times \hat{e}_F^t := n_{K,F} \times e_K + (\pi_{K,F}^t(h_K) + n_{K,F} \times s_H) \\
n_{K,F} \times \hat{h}_F^t := s_H
\end{array}
\right. \qquad \text{if } F \subset \Gamma_H, $$

$$ (2.4d) \qquad \left\{ 
\begin{array}{l}
n_{K,F} \times \hat{e}_F^t := n_{K,F} \times (\pi_{K,F}^t(e_K) - n_{K,F} \times h_K + s_I)/2 \\
n_{K,F} \times \hat{h}_F^t := -(\pi_{K,F}^t(e_K) - n_{K,F} \times h_K - s_I)/2
\end{array}
\right. \qquad \text{if } F \subset \Gamma_I. $$

These numerical fluxes satisfy

$$ 
\begin{align*}
\pi_{K,F}^t(\hat{e}_F^t) - n_{K,F} \times \hat{h}_F^t &= \pi_{K,F}^t(e_K) - n_{K,F} \times h_K, \\
n_{K,F} \times \hat{e}_F^t + \pi_{K,F}^t(\hat{h}_F^t) &= n_{K,F} \times e_K + \pi_{K,F}^t(h_K),
\end{align*}
$$

regardless of whether $F$ is an interior or boundary face.
<figure><img src="image_2.png" /><figcaption>Fig. 2.1: Notation for the outgoing and incoming transmission variables at the face F shared by an element K and a neighboring element K'.</figcaption></figure>

**2.3. Transmission variables.** The definition of the upwind fluxes is related to the characteristic analysis of the system, which exhibits quantities transported through the faces of the elements.

We consider the time-dependent Maxwell equations

$$ egin{cases} \frac{\partial e}{\partial t} - \nabla \times \mathbf{h} = \mathbf{0}, \\ \frac{\partial \mathbf{h}}{\partial t} + \nabla \times \mathbf{e} = \mathbf{0}. \end{cases} $$

Assuming only the variation with respect to a constant direction $\mathbf{n}$ leads to system

$$ egin{cases} \frac{\partial e}{\partial t} - \mathbf{n} \times \frac{\partial \mathbf{h}}{\partial n} = \mathbf{0}, \\ \frac{\partial \mathbf{h}}{\partial t} + \mathbf{n} \times \frac{\partial e}{\partial n} = \mathbf{0}. \end{cases} $$

Applying $-\mathbf{n} \times (\mathbf{n} \times \mathbf{e})$ and $-(\mathbf{n} \times \mathbf{e})$ to the first and second equations, respectively, and combining the resulting equations yields the following decoupled transport equations

$$ \frac{\partial \mathbf{q}^{\pm}}{\partial t} \pm \frac{\partial \mathbf{q}^{\pm}}{\partial n} = \mathbf{0}, $$

where the transported quantities are $\mathbf{q}^{\pm} := -\mathbf{n} \times (\mathbf{n} \times \mathbf{e}) \mp \mathbf{n} \times \mathbf{h}$, which are frequently called the *characteristic variables*.

The upwind fluxes can be rewritten by introducing transmission variables related to the quantities transported through the faces of the elements. At each face F of an element K, we define the *outgoing transmission variable* as

$$ (2.5a) \qquad \mathbf{g}_{K,F}^{\oplus} := \pi_{K,F}^{t}(e_K) - \mathbf{n}_{K,F} \times \mathbf{h}_K, $$

and the *incoming transmission variable* as

$$ (2.5b) \qquad \mathbf{g}_{K,F}^{\ominus} := \begin{cases} \mathbf{g}_{K',F}^{\oplus}, & \text{if } F \not\subset \partial\Omega \text{ is shared by } K \text{ and } K', \\ -\mathbf{g}_{K,F}^{\oplus} - 2\mathbf{n}_{K,F} \times \mathbf{s}_E, & \text{if } F \subset \Gamma_E, \\ \mathbf{g}_{K,F}^{\oplus} + 2\mathbf{s}_H, & \text{if } F \subset \Gamma_H, \\ \mathbf{s}_I, & \text{if } F \subset \Gamma_I. \end{cases} $$

This notation is illustrated in Figure 2.1. By using the transmission variables (2.5a)-(2.5b), the numerical fluxes (2.4a)-(2.4d) can be rewritten as

$$ (2.6) \qquad \begin{cases} \mathbf{n}_{K,F} \times \hat{\mathbf{e}}_F^t = \mathbf{n}_{K,F} \times (\mathbf{g}_{K,F}^{\oplus} + \mathbf{g}_{K,F}^{\ominus})/2, \\ \mathbf{n}_{K,F} \times \hat{\mathbf{h}}_F^t = -(\mathbf{g}_{K,F}^{\oplus} - \mathbf{g}_{K,F}^{\ominus})/2, \end{cases} $$

for any face $F \in \mathcal{F}_K$.

For an interior face ($F \not\subset \partial\Omega$), the outgoing variable of one side corresponds to the incoming one of the other side, i.e. $\mathbf{g}_{K,F}^{\oplus} = \mathbf{g}_{K',F}^{\ominus}$ and $\mathbf{g}_{K,F}^{\ominus} = \mathbf{g}_{K',F}^{\oplus}$. For a boundary face, the boundary condition is prescribed by a specific definition of the incoming variable. For example, for a condition on the tangential trace of a field ($F \subset \Gamma_E$ or $\Gamma_H$), the incoming variable is such that the numerical flux corresponding to the field is equal to the surface source.
**2.4. Hybridization.** In the standard hybridization approach, a new variable corresponding to the numerical trace associated to one of the field is introduced, see e.g. [18, 19] and references therein. This variable, denoted by $\zeta_h$, belongs to the space of the tangential vector-valued polynomials defined on the faces of the mesh,

$$ \hat{V}_h^t := \prod_{F \in \mathcal{F}_h} \mathcal{P}_p^t(F). $$

It can be defined such that $n_{K,F} \times \zeta_F := n_{K,F} \times \hat{h}_F^t$ for each face $F$ of each element $K$. The unknowns associated with the physical fields can then be eliminated from the system, resulting in a hybridized system defined on the skeleton.

As in [22], we propose an alternative approach where the new variable, denoted by $g_h^{\ominus}$, corresponds to the incoming transmission variable. It belongs to the space

$$ (2.7) \qquad G_h^t := \prod_{K \in \mathcal{T}_h} \prod_{F \in \mathcal{F}_K} \mathcal{P}_p^t(F). $$

By contrast with the standard approach, there are two tangential vector values for each interior face instead of one, i.e. $g_{K,F}^{\ominus}$ and $g_{K',F}^{\ominus}$ for a face $F$ shared by neighboring elements $K$ and $K'$. There is still one value for each boundary face.

The CHDG formulation is obtained by introducing the new variable $g_h^{\ominus}$ in Problem 2.1, and by using (2.6). To write the formulation, we define the global exchange operator $\Pi: G_h^t \to G_h^t$ given, for each face $F$ of each element $K$, by

$$ (2.8) \qquad \Pi(g^{\oplus})|_{K,F} := \begin{cases} g_{K',F}^{\oplus} & \text{if } F \not\subset \partial\Omega \text{ is shared by } K \text{ and } K', \\ -g_{K,F}^{\oplus} & \text{if } F \subset \Gamma_E, \\ g_{K,F}^{\oplus} & \text{if } F \subset \Gamma_H, \\ 0 & \text{if } F \subset \Gamma_I, \end{cases} $$

for any $g^{\oplus} \in G_h^t$. We also define the global right-hand side $b$ given, for each face $F$ of each element $K$, by

$$ b|_{K,F} := \begin{cases} 0 & \text{if } F \not\subset \partial\Omega, \\ -2n_{K,F} \times s_E & \text{if } F \subset \Gamma_E, \\ 2s_H & \text{if } F \subset \Gamma_H, \\ s_I & \text{if } F \subset \Gamma_I. \end{cases} $$

The CHDG formulation then reads as follows.

**PROBLEM 2.2 (CHDG formulation).** Find $(e_h, h_h, g_h^{\ominus}) \in V_h \times V_h \times G_h^t$ such that, for all $(v_h, w_h, \xi_h) \in V_h \times V_h \times G_h^t$,

$$ \left\{ \begin{array}{l} \iota\kappa(e_h, v_h)_{\mathcal{T}_h} - (h_h, \nabla \times v_h)_{\mathcal{T}_h} - \langle -\frac{1}{2}(g^{\oplus}(e_h, h_h) - g_h^{\ominus}), \pi^t v_h \rangle_{\partial \mathcal{T}_h} = 0, \\ \iota\kappa(h_h, w_h)_{\mathcal{T}_h} + (e_h, \nabla \times w_h)_{\mathcal{T}_h} + \langle \frac{1}{2}n \times (g^{\oplus}(e_h, h_h) + g_h^{\ominus}), \pi^t w_h \rangle_{\partial \mathcal{T}_h} = 0, \end{array} \right. $$

and

$$ (2.9) \qquad \langle g_h^{\ominus} - \Pi(g^{\oplus}(e_h, h_h)), \xi_h \rangle_{\partial \mathcal{T}_h} = \langle b, \xi_h \rangle_{\partial \mathcal{T}_h}, $$

with $g^{\oplus}(e_h, h_h) := \pi^t(e_h) - n \times h_h$.

After eliminating the physical variables, we are left with solving (2.9) on the skeleton with only the variable $g_h^{\ominus}$ as unknown. The elimination implicitly requires the solution of local element-wise problems, where the incoming transmission variable is considered data. The local problem for an element $K$ is given by
**PROBLEM 2.3 (Local problem).** Find $(e_K, h_K) \in \mathcal{P}_p(K) \times \mathcal{P}_p(K)$ such that, for all $(v_K, w_K) \in \mathcal{P}_p(K) \times \mathcal{P}_p(K)$,

$$ (2.10) \qquad \left\{ \begin{array}{l} \iota\kappa(e_K, v_K)_K - (h_K, \nabla \times v_K)_K + \sum_{F \in \mathcal{F}_K} \langle \frac{1}{2} g_{K,F}^{\oplus}(e_K, h_K), \pi_{K,F}^t(v_K) \rangle_F \\ \qquad = \sum_{F \in \mathcal{F}_K} \langle \frac{1}{2} g_{K,F}^{\ominus}, \pi_{K,F}^t(v_K) \rangle_F, \\[2ex] \iota\kappa(h_h, w_K)_K + (e_K, \nabla \times w_K)_K + \sum_{F \in \mathcal{F}_K} \langle \frac{1}{2} n_{K,F} \times g_{K,F}^{\oplus}(e_K, h_K), \pi_{K,F}^t(w_K) \rangle_F \\ \qquad = - \sum_{F \in \mathcal{F}_K} \langle \frac{1}{2} n_{K,F} \times g_{K,F}^{\ominus}, \pi_{K,F}^t(w_K) \rangle_F, \end{array} \right. $$

with $g_{K,F}^{\oplus}(e_K, h_K) := \pi_{K,F}^t(e_K) - n_{K,F} \times h_K$, for given data $g_{K,F}^{\ominus} \in \mathcal{P}_p^t(F)$, for all $F \in \mathcal{F}_K$.

The global problem resulting from the elimination of the physical variables is called the *hybridized problem*. This problem can be written in a convenient abstract form by introducing the global scattering operator $\mathbf{S}: G_h^t \to G_h^t$ satisfying, for each face $F$ of each element $K$,

$$ (2.11) \qquad \mathbf{S}(g_h^{\ominus})|_{K,F} := \pi_{K,F}^t(e_K(g_h^{\ominus})) - n_{K,F} \times h_K(g_h^{\ominus}), $$

where $(e_K, h_K)$ is a solution of the local problem 2.3 with $(g_h^{\ominus}|_{K,F})_{F \in \mathcal{F}_K}$ as surface data. The application of $\mathbf{S}$ involves the solution of all the local problems. The hybridized problem can then be written as

**PROBLEM 2.4 (Hybridized problem I).** Find $g_h^{\ominus} \in G_h^t$ such that, for all $\xi_h \in G_h^t$,

$$ \langle g_h^{\ominus} - \Pi(\mathbf{S}(g_h^{\ominus})), \xi_h \rangle_{\mathcal{T}_h} = \langle \mathbf{b}, \xi_h \rangle_{\mathcal{T}_h}. $$

The problem can be written in operator form by introducing the global projected right-hand side $\mathbf{b}_h$ is the $L^2$-orthogonal projection of $\mathbf{b}$ satisfying $\langle \mathbf{b}_h, \xi_h \rangle_{\partial \mathcal{T}_h} = \langle \mathbf{b}, \xi_h \rangle_{\partial \mathcal{T}_h}$ for all $\xi_h \in G_h^t$. We then obtain

**PROBLEM 2.5 (Hybridized problem II).** Find $g_h^{\ominus} \in G_h^t$ such that

$$ (2.12) \qquad (\mathbf{I} - \Pi\mathbf{S})g_h^{\ominus} = \mathbf{b}_h. $$

Problems 2.4 and 2.5 are equivalent to Problem 2.2 because the element-wise problems are well-posed, see below.

**3. Analysis of the CHDG problems.** In this section, we extend the results obtained in [22] for the Helmholtz equation to the time-harmonic Maxwell equations.

**3.1. Local problem.** We show that the local problems (Problem 2.3) are well-posed, which authorizes the proposed hybridization.

**THEOREM 3.1.** *Problem 2.3 is well-posed.*

*Proof.* It is sufficient to show that $g_{K,F}^{\ominus} = \mathbf{0}$ implies $e_K = h_K = \mathbf{0}$. For brevity, we drop the subscripts $K$ and $F$. Testing system (2.10) with $v = e$ and $w = h$, and avoiding the source terms yields

$$ \begin{aligned} \iota\kappa(e, e)_K - (h, \nabla \times e)_K + \langle \frac{1}{2}(\pi^t e - n \times h), \pi^t e \rangle_{\partial K} &= 0, \\ \iota\kappa(h, h)_K + (e, \nabla \times h)_K + \langle \frac{1}{2}n \times (\pi^t e - n \times h), \pi^t h \rangle_{\partial K} &= 0, \end{aligned} $$

which can be rewritten as

$$ \begin{aligned} \iota\kappa(e, e)_K - (h, \nabla \times e)_K - \frac{1}{2} \langle n \times h, \pi^t e \rangle_{\partial K} + \frac{1}{2} \langle \pi^t e, \pi^t e \rangle_{\partial K} &= 0, \\ \iota\kappa(h, h)_K + (e, \nabla \times h)_K + \frac{1}{2} \langle n \times e, \pi^t h \rangle_{\partial K} + \frac{1}{2} \langle \pi^t h, \pi^t h \rangle_{\partial K} &= 0, \end{aligned} $$

Taking the conjugate of both equations, integrating by part and using surface identities gives

$$ \begin{aligned} -\iota\kappa(e, e)_K - (e, \nabla \times h)_K - \frac{1}{2} \langle n \times e, \pi^t h \rangle_{\partial K} + \frac{1}{2} \langle \pi^t e, \pi^t e \rangle_{\partial K} &= 0, \\ -\iota\kappa(h, h)_K + (h, \nabla \times e)_K + \frac{1}{2} \langle n \times h, \pi^t e \rangle_{\partial K} + \frac{1}{2} \langle \pi^t h, \pi^t h \rangle_{\partial K} &= 0, \end{aligned} $$
Adding the four previous equations yields

$$ \|\pi^t e\|_{\partial K}^2 + \|\pi^t h\|_{\partial K}^2 = 0, $$

which gives $\pi^t e = \pi^t h = 0$ on $\partial K$. Therefore, Problem 2.3 combined with these boundary conditions implies the pair $(e, h)$ is a solution to the strong Maxwell system (2.1) in $K$ with homogeneous boundary conditions. Then, the stability of the Maxwell system implies that the pair $(e, h)$ must also be identically zero in $K$, as desired. □

**3.2. Hybridized problem.** We show that the operator $\Pi S$ is a strict contraction. As a consequence of the Banach fixed-point theorem, the hybridized problem (Problem 2.5) is well-posed, and it can be solved with the fixed-point iteration. We begin with the following lemma pertaining to the scattering operator $S$.

**LEMMA 3.2.** (i) The solution of the local problem 2.3 verifies

$$ (3.1) \quad \sum_{F \in \mathcal{F}_K} \| \pi_{K,F}^t (e_K) - n_{K,F} \times h_K \|_F^2 \\ \qquad + \sum_{F \in \mathcal{F}_K} \| g_{K,F}^{\ominus} - (\pi_{K,F}^t (e_K) + n_{K,F} \times h_K) \|_F^2 = \sum_{F \in \mathcal{F}_K} \| g_{K,F}^{\ominus} \|_F^2. $$

(ii) Moreover, the second term of the left-hand side of (3.1) vanishes if and only if $g_{K,F}^{\ominus} = 0$ for all $F \in \mathcal{F}_K$.

*Proof.* We suppress the subscripts $K, F$ for simplicity.

(i) Taking the test functions $v = e$ and $w = h$, Problem 2.3 reduces to

$$ \begin{aligned} \iota\kappa(e, e)_K - (h, \nabla \times e)_K - \frac{1}{2}\langle n \times h, \pi^t e \rangle_{\partial K} + \frac{1}{2}\langle \pi^t e, \pi^t e \rangle_{\partial K} &= \frac{1}{2}\langle g^{\ominus}, \pi^t e \rangle_{\partial K}, \\ \iota\kappa(h, h)_K + (e, \nabla \times h)_K + \frac{1}{2}\langle n \times e, \pi^t h \rangle_{\partial K} + \frac{1}{2}\langle \pi^t h, \pi^t h \rangle_{\partial K} &= \frac{1}{2}\langle g^{\ominus}, n \times h \rangle_{\partial K}. \end{aligned} $$

Taking the conjugate of both equations, integrating by part and using surface identities gives

$$ \begin{aligned} -\iota\kappa(e, e)_K - (e, \nabla \times h)_K - \frac{1}{2}\langle n \times e, \pi^t h \rangle_{\partial K} + \frac{1}{2}\langle \pi^t e, \pi^t e \rangle_{\partial K} &= \frac{1}{2}\langle \pi^t e, g^{\ominus} \rangle_{\partial K}, \\ -\iota\kappa(h, h)_K + (h, \nabla \times e)_K + \frac{1}{2}\langle n \times h, \pi^t e \rangle_{\partial K} + \frac{1}{2}\langle \pi^t h, \pi^t h \rangle_{\partial K} &= \frac{1}{2}\langle n \times h, g^{\ominus} \rangle_{\partial K}. \end{aligned} $$

Adding the four previous equations yields

$$ \langle \pi^t e, \pi^t e \rangle_{\partial K} + \langle \pi^t h, \pi^t h \rangle_{\partial K} = \frac{1}{2} \langle g^{\ominus}, \pi^t e + n \times h \rangle_{\partial K} + \frac{1}{2} \langle \pi^t e + n \times h, g^{\ominus} \rangle_{\partial K}, $$

and then

$$ 2 \| \pi^t e \|_{\partial K}^2 + 2 \| n \times h \|_{\partial K}^2 = \langle g^{\ominus}, \pi^t e + n \times h \rangle_{\partial K} + \langle \pi^t e + n \times h, g^{\ominus} \rangle_{\partial K}. $$

The left-hand side of this equation is equal to

$$ \| \pi^t e + n \times h \|_{\partial K}^2 + \| \pi^t e - n \times h \|_{\partial K}^2, $$

and the right-hand side is equal to

$$ \| g^{\ominus} \|_{\partial K}^2 + \| \pi^t e + n \times h \|_{\partial K}^2 - \| g^{\ominus} - (\pi^t e + n \times h) \|_{\partial K}^2. $$

We then have

$$ \| \pi^t e - n \times h \|_{\partial K}^2 = \| g^{\ominus} \|_{\partial K}^2 - \| g^{\ominus} - (\pi^t e + n \times h) \|_{\partial K}^2, $$

which yields (3.1).

(ii) If the second term of (3.1) is zero, we have that

$$ g^{\ominus} = \pi^t e + n \times h \quad \text{on } \partial K. $$
Applying this in Problem 2.3 yields that the pair $(e, h)$ solves the following system

$$ \begin{aligned} \iota\kappa(e, v)_K - (h, \nabla \times v)_K - \langle n \times h, \pi^t(v) \rangle_{\partial K} &= 0, \\ \iota\kappa(h, w)_K + (e, \nabla \times w)_K + \langle n \times e, \pi^t(w) \rangle_{\partial K} &= 0, \end{aligned} $$

for all $v, w \in \mathcal{P}_p(K)$. By using integration by parts in both equations, we see that the pair $(e, h)$ satisfies the problem in the strong form. However, there cannot be a nontrivial polynomial solution to the time-harmonic Maxwell equations. Hence, $e = h = 0$, and then $g^{\ominus} = 0$. The converse follows from the well-posedness of the problem (Theorem 3.1). □

We are now in a position to prove the following.

**THEOREM 3.3.** *The scattering operator S is a strict contraction, i.e.*

$$ \|\mathbf{S}(g_h^{\ominus})\|_{\partial T_h} < \|g_h^{\ominus}\|_{\partial T_h}, \quad \forall g_h^{\ominus} \in G_h^t \setminus \{0\}. $$

*Proof.* Let $g_h^{\ominus} \in G_h^t \setminus \{0\}$ be arbitrary. We set $g_{K,F}^{\ominus} := g_h^{\ominus}|_{K,F}$ and $e_K, h_K$ as the corresponding solutions to Problem 2.3. By Lemma 3.2, we have that

$$ (3.2) \qquad \sum_{F \in \mathcal{F}_K} \| \pi_{K,F}^t(e_K) - n_{K,F} \times h_K \|_F^2 \le \sum_{F \in \mathcal{F}_K} \| g_{K,F}^{\ominus} \|_F^2, $$

with equality occurring only if $g_{K,F}^{\ominus} = 0$ for all $F \in \mathcal{F}_K$. Next, by definition of the scattering operator (2.11), (3.2) is equivalent to

$$ \sum_{F \in \mathcal{F}_K} \| \mathbf{S}(g_{K,F}^{\ominus}) \|_F^2 \le \sum_{F \in \mathcal{F}_K} \| g_{K,F}^{\ominus} \|_F^2, $$

again with equality occurring only if $g_{K,F}^{\ominus} = 0$ for all $F \in \mathcal{F}_K$. Since $g_h^{\ominus} \ne 0$, summing this last estimate over all $K \in T_h$ yields the result. □

Now, we consider the contraction property for the exchange operator $\mathbf{\Pi}$.

**THEOREM 3.4.** *The exchange operator $\mathbf{\Pi}$ is a contraction, i.e.*

$$ \|\mathbf{\Pi}(g_h^{\ominus})\|_{\partial T_h} \le \|g_h^{\ominus}\|_{\partial T_h}, \quad \forall g_h^{\ominus} \in G_h^t. $$

*In addition, if $\Gamma_I = \emptyset$, $\mathbf{\Pi}$ is an involution and an isometry.*

*Proof.* These results are straightforward consequences of the definition (2.8). □

Putting together the previous two theorems, we arrive at the following.

**COROLLARY 3.5.** *The operator $\mathbf{\Pi}\mathbf{S}$ is a strict contraction, i.e.*

$$ \|\mathbf{\Pi}\mathbf{S}(g_h^{\ominus})\|_{\partial T_h} < \|g_h^{\ominus}\|_{\partial T_h}, \quad \forall g_h^{\ominus} \in G_h^t \setminus \{0\}. $$

**4. Discrete systems and iterative solvers.** The solution procedure consists in applying an iterative method to the hybridized problem 2.4, which involves solving the local problem 2.3 on all elements for a given incoming transmission variable. In this section, we derive the discrete systems and describe the corresponding iterative procedures. Throughout this section, for vector-valued functions, we use the non-bold counterpart of a symbol to denote a single component (e.g., $n_{K,F,d}$ for the $d^{th}$ component of $n_{K,F}$), and the notation $[\cdot]_i$ refers to the $i^{th}$ entry of an arbitrary vector.

**4.1. Representation of the numerical fields.** The numerical fields are represented with Lagrange basis functions of maximum degree $p$. The discrete unknowns correspond to the values of the Cartesian components of the fields at finite element nodes, which simplifies the implementation, see e.g. [12, 13]. The spatial distribution of the nodes in the reference tetrahedron is defined using a technique described in [31], and the nodes in the reference triangle are chosen to match the face nodes of the tetrahedron. The number of nodes per tetrahedron and per triangle is given by $N_p = (p+1)(p+2)(p+3)/6$ and $N_{fp} = (p+1)(p+2)/2$, respectively.
Denoting the sets of Lagrange basis functions by {$ℓ_{K,i}(\boldsymbol{x})$}$_{i=1...N_p}$ on each element K and {$ℓ_{K,F,i}(\boldsymbol{x})$}$_{i=1...N_{fp}}$ on each face F ∈ $\mathcal{F}_K$, the d-th components of the vector-valued polynomials $\boldsymbol{e}_K$, $\boldsymbol{h}_K$, $\boldsymbol{g}_{K,F}$ are written as

$$ (4.1a) \qquad e_{K,d}(\boldsymbol{x}) = \sum_{i=1}^{N_p} [\boldsymbol{e}_{K,d}]_i \ell_{K,i}(\boldsymbol{x}), \quad \text{on each } K, $$

$$ (4.1b) \qquad h_{K,d}(\boldsymbol{x}) = \sum_{i=1}^{N_p} [\boldsymbol{h}_{K,d}]_i \ell_{K,i}(\boldsymbol{x}), \quad \text{on each } K, $$

$$ (4.1c) \qquad g_{K,F,d}^{\ominus}(\boldsymbol{x}) = \sum_{i=1}^{N_{fp}} [\boldsymbol{g}_{K,F,d}^{\ominus}]_i \ell_{K,F,i}(\boldsymbol{x}), \quad \text{on each } F \in \mathcal{F}_K, $$

for $d = 1, 2, 3$, where $\boldsymbol{e}_{K,d}, \boldsymbol{h}_{K,d} \in \mathbb{C}^{N_p}$ and $\boldsymbol{g}_{K,F,d}^{\ominus} \in \mathbb{C}^{N_{fp}}$ are the corresponding unknown vectors of DOFs.

Since the transmission variables belong to the tangential trace space (2.3) on each face, they can be represented with only two components in local frames associated with the faces. However, the three Cartesian components are stored to simplify the implementation. Indeed, storing only the tangential components would correspond to a 33% reduction in DOFs ($2 \times N_{fp}$ instead of $3 \times N_{fp}$ per face) but would necessitate additional operations and storage to handle the transformation between the tangent space and the full space.

**4.2. Local systems.** To derive the discrete version of the local problem 2.3, we first write it in Cartesian components. First, for a given a vector-valued function $\boldsymbol{v}: \Omega \to \mathbb{C}^3$, we recall the formulas for the d-th component of the various operators involved:

$$ [\boldsymbol{\nabla} \times \boldsymbol{v}]_d = \sum_{e,f} \varepsilon_{def} \partial_{x_e} v_f, \quad [\boldsymbol{n} \times \boldsymbol{v}]_d = \sum_{e,f} \varepsilon_{def} n_e v_f, \quad [\boldsymbol{\pi}^t(\boldsymbol{v})]_d = v_d - n_d \sum_e n_e v_e, $$

for $d = 1, 2, 3$, where $\varepsilon_{def}$ is the Levi-Civita symbol. Simple calculations then yield the following problem.

**PROBLEM 4.1.** Find ($[e_{K,1}, e_{K,2}, e_{K,3}], [h_{K,1}, h_{K,2}, h_{K,3}] \in [\mathcal{P}_p(K)]^3 \times [\mathcal{P}_p(K)]^3$ such that, for all ($[v_{K,1}, v_{K,2}, v_{K,3}], [w_{K,1}, w_{K,2}, w_{K,3}] \in [\mathcal{P}_p(K)]^3 \times [\mathcal{P}_p(K)]^3$,

$$ \left\{ \begin{array}{l} i\kappa(e_{K,d}, v_{K,d})_K - \sum_{e,f} \varepsilon_{def}(h_{K,d}, \partial_{x_e} v_{K,f})_K + \sum_{F \in \mathcal{F}_K} \langle \frac{1}{2} [\boldsymbol{g}_{K,F}^{\oplus}]_d, v_{K,d} \rangle_F \\ \\ = \sum_{F \in \mathcal{F}_K} \langle \frac{1}{2} \boldsymbol{g}_{K,F,d}^{\ominus}, v_{K,d} \rangle_F, \\[2ex] i\kappa(h_{K,d}, w_{K,d})_K + \sum_{e,f} \varepsilon_{def}(e_{K,d}, \partial_{x_e} w_{K,f})_K + \sum_{F \in \mathcal{F}_K} \langle \frac{1}{2} [\boldsymbol{n}_{K,F} \times \boldsymbol{g}_{K,F}^{\oplus}]_d, w_{K,d} \rangle_F \\ \\ = - \sum_{F \in \mathcal{F}_K} \sum_{e,f} \langle \frac{1}{2} \varepsilon_{def} n_{K,F,e} \boldsymbol{g}_{K,F,f}^{\ominus}, w_{K,d} \rangle_F, \end{array} \right. $$

with

$$ \begin{aligned} [\boldsymbol{g}_{K,F}^{\oplus}]_d &= (e_{K,d} - n_{K,F,d} \sum_e n_{K,F,e} e_{K,e}) - \sum_{e,f} \varepsilon_{def} n_{K,F,e} h_{K,f}, \\ [\boldsymbol{n}_{K,F} \times \boldsymbol{g}_{K,F}^{\oplus}]_d &= \sum_{e,f} \varepsilon_{def} n_{K,F,e} e_{K,f} + (h_{K,d} - n_{K,F,d} \sum_e n_{K,F,e} h_{K,e}), \end{aligned} $$

for $d = 1, 2, 3$.

The discrete system is obtained by using (4.1) in the previous problem and by taking the Lagrange functions as test functions. We introduce the local mass and stiffness matrices

$$ \begin{aligned} [\mathbf{M}_K]_{ij} &= (\ell_{K,j}, \ell_{K,i})_K && \text{for } i, j = 1 \dots N_p, \\ [\mathbf{B}_K, \mathbf{F}]_{ij} &= \langle \ell_{K,F,j}, \ell_{K,i} \rangle_F && \text{for } i = 1 \dots N_p \text{ and } j = 1 \dots N_{fp}, \\ [\mathbf{S}_K, d]_{ij} &= (\partial_{x_d} \ell_{K,j}, \ell_{K,i})_K && \text{for } i, j = 1 \dots N_p \text{ and } d = 1, 2, 3, \end{aligned} $$
and the $N_{fp} \times N_p$ restriction matrix $\mathbf{R}_{K,F}$ that maps the nodes associates to an element $K$ to the nodes associated to a face $F \in \mathcal{F}_K$. Taking this notation into account, we arrive at the following version of the local problem.

**PROBLEM 4.2 (Nodal local system).** Find $([\mathbf{e}_{K,1}, \mathbf{e}_{K,2}, \mathbf{e}_{K,3}], [\mathbf{h}_{K,1}, \mathbf{h}_{K,2}, \mathbf{h}_{K,3}]) \in [\mathbb{C}^{N_p}]^3 \times [\mathbb{C}^{N_p}]^3$ such that

$$ \left\{ \begin{array}{l} i\kappa \mathbf{M}_K \mathbf{e}_{K,d} + \sum_{e,f} \varepsilon_{def} \mathbf{S}_{K,e}^{\top} \mathbf{h}_{K,f} + \sum_{F \in \mathcal{F}_K} \frac{1}{2} \mathbf{B}_{K,F} \mathbf{f}_{K,d}^{\mathrm{E}} = \sum_{F \in \mathcal{F}_K} \frac{1}{2} \mathbf{B}_{K,F} \mathbf{s}_{K,d}^{\mathrm{E}}, \\ i\kappa \mathbf{M}_K \mathbf{h}_{K,d} - \sum_{e,f} \varepsilon_{def} \mathbf{S}_{K,e}^{\top} \mathbf{e}_{K,f} + \sum_{F \in \mathcal{F}_K} \frac{1}{2} \mathbf{B}_{K,F} \mathbf{f}_{K,d}^{\mathrm{H}} = - \sum_{F \in \mathcal{F}_K} \frac{1}{2} \mathbf{B}_{K,F} \mathbf{s}_{K,d}^{\mathrm{H}}, \end{array} \right. $$

with

$$ \mathbf{f}_{K,d}^{\mathrm{E}} := \mathbf{R}_{K,F} \left[ \left( \mathbf{e}_{K,d} - n_{K,F,d} \sum_{e} n_{K,F,e} \mathbf{e}_{K,e} \right) - \sum_{e,f} \varepsilon_{def} n_{K,F,e} \mathbf{h}_{K,f} \right], $$

$$ \mathbf{f}_{K,d}^{\mathrm{H}} := \mathbf{R}_{K,F} \left[ \sum_{e,f} \varepsilon_{def} n_{K,F,e} \mathbf{e}_{K,f} + \left( \mathbf{h}_{K,d} - n_{K,F,d} \sum_{e} n_{K,F,e} \mathbf{h}_{K,e} \right) \right], $$

$$ \mathbf{s}_{K,d}^{\mathrm{E}} := \left( \mathbf{g}_{K,F,d}^{\ominus} - n_{K,F,d} \sum_{e} n_{K,F,e} \mathbf{g}_{K,F,e}^{\ominus} \right), $$

$$ \mathbf{s}_{K,d}^{\mathrm{H}} := \sum_{e,f} \varepsilon_{def} n_{K,F,e} \mathbf{g}_{K,F,f}^{\ominus}, $$

for $d = 1, 2, 3$, for given data $\{[\mathbf{g}_{K,F,1}^{\ominus}, \mathbf{g}_{K,F,2}^{\ominus}, \mathbf{g}_{K,F,3}^{\ominus}]\}_{F \in \mathcal{F}_K}$.

In practice, both equations can be multiplied by the inverse of the local mass matrix, and the local matrices associated to any element $K$ can be easily computed to those associated with the reference element, see e.g. [12, 13].

**4.3. Hybridized systems.** A discrete version of the hybridized problem 2.4 is directly obtained by using the representation (4.1c) and by taking the Lagrange basis functions as test functions. By recasting all the discrete values of $\mathbf{g}_h^{\ominus}$ and $\mathbf{b}_h$ into the global vectors $\mathbf{g}$ and $\mathbf{b}$, respectively, we obtain the system

$$ (4.2) \qquad \tilde{\mathbf{A}}\mathbf{g} := \mathbf{M}(\mathbf{I} - \mathbf{\Pi}\mathbf{S})\mathbf{g} = \mathbf{M}\mathbf{b}, $$

where $\mathbf{S}$ and $\mathbf{\Pi}$ are the algebraic counterparts of the scattering and exchange operators in the nodal basis, $\mathbf{I}$ is the identity matrix, $\mathbf{M}$ is a block diagonal matrix with the mass matrices associated to the faces of all the elements.

In our implementation, we have considered two preconditioned versions of System (4.2). For the first one, we precondition the system on the left with the mass matrix.

**PROBLEM 4.3 (Nodal hybridized system).** Find $\mathbf{g} \in \mathbb{C}^{N_{\mathrm{DOF}}}$ satisfying

$$ (4.3) \qquad \mathbf{A}\mathbf{g} := (\mathbf{I} - \mathbf{\Pi}\mathbf{S})\mathbf{g} = \mathbf{b}. $$

For the second version, we have used a symmetric preconditioning with the mass matrix, which corresponds to using modal basis functions that are orthonormal with respect to the inner product $\langle \cdot, \cdot \rangle_{\partial T_h}$ instead of nodal basis functions. Following this interpretation, we introduce the modal versions of the discrete matrices and vectors denoted by $\hat{\mathbf{A}}, \hat{\mathbf{\Pi}}, \hat{\mathbf{S}}, \hat{\mathbf{g}}$ and $\hat{\mathbf{b}}$. In the modal basis, the system then reads as follows.

**PROBLEM 4.4 (Modal hybridized system).** Find $\hat{\mathbf{g}} \in \mathbb{C}^{N_{\mathrm{DOF}}}$ satisfying

$$ (4.4) \qquad \hat{\mathbf{A}}\hat{\mathbf{g}} := (\mathbf{I} - \hat{\mathbf{\Pi}}\hat{\mathbf{S}})\hat{\mathbf{g}} = \hat{\mathbf{b}}. $$

In practice, this system is rewritten with the nodal versions of the matrices and vectors by using a change of basis. Let introduce the change of basis matrix $\mathbf{V}$ that maps the modal components of any vector to the nodal ones. Each column of $\mathbf{V}$ contains the values of one modal basis function
at the finite element nodes, see e.g. [13]. Using the relations $\mathbf{g} = \mathbf{V}\hat{\mathbf{g}}$, $\mathbf{b} = \mathbf{V}\hat{\mathbf{b}}$, $\Pi = \mathbf{V}\hat{\Pi}\mathbf{V}^{-1}$ and $\mathbf{S} = \mathbf{V}\hat{\mathbf{S}}\mathbf{V}^{-1}$ in System (4.4) gives

$$(\mathbf{V}^{-1}\mathbf{A}\mathbf{V})(\mathbf{V}^{-1}\mathbf{g}) = (\mathbf{V}^{-1}\mathbf{b}).$$

Because the mass matrix verifies $\mathbf{M} = \mathbf{V}^{-\top}\mathbf{V}^{-1}$, we also have

$$(\mathbf{V}^{\top}\tilde{\mathbf{A}}\mathbf{V})(\mathbf{V}^{-1}\mathbf{g}) = \mathbf{V}^{\top}(\mathbf{M}\mathbf{b}),$$

which shows the link with the preconditioning technique.

In the next section, we shall use the $L^2$-norm of fields of $G_h^t$, which can be easily computed in both nodal and modal bases. For any $g_h \in G_h^t$, we indeed have

$$\|g_h\|_{\partial\mathcal{T}_h} = \sqrt{\langle g_h, g_h \rangle_{\partial\mathcal{T}_h}} = \|\hat{g}\|_2 = \|g\|_{\mathbf{M}},$$

with the algebraic norms $\|\hat{g}\|_2 := \sqrt{\hat{g}^{\top}\hat{g}}$ and $\|g\|_{\mathbf{M}} := \sqrt{g^{\top}\mathbf{M}g}$. Because the modal basis is orthonormal with respect to $\langle \cdot, \cdot \rangle_{\partial\mathcal{T}_h}$, the modal mass matrix is the identity.

**4.4. Iterative linear solvers.** Following [22], we consider iterative solvers based on the fixed-point, CGNR and GMRES iterations.

*Fixed-point iteration.* The exchange and scattering matrices inherit from the properties of the corresponding operators, which are proved in §3.2. By Corollary 3.5, the matrices $\hat{\Pi}\hat{\mathbf{S}}$ and $\Pi\mathbf{S}$ are strict contractions in the norms $\|\cdot\|_2$ and $\|\cdot\|_{\mathbf{M}}$, respectively. In both cases, the spectral radius of the matrix is strictly smaller than 1, and the system can be solved with the fixed-point iteration. Writing the resulting algorithm with the nodal matrices and vectors, we obtain, for an initial guess $\mathbf{g}^{(0)}$,

$$\mathbf{g}^{(\ell+1)} = \Pi\mathbf{S}\mathbf{g}^{(\ell)} + \mathbf{b}, \quad \ell = 0, 1, \dots$$

This iteration is completely equivalent to the corresponding modal version which can be seen via the change of basis:

$$\underbrace{\mathbf{V}^{-1}\mathbf{g}^{(\ell+1)}}_{\hat{\mathbf{g}}^{(\ell+1)}} = \underbrace{\mathbf{V}^{-1}\Pi\mathbf{S}\mathbf{V}}_{\hat{\Pi}\hat{\mathbf{S}}} \underbrace{(\mathbf{V}^{-1}\mathbf{g}^{(\ell)})}_{\hat{\mathbf{g}}^{(\ell)}} + \underbrace{\mathbf{V}^{-1}\mathbf{b}}_{\mathbf{b}}, \quad \ell = 0, 1, \dots$$

*CGNR iterations.* The conjugate gradient normal method (CGNR) is applied to the nodal system $\mathbf{A}\mathbf{g} = \mathbf{b}$, which is equivalent to applying the conjugate gradient algorithm to the normal equation $\mathbf{A}^*\mathbf{A}\mathbf{g} = \mathbf{A}^*\mathbf{b}$. The CGNR iteration is also applied to the modal system, where the normal equation $\hat{\mathbf{A}}^*\hat{\mathbf{A}}\hat{\mathbf{g}} = \hat{\mathbf{A}}^*\hat{\mathbf{b}}$ is rewritten with the nodal matrices and vectors. These approaches are referred as the *nodal CGNR* and *modal CGNR* solvers, respectively.

These two solvers exhibit qualitatively different numerical properties. With the nodal solver, the iterate $\mathbf{g}^{(\ell)}$ minimizes $\mathbf{g} \mapsto \|\mathbf{A}\mathbf{g} - \mathbf{b}\|_2$ in the affine subspace $\mathbf{g}^{(0)} + \mathcal{K}_\ell(\mathbf{A}^*\mathbf{A}, \mathbf{A}^*\mathbf{r}^{(0)})$ with the initial residual $\mathbf{r}^{(0)} := \mathbf{A}\mathbf{g}^{(0)} - \mathbf{b}$ for an initial guess $\mathbf{g}^{(0)}$, see e.g. [29]. With the modal solver, the iterate $\mathbf{g}^{(\ell)}$ minimizes $\mathbf{g} \mapsto \|\mathbf{A}\mathbf{g} - \mathbf{b}\|_{\mathbf{M}}$ in $\mathbf{g}^{(0)} + \mathcal{K}_\ell(\mathbf{M}^{-1}\mathbf{A}^*\mathbf{M}\mathbf{A}, \mathbf{M}^{-1}\mathbf{A}^*\mathbf{M}\mathbf{r}^{(0)})$.

The modal CGNR solver is given in Algorithm 4.1. Compared to the nodal one, it requires one application of the mass matrix and its inverse at each iteration, as well as the use of the $\mathbf{M}$-norm instead of the 2-norm. These operations are cheap and easy to parallelize because $\mathbf{M}$ is block diagonal with small blocks of size $N_{\text{fp}} \times N_{\text{fp}}$.

*GMRES iterations.* The generalized minimal residual method (GMRES) is applied to both nodal and modal systems, which we refer to as the *nodal GMRES* and *modal GMRES* solvers, respectively. The iterate $\mathbf{g}^{(\ell)}$ minimizes $\mathbf{g} \mapsto \|\mathbf{A}\mathbf{g} - \mathbf{b}\|_2$ and $\mathbf{g} \mapsto \|\mathbf{A}\mathbf{g} - \mathbf{b}\|_{\mathbf{M}}$, respectively, in the affine subspace $\mathbf{g}^{(0)} + \mathcal{K}_\ell(\mathbf{A}, \mathbf{r}^{(0)})$. Again, the cost of a modal iteration is slightly higher than for a nodal iteration because it requires the application of the mass matrix and its inverse, as well as the use of the $\mathbf{M}$ norm instead of the 2-norm. We will also use the standard practice of restarting the GMRES solver and the notation $\text{GMRES}(n)$ means that we restart every $n$ iteration. Finally, we explicitly compute the norm of the residual $\|\mathbf{b} - \mathbf{A}\mathbf{g}\|_2$ instead of using the norm of the residual that is provided automatically during the Householder version of a GMRES iteration. Indeed, for the modal version, this residual is naturally computed in the $\|\cdot\|_{\mathbf{M}}$-norm. In our experiments, we did not remark a major difference between the two so we only present the $\|\cdot\|_2$-norm version.
<caption>**Algorithm 4.1 Modal CGNR iterative solver**</caption>

Initial guess: $\mathbf{g}^{(0)} \in \mathbb{R}^n$
Initial residual: $\mathbf{r}^{(0)} = \mathbf{b} - \mathbf{A}\mathbf{g}^{(0)}$
$\mathbf{z}^{(0)} = \mathbf{M}^{-1}\mathbf{A}^*\mathbf{M}\mathbf{r}^{(0)}$
$\mathbf{p}^{(0)} = \mathbf{z}^{(0)}$
for $\ell = 0, 1, \dots$ until convergence do
&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{q}^{(\ell)} = \mathbf{A}\mathbf{p}^{(\ell)}$
&nbsp;&nbsp;&nbsp;&nbsp;$\alpha^{(\ell)} = \|\mathbf{z}^{(\ell)}\|_{\mathbf{M}}^2/\|\mathbf{q}^{(\ell)}\|_{\mathbf{M}}^2$
&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{g}^{(\ell+1)} = \mathbf{g}^{(\ell)} + \alpha^{(\ell)}\mathbf{p}^{(\ell)}$
&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{r}^{(\ell+1)} = \mathbf{r}^{(\ell)} - \alpha^{(\ell)}\mathbf{q}^{(\ell)}$
&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{z}^{(\ell+1)} = \mathbf{M}^{-1}\mathbf{A}^*\mathbf{M}\mathbf{r}^{(\ell+1)}$
&nbsp;&nbsp;&nbsp;&nbsp;$\beta^{(\ell)} = \|\mathbf{z}^{(\ell+1)}\|_{\mathbf{M}}^2/\|\mathbf{z}^{(\ell)}\|_{\mathbf{M}}^2$
&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{p}^{(\ell+1)} = \mathbf{z}^{(\ell+1)} + \beta^{(\ell)}\mathbf{p}^{(\ell)}$
end for

**5. Numerical results.** In this section, we present numerical results to validate our implementation, to study and compare the different iterative solvers, and to illustrate the applicability of the approach on a realistic case. We have implemented the CHDG method with the different iterative linear solver in a dedicated C++ code already used for the acoustic case in [21]. This code has parallel capabilities with both MPI and OpenMP. The meshes generation and visualization are performed with gmsh [10]. In what follows, the mesh size parameter h is controlled via MeshSizeMin/MeshSizeMax in gmsh.

**5.1. Validation and comparison.** We consider two reference benchmarks corresponding to a free space case (with an analytic solution) and a cavity case (with a quasi-analytic solution) defined the unit cube $\Omega = (0, 1)^3 \subset \mathbb{R}^3$. For both benchmarks, we use a uniform tetrahedral mesh with $h = 0.4$ and polynomial degree $p = 4$. The wavenumber is set to $\kappa = 2.1\pi$.

**Benchmark 1 (Free space).** The reference solution of the first case corresponds to the plane-wave fields given by $e_{\text{ref}}(\mathbf{x}) = e_0 e^{i\kappa_r (d \cdot \mathbf{x})}$ and $h_{\text{ref}}(\mathbf{x}) = -\mathbf{d} \times e_{\text{ref}}(\mathbf{x})$ with the direction vector $\mathbf{d} := (1, 1, 1)^{\top}/\sqrt{3}$ and the amplitude vector $e_0 = (0, 1, -1)^{\top}/\sqrt{2}$. On the boundary we set $\partial\Omega := \Gamma_I$, and impose the exact solution via $s_I := -\mathbf{n} \times (\mathbf{n} \times e_{\text{ref}}) + \mathbf{n} \times h_{\text{ref}}$.

**Benchmark 2 (Cavity).** In the second case, we impose perfect electric conductor (PEC) boundary conditions, i.e. $\partial\Omega := \Gamma_E$ with $s_E = \mathbf{0}$, and the system is driven by an external (volumetric) electric current density $\mathbf{J}_e = (-\mathbf{i}/\kappa, 0, 0)^{\top}$. This problem admits a semi-analytic solution by considering an expansion of the eigenmode solutions for the electric field, i.e.

$$e_{\text{ref}}(\mathbf{x}) = \sum_{k_1, k_2, k_3 \in \mathbb{N}} \alpha_{k_1, k_2, k_3} \begin{pmatrix} a_1 \cos(k_1 \pi x) \sin(k_2 \pi y) \sin(k_3 \pi z) \\ a_2 \sin(k_1 \pi x) \cos(k_2 \pi y) \sin(k_3 \pi z) \\ a_3 \sin(k_1 \pi x) \sin(k_2 \pi y) \cos(k_3 \pi z) \end{pmatrix}$$

with $(a_1, a_2, a_3) \cdot (k_1, k_2, k_3) = 0$ in order to ensure $\nabla \cdot e_{\text{ref}} = 0$. For the chosen current density $\mathbf{J}_e$, the eigenmodes simplify substantially for $k_1 \equiv 0$ and $k_2, k_3$ odd, and the coefficients of the expansion can be calculated explicitly, leading to

$$e_{\text{ref}}(\mathbf{x}) = \sum_{k_2, k_3 \in \mathbb{N} \text{ odd}} \frac{16}{\pi^2 k_2 k_3 (\pi^2 (k_2^2 + k_3^2) - \kappa_r^2)} \begin{pmatrix} \sin(k_2 y) \sin(k_3 z) \\ 0 \\ 0 \end{pmatrix}.$$

In practice, we retain up to $k_2, k_3 = 25$ in the series to represent the reference solution.

**Relative error and residual.** At each iteration $\ell$, we consider the error of the physical fields for $e_h^{(\ell)}, h_h^{(\ell)} \in V_h$ corresponding to the global vector $\mathbf{g}^{(\ell)}$

$$\text{Relative error} := \sqrt{\frac{\|e_h^{(\ell)} - e_{\text{ref}}\|_{L^2(\Omega)}^2 + \|h_h^{(\ell)} - h_{\text{ref}}\|_{L^2(\Omega)}^2}{\|e_{\text{ref}}\|_{L^2(\Omega)}^2 + \|h_{\text{ref}}\|_{L^2(\Omega)}^2}}.$$
<figure><img src="image_3.png" /><figcaption>Fig. 5.1: Plots of Re($e_{ref}$) for the reference solutions. The coloring corresponds to the amplitude of the vector.</figcaption></figure>

and the relative residual in 2-norm

$$ \text{Relative residual} := \frac{\|\mathbf{b} - \mathbf{Ag}^{(\ell)}\|_2}{\|\mathbf{b} - \mathbf{Ag}^{(0)}\|_2}. $$

These quantities are plotted as a function of iterations in Figure 5.2 for both benchmarks.

The $L^2$ relative projection error, corresponding to the relative error between $e_{ref}$, $h_{ref}$ and their $L^2$ projections on the finite element space, is represented by the horizontal dashed lines. We observe that none of the methods reduce the error to the level of the projection error. This is expected since the numerical solution we produce is rigorously equivalent to the DG solution which in turn is not equal to the $L^2$ projection of the true solution. However, the final error produced the CHDG discretization is about twice of the relative projection error.

*Comparison of the iterative solvers.* We now compare the performance of different iterative solvers. Several remarks are in order with respect to the results:

* The fixed-point iteration converges for both test cases as guaranteed by Corollary 3.5. However, we see that this convergence is much slower for the cavity benchmark, which can be partially explained by the fact that the proximity of resonance modes results in worse conditioning at the discrete level.

* Both versions of the CGNR method outperform the fixed-point iteration on both benchmarks. This is more pronounced for the cavity benchmark due to the previous remark. Moreover, the modal CGNR as in (4.1) is more efficient for both benchmarks. As discussed in §4.4, we attribute this to the fact the CGNR is minimizing the residual in the same norm that we have the contraction.

* The sequence of iterates produced by GMRES(30) converges more quickly than both the fixed-point and nodal GGNR variant. In the case of the free space benchmark, the behavior of GMRES(30) and the modal CGNR variant is almost identical. On the other hand, for the cavity benchmark, the modal CGNR variant is more efficient. We also consider a shorter restart interval given by GMRES(5) and observe a clear deterioration in the speed of convergence for both benchmarks. We tested an non-restarted version
<figure><img src="image_4.png" /><figcaption>Fig. 5.2: Results of the different iterative solvers for the Benchmarks 1 (free space) and 2 (cavity).</figcaption></figure>

of GMRES as well which performed practically identically to GMRES(30) for the plane wave benchmark, and was faster than CGNR nodal but slower than CGNR modal for the cavity benchmark.

These results are consistent with those obtained in the acoustic case in [22].

**5.2. Solver results for the realistic test case.** We now present results with a benchmark designed and measured by EADS Aerospatiale Matra Missiles for Workshop EM-JINA 98. Numerical tests have been performed for the same setup in, e.g. [20, 6, 30, 1].

The domain $\Omega$ consists of five continuous segments. The first segment is a straight rectangular section with a length of 10 mm, followed by a circular bend of 35° with a radius of 186 mm. The third segment is another straight section measuring 80 mm in length, which is then followed by a second circular bend, also 35° with a radius of 186 mm. Finally, the domain ends with a straight segment of 100 mm. The domain has a uniform cross-section measuring 84 mm by 110 mm. We impose PEC boundary conditions $s_E = 0$ on all the boundaries except the aperture as shown in Figure 5.3. At the aperture we prescribe an incident free space with a
<figure><img src="image_5.png" /><figcaption>Fig. 5.3: Geometry of the COBRA domain with the aperture in gray (left) and the subdomain partition with 160 subdomains (right).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5.4: Snapshots of $|\operatorname{Re}(e_h)|$ for the COBRA benchmark at four different points during the solution process using modal CGNR, cf. §5.2.</figcaption></figure>

frequency of 10 GHz corresponding to $\kappa_r = 209 \text{ m}^{-1}$. We use a uniform mesh with $h = 6 \times 10^{-3}$ and polynomial degree $p = 3$ which respects approximately 5 cells per wavelength, resulting in 13, 308, 240 DOFs for the nodal eliminated system (4.2).

In order to render this computation tractable, we use a **METIS** [16] mesh partition with 160
<figure><img src="image_7.png" /><figcaption>Fig. 5.5: Relative residual history for the COBRA benchmark.</figcaption></figure>

subdomains as depicted in Figure 5.3. This partition was then used on a local cluster with four nodes, each equipped with 40 Intel Xeon CPU Gold 6230 cores at 2.1 GHz for a total of 160 MPI ranks.

In Figure 5.4 we plot the magnitude of the real part of the electric field for the solution obtained at intermediate solver iterations for the modal CGNR method. After 2000 iterations the solution stabilizes and we obtain a profile similar to the one obtained e.g. in [6], with higher magnitude near the first bent region after the aperture and near the “closed” (PEC) side of the cavity.

The residual history of the different methods is plotted in Figure 5.5. We remark a much slower convergence as compared to the previous test cases, which is reasonable since the linear system size is also much greater and the configuration is more complex. We impose a tolerance of $10^{-3}$ on the relative residual for the stopping criteria. The fixed point iteration failed to converge after 15 hours of wall time corresponding to $4.5 \times 10^4$ iterations and a final relative residual of $2.1 \times 10^{-3}$.

We give more details of the solvers performance in Table 5.1. In terms of iteration count, CGNR modal was the fastest to converge, taking approximately 66% of the iterations of the next faster method (GMRES(30) nodal). However, we see that the average time per iteration is higher for both the CGNR variants, practically double that of GMRES(30) and fixed point. This is reasonable since CGNR requires the action of $\mathbf{A}$ and $\mathbf{A}^*$ at each iteration whereas the other methods only require the action of $\mathbf{A}$. There is also a cost associated with the modal versions compared with their nodal counterparts, which turned out to be around 4% of the overall cost of an iteration for both CGNR and GMRES(30). We emphasize that these results are preliminary and that further optimizations to the code are underway.

<table><caption>Table 5.1: Runtime analysis for the COBRA benchmark of §5.2 to achieve a relative residual of $10^{-3}$. The fixed point iteration failed to achieve this criteria.</caption><thead><tr><th>Method</th><th># iterations</th><th>Total runtime</th><th>Average time/iteration</th></tr></thead><tbody><tr><td>GMRES(30) nodal</td><td>7379</td><td>2.74 hours</td><td>1.34 s</td></tr><tr><td>GMRES(30) modal</td><td>7818</td><td>3.02 hours</td><td>1.39 s</td></tr><tr><td>CGNR nodal</td><td>8838</td><td>6.26 hours</td><td>2.55 s</td></tr><tr><td>CGNR modal</td><td>4876</td><td>3.59 hours</td><td>2.65 s</td></tr><tr><td>Fixed point</td><td>X</td><td>X</td><td>1.29 s</td></tr></tbody></table>
**6. Conclusion.** In this paper, we developed and analyzed an extension of the CHDG method, originally studied in [22] for the Helmholtz equation, to the time-harmonic Maxwell equations. The main mathematical novelty is the introduction of a vector-valued transmission variable that lives in the tangential trace space of element faces. The same types of mathematical arguments lead to a hybridized equation of the form $(\mathbf{I} - \mathbf{\Pi}\mathbf{S})g = \mathbf{b}$ (Problem 2.5). We prove using similar techniques that the operator $\mathbf{\Pi}\mathbf{S}$ is a strict contraction, and that the local sub-problems arising from the hybridization are well-posed. A fixed-point iteration can be used to solve the hybridized problem, which is illustrated in our numerical experiments.

Practically, we consider 3D vector cases, while in [22] only 2D scalar cases were presented. In order to solve the hybridized system efficiently, careful considerations have to be made from an implementation point of view. In this paper we present a nodal implementation, using ideas borrowed from [13] for discretization and mapping between nodal and modal polynomial bases. We validate and study variants of our implementation by solving two reference benchmarks, and we also present a large-scale benchmark with a parallel C++ code.

We compare iterative strategies with CGNR and (restarted) GMRES iterations based on both nodal and modal bases, but implemented only with the nodal basis for computational performance as well as ease of implementation. The modal and nodal versions of the GMRES iterations yield similar results, which are similar or slightly better than those obtained with the modal CGNR iteration. However, when using the restarted version of GMRES, modal CGNR slightly outperforms GMRES. We believe that CGNR can be a competitive approach for practical cases, notably due to its low memory footprint.

In light of the recent work [28] where an extension to heterogeneous media for Helmholtz equation is presented, an obvious extension would be to consider heterogeneous media in the context of Maxwell's equations. We are also working to implement a version of the code that is accelerated by GPUs which could provide further performance gains.

## REFERENCES

[1] N. BOOTLAND, V. DOLEAN, P. JOLIVET, AND P.-H. TOURNIER, *A comparison of coarse spaces for Helmholtz problems in the high frequency regime*, Comput. Math. Appl., 98 (2021), pp. 239–253.

[2] A. BUFFA AND I. PERUGIA, *Discontinuous Galerkin approximation of the Maxwell eigenproblem*, SIAM J. Numer. Anal., 44 (2006), pp. 2198–2226.

[3] L. CAMARGO, B. LOPEZ-RODRIGUEZ, M. OSORIO, AND M. SOLANO, *An HDG method for Maxwell's equations in heterogeneous media*, Comput. Method. Appl. M., 368 (2020), p. 113178.

[4] O. CESSENAT AND B. DESPRES, *Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem*, SIAM J. Numer. Anal., 35 (1998), pp. 255–299.

[5] T. CHAUMONT-FRELET AND A. ERN, *A priori and a posteriori analysis of the discontinuous Galerkin approximation of the time-harmonic Maxwell's equations under minimal regularity assumptions*. HAL preprint 04589791, May 2024.

[6] V. DOLEAN, M. J. GANDER, S. LANTERI, J.-F. LEE, AND Z. PENG, *Effective transmission conditions for domain decomposition methods applied to the time-harmonic curl-curl Maxwell's equations*, J. Comput. Phys., 280 (2015), pp. 232–247.

[7] V. DOLEAN, S. LANTERI, AND R. PERRUSSEL, *A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods*, J. Comput. Phys., 227 (2008), pp. 2044–2072.

[8] M. EL BOUAJAJI, V. DOLEAN, M. J. GANDER, AND S. LANTERI, *Optimized Schwarz methods for the time-harmonic Maxwell equations with damping*, SIAM J. Sci. Comput., 34 (2012), pp. A2048–A2071.

[9] H. S. FURE, S. PERNET, M. SIRDEY, AND S. TORDEUX, *A discontinuous Galerkin Trefftz type method for solving the two dimensional Maxwell equations*, SN Partial Differential Equations and Applications, 1 (2020), pp. 1–25.

[10] C. GEUZAINE AND J.-F. C. REMACLE, *Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities*, Internat. J. Numer. Methods Engrg., 79 (2009), pp. 1309–1331.

[11] Y.-X. HE, L. LI, S. LANTERI, AND T.-Z. HUANG, *Optimized Schwarz algorithms for solving time-harmonic Maxwell's equations discretized by a hybridizable discontinuous Galerkin method*, Computer Physics Communications, 200 (2016), pp. 176–181.

[12] J. S. HESTHAVEN AND T. WARBURTON, *Nodal high-order methods on unstructured grids: I. Time-domain solution of Maxwell's equations*, J. Comput. Phys., 181 (2002), pp. 186–221.

[13] J. S. HESTHAVEN AND T. WARBURTON, *Nodal discontinuous Galerkin methods: algorithms, analysis, and applications*, Springer Science & Business Media, 2007.

[14] R. HIPTMAIR, A. MOIOLA, AND I. PERUGIA, *Error analysis of Trefftz-discontinuous Galerkin methods for the time-harmonic Maxwell equations*, Math. Comp., 82 (2013), pp. 247–268.

[15] T. HUTTUNEN, M. MALINEN, AND P. MONK, *Solving Maxwell's equations using the ultra weak variational*
formulation, J. Comput. Phys., 223 (2007), pp. 731–758.

[16] G. KARYPIS AND V. KUMAR, *METIS: A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices*, September 1998.

[17] A. KLÖCKNER, T. WARBURTON, J. BRIDGE, AND J. S. HESTHAVEN, *Nodal discontinuous Galerkin methods on graphics processors*, J. Comput. Phys., 228 (2009), pp. 7863–7882.

[18] L. LI, S. LANTERI, AND R. PERRUSSEL, *Numerical investigation of a high order hybridizable discontinuous Galerkin method for 2D time-harmonic Maxwell's equations*, COMPEL, 32 (2013), pp. 1112–1138.

[19] 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 equation*, J. Comput. Phys., 256 (2014), pp. 563–581.

[20] J. LIU AND J. JIN, *Scattering analysis of a large body with deep cavities*, IEEE Trans. Antennas Propag., 51 (2003), pp. 1157–1167.

[21] A. MODAVE, *Accelerated iterative DG finite element solvers for large-scale time-harmonic acoustic problems*, in *INTER-NOISE and NOISE-CON Congress and Conference Proceedings*, vol. 270, Institute of Noise Control Engineering, 2024, pp. 1301–1307.

[22] A. MODAVE AND T. CHAUMONT-FRELET, *A hybridizable discontinuous Galerkin method with characteristic variables for Helmholtz problems*, J. Comput. Phys., 493 (2023), p. 112459.

[23] N. C. NGUYEN, J. PERAIRE, AND B. COCKBURN, *Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations*, J. Comput. Phys., 230 (2011), pp. 7151–7175.

[24] S. PERNET, N. SERDIUK, M. SIRDEY, AND S. TORDEUX, *Discontinuous Galerkin method based on Riemann fluxes for the time domain Maxwell system*, tech. report, INRIA Bordeaux- Sud-Ouest, 2021.

[25] S. PERNET, M. SIRDEY, AND S. TORDEUX, *Ultra-weak variational formulation for heterogeneous Maxwell problem in the context of high performance computing*, ESAIM: Proceedings and Surveys, 75 (2023), pp. 96–121.

[26] I. PERUGIA AND D. SCHÖTZAU, *The hp-local discontinuous Galerkin method for low-frequency time-harmonic Maxwell equations*, Math. Comp., 72 (2003), pp. 1179–1214.

[27] I. PERUGIA, D. SCHÖTZAU, AND P. MONK, *Stabilized interior penalty methods for the time-harmonic Maxwell equations*, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 4675–4697.

[28] S. PESCUMA, G. GABARD, T. CHAUMONT-FRELET, AND A. MODAVE, *A hybridizable discontinuous Galerkin method with transmission variables for time-harmonic acoustic problems in heterogeneous media*, J. Comput. Phys., 534 (2025), p. 114009.

[29] Y. SAAD, *Iterative methods for sparse linear systems*, Society for Industrial and Applied Mathematics, Philadelphia, PA, second ed., 2003.

[30] A. VION AND C. GEUZAINE, *Improved sweeping preconditioners for domain decomposition algorithms applied to time-harmonic Helmholtz and Maxwell problems*, vol. 61 of ESAIM Proc. Surveys, 2018, pp. 93–111.

[31] T. WARBURTON, *An explicit construction of interpolation nodes on the simplex*, Journal of engineering mathematics, 56 (2006), pp. 247–262.