# Efficient approximation of high-frequency Helmholtz solutions by Gaussian coherent states

Théophile Chaumont-Frelet, Victorita Dolean, Maxime Ingremeau

► To cite this version:

Théophile Chaumont-Frelet, Victorita Dolean, Maxime Ingremeau. Efficient approximation of high-frequency Helmholtz solutions by Gaussian coherent states. CMAM 2022, Aug 2022, Vienna, Austria. hal-03940689

HAL Id: hal-03940689

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

Submitted on 16 Jan 2023

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Efficient approximation of high-frequency Helmholtz solutions by Gaussian coherent states

Théophile Chaumont-Frelet, Victorita Dolean, Maxime Ingremeau

► To cite this version:

Théophile Chaumont-Frelet, Victorita Dolean, Maxime Ingremeau. Efficient approximation of high-frequency Helmholtz solutions by Gaussian coherent states. CMAM 2022, Aug 2022, Vienna, Austria. hal-03940689v2

HAL Id: hal-03940689

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

Submitted on 2 Feb 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Efficient approximation of high-frequency Helmholtz solutions
by Gaussian coherent states

T. Chaumont-Frelet*,†, V. Dolean†,‡ and M. Ingremeau‡

CMAM 2022

MS11: Numerical methods for wave propagation problems

TU Wien, August 29

*Inria Université Côte d'Azur, Atlantis project-team

†Laboratoire J.A. Dieudonné

‡University of Strathclyde, department of mathematics and statistics
EFFICIENT APPROXIMATION OF HIGH-FREQUENCY HELMHOLTZ
SOLUTIONS BY GAUSSIAN COHERENT STATES

T. CHAUMONT-FRELET†, V. DOLEAN‡, AND M. INGREMEAU‡

ABSTRACT. We introduce new finite-dimensional spaces specifically designed to approximate the solutions to high-frequency Helmholtz problems with smooth variable coefficients in dimension d. These discretization spaces are spanned by Gaussian coherent states, that have the key property to be localised in phase space. We carefully select the Gaussian coherent states spanning the approximation space by exploiting the (known) micro-localisation properties of the solution. For a large class of source terms (including plane-wave scattering problems), this choice leads to discrete spaces that provide a uniform approximation error for all wavenumber k with a number of degrees of freedom scaling as $k^{d-1/2}$, which we rigorously establish. In comparison, for discretization spaces based on (piecewise) polynomials, the number of degrees of freedom has to scale at least as $k^d$ to achieve the same property. These theoretical results are illustrated by one-dimensional numerical examples, where the proposed discretization spaces are coupled with a least-squares variational formulation.

KEY WORDS. Gabor frames, Helmholtz equation, high-frequency problems

1. INTRODUCTION

Time-harmonic wave propagation is a mechanism at the center of a large number of physical and industrial applications. We may cite, among many, radar imaging [9], or seismic prospection [42]. In practice, numerical methods are required to approximately simulate the propagation of waves, and although several methods are available, it is still very challenging to compute accurate approximations in the high-frequency regime.

Here, we consider the scalar Helmholtz equation, which is probably the simplest model for this kind of problems. Specifically, given a compactly supported right-hand side $f: \mathbb{R}^d \to \mathbb{C}$, our model problem is to find $u: \mathbb{R}^d \to \mathbb{C}$ such that

$$ (1.1a) \qquad -k^2 \mu u - \nabla \cdot (A \nabla u) = f \quad \text{in } \mathbb{R}^d, $$

where $\mu$ and $A$ are (given) smooth coefficients that are respectively equal to 1 and $I$ outside a ball of radius $R > 0$, and $k > 0$ is the (given) wavenumber. This equation is supplemented with the Sommerfeld radiation condition at infinity. Namely, we require that

$$ (1.1b) \qquad \frac{\partial u}{\partial |\mathbf{x}|} (\mathbf{x}) - iku(\mathbf{x}) = o\left(|\mathbf{x}|^{-(d+1)/2}\right) \quad \text{as } |\mathbf{x}| \to +\infty. $$

A particularly important scenario covered by (1.1) is the scattering of a plane-wave, where the right-hand side takes the form

$$ (1.2) \qquad f := (k^2 \mu + \nabla \cdot (A \nabla \cdot)) e^{ikd \cdot x}, $$

where $d \in \mathbb{R}^d$, $|d| = 1$ is the incidence direction (such right-hand sides are compactly supported due to the assumptions on $\mu$ and $A$).

As we propose a “volumic” method, we will actually replace the Sommerfeld radiation condition (1.1b) by a Perfectly Matched Layer (PML). This approach is entirely standard [2, 7, 15], and amounts to slightly modifying the coefficients $\mu$ and $A$ in (1.1a). This process is detailed in Section 4.2. We therefore do not directly discretize (1.1), but rather, an alternative version where the coefficients have been modified away from the origin in (1.1a) and the Sommerfeld radiation

<footnote>†University Côte d'Azur, Inria, CNRS, LJAD</footnote>

<footnote>‡University Côte d'Azur, CNRS, LJAD and University of Strathclyde</footnote>

<footnote>‡University Côte d'Azur, CNRS, LJAD</footnote>
Given $f : \mathbb{R}^d \to \mathbb{C}$, find $u : \mathbb{R}^d \to \mathbb{C}$ such that

$$ -k^2 \mu u - \nabla \cdot (A \nabla u) = f, $$

where $\mu$ and $A$ are smooth (given) coefficients, and $k$ is the (large) frequency.
condition (1.1b) replaced by the requirement that the solution sits in $H^1(\mathbb{R}^d)$. Crucially, the solutions to both problem coincide close to the origin. In what follows, we work with modified coefficients $\mu$ and $\mathbf{A}$ that incorporate a PML, but keep the same notation for simplicity.

In this work, we investigate the use of discretization spaces based on Gaussian coherent states (GCS), that is, functions of the form

$$ \Phi_{k, x_0, \xi_0}(\mathbf{x}) := \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k}{2}|\mathbf{x}-x_0|^2} e^{-ik\xi_0 \cdot (\mathbf{x}-x_0)}, $$

where $x_0, \xi_0 \in \mathbb{R}^d$ are user-selected parameters. The idea of decomposing a function as a discrete sum of Gaussian coherent states goes back to [14]. Here, following [5, 8], we focus on a lattice of phase-space points $[\mathbf{x}^{k,m}, \xi^{k,n}] := \sqrt{(kR)^{-1}\pi}[\mathbf{m}, \mathbf{n}], [\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d}$, that are spaced by $\sim (kR)^{-1/2}$, and we write

$$ \Psi_{k,m,n} := \Phi_{k,x^{k,m},\xi^{k,n}}. $$

Then, our a discretization space is of the form

$$ W_{\Lambda} := \mathrm{Vect}\{\Psi_{k,m,n}; [\mathbf{m}, \mathbf{n}] \in \Lambda\}, $$

where $\Lambda \subset \mathbb{Z}^{2d}$ is a carefully chosen set of indices.

For the sake of simplicity, we will assume in the introduction that the domain is non-trapping. Our first result is that if $\Lambda$ is chosen as

$$ \Lambda := \{[\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d}; |\mathbf{x}^{k,m}|^2 + |\xi^{k,n}|^2 \le \rho k R\}, $$

for $\rho > 0$, we have

$$ (1.3) \quad \dim W_{\Lambda} \simeq \rho^d (kR)^d \quad \text{and} \quad k^2 \min_{w \in W_{\Lambda}} \|u-w\|_{\hat{H}^1(\mathbb{R}^d)} \le C \rho^{-1/2} \|f\|_{L^2(\mathbb{R}^d)}, $$

for a general $f \in L^2(\mathbb{R}^d)$ with $\mathrm{supp} f \subset B(0, R)$, where $\|\cdot\|_{\hat{H}^1(\mathbb{R}^d)}$ is a $H^1(\mathbb{R}^d)$-norm including a weight at infinity (see (4.1) below). As we describe in more length afterwards, this result is not very impressive on its own. Specifically, it is a standard approximation result similar to polynomial approximations: we need a fixed number of points per wavelength to achieve a constant accuracy. Our second result, which is key, deals with the case where $f$ takes the particular form (1.2). In this case, we select

$$ \Lambda := \{[\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d}; |p(\mathbf{x}^{k,m}, \xi^{k,n})| \le (kR)^{-1/2+\epsilon}\} $$

where $\epsilon > 0$ can be selected arbitrarily small and

$$ (1.4) \quad p(\mathbf{x}, \xi) := \mathbf{A}(\mathbf{x})\xi \cdot \xi - \mu(\mathbf{x}) \quad \forall \mathbf{x}, \xi \in \mathbb{R}^d, $$

is the principal symbol associated with the differential operator in (1.1), where the coefficients include a PML. We then have

$$ (1.5) \quad \dim W_{\Lambda} \simeq (kR)^{d-1/2+\epsilon} \quad \text{and} \quad \min_{w \in W_{\Lambda}} \|u-w\|_{\hat{H}^1(\mathbb{R}^d)} \le C_{\epsilon,m}(kR)^{-m} \quad \forall m \in \mathbb{N}. $$

It means that for right-hand sides corresponding to scattering problems (and actually, a wider family of right-hand sides), Gaussian coherent states provide an accurate solution with $O((kR)^{d-1/2+\epsilon})$ DOFs. In fact, the convergence is even super-algebraic as the frequency increases.

To put (1.3) and (1.5) into perspective, we compare them with other standard methods. Actually, there are several options to numerically solve (1.1) (either with the Sommerfeld condition (1.1b) or with a PML approximation), that we review below.

The most versatile approach is probably the finite element method (FEM). The method hinges on a triangulation of the domain into elements of size $h$, and piecewise polynomial basis functions of degree $p$. It can be shown that if $p$ grows logarithmically with $k$, then the condition that $kh/p$ is constant provides (at least) a constant accuracy as $k$ increases [30, 34, 35]. As a result, high-order FEM essentially requires $O((kR)^d)$ degrees of freedom (DOFs) to achieve a constant accuracy. The resulting matrix is sparse.

Trefftz-like methods are similar to FEM in that they also rely on a mesh of the domain, but the polynomial shape functions are replaced by local solutions to the Helmholtz problem, such as
Given $f: \mathbb{R}^d \to \mathbb{C}$, find $u: \mathbb{R}^d \to \mathbb{C}$ such that

$$ -k^2 \mu u - \nabla \cdot (A \nabla u) = f, $$

where $\mu$ and $A$ are smooth (given) coefficients, and $k$ is the (large) frequency.

A central ingredient of numerical methods is to design a finite-dimensional space $W$ such that

$$ \min_{w \in W} \|u - w\| $$

is small in a suitable norm $\|\cdot\|$.
<table><caption>TABLE 1. Comparison of commonly used discretization techniques</caption><thead><tr><th></th><th>FEM</th><th>GCS</th><th>BEM</th><th>asymptotic</th></tr></thead><tbody><tr><td>Cost</td><td>(kR)<sup>d</sup></td><td>(kR)<sup>d-1/2</sup></td><td>(kR)<sup>d-1</sup></td><td>1</td></tr><tr><td>Heterogeneous media</td><td>yes</td><td>yes</td><td>no</td><td>yes</td></tr><tr><td>Convergence at fixed k</td><td>yes</td><td>yes</td><td>yes</td><td>no</td></tr></tbody></table>

plane-waves [25], or generalised plane-waves [27, 28]. There are many ways to “glue” these local solution together, including partition of unity methods [32], least squares methods [36], the ultra weak variational method [3], the discontinuous enrichment method [13] or the variational theory of complex rays [39]. While these methods typically induce a large reduction of the number of DOFs as compared to FEM, they usually still need at least $O((kR)^d)$ DOFs, see, e.g., [6, 18, 24]. Similar to FEM, the resulting matrix is sparse.

The next family of methods we want to mention are boundary element methods (BEM) [40]. These methods rely on boundary integral equations which, strictly speaking, are not available for smoothly varying coefficients, since the expression of Green’s function must be available. It is nevertheless interesting to include them in the comparison. These methods typically provide a constant accuracy with only $O((kR)^{d-1})$ DOFs [16]. However, the resulting matrix is dense and its entries are costly to compute. These issues can be mitigated using compression techniques, such as the fast multi-pole method [21] or hierarchical matrices [23].

Finally, asymptotic methods rely on the fact that when the frequency is very large, it is sometimes possible to simplify the search for the solution of the Helmholtz equation and the properties of the solution itself to computations involving only the underlying classical dynamics [11]. This is done using tools of semi-classical analysis, such as the WKB method and can lead to discrete problems with a number of DOFs independent of $k$. The main drawback of these approaches is that they are only asymptotically valid: they do not converge for fixed value of $k$. Besides, it is not always clear from which range of $k$ they are relevant.

As compared to FEM, the proposed GCS method thus gains “half a dimension” at high-frequencies, but it is still half a dimension higher than BEM. As compared to BEM however, our methodology has the advantage to apply in a generic framework where the Green’s function is not available. Another important comment is that in the (very) high-frequency regime, our method is more expensive than asymptotic methods. However, asymptotic methods cannot converge at fixed $k$, which our method does. This is summarized in Table 1.

In addition to the approximability results (1.3) and (1.5), we also present a least-squares method based on Gaussian coherent states for Problem (1.1). As we show, the convergence of the method is easily established. Besides, although the matrix is dense, we show that the entries decay super-algebraically away from the diagonal. As a result, the matrix is essentially banded, and we believe that efficient linear solvers can be devised. This will be analysed in more depth in future works.

We finally present a set of one-dimensional numerical experiments using the proposed least-squares method. Although the setting is rather simple, the results perfectly fit the theoretical analysis and readily shows that proposed approach allows for a drastic reduction of the number of DOFs in the high-frequency regime.

To the best of our knowledge, our micro-locally adapted spaces of Gaussian coherent states appear to be entirely original, but we would like to mention that similar basis functions have already been employed to discretize PDE problems. In particular, generalised coherent states like Hagedorn wavepackets were used to describe the solution of time-dependent Schrödinger equation in [12, 19, 20, 31].

We would like to emphasize that the main goal of our work is to understand how many DOFs are required to maintain a constant accuracy as $k$ increases, rather than deriving precise convergence rates in terms of DOFs for a fixed frequency. For the sake of completeness, we nevertheless provide a quick comparison here. Assuming that $k$ is fixed and that the right-hand side $f$ is smooth, the results in [5] (see also [22, Chapter 11] and Corollary 4.7 below) show that the convergence to the solution will be super-algebraic in terms of the number of DOFs. Without further assumption
Given $f: \mathbb{R}^d \to \mathbb{C}$, find $u: \mathbb{R}^d \to \mathbb{C}$ such that

$$ -k^2 \mu u - \nabla \cdot (A \nabla u) = f, $$

where $\mu$ and $A$ are smooth (given) coefficients, and $k$ is the (large) frequency.

A central ingredient of numerical methods is to design a finite-dimensional space $W$ such that

$$ \min_{w \in W} \|u - w\| $$

is small in a suitable norm $\|\cdot\|$.

Here, we would especially like to design a sequence of spaces $W(k)$ so that

$$ \min_{w \in W(k)} \|u - w\| = Cst $$

uniformly in $k$ while minimizing the dimension of $W(k)$.
on $f$, such rates are equivalent to hp finite element and Trefftz methods. Under the additional assumption that the solution $u$ is (piecewise) analytic, hp finite element and Trefftz methods can achieve an exponential convergence rate in terms of DOFs [41], with an improved rate for Trefftz methods as compared to finite element methods [24, 33]. To the best our knowledge, it is an open question whether such a result holds true for the present method.

The remainder of our work is organised as follows. In Section 2, we precise the setting and state our key approximation result (1.5) in its most general form. Section 3 contains the proof of our findings. In Section 4, we apply the general theory of Sections 2 and 3 to our scattering model problem with PML. Numerical examples are reported in Section 5. Finally, Appendix A collects technical results concerning Gaussian coherent states.

## 2. SETTING AND MAIN RESULTS

### 2.1. Notations.
Throughout this work $\hbar \in \mathcal{H} \subset (0, 1]$ will denote a small parameter. When applying our general results to the Helmholtz equation, we will have $\hbar \sim (kR)^{-1}$, so that considering the set $(0, 1]$ amounts to ignoring low frequencies, and focusing on high frequencies when $\hbar \to 0$. For the sake of generality, we restrict our analysis to a subset $\mathcal{H} \subset (0, 1]$ for reasons that will become apparent in Section 4. Notice that the case $\mathcal{H} = (0, 1]$ is not excluded.

#### 2.1.1. Basic notation.
The canonical basis of $\mathbb{R}^d$ or of $\mathbb{C}^d$ will be denoted by $(e_1, ..., e_d)$. If $\mathbf{x}, \mathbf{y} \in \mathbb{C}^d$, we write

$$ \mathbf{x} \cdot \mathbf{y} := \sum_{j=1}^{d} x_j y_j $$

without complex conjugation on the second argument, and $|\mathbf{x}| = (\mathbf{x} \cdot \bar{\mathbf{x}})^{1/2}$ is the usual Euclidean norm.

For a multi-index $\alpha \in \mathbb{N}^d$, $[\alpha] := \alpha_1 + \cdots + \alpha_d$ denotes its usual $\ell_1$ norm. If $v : \mathbb{R}^d \to \mathbb{C}$, the notation

$$ \partial^{\alpha} v := \frac{\partial^{\alpha_1}}{\partial x_1} \cdots \frac{\partial^{\alpha_d}}{\partial x_d} v $$

is employed for the partial derivatives in the sense of distributions, whereas $\mathbf{x}^{\alpha} := x_1^{\alpha_1} \cdots x_d^{\alpha_d}$. Finally, if $\beta \in \mathbb{N}^d$ is another multi-index, we will sometimes need the notation

$$ \left( \begin{array}{c} \alpha \\ \beta \end{array} \right) = \prod_{j=1}^{d} \left( \begin{array}{c} \alpha_j \\ \beta_j \end{array} \right), $$

and the notation $\alpha \le \beta$ means that $\alpha_j \le \beta_j$ for all $j \in \{1, \dots, d\}$.

If $\mathbf{n} \in \mathbb{Z}^d$, we employ the notation $|\mathbf{n}|^2 := n_1^2 + \cdots + n_d^2$ for its $\ell_2$ norm. Finally, if $\Lambda \subset \mathbb{Z}^{2d}$, $\ell^2(\Lambda)$ has its usual definition, and we denote by $\| \cdot \|_{\ell^2(\Lambda)}$ its usual norm.

#### 2.1.2. Key functional spaces.
In what follows, $L^2(\mathbb{R}^d)$ is the usual Lebesgue space of complex-valued square integrable functions over $\mathbb{R}^d$. The usual norm and inner products of $L^2(\mathbb{R}^d)$ are respectively denoted by $\| \cdot \|_{L^2(\mathbb{R}^d)}$ and $(\cdot, \cdot)$.

Since we are dealing with the (unbounded) $\mathbb{R}^d$ space, following [5], our analysis will require the weighted Sobolev spaces

$$ \hat{H}^p(\mathbb{R}^d) := \{ v \in L^2(\mathbb{R}^d) \mid \mathbf{x}^{\alpha} \partial^{\beta} v \in L^2(\mathbb{R}^d) \quad \forall \alpha, \beta \in \mathbb{N}^d; [\alpha], [\beta] \le p \}, $$

that we equip with the family of equivalent $\hbar$-weighted norms given by

$$ \|v\|_{\hat{H}_{\hbar}^p(\mathbb{R}^d)}^2 := \sum_{[\alpha] \le p} \sum_{q \le p - [\alpha]} \hbar^{2[\alpha]} \|\mathbf{x}|^q \partial^{\alpha} v\|_{L^2(\mathbb{R}^d)}^2 $$

for all $p \in \mathbb{N}$.
Finite element discretization with a mesh of size $h \sim k^{-1}$

$\dim W(k) \sim k^d$, sparse matrix.
$C^0(\mathbb{R}^d)$ is the set of complex-valued continuous functions defined over $\mathbb{R}^d$, and $C^\ell(\mathbb{R}^d)$ is the set of functions $v: \mathbb{R}^d \to \mathbb{C}$ such that $\partial^\alpha v \in C^0(\mathbb{R}^d)$ for all $\alpha \in \mathbb{N}^d$ with $[\alpha] \le \ell$. We introduce the notation

$$ \|v\|_{C^\ell(\mathbb{R}^d)} := \max_{\substack{\alpha \in \mathbb{N}^d \\ [\alpha] \le \ell}} \max_{\mathbf{x} \in \mathbb{R}^d} |(\partial^\alpha v)(\mathbf{x})| \quad \forall v \in C^\ell(\mathbb{R}^d) $$

and $C_b^\ell(\mathbb{R}^d)$ is the subset of functions $v \in C^\ell(\mathbb{R}^d)$ such that $\|v\|_{C^\ell(\mathbb{R}^d)} < +\infty$. We also set

$$ C_b^\infty(\mathbb{R}^d) := \bigcap_{\ell \in \mathbb{N}} C_b^\ell(\mathbb{R}^d). $$

Finally, if $\Omega \subset \mathbb{R}^d$ is an open set, we denote by $C_c^\infty(\Omega)$ the set of smooth functions whose support is a compact subset of $\Omega$.

## 2.2. The frame of Gaussian coherent states.
The goal of this work is to efficiently approximate the solution $u_\hbar$ to the equation $P_\hbar u_\hbar = f$ with a finite span of Gaussian coherent states. For $[m, n] \in \mathbb{Z}^{2d}$, we thus consider the Gaussian state

$$ \Psi_{\hbar, m, n}(\mathbf{x}) := (\pi \hbar)^{-d/4} e^{-\frac{1}{2\hbar} |\mathbf{x} - \mathbf{x}^{\hbar, m}|^2} e^{\frac{i}{\hbar} \boldsymbol{\xi}^{\hbar, n} \cdot (\mathbf{x} - \mathbf{x}^{\hbar, m})}, $$

where $\mathbf{x}^{\hbar, m} := \sqrt{\pi \hbar} m$ and $\boldsymbol{\xi}^{\hbar, n} := \sqrt{\pi \hbar} n$. The family of Gaussian coherent states $(\Psi_{\hbar, m, n})_{[m, n] \in \mathbb{Z}^{2d}}$ actually forms a frame over $L^2(\mathbb{R}^d)$, meaning there exists two constants $0 < \alpha < \beta < +\infty$ solely depending on $d$ such that

$$ \alpha \|v\|_{L^2(\mathbb{R}^d)}^2 \le \sum_{[m,n] \in \mathbb{Z}^{2d}} |(v, \Psi_{\hbar,m,n})|^2 \le \beta \|v\|_{L^2(\mathbb{R}^d)}^2 \quad \forall v \in L^2(\mathbb{R}^d). $$

This result was first proved in [8], but the idea of decomposing a function as a discrete sum of Gaussian states goes back to [14], where it was proved that the span of $(\Psi_{\hbar, m, n})_{[m, n] \in \mathbb{Z}^{2d}}$ is dense in $L^2(\mathbb{R}^d)$.

Actually, the frame property implies that there exists another family of functions $(\Psi_{\hbar, m, n}^*)_{[m, n] \in \mathbb{Z}^{2d}}$ called the dual frame such that

$$ (2.1) \qquad v = \sum_{[m,n] \in \mathbb{Z}^{2d}} (v, \Psi_{\hbar,m,n}^*) \Psi_{\hbar,m,n} $$

for all $v \in L^2(\mathbb{R}^d)$. It is thus clear that any $v \in L^2(\mathbb{R}^d)$ may be well-approximated by (a large number of) Gaussian states. As we are going to develop hereafter, when considering the solution to a high-frequency PDE problem, a good approximation may be obtained with few Gaussian states, by carefully selecting the indices $[m, n]$ in (2.1).

**Remark 2.1** (General expansions in the Gaussian frame). The family $(\Psi_{\hbar, m, n})$ is not a Riesz basis, so that the expansion (2.1) of $v$ as a sum of $\Psi_{\hbar, m, n}$ is not unique. However, a crucial property of (2.1) is that this expansion is stable, in the sense that

$$ \sum_{[m,n] \in \mathbb{Z}^{2d}} |(v, \Psi_{\hbar,m,n}^*)|^2 \le \gamma \|v\|_{L^2(\mathbb{R}^d)}^2, $$

where $\gamma$ only depends on $d$, the dual frame being itself a frame. This is especially important at the numerical level in the presence of round-off errors [1].

## 2.3. Settings and key assumptions.
Throughout this work, we consider a second order differential operator on $\mathbb{R}^d$ depending on $\hbar$, and taking the form

$$ (2.2) \qquad (P_\hbar v)(\mathbf{x}) = \hbar^2 \sum_{j, \ell=1}^{d} a_{j\ell}^\hbar(\mathbf{x}) \frac{\partial^2 v}{\partial x_j \partial x_\ell}(\mathbf{x}) + i\hbar \sum_{j=1}^{d} b_j^\hbar(\mathbf{x}) \frac{\partial v}{\partial x_j}(\mathbf{x}) + c^\hbar(\mathbf{x})v(\mathbf{x}), $$

where $a_{j\ell}^\hbar, b_j^\hbar, c^\hbar \in C_b^\infty(\mathbb{R}^d)$ for $1 \le j, \ell \le d$. For the sake of simplicity, we introduce

$$ C_{\text{coef},p} := \sup_{\hbar \in \mathcal{H}} \left( \sum_{j, \ell=1}^{d} \|a_{j, \ell}^{\hbar}\|_{C^p(\mathbb{R}^d)} + \sum_{j=1}^{d} \|b_j^{\hbar}\|_{C^p(\mathbb{R}^d)} + \|c^{\hbar}\|_{C^p(\mathbb{R}^d)} \right) \quad \forall p \in \mathbb{N}, $$
Finite element discretization with a mesh of size $h \sim k^{-1}$

$\dim W(k) \sim k^d$, sparse matrix.

Boundary elements with a surface mesh of size $h \sim k^{-1}$

$\dim W(k) \sim k^{d-1}$, dense matrix.
and assume that $C_{\text{coef},p} < +\infty$ for all $p \in \mathbb{N}$.

The principal symbol of $P_h$ is the function $p_h \in C^\infty(\mathbb{R}^{2d})$ defined by

$$ (2.3) \qquad p_h(\mathbf{x}, \boldsymbol{\xi}) := \sum_{j, \ell=1}^{d} a_{j, \ell}^{\hbar}(\mathbf{x}) \xi_i \xi_j + \sum_{j=1}^{d} b_j^{\hbar}(\mathbf{x}) \xi_j + c^{\hbar}(\mathbf{x}). $$

For the sake of shortness, we will often write $p_h(\mathbf{m}, \mathbf{n}) := p_h(\mathbf{x}^{\hbar, m}, \boldsymbol{\xi}^{\hbar, n})$ for $[\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d}$.

**Remark 2.2 (h-dependent symbol).** When considering the standard Helmholtz differential operator $P_h v := -v - \hbar^2 \Delta v$ the symbol simply reads $p(\mathbf{x}, \boldsymbol{\xi}) = 1 - |\boldsymbol{\xi}|^2$ and in particular, it does not depend on $\hbar$. It is however interesting to allow for a mild dependence on $\hbar$ in the symbol. This is for instance the case when considering Helmholtz problems in the presence of dissipative materials. In this case, the differential operator reads $P_h v := -v - i\gamma\hbar v - \hbar^2 \Delta v$, where the function $\gamma \ge 0$ models the dissipation. One readily see in this case that the symbol is $p_h(\mathbf{x}, \boldsymbol{\xi}) = 1 - i\gamma(\mathbf{x})\hbar - |\boldsymbol{\xi}|^2$.

**Remark 2.3 (non-divergence form).** In contrast to (1.1) in the introduction, we present our model PDE in non-divergence form in (2.2). The key reason behind this choice is that it makes the link between the PDE operator and its symbol simpler. Notice that since we are only considering smooth coefficients in this work, this choice is absolutely not restrictive, and (1.1) can be easily recast into (2.2).

Along with the smoothness of the coefficients, we make two key assumptions. First, we assume that $P_h : \hat{H}^p(\mathbb{R}^d) \to \hat{H}^p(\mathbb{R}^d)$ is invertible with the norm of $P_h^{-1}$ being polynomially bounded in $\hbar^{-1}$. Specifically, we assume that for all $f \in L^2(\mathbb{R}^d)$ there exists a unique $u_h \in L^2(\mathbb{R}^d)$ such that $P_h u_h = f$. In addition, there exists $N \in \mathbb{N}$ such that for all $p \in \mathbb{N}$, if $f \in \hat{H}^p(\mathbb{R}^d)$, then $u_h \in \hat{H}^p(\mathbb{R}^d)$ with

$$ (2.4) \qquad \|u_h\|_{\hat{H}_0^p(\mathbb{R}^d)} \le C_{\text{sol},p} \hbar^{-N} \|f\|_{\hat{H}_h^p(\mathbb{R}^d)} \quad \forall \hbar \in \mathcal{H} $$

for some constant $C_{\text{sol},p}$ independent of $\hbar$.

In the context of high-frequency scattering, these assumptions are reasonable and hold in a variety of situations. We refer the reader to Remark 4.1 in Section 4.2 below where we expand on that aspect.

Our second assumption is that there exists a value $\delta_0 > 0$ such that

$$ (2.5) \qquad \exists D_0 > 0 \text{ such that } \forall \hbar \in \mathcal{H}, \{(\mathbf{x}, \boldsymbol{\xi}) \in \mathbb{R}^{2d}; |p_h(\mathbf{x}, \boldsymbol{\xi})| < \delta_0\} \subset B(0, D_0). $$

In the remainder of this work, we allow generic constants $C$ to depend on $\{C_{\text{coef},p}\}_{p \in \mathbb{N}}$, $\{C_{\text{sol},p}\}_{p \in \mathbb{N}}$, $N$ and $D_0$. We also employ the notation $C_{\alpha,\beta,...}$ if the constant $C$ is additionally allowed to depend on other previously introduced quantity $\alpha, \beta, ...$.

**2.4. Statement of the approximability result.** Our main result is that, if $f$ is micro-localised near the set $\{(\mathbf{x}, \boldsymbol{\xi}) \in \mathbb{R}^{2d}; p_h(\mathbf{x}, \boldsymbol{\xi}) = 0\}$, then so is the solution $u_h$ to $P_h u_h = f$. This is a standard result when micro-localisation is understood in terms of pseudo-differential operators (see for instance [43, Theorem 6.4]), but here, by micro-localisation properties, we mean that $f$ can be approached by a linear combination of $\Psi_{h,m,n}$ with $[\mathbf{m}, \mathbf{n}]$ in a certain region of $\mathbb{Z}^{2d}$. Hence, our results may not be easily recovered from standard results in semiclassical analysis (see e.g. [43, Theorem 6.4]).

**Theorem 2.4 (Approximability for Gaussian state right-hand sides).** Let $0 < \varepsilon < 1/2$ and $0 \le \alpha < \delta_0/2$. For $\hbar \in \mathcal{H}$, consider the right-hand side

$$ f_h := \sum_{[\mathbf{m}, \mathbf{n}] \in \Lambda_{h, \text{rhs}}} F_{\mathbf{m}, \mathbf{n}}^{\hbar} \Psi_{h, \mathbf{m}, \mathbf{n}}. $$

where $F^\hbar \in \Lambda_{h, \text{rhs}}$ with

$$ (2.6) \qquad \Lambda_{h, \text{rhs}} := \{[\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d} \mid |p_h(\mathbf{x}^{\hbar, m}, \boldsymbol{\xi}^{\hbar, n})| \le \alpha + 2\hbar^{1/2-\varepsilon}\}. $$
Finite element discretization with a mesh of size $h \sim k^{-1}$

dim $W(k) \sim k^d$, sparse matrix.

Boundary elements with a surface mesh of size $h \sim k^{-1}$

dim $W(k) \sim k^{d-1}$, dense matrix.

Our goal is to propose a new family of spaces such that

dim $W(k) \sim k^{d-1/2}$, sparse matrix.
Then, if $u_h$ is the solution to $P_h u_h = f_h$, we have

$$ \left\| u_h - \sum_{[m,n] \in \Lambda_{h,\text{sol}}} (u_h, \Psi_{h,m,n}^* \Psi_{h,m,n}) \right\|_{\hat{H}_h^p(\mathbb{R}^d)} \le C_{\varepsilon,p,m} \hbar^m \|F_h\|_{\ell^2(\mathbb{Z}^{2d})} \quad \forall m \in \mathbb{N}. $$

for all $p \in \mathbb{N}$ with

$$ \Lambda_{h,\text{sol}} := \{[m,n] \in \mathbb{Z}^{2d} \mid |p_h(x^{h,m}, \xi^{h,n})| \le \alpha + 4\hbar^{1/2-\varepsilon}\}. $$

Notice that the index sets $\Lambda_{h,\text{rhs}}$ and $\Lambda_{h,\text{sol}}$ depend of $\varepsilon$, which explains why the constants in Theorem 2.4 and Corollary 2.5 depend on $\varepsilon$.

In practice, the right-hand side of the problem is not a finite linear combination of Gaussian coherent states. However, many right-hand sides of interest become well-approximated by such combination in the high-frequency regime. This is in particular the case when considering scattering by a plane-wave (see Lemma 4.4 below).

**Corollary 2.5** (Approximability for micro-localised right-hand sides). Let $p \ge 0$. Consider a set of right-hand sides $(f_h)_{h \in \mathcal{H}} \subset \hat{H}^p(\mathbb{R}^d)$ and assume that there exists a set of sequences $(F^h)_{h \in \mathcal{H}} \subset \ell^2(\mathbb{Z}^{2d})$ such that

$$ \begin{gathered} \|F^h\|_{\ell^2(\mathbb{Z}^{2d})} \le C \\ \|f_h - \sum_{[m,n] \in \Lambda_{h,\text{rhs}}} F_{m,n}^h \Psi_{h,m,n} \|_{\hat{H}_h^p(\mathbb{R}^d)} \le C_{\varepsilon,m} \hbar^m \quad \forall m \in \mathbb{N} \end{gathered} $$

for all $\hbar \in \mathcal{H}$. Then, we have

$$ \left\| u_h - \sum_{[m,n] \in \Lambda_{h,\text{sol}}} (u_h, \Psi_{h,m,n}^* \Psi_{k,m,n}) \right\|_{\hat{H}_h^{p+2}(\mathbb{R}^d)} \le C_{\varepsilon,p,m} \hbar^m \quad \forall m \in \mathbb{N}. $$

### 3. PROOF OF THEOREM 2.4

This section is devoted to the detailed proofs of Theorem 2.4 and Corollary 2.5.

**3.1. Preliminary results on Gaussian states.** We start by stating that the following bound

$$ (3.1) \qquad \|\Psi_{h,m,n}\|_{\hat{H}_h^s(\mathbb{R}^d)} \le C_s (1 + (\hbar^{1/2} |[m,n]|)^s) $$

holds true for all $s \in \mathbb{N}$, see [5, Lemma C.1]. We will also need the following expansion result.

**Proposition 3.1** (Tight expansion). For all $[m,n] \in \mathbb{Z}^{2d}$, there exists a sequence of coefficients $U^{m,n} \subset \ell^1(\mathbb{Z}^{2d})$ such that

$$ \|U^{m,n}\|_{\ell^p(\mathbb{Z}^{2d})} \le C_p \quad \forall p \in [1, +\infty], $$

and for all $\varepsilon \in (0, 1/2)$ and $s \in \mathbb{N}$, we have

$$ (3.2) \qquad \left\| \Psi_{h,m,n}^* - \sum_{\substack{[m',n'] \in \mathbb{Z}^{2d} \\ |[m,n] - [m',n']| \le \hbar^{-\varepsilon}}} U_{m',n'}^{m,n'} \Psi_{h,m',n'} \right\|_{\hat{H}_h^s(\mathbb{R}^d)} \le C_{\varepsilon,s,m} \hbar^m \quad \forall m \in \mathbb{N}. $$

*Proof.* We start by applying (2.1) to $v = \Psi_{h,m,n}^*$, leading to

$$ \Psi_{h,m,n}^* = \sum_{[m',n'] \in \mathbb{Z}^{2d}} (\Psi_{h,m,n}^*, \Psi_{h,m',n'}^*) \Psi_{k,m',n'}. $$

Next, we recall from [5, Proposition 4.2] that

$$ |(\Psi_{h,m,n}^*, \Psi_{h,m',n'}^*)| \le C e^{-|[m,n]-[m',n']|^{1/2}}. $$
<table><thead><tr><th></th><th>BEM</th><th>FEM</th></tr></thead><tbody><tr><td>Cost</td><td>$k^{d-1}$</td><td>$k^d$</td></tr><tr><td>Heterogeneous media</td><td>$\times$</td><td>✓</td></tr><tr><td>Sparse matrix</td><td>$\times$</td><td>✓</td></tr></tbody></table>

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

Galkowski and Spence, 2022

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

Melenk and Sauter, 2010
As a result, we define $U_{m',n'}^{m,n} := (\Psi_{\hbar,m,n}^*, \Psi_{\hbar,m',n'})$, so that $U^{m,n}$ indeed belongs to $\ell^p(\mathbb{Z}^{2d})$ for $1 \le p \le +\infty$, and

$$ \mathcal{E} := \Psi_{\hbar,m,n}^* - \sum_{\substack{[m',n'] \in \mathbb{Z}^{2d} \\ |[m,n]-[m',n']| \le \hbar^{-\epsilon}}} U_{m',n'}^{m,n} \Psi_{\hbar,m',n'} = \sum_{\substack{[m',n'] \in \mathbb{Z}^{2d} \\ |[m,n]-[m',n']| > \hbar^{-\epsilon}}} U_{m',n'}^{m,n} \Psi_{\hbar,m',n'}. $$

We then observe that

$$ \begin{aligned} \|\mathcal{E}\|_{\hat{H}_{\hbar}^s(\mathbb{R}^d)} & \le \sum_{\substack{[m',n'] \in \mathbb{Z}^{2d} \\ |[m,n]-[m',n']| > \hbar^{-\epsilon}}} |U_{m',n'}^{m,n}| \cdot \|\Psi_{\hbar,m',n'}\|_{\hat{H}_{\hbar}^s(\mathbb{R}^d)} \\ & \le C \sum_{\substack{[m',n'] \in \mathbb{Z}^{2d} \\ |[m,n]-[m',n']| > \hbar^{-\epsilon}}} (1 + |[m',n']|)^s e^{-|[m,n]-[m',n']|^{1/2}}, \end{aligned} $$

due to (3.1) and the fact that $\hbar \le 1$. Estimate (3.2) follows since

$$ \sum_{\substack{[m',n'] \in \mathbb{Z}^{2d} \\ |[m,n]-[m',n']| > \hbar^{-\epsilon}}} (1 + |[m',n']|)^s e^{-|[m,n]-[m',n']|^{1/2}} \le C_{\epsilon,s,m} \hbar^m. \quad \square $$

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

We close this section with two technical results. As we believe they are of independent interest, and because their proof require tedious computations, they are presented later in Appendix A.

**Proposition 3.2 (Quasi orthogonality).** Consider two sets of indices $\Lambda, \Lambda' \subset \mathbb{Z}^{2d}$ with

$$ \rho := \mathrm{dist}(\Lambda, \Lambda') > 0, $$

and such that there exists $\mu > 0$ with

$$ \Lambda \subset B(0, \mu). $$

Consider $L \in \mathbb{N}$, smooth coefficients $(A_a)_{a \in \mathbb{N}^d} \subset C_b^\infty(\mathbb{R}^d)$ and the differential operator

$$ P_{\hbar,L,A} := \sum_{\substack{\alpha \in \mathbb{N}^d \\ [\alpha] \le L}} \hbar^{[\alpha]} A_\alpha \partial^\alpha. $$

Then, for all $q > 0$ and $\alpha \in \mathbb{N}^d$, we have

$$ \begin{aligned} \sum_{[m,n] \in \Lambda} \sum_{[m',n'] \in \Lambda'} |(x^\alpha P_{\hbar,L,A} \Psi_{\hbar,m,n}, \Psi_{\hbar,m',n'})|^q & \le C_{L,A,q,m,\alpha} (1 + (\mu\hbar^{-1/2})^{([\alpha]+L)q}) |\Lambda|(1+\rho)^{-m} \\ & \le C'_{L,A,q,m,\alpha} (1 + (\mu\hbar^{-1/2})^{([\alpha]+L)q+2d}) (1+\rho)^{-m} \end{aligned} $$

for all $m \in \mathbb{N}$.

**Proposition 3.3 (Control of $(P_\hbar - p_\hbar)^L$).** For all $[m, n] \in \mathbb{Z}^{2d}$, we have

$$ (3.4a) \qquad \left\| (P_{\hbar} - p_{\hbar}(x^{k,m}, \xi^{k,n}))^{L} \Psi_{\hbar,m,n} \right\|_{L^2(\mathbb{R}^d)} \le C_L (1 + \hbar|n|^2)^{L/2} \hbar^{L/2}, $$

and

$$ (3.4b) \qquad \left\| (P_{\hbar}^{*} - \bar{p}_{\hbar}(x^{k,m}, \xi^{k,n}))^{L} \Psi_{\hbar,m,n} \right\|_{L^2(\mathbb{R}^d)} \le C_L (1 + \hbar|n|^2)^{L/2} \hbar^{L/2}. $$

The first part of Proposition 3.3 directly follows from Proposition A.6. In the second part, $P_\hbar^*$ denotes the formal adjoint of $P_\hbar$, which is also of the form (A.6). A straightforward computation shows that its symbol $p_\hbar^*$ satisfies $|p_\hbar(x, \xi) - p_\hbar^*(x, \xi)| \le C\hbar^{1/2}|\xi|^2$ uniformly in $x \in \mathbb{R}^d$, so that the second part of Proposition 3.3 also follows from Proposition A.6.
<table><thead><tr><th></th><th>BEM</th><th>Gaussian states</th><th>FEM</th></tr></thead><tbody><tr><td>Cost</td><td>$k^{d-1}$</td><td>$k^{d-1/2}$</td><td>$k^d$</td></tr><tr><td>Heterogeneous media</td><td>$\times$</td><td>✓</td><td>✓</td></tr><tr><td>Sparse matrix</td><td>$\times$</td><td>✓</td><td>✓</td></tr></tbody></table>

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

Galkowski and Spence, 2022

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

Melenk and Sauter, 2010

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

T. Chaumont-Frelet and M. Ingremeau, 2022

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

T. Chaumont-Frelet, V. Dolean and M. Ingremeau, 2022
3.2. Main proof. We focus here on the proof of Theorem 2.4. The notation

$$
\Lambda(a, b) := \{[\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d} \mid a \le |p(\mathbf{x}^{\mathbf{h}}, \mathbf{m}, \boldsymbol{\xi}^{\mathbf{h}, \mathbf{n}})| \le b\},
$$

where $0 \le a \le b \le +\infty$, will be useful. In what follows, we fix an $\varepsilon \in (0, 1/2)$, and consider a
right-hand side micro-localised near $\{p_\hbar = 0\}$. Specifically, we will assume that

$$
f_h := \sum_{[m,n] \in \Lambda_{h,rhs}} F_{m,n}^h \Psi_{h,m,n},
$$

where $F^h \in \ell^2(\Lambda_{\text{rhs}})$ with $\Lambda_{\text{rhs}} := \Lambda(0, \alpha + 2\hbar^{1/2-\epsilon})$. Our goal is to show that the associated solution $u_h$ is essentially micro-localised near $\{p_h = 0\}$ as well. Specifically, setting $\Lambda_{\text{h,near}} := \Lambda(0, \alpha + 4\hbar^{1/2-\epsilon})$, our goal will be to show that

$$
u_h^{\text{near}} := \sum_{[m,n] \in \Lambda_{h,\text{near}}} (u_h, \Psi_{h,m,n}^*) \Psi_{h,m,n}
$$

is “close” to $u_h$.

The key idea is to separate the set of indices of $u_h$ into $\Lambda_{h,\text{near}}$, $\Lambda_{h,\text{mid}} := \Lambda(\alpha + 4\hbar^{1/2-\epsilon}, \alpha + 6\hbar^{1/2-\epsilon})$ and $\Lambda_{h,\text{far}} := \Lambda(\alpha + 6\hbar^{1/2-\epsilon}, +\infty)$. We shall also need the “enlarged” sets

$$
\Lambda_{h, \text{near}}^* := \Lambda(0, \alpha + 5\hbar^{1/2-\epsilon}), \quad \Lambda_{h, \text{mid}}^* := \Lambda(\alpha + 3\hbar^{1/2-\epsilon}, \alpha + 7\hbar^{1/2-\epsilon}), \quad \Lambda_{h, \text{far}}^* := \Lambda(\alpha + 5\hbar^{1/2-\epsilon}, +\infty), \\
\text{for the test functions.}
$$

We first state some elementary properties of these sets of indices. We do not report the (straight-
forward) proofs for the sake of shortness.

**Lemma 3.4 (Index sets).** *Assume that *α* + 7*ħ*¹/²⁻*ε* ≤ *δ*₀, then we have*

$$
(3.5) \quad \begin{array}{@{}l@{\quad}l@{}}
\operatorname{dist}(\Lambda_{\bar{h}, \text{far}}, \Lambda_{\bar{h}, \text{near}}^*) & \ge C \bar{h}^{-\epsilon} \\
\operatorname{dist}(\Lambda_{\bar{h}, \text{rhs}}, \Lambda_{\bar{h}, \text{mid}}^*) & \ge C \bar{h}^{-\epsilon} \\
\operatorname{dist}(\Lambda_{\bar{h}, \text{near}}, \Lambda_{\bar{h}, \text{far}}^*) & \ge C \bar{h}^{-\epsilon}
\end{array}
\quad \text{and there exists } C > 0 \text{ such that}
$$

$$
(3.6) \qquad \Lambda_{\bar{h}, \text{rhs}}^* \cup \Lambda_{\bar{h}, \text{near}}^* \cup \Lambda_{\bar{h}, \text{mid}}^* \subset B(0, C\bar{h}^{-1/2}).
$$

In addition, if $\bar{h}$ is small enough, we have

$$
(3.7) \quad \left\{ [\mathbf{m}', \mathbf{n}'] \in \mathbb{Z}^{2d} \mid \exists[\mathbf{m}, \mathbf{n}] \in \Lambda_{\bar{h}, \text{mid}}; \; \|[\mathbf{m}, \mathbf{n}] - [\mathbf{m}', \mathbf{n}']| \le \bar{h}^{-\epsilon/2} \right\} \subset \Lambda_{\bar{h}, \text{mid}}^*.
$$

**Lemma 3.5 (Quasi orthogonality away from RHS micro-support).** For $F^h \in \ell^2(\Lambda_{h,rhs})$, consider the right-hand side

$$
f_h = \sum_{[m,n] \in \Lambda_{h,rhs}} F_{m,n}^h \Psi_{h,m,n}
$$

and the associated solution $u_h$. Then, if $[\mathbf{m}, \mathbf{n}] \in \Lambda_{h,\text{mid}}^*$, we have

$$
| (u_h, \Psi_{h,m,n}) | \leq C_{\epsilon,m} \hbar^m \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}
$$

for all $m \in \mathbb{N}$.

Proof. Throughout the proof, we fix a pair of indices $[\mathbf{m}, \mathbf{n}] \in \Lambda_{\hbar, \text{mid}}^*$. By definition of $\Lambda_{\hbar, \text{mid}}^*$, the assumption that $\alpha + 7\hbar^{1/2-\epsilon} \le \delta_0$ and (2.5), we have

$$
(3.8) \qquad c\hbar^{1/2-\epsilon} \le |p_h(\mathbf{m}, \mathbf{n})| \le C.
$$

In particular, we can write

$$
u_h = \frac{1}{p_h(\mathbf{m}, \mathbf{n})} f_h + \frac{1}{p_h(\mathbf{m}, \mathbf{n})} (P_h - p_h(\mathbf{m}, \mathbf{n})) u_h.
$$

Since $f_h$ is smooth, so is $u_h$, and we can iterate this relation, leading to

$$
u_h = \sum_{\ell=1}^{r} \frac{1}{p_h(\mathbf{m}, \mathbf{n})^\ell} (P_h - p_h(\mathbf{m}, \mathbf{n}))^{\ell-1} f_h + \frac{1}{p_h(\mathbf{m}, \mathbf{n})^r} (P_h - p_h(\mathbf{m}, \mathbf{n}))^r u_h,
$$
1 The model scattering problem

2 Convergence for general right-hand sides

3 Convergence for scattering problems

4 Numerical examples

5 Concluding remarks
and

$$ (u_h, \Psi_{h,m,n}) = \sum_{\ell=1}^{r} \frac{1}{\mathrm{p}_h(m,n)^\ell} (f_h, (P_h^* - \bar{p}_h(m,n))^{\ell-1} \Psi_{h,m,n}) \\ + \frac{1}{\mathrm{p}_h(m,n)^r} (u_h, (P_h^* - \bar{p}_h(m,n))^r \Psi_{h,m,n}) \tag{3.9} $$

for all $r \in \mathbb{N}$.

Then, if $[m', n'] \in \Lambda_{h,rhs}$, using the upper-bound in (3.8), we have

$$ |(\Psi_{h,m',n'}, (P_h^* - \bar{p}_h(m,n))^{\ell-1} \Psi_{h,m,n})| = \left| \sum_{j=0}^{\ell-1} \binom{\ell-1}{j} \bar{p}_h(m,n)^{\ell-1-j} (\Psi_{h,m',n'}, (P_h^*)^j \Psi_{h,m,n}) \right| \\ \le C_\ell \sum_{j=0}^{\ell-1} |(\Psi_{h,m',n'}, (P_h^*)^j \Psi_{h,m,n})|. $$

Recalling (3.6), applying (3.3) gives

$$ |(\Psi_{h,m',n'}, (P_h^* - \bar{p}_h(m,n))^{\ell-1} \Psi_{h,m,n})| \le C_{\ell,n} (1 + |[m,n] - [m',n']|)^{-n} \le C_{\ell,n} \hbar^{\epsilon n}, $$

for all $n \in \mathbb{N}$, since $|[m,n] - [m',n']| \ge \hbar^{-\epsilon}$ due to (3.5). The case $\ell = 1$ also easily follows by Proposition 3.2 We then write that

$$ |(f_h, (P_h^* - \bar{p}_h(m,n))^{\ell-1} \Psi_{h,m,n})| \le \sum_{[m',n'] \in \Lambda_{h,rhs}} |F_{m,n}^h|(\Psi_{h,m',n'}, (P_h^* - \bar{p}_h(m,n))^{\ell-1} \Psi_{h,m,n})| \\ \le C_{\ell,n} \hbar^{\epsilon n} \sum_{[m',n'] \in \Lambda_{h,rhs}} |F_{m,n}^h| \\ \le C_{\ell,n} \hbar^{\epsilon n - d/2} \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}, $$

where we used the Cauchy-Schwarz inequality and the fact that $|\Lambda_{h,rhs}| \le C\hbar^{-d}$ due to (3.6). As a result, using the lower-bound in (3.8), we have

$$ \frac{1}{|\mathrm{p}_h(m,n)|^\ell} |(f_h, (P_h^* - \bar{p}_h(m,n))^{\ell-1} \Psi_{h,m,n})| \le C_{\ell,n} \hbar^{\epsilon n - d/2 - \ell/2} \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}, $$

and

$$ \sum_{\ell=1}^{r} \frac{1}{|\mathrm{p}_h(m,n)|^\ell} |(f_h, (P_h^* - \bar{p}_h(m,n))^{\ell-1} \Psi_{h,m,n})| \le C_{m,n} \hbar^{\epsilon n - d/2 - r/2} \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}, $$

for all $n \in \mathbb{N}$. Thus, for any $m \in \mathbb{N}$, we can select $n = n(m, d, r, \epsilon)$ such that $\epsilon n - (d + r)/2 \ge m$, leading to

$$ \sum_{\ell=1}^{r} \frac{1}{|\mathrm{p}_h(m,n)|^\ell} |(f_h, (P_h^* - \bar{p}_h(m,n))^{\ell-1} \Psi_{h,m,n})| \le C_{\epsilon,m,r} \hbar^m \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}. $$

On the other hand, using again (3.6) and applying (3.4b), we have

$$ \|(P_h^* - \bar{p}_h(m,n))^r \Psi_{h,m,n}\|_{L^2(\mathbb{R}^d)} \le C_r \hbar^{r/2}, $$

and the lower bound in (3.8) shows that

$$ \frac{1}{|\mathrm{p}_h(m,n)|^r} \|(P_h^* - \bar{p}_h(m,n))^r \Psi_{h,m,n}\|_{L^2(\mathbb{R}^d)} \le C_r \hbar^{\epsilon r}. $$

We then write that

$$ \frac{1}{|\mathrm{p}_h(m,n)|^r} |(u_h, (P_h^* - \bar{p}_h(m,n))^r \Psi_{h,m,n})| \le C_r \hbar^{\epsilon r} \|u_h\|_{L^2(\mathbb{R}^d)} \le \\ C_r \hbar^{\epsilon r - N} \|f_h\|_{L^2(\mathbb{R}^d)} \le C_r \hbar^{\epsilon r - N} \|F^h\|_{\ell^2(\mathbb{Z}^{2d})} \le C_{\epsilon,m} \hbar^m \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}, $$

up to picking $r$ such that $\epsilon r - N \ge m$. $\square$

$\square$
The model scattering problem
Proof of Theorem 2.4. We expand $u_h$ in the frame $(\Psi_{h,m,n})_{[m,n]\in\mathbb{Z}^{2d}}$ as

$$u_h = u_h^{\text{near}} + u_h^{\text{mid}} + u_h^{\text{far}}$$

where

$$u_h^{\text{near}} := \sum_{[m,n] \in \Lambda_{h,\text{near}}} (u_h, \Psi_{h,m,n}^*) \Psi_{h,m,n}$$

$$u_h^{\text{mid}} := \sum_{[m,n] \in \Lambda_{h,\text{mid}}} (u_h, \Psi_{h,m,n}^*) \Psi_{h,m,n}$$

$$u_h^{\text{far}} := \sum_{[m,n] \in \Lambda_{h,\text{far}}} (u_h, \Psi_{h,m,n}^*) \Psi_{h,m,n}.$$

The proof then consists in showing that $u_h^{\text{mid}}$ and $u_h^{\text{far}}$ are small.

**Step 1.** We first treat the $u_h^{\text{mid}}$ term. To do so, we start by introducing the approximation

$$(3.10) \qquad \tilde{u}_h^{\text{mid}} := \sum_{[m,n] \in \Lambda_{h,\text{mid}}} \sum_{[m',n'] \in \mathbb{Z}^{2d} \atop |[m,n]-[m',n']| \le \hbar^{-\varepsilon/2}} U_{m',n'}^{m,n} (u_h, \Psi_{h,m',n'}) \Psi_{h,m,n}.$$

Recalling (3.7), all the $[m', n']$ indices in the sum belong to the enlarged set $\Lambda_{h,\text{mid}}^*$, so that

$$(3.11) \qquad |(u_h, \Psi_{h,m',n'})| \le C_{\varepsilon,m} \hbar^m \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}$$

by Lemma 3.5. Recalling (3.6), $|[m, n]| \le C\hbar^{-1/2}$ for all $[m, n] \in \Lambda_{h,\text{mid}}$, and we have from (3.1)

$$(3.12) \qquad \| \Psi_{h,m,n} \|_{\hat{H}_h^p(\mathbb{R}^d)} \le C_p.$$

Thus, plugging (3.11) and (3.12) into (3.10), we have

$$\|\tilde{u}_h^{\text{mid}}\|_{L^2(\mathbb{R}^d)} \le C_{p,\varepsilon,m} \hbar^m \|F^h\|_{\ell^2(\mathbb{Z}^{2d})} \sum_{[m,n] \in \Lambda_{h,\text{mid}}} \sum_{[m',n'] \in \mathbb{Z}^{2d} \atop |[m,n]-[m',n']| \le \hbar^{-\varepsilon/2}} 1$$

$$(3.13) \qquad \le C_{p,\varepsilon,m} \hbar^{m-\varepsilon d^2/2} \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}.$$

We now estimate the difference between $u_h^{\text{mid}}$ and $\tilde{u}_h^{\text{mid}}$

$$u_h^{\text{mid}} - \tilde{u}_h^{\text{mid}} = \sum_{[m,n] \in \Lambda_{h,\text{mid}}} \left( u_h, \Psi_{h,m,n}^* - \sum_{[m',n'] \in \mathbb{Z}^{2d} \atop |[m,n]-[m',n']| \le \hbar^{-\varepsilon/2}} U_{m',n'}^{m,n} \Psi_{h,m',n'} \right) \Psi_{h,m,n}$$

so that

$$\|u_h^{\text{mid}} - \tilde{u}_h^{\text{mid}}\|_{\hat{H}_h^p(\mathbb{R}^d)} \le C_p \|u_h\|_{L^2(\mathbb{R}^d)} \sum_{[m,n] \in \Lambda_{h,\text{mid}}} \left\| \Psi_{h,m,n}^* - \sum_{[m',n'] \in \mathbb{Z}^{2d} \atop |[m,n]-[m',n']| \le \hbar^{-\varepsilon/2}} \Psi_{h,m',n'} \right\|_{L^2(\mathbb{R}^d)}$$

and it follows from Proposition 3.1 that

$$(3.14) \qquad \|u_h^{\text{mid}} - \tilde{u}_h^{\text{mid}}\|_{\hat{H}_h^p(\mathbb{R}^d)} \le C_{p,\varepsilon,m} \hbar^{m-d} \|u_h\|_{\hat{H}_h^p(\mathbb{R}^d)} \le C_{p,\varepsilon,m} \hbar^{m-d-N} \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}.$$

Then, it follows from (3.13) and (3.14) that

$$\|u_h^{\text{mid}}\|_{\hat{H}_h^p(\mathbb{R}^d)} \le C_{p,\varepsilon,m} \hbar^m \|F^h\|_{\ell^2(\mathbb{Z}^{2d})} \quad \forall m \in \mathbb{N},$$

up to redefining $m$.
# The model scattering problem

## Standard formulation of the scattering problem
**Step 2.** We then turn our attention to $u_h^{\text{far}}$. On the one hand, we can apply Proposition 3.2 with $\mathcal{P}_{h,L,A} := \hbar^{[\beta]} \partial^{\beta} \circ P_h$, $\Lambda = \{[\mathbf{m}, \mathbf{n}]\} \subset \Lambda_{h,\text{near}}^*$ and $\Lambda' = \Lambda_{h,\text{far}}$. Using Lemma 3.4 and the Cauchy-Schwarz inequality, this gives

$$ \begin{aligned} |(\mathbf{x}^{\alpha} \partial^{\beta} (P_h u_h^{\text{far}}), \Psi_{h,m,n})| &= \left| \sum_{[m',n'] \in \Lambda_{h,\text{far}}} (u_h, \Psi_{h,m',n'}^{*})(\mathbf{x}^{\alpha} \mathcal{P}_{h,L,A} \Psi_{h,m',n'}, \Psi_{h,m,n}) \right| \\ &\le C_{\alpha,\beta,r} \|u_h\|_{L^2(\mathbb{R}^d)} \hbar^r \quad \forall r \in \mathbb{N}. \end{aligned} $$

As a result, since $|\Lambda_{h,\text{near}}^*| \le C\hbar^{-d}$ due to (3.6), we have

$$ \sum_{[m,n] \in \Lambda_{h,\text{near}}^*} |(\mathbf{x}^{\alpha} \partial^{\beta} (P_h u_h^{\text{far}}), \Psi_{h,m,n})|^2 \le C_{\alpha,\beta,r} \hbar^{2r-d} \|u_h\|_{L^2(\mathbb{R}^d)}^2 $$

which we rewrite as

$$ \sum_{[m,n] \in \Lambda_{h,\text{near}}^*} |(\mathbf{x}^{\alpha} \partial^{\beta} (P_h u_h^{\text{far}}), \Psi_{h,m,n})|^2 \le C_{\alpha,\beta,m} \hbar^{2m} \|u_h\|_{L^2(\mathbb{R}^d)}^2 \quad \forall m \in \mathbb{N} $$

after changing variables.

On the other hand, if $[\mathbf{m}, \mathbf{n}] \in \Lambda_{h,\text{far}}^*$, we write that

$$ P_h u_h^{\text{far}} = f_h - P_h u_h^{\text{near}} - P_h u_h^{\text{mid}}, $$

so that

$$ |(\mathbf{x}^{\alpha} \partial^{\beta} (P_h u_h^{\text{far}}), \Psi_{h,m,n})|^2 \le C \left( |(\mathbf{x}^{\alpha} \partial^{\beta} f_h, \Psi_{h,m,n})|^2 + |(\mathbf{x}^{\alpha} \partial^{\beta} (P_h u_h^{\text{near}}), \Psi_{h,m,n})|^2 + |(\mathbf{x}^{\alpha} \partial^{\beta} (P_h u_h^{\text{mid}}), \Psi_{h,m,n})|^2 \right). $$

We then have

$$ \sum_{[m,n] \in \Lambda_{h,\text{far}}^*} |(\mathbf{x}^{\alpha} \partial^{\beta} (P_h u_h^{\text{mid}}), \Psi_{h,m,n})|^2 \le C \|\mathbf{x}^{\alpha} \partial^{\beta} (P_h u_h^{\text{mid}})\|_{L^2(\mathbb{R}^d)}^2 \le C \|u_h^{\text{mid}}\|_{\hat{H}_h^{[\beta]+2}(\mathbb{R}^d)}^2 \le C_{\beta,m} \hbar^m \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}^2 $$

due to **Step 1**. Next, thanks to Proposition 3.2,

$$ \begin{aligned} \sum_{[m,n] \in \Lambda_{h,\text{far}}^*} |(\mathbf{x}^{\alpha} \partial^{\beta} f_h, \Psi_{h,m,n})|^2 &\le \sum_{[m,n] \in \Lambda_{h,\text{far}}^*} \sum_{[m',n'] \in \Lambda_{h,\text{rms}}} |F_{m,n}^{\hbar}|^2 |(\mathbf{x}^{\alpha} \partial^{\beta} \Psi_{h,m',n'}, \Psi_{h,m,n})|^2 \\ &\le C_{\alpha,\beta,m}^2 \hbar^{2m} \|F^{\hbar}\|_{\ell^2(\mathbb{Z}^{2d})}^2. \end{aligned} $$

Finally, again by Proposition 3.2, we have

$$ \begin{aligned} \sum_{[m,n] \in \Lambda_{h,\text{far}}^*} |(\mathbf{x}^{\alpha} \partial^{\beta} (P_h u_h^{\text{near}}), \Psi_{h,m,n})|^2 &\le \sum_{[m,n] \in \Lambda_{h,\text{far}}^*} \sum_{[m',n'] \in \Lambda_{h,\text{near}}} |(u, \Psi_{k,m',n'}^*)(\mathbf{x}^{\alpha} \partial^{\beta} (P_h \Psi_{h,m',n'}), \Psi_{h,m,n})|^2 \\ &\le C_{\alpha,\beta,r} \hbar^{2r} \|u_h\|_{L^2(\mathbb{R}^d)}^2 \le C_{\alpha,\beta,m} \hbar^{2m} \|F^h\|_{\ell^2(\mathbb{Z}^{2d})}^2. \end{aligned} $$

Since $\Lambda_{h,\text{near}}^*$ and $\Lambda_{h,\text{far}}^*$ form a partition of $\mathbb{Z}^{2d}$, we have thus shown that

$$ (3.15) \qquad \|\mathbf{x}^{\alpha} \partial^{\beta} (P_h u_h^{\text{far}})\|_{L^2(\mathbb{R}^d)} \le C_{\alpha,\beta,m} \hbar^m \|F^h\|_{\ell^2(\mathbb{Z}^{2d})} $$

Letting $f_h^{\text{far}} = P_h u_h^{\text{far}}$, we see from (3.15) that $u_h^{\text{far}}$ solves $P_h u_h^{\text{far}} = f_h^{\text{far}}$ with a right-hand side $f_h^{\text{far}} \in \hat{H}_h^p(\mathbb{R}^d)$ such that

$$ \|f_h^{\text{far}}\|_{\hat{H}_h^p(\mathbb{R}^d)} \le C_{p,m} \hbar^m \|F^h\|_{\ell^2(\mathbb{Z}^{2d})} \quad \forall m \in \mathbb{N}. $$

Then, we conclude the proof with (2.4).

## 4. APPLICATION TO THE HELMHOLTZ EQUATION

We now turn our attention to the model problem of the Helmholtz equation (1.1), with a particular focus on the case of plane-wave scattering.
The scattering problem by a planewave is usually formulated that way:

<figure><img src="image_7.png" /></figure>
**4.1. Notation.** In this section $k$ will denote the wavenumber in the Helmholtz problem, For the sake of simplicity, we assume that $kR \ge 1$. We will apply the results of Section 2 with $\hbar \sim (kR)^{-1}$. As a result, the norms

$$ (4.1) \qquad \|v\|_{\hat{H}_k^p(\mathbb{R}^d)}^2 := \sum_{\substack{\alpha \in \mathbb{N}^d \\ [\alpha] \le p}} \sum_{q \le p - [\alpha]} k^{-2[\alpha]} \left\| \frac{\mathbf{x}}{R} \right\|^q \partial^\alpha v \Big|_{L^2(\mathbb{R}^d)}^2 $$

will be convenient. Notice that if

$$ \mathcal{F}(v)(\xi) := (2\pi)^{-d/2} \int_{\mathbb{R}^d} v(\mathbf{x})e^{-i\mathbf{x} \cdot \xi} d\mathbf{x}, \text{ a.e. } \xi \in \mathbb{R}^d. $$

is the Fourier transform of $v \in L^2(\mathbb{R}^d)$, and if we define the “reverse” norm by

$$ (4.2) \qquad \|v\|_{\hat{H}_k^p(\mathbb{R}^d)}^2 := \sum_{\substack{\alpha \in \mathbb{N}^d \\ [\alpha] \le p}} \sum_{q \le p - [\alpha]} k^{-2q} \left\| \frac{\mathbf{x}}{R} \right\|^q \partial^\alpha v \Big|_{L^2(\mathbb{R}^d)}^2, $$

then we have

$$ (4.3) \qquad c(R) \|\mathcal{F}v\|_{\hat{H}_k^p(\mathbb{R}^d)} \le \|v\|_{\hat{H}_k^p(\mathbb{R}^d)} \le C(R) \|\mathcal{F}v\|_{\hat{H}_k^p(\mathbb{R}^d)}. $$

We will also use the following “standard” norm

$$ \|v\|_{\hat{H}_k^p(\mathbb{R}^d)}^2 := \sum_{\substack{\alpha \in \mathbb{N}^d \\ [\alpha] \le p}} k^{-2[\alpha]} \|\partial^\alpha v\|_{L^2(\mathbb{R}^d)}^2. $$

**4.2. Model problem.** We consider smooth coefficients $\mu, A \in C_b^\infty(\mathbb{R}^d)$ that are respectively equal to 1 and $I$ outside $B(0, R)$. Given $f: \Omega \to \mathbb{C}$ our model problem is to find $u: \Omega \to \mathbb{C}$ such that

$$ (4.4a) \qquad -k^2 \mu u - \nabla \cdot (A \nabla u) = f \text{ in } \mathbb{R}^d $$

and

$$ (4.4b) \qquad \frac{\partial u}{\partial |\mathbf{x}|}(\mathbf{x}) - iku(\mathbf{x}) = o\left(|\mathbf{x}|^{-(d-1)/2}\right) \text{ as } |\mathbf{x}| \to +\infty. $$

Problem (4.4) is well-posed in the sense that for all $f \in L_{\text{comp}}^2(\mathbb{R}^d)$, there exists a unique $R_k f := u \in L_{\text{loc}}^2(\mathbb{R}^d)$ such that (4.4a) and (4.4b) hold true.

Following the assumption on $\hbar$ in Section 2.1 and assumption (2.4), we focus on wave numbers such that $k \in \mathcal{K} \subset [1/R, +\infty)$, and assume that there exists $N \ge 1$ such that

$$ (4.5) \qquad k^2 \|\chi R_k \chi\|_{L^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)} \le C(kR)^N \quad \forall k \in \mathcal{K} $$

where $\chi$ is a smooth cut-off function that takes the value 1 on $B(0, R)$ and 0 outside $B(0, 2R)$.

Notice that in view of the term $-k^2\mu u$, in (4.4a), the factor $k^2$ in (4.5) ensures that the right-hand side does not bear any physical dimension. We have explicitly required that $N \ge 1$ for the sake of clarity, but in fact, it can be shown that $N = 1$ is the best possible exponent for which (4.5) can hold true, see e.g. [4, Lemma 4.1].

**Remark 4.1** (When does the polynomial bound actually hold?). The bound (4.5) is known to hold in several situations:

* When the dynamics induced by the Hamiltonian $p$ has no trapped trajectory, i.e., when every trajectory leaves any compact set in finite time, the assumption holds with $\mathcal{K} = [1/R, +\infty)$. See for instance [17]. This situation is often referred to as “non-trapping”.

* When the dynamics induced by $p$ has a trapped set, and the dynamics is hyperbolic, close to this trapped set, it has been conjectured in [44] that (2.4) always holds with $\mathcal{K} = [1/R, +\infty)$. Actually, this is already known when the trapped set is “filamentary enough”, see [37, 38].
The scattering problem by a planewave is usually formulated that way:

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

$\mu$ and $A$ are smooth coefficients, with $\mu(x) = 1$ and $A(x) = I$ if $|x| \ge R$.
* Without any assumption on the dynamics, (2.4) holds taking $\mathcal{K}$ to be $[1/R, +\infty)$ from which we exclude a set of frequencies $k$ whose intersection with $\{k \ge k_0\}$ has a length going to zero as $k_0 \to +\infty$. We refer the reader to [29] for more details.

### 4.3. Perfectly matched layers.
As advertised in the introduction, the formulation (4.4) is not suited for immediate discretization by “volume” methods, as the radiation condition is hard to take into account. We will thus rely on an equivalent formulation that uses perfectly matched layers (PML).

Specifically, given $f: \Omega \to \mathbb{C}$, we consider the problem to find $u: \Omega \to \mathbb{C}$ such that $P_k u = f$ where

$$ (4.6) \qquad P_k u := -\frac{1}{k^2} ((\boldsymbol{I} + i\boldsymbol{M})^{-1}\boldsymbol{\nabla}) \cdot ((\boldsymbol{I} + i\boldsymbol{M})^{-1}\boldsymbol{A}\boldsymbol{\nabla}u) - \mu u. $$

In (4.6), the (SPD) matrix function $\boldsymbol{M}$ is given by

$$ \boldsymbol{M}(\boldsymbol{x}) := \frac{g(|\boldsymbol{x}|)}{|\boldsymbol{x}|^3} (|\boldsymbol{x}|^2 \mathrm{Id} - \boldsymbol{x} \otimes \boldsymbol{x}) + \frac{g'(|\boldsymbol{x}|)}{|\boldsymbol{x}|^2} \boldsymbol{x} \otimes \boldsymbol{x}, $$

where $g: \mathbb{R} \to \mathbb{R}$ is a user-defined function such that $g(r) = 0$ if $r \le R$, $g(r) = r$ if $r \ge R_0 > R$, and $g'(r) \ge 0$. In what follows, we will assume that $g$ is a smooth function to satisfy the assumptions of Section 2, but many results about PML still hold with less regular $g$ (see, e.g., [15]).

Notice that, if $|\boldsymbol{x}| \le R$, $\boldsymbol{M} = \boldsymbol{0}$, so that the original operator is not modified on the support of $\mu - 1$ and $\boldsymbol{A} - \boldsymbol{I}$. On the other hand, $\boldsymbol{M} = \boldsymbol{I}$ if $|\boldsymbol{x}| \ge R_0$, so that dissipation is introduced away from the origin, where the operator simply reads

$$ P_k v = \frac{i}{2k^2} \Delta v - v \quad \text{whenever} \quad \mathrm{supp} v \cap B(0, 2R) = \emptyset. $$

This transformation can be naturally interpreted as a complex deformation of coordinates. It is also often called the complex scaling technique. Crucially, the PML is designed in such way that

$$ (4.7) \qquad ((P_k)^{-1} f)|_{B(0,R)} = (R_k f)|_{B(0,R)} $$

whenever $\mathrm{supp} f \subset B(0, R)$. We refer the reader to [10, §4.5] or [15] for more information.

### 4.4. Abstract setting.
We now verify that the Helmholtz problem formulated with PML indeed fits the abstract setting of Section 2.3. The only non-trivial facts to establish are the polynomial resolvent estimates in $\hat{H}^q(\mathbb{R}^d)$ in (2.4) and the boundedness of the energy layer (2.5).

**Lemma 4.2 (Resolvent estimates).** Let $q \in \mathbb{N}$. For all $f \in \hat{H}^q(\mathbb{R}^d)$, there exists a unique $u \in L^2(\mathbb{R}^d)$ such that $P_k u = f$. In addition, we have $u \in \hat{H}^q(\mathbb{R}^d)$ with

$$ k^2 \|u\|_{\hat{H}_k^q(\mathbb{R}^d)} \le C(kR)^N \|f\|_{\hat{H}_k^q(\mathbb{R}^d)}. $$

Furthermore, if $\mathrm{supp} f \subset B(0, R)$, then $u \in \hat{H}^{q+2}(\mathbb{R}^d)$ with

$$ (4.8) \qquad k^2 \|u\|_{\hat{H}_k^{q+2}(\mathbb{R}^d)} \le C(kR)^N \|f\|_{\hat{H}_k^q(\mathbb{R}^d)}. $$

**Proof.** We first invoke Theorem 1.6 of [15] which states that

$$ \|P_k^{-1}\|_{L^2(\mathbb{R}^d) \to H_k^2(\mathbb{R}^d)} \le C \|\chi R_k \chi\|_{L^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)}. $$

Then, using a usual bootstrap technique, we easily show that

$$ (4.9) \qquad \|P_k^{-1}\|_{H_k^q(\mathbb{R}^d) \to H_k^{q+2}(\mathbb{R}^d)} \le C \|\chi R_k \chi\|_{L^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)} $$

for all $q \in \mathbb{N}$.

We then need to take care of the weights in the $\hat{H}^q(\mathbb{R}^d)$ norms. To do so, we observe that if $P_k u = f$, then we may write

$$ 2ik^2 u + \Delta u = g, $$

with

$$ g := -2ik^2 P_k u + (2ik^2 P_k u + \Delta u + 2ik^2)u = -2ik^2(f + Q_k u), $$
The scattering problem by a planewave is usually formulated that way:

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

$\mu$ and $A$ are smooth coefficients, with $\mu(x) = 1$ and $A(x) = I$ if $|x| \ge R$.

Problem:

we do not have $u^{\text{sca}} \in H^1(\mathbb{R}^d)$, and the right-hand side is not convenient.
where $Q_k$ is differential operator of degree 2 with smooth coefficients supported in $B(0, R_0)$. Let $\alpha, \beta \in \mathbb{N}^d$ with $|\alpha|, |\beta| \le 2$. Since $Q_k$ is compactly supported, we have

$$ (4.10) \qquad \left\| \left( \frac{\mathbf{x}}{R} \right)^{\alpha} g \right\|_{\hat{H}_k^{q}(\mathbb{R}^d)} \le C_{\alpha} k^2 \left( \|u\|_{H_k^{q+2}(\mathbb{R}^d)} + \|f\|_{\hat{H}_k^{q+2}(\mathbb{R}^d)} \right). $$

Now, we note that

$$ (4.11) \qquad \left( \frac{\mathbf{x}^{\alpha}}{R} \right) \partial^{\beta} u = i^{\alpha-\beta} \mathcal{F}^{-1} \left( \left( \frac{1}{R} \partial \right)^{\alpha} \left( \frac{\xi^{\beta}}{-|\xi|^2 + 2ik^2} \mathcal{F}(g) \right) \right), $$

Let us write

$$ \frac{\xi^{\beta}}{-|\xi|^2 + 2ik^2} = k^{[\beta]-2} \frac{(\xi/k)^{\beta}}{-|\xi/k|^2 + 2i}, $$

so that the map

$$ \xi \mapsto \frac{\xi^{\beta}}{-|\xi|^2 + 2ik^2} $$

has $C^\ell$ norm bounded by $C_\ell k^{[\beta]-2-\ell} \le C_\ell k^{[\beta]-2}$, since $k \ge 1$.

We deduce from (4.3) and (4.10) that

$$ \left\| \left( \frac{1}{R} \partial \right)^{\alpha} \left( \frac{\xi^{\beta}}{-|\xi|^2 + 2ik^2} \mathcal{F}(g) \right) \right\|_{\hat{H}_k^{q}(\mathbb{R}^d)} \le C k^{[\beta]} \left( \|u\|_{H_k^{q+2}(\mathbb{R}^d)} + \|f\|_{\hat{H}_k^{q+2}(\mathbb{R}^d)} \right) $$

and hence, thanks to (4.11),

$$ \begin{aligned} \|u\|_{\hat{H}_k^{q+2}(\mathbb{R}^d)} &\le C \sum_{[\alpha]+[\beta] \le 2} k^{-[\beta]} \left\| \left(\frac{\mathbf{x}}{R}\right)^{\alpha} \partial^{\beta} u \right\|_{\hat{H}_k^p(\mathbb{R}^d)} \\ &\le C \sum_{[\alpha]+[\beta] \le 2} k^{-[\beta]} \left\| \mathcal{F} \left( \left(\frac{\mathbf{x}}{R}\right)^{\alpha} \partial^{\beta} u \right) \right\|_{\hat{H}_k^p(\mathbb{R}^d)} \\ &\le C \sum_{[\alpha]+[\beta] \le 2} \left( \|u\|_{H_k^{q+2}(\mathbb{R}^d)} + \|f\|_{\hat{H}_k^{q+2}(\mathbb{R}^d)} \right) \\ &\le C \|\chi R_k \chi\|_{L^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)} \|f\|_{\hat{H}_k^{q+2}(\mathbb{R}^d)}, \end{aligned} $$

thanks to (4.9). The first part of the result follows from (4.5).

Now, if we further assume that $f$ is supported in $B(0, R)$, $g$ is compactly supported and equation (4.10) may be replaced by

$$ \left\| \left( \frac{\mathbf{x}}{R} \right)^{\alpha} g \right\|_{\hat{H}_k^q(\mathbb{R}^d)} \le C_{\alpha} k^2 \left( \|f\|_{H_k^q(\mathbb{R}^d)} + \|u\|_{H_k^{q+2}(\mathbb{R}^d)} \right), $$

and the same reasoning as above leads to

$$ \|u\|_{\hat{H}_k^{q+2}(\mathbb{R}^d)} \le C \left( \|u\|_{H_k^{q+2}(\mathbb{R}^d)} + \|f\|_{H_k^q(\mathbb{R}^d)} \right) \le C \|\chi R_k \chi\|_{L^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)} \|f\|_{\hat{H}_k^q(\mathbb{R}^d)}, $$

as announced.

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

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

**Lemma 4.3** (Boundedness of the energy layer). For all $\delta \in (0, 1)$, the set

$$ \{[\mathbf{x}, \boldsymbol{\xi}] \in \mathbb{R}^{2d} \mid |p(\mathbf{x}, \boldsymbol{\xi})| \le \delta\} $$

is bounded.

*Proof.* Fix $0 < \delta < 1$, and let $U := \{[\mathbf{x}, \boldsymbol{\xi}] \in \mathbb{R}^{2d}; |p(\mathbf{x}, \boldsymbol{\xi})| \le \delta\}$. Let $\mathbf{x} \in \mathbb{R}^d$. We first assume that $|\mathbf{x}| \le R_0$. Since we know that

$$ C|\boldsymbol{\xi}|^2 - C' \le |p(\mathbf{x}, \boldsymbol{\xi})| \le \delta, $$

we see that $\{\boldsymbol{\xi} \in \mathbb{R}^d | p(\mathbf{x}, \boldsymbol{\xi}) \le \delta\}$ is bounded, and thus $U \cap (B_\mathbf{x}(0, R_0) \times \mathbb{R}^d)$ is bounded. On the other hand, if $|\mathbf{x}| \ge R_0$, we have

$$ p(\mathbf{x}, \boldsymbol{\xi}) = (1 + i)^2 |\boldsymbol{\xi}|^2 - 1 = 2i|\boldsymbol{\xi}|^2 - 1 $$

so that

$$ |p(\mathbf{x}, \boldsymbol{\xi})|^2 = 4|\boldsymbol{\xi}|^4 + 1 \ge 1 > \delta, $$

which implies that $U \setminus (B_\mathbf{x}(0, R_0) \times \mathbb{R}^d) = \emptyset$.

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

<figure><img src="image_6.png" /></figure>
The model scattering problem

Volume right-hand side
We finally show that the right-hand side associated with plane-wave scattering are indeed well-approximated by Gaussian coherent states in order to apply Corollary 2.5 later on.

**Lemma 4.4** (Approximability of plane-wave right-hand sides). For $k > 0$, consider the right-hand side $f_k := (-k^2 - \Delta)(\chi e^{ikd \cdot x})$ where $\chi \in C_c^\infty(B_R)$ and $d \in \mathbb{R}^d$ with $|d| = 1$. Then, there exists $F^k \in \ell^2(\mathbb{Z}^{2d})$ such that

$$ (4.12) \qquad R^2 \left\| f_k - \sum_{[m,n] \in \Lambda_{k,\text{rhs}}} F_{m,n}^k \Psi_{k,m,n} \right\|_{\hat{H}_k^p(\mathbb{R}^d)} \le C_{\varepsilon,p,m} (kR)^{-m}, $$

for all $m \in \mathbb{N}$, where

$$ \Lambda'_{k,\text{rhs}} := \{[m,n] \in \mathbb{Z}^{2d} \mid \text{dist}(x^{k,m}, \text{supp}(\chi)) \le (kR)^{-1/2+\varepsilon} \text{ and } |\xi^{k,n}| - 1| \le (kR)^{-1/2+\varepsilon}\}. $$

In particular, if $p$ is the symbol of the function appearing in section 4.3 and if $\chi$ is supported in the region where $p(x, \xi) = |\xi|^2 - 1$, then $\Lambda'_{k,\text{rhs}}$ is of the same form as $\Lambda_{k,\text{rhs}}$ introduced in (2.6) with $\alpha = 0$.

*Proof.* We start by observing that

$$ f_k = -(\Delta\chi + 2ikd \cdot \nabla\chi)e^{ikd \cdot x} = -R^{-2}(R^2\Delta\chi + 2ikR(Rd \cdot \nabla\chi))e^{ikd \cdot x}, $$

so that the result holds true if we can show (4.12) for $\tilde{f}_k := \tilde{\chi}e^{ikd \cdot x}$ with $\tilde{\chi}$ smooth and compactly supported (without the factor $R^2$ in front of the left-hand side).

From now on, we fix $\tilde{\chi}$, $p$, $m$ and $\varepsilon$, and consider $\tilde{f}_k$ as above. We will first show that, if $[m, n] \notin \Lambda'_{k, \text{rhs}}$, then we have

$$ (4.13) \qquad \left| (\tilde{f}_k, \Psi_{k,m,n}) \right| \le C_m (kR)^{-m} \quad \forall m \in \mathbb{N}. $$

The quantity $(\tilde{f}_k, \Psi_{k,m,n})$ is of the form

$$ \int_{\mathbb{R}^d} \tilde{\chi}(x) e^{ik\varphi_{m,n}(x)} dx, $$

with

$$ \varphi_{m,n}(x) = x \cdot (\xi - \xi^{k,n}) + \frac{i}{2} |x - x^{k,m}|^2, $$

so that

$$ \nabla \varphi_{m,n}(x) = \xi - \xi^{k,n} + i(x - x^{m,k}). $$

In particular, we have

$$ |\nabla \varphi_{m,n}(x)| \ge |\xi^{k,n}| - 1 + \text{dist}(x^{k,m}, \text{supp}(\chi)) \ge \frac{1}{2}(kR)^{-\frac{1}{2}+\varepsilon}. $$

we may use the method of non-stationary phase (i.e., integrate by parts several times, as in [43, Lemma 3.14]) to deduce (4.13).

Now, combining (4.13) with Proposition 3.1, we see that there exists $C$ such that, for any $[m, n] \in \Lambda_{k, \text{rhs}}$, we have

$$ (4.14) \qquad \left| (\tilde{f}_k, \Psi_{k,m,n}^*) \right| \le C_m (kR)^{-m} \quad \forall m \in \mathbb{N}. $$

On the other hand, it follows from [5, Theorem 3.1] that we have

$$ (4.15) \qquad \left\| \tilde{f}_k - \sum_{\substack{[m,n] \in \mathbb{Z}^{2d} \\ |m,n| \le k}} (f_k, \Psi_{k,m,n}^*) \Psi_{k,m,n} \right\|_{\hat{H}_k^p(\mathbb{R}^d)} \le C_{p,m} (kR)^{-m}. $$
Fix a smooth cutoff $\chi$ s.t. $\chi(x) = 1$ if $|x| \le R$, and $\chi(x) = 0$ if $|x| \ge R'$. Instead of solving for $u^{\text{tot}}$, we look for $u = \chi u^{\text{tot}} + (1-\chi)u^{\text{sca}}$. Then

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

$$
\tilde{f}_k = \sum_{[m,n] \in \Lambda_{k,\text{rhs}}} (\tilde{f}_k, \Psi_{k,m,n}^*) \Psi_{k,m,n} + \sum_{\substack{[m,n] \in \mathbb{Z}^{2d} \setminus \Lambda_{k,\text{rhs}} \\ |m,n| \le k}} (\tilde{f}_k, \Psi_{k,m,n}^*) \Psi_{k,m,n} \\
+ \left( \tilde{f}_k - \sum_{\substack{[m,n] \in \mathbb{Z}^{2d} \\ |m,n| \le k}} (f_k, \Psi_{k,m,n}^*) \Psi_{k,m,n} \right),
$$

we deduce from (4.14) and (4.15) that the last two terms have a $\hat{H}_k^p$ norm bounded by $C_{\epsilon,p,m}(kR)^{-m}$,
and the result follows.
$\square$

4.5. **Approximability estimates.** We are now ready to present our approximability estimates for the Helmholtz problem. We start with an approximation result for general right-hand sides that does not hinge on Section 2, but rather on standard results on Gabor frames and modulation spaces, see e.g. [22, Chapter 11]. We also refer the reader to the research note [5] where elementary proofs are presented.

Recall that $N \ge 1$ was introduced in (4.5) and that $kR \ge 1$.

**Theorem 4.5** (Approximability at a fixed frequency). Assume $f \in \hat{H}^p(\mathbb{R}^d)$ with supp $f \subset B(0, R)$. If

$$
\Lambda := \{[m, n] \in \mathbb{Z}^{2d} \mid \|[m, n]| \le \sqrt{\rho k R}\},
$$

for some $\rho > 0$, then we have

$$
k^2 \left\| u - \sum_{[m,n] \in \Lambda} (u, \Psi_{k,m,n}^*) \Psi_{k,m,n} \right\|_{\hat{H}_k^{q}(\mathbb{R}^d)} \le C_p (kR)^{N-(2+p-q)} \rho^{-(2+p-q)/2} \|f\|_{\hat{H}_k^{p}(\mathbb{R}^d)}
$$

for all $p, q \in \mathbb{N}$ with $q \le p$.

Proof. Let us set $D := (kR/\pi)^{1/2}(\rho kR)^{1/2} = \sqrt{\rho/\pi}(kR)$. We start with [5, Theorem 3.1], showing
that

$$
\left\| u - \sum_{[m,n] \in \Lambda} (u, \Psi_{k,m,n}^*) \Psi_{k,m,n} \right\|_{\hat{H}_k^q(\mathbb{R}^d)} \le C_p D^{q-p-2} \|u\|_{\hat{H}_k^{p+2}(\mathbb{R}^d)}.
$$

The result then follows using (4.8), since

$$
\begin{align*}
D^{q-p-2} k^2 \|u\|_{\hat{H}_k^{p+2}(\mathbb{R}^d)} &\le C_p (\sqrt{\rho k R})^{q-p-2} k^2 \|u\|_{\hat{H}_k^{p+2}(\mathbb{R}^d)} \\
&\le C_p (\sqrt{\rho k R})^{q-p-2} (kR)^N \|f\|_{\hat{H}_k^p(\mathbb{R}^d)} \\
&\le C_p (kR)^{N+q-p-2} \rho^{(q-p-2)/2} \|f\|_{\hat{H}_k^p(\mathbb{R}^d)}. \quad \square
\end{align*}
$$

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

Our second approximability estimate applies specifically to high-frequency scattering problems.
It is a direct consequence of Theorem 2.4 and Lemma 4.4.

**Theorem 4.6** (Approximability in the high-frequency regime). Fix $0 < \varepsilon < 1/2$ and consider the index set

$$
\Lambda := \{[m, n] \in \mathbb{Z}^{2d} \mid |p(m, n)| \le (kR)^{-1/2+\varepsilon}\}.
$$

Then, if the right-hand side is of the form $f_k := (-k^2 - \Delta)(\chi e^{ikd \cdot x})$ where $d \in \mathbb{R}^d$ with $|d| = 1$ and $\chi \in C_c^\infty(B_R)$ is supported in the region where $\mu = 1$ and $A = I$, for all $q \in \mathbb{N}$, we have

$$
\left\| u_k - \sum_{[m,n] \in \Lambda} (u_k, \Psi_{k,m,n}^*) \Psi_{k,m,n} \right\|_{\hat{H}_k^q(\mathbb{R}^d)} \le C_{\epsilon,q,m} (kR)^{-m} \quad \forall m \in \mathbb{N}.
$$
Fix a smooth cutoff $\chi$ s.t. $\chi(x) = 1$ if $|x| \le R$, and $\chi(x) = 0$ if $|x| \ge R'$. Instead of solving for $u^{\text{tot}}$, we look for $u = \chi u^{\text{tot}} + (1-\chi)u^{\text{sca}}$. Then

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

where the right-hand side has the specific form

$$f := (2ikd \cdot \nabla\chi + \Delta\chi)e^{ikd \cdot x} = \phi e^{ikd \cdot x}$$

with $\phi$ smooth and compactly supported.
4.6. **A least-squares method.** In this section, we introduce a least-squares method based on Gaussian coherent states. For a finite set $\Lambda \subset \mathbb{Z}^{2d}$, we consider the discretization space

$$W_{\Lambda} := \mathrm{Vect}\{\Psi_{k,m,n}; [m, n] \in \Lambda\}.$$

Then, the least squares method consists in finding $u_{\Lambda} \in W_{\Lambda}$ such that

$$(4.16) \qquad (P_k u_{\Lambda}, P_k w_{\Lambda}) = (f, P_k w_{\Lambda})$$

for all $w_{\Lambda} \in W_{\Lambda}$. Simple manipulations reveal that (4.16) is the Euler-Lagrange equation corresponding to the minimization problem

$$(4.17) \qquad \|P_k(u - u_{\Lambda})\|_{L^2(\mathbb{R}^d)} = \min_{w_{\Lambda} \in W_{\Lambda}} \|P_k(u - w_{\Lambda})\|_{L^2(\mathbb{R}^d)}.$$

Since $W_{\Lambda}$ is finite dimensional, the cost functional is quadratic, and $P_k$ is an isomorphism, $u_{\Lambda}$ is uniquely defined.

Using Theorem 4.5 we can easily show that the method converges for any fixed frequency, if sufficiently many DOFs per wavelength are employed, as stated in the following corollary:

**Corollary 4.7 (Convergence at a fixed frequency).** Consider the set

$$\Lambda_{\rho} := \{[m, n] \in \mathbb{Z}^{2d} \mid |m|^2 + |n|^2 \le \rho k R\},$$

then the following error estimates holds true:

$$k^2 \|u - u_{\Lambda}\|_{L^2(\mathbb{R}^d)} \le C_p (kR)^{N-p} \rho^{-p/2} \|f\|_{\hat{H}^p(\mathbb{R}^d)} \quad \forall p \in \mathbb{N}.$$

Relying on Theorem 4.6 instead, we can show a refined error estimate in the high-frequency regime.

**Corollary 4.8 (Convergence in the high-frequency regime).** Fix $0 < \varepsilon < 1/2$,

$$\Lambda := \{[m, n] \in \mathbb{Z}^{2d} \mid |p(\mathbf{x}^{k,m}, \xi^{k,n})| \le (kR)^{-1/2+\varepsilon}\},$$

and let $f_k := (-k^2 - \Delta)(\chi e^{ikd \cdot x})$ with $\chi \in C_c^{\infty}(B_R)$ supported in the region where $\mu = 1$ and $\mathbf{A} = \mathbf{I}$. Then, if $kR$ is sufficiently large, we have

$$k^2 \|u - u_{\Lambda}\|_{L^2(\mathbb{R}^d)} \le C_{\varepsilon,m} (kR)^{-m} \quad \forall m \in \mathbb{N}.$$

## 5. NUMERICAL RESULTS

In this section, we provide numerical illustrations of the above theory in the one-dimensional case. The purpose of these examples is simply to illustrate and support our theoretical findings. Extension to higher dimensions and discussion about efficient linear system assembly and solve will be reported elsewhere.

### 5.1. Setting.
We consider the one dimensional case, where the differential operator in (4.6) simplifies to

$$(5.1) \qquad P_k v = -\mu \nu u - k^{-2}(\alpha \nu^{-1} v')'$$

where $\mu, \alpha > 0$ are smooth “physical” coefficients, and $\nu$ is the PML scaling. For the sake of simplicity we take $R := 1$, meaning that $\mu = \alpha = 1$ outside of $(-1, 1)$ and that our right-hand sides $f$ will be supported in $(-1, 1)$. Here, we select $\nu := 1 + i\sigma$ where $\sigma$ is a stretching function defined as

$$\sigma(x) := a(\mathbf{1}_{(1, +\infty)}(x-1)^r + \mathbf{1}_{(-\infty, -1)}(-1-x)^r)$$

with $a := 1/10$ and $r := 4$. Notice that this choice slightly departs from our theoretical framework as the coefficients are only $C^3$ and not $C^{\infty}$.

In all our experiments, the grid in phase space is chosen as

$$x^{k,m} = \sqrt{k^{-1}\pi m}, \quad \xi^{k,n} = \sqrt{k^{-1}\pi n}, \quad n, m \in \mathbb{Z}$$

and for a given values of $k$ and $\delta$, the set of indices is taken to be

$$\Lambda := \{(m, n) \in \mathbb{Z}^2 \mid |p(x^{k,m}, \xi^{k,n})| < \delta\}.$$
The model scattering problem

Perfectly matched layers
<figure><img src="image_8.png" /><figcaption>FIGURE 1. Representation of $\Lambda$ for $k = 400$ in the homogeneous example. Top-left: $\delta = 0.1$. Top-right: $\delta = 0.4$. Bottom-left: $\delta = 0.6$. The horizontal and vertical distances between dots are $\sqrt{\pi/k}$.</figcaption></figure>

5.2. **Homogeneous medium with analytical solution.** Our first example concerns the case where $\alpha = \mu = 1$. The right-hand side is given by

$$f_k := P_k(\phi e^{ikx}) = (\phi'' + 2ik\phi')e^{ikx},$$

where $\phi$ is the only even function in $C^3(\mathbb{R}^d)$ such that $\phi = 0$ on $(-\infty, -3/4)$, $\phi = 1$ on $(-1/2, 0)$, and $\phi$ is a polynomial of degree 7 in $(-3/4, -1/2)$, and the associated solution is $u_k(x) := \phi(x)e^{ikx}$. Results are reported in Table 2, where we have chosen $\delta$ to maintain a constant accuracy for different values of $k$. Figure 1 represents the points $(x^{k,m}, \xi^{k,n})$ included in the space $\Lambda$ for different values of $\delta$ for the case $k = 400$.

<table><caption>TABLE 2. Numerical results for different $k$ while varying $\delta$</caption><thead><tr><th>k</th><th colspan="4">Relative H<sup>1</sup><sub>k</sub> error (δ) [N<sub>dofs</sub>]</th></tr></thead><tbody><tr><td>20</td><td>3.3539e-01</td><td>(0.5)</td><td>[60]</td><td>2.2356e-05</td><td>(2.0)</td><td>[177]</td><td>6.6722e-06</td><td>(4.0)</td><td>[249]</td></tr><tr><td>50</td><td>2.3180e-02</td><td>(0.4)</td><td>[38]</td><td>1.7831e-05</td><td>(1.0)</td><td>[229]</td><td>4.4104e-06</td><td>(3.0)</td><td>[419]</td></tr><tr><td>100</td><td>4.4322e-01</td><td>(0.3)</td><td>[96]</td><td>3.9496e-05</td><td>(0.8)</td><td>[278]</td><td>2.5610e-06</td><td>(2.0)</td><td>[645]</td></tr><tr><td>200</td><td>3.9807e-02</td><td>(0.2)</td><td>[58]</td><td>1.8555e-05</td><td>(0.6)</td><td>[334]</td><td>4.4586e-06</td><td>(1.0)</td><td>[825]</td></tr><tr><td>400</td><td>3.6922e-01</td><td>(0.1)</td><td>[74]</td><td>3.6289e-05</td><td>(0.4)</td><td>[532]</td><td>2.3452e-05</td><td>(0.6)</td><td>[646]</td></tr></tbody></table>

On Figure 2, we present the convergence history of the method as $\delta$ is increased for different values of $k$. These curves illustrate the fact that the method converges for any fixed $k$ when increasing the number of coherent Gaussian states, as predicted by our theory. Besides, on Figure 3 we represent the values of $\delta$ and corresponding number of freedom $N_{dofs}$ required to achieve a value of about $2 \times 10^{-5}$ (second column of Table 2) for different frequencies. Interestingly, the expected rates, namely, $\delta \sim k^{-1/2}$ and $N_{dofs} \sim k^{1/2}$ are observed.
Finally, we employ a perfectly matched layer (PML).

In essence, we artificially modify the coefficients for $|x| > R'$ to add absorption.

This approach is standard, and do not modify the value of $u$ for $|x| < R'$.

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

We end up with a problem set over $H^1(\mathbb{R}^d)$ with smooth coefficients, no boundary conditions, and $f \in L^2(\mathbb{R}^d)$ compactly supported.
<figure><img src="image_9.png" /><figcaption>FIGURE 2. Convergence histories in the homogeneous example</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 3. Computational cost for a fixed accuracy in the homogeneous example</figcaption></figure>

5.3. **Scattering in an heterogeneous medium.** We now focus on the case where $\alpha = 1$ and $\mu$ is the only even $C^3$ functions that equals 2 on $[0, 0.7]$, equals 1 on $[0.8, +\infty]$ and is a polynomial of degree 7 on $[0.7, 0.8]$. We select the right-hand side $f_k(x) := k^2(\mu - 1)e^{ikx}$. Here, the analytical solution is not available, and instead, we rely on a reference solution computed by a Lagrange finite element method of order 4 on a grid with $h = 0.02 \cdot k^{-9/8}$ in order to avoid the pollution effect (see, e.g., [26]). The results are listed on Table 3, and Figure 6 represents the phase-space points included in the space $\Lambda$ for the case $k = 200$.

We provide the same figures than in the previous experiment, and arrive at similar observation. First, Figure 4 illustrates the convergence of the method for fixed values of $k$ as the number of degrees of freedom is increased. Second, we compare the values of $\delta$ and $N_{\text{dofs}}$ to achieve an accuracy about $10^{-4}$ (last column of Table 3) for different frequencies. In this example too, the expected rates are numerically observed.
# The model scattering problem

## Key assumptions
<table><caption>TABLE 3. Numerical results for $k = 20, 50, 100, 200$ while varying $\delta$</caption><thead><tr><th>k</th><th colspan="6">Relative H<sup>1</sup><sub>k</sub> error (δ) [N<sub>dofs</sub>]</th></tr></thead><tbody><tr><td>20</td><td>1.7861e-01 (2.0) [183]</td><td>4.1982e-02 (4.0) [255]</td><td>1.6388e-03 (8.0) [399]</td><td>1.3463e-04 (12.0) [569]</td><td></td><td></td></tr><tr><td>50</td><td>3.0122e-01 (2.0) [337]</td><td>2.5263e-02 (4.0) [511]</td><td>1.0504e-03 (6.0) [695]</td><td>4.5593e-04 (8.0) [841]</td><td></td><td></td></tr><tr><td>100</td><td>2.0614e-01 (1.0) [376]</td><td>7.7805e-03 (2.0) [665]</td><td>1.7647e-04 (4.0) [995]</td><td>1.3300e-04 (6.0) [1309]</td><td></td><td></td></tr><tr><td>200</td><td>3.8799e-01 (0.5) [250]</td><td>8.5455e-03 (1.0) [748]</td><td>6.9635e-04 (2.0) [1233]</td><td>1.2300e-04 (4.0) [1879]</td><td></td><td></td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>FIGURE 4. Convergence histories in the heterogeneous example</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 5. Computational cost for a fixed accuracy in the heterogeneous example</figcaption></figure>

## 6. CONCLUSION

We propose a new family of finite-dimensional spaces to approximate the solutions to high-frequency Helmholtz problems. These discretization spaces are spanned by Gaussian coherent
Our model problem is to find $u \in H^1(\mathbb{R}^d)$ such that

$$ -k^2 \mu u - \nabla \cdot (A \nabla u) = f $$

where

$f \in L^2(\mathbb{R}^d)$ is compactly supported, $\mu$ and $A$ are $C^\infty(\mathbb{R}^d, \mathbb{C})$ and bounded along with all their derivatives.
<figure><img src="image_13.png" /><figcaption>FIGURE 6. Representation of $\Lambda$ for $k = 200$ in the heterogeneous example. Top-left: $\delta = 0.5$. Top-right: $\delta = 1$. Bottom-left: $\delta = 2$. The horizontal and vertical distances between dots are $\sqrt{\pi}/k$.</figcaption></figure>

states, which have the key property to be micro-localised in phase space. This unique feature allows to carefully select which Gaussian coherent states are included in the discretization, leading to a frequency-aware discretization space specifically tailored to approximate the solutions to scattering problems efficiently.

Our key findings correspond to two types of approximability results. First, assuming for simplicity that the problem is non-trapping, for general $L^2$ right-hand sides, we show that the Gaussian state approximation converges and provides a uniform error for all frequencies if the number of degrees of freedom grows as $(kR)^d$. This result is similar to approximation results available for finite element discretizations. Our second result applies when the right-hand side corresponds to a plane-wave scattering problem. In this case, we show that it is sufficient for the number of degrees of freedom to grow only as $(kR)^{d-1/2}$ to achieve a constant accuracy for increasing frequencies. To the best of the authors knowledge, this last estimate suggests that the proposed discretization space requires substantially less degrees of freedom than any available method in the literature for high-frequency scattering problems in general smooth heterogeneous media.

We also present a set of numerical examples where our Gaussian state spaces are coupled with a least-squares variational formulation. Although the setting is elementary, these examples successfully illustrate the key features of our abstract analysis.

While we believe that the proposed results are very encouraging for a further development of the proposed method, there still remain several challenges that we would like to address in future works. First (i), we have chosen to focus on a least-squares method, because it is simpler to analyse than a Galerkin formulation. However, least-squares methods are typically poorly-conditioned as compared to their Galerkin counterparts. While there is no reason to believe that the proposed discrete space would not work with a Galerkin variational formulation, the analysis appears to be substantially more complex. Second (ii), our convergence analysis is currently split into two distinct cases: arbitrary precision for a fixed frequency (with about $(kR)^d$ DOFs) and or fixed
Our model problem is to find $u \in H^1(\mathbb{R}^d)$ such that

$$ -k^2 \mu u - \nabla \cdot (A \nabla u) = f $$

where

$f \in L^2(\mathbb{R}^d)$ is compactly supported, $\mu$ and $A$ are $C^\infty(\mathbb{R}^d, \mathbb{C})$ and bounded along with all their derivatives.

In the second part of the talk, we will also specifically use the assumption that

$$ f = \phi e^{ik \cdot x} $$

for a smooth compactly supported function $\phi$.
precision for large frequencies (where the number DOFs scales as $(kR)^{d-1/2}$). In future works, we will focus on unifying the convergence analysis over the whole frequency range. Third (iii), especially for three-dimensional problems, computing the entries of the matrices associated with Gaussian coherent states is not a trivial task, and specific algorithms should be employed. Besides (iv), although the entries of the discrete matrices decrease super-algebraically away from their diagonal, they are still dense in principle. While we are convinced that efficient truncation or compression methods can be employed, this is still to be analysed. In addition (v), the resulting linear systems (after possible truncation and/or compression) will likely have different properties than usual volume methods. As a result, the design of specific preconditioners will probably be required. Finally (vi), although we do not see a convenient way to handle boundary conditions with the proposed Gaussian coherent states, a promising solution would be to replace the Gaussian by another smooth, but compactly supported, function in the definition of the Gabor frame.

## APPENDIX A. INTERPLAY BETWEEN GAUSSIAN STATES AND DIFFERENTIAL OPERATORS

**Lemma A.1 (Action of a differential operator).** Given $\beta \in \mathbb{N}^d$ and smooth coefficients $(A_\alpha)_{\alpha \in \mathbb{N}^d}$ that are bounded along with all their derivatives, consider the differential operator

$$ (A.1) \qquad \mathcal{P}_{\hbar, \beta, A} := \sum_{\substack{\alpha \in \mathbb{N}^d \\ \alpha \le \beta}} \hbar^{[\alpha]} A_\alpha \partial^\alpha. $$

Then for all $[x_0, \xi_0] \in \mathbb{Z}^{2d}$ there exists a smooth bounded function $g_{\hbar, \beta, A, x_0, \xi_0}$ such that

$$ (A.2) \qquad (\mathcal{P}_{\hbar, \beta, A} \Phi_{\hbar, x_0, \xi_0})(\mathbf{x}) = g_{\hbar, \beta, A, x_0, \xi_0}(\mathbf{x}) \Phi_{\hbar, x_0, \xi_0}(\mathbf{x}) \quad \forall \mathbf{x} \in \mathbb{R}^d. $$

In addition, we have

$$ (A.3) \qquad \left| (\partial^\gamma g_{\hbar, \beta, A, x_0, \xi_0})(\mathbf{x}) \right| \le C_{\beta, A, \gamma} \left( 1 + |\mathbf{x} - \mathbf{x}_0|^{[\beta]} + |\xi_0|^{[\beta]} \right) \quad \forall \gamma \in \mathbb{N}^d. $$

In (A.3), the constant $C_{\beta, A, \gamma}$ depends on $A$ only through bounds on a finite number of derivatives of the functions $A_\alpha$ for $\alpha \le \beta$, where this number of derivatives depends on $\beta$ and $\gamma$. We will use similar notations for constants in the rest of the appendix.

*Proof.* Recalling [5, Lemma A.2], for $\alpha \in \mathbb{N}$, we have

$$ (\partial^\alpha \Phi_{\hbar, x_0, \xi_0})(\mathbf{x}) = \hbar^{-[\alpha]/2} q_\alpha(\hbar^{-1/2}(\mathbf{x} - \mathbf{x}_0 + i\xi_0)) \Phi_{\hbar, x_0, \xi_0}(\mathbf{x}), $$

where $q_\alpha$ is a polynomial of degree less than or equal to $[\alpha]$. Hence, we readily see that (A.2) holds true with

$$ g_{\hbar, \beta, A, x_0, \xi_0}(\mathbf{x}) := \sum_{\substack{\alpha \in \mathbb{N}^d \\ \alpha \le \beta}} \hbar^{[\alpha]/2} A_\alpha(\mathbf{x}) q_\alpha(\hbar^{-1/2}(\mathbf{x} - \mathbf{x}_0 + i\xi_0)). $$

To establish (A.3), we start with Leibniz' rule

$$ (\partial^\gamma g_{\hbar, \beta, A, x_0, \xi_0})(\mathbf{x}) := \sum_{\substack{\alpha \in \mathbb{N}^d \\ \alpha \le \beta}} \hbar^{[\alpha]/2} \sum_{\substack{\gamma \in \mathbb{N}^d \\ \delta \le \gamma}} \left( \begin{array}{cc} \delta & \gamma \end{array} \right) \hbar^{-[\delta]/2} (\partial^\gamma - \delta A_\alpha)(\mathbf{y}) (\partial^\delta q_\alpha)(\hbar^{-1/2}(\mathbf{x} - \mathbf{x}_0 + i\xi_0)), $$

so that

$$ \left| (\partial^\gamma g_{\hbar, \beta, A, x_0, \xi_0})(\mathbf{y}) \right| \le C_{\beta, A, \gamma} \sum_{\substack{\alpha \in \mathbb{N}^d \\ \alpha \le \beta}} \hbar^{[\alpha]/2} \sum_{\substack{\delta \in \mathbb{N}^d \\ \delta \le \gamma}} \hbar^{-[\delta]/2} \left| (\partial^\delta q_\alpha)(\hbar^{-1/2}(\mathbf{x} - \mathbf{x}_0 + i\xi_0)) \right|. $$

Next, we employ the fact that $\partial^\delta q_\alpha = 0$ whenever $\delta > \alpha$, and the estimate

$$ \left| (\partial^\delta q_\alpha)(\hbar^{-1/2}(\mathbf{x} - \mathbf{x}_0 + i\xi_0)) \right| \le C_{\alpha, \delta} \left( 1 + \hbar^{-([\alpha]-[\delta])/2} |\mathbf{x} - \mathbf{x}_0 + i\xi_0|^{[\alpha]-[\delta]} \right), $$
Our model problem is to find $u \in H^1(\mathbb{R}^d)$ such that

$$ -k^2 \mu u - \nabla \cdot (A \nabla u) = f $$

where

$f \in L^2(\mathbb{R}^d)$ is compactly supported, $\mu$ and $A$ are $C^\infty(\mathbb{R}^d, \mathbb{C})$ and bounded along with all their derivatives.

In the second part of the talk, we will also specifically use the assumption that

$$ f = \phi e^{ik \cdot x} $$

for a smooth compactly supported function $\phi$.

For the sake of simplicity, we assume that the coefficients are non-trapping:

$$ \|\nabla u\|_{L^2(\mathbb{R}^d)} \le C \|f\|_{L^2(\mathbb{R}^d)} \quad \forall k \ge 1. $$
leading to

$$
\begin{align*}
| (\partial^\gamma g_{h,\beta,A,x_0,\xi_0})(\mathbf{x}) | &\le C_{\beta,A,\gamma} \sum_{\substack{\alpha \in \mathbb{N}^d \\ \alpha \le \beta}} \hbar^{[\alpha]/2} \sum_{\substack{\gamma \in \mathbb{N}^d \\ \delta \le \gamma}} (1 + (\hbar^{-1/2} |\mathbf{x} - \mathbf{x}_0 + i\xi_0|)^{[\alpha]-[\delta]}) \\
&\le C_{\beta,A,\gamma} \sum_{\substack{\alpha \in \mathbb{N}^d \\ \alpha \le \beta}} \hbar^{[\alpha]/2} (1 + (\hbar^{-1/2} |\mathbf{x} - \mathbf{x}_0 + i\xi_0|)^{[\alpha]}) \\
&\le C_{\beta,A,\gamma} \sum_{\substack{\alpha \in \mathbb{N}^d \\ \alpha \le \beta}} (1 + |\mathbf{x} - \mathbf{x}_0 + i\xi_0|^{[\alpha]}) \\
&\le C_{\beta,A,\gamma} (1 + |\mathbf{x} - \mathbf{x}_0 + i\xi_0|^{[\beta]})
\end{align*}
$$

Then, (A.3) follows since

$$
|\mathbf{x} - \mathbf{x}_0 + i\boldsymbol{\xi}_0|^{[\beta]} \le C_{\beta} (|\mathbf{x} - \mathbf{x}_0|^{[\beta]} + |\boldsymbol{\xi}_0|^{[\beta]}). \quad \square
$$

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

For the reader's convenience, we recall the following basic fact.

**Proposition A.2** (Moments of a Gaussian). *Consider the Gaussian function*

$$
(A.4) \qquad G_{\hbar}(\mathbf{y}) := \hbar^{-d/2} e^{-\frac{1}{\hbar} |\mathbf{y}|^2}.
$$

Then, for $\beta, \gamma \in \mathbb{N}^d$, we have

$$
(A.5) \quad \int_{\mathbb{R}^d} |\mathbf{x}^{\beta} (\partial^{\gamma} G_{\hbar}) (\mathbf{x} - \mathbf{x}_0)| d\mathbf{x} \le C_{\beta, \gamma} \hbar^{-[\gamma]/2} (\hbar^{[\beta]/2} + |\mathbf{x}_0|^{[\beta]}) \quad \forall \mathbf{x}_0 \in \mathbb{R}^d.
$$

We now show that Gaussian states are localised in phase-space.

**Lemma A.3** (Quasi orthogonality). Let $\mathcal{P}_{h,\beta,A}$ be as in (A.1). Then, for all $\gamma \in \mathbb{N}^d$, $[\mathbf{x}_0, \boldsymbol{\xi}_0] \in \mathbb{R}^{2d}$ and $[\mathbf{x}'_0, \boldsymbol{\xi}'_0] \in \mathbb{R}^{2d}$, we have

$$
| (\mathbf{x}^{\gamma} \mathcal{P}_{h,\beta,A} \Phi_{h,x_0,\boldsymbol{\xi}_0}, \Phi_{h,x'_0,\boldsymbol{\xi}'_0}) | \le C_{\gamma,\beta,A,m} \hbar^{m/2} (1 + |\mathbf{x}_0 + \mathbf{x}'_0|^{[\gamma]}) \frac{1 + |\boldsymbol{\xi}_0|^{[\beta]} + |\mathbf{x}_0 - \mathbf{x}'_0|^{[\beta]}}{1 + |\boldsymbol{\xi}_0 - \boldsymbol{\xi}'_0|^m} e^{-\frac{1}{4\hbar} |\mathbf{x}_0 - \mathbf{x}'_0|^2} \quad \forall m \in \mathbb{N}.
$$

In particular,

$$
|(\mathbf{x}^{\gamma} \mathcal{P}_{h,\beta,A} \Phi_{h,x_0,\boldsymbol{\xi}_0}, \Phi_{h,x'_0,\boldsymbol{\xi}'_0})| \le C_{\gamma,\beta,A,m} \hbar^{m/2} (1 + |\mathbf{x}_0 + \mathbf{x}'_0|^{[\gamma]}) \frac{1 + |\boldsymbol{\xi}_0|^{[\beta]}}{1 + |[\mathbf{x}_0, \boldsymbol{\xi}_0] - [\mathbf{x}'_0, \boldsymbol{\xi}'_0]|^m} \quad \forall m \in \mathbb{N}.
$$

*Proof.* We have

$$
\Phi_{h,x_0,\xi_0}(\mathbf{x}) \overline{\Phi_{h,x'_0,\xi'_0}(\mathbf{x})} = \theta \hbar^{-d/2} e^{-\frac{1}{4\hbar} |\mathbf{x}_0 - \mathbf{x}'_0|^2} e^{-\frac{1}{\hbar} |\mathbf{x} - \frac{1}{2}(\mathbf{x}_0 + \mathbf{x}'_0)|^2} e^{\frac{i}{\hbar} (\boldsymbol{\xi}_0 - \boldsymbol{\xi}'_0) \cdot \mathbf{x}}
$$

where $\theta \in \mathbb{C}$ with $|\theta| = \pi^{-d/2}$. Recalling (A.2) and the notation (A.4) from Proposition A.2, we
have

$$
\mathbf{x}^{\gamma}(\mathcal{P}_{h,\beta,A} \Phi_{h,x_0,\boldsymbol{\xi}_0})(\mathbf{x}) \overline{\Phi_{h,x'_0,\boldsymbol{\xi}'_0}(\mathbf{x})} = \theta \mathbf{x}^{\gamma} g_{h,\beta,A,x_0,\boldsymbol{\xi}_0}(\mathbf{x}) e^{-\frac{1}{4\hbar} |\mathbf{x}_0 - \mathbf{x}'_0|^2} G_{\hbar} \left( \mathbf{x} - \frac{1}{2}(\mathbf{x}_0 + \mathbf{x}'_0) \right) e^{\frac{i}{\hbar} (\boldsymbol{\xi}_0 - \boldsymbol{\xi}'_0) \cdot \mathbf{x}}
$$

so that

$$
| (\mathbf{x}^{\gamma} \mathcal{P}_{h,\beta,A} \Phi_{h,x_0,\boldsymbol{\xi}_0}, \Phi_{h,x'_0,\boldsymbol{\xi}'_0}) | \le C e^{-\frac{1}{4\hbar} |\mathbf{x}_0 - \mathbf{x}'_0|^2} |I|,
$$

with

$$
I := \hbar^{-d/2} \int_{\mathbb{R}^d} \mathbf{x}^{\gamma} g_{h,\beta,A,x_0,\xi_0}(\mathbf{x}) G_{h,x_0,x'_0}(\mathbf{x}) e^{\frac{i}{\hbar}(\boldsymbol{\xi}_0 - \boldsymbol{\xi}'_0) \cdot \mathbf{x}} d\mathbf{x},
$$

where $G_{h,x_0,x'_0}(\mathbf{x}) := G_{\hbar}(\mathbf{x} - \frac{1}{2}(\mathbf{x}_0 + \mathbf{x}'_0))$.
Convergence for general right-hand sides
If $\xi_0 = \xi'_0$, this concludes the proof since

$$
\begin{align*}
|I| &\le \hbar^{-d/2} \int_{\mathbb{R}^d} |\mathbf{x}^{\gamma} g_{\hbar, \beta, A, \mathbf{x}_0, \xi_0}(\mathbf{x})| G_{\hbar, \mathbf{x}_0, \mathbf{x}'_0}(\mathbf{x}) \mathrm{d}\mathbf{x} \\
&\le C_{\beta} \int_{\mathbb{R}} |\mathbf{x}|^{[\gamma]} (1 + |\mathbf{x} - \mathbf{x}_0|^{[\beta]}) G_{\hbar, \mathbf{x}_0, \mathbf{x}'_0}(\mathbf{x}) \mathrm{d}\mathbf{x} \\
&= C_{\beta} \int_{\mathbb{R}} \left|\mathbf{y} + \frac{1}{2}(\mathbf{x}_0 + \mathbf{x}'_0)\right|^{[\gamma]} \left(1 + \left|\mathbf{y} + \frac{1}{2}(\mathbf{x}'_0 - \mathbf{x}_0)\right|^{[\beta]}\right) G_{\hbar}(\mathbf{y}) \mathrm{d}\mathbf{y} \\
&\le C_{\beta, \gamma} (1 + |\mathbf{x}_0 + \mathbf{x}'_0|^{[\gamma]})(1 + |\mathbf{x}_0 - \mathbf{x}'_0|^{[\beta]}).
\end{align*}
$$

If $\xi_0 \neq \xi'_0$, then there exists a $j \in \{1, \dots, d\}$ such that $|\langle \xi_0 - \xi'_0 \rangle_j| \ge 1/\sqrt{d}|\xi_0 - \xi'_0|$. We integrate $m$ times by part with respect to $\mathbf{x}_j$, leading to

$$
\begin{align*}
I &= \hbar^{-d/2} \left( \frac{i}{\hbar} (\boldsymbol{\xi}_0 - \boldsymbol{\xi}'_0)_j \right)^{-m} \int_{\mathbb{R}^d} \partial_j^m (\boldsymbol{x}^\gamma g_{\hbar, \beta, A, \boldsymbol{x}_0, \boldsymbol{\xi}_0} G_{\hbar, \boldsymbol{x}_0, \boldsymbol{x}'_0}) (\boldsymbol{x}) e^{\frac{i}{\hbar} (\boldsymbol{\xi}_0 - \boldsymbol{\xi}'_0) \cdot \boldsymbol{x}} \mathrm{d}\boldsymbol{x} \\
  &\le C \hbar^m |\boldsymbol{\xi}_0 - \boldsymbol{\xi}'_0|^{-m} I'.
\end{align*}
$$

where

$$
\begin{align*}
|I'| &:= \left| \hbar^{-d/2} \int_{\mathbb{R}^d} \partial_j^m (\mathbf{x}^\gamma g_{\hbar, \beta, A, \mathbf{x}_0, \xi_0} G_{\hbar, \mathbf{x}_0, \mathbf{x}'_0}) (\mathbf{x}) e^{\frac{i}{\hbar} (\boldsymbol{\xi}_0 - \boldsymbol{\xi}'_0) \cdot \mathbf{x}} \mathrm{d}\mathbf{x} \right| \\
&= \left| \hbar^{-d/2} \sum_{\ell+\ell'+\ell''=m} \frac{m!}{\ell! \cdot \ell'! \ell''!} \int_{\mathbb{R}^d} \partial_j^{\ell} \mathbf{x}^{\gamma} (\partial_j^{\ell'} g_{\hbar, \beta, A, \mathbf{x}_0, \xi_0}) (\mathbf{x}) (\partial_j^{\ell''} G_{\hbar, \mathbf{x}_0, \mathbf{x}'_0}) (\mathbf{x}) e^{\frac{i}{\hbar} (\boldsymbol{\xi}_0 - \boldsymbol{\xi}'_0) \cdot \mathbf{x}} \mathrm{d}\mathbf{x} \right| \\
&\le C_m \sum_{\ell+\ell'+\ell''=m} I_{\ell, \ell', \ell''}.
\end{align*}
$$

with

$$
I_{\ell, \ell', \ell''} := \hbar^{-d/2} \int_{\mathbb{R}^d} |\partial_j^\ell \mathbf{x}^\gamma| |(\partial_j^{\ell'} g_{\hbar, \beta, A, \mathbf{x}_0, \xi_0})(\mathbf{x})| |\partial_j^{\ell''} G_{\hbar, \mathbf{x}_0, \mathbf{x}'_0}(\mathbf{x})| \mathrm{d}\mathbf{x}.
$$

We then employ (A.3) and (A.5), showing that

$$
\begin{align*}
|I_{\ell, \ell', \ell''}| &\le C_{\beta, A, m, \gamma} \hbar^{-d/2} \int_{\mathbb{R}^d} (1 + |\mathbf{x}|^{[\gamma]})(1 + |\mathbf{x} - \mathbf{x}_0|^{[\beta]} + |\boldsymbol{\xi}_0|^{[\beta]}) |\partial_j^{\ell''} G_{\hbar, \mathbf{x}_0, \mathbf{x}'_0}(\mathbf{x})| \mathrm{d}\mathbf{x} \\
&= C_{\beta, A, m, \gamma} \hbar^{-d/2} \int_{\mathbb{R}^d} (1 + |\mathbf{x} + \mathbf{x}_0|^{[\gamma]})(1 + |\mathbf{x}|^{[\beta]} + |\boldsymbol{\xi}_0|^{[\beta]}) \left|\partial_j^{\ell''} (G_{\hbar}) \left(\mathbf{x} + \frac{1}{2}(\mathbf{x}_0 - \mathbf{x}'_0)\right)\right| \mathrm{d}\mathbf{x} \\
&\le C_{\beta, A, m, \gamma} \hbar^{-\ell/2} (1 + |\mathbf{x}_0 + \mathbf{x}'_0|^{[\gamma]}) \left( (1 + |\boldsymbol{\xi}_0|^{[\beta]}) + \left|\frac{1}{2}(\mathbf{x}_0 - \mathbf{x}'_0)\right|^{[\beta]} \right) \\
&\le C_{\beta, A, m, \gamma} \hbar^{-\ell/2} (1 + |\mathbf{x}_0 + \mathbf{x}'_0|^{[\gamma]})(1 + |\boldsymbol{\xi}_0|^{[\beta]} + |\mathbf{x}_0 - \mathbf{x}'_0|^{[\beta]}).
\end{align*}
$$

□

**Lemma A.4 (Action of $(P_h - p)$ on Gaussian states).** Consider the second-order differential operator

$$
(A.6) \qquad P_h := \hbar^2 \sum_{k,l=1}^{d} (A_h)_{kl} \partial_{kl} + i\hbar b_h \cdot \nabla + c_h,
$$

together with its symbol

$$
p_h(\mathbf{x}, \boldsymbol{\xi}) = A_h(\mathbf{x})\boldsymbol{\xi} \cdot \boldsymbol{\xi} + b_h(\mathbf{x}) \cdot \boldsymbol{\xi} + c_h(\mathbf{x}),
$$

where $A_h, b_h$ and $c_h$ are smooth functions that are bounded, along with all their derivatives. For $\mathbf{x}_0, \boldsymbol{\xi}_0 \in \mathbb{R}^d$, we have

$$
(A.7) \qquad (P_h - p_h(\mathbf{x}_0, \boldsymbol{\xi}_0)) \Phi_{\hbar, \mathbf{x}_0, \boldsymbol{\xi}_0} = r_{\hbar, \mathbf{x}_0, \boldsymbol{\xi}_0} \Phi_{\hbar, \mathbf{x}_0, \boldsymbol{\xi}_0}
$$
Convergence for general right-hand sides

The Gabor frame
with

$$
\begin{align*}
r_{h,x_0,\xi_0}(\mathbf{x}) &:= \hbar \operatorname{tr} \mathbf{A}_h(\mathbf{x}) \\
&\quad + i[(\mathbf{A}_h + \mathbf{A}_h^T)(\mathbf{x})\xi_0 + \mathbf{b}_h(\mathbf{x})] \cdot (\mathbf{x} - \mathbf{x}_0) \\
&\quad + (\mathbf{A}_h(\mathbf{x}) - \mathbf{A}_h(\mathbf{x}_0))\xi_0 \cdot \xi_0 + i(\mathbf{b}_h(\mathbf{x}) - \mathbf{b}_h(\mathbf{x}_0)) \cdot \xi_0 + c_h(\mathbf{x}) - c_h(\mathbf{x}_0) \\
&\quad + \mathbf{A}_h(\mathbf{x})(\mathbf{x} - \mathbf{x}_0) \cdot (\mathbf{x} - \mathbf{x}_0).
\end{align*}
$$

In particular, there exists smooth functions $\alpha_{h,x_0,\xi_0,\beta}$ such that

$$
(A.8) \qquad r_{h,x_0,\xi_0}(\mathbf{x}) = \sum_{\substack{\beta \in \mathbb{N}^d \\ |\beta| \le 2}} \hbar^{(1-|\beta|)/2} \alpha_{h,x_0,\xi_0,\beta}(\mathbf{x}) (\mathbf{x} - \mathbf{x}_0)^{\beta}
$$

and

$$
(A.9) \qquad \| \partial^\gamma \alpha_{h, x_0, \xi_0, \beta} \|_{L^\infty(\mathbb{R}^d)} \le C_{\beta, \gamma} (1 + |\xi_0|^2).
$$

Proof. We obtain (A.7) by straightforward computations. We easily identify that

$$
\alpha_{h,x_0,\xi_0,0} = \hbar^{1/2} \operatorname{tr} \mathbf{A}_h \qquad \alpha_{h,x_0,\xi_0,e_j+e_l} = \hbar^{1/2} (\mathbf{A}_h)_{jl}.
$$

For the terms with $[\beta] = 1$ we further write that for $\varphi \in \{(\mathbf{A}_h)_{j\ell}, (\mathbf{b}_h)_{\ell}, c_h\}$, there exist smooth functions $\varphi^{x_0,k}$ such that

$$
\varphi(\mathbf{x}) - \varphi(\mathbf{x}_0) = \sum_{k=1}^{d} \varphi^{x_0,k}(\mathbf{x}) \cdot (\mathbf{x} - \mathbf{x}_0),
$$

so that

$$
\alpha_{x_0, \xi_0, e_k} = \mathbf{A}_{\bar{h}}^{x_0, k} \xi_0 \cdot \xi_0 + i \mathbf{b}_{\bar{h}}^{x_0, k} \cdot \xi_0 + c_{\bar{h}}^{x_0, k} + i \left( \sum_{j=1}^{d} \left( (\mathbf{A}_{\bar{h}})_{kj} + (\mathbf{A}_{\bar{h}})_{jk} \right) \xi_{0,j} + (\mathbf{b}_{\bar{h}})_k \right)
$$

for $k \in \{1, \dots, d\}$.

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

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

We will now explore further the action of $P_h$ on $\Phi_{h,x_0,\xi_0}$, and show that it is close to $p(x_0, \xi_0)\Phi_{h,x_0,\xi_0}$

**Lemma A.5 (Action of $(P_h - p_h)^L$ on Gaussian states).** For $L \in \mathbb{N}$, we have

$$
(P_h - p_h(\mathbf{x}_0, \boldsymbol{\xi}_0))^L \Phi_{h, \mathbf{x}_0, \boldsymbol{\xi}_0} = r_{h, \mathbf{x}_0, \boldsymbol{\xi}_0, L} \Phi_{h, \mathbf{x}_0, \boldsymbol{\xi}_0},
$$

where $r_{h,x_0,\xi_0,L}$ can be written in the form

$$
r_{h,x_0,\xi_0,L}(\mathbf{x}) = \sum_{\substack{\beta \in \mathbb{Z}^d \\ |\beta| \le 2L}} \hbar^{(L-|\beta|)/2} \alpha_{h,x_0\xi_0,L,\beta}(\mathbf{x}) (\mathbf{x} - \mathbf{x}_0)^{\beta},
$$

where the functions $\alpha_{h,x_0,\xi_0,\beta}$ are smooth and satisfy

$$
(A.10) \qquad \| \partial^\gamma \alpha_{h, x_0, \xi_0, L, \beta} \|_{L^\infty(\mathbb{R}^d)} \le C_{L, \beta, \gamma} (1 + |\xi_0|^2)^L.
$$

Proof. We shall prove this result by induction. The case $L = 0$ trivially holds with $\alpha_{x_0,\xi_0,0,0} = 1$.
The case $L = 1$ is treated in Lemma A.4. Hence, let us assume that $L \in \mathbb{N}$ is such that (A.10)
holds for all $L' \le L$. We have

$$
\begin{align*}
(P_h - p_h(\mathbf{x}_0, \boldsymbol{\xi}_0))^{L+1} \Phi_{h, \mathbf{x}_0, \boldsymbol{\xi}_0} &= (P_h - p_h(\mathbf{x}_0, \boldsymbol{\xi}_0)) (r_{h, \mathbf{x}_0, \boldsymbol{\xi}_0, L} \Phi_{h, \mathbf{x}_0, \boldsymbol{\xi}_0}) \\
&= [P_h - p_h(\mathbf{x}_0, \boldsymbol{\xi}_0), r_{h, \mathbf{x}_0, \boldsymbol{\xi}_0, L}] \Phi_{h, \mathbf{x}_0, \boldsymbol{\xi}_0} + r_{h, \mathbf{x}_0, \boldsymbol{\xi}_0, L} r_{h, \mathbf{x}_0, \boldsymbol{\xi}_0, 1} \Phi_{h, \mathbf{x}_0, \boldsymbol{\xi}_0},
\end{align*}
$$
Start with the Gaussian distribution

$$ \pi^{-d/4} e^{-\frac{1}{2}|y|^2} $$
Writing the operator $P_h$ in the form (2.2), we deduce that

$$ [P_h - p_h(x_0, \xi_0), r_h(x_0, \xi_0, L)] \Phi_{h, x_0, \xi_0} = \hbar^2 \sum_{j, \ell=1}^{d} a_{j\ell} \left( \frac{\partial^2 r_{h, x_0, \xi_0, L}}{\partial x_j \partial x_\ell} \Phi_{h, x_0, \xi_0} + \frac{\partial r_{h, x_0, \xi_0, L}}{\partial x_j} \frac{\partial \Phi_{h, x_0, \xi_0}}{\partial x_\ell} \right) \\ + i\hbar \sum_{j=1}^{d} b_j \frac{\partial r_{h, x_0, \xi_0, L}}{\partial x_j} \Phi_{h, x_0, \xi_0} \\ = \tilde{r}_{k, x_0, \xi_0, L+1} \Phi_{h, x_0, \xi_0}, $$

where

$$ \tilde{r}_{h, x_0, \xi_0, L+1} := \sum_{j, \ell=1}^{d} \left( \hbar^2 a_{j\ell} \frac{\partial^2 r_{k, x_0, \xi_0, L}}{\partial x_j \partial x_\ell} + \hbar \frac{\partial r_{h, x_0, \xi_0, L}}{\partial x_j} (i\xi_{0,\ell} - (x_\ell - x_{0,\ell})) \right) + i\hbar \sum_{j=1}^{d} b_j \frac{\partial r_{h, x_0, \xi_0, L}}{\partial x_j} $$

Using the induction hypothesis, we have

$$ \frac{\partial r_{h, x_0, \xi_0, L}}{\partial x_j}(x) = \sum_{\substack{\beta \in \mathbb{N}^d \\ [\beta] \le 2L}} \hbar^{(L-[\beta])/2} \frac{\partial \alpha_{h, x_0, \xi_0, L, \beta}}{\partial x_j}(x)(x-x_0)^\beta \\ + \sum_{\substack{\beta \in \mathbb{N}^d \\ [\beta] \le 2L}} \hbar^{(L-[\beta])/2} \beta_j \alpha_{h, x_0, \xi_0, L, \beta}(x)(x-x_0)^{\beta-e_j} \\ = \sum_{\substack{\beta \in \mathbb{N}^d \\ [\beta] \le 2L}} \hbar^{(L-[\beta])/2} \frac{\partial \alpha_{h, x_0, \xi_0, L, \beta}}{\partial x_j}(x)(x-x_0)^\beta \\ + \sum_{\substack{\beta \in \mathbb{N}^d \\ [\beta] \le 2L-1}} \hbar^{(L-1-[\beta])/2} (\beta_j+1) \alpha_{h, x_0, \xi_0, L, \beta+e_j}(x)(x-x_0)^\beta. $$

Setting $\alpha_{h,x_0,\xi_0,\beta} = 0$ if $[\beta] > 2L$, this gives us

$$ \hbar \frac{\partial r_{h, x_0, \xi_0, L}}{\partial x_j}(x) = \sum_{\substack{\beta \in \mathbb{N}^d \\ [\beta] \le 2L}} \hbar^{(L+1-[\beta])/2} \left( \hbar^{1/2} \frac{\partial \alpha_{h, x_0, \xi_0, L, \beta}}{\partial x_j}(x) + (\beta_j + 1) \alpha_{h, x_0, \xi_0, L, \beta+e_j}(x) \right) (x-x_0)^\beta, $$

these terms are multiplied either by a smooth function, a term of the order $|\xi_0|$ or a power of $(x - x_0)$, which always enters the induction for $L + 1$. The term with the second derivative is treated similarly.

For the remaining term, we simply write that

$$ r_{h, x_0, \xi_0, L} r_{h, x_0, \xi_0, 1} = \left( \sum_{\substack{\beta' \in \mathbb{Z}^{2d} \\ [\beta'] \le 2}} \hbar^{(1-[\beta'])/2} \alpha_{h, x_0, \xi_0, 1, \beta'} (x-x_0)^{\beta'} \right) \left( \sum_{\substack{\beta \in \mathbb{Z}^{2d} \\ [\beta] \le 2L}} \hbar^{(L-[\beta])/2} \alpha_{h, x_0, \xi_0, L, \beta} (x-x_0)^{\beta} \right) \\ = \sum_{\substack{\beta \in \mathbb{Z}^{2d} \\ [\beta'] \le 2}} \sum_{\substack{\beta' \in \mathbb{Z}^{2d} \\ [\beta'] \le 2L}} \hbar^{(L-[\beta])/2 + (1-[\beta'])/2} \alpha_{h, x_0, \xi_0, L, \beta} \alpha_{h, x_0, \xi_0, 1, \beta'} (x-x_0)^{\beta+\beta'}, $$

and the result follows.

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

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

**Proposition A.6 (Control of $(P_h - p_h)^L)$. We have**

$$ \| (P_h - p_h(x_0, \xi_0))^L \Phi_{h, x_0, \xi_0} \|^2 \le C_{A,b,c} (1 + |\xi_0|^2)^L \hbar^L $$

for all $[x_0, \xi_0] \in \mathbb{R}^{2d}$.
Start with the Gaussian distribution

$$ \pi^{-d/4} e^{-\frac{1}{2}|y|^2}. $$

Shift in position by $m\sqrt{\pi}$, with $m \in \mathbb{Z}^d$

$$ \pi^{-d/4} e^{-\frac{1}{2}|y-m\sqrt{\pi}|^2}. $$
Proof. We deduce from the previous lemma that

$$
\begin{align*}
\|(P_h - p_h(\mathbf{x}_0, \boldsymbol{\xi}_0))^L \Phi_{\hbar, \mathbf{x}_0, \boldsymbol{\xi}_0}\|^2 &\le C(1 + |\boldsymbol{\xi}_0|^2)^L (\pi\hbar)^{-d/2} \int_{\mathbb{R}^d} (|\mathbf{x} - \mathbf{x}_0|^2 + \hbar)^L e^{-\frac{1}{\hbar}|\mathbf{x} - \mathbf{x}_0|^2} d\mathbf{x} \\
&= C(1 + |\boldsymbol{\xi}_0|^2)^L \pi^{-d/2} \int_{\mathbb{R}^d} (\hbar|\mathbf{y}|^2 + \hbar)^L e^{-|\mathbf{y}|^2} d\mathbf{y} \\
&\le C(1 + |\boldsymbol{\xi}_0|^2)^L \hbar^L. \quad \square
\end{align*}
$$

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

REFERENCES

1. B. Adcock and D. Huybrechs, *Frames and numerical approximation*, SIAM Review **61** (2019), no. 3, 443–473.

2. J. P. Bérenger, *A perfectly matched layer for the absorption of electromagnetics waves*, J. Comput. Phys. **114** (1994), 185–200.

3. O. Cessenat and B. Després, *Using plane waves as base functions for solving time harmonic equations with the ultra weak variational formulation*, J. Comput. Acous. **11** (2003), no. 02, 227–238.

4. T. Chaumont-Frelet, D. Gallistl, S. Nicaise, and J. Tomezyk, *Wavenumber explicit convergence analysis for finite element discretizations of time-harmonic wave propagation problems with perfectly matched layers*, submitted, preprint hal-01887267, 2018.

5. T. Chaumont-Frelet and M. Ingremeau, *Decay of coefficients and approximation rates in Gabor Gaussian frames*, hal-03746979, 2022.

6. T. Chaumont-Frelet and F. Valentin, *A multiscale hybrid-mixed method for the Helmholtz equation in heterogeneous domains*, SIAM J. Numer. Anal. **58** (2020), no. 2, 1029–1067.

7. F. Collino and P. Monk, *The perfectly matched layer in curvilinear coordinates*, SIAM J. Sci. Comp. **19** (1998), no. 6, 2061–2090.

8. I. Daubechies, A. Grossmann, and Y. Meyer, *Painless nonorthogonal expansions*, J. Math. Phys. **27** (1986), no. 5, 1271–1283.

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

10. S. Dyatlov and M. Zworski, *Mathematical theory of scattering resonances*, vol. 200, American Mathematical Soc., 2019.

11. B. Engquist and O. Runborg, *Computational high frequency wave propagation*, Acta Numerica **2003** (2003), 181–266.

12. E. Faou, V. Gradinaru, and C. Lubich, *Computing semiclassical quantum dynamics with hagedorn wavepackets*, SIAM J. Sci. Comput. **31** (2009), no. 4, 3027–3041.

13. C. Farhat, I. Harari, and L. P. Franca, *The discontinuous enrichment method*, Comput. Methods Appl. Mech. Engrg. **190** (2001), no. 48, 6455–6479.

14. D. Gabor, *Theory of communication. part 1: The analysis of information*, Journal of the Institution of Electrical Engineers-Part III: Radio and Communication Engineering **93** (1946), no. 26, 429–441.

15. J. Galkowski, D. Lafontaine, and E. A. Spence, *Perfectly-matched-layer truncation is exponentially accurate at high frequency*, arXiv preprint arXiv:2105.07737 (2021).

16. J. Galkowski and E.A. Spence, *Does the Helmholtz boundary element method suffer from the pollution effect?*, arXiv:2201.09721, 2022.

17. J. Galkowski, E.A. Spence, and J. Wunsch, *Optimal constants in nontrapping resolvent estimates and applications in numerical analysis*, Pure and Applied Analysis **2** (2019), no. 1, 157–202.

18. C.J. Gittelson, R. Hiptmair, and I. Perugia, *Plane wave discontinuous Galerkin methods: analysis of the h-version*, ESAIM Math. Model. Numer. Anal. **43** (2009), 297–331.

19. V. Gradinaru and G.A. Hagedorn, *Convergence of a semiclassical wavepacket based time-splitting for the Schrödinger equation*, Numer. Math. **126** (2014), no. 1, 53–73.

20. V. Gradinaru and O. Rietmann, *Hagedorn wavepackets and Schrödinger equation with time-dependent, homogeneous magnetic field*, J. Comput. Phys. **445** (2021), 110581.

21. L. Greengard and V. Rokhlin, *A fast algorithm for particle simulations*, J. Comput. Phys. **73** (1987), no. 2, 325–348.

22. K. Gröchenig, *Foundations of time-frequency analysis*, Springer Science & Business Media, 2001.

23. W. Hackbusch, *Hierarchical matrices: algorithms and analysis*, vol. 49, Springer, 2015.

24. R. Hiptmair, A. Moiola, and I. Perugia, *Plane wave discontinuous Galerkin methods for the 2d Helmholtz equation: analysis of the p-version*, SIAM J. Numer. Anal. **49** (2011), no. 1, 264–284.

25. __________, *A survey of Trefftz methods for the Helmholtz equation*, pp. 237–279, Springer International Publishing, 2016.

26. F. Ihlenburg and I. Babuška, *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.

27. L.-M. Imbert-Gérard and B. Després, *A generalized plane-wave numerical method for smooth nonconstant coefficients*, IMA Journal of Numerical Analysis **34** (2014), no. 3, 1072–1103.
Start with the Gaussian distribution

$$ \pi^{-d/4} e^{-\frac{1}{2}|y|^2}. $$

Shift in position by $m\sqrt{\pi}$, with $m \in \mathbb{Z}^d$

$$ \pi^{-d/4} e^{-\frac{1}{2}|y-m\sqrt{\pi}|^2}. $$

Shift in momentum by $n\sqrt{\pi}$, with $n \in \mathbb{Z}^d$

$$ \Psi_{m,n} := \pi^{-d/4} e^{-\frac{1}{2}|y-m\sqrt{\pi}|^2} e^{in\sqrt{\pi} \cdot y}. $$
28. L.M. Imbert-Gerard, *Amplitude-based generalized plane waves: New quasi-trefftz functions for scalar equations in two dimensions*, SIAM Journal on Numerical Analysis **59** (2021), no. 3, 1663–1686.

29. D. Lafontaine, E.A. Spence, and J. Wunsch, *For most frequencies, strong tapping has a weak effect in frequency-domain scattering*, preprint arXiv:1903.12172, 2019.

30. ________, *Wavenumber-explicit convergence of the hp-fem for the full-space heterogeneous Helmholtz equation with smooth coefficients*, Comput. Math. Appl. **113** (2022), 59–69.

31. C. Lasser and C. Lubich, *Computing quantum dynamics in the semiclassical regime*, Acta Numerica **29** (2020), 229–401.

32. J. M. Melenk and I. Babuška, *The partition of unity finite element method: basic theory and applications*, Comput. Methods Appl. Mech. Engrg. **139** (1996), no. 1-4, 289–314.

33. J.M. Melenk, *Operator adapted spectral element methods I: harmonic and generalized harmonic polynomials*, Numer. Math. **84** (1999), 35–69.

34. J.M. Melenk and S. Sauter, *Convergence analysis for finite element discretizations of the Helmholtz equation with Dirichlet-to-Neumann boundary conditions*, Math. Comp. **79** (2010), no. 272, 1871–1914.

35. ________, *Wavenumber explicit convergence analysis for Galerkin discretizations of the Helmholtz equation*, SIAM J. Numer. Anal. **49** (2011), no. 3, 1210–1243.

36. P. Monk and D.-Q. Wang, *A least-squares method for the Helmholtz equation*, Comput. Methods Appl. Mech. Engrg. **175** (1999), no. 1-2, 121–136.

37. S. Nonnenmacher, *Spectral problems in open quantum chaos*, Nonlinearity **24** (2011), no. 12, R123.

38. S. Nonnenmacher and M. Zworski, *Quantum decay rates in chaotic scattering*, Acta mathematica **203** (2009), no. 2, 149–233.

39. H. Riou, P. Ladeveze, and B. Sourcis, *The multiscale VTCR approach applied to acoustics problems*, J. Comput. Acous. **16** (2008), no. 04, 487–505.

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

41. C. Schwab, *p- and hp- finite element methods*, Clarendon Press, 1998.

42. A. Tarantola, *Inversion of seismic reflection data in the acoustic approximation*, Geophysics **49** (1984), no. 8, 1259–1266.

43. M. Zworski, *Semiclassical analysis*, vol. 138, American Mathematical Soc., 2012.

44. ________, *Mathematical study of scattering resonances*, Bulletin of Mathematical Sciences **7** (2017), no. 1, 1–85.
Start with the Gaussian distribution

$$ \pi^{-d/4} e^{-\frac{1}{2}|y|^2}. $$

Shift in position by $m\sqrt{\pi}$, with $m \in \mathbb{Z}^d$

$$ \pi^{-d/4} e^{-\frac{1}{2}|y-m\sqrt{\pi}|^2}. $$

Shift in momentum by $n\sqrt{\pi}$, with $n \in \mathbb{Z}^d$

$$ \Psi_{m,n} := \pi^{-d/4} e^{-\frac{1}{2}|y-m\sqrt{\pi}|^2} e^{in\sqrt{\pi} \cdot y}. $$

The resulting family $\{\Psi_{m,n}\}_{[m,n] \in \mathbb{Z}^{2d}}$ is the Gabor frame. It is dense in $L^2(\mathbb{R}^d)$.

<figure><img src="image_13.png" /><figcaption>Gabor, 1946.</figcaption></figure>
We can scale the family with the change of variable $x = k^{1/2}y$, leading to

$$ \Psi_{k,m,n}(x) := \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k^2}{2}|x|} \sqrt{\frac{\pi}{km}} e^{ik(\sqrt{\pi/k}n) \cdot x} $$
We can scale the family with the change of variable $x = k^{1/2}y$, leading to

$$ \Psi_{k,m,n}(x) := \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k}{2}|x| - \sqrt{\pi/km}|x|^2} e^{ik(\sqrt{\pi/kn}) \cdot x} $$

Introducing the lattice $x^m := \sqrt{\pi/km}$ and $\xi^n := \sqrt{\pi/kn}$, we have

$$ \Psi_{k,m,n}(x) = \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k}{2}|x-x^m|^2} e^{ik\xi^n \cdot x} $$

<figure><img src="image_14.png" /></figure>
We can scale the family with the change of variable $x = k^{1/2}y$, leading to

$$ \Psi_{k,m,n}(x) := \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k}{2}|x| - \sqrt{\pi/km}|x|^2} e^{ik(\sqrt{\pi/kn}) \cdot x} $$

Introducing the lattice $x^m := \sqrt{\pi/km}$ and $\xi^n := \sqrt{\pi/kn}$, we have

$$ \Psi_{k,m,n}(x) = \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k}{2}|x-x^m|^2} e^{ik\xi^n \cdot x} $$

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

This is still a dense set in $L^2(\mathbb{R}^d)$.
We can scale the family with the change of variable $x = k^{1/2}y$, leading to

$$ \Psi_{k,m,n}(x) := \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k}{2}|x| - \sqrt{\pi/km}|x|^2} e^{ik(\sqrt{\pi/kn}) \cdot x} $$

Introducing the lattice $x^m := \sqrt{\pi/km}$ and $\xi^n := \sqrt{\pi/kn}$, we have

$$ \Psi_{k,m,n}(x) = \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k}{2}|x-x^m|^2} e^{ik\xi^n \cdot x} $$

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

This is still a dense set in $L^2(\mathbb{R}^d)$.

We have a grid $[x^m, \xi^n]$ on $\mathbb{R}^{2d}$ instead of $\mathbb{R}^d$, spaced by $k^{-1/2}$ instead of $k^{-1}$.
The family {$\Psi_{k,m,n}$}$_{[m,n] \in \mathbb{Z}^{2d}}$ is dense in $L^2(\mathbb{R}^d)$.
The family {$\Psi_{k,m,n}$}$_{[m,n]\in\mathbb{Z}^{2d}}$ is dense in $L^2(\mathbb{R}^d)$.

It also forms a *frame* on $L^2(\mathbb{R}^d)$. It means that

$$ \alpha \|v\|_{L^2(\mathbb{R}^d)}^2 \le \sum_{[m,n] \in \mathbb{Z}^{2d}} |(v, \psi_{k,m,n})_{L^2(\mathbb{R}^d)}|^2 \le \beta \|v\|_{L^2(\mathbb{R}^d)}^2 \quad \forall v \in L^2(\mathbb{R}^d), $$

with $0 < \alpha \le \beta < +\infty$, two constants independent of $k$.
The family {$\Psi_{k,m,n}$}$_{[m,n]\in\mathbb{Z}^{2d}}$ is dense in $L^2(\mathbb{R}^d)$.

It also forms a *frame* on $L^2(\mathbb{R}^d)$. It means that

$$ \alpha \|v\|_{L^2(\mathbb{R}^d)}^2 \le \sum_{[m,n] \in \mathbb{Z}^{2d}} |(v, \Psi_{k,m,n})_{L^2(\mathbb{R}^d)}|^2 \le \beta \|v\|_{L^2(\mathbb{R}^d)}^2 \quad \forall v \in L^2(\mathbb{R}^d), $$

with $0 < \alpha \le \beta < +\infty$, two constants independent of $k$.

This property is a generalization of the concept of Hilbertian basis.

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

Daubechies, Grossman and Meyer, 1986.
Convergence for general right-hand sides

Approximating smooth functions in the Gabor frame
Think about $L^2(-\pi, \pi)$. Then the family $\{e^{inx}\}_{n\in\mathbb{Z}}$ is a Hilbertian basis.
Think about $L^2(-\pi, \pi)$. Then the family $\{e^{inx}\}_{n \in \mathbb{Z}}$ is a Hilbertian basis.

In particular, we have

$$ v = \sum_{n \in \mathbb{Z}} (v, e^{inx})_{L^2(-\pi, \pi)} e^{inx} \quad \forall v \in L^2(-\pi, \pi). $$
Think about $L^2(-\pi, \pi)$. Then the family $\{e^{inx}\}_{n \in \mathbb{Z}}$ is a Hilbertian basis.

In particular, we have

$$ v = \sum_{n \in \mathbb{Z}} (v, e^{inx})_{L^2(-\pi, \pi)} e^{inx} \quad \forall v \in L^2(-\pi, \pi). $$

If $v \in H^m(-\pi, \pi)$, we can actually show that

$$ |(v, e^{inx})_{L^2(-\pi, \pi)}| \le \frac{C}{|n|^m} |(v^{(m)}, e^{inx})_{L^2(-\pi, \pi)}| $$
Think about $L^2(-\pi, \pi)$. Then the family $\{e^{inx}\}_{n \in \mathbb{Z}}$ is a Hilbertian basis.

In particular, we have

$$v = \sum_{n \in \mathbb{Z}} (v, e^{inx})_{L^2(-\pi, \pi)} e^{inx} \quad \forall v \in L^2(-\pi, \pi).$$

If $v \in H^m(-\pi, \pi)$, we can actually show that

$$|(v, e^{inx})_{L^2(-\pi, \pi)}| \le \frac{C}{|n|^m} |(v^{(m)}, e^{inx})_{L^2(-\pi, \pi)}|$$

An easy consequence is that

$$\|v - v_N\|_{L^2(-\pi, \pi)} \le C N^{-m} |v|_{H^m(-\pi, \pi)}, \quad v_N := \sum_{\substack{n \in \mathbb{Z} \\ |n| < N}} (v, e^{inx})_{L^2(-\pi, \pi)} e^{inx}.$$
A general result on frames is the existence of a dual frame $\{\Psi_{k,m,n}^*\}_{[m,n]\in\mathbb{Z}^{2d}}$ s.t.

$$ v = \sum_{[m,n]\in\mathbb{Z}^{2d}} (v, \Psi_{k,m,n}^*)_{L^2(\mathbb{R}^d)} \Psi_{k,m,n}. $$

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

Christensen, 2003.
A general result on frames is the existence of a dual frame $\{\Psi_{k,m,n}^*\}_{[m,n]\in\mathbb{Z}^{2d}}$ s.t.

$$ v = \sum_{[m,n]\in\mathbb{Z}^{2d}} (v, \Psi_{k,m,n}^*)_{L^2(\mathbb{R}^d)} \Psi_{k,m,n}. $$

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

Christensen, 2003.

Our first contribution is to show that

$$ \sum_{\substack{[m,n]\in\mathbb{Z}^{2d} \\ |x^m|^2+|\xi^n|^2\ge D^2}} |(v, \Psi_{k,m,n}^*)_{L^2(\mathbb{R}^d)}|^2 \le C D^{-2m} \|v\|_{\dot{H}_k^m(\mathbb{R}^d)}^2 \quad \forall v \in \dot{H}^m(\mathbb{R}^d) $$

where $\dot{H}^m(\mathbb{R}^d)$ is a weighted Sobolev space:

$$ \|v\|_{\dot{H}_k^m(\mathbb{R}^d)}^2 := \sum_{[a] \le m} \sum_{q=0}^{p-[a]} k^{-2[a]} \|x\|^q \partial^q v \|_{L^2(\mathbb{R}^d)}^2. $$

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

Chaumont-Frelet and Ingremeau, 2022.
For $v \in L^2(\mathbb{R}^d)$ let us set

$$ \Pi_D v := \sum_{\substack{[m,n] \in \mathbb{Z}^{2d} \\ |x^m|^2 + |\xi^n|^2 < D^2}} (v, \Psi_{k,m,n}^*) \Psi_{k,m,n}. $$

<figure><img src="image_21.png" /></figure>
For $v \in L^2(\mathbb{R}^d)$ let us set

$$ \Pi_D v := \sum_{|m,n| \in \mathbb{Z}^{2d} \atop |x^m|^2 + |\xi^n|^2 < D^2} (v, \Psi_{k,m,n}^*) \Psi_{k,m,n}. $$

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

Then, we have

$$ \|v - \Pi_D v\|_{\dot{H}^p(\mathbb{R}^d)} \leq C D^{-r} \|v\|_{\dot{H}^{p+r}(\mathbb{R}^d)} $$

for all $v \in \dot{H}^{p+r}(\mathbb{R}^d)$ and $0 \le p \le r$.

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

Chaumont-Frelet and Ingremeau, 2022.
Convergence for general right-hand sides

Application to the Helmholtz equation
By combining standard arguments and decay in the PML, we can show that

$$ \|u\|_{H_k^2(\mathbb{R}^d)} \le C \|f\|_{L^2(\mathbb{R}^d)}. $$
By combining standard arguments and decay in the PML, we can show that

$$ \|u\|_{H_k^2(\mathbb{R}^d)} \le C \|f\|_{L^2(\mathbb{R}^d)}. $$

An easy consequence is

$$ \|u - \Pi_D u\|_{H_k^1(\mathbb{R}^d)} \le C D^{-1} \|f\|_{L^2(\mathbb{R}^d)} \quad \forall k \ge 1. $$
Let

$$W_{\Lambda} := \mathrm{Vect} \{ \Psi_{k,m,n}; \ [m, n] \in \Lambda \}$$

where

$$\Lambda := \{ [m, n] \in \mathbb{Z}^{2d}; \ |x^m|^2 + |\xi^n|^2 \le D^2 \}.$$

<figure><img src="image_24.png" /></figure>
Let

$$
W_{\Lambda} := \mathrm{Vect} \{ \Psi_{k,m,n}; [\mathbf{m}, \mathbf{n}] \in \Lambda \}
$$

where

$$
\Lambda := \left\{ [\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d} ; \; |\mathbf{x}^m|^2 + |\xi^n|^2 \le D^2 \right\}.
$$

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

Then, we have

$$
\dim W_{\Lambda} \sim D^{2d} k^d, \quad \text{and} \quad \min_{w \in W_{\Lambda}} \|u - w\|_{H_k^1(\mathbb{R}^d)} \le C D^{-1} \|f\|_{L^2(\mathbb{R}^d)}.
$$
Let

$$
W_{\Lambda} := \operatorname{Vect} \{ \Psi_{k,m,n}; [\mathbf{m}, \mathbf{n}] \in \Lambda \}
$$

where

$$
\Lambda := \left\{ [\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d} ; \; |\mathbf{x}^m|^2 + |\boldsymbol{\xi}^n|^2 \le D^2 \right\}.
$$

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

Then, we have

$$
\dim W_{\Lambda} \sim D^{2d} k^d, \text{ and } \min_{w \in W_{\Lambda}} \|u - w\|_{H_k^1(\mathbb{R}^d)} \le C D^{-1} \|f\|_{L^2(\mathbb{R}^d)}.
$$

For general $L^2(\mathbb{R}^d)$ rhs, same number of dofs than standard FEM: $k^d$.
Let

$$
W_{\Lambda} := \mathrm{Vect} \{ \Psi_{k,m,n}; [\mathbf{m}, \mathbf{n}] \in \Lambda \}
$$

where

$$
\Lambda := \left\{ [\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d} ; \; |\mathbf{x}^m|^2 + |\boldsymbol{\xi}^n|^2 \le D^2 \right\}.
$$

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

Then, we have

$$
\dim W_{\Lambda} \sim D^{2d} k^d, \quad \text{and} \quad \min_{w \in W_{\Lambda}} \|u-w\|_{H_k^1(\mathbb{R}^d)} \le C D^{-1} \|f\|_{L^2(\mathbb{R}^d)}.
$$

For general $L^2(\mathbb{R}^d)$ rhs, same number of dofs than standard FEM: $k^d$.

This is expected: for a general $f \in L^2(\mathbb{R}^d)$, $u$ is a general $H^2(\mathbb{R}^d)$ function.
Convergence for scattering problems
Convergence for scattering problems

Semi-classical concepts: symbols and micro-localization
Assuming for simplicity that **A** is constant, we have:

$$ \{-k^2 \mu u + \nabla \cdot (\mathbf{A}\nabla u)\} (x) = C \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} e^{ik(x-y)\cdot\xi} p(x, \xi) u(y) dy d\xi, $$

where $p : \mathbb{R}^{2d} \to \mathbb{C}$ with

$$ p(x, \xi) := \mathbf{A}(x)\xi \cdot \xi - \mu(x) $$

is the symbol associated with the Helmholtz differential operator.
Assuming for simplicity that **A** is constant, we have:

$$ \{-k^2 \mu u + \nabla \cdot (\mathbf{A}\nabla u)\} (x) = C \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} e^{ik(x-y)\cdot\xi} p(x, \xi) u(y) dy d\xi, $$

where $p: \mathbb{R}^{2d} \to \mathbb{C}$ with

$$ p(x, \xi) := \mathbf{A}(x)\xi \cdot \xi - \mu(x) $$

is the symbol associated with the Helmholtz differential operator.

We can see the action of operators through a function defined over $\mathbb{R}^{2d}$.
Assuming for simplicity that **A** is constant, we have:

$$ \{-k^2 \mu u + \nabla \cdot (\mathbf{A}\nabla u)\} (x) = C \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} e^{ik(x-y)\cdot\xi} p(x, \xi) u(y) dy d\xi, $$

where $p: \mathbb{R}^{2d} \to \mathbb{C}$ with

$$ p(x, \xi) := \mathbf{A}(x)\xi \cdot \xi - \mu(x) $$

is the symbol associated with the Helmholtz differential operator.

We can see the action of operators through a function defined over $\mathbb{R}^{2d}$.

We can define the micro-local support of a function $v \in L^2(\mathbb{R}^d)$ by duality

$$ \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} e^{ik(x-y)\cdot\xi} a(x, \xi) v(y) dy d\xi $$

as a subset of $\mathbb{R}^{2d}$.
Convergence for scattering problems

Micro-localization of high-frequency solutions
We can show that

$$ (a^W u, u)_{L^2(\mathbb{R}^d)} = \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} e^{ik(x-y)\cdot\xi} a\left(\frac{x+y}{2}, \xi\right) u(y) \bar{u}(x) dy d\xi dx \\ \xrightarrow{k \to \infty} \\ \int_{\mathbb{R}^{2d}} a(x, \xi) d\lambda(x, \xi) $$

for all "smooth" symbols a, where $\lambda$ is a measure over $\mathbb{R}^{2d}$.
We can show that

$$ (a^W u, u)_{L^2(\mathbb{R}^d)} = \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} e^{ik(x-y)\cdot\xi} a\left(\frac{x+y}{2}, \xi\right) u(y) \bar{u}(x) dy d\xi dx \\ \xrightarrow{k \to \infty} \\ \int_{\mathbb{R}^{2d}} a(x, \xi) d\lambda(x, \xi) $$

for all "smooth" symbols $a$, where $\lambda$ is a measure over $\mathbb{R}^{2d}$.

It is called the "defect measure", and we have

$$ \text{supp } \lambda \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\} \cup \text{supp } f. $$
We can show that

$$ (a^W u, u)_{L^2(\mathbb{R}^d)} = \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} e^{ik(x-y)\cdot\xi} a\left(\frac{x+y}{2}, \xi\right) u(y) \bar{u}(x) dy d\xi dx \\ \xrightarrow{k \to \infty} \\ \int_{\mathbb{R}^{2d}} a(x, \xi) d\lambda(x, \xi) $$

for all "smooth" symbols $a$, where $\lambda$ is a measure over $\mathbb{R}^{2d}$.

It is called the "defect measure", and we have

$$ \text{supp } \lambda \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\} \cup \text{supp } f. $$

More generally, $u$ is micro-localized close to $\text{supp } \lambda$ for large frequencies.
We can show that

$$ (a^W u, u)_{L^2(\mathbb{R}^d)} = \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} e^{ik(x-y)\cdot\xi} a\left(\frac{x+y}{2}, \xi\right) u(y) \bar{u}(x) dy d\xi dx \\ \xrightarrow{k \to \infty} \\ \int_{\mathbb{R}^{2d}} a(x, \xi) d\lambda(x, \xi) $$

for all "smooth" symbols $a$, where $\lambda$ is a measure over $\mathbb{R}^{2d}$.

It is called the "defect measure", and we have

$$ \text{supp } \lambda \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\} \cup \text{supp } f. $$

More generally, $u$ is micro-localized close to $\text{supp } \lambda$ for large frequencies.

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

Zworski, 2012
In an appropriate sense, we have

$$u \xrightarrow{k \to \infty} \lambda,$$

where $\lambda$ is a measure on $\mathbb{R}^{2d}$ such that

$$\text{supp } \lambda \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\} \cup \text{supp } f.$$
In an appropriate sense, we have

$$u \xrightarrow{k \to \infty} \lambda,$$

where $\lambda$ is a measure on $\mathbb{R}^{2d}$ such that

$$\text{supp } \lambda \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\} \cup \text{supp } f.$$

In general, this information is not that useful, because $\text{supp } f$ can be anything.
In an appropriate sense, we have

$$u \underset{k \to \infty}{\longrightarrow} \lambda,$$

where $\lambda$ is a measure on $\mathbb{R}^{2d}$ such that

$$\text{supp } \lambda \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\} \cup \text{supp } f.$$

In general, this information is not that useful, because $\text{supp } f$ can be anything.

For scattering problems, however, $\text{supp } f \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\}$, and

$$\text{supp } \lambda \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\}.$$
In an appropriate sense, we have

$$u \xrightarrow{k \to \infty} \lambda,$$

where $\lambda$ is a measure on $\mathbb{R}^{2d}$ such that

$$\text{supp } \lambda \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\} \cup \text{supp } f.$$

In general, this information is not that useful, because $\text{supp } f$ can be anything.

For scattering problems, however, $\text{supp } f \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\}$, and

$$\text{supp } \lambda \subset \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\}.$$

In the high-frequency regime, we know "where" to look for $u$!
<figure><img src="image_29.png" /></figure>

$$p(x, \xi) = 0 \iff \mu(x) - \alpha(x)\xi^2 = 0 \iff \xi^2 = \vartheta(x)^2$$
<figure><img src="image_30.png" /></figure>
Convergence for scattering problems

Micro-localization properties of Gaussian states
Consider the Gaussian distribution

$$ \Psi(x) := \pi^{-d/4} e^{-\frac{1}{2}|x|^2}. $$

It is localized around 0 in position.
Consider the Gaussian distribution

$$ \Psi(x) := \pi^{-d/4} e^{-\frac{1}{2}|x|^2}. $$

It is localized around 0 in position. But we also have

$$ \hat{\Psi}(\xi) := \int_{\mathbb{R}^d} \Psi(x) e^{i\xi \cdot x} dx = (2\pi)^{d/2} \pi^{-d/4} e^{-\frac{1}{2}|\xi|^2}, $$

so that it is also localized around 0 in momentum.
Consider the Gaussian distribution

$$ \Psi(x) := \pi^{-d/4} e^{-\frac{1}{2}|x|^2}. $$

It is localized around 0 in position. But we also have

$$ \hat{\Psi}(\xi) := \int_{\mathbb{R}^d} \Psi(x) e^{i\xi \cdot x} dx = (2\pi)^{d/2} \pi^{-d/4} e^{-\frac{1}{2}|\xi|^2}, $$

so that it is also localized around 0 in momentum.

Actually, it has the “best” micro-localization property, as it saturates Heisenberg’s uncertainty principle:

$$ \|x\Psi\|_{L^2(\mathbb{R}^d)} \|\xi\hat{\Psi}\|_{L^2(\mathbb{R}^d)} \underset{(\geq)}{=} \frac{1}{2} \|\Psi\|_{L^2(\mathbb{R}^d)}^2 \|\hat{\Psi}\|_{L^2(\mathbb{R}^d)}^2. $$
Recall that our discretization space is based on the *k*-scaled version

$$ \Psi_{k,m,n}(x) := \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k}{2}|x-x^m|^2} e^{ikx \cdot \xi^n}, $$

where $[x^m, \xi^n] := \sqrt{\pi/k}[m, n]$.
Recall that our discretization space is based on the *k*-scaled version

$$ \Psi_{k,m,n}(x) := \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k}{2}|x-x^m|^2} e^{ikx \cdot \xi^n}, $$

$$ \text{where } [x^m, \xi^n] := \sqrt{\pi/k}[m, n]. $$

The scaling is adapted to the scaled Fourier transform employed in symbols.
Recall that our discretization space is based on the *k*-scaled version

$$ \Psi_{k,m,n}(x) := \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k^2}{2}|x-x^m|^2} e^{ikx \cdot \xi^n}, $$

$$ \text{where } [x^m, \xi^n] := \sqrt{\pi/k}[m, n]. $$

The scaling is adapted to the scaled Fourier transform employed in symbols.

In fact, we have

$$ \Psi_{k,x_0,\xi_0}(x) := \left(\frac{k}{\pi}\right)^{d/4} e^{-\frac{k^2}{2}|x-x_0|^2} e^{ikx \cdot \xi_0} \underset{k \to \infty}{\longrightarrow} \delta_{[x_0, \xi_0]}, $$

in the sense that

$$ (a^W \Psi_{k,x_0,\xi_0}, \Psi_{k,x_0,\xi_0})_{L^2(\mathbb{R}^d)} \underset{k \to \infty}{\longrightarrow} a(x_0, \xi_0) $$

for all "smooth" symbols a.
<figure><img src="image_31.png" /></figure>
<figure><img src="image_32.png" /></figure>
<figure><img src="image_33.png" /></figure>
<figure><img src="image_34.png" /></figure>
Convergence for scattering problems

Approximation of high-frequency solutions
Fix $\varepsilon > 0$, and consider the set of indices

$$ \Lambda := \{[\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d} \mid |p(x^\mathbf{m}, \xi^\mathbf{n})| \le k^{-1/2+\varepsilon}\}, $$

and the associated space

$$ W_{\Lambda} := \mathrm{Vect}\{\Psi_{k,m,n}; \ [\mathbf{m}, \mathbf{n}] \in \Lambda\}. $$
Fix $\varepsilon > 0$, and consider the set of indices

$$ \Lambda := \{[\mathbf{m}, \mathbf{n}] \in \mathbb{Z}^{2d} \mid |p(x^m, \xi^n)| \le k^{-1/2+\varepsilon}\}, $$

and the associated space

$$ W_{\Lambda} := \mathrm{Vect} \{ \Psi_{k,m,n}; [\mathbf{m}, \mathbf{n}] \in \Lambda \}. $$

Then, we have

$$ \dim W_{\Lambda} \sim k^{d-1/2+\varepsilon}, $$

and

$$ \min_{w \in W_{\Lambda}} \|u - w\|_{H_k^1(\Omega)} \le C_{\varepsilon,m} k^{-m} \quad \forall m \in \mathbb{N}. $$

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

T. Chaumont-Frelet, V. Dolean and M. Ingremau, 2022
Convergence for scattering problems
Takeaways
We consider an hypersurface in $\mathbb{R}^{2d}$. Specifically, the $2d-1$ manifold

$$ \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\}. $$
We consider an hypersurface in $\mathbb{R}^{2d}$. Specifically, the $2d-1$ manifold

$$ \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\}. $$

Since $\mathbb{R}^{2d}$ is paved with a lattice spaced in $k^{-1/2}$, we have about

$$ (k^{1/2})^{2d-1} = k^{d-1/2} $$

dofs in our discretization space.
We consider an hypersurface in $\mathbb{R}^{2d}$. Specifically, the $2d-1$ manifold

$$ \{[x, \xi] \in \mathbb{R}^{2d}; p(x, \xi) = 0\}. $$

Since $\mathbb{R}^{2d}$ is paved with a lattice spaced in $k^{-1/2}$, we have about

$$ (k^{1/2})^{2d-1} = k^{d-1/2} $$

dofs in our discretization space.

Micro-localization properties of high-frequency solutions and of Gaussian states are the key to allow this “sparse” selection of shape functions.
Numerical examples
Numerical examples

A least-squares formulation
Define $P_k: H^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)$ as $P_k v := -k^2 \mu v - \nabla \cdot (A \nabla v)$.
Define $P_k: H^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)$ as $P_k v := -k^2 \mu v - \nabla \cdot (A \nabla v)$.

The least-squares formulation is to find $u_\Lambda \in W_\Lambda$ such that

$$ (P_k u_\Lambda, P_k v)_{L^2(\mathbb{R}^d)} = (f, P_k v)_{L^2(\mathbb{R}^d)} $$

for all $v \in W_\Lambda$.
Define $P_k: H^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)$ as $P_k v := -k^2 \mu v - \nabla \cdot (A \nabla v)$.

The least-squares formulation is to find $u_\Lambda \in W_\Lambda$ such that

$$ (P_k u_\Lambda, P_k v)_{L^2(\mathbb{R}^d)} = (f, P_k v)_{L^2(\mathbb{R}^d)} $$

for all $v \in W_\Lambda$.

We can show that there exists a unique discrete solution, and that

$$ \|u - u_\Lambda\|_{H_k^1(\mathbb{R}^d)} \le C_{\varepsilon,m} k^{-m} \quad \forall m \in \mathbb{N}. $$
Define $P_k: H^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)$ as $P_k v := -k^2 \mu v - \nabla \cdot (A \nabla v)$.

The least-squares formulation is to find $u_\Lambda \in W_\Lambda$ such that

$$ (P_k u_\Lambda, P_k v)_{L^2(\mathbb{R}^d)} = (f, P_k v)_{L^2(\mathbb{R}^d)} $$

for all $v \in W_\Lambda$.

We can show that there exists a unique discrete solution, and that

$$ \|u - u_\Lambda\|_{H_k^1(\mathbb{R}^d)} \le C_{\varepsilon,m} k^{-m} \quad \forall m \in \mathbb{N}. $$

In addition, we can also show that

$$ |(P_k \Psi_{k,m,n}, P_k \Psi_{k,m',n'})|_{L^2(\mathbb{R}^d)} \le C_m (|m - m'| + |n - n'|)^{-m} \quad \forall m \in \mathbb{N}, $$

so that the matrix is essentially sparse.
We solve the equation

$$-k^2 \mu \nu u - (\nu^{-1} u')' = f,$$

where $\nu$ is the PML scaling function.
We solve the equation

$$ -k^2 \mu \nu u - (\nu^{-1} u')' = f, $$

where $\nu$ is the PML scaling function.

We only consider scattering problems, hence

$$ f(x) := \phi e^{ikx}, $$

with a smooth $\phi$ supported where $\mu = 1$.
We solve the equation

$$ -k^2 \mu \nu u - (\nu^{-1} u')' = f, $$

where $\nu$ is the PML scaling function.

We only consider scattering problems, hence

$$ f(x) := \phi e^{ikx}, $$

with a smooth $\phi$ supported where $\mu = 1$.

In the first example, $\mu \equiv 1$.

In the second $\mu = 2$ on $(-1/2, 1/2)$ and $\mu = 1$ outside $(-1, 1)$.
We solve the equation

$$ -k^2 \mu \nu u - (\nu^{-1} u')' = f, $$

where $\nu$ is the PML scaling function.

We only consider scattering problems, hence

$$ f(x) := \phi e^{ikx}, $$

with a smooth $\phi$ supported where $\mu = 1$.

In the first example, $\mu \equiv 1$.

In the second $\mu = 2$ on $(-1/2, 1/2)$ and $\mu = 1$ outside $(-1, 1)$.

We solve the problem for different values of $k$ and $\delta$ where

$$ \Lambda := \{ [x^m, \xi^n] \in \mathbb{R}^2; |p(x^m, \xi^n)| \le \delta \}. $$
Numerical examples

Homogeneous test case
<figure><img src="image_36.png" /></figure>
8

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

Computational cost for a 10⁻⁵ relative error.
Numerical examples

Heterogeneous test case
<figure><img src="image_38.png" /></figure>
<figure><img src="image_39.png" /></figure>

Computational cost for a $10^{-5}$ relative error
Concluding remarks
Concluding remarks

Open questions and challenges
We currently only propose approximability estimates,
which we only know how to apply to least-squares formulations.
We currently only propose approximability estimates,
which we only know how to apply to least-squares formulations.
How can we cover the convergence of Galerkin discretizations?
We currently only propose approximability estimates, which we only know how to apply to least-squares formulations. How can we cover the convergence of Galerkin discretizations?

I lied: the matrices are *not* sparse!
But... The entries decay super-algebraically away from the diagonal.
We currently only propose approximability estimates, which we only know how to apply to least-squares formulations. How can we cover the convergence of Galerkin discretizations?

I lied: the matrices are *not* sparse!
But... The entries decay super-algebraically away from the diagonal. How should we properly sparsify and/or compress the matrices?
We currently only propose approximability estimates, which we only know how to apply to least-squares formulations. How can we cover the convergence of Galerkin discretizations?

I lied: the matrices are *not* sparse!
But... The entries decay super-algebraically away from the diagonal. How should we properly sparsify and/or compress the matrices?

Actually, even the computation of the entries is challenging.
We currently only propose approximability estimates, which we only know how to apply to least-squares formulations. How can we cover the convergence of Galerkin discretizations?

I lied: the matrices are *not* sparse!
But... The entries decay super-algebraically away from the diagonal. How should we properly sparsify and/or compress the matrices?

Actually, even the computation of the entries is challenging. How can we derive suited quadrature rules for Gabor wavelets?
We currently only propose approximability estimates, which we only know how to apply to least-squares formulations. How can we cover the convergence of Galerkin discretizations?

I lied: the matrices are *not* sparse!
But... The entries decay super-algebraically away from the diagonal. How should we properly sparsify and/or compress the matrices?

Actually, even the computation of the entries is challenging. How can we derive suited quadrature rules for Gabor wavelets?

The frame properties of Gabor wavelets are established in the whole $\mathbb{R}^d$. This enables us to treat (globally) smooth coefficients (including PML).
We currently only propose approximability estimates, which we only know how to apply to least-squares formulations. How can we cover the convergence of Galerkin discretizations?

I lied: the matrices are *not* sparse!
But... The entries decay super-algebraically away from the diagonal. How should we properly sparsify and/or compress the matrices?

Actually, even the computation of the entries is challenging. How can we derive suited quadrature rules for Gabor wavelets?

The frame properties of Gabor wavelets are established in the whole $\mathbb{R}^d$. This enables us to treat (globally) smooth coefficients (including PML). How can we incorporate boundary and/or interface conditions?
Concluding remarks

Main takeaways
Gabor wavelets form a frame over $L^2(\mathbb{R}^d)$.

They actually correspond to a $k^{-1/2}$ lattice of $\mathbb{R}^{2d}$.
Gabor wavelets form a frame over $L^2(\mathbb{R}^d)$.

They actually correspond to a $k^{-1/2}$ lattice of $\mathbb{R}^{2d}$.

For generic right-hand sides, they require $k^d$ dofs.
Gabor wavelets form a frame over $L^2(\mathbb{R}^d)$.
They actually correspond to a $k^{-1/2}$ lattice of $\mathbb{R}^{2d}$.

For generic right-hand sides, they require $k^d$ dofs.

A key aspect of Gaussian states is that they are localized in phase-space.
For scattering problems, we only need to "pave" the hypersurface $p(x, \xi) = 0$.
Gabor wavelets form a frame over $L^2(\mathbb{R}^d)$.
They actually correspond to a $k^{-1/2}$ lattice of $\mathbb{R}^{2d}$.

For generic right-hand sides, they require $k^d$ dofs.

A key aspect of Gaussian states is that they are localized in phase-space.
For scattering problems, we only need to "pave" the hypersurface $p(x, \xi) = 0$.

For scattering problems, we only need $k^{d-1/2}$ dofs with sparse matrices.