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

HAL
open science

# A multiscale hybrid-mixed method for the Helmholtz equation in heterogeneous domains

Théophile Chaumont-Frelet, Frédéric Valentin

► To cite this version:

Théophile Chaumont-Frelet, Frédéric Valentin. A multiscale hybrid-mixed method for the Helmholtz equation in heterogeneous domains. 2019. hal-01698914v2

HAL Id: hal-01698914

https://inria.hal.science/hal-01698914v2

Preprint submitted on 16 Apr 2019

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.
# A MULTISCALE HYBRID-MIXED METHOD FOR THE HELMHOLTZ EQUATION IN HETEROGENEOUS DOMAINS

THÉOPHILE CHAUMONT-FRELET * AND FRÉDÉRIC VALENTIN †

**Abstract.** This work proposes a novel multiscale finite element method for acoustic wave propagation in highly heterogeneous media which is accurate on coarse meshes. It originates from the primal hybridization of the Helmholtz equation at the continuous level, which relaxes the continuity of the unknown on the skeleton of a partition. As a result, face-based degrees of freedom drive the approximation on the faces, and independent local problems respond for the multiscale basis function computation. We show how to recover other well-established numerical methods when the basis functions are promptly available. A numerical analysis of the method establishes its well-posedness and proves the quasi-optimality for the numerical solution. Also, we demonstrate that the MHM method is super-convergent in natural norms. We assess theoretical results, as well as the performance of the method on heterogeneous domains, through a sequence of numerical tests.

**Key words.** Helmholtz equation, hybrid methods, multiscale methods, finite element methods

**AMS subject classifications.** 65N12, 65N15, 65N30

**1. Introduction.** Numerical approximation of waves deserves particular attention by the scientific community. Indeed, though wave propagation problems arise in a wide range of applications, solving realistic 3D problems in the high-frequency regime is still either very costly or impossible. Frequency-domain discretization methods for wave problems face the so-called “pollution effect”: the number of discretization points per wavelength must increase when the frequency is high. The pollution effect leads to restrictive conditions on the mesh size in the high-frequency regime, making the discretization of high-frequency problems computationally intensive.

The pollution effect has been vastly studied and is well-understood when the medium of propagation is homogeneous [1, 4, 26, 27, 31]. In order to reduce the pollution effect, the idea to capture *a priori* knowledge of the solution in the basis functions has recently emerged. The critical ingredient is to use local solutions to the Helmholtz equation, usually plane waves, as basis functions. Because plane waves are less flexible than polynomials, special techniques are required to construct a conforming discretization space or to stabilize the method. These techniques include the partition of unity [40], least squares methods [33], ultra weak variational formulations [9], the variational theory of complex rays [38] and the discontinuous enrichment method (DEM) [15].

The aforementioned methods assume the medium of propagation is homogeneous. When the propagation medium is highly heterogeneous, the wave speed might change inside the mesh elements, so that local analytical solutions are no longer available. Extensions of plane wave methods to heterogeneous media have been the subject of intense investigation. For instance, when the propagation medium varies smoothly, the so-called generalized plane wave methods can be used [28, 41]. Generalized plane waves assume that a polynomial function can represent the wave speed, so that a close formula is available for the basis functions. On the other hand, when the wave speed is piecewise constant, plane wave methods can still apply, but it requires that the wave speed be constant in each element.

Unfortunately, some applications, like seismic wave propagation, feature highly heterogeneous media with jumps in the wave speed. In this context, it is not clear that generalized plane waves can be adopted. Besides, restricting the mesh so that the wave speed is constant in each element is not always possible. T. Chaumont-Frelet and collaborators have developed a multiscale strategy based on high order polynomial basis functions [5, 16, 17]: the Multiscale Medium Approximation Method (MMAM). When applied to geophysical benchmarks, it yields excellent results for the constant density medium case. Unfortunately, the cases of acoustic medium with non-constant density and

<footnote>* theophile.chaumont@inria.fr. Inria Sophia Antipolis-Méditerranée Research Centre, Nachos Project Team and University of Nice-Sophia Antipolis, 2004 route des Lucioles BP 93, 06902 Sophia Antipolis Cedex, France - J. A. Dieudonné Mathematics Laboratory (UMR CNRS 6621), Parc Valrose 06108 Nice, Cedex 02, France</footnote>

<footnote>† valentin@lncc.br and frederic.valentin@inria.fr. Department of Numerical and Computational Methods, LNCC - National Laboratory for Scientific Computing, Av. Getúlio Vargas, 333, 25651-070 Petrópolis - RJ, Brazil, and Université Côte d'Azur, Inria, CNRS, LJAD, 06902 Sophia Antipolis Cedex, France</footnote>
elastic medium are more tricky to discretize. Though the MMAM outperforms standard finite element discretization in these cases, the results are not fully satisfactory, especially in the case of elastic problems in the high-frequency regime.

In this paper, we propose a multiscale method for the acoustic Helmholtz equation in highly heterogeneous media featuring jumps in the wave speed and density. Our multiscale strategy is an adaptation of the Multiscale Hybrid Mixed (MHM) method initially developed for Darcy flows in highly heterogeneous media [3,22,35], and further extended to other operators in [21, 23, 24, 29]. The key idea is to solve element-level heterogeneous Helmholtz problems to provide basis functions, leading to a multiscale strategy.

Like the DEM [15], the MHM method stems from the primal hybrid formulation of the Helmholtz equation. As a result, we retrieve the lowest-order DEM elements as a particular case of the MHM method when the medium of propagation is homogeneous. In fact, it turns out that the MHM method has a lot in common with the DEM when the medium is homogeneous. To the best of our knowledge, convergence analysis of the DEM is currently limited to the lowest-order elements [2] and, then, the material presented hereafter might be helpful for the analysis of the DEM as well.

We note that recently, other multiscale strategies for highly heterogeneous Helmholtz problems appeared [8,19,34]. In [8,19], the authors introduce a Petrov-Galerkin stabilization technique that ensures the quasi-optimality of the standard piecewise linear finite element space, using multiscale test functions. The key asset of this method is that it is stable under the sole assumption that there are sufficiently many degrees of freedom per wavelength. This stability arises at the price of solving local problems on patches whose size increases (logarithmically) with the frequency. The method is not designed to allow variations of the coefficients inside the mesh cells as it adopts standard polynomials for approximability. In contrast, the proposed MHM method requires more elements per wavelength to be stable, but the base functions are designed to capture small-scale heterogeneities. Also, the local problems are element-wise and thus, less costly and easier to implement. On the other hand, a Heterogeneous Multiscale Method (HMM) is introduced in [34]. In such a work, the authors propose an efficient method to capture the solution on coarse meshes, even when the coefficient strongly vary inside mesh elements. The method is stable under a condition on the mesh size that is similar to the one obtained for the present MHM method. While the HMM proposed in [34] is limited to first-order scheme, we propose a general method to compute high order multiscale basis functions. Our work thus appear as a novel approach to obtain high order basis functions in the presence of strongly oscillating coefficients inside mesh elements.

The aim of this paper is twofold. First, we introduce a novel MHM method for the Helmholtz equation with prescribed mixed boundary conditions and verify that it can cope with the solution of highly heterogeneous wave problems. Second, in the context of homogeneous media, we give a full frequency-explicit convergence analysis of the MHM method. Notably, we prove that the MHM's solution achieves super-convergence in natural norms and is high order accurate under usual regularity assumptions. Also, we analyze the use of curved elements, which is a novelty in the context of MHM methods.

The outline of the paper is as follows: In Section 2, we state the problem and introduce the MHM method. The global-local formulation, from which the MHM method emerges, is proved to be well-posed in Section 3. Well-posedness and best approximation results for the MHM method is the subject of Section 4. Convergence results are then established assuming polynomial shape functions on faces in Section 5, followed by numerical validation in Section 6. Section 7 presents conclusions, and the appendix section some technical results.

## 2. The MHM method.

### 2.1. Statement and preliminaries.
In this work, we focus on the acoustic Helmholtz equation set in a heterogeneous medium. Let $\Omega \subset \mathbb{R}^d$, $d \in \{2, 3\}$, be a bounded Lipschitz domain with boundary $\partial\Omega$, characterized by its bulk modulus $\kappa \in L^\infty(\Omega, \mathbb{R})$ and its density $\rho \in L^\infty(\Omega, \mathbb{R})$. Considering an angular frequency $\omega \in \mathbb{R}_+^*$ and a load term $f \in L^2(\Omega)$, the pressure
$u \in H^1(\Omega)$ satisfies

$$-\frac{\omega^2}{\kappa}u - \nabla \cdot \left(\frac{1}{\rho}\nabla u\right) = f \text{ in } \Omega.$$

Having in mind seismic wave propagation, we divide the boundary $\partial\Omega$ into two subsets $\Gamma_D$ and $\Gamma_A$. $\Gamma_D$ represents the Earth surface. We impose a Dirichlet boundary condition $u = 0$ on $\Gamma_D$. The Dirichlet boundary condition corresponds to a free surface condition. A “transparent” boundary condition is prescribed on $\Gamma_A$ in order to simulate a semi-infinite propagation medium. For the sake of simplicity, we consider a first-order absorbing boundary condition on $\Gamma_A$ [14]:

$$\frac{1}{\rho}\nabla u \cdot n - \frac{1}{\sqrt{\kappa\rho}}u = 0,$$

where $n$ is the unit outward normal to $\partial\Omega$. Other transparent conditions, like high order absorbing conditions [20] or perfectly matched layers [6], can be handled by the MHM method with minor modifications. The latter option is discussed in Section 3.3.

In the remaining, we focus on the boundary value problem to find $u \in H^1(\Omega)$ solution to

$$(2.1) \qquad \left\{ \begin{array}{ll} -\frac{\omega^2}{\kappa}u - \nabla \cdot \left(\frac{1}{\rho}\nabla u\right) & = f \quad \text{in } \Omega, \\ u & = 0 \quad \text{on } \Gamma_D, \\ \frac{1}{\rho}\nabla u \cdot n - \frac{1}{\sqrt{\kappa\rho}}u & = g \quad \text{on } \Gamma_A, \end{array} \right.$$

where $g \in L^2(\Gamma_A)$ is typically employed to model an incoming plane wave.

Like the DEM [15], the MHM method is based on the primal hybrid formulation of the Helmholtz equation. In this formulation, the pressure is sought in a “broken space” of piecewise continuous functions. Thus, we introduce a mesh $\mathcal{T}_H$ of $\Omega$. Without loss of generality, we shall use here the terminology usually employed for three-dimensional domains. We assume that $\mathcal{T}_H$ perfectly partitions $\Omega$ and that each element $K \in \mathcal{T}_H$ is a (possibly curved) tetrahedron such that $\text{diam}(K) \le H$. We also assume that the family of meshes $\{\mathcal{T}_H\}_{H>0}$ is regular in the sense of [7]. By $\mathcal{F}_H$ we denote the set of the faces $F$ of $\mathcal{T}_H$.

**REMARK 1.** Exact triangulations are seldom used in practice, but isoparametric elements are usually employed instead. In this work, we consider exact triangulations only for the sake of simplicity. The authors believe that similar results can be obtained for isoparametric elements, at the price of a greater complexity in the proof to handle the “variational crime” of approximately meshing the boundary.

Consider the space

$$(2.2) \qquad V := \{v \in L^2(\Omega) \mid v|_K \in H^1(K); \forall K \in \mathcal{T}_H\}.$$

The continuity of functions in $V$ is not ensured across adjacent elements. In order to weakly enforce the continuity of the solution, we consider a space $\Lambda$ of Lagrange multipliers. This space is defined by

$$(2.3) \qquad \Lambda := \left\{ \mu \in \prod_{K \in \mathcal{T}_H} H^{-1/2}(\partial K) \mid \begin{array}{l} \exists \sigma \in H(\text{div}, \Omega); \\ \mu|_{\partial K} = \sigma \cdot n_K \quad \forall K \in \mathcal{T}_H \\ \mu|_{\Gamma_A} = 0 \end{array} \right\},$$

where $n_K$ stands for the unit outward normal to $\partial K$. We introduce the sesquilinear forms $a : V \times V \to \mathbb{C}$ and $b : \Lambda \times V \to \mathbb{C}$ by

$$(2.4) \qquad a(u, v) := \sum_{K \in \mathcal{T}_H} \left\{ -\omega^2 \int_K \frac{1}{\kappa} u \bar{v} - i\omega \int_{\partial K \cap \Gamma_A} \frac{1}{\sqrt{\kappa\rho}} u \bar{v} + \int_K \frac{1}{\rho} \nabla u \cdot \bar{\nabla} v \right\},$$
for all $u, v \in V$ and

$$ (2.5) \qquad b(\mu, v) := \sum_{K \in \mathcal{T}_h} \langle \mu, \bar{v} \rangle_{\partial K}, $$

for $\mu \in \Lambda$ and $v \in V$. Here $\langle \cdot, \cdot \rangle_{\partial K}$ is understood as the duality product between $H^{-1/2}(\partial K)$ and $H^{1/2}(\partial K)$. We also denote by $(\cdot, \cdot)_D$ the $L^2(D)$-inner product (we don't make a distinction between vector-valued and scalar-valued functions).

The primal hybrid formulation consists of finding a couple $(u, \lambda) \in V \times \Lambda$ solution to

$$ (2.6) \qquad \left\{ \begin{array}{ll} a(u, v) + b(\lambda, v) & = (f, v)_{\Omega} + (g, v)_{\Gamma_A} & \text{for all } v \in V, \\ b(\mu, u) & = 0 & \text{for all } \mu \in \Lambda. \end{array} \right. $$

The first equation of (2.6) is obtained by piecewise integration by parts of (2.1). It also corresponds to an optimality condition where $\lambda$ plays the role of the Lagrange multiplier. On the other hand, the second equation of (2.6) enforces the continuity of $u$ and is equivalent to $u \in H^1(\Omega)$ and $u = 0$ on $\Gamma_D$. The couple $(u, \lambda) \in V \times \Lambda$ solves (2.6), if and only if $u$ is a weak solution to (2.1) and $\lambda = -(1/\rho)\nabla u \cdot n_K$ on $\partial K \setminus \Gamma_A$ for all $K \in \mathcal{T}_H$. We refer the reader to [15,37] for more details about the primal hybrid formulation.

At this point, it is possible to discretize problem (2.6) directly by selecting finite-dimensional subspaces of $V$ and $\Lambda$. This approach is taken in the reference work of Raviart and Thomas [37] for the discretization of the Laplace problem. In their work, the authors derive a family of stable pairs of polynomial discretization spaces. The DEM [2, 15] is also directly based on formulation (2.6). In the standard version of DEM, the Lagrange multiplier is discretized with piecewise polynomial functions and the pressure is discretized using plane waves.

**2.2. The one-level MHM method.** The novelty of the MHM method as compared to the DEM, is to substitute the pressure $u$ for the Lagrange multiplier $\lambda$ at the continuous level. Indeed, rewriting the first equation of (2.6) as

$$ a(u, v) = (f, v) + (g, v)_{\Gamma_A} - b(\lambda, v) \quad \text{for all } v \in V, $$

we observe that the pressure can be expressed as

$$ (2.7) \qquad u = T\lambda + \hat{T}f + \tilde{T}g, $$

where $T: \Lambda \to V$, $\hat{T}: L^2(\Omega) \to V$ and $\tilde{T}: L^2(\Gamma_A) \to V$ are three linear bounded operators defined by

$$ (2.8) \qquad a(T\mu, v) = -b(\mu, v), \quad a(\hat{T}f, v) = (f, v)_{\Omega}, \quad a(\tilde{T}g, v) = (g, v)_{\Gamma_A} $$

for all $\mu \in \Lambda$, $f \in L^2(\Omega)$, $g \in L^2(\Gamma_A)$ and $v \in V$.

The analytical expression of these operators is not available in general. However, since there is no build-in compatibility condition in $V$ for $u$, the operators $T$, $\hat{T}$ and $\tilde{T}$ can be defined thanks to local boundary value problems in each element. Indeed, considering an element $K \in \mathcal{T}_H$ and test functions $v \in H^1(K)$ in (2.8), we see that

$$ (2.9) \qquad \left\{ \begin{array}{ll} -\frac{\omega^2}{\kappa} (T\mu) - \nabla \cdot \left( \frac{1}{\rho} \nabla (T\mu) \right) & = 0 \quad \text{in } K, \\ \frac{1}{\rho} \nabla (T\mu) \cdot n_K & = -\mu \quad \text{on } \partial K \setminus \Gamma_A, \\ \frac{1}{\rho} \nabla (T\mu) \cdot n & = -\frac{1}{\sqrt{\kappa\rho}} (T\mu) \quad \text{on } \partial K \cap \Gamma_A, \end{array} \right. $$
$$ (2.10) \qquad \left\{ \begin{array}{l} -\frac{\omega^2}{\kappa} (\hat{T}f|_K) - \nabla \cdot \left( \frac{1}{\rho} \nabla (\hat{T}f|_K) \right) = f & \text{in } K, \\ \frac{1}{\rho} \nabla (\hat{T}f) \cdot \mathbf{n}_K = 0 & \text{on } \partial K \setminus \Gamma_A, \\ \frac{1}{\rho} \nabla (\hat{T}f) \cdot \mathbf{n} - \frac{1}{\sqrt{\kappa \rho}} (\hat{T}f) = 0 & \text{on } \partial K \cap \Gamma_A. \end{array} \right. $$

and

$$ (2.11) \qquad \left\{ \begin{array}{l} -\frac{\omega^2}{\kappa} (\tilde{T}g|_K) - \nabla \cdot \left( \frac{1}{\rho} \nabla (\tilde{T}g|_K) \right) = 0 & \text{in } K, \\ \frac{1}{\rho} \nabla (\tilde{T}g) \cdot \mathbf{n}_K = 0 & \text{on } \partial K \setminus \Gamma_A, \\ \frac{1}{\rho} \nabla (\tilde{T}g) \cdot \mathbf{n} - \frac{1}{\sqrt{\kappa \rho}} (\tilde{T}g) = g & \text{on } \partial K \cap \Gamma_A. \end{array} \right. $$

Thus, the operators $T$, $\hat{T}$ and $\tilde{T}$ are defined locally in each element $K$ as the solutions to local boundary value problems. Assuming that the operators $T$, $\hat{T}$ and $\tilde{T}$ are available, we can substitute $u$ for $\lambda$ by plugging (2.7) in the second equation of (2.6). We obtain

$$ (2.12) \qquad b(\mu, T\lambda) = -b(\mu, \hat{T}f + \tilde{T}g) \quad \text{for all } \mu \in \Lambda. $$

In contrast to (2.9), (2.10) and (2.11), problem (2.12) is global. The one-level MHM method is obtained by introducing a finite-dimensional subspace $\Lambda_H \subset \Lambda$. The MHM formulation then consists of finding $\lambda_H \in \Lambda_H$ such that

$$ (2.13) \qquad b(\mu_H, T\lambda_H) = -b(\mu_H, \hat{T}f + \tilde{T}g) \quad \text{for all } \mu_H \in \Lambda_H, $$

and to calculate the discrete pressure by

$$ (2.14) \qquad u_H = T\lambda_H + \hat{T}f + \tilde{T}g. $$

**REMARK 2.** It is worth to mentioning that the DEM may be seen as a particular case of the MHM method. To see this clearer, assume the propagation medium is homogeneous in a two-dimensional case, and let $T_H$ be a square element mesh compound of elements $K = [0, H] \times [0, H]$ (up to a translation), and set $\Lambda_H$ as the space of functions that are constant on each face of the mesh. Then, by replacing the right-hand side of (2.9) by the corresponding basis of $\Lambda_H$, we observe that the local problem (2.9) drives the multiscale basis functions. Particularly, if one denotes those base functions by $\eta_j$ at the local level, we arrive that $\{\eta_j\}_{j=1}^4$ satisfy the following problem when it is set in the interior of the domain (i.e. $\partial K \cap \Gamma_A = \{0\}$)

$$ \left\{ \begin{array}{l} -\frac{\omega^2}{\kappa} \eta_j - \nabla \cdot \left( \frac{1}{\rho} \nabla \eta_j \right) = 0 & \text{in } K, \\ \frac{1}{\rho} \nabla \eta_j \cdot \mathbf{n}_K = -\delta_{ij} & \text{on } E_i \subset \partial K \setminus \Gamma_A, \end{array} \right. $$

where $\delta_{ij}$ is the Kronecker delta. Interestingly, a close formula to the multiscale functions $\eta_j = T\delta_{ij}$ is available for all faces $F_i \subset \partial K$. They correspond to the set of plane waves used in the lowest-order DEM to approximate the primal function $u$ at the local level, while the Lagrange multiplier $\lambda$ in the DEM method is also approximated by piecewise constant functions on the faces as in the MHM method.

**2.3. The two-level MHM method.** In the one-level MHM formulation (2.13), we assume that the operators $T$, $\hat{T}$ and $\tilde{T}$ are available and exactly computed. In practice, these operators are only analytically available in some exceptional cases. The key ingredient of the two-level MHM method is to approximate these operators using a second-level discretization scheme (for instance, Lagrange finite elements [12]). Since all computations are local, the second-level approximation is naturally parallelized and corresponds to a pre-processing step before solving the global problem.
It is thus of interest to introduce a two-level MHM method in which the operators $T, \hat{T}$ and $\tilde{T}$ are replaced by their discrete counterparts $T_h, \hat{T}_h$ and $\tilde{T}_h$ obtained thanks to a second-level method. In this case, a discretization space $V_h \subset V$ is introduced, and the discrete operators are defined by

$$a(T_h\mu, v_h) = -b(\mu, v_h), \quad a(\hat{T}_h f, v_h) = (f, v_h)_{\Omega}, \quad a(\tilde{T}_h g, v) = (g, v_h)_{\Gamma_A},$$

for all $\mu \in \Lambda$, $f \in L^2(\Omega)$ and $v_h \in V_h$. The underlying two-level MHM method corresponds to solving problem (2.13) wherein $T, \hat{T}$ and $\tilde{T}g$ are replaced by $T_h, \hat{T}_h$ and $\tilde{T}_h g$, respectively. It results in the second-level MHM solution which reads

$$(2.15) \qquad u_{H,h} := T_h \lambda_H + \hat{T}_h f + \tilde{T}_h g.$$

The definition of the operators $T_h, \hat{T}_h$ and $\tilde{T}_h$ decouples over each element $K \in \mathcal{T}_H$. As a result, evaluating $T_h\mu, \hat{T}_h f$ and $\tilde{T}_h g$ for particular arguments $\mu \in \Lambda, f \in L^2(\Omega)$ and $g \in L^2(\Gamma_A)$ amounts to solve a sequence of local element-wise discrete Galerkin problems. Owing to the flexibility of the MHM strategy, other numerical methods can replace the standard Galerkin method as the second-level solver.

**REMARK 3.** The DEM can be seen as a particular two-level MHM method. In this case, the second-level discretization space $V_h$ consists of a linear combination of plane waves in each element $K \in \mathcal{T}_H$.

**REMARK 4.** The MHM method relies on a mesh characterized by the step $H$. One of the aims of the MHM method is to account for highly heterogeneous media. As depicted in [3, 22], one asset of the MHM method is to incorporate heterogeneities of characteristic length $\epsilon \ll H$ on a coarse mesh of characteristic size $H$. In this situation, it is usually required to use a mesh of size $\epsilon$ for the second-level discretization. However, since the second-level computations are local, using a fine mesh which fits the heterogeneities (hence $h \simeq \epsilon$) is affordable.

**2.4. Practical implementation of the two-level MHM method.** We briefly sketch the main steps of the two-level MHM algorithm from a computational point of view.

1. We consider a mesh $\mathcal{T}_H$ and a finite-dimensional discretization space $\Lambda_H \subset \Lambda$ for $\lambda$ spanned by functions $(\mu_k)_{k=1}^n$, and $V_h \subset V$.

2. We compute the images of the basis functions $\mu_k$, $k \in \{1, \dots, n\}$, by the operator $T_h$ and the image of $f$ and $g$ by $\hat{T}_h$ and $\tilde{T}_h$ using $V_h$, respectively. This is done by searching functions $\eta_k = T_h\mu_k \in V_h$ for $i \in \{1, \dots, n\}$, $\eta_f = \hat{T}_h f \in V_h$ and $\eta_g = \tilde{T}_h g \in V_h$ solutions to

$$a(\eta_k, v_h) = -b(\mu_k, v_h), \quad a(\eta_f, v_h) = (f, v_h)_{\Omega}, \quad a(\eta_g, v_h) = (g, v_h)_{\Gamma_A},$$

for all $v_h \in V_h$. We emphasize that all the computations are local to each element. Such computations correspond to local (element-wise) Helmholtz problems (2.9), (2.10) and (2.11) and are a pre-processing step before solving the global problem in the next stage.

3. The MHM approximation $\lambda_H$ of $\lambda$ is defined by

$$\lambda_H = \sum_{k=1}^{n} c_k \mu_k,$$

where $c := \{c_k\}_{k=1}^n \in \mathbb{C}^n$ is a vector of unknown coefficients. We obtain $c$ by solving the $n \times n$ linear system stemming from (2.13):

$$\sum_{k=1}^{n} b(\mu_{\ell}, \eta_k) \bar{c}_k = -b(\mu_{\ell}, \eta_f + \eta_g) \quad \text{for all } \ell \in \{1, \dots, n\}.$$

The coefficient of this system simply involve boundary integral of the basis functions, that can be easily computed from their underlying second-level Galerkin's representation.
4. Compute the two-level approximate solution $u_{H,h}$ from

$$u_{H,h} = \sum_{k=1}^{n} c_k \eta_k + \eta_f + \eta_g.$$

**3. Well-posedness of the continuous MHM formulation.** In this section, we establish the well-posedness of the continuous MHM formulation (2.12). As our analysis focuses on the high-frequency case, we assume that $\omega \ge \omega_0 > 0$. We can select $\omega_0$ as small as we wish in our analysis, but the constants may “blow up” as $\omega_0 \to 0$. In order to treat general settings, we assume that (2.1) is uniquely solvable for all $\omega \ge \omega_0$, and that the stability constant growth with respect to the frequency is known. Notably, we shall assume that the following statement holds.

**ASSUMPTION 1.** For all $\omega \ge \omega_0$, $f \in L^2(\Omega)$ and $g \in L^2(\Gamma_A)$, there exists a unique $u \in H^1(\Omega)$ solution to problem (2.1). Moreover, $u$ satisfies

$$(3.1) \qquad \omega \|u\|_{0,\Omega} + |u|_{1,\Omega} \lesssim \frac{C_s(\omega)}{\omega} (\|f\|_{0,\Omega} + \|g\|_{0,\Gamma_A})$$

where $C_s(\omega)$ is a positive constant such that $C_s(\omega) \ge \omega$.

**REMARK 5.** Stability assumption (3.1) is equivalent to the inf-sup condition

$$\inf_{u \in H_{\Gamma_D}^1(\Omega)} \sup_{v \in H_{\Gamma_D}^1(\Omega)} \frac{a(u, v)}{\|u\|_{V,\omega} \|v\|_{V,\omega}} \gtrsim \frac{1}{C_s(\omega)}.$$

As shown for instance in [10, 16, 18, 25], we have $C_s(\omega) = \omega$ in the case of non-trapping geometries, but $C_s(\omega)$ can have a less favourable behaviour in other cases [32].

Above and hereafter, we lighten notation and understand the infimum (supremum) to be taken over sets excluding the zero function, even though this is not specifically indicated. Also, given two positive real numbers $A, B \ge 0$, we will employ the notation $A \lesssim B$, if there exists a constant $C > 0$ that is independent of $A, B, T_H$ and $\omega$, but that possibly depends on $\Omega, \kappa, \rho$ and $\delta$ (introduced hereafter) such that $A \le CB$. We also write $A \gtrsim B$ when $B \lesssim A$.

Next, we precise the requirements on the mesh $T_H$. Each element $K \in T_H$ is a (possibly curved) tetrahedron, and we write $H_K = \text{diam} K$. $\mathcal{F}_H$ is the set of faces of $T_H$. Also, we assume that each element $K \in T_H$ and it has at most one curved face $F$. We will establish the well-posedness of the MHM formulation under the assumption that the mesh is sufficiently refined. Specifically, we assume that there exists a constant $\delta > 0$ such that for each element $K \in T_H$, we have

$$(3.2) \qquad \frac{\omega H_K}{c_K} \le (1 - \delta)^{1/2} \pi, \quad c_K = \frac{\inf_K \kappa}{\sup_K \rho}.$$

Mesh refinement condition (3.2) states that there are at least two elements per wavelength. $c_K$ is a lower bound of the wavespeed over $K$. Actually, it is the minimum value of the wavespeed over $K$ if the inf and sup are attained at the same point in (3.2). Then $\ell_K = 2\pi c_K / \omega$ is the minimum wavelength over $K$ and we see that (3.2) means that $H_K < \ell_K / 2$.

We shall consider non-convex elements $K$, for which additional assumptions are required. The aim of the present work is not to focus on general element shapes. However, non-convex curved elements appear naturally when meshing both sides of a curved interface. In this context, let $K$ be a non-convex element such that $\partial K \cap \Gamma_A = \emptyset$, $\tilde{K}$ be the simplex formed by the convex hull of the vertices of $K \in T_H$, and $\tilde{M}: \tilde{K} \to K$ be a bijective map. We can think of $\tilde{M}$ as the “curved deformation” of $K$ from the corresponding straight element $\tilde{K}$. Then, in addition to (3.2), we assume that

$$(3.3) \qquad \frac{J_+}{J_-} B_- \ge 1 - \frac{\delta}{2} \quad \text{and} \quad B_+ \le 1 + \delta,$$
where

$$J_+ = \sup_{\tilde{K}} J, \quad J_- = \inf_{\tilde{K}} J, \quad B_+ = \sup_{\tilde{K}} \beta, \quad B_- = \inf_{\tilde{K}} \alpha,$$

with $J$ standing for the Jacobian matrix of $\tilde{M}$, $B = JJ^T$, $J = |\det J|$, and $\alpha(\tilde{x})$ and $\beta(\tilde{x})$ the smallest and highest eigenvalues of $B(\tilde{x})$, for each $\tilde{x} \in \tilde{K}$. Observe that if the mesh size is small compared to the wavelength then $\delta$ is close to one. This fact allows general curved elements to be used since restrictions on $\tilde{M}$ are mild. On the other hand, if the mesh size is close to the critical value (with respect to (3.2)), then the “non-convexity” of the element must be small since $\delta \ll 1$.

We thus assume conditions (3.2) and (3.3) for those elements $K$ that are not in contact with the absorbing boundary. On the other hand, if $\partial K \cap \Gamma_A \neq \emptyset$, then we distinguish two scenarios. First, if $\kappa$ and $\rho$ are constant in $K$, we only assume that (3.2) holds. On the contrary, when $\kappa$ and $\rho$ are not constant inside $K$, we need the more restrictive condition that

$$(3.4) \qquad \frac{\omega H_K}{c_K} \le (1 - \delta)^{1/2} \sqrt{d-1}.$$

Some comments are mandatory at that point. We observe that while the MHM formulation can become ill-posed if condition (3.2) is not fulfilled, condition (3.4) is only a technical assumption that is mandatory for our convergence analysis when the coefficients are varying in boundary cells, but not a requirement to actually apply the MHM formulation. We also stress that in the case where $\kappa$ and $\rho$ are constant in the cells in contact with $\Gamma_A$ (which is usually the case in scattering applications), our analysis fully holds under assumption (3.2) that there are at least two elements per wavelength (and assumption (3.3) that the curve deformation of the elements is not too large, if curved interfaces need to be exactly meshed).

Also, condition (3.2) is not an important constraint for low order methods, since they usually require several element per wavelength to correctly capture the solution. On the other hand, it is a severe limitation for high order methods, especially in the low- and mid-frequency regimes. In contrast, it is shown in [31] that finite element methods is well-posed under the condition that $\omega H/p$ is small enough (and the additional condition that $\log p \ge C\omega$). In particular, less than two elements per wavelength can be employed if $p$ is large enough. In the end of Section 3.1, we will describe the difficulties faced by the MHM method when (3.2) is not satisfied and propose solutions to remedy such difficulties.

In the remaining of this work, we equip the spaces $V$ and $\Lambda$ with the norms

$$\|v\|_{V,\omega}^2 := \sum_{K \in \mathcal{T}_H} \{\omega^2 \|v\|_{0,K}^2 + \omega \|v\|_{0,\partial K \cap \Gamma_A}^2 + |v|_{1,K}^2\} \quad \text{for all } v \in V,$$

and

$$(3.5) \qquad \|\mu\|_{\Lambda,\omega} := \sup_{v \in V} \frac{\operatorname{Re} b(\mu, v)}{\|v\|_{V,\omega}} \quad \text{for all } \mu \in \Lambda.$$

These norms are equivalent to the natural ones [3,37], so that $V$ and $\Lambda$ are Hilbert spaces equipped with $\|\cdot\|_{V,\omega}$ and $\|\cdot\|_{\Lambda,\omega}$. In view of (3.1), the weighted norm $\|\cdot\|_{V,\omega}$ is convenient as it “balances” the $\|\cdot\|_{0,\Omega}$ and $|\cdot|_{1,\mathcal{T}_H}$ terms of the $\|\cdot\|_{1,\mathcal{T}_H}$ norm of solutions to the Helmholtz equation. Also, the following local norm will be useful in the sequel

$$(3.6) \qquad \|v\|_{V(K),\omega}^2 := \omega^2 \|v\|_{0,K}^2 + \omega \|v\|_{0,\partial K \cap \Gamma_A}^2 + |v|_{1,K}^2 \quad \text{for all } v \in V(K),$$

where $V(K)$ stands for the space of functions in $V$ restricted to element $K \in \mathcal{T}_H$.

**3.1. Well-posedness of the local problems.** We start by proving that the local problems defining the operators $T$, $\hat{T}$ and $\tilde{T}$ are well-posed. This is classically done by showing that the sesquilinear form $a(\cdot, \cdot)$ is continuous and satisfies an inf-sup condition. The continuity of $a(\cdot, \cdot)$ is rather straightforward, and is recorded in Lemma 3.1.
LEMMA 3.1. For all $u, v \in V$, it holds that

$$|a(u, v)| \lesssim \|u\|_{V, \omega} \|v\|_{V, \omega}.$$

We require more involved arguments to demonstrate the inf-sup condition. In particular, we make use of mesh refinement conditions (3.2), (3.3) and (3.4). The inf-sup condition is constructed element-wise. Roughly speaking, the proofs rely on the fact that the local problems are coercive if $\omega H$ is small enough. It corresponds to the condition that $\omega^2$ is smaller than the Poincaré constant of the element $K$. While the key concepts are simple, we leave the proof in Appendix A as tedious computations are required to exactly obtain mesh refinements conditions (3.2), (3.3) and (3.4).

In Lemma 3.2, we establish an inf-sup condition for the sesquilinear form $a(\cdot, \cdot)$ as an easy consequence of the propositions established in Appendix A.

LEMMA 3.2. Given $u \in V$, it holds

$$\sup_{v \in V} \frac{\mathrm{Re} a(u, v)}{\|v\|_{V, \omega}} \gtrsim \|u\|_{V, \omega}.$$

*Proof.* For each element $K \in \mathcal{T}_H$, let us pick $u^* \in H^1(K)$ like in the propositions of Appendix A. We define $v \in V$ such that $v|_K = u^*$ for all $K \in \mathcal{T}_H$. Then, we have

$$\mathrm{Re} a(u, v) = \sum_{K \in \mathcal{T}_H} \mathrm{Re} a(u, v) = \sum_{K \in \mathcal{T}_H} \mathrm{Re} a(u, u^*) \gtrsim \sum_{K \in \mathcal{T}_H} \|u\|_{V, \omega}^2 = \|u\|_{V, \omega}^2,$$

□

and the result follows, since one easily sees that $\|v\|_{V, \omega} \lesssim \|u\|_{V, \omega}$.

We are now ready to establish the well-posedness of the local problems in Theorem 3.3.

THEOREM 3.3. For all $\mu \in \Lambda$, $f \in L^2(\Omega)$ and $g \in L^2(\Gamma_A)$, there exist unique elements $T\mu, \hat{T}f, \tilde{T}g \in V$ such that

$$a(T\mu, v) = -b(\mu, v), \quad a(\hat{T}f, v) = (f, v)_{\Omega}, \quad a(\tilde{T}g, v) = (g, v)_{\Gamma_A} \quad \text{for all } v \in V.$$

Furthermore, we have

$$(3.7) \qquad \|T\mu\|_{V, \omega} \lesssim \|\mu\|_{\Lambda, \omega}, \qquad \|\hat{T}f\|_{V, \omega} \lesssim \omega^{-1} \|f\|_{0, \Omega}, \qquad \|\tilde{T}g\|_{V, \omega} \lesssim \omega^{-1/2} \|g\|_{0, \Gamma_A}.$$

*Proof.* The continuity of the sesquilinear form $a(\cdot, \cdot)$ established in Lemma 3.1 together with the inf-sup condition of Lemma 3.2 ensure the existence and uniqueness of $T\lambda$ and $\hat{T}f$ for all $\lambda \in \Lambda$ and $f \in L^2(\Omega)$.

Let us show (3.7). Consider $\lambda \in \Lambda$. Then for $v \in V \setminus \{0\}$, we have

$$\frac{\mathrm{Re} a(T\lambda, v)}{\|v\|_{V, \omega}} = -\frac{\mathrm{Re} b(\lambda, v)}{\|v\|_{V, \omega}} \le \|\lambda\|_{\Lambda, \omega},$$

and the first estimate of (3.7) follows from Lemma 3.2. Now, if $f \in L^2(\Omega)$, we have

$$\mathrm{Re} a(\hat{T}f, v) = \mathrm{Re}(f, v)_{\Omega} \le \|f\|_{0, \Omega} \|v\|_{0, \Omega} \le \omega^{-1} \|f\|_{0, \Omega} \|v\|_{V, \omega} \quad \text{for all } v \in V$$

so that

$$\|\hat{T}f\|_{V, \omega} \lesssim \sup_{v \in V} \frac{\mathrm{Re} a(\hat{T}f, v)}{\|v\|_{V, \omega}} \le \omega^{-1} \|f\|_{0, \Omega},$$

and the second estimate of (3.7) follows. Likewise, for $g \in L^2(\Gamma_A)$, we have

$$\mathrm{Re} a(\tilde{T}g, v) = \mathrm{Re}(g, v)_{\Gamma_A} \le \|g\|_{0, \Gamma_A} \|v\|_{0, \Gamma_A} \le \omega^{-1/2} \|g\|_{0, \Gamma_A} \|v\|_{V, \omega} \quad \text{for all } v \in V,$$
thus

$$ \|\tilde{T}g\|_{V,\omega} \lesssim \sup_{v \in V} \frac{\mathrm{Re} a(\tilde{T}g, v)}{\|v\|_{V,\omega}} \lesssim \omega^{-1/2} \|g\|_{0,\Gamma_A}, $$

and we obtain the last estimate of (3.7).

COROLLARY 3.4. The applications $\| \cdot \|_{\Lambda, \omega}$ and $\|T \cdot \|_{V, \omega}$ define equivalent norms on $\Lambda$, i.e.,
$$ (3.8) \qquad \|\mu\|_{\Lambda, \omega} \lesssim \|T\mu\|_{V, \omega} \lesssim \|\mu\|_{\Lambda, \omega} \quad \text{for all } \mu \in \Lambda. $$

*Proof.* From Theorem 3.3, the upper bound of (3.8) is already established. To prove the lower bound, we only need to observe that

$$ \|\mu\|_{\Lambda, \omega} = \sup_{v \in V} \frac{\mathrm{Re} b(\mu, v)}{\|v\|_{V, \omega}} = \sup_{v \in V} \frac{\mathrm{Re} a(T\mu, v)}{\|v\|_{V, \omega}} \lesssim \|T\mu\|_{V, \omega}, $$

by Lemma 3.1.

We have employed mesh refinement conditions (3.2), (3.3) and (3.4) to establish the inf-sup condition for the sesquilinear form $a(\cdot, \cdot)$, and show that the operators $T$, $\tilde{T}$ and $\tilde{T}$ are well-defined. Here, we give more insight into what happens if these conditions are not satisfied. First, the local problems associated with the outer element $(\partial K \cap \Gamma_A \neq \emptyset)$ are always well-posed, since they correspond to Helmholtz problems with an absorbing boundary. However, the inf-sup constant (and thus, the operator norms of $T$, $\tilde{T}$ and $\tilde{T}$) will depend on the number of wavelength inside $K$, and thus on $H_K$ and $\omega$.

The case of interior elements $(\partial K \cap \Gamma_A = \emptyset)$ is more subtle. In this case, the local problem corresponds to a Helmholtz problem without absorption, and there is a discrete set of “resonance frequencies” for which the problem is not well-posed. Such resonance frequencies correspond to the eigenvalues $\tilde{\omega}^2$ of

$$ (\rho^{-1}\nabla u, \nabla v)_K = \tilde{\omega}^2 (\kappa^{-1}u, v)_K. $$

To give more details, we assume that the medium is homogeneous $(\kappa = \rho = 1)$. In this case, for a fixed element shape, we can show using a scaling argument that the discrete set of “resonance frequencies” corresponds to a discrete set of mesh sizes $H_K$ for which the local problem is not well-posed. Assuming that $K$ is square of size $H_K$, we can be more specific: the eigenvalues are the set $(j^2\pi^2H^{-2})_{j \in \mathbb{N}}$. Thus, the local problem is well-posed if and only if $\omega H \neq j\pi$, $j \in \mathbb{N}$. Condition (3.2) then states that $\omega H$ is between the first two forbidden values $0$ and $\pi$. We then see that the MHM formulation is well-posed as soon as $\omega H \neq j\pi$, but in this case, the operator norms of $T$, $\tilde{T}$ and $\tilde{T}$ will depend on the distance between $\omega H/\pi$ and the closest integer. Also, though constructing meshes that avoids “forbidden sizes” might be feasible for square elements in homogeneous media, it is probably not practically possible for general simplicial elements or heterogeneous media.

We propose three different ways to overcome the difficulty mentioned above, although we do not pursue their analysis in this work. The first idea uses concepts adopted in MHM formulations to other PDEs. The key solution is to “extract the kernel” of the local problems, as shown in [3] for the Laplace problem (that correspond to the particular case $\omega = 0$ here). For the particular case of the Laplace operator, this kernel corresponds to constant functions. Such constants are “removed” from the local problems and computed using the global problem. For the Helmholtz equation, the similar idea would be to extract the eigenfunctions corresponding to the $m$-first eigenvalues of the local problem.

The second technique consists of introducing artificial dissipation in the local problem. Specifically, we replace the operator $T$ by a perturbed operator $T_\epsilon$ define as the solution to

$$ (3.9) \qquad \begin{cases} -\omega^2 T \mu - i\epsilon T \mu - \Delta \mu = 0 & \text{in } K, \\ \nabla(T \mu) \cdot n_K = \mu & \text{on } \partial K. \end{cases} $$
A similar perturbation is applied to $\hat{T}$ and $\tilde{T}$. We can easily show that (3.9) is always well-posed when $\varepsilon > 0$. Also, if the problem is well-posed for $\varepsilon = 0$, then the operator norms of $T_\varepsilon$, $\hat{T}_\varepsilon$ and $\tilde{T}_\varepsilon$ converge to the norms of $T$, $\hat{T}$ and $\tilde{T}$ as $\varepsilon \to 0$. On the other hand, if the original problem is not well-posed the operator norms of $T_\varepsilon$, $\hat{T}_\varepsilon$ and $\tilde{T}_\varepsilon$ “blow up” as $\varepsilon \to 0$, and the local problems become very difficult to solve numerically. Nevertheless, we can show that the solution to the perturbed MHM formulation corresponds to the solution $u_\varepsilon$ to

$$-\omega^2 u_\varepsilon - i\varepsilon u_\varepsilon - \Delta u_\varepsilon = f \text{ in } \Omega,$$

that converges to $u$. Hence, selecting $\varepsilon$ small enough (typically, in $O(\omega H)$), the MHM formulation is well-posed for any mesh and we obtain an accurate solution.

The third method employs a modified Lagrange multiplier $\lambda \in \Lambda$. The idea is to select one face $F = \partial K_+ \cap \partial K_-$ per couple of elements $K_\pm \in \mathcal{T}_H$ and change the definition of $\Lambda$ and change the original definition of $\lambda$ ($\lambda = -\frac{1}{\rho} \nabla u \cdot n_{K_\pm}$ on each side of $F$) by $\lambda = -\frac{1}{\rho} \nabla u \cdot n_{K_\pm} \pm i\omega u$ on each side of $F$. This modification of $\lambda$ introduces an absorbing boundary condition on one face of each local problem so that they are always well-posed. The main drawback of the method is that a particular partition of the mesh must be constructed to select the faces on which the definition of $\lambda$ should be modified. Furthermore, though we can easily construct such partition for particular mesh topologies, it is not clear that it is possible for general unstructured meshes.

**3.2. Well-posedness of the global MHM problem.** We now prove that the MHM global formulation is well-posed. We start by showing existence and uniqueness of the solution in the following theorem.

**THEOREM 3.5.** There exists a unique solution $\lambda \in \Lambda$ solution to (2.12). Furthermore, it holds that

$$ (3.10) \qquad \|\lambda\|_{\Lambda, \omega} \lesssim \frac{C_s(\omega)}{\omega} (\|f\|_{0, \Omega} + \|g\|_{0, \Gamma_A}). $$

*Proof.* From Assumption 1, we know that there exists a couple $(u, \lambda) \in V \times \Lambda$ such that

$$ (3.11) \qquad \left\{ \begin{array}{ll} a(u, v) + b(\lambda, v) & = (f, v)_\Omega + (g, v)_{\Gamma_A} & \text{for all } v \in V, \\ b(\mu, u) & = 0 & \text{for all } \mu \in \Lambda, \end{array} \right. $$

where $u \in H^1(\Omega)$ is the usual solution to the Helmholtz equation and $\lambda$ is defined as the normal derivative of $u$ on the boundary of each $K \in \mathcal{T}_H$. On the other hand, Theorem 3.3 states that the operators $T$, $\hat{T}$ and $\tilde{T}$ are well-defined and invertible. As a result, the first equation of (3.11) shows that

$$ (3.12) \qquad u = T\lambda + \hat{T}f + \tilde{T}g. $$

Injecting (3.12) into the second equation of (3.11) shows that $\lambda$ is the solution to the continuous MHM formulation (2.12), and existence follows. Also, we have uniqueness, as the couple $(u, \lambda)$ is unique.

As for (3.10), observe that

$$ \begin{aligned} \|T\lambda\|_{V, \omega} &\le \|u\|_{V, \omega} + \|\hat{T}f\|_{V, \omega} + \|\tilde{T}g\|_{V, \omega} \\ &\lesssim \left( \frac{C_s(\omega)}{\omega} + \omega^{-1} + \omega^{-1/2} \right) (\|f\|_{0, \Omega} + \|g\|_{0, \Gamma_A}) \\ &\lesssim \frac{C_s(\omega)}{\omega} (\|f\|_{0, \Omega} + \|g\|_{0, \Gamma_A}), \quad \Box \end{aligned} $$

where we have used the fact that $C_s(\omega) \ge \omega$ and $\omega \gtrsim 1$. The result follows from the norm equivalence (3.8).

The next results are devoted to the analysis of the MHM sesquilinear form $\Lambda \ni \mu, \lambda \to b(\mu, T\lambda)$. The next proposition addresses a symmetry result.
PROPOSITION 3.6. For all $\mu, \lambda \in \Lambda$, we have

$$ (3.13) \qquad b(\mu, T\lambda) = b(\bar{\lambda}, T\bar{\mu}). $$

*Proof.* Let $\mu, \lambda \in \Lambda$. We have

$$ \overline{b(\mu, T\lambda)} = b(\bar{\mu}, \overline{T\lambda}) = -a(T\bar{\mu}, \overline{T\lambda}) = -a(T\lambda, \overline{T\bar{\mu}}) = b(\lambda, \overline{T\bar{\mu}}), $$

and (3.13) follows by taking the complex conjugate.

The next lemma associates to $\lambda \in \Lambda$ an element $\eta_\lambda \in \Lambda$ that plays a crucial role in the derivation of the inf-sup condition for the sesquilinear form $b(\cdot, \cdot)$.

LEMMA 3.7. For $\lambda \in \Lambda$, define $\eta_\lambda \in \Lambda$ as the unique solution to

$$ (3.14) \qquad b(\mu, T\overline{\eta_\lambda}) = -b(\mu, \hat{T}(\overline{T\lambda})) \quad \text{for all } \mu \in \Lambda. $$

Then, we have

$$ (3.15) \qquad b(\eta_\lambda, T\lambda) = \|T\lambda\|_{0,\Omega}^2, $$

and

$$ (3.16) \qquad \|T\eta_\lambda\|_{V,\omega} \lesssim \frac{C_s(\omega)}{\omega} \|T\lambda\|_{0,\Omega}. $$

*Proof.* Consider $\lambda \in \Lambda$. The existence and uniqueness result of Theorem 3.5 ensures that definition (3.14) makes sense for $\eta_\lambda$. Then, from Proposition 3.6, it holds that

$$ (3.17) \qquad b(\eta_\lambda, T\lambda) = b(\bar{\lambda}, T\overline{\eta_\lambda}) = -b(\bar{\lambda}, \hat{T}(\overline{T\lambda})). $$

Next, using the definitions of $T$ and $\hat{T}$, we show that

$$ (3.18) \qquad b(\bar{\lambda}, \hat{T}(\overline{T\lambda})) = -b(\bar{\lambda}, \overline{T(\overline{T\lambda})}) = a(T\lambda, \overline{T(\overline{T\lambda})}) = a(\hat{T}(\overline{T\lambda}), \overline{T\lambda}) = (\overline{T\lambda}, \overline{T\lambda})_\Omega = \|\overline{T\lambda}\|_{0,\Omega}^2. $$

Taking the complex conjugate of (3.18), we conclude the demonstration of (3.14). Indeed, since the right-hand side of (3.18) is real, we have

$$ -b(\bar{\lambda}, \hat{T}(\overline{T\lambda})) = \overline{\|\overline{T\lambda}\|_{0,\Omega}^2} = \|\overline{T\lambda}\|_{0,\Omega}^2 = \|T\lambda\|_{0,\Omega}^2. $$

Finally, (3.16) follows from the definition of $\eta_\lambda$ and Theorem 3.5. $\square$

We close this section with an inf-sup condition for the MHM sesquilinear form $b(\cdot, \cdot)$. Essentially, Theorem 3.8 establishes that the inf-sup condition of the MHM formulation has the same frequency behavior as the one of the original problem.

THEOREM 3.8. For all $\lambda \in \Lambda$, it holds

$$ (3.19) \qquad \sup_{\mu \in \Lambda} \frac{\operatorname{Re} b(\mu, T\lambda)}{\|\mu\|_{\Lambda, \omega}} \ge \frac{1}{C_s(\omega)} \|\lambda\|_{\Lambda, \omega}. $$

*Proof.* Fix $\lambda \in V$. We have

$$ -\operatorname{Re} b(\lambda, T\lambda) = \operatorname{Re} a(T\lambda, T\lambda) \ge A|T\lambda|_{1, T_H}^2 - B\omega^2 \|T\lambda\|_{0, \Omega}^2, $$

and

$$ -\operatorname{Re} b(i\lambda, T\lambda) = \operatorname{Re} a(T(i\lambda), T\lambda) = \operatorname{Re} ia(T\lambda, T\lambda) \ge C\omega \|T\lambda\|_{0, \Gamma_A}^2, $$
where the constants A, B and C only depend on $κ$ and $ρ$. As a result, we have

$$-\mathrm{Re} b((1+i)\lambda, T\lambda) \geq A' \|T\lambda\|_{V,\omega}^2 - B' \omega^2 \|T\lambda\|_{0,\Omega}^2,$$

where $A'$ and $B'$ only depends on $κ$ and $ρ$. Next, we define $η_λ \in Λ$ as in Lemma 3.7, and we have

$$-\mathrm{Re} b(\eta_λ, T\lambda) = \|T\lambda\|_{0,\omega}^2.$$

We can thus define $μ = (1+i)λ + B'ω^2η_λ$, and obtain

$$-\mathrm{Re} b(μ, Tλ) \geq A' \|Tλ\|_{V,ω}^2 \gtrsim \|Tλ\|_{V,ω}^2.$$

Hence, it remains to show that $\|Tμ\|_{V,ω} \lesssim C_s(ω)\|Tλ\|_{V,ω}^2$. But in view of (3.16), we have

$$\|Tμ\|_{V,ω} \lesssim \|Tλ\|_{V,ω} + ω^2 \|Tη_λ\|_{V,ω} \lesssim \|Tλ\|_{V,ω} + ω^2 \frac{C_s(ω)}{ω} \|Tλ\|_{0,ω} \\ \lesssim \|Tλ\|_{V,ω} + C_s(ω) \|Tλ\|_{V,ω} \lesssim (1 + C_s(ω)) \|Tλ\|_{V,ω}$$

and (3.19) follows since $C_s(ω) \geq ω \geq 1$.

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

**3.3. Well-posedness with perfectly matched layers.** For the sake of simplicity, the present work focuses on absorbing boundary conditions. However, the proposed MHM method applies to perfectly matched layers (PML) with very minor modifications. Indeed, in this case, the governing equation reads

$$\left\{ \begin{array}{ll} -\frac{\omega^2}{\kappa} du - \nabla \cdot \left( \frac{1}{\rho} D \nabla u \right) & = f \quad \text{in } \Omega, \\ u & = 0 \quad \text{on } \partial \Omega, \end{array} \right.$$

where the scalar and matrix functions $d$ and $D$ equal to 1 and $I$, respectively, in the “region of interest” $\Omega_0$ and take complex values in the PML region $\Omega \setminus \Omega_0$ (see [6] for details).

As a result, the derivation of the MHM formulation follows the same steps as in the absorbing boundary case. Here, we face two types the local problems. From one side, the “interior” local problems for elements $K \subset \Omega_0$ are the same as in the case of absorbing boundary conditions. On the other hand, the non-zero imaginary parts of $D$ and $d$ in the “PML” local problems (when $K \subset \Omega \setminus \Omega_0$) yield the coercivity of the associated weak form, which implies the well-posedness of the local problems.

**4. Well-posedness and convergence of the MHM method.** In this section, we assume that the (local) second-level approximations involved in the two-level MHM method (see Section 2.3) are “accurate enough” such that their underlying errors may be disregarded. As such, we analyze the one-level MHM method proposed in Section 2.2, where the space $Λ$ is replaced by an internal approximation subspace $Λ_H \subset Λ$. We show that if the elements of $Λ_H$ can accurately reproduce continuous solutions to (2.12) then discrete problem (2.13) is well-posed, and the discrete solution is quasi-optimal.

We derive a general theory under the mere assumption that $Λ_H$ is a finite-dimensional subspace of $Λ$. An application to the particular case of spaces $Λ_H$ build on polynomials is discussed later in Section 5. Our analysis is based on the so-called “Shatz argument” [39]. Specifically, following [31], we introduce the real number $α_{ω,H}$ in order to characterize the approximation properties of $Λ_H$. It is defined by

$$ (4.1) \qquad \alpha_{\omega,H} = \sup_{f \in L^2(\Omega)} \inf_{\mu_H \in \Lambda_H} \frac{\|\lambda_f - \mu_H\|_{\Lambda,\omega}}{\|f\|_{0,\Omega}}, $$

where $λ_f$ is the unique element of $Λ$ such that $b(μ, Tλ_f) = -b(μ, \hat{T}f)$ for all $μ \in Λ$. Since $Λ_H$ is finite-dimensional, as a direct consequence of (4.1), for all $f \in L^2(\Omega)$, there exists an element $μ_H \in Λ_H$ such that

$$ (4.2) \qquad \|\lambda_f - \mu_H\|_{\Lambda,\omega} \leq \alpha_{\omega,H} \|f\|_{0,\Omega}. $$
We follow the path of [31] to derive an inf-sup condition for the discrete problem as well as error estimates. Though the key ideas are the same as in [31], several slight modifications are required to adapt the proof (performed having standard finite element methods in mind) to the MHM method. We recall that the results of this section are established under assumption (3.2) that the mesh contains at least two elements per wavelength, and the additional assumption that (3.3) and (3.4) respectively hold for non-convex elements and for the boundary elements if $\kappa$ and $\rho$ are not constant close to $\Gamma_A$.

**THEOREM 4.1.** Assume that $\omega\alpha_{\omega,H}$ is small enough. Then, it holds that

$$ (4.3) \qquad \sup_{\mu_H \in \Lambda_H} \frac{\mathrm{Re} b(\mu_H, T\lambda_H)}{\|\mu_H\|_{\Lambda,\omega}} \gtrsim \frac{1}{C_s(\omega)} \|\lambda_H\|_{\Lambda,\omega}, $$

for all $\lambda_H \in \Lambda_H$.

*Proof.* We fix an arbitrary $\lambda_H \in \Lambda_H$. Recalling the proof of Theorem 3.8, we have

$$ -\mathrm{Re} b(\mu, \lambda_H) \gtrsim \|T\lambda_H\|_{V,\omega}^2, $$

with $\mu = (1+i)\lambda_H + B\omega^2\eta_{\lambda_H}$ where the constant $B$ only depends on $\kappa$ and $\rho$. Then, we define $\mu_H \in \Lambda_H$ as

$$ \mu_H = (1+i)\lambda_H + B\omega^2\eta_H, $$

where $\eta_H$ is the best approximation of $\eta_{\lambda_H}$. It follows that

$$ \mu - \mu_H = B\omega^2(\eta_{\lambda_H} - \eta_H) $$

and recalling the definitions of $\eta_{\lambda_H}$ from Lemma 3.7, and the definition of $\alpha_{\omega,H}$ we see that

$$ (4.4) \qquad \|T(\mu - \mu_H)\|_{V,\omega} \lesssim \omega^2 \|T(\eta_{\lambda_H} - \eta_H)\|_{V,\omega} \lesssim \omega^2 \alpha_{\omega,H} \|T\lambda_H\|_{0,\Omega} \lesssim \omega\alpha_{\omega,H} \|T\lambda_H\|_{V,\omega}. $$

Thus, we have

$$ -\mathrm{Re} b(\mu_H, T\lambda_H) = -\mathrm{Re} b(\mu, T\lambda_H) + \mathrm{Re} b(\mu - \mu_H, T\lambda_H) \\ \gtrsim \|T\lambda_H\|_{V,\omega}^2 - \|T(\mu - \mu_H)\|_{V,\omega} \|T\lambda_H\|_{V,\omega} \gtrsim (1 - \omega\alpha_{\omega,H}) \|T\lambda_H\|_{V,\omega}^2 \gtrsim \|T\lambda_H\|_{V,\omega}^2 $$

under the assumption that $\omega\alpha_{\omega,H}$ is small enough. Thus, it remains to show that

$$ \|T\mu_H\|_{V,\omega} \lesssim C_s(\omega) \|T\lambda_H\|_{V,\omega}. $$

But, recalling Lemma 3.7 and (4.4), we have

$$ \|T\mu_H\|_{V,\omega} \lesssim \|T\mu\|_{V,\omega} + \|T(\mu - \mu_H)\|_{V,\omega} \lesssim (\omega\alpha_{\omega,H} + C_s(\omega)) \|T\lambda_H\|_{V,\omega} \lesssim C_s(\omega) \|T\lambda_H\|_{V,\omega}, $$

assuming $\omega\alpha_{\omega,H}$ is small enough, and (4.3) follows. $\square$

Theorem 4.1 shows that the one-level MHM method is well-posed as soon as the discretization subspace $\Lambda_H$ reproduces solutions to (2.13) accurately. This fact corresponds to the condition that $\omega\alpha_{\omega,H}$ must be sufficiently small. Under this condition, there exists a unique $\lambda_H \in \Lambda_H$ solution to (2.13).

In the remaining of this section, we derive an error-estimate for the discrete solution. Though the error-estimate can arise from the inf-sup condition in Theorem 4.1, we perform an additional analysis to obtain a sharper bound. We start by showing an Aubin-Nitsche type inequality in Lemma 4.2.

**LEMMA 4.2.** Let $\lambda \in \Lambda$ solve (2.12) and $\lambda_H \in \Lambda_H$ satisfy (2.13). It holds that

$$ (4.5) \qquad \|T(\lambda - \lambda_H)\|_{0,\Omega} \lesssim \alpha_{\omega,H} \|T(\lambda - \lambda_H)\|_{V,\omega}. $$
*Proof.* We use again Lemma 3.7, and define $\eta \in \Lambda$ as

$$b(\mu, T\bar{\eta}) = -b(\mu, \hat{T}(\bar{T}(\lambda - \lambda_H)))$$

so that we have

$$b(\eta, T(\lambda - \lambda_H)) = \|T(\lambda - \lambda_H)\|_{0,\Omega}^2.$$

Then, by Galerkin's orthogonality, the definition of the norm $\| \cdot \|_{\Lambda, \omega}$, and the continuity of the application $a(\cdot, \cdot)$ in Lemma 3.1, it is clear that

$$(4.6) \quad \begin{aligned} \|T(\lambda - \lambda_H)\|_{0,\Omega}^2 &= b(\eta, T(\lambda - \lambda_H)) = b(\eta - \bar{\eta}_H, T(\lambda - \lambda_H)) = b(\bar{\eta} - \eta_H, \bar{T}(\lambda - \lambda_H)) \\ &= -a(T(\bar{\eta} - \eta_H), \bar{T}(\lambda - \lambda_H)) \lesssim \|T(\bar{\eta} - \eta_H)\|_{V,\omega} \|T(\lambda - \lambda_H)\|_{V,\omega}, \end{aligned}$$

for all $\eta_H \in \Lambda_H$. By definition of $\alpha_{\omega,H}$ in (4.1), recalling (4.2), there exists a $\eta_H \in \Lambda_H$ such that

$$(4.7) \quad \|T(\bar{\eta} - \eta_H)\|_{V,\omega} \lesssim \alpha_{\omega,H} \|\bar{T}(\lambda - \lambda_H)\|_{0,\Omega} = \alpha_{\omega,H} \|T(\lambda - \lambda_H)\|_{0,\Omega},$$

and (4.5) follows from (4.6) and (4.7). $\square$

We are now ready to establish our main convergence result for the one-level MHM method.

**THEOREM 4.3.** Assume that $\omega\alpha_{\omega,H}$ is small enough. Then, there exists a unique $\lambda_H \in \Lambda_H$ solution to (2.13), and we have

$$(4.8) \quad \|\lambda - \lambda_H\|_{\Lambda, \omega} \lesssim \inf_{\mu_H \in \Lambda_H} \|\lambda - \mu_H\|_{\Lambda, \omega}.$$

*Proof.* By definition of $a(\cdot, \cdot)$ and $b(\cdot, \cdot)$, we have that

$$\begin{aligned} - \operatorname{Re} b((\lambda - \lambda_H), (1+i)T(\lambda - \lambda_H)) &= \operatorname{Re} a(T(\lambda - \lambda_H), T(\lambda - \lambda_H)) - \operatorname{Im} a(T(\lambda - \lambda_H), T(\lambda - \lambda_H)) \\ &\gtrsim \|T(\lambda - \lambda_H)\|_{V,\omega}^2 - \omega^2 \|T(\lambda - \lambda_H)\|_{0,\Omega}^2. \end{aligned}$$

Recalling Lemma 4.2, we have

$$\omega^2 \|T(\lambda - \lambda_H)\|_{0,\Omega}^2 \lesssim \omega^2 \alpha_{\omega,H}^2 \|T(\lambda - \lambda_H)\|_{V,\omega}^2,$$

and it follows that

$$- \operatorname{Re} b(\lambda - \lambda_H, T(\lambda - \lambda_H)) \gtrsim (1 - \omega^2 \alpha_{\omega,H}^2) \|T(\lambda - \lambda_H)\|_{V,\omega}^2.$$

Assuming that $\omega\alpha_{\omega,H}$ is small enough, we obtain that

$$(4.9) \quad - \operatorname{Re} b(\lambda - \lambda_H, T(\lambda - \lambda_H)) \gtrsim \|T(\lambda - \lambda_H)\|_{V,\omega}^2.$$

We can now end the proof of error-estimate (4.8) using the Galerkin's orthogonality. Indeed, using Galerkin's orthogonality in (4.9), it holds that

$$\|T\lambda - T\lambda_H\|_{V,\omega}^2 \lesssim |b(\lambda - \lambda_H, T(\lambda - \lambda_H))| \lesssim |b(\lambda - \mu_H, T(\lambda - \lambda_H))| \lesssim \|\lambda - \mu_H\|_{\Lambda,\omega} \|T(\lambda - \lambda_H)\|_{V,\omega},$$

for all $\mu_H \in \Lambda_H$, and (4.8) follows from norm equivalence (3.8). $\square$

**5. Convergence with polynomial discretizations.** In this section, we assume that $\Gamma_D$ and $\Gamma_A$ are disjoint, and we assume that $\kappa$ and $\rho$ are piecewise smooth. Specifically, as depicted in Figure 5.1, we assume that there exists a partition $\mathcal{P} = \{\Omega_p\}_{p=1}^P$ of $\Omega$ such that for each $p \in \{1, \dots, P\}$, $\Omega_p$ has a smooth boundary $\partial\Omega_p$ of class $C^{m+1,1}$, and $\kappa|_{\Omega_p}, \rho|_{\Omega_p} \in C^{m,1}(\bar{\Omega}_p)$, where $m \ge 0$ is a positive integer. The following regularity theorem is established in [10].
<figure><img src="image_3.png" /><figcaption>FIG. 5.1. Example of a partition $\mathcal{P}$.</figcaption></figure>

**THEOREM 5.1.** For every $\omega \ge \omega_0$, there exists a unique $u \in H^1(\Omega)$ solution to (2.1) with $g=0$. Furthermore, for every $\ell \le m$, the exact solution $u$ admits the following splitting<fnref n="1" />

$$ (5.1) \qquad u = \sum_{j=0}^{\ell-1} \omega^j u_j + r_\ell, $$

where $u_j, r_\ell \in H_{\Gamma_D}^1(\Omega)$. In addition, for each $p \in \{1, \dots, P\}$, we have $u_j \in H^{j+2}(\Omega_p)$, $r_\ell \in H^{\ell+2}(\Omega_p)$ and

$$ \|u_j\|_{j+2, \Omega_p} \lesssim \|f\|_{0, \Omega}, \quad \|r_\ell\|_{\ell+2, \Omega_p} \lesssim C_s(\omega)\omega^\ell \|f\|_{0, \Omega}. $$

In this section, we consider discretization spaces that are built on polynomials. Such spaces are extensively described in [37] for the case of straight elements. Here, as we consider smooth boundaries of class $C^{m+1,1}$, we introduce slight modifications to consider curved elements. Specifically, we allow curved faces on the boundaries $\partial\Omega_p$ ($1 \le p \le P$). The other faces of the mesh are straight. For the sake of simplicity, we assume that $\Omega$ is exactly triangulated by a mesh $\mathcal{T}_H$ that is regular of class $C^m$ in the sense of [7]. We also assume that mesh $\mathcal{T}_H$ is compatible with the partition $\{\Omega_p\}_{p=1}^P$. Specifically, we assume that for each $K \in \mathcal{T}_H$, there exists a $p \in \{1, \dots, P\}$ such that $K \subset \bar{\Omega}_p$.

**REMARK 6.** It is also possible to consider polygonal domains and straight elements, if the mesh is properly refined close to the corners and edges of the geometry, where the solution can exhibit a singular behaviour [10]. For the sake of simplicity though, we do not consider such configuration here.

Each element $K \in \mathcal{T}_H$ is obtained from a single reference tetrahedron $\hat{K}$ through a smooth invertible mapping $M_K$. We employ a usual definition for the discretization space, namely,

$$ (5.2) \qquad \Lambda_H := \{ \mu_H \in \Lambda : \mu_H = \hat{p} \circ M_K^{-1} \text{ on } F; \hat{p} \in \hat{P}_\ell(\hat{K}); \forall K \in \mathcal{T}_H; \forall F \in \mathcal{F}_H(K) \}, $$

where $\hat{P}_\ell(D)$ stands for the space of polynomial functions of degree less or equal than $\ell$ on an open set $D$ and $0 \le \ell \le m$ is a fixed integer, and $\mathcal{F}_H(K)$ stands for the faces of the element $K$. We stress that since the condition $\mu_H = \hat{p} \circ \mathcal{F}_K^{-1}$ is only imposed face-wise (and not over the whole boundary $\partial K$), functions of $\Lambda_H$ are continuous on faces, but do not satisfy any compatibility conditions at the nodes and edges of the mesh. Furthermore, in the case of straight elements (with affine mappings $\mathcal{F}_K$), $\mu_H$ is simply a polynomial function over each face $F$ in the two-dimensional coordinate system of $F$.

For each face $\hat{F}$ of the reference element $\hat{K}$, we define the face interpolant $\pi_{\hat{F}}\hat{\mu} \in L^2(\hat{F})$ as the unique element of $\hat{P}_\ell(\hat{F})$ such that

$$ \int_{\hat{F}} \pi_{\hat{F}} \hat{\mu} q = \int_{\hat{F}} \hat{\mu} q \quad \text{for all } q \in \hat{P}_\ell(\hat{F}). $$

<footnote n="1">¹We take the convention that the $u = r_0$ when $\ell = 0$.</footnote>
If $\mu \in \Lambda$ satisfies $\mu|_{\partial K} \in L^2(\partial K)$ for all $K \in \mathcal{T}_H$, we define its interpolant $\pi_H \mu \in \Lambda_H$ on each face $F \in \mathcal{F}_H$ through $L^2(F)$-projection. Specifically, we set

$$ (\tilde{\pi}_H \mu)|_F := \pi_F \hat{\mu}. $$

Above and hereafter, we employ the notation $\hat{v} := v \circ \mathcal{M}_K$ for any function $v$. Lemma 5.2 states that the interpolant converges with the expected rate if the interpolated function is sufficiently smooth. The result is standard for straight elements [37], and is extended here for the curved case. The proof is given in the appendix, as it is not difficult, but technical. Although we only consider the two-dimensional case here, three-dimensional results can be obtained at the price of more technicalities in the proof.

**LEMMA 5.2.** Let $\phi \in H^{j+2}(\Omega_p)$ for all $p \in \{1, \dots, P\}$ for some $j \in \mathbb{N}$, and define $\mu \in \Lambda$ as

$$ \mu|_{\partial K \setminus \Gamma_A} := \frac{1}{\rho} \nabla \phi \cdot n_K, \quad \mu|_{\partial K \cap \Gamma_A} := 0, $$

for all $K \in \mathcal{T}_H$. Then, we have

$$ \|\mu - \pi_H \mu\|_{\Lambda, \omega} \lesssim H^{r+1} \sum_{p=1}^{P} \|\phi\|_{r+2, \Omega_p}, $$

where $r = \min(j, \ell)$.

Next, Lemma 5.2 is used to estimate $\alpha_{\omega,H}$ given in (4.1).

**LEMMA 5.3.** *From the definition (4.1), it holds*

$$ (5.3) \qquad \alpha_{\omega,H} \lesssim H + \frac{C_s(\omega)}{\omega} (\omega H)^{\ell+1}. $$

*Proof.* Let $f \in L^2(\Omega)$, and let $\lambda^f \in \Lambda$ be the solution of $b(\mu, T\lambda^f) = -b(\mu, \hat{T}f)$ for all $\mu \in \Lambda$. We know that if $u \in H^1(\Omega)$ solves (2.1), then for each element $K \in \mathcal{T}_H$, we have $\mu|_{\partial K \setminus \Gamma_A} = -(1/\rho)\nabla u \cdot n_K$ and $\mu|_{\partial K \cap \Gamma_A} = 0$. Furthermore, according to Theorem 5.1, we know that $u$ admits decomposition (5.1). Similarly, $\lambda^f$ admits the decomposition

$$ \lambda^f = \sum_{j=0}^{\ell-1} \omega^j \lambda_j + \eta_\ell, $$

where $\lambda_j|_{\partial K \setminus \Gamma_A} = -(1/\rho)\nabla u_j \cdot n_K$, $\eta_\ell|_{\partial K \setminus \Gamma_A} = -(1/\rho)\nabla r_\ell \cdot n_K$, for all $K \in \mathcal{T}_H$ and $\lambda_j = \eta_\ell = 0$ on $\Gamma_A$.

Then, Lemma 5.2 and Theorem 5.1 state that

$$ \|\lambda_j - \pi_H \lambda_j\|_{\Lambda, \omega} \lesssim H^{j+1} \|u_j\|_{j+2} \lesssim H^{j+1} \|f\|_{0, \Omega}, \quad \|\eta_\ell - \pi_H \eta_\ell\|_{\Lambda, \omega} \lesssim H^{\ell+1} \|r_\ell\|_{l+2} \lesssim \frac{C_s(\omega)}{\omega} \omega^{\ell+1} H^{\ell+1} \|f\|_{0, \Omega}. $$

It follows that

$$ \begin{aligned} \|\lambda^f - \pi_H \lambda^f\|_{\Lambda, \omega} &\lesssim \sum_{j=0}^{\ell-1} \omega^j \|\lambda_j - \pi_H \lambda_j\|_{\Lambda, \omega} + \|\eta_\ell - \pi_H \eta_\ell\|_{\Lambda, \omega} \\ &\lesssim \sum_{j=0}^{\ell-1} \omega^j H^{j+1} \|f\|_{0, \Omega} + \frac{C_s(\omega)}{\omega} \omega^{\ell+1} H^{\ell+1} \|f\|_{0, \Omega} \\ &\lesssim \left( H \left( \sum_{j=0}^{\ell-1} (\omega H)^j \right) + \frac{C_s(\omega)}{\omega} \omega^{\ell+1} H^{\ell+1} \right) \|f\|_{0, \Omega}, \end{aligned} $$
and we conclude that

$$
(5.4) \quad \inf_{\mu_H \in \Lambda_H} \|\lambda^f - \mu_H\|_{\Lambda, \omega} \le \|\lambda^f - \pi_H \lambda^f\|_{\Lambda, \omega} \lesssim \left(H + \frac{C_s(\omega)}{\omega} (\omega H)^{\ell+1}\right) \|f\|_{0, \Omega},
$$

since $\omega H \lesssim 1$ by assumption (3.2). As a result, (5.3) follows directly from (5.4) and Definition (4.1) of $\alpha_{\omega,H}$.

The next result is a direct consequence of Theorems 4.3 and 5.3, and establishes the conver-
gence rates of the one-level MHM method with polynomial interpolation on faces.

THEOREM 5.4. Assume that $\omega H$ and $C_s(\omega)(\omega H)^{\ell+1}$ are small enough. Then, problems (2.12) and (2.13) admit unique solutions $\lambda \in \Lambda$ and $\lambda_H \in \Lambda_H$, and

$$
(5.5) \qquad \|\lambda - \lambda_H\|_{\Lambda, \omega} \lesssim \inf_{\mu_H \in \Lambda_H} \|\lambda - \mu_H\|_{\Lambda, \omega}.
$$

In addition, if $f \in H^j(\Omega)$ and $g \in H^{j+1/2}(\Gamma_A)$ for some $0 \le j \le m$, then we have

$$
(5.6) \quad \|\lambda - \lambda_H\|_{\Lambda, \omega} \lesssim \frac{C_s(\omega)}{\omega} (\omega H)^{r+1} (\|f\|_{r,\Omega} + \|g\|_{r+1/2,\Gamma_A}),
$$

where $r = \min(\ell, j)$. Furthermore, it holds that

$$
(5.7) \quad \|u - u_H\|_{V,\omega} \lesssim \frac{C_s(\omega)}{\omega} (\omega H)^{r+1} (\|f\|_{r,\Omega} + \|g\|_{r+1/2,\Gamma_A}),
$$

$$
(5.8) \quad \omega \|u - u_H\|_{0,\Omega} \lesssim \frac{C_s(\omega)}{\omega} \left( (\omega H)^{r+2} + \frac{C_s(\omega)}{\omega} (\omega H)^{r+\ell+2} \right) (\|f\|_{r,\Omega} + \|g\|_{r+1/2,\Gamma_A}),
$$

where $u = T\lambda + \hat{T}f + \tilde{T}g$, and $u_H = T\lambda_H + \hat{T}f + \tilde{T}g$.

Proof. As established in the proof of Theorem 2.6 of [10], we have<fnref n="2" /> $u \in H^{j+2}(\Omega_p)$ for all $p \in \{1, \dots, P\}$, and

$$
(5.9) \quad \|u\|_{j+2, \Omega_p} \lesssim \frac{C_s(\omega)}{\omega} \omega^{j+1} \sum_{p=1}^{P} \left( \|f\|_{j, \Omega_p}^2 + \|g\|_{j+1/2, \Gamma_A \cap \partial \Omega_p}^2 \right).
$$

Then, from (5.9) and Lemma 5.2, it holds that

$$
\|\lambda - \pi_H \lambda\|_{\Lambda, \omega} \lesssim \frac{C_s(\omega)}{\omega} (\omega H)^{r+1} (\|f\|_{r,\Omega} + \|g\|_{r+1/2,\Gamma_A}),
$$

and (5.6) and (5.7) directly follow from (5.5). Then, recalling (4.5) and (5.3), we have

$$
\begin{align*}
\omega \|u - u_H\|_{0,\Omega} &= \omega \|T(\lambda - \lambda_H)\|_{0,\Omega} \lesssim \omega \alpha_{\omega,H} \|T(\lambda - \lambda_H)\|_{V,\omega} \\
&\lesssim (\omega H + C_s(\omega) \omega^{\ell+1} H^{\ell+1}) \|T(\lambda - \lambda_H)\|_{V,\omega}
\end{align*}
$$

and (5.8) follows from (5.7).

**6. Numerical experiments.** In this section, we present numerical experiments illustrating our key results and demonstrating the performance of the proposed MHM method. We employ Cartesian meshes that are made of square elements. We use two different spaces $\Lambda_H$. On the one hand, we employ the polynomial spaces introduced in Section 5, and we will denote by $\Lambda_{H,\ell}$ the space spanned by piecewise polynomial of degree $\ell$. On the other hand, for each $\ell \ge 0$, we introduce the space $\tilde{\Lambda}_{H,\ell} := \Lambda_{H,0} \oplus S_{[\ell/2]}$ if $\ell$ is even and $\tilde{\Lambda}_{H,\ell} := \Lambda_{H,1} \oplus S_{[(\ell-1)/2]}$ if $\ell$ is odd, where

$$
S_k := \{ \lambda \in \Lambda_H \mid \lambda|_F \in \tilde{S}_k, \forall F \in \mathcal{F}_H \},
$$

<footnote n="2">² The proof is actually carried out for *g* = 0, but easily carries over the general case.</footnote>
with

$$ \tilde{S}_k := \operatorname{span}_{n \in \{1, \dots, k\}} \{e^{i\omega\nu_n t}, e^{-i\omega\nu_n t}\}, $$

and $\nu_n := \cos(\frac{n}{k+1}\frac{\pi}{2})$. If $\ell \le 1$, then $\tilde{\Lambda}_{H,\ell} = \Lambda_{H,\ell}$. For $\ell \ge 2$ however, oscillating basis functions are added instead of polynomial shape functions. Specifically, on a Cartesian mesh, one easily sees that the function

$$ (6.1) \qquad \mu|_{\partial K \setminus \Gamma_A} = \nabla e^{i\omega d \cdot x} \cdot n_K|_{\partial K \setminus \Gamma_A}, $$

belongs to the space $\tilde{\Lambda}_{H,\ell}$ if $d = (\cos \theta, \sin \theta)$ with $\theta = (n\pi)/(2k+2)$, with $n \in \{0, \dots, 4k+4\}$. As a result, we expect the resulting scheme to exactly reproduce plane wave propagating in a discrete set of directions $d$. Actually, we see that if $d = (1,0)$ or $d = (0,1)$, then $\mu$ in (6.1) restricted to $F \subset \partial K$ belongs to $\mathcal{P}_0(F)$ on a Cartesian mesh. Hence, we expect all MHM schemes to exactly reproduce plane wave travelling in the $x_1$ and $x_2$ directions.

**6.1. Stability analysis.** This experiment assesses the main theoretical results concerning the one-level MHM method, stated in Theorem 4.3 and Corollary 5.4. We consider the problem to find $u$ such that

$$ (6.2) \qquad \left\{ \begin{array}{ll} -\omega^2 u - \Delta u = 0 & \text{in } \Omega, \\ \nabla u \cdot n - i\omega u = g & \text{on } \partial\Omega, \end{array} \right. $$

where $g \in L^2(\partial\Omega)$ is selected so that

$$ u(x) = J(\omega|x - y|) + iY(\omega|x - y|), $$

with $y = (1.5, 0.5)$. This $u$ corresponds to the fundamental solution of the Helmholtz problem with Dirac mass at $y$. It is representative of applications, since a wave scattered by a small obstacle centered at $y$ behaves like $u$ far enough from the scatterer. Figure 6.1 depicts the solution for the frequency $\omega = 40\pi$.

<figure><img src="image_4.png" /><figcaption>FIG. 6.1. Solution $u$ to (6.2) for $\omega = 40\pi$. Real part (left) and imaginary part (right).</figcaption></figure>

We adopt the following strategy to verify the theoretical results presented in Theorem 4.3 and Corollary 5.4. For a fixed $\ell$, we solve (6.2) for different frequencies. Then, for each frequency, we fix a sequence of mesh sizes $H$, and we construct the convergence curves

$$ \mathcal{E}_{\text{MHM}}(H) = |u - T\lambda_H - \tilde{T}g|_{1,\tau_H}, \quad \mathcal{E}_{\text{INT}}(H) = |u - T(\Pi_H \lambda) - \tilde{T}g|_{1,\tau_H}. $$
<figure><img src="image_5.png" /><figcaption>FIG. 6.2. Convergence curves for $\ell = 0$ with $\omega = 20\pi$ (left) and $\omega = 40\pi$ (right).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 6.3. Convergence curves for $\ell = 1$ with $\omega = 30\pi$ (left) and $\omega = 150\pi$ (right).</figcaption></figure>

Here, $\lambda_H$ is the MHM solution and $\Pi_H \lambda$ an interpolant of $\lambda$. Such an interpolant is constructed edge-wise, and since in the experiment $u \in C^\infty(\bar{\Omega})$, then we have $\lambda|_F \in L^2(F)$ for all $F \in \mathcal{F}_H$. Thereby, $(\Pi_H \lambda)|_F$ is defined by $L^2(E)$-projection of $\lambda|_F$

We observe on Figures 6.2 and 6.3 the super-convergence predicted in (5.7), i.e.,

$$ \mathcal{E}_{\text{MHM}}(H) \lesssim H^{\ell+1}, \quad \mathcal{E}_{\text{INT}}(H) \lesssim H^{\ell+1}, $$

for $\ell = 0$ and $\ell = 1$. The pollution effect is visible on these two figures, as we can observe the gap between the interpolation error and the MHM error. Also, this gap is more important for higher frequencies, and less important when $\ell$ is increased.

Finally, we verify that the condition $\omega^{\ell+2}H^{\ell+1} \lesssim 1$ is necessary to ensure quasi-optimality. To this end, for each selected frequency $\omega$, we compute $H^*(\omega)$ as the largest value of $H$ (in the selected sequence of mesh sizes) such that

$$ (6.3) \qquad \mathcal{E}_{\text{MHM}}(H) \le 3 \mathcal{E}_{\text{INT}}(H) \quad \text{for all } H \le H^*(\omega). $$

The purpose of definition (6.3) is that if $H \le H^*(\omega)$, then we have

$$ (6.4) \qquad \mathcal{E}_{\text{MHM}}(H) \lesssim \mathcal{E}_{\text{INT}}(H), $$

with a constant independent of $\omega$ and $H$ (here, we have arbitrarily selected the constant 3). In Corollary 5.4, we showed that if $\omega^{\ell+2}H^{\ell+1} \lesssim 1$, then (6.4) holds. As illustrated by Figure 6.4, we have $H^*(\omega) \simeq \omega^{-1-1/(\ell+1)}$ for $\ell = 0$ and $\ell = 1$. Then, the necessary condition that $H \le H^*(\omega)$ is equivalent to $\omega^{\ell+2}H^{\ell+1} \lesssim 1$ which proves that Corollary 5.4 is sharp.
<figure><img src="image_7.png" /><figcaption>FIG. 6.4. $H^*(\omega)$ for $\ell = 0$ (left) and $\ell = 1$ (right).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 6.5. Anisotropy study for $\ell = 2$ (left) and $\ell = 4$ (right).</figcaption></figure>

**6.2. Anisotropy analysis.** We consider again test-case (6.2), but this time $g \in L^2(\Omega)$ is chosen so that

$$u(x) = e^{i\omega d \cdot x},$$

where $d = (\cos\theta, \sin\theta)$ is a unit vector representing a direction of propagation. We fix a mesh size $H$ and a frequency $\omega$ and solve (6.2) for 256 values of $\theta$ ranging from 0 to $\pi/2$ (as the mesh is made is squares, we obtain the remaining angles by symmetry). Figure 6.5 presents the results. For $\ell = 2$, we have chosen $\omega = 20\pi$ and $H = 1/11$, while the values $\omega = 40\pi$ and $H = 1/21$ have been selected for $\ell = 4$.

We employ both polynomial and oscillating shape functions. We see that as announced, the MHM solution is exact (up to the second-level accuracy) for some direction of propagation. When, using polynomial basis functions these “exact directions” are orthogonal to the mesh faces. On the other hand, oscillating basis functions permit to increase the number of “exact directions”. We also observe that overall in the experiment, oscillating basis functions provide a better accuracy than polynomials for the same number of degrees of freedom.

**6.3. Convergence study.** We consider the same test case as before with $\theta = \pi/13$ (so that the MHM method does not give the exact solution). Then, we fix the frequency $\omega = 10\pi$ and let $H \to 0$ in order to analyze the convergence rates. Figures 6.6 and 6.7 depict the results for $\ell = 2$ and $\ell = 3$.

We observe that the convergence rates are $O(H^{\ell+2})$ and $O(H^{\ell+1})$ in the $L^2$ and $H^1$ norms,
<figure><img src="image_9.png" /><figcaption>FIG. 6.6. Convergence study for $l = 2$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 6.7. Convergence study for $l = 3$</figcaption></figure>

respectively, as proved for polynomial basis functions. The convergence rates are the same for oscillating basis functions, though we did not prove this result. Also, we again observe that the oscillating basis functions provide more accurate results than polynomials.

**6.4. A multiscale test-case.** We consider the Marmousi II synthetical model [30]. The domain of propagation is 10,240 m large and 2,560 m deep. The P-wave velocity $c_p$ and the density $\rho$ are given on a 2,048 × 512 grid (see Figures 6.8 and 6.9). The bulk modulus is defined from the P-wave velocity and the density through the formula $\kappa = \rho c_p^2$.

A Dirichlet boundary condition is imposed on the top of the domain and a first-order absorbing boundary condition is applied on the rest of the boundary to simulate a semi-infinite propagation medium. The seismic source is represented by a Dirac right-hand side $\phi \in D'(\Omega)$ located at $x = 5000$ m and $z = 50$ m:

$$ \langle \phi, v \rangle = \overline{v(5000, 50)} \quad \text{for all } v \in D(\Omega). $$

There is no analytical solution for this benchmark. Hence, we use a finite element solution $u_{ref}$ computed on a very fine mesh as a reference. More precisely, we compute this solution with triangular Lagrangian elements of degree 4. The mesh is a 2048 × 512 cartesian grid, each square of the grid being subdivided into two triangles. It is worth noting that the mesh coincides with the same Cartesian grid as the media parameters. Hence, these parameters are constant on each cell, and we can use a standard finite element method to approximate the solution.

We evaluate the MHM solutions on a 513 × 129 cartesian grid which is used to compute relative
<figure><img src="image_11.png" /><figcaption>FIG. 6.8. Marmousi II: velocity model</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 6.9. Marmousi II: density model</figcaption></figure>

$L^2$ errors. The relative error of a solution $u_{mhm}$ computed with the MHM method is given by

$$E(u_{mhm}) = \left( \frac{\sum_{i=1}^{513} \sum_{j=1}^{129} |u_{ref}(x_i, z_j) - u_{mhm}(x_i, z_j)|^2}{\sum_{i=1}^{513} \sum_{j=1}^{129} |u_{ref}(x_i, z_j)|^2} \right)^{1/2},$$

where $x_i$ and $z_j$ correspond to the offset and depth of the evaluation grid lines.

We tabulate the error of the MHM solutions for different choices of $H$ and $\ell$. We use 3 different values of $H$: 20, 40 and 80 m, which corresponds to Cartesian grids of size 512 × 128, 256 × 64 and 128 × 32 respectively. Thus, it is clear that the original medium parameters grid is coarsened by a factor 4, 8 or 16.

For the second level methodology, we use square Lagrange finite elements of degree 3 based on a 8 × 8, 16 × 16 or 32 × 32 cartesian grid for $h = 20, 40$ or $80$ m, respectively. This choice ensures that the second-level computational scale is twice smaller than the medium parameters grid. In particular, it is clear that the second-level mesh size matches the scale of the heterogeneities.

In order to give a comparison with more standard finite element methods, we also compute the solutions with the MMAM [5, 16, 17]. The method is based on triangular Lagrangian elements so that we subdivide each square of the MHM grids into two triangles to apply the MMAM and we use enough subcells to take into account the medium parameters exactly.

We solve the problem for the frequency $f = 20$ Hz, the angular frequency being defined as $\omega = 2\pi f$. The results are in Table 6.1.

<table><thead><tr><th colspan="6">MHM method</th></tr><tr><th><i>h</i></th><th><i>l</i> = 0</th><th><i>l</i> = 1</th><th><i>l</i> = 2</th><th><i>l</i> = 3</th><th><i>l</i> = 4</th></tr></thead><tbody><tr><td>20</td><td>138</td><td>2.90</td><td>0.13</td><td>0.11</td><td>0.03</td></tr><tr><td>40</td><td>225</td><td>45.1</td><td>1.47</td><td>0.37</td><td>0.20</td></tr><tr><td>80</td><td>161</td><td>154</td><td>44.8</td><td>4.12</td><td>0.33</td></tr></tbody></table>

<table><caption>TABLE 6.1</caption><thead><tr><th colspan="6">MMAm</th></tr><tr><th><i>h</i></th><th><i>p</i> = 1</th><th><i>p</i> = 2</th><th><i>p</i> = 3</th><th><i>p</i> = 4</th><th><i>p</i> = 5</th></tr></thead><tbody><tr><td>20</td><td>127</td><td>90.2</td><td>2.41</td><td>0.85</td><td>0.28</td></tr><tr><td>40</td><td>121</td><td>134</td><td>80.3</td><td>5.23</td><td>1.34</td></tr><tr><td>80</td><td>101</td><td>124</td><td>132</td><td>124</td><td>57.4</td></tr></tbody></table>

Relative error (%) in the Marmousi II model.

From Table 6.1, it is clear that the MHM method can produce accurate solutions on coarse meshes. Furthermore, when considering the same mesh and the same order of discretization, the
MHM method outperforms the MMAM regarding accuracy. This fact is interesting because the global linear system has approximately the same size and filling in both cases. The drawback for the MHM method is that it is required to solve local problems to construct the global linear system. On the other hand, such expensive computations are independent to one another, and therefore, they can take advantage of parallel facilities.

**7. Conclusion.** We presented a new approach to solve the heterogeneous Helmholtz equation: the Multiscale Hybrid Mixed (MHM) method. The MHM method has many similarities with the DEM of Farhat et al. [15] since both methods rely on the primal hybrid formulation of the Helmholtz equation. The difference lies in the fact that the MHM base functions are local solutions of the Helmholtz equation which need to be computed with a second-level strategy, while the DEM uses plane waves.

The MHM and DEM are similar in homogeneous propagation media. Notably, the lowest-order MHM method recovers DEM elements in some cases. However, the MHM method behaves differently when the propagation medium is heterogeneous. In this context, the MHM method can leverage sub-element variations of the media, whereas coefficient parameters need to be constant in each mesh element for the DEM. Our numerical experiments illustrate the robustness of the MHM regarding small-scale heterogeneities.

We proposed a convergence analysis of the MHM method for elements of arbitrary orders with curved boundaries. We proved that the method is super-convergent under usual constraints on the mesh. Also, the error analysis generalizes the convergence results for the lowest-order DEM elements presented in [2].

Numerical experiments illustrated the accuracy of the MHM method on analytical solutions and geophysical benchmarks. The analytical experiments validated our convergence analysis. On the other, our experiments on the Marmousi II geophysical benchmark [30] showed that the method is very efficient for geophysical applications. In particular, we highlighted the superiority of the MHM method over polynomial Lagrangian elements.

**Acknowledgements.** The second author was supported by INRIA/France, CNPq/Brazil under the Project No. 301576/2013-0, EU H2020 Program and from MCTI/RNP-Brazil under the HPC4E Project, Grant Agreement No. 689772, and by CAPES/Brazil No. 88881.143295/2017-01 under the PHOTOM Project.

**Appendix A. Stability of local problems.** In this section, we establish inf-sup conditions for the local problems. Specifically, given $K \in \mathcal{T}_H$, we show that for all $u \in H^1(K)$, there exists an element $u^* \in H^1(K)$ such that

$$ (A.1) \qquad \mathrm{Re} a(u, u^*) \gtrsim \|u\|_{V(K), \omega}^2 \quad \text{and} \quad \|u^*\|_{V(K), \omega} \lesssim \|u\|_{V(K), \omega}. $$

We give four different proofs, depending on whether $K$ is convex, $\partial K \cap \Gamma_A = \emptyset$, and if the coefficient $\kappa$ and $\rho$ are constant functions in $K$.

**PROPOSITION A.1.** Let $K \in \mathcal{T}_H$ be a convex element such that $\partial K \cap \Gamma_A = \emptyset$. If (3.2) is satisfied, then (A.1) holds.

*Proof.* Let $u \in H^1(K)$, and denote

$$ u_0 := \frac{1}{|K|} \int_K u \quad \text{and} \quad u^{\perp} := u - u_0, $$

as well as $\kappa_* := \inf_K \kappa$, $\kappa^* := \sup_K \kappa$ and $\rho_* := \inf_K \rho$, $\rho^* = \sup_K \rho$.

We first establish that

$$ (A.2) \qquad \mathrm{Re} a(u, v^*) \geq \frac{\delta}{\rho^*} |u|_{1,K}^2, $$

where $v^* := u - 2u_0$. We observe that as $\partial K \cap \Gamma_A = \emptyset$, we have

$$ \mathrm{Re} a(u, v^*) = \mathrm{Re} a(u^{\perp} + u_0, u^{\perp} - u_0) = \mathrm{Re} a(u^{\perp}, u^{\perp}) - \mathrm{Re} a(u_0, u_0). $$
Using that $\nabla u_0 = 0$, we easily see that

$$-a(u_0, u_0) = \omega^2 \int_K \frac{1}{\kappa} |u_0|^2 \geq 0,$$

so that

$$(A.3) \qquad \mathrm{Re} a(u, v^*) \geq \mathrm{Re} a(u^{\perp}, u^{\perp}).$$

Next, using that $u^{\perp}$ has zero mean value and $K$ is convex with $\mathrm{diam} K \leq H$, we apply the Poincaré-Wirtinger inequality (c.f. [36]) to get

$$\|u^{\perp}\|_{0,K} \leq \frac{H}{\pi} |u^{\perp}|_{1,K},$$

from which we deduce that

$$a(u^{\perp}, u^{\perp}) \geq -\frac{\omega^2}{\kappa_{\star}} \|u^{\perp}\|_{0,K}^2 + \frac{1}{\rho^{\star}} |u^{\perp}|_{1,K}^2 \geq \left( \frac{1}{\rho^{\star}} - \frac{\omega^2 H^2}{\pi^2 \kappa_{\star}} \right) |u^{\perp}|_{1,K}^2 \geq \frac{1}{\rho^{\star}} \left( 1 - \frac{\rho^{\star}}{\kappa_{\star}} \frac{\omega^2 H^2}{\pi^2} \right) |u^{\perp}|_{1,K}^2.$$

We recall that by assumption (3.2), we have

$$\frac{\rho^{\star}}{\kappa_{\star}} \frac{\omega^2 H^2}{\pi^2} \leq 1 - \delta,$$

and therefore

$$(A.4) \qquad \mathrm{Re} a(u^{\perp}, u^{\perp}) \geq \frac{\delta}{\rho^{\star}} |u^{\perp}|_{1,K}^2 = \frac{\delta}{\rho^{\star}} |u|_{1,K}^2.$$

where we used $\nabla u_0 = 0$. At that point, (A.2) follows from (A.3) and (A.4).

Now, that (A.2) is established, we observe that

$$\mathrm{Re} a(u, -u) \geq \frac{\omega^2}{\kappa^*} \|u\|_{0,K}^2 - \frac{1}{\rho^*} |u|_{1,K}^2,$$

and, as a result, it holds

$$\mathrm{Re} a(u, 2(\rho^*/\rho_\star)\delta^{-1}v^* - u) \geq \frac{\omega^2}{\kappa^*} \|u\|_{0,K}^2 + \frac{1}{\rho_\star} |u|_{1,K}^2 \gtrsim \|u\|_{V(K),\omega}^2,$$

and (A.1) follows taking $u^* := 2(\rho^*/\rho_\star)\delta^{-1}v^* - u$. $\square$

**PROPOSITION A.2.** Let $K \in \mathcal{T}_H$ be such that $\partial K \cap \Gamma_A = \emptyset$, and assume that there exists a bijective mapping $\tilde{M} : \tilde{K} \to K$, where $\tilde{K}$ is the straight simplex having the same vertices of $K$. Moreover, assume that

$$(A.5) \qquad \frac{J_+}{J_-} B_- \geq 1 - \frac{\delta}{2} \quad \text{and} \quad B_+ \leq 1 + \delta,$$

where

$$J_+ := \sup_{\tilde{K}} J, \quad J_- := \inf_{\tilde{K}} J, \quad B_+ := \sup_{\tilde{K}} \beta, \quad B_- := \inf_{\tilde{K}} \alpha,$$

and $\mathbf{J}$ is the jacobian matrix of $\tilde{M}$, $\mathbf{B} := \mathbf{J}\mathbf{J}^T$, $J := |\det \mathbf{J}|$, and $\alpha(\tilde{\mathbf{x}})$ and $\beta(\tilde{\mathbf{x}})$ are the smallest and highest eigenvalue of $\mathbf{B}(\tilde{\mathbf{x}})$, for each $\tilde{\mathbf{x}} \in \tilde{K}$. Then, (A.1) holds.

*Proof.* Given $v \in H^1(K)$, we denote $\tilde{v} := v \circ \tilde{M}$. Observe that

$$J_- \|\tilde{v}\|_{0, \tilde{K}}^2 \leq \|v\|_{0, K}^2 \leq J_+ \|\tilde{v}\|_{0, \tilde{K}}^2,$$
and

$$J_{-}B_{-}|\tilde{v}|_{1,\tilde{K}}^{2} \leq |v|_{1,K}^{2} \leq J_{+}B_{+}|\tilde{v}|_{1,\tilde{K}}^{2}.$$

Next, define $u_0 \in H^1(K)$ as follows

$$u_0 = \frac{1}{|\tilde{K}|} \int_{\tilde{K}} \tilde{u} \quad \text{for } u \in H^1(K),$$

and remark that $\tilde{u}_0 := u_0$ since $u_0$ is a constant function. Also, set $u^{\perp} = u - u_0$ and $v^{\star} = u - 2u_0$, and then, arguing as in the convex case, it holds

$$\mathrm{Re} a(u, v^{\star}) \geq a(u^{\perp}, u^{\perp}).$$

We further remark that

$$\int_{\tilde{K}} \tilde{u}^{\perp} = \int_{\tilde{K}} (\tilde{u} - u_0) = 0.$$

Then, we have

$$\begin{aligned} a(u^{\perp}, u^{\perp}) &\geq -\frac{\omega^2}{\kappa_{\star}} \|u^{\perp}\|_{0,K}^2 + \frac{1}{\rho^{\star}} |u^{\perp}|_{1,K}^2 \\ &\geq -\frac{\omega^2}{\kappa_{\star}} J_{+} \|\tilde{u}^{\perp}\|_{0,\tilde{K}}^2 + \frac{1}{\rho^{\star}} J_{-} B_{-} |\tilde{u}^{\perp}|_{1,\tilde{K}}^2 \\ &\geq \left( \frac{1}{\rho^{\star}} J_{-} B_{-} - \frac{\omega^2 H^2}{\pi^2 \kappa_{\star}} J_{+} \right) |\tilde{u}^{\perp}|_{1,\tilde{K}}^2 \\ &\geq \frac{1}{\rho^{\star}} J_{+} \left( \frac{J_{+}}{J_{-}} B_{-} - \frac{\rho^{\star} \omega^2 H^2}{\kappa_{\star} \pi^2} \right) |\tilde{u}^{\perp}|_{1,\tilde{K}}^2 \end{aligned}$$

since $\tilde{K}$ is convex and $\mathrm{diam} \tilde{K} \leq H$. Next, we recall that

$$\frac{\rho^{\star} \omega^2 H^2}{\kappa_{\star} \pi^2} \leq 1 - \delta$$

by assumption (3.2), so that

$$a(u, v^*) \geq \frac{1}{\rho^*} J_+ \left( \frac{J_+}{J_-} B_- - 1 + \delta \right) |\tilde{u}^{\perp}|_{1, \tilde{K}}^2 \geq \frac{1}{\rho^*} J_+ \frac{\delta}{2} |\tilde{u}^{\perp}|_{1, \tilde{K}}^2$$

since $|(J_+/J_-)B_- - 1| \geq -\delta/2$ by assumption (A.5). Then, we have

$$a(u, v^*) \geq \frac{1}{\rho^*} \frac{1}{B_+} \frac{\delta}{2} |u^{\perp}|_{1,K}^2 \geq \frac{1}{\rho^*} \frac{\delta}{2(1+\delta)} |u^{\perp}|_{1,K}^2,$$

and we conclude following the convex case proof. □

**PROPOSITION A.3.** Let $K \in \mathcal{T}_H$ such that $\partial K \cap \Gamma_A \neq \emptyset$. If (3.4) is satisfied, then (A.1) holds.

*Proof.* As before, we denote $\kappa_{\star} := \inf_K \kappa$, $\kappa^{\star} := \sup_K \kappa$, $\rho_{\star} := \inf_K \rho$ and $\rho^{\star} = \sup_K \rho$. We also set $c_{\star} := \sqrt{\kappa_{\star}/\rho^{\star}}$.

We denote by $F \subset \Gamma_A$ the (possibly) curved face of $K$ that belongs to the boundary of $\Omega$. We first establish that

$$(A.6) \qquad \omega^2 \|v\|_{0,K}^2 \leq \frac{1}{d-1} \left( \omega H \left( \omega \|v\|_{0,E}^2 \right) + \omega^2 H^2 |v|_{1,K}^2 \right) \quad \text{for all } v \in H^1(K).$$
We denote by $b \in \mathbb{R}^d$ the vertex of $K$ opposite to $F$ and define $\sigma(x) := x - b$. We note that $|\sigma|_{0,\infty,K} \le H$ and $\nabla \cdot \sigma = d$. In addition, we have $\sigma \cdot n_K = 0$ on $\partial K \setminus E$. Thus, for every $v \in H^1(K)$, we have

$$ \int_{\partial K} |v|^2 \sigma \cdot n_K = \int_E |v|^2 \sigma \cdot n_K \le H \|v\|_{0,E}^2. $$

On the other hand, from the Stokes' formula we get

$$ \int_{\partial K} |v|^2 \sigma \cdot n_K = \int_K \nabla \cdot (|v|^2 \sigma) = \int_K \nabla \cdot \sigma |v|^2 + \int_K \sigma \cdot \nabla |v|^2 = d \|v\|_{0,K}^2 + \int_K \sigma \cdot \nabla |v|^2, $$

and we obtain

$$ (A.7) \qquad d\|v\|_{0,K}^2 \le H\|v\|_{0,E}^2 - \int_K \sigma \cdot \nabla |v|^2. $$

We focus on the last term in (A.7). Recalling that $\nabla|v|^2 = 2\operatorname{Re}\bar{v}\nabla v$, we have

$$ \int_K \sigma \cdot \nabla |v|^2 = 2 \operatorname{Re} \int_K \bar{v} \sigma \cdot \nabla v. $$

Hence, using Young's inequality and recalling that $|\sigma|_{0,\infty,K} \le H$, we have

$$ (A.8) \qquad \left| \int_K \sigma \cdot \nabla |v|^2 \right| = 2H \operatorname{Re} \int_K |v| |\nabla v| \le \|v\|_{0,K}^2 + H^2 |v|_{1,K}^2. $$

At that point, we plug (A.8) into (A.7) to obtain

$$ (d-1) \|v\|_{0,K}^2 \le H \|v\|_{0,E}^2 + H^2 |v|_{1,K}^2, $$

and (A.6) follows by multiplying by $\omega^2/(d-1)$.

Now, let $u \in H^1(K)$. We have

$$ \operatorname{Re} a(u, iu) = \omega \int_{\partial K \cap \Gamma_A} \frac{1}{\sqrt{\kappa \rho}} |u|^2 = \omega \int_E \frac{1}{\sqrt{\kappa \rho}} |u|^2 \ge \frac{\omega}{\sqrt{\kappa^* \rho^*}} \|u\|_{0,E}^2 $$

and

$$ \operatorname{Re} a(u, u) \ge -\frac{\omega^2}{\kappa_*} \|u\|_{0,K}^2 + \frac{1}{\rho^*} |u|_{1,K}^2. $$

Then, using (A.6), it follows that

$$ \begin{align*} \operatorname{Re} a(u, (1+i\beta)u) &\ge -\frac{\omega^2}{\kappa_*} \|u\|_{0,K}^2 + \frac{1}{\rho^*} |u|_{1,K}^2 + \frac{\beta\omega}{\sqrt{\kappa^*\rho^*}} \|u\|_{0,E}^2 \\ &\ge -\omega H \frac{\omega}{(d-1)\kappa_*} \|u\|_{0,E}^2 - \frac{\omega^2 H^2}{(d-1)\kappa_*} |u|_{1,K}^2 + \frac{1}{\rho^*} |u|_{1,K}^2 + \frac{\beta\omega}{\sqrt{\kappa^*\rho^*}} \|u\|_{0,E}^2 \\ &\ge \omega \left( \frac{\beta}{\kappa^*\rho^*} - \frac{\omega H}{(d-1)\kappa_*} \right) \|u\|_{0,E}^2 + \left( \frac{1}{\rho^*} - \frac{\omega^2 H^2}{(d-1)\kappa_*} \right) |u|_{1,K}^2 \\ &\ge \frac{\omega}{\rho^*} \left( \frac{\beta}{\kappa^*} - \frac{1}{d-1} \frac{\omega H}{c_*^2} \right) \|u\|_{0,E}^2 + \frac{1}{(d-1)\rho^*} \left( (d-1) - \frac{\omega^2 H^2}{c_*^2} \right) |u|_{1,K}^2 \\ &\ge \frac{\omega}{(d-1)\rho^* c_*} \left( \beta \frac{c_*}{\kappa^*} (d-1) - \frac{\omega H}{c_*} \right) \|u\|_{0,E}^2 + \frac{1}{(d-1)\rho^*} \left( (d-1) - \frac{\omega^2 H^2}{c_*^2} \right) |u|_{1,K}^2. \end{align*} $$

We select $\beta = \kappa^*/c_*(d-1)^{-1/2}$ and from (3.2), it holds

$$ \begin{align*} \operatorname{Re} a(u, (1+i\beta)u) &\ge \frac{\omega}{\rho^* c_*} \left( \sqrt{d-1} - \frac{\omega H}{c_*} \right) \|u\|_{0,E}^2 + \left( (d-1) - \frac{\omega^2 H^2}{c_*^2} \right) |u|_{1,K}^2 \\ &\ge \frac{\omega}{(d-1)\rho^* c_*} (1-\delta)^{1/2} \|u\|_{0,E}^2 + \frac{1}{(d-1)\rho^*} (1-\delta) |u|_{1,K}^2. \end{align*} $$
Then, recalling that

$$ \mathrm{Re} a(u, -u) \geq \frac{\omega^2}{\kappa^*} \|u\|_{0,K}^2 - \frac{1}{\rho_*} |u|_{1,K}^2, $$

we obtain

$$ \mathrm{Re} a(u, 2(d-1))\rho^*/\rho_*(1-\delta)^{-1}(1+i\beta)u - u \gtrsim \|u\|_{V(K),\omega}^2 $$

and since we have

$$ \|2(d-1)(\rho^*/\rho_*)(1-\delta)^{-1}(1+i\beta)u - u\|_{V(K),\omega} \lesssim \|u\|_{V(K),\omega} $$

the result follows.

PROPOSITION A.4. Let $K \in \mathcal{T}_H$ such that $\partial K \cap \Gamma_A \neq \emptyset$ and that $\kappa|_K = \kappa_K \in \mathbb{R}$, $\rho|_K = \rho_K \in \mathbb{R}$. Then, (A.1) holds.

*Proof.* We conserve the notations of the previous propositions in this proof. We have

$$ \mathrm{Re} a(u, (1+i)u) = \frac{\omega^2}{\kappa_K} \|u\|_{0,K}^2 + \frac{\omega}{\sqrt{\kappa_K \rho_K}} \|u\|_{0,F}^2 + \frac{1}{\rho_K} |u|_{1,K}^2 - \frac{2\omega^2}{\kappa_K} \|u\|_{0,K}^2. $$

Thus, defining $z$ as the unique element of $H^1(K)$ such that $a(w, z) = (u, w)$ for all $w \in H^1(K)$, and as a result it holds $a(u, z) = \|u\|_{0,K}^2$, and

$$ \mathrm{Re} a \left( u, (1+i)u + \frac{2\omega^2}{\kappa_K} z \right) = \frac{\omega^2}{\kappa_K} \|u\|_{0,K}^2 + \frac{\omega}{\sqrt{\kappa_K \rho_K}} \|u\|_{0,F}^2 + \frac{1}{\rho_K} |u|_{1,K}^2 \gtrsim \|u\|_{V(K),\omega}^2. $$

Hence, it remains to show that $\|z\|_{V(K),\omega} \lesssim \|u\|_{V(K),\omega}$. To do so, we proceed as in [25] and pick the test function $w = \sigma \cdot \nabla z$ in the definition of $z$. Standard computations show that

$$ 2 \mathrm{Re} a(\sigma \cdot \nabla z, z) = d \frac{\omega^2}{\kappa_K} \|z\|_{0,K}^2 + \frac{1}{\rho_K} \int_{\partial K} |\nabla z| \sigma \cdot n - \frac{d-2}{\rho_K} |z|_{1,K}^2 - \frac{\omega^2}{\kappa_K} \int_{\partial K} |z|^2 \sigma \cdot n - 2 \mathrm{Re} \frac{i\omega}{\sqrt{\kappa_K \rho_K}} \int_{\Gamma_A} z \sigma \cdot \nabla z, $$

so that

$$ d \frac{\omega^2}{\kappa_K} \|z\|_{0,K}^2 + \frac{1}{\rho_K} \int_{\partial K} |\nabla z| \sigma \cdot n = 2 \mathrm{Re}(u, \sigma \cdot \nabla z) + \frac{d-2}{\rho_K} |z|_{1,K}^2 + \frac{\omega^2}{\kappa_K} \int_{\partial K} |z|^2 \sigma \cdot n + 2 \mathrm{Re} \frac{i\omega}{\sqrt{\kappa_K \rho_K}} \int_{\Gamma_A} z \sigma \cdot \nabla z, $$

and

$$ d \frac{\omega^2}{\kappa_K} \|z\|_{0,K}^2 + \frac{\gamma_K H_K}{\rho_K} \|\nabla z\|_{0,\Gamma_A}^2 \leq 2H_K \|u\|_{0,K} \|z\|_{1,K} + \frac{d-2}{\rho_K} \|z\|_{1,K}^2 + \frac{\omega^2 H_K}{\kappa_K} \|z\|_{0,\Gamma_A}^2 + 2 \frac{\omega H_K}{\sqrt{\kappa_K \rho_K}} \|z\|_{0,\Gamma_A} \|\nabla z\|_{0,\Gamma_A}. $$

Then, since

$$ \begin{aligned} 2 \frac{\omega H_K}{\sqrt{\kappa_K \rho_K}} \|z\|_{0,\Gamma_A} \|\nabla z\|_{0,\Gamma_A} &\leq \frac{\omega^2 H_K}{\sqrt{\kappa_K \rho_K}} \frac{\rho_K}{\gamma_K} \|z\|_{0,\Gamma_A}^2 + \frac{\gamma_K H_K}{\rho_K} \|\nabla z\|_{0,\Gamma_A}^2 \\ &\leq \frac{\omega^2 H_K}{\gamma_K c_K} \|z\|_{0,\Gamma_A}^2 + \frac{\gamma_K H_K}{\rho_K} \|\nabla z\|_{0,\Gamma_A}^2, \end{aligned} $$
for a positive constant $\gamma_K$, it holds

$$ 
\begin{aligned}
d_{\kappa K} \frac{\omega^2}{\kappa_K} \|z\|_{0,K}^2 &\le 2H_K \|u\|_{0,K} |z|_{1,K} + \frac{d-2}{\rho_K} |z|_{1,K}^2 + \left( \frac{\omega^2 H_K}{\kappa_K} + \frac{\omega^2 H_K}{\gamma_K c_K} \right) \|z\|_{0,\Gamma_A}^2 \\
&\le 2H_K \|u\|_{0,K} |z|_{1,K} + \frac{d-2}{\rho_K} |z|_{1,K}^2 + \omega H_K \left( \frac{1}{\kappa_K} + \frac{1}{\gamma_K c_K} \right) \omega \|z\|_{0,\Gamma_A}^2 \\
(A.9) \qquad &\le \rho_K H_K^2 \|u\|_{0,K}^2 + \frac{d-1}{\rho_K} |z|_{1,K}^2 + \omega H_K \left( \frac{1}{\kappa_K} + \frac{1}{\gamma_K c_K} \right) \omega \|z\|_{0,\Gamma_A}^2,
\end{aligned}
$$

where $\rho_K$ is a positive constant. Then, we have

$$ (A.10) \qquad \frac{1}{\rho_K} |z|_{1,K}^2 = \operatorname{Re} a(z, z) + \frac{\omega^2}{\kappa_K} \|z\|_{0,K}^2 \le \|u\|_{0,K} \|z\|_{0,K} + \frac{\omega^2}{\kappa_K} \|z\|_{0,K}^2 $$

so that

$$ (A.11) \qquad \frac{d-1}{\rho_K} |z|_{1,K}^2 \le (d-1) \|u\|_{0,K} \|z\|_{0,K} + \frac{\omega^2(d-1)}{\kappa_K} \|z\|_{0,K}^2 $$

and using (A.9) and (A.11), we get

$$ (A.12) \qquad \frac{\omega^2}{\kappa_K} \|z\|_{0,K}^2 \le \rho_K H_K^2 \|u\|_{0,K}^2 + (d-1) \|u\|_{0,K} \|z\|_{0,K} + \omega H_K \left( \frac{1}{\kappa_K} + \frac{1}{\gamma_K c_K} \right) \omega \|z\|_{0,\Gamma_A}^2 $$

We also have

$$ (A.13) \qquad \frac{\omega}{\sqrt{\kappa_K \rho_K}} \|z\|_{0,\Gamma_A}^2 = - \operatorname{Im} a(z, z) = - \operatorname{Im}(u, z)_K \le \|u\|_{0,K} \|z\|_{0,K} $$

and therefore, since $\omega H_K \lesssim 1$, from (A.12) and (A.13) it holds

$$ 
\begin{aligned}
\frac{\omega^2}{\kappa_K} \|z\|_{0,K}^2 &\le \rho_K H_K^2 \|u\|_{0,K}^2 + (d-1) \|u\|_{0,K} \|z\|_{0,K} + \omega H_K \left( \frac{1}{\kappa_K} + \frac{1}{\gamma_K c_K} \right) \sqrt{\kappa_K \rho_K} \|u\|_{0,K} \|z\|_{0,K} \\
&\le \rho_K H_K^2 \|u\|_{0,K}^2 + A \|u\|_{0,K} \|z\|_{0,K} \\
&\le \left( \rho_K H_K^2 + \frac{1}{2} \frac{\kappa_K A^2}{\omega^2} \right) \|u\|_{0,K}^2 + \frac{1}{2} \omega^2 \kappa_K \|z\|_{0,K}^2,
\end{aligned}
$$

where $A := \max\{(d-1), \omega H_K \left( \frac{1}{\kappa_K} + \frac{1}{\gamma_K c_K} \right) \sqrt{\kappa_K \rho_K} \}$, and then

$$ \frac{\omega^2}{\kappa_K} \|z\|_{0,K}^2 \le 2 \left( \rho_K H_K^2 + \frac{1}{2} \frac{\kappa_K A^2}{\omega^2} \right) \|u\|_{0,K}^2 \lesssim (H_K^2 + \omega^{-2}) \|u\|_{0,K}^2 \lesssim \omega^{-2} (1 + \omega^2 H_K^2) \|u\|_{0,K}^2 \lesssim \omega^{-2} \|u\|_{0,K}^2, $$

and therefore $\omega \|z\|_{0,K} \lesssim \omega^{-1} \|u\|_{0,K}$. Recalling (A.10) and (A.13), we arrive that

$$ \omega \|z\|_{0,\Gamma_A}^2 \lesssim \|u\|_{0,K} \|z\|_{0,K} \lesssim \omega^{-2} \|u\|_{0,K}^2, \quad |z|_{1,K}^2 \lesssim \|u\|_{0,K} \|z\|_{0,K} + \omega^2 \|z\|_{0,K} \lesssim \omega^{-2} \|u\|_{0,K}^2, $$

which leads to

$$ \|z\|_{V(K),\omega} \lesssim \omega^{-1} \|u\|_{0,K} \lesssim \omega^{-2} \|u\|_{V(K),\omega}, $$

and the result follows.

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

REMARK 7. It is possible to extend the above results to allow some variations of $\kappa$ and $\rho$ inside $K$, as shown for instance in Theorem B.2 of [11]. Such requirements on $\kappa$ and $\rho$ would correspond to non-trapping geometries, as discussed in [18]. However, we focus on the case where they are constant functions to avoid unnecessary technicalities.
<figure><img src="image_14.png" /><figcaption>FIG. B.1. Local change of coordinates.</figcaption></figure>

**Appendix B. Curved elements.** We gather here results concerning interpolation on curved elements. We refer the reader to Figure B.1 for an illustration of the adopted notations.

**PROPOSITION B.1.** There exists a constant $H_{\mathcal{P}}$ that only depends on the partition $\mathcal{P}$ of $\Omega$ such that, if $H \le H_{\mathcal{P}}$ and if $K \in \mathcal{T}_H$ is such that $F = \partial K \cap \partial \Omega_p \ne \emptyset$ for some $\Omega_p \in \mathcal{P}$, then

• There exists a vector field $\tilde{n} \in C^m(K, \mathbb{R}^2)$ with $\|\tilde{n}\|_{m,\infty,K} \lesssim 1$, such that $\tilde{n} = n$ on $F$.

• There exists a $C^1$ diffeomorphism $\phi_F : (0,1) \to F$ such that $\sup_{s \in (0,1)} |\phi'_F(s)| \lesssim H$.

*Proof.* Since $\partial \Omega_p$ is of class $C^{m,1}$, there exists a finite number of open sets $\{U_j\}_{j=1}^{N_p}$ such that $\partial \Omega_p \subset \cup_{j=1}^{N_p} U_j$ and for $j = 1, \dots, N$ there exists a $C^{m,1}$ diffeomorphism $\psi_j : (-1,1)^2 \to U_j$ such that

$$U_j \cap \Omega = \psi_j((-1,1) \times (-1,0)), \quad U_j \setminus \Omega = \psi_j((-1,1) \times (1,1)), \quad U_j \cap \partial \Omega = \psi_j((-1,1) \times \{0\}).$$

In addition, we have

$$(B.1) \qquad \mathbf{n}(\mathbf{x}) = \frac{1}{\sqrt{1 + |\partial_1 \psi_j(\mathbf{y})|}} \begin{pmatrix} -\partial_1 \psi_j(\mathbf{y}) \\ 1 \end{pmatrix}$$

for all $\mathbf{x} \in \partial \Omega \cap U_j$, where $\mathbf{y} = \psi_j^{-1}(\mathbf{x})$.

We define

$$H_{\Omega_p} := \min_{j \in \{1, \dots, N_p\}} \mathrm{diam} U_j, \quad H_{\mathcal{P}} := \min_{p \in \{1, \dots, P\}} H_{\Omega_p}$$

As a result, under the assumption that $H \le H_{\mathcal{P}}$ for each $K \in \mathcal{T}_H$ such that $K \cap \partial \Omega_p \ne \emptyset$ for some $p \in \{1, \dots, P\}$, we can select $j \in \{1, \dots, N_p\}$ such that $K \subset U_j$. To simplify, we fix one element $K$ and we ommit the index $j$ in the remaining of the proof.

On the one hand, we see from (B.1) that we can define $\tilde{n}$ as

$$\tilde{n}(\mathbf{x}) := \frac{1}{\sqrt{1 + |\partial_1 \psi(\psi^{-1}(\mathbf{x}))|}} \begin{pmatrix} -\partial_1 \psi(\psi^{-1}(\mathbf{x})) \\ 1 \end{pmatrix}$$

for all $\mathbf{x} \in K \subset \partial U_j$. Since $\psi_j$ is a $C^m$ diffeomorphism whose norm depends only on $\Omega$, we see that $\|\tilde{n}\|_{m,\infty,K} \lesssim 1$ where the hidden constant only depends on $\Omega$.

On the other hand, we see that $\psi_{\partial \Omega}(t) := \psi(t,0)$ is a $C^m$ diffeomorphism from $(-1,1)$ to $\partial \Omega \cap U$, and we have $\sup_{t \in (0,1)} |\psi'_{\partial \Omega}(t)| \gtrsim 1$ where the hidden constant only depends on $U$. Then, the measure of $F$ is at most $H$, and a simple scaling argument shows that $F = \psi_{\partial \Omega}((t',t''))$ with $t'' - t' \lesssim H/H_{\Omega}$. At that point, we define

$$\phi_E(s) := \psi_{\partial \Omega}(t' + (t'' - t')s),$$

where $s \in (0,1)$, which satisfies the statement of the proposition. □
Now, let us assume that $H \le H_\Omega$, so that Proposition B.1 applies. Also, in the following, $\Lambda_H$ is defined according to (5.2) with $\ell \le m$, where $m$ is the regularity of the vector field $\tilde{n}$ in Proposition B.1.

*Proof of Lemma 5.2.* We note that for all $v \in V$, it holds

$$ \left| \sum_{K \in \mathcal{T}_H} \int_{\partial K} (\mu - \pi_h \mu) v \right| \le \sum_{K \in \mathcal{T}_H} \sum_{E \subset \partial K \setminus \Gamma_A} \int_E |(\mu - \pi_h \mu) v|, $$

where the integrals are usual $L^2$ integrals, since $\mu \in L^2(E)$ for all $F \in \mathcal{F}_H$. In the following, we use the notation $\chi := (1/\rho)\nabla\phi \cdot \tilde{n}_F$, where $\tilde{n}_F$ is an extension of $n_F$ to $K$. We note if the face $F$ is straight, such an extension trivially exists, as $n_F$ is a constant function. Otherwise, if $F \subset \partial\Omega_p$ for some $p \in \{1, \dots, P\}$, the existence of $\tilde{n}_F$ is ensured by Proposition B.1. In both case, we have $\tilde{n}_F \in C^m(\bar{K}, \mathbb{R}^2)$ with $\|n_F\|_{m,\infty,K} \lesssim 1$, so that $\|\chi\|_{j+1,K} \lesssim \|\phi\|_{j+2,K}$. Now, using Proposition B.1 again, we have

$$ \int_F |(\chi - \pi_H \chi)v| = \int_{\hat{F}} |\phi'_F| |(\hat{\chi} - \pi_{\hat{F}} \hat{\chi})\hat{v}| = \int_{\hat{F}} |\phi'_F| |(\hat{\chi} - \pi_{\hat{F}} \hat{\chi})\hat{v}| \\ \le \sup_{s \in (0,1)} |\phi'_F| \int_{\hat{F}} |(\hat{\chi} - \pi_{\hat{F}} \hat{\chi})\hat{v}| \lesssim H \int_{\hat{F}} |(\hat{\chi} - \pi_{\hat{F}} \hat{\chi})\hat{v}|. $$

Then, arguing as in Lemma 3 of [13], we obtain that

$$ \int_{\hat{E}} |(\hat{\chi} - \pi_{\hat{E}} \hat{\chi})\hat{v}| \lesssim |\hat{\chi}|_{j+1, \hat{K}} |\hat{v}|_{1, K}, $$

and Lemma 2.3 of [7] shows that $|\hat{\chi}|_{j+1, \hat{K}} \lesssim H^j \|\chi\|_{j+1, K}$ and $|\hat{v}|_{1, \hat{K}} \lesssim \|v\|_{1, K}$, so that

$$ \int_E |(\mu - \pi_h \mu)v| \lesssim H^{j+1} \|\chi\|_{j+1,K} \|v\|_{1,K} \lesssim H^{j+1} \|\phi\|_{j+2,K} \|v\|_{V(K),\omega}, $$

and the result follows by summation, and the definition of norm $\| \cdot \|_{\Lambda, \omega}$ in (3.5). $\square$

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

## REFERENCES

[1] M. Ainsworth, Discrete dipersion relation for hp-version finite element approximation at high wave number, SIAM J. Numer. Anal. **42** (2004), no. 2, 553–575.

[2] M. Amara, R. Djellouli, and C. Farhat, Convergence analysis of a discontinuous galerkin method with plane waves and lagrange multipliers for the solution of Helmholtz problems, SIAM J. Numer. Anal. **47** (2009), 1038–1066.

[3] R. Araya, C. Harder, D. Paredes, and F. Valentin, Multiscale hybrid-mixed method, SIAM J. Numer. Anal. **51** (2013), no. 6, 3505–3531.

[4] I. Babuška and S. Sauter, Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers?, SIAM REVIEW **42** (2000), no. 3, 451–484.

[5] H. Barucq, T. Chaumont Frelet, and C. Gout, Stability analysis of heterogeneous Helmholtz problems and finite element solution based on propagation media approximation, Math. Comp. **86** (2016), no. 307, 2129–2157.

[6] J.P. Bérenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comp. Phys. **114** (1994), 185–200.

[7] C. Bernardi, Optimal finite-element interpolation on curved domains, SIAM J. Numer. Anal. **26** (1989), no. 5, 1212–1240.

[8] D.L. Brown, D.Gallistl, and D. Peterseim, Multiscale Petrov-Galerkin method for high-frequency heterogeneous Helmholtz equations, Meshfree Methods for Partial Differential Equations VIII (2017), 85–115.

[9] O. Cessenat and B. Despres, Application of an ultra weak variational formulation of elliptic pdes to the two-dimensional Helmholtz problem, SIAM J. Numer. Anal. **35** (1998), no. 1, 255–299.

[10] T. Chaumont-Frelet and S. Nicaise, Wavenumber explicit convergence analysis for finite element discretizations of general wave propagation problems, submitted to IMA J. Numer. Anal.

[11] ———, High-frequency behaviour of corner singularities in Helmholtz problems, ESAIM: Math. Model. Numer. Anal. **52** (2018), 1803–1845.

[12] P.G. Ciarlet, The finite element method for elliptic problems, SIAM, 1978.

[13] M. Crouzeix and P.A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations I, R.A.I.R.O **7** (1973), no. R3, 33–76.
[14] B. Engquist and A. Majda, *Absorbing boundary conditions for numerical simulation of waves*, Proc. Natl. Acad. Sci. USA **74** (1977), no. 5, 1765-1766.

[15] C. Farhat, I. Harari, and L.P. Franca, *The discontinuous enrichment method*, Comput. Methods Appl. Mech. Engrg. **190** (2001), 6455-6479.

[16] T. Chaumont Frelet, *Finite element approximation of Helmholtz problems and application to seismic wave propagation*, Ph.D. thesis, INSA Rouen, 2015.

[17] ——, *On high order methods for the Helmholtz equation in highly heterogeneous media*, Comput. Math. App. **76** (2016), 2203-2225.

[18] J. Galkowski, E.A. Spence, and J. Wunsch, *Optimal constants in nontrapping resolvent estimates and applications in numerical analysis*, arXiv:1810.13426v2, 2018.

[19] D. Gallistl and D. Peterseim, *Stable multiscale Petrov-Galerkin finite element method for high frequency acoustic scattering*, Comput. Methods Appl. Engrg. **295** (2015), 1-17.

[20] D. Givoli, *High-order local non-reflection boundary conditions: a review*, Wave Motion **39** (2004), 319-326.

[21] C. Harder, A.L. Madureira, and F. Valentin, *A hybrid-mixed method for elasticity*, ESAIM: Math. Model. Num. Anal. **50** (2016), no. 2, 311-336.

[22] C. Harder, D. Paredes, and F. Valentin, *A family of multiscale hybrid-mixed finite element methods for the darcy equation with rough coefficients*, J. of Comput. Phys. **245** (2013), 107-130.

[23] C. Harder, D. Paredes, and F. Valentin, *On a multiscale hybrid-mixed method for advective-reactive dominated problems with heterogeneous coefficients*, SIAM Multiscale Model. and Simul. **13** (2015), no. 2, 491-518.

[24] C. Harder and F. Valentin, *Foundations of the MHM method*, Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations (G. R. Barrenechea, F. Brezzi, A. Cangiani, and E. H. Georgoulis, eds.), Lecture Notes in Computational Science and Engineering, Springer, 2016.

[25] U. Hetmaniuk, *Stability estimates for a class of helmholtz problems*, Commun. Math. Sci. **5** (2007), no. 3, 665-678.

[26] F. Ihlenburg and I. Babuška, *Finite element solution of the helmholtz equation with high wave number part i: the h-version of the fem*, Comp. Math. Appl. **30** (1995), no. 9, 9-37.

[27] ——, *Finite element solution of the helmholtz equation with high wave number part ii: the h-p version of the fem*, SIAM J. Numer. Anal. **34** (1997), no. 1, 315-358.

[28] L.-M. Imbert-Gérard and B. Després, *A generalized plane wave numerical method for smooth non constant coefficients*, IMA J. Numer. Anal. **34** (2014), 1072-1103.

[29] S. Lanteri, D. Paredes, C. Scheid, and F. Valentin, *The multiscale hybrid-mixed method for the maxwell equations in heterogeneous media*, SIAM Multiscale Model. Simul. **16** (2018), no. 4, 1648-1683.

[30] G.S. Martin, R. Wiley, and J. Manfurt, *Marmousi2: An elastic upgrade for marmousi*, The Leading Edge **25** (2006), no. 2, 156-166.

[31] J.M. Melenk and S. Sauter, *Wavenumber explicit convergence analysis for galerkin discretizations of the Helmholtz equation*, SIAM J. Numer. Anal. **49** (2011), no. 3, 1210-1243.

[32] A. Moiola and E.A. Spence, *Acoustic transmission problems: wavenumber-explicit bounds and resonance-free regoins*, Math. Models Meth. Appl. Sci., in press, 2018.

[33] P. Monk and D.Q. Wang, *A least-squares method for the Helmholtz equation*, Comput. Meth. Appl. Mech. Engrg. (1999), no. 175, 121-136.

[34] M. Ohlberger and B. Verfürth, *A new heterogeneous multiscale method for the Helmholtz equation with high contrast*, Multiscale Model. Simul. **16** (2018), no. 1, 385-411.

[35] D. Paredes, F. Valentin, and H. M. Versieux, *On the robustness of multiscale hybrid-mixed methods*, Math. Comp. **86** (2017), no. 304, 525-548.

[36] L.E. Payne and H.F. Weinberger, *An optimal poincaré inequality for convex domains*, Rational Mechanics and Analysis **5** (1960), no. 1, 286-192.

[37] P.A. Raviart and J.M. Thomas, *Primal hybrid finite element methods for 2nd order elliptic equations*, Math. Comp. **31** (1977), no. 138, 391-413.

[38] H. Riou, P. Ladeveze, and B. Sourcis, *The multiscale vtcr approach applied to acoustics problems*, J. Comput. Acous. **16** (2008), no. 4, 487-505.

[39] A.H. Schatz, *An observation concerning ritz-galerkin methods with indefinite bilinear forms*, Math. Comp. **28** (1974), no. 128, 959-962.

[40] T. Strouboulis, I. Babuška, and R. Hidakat, *The generalized finite element method for Helmholtz equation: Theory, computation, and open problems*, Comput. Methods Appl. Mech. Engrg. **195** (2006), 4711-4731.

[41] R. Tezaur, I. Kalashnikova, and C. Farhat, *The discontinuous enrichment method for medium-frequency Helmholtz problems with a spatially variable wavenumber*, Comput. Methods Appl. Mech. Engrg. **268** (2013), 126-140.