# A hybridizable discontinuous Galerkin method with transmission variables for time-harmonic acoustic problems in heterogeneous media

Simone Pescuma, Gwenael Gabard, Théophile Chaumont-Frelet, Axel Modave

► To cite this version:

Simone Pescuma, Gwenael Gabard, Théophile Chaumont-Frelet, Axel Modave. A hybridizable discontinuous Galerkin method with transmission variables for time-harmonic acoustic problems in heterogeneous media. Journal of Computational Physics, 2025, 534, pp.114009. 10.1016/j.jcp.2025.114009. hal-04821539v2

HAL Id: hal-04821539

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

Submitted on 23 Apr 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 acoustic problems in heterogeneous media*

S. Pescuma¹, G. Gabard², T. Chaumont-Frelet³, and A. Modave¹

¹POEMS, CNRS, Inria, ENSTA, Institut Polytechnique de Paris, 91120 Palaiseau, France,
simone.pescuma@ensta.fr, axel.modave@ensta.fr

²Laboratoire d'Acoustique de l'Université du Mans (LAUM), IA-GS, CNRS, 72085 Le Mans, France
gwenael.gabard@univ-lemans.fr

³Inria, Laboratoire Paul Painlevé, Université de Lille, 59655 Villeneuve-d'Ascq, France
theophile.chaumont@inria.fr

## Abstract

We consider the finite element solution of time-harmonic wave propagation problems in heterogeneous media with hybridizable discontinuous Galerkin (HDG) methods. In the case of homogeneous media, it has been observed that the iterative solution of the linear system can be accelerated by hybridizing with transmission variables instead of numerical traces, as performed in standard approaches. In this work, we extend the HDG method with transmission variables, which is called the CHDG method, to the heterogeneous case with piecewise constant physical coefficients. In particular, we consider formulations with standard upwind and general symmetric fluxes. The CHDG hybridized system can be written as a fixed-point problem, which can be solved with stationary iterative schemes for a class of symmetric fluxes. The standard HDG and CHDG methods are systematically studied with the different numerical fluxes by considering a series of 2D numerical benchmarks. The convergence of standard iterative schemes is always faster with the extended CHDG method than with the standard HDG methods, with upwind and scalar symmetric fluxes.

## 1 Introduction

Discontinuous Galerkin (DG) methods are widely used for solving boundary value problems because of their ability to provide high-fidelity numerical solutions for complicated geometric and physical configurations, see e.g. [7, 11, 31, 33, 43]. In the context of time-harmonic acoustic problems, they allow the use of high-order polynomial basis functions, which limits the dispersion error that occurs when considering high-frequency cases, e.g. [6, 10, 16, 32]. However, from a computational point of view, these methods require the solution of large sparse unstructured linear systems. On the one hand, direct solvers require huge amounts of computation and are complicated to run efficiently on parallel computers. On the other hand, iterative solvers require much less memory and allow efficient parallel implementations, but they may exhibit slow convergence due to intrinsic properties of the time-harmonic problems [21]. To speed up convergence, preconditioning techniques and domain decomposition methods (DDM) have been and are being intensively studied, see e.g. [5, 19, 20, 24].

<footnote>*Published in Journal of Computational Physics (doi: 10.1016/j.jcp.2025.114009). Distributed under Creative Commons CC BY 4.0 license.</footnote>
Here, we consider the iterative solution of time-harmonic acoustic models with Hybridizable Dis-
continuous Galerkin (HDG) methods. The principle of HDG methods is to introduce a hybrid
variable on the mesh skeleton in order to decouple the physical variables at the interface between
the elements. Eliminating the physical variables within every element then leads to a hybridized
system where the only unknown is the hybrid variable itself, see e.g. [2, 3, 9, 13, 26, 37, 39]. This
approach modifies the size, structure, and conditioning of the global system that is actually solved,
potentially allowing for more efficient solution procedures. Once the hybrid variables are calcu-
lated, the physical variables inside each element can be recovered with a simple and inexpensive
post-processing step.

In standard HDG methods, the so-called numerical trace is taken as the hybrid variable. Here,
we consider an alternative hybridization strategy where the hybrid variable corresponds to a
characteristic variable, which can be interpreted as a transmission variable or a Robin trace. This
approach, which we call CHDG, was recently introduced by two of the authors in [34]. It is related
to some types of non-overlapping DDM [15, 17, 22, 35], Ultra Weak Variational Formulations
(UWVF) [1, 4, 8, 23, 29, 38] and a hybridization technique studied in [36].

In Ref. [34], the performance of the DG, standard HDG and CHDG methods for fast iterative procedures has been compared in the cases of time-harmonic scalar wave propagation problems set in a homogeneous medium. In particular, it has been observed that the convergence of classical iterative schemes is faster with the CHDG method than with the standard DG and HDG methods.

In this work, we investigate extensions of the CHDG approach to scalar time-harmonic problems in heterogeneous media with material coefficients that are (possibly) discontinuous at the interfaces between the elements. We consider DG schemes using either standard upwind fluxes, or a general family of symmetric fluxes based on differential operators defined on the interfaces between elements. The corresponding transmission variables are used as hybrid variable in the hybridization process. By extending the results in [34], we show that the hybridized system is well-posed with both types of fluxes, and that it can be written in the form $(I - \Pi S)g = b$, with the set of all hybrid variables $g$, an exchange operator $\Pi$, a scattering operator $S$, and a global right-hand side $b$. In the case of symmetric fluxes, we prove under reasonable assumptions that the operator $\Pi S$ is a strict contraction, and that the hybridized system can be solved with fixed-point iterations.

In practice, Krylov-type iterative solvers are widely used to solve time-harmonic problems. Here, by using a set of academic benchmarks, we study and compare the performance of the GMRES iteration and the CGNR iteration (i.e. the conjugate gradient applied to the normal equation) for solving the hybridized systems. More specifically, we consider hybridized systems corresponding to the standard HDG and CHDG methods with different numerical fluxes. We observe that, as in [34], CHDG always requires fewer iterations than the standard HDG to reach a given accuracy with the GMRES and CGNR iterations. Moreover, the convergence of standard upwind fluxes and symmetric fluxes is similar, and using symmetric fluxes with higher-order derivatives does not seem to speed up the convergence.

The remainder of this paper is structured as follows. In Section 2, we introduce the problem, the standard DG method, and the upwind and general symmetric numerical fluxes. The standard HDG and CHDG methods are described and studied in Sections 3 and 4. In particular, in the latter case, the strict contraction of the operator $\Pi S$ in the CHDG system is proven for a class of symmetric numerical fluxes. In Section 5, the hybridizable methods and the numerical fluxes are systematically assessed by considering a set of 2D reference numerical benchmarks with different iterative schemes. Conclusions and perspectives are proposed in Section 6.
## 2 Discontinuous Galerkin methods

Let $\Omega \subset \mathbb{R}^d$, with $d = 2$ or 3, be a Lipschitz polytopal domain. The boundary $\partial\Omega$ of the domain is partitioned into three non-overlapping polytopal Lipschitz subsets $\Gamma_D$, $\Gamma_N$ and $\Gamma_R$, corresponding to the Dirichlet, Neumann and Robin boundaries, respectively. We consider the following time-harmonic wave propagation problem:

$$ \left\{ \begin{array}{ll} -i\kappa\eta^{-1}p + \nabla \cdot \mathbf{u} = 0, & \text{in } \Omega, \\ -i\kappa\eta\mathbf{u} + \nabla p = \mathbf{0}, & \text{in } \Omega, \\ p = s_D, & \text{on } \Gamma_D, \\ \mathbf{n} \cdot \mathbf{u} = s_N, & \text{on } \Gamma_N, \\ p - \eta\mathbf{n} \cdot \mathbf{u} = s_R, & \text{on } \Gamma_R, \end{array} \right. \qquad (2.1) $$

written for the unknown scalar field $p : \Omega \to \mathbb{C}$ and vector field $\mathbf{u} : \Omega \to \mathbb{C}^d$ corresponding to the acoustic pressure and velocity. On the boundary, we have introduced the source terms $s_D : \Gamma_D \to \mathbb{C}$, $s_N : \Gamma_N \to \mathbb{C}$ and $s_R : \Gamma_R \to \mathbb{C}$. The propagation medium is characterised by its density $\rho : \Omega \to \mathbb{R}$ and sound speed $c : \Omega \to \mathbb{R}$. These two quantities are strictly positive and assumed to be piecewise constant. The problem (2.1) is written in the frequency domain, assuming a time dependence $e^{-i\omega t}$ for the solution, where $\omega$ is the real-valued angular frequency. We also introduced the wavenumber $\kappa : \Omega \to \mathbb{R}$ and the acoustic impedance $\eta : \Omega \to \mathbb{R}$ such that $\kappa = \omega/c$ and $\eta = \rho c$. The vector field $\mathbf{n} : \partial\Omega \to \mathbb{R}^d$ is the unit outward normal to $\Omega$. For brevity, we do not consider source terms on the right-hand sides of the first two equations in (2.1), but these could be included without difficulty.

### 2.1 Mesh, approximation spaces and inner products

To solve problem (2.1), we use 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\}$, and the collection of faces is denoted by $\mathcal{F}_h$. The collection of faces of an element $K$ is denoted by $\mathcal{F}_K$.

For simplicity, we fix a polynomial degree $p \ge 0$ and introduce the following approximation spaces for the scalar and vector fields:

$$ V_h := \prod_{K \in \mathcal{T}_h} \mathcal{P}_p(K) \quad \text{and} \quad \mathbf{V}_h := \prod_{K \in \mathcal{T}_h} \mathcal{P}_p(K), $$

where $\mathcal{P}_p(\cdot)$ and $\mathcal{P}_p(\cdot)$ denote spaces of scalar and vector complex-valued polynomial of degree smaller or equal to $p$. The restrictions of $u_h \in V_h$ and $\mathbf{u}_h \in \mathbf{V}_h$ on $K$ are denoted by $u_K$ and $\mathbf{u}_K$, respectively. The coefficients $\kappa$ and $\eta$ are assumed to be constant on each element. Their values on $K$ are denoted by $\kappa_K$ and $\eta_K$, respectively. For a face $F$ of $K$, $\mathbf{n}_{K,F}$ is the unit normal on $F$ pointing outside $K$.

We introduce the following sesquilinear forms

$$ (u, v)_K := \int_K u \bar{v} \, d\mathbf{x}, \qquad (\mathbf{u}, \mathbf{v})_K := \int_K \mathbf{u} \cdot \bar{\mathbf{v}} \, d\mathbf{x}, \qquad \langle u, v \rangle_{\partial K} := \sum_{F \in \mathcal{F}_K} \int_F u \bar{v} \, d\sigma(\mathbf{x}), $$

$$ (u, v)_{\mathcal{T}_h} := \sum_{K \in \mathcal{T}_h} (u, v)_K, \qquad (\mathbf{u}, \mathbf{v})_{\mathcal{T}_h} := \sum_{K \in \mathcal{T}_h} (\mathbf{u}, \mathbf{v})_K, \qquad \langle u, v \rangle_{\partial \mathcal{T}_h} := \sum_{K \in \mathcal{T}_h} \langle u, v \rangle_{\partial K}, $$
where the quantities used in the surface integral $\langle \cdot, \cdot \rangle_{\partial K}$ correspond to the restriction of fields or coefficients defined on $K$ (e.g. $v_K$, $v_K$, $\kappa_K$ and $\eta_K$) or quantities associated with the faces of $K$ (e.g. $n_{K,F}$ with $F \in \mathcal{F}_K$), unless explicitly specified.

## 2.2 Variational formulation of the problem

The general discontinuous Galerkin (DG) formulation of system (2.1) reads:

**Problem 2.1.** Find $(p_h, u_h) \in V_h \times V_h$ such that, for all $(q_h, v_h) \in V_h \times V_h$,

$$ \begin{cases} -\iota(\kappa\eta^{-1}p_h, q_h)_{\mathcal{T}_h} - (u_h, \nabla q_h)_{\mathcal{T}_h} + \langle \mathbf{n} \cdot \hat{\mathbf{u}}, q_h \rangle_{\partial\mathcal{T}_h} = 0, \\ -\iota(\kappa\eta u_h, v_h)_{\mathcal{T}_h} - (p_h, \nabla \cdot v_h)_{\mathcal{T}_h} + \langle \hat{p}, \mathbf{n} \cdot v_h \rangle_{\partial\mathcal{T}_h} = 0, \end{cases} $$

with the numerical fluxes $\hat{p}$ and $\mathbf{n} \cdot \hat{\mathbf{u}}$ to be defined.

The numerical fluxes are defined face by face. At each *interior face* (i.e. $F \not\subset \partial\Omega$), they depend on the values of the physical fields and coefficients associated to both neighboring elements. At each *boundary face* (i.e. $F \subset \partial\Omega$), a specific definition is used to prescribe a boundary condition. The numerical fluxes then depend on the boundary data. Hereafter, two kinds of numerical fluxes are considered: standard upwind fluxes and general symmetric fluxes.

By convention, if $F$ is an *interior face* of an element $K$, then $K'$ is the neighboring element of $K$ sharing this face. If $F$ is a *boundary face*, we define $\eta_{K'} := \eta_K$ to simplify the presentation, even if there is no neighboring element $K'$.

## 2.3 Upwind numerical fluxes

Standard upwind numerical fluxes are obtained by solving a Riemann problem associated to each interface between two elements, see e.g. [30, 46]. For the wave propagation problem (2.1) with discontinuous coefficients between elements, the fluxes are given by Equations (9.58) at the end of Section 9.9 of Ref. [30]. For a given face $F$, the *upwind numerical fluxes* are defined as

$$ \begin{cases} \hat{p}_F := \frac{1}{\eta_K + \eta_{K'}} \left( \eta_{K'} g_{K,F}^{\oplus} + \eta_K g_{K,F}^{\ominus} \right), \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F := \frac{1}{\eta_K + \eta_{K'}} \left( g_{K,F}^{\oplus} - g_{K,F}^{\ominus} \right), \end{cases} \qquad (2.2) $$

where $g_{K,F}^{\oplus}$ and $g_{K,F}^{\ominus}$ are the outgoing and incoming transmission variables, respectively. In the standard terminology of Riemann solvers, the transmission variables are generally called characteristic variables, see e.g. [28, 44]. The *outgoing transmission variable* depends on the physical variables and coefficients of the considered element $K$. It is defined as

$$ g_{K,F}^{\oplus} := p_K + \eta_K \mathbf{n}_{K,F} \cdot \mathbf{u}_K. $$
The incoming transmission variable depends on whether there is a neighboring element $K'$ or not. In the latter case, it depends on the boundary condition. It is defined as

$$ g_{K,F}^{\ominus} := \begin{cases} p_{K'} - \eta_{K'} \mathbf{n}_{K,F} \cdot \mathbf{u}_{K'} = g_{K',F}^{\oplus}, & \text{if } F \not\subset \partial\Omega, \\ 2s_{\mathrm{D}} - g_{K,F}^{\oplus}, & \text{if } F \subset \Gamma_{\mathrm{D}}, \\ g_{K,F}^{\oplus} - 2\eta_{K}s_{\mathrm{N}}, & \text{if } F \subset \Gamma_{\mathrm{N}}, \\ s_{\mathrm{R}}, & \text{if } F \subset \Gamma_{\mathrm{R}}. \end{cases} \quad (2.3) $$

The upwind numerical fluxes can be written more explicitly as

$$ \left\{ \begin{aligned} \hat{p}_F &= \frac{\eta_K \eta_{K'}}{\eta_K + \eta_{K'}} \left( \frac{1}{\eta_K} p_K + \frac{1}{\eta_{K'}} p_{K'} + \mathbf{n}_{K,F} \cdot (\mathbf{u}_K - \mathbf{u}_{K'}) \right), & \text{if } F \not\subset \partial\Omega, \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F &= \frac{1}{\eta_K + \eta_{K'}} (\mathbf{n}_{K,F} \cdot (\eta_K \mathbf{u}_K + \eta_{K'} \mathbf{u}_{K'}) + p_K - p_{K'}), & \\ \hat{p}_F &= s_{\mathrm{D}}, & \text{if } F \subset \Gamma_{\mathrm{D}}, \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F &= \mathbf{n}_{K,F} \cdot \mathbf{u}_K + \frac{1}{\eta_K} (p_K - s_{\mathrm{D}}), & \\ \hat{p}_F &= p_K + \eta_K (\mathbf{n}_{K,F} \cdot \mathbf{u}_K - s_{\mathrm{N}}), & \text{if } F \subset \Gamma_{\mathrm{N}}, \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F &= s_{\mathrm{N}}, & \\ \hat{p}_F &= \frac{1}{2} (p_K + \eta_K \mathbf{n}_{K,F} \cdot \mathbf{u}_K + s_{\mathrm{R}}), & \text{if } F \subset \Gamma_{\mathrm{R}}, \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F &= \frac{1}{2\eta_K} (p_K + \eta_K \mathbf{n}_{K,F} \cdot \mathbf{u}_K - s_{\mathrm{R}}), & \end{aligned} \right. $$

These fluxes satisfy the relation

$$ \hat{p}_F + \eta_K \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F = p_K + \eta_K \mathbf{n}_{K,F} \cdot \mathbf{u}_K, \quad (2.4) $$

which will be useful to derive the standard HDG formulation.

## 2.4 Symmetric numerical fluxes

We now consider a general family of numerical fluxes where the coefficients do not depend on the elements, but only on the faces. For each face $F$, we introduce a general operator $\mathcal{A}_F$ that verifies the following assumption where $\mathcal{P}_p(F)$ denotes the space of polynomials of degree at most $p$ defined on $F$.

**Assumption 2.2.** For each face $F$, the linear operator $\mathcal{A}_F : \mathcal{P}_p(F) \to \mathcal{P}_p(F)$ is positive and self-adjoint.

For a given face $F$, the symmetric numerical fluxes are defined as

$$ \left\{ \begin{aligned} \hat{p}_F &:= \frac{1}{2} \mathcal{A}_F (g_{K,F}^{\oplus} + g_{K,F}^{\ominus}), \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F &:= \frac{1}{2} (g_{K,F}^{\oplus} - g_{K,F}^{\ominus}), \end{aligned} \right. \quad (2.5) $$

where the outgoing transmission variable is defined as

$$ g_{K,F}^{\oplus} := \mathcal{A}_F^{-1} p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K, $$
and the incoming transmission variable is defined as

$$ g_{K,F}^{\ominus} := \begin{cases} \mathcal{A}_F^{-1} p_{K'} - \mathbf{n}_{K,F} \cdot \mathbf{u}_{K'} = g_{K',F}^{\oplus}, & \text{if } F \not\subset \partial\Omega, \\ 2\mathcal{A}_F^{-1} s_D - g_{K,F}^{\oplus}, & \text{if } F \subset \Gamma_D, \\ g_{K,F}^{\oplus} - 2s_N, & \text{if } F \subset \Gamma_N, \\ \mathcal{B}_{K,F,+}^{-1} (\mathcal{B}_{K,F,-} g_{K,F}^{\oplus} + 2\eta_K^{-1} s_R), & \text{if } F \subset \Gamma_R, \end{cases} \quad (2.6) $$

with $\mathcal{B}_{K,F,\pm} := 1 \pm \eta_K^{-1} \mathcal{A}_F$. The symmetric numerical fluxes can be written more explicitly as

$$ \left\{ \begin{aligned} \hat{p}_F &= \frac{p_K + p_{K'}}{2} + \mathcal{A}_F \left( \mathbf{n}_{K,F} \cdot \frac{\mathbf{u}_K - \mathbf{u}_{K'}}{2} \right), & \text{if } F \not\subset \partial\Omega, \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F &= \mathbf{n}_{K,F} \cdot \frac{\mathbf{u}_K + \mathbf{u}_{K'}}{2} + \mathcal{A}_F^{-1} \left( \frac{p_K - p_{K'}}{2} \right), \\ \hat{p}_F &= s_D, & \text{if } F \subset \Gamma_D, \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F &= \mathbf{n}_{K,F} \cdot \mathbf{u}_K + \mathcal{A}_F^{-1} (p_K - s_D), \\ \hat{p}_F &= p_K + \mathcal{A}_F (\mathbf{n}_{K,F} \cdot \mathbf{u}_K - s_N), & \text{if } F \subset \Gamma_N, \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F &= s_N, \\ \hat{p}_F &= \mathcal{B}_{K,F,+}^{-1} (p_K + \mathcal{A}_F (\mathbf{n}_{K,F} \cdot \mathbf{u}_K) + \eta_K^{-1} \mathcal{A}_F s_R), & \text{if } F \subset \Gamma_R. \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F &= \eta_K^{-1} \mathcal{B}_{K,F,+}^{-1} (p_K + \mathcal{A}_F (\mathbf{n}_{K,F} \cdot \mathbf{u}_K) - s_R), \end{aligned} \right. $$

Similar to the upwind fluxes, the following relation holds

$$ \hat{p}_F + \mathcal{A}_F (\mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F) = p_K + \mathcal{A}_F (\mathbf{n}_{K,F} \cdot \mathbf{u}_K). \quad (2.7) $$

It will be used in the derivation of the standard HDG formulation.

## Zeroth-order symmetric fluxes

If $\mathcal{A}_F = \mu_F$ is a strictly positive real parameter, the numerical fluxes simply read

$$ \left\{ \begin{aligned} \hat{p}_F &:= \frac{1}{2} \mu_F (g_{K,F}^{\oplus} + g_{K,F}^{\ominus}), \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F &:= \frac{1}{2} (g_{K,F}^{\oplus} - g_{K,F}^{\ominus}), \end{aligned} \right. $$

with

$$ \begin{aligned} g_{K,F}^{\oplus} &:= \frac{1}{\mu_F} p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K, \\ g_{K,F}^{\ominus} &:= \frac{1}{\mu_F} p_{K'} - \mathbf{n}_{K,F} \cdot \mathbf{u}_{K'}, & \text{if } F \not\subset \partial\Omega. \end{aligned} $$

The definitions of the incoming transmission variables when $F \subset \partial\Omega$ are adapted accordingly. Note that these fluxes are equivalent to the standard upwind fluxes if the impedance at the interfaces between elements is continuous, and $\mu_F = \eta_K = \eta_{K'}$ for neighboring elements K and K', regardless of whether the wavenumber is continuous or not. We also point out that these transmission variables correspond to the lowest-order absorbing boundary conditions approximating Dirichlet-to-Neumann (DtN) maps.
## Second-order symmetric fluxes

By analogy with non-overlapping domain decomposition methods (e.g. [35]), the operator $\mathcal{A}_F$ can be defined by using transmission operators based on more involved domain truncation techniques, such as high-order absorbing boundary conditions. We consider here an operator with second-order partial derivatives used in a second-order absorbing boundary conditions in [18]. This operator is implicitly defined by

$$ \langle \mathcal{A}_F \phi_F, v_F \rangle_F + \frac{1}{2\kappa_F^2} \langle \nabla_F (\mathcal{A}_F \phi_F), \nabla_F v_F \rangle_F = \mu_F \langle \phi_F, v_F \rangle_F, \quad \forall \phi_F, v_F \in \mathcal{P}_p(F), \quad (2.8) $$

where the real parameter $\mu_F$ is strictly positive and $\nabla_F$ is the gradient operator on $F$. The operator $\mathcal{A}_F$ is then a discrete realization of

$$ \mathcal{A}_F := \mu_F \left( 1 - \frac{1}{2\kappa_F^2} \Delta_F \right)^{-1}, $$

where $\Delta_F$ is the Laplace-Beltrami operator on $F$, completed by homogeneous Neumann boundary conditions at the extremities of $F$. It satisfies Assumption 2.2.

**Proposition 2.3** (Positivity and self-adjointness of the second-order operator). The operator $\mathcal{A}_F$ (2.8) defined on an edge $F$ is positive and self-adjoint.

*Proof.* Because we are in a finite dimensional framework, it is sufficient to prove that $\mathcal{A}_F^{-1}$ is positive and self-adjoint. First, we prove that $\mathcal{A}_F^{-1}$ is positive. To do so, given $\phi_F \in \mathcal{P}_p(F)$, we let $u_F := \mathcal{A}_F \phi_F$, and pick $v_F = u_F$ in (2.8). This gives

$$ 0 \le \langle u_F, u_F \rangle_F + \frac{1}{2\kappa_F^2} \langle \nabla_F u_F, \nabla_F u_F \rangle_F = \mu_F \langle \mathcal{A}_F^{-1} u_F, u_F \rangle_F, $$

which is the desired inequality, as $\mu_F > 0$. Then, we prove that $\mathcal{A}_F^{-1}$ is self-adjoint, which follows from the fact that the left-hand side in (2.8) is an Hermitian form. Indeed,

$$ \begin{aligned} \langle \mathcal{A}_F^{-1} u_F, v_F \rangle_F &= \frac{1}{\mu_F} \langle u_F, v_F \rangle_F + \frac{1}{2\kappa_F^2 \mu_F} \langle \nabla_F u_F, \nabla_F v_F \rangle_F, \\ &= \frac{1}{\mu_F} \langle v_F, u_F \rangle_F + \frac{1}{2\kappa_F^2 \mu_F} \langle \nabla_F v_F, \nabla_F u_F \rangle_F, \\ &= \langle \mathcal{A}_F^{-1} v_F, u_F \rangle_F = \langle u_F, \mathcal{A}_F^{-1} v_F \rangle_F \end{aligned} $$

for all $u_F, v_F \in \mathcal{P}_p(F)$. $\square$

## 3 Hybridization with numerical trace

In HDG methods, the hybridization consists in introducing an additional variable on the mesh skeleton to decouple the physical unknowns defined within neighboring elements. The physical unknowns are then eliminated to obtain a hybridized system where the unknowns are associated only with the additional variable, called the *hybrid variable*. This step requires solving local element-wise problems.

We emphasize that this process does not affect the accuracy of the numerical solution, since it can be interpreted as an algebraic manipulation leading to an equivalent linear system. When this
linear system is solved with a direct solver, the physical unknowns that can be reconstructed are those of the original DG system. However, the hybridized system has different algebraic properties that affect the iterative numerical solution.

## 3.1 Standard HDG formulations

In the standard HDG methods, the hybrid variable, denoted $\hat{p}_h$, corresponds to the numerical flux $\hat{p}$, frequently called the numerical trace in the literature, see e.g. [27]. It belongs to the space of polynomials of order p on each face F of the mesh, denoted by

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

For any field $\hat{p}_h \in \hat{V}_h$, there is one set of scalar unknowns associated to each face of the mesh.

For the case with symmetric numerical fluxes, the variable $\hat{p}_h$ is defined by equation (2.5). Using equation (2.7) in Problem 2.1 then leads to the following HDG formulation.

**Problem 3.1 (HDG formulation with symmetric fluxes).** Find $(p_h, \mathbf{u}_h, \hat{p}_h) \in V_h \times \mathbf{V}_h \times \hat{V}_h$ such that, for all $(q_h, \mathbf{v}_h, \hat{q}_h) \in V_h \times \mathbf{V}_h \times \hat{V}_h$,

$$ \left\{ \begin{array}{l} -i(\kappa \eta^{-1} p_h, q_h)_{\mathcal{T}_h} - (\mathbf{u}_h, \nabla q_h)_{\mathcal{T}_h} + \langle \mathbf{n} \cdot \mathbf{u}_h + \mathcal{A}^{-1}(p_h - \hat{p}_h), q_h \rangle_{\partial \mathcal{T}_h} = 0, \\ -i(\kappa \eta \mathbf{u}_h, \mathbf{v}_h)_{\mathcal{T}_h} - (p_h, \nabla \cdot \mathbf{v}_h)_{\mathcal{T}_h} + \langle \hat{p}_h, \mathbf{n} \cdot \mathbf{v}_h \rangle_{\partial \mathcal{T}_h} = 0, \end{array} \right. $$

and

$$ \langle \hat{p}_h, \hat{q}_h \rangle_{\mathcal{F}_h} - \langle \frac{1}{2}(p_h + \mathcal{A}(\mathbf{n} \cdot \mathbf{u}_h)), \hat{q}_h \rangle_{\partial \mathcal{T}_h \setminus \partial \Omega} \\ - \langle p_h + \mathcal{A}(\mathbf{n} \cdot \mathbf{u}_h), \hat{q}_h \rangle_{\Gamma_N} - \langle (1 + \eta^{-1} \mathcal{A})^{-1}(p_h + \mathcal{A}(\mathbf{n} \cdot \mathbf{u}_h)), \hat{q}_h \rangle_{\Gamma_R} \\ = \langle s_D, \hat{q}_h \rangle_{\Gamma_D} - \langle \mathcal{A}s_N, \hat{q}_h \rangle_{\Gamma_N} + \langle (1 + \eta^{-1} \mathcal{A})^{-1}(\eta^{-1} \mathcal{A}s_R), \hat{q}_h \rangle_{\Gamma_R}. $$

*Remark 3.2 (Case with upwind fluxes).* For the case with the upwind numerical fluxes, the variable $\hat{p}_h$ is defined by equation (2.2), and the HDG formulation is obtained similarly by using equation (2.4) in Problem 2.1. We refer to [12, 14] for the complete description of the HDG formulation with upwind fluxes.

## 3.2 Local element-wise discrete problems

In the hybridization procedure, the physical fields $p_h$ and $\mathbf{u}_h$ are eliminated by solving local element-wise problems, where the numerical trace $\hat{p}_h$ is considered as a given datum.

For each element K, the local problem corresponding to the symmetric fluxes reads as follows.

**Problem 3.3.** Find $(p_K, \mathbf{u}_K) \in \mathcal{P}_p(K) \times \mathcal{P}_p(K)$ such that, for all $(q_K, \mathbf{v}_K) \in \mathcal{P}_p(K) \times \mathcal{P}_p(K)$,

$$ \left\{ \begin{array}{l} -i(\kappa_K \eta_K^{-1} p_K, q_K)_K - (\mathbf{u}_K, \nabla q_K)_K + \sum_{F \in \mathcal{F}_K} \langle (\mathbf{n}_{K,F} \cdot \mathbf{u}_K + \mathcal{A}_F^{-1} p_K), q_K \rangle_F \\ = \sum_{F \in \mathcal{F}_K} \langle \mathcal{A}_F^{-1} \hat{p}_F, q_K \rangle_F, \\ -i(\kappa_K \eta_K \mathbf{u}_K, \mathbf{v}_K)_K - (p_K, \nabla \cdot \mathbf{v}_K)_K = - \sum_{F \in \mathcal{F}_K} \langle \hat{p}_F, \mathbf{n}_{K,F} \cdot \mathbf{v}_K \rangle_F, \end{array} \right. $$
for a given surface datum $\hat{p}_F \in \mathcal{P}_p(F)$ for all $F \in \mathcal{F}_K$.

This local problem can be interpreted as a discretized Helmholtz problem defined on $K$ with a non-homogeneous Dirichlet boundary condition on $\partial K$. The discrete problem is well-posed without any condition, as stated and proved below. It is worth pointing out that at the continuous level, Helmholtz problems set in $K$ with Dirichlet boundary conditions might not be well-posed, due to possible resonance frequencies. Well-posedness at the discrete level follows from the fact that the solutions for the resonances do not belong to $\mathcal{P}_p(K)$ and $\mathcal{P}_p(K)$, as it can be seen at the end of the proof.

**Theorem 3.4** (Well-posedness of the local discrete problem). *Problem 3.3 is well-posed.*

*Proof.* We have to prove that, if $\hat{p}_F = 0$ for all $F \in \mathcal{F}_K$, the unique solution of Problem 3.3 is $p_K = 0$ and $\mathbf{u}_K = \mathbf{0}$. For brevity, the subscripts $K$ and $F$ are omitted for the local fields, the test functions, the outgoing unit normal and the coefficients. Taking both equations of Problem 3.3 with $q = p$ and $\mathbf{v} = \mathbf{u}$ gives

$$-i(\kappa\eta^{-1}p, p)_K - (\mathbf{u}, \nabla p)_K + \langle (\mathbf{n} \cdot \mathbf{u}) + \mathcal{A}^{-1}p, p \rangle_{\partial K} = 0, \\ -i(\kappa\eta\mathbf{u}, \mathbf{u})_K - (p, \nabla \cdot \mathbf{u})_K = 0.$$

Integrating by parts in both equations and taking the complex conjugate lead to

$$\begin{aligned} &i(\kappa\eta^{-1}p, p)_K + (p, \nabla \cdot \mathbf{u})_K + \langle p, \mathcal{A}^{-1}p \rangle_{\partial K} = 0, \\ &i(\kappa\eta\mathbf{u}, \mathbf{u})_K + (\mathbf{u}, \nabla p)_K - \langle \mathbf{n} \cdot \mathbf{u}, p \rangle_{\partial K} = 0. \end{aligned}$$

Adding the four previous equations yields $\langle \mathcal{A}^{-1}p, p \rangle_{\partial K} = 0$. Combined with Assumption 2.2, this implies that $p = 0$ on $\partial K$. By using this result in Problem 3.3, one has

$$\left\{ \begin{array}{l} -i(\kappa\eta^{-1}p, q)_K + (\nabla \cdot \mathbf{u}, q)_K = 0, \\ -i(\kappa\eta\mathbf{u}, \mathbf{v})_K + (\nabla p, \mathbf{v})_K = 0, \end{array} \right.$$

for all $(q, \mathbf{v}) \in \mathcal{P}_p(K) \times \mathcal{P}_p(K)$. We conclude that

$$\begin{aligned} &-i\kappa\eta^{-1}p + \nabla \cdot \mathbf{u} = 0, \\ &-i\kappa\eta\mathbf{u} + \nabla p = 0, \end{aligned}$$

in a strong sense. Because there is no non-zero polynomial solution to the previous equations, this yields the result. $\square$

*Remark 3.5* (Case with upwind fluxes). The well-posedness of the local problem in the case of upwind fluxes can be proved in a very similar way by following the same steps, see e.g. [12, 14].

# 4 Hybridization with transmission variables

Following the approach proposed in [34], the hybridization is performed by taking the incoming transmission variable as hybrid variable. As with standard hybridization, this approach does not alter the accuracy of the numerical solution when a direct solver is used, but it changes the properties of the linear system that needs to be solved.
## 4.1 CHDG formulations

An additional variable, denoted $g_h^{\ominus}$, corresponding to the incoming transmission variable is introduced at each face of each element. This variable belongs to the space

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

Each field of this space has one value associated to each face of each element. Therefore, there are two values per interior face of the mesh, and one value per boundary face.

For the case with upwind fluxes, the hybrid variable is defined by equation (2.3), and the numerical fluxes are defined by equation (2.5). For the case with symmetric fluxes, equations (2.6) and (2.2) are used instead. The following formulation is obtained in both cases.

**Problem 4.1** (CHDG formulation with upwind/symmetric fluxes). Find $(p_h, u_h, g_h^{\ominus}) \in V_h \times V_h \times G_h$ such that, for all $(q_h, v_h, \xi_h) \in V_h \times V_h \times G_h$,

$$\begin{cases} -i(\kappa\eta^{-1}p_h, q_h)_{\mathcal{T}_h} - (u_h, \nabla q_h)_{\mathcal{T}_h} + \langle \mathbf{n} \cdot \hat{\mathbf{u}}(g^{\oplus}(p_h, u_h), g_h^{\ominus}), q_h \rangle_{\partial \mathcal{T}_h} = 0, \\ -i(\kappa\eta\mathbf{u}_h, v_h)_{\mathcal{T}_h} - (p_h, \nabla \cdot v_h)_{\mathcal{T}_h} + \langle \hat{p}(g^{\oplus}(p_h, u_h), g_h^{\ominus}), \mathbf{n} \cdot v_h \rangle_{\partial \mathcal{T}_h} = 0, \end{cases}$$

and

$$\langle g_h^{\ominus} - \Pi(g^{\oplus}(p_h, u_h)), \xi_h \rangle_{\partial \mathcal{T}_h} = \langle b, \xi_h \rangle_{\partial \mathcal{T}_h},$$

with $g^{\oplus}(p_h, u_h)|_{K,F} := p_K + \eta_K \mathbf{n}_{K,F} \cdot \mathbf{u}_K$ in the upwind case and $g^{\oplus}(p_h, u_h)|_{K,F} := A_F^{-1} p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K$ in the symmetric case.

To simplify the presentation, we have introduced the global exchange operator $\Pi: G_h \to G_h$ and the global right-hand side $b$, whose definitions depend on the choice of numerical fluxes. For each face $F$ of each element $K$, and for any $g^{\oplus} \in G_h$, they are defined as

$$\Pi(g^{\oplus})|_{K,F} = \left\{ \begin{array}{ll} g_{K',F}^{\oplus}, & \text{if } F \not\subset \partial\Omega, \\ -g_{K,F}^{\oplus}, & \text{if } F \subset \Gamma_D, \\ g_{K,F}^{\oplus}, & \text{if } F \subset \Gamma_N, \\ 0, & \text{if } F \subset \Gamma_R, \end{array} \right. \quad \text{and} \quad b|_{K,F} = \left\{ \begin{array}{ll} 0, & \text{if } F \not\subset \partial\Omega, \\ 2s_D, & \text{if } F \subset \Gamma_D, \\ -2\eta_K s_N, & \text{if } F \subset \Gamma_N, \\ s_R, & \text{if } F \subset \Gamma_R, \end{array} \right.$$

for the upwind fluxes, and as

$$\Pi(g^{\oplus})|_{K,F} = \left\{ \begin{array}{ll} g_{K',F}^{\oplus}, & \text{if } F \not\subset \partial\Omega, \\ -g_{K,F}^{\oplus}, & \text{if } F \subset \Gamma_D, \\ g_{K,F}^{\oplus}, & \text{if } F \subset \Gamma_N, \\ \mathcal{B}_{K,F,+}^{-1} \mathcal{B}_{K,F,-} g_{K,F}^{\oplus}, & \text{if } F \subset \Gamma_R, \end{array} \right. \quad \text{and} \quad b|_{K,F} = \left\{ \begin{array}{ll} 0, & \text{if } F \not\subset \partial\Omega, \\ 2\mathcal{A}_F^{-1}s_D, & \text{if } F \subset \Gamma_D, \\ -2s_N, & \text{if } F \subset \Gamma_N, \\ 2\eta_K^{-1}\mathcal{B}_{K,F,+}^{-1}s_R, & \text{if } F \subset \Gamma_R, \end{array} \right.$$

for the symmetric fluxes.

## 4.2 Local element-wise discrete problems

The physical fields $p_h$ and $u_h$ are eliminated by solving local element-wise problems where the incoming transmission variable $g_h^{\ominus}$ is considered as a given datum. In the symmetric case, for each
element K, the local problem reads:

**Problem 4.2.** Find $(p_K, \mathbf{u}_K) \in \mathcal{P}_p(K) \times \mathcal{P}_p(K)$ such that, for all $(q_K, \mathbf{v}_K) \in \mathcal{P}_p(K) \times \mathcal{P}_p(K)$,

$$ 
\left\{ 
\begin{aligned}
& -\iota(\kappa_K \eta_K^{-1} p_K, q_K)_K - (\mathbf{u}_K, \nabla q_K)_K + \sum_{F \in \mathcal{F}_K} \frac{1}{2} \langle (\mathcal{A}_F^{-1} p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K), q_K \rangle_F \\
& = \sum_{F \in \mathcal{F}_K} \frac{1}{2} \langle g_{K,F}^{\ominus}, q_K \rangle_F, \\
& -\iota(\kappa_K \eta_K \mathbf{u}_K, \mathbf{v}_K)_K - (p_K, \nabla \cdot \mathbf{v}_K)_K + \sum_{F \in \mathcal{F}_K} \frac{1}{2} \langle (p_K + \mathcal{A}_F(\mathbf{n}_{K,F} \cdot \mathbf{u}_K)), \mathbf{n}_{K,F} \cdot \mathbf{v}_K \rangle_F \\
& = -\sum_{F \in \mathcal{F}_K} \frac{1}{2} \langle \mathcal{A}_F g_{K,F}^{\ominus}, \mathbf{n}_{K,F} \cdot \mathbf{v}_K \rangle_F,
\end{aligned}
\right. 
$$

for a given surface datum $g_{K,F}^{\ominus} \in \mathcal{P}_p(F)$ for all $F \in \mathcal{F}_K$.

This local discrete problem corresponds to a discretized Helmholtz problem defined on K with a non-homogeneous Robin boundary condition on $\partial K$. This problem is well-posed without any condition, as stated and proved hereafter. In contrast to the standard HDG framework, these local problems are also well-posed at the continuous level.

**Theorem 4.3 (Well-posedness of the local discrete problem).** *Problem 4.2 is well-posed.*

*Proof.* We simply have to prove that, if $g_{K,F}^{\ominus} = 0$ for all $F \in \mathcal{F}_K$, the unique solution of Problem 4.2 is $p_K = 0$ and $\mathbf{u}_K = \mathbf{0}$. For brevity, the subscripts K and F are omitted for the local fields, the test functions, the outgoing unit normal and the coefficients. Taking both equations of Problem 4.2 with $q = p$ and $\mathbf{v} = \mathbf{u}$ gives

$$ 
\begin{aligned}
& -\iota(\kappa \eta^{-1} p, p)_K - (\mathbf{u}, \nabla p)_K + \sum_{F \in \mathcal{F}_K} \frac{1}{2} \langle (\mathcal{A}^{-1} p + \mathbf{n} \cdot \mathbf{u}), p \rangle_F = 0, \\
& -\iota(\kappa \eta \mathbf{u}, \mathbf{u})_K - (p, \nabla \cdot \mathbf{u})_K + \sum_{F \in \mathcal{F}_K} \frac{1}{2} \langle (p + \mathcal{A}(\mathbf{n} \cdot \mathbf{u})), \mathbf{n} \cdot \mathbf{u} \rangle_F = 0.
\end{aligned}
$$

Integrating by parts in both equations and taking the complex conjugate lead to

$$ 
\begin{aligned}
& \iota(\kappa \eta^{-1} p, p)_K + (p, \nabla \cdot \mathbf{u})_K + \sum_{F \in \mathcal{F}_K} \frac{1}{2} \langle p, (\mathcal{A}^{-1} p - \mathbf{n} \cdot \mathbf{u}) \rangle_F = 0, \\
& \iota(\kappa \eta \mathbf{u}, \mathbf{u})_K + (\mathbf{u}, \nabla p)_K + \sum_{F \in \mathcal{F}_K} \frac{1}{2} \langle \mathbf{n} \cdot \mathbf{u}, (-p + \mathcal{A}(\mathbf{n} \cdot \mathbf{u})) \rangle_F = 0.
\end{aligned}
$$

Adding the four previous equations yields

$$ \sum_{F \in \mathcal{F}_K} \left( \frac{1}{2} \langle \mathcal{A}^{-1} p, p \rangle_F + \frac{1}{2} \langle \mathcal{A}(\mathbf{n} \cdot \mathbf{u}), \mathbf{n} \cdot \mathbf{u} \rangle_F + \frac{1}{2} \langle p, \mathcal{A}^{-1} p \rangle_F + \frac{1}{2} \langle \mathbf{n} \cdot \mathbf{u}, \mathcal{A}(\mathbf{n} \cdot \mathbf{u}) \rangle_F \right) = 0. $$

Since $\mathcal{A}$ and $\mathcal{A}^{-1}$ are self-adjoint, we have

$$ \langle \mathcal{A}^{-1} p, p \rangle_{\partial K} + \langle \mathcal{A}(\mathbf{n} \cdot \mathbf{u}), \mathbf{n} \cdot \mathbf{u} \rangle_{\partial K} = 0, $$

which gives $p = 0$ and $\mathbf{n} \cdot \mathbf{u} = 0$ on $\partial K$, thanks to the positivity of $\mathcal{A}$ and $\mathcal{A}^{-1}$. By using these boundary conditions in Problem 4.2, we have that the fields should be a solution of the strong
problem. Because there is no solution with both homogeneous Neumann and Dirichlet boundary conditions, this yields the result. □

Remark 4.4 (Case with upwind fluxes). It is possible to write a local element-wise discrete problem similar to Problem 4.2 in the case of upwind fluxes and prove its well-posedness by following the steps of the proof of Theorem 4.3. The task can be achieved with a similar reasoning. For brevity, the definition of the problem and the proof of its well-posedness are omitted.

## 4.3 Abstract form of the hybridized system

The hybridized CHDG problem can be written in a convenient abstract form by introducing the global scattering operator S: $G_h \to G_h$ defined such that, for each face F of each element K,

$$S(g_h^{\ominus})|_{K,F} := A_F^{-1} p_K(g_h^{\ominus}) + n_{K,F} \cdot u_K(g_h^{\ominus}), \quad (4.1)$$

where $(p_K(g_h^{\ominus}), u_K(g_h^{\ominus}))$ is the solution of Problem 4.2 with the incoming transmission variable $(g_{K,F}^{\ominus})_{F \in \mathcal{F}_K}$ contained in $g_h^{\ominus}$ as a given surface datum. This operator can be interpreted as an “incoming transmission variable to outgoing transmission variable” operator.

By using the operator S, Problem 4.2 is rewritten as:

**Problem 4.5.** Find $g_h^{\ominus} \in G_h$ such that, for all $\xi_h \in G_h$,

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

We introduce the global projected right-hand side $b_h := P_h b \in G_h$, where $P_h: L^2(\partial \mathcal{T}_h) \to G_h$ is the projection operator defined such that $\langle P_h b, \xi_h \rangle_{\partial \mathcal{T}_h} = \langle b, \xi_h \rangle_{\partial \mathcal{T}_h}$ for all $\xi_h \in G_h$. Problem 4.5 can then be rewritten as:

**Problem 4.6.** Find $g_h^{\ominus} \in G_h$ such that

$$(I - \Pi S)g_h^{\ominus} = b_h,$$

where I is the identity operator on $G_h$. Problem 4.5 and Problem 4.6 are equivalent to Problem 4.1 because the element-wise local problems are well-posed.

## 4.4 Strict contraction of the operator ΠS (case with symmetric fluxes)

An interesting property of the CHDG formulation when using symmetric fluxes satisfying Assumption 2.2 is that the operator ΠS is a strict contraction. As a consequence, Problem 4.6 is always well-posed, and it can be solved with the fixed-point iteration without relaxation.

Properties of S and Π are proved by using norms associated to $G_h$ and $\bigoplus_{F \in \mathcal{F}_K} \mathcal{P}_P(F)$ defined as

$$\|g_h^{\ominus}\| := \sqrt{\sum_{K \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_K} \|g_{K,F}^{\ominus}\|_F^2} \quad \text{and} \quad \|u\|_{\partial K} := \sqrt{\sum_{F \in \mathcal{F}_K} \|u\|_F^2},$$

where $\| \cdot \|_F$ is the norm of $L^2(F)$ induced by $\mathcal{A}_F$, namely

$$\|u\|_F := \sqrt{\langle \mathcal{A}_F u, u \rangle_F}.$$
**Lemma 4.7.** (i) The solution of Problem 4.2 verifies

$$ \sum_{F \in \mathcal{F}_K} \| \mathcal{A}_F^{-1} p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K \|_F^2 + \sum_{F \in \mathcal{F}_K} \| \mathcal{A}_F^{-1} p_K - \mathbf{n}_{K,F} \cdot \mathbf{u}_K - g_{K,F}^{\ominus} \|_F^2 = \sum_{F \in \mathcal{F}_K} \| g_{K,F}^{\ominus} \|_F^2. \quad (4.2) $$

(ii) The second term on the left-hand side of (4.2) vanishes if and only if $g_{K,F}^{\ominus} = 0$.

*Proof.* For brevity, the subscripts K and F are omitted for the local fields, the test functions, the unit outgoing normal, the surface data and the coefficients.

(i) Taking both equations of Problem 4.2 with $q = p$ and $v = u$ gives

$$ \begin{aligned} -i(\kappa\eta^{-1}p, p)_K - (\mathbf{u}, \nabla p)_K + \frac{1}{2}\langle \mathcal{A}(\mathcal{A}^{-1}p + \mathbf{n} \cdot \mathbf{u}), \mathcal{A}^{-1}p \rangle_{\partial K} &= \frac{1}{2}\langle \mathcal{A}g^{\ominus}, \mathcal{A}^{-1}p \rangle_{\partial K}, \\ -i(\kappa\eta\mathbf{u}, \mathbf{u})_K - (p, \nabla \cdot \mathbf{u})_K + \frac{1}{2}\langle \mathcal{A}(\mathcal{A}^{-1}p + \mathbf{n} \cdot \mathbf{u}), \mathbf{n} \cdot \mathbf{u} \rangle_{\partial K} &= -\frac{1}{2}\langle \mathcal{A}g^{\ominus}, \mathbf{n} \cdot \mathbf{u} \rangle_{\partial K}. \end{aligned} $$

Integrating by parts in both equations and taking the complex conjugate lead to

$$ \begin{aligned} i(\kappa\eta^{-1}p, p)_K + (p, \nabla \cdot \mathbf{u})_K + \frac{1}{2}\langle \mathcal{A}^{-1}p, \mathcal{A}(\mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u}) \rangle_{\partial K} &= \frac{1}{2}\langle \mathcal{A}^{-1}p, \mathcal{A}g^{\ominus} \rangle_{\partial K}, \\ i(\kappa\eta\mathbf{u}, \mathbf{u})_K + (\mathbf{u}, \nabla p)_K - \frac{1}{2}\langle \mathbf{n} \cdot \mathbf{u}, \mathcal{A}(\mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u}) \rangle_{\partial K} &= -\frac{1}{2}\langle \mathbf{n} \cdot \mathbf{u}, \mathcal{A}g^{\ominus} \rangle_{\partial K}. \end{aligned} $$

Adding the four previous equations and multiplying by two yield

$$ \begin{aligned} \langle \mathcal{A}(\mathcal{A}^{-1}p + \mathbf{n} \cdot \mathbf{u}), \mathcal{A}^{-1}p + \mathbf{n} \cdot \mathbf{u} \rangle_{\partial K} &+ \langle \mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u}, \mathcal{A}(\mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u}) \rangle_{\partial K} \\ &= \langle \mathcal{A}g^{\ominus}, \mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u} \rangle_{\partial K} + \langle (\mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u}), \mathcal{A}g^{\ominus} \rangle_{\partial K}. \end{aligned} \quad (4.3) $$

Using the following identities for rewriting the right-hand side,

$$ \begin{aligned} \langle \mathcal{A}g^{\ominus}, \mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u} \rangle_{\partial K} &+ \langle (\mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u}), \mathcal{A}g^{\ominus} \rangle_{\partial K} \\ &= \langle \mathcal{A}(\mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u}), \mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u} \rangle_{\partial K} + \langle \mathcal{A}g^{\ominus}, g^{\ominus} \rangle_{\partial K} \\ &\quad - \langle \mathcal{A}(\mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u} - g^{\ominus}), \mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u} - g^{\ominus} \rangle_{\partial K} \\ &= \|\mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u}\|_{\partial K}^2 - \|\mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u} - g^{\ominus}\|_{\partial K}^2 + \|g^{\ominus}\|_{\partial K}^2, \end{aligned} $$

equation (4.3) becomes

$$ \begin{aligned} \| \mathcal{A}^{-1}p + \mathbf{n} \cdot \mathbf{u} \|_{\partial K}^2 + \| \mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u} \|_{\partial K}^2 \\ &= \| \mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u} \|_{\partial K}^2 - \| \mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u} - g^{\ominus} \|_{\partial K}^2 + \| g^{\ominus} \|_{\partial K}^2, \end{aligned} $$

which gives the result (4.2).

(ii) If the second term on the left-hand side of (4.2) vanishes, then $g^{\ominus} = \mathcal{A}^{-1}p - \mathbf{n} \cdot \mathbf{u}$ on $\partial K$. Using this relation in Problem 4.2, we see that $p$ and $u$ must satisfy

$$ \begin{aligned} -i(\kappa\eta^{-1}p, q)_K - (\mathbf{u}, \nabla q)_K + \langle \mathbf{n} \cdot \mathbf{u}, q \rangle_{\partial K} &= 0, \\ -i(\kappa\eta\mathbf{u}, \mathbf{v})_K - (p, \nabla \cdot \mathbf{v})_K + \langle p, \mathbf{n} \cdot \mathbf{v} \rangle_{\partial K} &= 0, \end{aligned} $$

for all $q \in \mathcal{P}_p(K)$ and $v \in \mathcal{P}_p(K)$, and integration by parts shows that $p$ and $u$ solve the Helmholtz equation in strong form. Because there is no non-zero polynomial solution to the previous equations, meaning that $p = 0$ and $u = 0$, and then $g^{\ominus} = 0$, this yields the result. The
converse statement is direct, because the local problem is well-posed.

☐

**Theorem 4.8.** The scattering operator S is a strict contraction, i.e.

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

*Proof.* Let $g_h^{\ominus} \in G_h \setminus \{0\}$. By Lemma 4.7, one has

$$ \sum_{F \in \mathcal{F}_K} \|A_F^{-1} p_K + n_{K,F} \cdot u_K\|_F^2 < \sum_{F \in \mathcal{F}_K} \|g_{K,F}^{\ominus}\|_F^2. $$

Note that the equality cannot happen because $g_h^{\ominus} \neq 0$. Then, by using the definition of S in equation (4.1), one has

$$ \sum_{F \in \mathcal{F}_K} \|S(g_h^{\ominus})|_{K,F}\|_F^2 < \sum_{F \in \mathcal{F}_K} \|g_{K,F}^{\ominus}\|_F^2. $$

Summing this estimate over all $K \in \mathcal{T}_h$ gives the result.

☐

**Theorem 4.9.** (i) If $\Gamma_R = \emptyset$, $\Pi$ is an involution, i.e. $\Pi^2 = I$, and an isometry, i.e.

$$ \|\Pi(g_h^{\ominus})\| = \|g_h^{\ominus}\|, \quad \forall g_h^{\ominus} \in G_h. $$

(ii) If $\Gamma_R \neq \emptyset$, the exchange operator $\Pi$ is a contraction, i.e.

$$ \|\Pi(g_h^{\ominus})\| \le \|g_h^{\ominus}\|, \quad \forall g_h^{\ominus} \in G_h. $$

The inequality is strict for all $g_h^{\ominus} \in G_h$ such that there is at least one face $F \subset \Gamma_R$ where $g_{K,F}^{\ominus}$ is non-zero.

*Proof.* (i) The result is a straightforward consequence of the definition of $\Pi$ for interior faces and for boundary faces belonging to $\Gamma_D$ or $\Gamma_N$.

(ii) For every boundary face $F \subset \Gamma_R$ belonging to an element $K$, we have to prove that

$$ \|(1 + \eta_K^{-1} A_F)^{-1} (1 - \eta_K^{-1} A_F) g\|_F \le \|g\|_F, \quad \forall g \in \mathcal{P}_p(F). $$

For any given $g \in \mathcal{P}_p(F)$, the inequality holds if and only if

$$ \|(1 - \eta_K^{-1} A_F) \xi\|_F \le \|(1 + \eta_K^{-1} A_F) \xi\|_F, $$

with $\xi := (1 + \eta_K^{-1} A_F)^{-1} g$. Because

$$ \|(1 + \eta_K^{-1} A_F) \xi\|_F^2 = \|\xi\|_F^2 + \|(\eta_K^{-1} A_F) \xi\|_F^2 + 2 \langle \eta_K^{-1} A_F^2 \xi, \xi \rangle_F, $$

$$ \|(1 - \eta_K^{-1} A_F) \xi\|_F^2 = \|\xi\|_F^2 + \|(\eta_K^{-1} A_F) \xi\|_F^2 - 2 \langle \eta_K^{-1} A_F^2 \xi, \xi \rangle_F, $$

and $\langle \eta_K^{-1} A_F^2 \xi, \xi \rangle_F \ge 0$, the result holds true. In addition, if $g \neq 0$, then $\xi \neq 0$ and the inequality is strict.

☐

**Corollary 4.10.** The operator $\Pi S$ is a strict contraction, i.e.

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

*Proof.* Corollary 4.10 is a direct consequence of Theorems 4.8 and 4.9.

☐
Remark 4.11 (Interpretation). The strict contraction of S is a numerical effect: we have used the fact that the Helmholtz equation has no homogeneous polynomial solution in the proof of Lemma 4.7. The version of equation (4.2) for the continuous case corresponds to a conservation of energy. By contrast, the properties of Π are related to the boundary conditions, and would be preserved at the continuous level. This operator is contracting only with a boundary condition that does not preserve energy, i.e. the Robin boundary condition. More precisely, it is strictly contracting only for fields that are different from zero on Γ<sub>R</sub>.

# 5 Numerical comparison of the methods

In this section, iterative solution procedures for the HDG and CHDG methods are studied and compared by using two simple benchmarks and a more realistic application. For the first two benchmarks, we consider configurations with homogeneous and heterogeneous media in Sections 5.2 and 5.3, respectively. For the last benchmark (Section 5.4), the complex heterogeneous medium illustrates a realistic subsurface medium for seismic wave propagation.

## 5.1 Benchmark and numerical setting

The numerical simulations are performed with a dedicated MATLAB code already used in Ref. [34]. The mesh generation and the visualization are performed with gmsh [25]. In all the cases, third-order polynomial Lobatto basis functions (p = 3) are used, and h is the element size specified in gmsh. For the symmetric fluxes, we take $\mu_F = \sqrt{\eta_K \eta_{K'}}$ and $\kappa_F = \sqrt{\kappa_K \kappa_{K'}}$ for an interior face F shared by two elements K and K'.

**Benchmark 1 (Plane wave).** The first benchmark represents the reflection of an incident plane wave at the interface between two media. The problem is defined on a square domain $\Omega = (0, 1) \times (0, 1)$ partitioned into two rectangular regions $\Omega_1 = (0, 1/2) \times (0, 1)$ and $\Omega_2 = (1/2, 1) \times (0, 1)$ corresponding to the two media. The wavenumber and impedance, which are constant in each region, are denoted by $\kappa_1$ and $\eta_1$ in $\Omega_1$, and by $\kappa_2$ and $\eta_2$ in $\Omega_2$. The exact solution in region $\Omega_1$ is the sum of the incident and reflected waves, while in region $\Omega_2$ it is the transmitted wave. The solution satisfies the continuity of $p$ and $n \cdot u$ across the interface ($x = 1/2$) between the two regions $\Omega_1$ and $\Omega_2$. It can be written as

$$p_{\text{ref}}(x, y) = \begin{cases} e^{i\kappa_1}(x \cos \theta_I + y \sin \theta_I) + Re^{i\kappa_1}(-x \cos \theta_I + y \sin \theta_I), & \text{in } \Omega_1, \\ Te^{i\kappa_2}(x \cos \theta_T + y \sin \theta_T), & \text{in } \Omega_2, \end{cases}$$

with the reflection coefficient R and transmission coefficient T given by

$$R = \frac{\eta_2 \cos \theta_I - \eta_1 \cos \theta_T}{\eta_1 \cos \theta_T + \eta_2 \cos \theta_I} e^{i\kappa_1 \cos \theta_I} \quad \text{and} \quad T = \frac{2\eta_2 \cos \theta_I}{\eta_1 \cos \theta_T + \eta_2 \cos \theta_I} e^{i(\kappa_1 \cos \theta_I - \kappa_2 \cos \theta_T)/2}.$$

The angle $\theta_I$ of the incident wave is set to $\theta_I = \pi/4$ in the following. The angle of the transmitted wave is $\theta_T = \arcsin(\frac{\kappa_1}{\kappa_2} \sin \theta_I)$. The reference solution is shown for cases with homogeneous and heterogeneous media in Figures 1a and 1b, respectively. In the numerical model, these solutions are enforced by using non-homogeneous Robin boundary conditions on the boundary of $\Omega$ with boundary data specified by the reference solution, i.e. $s_R = p_{\text{ref}} - \eta n \cdot u_{\text{ref}}$ on $\partial\Omega = \Gamma_R$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Real part of the reference solutions $p_{\text{ref}}$ for both benchmark problems with the default parameters corresponding to the first cases of each benchmark in Tables 1 and 2.</figcaption></figure>

**Benchmark 2 (Cavity).** We consider a circular domain $\Omega$ with radius 1/2 partitioned into the circular region $\Omega_1$ with radius 1/4 and the annulus $\Omega_2$ with the radial coordinate $r \in (1/4, 1/2)$, all centered at the origin. The physical coefficients are constant in each region, with the same notation as in the first benchmark. A homogeneous Dirichlet boundary condition is prescribed on the boundary of the domain, i.e. $\partial\Omega = \Gamma_D$ and $s_D = 0$. The constant volume source term $f = -1/(i\kappa\eta)$ is introduced on the right-hand side of the first equation in (2.1).

If the physical coefficients do not correspond to a resonance mode, the solution of this benchmark is unique and real. In this case, the reference solution, which depends only on the radial coordinate $r$, is given by

$$p_{\text{ref},j}(r) = A_j J_0(\kappa_j r) + B_j Y_0(\kappa_j r) - \kappa_j^{-2}, \quad \text{in } \Omega_j, \quad j = 1, 2,$$

where $A_j$ and $B_j$ are constant coefficients, $J_0$ is the Bessel function of the first kind of order 0 and $Y_0$ is the Bessel function of the second kind of order 0. To preclude singular solutions at the origin, one has to set $B_1 = 0$. The other three coefficients can be determined by enforcing the continuity of $p$ and $n \cdot u$ across the interface and the homogeneous Dirichlet boundary condition. This leads to the following linear system:

$$ \begin{bmatrix} 0 & J_0(\kappa_2 R_2) & Y_0(\kappa_2 R_2) \\ J_0(\kappa_1 R_1) & -J_0(\kappa_2 R_1) & -Y_0(\kappa_2 R_1) \\ J_1(\kappa_1 R_1)/\eta_1 & -J_1(\kappa_2 R_1)/\eta_2 & -Y_1(\kappa_2 R_1)/\eta_2 \end{bmatrix} \begin{bmatrix} A_1 \\ A_2 \\ B_2 \end{bmatrix} = \begin{bmatrix} \kappa_2^{-2} \\ \kappa_1^{-2} - \kappa_2^{-2} \\ 0 \end{bmatrix}, $$

with $R_1 = 1/4$ and $R_2 = 1/2$. Resonance phenomena occur whenever the values of $\kappa_1 R_1$ and $\kappa_2 R_2$ are such that the determinant of the matrix associated to the linear system is zero, namely the matrix is not invertible. The reference solutions are shown for cases with homogeneous and heterogeneous media in Figures 1c and 1d, respectively.

**Iterative solvers.** Three standard iterations are considered: the fixed-point iteration (only for the CHDG system), the CGNR (conjugate gradient normal residual) iteration, and the GMRES (generalized minimal residual) iteration without restart. For a given system $Ag = b$, the second one corresponds to the conjugate gradient iteration applied to the normal system $A^*Ag = A^*b$, see e.g. [40, 42].

Similarly to the approach used in [34], we have used a symmetric preconditioning with the mass
matrix **M** associated to the faces of the elements. Denoting the Cholesky factorization of the mass matrix with **M** = **L****L**$^T$, this leads to the system $\tilde{\mathbf{A}}\tilde{\mathbf{g}} = \tilde{\mathbf{b}}$ with $\tilde{\mathbf{A}} := \mathbf{L}^{-1}\mathbf{A}\mathbf{L}^{-1}$, $\tilde{\mathbf{g}} := \mathbf{L}^T\mathbf{g}$ and $\tilde{\mathbf{b}} := \mathbf{L}^{-1}\mathbf{b}$. The preconditioned system corresponds to the one which would be used if orthonormal basis functions would be used on the face, see [34]. Then, the 2-norm of an algebraic vector $\tilde{\mathbf{g}}$ is equal to the $L^2$-norm of the corresponding field, i.e. $\|\tilde{\mathbf{g}}\|_2 = \sqrt{\tilde{\mathbf{g}}^*\tilde{\mathbf{g}}} = \sqrt{\tilde{\mathbf{g}}^*\mathbf{M}\tilde{\mathbf{g}}} = \|g_h\|$.

**Numerical error.** In this study, we consider a relative numerical error based on the energy norm of the physical fields. It is defined as

$$ \text{relative error} := \frac{\|p_h - p_{\text{ref}}, \mathbf{u}_h - \mathbf{u}_{\text{ref}}\|_E}{\|p_{\text{ref}}, \mathbf{u}_{\text{ref}}\|_E}, \quad \text{with } \|p, \mathbf{u}\|_E^2 := \sum_{K \in \mathcal{T}_h} \frac{\|p_K\|_{L^2(K)}^2}{2\rho_K c_K^2} + \frac{1}{2}\rho_K \|\mathbf{u}_K\|_{L^2(K)}^2, $$

where $(p_h, \mathbf{u}_h)$ is the numerical solution, $(p_{\text{ref}}, \mathbf{u}_{\text{ref}})$ is the reference solution, and $\|\cdot\|_E$ denotes the energy norm.

## 5.2 Comparison for homogeneous media

The benchmark problems are solved for constant physical parameters corresponding to two levels of difficulty: low- and high-frequency cases for the plane-wave benchmark, and wavenumbers close and very close to the resonance mode $\kappa_3 \approx 17.3075$ for the cavity benchmark. For the circular cavity with homogeneous medium, a resonance corresponds to $\kappa_j = 2x_j$, where $\{x_j\}_{j \in \mathbb{N}}$ are the zeros of $J_0(x)$. In all the cases, the mesh sizes are chosen to achieve a relative error close to $10^{-2}$ when a direct solver is employed. The parameters are given in Table 1.

<table><caption>(a) Benchmark 1 (plane wave)</caption><thead><tr><th>Case</th><th>&kappa;</th><th>h</th><th>Numerical flux</th><th>&rho;(&Pi;S)</th><th>Relative error</th></tr></thead><tbody><tr><td rowspan="2">1</td><td rowspan="2">15&pi;</td><td rowspan="2">1/16</td><td>Sym-0</td><td>1 &minus; 4.18 &middot; 10<sup>-3</sup></td><td>1.44 &middot; 10<sup>-2</sup></td></tr><tr><td>Sym-2</td><td>1 &minus; 7.97 &middot; 10<sup>-3</sup></td><td>1.70 &middot; 10<sup>-2</sup></td></tr><tr><td rowspan="2">2</td><td rowspan="2">30&pi;</td><td rowspan="2">1/34</td><td>Sym-0</td><td>1 &minus; 1.88 &middot; 10<sup>-3</sup></td><td>1.37 &middot; 10<sup>-2</sup></td></tr><tr><td>Sym-2</td><td>1 &minus; 3.48 &middot; 10<sup>-3</sup></td><td>1.59 &middot; 10<sup>-2</sup></td></tr></tbody></table>

<table><caption>(b) Benchmark 2 (cavity)</caption><thead><tr><th>Case</th><th>&kappa;</th><th>h</th><th>Numerical flux</th><th>&rho;(&Pi;S)</th><th>Relative error</th></tr></thead><tbody><tr><td rowspan="2">1</td><td rowspan="2">16.5</td><td rowspan="2">0.04</td><td>Sym-0</td><td>1 &minus; 1.22 &middot; 10<sup>-8</sup></td><td>1.07 &middot; 10<sup>-2</sup></td></tr><tr><td>Sym-2</td><td>1 &minus; 1.71 &middot; 10<sup>-6</sup></td><td>1.07 &middot; 10<sup>-2</sup></td></tr><tr><td rowspan="2">2</td><td rowspan="2">17</td><td rowspan="2">0.025</td><td>Sym-0</td><td>1 &minus; 1.15 &middot; 10<sup>-9</sup></td><td>1.16 &middot; 10<sup>-2</sup></td></tr><tr><td>Sym-2</td><td>1 &minus; 1.44 &middot; 10<sup>-7</sup></td><td>1.16 &middot; 10<sup>-2</sup></td></tr></tbody></table>

<caption>Table 1: Parameters for the benchmark problems with homogeneous media, including the spectral radius $\rho(\Pi S)$ of the matrix of the CHDG hybridized system, and the relative numerical error. In these cases, $\omega = \kappa$, $c = 1$, $\rho = 1$ and $\eta = 1$.</caption>

The history of relative error is plotted in Figure 2 during the fixed-point iterations (lines with marker $\times$, only for CHDG), the CGNR iterations (lines with marker $\circ$) and the GMRES iterations (lines with marker $\bullet$) applied to the HDG and CHDG hybridized systems. The 0th-order symmetric flux (Sym-0) is considered for both methods, and the 2nd-order symmetric flux (Sym-2) is used for CHDG. The relative errors obtained with a direct solver is indicated by the horizontal dashed lines for both numerical fluxes. Let us note that, since the impedance $\eta$ is constant, the upwind flux is identical to the 0th-order symmetric flux with $A_F = \eta$.

The following observations can be made:
<figure><img src="image_3.png" /></figure>

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

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

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

<figure><img src="image_7.png" /><figcaption>Figure 2: Results for the benchmark problems with homogeneous media. Error history with different iterative schemes and different DG methods. The dashed horizontal lines correspond to the relative numerical errors obtained with a direct solver for the different numerical fluxes.</figcaption></figure>
* The fixed-point iteration applied to the CHDG system converges in all the cases. The convergence is fast for the plane-wave benchmark, and very slow for the cavity benchmark. In both cases, the convergence is slower for the parameters corresponding to the more difficult configuration (i.e. high frequency case or close to a resonance frequency). These results are in accordance with the spectral radii given in Table 1: the convergence is slower for radii closer to 1.

* Comparing the numerical fluxes used in the CHDG method, we observe that the convergence of the fixed point iteration is comparable with 'Sym-0' and 'Sym-2'. The results are similar for the GMRES iteration. In contrast, the CGNR iteration with 'Sym-2' is slower than with 'Sym-0', sometimes significantly so.

* The convergence of the CGNR and GMRES iterations is slower with HDG than with CHDG for the plane-wave benchmark. The difference is more important for the high-frequency case. For the cavity benchmark, the convergence of GMRES is similar with both methods, while the convergence of CGNR is slower with HDG than with CHDG.

In a nutshell, CHDG with 'Sym-0' is one of the best solution in all the cases. The fixed-point iteration applied to the CHDG system converges, but its performance is strongly influenced by the physical parameters, and it becomes unusable for cavities. The second-order numerical flux does not accelerate the convergence of the iterative procedures.

## 5.3 Comparison for heterogeneous media

The benchmark problems are now tested with two sets of parameters. For both benchmarks, the impedance $\eta$ is constant in the first set, it is not constant in the second set, and the wavenumber $\kappa$ is not constant in all the sets. The parameters are given in Table 2. For all the cases, the mesh size is chosen to obtain a relative error close to $10^{-2}$ when a direct solver is employed.

<table><caption>(a) Benchmark 1 (plane wave)</caption><thead><tr><th>Case</th><th>ω</th><th>c</th><th>ρ</th><th>κ</th><th>η</th><th>h</th><th>Numerical flux</th><th>ρ(ΠS)</th><th>Relative error</th></tr></thead><tbody><tr><td rowspan="3">1</td><td rowspan="3">15π</td><td>1</td><td>1</td><td>15π</td><td>1</td><td>1/16</td><td>Upw</td><td>1 - 1.63 · 10<sup>-3</sup></td><td>1.01 · 10<sup>-2</sup></td></tr><tr><td>1/2</td><td>2</td><td>30π</td><td>1</td><td>1/34</td><td>Sym-0</td><td>1 - 1.63 · 10<sup>-3</sup></td><td>1.01 · 10<sup>-2</sup></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>Sym-2</td><td>1 - 4.32 · 10<sup>-3</sup></td><td>1.20 · 10<sup>-2</sup></td></tr><tr><td rowspan="3">2</td><td rowspan="3">15π</td><td>1</td><td>1</td><td>15π</td><td>1</td><td>1/16</td><td>Upw</td><td>1 + 2.38 · 10<sup>-2</sup></td><td>9.75 · 10<sup>-3</sup></td></tr><tr><td>1/2</td><td>1</td><td>30π</td><td>1/2</td><td>1/34</td><td>Sym-0</td><td>1 - 1.55 · 10<sup>-2</sup></td><td>9.74 · 10<sup>-3</sup></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>Sym-2</td><td>1 - 4.31 · 10<sup>-3</sup></td><td>1.16 · 10<sup>-2</sup></td></tr></tbody></table>

<table><caption>(b) Benchmark 2 (cavity)</caption><thead><tr><th>Case</th><th>ω</th><th>c</th><th>ρ</th><th>κ</th><th>η</th><th>h</th><th>Numerical flux</th><th>ρ(ΠS)</th><th>Relative error</th></tr></thead><tbody><tr><td rowspan="3">1</td><td rowspan="3">10π</td><td>1</td><td>1</td><td>10π</td><td>1</td><td>1/12</td><td>Upw</td><td>1 - 5.19 · 10<sup>-5</sup></td><td>1.57 · 10<sup>-2</sup></td></tr><tr><td>2/3</td><td>3/2</td><td>15π</td><td>1</td><td>1/16</td><td>Sym-0</td><td>1 - 5.19 · 10<sup>-5</sup></td><td>1.57 · 10<sup>-2</sup></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>Sym-2</td><td>1 - 5.28 · 10<sup>-5</sup></td><td>1.67 · 10<sup>-2</sup></td></tr><tr><td rowspan="3">2</td><td rowspan="3">10π</td><td>1</td><td>1</td><td>10π</td><td>1</td><td>1/12</td><td>Upw</td><td>1 + 3.45 · 10<sup>-4</sup></td><td>8.59 · 10<sup>-3</sup></td></tr><tr><td>2/3</td><td>1</td><td>15π</td><td>2/3</td><td>1/16</td><td>Sym-0</td><td>1 - 4.99 · 10<sup>-5</sup></td><td>8.61 · 10<sup>-3</sup></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>Sym-2</td><td>1 - 5.29 · 10<sup>-5</sup></td><td>1.01 · 10<sup>-2</sup></td></tr></tbody></table>

<caption>Table 2: Parameters of the benchmark problems with heterogeneous media, including spectral radius $ρ(ΠS)$ of the matrix of the CHDG hybridized system, and the relative numerical error.</caption>

The history of relative error is plotted in Figure 3 for different combinations of hybridized methods, numerical fluxes and iterative schemes. For CHDG, we have considered the upwind fluxes (Upw) and both zeroth- and second-order symmetric fluxes (Sym-0 and Sym-2). The numerical fluxes
<figure><img src="image_8.png" /></figure>

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

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

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

<figure><img src="image_12.png" /><figcaption>Figure 3: Results of the benchmark problems with heterogeneous media. Error history with different iterative schemes and different DG methods. The dashed horizontal lines correspond to the relative numerical errors obtained with a direct solver for the different numerical fluxes. In most graphs, the green lines (CHDG ‘Upw’) are hidden by the corresponding blue lines (CHDG ‘Sym-0’).</figcaption></figure>
‘Upw’ and ‘Sym-0’ are identical when the impedance is constant. Nonetheless, the results are very similar when the impedance is not constant. For HDG, both fluxes lead to the same results, even when the impedance is discontinuous.

We observe that the fixed-point iteration applied to the CHDG systems with the symmetric fluxes converges in all cases, while it fails to converge with the upwind fluxes when the impedance is discontinuous. The divergence can be seen in Figure 3b, and it also appears when increasing the number of iterations for the case corresponding to Figure 3d. These results are consistent with the theory and the spectral radii in Table 2: the contraction property of the operator IIS was proved in Corollary 4.10 only for symmetric fluxes, and the convergence occurs numerically for radii strictly less than 1. The convergence is quite fast for the plane-wave benchmark, while it slows down dramatically for the cavity benchmark.

Similarly to the cases with homogeneous media, when comparing the symmetric numerical fluxes used in the CHDG system, we observe that the convergence of CGNR and GMRES with ‘Sym-2’ is always either comparable or slower than with ‘Sym-0’ for both benchmarks. Once again, the convergence of the iterative schemes is much slower with HDG than with CHDG, regardless of the iterative procedure that is used.

In a nutshell, the CGNR and GMRES iterations combined with CHDG and either ‘Upw’ or ‘Sym-0’ are always the best options. The fixed-point iteration applied to CHDG with ‘Sym-0’ always converges, whereas it does not with HDG or CHDG with ‘Upw’.

## 5.4 Realistic application

To compare the methods on a more illustrative case, we consider the well-known Marmousi benchmark problem, which consists of a realistic underground structure commonly used to assess numerical methods, see e.g. [45]. The variations of the velocity $c(\mathbf{x})$ and the density $\rho(\mathbf{x})$ are shown in Figure 4.

An acoustic field is generated by a source point at the coordinates (4585m, -10m), near the top of the rectangular domain $\Omega = ]0m, 9192m[ \times ]-2094m, 0m[$. The mesh is generated with `gmsh` [25], and includes 60 761 elements. The element size is adjusted following the empirical rule $h \approx \lambda(\mathbf{x})/n_\lambda$, with the wavelength $\lambda(\mathbf{x}) = c(\mathbf{x})/f$, the frequency $f = 30$ Hz, the angular frequency $\omega = 2\pi f$, the mesh density $n_\lambda = 10/(p+1)$, and the polynomial degree $p = 3$. The physical coefficients are constant in each cell. The Dirichlet boundary condition $p = 0$ is applied at the top of the domain, and the absorbing boundary condition $p - \eta \mathbf{n} \cdot \mathbf{u} = 0$ is prescribed on the other sides of the domain. The mesh and the solution are shown in Figure 4.

The problem is solved using HDG and CHDG with the zeroth-order symmetric flux (Sym-0) and the upwind flux (Upw). We consider the fixed-point iteration (only for CHDG with ‘Sym-0’), the CGNR iteration, and the GMRES iterations (with a restart after every 10 iterations). The restart strategy, which is widely used in practice, helps to reduce computation and memory usage.

Figure 5 shows the histories of residual and error for each combination of methods. For each case, we show the relative residual associated with the preconditioned hybridized system and the one corresponding to the physical system. In contrast to the previous sections, the relative error is computed by comparing the numerical solution obtained at each iteration with that obtained with a direct solver.

We observe that the iterative solution of HDG is considerably slower compared to that of CHDG, regardless of whether CGNR or GMRES is used, both in terms of residuals and errors. The choice
<figure><img src="image_13.png" /><figcaption>Figure 4: Marmousi benchmark problem. Profile of velocity $c(\mathbf{x})$ (a), of density $\rho(\mathbf{x})$(b), mesh (c) and real part of the reference pressure field obtained with a direct solver (d).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 5: Marmousi benchmark problem. Relative residual (left) and relative error history (right) with different iterative schemes and DG methods. For the relative residual history, we consider for each case both the residual associated with the hybridized system (in $\times$, plain lines) and the one associated with the physical system (in physical norm, dashed lines). The relative error is computed by comparing the numerical physical solution obtained at each iteration with the reference numerical physical solution obtained with a direct solver.</figcaption></figure>

of numerical flux has a negligible effect on the convergence of the iterations with HDG (only one curve is shown for clarity), while the iterations with 'Upw' are slightly slower than with 'Sym-0'
when considering CHDG.

Comparing the different iterative strategies for solving the CHDG system with ‘Sym-0’, we observe
that the decay of both error and residual is faster with GMRES than with CGNR. The slopes of
the curves of error and residual are nonetheless similar. The fixed-point iteration performs well
until approximately 4 000 iterations, where the decay of error slows down.

One should note that these results do not give a final conclusion on which iterative scheme should be used in practice, since the runtime per iteration depends on the iterative scheme. The fixed-point iteration requires a single multiplication by the matrix of the system, while CGNR involves the multiplications by both the matrix and its adjoint at each iteration. The cost of GMRES increases at each iteration until the restart, and it is, overall, the most expensive method per iteration. Because these results are obtained with a MATLAB code that is not fully optimized, the observed runtimes are not reported here as they do not provide a representative comparison of efficiency. A performance analysis will be performed on a 3D parallel code under development.

6 Conclusion

In this work, we have extended and studied the CHDG method, originally proposed for Helmholtz problems in [34], in order to address wave propagation in heterogeneous media with piecewise constant physical coefficients. In particular, we considered formulations with standard upwind numerical fluxes, which are widely used in the literature, and symmetric numerical fluxes involving a general impedance operator $\mathcal{A}$. For both numerical fluxes, transmission variables are introduced at the interfaces between the elements in the hybridization procedure. In contrast to [34], the definition of the transmission variables involves the physical coefficients and the operator $\mathcal{A}$.

The CHDG hybridized system can be written in the form $(I - \Pi_S)g = b$. In the case of symmetric fluxes, if $\mathcal{A}$ is a positive and self-adjoint operator, it is proven that the operator $\Pi_S$ is a strict contraction, which makes it possible to solve the system with a fixed-point iteration. This result also holds for the standard upwind fluxes if the impedance $\eta$ is constant. Otherwise, $\Pi_S$ is not a contraction and the fixed-point iteration does not converge, as confirmed by the numerical results.

We have systematically studied and compared the different methods, numerical fluxes and iterative
schemes by considering 2D benchmark problems. Most of the observations that were made in the
homogeneous case in [34] remain valid. The fixed-point iteration applied to the CHDG system
converges for the cases where it is proven that the operator $\Pi_S$ is a strict contraction. While
the convergence of the fixed-point iteration is in general very fast for benchmarks with physical
dissipation, it is very slow for cavity cases. Considering the CGNR and GMRES iterations,
we observe that their convergence is nearly always faster with CHDG than with the standard
HDG method. The difference is very significant for the benchmarks with impedance boundary
conditions, and it is drastic for the Marmousi benchmark.

In this work, we have also investigated the use of symmetric fluxes with a second-order differential operator $\mathcal{A}$, inspired by recent work on domain decomposition methods. While this operator satisfies the assumptions of the theoretical framework, it always results in a slower convergence of the iterative schemes compared to the symmetric fluxes with a scalar operator. Comparing the upwind fluxes and the scalar symmetric fluxes, we have observed that, in the cases where they are not identical (i.e. when the impedance $\eta$ is discontinuous at the interfaces), they provide similar results, except for the fixed-point iteration.

This approach can be extended to 3D cases and to other physical problems, such as aeroacoustic,
electromagnetic [41] and elastic wave problems. We are currently making effort in these directions. We are investigating in more detail the computational aspects, as well as combinations with preconditioning and domain decomposition methods, to further accelerate the convergence of the iterations.

**Acknowledgments.** This work was supported in part by the ANR JCJC project WavesDG (research grant ANR-21-CE46-0010).

## References

[1] H. Barucq, A. Bendali, J. Diaz, and S. Tordeux. Local strategies for improving the conditioning of the plane-wave Ultra-Weak Variational Formulation. *Journal of Computational Physics*, 441:110449, 2021.

[2] H. Barucq, J. Diaz, R.-C. Meyer, and H. Pham. Implementation of hybridizable discontinuous Galerkin method for time-harmonic anisotropic poroelasticity in two dimensions. *International Journal for Numerical Methods in Engineering*, 122(12):3015-3043, 2021.

[3] H. Barucq, N. Rouxelin, and S. Tordeux. Construction and analysis of a HDG solution for the total-flux formulation of the convected Helmholtz equation. *Mathematics of Computation*, 92(343): 2097-2131, 2023.

[4] H. Barucq, A. Bendali, J. Diaz, and S. Tordeux. Revisiting the Plane-Wave Ultra-Weak Variational Formulation. *HAL preprint hal-04591459*, 2024.

[5] A. Bayliss, C. I. Goldstein, and E. Turkel. An iterative method for the Helmholtz equation. *Journal of Computational Physics*, 49:443-457, 1983.

[6] H. Bériot, A. Prinn, and G. Gabard. Efficient implementation of high-order finite elements for Helmholtz problems. *International Journal for Numerical Methods in Engineering*, 106(3):213-240, 2015.

[7] C. Carstensen, L. Demkowicz, and J. Gopalakrishnan. Breaking spaces and forms for the DPG method and applications including Maxwell equations. *Computers & Mathematics with Applications*, 72(3):494-522, 2016.

[8] O. Cessenat and B. Després. Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem. *SIAM Journal on Numerical Analysis*, 35(1):255-299, 1998.

[9] H. Chen, P. Lu, and X. Xu. A hybridizable discontinuous Galerkin method for the Helmholtz equation with high wave number. *SIAM Journal on Numerical Analysis*, 51(4):2166-2188, 2013.

[10] K. Christodoulou, O. Laghouche, M. S. Mohamed, and J. Trevelyan. High-order finite elements for the solution of Helmholtz problems. *Computers & Structures*, 191:129-139, 2017.

[11] E. T. Chung and P. Ciarlet. A staggered discontinuous Galerkin method for wave propagation in media with dielectrics and meta-materials. *Journal of Computational and Applied Mathematics*, 239: 189-207, 2013.

[12] B. Cockburn. *Static Condensation, Hybridization, and the Devising of the HDG Methods*, pages 129-177. Springer International Publishing, 2016.

[13] B. Cockburn, B. Dong, and J. Guzmán. A superconvergent LDG-hybridizable Galerkin method for second-order elliptic problems. *Mathematics of Computation*, 77(264):1887-1916, 2008.

[14] B. Cockburn, J. Gopalakrishnan, and R. Lazarov. Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems. *SIAM Journal on Numerical Analysis*, 47(2):1319-1365, 2009.

[15] F. Collino, P. Joly, and M. Lecouvez. Exponentially convergent non overlapping domain decomposition methods for the Helmholtz equation. *ESAIM: Mathematical Modelling and Numerical Analysis*, 54(3):775-810, 2020.

[16] S. Congreve, J. Gedicke, and I. Perugia. Robust adaptive hp-discontinuous Galerkin finite element methods for the Helmholtz equation. *SIAM Journal on Scientific Computing*, 41(2):A1121-A1147, 2019.
[17] B. Després. Méthodes de décomposition de domaine pour la propagation d'ondes en régime harmonique. Le théorème de Borg pour l'équation de Hill vectorielle. PhD thesis, 1991. Université de Paris-IX.

[18] B. Després, A. Nicolopoulos, and B. Thierry. Corners and stable optimized domain decomposition methods for the Helmholtz problem. *Numerische Mathematik*, 149:779–818, 2021.

[19] B. Engquist and L. Ying. Sweeping preconditioner for the Helmholtz equation: Moving Perfectly Matched Layers. *Multiscale Modeling & Simulation*, 9(2):686–710, 2011.

[20] Y. A. Erlangga, C. Vuik, and C. W. Oosterlee. On a class of preconditioners for solving the Helmholtz equation. *Applied Numerical Mathematics*, 50(3):409–425, 2004.

[21] O. G. Ernst and M. J. Gander. Why it is difficult to solve Helmholtz problems with classical iterative methods. *Numerical analysis of multiscale problems*, pages 325–363, 2011.

[22] C. Farhat, R. Tezaur, and J. Toivanen. A domain decomposition method for discontinuous Galerkin discretizations of Helmholtz problems with plane waves and Lagrange multipliers. *International journal for numerical methods in engineering*, 78(13):1513–1531, 2009.

[23] G. Gabard. Discontinuous Galerkin methods with plane waves for time-harmonic problems. *Journal of Computational Physics*, 225(2):1961–1984, 2007.

[24] M. J. Gander and H. Zhang. A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods. *SIAM Review*, 61(1):3–76, 2019.

[25] C. Geuzaine and J.-F. Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. *International Journal for Numerical Methods in Engineering*, 79(11): 1309–1331, 2009.

[26] G. Giorgiani, S. Fernández-Méndez, and A. Huerta. Hybridizable discontinuous Galerkin *p*-adaptivity for wave propagation problems. *International Journal for Numerical Methods in Fluids*, 72(12):1244–1262, 2013.

[27] M. Griesmaier and P. Roland. Error analysis for a hybridizable discontinuous Galerkin method for the Helmholtz equation. *Journal of Scientific Computing*, 49(3):291–310, 2011.

[28] J. S. Hesthaven and T. Warburton. *Nodal discontinuous Galerkin methods: algorithms, analysis, and applications*. Springer Science & Business Media, 2007.

[29] T. Huttunen, P. Monk, and J. P. Kaipio. Computational aspects of the ultra-weak variational formulation. *Journal of Computational Physics*, 182(1):27–46, 2002.

[30] R. J. LeVeque. *Finite Volume Methods for Hyperbolic Problems*. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2002.

[31] J. Li and J. S. Hesthaven. Analysis and application of the nodal discontinuous Galerkin method for wave propagation in metamaterials. *Journal of Computational Physics*, 258:915–930, 2014.

[32] A. Lieu, G. Gabard, and H. Bériot. A comparison of high-order polynomial and wave-based methods for Helmholtz problems. *Journal of Computational Physics*, 321:105–125, 2016.

[33] T. Lu, P. Zhang, and W. Cai. Discontinuous Galerkin methods for dispersive and lossy Maxwell's equations and pml boundary conditions. *Journal of Computational Physics*, 200(2):549–580, 2004.

[34] A. Modave and T. Chaumont-Frelet. A hybridizable discontinuous Galerkin method with characteristic variables for Helmholtz problems. *Journal of Computational Physics*, 493:112459, 2023.

[35] A. Modave, A. Royer, X. Antoine, and C. Geuzaine. A non-overlapping domain decomposition method with high-order transmission conditions and cross-point treatment for Helmholtz problems. *Computer Methods in Applied Mechanics and Engineering*, 368:113162, 2020.

[36] P. Monk, J. Schöberl, and A. Sinwel. Hybridizing Raviart-Thomas elements for the Helmholtz equation. *Electromagnetics*, 30(1-2):149–176, 2010.

[37] N. C. Nguyen, J. Peraire, and B. Cockburn. High-order implicit hybridizable discontinuous Galerkin methods for acoustics and elastodynamics. *Journal of Computational Physics*, 230(10):3695–3718, 2011.

[38] E. Parolin, D. Huybrechs, and A. Moiola. Stable approximation of Helmholtz solutions in the disk by evanescent plane waves. *ESAIM: M2AN*, 57(6):3499–3536, 2023.

[39] H. Pham, F. Faucher, and H. Barucq. Numerical investigation of stabilization in the Hybridizable Discontinuous Galerkin method for linear anisotropic elastic equation. *Computer Methods in Applied Mechanics and Engineering*, 428:117080, 2024.
[40] A. Quarteroni. *Numerical Models for Differential Problems*. Springer Publishing Company, Incorporated, 2nd edition, 2013.

[41] A. E. Rappaport, T. Chaumont-Frelet, and A. Modave. A hybridizable discontinuous Galerkin method with transmission variables for time-harmonic electromagnetic problems. *HAL preprint hal-05016261*, 2025.

[42] Y. Saad. *Iterative methods for sparse linear systems*. SIAM, 2003.

[43] C. Shi, J. Li, and C.-W. Shu. Discontinuous Galerkin methods for Maxwell's equations in Drude metamaterials on unstructured meshes. *Journal of Computational and Applied Mathematics*, 342: 147-163, 2018.

[44] E. F. Toro. *Riemann solvers and numerical methods for fluid dynamics: a practical introduction*. Springer Science & Business Media, 2013.

[45] A. Vion and C. Geuzaine. Improved sweeping preconditioners for domain decomposition algorithms applied to time-harmonic Helmholtz and Maxwell problems. *ESAIM: Proceedings and Surveys*, 61: 93-111, 2018.

[46] L. C. Wilcox, G. Stadler, C. Burstedde, and O. Ghattas. A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media. *Journal of Computational Physics*, 229 (24):9373-9396, 2010.