{'fraction_non_alphanumeric': 0.047147210010932014, 'fraction_numerical': 0.015770009071244154, 'mean_word_length': 3.8124020595477948, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interpolation properties of generalized plane waves

This paper aims at developing new shape functions adapted to smooth vanishing coefficients for scalar wave equation. It proposes the numerical analysis of their interpolation properties. The interpolation is local but high order convergence is shown with respect to the size of the domain considered. The new basis functions are then implemented in a numerical method to solve a scalar wave equation problem with a mixed boundary condition. The order of convergence of the method varies linearly with the one of the interpolation.


Introduction
This paper focuses designing Generalized Plane Waves (GPW) to approximate smooth solutions u ∈ C ∞ (Ω) of the model problem where β is in C ∞ (Ω). This time-harmonic equation, generally called the scalar wave equation, models for instance the acoustic pressure describing the behavior of sound in matter or a polarized electromagnetic wave propagating in an isotropic medium. If β = −ω 2 , ω ∈ R, the equation is classically named the Helmholtz equation and is still the subject of recent research, see for instance [24]. If β < 0 is non constant, this is a simple model of wave propagation in an inhomogeneous medium. If β > 0, it models an absorbing medium and the partial differential equation is coercive. The applications considered here include both propagative and absorbing medium, as well as smooth transitions in between them, i.e. respectively β > 0, β < 0 and β = 0. Several types of numerical methods are used for the simulation of wave propagation. Classical finite element methods applied to such problems are known to be polluted by dispersion, see [1]. An alternative is to consider approximation methods based on shape functions that are local solutions of the homogeneous equation: this justifies the development of Trefftz-based methods, first introduced in [28], that rely on solutions of the homogeneous governing domain equation: information about the problem is embedded in the finite dimension basis functions set. The present work originated from the idea to apply such a method to a problem modeled by (1) in which the coefficient is likely to vanish: shape functions adapted to this problem are here designed and studied. See the previous work [20] for the physical motivation of the problem. Refer to [26] and references therein for more recent developments of these Trefftz-based methods, and to [11,16] for applications linked to one specific method, the so-called Ultra Weak Variational Formulation (UWVF). The method coupling the latter to the adapted shape functions is the topic of [20], and the present work includes numerical results for the h convergence of the coupled method.
The novelty in the present paper lies in the smooth feature of the coefficient β of the governing domain equation (1) and on the explicit procedure proposed to design the corresponding shape functions. This work can be compared to recent works that focus on non polynomial methods for smooth varying coefficients, see for instance [3,27]. The design of shape functions adapted to smooth and possibly vanishing coefficients that is the core of this work starts from mimicking the equation that shows that classical plane waves functions e Iω − → k · − → x are exact solutions of (1) when β = −ω 2 is constant and negative. The case of a piecewise constant coefficient is addressed for example in [5,10], and the more general case of a smooth coefficient is generally approximated by a piecewise constant coefficient on each cell of the mesh. A very simple extension of classical plane waves for a positive or negative constant coefficient would be to consider at a point G = (x G , y G ) the shape function ϕ(x, y) = exp sgn(β(G)) |β(G)| ((x − x G ) cos θ + (y − y G ) sin θ) , where the parameter θ represents the direction of the plane wave. Indeed, the case β(G) < 0 corresponds to the classical plane wave whereas the case β(G) > 0 corresponds to a complex wavenumber. This choice will provide a tool to extend the interpolation results cited previously. Remark that if the coefficient β(G) < 0 goes to zero then the corresponding classical plane waves functions, generates by equi-spaced directions θ, tend not to be independent anymore. Consider here the case of a general smooth coefficient β. Typically, in the case β(x, y) = x the Airy functions Ai and Bi are solutions ; however, in the general case there is no exact analytic solution known. Indeed, as was explained in Section 2.1 of [20], no exponential of a polynomial can solve a generic scalar wave equation. So the idea is to generalize the classical plane wave function as an approximated solution of the initial equation, in the following sense : design ϕ = e P , with a complex polynomial P (x, y) = dP i=0 dP −i j=0 λ i,j (x − x G ) i (y − y G ) j such that (−∆ + β)e P = − (∂ 2 x P + (∂ x P ) 2 + ∂ 2 y P + (∂ y P ) 2 ) + β e P is locally small. More precisely, the generalized plane wave described in this paper will be designed to satisfy locally up to a given order q ∈ N, satisfying q ≥ 1. In this process, the degree and coefficients of P will be chosen to satisfy (3), and this approximation identity is equivalent to canceling the q(q + 1)/2 coefficients of lower degree terms in the Taylor expansion of its left hand side. It provides a system which unknowns are the coefficients of P and which size does depend on q. This resulting system can be either underdetermined or overdetermined, depending on the degree of P , denoted dP , with respect to the value of q. A specific procedure will be described in order to obtain a square invertible system. Its main feature is based on the idea of generalizing the classical plane wave function as displayed for β(G) < 0 in (2), by setting P (x, y) = β(G)((x − x G ) cos θ + (y − y G ) sin θ)+ higher order terms. As a result the reasoning leading to interpolation properties of such new functions is built on the simpler case of classical plane waves.
In addition to the approximation order q, the general design procedure proposed in this paper involves two parameters: a parameter θ corresponding to the direction of a classical plane wave, -a parameter N = 0 , where N/I can be interpreted as the local wave number of a classical plane wave.
These parameters are used to set (λ 1,0 , λ 0,1 ) = N (cos θ, sin θ). It justifies the name given to the new shape functions: generalized plane waves. Additional constraints greatly simplify both the numerical computations and the analysis of the method. See Definition 3. For a given value of (N, θ), formula (5) together with the Definitions 4 and 5 provide an explicit function ϕ = e P , which satisfies the approximation identity (3). Varying θ then provides different functions ϕ, only as long as N = 0. This condition N = 0 is mandatory to define a set of linearly independent shape functions: they then form a basis of an approximation space E(G, N, p, q).
The set of p shape functions denoted E(G, N, p, q) is defined by {ϕ l } l∈[ [1,p]] . This set of p basis functions is meant to approximate the solution of scalar wave equation. The goal of the theoretical part of this paper is to prove high order approximation properties on such sets of basis functions, provided that enough basis functions are used with respect to the approximation parameter q. Note that the design process does not involve the number p of basis functions. See the hypothesis of the following claim to quantify the relation between the parameters p and q.
The parameter N is then the main degree of freedom to be fixed to define explicitly the approximation space E(G, N, p, q). As will be detailed in Section 2, two different choices will be considered in this paper. A first choice is comprised of setting N = β(G), see Definition 4. It gives a direct generalization of a classical plane waves, since in this case β(G) < 0, so that −β(G) is the local wave number. However, this choice is local since N does depend on G ∈ Ω, and cannot be used if β(G) = 0: it is a classical problem in low frequency regime since the linear independence of such shape functions is damaged, see [14,19].
To overcome this limitation and consider the stationary limit case, a second possibility is to choose one constant and non zero value for N : it will not depend on G anymore. The choice N = I = √ −1 is proposed to address the case β(G) = 0 and ensure the desired interpolation property.
Claim Denote by n ∈ N an interpolation parameter and by G ∈ Ω a point in R 2 . Consider a smooth solution u ∈ C n+1 (Ω) of scalar wave equation (1). Set then -N = 0 the degree of freedom in the design process, -q ≥ n + 1 the order of approximation in (3), -p = 2n + 1 the number of basis functions in E(G, N, p, q).
There exists an approximation u a ∈ Span E(G, N, p, q) of order n + 1 of u in the following sense: The behavior of the constant C N,Ω as N goes to zero is commented in Subsection 3.2. It suggests the need for a parameter N that is bounded away from zero, see Subsections 2.4 and 2.5.
There are two main streams in proving such interpolation results for Helmholtz equation that have been developed in the literature. One of them is based on Vekua theory, which was first translated into English in [13] for functions in R 2 . A more recent introduction to the topic can be found in [2]. Theoretical studies based on this technical tool can be found in [23], and more recently in [25]. In the latter, the case of Helmholtz equation with constant coefficient is explicitly studied and interpolation properties are obtained with explicit dependence with respect to the parameters. However, even if this theory is powerful, in the case of a smooth coefficient it gives no explicit estimates with respect to the different parameters. On the other hand, another method using Taylor expansions was proposed in [5]. Since the design of solutions developed in this paper is based on Taylor expansions as well, this second method will be used here. Section 2 describes precisely the design process, and develops some properties of the resulting approximated functions. It defines two different ways of defining the new functions, called normalizations. Section 3 focuses on the proof of Theorem 1, considering these two normalizations as well. A last section presents a numerical application with a method based on the generalized plane wave basis functions and some numerical results. The numerical test cases are chosen to consider problems linked with reflectometry, a radar diagnostic technique for fusion plasma, see [20] for more details.
Notation. The symbol ∂ z represents the partial derivative with respect to the variable z. The symbol I represents the complex number I = √ −1 to avoid any confusion with the summation index i.

Design and properties of a shape function
This section concerns the design of shape functions that are locally approximated solutions of the scalar wave equation (1). The point G = (x G , y G ) ∈ Ω is fixed, and the design process, based on Taylor expansions, depends on that point G, the degree of freedom N , the number p of shape functions and the order of approximation q. This order of approximation will satisfy q ≥ 1. The case q = 1 corresponds to the simplest generalization of plane waves described previously in (2). The design of the polynomial P starts with the choice of its degree, and then focuses on computing its coefficients to satisfy (3).
Some properties of two different types of shape functions follow. They are meant to be used in the proof of Theorem 1.

Design procedure
Definition 2 Denote by P a bivariate polynomial. The polynomial (∆e P )/e P will be denoted P ∆ , so that P ∆ = ∂ 2 x P + (∂ x P ) 2 + ∂ 2 y P + (∂ y P ) 2 In order to satisfy the local approximation (3), the design is based on a non linear system on the coefficients of P that arises from considering the Taylor expansions in scalar wave equation of β − P ∆ up to the order q. Thanks to Definition 2 it reads The procedure includes choosing the degree of the polynomial and giving an explicit expression to compute the coefficients of the polynomial. These two choices are not independent. A precise analysis of equation (4) leads to choosing the degree of P such that the computation of the coefficients appears to be straightforward.

Remark 1
Since the constant coefficient λ 0,0 does not appear in (4), it is set to zero. This will simplify all the upcoming computations. Moreover, the fact that it does not depend either on G or on β prevents any blow up of the corresponding shape function since then ϕ(G) = e λ0,0 is constant.
The system to be solved to ensure that equation (4) holds has: -N un = (dP +1)(dP +2) 2 − 1 unknowns, namely the coefficients of P except λ 0,0 , equations, corresponding to the cancellation of the terms of degree lower than q in the Taylor expansion of β − P ∆ .
As a result the system is overdetermined if dP < q, and in such a case the existence of a solution is not guaranteed. The idea is then to find the smallest value of dP ≥ q that would provide an invertible system.
The case dP = q is more intricate than the next one, since the q equation stemming from the terms of degree q − 1 have no linear term. It does not -in general -lead to a convenient invertible system. Indeed, in a such case, the system is underdetermined however there is no straightforward way to obtain an invertible system, because of the nonlinearity.
For a given (i 0 , j 0 ), the left part of the figure shows the contributions from P ∆ to the x i 0 y j 0 term in β − P ∆ . The right part shows that λ i 0 +2,j 0 can be explicitly expressed as long as λ k,l are known for all k ≤ i 0 + 1 and l ≤ dP − 2 − k.
As for the case dP = q + 1, the system is underdetermined and the number of additional equations to be imposed to get a square system is N un − N eq = 2q + 2. Moreover, since then the N eq equations of the system that come from (4) actually reads As a consequence, to obtain an invertible system the choice proposed in this paper is to fix the set of coefficients {λ i,j , i ∈ {0, 1}, j ∈ [[0, q + 1 − i]]}. Thus this choice corresponds to the 2q + 3 additional constraints that, together with equations (5), form a square system.

Proposition 1
The system described by (5) together with the additional constraints of fixing the elements of } has a unique solution, given by Proof For any given set of coefficients , the existence and uniqueness of a solution of (5) stems directly from the induction relation (6). See Figure 1.
The last item is the simplest possible choice and is meant to simplify both the numerical computations -by a substantial decrease of basic operations necessary to evaluate a shape function -and the analysis of the method.
Remark 2 (Other possible choices) Other choices to obtain an invertible system would give the same theoretical results. For instance choosing to fix ]} is possible as well. But numerically, as will be seen later on, there is no evidence of the lack of symmetry with respect to the two space variables.

A fundamental property of a generalized pane wave
Since the design and the interpolation study are based on different Taylor expansions, the derivatives of the shape function ϕ are important quantities. Both • the coefficients λ i,j s defining a shape function ϕ • the derivatives of ϕ are here expressed as polynomials with two variables with respect to (λ 1,0 , λ 0,1 ). The following Lemma 1 and Proposition 2 give a description of these quantities with respect to the only non zero coefficients fixed as constraints, namely (λ 1,0 , λ 0,1 ).

Lemma 1 For a given set of coefficients
that are the unique solution of (5) from Lemma 1 can be described as polynomials with two variables in (λ 1,0 , λ 0,1 ) as follows.
The following proof relies on a close examination of the induction formula (5), considered as polynomial with two variables, namely (λ 1,0 , λ 0,1 ). The idea is to track the terms with higher degree.
As a consequence the terms with higher degree appearing in the expression of λ i+2,j have degree at most equal to i. It completes the proof of (7) for i > 2 by induction.
Proposition 2 Suppose θ ∈ R and N ∈ C is such that N = 0. Consider a shape function ϕ = e P constructed in Subsection 2.1. Then for all (i, j) ∈ N 2 such that i + j ≤ q + 1 there is a complex polynomial R i,j such that its total degree satisfies dR i,j ≤ i − 2 and such that The coefficients of R i,j only depend on N and on the derivatives of β.
Remark 3 Since (λ 1,0 ) 2 + (λ 0,1 ) 2 is fixed, none of the polynomial expressions that are at stake can be unique. For instance, any occurrence of (λ 1,0 ) 2 could be replaced by N 2 − (λ 0,1 ) 2 which would change the term of higher degree. This is the reason why R i,j is not unique: see Subsection 2.3 for a different point of view. However, formula (6) from Proposition 1 gives an explicit procedure for the computation of all λ i,j s: this is the crucial point that will be used for practical implementation.
One could have expected the degree of R i,j to be smaller than i + j − 1. The fact that it does actually not depend on j is due to the choice of {λ i,j , i ∈ {0, 1}, i + j > 1} to be zero. The fact that it is smaller than i − 2 is due to the fact that the degree of λ 2,j is 0, since (λ 1,0 ) 2 + (λ 0,1 ) 2 = N 2 is constant with respect to λ 0,1 and λ 1,0 . See Definition 3.
Proof Applying the chain rule introduced Appendix A.2 to ϕ = e P one gets for all (i, j) ∈ N 2 , where p s ((i, j), µ) is the set of partitions of (i, j) with length µ: Now consider such a partition to be given and focus on the degree of the corresponding product term, namely s l=1 (λ i l ,j l ) k l . Thanks to Lemma 1 one can split this product into different terms regarding their degree as polynomials with respect to (λ 1,0 λ 0,1 ). As a result, since Deg k l Deg λ i l ,j l , this quantity is also at most equal to where the two first sums contain at most one term each.
Obviously the leading term in 0). Indeed, as long as a partition contains at least one term such that i l ≥ 2, the resulting degree computed from (10) will contain at least one term k l · 0 or k l (i l − 2), and any of them is at most k l (i l + j l ) − 2; as a consequence the degree computed in (10) is then strictly lower Since the product term corresponding to the partition j(0, 1) + i(1, 0) is (λ 0,1 ) j (λ 1,0 ) i /(j!i!) it completes the proof.

A more algebraic viewpoint
This paragraph presents a more algebraic point of view on Remark 3. Suppose N ∈ C is such that N = 0. The value of (λ 1,0 , λ 0,1 ) gives that P N = (λ 1,0 ) 2 + (λ 0,1 ) 2 − N 2 satisfies P N = 0 for the p different functions of E(G, N, p, q). From then on, considering other quantities as polynomials with two variables in (λ 0,1 , λ 1,0 ) is in fact computing in the quotient ring C[λ 1,0 , λ 0,1 ]/(P N ) of C[λ 1,0 , λ 0,1 ] modulo the ideal generated by P N . For instance, the system (8) reads Of course in this quotient ring, each equivalence class has an infinite number of elements, and all the computations of the previous subsection are performed on elements of these classes. Thus any equality applies to all the elements of the same class. Note that since the ring considered here is the ring of polynomials with two variables, there is no such thing as the Euclidean division. As a result there is nothing like a canonical element of a class used for computations. One can easily see that for q ≥ 4 which gives two possible R 4,1 ∈ C[λ 1,0 , λ 0,1 ] satisfying (9) in Proposition 2.

First family of generalized plane waves
The first type of shape functions corresponds to a local choice since it does depend on G ∈ R 2 .

Definition 4
The β-normalization is defined by choosing N = β(G) in Definition 3, which means setting Remark 4 (Back to classical Plane Waves from the β-normalization) The fact that the quantity (λ 1,0 ) 2 + (λ 0,1 ) 2 is equal to β(G) however gives that the value of β(G) does actually never appear in the expression of the other coefficients explicitley, but only in product terms involving λ 1,0 or λ 0,1 . One can easily check by induction that all the terms appearing in formula (6) are then linear combinations of the derivatives of β. As a consequence, for β = −ω 2 < 0 and for any q ≥ 1, all the coefficients λ i,j such that i > 1 are actually zero, which means that the corresponding function ϕ = e P is nothing more than a classical plane wave. As remarked in the introduction, it is also obvious that for q = 1 this new shape function is again nothing more than a classical plane wave as long as β < 0. This case q = 1 corresponds to the classical fact of approximating a smooth coefficient by its piecewise constant value at the center of the cells.

Second possibility
The fact is that since the terms (λ 1,0 , λ 0,1 ) of the β-normalization are proportional to the square root of β, they will tend to zero with β. A theoretical estimate displayed in Subseciton 3.2 will justify the need for a second normalization. Moreover numerical results show that it causes severe damaging to the conditioning of the discrete UWVF problem if β tends to zero. As a consequence, a second normalization is considered, with a global choice of N independent from G ∈ R 2 .

Definition 5
The constant-normalization is defined by choosing N = i in Definition 3, which means setting It corresponds to N = i.

Remark 5 (Classical Plane Waves and the constant-normalization)
In order to balance Remark 4, note that since the constant-normalization does not depend on β it arises that for β = −ω 2 = −1 the term β(G) appears in higher order terms. For instance it is clear that As a consequence, neither when β( = −1) < 0 is constant nor when q = 1 the shape function designed from the constant-normalization can be a classical plane wave.

Interpolation
The interpolation properties of the set E(G, N, p, q) are defined in Theorem 1. This section is devoted to the proof of this result, which states that, in order to approximate to a given order n + 1 the solution of the scalar wave equation (1) around a point G, a sufficient number p = 2n + 1 of basis functions together with a sufficient approximation parameter q = n + 1 are required. The gradient of the solution is then approximated to the order n.

Theoretical result
This subsection focuses on the interpolation property of the set of basis functions E(G, N, p, q). The sketch of the proof is inspired by the one developed by Cessenat in [5], but it is adapted to the generalized plane wave basis functions. Note that the application to the UWVF is postponed to Section 4.

Definition 6
Suppose that N ∈ C is such that N = 0, n ∈ N * and G ∈ R 2 . For all l ∈ N such that 1 ≤ l ≤ n consider the direction θ l = 2π(l − 1)/(2n + 1), the generalized plane wave ϕ l , κ = −iN ∈ C * and the function e l (x, y) = e Iκ((x−xG) cos θ l +(y−yG) sin θ l ) , which is a classical plane wave if N ∈ iR. The (n + 1)(n + 2)/2 × (2n + 1) matrices M C n and M n are defined as follows: for all (k 1 , k 2 ) ∈ N 2 , such that Their lth columns contain respectively the Taylor expansion coefficients of the functions e l and ϕ l .
The rank of the matrix M C n is computed in Lemma 2, which profits from the fact that the result proved by Cessenat and Després in [5] for κ > 0 is actually still valid for κ ∈ C * . The proof of Theorem 1 relies on Lemma 3 that explicits the link between the matrix M C n and the corresponding matrix M n built with the generalized plane waves. Lemma 2 Suppose that N ∈ C is such that N = 0, n ∈ N * and G ∈ R 2 . There are two matrices: a rectangle matrix P n only depending on β(G) and a square invertible matrix S n only depending on the directions θ l such that S n = P n · M C n . Moreover rk(M C n ) = 2n + 1.
Proof Consider M C n be the matrix introduced in Definition 6 so that for all (k 1 , k 2 ) ∈ N 2 , such that Thanks to the definition of M C n one can check that so that S n is a (2n + 1) × (2n + 1) matrix that is a linear transform of M C n . More precisely, define P n as a (2n + 1) × (n+1)(n+2) 2 matrix such that Then S n = P n · M C n . As a consequence, rk(M C n ) ≥ rk(S n ). The rank of S n is now to be evaluated thanks to the definition of the plane waves e l . Since e l (x, y) = e (iκ)((x−xG) cos θ l +(y−yG) sin θ l ) then (∂ x ± I∂ y ) k e l = (iκ) k (cos θ l ± I sin θ l ) k e l .
Thus S n 's columns are proportional to the one of a VanDerMonde matrix and From the choice of θ l s, for all i = j: z i = z j so that S n is invertible and rk(M C n ) ≥ rk(S n ) = 2n + 1. Since rk(M C n ) ≤ min 2n + 1, (n + 1)(n + 2) 2 = 2n + 1 the proof is then completed.
Lemma 3 Suppose that N ∈ C is such that N = 0, n ∈ N * and G ∈ R 2 . Consider E(G, N, p, q) introduced in Definition 1, together with M n and M C n introduced in Definition 6. Then there is a lower triangular matrix L n , which diagonal coefficients are all equal to 1 and which other coefficients are linear combinations of the derivatives of β evaluated at G, such that M n = L n · M C n .
As a consequence rk(M n ) = rk(M C n ) and both L n and (L n ) −1 are bounded by a constant only depending on β.
The following proof is straightforward considering the feature of the derivatives of ϕ l described in Proposition 2.
Proof From (9) there exists a polynomial R i,j ∈ C[X, Y ] with Deg R i,j ≤ i − 2 such that The coefficients of R i,j do not depend on the basis function considered, but only depends on β and its derivatives evaluated at G. By construction of the classical plane wave e l , one has The numbering of the rows in matrices M C n and M n is set up such that the derivatives of smaller order appear higher in the matrix, which proves (11). Indeed (12) shows that any coefficient of M n is the sum of the corresponding coefficient in M C n plus a linear combination -which coefficients do not depend on the column that is considered but only on β and its derivatives evaluated at G -of terms that appear higher in the corresponding column of M n .
The rank of M n is then equal to the rank of M C n , and L n and (L n ) −1 do only depend on the coefficients of R i,j . As a result they do not depend on the basis functions but only on the coefficient β and its derivatives at G.
Theorem 1 Suppose that n ∈ N and that u is a solution of scalar wave equation (15) belongs to C n+1 . Consider then q ≥ n + 1, p = 2n + 1, and E(G, N, p, q) introduced in Definition 1. Then there are a function u a ∈ Span E(G, N, p, q) depending on β and n, and a constant C N,Ω depending on N , β and n such that for all M ∈ R 2 Proof The idea of the proof is to look for u a = where for the sake of simplicity M l k1k2 stands for the coefficient of M n that corresponds to the term ∂ k1 x ∂ k2 y ϕ l /(k 1 !k 2 !), namely the coefficient (M n ) (k 2 +k 1 )(k 2 +k 1 +1) 2 +k2+1,l , and in the same way B k1,k2 stands for the coefficient (B n ) (k 2 +k 1 )(k 2 +k 1 +1) 2 +k2+1 . The system to be solved is then In order to study the system's matrix, the equations depending on (k 1 , k 2 ) have to be numbered: they will be considered with increasing m = k 1 + k 2 , and with decreasing k 1 for a fixed value of m. Defining the corresponding vector B n ∈ C (n+1)(n+2) 2 , together with the unknown X n = (x 1 , x 2 , · · · , x 2n+1 ) ∈ C 2n+1 , the system now reads Find X n ∈ C 2n+1 such that M n · X n = B n where M n ∈ C (n+1)(n+2) 2 ×(2n+1) is the matrix from Definition 6. Since the system is not square, there is a solution if and only if B n ∈ Im(M n ).
i) The technical point is to prove that rk(M n ) = 2n + 1. It is straightforward from Lemmas 3 and 2.
ii) There exists a subset K ⊂ C (n+1)(n+2) 2 such that Im(M n ) ⊂ K and B n ∈ K. This subspace K is built from the fact that the basis functions are designed to fit the Taylor expansion of the scalar wave equation: All basis functions ϕ l , l ∈ [[1, 2n + 1]], satisfy (−∆ + β)ϕ l = (−P ∆,l + β)ϕ l . From the equation (4) with q ≥ n + 1, it is then straightforward to see that Im(M n ) ⊂ K. The fact that B n ∈ K simply stems from plugging the Taylor expansions of u and β into scalar wave equation.
iii) The dimension of K defined by (14) is 2n + 1. Indeed, one can check -using the same numbering as previously for the equations -that K is defined by n(n + 1)/2 linearly independent relations on C (n+1)(n+2) 2 , so that its dimension is (n + 1)(n + 2)/2 − n(n + 1)/2. As a consequence, from the solution to the system M n · X n = B n that now is known to exist, one can define u a = 2n+1 l=1 x l ϕ l . Thanks to that definition and to the Taylor expansions of u and the ϕ l s it yields Moreover one has the identity X n = (S C n ) −1 P C n (L n ) −1 B n , where (S C n ) −1 P C n is bounded from above by sup e l C n+1 , see Lemma 2, (L n ) −1 is bounded from above by a constant depending only on β and its derivatives from Lemma 3, and B n is bounded by u C n+1 . Since for all l ∈ [[1, 2n + 1]] it yields |x l | ≤ C N,Ω u C n+1 , it turns out to be the first part of (13): At last, the second part of (13) stems from taking the Taylor Lagrange formula of the gradient of u − u a , up to the order n, since n m=0 k1+k2=m That is: for all M = (x, y) ∈ Ω there are ζ 1 , ζ 2 in R 2 on the segment line between M and G such that which indeed leads to the desired inequality.
Remark 6 Some comments on the hypothesis on q are to be found in the next paragraph.
3.2 Estimate of C N,Ω with respect to N → 0 Because the Taylor expansion actually reads blows up when N goes to zero, then so does C N,Ω . As displayed in Proposition 2, each ∂ j x ∂ q−j y ϕ l (G) term is a polynomial with respect to (λ 1,0 , λ 0,1 ) which higher degree term is λ j 1,0 λ q−j 0,1 , that is to say that each derivative term from C N,Ω is a polynomial with respect to N which higher degree term is cos(θ l ) j (i sin(θ l )) q−j N q . As a result these terms are not zero and tend to zero as N tends to zero, at most as N q .
On the other hand since X n satisfies M C n · X n = (L n ) −1 B n , one can describe the asymptotic behavior of the x l s with respect to N as follows. Consider the (2n + 1) × (2n + 1) square system defined as: for all k ∈ N such that 0 The right hand side of this system is independent on N while its determinant is N n(n+1) det(S n ). As a result the coefficients x l behave (as long as it is non zero) as C N,Ω /N n(n+1) as N tends to zero. As a consequence, C N,Ω as well as C N,Ω blow up at least as 1/N (n−1)(n+1) as N goes to zero.

Numerical validation
In order to validate the Theorem 1, each of the numerical validation case is computed, for a given value of n, setting q = n + 1 and p = 2n + 1. The test case considered is β(x, y) = x − 1, to approximate the exact solution u e (x, y) = Airy(x)e iy . See [20] for the physical motivation of this test case: its main interest is that the coefficient vanishes along the line x = 1, which represents a plasma cut-off that reflects incoming waves. Of course since the theoretical results give local approximation properties, the validation procedure itself will be local as well. As stated in the theorem u e can be approximated by a function u a that belongs to the approximation space E (G, N, p, q), space that is built with either the βnormalization or the constant-normalization.
The idea is to follow the error max |u − u a | on disks with decreasing radius h in order to observe the order of convergence with respect to h. Several different cases are proposed to validate the theoretical order of convergence, and an additional case concerns the behavior of the basis functions designed with the β-normalization as the approximation point gets closer to the cut-off.

In the propagative zone
The point G = [−3, 1] is in the propagative zone. Then concentric disks are centered on G with radius h = 1/2 k , increasing the value of k. Following the theorem, the expected order of convergence is n + 1. Figure 2 displays computed convergence results that fit perfectly the theoretical result. A set of p = 11 classical Plane Waves is used as a control case, since p = 11 is the highest number of basis functions used in the different cases with the Generalized Plane Waves. Note that machine precision is reached in some cases.

In the non propagative zone
The point G = [2,1] is in the non propagative zone. Again concentric disks are centered on G with radius h = 1/2 k , increasing the value of k, and the expected order of convergence is n + 1. There is no classical Plane Wave that can be computed here since β(G) > 0. Figure 3 displays computed convergence results that fit perfectly the theoretical result as well. Again machine precision is reached in some cases.  Figure 4 show that the β normalized Generalized plane waves give a high order approximation of u even getting closer to the vanishing line x = 1, as long as h is not too small. Note that there does not seem to be a significant difference between the two type of functions. This is observed on numerical results even if there is no corresponding theoretical result.
Another possibility is to compare the influence of two parameters : the size of the disk h and the distance d between G and the line x = 1. The error e depends on both parameters, so one can write e(h, d). Figure 5 displays the error computed for h and d convergence with the β-normalization. The h convergence is clearly damaged for decreasing values of d. This is linked to the low frequency   As to the Generalized Plane waves with constant-normalization, the theoretical results show that their interpolation property holds along the cut-off as well as anywhere else in the domain.
As Figures 2 and 3, Figure 6 displays results that fit perfectly the theoretical result. It is an example of efficient approximation of the exact solution u e along the cut-off.

A generalized plane wave numerical method for smooth non constant coefficients
Associating to the scalar wave equation a boundary condition, consider now the following problem.  inscribed in Ω k . The mesh is such that ∃σ such that h k ≤ σρ k . The refinement parameter or mesh size parameter h is then defined by h = max h k . The terminology "regular mesh" comes from [6]. On a given element of the mesh Ω k , the center of gravity is denoted G k , and p(k) is the number of basis functions on Ω k . These basis functions {ϕ l k } l∈[ [1,p(k)]] are set to be zero on Ω\Ω k . This process defines a set of basis functions on Ω, namely The function space for the UWVF is denoted V as It defines a norm: X = (X, X). In particular for any operator A ∈ L(V ), the norm is As a consequence any element of V is actually defined on the edges of the mesh elements.

Adapted Ultra Weak Variational Formulation
Let us start from the classical ultra weak Variational formulation. The test functions space is defined by Theorem 2 Let u ∈ H 1 (Ω) be a solution of problem (15) such that ∂ ν k u ∈ L 2 (∂Ω k ) for any k. Let γ > 0 be a given real number. Then X ∈ V defined by for any e = (e k ) k∈[[1,N h ]] ∈ H. Conversely, if X ∈ V is solution of (17) then the function u defined locally by is the unique solution of the problem (15).
This result is classical in the context of UWVF. We refer to [5,4,15,17,18]. Even though the formalism and the use of approximated solutions as test functions has been described in [20], the compact formulation is once again introduced hereafter, together with its adaptation to the generalized plane waves.
Definition 7 For any f ∈ L 2 (Ω), let E f be the extension mapping defined by : Also define E which is the homogeneous extension operator with vanishing right hand side, namely E = E 0 .
Let F and Π be the mapping defined by Π : Then, see [5], the problem (17) is equivalent to where the right hand side b ∈ V is given by the Riesz theorem The classical discretization process would be to consider V h ⊂ V of finite dimension and solve the problem (19) on V h . However, since the generalized plane waves do not belong to H, this formulation (19) has to be adapted to the new basis functions. Indeed, there is no such thing as a unique equation satisfied by the basis functions to define an extension mapping from V to H. The extension mapping is then defined as a one-to-one function from a subset of V to the set of basis functions E.
Some additional notations will be useful. The local discrete space is It is therefore convenient to define the trace v l k ∈ V by V l k = (−∂ ν + iγ)ϕ l k on L 2 (∂Ω k ), and V l k = 0 on L 2 (∂Ω k ′ ) k ′ = k. An equivalent way to define W k and V q could be Next define what are the generalizations of operators E and F in this context. Let E q belong to L(V q , where P ∆,l = ∂ 2 x P l + ∂ 2 y P l + (∂ x P l ) 2 + (∂ y P l ) 2 . Each of these equations are different, but they satisfy for each l: P ∆,l − β = O(h q ).
With these notations and definitions, the abstract UWVF with generalized plane waves is defined as follows.

Definition 8 (UWVF method with generalized plane waves)
with the right hand side given by

Interpolation interpretation
P h denotes the orthogonal projection in V on V q .
Proposition 3 Let u be a solution of a homogeneous scalar wave equation problem. Assume that u is of class C n+1 with n ≥ 1. Let X ∈ V satisfy X |∂Ω k = (−∂ ν k + iω)u |∂Ω k . The number of basis functions Z k,l = (−∂ ν + iγ)ϕ l k per element Ω k is fixed p = 2n + 1. Let us assume, for the sake of the simplicity of the proof, that their directions are fixed once for all. Then ∃C > 0 depending on n and the problem's data γ such that Proof Applying Theorem 1 on every element of the mesh, one gets ∃u a ∈ E, u a = k,l x l k ϕ l k such that ∀ − → x ∈ Ω: so that for h small enough The results then stems from the fact that, for a regular mesh, the total number of elements can be bounded by C/h 2 .

Numerical results
The numerical method presented in Subsection 4.1.2 has been implemented in two dimensions.
The main difference with a classical UWVF code lies in the need for a quadrature formula. Indeed, there is no exact integration formula when the basis functions are no more classical Plane Waves. Results corresponding to two Newton-Cotes formulas will be compared: a Boole formula with five points, a Weddle formula with seven points and a formula with ten points.

Comparing h convergence results on a test case
Consider the system −∆u + (x − 1)u = 0,  Figure 7. Figure 8 evidences the limitation of low order quadrature formulas. It compares results obtained with Boole (5 points), Weddle (7 points) and the 10 points quadrature formulas. All the results are computed with p = 2n + 1 basis functions per element and with the order of approximation q = n + 1. The results obtained with n = 3 are plotted with different cross-marks, whereas those obtained with n = 4 are plotted with different square-and disk-marks. For n = 3 the choice of quadrature makes obviously no difference. However, to obtain a higher order convergence, there is a price to pay regarding quadrature formula. Indeed, for n = 4 the situation is different : the computations performed with Boole formula do not reach the high convergence rate obtained with both Weddle and the 10 points formulas. This is due to the pollution introduced by the low order quadrature formula. Figure 9 compares convergence results with respect to the mesh size h, all the integral being computed with a Weddle quadrature formula. Results obtained with the β-normalization and with the constant-normalization are displayed respectively on the left and on the right, at the same scale. As suggested by the interpolation theorem, all the results are computed with q = n + 1 and p = 2n + 1 where n ∈ N. The left hand part of the figure displays results computed with the β-normalization. For p = 11, the ill-conditioning of the system matrix prevent the method from converging. It is remarkable that even though the basis functions with the β-normalization do not exist along the cut-off, this does not prevent the method from converging. The right hand side of the figure displays results computed with the constant-normalization. It shows that even if the 10 100 L2 error at the center of the mesh cells 1/h q=2 and p=3 q=3 and p=5 q=4 and p=7 q=5 and p=9 q=6 and p=11 Fig. 9 h-convergence results displayed for β-normalization on the left and constant-normalization on the right, both using Weddle quadrature formula to compute the integrals. The relative discrete L 2 error computed at the centers of the mesh cells is represented depending on the inverse of the mesh size h. The results are computed with p = 2n + 1 and q = n + 1, for n between 1 and 5.
constant-normalization ensures high order interpolation along the cut-off, the numerical method obtained by described in this section provides much smaller error when used with the β-normalized basis functions than with the constant-normalized basis functions. It also seems that the method with the constant-normalization is inadequate when the number p of basis functions per element increases.

A first physical test case
This case was proposed by Stéphane Heuraux as a second step toward reflectometry applications, for which the presence of a plasma cut-off is crucial. It models a wave sent in a plasma by an antenna from the wall of a fusion reactor, see [21]. The antenna is represented by a wave guide added outside the reactor on a wall plus a horn inside the reactor.
The geometry is described in Figure 10. The domain Ω is a L × L square, the width and length of the waveguide are l 0 and 4l 0 . The size of the domain is set to be L = 50l 0 , where l 0 is the wavelength of the incoming signal in the horn. The cut-off is set at x = 40l 0 . The heterogeneous  medium is modeled by the coefficient One gets a wave propagating from the wave guide through the horn toward the right end of the domain, reflected by the cut off situated at x = 4. See Figure 11. This result was computed with the normalization N = √ β.

Conclusion
A procedure to design a set of generalized plane waves that are locally approximated solution of the scalar wave equation with smooth non constant coefficient has been successfully developed. Both theoretical and numerical results evidence the high order approximation. This procedure can easily be generalized to many differential operators, as described in [21]. Moreover a natural idea would be to extend the generalization process from the phase to the amplitude of plane waves, by considering a looking for a shape function as ϕ = Qe P where P and Q are two polynomials. These generalized plane waves have been coupled with the UWVF to obtain a numerical method adapted to problems with smooth vanishing coefficients. This resulting numerical method has been described, and the first numerical results are promising but no theoretical convergence is available yet. The behavior of the method along a cut-off requires further investigation, since the comparison between numerical results obtained with the two normalizations of the basis functions is very not yet explained theoretically.