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

HAL
open science

# Construction and analysis of a HDG+ method for the diffusive-flux formulation of the convected Helmholtz equation

Hélène Barucq, Nathan Rouxelin, Sébastien Tordeux

► To cite this version:

Hélène Barucq, Nathan Rouxelin, Sébastien Tordeux. Construction and analysis of a HDG+ method for the diffusive-flux formulation of the convected Helmholtz equation. 2023. hal-04160700

HAL Id: hal-04160700

https://inria.hal.science/hal-04160700

Preprint submitted on 12 Jul 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
CONSTRUCTION AND ANALYSIS OF A HDG+ METHOD FOR THE
DIFFUSIVE-FLUX FORMULATION OF THE CONVECTED HELMHOLTZ
EQUATION *

HÉLÈNE BARUCQ<fnref n="1" />, NATHAN ROUXELIN<fnref n="1" />,<fnref n="2" /> AND SÉBASTIEN TORDEUX<fnref n="1" />

**Abstract.** We construct HDG methods based on the diffusive flux formulation of the convected acoustic wave equation. We mostly describe the HDG+ method which involves different polynomial degrees for approximating the unknowns, hence leading to a more efficient method with a super-convergence-like behaviour. A detailed analysis of the methods including local and global well-posedness, as well as convergence estimates is carried out. The HDG+ method is also compared with a more conventional HDG method to demonstrate its effectiveness.

**2020 Mathematics Subject Classification.** Primary 65N12, 65N30.

*This paper is dedicated to the memory of Prof. Francisco-Javier Sayas.*

INTRODUCTION

The numerical simulation of aeroacoustic waves has become an important problem for various applications ranging from *aeronautic noise* [LMG+20, BCD+15] to *solar physics* [GBD+17]. In solar physics, aeroacoustic waves model the solar oscillations that can be measured on the surface of the Sun and used to characterize the solar interior, see [Chr04]. This work is part of a scientific project aiming at constructing an efficient computational framework that can reproduce the observed data from the Sun. A longer term goal is to solve *inverse problems* to reconstruct the solar medium from those surface observations. With this in mind, we have chosen to work with the *hawen* simulation code [Fau21] which solves wave propagation forward and inverse problems based upon high-order finite element methods. Our aim is to equip *hawen* with a direct solver for aeroacoustics based on discontinuous finite elements, whose effectiveness has been proven in [FS20] for inverting acoustic and elastic data with Full Waveform Inversion (FWI) from terrestrial measurements. One of the reasons for using discontinuous finite element spaces is that, by not imposing continuity of shape functions between elements, discontinuous Galerkin methods are much more resistant to the effect of numerical pollution, which has been studied by many authors at the initiative of the seminal paper [BS97]. This effect must be

<footnote>*Keywords and phrases:* Hybridizable Discontinuous Galerkin Method (HDG), aeroacoustics, convected Helmholtz equation, harmonic regime, error analysis</footnote>

<footnote>* Nathan Rouxelin acknowledges financial support from e2s-UPPA, Maison Normande des Sciences du Numérique and Région Normandie.</footnote>

<footnote n="1">¹ Makutu, Inria, Université de Pau et des Pays de l'Adour, TotalEnergies, CNRS UMR 5142, France ; e-mail: helene.barucq@inria.fr & sebastien.tordeux@univ-pau.fr</footnote>

<footnote n="2">² Laboratoire de Mathématiques de l'INSA – UR 3226, INSA de Rouen-Normandie, Normandie Université, France ; e-mail: nathan.rouxelin@insa-rouen.fr</footnote>
taken into account all the more when the propagation domain is very large, in the sense that it contains a
very large number of characteristic wavelengths. Moreover, hp-adaptivity can be fully exploited which makes
the inverse problem easier to solve as it allows the use of a local degree of approximation on an unstructured
mesh. This allows variable parameters to be taken into account, improving the updating of the physical model
during FWI iterations. However, DG methods are known to be more expensive than continuous ones as they
use increased number of degrees of freedom. The reduction of the number of degrees of freedom is then
critical to build inversion software based on those methods. The *Hybridizable Discontinuous Galerkin* (HDG)
methods are attractive for this purpose. Indeed, HDG methods can be efficiently implemented by using a *static
condensation* process [Coc14] which eliminates the volumetric degrees of freedom by expressing them in terms
of a new unknown, the *numerical trace*, defined only on the skeleton of the mesh. Once the *global problem* for
the numerical trace has been solved, the full unknown can be reconstructed by solving small *local problems* in
parallel.

HDG methods have been used and validated for solving problems governed by elliptic equations in [CGL09, CDG+09,CC12,CC14]. As far as wave propagation is concerned, acoustic wave propagation has been considered in [GM11,GSV18,NPRC15], elastic wave propagation in [HPS17,BDMP21,CS13,FCS15,BCDL15], and Maxwell equations in [CQSS17, CQS18, CLOS20]. HDG method construction is rather similar to the one of mixed finite elements and the actual connection was first established by Cockburn and his coworkers in [CGS10]. For a self-contained introduction to the theory of HDG methods, we refer to [DS19]. For a comparison between HDG and Continuous Galerkin (CG) methods, we refer to [KSC12, YMKS16] where it is demonstrated that even the CG static condensation system rank is smaller than the DG one, its bandwith is really large, making HDG method clearly more efficient at high order. The relationship between HDG and Hybrid High-Order methods has been studied in [CDPE16]. It is shown that there exists a numerical flux associated with HHO formulation such that both methods can be compared.

It is worth noticing that other ways to reduce the computational cost of aeroacoustic simulations have been considered in the literature, such as *Trefftz methods* base on Generalized Plane-Waves (see [IGS22] and the references therein). For our purpose, we focus on finite element methods based on discontinuous polynomial approximations and the idea of using Trefftz methods will be addressed in a former work.

This work follows [BRT23], dealing with the *total-flux formulation* of the convected Helmholtz equation in which the vector unknown includes both diffusive and convective effects. It was proven that the standard HDG method associated with this first-order formulation is super-convergent. The well-posedness proofs rely on a discrete compactness argument, and it was proven that the method is well-posed if the mesh is sufficiently refined, but it was not possible to characterize the level of refinement needed. This work is motivated by the observation that those well-posedness proofs become more natural when considering the *diffusive-flux formulation* in which the vector unknown only encompasses diffusion effects. In particular, it was noted that the standard HDG method for the diffusive-flux formulation is suboptimal. To circumvent this difficulty, we consider the HDG+ method, also called *reduced-stabilization* HDG method, which was introduced in [Leh10]. HDG+ method mainly differs from HDG one by using a different polynomial degree of approximation for the volume unknown and for the unknown defined on the mesh skeleton. By this way, HDG+ solution to the diffusive-flux formulation of the convected Helmholtz equation is *super-convergent*. It is worth noting that due to the particular nature of the HDG+ method and to the convected nature of the problem, the convergence analysis is carried out using $L^2$-orthogonal projections instead of the tailored HDG projections. The analysis is therefore more technical, but is valid on general meshes and not only on simplicial ones. HDG+ methods have been considered for various applications in [CQSS17, Oik14, Oik16, Oik18, QSS16, QS16a, QS16b, Hun19] and to the best of our knowledge, the case of the convected Helmholtz equation has not been addressed yet.

The organization of this paper is as follows. **SECTION 1** deals with the convected Helmholtz, followed by useful notations and approximation settings are introduced in **SECTION 2**. The HDG+ formulation is constructed in **SECTION 3**. The well-posedness of the local problems is proven in **SECTION 4**, the error analysis is carried out
in SECTION 5 and the well-posedness of the global problem is proven in SECTION 6. Finally, some numerical results illustrating the theoretical properties of the method are presented in SECTION 7.

# 1. MODEL PROBLEM

As a model problem we consider the so-called *convected Helmholtz equation*

$$ \rho_0 (-\omega^2 p - 2i\omega\mathbf{v}_0 \cdot \nabla p + \mathbf{v}_0 \cdot \nabla(\mathbf{v}_0 \cdot \nabla p)) - \operatorname{div}(\rho_0 c_0^2 \nabla p) = s \quad (1.1) $$

where $\omega$ is the angular frequency, $\rho_0$ is the density of the fluid, $\mathbf{v}_0$ is the velocity of the fluid, $c_0$ is the adiabatic sound speed, and $s$ is the acoustic source.

Equation (1.1) is the simplest aeroacoustic models and therefore has a limited validity. This equation can be used for

* a *uniform background flow*, in this case the unknown $p$ can be interpreted as a pressure perturbation,

* a *potential background flow*, in this case the unknown $p$ should be interpreted as an *acoustic potential* and the physical quantities can be retrieved using the following identities

Pressure perturbation: $p' = -\rho_0 c_0 (-i\omega + \mathbf{v}_0 \cdot \nabla)p,$

Velocity perturbation: $\mathbf{v}' = -c_0 \nabla p,$

see [Pie90, Sec. II.].

In this paper, we only consider finite computational domains, which we denote by $\mathcal{O}$ and whose boundary is denoted by $\Gamma$. More precisely $\mathcal{O}$ is a bounded open subset of $\mathbb{R}^n$ with $n = 2$ or 3.

We will assume that the background flow is incompressible which leads to the following local mass conservation equation

$$ \operatorname{div}(\rho_0 \mathbf{v}_0) = 0. \quad (1.2) $$

Furthermore, we require some additional regularity for the velocity field $\rho_0 \mathbf{v}_0$ assuming that it is Lipschitz continuous, i.e. $\rho_0 \mathbf{v}_0 \in W^{1,\infty}(\mathcal{O})$. This will be useful to derive convergence estimates for the method as it allows us to estimate the difference between $\rho_0 \mathbf{v}_0$ and its average on a mesh element.

To get a mixed Discontinuous Galerkin approximation of the convected Helmholtz equation, we rewrite (1.3) as a system involving only first-order in space derivatives. We first combine the Laplace operator and the second-order convection term to obtain an anisotropic Laplace operator, which can be naturally handled in a HDG formulation. Using the mass conservation assumption (1.2), we have

$$ \rho_0 \mathbf{v}_0 \cdot \nabla (\mathbf{v}_0 \cdot \nabla p) = \operatorname{div} (\rho_0 \mathbf{v}_0 \mathbf{v}_0^T \nabla p). $$

Introducing the anisotropy tensor $\mathbf{K}_0 := \rho_0 (c_0^2 \mathbf{Id} - \mathbf{v}_0 \mathbf{v}_0^T)$, we obtain

$$ \rho_0 (-\omega^2 p - 2i\omega\mathbf{v}_0 \cdot \nabla p) - \operatorname{div}(\mathbf{K}_0 \nabla p) = s. \quad (1.3) $$

To lighten the notations in the remaining of this paper, we introduce the following vector field

$$ \mathbf{b}_0 := \rho_0 \mathbf{v}_0, $$

that satisfies the following incompressibility equation

$$ \operatorname{div}(\mathbf{b}_0) = 0. $$

If the background flow is subsonic, i.e.

$$ \inf_{\mathcal{O}} (c_0^2 - |\mathbf{v}_0|^2) > 0, \quad (1.4) $$

then we have that the following lemma.
**Lemma 1.1.** The anisotropy tensor $K_0$ is symmetric positive-definite and its spectrum is

$$ \mathrm{Sp}(K_0) = \{ \rho_0 c_0^2, \rho_0 (c_0^2 - |\mathbf{v}_0|^2) \} $$

**Proof:** We have $K_0 \mathbf{v}_0 = \rho_0 (c_0^2 - |\mathbf{v}_0|^2) \mathbf{v}_0$ and $K_0 \mathbf{u} = \rho_0 c_0^2 \mathbf{u}$ for all $\mathbf{u} \in \mathbf{v}_0^\perp$.

We can conclude that $-\mathrm{div}(K_0 \nabla p)$ is an elliptic operator.

To obtain a well-posed problem, the equation (1.3) must be closed by adding some boundary conditions on $\Gamma$.

In this paper, we only consider Neumann boundary conditions

$$ (K_0 \nabla p + 2i\omega p \mathbf{b}_0) \cdot \mathbf{n} = -g_N, \quad \text{on } \Gamma, \qquad (1.5) $$

where $\mathbf{n}$ is the outward-facing unitary normal vector to $\Gamma$. In this case, the bilinear form associated with the convected Helmholtz equation has a *coercive + compact* structure and is therefore of Fredholm type. This implies that the system (1.3)-(1.5) has a unique solution, except for some frequencies $\omega$ for which a resonant phenomenon can occur. In this paper, we will always assume that $\omega$ is not a resonant frequency.

To construct a HDG method, we need to rewrite (1.3) as a first-order in space system. We therefore introduce the *diffusive flux*

$$ \mathbf{q} := -K_0 \nabla p. $$

The resulting first-order formulation for (1.3) supplemented with the Neumann boundary condition (1.5) reads

$$ W_0 \mathbf{q} + \nabla p = 0, \qquad \text{in } \mathcal{O}, \qquad (1.6a) $$

$$ -\rho_0 \omega^2 p - 2i\omega \mathbf{b}_0 \cdot \nabla p + \mathrm{div}(\mathbf{q}) = s, \qquad \text{in } \mathcal{O}, \qquad (1.6b) $$

$$ (\mathbf{q} + 2i\omega p \mathbf{b}_0) \cdot \mathbf{n} = g_N, \qquad \text{on } \Gamma, \qquad (1.6c) $$

here $W_0$ is the inverse of $K_0$. Note that $K_0$ is always invertible as

$$ \det K_0 = \rho_0 c_0^2 (c_0^2 - |\mathbf{v}_0|^2) \neq 0, $$

and its inverse can be expressed as

$$ W_0 := K_0^{-1} = \frac{1}{\rho_0 c_0^2} \left[ \mathrm{Id} + \frac{\mathbf{v}_0 \mathbf{v}_0^T}{c_0^2 - |\mathbf{v}_0|^2} \right]. $$

thanks to the *Sherman-Morrison formula* [SM50].

Notice that even if we have chosen to work with second-order in frequency reading (1.6a)-(1.6b)-(1.6c), the resulting method can easily be adapted to obtain a first-order in frequency one. Indeed, by introducing

$$ \tilde{p} := -i\omega p, \quad \text{and} \quad \tilde{q} := \frac{1}{i\omega} K_0 \nabla \tilde{p}, $$

we can rewrite (1.3) as

$$ -i\omega W_0 \tilde{q} + \nabla \tilde{p} = 0, \qquad \text{in } \mathcal{O}, $$

$$ -i\omega \rho_0 \tilde{p} + 2\mathbf{b}_0 \cdot \nabla \tilde{p} + \mathrm{div}(\tilde{q}) = s, \qquad \text{in } \mathcal{O}, $$

$$ (i\omega \tilde{q} + 2\tilde{p} \mathbf{b}_0) \cdot \mathbf{n} = g_N, \qquad \text{on } \Gamma, $$

and the HDG methods can easily be adapted to this system.
## 2. APPROXIMATION SETTINGS

In this section, we introduce the notations and approximation spaces that will be used to construct the HDG methods considered in this paper.

### 2.1. Approximation spaces

We consider a mesh $\mathcal{T}_h$ of the domain $\mathcal{O}$. For a given element $K \in \mathcal{T}_h$ we denote its diameter by $h_K$ and we set

$$h := \max_{K \in \mathcal{T}_h} h_K.$$

For an element $K \in \mathcal{T}_h$, we denote by $\mathcal{F}(K)$ the set of its edges. We also consider

* The set of all edges of $\mathcal{T}_h$:

$$\mathcal{F}_h := \bigcup_{K \in \mathcal{T}_h} \mathcal{F}(K).$$

* The set of boundary edges:

$$\mathcal{F}_h^b := \{e \in \mathcal{F}_h \mid e \subset \Gamma\}.$$

* The set of interior edges:

$$\mathcal{F}_h^i := \{e \in \mathcal{F}_h \mid \exists K_+, K_- \in \mathcal{T}_h, K_- \neq K_+, e = \partial K_+ \cap \partial K_-\}.$$

**Remark 2.1.** We will assume that the mesh $\mathcal{T}_h$ has the usual shape-regularity property, see [EG04, Def. 1.107].

For $K \in \mathcal{T}_h$, we denote by $\mathcal{P}_k(K)$ the space of polynomial functions of total degree at most $k$ defined on $K$. We will also use the space of vectorial polynomials $\mathcal{P}_k(K) = \mathcal{P}_k(K)^n$. Even if those spaces can be defined for $k \ge 0$, in this paper we will usually assume that $k > 2$ as HDG methods of lower order have limited interest from a computational point of view. Indeed the key step in HDG methods is a static condensation process which consists in eliminating the interior degrees of freedom. The later do not exist for polynomial approximation of degree 1 or 2. Furthermore it was noted in [KSC12] that HDG method have a cost similar to CG methods for polynomial approximation of degree 5 or higher. However the static condensation of lower-order HDG methods has one practical interest: it leads to a mixed DG method which has the same cost as a primal DG method.

On each element $K \in \mathcal{T}_h$, we introduce the following approximation spaces for the pressure and the flux

$$V_h(K) := \mathcal{P}_k(K) \quad \text{for the flux } q_h,$$

$$W_h(K) := \mathcal{P}_\ell(K) \quad \text{for the potential } p_h,$$

with $\ell = k$ or $\ell = k + 1$.

As the approximation spaces are discontinuous, we introduce the numerical fluxes $\hat{q}_h$ and $\hat{p}_h$ which are designed to approximate the traces of $q$ and $p$ on the boundary of the elements. Those numerical fluxes also include some stabilization terms that ensure the stability of the DG method. In the particular case of HDG methods, a static condensation process is used to express $q_h$, $p_h$ and $\hat{q}_h$ as a function of $\hat{p}_h$. This leads to a so-called global problem whose unknown is the numerical flux $\hat{p}_h$. To approximate $\hat{p}_h$ we introduce the following space for $e \in \mathcal{F}(K)$

$$M_h(e) := \mathcal{P}_k(e).$$

As those approximation spaces are discontinuous, we can construct the global approximation spaces from the local ones

$$V_h := \{\sigma \in L^2(\mathcal{O}) \mid \sigma|_K \in V_h(K), \forall K \in \mathcal{T}_h\} \quad \text{for the flux } q_h,$$

$$W_h := \{p \in L^2(\mathcal{O}) \mid p|_K \in W_h(K), \forall K \in \mathcal{T}_h\} \quad \text{for the potential } p_h,$$

$$M_h := \{\mu \in L^2(\mathcal{F}_h) \mid \mu|_e \in M_h(e), \forall e \in \mathcal{F}_h\} \quad \text{for the trace } \hat{p}_h.$$
In FIGURE 2.1, we have depicted the differences in the *degrees of freedom* for the continuous (CG), discontinuous (DG) and hybridizable discontinuous (HDG) Galerkin methods. The degrees of freedom of the HDG methods are the ones associated with the numerical trace $\hat{p}_h$. If a mixed DG method is used, there are three unknowns for each degree of freedom, thus rendering those methods even more expensive. It was demonstrated in [KSC12] that HDG methods have numerical cost similar to the one of CG methods when they are properly implemented. Despite being expensive from a computational point of view, DG methods have been known to have some attractive properties. In particular, they can naturally be implemented for an arbitrary high-order with hp-adaptativity and in a parallel way. Using HDG methods therefore allows to keep those advantages of the DG methods for a reduced numerical cost, as it is illustrated in FIGURE 2.2. As the numerical cost of the method is directly linked to the number of degrees of freedom, we can clearly see that the HDG method is less expensive than the DG method.

<figure><img src="image_3.png" /><figcaption>FIGURE 2.1. Degrees of freedom for polynomial interpolation of degree 3 in 2D.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2.2. Number of degrees of freedom for a mixed DG method, a primal DG method (based on the second-order formulation) and a HDG method with interpolation degree $k$ in 3D with $10^3$ elements.</figcaption></figure>
## 2.2. Hermitian products and norms

The complex conjugate of $z$ is denoted by $\bar{z}$. For an element $K \in \mathcal{T}_h$, we denote the $L^2$-inner product<fnref n="1" /> and its associated norm by

$$ (u, v)_K := \int_K u \cdot \bar{v} \, d\boldsymbol{x} \quad \text{and} \quad \|u\|_K^2 := (u, u)_K. $$

We then introduce the broken inner product

$$ (u, v)_{\mathcal{T}_h} = \sum_{K \in \mathcal{T}_h} (u, v)_K, $$

and we denote by $\|\cdot\|_{\mathcal{T}_h}$ the associated norm. On the boundary of an element $K$, we also introduce the local duality bracket, which reduces to the $L^2(\partial K)$ inner product when $u$ and $v$ are sufficiently regular,

$$ \langle u, v \rangle_{\partial K} := \sum_{e \in \mathcal{F}(K)} \int_e u \cdot \bar{v} \, d\sigma, $$

and the associated norm is denoted by $\|\cdot\|_{\partial K}$. The broken product is then defined as

$$ \langle u, v \rangle_{\partial \mathcal{T}_h} := \sum_{K \in \mathcal{T}_h} \langle u|_K, v|_K \rangle_{\partial K}, $$

and we denote by $\|\cdot\|_{\partial \mathcal{T}_h}$ the associated norm. We will sometimes need to work on the interior edges only and we define the following broken product

$$ \langle u, v \rangle_{\partial \mathcal{T}_h \setminus \Gamma} := \sum_{K \in \mathcal{T}_h} \sum_{e \in \mathcal{F}(K) \cap \mathcal{F}_h^i} \langle u|_K, v|_K \rangle_e. $$

Finally, we also introduce the following weighted norms

$$ \begin{align*} \|u\|_{\rho_0, K}^2 &:= (\rho_0 u, u)_K && \text{which satisfies} & \quad \|u\|_{\rho_0, K} &\le \|\rho_0\|_{L^\infty(K)}^{\frac{1}{2}} \|u\|_K \\ \|\boldsymbol{q}\|_{\boldsymbol{W}_0, K}^2 &:= (\boldsymbol{W}_0 \boldsymbol{q}, \boldsymbol{q})_K && \text{which satisfies} & \quad \|\boldsymbol{q}\|_{\boldsymbol{W}_0, K} &\le C_{\boldsymbol{W}_0, K} \|\boldsymbol{q}\|_K \end{align*} $$

where

$$ C_{\boldsymbol{W}_0, K} = \left( \max_K \frac{1}{\rho_0 (c_0^2 - |\boldsymbol{v}_0|^2)} \right)^{\frac{1}{2}} $$

is the largest eigenvalue of $\boldsymbol{W}_0$ in $K$, we recall that the spectrum of $\boldsymbol{K}_0$ (and therefore of $\boldsymbol{W}_0 := \boldsymbol{K}_0^{-1}$) was given in LEMMA 1.1.

## 2.3. Edges, jumps and averages

Discontinuity at the interface between elements distinguish DG formulations from the CG ones. For stability and implementation purposes, it is then required to define quantities related to the edges of the elements.

For an interior face $\mathcal{F}_h^i \ni e = \partial K_+ \cap \partial K_-$, we denote by $\boldsymbol{n}^+$ (resp. $\boldsymbol{n}^-$) the unitary outgoing normal vector of $\partial K_+$ (resp. $\partial K_-$). We will always assume that the flow $\boldsymbol{v}_0$ goes from $K_-$ to $K_+$, as depicted on FIGURE 2.3. If $e$ is a boundary edge, then $\boldsymbol{n}$ denotes the unitary normal vector outwardly directed to $\mathcal{O}$.

<footnote n="1">¹For vector fields, the $\mathbb{R}^n$ dot-product is used inside the integral as the conjugate is already applied.</footnote>
<figure><img src="image_5.png" /><figcaption>FIGURE 2.3. Normal vectors on an interior face</figcaption></figure>

We will often use the *average operator* defined on the interior and boundary faces by

$$ \text{On } \mathcal{F}_h^i \ni e = \partial K_+ \cap \partial K_-, \qquad \{\{\varphi\}\}_e := \frac{1}{2}(\varphi^+ + \varphi^-), $$

$$ \text{On } \mathcal{F}_h^b \ni e = \partial K \cap \Gamma, \qquad \{\{\varphi\}\}_e := \frac{1}{2}\varphi, $$

where $\varphi$ can either be a scalar or vectorial quantity. We will also make frequent use of the *jump operator* defined on the interior and boundary faces by

$$ \text{On } \mathcal{F}_h^i \ni e = \partial K_+ \cap \partial K_-, \qquad \llbracket q \rrbracket_e := q^+ \cdot n^+ + q^- \cdot n^-, $$

$$ \text{On } \mathcal{F}_h^b \ni e = \partial K \cap \Gamma, \qquad \llbracket q \rrbracket_e := q \cdot n, $$

for a vectorial quantity. Notice that with this definition, the jump operator only controls the normal part of the vector. For a scalar quantity, the *jump operator* is defined on the interior and boundary faces by

$$ \text{On } \mathcal{F}_h^i \ni e = \partial K_+ \cap \partial K_-, \qquad \llbracket p \rrbracket_e := p^+ n^+ + p^- n^-, $$

$$ \text{On } \mathcal{F}_h^b \ni e = \partial K \cap \Gamma, \qquad \llbracket p \rrbracket_e := pn, $$

for a scalar quantity. A sketch of those quantities is given in FIGURE 2.4.

<figure><img src="image_6.png" /><figcaption>FIGURE 2.4. 1D-sketch of the jump and average on an interior node</figcaption></figure>
## 3. CONSTRUCTION OF THE METHOD

In this paper we consider both HDG and HDG+ variants of the method. The construction is detailed for the HDG+ method as it is the most complicated case. This construction can then be adapted to the standard HDG method without any major difficulties.

On an element $K \in \mathcal{T}_h$, the construction of the HDG(+) methods relies on the following integration by parts formula

$$ \int_K \mathbf{W}_0 \mathbf{q} \cdot \bar{\mathbf{r}} \, d\mathbf{x} - \int_K p \operatorname{div}(\bar{\mathbf{r}}) \, d\mathbf{x} + \langle p, \mathbf{r} \cdot \mathbf{n} \rangle_{\partial K} = 0, \quad (3.1a) $$

$$ -\omega^2 \int_K \rho_0 p \bar{w} \, d\mathbf{x} - 2i\omega \int_K \mathbf{b}_0 \cdot \nabla p \bar{w} \, d\mathbf{x} + \int_K \operatorname{div}(\mathbf{q}) \bar{w} \, d\mathbf{x} = \int_K s \bar{w} \, d\mathbf{x}, \quad (3.1b) $$

which is valid for all $(\mathbf{q}, p), (\mathbf{r}, w) \in \mathbf{H}_{\operatorname{div}}(K) \times H^1(K)$.

### 3.1. Choice of approximation spaces

For the HDG+ method, the choice of approximation spaces is different from the choice made for the previous HDG method. We consider the following *local approximation spaces*

$$ \mathbf{V}_h(K) = \mathcal{P}_k(K), \quad \text{for the flux } \mathbf{q}_h, $$

$$ W_h(K) = \mathcal{P}_{k+1}(K), \quad \text{for the pressure } p_h, $$

where $k \ge 2$ is the degree of the method. The use of a higher polynomial degree for $p_h$ is the distinctive feature of the HDG+ method.

### 3.2. Introduction of the hybrid unknown

As we did before, we introduce the numerical trace $\widehat{p_h}$ which approximates $p$ on the skeleton $\mathcal{F}_h$ of the mesh. As before the boundary integral in (3.1a) will be discretized as

$$ \langle p, \mathbf{r} \cdot \mathbf{n} \rangle_{\partial K} \quad \text{becomes} \quad \langle \widehat{p_h}, \mathbf{r}_h \cdot \mathbf{n} \rangle_{\partial K}. $$

For the HDG+ method, we use the following approximation space for $\widehat{p_h}$

$$ M_h(e) = \mathcal{P}_k(e), \quad \forall e \in \mathcal{F}(K). $$

With this choice, $p_h$ and $\widehat{p_h}$ do not have the same polynomial degree and we therefore have two approximations of $p$ with different polynomial degrees on the skeleton of the mesh. We therefore need to change the *penalization term* to

$$ \tau P_M(p_h - \widehat{p_h}), \quad (3.2) $$

where $P_M$ is the $L^2$-orthogonal projection onto $M_h$. This is called the *reduced stabilization* or *Lehrenfeld-Schöberl* stabilization and it was introduced in [Lch10]. Even if the result was known by Lehrenfeld and Schöberl, the first published analysis of the HDG+ method was given in [Oik16]. It allows to get convergence rate of $k+2$ for $p_h$ for the cost of a method of degree $k$. A large penalization parameter $\tau \sim h_K^{-1}$ is needed to obtain optimal convergence as it will be detailed in SECTION 5.

*Remark 3.1.* At the discrete level, we can rewrite (3.2) as

$$ \tau P_M(p_h - \widehat{p_h}) = \tau(P_M p_h - \widehat{p_h}), $$

however the left-hand side is consistent, but the right-hand side is not.
### 3.3. Local problems

Assuming that the numerical trace $\hat{p}_h$ is known, we approximate the weak formulation (3.1a)-(3.1b) on an element $K \in \mathcal{T}_h$ leading to the so-called local problem : seek $(\mathbf{q}_h, p_h) \in \mathbf{V}_h(K) \times W_h(K)$ such that

$$ (\mathbf{W}_0 \mathbf{q}_h, \mathbf{r}_h)_K - (p_h, \mathrm{div}(\mathbf{r}_h))_K + \langle \hat{p}_h, \mathbf{r}_h \cdot \mathbf{n} \rangle_{\partial K} = 0, \quad (3.3a) $$

$$ \begin{aligned} & -\omega^2 (\rho_0 p_h, w_h)_K - 2i\omega (\mathbf{b}_0 \cdot \nabla p_h, w_h)_K + (\mathrm{div}(\mathbf{q}_h), w_h)_K \\ & \qquad + 2i\omega \langle \tau P_M (p_h - \hat{p}_h), w_h \rangle_{\partial K} = (s, w_h)_K, \quad (3.3b) \end{aligned} $$

for all $(\mathbf{r}_h, w_h) \in \mathbf{V}_h(K) \times W_h(K)$. An equation for $\hat{p}_h$ will be constructed in the next paragraph. On $K$, the system (3.3a)-(3.3b) allows to reconstruct the local fields $(\mathbf{q}_h^K, p_h^K)$ when the source term $s$ and the numerical trace $\hat{p}_h$ are known. We can therefore write those local fields as

$$ \mathbf{q}_h^K = \mathbf{Q}^K \hat{p}_h + \mathbf{Q}_{\mathrm{src}}^K s, $$

$$ p_h^K = \mathbf{P}^K \hat{p}_h + \mathbf{P}_{\mathrm{src}}^K s, $$

where $(\mathbf{Q}^K \hat{p}_h, \mathbf{P}^K \hat{p}_h)$ satisfy (3.3a)-(3.3b) with $s = 0$, and $(\mathbf{Q}_{\mathrm{src}}^K s, \mathbf{P}_{\mathrm{src}}^K s)$ satisfy (3.3a)-(3.3b) with $\hat{p}_h = 0$. The linear operators $\mathbf{Q}^K$ and $\mathbf{P}^K$ are called local solvers, they are well-defined as it will be proven in THEOREM 1. The local problems are therefore independent of each other, and the numerical trace $\hat{p}_h$ is the only quantity that couples the elements together.

Following [QS16a], we have introduced a second penalization parameter $\tau_{\mathrm{upw}}$ defined by

$$ \tau_{\mathrm{upw}} := \max(\mathbf{b}_0 \cdot \mathbf{n}, 0). $$

To understand why this second parameter is required, we recall that in HDG methods the penalization serves two purposes:

(1) it enforces the Dirichlet boundary condition for the local problems,

(2) it controls the stability of the method.

Here as $\mathbf{q}_h$ does not take the convection into account, the penalization term (3.2) with $\tau$ only stabilizes the diffusion. Therefore we had a second penalization to stabilize the convection. We denoted it $\tau_{\mathrm{upw}}$ as it leads to an upwinding behavior that will be detailed in the next paragraph.

*Remark 3.2.* The restriction $(\mathbf{q}|_K, p|_K, p|_{\partial K})$ of the exact solution $(\mathbf{q}, p)$ satisfies the local problem (3.3a)-(3.3b) on each element $K$ of the mesh.

### 3.4. Transmission condition

Due to the discontinuous nature of the approximation spaces, the local problems (3.3a)-(3.3b) are independent of each other. To obtain a coupling between those problems, we introduce the following numerical flux

$$ \hat{\mathbf{q}}_h \cdot \mathbf{n} := \mathbf{q}_h \cdot \mathbf{n} + 2i\omega\tau P_M (p_h - \hat{p}_h), $$

where $\tau = \mathcal{O}(h_K^{-1})$. The quantity $\hat{\mathbf{q}}_h \cdot \mathbf{n}$ only ensures the normal continuity of the diffusive flux between two elements. To ensure the normal continuity of the total flux we add a second numerical flux

$$ 2i\omega p_h \hat{\mathbf{b}}_0 \cdot \mathbf{n} := 2i\omega(\mathbf{b}_0 \cdot \mathbf{n})\hat{p}_h + 2i\omega\tau_{\mathrm{upw}}(p_h - \hat{p}_h). \quad (3.4) $$
It is important to notice that this flux has an upwind behavior. Let $e = \partial K_+ \cap \partial K_-$ be an interior edge with $b_0 \cdot n_- > 0$ on $\partial K_-$. We have

$$ \text{On } \partial K_-: \quad \tau_{\text{upw}} := \max(b_0 \cdot n, 0) = b_0 \cdot n, \quad \text{so} \quad 2i\omega \widehat{p_h b_0} \cdot n = 2i\omega(b_0 \cdot n)p_h, $$

$$ \text{On } \partial K_+: \quad \tau_{\text{upw}} := \max(b_0 \cdot n, 0) = 0, \quad \text{so} \quad 2i\omega \widehat{p_h b_0} \cdot n = 2i\omega(b_0 \cdot n)\widehat{p_h}. $$

So on the outflow boundary we use the interior value $p_h$, whereas on the inflow boundary we use the trace value $\widehat{p_h}$.

Finally we write the transmission condition as

$$ \left\langle \left( \widehat{q_h} - 2i\omega p_h \mathbf{b_0} \right) \cdot \mathbf{n}, \mu_h \right\rangle_{\partial \mathcal{T}_h \setminus \Gamma} + \langle \widehat{p_h} - g_D, \mu_h \rangle_{\Gamma_D} = \langle g_N, \mu_h \rangle_{\Gamma_N}. \quad (3.6) $$

This formulation enforces normal continuity of the total flux between the elements and the boundary conditions on $\Gamma_D$ and $\Gamma_N$.

*Remark 3.3.* To ensure the well-posedness of the local problems, the second penalization must be

$$ \tau_{\text{upw}}(p_h - \widehat{p_h}), \quad \text{rather than} \quad \tau_{\text{upw}}(P_M p_h - \widehat{p_h}), $$

as it will be proven later in **THEOREM 1**.

*Remark 3.4.* We would like to point out the main theoretical difficulty of this method : when the background flow is not constant, the second flux (3.4) leads to non-polynomial terms on the skeleton. This is usually avoided as much as possible in HDG methods.

### 3.5. Compact formulation of the methods:

HDG methods are usually stated in a compact way that can be obtained by summing the local problems (3.3a)-(3.3b) over the mesh elements and by adding the transmission condition (3.6). This formulation reads : seek $(q_h, p_h, \widehat{p_h}) \in V_h \times W_h \times M_h$, such that

$$ (W_0 q_h, r_h)_{\mathcal{T}_h} - (p_h, \text{div}(r_h))_{\mathcal{T}_h} + \langle \widehat{p_h}, r_h \cdot n \rangle_{\partial \mathcal{T}_h} = 0, \quad (3.7a) $$

$$ -\omega^2 (\rho_0 p_h, w_h)_{\mathcal{T}_h} - 2i\omega (b_0 \cdot \nabla p_h, w_h)_{\mathcal{T}_h} + (\text{div}(q_h), w_h)_{\mathcal{T}_h} \quad (3.7b) $$

$$ +2i\omega \langle \tau P_M (p_h - \widehat{p_h}) - \tau_{\text{upw}}(p_h - \widehat{p_h}), w_h \rangle_{\partial \mathcal{T}_h} = (s, w_h)_{\mathcal{T}_h} $$

$$ \left\langle \left( \widehat{q_h} - 2i\omega p_h \mathbf{b_0} \right) \cdot \mathbf{n}, \mu_h \right\rangle_{\partial \mathcal{T}_h} = \langle g_N, \mu_h \rangle_{\Gamma_N}, \quad (3.7c) $$

for all $(r_h, w_h, \mu_h) \in V_h \times W_h \times M_h$.

### 3.6. Condensed variational formulation

The compact formulation (3.7a)-(3.7b)-(3.7c) cannot be used to efficiently implement the HDG method, indeed with this formulation it is not clear how the global problem for $\widehat{p_h}$ only can be obtained. To describe this process, we will now write a condensed variational formulation for $\widehat{p_h}$ only. On an element $K \in \mathcal{T}_h$, we write the local fields as

$$ q_h^K = Q^K \widehat{p_h} + Q_{\text{src}}^K s, $$

$$ p_h^K = P^K \widehat{p_h} + P_{\text{src}}^K s, $$

where $(Q^K \widehat{p_h}, P^K \widehat{p_h})$ satisfy (3.3a)-(3.3b) with $s = 0$, and $(Q_{\text{src}}^K s, P_{\text{src}}^K s)$ satisfy (3.3a)-(3.3b) with $\widehat{p_h} = 0$. The linear operators $Q^K$ and $P^K$ are called local solvers, they are well-defined as it will be proven in **THEOREM 1**.
We can therefore rewrite the transmission condition (3.7c) as

$$a_h(\widehat{p}_h, \mu_h) = \ell_h(\mu_h), \qquad (3.8)$$

where

$$a_h(\widehat{p}_h, \mu_h) := \langle \mathbf{Q}^K \widehat{p}_h \cdot \mathbf{n} + 2i\omega\tau P_M(\mathbf{P}^K \widehat{p}_h - \widehat{p}_h) - 2i\omega\tau_{\text{upw}}(\mathbf{P}^K \widehat{p}_h - \widehat{p}_h), \mu_h \rangle_{\partial \mathcal{T}_h}$$

and the terms involving $\mathbf{Q}^K_{\text{src}}$ and $\mathbf{P}^K_{\text{src}}$ are put into $\ell_h$. Equation (3.8) is the so-called global problem and is the main equation of the HDG method. As the local solvers satisfy the local problem (3.3a)-(3.3b), we can obtain the following characterization for $a_h$

$$a_h(\widehat{p}_h, \mu_h) = (\mathbf{W}_0 \mathbf{Q}^K \widehat{p}_h, \mathbf{Q}^K \mu_h)_{\mathcal{T}_h} - \omega^2 (\rho_0 \mathbf{P}^K \widehat{p}_h, \mathbf{P}^K \mu_h)_{\mathcal{T}_h} - 2i\omega (\mathbf{b}_0 \cdot \nabla \mathbf{P}^K \widehat{p}_h, \mathbf{P}^K \mu_h)_{\mathcal{T}_h} \\ + 2i\omega \langle \tau P_M(\mathbf{P}^K \widehat{p}_h - \widehat{p}_h), P_M(\mathbf{P}^K \mu_h - \mu_h) \rangle_{\partial \mathcal{T}_h} - 2i\omega \langle \tau_{\text{upw}}(\mathbf{P}^K \widehat{p}_h - \widehat{p}_h), \mathbf{P}^K \mu_h - \mu_h \rangle_{\partial \mathcal{T}_h},$$

by following [CGL09, Sec 2.3]. In particular, this characterization shows that the global problem of the HDG+ method has a structure which is similar to the one of the convected Helmholtz equation. In particular testing with $\mu_h = \widehat{p}_h$ leads to a global discrete Garding-like inequality. From a computational point of view, we proceed as described in ALGORITHM 1.

<table><caption>**Algorithm 1:** Solving HDG+</caption><tr><td>1</td><td>for K ∈ T_h do</td><td></td></tr><tr><td>2</td><td>/* Assembling step */</td><td></td></tr><tr><td>3</td><td>Construct the local solvers Q^K, and P^K</td><td></td></tr><tr><td>3</td><td>Add local contribution to the global problem (3.8)</td><td></td></tr><tr><td>4</td><td>Solve (3.8) for p_h using a direct solver</td><td>// Main linear system to solve</td></tr><tr><td>5</td><td>for K ∈ T_h do</td><td></td></tr><tr><td>6</td><td>/* Reconstruction step */</td><td></td></tr><tr><td>6</td><td>Reconstruct the local unknowns p_h^K = P^K \widehat{p}_h + P_{\text{src}}^K s and q_h^K = Q^K \widehat{p}_h + Q_{\text{src}}^K s</td><td></td></tr></table>

This algorithm is the blueprint of the practical implementation of the HDG method which is discussed in [Rou21].

## 3.7. Adaptation to a standard HDG method

It is also possible to apply a standard HDG method to this formulation and thus having the same polynomial degree for approximating $p_h$ and $q_h$, that is:

$$\begin{align*} V_h(K) &= \mathcal{P}_h(K), & \text{for the flux } q_h, \\ W_h(K) &= \mathcal{P}_h(K), & \text{for the pressure } p_h. \end{align*}$$

In this case, as $W_h$ and $M_h$ are spaces of polynomials having the same degree, the projection term becomes simpler, indeed

$$P_M p_h = p_h.$$

For this formulation, a large penalization parameter is no longer required, and we only need $\tau = O(1)$. However this HDG method is less efficient than the HDG+ method. In particular the upwinding term $\tau_{\text{upw}}$ is crucial to ensure the stability of the method. To illustrate this phenomenon, we have depicted numerical results obtained by using the parameters of SECTION 7 in FIGURE 3.1. When the upwinding term is not introduced, instabilities are clearly visible inside the domain. However adding this term to the formulation leads to a sub-optimal convergence rate for the vector unknown $q$ by following the error analysis of the HDG+ method carried out in SECTION 5. For the HDG+ method, this loss in the error rate is compensated by the higher interpolation degree used for the scalar unknown $p$.
<figure><img src="image_7.png" /><figcaption>FIGURE 3.1. Influence of the upwinding term for the standard HDG method</figcaption></figure>

## 4. LOCAL SOLVABILITY

It is worth remembering that HDG methods were originally developed for strongly elliptic problems and that harmonic wave equations are only quasi-elliptic. It is well-known that solving those equations with Dirichlet boundary conditions<fnref n="2" /> leads the resonance phenomenon. In this section we will show that the static condensation process is well-defined when the mesh is fine enough, ie. the local problem does not produce resonance. Before actually showing the local solvability, we need to prove the following lemma.

**Lemma 4.1.** If $p \in H^1(K)$ and $b_0 \in L^\infty(K) \cap C(O)$, where $C(O)$ is the space of vector functions continuous in the domain $O$, then the following identity holds

$$ \Re \langle (pb_0, \nabla p)_K, p, p \rangle_{\partial K} = \frac{1}{2} \langle (b_0 \cdot n)p, p \rangle_{\partial K}. $$

*Proof.* We use an integration by parts to obtain a relationship between $(pb_0, \nabla p)_K$ and its complex conjugate:

$$ \begin{aligned} 2\Re \langle (pb_0, \nabla p)_K &= (pb_0, \nabla p)_K + \overline{(pb_0, \nabla p)_K} \\ &= (pb_0, \nabla p)_K + (\nabla p, pb_0)_K \\ &= -(\operatorname{div}(pb_0), p)_K + \langle (b_0 \cdot n)p, p \rangle_{\partial K} + (\nabla p, pb_0)_K \\ (\operatorname{div}(b_0) = 0) &= -(\nabla p, pb_0)_K + \langle (b_0 \cdot n)p, p \rangle_{\partial K} + (\nabla p, pb_0)_K \\ &= \langle (b_0 \cdot n)p, p \rangle_{\partial K}. \end{aligned} $$

□

We can now state and prove the main result of this section.

**Theorem 1** (Local solvability). If the penalization parameter $\tau$ satisfies

$$ \forall e \in \mathcal{F}(K), \quad \tau|_e < 0, \qquad (4.1) $$

if the pulsation $\omega$ and the diameter $h_K$ of each element $K \in \mathcal{T}_h$ satisfy

$$ \omega h_K < \frac{-C_{W_0,K} \|b_0\|_{L^\infty(K)} + \left(C_{W_0,K}^2 \|b_0\|_{L^\infty(K)}^2 + \|\rho_0\|_{L^\infty(K)}\right)^{\frac{1}{2}}}{C_{W_0,K} C \|\rho_0\|_{L^\infty(K)}}, \qquad (4.2) $$

<footnote n="2">²Which is what the local solver does.</footnote>
where $C > 0$ is a constant that depends only on the shape regularity of $K$, and if the diameter $h_K$ of each element $K \in \mathcal{T}_h$ is sufficiently small for the sign of $b_0 \cdot n$ to be fixed on each edge of $K^3$, then the local solver

$$ (\widehat{p}_h, s) \mapsto (p_h^K, q_h^K) $$

is well-posed on each element $K \in \mathcal{T}_h$.

*Proof.* As the local problems have a finite dimension, we only need to prove uniqueness of the solution. We therefore assume that $\widehat{p}_h = s = 0$ and we need to prove that the system

$$ (W_0 q_h, r_h)_K - (p_h, \operatorname{div}(r_h))_K = 0, \quad (4.3a) $$

$$ -\omega^2 (\rho_0 p_h, w_h)_K - 2i\omega (b_0 \cdot \nabla p_h, w_h)_K $$

$$ + (\operatorname{div}(q_h), w_h)_K + 2i\omega \langle \tau P_M p_h - \tau_{\operatorname{upw}} p_h, w_h \rangle_{\partial K} = 0, \quad (4.3b) $$

has only one solution: $(p_h, q_h) = (0, 0)$.

We will prove the theorem by contradiction. We therefore assume that the system (4.3a)-(4.3b) has a non-zero solution $(p_h, q_h)$.

**Step 1:** An energy-like system

We begin by testing (4.3b) with $w_h = p_h$

$$ -\omega^2 \|p_h\|_{\rho_0, K}^2 + 2i\omega (p_h b_0, \nabla p_h)_K + (\operatorname{div}(q_h), p_h)_K + \langle 2i\omega\tau P_M p_h - 2i\omega\tau_{\operatorname{upw}} p_h, p_h \rangle_{\partial K} = 0 \quad (4.4) $$

Then, (4.3a) is tested with $r_h = q_h$ and conjugated:

$$ \|q_h\|_{W_0, K}^2 - (\operatorname{div}(q_h), p_h)_K = 0 \quad (4.5) $$

We now add (4.4) and (4.5) leading to

$$ \|q_h\|_{W_0, K}^2 - \omega^2 \|p_h\|_{\rho_0, K}^2 + 2i\omega (p_h b_0, \nabla p_h)_K + \langle 2i\omega\tau P_M p_h - 2i\omega\tau_{\operatorname{upw}} p_h, p_h \rangle_{\partial K} = 0 \quad (4.6) $$

We now obtain the following system by taking the real and imaginary parts of (4.6)

$$ \Re \mathcal{E}: \qquad \|q_h\|_{W_0, K}^2 - \omega^2 \|p_h\|_{\rho_0, K}^2 - 2\omega \Im (p_h b_0, \nabla p_h)_K = 0 \quad (4.7) $$

$$ \Im \mathcal{E}: \qquad \Re \mathcal{E} (p_h b_0, \nabla p_h)_K + \langle \tau P_M p_h, p_h \rangle_{\partial K} = \langle \tau_{\operatorname{upw}} p_h, p_h \rangle_{\partial K} \quad (4.8) $$

Indeed, as $P_M p_h \in M_h$ and $\tau$ is constant on each edge, one has

$$ \langle \tau P_M p_h, p_h \rangle_{\partial K} = \langle \tau P_M p_h, P_M p_h \rangle_{\partial K} \in \mathbb{R} $$

**Step 2:** We focus on (4.8) to express $p_h|_{\partial K}$.

By the LEMMA 4.1 we have

$$ \Re \mathcal{E} (p_h b_0, \nabla p_h)_K = \frac{1}{2} \langle (b_0 \cdot n) p_h, p_h \rangle_{\partial K} $$

and (4.8) becomes

$$ \frac{1}{2} \langle (b_0 \cdot n) p_h, p_h \rangle_{\partial K} + \langle \tau P_M p_h, p_h \rangle_{\partial K} = \langle \tau_{\operatorname{upw}} p_h, p_h \rangle_{\partial K} $$

For the sake of simplicity, we have assumed that the sign of $b_0 \cdot n$ is constant on each edge. For a given edge $e \in \mathcal{F}(K)$, the three following cases are exhaustive:

<footnote n="3">³This assumptions amounts to requiring than the inflow and outflow boundaries of $K$ are separated, which is a standard assumption to study problems with connvection.</footnote>
• Case 1: $b_0 \cdot n < 0$: therefore $\tau_{upw} := \max(b_0 \cdot n, 0) = 0$ and

$$ \underbrace{\langle \tau P_M p_h, P_M p_h \rangle_{\partial K}}_{\le 0 \text{ by (4.1) as } \tau < 0} - \frac{1}{2} \langle |b_0 \cdot n| p_h, p_h \rangle_{\partial K} = 0 $$

• Case 2: $b_0 \cdot n > 0$: therefore $\tau_{upw} := \max(b_0 \cdot n, 0) = b_0 \cdot n$ and

$$ \frac{1}{2} \langle |\mathbf{b}_0 \cdot \mathbf{n}| p_h, p_h \rangle_{\partial K} + \langle \tau P_M p_h, P_M p_h \rangle_{\partial K} = \langle |\mathbf{b}_0 \cdot \mathbf{n}| p_h, p_h \rangle_{\partial K} \\ \iff \underbrace{\langle \tau P_M p_h, P_M p_h \rangle_{\partial K}}_{\le 0 \text{ by (4.1) as } \tau < 0} - \frac{1}{2} \langle |\mathbf{b}_0 \cdot \mathbf{n}| p_h, p_h \rangle_{\partial K} = 0 $$

• Case 3: $b_0 \cdot n = 0$: in this case we only have

$$ \langle \tau P_M p_h, P_M p_h \rangle_{\partial K} = 0. $$

In the first two cases we have $p_h|_{\partial K} = P_M p_h = 0$ and in the third one we only have $P_M p_h = 0$. In particular, the following identity holds for all the three previous cases

$$ \int_{\partial K} p_h d\sigma = 0, $$

indeed as $P_M p_h$ is the $L^2$-orthogonal projection of $p_h$ onto $M_h$, we have

$$ \int_{\partial K} P_M p_h \overline{\mu_h} d\sigma = \int_{\partial K} p_h \overline{\mu_h} d\sigma, \quad \forall \mu_h \in M_h := \prod_{e \in \mathcal{F}(K)} \mathcal{P}_k(e), \quad (4.9) $$

and the previous identity is obtained by taking $\mu_h = 1$.

**Step 3: Contradiction**

As $p_h \in \mathcal{P}_{k+1}(K)$, we have $p_h \in H^1(K)$ and the following Poincaré-Friedrichs inequality holds<fnref n="4" />

$$ \|p_h\|_K \le C h_K \|\nabla p_h\|_K, $$

see [EG04, Lemma B.66] with $f(v) = \int_{\partial K} v d\sigma$. The constant $C$ is the same one as in (4.2).

Going back to (4.3a), integrating by parts and testing it with $r_h = \nabla p_h$ we have

$$ \begin{aligned} \|\nabla p_h\|_K^2 &= |(\mathbf{W}_0 \mathbf{q}_h, \nabla p_h)_K| \\ &\le C_{\mathbf{W}_0, K} \|\mathbf{q}_h\|_{\mathbf{W}_0, K} \|\nabla p_h\|_K \\ \|\nabla p_h\|_K &\le C_{\mathbf{W}_0, K} \|\mathbf{q}_h\|_{\mathbf{W}_0, K} \end{aligned} \quad (4.10) $$

On the other hand, from (4.7) we see that

$$ \begin{aligned} \|\mathbf{q}_h\|_{\mathbf{W}_0, K}^2 &= \omega^2 \|p_h\|_{\rho_0, K}^2 + 2\omega \mathfrak{Im} (p_h \mathbf{b}_0, \nabla p_h)_K \\ &\le \omega^2 \|\rho_0\|_{L^\infty(K)} \|p_h\|_K^2 + 2\omega \|\mathbf{b}_0\|_{L^\infty(K)} \|p_h\|_K \|\nabla p_h\|_K \\ \|\mathbf{q}_h\|_{\mathbf{W}_0, K}^2 &\le C^2 \|\rho_0\|_{L^\infty(K)} \omega^2 h_K^2 \|\nabla p_h\|_K^2 + 2C \|\mathbf{b}_0\|_{L^\infty(K)} \omega h_K \|\nabla p_h\|_K^2 \end{aligned} \quad (4.11) $$

<footnote n="4">⁴When $b_0 \cdot n \ne 0$ we can use the standard Poincaré inequality instead.-</footnote>
Combining (4.10) and (4.11) we have

$$ \|\nabla p_h\|_K^2 \le C_{\mathbf{W}_0, K}^2 [C^2 \|\rho_0\|_{L^\infty(K)} \omega^2 h_K^2 \|\nabla p_h\|_K^2 + 2C \|\mathbf{b}_0\|_{L^\infty(K)} \omega h_K \|\nabla p_h\|_K^2], $$

as we assumed $(\mathbf{q}_h, p_h) \ne (0, 0)$ we can divide by $\|\nabla p_h\|_K$ to obtain

$$ 1 \le C_{\mathbf{W}_0, K}^2 [C^2 \|\rho_0\|_{L^\infty(K)} \omega^2 h_K^2 + 2C \|\mathbf{b}_0\|_{L^\infty(K)} \omega h_K]. \quad (4.12) $$

We now define the function

$$ f: \alpha \mapsto C_{\mathbf{W}_0, K}^2 C^2 \|\rho_0\|_{L^\infty(K)} \alpha^2 + 2C_{\mathbf{W}_0, K}^2 C \|\mathbf{b}_0\|_{L^\infty(K)} \alpha - 1 $$

Rewriting (4.12) in terms of $f$ gives

$$ f(\omega h_K) \ge 0. $$

We notice that $f$ is a second-order polynomial whose roots are

$$ \alpha_{\pm} = \frac{-C_{\mathbf{W}_0, K} \|\mathbf{b}_0\|_{L^{\infty}(K)} \pm \left(C_{\mathbf{W}_0, K}^2 \|\mathbf{b}_0\|_{L^{\infty}(K)}^2 + \|\rho_0\|_{L^{\infty}(K)}\right)^{\frac{1}{2}}}{C_{\mathbf{W}_0, K} C \|\rho_0\|_{L^{\infty}(K)}} $$

As the leading coefficient of $f$ is positive, we know that

$$ \forall \alpha \in (\alpha_{-}, \alpha_{+}), \quad f(\alpha) < 0 $$

and it is obvious that $\alpha_{-} < 0$ and $\alpha_{+} > 0$.

Finally, we can see that the assumption on $\omega h_K$ (4.2) is exactly

$$ 0 < \omega h_K < \alpha_+ $$

which means

$$ f(\omega h_K) < 0. $$

This is the desired contradiction and concludes the proof, as we necessarily have $p_h \equiv 0$ and $\mathbf{q}_h \equiv \mathbf{0}$. $\square$

**Remark 4.1.** For triangular elements, the constant $C$ satisfies

$$ C < \frac{1}{\pi}. $$

**Remark 4.2.** when $\mathbf{b}_0 = \mathbf{0}$, the solvability assumption (4.1) becomes

$$ \omega h_K < \frac{1}{CC_{\mathbf{W}_0, K} \|\rho_0\|_{L^\infty(K)}^{\frac{1}{2}}} $$

which is similar to the ones given in [DS19, Prop. 3.9] and [Hun19, Prop. 3.4.2].

**Remark 4.3.** this proof is written for the HDG+ method, for the more standard HDG method only minor changes are needed : in Step 2, $P_M p$ should be replaced by $p$. Assumption (4.1) can therefore be replaced with

$$ \forall e \in \mathcal{F}(K), \quad \tau|_e < 0 \quad \text{or} \quad \tau|_e > \max_e \left(\frac{1}{2} |\mathbf{b}_0 \cdot \mathbf{n}|\right). $$
## 5. ERROR ANALYSIS OF THE HDG+ METHOD

In this section we will carry out a detailed error analysis of the HDG+ method.

The analysis of this section is carried out by mixing elements coming from the analysis of HDG discretization of the Helmholtz equation carried out in [DS19] and of the analysis of HDG+ discretization of elliptic convection-diffusion problems carried out in [QS16a].

We chose to use the orthogonal $L^2$ projections instead of the tailored HDG(+) projections that fit the numerical trace. As we study problems involving convection, the design of a new projection would be required as using the standard HDG(+) projection would not lead to cleaner error system. The design of such a projection seems very difficult when $b_0$ is not constant. This analysis is therefore valid for general polyhedral meshes, and not only simplicial ones.

We denote by $\pi_V$, $\pi_W$ and $P_M$ the $L^2$-orthogonal projections onto $V_h$, $W_h$ and $M_h$ respectively.

We recall the following estimates due to standard approximation theory for polynomials and trace inequalities which will be useful for our analysis, see eg. [EG04, Prop. 1.135]:

$$
\begin{align}
\|p - \pi_W p\|_{\mathcal{O}} &\lesssim h^s \|p\|_{s,\mathcal{O}}, & 0 \le s \le k+2, \nonumber \\
\|q - \pi_V q\|_{\mathcal{O}} &\lesssim h^t \|q\|_{t,\mathcal{O}}, & 0 \le t \le k+1, \nonumber \\
\|p - P_M p\|_{\partial T_h} &\lesssim h^{s-\frac{1}{2}} \|p\|_{s,\mathcal{O}}, & 1 \le s \le k+1, \tag{5.1a} \\
\|p - \pi_W p\|_{\partial K} &\lesssim h^{s-\frac{1}{2}} \|p\|_{s,K}, & 1 \le s \le k+2, \nonumber \\
\|q \cdot n - \pi_V q \cdot n\|_{\partial K} &\lesssim h^{t-\frac{1}{2}} \|q\|_{t,K}, & 1 \le t \le k+1, \nonumber
\end{align}
$$

where $a \lesssim b$ means that there exists a constant $C > 0$ independent of the mesh size and frequency such that $a \le Cb$, we also write $a \approx b$ when $a \lesssim b \lesssim a$.

We will also frequently use the following inverse inequality

$$ \|w\|_{\partial K} \lesssim h^{-\frac{1}{2}} \|w\|_{K}, \quad \forall w \in W_h. \qquad (5.2) $$

Contrary to the case of elliptic problems studied in [QS16a], the convergence of the HDG(+) methods for the convected Helmholtz equation cannot be directly obtained in an energy norm. The analysis of this paper therefore strongly relies on the *Aubin-Nitsche method* to obtain convergence directly in $L^2$-norm. For an introduction to this method, we refer to [EG04, Sec. 2.3.4], similar processes for HDG(+) methods in the context of wave equations have been carried out in [QSS16], [Hun19, Sec. 3.5] and [DS19, Sec. 3.5.2] and for coercive problems with convection in [QS16a]. We introduce the following auxiliary problem

$$
\begin{align}
W_0 \xi - \nabla \theta &= 0, & \text{in } \mathcal{O}, \tag{5.3a} \\
-\rho_0 \omega^2 \theta - 2i\omega b_0 \cdot \nabla \theta - \operatorname{div}(\xi) &= \varepsilon_h^p, & \text{in } \mathcal{O}, \tag{5.3b} \\
\xi \cdot n - 2i\omega(b_0 \cdot n)\theta &= 0, & \text{on } \partial \mathcal{O}. \tag{5.3c}
\end{align}
$$

To carry out the error analysis, we assume that the problem (5.3a)-(5.3b)-(5.3c) is well-posed and satisfies the elliptic regularity assumption

$$ \|\theta\|_{2,\mathcal{O}} + \|\xi\|_{1,\mathcal{O}} \le C_{\text{reg}} \|\varepsilon_h^p\|_{\mathcal{O}}. \qquad (5.4) $$

*Remark 5.1.* As the problem (5.3a)-(5.3b)-(5.3c) can be equivalently written as

$$ -\omega^2 \rho_0 \theta - 2i\omega b_0 \cdot \nabla \theta - \operatorname{div}(K_0 \nabla \theta) = \varepsilon_h^p, $$
it has the same coercive + compact structure as the original problem (1.3) if the background flow is subsonic (see (1.4)). This problem is therefore of Fredholm type, and well-posedness is equivalent to uniqueness of the solution. The well-posedness hypothesis thus means that we will consider frequencies $\omega$ that are not resonant.

We can now state the assumptions under which the error analysis of the HDG+ method is carried out.

**Assumption 1.** The background flow is subsonic.

**Assumption 2.** The frequency $\omega$ is not a resonant frequency of the continuous problem.

**Assumption 3.** The maximum diameter $h$ of the elements is sufficiently small.

**Assumption 4.** The domain $O$ is either a regular bounded open set or a convex polyhedron.

**Assumption 5.** The physical parameters $W_0$ and $b_0$ are Lipschitz continuous in $\bar{O}$.

**ASSUMPTION 1** and **ASSUMPTION 2** ensure that the continuous problem is well-posed. **ASSUMPTION 3** ensures that the local problems are well-posed through **THEOREM 1**. **ASSUMPTION 1**, **ASSUMPTION 2**, and **ASSUMPTION 3** imply that the global problem is well-posed as it will be detailed in **THEOREM 3**. **ASSUMPTION 4** and **ASSUMPTION 5** ensure that the elliptic regularity assumption (5.4) holds as $\varepsilon_h^p \in W_h(T_h) \subset L^2(O)$, on this topic see [Grill, Th. 2.2.2.3 & 4.3.1.4] for the regularity of the domain or [BC13, Sec. 7] for the regularity of the coefficients. Finally the amount of regularity of the exact solution $(q, p)$ is only limited by the regularity of the source term $s$.

To study the convergence of the method, we split the errors as

$$q - q_h = q - \pi_V q + \pi_V q - q_h = \varepsilon_h^q - \delta_h^q,$$

$$p - p_h = p - \pi_W p + \pi_W p - p_h = \varepsilon_h^p - \delta_h^p,$$

$$p - \widehat{p_h} = p - P_M p + P_M p - \widehat{p_h} = \widehat{\varepsilon_h^p} + \widehat{\delta_h^p},$$

where

$$\delta_h^q := \pi_V q - q; \quad \delta_h^p := \pi_W p - p; \quad \widehat{\delta_h^p} := p - P_M p$$

and

$$\varepsilon_h^q := \pi_V q - q_h \in V_h; \quad \varepsilon_h^p := \pi_W p - p_h \in W_h; \quad \widehat{\varepsilon_h^p} := P_M p - \widehat{p_h} \in M_h.$$

For the sake of conciseness, we define

$$\mathcal{E}_h^2 := \|\varepsilon_h^q\|_{\mathcal{T}_h}^2 + \left\|\tau\right\|^2 \left\|\left(P_M \varepsilon_h^p - \widehat{\varepsilon}_h^p\right)\right\|_{\partial \mathcal{T}_h}^2,$$

which should be interpreted as an error estimate in the energy semi-norm. Notice that thank to the Young's inequality, we have

$$\mathcal{E}_h \approx \|\varepsilon_h^q\|_{\mathcal{T}_h} + \left\|\tau\right\|^2 \left\|\left(P_M \varepsilon_h^p - \widehat{\varepsilon}_h^p\right)\right\|_{\partial \mathcal{T}_h}.$$

We also introduce the consistency error

$$\Delta_h := \|\delta_h^q\|_{\mathcal{T}_h} + h^{-1} \|\delta_h^p\|_{\mathcal{T}_h} + h^{\frac{1}{2}} \|\widehat{\delta}_h^p\|_{\partial \mathcal{T}_h}.$$

The approximation errors $\delta_h^q, \delta_h^p, \widehat{\delta}_h^p$ can be estimated using the approximation bounds in (5.1). In particular, if the exact solution $(q, p)$ is sufficiently smooth, we have

$$\Delta_h = \mathcal{O}(h^{k+1}).$$

We therefore only need to prove the following lemma to obtain a convergence result.
**Lemma 5.1.** For $k \ge 1$, and under ASSUMPTION 1, ASSUMPTION 2, ASSUMPTION 3, ASSUMPTION 4 and ASSUMPTION 5, the optimal error estimates hold

$$ \|\boldsymbol{\varepsilon}_h^{\boldsymbol{q}}\|_{\mathcal{T}_h} = \mathcal{O}(h^{k+1}), \quad \text{and} \quad \|\boldsymbol{\varepsilon}_h^{\boldsymbol{p}}\|_{\mathcal{T}_h} = \mathcal{O}(h^{k+2}). \qquad (5.5) $$

We can now state the main result of this section, which is a direct consequence of (5.1) and of LEMMA 5.1.

**Theorem 2** (Convergence of the HDG+ method). For $k \ge 1$, and under ASSUMPTION 1, ASSUMPTION 2, ASSUMPTION 3, ASSUMPTION 4 and ASSUMPTION 5, the HDG+ method (3.7a)-(3.7b)-(3.7c) converges with optimal rates, i.e.

$$ \|p - p_h\|_{\mathcal{T}_h} = \mathcal{O}(h^{k+2}), \quad \text{and} \quad \|\boldsymbol{q} - \boldsymbol{q}_h\|_{\mathcal{T}_h} = \mathcal{O}(h^{k+1}). $$

The end of this section is devoted to the proof LEMMA 5.1 whose sketch is given below.

**Proof of Lemma 5.1:**

The proof can be decomposed into two main steps:

(i) We derive an energy-like estimate of the form

$$ \mathcal{E}_h \lesssim \|\boldsymbol{\varepsilon}_h^{\boldsymbol{p}}\|_{\mathcal{T}_h} + \Delta_h. \qquad (5.6a) $$

(ii) We use the Aubin-Nitsche method to derive an estimate of the form

$$ \|\boldsymbol{\varepsilon}_h^{\boldsymbol{p}}\|_{\mathcal{T}_h} \lesssim h(\mathcal{E}_h + \Delta_h). \qquad (5.6b) $$

The estimates (5.5) are an immediate consequence of (5.6a)-(5.6b). Indeed, by inserting (5.6b) into (5.6a) and using the absorption technique of LEMMA 5.6, we obtain

$$ \mathcal{E}_h \lesssim \Delta_h. $$

Then, using this inequality into (5.6b) and LEMMA 5.6 again, we obtain

$$ \|\boldsymbol{\varepsilon}_h^{\boldsymbol{p}}\|_{\mathcal{T}_h} \lesssim h\Delta_h. $$

Finally, the approximation estimates (5.1) yield

$$ \Delta_h = \mathcal{O}(h^{k+1}), $$

when the exact solution $(\boldsymbol{q}, p)$ is sufficiently regular. This concludes the proof of LEMMA 5.1.

Remark 5.2. We would like to point out that the separation of this proof into two steps is specific of Helmholtz-like problems. Indeed, for coercive problems, the energy estimate (5.6a) is enough to establish the convergence of the method.

## 5.1. Preliminary results

Before actually working on Steps (i) and (ii), we first need to prove some technical results.

### 5.1.1. Error equations

We first establish the error equations satisfied by $\boldsymbol{\varepsilon}_h^{\boldsymbol{q}}$, $\boldsymbol{\varepsilon}_h^{\boldsymbol{p}}$ and $\hat{\boldsymbol{\varepsilon}}_h^{\boldsymbol{p}}$.
**Lemma 5.2.** The error quantities $(\boldsymbol{\varepsilon}_h^{\boldsymbol{q}}, \varepsilon_h^p, \hat{\varepsilon}_h^p) \in \mathbf{V}_h \times \mathbf{W}_h \times \mathbf{M}_h$ satisfy the following error equations:

$$ (\mathbf{W}_0 \boldsymbol{\varepsilon}_h^{\boldsymbol{q}}, \mathbf{r}_h)_{\mathcal{T}_h} - (\varepsilon_h^p, \operatorname{div}(\mathbf{r}_h))_{\mathcal{T}_h} + \langle \hat{\varepsilon}_h^p, \mathbf{r}_h \cdot \mathbf{n} \rangle_{\partial \mathcal{T}_h} = (\mathbf{W}_0 \boldsymbol{\delta}_h^{\boldsymbol{q}}, \mathbf{r}_h)_{\mathcal{T}_h} \qquad (5.7a) $$

$$ (\operatorname{div}(\boldsymbol{\varepsilon}_h^{\boldsymbol{q}}), w_h)_{\mathcal{T}_h} + 2i\omega \langle \tau(P_M \boldsymbol{\varepsilon}_h^p - \hat{\boldsymbol{\varepsilon}}_h^p), w_h \rangle_{\partial \mathcal{T}_h} = \omega^2 (\rho_0 \boldsymbol{\varepsilon}_h^p, w_h)_{\mathcal{T}_h} - 2i\omega (\boldsymbol{\varepsilon}_h^p \mathbf{b}_0, \nabla w_h)_{\mathcal{T}_h} + \langle \hat{\mathcal{Q}}_1 + \hat{\mathcal{Q}}_2 + \hat{\mathcal{Q}}_3, w_h \rangle_{\partial \mathcal{T}_h} $$

$$ - \omega^2 (\rho_0 \delta_h^p, w_h)_{\mathcal{T}_h} + 2i\omega (\delta_h^p \mathbf{b}_0, \nabla w_h)_{\mathcal{T}_h} \qquad (5.7b) $$

$$ - \langle \boldsymbol{\varepsilon}_h^{\boldsymbol{q}} \cdot \mathbf{n} + 2i\omega\tau(P_M \boldsymbol{\varepsilon}_h^p - \hat{\boldsymbol{\varepsilon}}_h^p), \mu_h \rangle_{\partial \mathcal{T}_h} = - \langle \hat{\mathcal{Q}}_1 + \hat{\mathcal{Q}}_2 + \hat{\mathcal{Q}}_3, \mu_h \rangle_{\partial \mathcal{T}_h} \qquad (5.7c) $$

for all $(\mathbf{r}_h, w_h, \mu_h) \in \mathbf{V}_h \times \mathbf{W}_h \times \mathbf{M}_h$, where

$$ \hat{\mathcal{Q}}_1 := 2i\omega(\mathbf{b}_0 \cdot \mathbf{n})\hat{\boldsymbol{\varepsilon}}_h^p + 2i\omega\tau_{upw}(\boldsymbol{\varepsilon}_h^p - \hat{\boldsymbol{\varepsilon}}_h^p), $$

$$ \hat{\mathcal{Q}}_2 := 2i\omega(\mathbf{b}_0 \cdot \mathbf{n})\hat{\delta}_h^p - 2i\omega\tau_{upw}(\delta_h^p - \hat{\delta}_h^p), $$

$$ \hat{\mathcal{Q}}_3 := \boldsymbol{\delta}_h^{\boldsymbol{q}} \cdot \mathbf{n} + 2i\omega\tau P_M \delta_h^p. $$

**Proof:** Notice that $(\mathbf{q}, p, \hat{p})$, where $\hat{p}$ is the trace of $p$ on the skeleton of the mesh, satisfy the equations (3.7a)–(3.7b)–(3.7c) for discrete test-functions, introduce the projections wherever possible and subtract the actual discrete equations. ■

### 5.1.2. An estimate for the gradient of the error

We use the following estimate, obtained following [QS16a, Lemma 3.2], for $\|\nabla \boldsymbol{\varepsilon}_h^p\|_{\partial \mathcal{T}_h}$ to carry out our analysis

**Lemma 5.3.** The following estimate holds

$$ \|\nabla \boldsymbol{\varepsilon}_h^p\|_{\mathcal{T}_h} \lesssim \mathcal{E}_h + \Delta_h. $$

**Proof:** On an element $K \in \mathcal{T}_h$, we notice that $\boldsymbol{\varepsilon}_h^p \in W_h(K)$, which implies $\nabla \boldsymbol{\varepsilon}_h^p \in \mathbf{V}_h(K)$ and $\nabla \boldsymbol{\varepsilon}_h^p \cdot \mathbf{n} \in \mathbf{M}_h(K)$. We can therefore go back to (5.7a), test with $\mathbf{r}_h = \nabla \boldsymbol{\varepsilon}_h^p$ and integrate by parts to obtain

$$ \|\nabla \boldsymbol{\varepsilon}_h^p\|_K^2 = (\mathbf{W}_0 \boldsymbol{\delta}_h^{\boldsymbol{q}}, \nabla \boldsymbol{\varepsilon}_h^p)_K - (\mathbf{W}_0 \boldsymbol{\varepsilon}_h^{\boldsymbol{q}}, \nabla \boldsymbol{\varepsilon}_h^p)_K + \langle \boldsymbol{\varepsilon}_h^p - \hat{\boldsymbol{\varepsilon}}_h^p, \nabla \boldsymbol{\varepsilon}_h^p \cdot \mathbf{n} \rangle_{\partial K} \qquad (5.8) $$

Noticing that $\nabla \boldsymbol{\varepsilon}_h^p \cdot \mathbf{n} \in \mathbf{M}_h(K)$ and resorting to (4.9), we have thanks to the Cauchy-Schwarz inequality

$$ |\langle \boldsymbol{\varepsilon}_h^p - \hat{\boldsymbol{\varepsilon}}_h^p, \nabla \boldsymbol{\varepsilon}_h^p \cdot \mathbf{n} \rangle_{\partial K}| = |\langle P_M \boldsymbol{\varepsilon}_h^p - \hat{\boldsymbol{\varepsilon}}_h^p, \nabla \boldsymbol{\varepsilon}_h^p \cdot \mathbf{n} \rangle_{\partial K}| \lesssim \|P_M \boldsymbol{\varepsilon}_h^p - \hat{\boldsymbol{\varepsilon}}_h^p\|_{\partial K} \|\nabla \boldsymbol{\varepsilon}_h^p\|_{\partial K}. $$

Since $\nabla \boldsymbol{\varepsilon}_h^p \in \mathbf{V}_h(K)$, the trace inequality (5.2) leads to

$$ |\langle \boldsymbol{\varepsilon}_h^p - \hat{\boldsymbol{\varepsilon}}_h^p, \nabla \boldsymbol{\varepsilon}_h^p \cdot \mathbf{n} \rangle_{\partial K}| \lesssim \|\hat{\boldsymbol{\varepsilon}}_h^p - P_M \boldsymbol{\varepsilon}_h^p\|_{\partial K} h^{-\frac{1}{2}} \|\nabla \boldsymbol{\varepsilon}_h^p\|_K $$

Incorporating $\tau|_K = \mathcal{O}(h_K^{-1})$, $h^{-\frac{1}{2}}$ in the boundary norm yields

$$ |\langle \boldsymbol{\varepsilon}_h^p - \hat{\boldsymbol{\varepsilon}}_h^p, \nabla \boldsymbol{\varepsilon}_h^p \cdot \mathbf{n} \rangle_{\partial K}| \lesssim \|\tau|^{\frac{1}{2}}(\hat{\boldsymbol{\varepsilon}}_h^p - P_M \boldsymbol{\varepsilon}_h^p)\|_{\partial K} \|\nabla \boldsymbol{\varepsilon}_h^p\|_K \qquad (5.9) $$

We can now use (5.9) into (5.8) and apply the Cauchy-Schwarz inequality to obtain

$$ \|\nabla \boldsymbol{\varepsilon}_h^p\|_K^2 \lesssim (\|\boldsymbol{\varepsilon}_h^{\boldsymbol{q}}\|_{\mathbf{W}_0,K} + \|\tau|^{\frac{1}{2}}(\hat{\boldsymbol{\varepsilon}}_h^p - P_M \boldsymbol{\varepsilon}_h^p)\|_{\partial K} + \|\boldsymbol{\delta}_h^{\boldsymbol{q}}\|_{\mathbf{W}_0,K}) \|\nabla \boldsymbol{\varepsilon}_h^p\|_K, $$
This allows to conclude by summing over the elements and using that

$$ \|\varepsilon_h^q\|_{\mathcal{T}_h} + \left\|\left|\tau\right|^{\frac{1}{2}}(P_M \varepsilon_h^p - \tilde{\varepsilon}_h^p)\right\|_{\partial \mathcal{T}_h} \lesssim \mathcal{E}_h, \quad \text{and} \quad \|\delta_h^q\|_{\mathcal{T}_h} \lesssim \Delta_h. $$

### 5.1.3. Using the Poincaré-Wirtinger inequality

In the error analysis, we will apply the Poincaré-Wirtinger inequality. We denote by {$u$} the $L^2$-projection of $u$ on $\mathcal{P}_0$, ie

$$ \forall K \in \mathcal{T}_h, \quad \{u\}|_K = \frac{1}{|K|} \int_K u d\boldsymbol{x} $$

We will need to subtract {$u$} from the equations to apply the Poincaré-Wirtinger inequality :

$$ \|u - \{u\}\|_{\mathcal{T}_h} \le Ch \|\nabla u\|_{\mathcal{T}_h}. \qquad (5.10) $$

We can do that thanks to the following property of the projections $\pi_W$ and $\pi_V$, indeed we have for $\boldsymbol{\xi}$ and $\boldsymbol{q}$

$$ (\pi_V \boldsymbol{q}, \{\boldsymbol{W}_0 \boldsymbol{\xi}\})_{\mathcal{T}_h} = (\boldsymbol{q}, \{\boldsymbol{W}_0 \boldsymbol{\xi}\})_{\mathcal{T}_h} \quad \text{as} \quad \{\boldsymbol{W}_0 \boldsymbol{\xi}\} \in \mathcal{P}_0 \subset \mathcal{P}_k $$

therefore

$$ (\delta_h^q, \boldsymbol{W}_0 \boldsymbol{\xi})_{\mathcal{T}_h} = (\boldsymbol{q} - \pi_V \boldsymbol{q}, \boldsymbol{W}_0 \boldsymbol{\xi})_{\mathcal{T}_h} = (\boldsymbol{q} - \pi_V \boldsymbol{q}, \boldsymbol{W}_0 \boldsymbol{\xi} - \{\boldsymbol{W}_0 \boldsymbol{\xi}\})_{\mathcal{T}_h}. $$

Similar results can be obtained in the same way for the other quantities.

### 5.1.4. Best approximation property of $P_M$

During the analysis, we will often need to compare quantities like $\|u - P_M u\|_{\partial K}$ and $\|u - \{u\}\|_{\partial K}$.

**Lemma 5.4.** For $u \in \mathcal{P}_{k+1}(\mathcal{T}_h)$, the following inequality holds

$$ \|u - P_M u\|_{\partial K} \lesssim \|u - \{u\}\|_{\partial K} $$

**Proof:**

We recall that $M_h := \prod_{e \in \mathcal{F}_h} \mathcal{P}_k(e)$ is a finite-dimensional vector subspace of $L^2(\partial \mathcal{T}_h)$. We recalled that functions in $M_h$ are piecewise polynomials of degree up to $k$ on the skeleton of the mesh. On an internal edge $e = \partial K_- \cap \partial K_+$, we define

$$ \{u\}_e := \begin{cases} \{u^-\} & \text{on } K_-, \\ \{u^+\} & \text{on } K_+, \end{cases} \quad \text{where } u^\pm = u|_{K^\pm}. $$

With this definition {$u$}_e is a bi-valued piecewise constant on the skeleton of the mesh, and therefore {$u$}_e $\in M_h$. As $P_M$ is the orthogonal projection onto $M_h$, we can use the Hilbert projection theorem to obtain

$$ \|u - P_M u\|_{\partial \mathcal{T}_h} \le \inf_{v \in M_h} \|u - v\|_{\partial \mathcal{T}_h}. $$

We can therefore conclude that

$$ \|u - P_M u\|_{\partial \mathcal{T}_h} \lesssim \|u - \{u\}_e\|_{\partial \mathcal{T}_h}. $$

When no confusions are possible, we will denote {$u$}_e by {$u$}.

This property will be referred to as the best approximation property of $P_M$.
### 5.1.5. A jump estimate

During the analysis, the quantity $||\epsilon_h^p - \hat{\epsilon}_h^p||_{\partial T_h}$ will be encountered. For standard HDG methods, this quantity is incorporated into the energy (semi)norm. For the particular case of HDG+ however, the natural jump term is $||\tau|^{\frac{1}{2}}(P_M \epsilon_h^p - \hat{\epsilon}_h^p)||_{\partial T_h}$ which is smaller than $\mathcal{E}_h$. Those two terms are linked by the following lemma.

**Lemma 5.5 (Jump estimate).** *The following inequality holds*

$$ ||\epsilon_h^p - \hat{\epsilon}_h^p||_{\partial T_h} \lesssim h^{\frac{1}{2}} (\mathcal{E}_h + \Delta_h). $$

**Proof:** We notice that

$$ ||\epsilon_h^p - \hat{\epsilon}_h^p||_{\partial T_h} \lesssim ||\epsilon_h^p - P_M \epsilon_h^p||_{\partial T_h} + ||P_M \epsilon_h^p - \hat{\epsilon}_h^p||_{\partial T_h} $$

we use the LEMMA 5.4 and we remark that $\tau = \mathcal{O}(h^{-1})$ to obtain

$$ \lesssim ||\epsilon_h^p - \{\epsilon_h^p\}||_{\partial T_h} + h^{\frac{1}{2}} ||\tau|^{\frac{1}{2}}(P_M \epsilon_h^p - \hat{\epsilon}_h^p)||_{\partial T_h}, $$

then, by (5.2), we also have

$$ \lesssim h^{-\frac{1}{2}} ||\epsilon_h^p - \{\epsilon_h^p\}||_{\partial T_h} + h^{\frac{1}{2}} ||\tau|^{\frac{1}{2}}(P_M \epsilon_h^p - \hat{\epsilon}_h^p)||_{\partial T_h}, $$

we then use (5.10) to obtain

$$ \lesssim h^{\frac{1}{2}} ||\nabla \epsilon_h^p||_{\partial T_h} + h^{\frac{1}{2}} ||\tau|^{\frac{1}{2}}(P_M \epsilon_h^p - \hat{\epsilon}_h^p)||_{\partial T_h}, $$

and the gradient estimate of LEMMA 5.3 leads to

$$ ||\epsilon_h^p - \hat{\epsilon}_h^p||_{\partial T_h} \lesssim h^{\frac{1}{2}} (\mathcal{E}_h + \Delta_h). $$

### 5.1.6. Absorption lemma

The subsequent error analysis heavily relies on the absorption technique which is summarised in the following lemma.

**Lemma 5.6 (Absorption).** Let a, b, h be positive real numbers and c be a real number. The two following propositions hold:

(1) if $a \lesssim ah + b$ and if $h$ is sufficiently small, then $a \lesssim b$,

(2) if $a^2 \lesssim ab + c$ then $a^2 \lesssim b^2 + c$.

## 5.2. Step (i): Energy-like estimate

We can now prove the energy-like estimate (5.6a).

Test (5.7a)-(5.7b)-(5.7c) with $(\epsilon_h^q, \epsilon_h^p, \hat{\epsilon}_h^p)$, conjugate (5.7a), and sum the resulting equations. The left-hand side of the sum is

$$ \begin{aligned} & ||\epsilon_h^q||_{W_0, T_h}^2 - (\text{div}(\epsilon_h^q), \epsilon_h^p)_{T_h} + \langle \epsilon_h^q \cdot n, \hat{\epsilon}_h^p \rangle_{\partial T_h} + (\text{div}(\epsilon_h^q), \epsilon_h^p)_{T_h} + 2i\omega \langle \tau(P_M \epsilon_h^p - \hat{\epsilon}_h^p), P_M \epsilon_h^p \rangle_{\partial T_h} \\ & - \langle \epsilon_h^q \cdot n + 2i\omega\tau(P_M \epsilon_h^p - \hat{\epsilon}_h^p), \hat{\epsilon}_h^p \rangle_{\partial T_h} = ||\epsilon_h^q||_{W_0, T_h}^2 + 2i\omega \langle \tau(P_M \epsilon_h^p - \hat{\epsilon}_h^p), P_M \epsilon_h^p - \hat{\epsilon}_h^p \rangle_{\partial T_h}, \\ & = ||\epsilon_h^q||_{W_0, T_h}^2 - 2i\omega ||\tau|^{\frac{1}{2}}(P_M \epsilon_h^p - \hat{\epsilon}_h^p)||_{\partial T_h}^2, \end{aligned} $$
as $\tau < 0$. Summing the right-hand sides, we obtain the discrete energy-like identity

$$ \|\boldsymbol{\varepsilon}_h^{\boldsymbol{q}}\|_{\boldsymbol{W}_0, \mathcal{T}_h}^2 - 2i\omega \left\|\tau^{\frac{1}{2}}(P_M \varepsilon_h^p - \widehat{\varepsilon}_h^p)\right\|_{\partial \mathcal{T}_h}^2 = \omega^2 \|\varepsilon_h^p\|_{\rho_0, \mathcal{T}_h}^2 + R_1 + R_2 + R_3 + R_4, \quad (5.11) $$

where

$$ R_1 := -2i\omega (\varepsilon_h^p \mathbf{b}_0, \nabla \varepsilon_h^p)_{\mathcal{T}_h} + \langle \hat{\mathcal{Q}}_1, \varepsilon_h^p - \widehat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h}, $$

$$ R_2 := \langle \hat{\mathcal{Q}}_2, \varepsilon_h^p - \widehat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h}, $$

$$ R_3 := \langle \hat{\mathcal{Q}}_3, \varepsilon_h^p - \widehat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h}, $$

$$ R_4 := -\omega^2 (\rho_0 \delta_h^p, \varepsilon_h^p)_{\mathcal{T}_h} + 2i\omega (\delta_h^p \mathbf{b}_0, \nabla \varepsilon_h^p)_{\mathcal{T}_h} + (\varepsilon_h^q, \boldsymbol{W}_0 \delta_h^q)_{\mathcal{T}_h}. $$

To prove (5.6a), we now need to separately estimate the terms in the right-hand side of (5.11).

### 5.2.1. Estimating $R_1$

We recall that $\hat{\mathcal{Q}}_1$ is defined as

$$ \hat{\mathcal{Q}}_1 := 2i\omega(\mathbf{b}_0 \cdot \mathbf{n})\widehat{\varepsilon}_h^p + 2i\omega\tau_{\text{upw}}(\varepsilon_h^p - \widehat{\varepsilon}_h^p). $$

To estimate this term, it is easier to separate its real and imaginary parts. Using LEMMA 4.1, we have

$$ \mathfrak{Re} (\varepsilon_h^p \mathbf{b}_0, \nabla \varepsilon_h^p)_{\mathcal{T}_h} = \frac{1}{2} \langle (\mathbf{b}_0 \cdot \mathbf{n}) \varepsilon_h^p, \varepsilon_h^p \rangle_{\partial \mathcal{T}_h}, $$

and we therefore have

$$ \begin{aligned} \frac{\mathfrak{Im} R_1}{2\omega} &= \mathfrak{Re} \langle (\mathbf{b}_0 \cdot \mathbf{n}) \widehat{\varepsilon}_h^p, \varepsilon_h^p - \widehat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h} + \langle \tau_{\text{upw}} (\varepsilon_h^p - \widehat{\varepsilon}_h^p), \varepsilon_h^p - \widehat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h} - \frac{1}{2} \langle (\mathbf{b}_0 \cdot \mathbf{n}) \varepsilon_h^p, \varepsilon_h^p \rangle_{\partial \mathcal{T}_h}, \\ &= \left\langle -\frac{1}{2} (\mathbf{b}_0 \cdot \mathbf{n}) (\varepsilon_h^p - \widehat{\varepsilon}_h^p), \varepsilon_h^p - \widehat{\varepsilon}_h^p \right\rangle_{\partial \mathcal{T}_h} - \frac{1}{2} \langle (\mathbf{b}_0 \cdot \mathbf{n}) \widehat{\varepsilon}_h^p, \widehat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h} + \langle \tau_{\text{upw}} (\varepsilon_h^p - \widehat{\varepsilon}_h^p), \varepsilon_h^p - \widehat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h}. \end{aligned} $$

Noticing that $\widehat{\varepsilon}_h^p$ is single-valued across the skeleton of the mesh we have

$$ \langle (\mathbf{b}_0 \cdot \mathbf{n}) \widehat{\varepsilon}_h^p, \widehat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h} = \langle [\widehat{\varepsilon}_h^p \mathbf{b}_0], \widehat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h} = 0, $$

We notice that $\tau_{\text{upw}} - \frac{1}{2}(\mathbf{b}_0 \cdot \mathbf{n}) = \frac{1}{2}|\mathbf{b}_0 \cdot \mathbf{n}| \ge 0$ and obtain

$$ \frac{\mathfrak{Im} R_1}{2\omega} = \left\langle \left( \tau_{\text{upw}} - \frac{1}{2}(\mathbf{b}_0 \cdot \mathbf{n}) \right) (\varepsilon_h^p - \widehat{\varepsilon}_h^p), \varepsilon_h^p - \widehat{\varepsilon}_h^p \right\rangle_{\partial \mathcal{T}_h} = \left\| \left( \frac{1}{2}|\mathbf{b}_0 \cdot \mathbf{n}| \right)^{\frac{1}{2}} (\varepsilon_h^p - \widehat{\varepsilon}_h^p) \right\|_{\partial \mathcal{T}_h}^2, $$

On the other hand, we also have

$$ \frac{\mathfrak{Re} R_1}{2\omega} = \mathfrak{Im} \langle (\mathbf{b}_0 \cdot \mathbf{n}) \widehat{\varepsilon}_h^p, \varepsilon_h^p \rangle_{\partial \mathcal{T}_h} - \mathfrak{Im} (\varepsilon_h^p \mathbf{b}_0, \nabla \varepsilon_h^p)_{\mathcal{T}_h}, $$

noticing that $\langle (\mathbf{b}_0 \cdot \mathbf{n}) \varepsilon_h^p, \varepsilon_h^p \rangle_{\partial \mathcal{T}_h} \in \mathbb{R}$, we have

$$ \frac{\mathfrak{Re} R_1}{2\omega} = \mathfrak{Im} \langle (\mathbf{b}_0 \cdot \mathbf{n}) (\varepsilon_h^p - \widehat{\varepsilon}_h^p), \varepsilon_h^p \rangle_{\partial \mathcal{T}_h} - \mathfrak{Im} (\varepsilon_h^p \mathbf{b}_0, \nabla \varepsilon_h^p)_{\mathcal{T}_h}. $$
We therefore obtain that

$$|R_1| \lesssim \left\| \left( \frac{1}{2} |\mathbf{b}_0 \cdot \mathbf{n}| \right)^{\frac{1}{2}} (\varepsilon_h^p - \hat{\varepsilon}_h^p) \right\|_{\partial \mathcal{T}_h}^2 + \left| \langle (\mathbf{b}_0 \cdot \mathbf{n})(\hat{\varepsilon}_h^p - \varepsilon_h^p), \varepsilon_h^p \rangle_{\partial \mathcal{T}_h} \right| + \left| (\varepsilon_h^p \mathbf{b}_0, \nabla \varepsilon_h^p)_{\mathcal{T}_h} \right|,$$

Cauchy Schwarz inequality leads to

$$|R_1| \lesssim \|\varepsilon_h^p - \hat{\varepsilon}_h^p\|_{\partial \mathcal{T}_h}^2 + \|\varepsilon_h^p - \hat{\varepsilon}_h^p\|_{\partial \mathcal{T}_h} \|\varepsilon_h^p\|_{\partial \mathcal{T}_h} + \|\varepsilon_h^p\|_{\mathcal{T}_h} \|\nabla \varepsilon_h^p\|_{\mathcal{T}_h}.$$

Using the jump estimate of LEMMA 5.5, the gradient estimate of LEMMA 5.3 and that $\tau = O(h^{-1})$, we obtain

$$|R_1| \lesssim (\mathcal{E}_h + \Delta_h) \|\varepsilon_h^p\|_{\mathcal{T}_h} + h (\mathcal{E}_h^2 + \Delta_h^2). \quad (5.12)$$

### 5.2.2. Estimating $R_2$

We recall that $\hat{Q}_2$ is defined as

$$\hat{Q}_2 := 2i\omega(\mathbf{b}_0 \cdot \mathbf{n})\hat{\delta}_h^p - 2i\omega\tau_{\text{upw}}(\delta_h^p - \hat{\delta}_h^p).$$

Cauchy-Schwarz inequality leads to

$$\left| \left\langle (\mathbf{b}_0 \cdot \mathbf{n}) \hat{\delta}_h^p, \varepsilon_h^p - \hat{\varepsilon}_h^p \right\rangle_{\partial \mathcal{T}_h} \right| \lesssim \left\| \hat{\delta}_h^p \right\|_{\partial \mathcal{T}_h} \|\varepsilon_h^p - \hat{\varepsilon}_h^p\|_{\partial \mathcal{T}_h},$$

Using the trace inequality (5.2) and the jump estimate of LEMMA 5.5, we have

$$\left| \left\langle (\mathbf{b}_0 \cdot \mathbf{n}) \hat{\delta}_h^p, \varepsilon_h^p - \hat{\varepsilon}_h^p \right\rangle_{\partial \mathcal{T}_h} \right| \lesssim h^{\frac{1}{2}} \left\| \hat{\delta}_h^p \right\|_{\mathcal{T}_h} (\mathcal{E}_h + \Delta_h) \lesssim \Delta_h (\mathcal{E}_h + \Delta_h),$$

$$\text{as } h^{\frac{1}{2}} \left\| \hat{\delta}_h^p \right\|_{\partial \mathcal{T}_h} \lesssim \Delta_h.$$

Using the trace inequality (5.2), we have

$$\left\| \delta_h^p - \hat{\delta}_h^p \right\|_{\partial \mathcal{T}_h} \lesssim h^{-\frac{1}{2}} \| \delta_h^p \|_{\mathcal{T}_h} + \left\| \hat{\delta}_h^p \right\|_{\partial \mathcal{T}_h} \lesssim h^{-\frac{1}{2}} (1+h) \Delta_h \lesssim h^{-\frac{1}{2}} \Delta_h,$$

and the jump estimate of LEMMA 5.5 yields

$$\left| \left\langle \tau_{\text{upw}}(\delta_h^p - \hat{\delta}_h^p), \varepsilon_h^p - \hat{\varepsilon}_h^p \right\rangle_{\partial \mathcal{T}_h} \right| \lesssim h^{-\frac{1}{2}} \Delta_h \|\varepsilon_h^p - \hat{\varepsilon}_h^p\|_{\partial \mathcal{T}_h} \lesssim \Delta_h (\mathcal{E}_h + \Delta_h).$$

So we finally have

$$|R_2| \lesssim \Delta_h (\mathcal{E}_h + \Delta_h). \quad (5.13)$$

### 5.2.3. Estimating $R_3$

We recall that $\hat{Q}_3$ is defined as

$$\hat{Q}_3 := \delta_h^q \cdot \mathbf{n} + 2i\omega\tau P_M \delta_h^p.$$
For the first term, using the jump estimate of LEMMA 5.5 and the trace inequality (5.2), we have

$$ \left| \langle \boldsymbol{\delta}_h^q \cdot \boldsymbol{n}, \varepsilon_h^p - \hat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h} \right| \lesssim h^{\frac{1}{2}} (\mathcal{E}_h + \Delta_h) \| \boldsymbol{\delta}_h^q \|_{\partial \mathcal{T}_h} \lesssim (\mathcal{E}_h + \Delta_h) \Delta_h, $$

as $\|\boldsymbol{\delta}_h^q\|_{\mathcal{T}_h} \lesssim \Delta_h$. For the second term, we have

$$ \begin{aligned} \left| \langle \tau P_M \delta_h^p, \varepsilon_h^p - \hat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h} \right| & \le \left| \langle |\tau| P_M \delta_h^p, P_M \varepsilon_h^p - \hat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h} \right| \\ & \le \left| \left\langle |\tau|^{\frac{1}{2}} \delta_h^p, |\tau|^{\frac{1}{2}} (P_M \varepsilon_h^p - \hat{\varepsilon}_h^p) \right\rangle_{\partial \mathcal{T}_h} \right| \end{aligned} $$

using that $\tau = \mathcal{O}(h^{-1})$, the trace inequality (5.2) and that $h^{-1} \|\delta_h^p\|_{\mathcal{T}_h} \lesssim \Delta_h$, we obtain

$$ \left| \langle \tau P_M \delta_h^p, \varepsilon_h^p - \hat{\varepsilon}_h^p \rangle_{\partial \mathcal{T}_h} \right| \lesssim \left\| |\tau|^{\frac{1}{2}} (P_M \varepsilon_h^p - \hat{\varepsilon}_h^p) \right\|_{\partial \mathcal{T}_h} h^{-\frac{1}{2}} \|\delta_h^p\|_{\partial \mathcal{T}_h} \lesssim h^{-1} \mathcal{E}_h \|\delta_h^p\|_{\mathcal{T}_h} \lesssim \mathcal{E}_h \Delta_h. $$

So we finally have

$$ |R_3| \lesssim \Delta_h (\mathcal{E}_h + \Delta_h). \qquad (5.14) $$

### 5.2.4. Estimating $R_4$

By similar computations using the Cauchy-Schwarz inequality, and the gradient estimate of LEMMA 5.3 we can show that

$$ \begin{aligned} & \left| (\rho_0 \delta_h^p, \varepsilon_h^p)_{\mathcal{T}_h} \right| \lesssim \|\delta_h^p\|_{\mathcal{T}_h} \|\varepsilon_h^p\|_{\mathcal{T}_h} \lesssim h \Delta_h \|\varepsilon_h^p\|_{\mathcal{T}_h}, \\ & \left| (\delta_h^p \mathbf{b}_0, \nabla \varepsilon_h^p)_{\mathcal{T}_h} \right| \lesssim \|\delta_h^p\|_{\mathcal{T}_h} (\mathcal{E}_h + \Delta_h) \lesssim h \Delta_h (\mathcal{E}_h + \Delta_h), \\ & \left| (\mathbf{W}_0 \varepsilon_h^q, \delta_h^q)_{\mathcal{T}_h} \right| \lesssim \|\varepsilon_h^q\|_{\mathbf{W}_0, \mathcal{T}_h} \|\delta_h^q\|_{\mathcal{T}_h} \lesssim \mathcal{E}_h \Delta_h. \end{aligned} $$

We therefore have

$$ |R_4| \lesssim h \Delta_h \left( \|\varepsilon_h^p\|_{\mathcal{T}_h} + \Delta_h \right) + \mathcal{E}_h \Delta_h. \qquad (5.15) $$

### 5.2.5. Discrete energy-like estimate

We first notice that

$$ \mathcal{E}_h^2 := \|\varepsilon_h^q\|_{\mathcal{T}_h}^2 + \| |\tau|^2 (P_M \varepsilon_h^p - \hat{\varepsilon}_h^p) \|_{\partial \mathcal{T}_h}^2 \lesssim \left| \|\varepsilon_h^q\|_{\mathbf{W}_0, \mathcal{T}_h}^2 - 2i\omega \| |\tau|^2 (P_M \varepsilon_h^p - \hat{\varepsilon}_h^p) \|_{\partial \mathcal{T}_h}^2 \right|. $$

By collecting the estimates (5.12), (5.13), (5.14) and (5.15) and by using them into the discrete energy-like identity (5.11), we obtain

$$ \mathcal{E}_h^2 \lesssim \|\varepsilon_h^p\|_{\mathcal{T}_h}^2 + (\mathcal{E}_h + \Delta_h) \|\varepsilon_h^p\|_{\mathcal{T}_h} + h \mathcal{E}_h^2 + \Delta_h \mathcal{E}_h + \Delta_h^2, $$

we then resort to the absorption technique described in LEMMA 5.6 to obtain

$$ \mathcal{E}_h^2 \lesssim \|\varepsilon_h^p\|_{\mathcal{T}_h}^2 + \Delta_h^2 \lesssim (\|\varepsilon_h^p\|_{\mathcal{T}_h} + \Delta_h)^2, $$

which yields

$$ \mathcal{E}_h \lesssim \|\varepsilon_h^p\|_{\mathcal{T}_h} + \Delta_h. $$
This is the desired estimate (5.6a). Moreover, we would like to point out that the consistency term $\Delta_h$ be estimated as $O(h^{k+1})$ by using (5.1) when the exact solution $(q, p)$ is sufficiently regular.

### 5.3. Step (ii): Aubin-Nitsche method

We now prove the Aubin-Nitsche estimate (5.6b). Let $(\xi, \theta) \in H^1(\mathcal{O}) \times H^2(\mathcal{O})$ be the solution to the auxiliary problem (5.3a)-(5.3b)-(5.3c), it satisfies the following weak formulation

$$ (\mathbf{W}_0\xi, \mathbf{r}_h)_{\mathcal{T}_h} + (\theta, \mathrm{div}(\mathbf{r}_h))_{\mathcal{T}_h} - \langle \theta, \mathbf{r}_h \cdot \mathbf{n} \rangle_{\partial \mathcal{T}_h} = 0, \quad (5.16a) $$

$$ -\omega^2 (\rho_0 \theta, w_h)_{\mathcal{T}_h} + (2i\omega\theta\mathbf{b}_0 + \xi, \nabla w_h)_{\mathcal{T}_h} + \langle \xi \cdot \mathbf{n} - 2i\omega(\mathbf{b}_0 \cdot \mathbf{n})\theta, w_h \rangle_{\partial \mathcal{T}_h} = (\varepsilon_h^p, w_h)_{\mathcal{T}_h}, \quad (5.16b) $$

$$ \langle \xi \cdot \mathbf{n}, \mu_h \rangle_{\partial \mathcal{T}_h \setminus \Gamma} = 0, \quad (5.16c) $$

for all discrete test-functions $(\mathbf{r}_h, w_h, \mu_h) \in V_h \times W_h \times M_h$. We recall that according to the assumptions, the solution $(\xi, \theta)$ exists in $H^1(\mathcal{O}) \times H^2(\mathcal{O})$ as it satisfies the elliptic regularity estimate (5.4). The last equation (5.16c) translates the continuity of $\xi \cdot \mathbf{n}$ between the elements and should be interpreted as a jump term. Indeed by the same argument as when we discussed weak continuity of $q_h \cdot \mathbf{n}$ in SUBSECTION 3, we have

$$ 0 = \sum_{K \in \mathcal{T}_h} \int_{\partial K} \xi \cdot n \overline{\mu_h} d\sigma = \sum_{e \in \mathcal{F}_h^i} \int_e \llbracket \xi \rrbracket \overline{\mu_h} d\sigma. $$

#### 5.3.1. Aubin-Nitsche identity

Introducing the projections in (5.16a)-(5.16b)-(5.16c), testing with $(\varepsilon_h^q, \varepsilon_h^p, \tilde{\varepsilon}_h^p)$, conjugating and summing the resulting equations, we obtain

$$ (\mathbf{W}_0 \varepsilon_h^q, \pi_V \xi)_{\mathcal{T}_h} + (\mathrm{div}(\varepsilon_h^p), \pi_W \theta)_{\mathcal{T}_h} - \langle \varepsilon_h^q \cdot \mathbf{n}, P_M \theta \rangle_{\partial \mathcal{T}_h} - \omega^2 (\rho_0 \varepsilon_h^p, \pi_W \theta)_{\mathcal{T}_h} - 2i\omega (\mathbf{b}_0 \cdot \nabla \varepsilon_h^p, \pi_W \theta)_{\mathcal{T}_h} \\ + (\nabla \varepsilon_h^p, \pi_V \xi)_{\mathcal{T}_h} + \langle \varepsilon_h^p, \pi_V \xi \cdot \mathbf{n} + 2i\omega (\mathbf{b}_0 \cdot \mathbf{n}) \pi_W \theta \rangle_{\partial \mathcal{T}_h} - \langle \tilde{\varepsilon}_h^p, \pi_V \xi \cdot \mathbf{n} \rangle_{\partial \mathcal{T}_h} = \ell_1, $$

on the other hand summing (5.7a)-(5.7b)-(5.7c) tested with $(\pi_V \xi, \pi_W \theta, P_M \theta)$ leads to

$$ (\mathbf{W}_0 \varepsilon_h^q, \pi_V \xi)_{\mathcal{T}_h} + (\nabla \varepsilon_h^p, \pi_V \xi)_{\mathcal{T}_h} + \langle \varepsilon_h^p - \tilde{\varepsilon}_h^p, \pi_V \xi \cdot \mathbf{n} \rangle_{\partial \mathcal{T}_h} - \omega^2 (\rho_0 \varepsilon_h^p, \pi_W \theta)_{\mathcal{T}_h} - 2i\omega (\mathbf{b}_0 \cdot \nabla \varepsilon_h^p, \pi_W \theta)_{\mathcal{T}_h} \\ + 2i\omega \langle (\mathbf{b}_0 \cdot \mathbf{n}) \varepsilon_h^p, \pi_W \theta \rangle_{\partial \mathcal{T}_h} + (\mathrm{div}(\varepsilon_h^q), \pi_W \theta)_{\mathcal{T}_h} - \langle \varepsilon_h^q \cdot \mathbf{n}, P_M \theta \rangle_{\partial \mathcal{T}_h} = \ell_2. $$

Hence $\ell_1 = \ell_2$ where

$$ \ell_1 = (\mathbf{W}_0 \varepsilon_h^q, \pi_V \xi - \xi)_{\mathcal{T}_h} + \|\varepsilon_h^p\|_{\rho_0, \mathcal{T}_h}^2 - \omega^2 (\rho_0 \varepsilon_h^p, \pi_W \theta - \theta)_{\mathcal{T}_h} - 2i\omega (\mathbf{b}_0 \cdot \nabla \varepsilon_h^p, \pi_W \theta - \theta)_{\mathcal{T}_h} \\ + \langle \varepsilon_h^p - \tilde{\varepsilon}_h^p, (\pi_V \xi - \xi) \cdot \mathbf{n} \rangle_{\partial \mathcal{T}_h} + 2i\omega \langle (\mathbf{b}_0 \cdot \mathbf{n}) \varepsilon_h^p, \pi_W \theta - \theta \rangle_{\partial \mathcal{T}_h}, $$

$$ \ell_2 = (\mathbf{W}_0 \delta_h^q, \pi_V \xi)_{\mathcal{T}_h} - \omega^2 (\rho_0 \delta_h^p, \pi_W \theta)_{\mathcal{T}_h} + 2i\omega (\delta_h^p \mathbf{b}_0, \nabla \pi_W \theta)_{\mathcal{T}_h} \\ - 2i\omega \langle \tau(P_M \varepsilon_h^p - \tilde{\varepsilon}_h^p), \pi_W \theta - P_M \theta \rangle_{\partial \mathcal{T}_h} + \langle \hat{Q}_1 + \hat{Q}_2 + \hat{Q}_3, \pi_W \theta - P_M \theta \rangle_{\partial \mathcal{T}_h}. $$

It follows that

$$ \|\varepsilon_h^p\|_{\rho_0, \mathcal{T}_h}^2 = m_1 + m_2 + m_3 + m_4, $$
with

$$m_1 = -(\mathbf{W}_0 \boldsymbol{\varepsilon}_h^q, \pi_V \boldsymbol{\xi} - \boldsymbol{\xi})_{\mathcal{T}_h} + (\mathbf{W}_0 \delta_h^q, \pi_V \boldsymbol{\xi})_{\mathcal{T}_h} + \omega^2 (\rho_0 \varepsilon_h^p, \pi_W \theta - \theta)_{\mathcal{T}_h} - \omega^2 (\rho_0 \delta_h^p, \pi_W \theta)_{\mathcal{T}_h}$$

$$m_2 = 2i\omega (\mathbf{b}_0 \cdot \nabla \varepsilon_h^p, \pi_W \theta - \theta)_{\mathcal{T}_h} - 2i\omega \langle (\mathbf{b}_0 \cdot \mathbf{n}) \varepsilon_h^p, \pi_W \theta - \theta \rangle_{\partial \mathcal{T}_h} + 2i\omega (\delta_h^p \mathbf{b}_0, \nabla \pi_W \theta)_{\mathcal{T}_h},$$

$$m_3 = -2i\omega \langle \tau(P_M \varepsilon_h^p - \varepsilon_h^p), \pi_W \theta - P_M \theta \rangle_{\partial \mathcal{T}_h} + \langle \hat{Q}_1 + \hat{Q}_2 + \hat{Q}_3, \pi_W \theta - P_M \theta \rangle_{\partial \mathcal{T}_h},$$

$$m_4 = - \langle \varepsilon_h^p - \varepsilon_h^p, (\pi_V \boldsymbol{\xi} - \boldsymbol{\xi}) \cdot \mathbf{n} \rangle_{\partial \mathcal{T}_h}.$$

We now estimate the terms $m_1, m_2, m_3$ and $m_4$ separately.

### 5.3.2. Estimating $m_1$

We follow the usual computation used in HDG convergence proof to combine the terms in $m_1$, see [DS19, Proof of Prop. 3.8]. We have

$$\begin{aligned} (\varepsilon_h^p, \rho_0(\pi_W \theta - \theta))_{\mathcal{T}_h} - (\delta_h^p, \rho_0 \pi_W \theta)_{\mathcal{T}_h} &= (\varepsilon_h^p - \delta_h^p, \rho_0(\pi_W \theta - \theta))_{\mathcal{T}_h} - (\delta_h^p, \rho_0 \theta)_{\mathcal{T}_h}, \quad (5.17) \\ &= (\varepsilon_h^p - \delta_h^p, \rho_0(\pi_W \theta - \theta))_{\mathcal{T}_h} - (\delta_h^p, \rho_0 \theta - \{\rho_0 \theta\})_{\mathcal{T}_h}, \end{aligned}$$

where $\{u\}|_K \in \mathcal{P}_0(K)$ is the average of $u$ on each element $K \in \mathcal{T}_h$, and as $\pi_W$ is the $L^2$-orthogonal projection onto $W_h(K)$, which implies that

$$(\delta_h^p, \{\rho_0 \theta\})_{\mathcal{T}_h} = (p - \pi_W p, \{\rho_0 \theta\})_{\mathcal{T}_h} = 0.$$

A similar computation can be performed for the terms involving $\varepsilon_h^q$ and $\delta_h^q$, yielding

$$\begin{aligned} |m_1| &\lesssim (\|\varepsilon_h^p\|_{\mathcal{T}_h} + \|\delta_h^p\|_{\mathcal{T}_h}) \|\pi_W \theta - \theta\|_{\mathcal{T}_h} + \|\delta_h^p\|_{\mathcal{T}_h} \|\rho_0 \theta - \{\rho_0 \theta\}\|_{\mathcal{T}_h} \\ &\quad + (\|\varepsilon_h^q\|_{\mathcal{T}_h} + \|\delta_h^q\|_{\mathcal{T}_h}) \|\pi_V \boldsymbol{\xi} - \boldsymbol{\xi}\|_{\mathcal{T}_h} + \|\delta_h^q\|_{\mathcal{T}_h} \|\mathbf{W}_0 \boldsymbol{\xi} - \{\mathbf{W}_0 \boldsymbol{\xi}\}\|_{\mathcal{T}_h}, \\ |m_1| &\lesssim (\|\varepsilon_h^p\|_{\mathcal{T}_h} + h\Delta_h) \|\pi_W \theta - \theta\|_{\mathcal{T}_h} + h\Delta_h \|\rho_0 \theta - \{\rho_0 \theta\}\|_{\mathcal{T}_h} \\ &\quad + (\mathcal{E}_h + \Delta_h) \|\pi_V \boldsymbol{\xi} - \boldsymbol{\xi}\|_{\mathcal{T}_h} + \Delta_h \|\mathbf{W}_0 \boldsymbol{\xi} - \{\mathbf{W}_0 \boldsymbol{\xi}\}\|_{\mathcal{T}_h}. \end{aligned} \quad (5.18)$$

### 5.3.3. Estimating $m_2$

An integration by parts yields

$$(\mathbf{b}_0 \cdot \nabla \varepsilon_h^p, \pi_W \theta - \theta)_{\mathcal{T}_h} - \langle (\mathbf{b}_0 \cdot \mathbf{n}) \varepsilon_h^p, \pi_W \theta - \theta \rangle_{\partial \mathcal{T}_h} = - (\varepsilon_h^p \mathbf{b}_0, \nabla[\pi_W \theta - \theta])_{\mathcal{T}_h},$$

and we therefore have

$$m_2 = -2i\omega (\varepsilon_h^p \mathbf{b}_0, \nabla[\pi_W \theta - \theta])_{\mathcal{T}_h} + 2i\omega (\delta_h^p \mathbf{b}_0, \nabla \pi_W \theta)_{\mathcal{T}_h}.$$

Following the same computation as the ones used to estimate $m_1$, we obtain

$$\begin{aligned} |m_2| &\lesssim (\|\varepsilon_h^p\|_{\mathcal{T}_h} + \|\delta_h^p\|_{\mathcal{T}_h}) \|\nabla[\pi_W \theta - \theta]\|_{\mathcal{T}_h} + \|\delta_h^p\|_{\mathcal{T}_h} \|\mathbf{b}_0 \cdot \nabla \theta - \{\mathbf{b}_0 \cdot \nabla \theta\}\|_{\mathcal{T}_h}, \\ |m_2| &\lesssim (\|\varepsilon_h^p\|_{\mathcal{T}_h} + h\Delta_h) \|\nabla[\pi_W \theta - \theta]\|_{\mathcal{T}_h} + h\Delta_h \|\mathbf{b}_0 \cdot \nabla \theta - \{\mathbf{b}_0 \cdot \nabla \theta\}\|_{\mathcal{T}_h}. \end{aligned} \quad (5.19)$$
### 5.3.4. Estimating $m_3$

We recall that

$$ \hat{Q}_1 := 2i\omega(\mathbf{b}_0 \cdot \mathbf{n})\hat{\varepsilon}_h^p + 2i\omega\tau_{\text{upw}}(\varepsilon_h^p - \hat{\varepsilon}_h^p), $$

$$ \hat{Q}_2 := 2i\omega(\mathbf{b}_0 \cdot \mathbf{n})\hat{\delta}_h^p - 2i\omega\tau_{\text{upw}}(\delta_h^p - \hat{\delta}_h^p), $$

$$ \hat{Q}_3 := \boldsymbol{\delta}_h^q \cdot \mathbf{n} + 2i\omega\tau P_M \delta_h^p, $$

and that we have established the following estimate in LEMMA 5.5,

$$ \|\varepsilon_h^p - \hat{\varepsilon}_h^p\|_{\partial \mathcal{T}_h} \lesssim h^{\frac{1}{2}} (\mathcal{E}_h + \Delta_h). \qquad (5.20) $$

We have

$$ \begin{aligned} \left| \langle (\mathbf{b}_0 \cdot \mathbf{n}) \hat{\varepsilon}_h^p, \pi_W \theta - P_M \theta \rangle_{\partial \mathcal{T}_h} \right| & \lesssim \left| \langle (\mathbf{b}_0 \cdot \mathbf{n})(\hat{\varepsilon}_h^p - \varepsilon_h^p), \pi_W \theta - P_M \theta \rangle_{\partial \mathcal{T}_h} \right| + \left| \langle (\mathbf{b}_0 \cdot \mathbf{n}) \varepsilon_h^p, \pi_W \theta - P_M \theta \rangle_{\partial \mathcal{T}_h} \right|, \\ & \lesssim \left[ h^{\frac{1}{2}} (\mathcal{E}_h + \Delta_h) + h^{-\frac{1}{2}} \|\varepsilon_h^p\|_{\mathcal{T}_h} \right] \|\pi_W \theta - P_M \theta\|_{\partial \mathcal{T}_h}, \end{aligned} $$

where we used the trace inequality (5.2). Using this trace inequality (5.2) again, that $\tau = \mathcal{O}(h^{-1})$ and keeping only the lowest power of $h$ as $h \to 0$, we obtain

$$ |m_3| \lesssim \|\pi_W \theta - P_M \theta\|_{\partial \mathcal{T}_h} \left[ h^{\frac{1}{2}} \|\boldsymbol{\varepsilon}_h^q\|_{\mathbf{W}_0, \mathcal{T}_h} + \left\| |\tau|^{\frac{1}{2}} (P_M \boldsymbol{\varepsilon}_h^p - \hat{\boldsymbol{\varepsilon}}_h^p) \right\|_{\partial \mathcal{T}_h} + h^{-\frac{1}{2}} \left( \|\boldsymbol{\varepsilon}_h^p\|_{\mathcal{T}_h} + \|\boldsymbol{\delta}_h^q\|_{\mathcal{T}_h} \right) + h^{-\frac{3}{2}} \|\boldsymbol{\delta}_h^p\|_{\mathcal{T}_h} + \left\| \hat{\boldsymbol{\delta}}_h^p \right\|_{\partial \mathcal{T}_h} \right], $$

which yields

$$ |m_3| \lesssim [\mathcal{E}_h + h^{-\frac{1}{2}} \Delta_h + h^{-\frac{1}{2}} \|\boldsymbol{\varepsilon}_h^p\|_{\mathcal{T}_h}] \|\pi_W \theta - P_M \theta\|_{\partial \mathcal{T}_h}. \qquad (5.21) $$

### 5.3.5. Estimating $m_4$

Using (5.20) and the trace inequality (5.2), we readily obtain

$$ |m_4| \lesssim (\mathcal{E}_h + \Delta_h) \|\boldsymbol{\pi}_V \boldsymbol{\xi} - \boldsymbol{\xi}\|_{\mathcal{T}_h}. \qquad (5.22) $$

### 5.3.6. Proof of the Aubin-Nitsche estimate

By collecting the estimates (5.18), (5.19), (5.21) and (5.22) and by using them into (5.17), we obtain

$$ \begin{aligned} \|\boldsymbol{\varepsilon}_h^p\|_{\mathcal{T}_h}^2 & \lesssim (\|\boldsymbol{\varepsilon}_h^p\|_{\mathcal{T}_h} + h\Delta_h) \|\pi_W \theta - \theta\|_{\mathcal{T}_h} + h\Delta_h \|\rho_0 \theta - \{\rho_0 \theta\}\|_{\mathcal{T}_h} \\ & + (\mathcal{E}_h + \Delta_h) \|\boldsymbol{\pi}_V \boldsymbol{\xi} - \boldsymbol{\xi}\|_{\mathcal{T}_h} + \Delta_h \|\mathbf{W}_0 \boldsymbol{\xi} - \{\mathbf{W}_0 \boldsymbol{\xi}\}\|_{\mathcal{T}_h} \\ & + (\|\boldsymbol{\varepsilon}_h^p\|_{\mathcal{T}_h} + h\Delta_h) \|\nabla[\pi_W \theta - \theta]\|_{\mathcal{T}_h} + h\Delta_h \|\mathbf{b}_0 \cdot \nabla \theta - \{\mathbf{b}_0 \cdot \nabla \theta\}\|_{\mathcal{T}_h} \\ & + [\mathcal{E}_h + h^{-\frac{1}{2}}\Delta_h + h^{-\frac{1}{2}}\|\boldsymbol{\varepsilon}_h^p\|_{\mathcal{T}_h}] \|\pi_W \theta - P_M \theta\|_{\partial \mathcal{T}_h} \\ & + (\mathcal{E}_h + \Delta_h) \|\boldsymbol{\pi}_V \boldsymbol{\xi} - \boldsymbol{\xi}\|_{\mathcal{T}_h} \end{aligned} $$
Using the elliptic regularity assumption (5.4) and the approximation property of the $L^2$-orthogonal projections (5.1) for $\xi \in H^1(\mathcal{O})$ and $\theta \in H^2(\mathcal{O})$, we have

$$
\begin{align*}
\|\pi_W \theta - \theta\|_{\rho_0, \mathcal{T}_h} &\lesssim h^2 \|\varepsilon_h^p\|_{\mathcal{T}_h} \lesssim h \|\varepsilon_h^p\|_{\mathcal{T}_h}, \\
\|\pi_W \xi - \xi\|_{\mathcal{T}_h} &\lesssim h \|\varepsilon_h^p\|_{\mathcal{T}_h}, \\
\|\rho_0 \theta - \{\rho_0 \theta\}\|_{\mathcal{T}_h} &\lesssim h \|\varepsilon_h^p\|_{\mathcal{T}_h}, \\
\|\mathbf{W}_0 \xi - \{\mathbf{W}_0 \xi\}\|_{\mathcal{T}_h} &\lesssim h \|\varepsilon_h^p\|_{\mathcal{T}_h}, \\
\|\nabla[\pi_W \theta - \theta]\|_{\mathcal{T}_h} &\lesssim h \|\theta\|_{2, \mathcal{O}} \lesssim h \|\varepsilon_h^p\|_{\mathcal{T}_h},
\end{align*}
$$

for the approximation properties of the projections in the $H^1$ semi-norm, we refer to [PE12, Lemma 1.58]. We also notice that

$$
\begin{align*}
\|\pi_W \theta - P_M \theta\|_{\partial \mathcal{T}_h} &\le \|\pi_W \theta - \theta\|_{\partial \mathcal{T}_h} + \|\theta - P_M \theta\|_{\partial \mathcal{T}_h}, \\
&\lesssim h^{-\frac{1}{2}} \|\pi_W \theta - \theta\|_{\partial \mathcal{T}_h} + \|\theta - P_M \theta\|_{\partial \mathcal{T}_h}, \\
&\lesssim h^{\frac{3}{2}} \|\varepsilon_h^p\|_{\mathcal{T}_h},
\end{align*}
$$

where we used the trace inequality (5.2) and the elliptic regularity assumption (5.4). We therefore obtain

$$
\|\varepsilon_h^p\|_{\rho_0, \mathcal{T}_h}^2 \lesssim h \|\varepsilon_h^p\|_{\mathcal{T}_h}^2 + h \|\varepsilon_h^p\|_{\mathcal{T}_h} (\mathcal{E}_h + \Delta_h).
$$

Finally the absorption technique of LEMMA 5.6 gives

$$
\| \varepsilon_h^p \|_{\mathcal{T}_h} \lesssim h (\varepsilon_h + \Delta_h),
$$

which is exactly the Aubin-Nitsche estimate (5.6b).

5.4. Additional results

We complete this section with two additional results. The first one is a convergence estimate for the numerical trace $\widehat{p_h}$ of the HDG+ method.

**Corollary 5.1.** Under the assumptions of **THEOREM 2**, the following error estimates for $\widehat{p_h}$ hold

$$
\|\widehat{\varepsilon}_h^p\|_{\partial \mathcal{T}_h} = \mathcal{O}(h^{k+\frac{3}{2}}) \quad \text{and} \quad \|p - \widehat{p}_h\|_{\partial \mathcal{T}_h} = \mathcal{O}(h^{k+\frac{1}{2}}).
$$

**Proof:** We have

$$
\|\widehat{\varepsilon}_h^p\|_{\partial \mathcal{T}_h} \le \|P_M \varepsilon_h^p\|_{\partial \mathcal{T}_h} + \|P_M \varepsilon_h^p - \widehat{\varepsilon}_h^p\|_{\partial \mathcal{T}_h}
$$

using that $\tau = \mathcal{O}(h^{-1})$, we get

$$
\lesssim \|P_M \varepsilon_h^p\|_{\partial \mathcal{T}_h} + h^{\frac{1}{2}} \left\| \tau^{\frac{1}{2}} (P_M \varepsilon_h^p - \widehat{\varepsilon}_h^p) \right\|_{\partial \mathcal{T}_h},
$$

the continuity of $P_M$ yields

$$
\lesssim \|\varepsilon_h^p\|_{\partial \mathcal{T}_h} + h^{\frac{1}{2}} \left\| \tau^{\frac{1}{2}} (P_M \varepsilon_h^p - \widehat{\varepsilon}_h^p) \right\|_{\partial \mathcal{T}_h},
$$
we then resort to the trace estimate (5.2) to obtain

$$ \lesssim h^{-\frac{1}{2}} \|\varepsilon_h^p\|_{\mathcal{T}_h} + h^{\frac{1}{2}} \left\|\tau|^{\frac{1}{2}}(P_M \varepsilon_h^p - \hat{\varepsilon}_h^p)\right\|_{\partial \mathcal{T}_h}, $$

finally, we use the error estimates of LEMMA 5.1 and we end up with

$$ \|\hat{\varepsilon}_h^p\|_{\partial \mathcal{T}_h} \lesssim \mathcal{O}(h^{-\frac{1}{2}}h^{k+2} + h^{\frac{1}{2}}h^{k+1}) = \mathcal{O}(h^{k+\frac{3}{2}}), $$

which is the first estimate. The second one comes from (5.1a) with $s = k + 1$. ■
The second result is the sub-optimal convergence of the standard HDG method.

**Proposition 5.1** (Convergence of the standard HDG method). Under the assumptions of THEOREM 2 and if $\tau = \mathcal{O}(1)$, the standard HDG method with $V_h = \mathcal{P}_k(\mathcal{T}_h)$, $W_h = \mathcal{P}_k(\mathcal{T}_h)$ and $M_h = \mathcal{P}_k(\mathcal{F}_h)$ satisfies

$$ \|\varepsilon_h^p\|_{\mathcal{T}_h} = \mathcal{O}(h^{k+\frac{3}{2}}), \quad \text{and} \quad \|\varepsilon_h^q\|_{\mathcal{T}_h} = \mathcal{O}(h^{k+\frac{1}{2}}). $$

Those error estimates can be obtained by a straightforward adaptation of the proof of THEOREM 2.

## 6. GLOBAL SOLVABILITY

The analysis that we have carried out in the previous subsection works for any solution $(q_h, p_h, \hat{p}_h)$ of the discrete system (3.7a)–(3.7b)–(3.7c) provided that such solution exists. We already discussed the well-posedness of the local problems in THEOREM 1, but we have not yet proved that the global problem (3.8) for $\hat{p}_h$ was well-posed. To do that we can either directly show the well-posedness of the global problem (3.8). Or we can choose take advantage of the error estimates of THEOREM 2 as we will describe below<fnref n="5" />.

We recall that the convected Helmholtz equation is of Fredholm type. It is therefore uniquely solvable except on a set of resonant frequencies. For those frequencies, there exist non-zero solutions to the homogeneous equation and unique solvability cannot be guaranteed.

We can now state and prove the main result of this section.

**Theorem 3** (Global solvability). Under the assumptions of THEOREM 1 and THEOREM 2 and if $\omega$ is not a resonant frequency of the convected Helmholtz equation (1.1) then the global problem is well-posed, ie $\hat{p}_h$ is uniquely defined by (3.8).

**Proof:** First we recall that (3.7a)–(3.7b)–(3.7c), or equivalently (3.8), is a square system of linear equations, we therefore only need to show the uniqueness of the solution of the homogeneous system (when $g_N = s = 0$). Assuming that $\omega$ is not a resonant frequency of (1.1), the exact solution is $p = 0$ and $\sigma = 0$, and therefore

$$ \|p\|_{s,\mathcal{O}} = 0 \quad \text{and} \quad \|q\|_{t,\mathcal{O}} = 0 $$

and

$$ \varepsilon_h^p = -p_h ; \quad \varepsilon_h^q = -q_h ; \quad \hat{\varepsilon}_h^p = -\hat{p}_h. $$

During the error analysis, we have proven that

$$ \|\varepsilon_h^p\|_{\mathcal{T}_h} \lesssim \|p\|_{s,\mathcal{O}} + \|q\|_{t,\mathcal{O}} = 0, \quad (6.1a) $$

$$ \|\varepsilon_h^q\|_{\mathcal{T}_h} \lesssim \|p\|_{s,\mathcal{O}} + \|q\|_{t,\mathcal{O}} = 0, \quad (6.1b) $$

$$ \left\| |\tau|^{\frac{1}{2}} (P_M \varepsilon_h^p - \hat{\varepsilon}_h^p) \right\|_{\partial \mathcal{T}_h} \lesssim \|p\|_{s,\mathcal{O}} + \|q\|_{t,\mathcal{O}} = 0, \quad (6.1c) $$

<footnote n="5">⁵In [DS19] this idea is attributed to B. Cockburn.</footnote>
when h is sufficiently small.

Notice that we have hidden the powers of h in ≲ as they do not play an important part here.

Therefore using (6.1a) and (6.1b) we have shown that

$$p_h \equiv 0 \quad \text{and} \quad \sigma_h \equiv 0$$

when h is small enough. We also notice that

$$\begin{aligned}
\| |\tau|^{\frac{1}{2}} (P_M \varepsilon_h^p - \widehat{\varepsilon}_h^p) \|_{\partial \mathcal{T}_h} &= \| |\tau|^{\frac{1}{2}} (\widehat{p}_h - P_M p_h) \|_{\partial \mathcal{T}_h} \\
&= \| |\tau|^{\frac{1}{2}} \widehat{p}_h \|_{\partial \mathcal{T}_h}, \\
&= 0,
\end{aligned}$$

by using (6.1c) when h is small enough as $\tau \neq 0$. We therefore obtain that

$$\widehat{p}_h \equiv 0.$$

## 7. NUMERICAL EXPERIMENTS

In this section, we present some numerical experiments to illustrate our theoretical results. This method has been implemented in the open-source software hawen, see [Fau21]. As most of the estimates obtained in our analysis involve projection errors of the form

$$\|p_h - \pi_W p\|_{\mathcal{T}_h} \quad \text{or} \quad \|q_h - \pi_V q\|_{\mathcal{T}_h},$$

we will therefore provide numerical errors involving those projections instead of the usual terms

$$\|p_h - p\|_{\mathcal{T}_h} \quad \text{or} \quad \|q_h - q\|_{\mathcal{T}_h}.$$

As the projections used are $L^2$-orthogonal projections, they can be computed by solving a linear system. This is done using lapack and a 91-point Gauss-Lobatto quadrature rule to evaluate the integrals. For the purpose of comparing numerical errors on different meshes, we introduce the relative $L^2$-errors

$$\mathcal{E}_u := \frac{\|u_h - \pi u\|_{\mathcal{T}_h}}{\|\pi u\|_{\mathcal{T}_h}}, \quad \text{for } u \in \{p, q_x, q_y\}.$$

Those quantities are plotted against $k/h$ which is proportional to the number of degrees of freedom per wavelength and a log-log scale is used.

### 7.1. Geometrical settings

As depicted on FIGURE 7.1 we consider a uniform directional flow $v_0 = Mc_0e_x$, where M is the Mach number.
<figure><img src="image_8.png" /><figcaption>FIGURE 7.1. Sketch of the geometrical configuration</figcaption></figure>

Unless stated otherwise, we will always use the following parameters for the convergence tests

$$ \mathcal{O} = (0, 2) \times (0, 1) ; \rho_0, c_0 \equiv 1 ; \omega = 5.55\pi, $$

and the choice of $M$ will be specified for each numerical experiment.

The duct modes are a family of analytic solutions of (1.1) in a waveguide, see [BBL03]. They are given by

$$ p_n^{\pm}(x, y) = e^{i\beta_n^{\pm}x} \varphi_n(y) $$

where

$$ \begin{aligned} n < N_0 : \qquad & \beta_n^{\pm} = \frac{-\kappa M \pm \sqrt{\kappa^2 - \frac{n^2\pi^2}{\ell^2}(1-M^2)}}{1-M^2} \\ n > N_0 : \qquad & \beta_n^{\pm} = \frac{-\kappa M \pm i\sqrt{\frac{n^2\pi^2}{\ell^2}(1-M^2) - \kappa^2}}{1-M^2} \end{aligned} $$

with

$$ \begin{gathered} \kappa = \frac{\omega}{c_0} \quad \text{and} \quad M = \frac{v_0}{c_0} \\ N_0 = \left\lfloor \frac{\kappa \ell}{\pi \sqrt{1-M^2}} \right\rfloor \end{gathered} $$

and

$$ \begin{aligned} \varphi_0(y) &:= \sqrt{\ell^{-1}} \\ \varphi_n(y) &:= \sqrt{2\ell^{-1}} \cos\left(\frac{n\pi y}{\ell}\right), \quad n \in \mathbb{N}^* \end{aligned} $$

The choice of $n$ will be specified for each numerical experiment.

## 7.2. Numerical experiments with a low Mach number

We then move to a flow with a low Mach number. In this case we have used the following parameters

$n = 3, \quad \text{and} \quad M = 0.2.$

On FIGURE 7.2 the convergence rate for the volumetric unknown $p_h$ for the HDG+ method is displayed. As expected the optimal convergence rate of $k+2$ is obtained. If we now move to the flux unknown, we can see on FIGURE 7.3 that $q_h$ converges with the optimal order $k+1$ as expected.

In SUBSECTION 3.7 and in PROPOSITION 5.1, we have mentioned that using a standard HDG method in the diffusive-flux formulation leads to a sub-optimal convergence rate for the vector unknown. This is illustrated in FIGURE 7.4 where a $k+1/2$ rate can be seen instead of the optimal $k+1$ rate. By reproducing the convergence
<figure><img src="image_9.png" /><figcaption>FIGURE 7.2. Low Mach convergence history for the volumetric unknown $p_h$ for the HDG+ method with $(q_h, p_h, \hat{p}_h) \in \mathcal{P}_k \times \mathcal{P}_{k+1} \times \mathcal{P}_k$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 7.3. Low Mach convergence history for the first component of volumetric unknown $q_h$ for the HDG+ method with $(q_h, p_h, \hat{p}_h) \in \mathcal{P}_k \times \mathcal{P}_{k+1} \times \mathcal{P}_k$</figcaption></figure>

proof of the HDG+ method, we obtain a convergence rate of $k + 3/2$ for $p_h$ instead of the super-convergent rate of $k + 2$ expected for HDG methods. Similarly, we get a $k + 1/2$ convergence rate for $q_h$ instead of the optimal rate of $k + 1$ that is expected for HDG methods. This can be seen in FIGURE 7.5. Even if this rate is higher than the optimal convergence rate, it cannot be exploited through a post-processing scheme to get a more accurate solution.

The convergence history for the volumetric unknown for the HDG+ method is displayed on FIGURE 7.6 and we can see that it still achieves the optimal convergence rate of $k + 2$. The convergence history for the volumetric
<figure><img src="image_11.png" /><figcaption>FIGURE 7.4. Convergence of the vectorial unknown $q_h$ with the standard HDG method with $(q_h, p_h, \hat{p}_h) \in \mathcal{P}_k \times \mathcal{P}_k \times \mathcal{P}_k$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 7.5. Convergence of the scalar unknown $p_h$ with the standard HDG method with $(q_h, p_h, \hat{p}_h) \in \mathcal{P}_k \times \mathcal{P}_k \times \mathcal{P}_k$</figcaption></figure>

flux unknown $q_h$ is depicted in FIGURE 7.7 for the HDG+ method. As in the low-Mach case, the HDG+ method has a convergence rate of $k + 1$.

CONCLUSION

In this paper we have introduced and analyzed a HDG+ method for the diffusive-flux formulation of the convected Helmholtz equation. The scalar unknown is approximated with polynomials of degree $k + 1$ whereas the volumetric vector unknown and the numerical trace are both approximated with polynomials of degree
<figure><img src="image_13.png" /><figcaption>FIGURE 7.6. Large Mach convergence history for the volumetric unknown $p_h$ for the HDG+ method with $(q_h, p_h, \hat{p}_h) \in \mathcal{P}_k \times \mathcal{P}_{k+1} \times \mathcal{P}_k$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 7.7. Large Mach convergence history for the second component of the volumetric unknown $q_h$ for the HDG+ method</figcaption></figure>

$k$. The introduction of this HDG+ method is motivated by the observation that the standard HDG method (where all the unknowns are approximated by polynomials of degree $k$) exhibits a sub-optimal convergence rate and can be unstable, see FIGURE 3.1 and FIGURE 7.4. As it was established in [BRT23], it is possible to devise a super-convergent HDG method when using a different first-order in space formulation of the convected Helmholtz equation.

We provide detailed theoretical results on well-posedness and convergence. In particular the HDG+ is super-convergent: $p_h$ converges to $p$ at order $k+2$, but the global problem is only of degree $k$. We therefore obtained
a convergence rate of $k + 2$ for the cost of a method of order $k + 1$. The choice of approximation spaces avoids applying post-processing schemes. Numerical experiments illustrating this property, as well as the sub-optimal convergence of the standard HDG method, are also provided.

## Acknowledgement

An open-source implementation of those three methods in the open-source *hawen* solver, see [Fau21], will be released soon. The authors would like to thank its main developer, Florian Faucher, for his help with the numerical implementation.

Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr).

Nathan Rouxelin acknowledges financial support from e2s-UPPA (see https://e2s-uppa.eu), Maison Normande des Sciences du Numérique (see http://www.criann.fr/mnsn/), and Région Normandie.

## REFERENCES

[BBL03] Eliane Bécache, Anne-Sophie Bonnet-Ben Dhia, and Guillaume Legendre. Perfectly matched layers for the convected Helmholtz equation. Technical report, 2003. Cited on page 32.

[BC13] Lucio Boccardo and Gisella Croce. *Elliptic Partial Differential Equations*. De Gruyter, Berlin, Boston, 2013. Cited on page 18.

[BCD+15] Nolwenn Balin, Fabien Casenave, François Dubois, Eric Duceau, Stefan Duprey, and Isabelle Terrasse. Boundary element and finite element coupling for aeroacoustics simulations. *Journal of Computational Physics*, 294:274–296, August 2015. Cited on page 1.

[BCDL15] Marie Bonnasse-Gahot, Henri Calandra, Julien Diaz, and Stéphane Lanteri. Hybridizable Discontinuous Galerkin method for the simulation of the propagation of the elastic wave equations in the frequency domain. Research Report RR-8990, INRIA Bordeaux ; INRIA Sophia Antipolis - Méditerranée, June 2015. Cited on page 2.

[BDMP21] Hélène Barucq, Julien Diaz, Rose-Cloé Meyer, and Ha 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. Cited on page 2.

[BRT23] Hélène Barucq, Nathan Rouxelin, and Sébastien 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, May 2023. Cited on pages 2 and 35.

[BS97] Ivo M Babuska and Stefan A Sauter. Is the pollution effect of the fem avoidable for the helmholtz equation considering high wave numbers? *SIAM Journal on numerical analysis*, 34(6):2392–2423, 1997. Cited on page 1.

[CC12] Yanlai Chen and Bernardo Cockburn. Analysis of variable-degree HDG methods for convection–diffusion equations. Part I: General nonconforming meshes. *IMA Journal of Numerical Analysis*, 32(4):1267–1293, October 2012. Cited on page 2.

[CC14] Yanlai Chen and Bernardo Cockburn. Analysis of variable-degree HDG methods for convection-diffusion equations. Part II: Semimatching nonconforming meshes. *Mathematics of Computation*, 83(285):87–111, January 2014. Cited on page 2.

[CDG+09] Bernardo Cockburn, Bo Dong, Johnny Guzmán, Marco Restelli, and Riccardo Sacco. A Hybridizable Discontinuous Galerkin Method for Steady-State Convection-Diffusion-Reaction Problems. *SIAM Journal on Scientific Computing*, 31(5):3827–3846, January 2009. Cited on page 2.

[CDPE16] Bernardo Cockburn, Daniele A. Di Pietro, and Alexandre Ern. Bridging the hybrid high-order and hybridizable discontinuous Galerkin methods. *ESAIM: Mathematical Modelling and Numerical Analysis*, 50(3):635–650, May 2016. Cited on page 2.

[CGL09] Bernardo Cockburn, Jayadeep Gopalakrishnan, and Raytcho Lazarov. Unified Hybridization of Discontinuous Galerkin, Mixed, and Continuous Galerkin Methods for Second Order Elliptic Problems. *SIAM J. Numer. Anal*, 47:1319–1365, August 2009. Cited on pages 2 and 12.

[CGS10] Bernardo Cockburn, Jayadeep Gopalakrishnan, and Francisco-Javier Sayas. A projection-based error analysis of HDG methods. *Mathematics of Computation*, 79(271):1351–1367, March 2010. Cited on page 2.

[Chr04] Jorgen Christensen-Dalsgaard. *Lecture Notes on Stellar Oscillations*. 2004. Cited on page 1.

[CLOS20] Liliana Camargo, Bibiana López-Rodríguez, Mauricio Osorio, and Manuel Solano. An HDG method for Maxwell's equations in heterogeneous media. *Computer Methods in Applied Mechanics and Engineering*, 368:113178, August 2020. Cited on page 2.
[Coc14] Bernardo Cockburn. Static Condensation, Hybridization, and the Devising of the HDG Methods. springerprofes-
sional.de, 2014. Cited on page 2.

[CQS18] Huangxin Chen, Weifeng Qiu, and Ke Shi. A priori and computable a posteriori error estimates for an HDG method for the coercive Maxwell equations. *Computer Methods in Applied Mechanics and Engineering*, 333:287–310, May 2018. Cited on page 2.

[CQSS17] Huangxin Chen, Weifeng Qiu, Ke Shi, and Manuel Solano. A Superconvergent HDG Method for the Maxwell Equations. *Journal of Scientific Computing*, 70(3):1010–1029, March 2017. Cited on page 2.

[CS13] Bernardo Cockburn and Ke Shi. Superconvergent HDG methods for linear elasticity with weakly symmetric stresses | IMA Journal of Numerical Analysis | Oxford Academic. *IMA Journal of Numerical Analysis*, 2013. Cited on page 2.

[DS19] Shukai Du and Francisco-Javier Sayas. *An Invitation to the Theory of the Hybridizable Discontinuous Galerkin Method: Projections, Estimates, Tools*. SpringerBriefs in Mathematics. Springer International Publishing, Cham, 2019. Cited on pages 2, 16, 17, 27, and 30.

[EG04] Alexandre Ern and Jean-Luc Guermond. *Theory and Practice of Finite Elements*. Applied Mathematical Sciences. Springer-Verlag, New York, 2004. Cited on pages 5, 15, and 17.

[Fau21] Florian Faucher. ‘hawen’: Time-harmonic wave modeling and inversion using hybridizable discontinuous Galerkin discretization. *Journal of Open Source Software*, 6(57):2699, January 2021. Cited on pages 1, 31, and 36.

[FCS15] Guosheng Fu, Bernardo Cockburn, and Henryk Stolarski. Analysis of an HDG method for linear elasticity. *International Journal for Numerical Methods in Engineering*, 102(3-4):551–575, 2015. Cited on page 2.

[FS20] Florian Faucher and Otmar Scherzer. Adjoint-state method for Hybridizable Discontinuous Galerkin discretization, application to the inverse acoustic wave problem. *Computer Methods in Applied Mechanics and Engineering*, 372:113406, December 2020. Cited on page 1.

[GBD+17] Laurent Gizon, Hélène Barucq, Marc Duruflé, Chris S. Hanson, Michael Leguèbe, Aaron C. Birch, Juliette Chabassier, Damien Fournier, Thorsten Hohage, and Emanuele Papini. Computational helioseismology in the frequency domain: acoustic waves in axisymmetric solar models with flows. *Astronomy & Astrophysics*, 600:A35, March 2017. Cited on page 1.

[GM11] Roland Griesmaier and Peter Monk. Error Analysis for a Hybridizable Discontinuous Galerkin Method for the Helmholtz Equation. *Journal of Scientific Computing*, 49(3):291–310, December 2011. Cited on page 2.

[Gri11] Pierre Grisvard. *Elliptic Problems in Nonsmooth Domains*. Society for Industrial and Applied Mathematics, 2011. Cited on page 18.

[GSV18] Jay Gopalakrishnan, Manuel Solano, and Felipe Vargas. Dispersion Analysis of HDG Methods. *Journal of Scientific Computing*, 77(3):1703–1735, December 2018. Cited on page 2.

[HPS17] Allan Hungria, Daniele Prada, and Francisco-Javier Sayas. HDG methods for elastodynamics. *Computers & Mathematics with Applications*, 74(11):2671–2690, December 2017. Cited on page 2.

[Hun19] Allan Hungria. *Using HDG+ to Compute Solutions of the 3D Linear Elastic and Poroelastic Wave Equations*. PhD thesis, University of Delaware, 2019. Cited on pages 2, 16, and 17.

[IGS22] Lise-Marie Imbert-Gérard and Guillaume Sylvand. Three types of quasi-trefftz functions for the 3d convected Helmholtz equation: construction and theoretical approximation properties. 2022. Cited on page 2.

[KSC12] Robert J. Kirby, Spencer J. Sherwin, and Bernardo Cockburn. To CG or to HDG: A Comparative Study. *Journal of Scientific Computing*, 51(1):183–212, April 2012. Cited on pages 2, 5, and 6.

[Leh10] Christoph Lehrenfeld. *Hybrid Discontinuous Galerkin Methods for Solving Incompressible Flow Problems*. PhD thesis, RWTH Aachen, 2010. Cited on pages 2 and 9.

[LMG+20] Alice Lieu, Philippe Marchner, Gwénaël Gabard, Hadrien Bériot, Xavier Antoine, and Christophe Geuzaine. A non-overlapping schwarz domain decomposition method with high-order finite elements for flow acoustics. *Computer Methods in Applied Mechanics and Engineering*, 369:113223, September 2020. Cited on page 1.

[NPRC15] Ngoc-Cuong Nguyen, Jaime Peraire, Fernando Reitich, and Bernardo Cockburn. A phase-based hybridizable discontinuous Galerkin method for the numerical solution of the Helmholtz equation. *Journal of Computational Physics*, 290:318–335, June 2015. Cited on page 2.

[Oik14] Issei Oikawa. A hybridized discontinuous Galerkin method with reduced stabilization. *arXiv:1405.2491 [math]*, November 2014. Cited on page 2.

[Oik16] Issei Oikawa. Analysis of a Reduced-Order HDG Method for the Stokes Equations. *Journal of Scientific Computing*, 67(2):475–492, May 2016. Cited on pages 2 and 9.

[Oik18] Issei Oikawa. An HDG Method with Orthogonal Projections in Facet Integrals. *Journal of Scientific Computing*, 76(2):1044–1054, August 2018. Cited on page 2.

[PE12] Daniele Antonio Di Pietro and Alexandre Ern. *Mathematical Aspects of Discontinuous Galerkin Methods*. Mathématiques et Applications. Springer-Verlag, Berlin Heidelberg, 2012. Cited on page 29.

[Pie90] Allan Pierce. Wave equation for sound in fluids with unsteady inhomogeneous flow. *The Journal of the Acoustical Society of America*, 87(6):2292–2299, June 1990. Cited on page 3.

[QS16a] Weifeng Qiu and Ke Shi. An HDG Method for Convection Diffusion Equation. *Journal of Scientific Computing*, 66(1):346–357, January 2016. Cited on pages 2, 10, 17, and 20.
[QS16b] Weifeng Qiu and Ke Shi. A superconvergent HDG method for the incompressible Navier-Stokes equations on general polyhedral meshes. *IMA Journal of Numerical Analysis*, 36(4):1943–1967, October 2016. Cited on page 2.

[QSS16] Weifeng Qiu, Jiguang Shen, and Ke Shi. An HDG method for linear elasticity with strong symmetric stresses. *arXiv:1312.1407 [math]*, February 2016. Cited on pages 2 and 17.

[Rou21] Nathan Rouxelin. *Condensed Mixed Numerical Methods for Convected Acoustics. Applications in Helioseismology*. PhD thesis, Université de Pau et des Pays de l'Adour, 2021. Cited on page 12.

[SM50] Jack Sherman and Winifried Morrison. Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix. *Annals of Mathematical Statistics*, 21(1):124–127, March 1950. Cited on page 4.

[YMKS16] Sergey Yakovlev, David Moxey, Robert M. Kirby, and Spencer J. Sherwin. To CG or to HDG: A Comparative Study in 3D. *Journal of Scientific Computing*, 67(1):192–220, April 2016. Cited on page 2.