Xavier Antoine

Christophe Geuzaine

Phase Reduction Models for Improving the Accuracy of the Finite Element Solution of Time-Harmonic Scattering Problems I: General Approach and Low-Order Models

Keywords: Helmholtz equation, acoustic scattering, short-wave problem, finite element method, pollution, accuracy, On-Surface Radiation Condition method PACS:

This paper introduces a new formulation of high-frequency time-harmonic scattering problems in view of a numerical finite element solution. It is well-known that pollution error causes inaccuracies in the finite element solution of short-wave problems. To partially avoid this precision problem, the strategy proposed here consists in firstly numerically computing at a low cost an approximate phase of the exact solution through asymptotic propagative models. Secondly, using this approximate phase, a slowly varying unknown envelope is introduced and is computed using coarser mesh grids. The global procedure is called Phase Reduction. In this first paper, the general theoretical procedure is developed and low-order propagative models are numerically investigated in detail. Improved solutions based on higher order models are discussed showing the potential of the method for further developments.

Introduction

The purpose of this paper is to investigate simple evolution models for exterior scattering problems around convex obstacles, and use these models to reduce the pollution error in the finite element solution of time-harmonic scattering problems.

The pollution error finds its origin in the oscillatory nature of the solution of such scattering problems, and is deeply linked to the loss of stability of the Helmholtz operator for high wavenumbers. As evidenced by Babuska et al. [START_REF] Babuška | Is the pollution effect of the fem avoidable for the helmholtz equation considering high wave numbers[END_REF][START_REF] Ihlenburg | Finite element solution of the Helmholtz equation with high wave number. I. The h-version of the FEM[END_REF][START_REF] Ihlenburg | Finite element solution of the Helmholtz equation with high wave number. II. The h-p version of the FEM[END_REF], the contribution of the pollution error to the overall error in the finite element solution of scattering problems thus becomes predominant at high wavenumbers k. This requires to increase the mesh density faster than the wavenumber to maintain a prescribed accuracy [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF], and leads to intractable computational costs at very high frequencies, i.e., when the wavelength λ = 2π/k is much smaller than the characteristic length of the geometrical structure.

Several techniques have been proposed during the last decade to overcome this problem: hp finite element methods, stabilized Galerkin finite element methods, multi-scales techniques, wave-based discretization techniques, etc. We cannot mention here the numerous papers related to all these improvements and instead refer the reader to the recent review by Thompson [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF] for further details and references.

A common point between all these techniques lies in the fact that information about the scattering problem is incorporated into the finite element procedure. One popular approach is to replace the standard polynomial basis functions by plane wave functions [START_REF] Laghrouche | Short wave modelling using special finite elements[END_REF][START_REF] Giladi | A hybrid numerical asymptotic method for scattering problems[END_REF]. Another viewpoint [START_REF] Baumeister | Analysis of sound propagation in ducts using the wave envelope concept[END_REF][START_REF] Baumeister | Finite-difference theory of sound propagation in ducts with variable axial dependence[END_REF][START_REF] Baumeister | Finite difference theory for sound propagation in lined duct with uniform flow using wave envelope concept[END_REF][START_REF] Baumeister | Time-dependent wave envelope finite difference analysis of sound propagation[END_REF][START_REF] Astley | Finite element formulations for acoustical radiation[END_REF][START_REF] Astley | Wave envelope and infinite elements for acoustical radiation[END_REF][START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case[END_REF][START_REF] Geuzaine | An amplitude formulation to reduce the pollution error in the finite element solution of time-harmonic scattering problems[END_REF][START_REF] Turkel | Improved accuracy for the helmholtz equation in unbounded domains[END_REF], which we adopt in this paper, is to approximate the phase of the solution and use this phase to reformulate the problem in terms of a slowly oscillatory envelope. If we denote the solution of the original scattering problem by u, this approach thus involves a two-fold solution process:

1) find an approximation φ of the phase φ of u in the whole computational domain; 2) use φ to solve the scattering problem in terms of a new, slowly varying unknown.

This has the advantage that the resulting formulation can be easily coded into a classical finite element solver and does not require any integration of new basis functions. Moreover, this technique is not restricted to the finite element method, and can be used in other numerical schemes like finite difference, integral or spectral methods.

In this paper, we propose an original approach for solving point 1) above, i.e., to find an approximation of the phase of the highly-oscillatory solution in the computational domain. Our approach is based on three main ingredients:

1.a) first, we build an initial condition on the surface of the scatterer (that we call the starter ); 1.b) then, we construct an approximate pseudo-differential evolution operator (that we call the propagator ) to extend the starter into the exterior computational domain; 1.c) and finally, to reduce the computational cost, we localize the propagator using purely differential operators.

We then solve an alternative variational problem, equivalent to the original one, by using the phase φ of the approximate solution. The new unknown is expected to oscillate less than the original one, which permits to significantly reduce the pollution error, or, alternatively, to obtain the same accuracy as the classical finite element method with coarser discretizations. We call the resulting technique the Phase Reduction procedure (or PR for short). It can be seen as a numerical extension of asymptotic approaches like the WKB approximation [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF], with the added advantage that, provided that the finite element mesh resolves the amplitude, the method results in a convergent solution of the original problem-even if the phase approximation is inaccurate.

To some extend, the PR procedure can also be seen as a general extension of the concept of Wave Envelope introduced by Baumeister [START_REF] Baumeister | Analysis of sound propagation in ducts using the wave envelope concept[END_REF][START_REF] Baumeister | Finite-difference theory of sound propagation in ducts with variable axial dependence[END_REF][START_REF] Baumeister | Finite difference theory for sound propagation in lined duct with uniform flow using wave envelope concept[END_REF][START_REF] Baumeister | Time-dependent wave envelope finite difference analysis of sound propagation[END_REF] for treating particular guided waves propagation problems, and later used by Astley and Everman [START_REF] Astley | Finite element formulations for acoustical radiation[END_REF][START_REF] Astley | Wave envelope and infinite elements for acoustical radiation[END_REF] for acoustical radiation.

To show the potential our approach, in this paper we mainly consider the lowest order (and computationally simplest) models associated with points 1.a), 1.b) and 1.c) above. Even with these low order approximations, we will show that significant improvements can be achieved on general scattering problems around convex obstacles. Improvements that can be expected at very high frequencies using higher order models are also briefly presented.

The paper is structured as follows: in Section 2 we outline our phase reduction strategy on a simple one-dimensional scattering problem. In Section 3 we present the usual variational formulation for higher-dimensional scattering problems, using a Bayliss-Gunzburger-Turkel-like radiation condition to truncate the infinite domain. The phase reduction procedure for higherdimensional scattering problems is then detailed in Section 4, and Section 5 provides numerical simulations to show the improvements related to the phase reduction approach. In Section 6 we show the improvements that can be expected at very high frequencies by using higher order models.

2 The one-dimensional case: a simple and explicit example Assume that we wish to solve the one-dimensional Neumann scattering problem

∂ 2 x u + k 2 u = 0, in Ω b =]0, 1[, ∂ n Γ u = ik, at Γ = {0}, ∂ x u -iku = 0, at Σ = {1}. (1)
The above boundary-value problem corresponds to the solution of the scattering of a one-dimensional incident plane wave by the left half-space, introducing an exact transparent boundary condition on the fictitious boundary Σ = {1}. This non-reflecting boundary condition is given through the Dirichlet-to-Neumann (DtN) operator Λ = ik on Σ by the relation:

∂ x u = Λu.

The pollution problem

It is well-known since the pioneering works of Babuska and Ihlenburg [START_REF] Ihlenburg | Finite element solution of the Helmholtz equation with high wave number. I. The h-version of the FEM[END_REF][START_REF] Ihlenburg | Finite element solution of the Helmholtz equation with high wave number. II. The h-p version of the FEM[END_REF][START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF][START_REF] Ihlenburg | Finite Element Analysis of Acoustic Scattering[END_REF] that the finite element approximation of the solution u(x) = e ikx of problem (1) suffers from numerical pollution, especially for large wavenumbers k. Essentially, this is related to the inadequacy of the polynomial finite element basis used to represent the unknown wavefield. To clarify this problem, we consider the covering of the computational bounded domain Ω b using N h uniform segments of length h, setting h = 1/N h . We report in Figure 1 the exact solution u and its linear finite element approximation u h for k = 40. We denote by n λ = λ/h the density of discretization points per wavelength λ. As we can see in Figure 1, a density n λ = 8 is too small to compute an accurate approximate solution. The usual rule of thumb of "10 points per wavelength" is generally not sufficient and a mesh refinement is required, resulting in an added computational cost. This loss of accuracy is known as pollution and translates the accumulation of numerical phase error over the computational domain inherent in the FEM. Classical error estimates can be derived (see e.g. [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF][START_REF] Ihlenburg | Finite Element Analysis of Acoustic Scattering[END_REF]) to quantify this pollution error, which becomes predominant compared to the approximation error for large wavenumbers.

A possible strategy to reduce the pollution error

To partially avoid the pollution error, we investigate in this paper the procedure based on the two following successive steps: 1) The first step consists in determining an approximate solution

u = A(x)e ik φ(x) (2)
in the whole computational domain Ω b as the solution of an evolution equation (the propagator) in the x-direction, with an initial condition to be determined on Γ (the starter); once the approximate field (2) is known in the volumetric computational domain, we extract its phase φ(x).

2) The second step consists in solving a new variational formulation for the approximate slowly varying complex-valued envelope A(x) of the true solution u(x) = A(x)e ik φ(x) , setting

A(x) = A(x)e ik(φ(x)-φ(x)) , (3)
where A and φ respectively designate the true real-valued amplitude and phase of u:

u(x) = A(x)e ikφ(x) . (4)
Following this strategy, the pollution error is expected to be reduced since an approximate a priori phase determination leads to the computation of a slowly varying field A.

Computing the approximate solution u

The solution of step 1) is split in two distinct problems:

1.a) the construction of the starter; 1.b) the obtention of the propagator.

Point 1.a) is solved through the following considerations. From the nature of the scattering problem (1), the Neumann boundary condition is given. To get an approximate starting field, we apply the transparent boundary condition directly at the physical boundary Γ. In other words, we consider the so-called On-Surface Radiation Condition (OSRC) [START_REF] Antoine | Advances in the On-Surface Radiation Condition Method: Theory, Numerics and Applications[END_REF][START_REF] Kriegsmann | A new formulation of electromagnetic wave scattering using the on-surface radiation condition method[END_REF] solution in the onedimensional case. In our simple situation, this condition is exact and gives the initial wavefield u 0 as

u 0 = u(0) = (ik) -1 ∂ n Γ u = 1, at Γ. ( 5
)
Interpreting x as a time variable, the equation ( 5) is an initial condition.

In the case of a Dirichlet boundary condition u = f on Γ, the initial condition is directly given by u 0 = f on Γ and hence does not require the help of the OSRC technique [START_REF] Geuzaine | An amplitude formulation to reduce the pollution error in the finite element solution of time-harmonic scattering problems[END_REF].

Let us now consider point 1.b). To construct the propagator, we rewrite the Helmholtz equation under its (exact) factorized form

(∂ 2 x + k 2 )u = (∂ x + iP -)(∂ x + iP + )u = 0, for x ≥ 0, (6)
setting P ± = ∓k. The solution of the evolution equation

∂ x u + iP + u = 0 (7)
corresponds to the forward propagative wave and is considered as our propagator. The solution u of [START_REF] Astley | Finite element formulations for acoustical radiation[END_REF] with the initial condition (5) can be directly obtained as u = e ikx u 0 = e ikx .

The procedure consisting in writing an approximate forward propagating model is closely related to the techniques employed in optics, underwater or electromagnetic wave propagation and better known under the denominations respectively of Beam Propagation Methods (BPM) [START_REF] Lu | Some techniques for computing wave propagation in optical waveguides[END_REF] or parabolic equations [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF]. As discussed next, the situation is more complicated in higher dimensions since the wavefield propagates from a (usually) curved surface Γ.

Computing the true solution u

The second step of the algorithm consists in solving an equation for A defined by: A = e -ik φ u. (8) Again, in the considered configuration, we know that u(x) = e ikx and consequently that A(x) = 1 since φ(x) = x. Therefore, A is a slowly varying envelope and the polynomial basis functions lead to an optimal numerical computation. However, to explain in detail the real strategy considered for higher dimensions, we rather proceed as follows. The variational formulation of problem (1) is given by: find u ∈ H 1 (Ω b ) such that

Ω b ∂ x u ∂ x v -k 2 uv dx -ik(uv)(1) = -(gv)(0), ∀v ∈ H 1 (Ω b ). ( 9
)
Since A is defined by [START_REF] Astley | Wave envelope and infinite elements for acoustical radiation[END_REF], with φ known from step 1), we replace u in the weak formulation (9) using u = e ik φ A, and we choose some test-functions v defined by v = e ik φ B. Using these relations, we get the new formulation: find

A ∈ H 1 (Ω b ) such that Ω b ∂ x A ∂ x B + ik∂ x φ(A∂ x B -B∂ x A) -k 2 (1 -|∂ x φ| 2 )A B dx -ik(A B)(1) = -(f B)(0), ∀B ∈ H 1 (Ω b ), ( 10
)
where f (0) = e -ik φ(x) |x=0 g(0) = ik. This formulation is nothing else than the weak formulation of the initial boundary-value problem (1) with A as the new unknown. The weak formulation ( 10) is easy to implement in a finite element code and only requires some elementary modifications of the initial formulation [START_REF] Baumeister | Analysis of sound propagation in ducts using the wave envelope concept[END_REF]. The FEM solution to formulation [START_REF] Baumeister | Finite-difference theory of sound propagation in ducts with variable axial dependence[END_REF] will be referred to as Phase-Reduction FEM (PR-FEM) in the sequel. We present in Figure 1 the computed discrete finite element solution based on A using [START_REF] Baumeister | Finite-difference theory of sound propagation in ducts with variable axial dependence[END_REF]. We observe that we avoid any pollution fundamentally because we compute a non oscillatory solution. The polynomial representation is therefore optimal.

The goal of the following sections is to extend this approach to scattering problems in higher dimensions, where the computation of u can only be realized in an approximate way. In this case, we cannot completely avoid the pollution error in the finite element solution, but we can expect to reduce it significantly at high wavenumbers. Furthermore, the extension to variable wavenumbers can also be considered. Indeed, the construction based on pseudodifferential operators and presented in the next section enable one to consider Helmholtz equations with variable coefficients.

3 The classical variational formulation for higher dimensions

We investigate the numerical solution of the time-harmonic acoustic scattering problem of a plane wave u inc (x) = e ikα•x , |α| = 1, by a sound-hard or a soundsoft convex obstacle Ω -⊂ R d , d > 1, with a C 1 closed boundary Γ. We set Ω + := R d \Ω -as the exterior domain of propagation associated with Ω -. The spatial variable is denoted by x = (x 1 , ..., x d ). The boundary value problem reads:

∆u + k 2 u = 0, in Ω + , ∂ n Γ u = g := -∂ n Γ u inc or u = f := -u inc , on Γ, lim |x|→+∞ |x| (d-1) (∇u • x |x| -iku) = 0. ( 11
)
The operator ∆ is the Laplacian operator and a • b designates the hermitian inner product of two d-dimensional complex-valued vector fields a and b. The outwardly directed unit normal vector to Ω -is n Γ . In a suitable functional setting, this problem is known to be well-posed [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF].

Let us begin by introducing the classical weak formulation and the finite element method to solve [START_REF] Baumeister | Finite difference theory for sound propagation in lined duct with uniform flow using wave envelope concept[END_REF]. We consider a smooth convex fictitious boundary Σ enclosing the scatterer Ω -and we set Ω b as the bounded computational domain delimited by Γ and Σ. Since our aim here is not to investigate all the available solutions to bound the computational domain, we only consider the second-order Bayliss-Gunzburger-Turkel-like (BGT2-like) Artificial Boundary Condition (ABC) derived in [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF]. Other truncation techniques like PML [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Harari | Analytical and numerical studies of a finite element PML for the Helmholtz equation[END_REF] could also be used and adapted as easily to the proposed strategy. The BGT2like ABC on Σ is given by

∂ n Σ u = Bu, on Σ, ( 12
)
where B is a second-order symmetrical tangential boundary operator and n Σ designates the outwardly directed unit normal vector to Σ. To simplify the notations, we keep on denoting by u the solution of the approximate model with the BGT2-like ABC (see system ( 16)) even if this function is obviously different from the solution of the initial problem [START_REF] Baumeister | Finite difference theory for sound propagation in lined duct with uniform flow using wave envelope concept[END_REF]. The boundary operator B is given by

Bu = div Σ (A∇ Σ u) -βu, on Σ, (13)
where the operator div Σ is the surface divergence of a tangential complexvalued vector field and ∇ Σ is the surface gradient operator of a complex-valued scalar surface field, all these quantities being defined over Σ. Furthermore, the operator A is the complex-valued tensor field given by

A = - 1 2ik (I + iR k ) -1 ( 14
)
and β is the complex-valued scalar function

-ik + H + i 2k (1 + 2H k ) -1 (K -H 2 ) - ∆ Γ H 4k 2 . ( 15
)
In the above notations, the operator I is the identity operator of the tangent plane and R is the curvature tensor. We denote by H and K respectively the mean and Gauss curvatures. The Laplace-Beltrami operator is ∆ Σ u = div Σ (∇ Σ u). It results that the truncated boundary-value problem is given by

∆u + k 2 u = 0, in Ω b , ∂ n Γ u = g or u = f, on Γ, ∂ n Σ u = Bu, on Σ. ( 16
)
For the Neumann problem, the variational formulation consists in computing [START_REF] Geuzaine | A Galerkin projection method for mixed finite elements[END_REF] and the linear form ℓ appearing in the right-hand side is

u ∈ H 1 (Ω b ) such that a(u, v) = ℓ(v), ( 17
) for any test-function v ∈ H 1 (Ω b ). The sesquilinear form a(•, •) is defined by a(u, v) = (∇u, ∇v) 0,Ω b -k 2 (u, v) 0,Ω b + (A∇ Σ u, ∇ Σ v) 0,Σ + (βu, v) 0,Σ
ℓ(v) = -(g, v) 0,Γ . (19)
We define (•, •) 0,D as the hermitian inner product of two complex-valued squareintegrable functions u and v of L 2 (D)

(u, v) 0,D = D uvdD, (20)
where D specifies the integration domain. In the case of two complex-valued vector fields u and v defined on D, the inner product is

(u, v) 0,D = D u • vdD. ( 21
)
The finite element solution consists in introducing a covering Ω h of Ω b using some tetrahedrons K: Ω h = ∪ K∈K h K, where K h designates a triangulation of the domain. The corresponding interpolated boundaries associated with Γ and Σ are respectively denoted by Γ h and Σ h . The p-finite element version of (17) yields the discrete formulation: find

u h ∈ V h such that a h (u h , v h ) = ℓ h (v h ), (22)
for any test-function v h of V h , setting

a h (u h , v h ) = (∇u h , ∇v h ) 0,Ω h -k 2 (u h , v h ) 0,Ω h + (A h ∇ Σ h u h , ∇ Σ h v h ) 0,Σ h + (β h u h , v h ) 0,Σ h (23) and ℓ h (v h ) = -(g h , v h ) 0,Γ h . ( 24
)
The classical finite element space of order p is given by

V h := {v h ∈ C 0 (Ω h )/v h |K ∈ P p (K), ∀K ∈ K h }, (25)
where P p denotes the space of polynomials of degree less than or equal to p. The approximate fields A h and β h are computed by some suitable schemes based on the surface mesh. We refer to [START_REF] Antoine | Fast approximate computation of a time-harmonic scattered field using the on-surface radiation condition method[END_REF] for implementation details. Finally, the solution of ( 22) leads to the solution of a linear system

[a h ]u h = b h , (26)
where [a h ] is a non-hermitian complex-valued matrix of size N h × N h , u h and b h are two complex-valued vectors of C N h , and N h is the number of degrees of freedom associated with the finite element approximation. The whole procedure is referred to as FEM in what follows.

Concerning the Dirichlet problem, a similar weak formulation can be obtained in a classical way. We do not detail here this point which is immediate.

The Phase Reduction formulation for higher dimensions

We now detail our strategy to solve the two successive steps outlined in the introduction for higher dimensions. Point 1) is developed in Section 4.1 (more specifically, Section 4.1.1 deals with the proposition of a starter, i.e., Point 1.a), while Section 4.1.2 concerns the construction of the propagator, i.e., Point 1.b)). For now we only consider the lowest order propagator and therefore Point 1.c) is direct-improvements will be considered later in Section 6. Point 2) is developed in Section 4.2.

We focus the presentation on the sound-hard case, i.e., on the problem with a Neumann boundary condition. Like in the one-dimensional case, the PR-FEM for the Dirichlet problem is obtained in a similar way, with the added simplification that the starter is obtained immediately from the boundary condition.

Preliminary results for the two-dimensional Dirichlet scattering problem were presented in the short paper [START_REF] Geuzaine | An amplitude formulation to reduce the pollution error in the finite element solution of time-harmonic scattering problems[END_REF].

Computing the approximate solution u

Construction of the initial condition u 0 (starter)

In the one-dimensional case, we have seen that the construction of the initial condition is essentially based on the knowledge of the Dirichlet-to-Neumann (DtN) operator Λ on Γ. This approach, called the On-Surface Radiation Condition method, has been formally introduced in the middle of the eighties by Kriegsmann, Taflove and Umashankar [START_REF] Kriegsmann | A new formulation of electromagnetic wave scattering using the on-surface radiation condition method[END_REF] for the computation of an electromagnetic wave by a perfectly conducting body. Since then, numerous im-provements and generalizations have been developed both in acoustics and electromagnetism. We restrict our presentation here to the use of the OSRC method applied to the construction of the initial condition at high frequencies. We refer to [START_REF] Antoine | Advances in the On-Surface Radiation Condition Method: Theory, Numerics and Applications[END_REF] for a recent presentation of the OSRC approach and its applications. An interesting aspect of the OSRC approach is that the construction of the approximate solution is not based on only taking into account the Helmholtz operator but rather on considering the complete scattering boundary value problem. In this sense, the OSRC method provides an approximate solution which has a behavior close to the true solution and can be seen as a numerical asymptotic solution for a general scattering problem.

The accuracy and efficiency of these techniques are deeply related to the way of constructing the approximation of the DtN map. Various choices exist and are in some way all related to a subjacent (local or microlocal) asymptotic analysis. Since we wish to solve a high frequency scattering problem, we choose the square-root OSRC proposed in [START_REF] Antoine | An improved surface radiation condition for high-frequency acoustic scattering problems[END_REF]. For a given normal derivative trace on Γ, this operator expresses an approximation ϕ of the exact trace of the solution u through the relation

∂ n Γ u = Λϕ, on Γ, ( 27
)
where Λ is the pseudo-local OSRC operator given by

Λ = ik √ 1 + X. (28)
In our case, we will set u 0 = ϕ. The symmetrical partial differential operator X is

X = div Γ (k -2 ε ∇ Γ •), (29)
defining the complex wavenumber k ε as

k ε = k + iε, with ε = 0.4k 1/3 H 2/3
Γ , where H Γ is the mean curvature on Γ. Since the operator is non-local like in the integral equation approach, its numerical inversion is computationally expensive. However, its simulation can be efficiently carried out by using a suitable complex Padé approximant of the square-root √ z of a complex number z. The principal determination of the square-root is considered with a branch-cut along the negative real axis. Following [START_REF] Antoine | An improved surface radiation condition for high-frequency acoustic scattering problems[END_REF], we use the approximation described in [29]

√ 1 + z ≈ e iθ/2 R Np (e -iθ (1 + z) -1) = A 0 + Np j=1 A j z 1 + B j z , (30)
setting

A j = e -iθ/2 a j (1 + b j (e -iθ -1)) 2 and B j = e -iθ b j 1 + b j (e -iθ -1) , (31)
and A 0 = e iθ/2 R Np (e -iθ -1). The function R Np designates the usual real Padé approximant of order N p of the square-root defined by

√ 1 + z ≈ R Np (z) = 1 + Np j=1 a j z 1 + b j z , (32)
with the coefficients

a j = 2 2N p + 1 sin 2 ( jπ 2N p + 1 ) and b j = cos 2 ( jπ 2N p + 1 ). ( 33
)
In the sequel, we take the optimal values N p = 8 and θ = π/3 considered in [START_REF] Antoine | An improved surface radiation condition for high-frequency acoustic scattering problems[END_REF]. The important feature in using a Padé interpolant is that the approximate and accurate application of the operator (28) to a given surface function ψ can be realized through the solution of (N p + 1) complex-valued coefficients Helmholtz-type partial differential equations over Γ. Indeed, we can reformulate the problem of evaluating Ψ = Λψ as the application of the Padé-type operator to a given datum ψ applying the following procedure: solve the N p surface partial differential equations with respect to the auxiliary unknowns ψ j through the variational formulation

Γ B j k 2 ε ∇ Γ ψ j • ∇ Γ ϕ -ψ j ϕdΓ + Γ ψϕdΓ = 0, 1 ≤ j ≤ N p , (34)
and next compute Ψ variationally using

Γ ΨϕdΓ = ikA 0 Γ ψϕdΓ -ik Np j=1 Γ A j k 2 ε ∇ Γ ψ j • ∇ Γ ϕdΓ. (35)
In the above relations, ϕ designates a suitable test-function defined on Γ. In our case, we adapt easily this procedure to Equation ( 27) by rewriting the pseudodifferential equation as

Γ ΛgϕdΓ = Γ k 2 k 2 ε ∇ Γ u 0 • ∇ Γ ϕ -k 2 u 0 ϕdΓ. (36)
This formulation only requires the solution to an additional surface partial differential equation over Γ. All these equations can be solved efficiently [START_REF] Antoine | An improved surface radiation condition for high-frequency acoustic scattering problems[END_REF] by using a discretization based on a surface finite element method on Γ h . Let us introduce N Γ h as the number of degrees of freedom arising in the surface finite element method. Then, the resulting linear systems defining the OSRC through the Padé approximants can be solved at a linear cost according to N Γ h by a preconditioned ILUT Krylov iterative solver [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. All this process yields the efficient construction of the approximate initial condition u 0 on the boundary Γ. The theory of OSRCs assumes that Ω -is a convex domain [START_REF] Antoine | Advances in the On-Surface Radiation Condition Method: Theory, Numerics and Applications[END_REF]. Some coupling procedures exist to extend its application range to non-convex bodies [START_REF] Antoine | An algorithm coupling the OSRC and FEM for the computation of an approximate scattered acoustic field by a non-convex body[END_REF], which we won't consider here. Finally, an explicit analytical starter u is obtained if one uses the lowest-order approximation called Sommerfeld

approximation u = -i∂ n Γ u/k, on Γ. ( 37
)
This approximation will be used in the numerical tests presented in Section 5.

We will see that it yields interesting but also limited accuracy improvements for high wavenumbers and/or coarse meshes. Section 6 will reveal the interest of numerically using high-order models like the Padé model introduced above if high-order propagators are considered at the same time.

Once the surface field u 0 is known on Γ, the corresponding surface phase function φ can be computed. This is of interest in the sequel since we use a low-order propagation model giving an explicit phase function φ over the computational domain once its trace is known (through Equation ( 50)). If one uses the solution which consists in writing the phase at a surface point using the relation

φ = 1 ik log( u 0 A ), for u 0 = 0, (38)
the resulting calculation yields a discontinuous determination of the phase since the exponential function is a periodic function. A test example of such a problem is given on Figure 2 in the simple case of a circular cylinder at k = 10 for the Neumann boundary condition. This difficulty is known as the phase unwrapping problem. This is an important point to deal with in our approach because the involvement of a discontinuous phase in the weak variational formulation on A (see for example Equations (51)-( 52)) would clearly lead to difficulties in the finite element solution. To solve this problem, we use the following strategy. Consider the known field u 0 and differentiate the relation e ik φ = u 0 / A assuming that the surface field is not equal to zero. Taking the real part of the equation, one gets

∇ Γ φ = F := ℜ( 1 ik A u 0 ∇ Γ ( u 0 A )), (39)
where the surface vector field F is known. If one fixes the value of φ at an arbitrary surface point x * where u 0 is non-zero and one takes the surface divergence of the equation (39), then, a continuous determination of the phase φ can be computed as the unique solution to the well-posed surface partial differential equation

       -∆ Γ φ = -div Γ F , on Γ, φ(x * ) = arg u 0 (x * ) A(x * ) . ( 40
)
The numerical solution can be easily obtained by a surface finite element method, similarly to the OSRC technique, by the formulation

Γ ∇ Γ φ • ∇ Γ ϕdΓ = - Γ div Γ (F )ϕdΓ, (41)
for some test-functions ϕ and fixing the value of φ at a surface point x * . This process yields a continuous phase φ over the surface Γ, as shown on Figure 2 where we see the unwrapped version of the phase for our previous example, computed using the high-order OSRC solution. We also see the accuracy improvement related to the choice of the OSRC.

Construction of the evolution equation (propagator): general approach and low-order approximation

The next step of the approach consists in developing an approximate forward propagating equation called the propagator. Unlike the one-dimensional case, the situation is much more complex because of the boundary. We propose here a general construction of an approximate DtN map using pseudodifferential operators theory and elements developed in [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF]. We will use its lowest order approximation for the PR-FEM computations performed in Section 5 but will show in Section 6 that it is crucial to consider the proposed high-order models to investigate improved solutions at high-frequency and for coarse meshes.

Let us consider a point x in Ω b (or more generally in Ω + ). Since Ω -is convex, any observation point x can be projected onto Γ as a unique point x 0 = π Γ x ∈ Γ, where π Γ is the projection operator onto Γ. Let us introduce r = x 0 x and r = r . The outwardly directed unit normal vector to Γ at point x 0 is given by n Γ (x 0 ) = r/r (see Fig. 3). Since Γ is a compact submanifold of R d , a local coordinates system at a point x 0 of Γ can be chosen. Let us introduce the tangential variable s = (s 1 , s 2 ) and the radial variable r along n Γ (x 0 ), setting r = 0 at Γ = Γ 0 and Γ r for the parallel surface. We choose an orthogonal coordinates system on Γ. The covariant basis (τ 1 , τ 2 ) of the tangent plane T x 0 (Γ) compatible with the orientation of n(x 0 ) is better known as the principal basis, where τ 1 and τ 2 are the principal directions of curvatures to the surface. Setting R Γ as the curvature tensor of the tangent plane at a given point of the surface, the principal curvatures κ 1 and κ 2 of Γ satisfy: R Γ τ β = κ β τ β , for β = 1, 2, and the mean curvature H Γ = trace(R Γ )/2. Introducing the functions h β = 1 + rκ β , β = 1, 2, we get the expression of the Helmholtz operator in generalized coordinates

e ikα•x r R 0 R x x 0 Γ Γ r
L(r, s, ∂ r , ∂ s )u = ∂ 2 r + 2H r ∂ r + h -1 1 h -1 2 ∂ s • (h 2 h -1 1 ∂ s 1 , h 1 h -1 2 ∂ s 2 ) + k 2 , ( 42
) setting H r = (h -2 1 κ 1 + h -2 2 κ 2 )
. Now adapting the techniques used in [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF], we prove the existence and uniqueness of two pseudodifferential operators P ± such that the equation

(∂ r + iP + ) u = 0, on Γ r , (43)
characterizes the forward propagating part of the wavefield. The other operator P -yields the reentrant part of the field to Γ r . Let us introduce the total symbols p ± (r, s, ξ) = σ P ± of P ± = P ± (r, s, ∂ s ), setting ξ = (ξ 1 , ξ 2 ) as the Fourier covariable of s. These two symbols admit a unique asymptotic expansion in homogeneous symbols p ± -j j≥-1 of order -j as

p ± ∼ j≥-1 p ± -j , (44)
where the functions p ± -j satisfy: p ± -j (r, s, λξ) = λ -j p ± -j (r, s, ξ), ∀λ > 0 (see [START_REF] Taylor | Pseudodifferential Operators[END_REF] for the definition of ∼). The uniqueness of the expansion is fixed through the condition on the principal symbol

p ± 1 = ∓ k 2 -||ξ|| 2 , ( 45
) setting ||ξ|| 2 = h -2 1 ξ 2 1 + h -2 2 ξ 2 2
. Now if one considers for example the approximation of (43) by the equation

(∂ r + iOp(p + 1 )) u = 0, on Γ r , (46)
we obtain a nonlocal propagation equation to solve (Op(σ) designates the pseudodifferential operator with symbol σ). Using the symbol (45) leads to a correct computation of propagative modes but can exhibit a loss of accuracy for high-order spatial frequencies ξ such that ||ξ|| ≈ k since it is singular. A regularization (like in the OSRC technique) of p + 1 must be considered. In the same spirit as in [START_REF] Antoine | An improved surface radiation condition for high-frequency acoustic scattering problems[END_REF], we take:

p + 1,η = -k 1 - ||ξ|| 2 k 2 η , (47)
setting

k η = k + iη, with η = 0.4k 1/3 H 2/3 r (where H r = (h -2 1 κ 1 + h -2 2 κ 2 )
is the local mean curvature of the front Γ r ), and the associated propagator

(∂ r + iOp(p + 1,η )) u = 0, on Γ r . (48)
Directly using (47)-( 48) would be computationally expensive since it is a pseudodifferential nonlocal equation. In this paper, we restrict our numerical implementation to the simplest low-order model considering a zeroth-order approximation of the principal symbol corresponding to ξ = 0. A higher order localization process based for example on a second-order Taylor expansion for small values of ||ξ|| would yield a second-order partial differential paraxial wave equation valid in a cone of aperture about 15 degrees in the normal direction (more details can be found in the litterature about parabolic equations models-see e.g. [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF]). Considering hence the rough approximation:

(∂ r -ik) u = 0, on Γ r , (49)
we get the solution: u(r, s) = e ikr u(0, s), or in other words u(x) = u(x 0 )e ik||x-x 0 || . Since the phase φ(x 0 ) is known at x 0 by using the OSRC approximation, the phase at point x of Ω + is approximated by

φ(x) = φ(x 0 ) + ||x -x 0 ||. (50)
Even if we only use a very simple approximation of (43), our a priori complex approach is constructive and can be improved following two directions. The first one consists in using some real parabolic equations using for example the first symbol p + 1 and a Padé expansion to localize the resulting pseudodifferential operator (see Section 6). The second direction consists in incorporating more terms in the asymptotic expansion, e.g. computing p + 0 (and next using some approximations, e.g. for ξ = 0). All the corrective symbols can be obtained through some adaptations of the computations performed in [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF]. However, this also implies some difficulties for the Phase Reduction approximation since special structured meshes must be used following the normal directions to solve the resulting equations. Another interesting point of our approach is that the previous construction can be adapted to variable coefficients problems like inhomogeneous Helmholtz equations. All these more advanced points are beyond the scope of the present paper and are currently being investigated.

Computing the true solution u

Since φ has now been computed, we follow the same path as in the onedimensional case. Replacing u by Ae ik φ in ( 18) and taking some test-functions v = Be ik φ , we obtain the alternative formulation: find

A ∈ H 1 (Ω b ) such that A(A, B) = L(B), (51)
for all B ∈ H 1 (Ω b ). The sesquilinear form is given by

A(A, B) = (∇A, ∇B) 0,Ω b + ik (A∇ φ, ∇B) 0,Ω b -(∇A, B∇ φ) 0,Ω b -k 2 ((1 -||∇ φ|| 2 )A, B)) 0,Ω b (A∇ Σ A, ∇ Σ B) 0,Σ + ik ((A∇ Σ φ)A, ∇ Σ B) 0,Σ -(A∇ Σ A, B∇ Σ φ) 0,Σ + k 2 ((A∇ Σ φ, ∇ Σ φ)A, B)) 0,Σ + (βA, B) 0,Σ (52)
and the linear form L by

L(B) = -(f, B) 0,Γ (53)
with f = ge -ik φ and φ given. Even if the formulation seems complicated at first sight, all the quantities can be easily computed by using assembling procedures available in most basic finite element codes.

Numerical performance of the lowest order model

To show the improvements that can be expected from our approach using the lowest order model based on Equations (37) and (50), we begin by detailing in Section 5.1 a modal study for the sound-hard disk. Full plane wave problems are considered for the Neumann problem in Section 5.2 and for the Dirichlet problem in Section 5.3. All results were obtained using linear (P 1 ) triangular finite elements and uniform meshes. The discretization density is defined as n λ = λ/h, where h measures the size of the triangles. The linear systems were solved with a sparse direct solver.

For a given approximated solution f calc (x), x ∈ Ω b , with exact counterpart f exact (x), we define the relative L 2 (Ω b )-error as:

ε 2 (f ) := Ω b |f calc (x) -f exact (x)| 2 dx 1/2 Ω b |f exact (x)| 2 dx 1/2 (54)
In all the following examples, f exact denotes the exact solution of the problem with the BGT2 ABC, so that the error is not influenced by the non-perfectly transparent boundary condition.

Scattering by a sound-hard circular cylinder: mode-by-mode analysis

A first example concerns the scattering problem by a circular cylinder Ω -= D 0 of radius R 0 centered at the origin. Its boundary Γ is therefore given by the circle Γ = C 0 . We also consider that the fictitious boundary Σ = C 1 is a larger circle of radius R 1 > R 0 , again centered at the origin. As a consequence, the bounded computational domain Ω b is the crown delimited by the boundaries C 0 and C 1 .

If we now consider an incident wavefield fixed by a mode of order m, that is

u inc m (x) = J m (kr)e imϕ , m ∈ Z, (55)
in the polar coordinates system (r, ϕ), then, the exact exterior modal solution to the scattering problem by the sound-hard circular cylinder D 0 is given by

u m (x) = - J ′ m (kR 0 ) H ′(1) m (kR 0 ) H (1) m (kr)e imϕ , r ≥ R 0 , m ∈ Z. (56)
In the case where we truncate the exterior problem using the BGT2-like radiation condition ( 13)-( 15), the analytical exact solution reads

u exact m (x) = (a m H (1) m (kr) + b m H (2) m (kr))e imϕ , r ≥ R 0 , m ∈ Z, (57)
since spurious reflection occurs at C 1 . The coefficients a m and b m are given by solving a 2 × 2 linear system of equations imposing the Neumann boundary condition at C 0 and the ABC ( 13)-( 15) at C 1 . This gives the following expressions

a m = - A 22 J ′ m (kR 0 ) D , b m = A 21 J ′ m (kR 0 ) D , ( 58
) setting D = A 11 A 22 -A 21 A 12 ,                        A 11 = H ′(1) m (kR 0 ), A 21 = kH ′(1) m (kR 1 ) -B m H (1) m (kR 1 ), A 12 = H ′(2) m (kR 0 ), A 22 = kH ′(2) m (kR 1 ) -B m H (2) m (kR 1 ), B m = -(α m m 2 R 2 1 + β m ), α m = - 1 2ik (1 + i kR 1 ) -1 , β m = -ik + 1 2R 1 + 1 8iR 1 (1 + kR 1 )
.

(59)

Hence, being able to compute both the exact solutions to the exterior and truncated problems allow us to separate in a numerical study i) the analytical error coming from the domain truncation using the ABC ( 13)-( 15) and ii) the error linked to the numerical approximation of the truncated problem e.g. by a finite element method. Therefore, the convergence of the numerical approximation is only visible from the exact solution u exact m given by ( 57)-(59) which is now considered as our reference solution for R 0 = 1 and R 1 = 2. Concerning the PR-FEM, we use the exact phase on Γ, i.e., φ(θ) = cmθ. The reason is that the unwrapping technique cannot be applied since φ is not a periodic function on [0, 2π] (unlike for the plane wave problem).

We present on Figure 4 the real parts of the wavefield u h and A h computed respectively by the FEM and the PR-FEM. We also report the absolute error compared to the exact analytical fields u = u exact m and A. The wavenumber is fixed to k = 25 and the density to n λ = 20. From these computations, we observe that the field A h oscillates much less than u h . As expected, the accuracy of the FEM for computing A h is higher than for u h since less pollution arises. The mode m is propagative because m ≤ k. This means that the rough approximation (50) is certainly not so bad. On Figure 5 we draw the evolution of the error ε 2 for k = 50 according to the mesh refinement for the three modes m = 0, m = 25 (both purely propagative) and m = 50 (which is at the transition between the region of propagative and evanescent modes). We see that the phase approximation gives an interesting improvement for the propagative part but the accuracy improvement is affected for higher harmonics where an approximate phase is much more complicated to obtain. The importance of this accuracy degradation is however moderated in practice by the fact that the solution of general scattering problems (e.g. under plane wave incidence) is a superposition of harmonics. We will see in Section 6 that higher order propagative models can capture correctly the oscillations linked to large modes even at high frequencies. We also see on Figure 5 that PR-FEM converges with the meshsize as expected for linear FEM. Second-order FEM have also been tested and lead to the expected convergence. Finally, Figure 6 shows the error behavior of the different methods for the modes m = 0 and m = k with respect to the wavenumber. We get a higher accuracy for the PR-FEM which is again penalized for higher order harmonics m.

ℜ(u h -u) 0.00397 0.0892 ℜ(A h ) -0.000149 -1.2e -05 ℜ(A h -A)

Scattering by a sound-hard circular cylinder: plane-wave incidence

We consider an incident plane wave e ikx . From now on, the PR-FEM uses the lowest-order (Sommerfeld) OSRC approximation (37). In this case, the phase φ of u is actually the same as the phase of u inc and is thus the same as the phase of the Dirichlet boundary condition. However, for completeness, the unwrapping phase technique is used to get the continuous extension of φ. The reason why we do not consider higher order OSRCs here is that the accuracy improvement of the starter would not be visible if high order exterior propagators are not used at the same time (this point will be analyzed in Section 6 on a simple analytical problem). Another point related to using higher order OSRCs is that we would have to refine the surface mesh just for doing the higher order OSRC (since it requires the solution of surface PDEs). It should be possible to relax this point by using two different grids-a refined (d -1)-dimensional grid for the OSRC computation, and a coarser d-dimensional grid for solving (51)-(53). Interpolation between the two grids could for example be performed using an L 2 -projection [START_REF] Geuzaine | A Galerkin projection method for mixed finite elements[END_REF].

We begin by considering the same two concentric circular boundaries of Section 5.1. For a plane wave, the exact reference solution u (and so A) is given as the linear superposition of the elementary reference solutions (using BGT2) of the previous subsection for each harmonics. We plot on Figures 7 and 8 repectively the fields u h and A h and the approximate phase φ. We see that the finite element error is less important for the proposed formulation than for the classical formulation. Figure 9 presents the error ε 2 with respect to the wavenumber k for two different discretization densities. This again shows that the error is maintained for a fixed mesh as the wavenumber increases. The rate of convergence of the error for three different wavenumbers is reported on Figure 10. An admissible error can be obtained for low discretization densities or equivalently a higher accuracy can be expected for a fixed density. Therefore, the pollution error is reduced significantly on this example. Finally, we draw on Figure 11 the RCS for k = 25 and n λ = 3. We observe a much better prediction of the far-field for these parameters when the PR-FEM is used. This is in particular clearly visible in the transition and shadow zones.

Scattering by sound-soft obstacles

Let us now consider the Dirichlet problem. In this case, the starter is given immediately from the trace of the wavefield over Γ.

The first test-case is again related to the previous two concentric boundaries. We present on the top of Figure 13 the computed approximate phases using the low-order approximation (50) at k = 25. Compared to the exact solution u h reported on Figure 12, we see that the complex amplitude A h oscillates less than u h . We report on Figure 14 the evolution of the error ε 2 with respect to k for two densities of discretization points. We can see that, compared to the classical approach, the new PR-FEM formulation even with a low-order model leads to significant accuracy improvements. Moreover, for both discretization densities n λ , the error appears to be almost constant, which means that the pollution error has been virtually eliminated. Figure 15 presents the behavior of the error according to n λ for three wavenumbers k = 10, k = 25 and k = 50. Again, the benefits of the new formulation are clearly visible: with 6 points per wavelength the relative error with the new formulation is comprised between 2.7% and 5% whereas the error with the original formulation varies between 27% and 100%.

We can expect similar gains in three dimensions. For example, the scattering problem by the unit sphere for a wavenumber k = 10 is considered on Figure 16, where the traces of the fields are plotted in the planes (Ox 1 x 2 ) and (Ox 1 x 3 ). We can see that, indeed, the field based on the slowly varying envelope A h oscillates less than with the physical scattered field u h .

To get a qualitative idea on how the PR-FEM method behaves on non-smooth geometries, we show on Figure 17 the results obtained for a thin plate under oblique plane wave incidence. The amplitude again clearly oscillates less than the original unknown, even if the reduction in the number of oscillations is less than in the case of smooth scatterers. This might be due to the low order nature of the propagator. However, as shown on Figure 18, the improvement is still substantial. For n λ = 4, the maximum error is reduced by one order of magnitude with PR-FEM. This is also true for the Neumann problem, where the low order approximation used for the starter also plays a role.

More elaborate techniques can be envisioned to approximate the phase φ(x) attains its limitation while a very satisfying error can be expected for a better approximation of the phase using e.g. the eikonal solution. The reason of the limitation of the low-order model is that high order harmonics must be considered to compute the solution, and we showed in the mode-by-mode study developed in Section 5.1 that, indeed, high order harmonics include larger errors for the low-order model, resulting in limited accuracy of the solution for coarse meshes. As we will see in Section 6, high-order models show great potential for computing an improved phase at very high frequencies for a low computational cost. Figure 20 shows how the error evolves for a fixed mesh and increasing k. Again, a more accurate phase computation shows that an improved accuracy can be expected. The advantages of the new formulation are clearly visible: for a prescribed tolerance, the PR-FEM allows to solve the problem on much coarser grids than the original FEM, the error control depending on the phase accuracy.

Toward higher order models

We have seen in the previous section that an interesting accuracy improvement is obtained via the PR-FEM with a low-order model even for relatively high wavenumbers. At the same time, we have also shown that an increased accuracy can be expected if a better approximation of the phase is at hand, using for example the solution of the eikonal equation.

We show in this section that higher order propagative models in conjunction with an improved OSRCs can provide the phase with high accuracy, both on the surface of the scatterer and in the computational domain, even for extremely large frequencies. An example of such high-order models is proposed in [START_REF] Antoine | Computing high-frequency scattered fields by beam propagation methods: A prospective study[END_REF] using BPM techniques. Here, we consider the microlocal approximation given in Section 4 and test it in the case of the circular cylinder D 0 with radius R 0 = 1 (see Section 5.1) where analytical formulae can be obtained using the mode-by-mode superposition of the solution to each problem. All the results are presented for k = 100 (see Figure 21); tests for a similar configuration but with k = 2000 for example lead to the same conclusions.

Let us introduce the following notations. For a given boundary condition (related to an OSRC operator Λ) and an exterior model (Helmholtz equation or propagator), we denote the corresponding solution by u BC,Model , where BC reflects the type of boundary condition and Model the equation solved in the exterior domain:

• BC = ik refers to the low-order OSRC (37) and BC = √ refers to the square-root OSRC (28); • Model = H refers to the Helmholtz equation, Model = LOP refers to the loworder propagator (49) and Model = HOP refers to high-order propagator (48).

Using the same notations, the slowly varying amplitudes are defined by

A BC,Model = u exact u BC,Model . (61)
First, we test the effect of using an accurate starter based on high-order OS-RCs. We report on Figure 22 the normalized real parts of the two solutions A ik,H (left) and A √ ,H (right). We can observe an interesting improvement if we use the OSRC based on the square-root operator compared to the Sommerferld condition. This is particularly clear near the boundary where the modeling of surface rays is crucial. The attenuated oscillations related to us- ing a high-order OSRC implies that a coarser mesh should be possible close to the boundary.

Second, we test the influence of the exterior propagation equation. We have seen in Section 4.1.2 that low and high-order approximate propagators can be considered. We report on Figure 23 the real parts of A √ ,LOP and A √ ,HOP .

We see that using a low-order model leads to a less oscillating solution than the initial field. However, visible oscillations remain, explaining the limitations met in Section 5 when higher frequencies and/or coarser grids were considered. Using a high-order propagator solves the problem and yields a very close solution to using the Helmholtz equation (see right part of Figure 22). We can observe that even the field behavior near the boundary is satisfactory. We do not present here the numerical results obtained by using the localized Padé version of the square-root operator, which are similar to those obtained by using the square-root propagator.

Investigating PR-FEM solutions based on these two high-order generalizations shows therefore some great potential in reducing pollution problems. This is beyond the scope of the present paper and will be the subject of a second 21). The behavior of the near-field is close to the one expected using the Helmholtz equation (see right pictures of Figure 22). The high-order OSRC is used for building the starter.

part.

Conclusion

We have presented a procedure, dubbed PR-FEM, to gain significant accuracy in the finite element solution of time-harmonic scattering problems at high wavenumbers. It is based on an a priori approximate determination of the phase of the scattered wave using a starter/propagator technique together with a variational formulation in terms of the resulting slowly varying envelope. Numerical tests on simple two-and three-dimensional convex scatterers using the lowest-order starter/propagator model show that the technique reduces the pollution error and allows use of much coarser grids than the standard FEM. Contrary to competing approaches, the proposed method does not require any new finite element basis functions and can thus be easily implemented in existing finite element codes. Finally, accurate phase computations based on high-order starters and propagators show that improved solutions for large frequencies can be expected. This last aspect will be analyzed in Part II of this paper.

xFig. 1 .

 1 Fig.1. Real part of the scattered field at a frequency k = 40 for a density of discretization points n λ = 8. The pollution effect can be observed on the classical FEM approach but is avoided by the PR-FEM formulation.

Fig. 2 .

 2 Fig. 2. Surface phase φ for the Neumann problem and the unit circular cylinder at k = 10: exact solution, Sommerfeld solution, high-order OSRC and unwrapped high-order OSRC solutions.

Fig. 3 .

 3 Fig. 3. Illustration of the notations used in the construction of the forward propagating model.

Fig. 4 .

 4 Fig. 4. Neumann problem: Real part of u h and A h and the associated errors for the mode m = 10. The wave number is k = 25 and the density is n λ = 20.

n λ ε 2 FEM 2 Fig. 5 .

 225 Fig. 5. Neumann problem: Evolution of the error ε 2 for several modes at k = 50 according to the density of discretization n λ .

k ε 2 FEMFig. 6 .

 26 Fig. 6. Neumann problem: Evolution of the error ε 2 for the modes 0 and k with respect to the wave number k, setting n λ = 10.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Neumann problem: Real part of u h and the associated errors. The wave number is k = 25 and the density is n λ = 20.

k ε 2 FEM (n λ = 10 )Fig. 9 .

 2109 Fig. 9. Neumann problem: Evolution of the error ε 2 with respect to the wave number k for two discretizations.

n λ ε 2 FEM 1 Fig. 10 .

 2110 Fig. 10. Neumann problem: Evolution of the error ε 2 with respect to the density of discretization points n λ for three wave numbers k.

Fig. 11 .

 11 Fig. 11. Neumann problem: Bistatic Radar Cross Section for k = 25 and n λ = 3.

Fig. 12 .

 12 Fig. 12. Dirichlet problem: real part of the original solution u h for k = 25.

Fig. 13 .

 13 Fig. 13. Dirichlet problem: two different approximations of the phase of the complexvalued scattering solution for k = 25. Left: schematic depiction of the trajectories for the simple evolution equation (top) and for the eikonal equation (bottom). Middle and right: approximate phase φ and real part of the resulting amplitude A h for the simple evolution equation (top) and for the eikonal equation (bottom).

k ε 2 FEM (n λ = 10 )Fig. 14 .

 21014 Fig. 14. Dirichlet problem: Evolution of the error ε 2 with respect to the wavenumber k for two densities of discretization points per wavelength n λ (n λ = 10 and n λ = 20).

n λ ε 2 FEM 1 Fig. 15 .

 2115 Fig. 15. Dirichlet problem: Evolution of the error ε 2 according to the density n λ for three wavenumbers (k = 10, k = 25 and k = 50).

Fig. 18 .

 18 Fig. 18. Dirichlet problem: Plate under 45 • plane wave incidence for k = 25 and n λ = 4. Left: Absolute error between the reference solution u (computed with n λ = 20) and the direct solution u h . Right: Absolute error between the reference solution u and the PR-FEM solution.

k ε 2 FEM (n λ = 2 ) 1 Fig. 19 .

 22119 Fig.[START_REF] Giladi | A hybrid numerical asymptotic method for scattering problems[END_REF]. Dirichlet problem: Evolution of the error ε 2 with respect to the wavenumber k for two densities of discretization points n λ = 2 and n λ = 4 using the classical FEM, the PR-FEM with the low-order propagator and the eikonal equation.

1 Fig. 20 .

 120 Fig.[START_REF] Harari | Analytical and numerical studies of a finite element PML for the Helmholtz equation[END_REF]. Dirichlet problem: Evolution of the error ε 2 with respect to the wavenumber k (k = 1 to k = 256) using a fixed uniform mesh of 7592 triangles. The results are presented for the classical FEM, the PR-FEM with the low-order propagator and the eikonal equation.

29

 29

-1 1 ℜFig. 21 .

 121 Fig. 21. Neumann problem: normalized exact solution at k = 100.

Fig. 22 .

 22 Fig. 22. Choice of the starter: Representation of ℜ( A ik,H |A ik,H | (x)) (left) and ℜ( A √ ,H |A √ ,H | (x)) (right) for the Neumann problem. We can observe the effect of using the low-order OSRC approximation BC = ik or the high-order model BC = √ . Zooming shows the boundary layer related to the incorporation of creeping waves in the OSRC model. This can affect the PR-FEM since a low-order OSRC would require a local remeshing to the surface for numerically capturing these oscillations.

Fig. 23 .

 23 Fig. 23. Choice of the propagator: Effect of using the low-order (49) (left) or the high-order (48) propagator (right) for the Neumann problem. Zooming shows the remaining oscillations close to the boundary for the low-order model. The solution is however much less oscillating than the true solution (see Figure21). The behavior of the near-field is close to the one expected using the Helmholtz equation (see right pictures of Figure22). The high-order OSRC is used for building the starter.

Acknowledgement

The authors wish to thank Jacob Bedrossian for the implementation of the Eikonal equation solver. This work was supported in part by NSF under grant DMS-0609824.

in the computational domain Ω. At the expense of costlier numerics, a prime candidate is to compute the phase by solving the eikonal equation

i.e., by using the asymptotic geometrical optics solution to the Helmholtz equation [START_REF] Keller | Asymptotic methods for partial differential equations: the reduced wave equation and Maxwell's equations[END_REF]. Solving Equation (60), we get a new approximate phase φ and slowly varying envelope A h . These two fields are presented on the bottom of Figure 13. We see that the phase is close to the one expected from u h (see Figure 12). The resulting field A h therefore oscillates less than using the low-order solution. To see the impact of improving the phase computation, we report on Figure 19 the evolution of ε 2 according to the wavenumber for two coarse grids. We see that for large wave numbers, the low-order model


Xavier Antoine
email: xavier.antoine@iecn.u-nancy.fr

Yuexia Huang

Ya Yan

Computing High-Frequency Scattered Fields by Beam Propagation Methods: A Prospective Study

Keywords: acoustic scattering, beam propagation methods, one-way equations

This paper presents some theoretical and numerical investigations concerning the fast computation of an exterior wave field to a scatterer by the Beam Propagation Method (BPM). Different models are presented and compared. It appears that the approach is able to correctly model the propagation of the propagative modes of the wave field while inaccuracies still remain for the evanescent and transition modes.

Introduction

The aim of this paper is to prospect the possible application of the Beam Propagation Method (BPM) for solving high-frequency scattering problems.

To focus our study on special features of the BPM, we rather restrict our developments to the sound-soft scattering problems but extensions could also be considered to more general problems like sound-hard or impedance problems.

Hence, the wave field u satisfies the Helmholtz equation in the unbounded domain exterior to the cylinder and can be decomposed as u = u (i) + u (s) , where u (i) is the given incident wave, u (s) is the unknown scattered wave solution to the boundary-value problem

         ∆u (s) + k 2 u (s) = 0, in Ω + , u (s) = -u (i) , on ∂Ω, lim r→∞ √ r ∂u (s) ∂r -iku (s) = 0,
where Ω is the cross-section of the cylinder, Ω + is the domain exterior to Ω, ∂Ω is the boundary of Ω, r is the radial variable in the polar coordinate system and k = 2π/λ is the wavenumber. The radiation condition at infinity ensures the uniqueness of the solution. A brief illustration of the problem is given on Various numerical methods have been developed over the past decades to solve the time-harmonic acoustic problem. One well-known possible approach is to truncate the unbounded domain Ω + with an Artificial Boundary Condition (ABC) and to use the Finite Element Method (FEM) in a bounded domain [START_REF] Thompson | A review of finite element methods for time-harmonic acoustics[END_REF]. This method can handle scatterers with complicated boundaries.

Another important method is the Boundary Integral Equation (BIE) method [START_REF] Burton | The application of integral equation methods to the numerical solution of some exterior boundary-value problems[END_REF][START_REF] Chen | Boundary Element Methods[END_REF][START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF]. The problem is then reformulated as an integral equation defined on ∂Ω and solved through techniques combining an iterative linear algebra solver [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] and compression algorithms like e.g. the Fast Multilevel Multipole (FMM) technique [START_REF] Rokhlin | Rapid solution of integral equations of scattering theory in two dimensions[END_REF]. In many applications, the angular frequency ω is extremely high, the wavenumber k = ω/c is large and the wavelength λ = 2π/k is much smaller than the characteristic length of the scatterer. Then, both the ABC-FEM and the BIE have difficulties in solving such large scale problems since, in particular, the associated linear systems have a very large size.

Moreover, these two methods require that an integral must be evaluated at any exterior point where the solution is desired. Both these two methods require more than O(N ) operations for each point where N is the number of grid points used to discretize the boundary. When N is large, this is quite expensive.

The Beam Propagation Method (BPM) is an efficient approximation method widely used in optical waveguide modelling. Under the assumption that the wave field is dominated by its forward component, the BPM gives rise to one-way equations that approximate the Helmholtz equation. It is useful for waveguides that change slowly in the propagation direction. The one-way equations have only a first order derivative in the propagation direction and can be efficiently solved as an initial value problem. Operator rational approximations are involved in solving the one-way models. There are various types of rational approximations to the one-way equations. For example, the BPMs can be solved by rational approximations to the square root operator [START_REF] Hadley | Wide-angle beam propagation using Padé approximant operators[END_REF][START_REF] Lu | A complex coefficient rational Padé approximaton of √ 1 + x[END_REF][START_REF] Milinazzo | Rational square-root approximations for parabolic equation algorithms[END_REF] or rational approximations to the exponential of the square root operator, i.e.

the propagator [START_REF] Lu | Beam propagation method using a [(p-1)/p] Padé approximant of the propagator[END_REF][START_REF] Yevick | The application of complex Padé approximants to vector field propagation[END_REF].

In this paper, an attempt is made to apply the BPM method to the scattering problem. We develop an approximate method which can efficiently solve the scattered wave field. Our method requires only O(1) operations for each point where the solution is desired.

The plan of the paper is the following. In Section 2, we use a curvilinear coordinate system to rewite the scattering problem following the parallel surfaces to the scatterer. We then explain in details in the third Section the approach by Beam Propagation Methods, developing different models and approximations as well as rational approximations for more efficiency. Section 4

proposes some numerical simulations and investigations to see how the method applies in terms of efficiency and accuracy. Finally, the last Section gives some conclusions.

Curvilinear coordinate system

The scattering problem is different from the propagation problem in waveguides and in particular, the geometry itself must be treated with special attention. In the following, we first develop a coordinate transformation to apply the BPM method in the new coordinate system.

We assume that the boundary of the scatterer ∂Ω is given by

x = f (s), y = g(s), 0 ≤ s ≤ l,
where s is the arclength and l is the total length of ∂Ω. Using the counterclockwise direction as the direction of increasing s, we have an outward unit normal vector given by

n(s) = (g ′ (s), -f ′ (s)) . (1)
As illustrated in Figure 2, any point (x, y) in the exterior domain Ω + (assuming that Ω is convex) can be written as

x = f (s) + ξg ′ (s) y = g(s) -ξf ′ (s), (2)
for 0 ≤ s ≤ l and 0 ≤ ξ ≤ +∞. Under the new coordinate system (ξ, s), the Helmholtz equation becomes

u ξξ + γ ξ γ u ξ + 1 γ ( 1 γ u s ) s + k 2 u = 0,
where γ = 1+ξκ(s) and κ is the curvature of the boundary ∂Ω. For u = γ -1/2 v, we have

v ξξ + 1 γ 2 v s s + k 2 + γ 2 ξ 4γ 2 - γ ss 2γ 3 + γ 2 s 4γ 4 v = 0. ( 3
)
For convenience, we define an operator L by

L(ξ) = ∂ ∂s 1 γ 2 ∂ ∂s + k 2 + γ 2 ξ 4γ 2 - γ ss 2γ 3 + γ 2 s 4γ 4 . (4)
It acts on functions of s and it depends on ξ as a parameter.

The BPM approach

In the new coordinate system, the scattered wave propagates outwards in the increasing ξ direction. This suggests that the scattered field may be approximated by one-way models that are first order in ξ. Following the BPM for optical waveguides, we first approximate equation ( 3) by one-way equations and then apply operator rational approximations.

There are a few different one-way models used in the BPM. One possibility is to approximate equation (3) by v ξ = i L(ξ)v using the square root operator L(ξ). Here, we adopt the Energy-Conserving model [START_REF] Collins | An energy-conserving parabolic equation for elastic media[END_REF][START_REF] Godin | Reciprocity and energy conservation within the parabolic approximation[END_REF][START_REF] Vassallo | Difficulty with vectorial BPM[END_REF] which gives improved accuracy for slowly varying waveguide. The energy-conserving model can be derived from the continuity of the power flux and it involves a transform using the fourth root of the operator L. We have

φ ξ = i √ Lφ, φ = 4 √ Lv, φ(0) = 4 √ Lv(0). (5)
In terms of the original function v, the energy-conserving one-way model is

v ξ = i √ L -L -1/4 d(L 1/4 ) dξ v.
Notice that the original boundary value problem of the Helmholtz equation for u is now approximated by an initial value problem for φ. As in the standard BPM, the square root (and the fourth root) of L must be approximated. For this purpose, we introduce the operator X by

L = k 2 0 (1 + X),
where k 0 is a reference wavenumber. Then, equation [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF] becomes

φ ξ = ik 0 √ 1 + Xφ φ(0) = -4 L(0)u (i) .
Formally, for the step from ξ j to ξ j+1 = ξ j + h, we have

φ j+1 = P φ j , P = e ik 0 h √ 1+X φ j , ( 6
)
where P is the one-way propagator and X is evaluated at the midpoint ξ =

ξ j + h/2.
The propagator P can be approximated by rational functions of X. For waveguide problem, the standard [p/p] Padé approximation provides a good approximation for the propagating modes but fails to suppress the evanescent modes, while the [(p-1)/p] Padé approximant can damp the evanescent modes but gives less accurate results for the propagating modes. For the current scattering problem, a proper treatment of the evanescent modes appears necessary.

The propagator-θ method [START_REF] Chui | A propagator-θ beam propagation method[END_REF] combines the [p/p] and [(p -1)/p] Padé approximants by a parameter θ. With a suitable choice of p and θ, the propagator-θ method gives a better balance for approximating both the propagating and evanescent modes.

The propagator-θ approximant of degree p takes the following form

P ≈ R p (θ) = 1 + α 1 X + . . . + α p-1 X p-1 + α p X p 1 + β 1 X + . . . + β p-1 X p-1 + β p X p (7)
= p k=1 1 + c k X 1 + b k X = 1 + p k=1 a k X 1 + b k X (8) = d 0 + p k=1 d k 1 + b k X . ( 9
)
It is related to the following [p/p] and [(p -1)/p] Padé approximants

R p-1,p = 1 + α (0) 1 X + . . . + α (0) p-1 X p-1 1 + β (0) 1 X + . . . + β (0) p-1 X p-1 + β (0) p X p (10) R p,p = 1 + α (1) 1 X + . . . + α (1) p-1 X p-1 + α (1) p X p 1 + β (1) 1 X + . . . + β (1) p-1 X p-1 + β (1) p X p . ( 11
)
The coefficients {α

(0) i }, {β (0)
i }, {α

i } and {β

i } can be solved from a linear system of equations assuming that the Taylor coefficients of the function e ihk 0 √ 1+X is available. The coefficients {α i } and {β i }, i = 1, 2, ..., p, can be calculated from

α i = θα (1) i + (1 -θ)α (0) i , β i = θβ (1) i + (1 -θ)β (0) i ,
where α (0) p = 0 and θ ∈ [0, 1]. The expression ( 9) is more convenient for numerical implementation than the standard form [START_REF] Collins | An energy-conserving parabolic equation for elastic media[END_REF]. The coefficients c k and b k are first obtained by factorizations of the dominator and numerator of [START_REF] Collins | An energy-conserving parabolic equation for elastic media[END_REF].

Coefficients a k and d k can be calculated from

a k = (c k -b k ) j =k c j -b k b j -b k , (12)
d k = - a k b k , d 0 = 1 - p k=1 d k . ( 13
)
With the rational approximation ( 9), the propagation step ( 6) is approximated by

φ j+1 = P φ j ≈ e ihk 0 d 0 + p k=1 d k 1 + b k X φ j = e ihk 0 d 0 φ j + p k=1 d k w k , (14)
where w k is solved from

(1 + b k X)w k = φ j , k = 1, 2, ..., p. (15)
Since X is related to L as in (4), the discretized form of ( 15) is a periodic tridiagonal system. If N is the number of grid points for discretizing the s variable, the required number of operations in each step is O(N ). Therefore, the number of operations required for each point is only of O(1).

The 4-th root operator 4 √ L = √ k 0 4 √ 1 + X must be evaluated when the unknown function v is transformed to φ, or vice versa. If the scattered wave field is required in the entire computational domain, this operator should be evaluated after each marching step of φ. In order to keep O(1) operations, a rational approximation for 4 √ 1 + X is used.

For a constant ν, the function q(x) = (1 + x) ν has a continued fraction expansion as

q(x) = 1 + σ 1 x 1 + σ 2 x 1 + σ 3 x 1 + . . . ( 16
)
where σ 1 = ν and

σ 2k = k -ν 2(2k -1) , σ 2k+1 = k + ν 2(2k + 1) , (17)
for k = 1, 2, • • • . Let q m (x) be the truncation of ( 16) up to the term σ m x, we have

q 2m (x) = R m,m (x), ( 18
)
where R m,m is the [m, m] Padé approximant of q(x). Furthermore, the sequence {q 2m (x)} converges to q(x) for all complex x, except for -∞ < x < -1 [START_REF] Lu | A complex coefficient rational Padé approximaton of √ 1 + x[END_REF].

For practical computation we rearrange q 2m (x) as

q 2m (x) = 1 + σ 1 x 1 + • • • 1 + σ 2m x 1 + iτ = m k=1 1 + ck x 1 + bk x = 1 + m k=1 ãk x 1 + bk x (19) = d0 + m k=1 dk 1 + bk x . (20)
The coefficients { bk } and {c k } can be calculated from a tridiagonal matrix related to the coefficients {σ k }, and we can use the same formulas ( 12) and ( 13) to calculate coefficients {ã k } and { dk }.

For our case, we set ν = 1/4. The initial value of φ can be approximately evaluated by

φ 0 = k 0 (1 + X) 1/4 v 0 ≈ k 0 ( d0 + m k=1 dk 1 + bk X )v 0 = k 0 ( d0 v 0 + m k=1 dk η k ), (21)
where the operator X is evaluated at ξ = 0 and w k is solved from

(1 + bk X)η k = v 0 , k = 1, 2, ..., m. (22)
The transformation from φ to v, i.e.

v j+1 = k -1/2 0 (1 + X) -1/4 φ j+1 (23)
can be implemented similarly. As in ( 14) and ( 15), the total number of operations needed in these approximations is also of O(N ). As far as the scattered wave field is concerned, the number of operations spent on each point is only of O(1).

Numerical results

We consider the scattering problem of a circular cylinder centered at the origin with radius a as shown in Figure 3. The boundary of the circular cross-section is given by x = a cos θ and y = a sin θ. The outward unit normal vector is n(s) = (cos θ, sin θ). In terms of ξ and s, the coordinate of a point (x, y) in Ω + are x = (a + ξ) cos s a , y = (a + ξ) sin s a .

Note that r = a + ξ for this case. The operator L(ξ) becomes

L(ξ) = a 2 (a + ξ) 2 ∂ 2 ∂s 2 + k 2 + 1 4(a + ξ) 2 . ( 24
)
The incident plane wave u (i) is given by

u (i) = -exp (ik(x cos ψ + y sin ψ)), ( 25
)
where ψ is the incident angle. In the following, we assume that a = 1, ψ = 0 In particular, the BPM solution has some undesirable oscillations.

For another example, we consider the scattering of an elliptical cylinder.

The incident plane wave has an incident angle ψ = 35 o (as shown in Figure 5).

The cross-section of the elliptical cylinder has semi-axis a = 1 and b = 0.25. approximations when the scattered wave is strong, but they produce incorrect oscillations when the scattered wave is weak. To understand the limitation of our method, we expand the plane incident wave u (i) in a series of Bessel functions as follows: where

u (i) = -exp(ikr cos(θ -ψ)) = - +∞ m=-∞ i m J m (kr) exp(im(θ -ψ)) = +∞ m=0 u (i) m ,
u (i) m = i m c m J m (kr) cos(m(θ -ψ))
with c 0 = 1 and c m = 2 for m = 1, 2, .... For this scattering problem, each term u (i) m gives rise to its scattered wave

u (s) m = c m J m (ka) H (1) m (ka) H (1) m (kr) cos(m(θ -ψ)),
where H

m is the Hankel function of the first kind. Therefore, the exact scattered wave for the plane incident wave can be written in the following Mie series:

u (s) = +∞ m=0 u (s) m . (26)
The above infinite sum can be truncated as follows:

u (s) ≈ u (s),M = M m=0 u (s) m .
In Figure 7, we plot the absolute values of u (s)

m for θ = 0, m = 0, ..., 60 and at r = 2 and r = 9. We observe that the first 40 (a quantity slightly larger than ka) terms are significant. Therefore, M should be at least 40 if u (s),M is to be a good approximation of u (s) . This suggests that to obtain a BPM solution with acceptable accuracy at least the same number of terms should be properly modelled. For this example, the E-C model can be solved analytically. In fact, the exact solution u m = u (s) m is separable and

L(ξ)u m = -m 2 a 2 r 2 + k 2 + 1 4r 2 u m = 1 -4m 2 4r 2 + k 2 u m .
That is, the operator L(ξ) becomes a scalar multiplicator on u m . Thus, its square root can be evaluated directly (without using rational approximations).

Following the decomposition of the exact solution (26), we can decompose φ

(of the E-C model) as φ = ∞ m=0 φ m .
Then, the propagation step from ξ j to ξ j+1 is reduced to:

φ m,j+1 ≈ exp ih L(ξ j+1/2 ) φ m,j = exp ih 1 -4m 2 4(a + ξ j+1/2 ) 2 + k 2 φ m,j .
The transformation between φ and v can be implemented in a similar manner.

Conclusion

In this paper, we develop and implement the Beam Propagation Method for 2-D scattering problems associated with acoustic sound-soft cylinders . This is an approximate method for scattering problems in the high frequency regime.

The method is very efficient, since the required number of operations is O(1) for each point where a solution is calculated. However, the accuracy of this method is limited. The BPM approach cannot accurately compute the scattered wave from a circular cylinder. Nevertheless, the methods may be useful to large scale scattering problems for which a more accurate solution is difficult to get.

Moreover, these new propagation models can be helpful in building efficiently an approximate exterior wave field used in the background of FEM for reducing pollution effects like e.g. in [START_REF] Geuzaine | An Amplitude Formulation to Reduce the Pollution Error in the Finite Element Solution of Time-Harmonic Scattering Problems[END_REF][START_REF] Antoine | Phase Reduction Models for Improving the Accuracy of the Finite Element Solution of Time-Harmonic Scattering Problems I: General Approach and Low-Order Models[END_REF]. In particular, access to high-order one-way exterior models would be valuable for accuracy improvement.

Figure 1

 1 Figure 1.

Figure 1 :

 1 Figure 1: Scattering by a two-dimensional obstacle.

Figure 2 :

 2 Figure 2: The curvilinear coordinate transform.

Figure 3 :

 3 Figure 3: The scattering problem of a circular cylinder.

Figure 4 :

 4 Figure 4: Magnitude of the scattered wave of a circular cylinder at r = 9.

Figure 5 :

 5 Figure 5: The scattering problem of an elliptic cylinder for plane incident wave with an incident angle ψ.

Figure 6 :

 6 Figure 6: Magnitude of the scattered wave of an elliptic cylinder for a plane incident wave with ψ = 35 o .

Figure 7 :

 7 Figure 7: The magnitude of u (s) m for m = 0, . . . , 60 at r = 2 (the upper one) and r = 9 (the lower one).

  Next, we compare the truncated analytic solution u (s),M with the analytic solution (E-C) of the E-C model and the fully numerical E-C BPM solution (RAtoEC). The solution E-C is also truncated to M terms and the numerical solution RAtoEC follows a starting field corresponding to the incident field truncated to M terms. These solutions are compared for M = 30 and M = 40.For M = 30, Figure8indicates that both E-C and RAtoEC coincide with the exact solution u (s),30 at r = 2 and r = 9. This, however, is not the case forM = 40as shown in Figure 9. It appears that the one-way E-C equation cannot accurately model the modes u (s) m for large m. To gain a better understanding, we compare the solutions for a single m. That is, we compare u (s) m with the corresponding analytic solution of the E-C model and the fully numerical E-C BPM solution (with a single mode incident wave). In Figure 10-12, we compare these solutions u (s) m for m = 32, 34, 35, 36, 37 and 40 for r ∈ [1, 9]. It is clear that the solutions of the E-C model are quite different from u (s) m for m = 34 and m = 35, and the rational approximations used with the E-C model give poor results for m ≥ 36. A possible explanation is that u (s) m has a strong evanescent behavior near the scatterer and this is difficult to model by the E-C equation and rational approximations. In conclusion, the accuracy of E-C one-way model is limited, since it fails to approximate all the modes that are important in the scattered wave. Nevertheless, the method does give a rough approximation very efficiently.

Figure 8 :

 8 Figure 8: Magnitude of the scattered wave of a circular cylinder for the incident field truncated to 30 terms.

Figure 9 :

 9 Figure 9: Magnitude of the scattered wave of a circular cylinder for the incident field truncated to 40 terms.

Figure 10 :

 10 Figure 10: Magnitude of the scattered waves u

  (s) m (along a certain direction) for a circular cylinder with m = 32 and 34.

Figure 11 :

 11 Figure 11: Magnitude of the scattered waves u

  (s) m (along a certain direction) for a circular cylinder with m = 35 and 36.

Figure 12 :

 12 Figure 12: Magnitude of the scattered waves u

  (s) m (along a certain direction) for a circular cylinder with m = 37 and 40.