Windowed Green Function method for the Helmholtz equation in presence of multiply layered media

# Windowed Green Function method for the Helmholtz equation in presence of multiply layered media

Oscar P. Bruno obruno@caltech.edu Computing & Mathematical Sciences, California Institute of Technology. Carlos Pérez-Arancibia Department of Mathematics, Massachusetts Institute of Technology.
July 12, 2019
###### Abstract

This paper presents a new methodology for the solution of problems of two- and three-dimensional acoustic scattering (and, in particular, two-dimensional electromagnetic scattering) by obstacles and defects in presence an arbitrary number of penetrable layers. Relying on use of certain slow-rise windowing functions, the proposed Windowed Green Function approach (WGF) efficiently evaluates oscillatory integrals over unbounded domains, with high accuracy, without recourse to the highly expensive Sommerfeld integrals that have typically been used to account for the effect of underlying planar multi-layer structures. The proposed methodology, whose theoretical basis was presented in the recent contribution (SIAM J. Appl. Math. 76(5), p. 1871, 2016), is fast, accurate, flexible, and easy to implement. Our numerical experiments demonstrate that the numerical errors resulting from the proposed approach decrease faster than any negative power of the window size. In a number of examples considered in this paper the proposed method is up to thousands of times faster, for a given accuracy, than corresponding methods based on use of Sommerfeld integrals.

## 1 Introduction

This paper presents a new methodology for the solution of problems of acoustic scattering by obstacles and defects in presence an arbitrary number of penetrable layers in two and three-dimensional space; naturally, the two-dimensional Helmholtz solvers also apply, by mathematical analogy, to corresponding two-dimensional electromagnetic scattering problems. This “Windowed Green Function” (WGF) method, whose theoretical basis was presented in the recent contribution [4], is based on use of smooth windowing functions and integral kernels that can be expressed directly in terms of the free-space Green function, and, importantly, it does not require use of expensive Sommerfeld integrals. The proposed methodology is fast, accurate, flexible, and easy to implement. Our experiments demonstrate that, as predicted by theory, the numerical errors resulting from the proposed approach decrease faster than any negative power of the window size. In a number of examples considered in this paper the proposed method is up to thousands of times faster, for a given accuracy, than corresponding methods based on use of Sommerfeld integrals.

The classical layer Green functions and associated Sommerfeld integrals automatically enforce the relevant transmission conditions on the unbounded flat surfaces and thus reduce the scattering problems to integral equations on the obstacles and/or defects (cf. [19, 24]). The Sommerfeld integrals amount to singular Fourier integrals [8, 25] whose evaluation is generally quite challenging. A wide range of approaches have been proposed for evaluation of these quantities [1, 6, 7, 12, 13, 18, 21, 22, 23] but, as is known, all of these methods entail significant computational costs [6, 18, 19, 21].

The WGF approach proceeds as follows. The integral equation formulations of the scattering problems under consideration, which are at first posed on the complete set of material interfaces (including all unbounded interfaces), are then smoothly truncated to produce an approximating integral-equation system posed over bounded integration domains that include the surface defects and relatively small portions of the flat interfaces. The integral-operator truncation is effected by means of a certain slow-rise smooth window function which, importantly, gives rise to solution errors which decrease faster than any negative power of the window size. In practice the proposed solution method is up to thousands of times faster, for a given accuracy, than corresponding methods [23] based on use of Sommerfeld integrals; the speedups in evaluations of near fields are even more significant, in view of the large computing times required for evaluation of Sommerfeld integrals near the planar interface.

This paper is organized as follows. Section 2 presents a description of the multilayer scattering problem. Section 3 then presents two types of direct multi-layer integral equations for a physical, which can be obtained by means of a generalized version of Green’s third identity (which is itself derived in Appendix A). The windowed integral equations are derived in Section 4 and the corresponding expressions for the field evaluation are presented in Section 5. Section 7, finally, presents a variety of numerical examples which demonstrate the super-algebraic convergence and the efficiency of the proposed approach.

## 2 Preliminaries

We consider the problem of scattering of an acoustic incoming wave by a two- or three-dimensional configuration such as the one depicted in Figure (b)b—in which an incoming wave is scattered by localized (bounded) surface defects and/or scattering objects in presence of a layered medium containing a number of layers. For notational simplicity our descriptions are presented in the two-dimensional case, but applications to three-dimensional configurations are presented in Section 7. The unperturbed configuration, which is shown in Figure (a)a for reference, consists of planar layers given by for . The planar boundary at the interface between the layers and is denoted by (). The corresponding perturbed layers and their boundaries will be denoted by , and , respectively; naturally, it is assumed that .

Letting and denote the angular frequency and the speed of sound in the layer , the wavenumber in that layer is given by . Assuming e.g. an incident plane wave of the form (where and where denotes the angle of incidence measured with respect to the horizontal) and letting denote the acoustic pressure, the restrictions of the total field to the domains () satisfy the homogeneous Helmholtz equation

 Δuj+k2juj=0inΩj,j=1,…,N, (1)

together with the transmission conditions

 uj=uj+1and∂uj∂n=νj∂uj+1∂nonΓj,  j=1,…,N−1, (2)

where where denotes the fluid density in . For definiteness, here and throughout this paper the unit normal for is assumed to point into .

As is well known, a closed form expression exists [8, 26] for the total field throughout space ( in , ), that results as a plane wave impinges on the planar layer medium . In detail, letting and , (where the complex square-root is defined in such a way that , which, noting that , requires as well), the planar-medium solution in is given by

 upj(x,y)=Ajeik1xx{e−ikjyy+˜Rj,j+1eikjy(y+2dj)}inDj, 1≤j≤N, (3)

in terms of certain generalized reflection coefficients and amplitudes . The amplitudes and the generalized reflection coefficients can be obtained recursively by means of the relations

 Aj = ⎧⎪ ⎪⎨⎪ ⎪⎩1ifj=1,Tj−1,jAj−1ei(kj−1,y−kjy)dj−11−Rj,j−1˜Rj,j+1e2ikjy(dj−dj−1)ifj=2,…,N, (4)

and

 ˜Rj,j+1 = ⎧⎪ ⎪⎨⎪ ⎪⎩0ifj=N,Rj,j+1+Tj+1,j˜Rj+1,j+2Tj,j+1e2ikj+1,y(dj+1−dj)1−Rj+1,j˜Rj+1,j+2e2ikj+1,y(dj+1−dj)ifj=N−1,…,1, (5)

in terms of the reflection and transmission coefficients

 Rj,j+1=kjy−νjkj+1,ykjy+νjkj+1,yandTj,j+1=2kjykjy+νjkj+1,y,

respectively.

## 3 Integral equation formulations

This section presents an integral equation for the unknown interface values of the total field and its normal derivative from below, at each one of the interfaces , . As in the contribution [4] we utilize the single- and double-layer potentials

 Stj[ϕ](r)=∫Γj−1Gkj(r,r′)ϕ(r′)dsr′,Dtj[ϕ](r)=∫Γj−1∂Gkj∂nr′(r,r′)ϕ(r′)dsr′,Sbj[ϕ](r)=∫ΓjGkj(r,r′)ϕ(r′)dsr′,Dbj[ϕ](r)=∫Γj∂Gkj∂nr′(r,r′)ϕ(r′)dsr′, (6)

which are defined for and are expressed in terms of improper integrals whose convergence is conditioned upon the oscillatory behavior of the integrand. Here we have called the free-space Green function for the Helmholtz equation with wavenumber . Additionally, we define the integral operators

 Ktj[ϕ](r)=∫Γj−1∂Gkj∂nr(r,r′)ϕ(r′)dsr′Ntj[ϕ](r)=∫Γj−1∂2Gkj∂nr∂nr′(r,r′)ϕ(r′)dsr′,Kbj[ϕ](r)=∫Γj∂Gkj∂nr(r,r′)ϕ(r′)dsr′Nbj[ϕ](r)=∫Γj∂2Gkj∂nr∂nr′(r,r′)ϕ(r′)dsr′, (7)

where the evaluation point belongs to either or .

In order to formulate an integral equation for the unknown interface values we define the unknown density functions and () by

 φj=uj+1andψj=∂uj+1∂nonΓj. (8)

Additionally we define the vector density functions

 ϕj=[φj,ψj]T,ϕinc=[uinc∣∣Γ1,∂uinc∂n∣∣Γ1]Tandϕ∥=⎡⎣u∥N∣∣ΓN−1,∂u∥N∂n∣∣ΓN−1⎤⎦T (9)

where is defined in (52), and the matrix operators

 Ej=[1001+νj2],Tj=[Dtj+1−Dbj−Stj+1+νjSbjNtj+1−Nbj−Ktj+1+νjKbj], Lj=[Dtj−StjNtj−Ktj]and Rj=[−Dbj+1νj+1Sbj+1−Nbj+1νj+1Kbj+1] (all operators % evaluated at observation points r on Γj). (10)

A general multi-layer integral formulation of the problem (1)–(2) can now be obtained in terms of these densities and operators. Indeed, as is shown in Appendix A, the fields within the layers admit the integral representations

 u1(r)=Db1[φ1](r)−ν1Sb1[ψ1](r)+uinc(r),uj(r)=Dbj[φj](r)−νjSbj[ψj](r)−Dtj[φj−1](r)+Stj[ψj−1](r),j=2,…,N−1,uN(r)=−DtN[φN−1](r)+StN[ψN−1](r)+u∥N(r), (11)

in terms of the interface values (8). Therefore, evaluating and on from the boundary values on of the expressions in (11) and their normal derivatives, and using the notations (9) and (10), we obtain the interface equation

 E1ϕ1+T1[ϕ1]+R1[ϕ2]=ϕinconΓ1. (12a) A similar procedure yields the integral equations Ejϕj+Lj[ϕj−1]+Tj[ϕj]+Rj[ϕj+1]=0onΓj,j=2,…,N−2 (12b) and EN−1ϕN−1+LN−1[ϕN−2]+TN−1[ϕN−1]=ϕ∥onΓN−1. (12c)

(Note that, of course, the calculations leading to equations (12) rely on the well-known jump relations for the single- and double-layer potentials and their normal derivatives [11].)

###### Remark 3.1.

In what follows equations (12) are expressed in terms of a single column vector function (defined on the Cartesian product of the curves ) whose -entry equals the density pair for . We may thus write

 ϕ=[ϕ1,ψ1,φ2,ψ2,⋯,φN−1,ψN−1]T:Γ→C2(N−1).

Similarly we define

 ϕinc=⎡⎣uinc|Γ1,∂uinc∂n∣∣Γ1,0,0,⋯,0,0,u∥N|ΓN−1,∂u∥N∂n∣∣ΓN−1⎤⎦T:Γ→C2(N−1).

With a slight notational abuse we will write . More generally, given arbitrary vectors for we will use the “block-vector” notation .

Using the operators

 E=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣E1E2⋱EN−1⎤⎥ ⎥ ⎥ ⎥ ⎥⎦andTΓ=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣T1R1L2T2R2L3⋱⋱⋱⋱RN−2LN−1TN−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (13)

together with the notations introduced in Remark 3.1, equations (12) can be expressed in the form

 (14)

## 4 Windowed integral equations

Following [4], in this section we introduce rapidly-convergent windowed versions of the integral formulation (14). In order to do so we utilize the block-diagonal matrix-valued window function given by

 WA(r1,r2,…,rN−1)=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣wA(x1)IwA(x2)I⋱wA(xN−1)I⎤⎥ ⎥ ⎥ ⎥ ⎥⎦,rj=(xj,yj)∈Γj, (15)

in terms of the two-by-two identity matrix and the smooth window function

 wA(x)=η(x/A;c,1), (16)

where and where

 η(t;t0,t1)=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩1,|t|≤t0,exp(2e−1/uu−1),t0<|t|t1. (17)

Clearly and are infinitely differentiable compactly-supported functions of and , respectively. The support of the window function as a function of equals the set . Note that the parameter , which controls the steepness of the rise of the window function , is not displayed as part of the notation .

(While different values of the window-size , , could in principle be used for the various layer interfaces and corresponding block entries in (15)—possibly utilizing smaller (resp. larger) values in higher (resp. lower) frequency layers, and therefore reducing the overall number of unknowns required for the WGF method to produce a given accuracy. For simplicity, however, throughout this paper a single window-size value is used for all the interfaces.)

In order to produce a windowed version of equation (14) we proceed in two stages. At first the integrand is multiplied by the window matrix and the equation is restricted to the windowed region —so that, moving the remainder of the windowed integral operator to the right-hand side and letting denote the identity matrix of dimension , the exact relation

 Eϕ+TΓ[WAϕ]=ϕinc−TΓ[(I−WA)ϕ]onΓA (18)

results. Note that defining we have

 ΓA=N−1∏j=1Γj,A. (19)

A successful implementation of the WGF idea requires use of an accurate substitute for the quantity throughout , which does not depend on knowledge of the unknown density (cf. [4]). In order to obtain such an approximation we introduce an operator which is defined just like in (13) but in terms of potentials (6) and operators (7) given by integrals on the flat interfaces depicted in Figure (a)a. Since vanishes wherever differs from , we clearly have . Additionally, we consider the aforementioned scalar densities and on () that are associated with the planarly layered medium . As shown in [4], letting (), substitution of by on the right-hand-side of (21) results in errors that decay super-algebraically fast as within the subset

 ˜ΓA=ΓA∩N−1∏j=1{(xj,yj)∈Γj:wA(xj)=1} (20)

of wherein the window function equals one. Indeed, even though may differ significantly from , the corresponding integrated terms result in super-algebraically small errors, as it may be checked via stationary phase analysis (see [4], for details). We thus obtain the super-algebraically-accurate windowed integral equation system

 Eϕw+TΓ[WAϕw]=ϕinc−TP[(I−WA)ϕp]onΓA (21)

which we re-express in the form

 Eϕw+TΓ[WAϕw]=ϕinc+TP[WAϕp]−TP[ϕp]onΓA. (22)

As shown in what follows, the right-hand term in (22) can be expressed in closed form, and thus, using numerical integration over the bounded domain to produce the term , the complete right-hand side can be efficiently evaluated for any given .

A closed-form expression for (cf. Remark 3.1) can be obtained via an application of Green’s formula: using (48) with (), equations (53) yield the desired relations:

 μj=δ1,jϕinc+δN−1,jϕ∥−⎧⎪⎨⎪⎩EjϕpjonΓj∩Pj,[up∇up⋅n]onΓj∩(Dj∪Dj+1), (23)

for where denotes the Kronecker delta symbol.

As demonstrated in Section 7 through a variety of numerical examples, the vector density function , which is the solution of the windowed integral equation (22), converges super-algebraically fast to the exact solution of (14) within as the window size increases. This observation can be justified via arguments analogous to those presented in [4].

###### Remark 4.1.

The difference of hypersingular operators that appears in the definition of the diagonal blocks of  is in fact a weakly singular integral operators (cf. [11, Sec. 3.8]).

## 5 Near-field evaluation

This section presents a super-algebraically accurate WGF approximation of the solution of (1)–(2) near the localized defects. In order to obtain this approximation we consider the “defect” field

 udj=uj−~upjinΩj(j=1,…,N−1), (24)

given by the difference between the total field and the planar-structure total field

 ~upj(x,y)=Ajeik1xx{e−ikjyy+˜Rj,j+1eikjy(y+2dj)}inΩj, 1≤j≤N. (25)

Note that is given in by the expressions on the right-hand side of equation (3).

Subtracting the integral representation

 ~up1(r)=Db1[~φp1+f1](r)−ν1Sb1[~ψp1+g1](r)+uinc(r),~upj(r)=Dbj[~φpj+fj](r)−νjSbj[~ψpj+gj](r)−Dtj[~φpj−1](r)+Stj[~ψpj−1](r),j=2,…,N−1,~upN(r)=−DtN[~φpN−1](r)+StN[~ψpN−1](r)+u∥N(r), (26)

—which follows as equation (53) is applied to —from the integral representation (11) we obtain the exact integral relations

 ud1(r)=Db1[φ1−~φp1−f1](r)−ν1Sb1[ψ1−~ψp1−g1](r),udj(r)=Dbj[φj−~φpj−fj](r)−νjSbj[ψj−~ψpj−gj](r)−Dtj[φj−1−~φpj−1](r)+Stj[ψj−1−~ψpj−1](r),j=2,…,N−1,udN(r)=−DtN[φN−1−~φpN−1](r)+StN[ψN−1−~ψpN−1](r) (27)

for the defect fields. These relations can be used to evaluate the defect fields in terms of the solution  of the integral equation (14) together with the planar-structure total fields

 ~ϕpj=[~φpj,~ψpj]T,where~φpj=~upj+1|Γjand~ψpj=∂~upj+1∂n∣∣Γj, (28)

and the jumps

 \uppsij=[fj,gj]T,wherefj=~upj−~upj+1andgj=1νj∂~upj∂n−∂~upj+1∂nonΓj. (29)

Note that, importantly, for each , , the functions and vanish outside the -th portion of the boundary of the localized defects.

Relying on the WGF solutions of equation (22) and applying a windowing procedure similar to the one used in the previous section, a highly-accurate approximation to the defect near-fields (27) results. In detail, substitution of by and by in (27) yields the approximate expressions

 ud,w1(r)=Db1[wA(φw1−~φp1)−f1]−ν1Sb1[wA(ψw1−~ψp1)−g1],ud,wj(r)=Dbj[wA(φwj−~φpj)−fj]−νjSbj[wA(ψwj−~ψpj)−gj]−Dtj[wA(φwj−1−~φpj−1)]+Stj[wA(ψwj−1−~ψpj−1)],j=2,…,N−1,ud,wN(r)=−DtN[wA(φwN−1−~φpN−1)]+StN[wA(ψwN−1−~ψpN−1)], (30)

for the defect field . The desired approximation for the total field then follows from (24):

 uwj=~upj+ud,wjinΩj(j=1,…,N). (31)

Formulae (31) provide super-algebraically accurate approximations of the total near-fields within the region

 ˜ΩA=N⋃j=1Ωj∩{r∈R2:wA(x)=1} (32)

containing the localized defects—with uniformly small errors, as , within every bounded subset of . A theoretical discussion in these regards (for the two-layer case) can be found in [4] (see e.g. Remark 4.1 in that reference).

## 6 Far-field evaluation

As indicated in the previous section, formulae (30)–(31) only provide uniformly accurate approximations within bounded subsets of . But, once accurate defect fields () have been obtained within , correspondingly accurate far-field values for the solution can be obtained by applying certain Green-type formulae on a bounding curve , such as the one depicted in Figure 2, which encloses all of the local defects, and which is contained within . In detail, defining the defect field to equal for (), use of a Green identity based on the -layer Green function  over the region exterior to leads to the integral representation [24, Lemma 4.2.6]

 ud(r)=∫S{∂H∂nr′(r,r′)ud(r′)−H(r,r′)∂ud∂n(r′)}dsr′, (33)

which is valid for everywhere outside . Note that the necessary values of and their normal derivatives on can be computed by means of (31)—since, by construction, lies inside the region where (31) provides an accurate approximation of the field .

The far-field approximation of the defect field as in any direction is then obtained by replacing the layer Green function and its normal derivative in (33) by the respective first-order asymptotic expansions and —which can be obtained for the -layer case (as illustrated in [2, 8, 24, 4, 3] for and below in this section for ) by means of the method of steepest descents. (The fact that the far field of the function coincides with can be verified by direct inspection of these two quantities.) The far field is thus given by

 uf(r)=∫S{∂Hf∂nr′(r,r′)ud(r′)−Hf(r,r′)∂ud∂n(r′)}dsr′. (34)

It is important to note that, unlike the layer Green function itself, the corresponding far-field and its normal derivative can be evaluated inexpensively by means of explicit expressions.

As an example we sketch here the calculation of the far-field for a slab—that is, a three-layer medium with wavenumbers , where —in two-dimensional space. We assume the case for which the slab can sustain guided modes that propagate along the -axis. In order to evaluate we first note that, for a source point and a target point (), the layer Green function is given by the contour integral [2, 8, 24, 3]

 Hj(r,r′) = 14π∫SCpj(ξ,r′)q(ξ)e|r|ϕ(ξ)dξ. (35)

Here, letting , , we have set

 ϕ(ξ) = iξcosθ−γ1(ξ)sinθ, (36a) p1(ξ,r′) = {R