Stable approximation of Helmholtz solutions in the disk by evanescent plane waves

Superpositions of plane waves are known to approximate well the solutions of the Helmholtz equation. Their use in discretizations is typical of Trefftz methods for Helmholtz problems, aiming to achieve high accuracy with a small number of degrees of freedom. However, Trefftz methods lead to ill-conditioned linear systems, and it is often impossible to obtain the desired accuracy in floating-point arithmetic. In this paper we show that a judicious choice of plane waves can ensure high-accuracy solutions in a numerically stable way, in spite of having to solve such ill-conditioned systems. Numerical accuracy of plane wave methods is linked not only to the approximation space, but also to the size of the coefficients in the plane wave expansion. We show that the use of plane waves can lead to exponentially large coefficients, regardless of the orientations and the number of plane waves, and this causes numerical instability. We prove that all Helmholtz fields are continuous superposition of evanescent plane waves, i.e., plane waves with complex propagation vectors associated with exponential decay, and show that this leads to bounded representations. We provide a constructive scheme to select a set of real and complex-valued propagation vectors numerically. This results in an explicit selection of plane waves and an associated Trefftz method that achieves accuracy and stability. The theoretical analysis is provided for a two-dimensional domain with circular shape. However, the principles are general and we conclude the paper with a numerical experiment demonstrating practical applicability also for polygonal domains.


Introduction
The space dependence of time-harmonic solutions U (x, t) = ℜ{e −ıωt u(x)} of the scalar wave equation 1 The wavenumber is κ := ω/c > 0, with c the wave speed and ω the time frequency.Solutions of boundary value problems for the Helmholtz equation are oscillatory, making their numerical approximation notoriously computationally expensive at high frequencies, namely when the wavelength λ := 2π/κ is much smaller than the characteristic length of the domain.
A well-studied way to efficiently represent Helmholtz solutions in a domain of R n is to approximate them with linear combinations of propagative plane waves x → e ıκd•x , which are particular solutions of (1.1) if the propagation direction d ∈ R n satisfies d • d = 1.Plane waves indeed offer better accuracy for fewer degrees of freedom compared to polynomial spaces, as supported by the theory developed in [30], building on previous results in [28,Sec. 8.4] and [9,Sec. 3.3.5].Approximation by plane waves has been extensively used in the context of Trefftz schemes for the Helmholtz equation, a class of methods that use trial and test functions satisfying (1.1) locally on each element of a mesh, see [22] for a comprehensive survey.The simple expression of plane waves allows for very cheap implementations; for instance, integrals of products of these functions can be computed in closed form with wavenumber-independent effort, see [22,Sec. 4.1].A second widespread use of plane wave approximation is the reconstruction of sound fields from point measurements (representing microphones) in experimental acoustics, see [11,25,34,20].
The computation of plane wave approximations is however known to be numerically unstable, imposing strong limits to the achievable accuracy [33,6].This issue is often understood as an effect of the ill-conditioning of the linear system that is solved [22,Sec. 4.3], which inevitably arises from the almost-linear dependence of plane waves with similar propagation directions.Different techniques have been proposed to overcome this instability, e.g.[3,15,6].A wellknown recommendation suggests using not more than a prescribed number of waves in elements of a given size, e.g.[23,Eq. (14)]: this keeps the instability at bay but limits the achievable accuracy.
The first purpose of this paper is to shed a new light on the numerical instability experienced with propagative plane waves and explain the fundamental mathematical reasons for their limitations as described above.The second objective is to propose a practical remedy, in the form of including evanescent plane waves, which may decay exponentially in one direction, and using which one can achieve arbitrary accuracy in a numerically stable way.The approach is substantiated by theoretical analysis in combination with numerical evidence.As a first step in this direction, we focus mainly on the model approximation problem of Helmholtz solutions in the unit disk, using the modal analysis tools described in Section 2.
A new point of view on plane wave instability.Recent advances in approximation theory, in particular based on the theory of frames and overcomplete bases [12], have shown that in the presence of ill-conditioning it is not sufficient to study best approximation errors in order to obtain accurate results in floating-point arithmetic [1,2].Rather, one is led to study the approximation error in relation to the coefficient norm, i.e., the norm of the coefficients in the expansion.The former depends solely on the approximation space, but the coefficient norm also depends on its chosen representation (i.e. on the spanning set used).We formalize this in Section 3 with the notion of stable approximation in Definition 3.1.The corresponding error analysis in Section 3.4, based largely on results in [1,2], allows us to conclude in Section 4 that the set of propagative plane waves does not yield stable approximations.That is, we can formally state that there exist Helmholtz solutions, with relatively high Fourier frequency components in the angular coordinate, that are well approximated in the approximation space, but are nevertheless not numerically computable, see Theorem 4.3.In the terminology of approximation theory, no countable set of propagative plane waves is a frame for the space of Helmholtz solutions.(We recall that a frame of a Hilbert space is a natural generalization of a basis that allows for redundancy, see [12,1].)In particular, it lacks a so-called lower frame bound which is the property that ensures that bounded functions can be represented with bounded coefficients.This point of view is reminiscent of a similar work in the context of the Method of Fundamental Solutions [5], which pre-dates the stability analysis from frame theory.
Unfortunately, while the theory in [1,2] allows to identify this problem, it offers no concrete suggestions as to how it can be remedied.If the approximation space remains unchanged, a lower frame bound can only be established through a change of basis, such as orthogonalization as in [3,15,8].However, that changes the representation: the solution would no longer be represented in the simple form of an expansion in plane waves, which is a key feature we would like to retain.Moreover, it may not be straightforward to ensure that the orthogonalization process itself is numerically stable.
The evanescent plane wave remedy.To obtain stable representations, we propose to enrich the approximation space with evanescent plane waves, i.e. plane waves whose direction vectors are complex-valued, d ∈ C n , as defined in Section 5.The Helmholtz equation is still satisfied provided d • d = 1 and, importantly, the expression remains simple and cheap to use in numerical schemes.Since their modulus decays exponentially in the direction ℑ[d], evanescent plane waves are localized in bounded physical domains but also in the Fourier domain, hence are natural candidates for the approximation of the high frequency Fourier content exhibited by certain Helmholtz solutions.This idea is already present in the Wave Based Method, a special class of Trefftz schemes, see e.g.[17] for a survey.Evanescent plane waves also proved particularly effective in the approximation of interface problems in Trefftz methods, e.g.[27], and the approximation of integral kernels in some versions of the Fast Multipole Method [10].
To support the use of evanescent plane wave, we prove in Section 6 our main positive result, Theorem 6.7, which states that any Helmholtz solution in the unit disk can be uniquely represented in the form of a continuous superposition of evanescent plane waves.This integral representation has the key property of being stable, i.e. it features a provably bounded density (in a weighted L 2 space), and it can be seen as a generalization of the classical Herglotz representation, see e.g.[14,35].This result implies that evanescent plane waves form a continuous frame for the space of Helmholtz solutions, see Theorem 6.10.While this is stated at the continuous level, such a property paves the way for successful stable discrete expansions.Indeed, from the stability of the representation one may expect that discretizations exist with bounded coefficient norms, thereby solving the main issue with propagative plane waves.
A practical numerical recipe.In view of practical implementations, we investigate the non-trivial task of identifying suitable sets of evanescent plane waves which deliver controllable accuracy in combination with stability.A heuristic choice for a set of complex directions d is suggested in [17,Sec. 3.2] (see also [22,Sec. 3.2]), but no mathematical justification is provided.
A first idea to obtain stable discrete representations (i.e. with bounded coefficients) would be to discretize the continous frame, but unfortunately, our setting does not fall within the assumptions of existing results (e.g. the boundedness assumption of [18,Th. 1.3] is not satisfied).Instead, the construction of approximation sets described in Section 7 is largely based on the optimal sampling procedure for weighted least-squares recently described by Cohen and Migliorati [13] (see also [21]) and subsequently used in [29], and it is illustrated with numerical experiments in Section 8.The strategy employed can be interpreted as the construction of a quadrature rule in the two-dimensional unbounded parametric domain of the integral representation.In practice, the recipe consists in drawing the quadrature points (i.e.select the directions of the plane waves) according to an explicit probability density function (7.3) which is a generalization to the multivariate setting of the Christoffel function density, The latter is sometimes called spectral function.While the rigorous numerical analysis of the above approach is thus far incomplete, we conjecture that such a construction provides stable discrete representations, see Conjecture 7.1.In fact, the experimental results in Section 8 show that the resulting approximations are both controllably accurate and numerically stable, provided one uses sufficient oversampling and regularization.
Although the recipe is derived from the analysis on the disk, we include numerical results on a triangular cell showing that it appears to be effective also for other shapes.Approximation and stability properties of evanescent plane waves in more general domains and their use in mesh-based Trefftz methods (e.g. the Trefftz-Discontinuous Galerkin method [22, Sec.2.2]) will be considered in future publications (see [19] for the extension of the theory of this paper to three-dimensional problems).

Helmholtz equation in circular geometry
We first present the setting of the paper and introduce notation.The proofs of the statements follow standard arguments and are collected in Appendix A.

Circular waves
In this paper we only consider circular two-dimensional geometries.Without loss of generality, we assume that the domain is the open unit disk, henceforth denoted The circular geometry enables modal analysis via separation of variables.The circular waves are the bounded solutions of the Helmholtz equation in the unit disk that are separable in polar coordinates.They are sometimes also referred to as Fourier-Bessel functions or as Generalized Harmonic Polynomials [28].The results of this paper are fomulated most concisely using the following κ-dependent scalar product and norm: for any u, v ∈ H 1 (B 1 ), In this definition, J p is the usual Bessel function of the first kind [31, Eq. (10.2.2)] and ı the imaginary unit ı 2 = −1.The space B is a strict subspace of H 1 (B 1 ), whose elements are solutions of the Helmholtz equation, see Lemma 2.3 below.A representation of the real part of some circular waves is given in Figure 1.We will refer to the circular waves with mode number |p| < κ as propagative modes.The 'energy' of such modes is distributed in the bulk of the domain.On the contrary, for |p| ≫ κ, the circular waves are termed evanescent.Their 'energy' is concentrated near the boundary of the domain.In between, the waves such that |p| ≈ κ are called grazing modes.
Lemma 2.2.The space (B, ∥ • ∥ B ) is a Hilbert space and the family {b p } p∈Z is a Hilbert basis (i.e. an orthonormal basis): The main reason for introducing circular waves is the possibility to use them to expand any Helmholtz solution on the disk, as we show in the next lemma.Related results for more general domains and different norms are available, see e.g.[22,Sec. 3 Circular and spherical waves have been used as basis functions in many Trefftz schemes, see [22,Sec. 3.1] and the references therein.An interesting feature of such waves is that the approximation sets are naturally hierarchical.

Asymptotics of normalization coefficients
The normalization coefficients β p in (2.2) grow super-exponentially with |p| after a pre-asymptotic regime up to |p| ≈ κ.The precise asymptotic behavior is given by the following lemma.

Stable numerical approximation
The purpose of this section is to explain the crucial notion of stable approximation, which we could informally call "approximation with small coefficients", and to clarify how it enables accurate numerical computations.Our approach builds on the results in [1,2] which highlight the importance for stability of having representations with bounded coefficients.We also describe the practical procedure, a regularized sampling method, that we use to investigate the existence of stable numerical approximations of Helmholtz solutions in this paper.An error bound is formulated in Proposition 3.2.

The notion of stable approximation
Let us consider a sequence of finite approximation sets in B and for each k, l, ϕ k,l ∈ B is a solution of the Helmholtz equation (1.1) in the unit disk.These sets need not be nested.When the {ϕ k,l } l are linearly independent, Φ k is a basis of the approximation space used for numerical computations.However, more generally, we also allow for linearly dependent sets.Associated to any set Φ k for some k ∈ N, we define the synthesis operator Here and in the following, we use the notation |X| to indicate the cardinality of the set X. We are now ready to define a notion of stable approximation, which is at the heart of this paper.
Definition 3.1 (Stable approximation).The sequence Φ of approximation sets (3.1) is said to be a stable approximation for B if, for any tolerance η > 0, there exist a stability exponent s ≥ 0 and a stability constant Having a sequence of stable approximation sets means that one can approximate any Helmholtz solution to a given accuracy in the form of a finite expansion T Φ k µ where the coefficients µ have bounded ℓ 2 -norm.This bound on the coefficients admits a polynomial growth in the number |Φ k | of terms in the expansion, but not an exponential growth.The stability exponent s ≥ 0 of a stable approximation sequence controls the growth of the coefficient norm ∥µ∥ ℓ 2 : the smaller s the more stable the sequence.This notion of stability is not related to a space but rather to a particular sequence of sets that are used to represent the numerical approximation.In practice the computation of approximations using stable sequences may lead to ill-conditioned linear systems if there is redundancy in the approximation sets.The rest of this section shows that, in spite of possible ill-conditioning, stable sequences lead to accurate approximations, thanks to the boundedness of the expansion coefficients.
The simplest stable approximation is provided by the truncation of any orthonormal basis of B, in which case s = 0 and C stb = 1, e.g. the circular waves Φ k = {b p } |p|≤k .However, in view of the application to Trefftz methods on polygonal meshes, we describe two examples of approximations sets of the type of (3.1): propagative plane waves (PPWs) in (4.2) and evanescent plane waves (EPWs) in (7.6).They exhibit different stability properties.In Theorem 4.3 we prove rigorously that PPWs are necessarily unstable.In contrast, numerical evidence from Section 8 indicates that the sets of EPWs constructed following the numerical recipe that we propose in Section 7.3 are stable.

Boundary sampling method
We explain how we compute the coefficients in practice, which builts on results in [24].All the numerical results obtained in this paper are obtained using the method described here.
Let us introduce a 'trace operator' γ, namely a (continuous) linear operator defined on H 1 (B 1 ) such that the following problem is well-posed: find u ∈ H 1 (B 1 ) such that We aim at reconstructing a solution u ∈ B having access to its trace γu on the boundary for such a 'good' trace operator γ.For simplicity, we use the Dirichlet trace operator and therefore assume that κ 2 is away from the eigenvalues of the Dirichlet Laplacian.Further we will assume that u ∈ B ∩ C 0 (B 1 ), so that it makes sense to consider point evaluations of the Dirichlet trace.
The reconstruction process is not the main subject of the paper and we stress that we make these two assumptions mainly for convenience and definiteness (in particular for the numerical experiments).One can consider alternative reconstruction procedures using other types of data, such as point evaluation in the bulk of the domain or by taking inner product of the solution with suitable test functions.See [11] for a more general discussion on the subject of reconstructing Helmholtz solutions from point evaluations.
Let u ∈ B ∩ C 0 (B 1 ) be the target of the approximation problem.We look for a set of coefficients ξ ∈ C |Φ k | for a given approximation set Φ k (introduced in (3.1)) such that T Φ k ξ ≈ u.We also assume that for any l, ϕ k,l ∈ B ∩ C 0 (B 1 ).Define the set of S ≥ |Φ k | sampling points {x s } S s=1 on the unit circle parametrized by the angle Let us introduce the matrix A = (A s,l ) s,l ∈ C S×|Φ k | and the vector b = (b s ) s ∈ C S such that The sampling method then consists in approximately solving the rectangular linear system Aξ = b. (3.5)

Regularization
It often happens that the matrix A is ill-conditioned (see Section 4.3).In finite precision arithmetic, severe ill-conditioning may prevent us from obtaining accurate approximations.However, the type of ill-conditioning encountered here is benign if it arises only from the redundancy of the approximating functions.In that case, ill-conditioning is associated with the numerical nonuniqueness of the solution of the linear system, yet all associated expansions may approximate the target to similar accuracy.If among those expansions there exist some with small coefficient norms, then it is possible to numerically compute an accurate approximation.To this aim, we rely on the combination of oversampling and regularization techniques developed in [1,2].Alternative techniques to curb ill-conditioning can be found in the literature, see [3] where a suitable change of basis is used that works well for circular geometries, [15] which uses orthogonalization, and [6,24] in the context of Trefftz methods.The first step is to compute the Singular Value Decomposition (SVD) of the matrix A, namely Let us denote by (σ m ) m for m = 1, . . ., |Φ k | the singular values of A, assumed to be sorted in descending order.For notational clarity, the largest singular value is renamed σ max := σ 1 .Then, the regularization amounts to trimming the tail of relatively small singular values, which are approximated by zero.Let ϵ ∈ (0, 1] be a chosen threshold, we denote by Σ ϵ the approximation of the diagonal matrix Σ where all diagonal elements σ m such that σ m < ϵσ max are replaced by zero.This leads to the approximate factorization of the matrix A. An approximate solution to (3.5) is then obtained by Here Σ † ϵ denotes the pseudo-inverse of the matrix Σ ϵ , namely the diagonal matrix with (Σ † ϵ ) j,j = (Σ j,j ) −1 if Σ j,j ≥ ϵσ max and (Σ † ϵ ) j,j = 0 otherwise.Robust computation of ξ S,ϵ requires to compute the right-hand-side of (3.7) from right to left, namely ξ S,ϵ := V Σ † ϵ (U * b) , in order to avoid mixing small and large values on the diagonal of Σ † ϵ .

Error estimates for the sampling method with regularization
With the regularization technique described above together with oversampling, i.e., S larger than |Φ k |, accurate approximations can be effectively computed, provided the set sequence is a stable approximation in the sense of Definition 3.1.This broad statement is the main message of [1, Th. 5.3] and [2, Th. 1.3 and 3.7], and is the starting point of our quest of stable approximation sets for Helmholtz solutions.More precisely, the following proposition is a rewording of [2, Th. 3.7] from the context of generalized sampling to our setting, with the notations just introduced.See Appendix B for the proof.
Proposition 3.2.Let γ be the Dirichlet trace operator and u ∈ B ∩ C 0 (B 1 ).Given some approximation set Φ k (k ∈ N fixed) such that for any l, ϕ k,l ∈ B ∩ C 0 (B 1 ); a sampling set of size S ∈ N as described in (3.3) and some regularization parameter ϵ ∈ (0, 1], we consider the approximate solution of the linear system Assume moreover that κ 2 is not an eigenvalue of the Dirichlet Laplacian in B 1 .Then, there exists a constant C err independent of u and Φ k , such that (3.9) Proposition 3.2 shows that having stable approximation sets in the sense of Definition 3.1 is a sufficient condition for the accurate reconstruction of a Helmholtz solution from its samples on the boundary of the disk, provided enough sampling points S and a sufficiently small regularization parameter ϵ are used.This is summed up in the following result, see Appendix B for its proof.
Corollary 3.3.Let δ > 0. We assume to have a sequence of approximation sets {Φ k } k∈N that is stable in the sense of Definition 3.1.Assume also that κ 2 is not a Dirichlet eigenvalue in ) . Moreover, we can take the regularization parameter ϵ as large as The point of Corollary 3.3 is not only that the solution of the regularized SVD problem provides an accurate approximation of u, but also that it is numerically computable.This is in contrast with the classical theory for the approximation by PPWs, e.g.[30], which provides rigorous best-approximation error bounds that often can not be attained numerically, precisely because accurate approximations require large coefficients and cancellation, so exact-arithmetic results cannot be reflected by floating-point computations.
The assumption on the eigenvalues in Corollary 3.3 can be lifted if in (3.10) the L 2 (B 1 ) norm is replaced by L 2 (∂B 1 ).Moreover, in this case, the constant C err at the right-hand side of (3.11) can be dropped.Finally, the largest singular value σ max of the matrix A appears in the above results: in our numerical experiments σ max is moderate, see Figure 10.
In the following, it will be convenient to measure the approximation error by the relative residual where ξ S,ϵ is the solution (3.7) of the regularized linear system.Arguing as in the proof of Proposition 3.2, (see Appendix B) for sufficiently large S, the quantity E satisfies (for a constant

Instability of propagative plane wave sets
The purpose of this section is to present the pitfalls encountered when using propagative plane waves (PPW) to approximate Helmholtz solutions in the unit disk.In particular, we show that PPWs with equispaced angles in general fail to yield stable approximations.This implies that problems can not be solved numerically to arbitrary accuracy or, in some cases, to any accuracy at all.The main effort is to show that PPW approximations lead to large expansion coefficients and that this problem can not be avoided using PPWs alone.

Propagative plane waves and Jacobi-Anger identity
We introduce the notion of a propagative plane wave.The adjective propagative is not customary in the literature, but serves to distinguish the following definition with the notion of evanescent plane waves (EPWs) that will be introduced in Definition 5.1.Propagative plane waves are a common choice in Trefftz schemes, see [22,Sec. 3.2] and the references therein.In 2D, isotropic approximations are obtained by using equispaced angles: for some M ∈ N, the approximation set is defined as In contrast to circular waves, the approximation sets based on such PPWs are in general not hierarchical.Plane waves spaces have been studied in the literature, in particular explicit hpestimates in suitable Sobolev semi-norms are available for general domains, see [30, Th. 5.2 and 5.3].These results ensure more than exponential convergence (with respect to the number of plane waves used) of the approximation of homogeneous Helmholtz solutions by a finite superposition of PPWs.Therefore, at least in principle, PPWs are well-suited for Trefftz approximations.
The Jacobi-Anger identity [31, Eq. ( 10.12.1)] provides a link between plane waves and circular waves and is ubiquitous in the analysis that follows:

Herglotz representation
We recall the so-called Herglotz functions.They are defined for any v ∈ L 2 ([0, 2π]) as see [14,Eq.(1.1)], [35, Eq. ( 6)] and [16,Def. 3.18].Such an expression is termed Herglotz representation.The function v is called Herglotz kernel or density.These functions u( are entire solutions of the Helmholtz equation and can be seen as a continuous superposition of PPWs, weighted according to v. ), which we write as a Fourier expansion for a sequence of coefficients (v p ) p∈Z ∈ ℓ 2 (Z).Plugging this expression into (4.4) and using the Jacobi-Anger expansion (4.3) together with the orthogonality of the complex exponentials {θ → e ıpθ } p∈Z , we obtain, for any thanks to the super-exponential growth of the coefficients {β p } p∈Z shown in Lemma 2.4.
While circular waves do have a Herglotz representation, their Herglotz densities are not bounded uniformly with respect to the index p.For any p ∈ Z and x = (r, θ) ∈ B 1 , using once again Jacobi-Anger expansion (4.3) together with the orthogonality of the complex exponentials, we have 2π 0 e ıpφ PW φ (x) dφ = 2π 0 e ıpφ q∈Z ı q J q (κr)e ıp(θ−φ) dφ = 2πı p J q (κr)e ıpθ .
Hence, we obtain the Herglotz representation of the circular waves, sometimes referred to as Bessel's first integral identity [33, Eq. ( 6)].The associated Herglotz density, φ → β p (2π) −1 ı −p e ıpφ , is clearly not bounded uniformly with respect to the mode number p, as a consequence of Lemma 2.4.As a result, the discretization of this exact integral representation (e.g. by the trapezoidal rule), cannot yield approximate discrete representations with bounded coefficients, as we establish next.Moreover, several solutions of the Helmholtz equation can not be represented in the form (4.4) for any v ∈ L 2 ([0, 2π]).For any sequence (û p ) p∈Z ∈ ℓ 2 (Z), the function u = p∈Z ûp b p belongs to B, because {b p } p∈Z is a Hilbert basis.If this u admits a Herglotz representation in the form (4.4) then the coefficients {v p } p∈Z of the Fourier expansion (4.5) of the density v satisfy the relation vp = ı −p β p ûp for all p ∈ Z.For v to belong to L 2 ([0, 2π]), these coefficients would need to belong to ℓ 2 (Z).This is only possible if the coefficients {û p } p∈Z decay superexponentially, to compensate for the growth of {β p } p∈Z , again by Lemma (2.4).For instance, the PPWs themselves are not Herglotz functions, because their Fourier coefficients do not decay sufficiently fast, as can be readily seen from the Jacobi-Anger identity (4.3) (in particular, for a PPW |ı −p β p ûp | = 1 for all p).In fact, the density v for a PPW would need to be a generalized function, the Dirac distribution.

Propagative plane waves do not give stable approximations
We investigate the approximation of a circular wave b p for some p ∈ Z by a generic sequence of approximation sets made of PPWs.It is shown that the two conditions in (3.2), namely accurate approximation and bounded coefficients, are mutually exclusive.Thus, stable approximations with PPWs are not possible.
Proof.Let M ∈ N and µ := {µ m } M m=1 ∈ C M .Using the Jacobi-Anger identity (4.3) we obtain at µ m e −ıqφm J q (κr)e ıqθ , so that T Φ M µ = q∈Z c q bp , where the coefficients c q := ı q / √ M 1≤m≤M µ m e −ıqφm satisfy To ensure that the approximation error p , which can be rewritten as (4.7), recalling that ∥b p ∥ B = 1.
This bound means that if one approximates circular waves b p in the form of PPW expansions T Φ M µ with a given accuracy (i.e.small η > 0), then the norms of the coefficients ∥µ∥ ℓ 2 need to increase at least like the normalization constant β p , i.e. super-exponentially fast in |p|, see Lemma 2.4.This is a clear example of accuracy and stability properties being opposite to each other.We state this important conclusion as a theorem to stress the message.Theorem 4.3.There does not exist a sequence of approximation sets made of PPWs that is a stable approximation for the space of Helmholtz solutions on the disk.
Proof.Lemma 4.2 exhibits a particular sequence, the sequence of circular waves {b p } p∈Z , for which any generic sequence of PPW approximation sets {Φ M } M ∈N does not provide stable approximations.Indeed, let p ∈ Z and suppose there exist M ∈ N and µ ∈ C M such that ∥b p − T Φ M µ∥ B ≤ η∥b p ∥ B for some η ∈ (0, 1).Then ∥µ∥ ℓ 2 ≥ (1 − η)β p ∥b p ∥ B , which implies that ∥µ∥ ℓ 2 cannot be bounded uniformly with respect to p in virtue of Lemma 2.4.The stability condition (3.2) is not satisfied and we conclude that any sequence of PPW approximation sets {Φ M } M ∈N is unstable in the sense of Definition 3.1.
More generally, this statement has implications for other Trefftz methods as well.It is not sufficient to study the best approximation error in a space spanned by Trefftz elements.If one is interested in numerical methods, one has to study approximation properties in relation to coefficient norm, and the latter depends not only on the approximation space but also on its chosen representation, i.e. the approximation set.
In the context of the Method of Fundamental Solutions (MFS), similar instability results (exponential growth of the coefficient size) are obtained if the analytic extension of the Helmholtz solution presents a singularity closer to the boundary than the MFS charge points [5,Th. 7].
Modal analysis of a propagative plane wave.Another point of view on the same issue is directly given by the Jacobi-Anger identity (4.3).This identity allows us to get quantitative insight into the modal content of PPWs.The modulus of the coefficients ı p e −ıpφ β −1 p in the expansion as a function of p can be directly deduced from Lemma 2.4 (for large |p|) and is reported in Figure 3 (left).This quantity does not depend on the propagation angle φ which parametrizes the PPW.
These coefficients decay super-exponentially fast in modulus in the evanescent regime |p| ≥ κ.Recalling Remark 2.5, the coefficients with respect to a normalization in alternative sensible norms (L 2 (B 1 ), L 2 (∂B 1 ) or L ∞ (∂B 1 ) for instance) modify the decay only by some moderate powers of |p|.This does not come as a surprise, since PPWs are entire functions.Yet, the modal content of any PPW is fixed and low-frequency.The direct implication is that they are not suited for approximating Helmholtz solutions with a high-frequency modal content (large |p|).

Evanescent plane waves
The goal of this section is to introduce evanescent plane waves (EPWs) with a complex-valued direction vector d ∈ C 2 , as opposed to propagative ones with d ∈ R 2 , and to provide intuitive reasons for their better stability properties.PPWs and EPWs are sometimes respectively called homogeneous and inhomogeneous plane waves, since only the former have constant amplitude.Combinations of PPWs and EPWs have already been used to approximate Helmholtz solutions, e.g. in the Wave Base Method [17], and Laplace eigenfunctions, e.g. in [4, §6.1.3].Since the angle is complex, the behavior of the "wave" might be unclear.Two more explicit expressions of EPWs are, for x = (r, θ) ∈ R 2 :
We see from these formulas that the wave oscillates with apparent wavenumber κ cosh ζ ≥ κ in the direction of d(φ) := (cos φ, sin φ), which was defined in (4.

Modal analysis of evanescent plane waves
The Jacobi-Anger expansion (4.(5.1) The modulus of the coefficients ı p e −ıpφ e pζ β −1 p in the modal expansion are reported in Figure 3 (right) as functions of p.On this graph, we have conveniently normalized the coefficients according to a normalization factor (depending only on ζ) which is described in the following sections, see (7.6).We see that by tuning the evanescence parameter ζ we are able to shift the modal content of the plane waves to higher-frequency regimes.As a result, we expect EPWs to be able to capture well the higher-frequency modes of Helmholtz solutions that are less regular.These may arise, for instance, in the presence of close-by singularities.The difficulty then is to properly choose suitable values for this new evanescence parameter ζ in order to build approximation spaces that are reasonable in size.This will be the main objective of the remainder of this paper.

Mapping Herglotz densities to Helmholtz solutions
In this section we introduce an integral transform between a space of functions defined on the parametric domain [0, 2π) × R and the space of Helmholtz solutions in the unit disk B.

Space of Herglotz densities
To shorten notations we denote the parametric domain as the cylinder We introduce a weighted L 2 space defined on Y .The weight function is (the square of) for some z ∈ R. In this section, the parameter z is temporarily not specified, although the following analysis shows that it cannot be chosen freely and should take the specific value z = 1/4, see (6.6).We stress that w z does not depend on the angle φ.The weighted scalar product and associated norm are then defined by: A := (u, u) A .
We now introduce a subspace of L 2 (Y ; w 2 z ) which we call space of Herglotz densities for reasons that will be clear in the following.
The wavenumber κ appears explicitly in the weight function w z .Therefore, each a p for p ∈ Z has an implicit dependence in the wavenumber κ through the normalization factor α p .Some functions a p , weighted by w 1/4 (see (6.6)), are represented in Figure 4.
For any p ∈ Z, the complex-valued function It follows that its real and imaginary parts are harmonic functions on the cylinder Y .Lemma 6.2.The space (A, ∥ • ∥ A ) is a Hilbert space and the family {a p } p∈Z is a Hilbert basis: The coefficients α p defined in (6.2) decay super-exponentially with |p| after a pre-asymptotic regime up to |p| ≈ κ.The precise asymptotic behavior is given by the following lemma.Lemma 6.3.For a constant c(κ) only depending on κ, we have and the last term is in fact equivalent to e 2z at infinity; the claimed result follows.
Using our definitions, the Jacobi-Anger expansion (5.1) takes the simple form where we introduced Formula (6.3) relates the basis {a p } p∈Z of the space A to EPWs EW y and circular waves b p on B 1 and is the key reason for introducing the space A. The behavior of |τ p | is of crucial importance in the following analysis and is given in Figure 5 for various wavenumber κ.From the asymptotics given in Lemma 2.4 and Lemma 6.3 we deduce the following result.
Lemma 6.4.We have where the constant c(κ) only depends on κ.Hence, choosing z = 1/4, we get It is clear that the uniform bounds for |τ p | are possible only for a precise pair of norms for the space of Helmholtz solutions and the space of Herglotz densities.The bounds τ ± depend implicitly on the wavenumber κ, see Figure 5.
The uniform boundedness of τ p is the key to the following analysis.In the remainder of the paper, we set z = 1/4 in (6.1) and we let w := w 1/4 .(6.6) We conclude this subsection with a lemma that will be useful in the following.Lemma 6.5.For any x ∈ B 1 , y → EW y (x) ∈ A.
Proof.Let x ∈ B 1 and define v x : y → EW y (x).The Jacobi-Anger identity (6.3) reads v x (y) = p∈Z τ p b p (x) a p (y) for all y ∈ Y .Since {a p } p∈Z is a Hilbert basis for A, if we write x = (r, θ) ∈ [0, 1) × [0, 2π), we get Using the estimates (2.3) and (A.5) from the proof of Lemma 2.4, we get from which we conclude that ∥v x ∥ A < ∞.
If x ∈ ∂B 1 , so that r = |x| = 1, then y → EW y (x) does not belong to A, as is readily seen from the proof of Lemma 6.5.

Herglotz transform
We introduce an integral operator T that allows to write every Helmholtz solution in B as a continuous linear combination of EPWs weighted by an element of A. We also describe its adjoint operator T * , the corresponding frame and Gram operators S and G, and prove some of their properties.The terminology of this section is borrowed from Frame Theory, see [12] for a reference on this field.
Synthesis operator.The first and most important definition concerns the transform that maps Herglotz densities to Helmholtz solutions as we prove next.Definition 6.6.Using the weight (6.6), we introduce the Herglotz transform T : for any v ∈ A, This operator is well-defined on A thanks to Lemma 6.5.In the setting of continuous-frame theory, see e.g.[12,Eq. (5.27)], this operator is called synthesis operator.
The Herglotz transform T is bounded and invertible between the space of Herglotz densities A and the space of Helmholtz solutions B. Theorem 6.7.The operator T is bounded and invertible from A to B: Moreover, T a p = τ p b p for all p ∈ Z.
Proof.Using the Jacobi-Anger formula (6.3), for any v ∈ A and x ∈ B 1 we get Hence, from Lemma 2.2, , and the result (6.8) follows from Lemma 6.2 and Lemma 6.4.It is readily checked that the inverse is given, for any u ∈ B, by From (6.9), the inverse operator T −1 can also be written as an integral operator: for u ∈ B, The integral representation T v in (6.7) is similar to the Herglotz representation (4.4).This is the reason why we refer to elements of A as Herglotz densities.For any p ∈ Z, the Herglotz densities τ −1 p a p of the circular waves b p are bounded in the A-norm by τ −1 − , hence uniformly with respect to the index p.This should be contrasted with the standard Herglotz representation (4.6) using only PPWs, where the associated Herglotz densities cannot be bounded uniformly with respect to the index p in L 2 ([0, 2π]).As we explained in Section 4.2, not all Helmholtz solutions admit a bounded Herglotz representation that uses only PPWs (4.4) (with density v ∈ L 2 ([0, 2π])).In contrast, using EPWs the generalized Herglotz representation (6.7) can represent any Helmholtz solution.Indeed, since T is an isomorphism between A and B, for any u ∈ B, there exists a unique v ∈ A such that u = T v.The price to pay for this result is the need for a two-dimensional parameter domain, the cylinder Y , in place of a one-dimensional one, the interval [0, 2π), and thus of a double integral; the added dimension corresponds to the evanescence parameter ζ.Theorem 6.7 is a stability result stated at the continuous level.Next, we aim to obtain a discrete version of this integral representation.
Analysis operator.In the continuous-frame setting, see [12, Eq. (5.28)], the adjoint operator T * of T is called analysis operator.
Lemma 6.8.The adjoint T * of T is given for any u ∈ B by (T * u)(y) := (u, EW y ) B , ∀y ∈ Y .The operator T * is bounded and invertible on B: Proof.We have, for any v ∈ A and u ∈ B In addition, using the Jacobi-Anger formula (6.Frame and Gram operators.We introduce two other important operators in Frame Theory.
Corollary 6.9.The frame operator S := T T * and the Gram operator G := T * T are bounded, invertible, self-adjoint and positive operators: Proof.This result stems directly from Theorem 6.7 and Lemma 6.8.
The frame operator admits the more explicit formula: for any u ∈ B, A continuous frame result.We are now ready to prove that EPWs form a continuous frame for the space of Helmholtz solutions in the unit disk.We recall [12, Def.5.6.1]:given a complex Hilbert space H and a measure space M with positive measure µ, a family Theorem 6.10.The family {EW y } y∈Y is a continuous frame for B. Besides, the optimal frame bounds are A = τ 2 − and B = τ 2 + .Proof.We need to verify the definition of a continuous frame, see [12,Def. 5.6.1].For any u ∈ B, the measurability of y → (u, EW y ) B = (T * u)(y), stems from T * u ∈ A according to Lemma 6.8 and A ⊂ L 2 (Y ; w 2 ).The frame condition, namely for some constants A and B is a consequence of the boundedness and positivity of the frame operator S which was established in Corollary 6.9.Indeed, for any u ∈ B, we have which also establishes the optimality of the claimed frame bounds.

The reproducing kernel property
The continuous frame result implies additional structure on the Herglotz density space A, which then allows to characterize the preimages of the EPWs under the integral transform T .For a general reference on Reproducing Kernel Hilbert Spaces (RKHS), we refer to [32].
Lemma 6.11.The range of the analysis operator T * , i.e. the space A defined in (6.2), has the reproducing kernel property.The reproducing kernel is given by with pointwise convergence of the series and where K y ∈ A is the (unique) Riesz representation of the evaluation functional at y ∈ Y , namely Proof.Take any v ∈ A and let u ∈ B such that v = T * u, which exists thanks to Lemma 6.8.From Corollary 6.9, we have Then we obtain the reproducing identity, for any where we introduced (the Riesz representation of) the evaluation functional at the point y defined as K y (z) := EW y , S −1 EW z B , ∀z ∈ Y.It is a direct consequence of Corollary 6.9 that the kernel admits the series representation (6.10).Alternatively, we refer to [32,Th. 2.4] for a direct proof of this result (valid in the general setting), since {a p } p∈Z is an orthonormal basis for A.
Examples of (normalized) evaluation functionals are given in Figure 6.
The interest in introducing the reproducing kernel property stems from the following result, which is a direct consequence of Lemma 6.11, Theorem 6.7, and the Jacobi-Anger identity (6.3).Corollary 6.12.The EPWs are the images under T of the Riesz representation of the evaluation functionals, namely As a consequence, the construction of an approximation of a Helmholtz solution u ∈ B as an expansion of EPWs is, up to the isomorphism T , equivalent to the approximation of its Herglotz density v := T −1 u ∈ A as an expansion of evaluation functionals, i.e.
for some set of coefficients µ = {µ m } M m=1 .This remark justifies the use of the sampling techniques described in the next section to discretize the integral representation in (6.7).Section 8 provides numerical evidence that such approximations can be built, for a suitable normalization of the sets {K ym } m and {EW ym } m .

A concrete evanescent plane wave approximation set
We describe a method for the numerical approximation of a general Helmholtz solution in the unit disk by EPWs.We exploit the equivalence of this approximation problem with the approximation problem of the corresponding Herglotz density, see (6.11).The main idea is to adapt the sampling procedure of [21,13,29] (sometimes called coherence-optimal sampling) to our case, in order to generate a distribution of sampling nodes in the cylinder Y that will be used to reconstruct the Herglotz density.While the numerical recipe that we describe below is found to be numerically very effective, see Section 8, our theoretical analysis still lacks a formal proof of the accuracy and stability of the approximation of Helmholtz solutions using EPWs.
Let u ∈ B be the Helmholtz solution, target of the approximation problem, and let v := T −1 u ∈ A be its associated Herglotz density.Let also some tolerance η > 0 be given.

Truncation of the modal expansion
Since u (resp.v) a priori lives in an infinite dimensional space B (resp.A), the idea behind the construction of finite dimensional approximation sets is to exploit the natural hierarchy of finite dimensional subspaces constructed by truncation of the Hilbert basis {b p } p∈Z .
Truncation in the Helmholtz solution space.For any P ∈ N, we define and Here Π P is the orthogonal projection from B onto the finite dimensional subspace B P .A natural approach to compute an approximation of u ∈ B is to approximate its projection onto B P , namely for some P large enough.It is immediate that the sequence of projections {u P } P ∈N converges to u in B. In particular, we can define for any η > 0 Unfortunately, it is not possible to compute such a P * in most practical configurations.It may be possible though to give estimates on P * , based on some regularity assumption on u and the decay of its coefficients in its modal expansion.For instance, it might be physically realistic to assume that all coefficients of the propagative modes |p| ≤ κ are O(1) and the coefficients associated to the subsequent evanescent modes |p| ≥ κ decay in modulus with a given algebraic or exponential rate.
Truncation in the Herglotz density space.Similarly, for any P ∈ N, we define and v P := T −1 u P ∈ A P , ∀P ∈ N.
Theorem 6.7 implies that the sequence {v P } P ∈N converges to v in A. In particular, for any P ≥ P * , where P * was defined in (7.1), we have (up to some normalization factor) is expected to provide a good approximation of v P .The approximation set for u P will then be given by the EPWs {EW ym } m (up to some normalization factor).We denote the dimension of both spaces A P and B P by The probability density function ρ P is defined (up to normalization) as the reciprocal of the N P -term Christoffel function µ P in the spirit of [13, Eq. (2.6)]: Observe that ρ P and µ P are well-defined since 0 < µ P ≤ µ 0 < ∞ from the fact that a 0 is just a non-vanishing constant.The density function ρ P is a univariate function on Y since it is independent of the angle φ.We point out that 1/µ P corresponds to the truncated series expansion of the diagonal of the reproducing kernel K, which amounts to taking z = y and truncating at P the series in (6.10).The numerical recipe consists, for each P ∈ N, in generating a sequence of sampling node sets in the parametric domain Y where using one's preferred sampling strategy such that |Y P,M | = M for all M ∈ N and the sequence Y P converges (in a suitable sense) to the density ρ P defined in (7. 3) as M tends to infinity.
The sampling method could be a deterministic, a random or even a quasi-random strategy, see Section 8.The sets are not assumed to be nested.This choice of EPW parameters is a major difference from the heuristic choice described in [24, Eq. ( 5)] where the parameters are chosen in order to approximate solutions defined in a rectangle containing the physical domain of interest (B 1 in our case).

Evanescent plane wave approximation sets
From the sampling node sets (7.4) we can construct two approximations sets: one set of sampling functionals in A and one set of EPWs in B.
Approximation sets in the Herglotz density space.Associated to the sampling node sets (7.4), we introduce a sequence of finite sets in A The normalization of K ym in (7.5) is crucial for the stable approximation property (3.2).In the approximation sets, each sampling functional K ym has been normalized by the real constant µ P (y m )/M which is (numerically) close to ∥K ym ∥ −1 A / √ M .More precisely, we have Approximation sets in the Helmholtz solution space.Associated to the sampling set sequences (7.4) and approximation set sequences (7.5) in A, we define the sequence of approximation sets of (normalized) EPWs in B as follows Following Corollary 6.12, the sequence of sets (7.6) is the image of the sequence of sets (7.5) by the Herglotz transform operator T .
Discussion on the parameters.Our numerical recipe for building the approximation sets Φ P,M is based on only two parameters, P and M , whose tuning is intuitive: 1.The first one is the Fourier truncation parameter P .Increasing P will improve the accuracy of the approximation of u (resp.v = T −1 u) by u P = Π P u (resp.v P = T −1 u P ).The appropriate value for P ≥ P * will solely depend on the decay of the coefficients in the modal expansion, which is intimately linked to the regularity of the Helmholtz solution.
2. The second one is the dimension M of the EPW approximation space, which is also the number of sampling points in the parameter cylinder Y .For a fixed P , increasing M should allow to control the accuracy of the approximation of u P (resp.v P = T −1 u P ) by T Φ P,M ξ (resp.T Ψ P,M ξ) for some bounded coefficients ξ ∈ C M .The numerical results presented below corroborate this conjecture and show experimentally that M should scale linearly with P , with a moderate proportionality constant (see Section 8.5).
For a fixed DOF budget M , the numerical experiments in Section 8.5 suggests that using a Fourier truncation parameter P = max (⌈κ⌉, ⌊M/4⌋) gives accurate and reliable approximations.
Once the approximation sets Φ P,M are chosen, our concrete implementation (see Section 3.3) to compute a particular set of coefficients ξ S,ϵ includes two additional parameters, S and ϵ: 1.The first parameter S is the number of sampling points on the boundary of the physical domain B 1 .According to (B.1) and following [1,2], sufficient oversampling should be used.
In practice, we chose for simplicity an oversampling ratio of 2, namely S = 2M .This amount of oversampling may not be necessary and further numerical experiments could investigate a reduction of the oversampling ratio S/M to reduce the computational cost.
2. The second parameter ϵ is the regularization parameter, i.e. the truncation threshold of the singular values.We set this parameter to ϵ = 10 −14 in the numerical experiments presented below.If one is interested in less accurate approximations than ours, this parameter could be set to larger values.
We stress that the construction of the approximation sets Φ P,M , together with their accuracy and stability, are not influenced by the choice of the reconstruction strategy made in Section 3.2.
Although we focus on the simple method of boundary sampling together with regularized SVD, alternative reconstruction strategies (such as sampling in the bulk of the domain or taking inner product with elements of other types of test spaces, for instance) and other regularization techniques (such as Tikhonov regularization) can also be successfully used in practice.Irrespective of the strategy, sufficient oversampling and regularization need to be used.
Relation with the literature.As we have already alluded to, our construction is based on similar ideas that pre-exist in the literature but in a different context.Indeed, sampling node sets similar to the ones we propose here can be found in [21,13,29].The context of these works is the reconstruction of elements of finite-dimensional subspaces (with explicit orthonormal basis) in weighted L 2 spaces from sampling [13] and it was subsequently used to construct random cubature rules [29].The underlying idea is that the information gathered from sampling at these nodes is enough to allow accurate reconstruction as an expansion in the (truncated) orthonormal basis.Translated into our setting, the results available in the literature say that to reconstruct an element v P = Π P v of the finite dimensional subspace A P , it is enough to sample at the nodes Ψ P,M for some sufficiently large M .In contrast, the numerical recipe described above seeks to construct an approximation of the element v P = Π P v ∈ A P as an expansion in the set of evaluation functionals Ψ P,M for some sufficiently large M .In other words, the approximation we are looking for belongs to the span of the evaluation functionals, span Ψ P,M , which has trivial intersection with A P .By Corollary 6.12, applying the Herglotz transform T to this approximation in span Ψ P,M yields an element in span Φ P,M (i.e. a finite superposition of EPWs) that approximates u P = T v P ∈ B P .
Unfortunately, besides the links with these works, we are not yet able to prove a rigorous theoretical analysis to support our numerical recipe.Yet, extensive numerical experiments in Section 8 illustrate the excellent approximation and stability properties of the sets Φ P,M .

A conjectural stable approximation result
We formalize below our speculations, which are hinted by the numerical experiments given in the next section.First, we state our main conjecture.
Conjecture 7.1.The sequence of approximation sets Ψ P defined in (7.5) is a stable approximation for A P , in the following sense: there exist s ≥ 0 and C > 0 such that, for all P ∈ N, there exists M * = M (P, η) such that In the following we assume for simplicity that all M ≥ M * satisfy the two inequalities appearing in (7.7) (otherwise the proofs can be easily adapted).This holds true if the sets are hierarchical, for instance, but this is not necessary.
Provided the above conjecture holds, the stability of the approximation sets of EPWs constructed above would follow as we prove next.Proposition 7.2.Let δ > 0. If Conjecture 7.1 holds then the sequence of approximation sets (7.6) provides a stable approximation for B. In particular, if κ 2 is not a Dirichlet eigenvalue on B 1 , where ξ S,ϵ ∈ C |Φ P,M | is computed with the regularization procedure in (3.7).The SVD regularization parameter ϵ can be chosen as (3.11).
Proof.We need to prove the stability of the sequence of approximation sets, namely that for any η > 0, there exists s ≥ 0 and C > 0 such that Assuming that Conjecture 7.1 holds, there exist s and C (both independent of P ) such that, for any M ≥ M * (P * , η), there exists a set of coefficients µ ∈ C M such that ∥v P − T Ψ P,M µ∥ A ≤ η∥v P ∥ A , and ∥µ∥ ℓ 2 ≤ CM s ∥v P ∥ A .
The properties of the isomorphism T given in Theorem 6.7 imply that For any P ≥ P * (u, η) and M ≥ M * (P * , η), the total approximation error for the Herglotz density v can be estimated as

Numerical results
We provide numerical evidence that the procedure described above allows to compute controllably accurate approximations of Helmholtz solutions in the unit disk and in other domains1 .

Probability densities and samples
Probability density and cumulative distributions functions.We represent the probability density function ρ P (see (7.3)) as a function of the evanescence parameter ζ on the left in Figure 7.Here P denotes the truncation parameter, meaning that the sampling is performed to approximate elements of A P , which has dimension N P .The associated cumulative distribution function with respect to the evanescence parameter ζ is defined as It is represented in the right of Figure 7. Recall that while ρ P is a bi-variate function on the cylinder Y , it is constant with respect to the angle φ.As a result, the cumulative distribution with respect to this variable φ is a linear function.This is why we represent these two functions ρ P and Υ P only with respect to the evanescence parameter ζ.
We observe that the probability density ρ P is a symmetric even function and exhibits a main mode at ζ = 0 which corresponds to purely PPWs.Moreover, the ϵ-support of this density is rather tight and the probability eventually tends to zero exponentially as |ζ| gets large enough.When P ≤ κ the density is a unimodal distribution whereas for P ≫ κ (e.g.P = 4κ) the density is a multimodal distribution.Indeed, in the latter case, there are two symmetric modes for relatively large evanescence parameter, in addition to the main mode at ζ = 0.The cumulative distribution function Υ P is close to a step function in the case where A P contains only elements associated to the propagative regime P ≤ κ.In contrast, for P > κ the distribution is non-trivial for moderate values of the evanescence parameter ζ.This means that for P ≤ κ one can safely choose only PPWs, while for P > κ EPWs are needed and their choice is non-trivial.The fact that the density function is constant with respect to φ considerably simplifies the generation of the samples.The inversion Υ −1 P can be performed using elementary root-finding techniques, our implementation resorts to the bisection method.
In our numerical experiments we tested three types of sampling methods, which differ by how we generate the first sampling distribution {z m } m in the unit square: 1. deterministic sampling: the initial samples in the unit square are a Cartesian product of two sets of equispaced points with the same number of points in both directions (all numerical results presented are obtained by using as approximation set dimension the smallest square integer larger than or equal to M ); 2. Sobol sampling: the initial samples in the unit square corresponds to Sobol sequences which are quasi-random low-discrepancy sequences2 ; 3. random sampling: the initial samples in the unit square are drawn randomly according to the product of two uniform distributions Some examples of sampling sets corresponding to the probability density function ρ P for κ = 16 are reported in Figure 8.For these examples the number of sampling nodes is set to M = νN P with ν = 4, for the three types of sampling considered.As expected, the sampling points cluster near the line ζ = 0 for smaller P .This is the (propagative) regime for which PPWs alone provide a good approximation.When P > κ the evanescence parameter ζ spreads in a wider domain, with some clustering at the secondary modes of the distribution, in agreement with Figure 7.

Propagative plane waves are unstable
Before presenting EPW approximations, we report some numerical experiments dedicated to verifying numerically the instability result of Lemma 4.2 when using PPWs.These will also serve as a reference point to compare with the results obtained using our EPW recipe.
Let us consider the approximation problem of Section 4.3, namely the approximation of the circular wave b p for some p ∈ Z by an approximation set Φ M of M ∈ N PPWs defined in (4.2).The sampling matrix A was defined in (3.4), using M PPWs with equispaced angles and S := max(2M, 2|p|) sampling points (we impose S ≥ 2|p| to avoid spurious results due to aliasing).The entries of the matrix A are immediately computed as A s,m = e ıκ cos(2π( s S − m M )) for s = 1, . . ., S, m = 1, . . .M .The right-hand side b is defined as in (3.4) for b p in place of u; we recall that we use Dirichlet data in all our numerical experiments.
The matrix A is notoriously ill-conditioned (see Figure 10a): its condition number grows exponentially with respect to the number of plane waves M in the approximation set Φ M .This is well-known, see for instance the numerical experiments in [33,Sec. 2.3] for the circular geometry and S = M .This is not a feature of the sampling method: we refer to similar experiments in [22,Sec. 4.3] for the mass matrix of a Galerkin formulation in a Cartesian geometry, again for S = M .The least-squares formulation suffers from an even worse condition number: proportional to the 0 square of the condition number of the sampling method, see e.g.[33,Eq. (30)].We apply the regularization procedure described in Section 3.3 with threshold parameter ϵ = 10 −14 .
The numerical results are reported in Figure 9a.On the left panel we report the relative residual E defined in (3.12) as a measure of the accuracy of the approximation.On the right panel we report the size of the coefficients ∥ξ S,ϵ ∥ ℓ 2 as a measure of the stability of the approximation.Relative residuals and coefficient norms were already used in [24] to assess the stability of the approximations.
We observe three regimes.First, for the propagative modes, i.e. the circular waves with mode number |p| ≤ κ, the approximation is accurate (E < 10 −13 ) and the size of the coefficients is moderate (∥ξ S,ϵ ∥ < 10).Second, for mode numbers |p| roughly larger than the wavenumber κ, the norms of the coefficients of the computed approximations blow up exponentially.The accuracy is spoiled proportionally.Third, for evanescent modes with |p| larger than about 2κ or 3κ, the size of the coefficients completely destroys the stability of the approximation, and we cannot approximate the target b p with any decent accuracy.Of course, for a relative error at O(1), the coefficient norm reported is not meaningful, and taking ξ S,ϵ identically zero would provide a similar error.
Increasing the number of plane waves M has no effect on the accuracy beyond a certain point.Indeed, Figure 10a shows that the ϵ-rank (i.e. the number of singular values larger than ϵ) of the matrix A does not increase when M is raised.Although increasing M does not bring any higher accuracy, it does not increase any further the numerical instability.For a fixed M , the same matrix A is used to approximate all the b p 's for any mode number p (i.e. to compute all markers of the same color in Figure 9a).Even when the matrix A is extremely ill-conditioned (say M = 32κ in the numerical experiments presented here), we get at the same time almost machine-precision accuracy for all propagative modes |p| ≤ κ while having O(1) error for evanescent modes with larger mode number |p| ≥ 3κ.It is the simple regularization procedure described in Section 3.3 that allows us to obtain such results.No other technique can overcome the inherent instability of PPWs.In particular, even with regularization, accuracy in the approximation of the evanescent modes remains out of reach for a given floating-point precision.
Analoguous numerical results are also observed in the context of the MFS, see [5, Fig. 3]. Figure 11: Accuracy E, as defined in (3.12), (left) and stability ∥ξ S,ϵ ∥ ℓ 2 /∥u∥ B (right) of the approximation by M EPWs (constructed using Sobol sampling) of a solution u in the form (8.2) that belong to the space B P of dimension N P = 2P + 1.The horizontal axis represents the ratio M/N P .Wavenumber κ = 16.The number M of EPWs necessary to approximate elements of the space B P seems to scale linearly with the space dimension N P .approximation spaces made of PPWs are tuned for propagative modes, which span a space of dimension 2κ + 1.In contrast, the approximation spaces made of EPWs target a larger number of modes, including some evanescent modes, which span a space of dimension N P = 2P + 1 with P = 4κ in this numerical experiment.For a general target solution containing evanescent modes, one does not expect any advantage in using PPWs only.

Approximation of random-expansion solutions
We test the procedure described so far by reconstructing a solution of the form in which ûp are normally-distributed random numbers with mean 0 and standard deviation 1.
The coefficients of any element of B decay in modulus as o(|p| −1/2 ) for |p| → ∞; this is therefore a rather difficult scenario for an approximation problem.
We then apply the procedure described above for the three types of sampling strategies considered.The sampling points are constructed knowing that T −1 u is an element of A P .In other words, the optimal modal truncation parameter P * = P (where P appears in (8.2)) is assumed to be known in this numerical experiment.The main purpose is to investigate the validity of Conjecture 7.1.We study here the convergence of the error with respect to the dimension of the approximation space M .The number of sampling points on the boundary of the disk is set to S = 2M .The numerical results are given in Figure 11 for the Sobol sampling strategy only.On the left panel we report the relative residual E, defined in (3.12), as a measure of the accuracy of the approximation.On the right panel we report the size of the coefficients, namely ∥ξ S,ϵ ∥ ℓ 2 /∥u∥ B , as a measure of the stability of the approximation.
The main observation is that the error quickly decays with respect to the ratio M/N P = M/(2P + 1), which represents the ratio of the dimension of the approximation set M over the dimension of the space B P the solution (8.2) lives in.When P is large enough (say P ≥ 2κ which remains moderate), the decay is relatively independent of P .The second observation is that the norm of the coefficients ∥ξ S,ϵ ∥ ℓ 2 /∥u∥ B in the expansions is a decreasing function of the size M of the approximation space.We see once more that one gets accurate and stable approximations.The values of ∥ξ S,ϵ ∥ ℓ 2 /∥u∥ B reported for small values of M/N P , and in particular the increase at the start, are not significant since they correspond to inaccurate approximations.We report in Figure 12 the plots of a solution (8.2) for a larger frequency κ = 64 and truncation parameter P = 3κ = 192.The approximation error when using M = 3(2P + 1) = 1155 PPWs or EPWs is also given, with points in Y sampled as a Sobol sequence.In the first case the absolute error in the disk is much larger, more than 12 orders of magnitude larger if measured in L ∞ (B 1 ) norm, and concentrated near the boundary.The number of degrees of freedom per wavelength λ = 2π/κ used in each direction can be estimated by λ M/π ≈ 1.9.Note that π here represents the area of the unit disk.For low-order methods, a common rule of thumb is to use around 6 ∼ 10 degrees of freedom per wavelength to have 1 or 2 digits of accuracy.We obtain 12 digits of accuracy for only a fraction of this number.For M = 2(2P + 1) = 770, the maximum absolute error reached is measured to 1.3 • 10 −10 (not plotted).
Overall, the numerical results are perfectly consistent with Conjecture 7.1.

Numerical evidence of quasi-optimality
An important question regarding the efficiency of the proposed method concerns how the size of the approximation set M should vary with respect to the truncation parameter P .Fixing P amounts to looking at the finite dimensional subspace B P which contains the first N P = 2P + 1 modes.Since N P is the dimension of B P there is no hope to have approximation spaces with dimension M < N P that are able to approximate all elements of this space.An optimal approximation set would therefore achieve this with M = N P elements at best.We show  The number of EPWs necessary to approximate elements of the space B P to relative accuracy σ seems to scale linearly with the space dimension N P .
numerical evidence that we achieve quasi-optimality, in the sense that the approximation spaces Φ P,M defined in (7.6)only need M = O(N P ) with a moderate proportionality constant to approximate the N P circular modes with reasonable accuracy.We investigate numerically the linearity of the relation P → M * (P, η), where M * (P, η) was defined in Conjecture 7.1 (for a fixed η), namely the validity of a law of the form M * (P, η) ≈ νN P = ν(2P + 1) for some ν = ν(η) > 0. To that end, for some σ > 0, we vary P and compute M * = M * (P, σ) := min M ∈ N | E(b p , Φ P,M , S, ϵ) ≤ σ, ∀|p| ≤ P , where E was defined in (3.12).The quantity M * is expected to be a good estimate of M * (P, η).The number of sampling points on the boundary of the disk is set to S = 2M .
The numerical results are given in Figure 13 for the accuracy level σ = 10 −12 .We represent here the variation of the ratio M * (P, σ)/N P with respect to the truncation parameter P .If the optimal law for M * (P, σ) was linear with respect to P , we would expect constant values.Regardless of the type of sampling, we observe decreasing curves that converge to some asymptotic value for ν that falls within the rather moderate range [3,6].This means that the first N P circular modes (propagative and evanescent) can be stably approximated with uniform relative error ≤ 10 −12 using roughly 3N P to 6N P EPWs.Moreover, this asymptotic behavior seems to be robust with respect to the wavenumber κ.These more systematic results confirm what was already observed in Section 8.4.The behavior of the optimal asymptotic M * with respect to N P seems indeed to be linear or even sub-linear.

Triangular domain
We conclude this section with some numerical results on a triangular geometry.Our purpose is to show that the approximation sets that we constructed also exhibit good approximation properties on other shapes, despite being built following the analysis for the disk.
We study the convergence of the approximation by plane waves for increasing size of the approximation set M .The approximation is constructed as indicated in Section 3.2-3.3from Dirichlet data at equispaced points on the boundary of the triangle and by solving the oversampled linear systems using a regularized SVD.The plane waves used in the approximation sets are either propagative, with uniformly spaced angles as described in (4.2), or evanescent, as Figure 16: Point-wise error in the triangle between the target of the approximation problem (see Figure 14) and the approximation using propagative (top) and evanescent (bottom) plane waves.The singularity in the solution is either close to the edge (left) or close to the vertex (right).Wavenumber κ = 16 and M = 300.

Conclusions
Ill-conditioning is inherent in plane-wave based Trefftz schemes but can be overcome if there exist accurate approximations that are moreover stable, in the sense of having expansions with bounded coefficients.To approximate Helmholtz solutions, PPWs are known to provide accurate approximations.However, the associated expansions are necessarily unstable: the norm of the coefficients blow up for solutions with high-frequency Fourier modes.In contrast, EPWs, which contain high-frequency content, give accurate as well as stable results.To construct stable sets of EPWs, we show numerically that an effective strategy is to sample the parametric domain according to a fully explicit probability measure.
This paper is only the first step towards stable and accurate approximation schemes based on EPWs.A theoretical problem that we have left open is the analysis of the approximation properties of the sets of EPWs constructed using our numerical recipe.Next steps include the extensions to more general geometries, three-dimensional problems (see [19]), time-harmonic Maxwell and elastic wave equations, the application to Trefftz schemes and to sound-field reconstruction algorithms.Preliminary experiments show that the proposed numerical recipe performs well for convex polygons and in Trefftz-Discontinuous Galerkin schemes with several cells, and provides a considerable improvement over standard PPW schemes.

Figure 1 :
Figure 1: Real part of the circular waves bp for three different modes (wavenumber κ = 16).

1 )
and is parallel to ℜ[d(y)].In addition, the wave decays exponentially with rate κ sinh ζ in the orthogonal direction d(φ) ⊥ , which is parallel to ℑ[d(y)].This justifies naming the new parameter ζ ∈ R, which controls the imaginary part of the angle, the evanescence parameter.

2 Figure 3 :
Figure 3: Modal analysis computed using Jacobi-Anger identity (5.1) of PW φ (left) and EW φ,ζ after normalization (right).In both cases, the absolute values of the coefficients of the expansion of the plane wave in the basis {b p } p∈Z is represented against the mode number p. Wavenumber κ = 16.Modifying φ has no influence, modifying ζ shifts the modal content in the Fourier space.

Figure 6 :
Figure 6: Representation of normalized evaluation functionals |wK y |/∥K y ∥ A in the cylinder Y for wavenumber κ = 16.

Figure 12 :
Figure 12: Solution u, target of the approximation, defined in (8.2) with P = 3κ = 192 (top) and associated absolute errors when approximated by M = 3(2P + 1) = 1155 plane waves, either propagative ones Φ M from (4.2) (bottom left) or evanescent ones Φ P,M from (7.6), whose parameters are constructed using a Sobol type sampling (bottom right).The colormaps associated to absolute errors are logarithmic for better visualization.Wavenumber κ = 64.Note the different color scales, which shows a factor-10 12 improvement in using EPWs instead of PPWs.

Figure 13 :
Figure13: Ratio M * (P, σ)/N P with respect to the truncation parameter P for various types of sampling method and σ = 10 −12 .The number of EPWs necessary to approximate elements of the space B P to relative accuracy σ seems to scale linearly with the space dimension N P .
in B 1 , and γu = g, on ∂B 1 , for some suitable boundary data g.Examples of such a trace operator γ are: the Dirichlet trace operator, extension to H 1 (B 1 ) of u → u| ∂B1 , when κ 2 is not an eigenvalue of the Dirichlet Laplacian; the Neumann trace operator, extension to H 1 (B 1 ) of u → ∂ n u, when κ 2 is not an eigenvalue of the Neumann Laplacian; the Robin trace operator, extension to H 1 (B 1 ) of u → ∂ n u − ıκu| ∂B1 (without assumptions on the wavenumber κ).
Sec. 2.1], [13, Sec.2.2] and [29, Sec.2].Despite having an unbounded parametric domain Y , the finite integrability of the weight function w 2 allows to sample Y on a bounded region only.The associated set of sampling functionals {K ym } m .2)7.2 Parameter sampling in the cylinder YOur objective is to approximate the truncated Fourier series u P = Π P u ∈ B P for some P ∈ N, instead of u.Up to the Herglotz transform, this problem is equivalent to the approximation of v P = T −1 u P ∈ A P .In this subsection let us fix a P ∈ N, not necessarily equal to P * .We propose to build approximations of elements of A P (resp.B P ) by constructing a finite set of sampling nodes {y m } m in the cylinder Y according to the distribution advocated in[21,