## 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 the presence of an arbitrary number of penetrable layers. Relying on the use of certain slow-rise windowing functions, the proposed windowed Green function approach 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 multilayer structures. The proposed methodology, whose theoretical basis was presented in the recent contribution (Bruno *et al.* 2016 *SIAM J. Appl. Math.* **76**, 1871–1898. (doi:10.1137/15M1033782)), 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 the 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 the presence of 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 [1], is based on the 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 the 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 the 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. [2,3]). The Sommerfeld integrals amount to singular Fourier integrals [4,5] whose evaluation is generally quite challenging. A wide range of approaches have been proposed for evaluation of these quantities [6–14] but, as is known, all of these methods entail significant computational costs [2,7,11,12].

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 [14] based on the use of Sommerfeld integrals; the speed-ups 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 a direct multilayer integral equation which is 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 §4 and the corresponding expressions for the field evaluation are presented in §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 1*b*—in which an incoming wave is scattered by localized (bounded) surface defects and/or scattering objects in the presence of a layered medium containing a number *N*>1 of layers. For notational simplicity, our descriptions are presented in the two-dimensional case, but applications to three-dimensional configurations are presented in §7. The unperturbed configuration, which is shown in figure 1*a* for reference, consists of *N* planar layers given by *j*=1,…,*N*. The planar boundary at the interface between the layers *D*_{j} and *D*_{j+1} is denoted by *j*=1,…,*N*−1). The corresponding perturbed layers and their boundaries will be denoted by *Ω*_{j}, *j*=1,…,*N* and *Γ*_{j}, *j*=1,…,*N*−1, respectively; naturally, it is assumed that *Γ*_{j}∩*Γ*_{i}=∅.

Letting *ω*>0 and *c*_{j}>0 denote the angular frequency and the speed of sound in the layer *Ω*_{j}, the wavenumber *k*_{j} in that layer is given by *k*_{j}=*ω*/*c*_{j}. Assuming, for example, an incident plane wave of the form ** r**=(

*x*,

*y*) and where −

*π*<

*α*<0 denotes the angle of incidence measured with respect to the horizontal) and letting

*u*denote the acoustic pressure, the restrictions

*u*

_{j}=

*u*|

_{Ωj}of the total field

*u*to the domains

*Ω*

_{j}(

*j*=1,…,

*N*) satisfy the homogeneous Helmholtz equation

*ν*

_{j}denotes certain ratios of material constants. In detail, letting

*ϱ*

_{j},

*ε*

_{j}and

*μ*

_{j}denote the mass density, dielectric constant and magnetic permeability of the

*j*th layer, respectively, we have set

*ν*

_{j}=

*ϱ*

_{j}/

*ϱ*

_{j+1}for the two- and three-dimensional acoustic problems and, in the two-dimensional electromagnetic case,

*ν*

_{j}=

*ε*

_{j}/

*ε*

_{j+1}for transverse-magnetic polarization and

*ν*

_{j}=

*μ*

_{j}/

*μ*

_{j+1}for transverse-electric polarization. For definiteness, here and throughout this paper the unit normal

*n*=

*n*(

**) for**

*r***∈**

*r**Γ*

_{j}is assumed to point into

*Ω*

_{j}.

As is well known, a closed-form expression exists [4,15] for the total field *u*^{p} throughout space (*D*_{j}, *j*=1,…,*N*) that results as a plane wave *u*^{inc} impinges on the planar layer medium *j*=1,…,*N* (where the complex square root is defined in such a way that Im *k*_{jy}≥0, which, noting that *k*_{jy}≥0), the planar-medium solution *D*_{j} is given by
*A*_{j}. The amplitudes and the generalized reflection coefficients can be obtained recursively by means of the relations

## 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 *Γ*_{j}, *j*=1,…,*N*−1. As in the contribution [1], we utilize the single- and double-layer potentials
*k*_{j}. Additionally, we define the integral operators
** r** belongs to either

*Γ*

_{j}or

*Γ*

_{j−1}.

To formulate an integral equation for the unknown interface values, we define the unknown density functions *j*=1,…,*N*−1) by

A general multilayer integral formulation of the problems (2.1)–(2.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
*u*_{1}+*u*_{2} and ∂*u*_{1}/∂*n*+∂*u*_{2}/∂*n* on *Γ*_{1} from the boundary values on *Γ*_{1} of the expressions in (3.6) and their normal derivatives, and using the notations (3.4) and (3.5), we obtain the *j*=1 interface equation

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

### Remark 3.1.

In what follows equations (3.7) are expressed in terms of a single column vector function ** ϕ** (defined on the Cartesian product

*Γ*

_{j}) whose

*j*-entry equals the density pair

*j*=1,…,

*N*−1. We may thus write

With a slight notational abuse, we will write ** ϕ**=[

*ϕ*_{1},

*ϕ*_{2},…,

*ϕ*_{N−1}]

^{T}=[

*φ*

_{1},

*ψ*

_{1},

*φ*

_{2},

*ψ*

_{2},…,

*φ*

_{N−1},

*ψ*

_{N−1}]

^{T}. More generally, given arbitrary vectors

*j*=1,…,

*N*−1, we will use the ‘block-vector’ notation

Using the operators

## 4. Windowed integral equations

Following [1], in this section, we introduce rapidly convergent windowed versions of the integral formulation (3.9). To do so, we utilize the (*N*−1)×(*N*−1) block-diagonal matrix-valued window function *I* and the smooth window function
*c*<1 and where

Clearly, *η* and *w*_{A} are infinitely differentiable compactly supported functions of *x* and *t*, respectively. The support of the window function *w*_{A}=*w*_{A}(*x*) as a function of *c*, which controls the steepness of the rise of the window function *w*_{A}, is not displayed as part of the notation *w*_{A}.

(Different values *A*_{j} of the window-size *A*, *j*=1,…,*N*−1, could in principle be used for the various layer interfaces and corresponding block entries in (4.1)—possibly utilizing smaller (respectively, larger) values *A*_{j} in higher (respectively, 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 *A* is used for all the interfaces.)

To produce a windowed version of equation (3.9), the integrand is multiplied by the window matrix *Γ*_{A}={(*r*_{1},…,*r*_{N−1})∈*Γ*:*r*_{j}=(*x*_{j},*y*_{j}) and *N*−1)×2(*N*−1), the exact relation

A successful implementation of the WGF idea requires use of an accurate substitute for the quantity *Γ*_{A}, which does not depend on knowledge of the unknown density ** ϕ** on the complement of

*P*

_{j}depicted in figure 1

*a*. Since

*Γ*differs from

*b*), we clearly have

*P*

_{j}(

*j*=1,…,

*N*−1) that are associated with the planarly layered medium

*P*. As shown in [1], letting

*j*=1,…,

*N*−1), substitution of

**by**

*ϕ*

*ϕ*^{p}on the right-hand side of (4.4) results in errors that decay super-algebraically fast as

*Γ*

_{A}wherein the window function

*w*

_{A}equals one. Indeed, even though

**may differ significantly from**

*ϕ*

*ϕ*^{p}, the corresponding integrated terms result in super-algebraically small errors, as may be checked via stationary phase analysis ([1] for details). We thus obtain the super-algebraically accurate windowed integral equation system

*ϕ*^{w}, which we re-express in the form

*Γ*

_{A}to produce the term

**∈**

*x**Γ*

_{A}.

A closed-form expression for *C*=*P*_{j} (*j*=1,…,*N*−1), equations (A 15) yield the desired relations
*j*=1,…,*N*−1, where *δ*_{i,j} denotes the Kronecker delta symbol.

As demonstrated in §7 through a variety of numerical examples, the vector density function *ϕ*^{w}, which is the solution of the windowed integral equation (4.8), converges super-algebraically fast to the exact solution ** ϕ** of (3.9) within

*A*>0 increases. This observation can be justified via arguments analogous to those presented in [1].

### Remark 4.1.

The difference *T*_{j} of

## 5. Near-field evaluation

This section presents a super-algebraically accurate WGF approximation *u*^{w} of the solution *u* of (2.1)–(2.2) near the localized defects. To obtain this approximation, we consider the ‘defect’ field
*u*_{j} and the planar-structure total field
*Ω*_{j} by the expressions on the right-hand side of equation (2.3).

Subtracting the integral representation
*ϕ*_{j}=[*φ*_{j},*ψ*_{j}]^{T} of the integral equation (3.9) together with the planar-structure total fields
*j*, 1≤*j*≤*N*, the functions *f*_{j} and *g*_{j} vanish outside the *j*th portion *Γ*_{j}∖*Π*_{j} of the boundary of the localized defects.

Relying on the WGF solutions *ϕ*^{w} of equation (4.8) and applying a windowing procedure similar to the one used in the previous section, a highly accurate approximation to the defect near-field (5.4) results. In detail, substitution of *φ*_{j} by *ψ*_{j} by *u*_{j} then follows from (5.1):

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

## 6. Far-field evaluation

As indicated in the previous section, formulae (5.7)–(5.8) only provide uniformly accurate approximations within bounded subsets of *j*=1,…,*N*) have been obtained within *u* can be obtained by applying certain Green-type formulae on a bounding curve *S*, such as the one depicted in figure 2, which encloses all of the local defects, and which is contained within *u*^{d}=*u*^{d}(** r**) to equal

**∈**

*r**Ω*

_{j}(

*j*=1,…,

*N*), use of a Green identity based on the

*N*-layer Green function

*H*over the region exterior to

*S*leads to the integral representation [3], Lemma 4.2.6

**everywhere outside**

*r**S*. Note that the necessary values of

*S*can be computed by means of (5.8)—since, by construction,

*S*lies inside the region where (5.8) provides an accurate approximation of the field

The far-field approximation *u*^{f} of the defect field *u*^{d} as *H* and its normal derivative ∂*H*/∂*n*_{r′} in (6.1) by the respective first-order *H*^{f} and ∂*H*^{f}/∂*n*_{r′}—which can be obtained for the *N*-layer case (as illustrated in [1,3,4,17,18] for *N*=2 and below in this section for *N*=3) by means of the method of steepest descents. (The fact that the far field of the function ∂*H*/∂*n*_{r′} coincides with ∂*H*^{f}/∂*n*_{r′} can be verified by direct inspection of these two quantities.) The far-field *u*^{f} is thus given by
*H* itself, the corresponding far-field *H*^{f} 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 *H*^{f} for a slab—that is, a three-layer medium with wavenumbers *k*_{j}, *j*=1,2,3 where *k*_{1}=*k*_{3}—in two-dimensional space. We assume the case *k*_{2}>*k*_{1} for which the slab can sustain guided modes that propagate along the *x*-axis. To evaluate *H*^{f}, we first note that, for a source point ** r**′=(

*x*′,

*y*′)∈

*D*

_{j}and a target point

*θ*∈[0,

*π*]), the layer Green function

*H*is given by the contour integral [3,4,17,18]

*j*=1,2,3, we have set

The determination of physically admissible branches of the functions *ξ*=*k*_{1}=*k*_{3} is in fact the only branch cut in the domain of definition of the functions *p*_{j}(*ξ*,** r**′)/

*f*(

*ξ*) (

*j*=1,2,3), as it can be shown that these are even functions of

*γ*

_{2}[4]. The branch cuts and Sommerfeld contour SC utilized in (6.3) are depicted in figure 3.

As suggested above, in order to obtain the far-field form of the layer Green function *H*^{f}, we resort to the method of steepest descents [17]. Analysis of the phase function *ϕ* (6.4a) readily shows that there is only one saddle point on the real axis at *ϕ*(*ξ*)=*k*_{1}, also intersects the real axis at *γ*_{1} it can shown that
*ξ* on SD. This information suffices to sketch the paths of steepest descent that are displayed in figure 3.

To produce asymptotic expansions of the integrals (6.3), we then proceed to deform the Sommerfeld contour (SC) to the steepest descent contour SD (figure 3). Considering the saddle point at *ξ*_{0} and taking into account the poles of the integrand *p*_{j}(*ξ*,** r**′)/

*q*(

*ξ*) at the points

*ξ*

_{p}which are enclosed by the curves SD and SC, we obtain the following expression for the far-field form of the Green function for the two-dimensional slab:

Clearly, as indicated above, the far-field asymptotics *H*^{f} of the layer Green function *H*, and thus its normal derivative, can be evaluated inexpensively by means of a simple explicit expression.

## 7. Numerical examples

This section presents a set of two- and three-dimensional numerical examples that demonstrate the character of the proposed multilayer WGF methodology. For the sake of definiteness, a window function *w*_{A} (4.2) with *c*=0.7 was used in all cases. Numerical errors were evaluated by resorting to numerical-convergence studies and/or increases in the window size *A*. As additional references, in some cases adequately accurate solutions obtained by the Sommerfeld layer–Green-function (LGF) method [3,14] (with accuracy evaluated by means of convergence studies) were used to evaluate the accuracy of the WGF approach. Brief indications will be provided when necessary to indicate which method is being used in each case. The two-dimensional results were obtained via solution of the integral equation system (4.8) by means of the Nyström method described in [19], §3.5. The three-dimensional solutions, in turn, were obtained by means of the algorithm presented in [20].

Our first example concerns the structure depicted in figure 4*a*, in which semicircular defects of radii *a*=1 are placed at the planar interfaces *k*_{1}=10, *k*_{2}=20 and *k*_{3}=30. Figure 4*b* displays the maximum relative errors (in log–log scale) in the total field produced by the WGF method on the surface of the semicircular defects (the curves marked in red in figure 4*a*) for various windows sizes *A*>0 and incidences *α*. The number of quadrature points was selected in such a way that, for any given *A*>0, the Nyström discretization error in the integral equation solution is not larger than 10^{−9}. The WGF solution obtained for *A*=32*λ* is utilized as the reference for the error estimation. As it can be inferred from the error curves displayed in figure 4, super-algebraic convergence is observed as *A* increases. In particular, these results demonstrate the uniformly fast convergence exhibited by the WGF method as the incidence angles approach grazing.

To compare the computational cost of the LGF method [14] and proposed WGF method for a given accuracy, we consider a planar three-layer structure similar to those considered previously, but now containing only one surface defect: a semicircular cavity of radius *a*=1 at *P*_{1}. (The use of a single defect reduces somewhat the LGF cost, which seemed inordinately large for the two-defect problem.) A plane-wave *u*^{inc} with *α*=−*π*/6 illuminates the structure. Five sets of wavenumbers given by *k*_{1}=*κ*, *k*_{2}=2*κ* and *k*_{3}=3*κ* with *κ*=2^{j}, *j*=1,…,5 are considered. The resulting problems of scattering are then solved by employing a Nyström discretization of the WGF equations (4.8), and a numerical version of (5.8) is used to evaluate near fields. The same problem of scattering is then solved, with a relative error not larger than 10^{−4}, by means of a generalization to the present three-layer case, of the two-layer LGF method presented in [14] (see also [3]). The reference solution used to estimate the accuracy of the LGF solution is obtained by solving the resulting LGF integral equation with an error not larger 10^{−9} (this accuracy is achieved by utilizing a large number of Nyström quadrature points and evaluating the layer Green function with an error not larger than 10^{−10}).

Table 1 displays the computing times needed by both methods to construct the system matrices. To allow for a fair comparison of the computing times and the field values on the surface defect, the same set of quadrature points is utilized to discretize the currents on the surface of the cavity in each case. The number of quadrature points was increased in direct proportion to the value of *κ*. The maximum of the absolute value of the difference between the LGF and WGF solutions (using *A*=8*λ*) on the surface of the defect is no larger than 10^{−4} in all the examples considered. Remarkably, in the *κ*=32 case the proposed WGF method is 260 times faster than the LGF method.

Figure 5 presents a comparison of the near fields obtained by means of the WGF and LGF methods for some of the test cases considered in table 1. The first and second columns in figure 5 display the real part of the total near fields produced by the WGF method (first column) and by the LGF method (second column), respectively, for *κ*=32. The fields are evaluated in the rectangular region *A*=16*λ* had to be used to produce accurate near fields throughout the plotted region.

Figure 6, in turn, compares the far-field patterns
*k*_{1}=*k*_{3}=10 and *k*_{2}=15. The WGF far-field pattern (blue solid line) was obtained by letting *u*=*u*^{f} in (7.1), where *u*^{f} is given by (6.2) with WGF defect fields *Ω*_{j}, *j*=1,…,*N* (equation (5.8)). The corresponding LGF far-field pattern (red dots) was obtained on the basis of a highly accurate LGF solution together with the far-field asymptotics of the layer Green function [14]. We have verified that, as expected, the accuracy of the WGF far-field patterns is comparable to the accuracy of the corresponding defect fields

Figure 7 displays near-fields resulting from the WGF method, with window size *A*=12*λ*, for a structure consisting of nested circular surface defects in a nine-layer medium with planar interfaces *j*=1,…,8. The corresponding wavenumbers are *k*_{2j−1}=15 for *j*=1,…,5 and *k*_{2j}=30 for *j*=1,…4. The structure is illuminated by plane waves with two different incidence angles. A 112-s overall computing time sufficed to evaluate each one of the two near-fields displayed. Note the resonance that takes place in the third upper and lower rings in figure 7*b*.

Figure 8, finally, presents applications of the WGF methodology to the problem of scattering by three-dimensional structures in the presence of layer media. The two-dimensional descriptions presented in §§2 through 6 extend directly to the present three-dimensional context.

## Authors' contributions

O.P.B. and C.P.-A. contributed to all aspects of this article.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by NSF and AFOSR through contracts DMS-1411876 and FA9550-15-1-0043, and by the NSSEFF Vannevar Bush Fellowship under contract number N00014-16-1-2808.

## Appendix A. Integral representation based on non-windowed free-space Green functions

This section presents an integral representation formula, based on the *free-space Green function*, for fields of the form *v*(** r**)=

*v*

_{j}(

**) for**

*r***∈**

*r**Ω*

_{j},

*j*=1,…,

*N*, where, letting

*N*-layer structures is straightforward.

Our derivations utilize three local polar-coordinate systems, each one of which is associated with one of the layers *Ω*_{j}. These coordinate systems are centred at (0,−*d*_{1}), (0,−(*d*_{1}+*d*_{2})/2) and (0,−*d*_{2}) and, thus, the radial variables are given by
*x* and *y*, as illustrated in figure 9.

Additionally, some of the subsequent derivations utilize the decomposition
*q*_{j}=*A*_{j} are expressed in terms of the amplitudes *A*_{j} and the generalized reflection coefficients *β*_{2}=1 and *β*_{j}=2/3 for *j*=1,3 (see (A 2)), verify
*m*=1,…,*M*_{j}. The propagation constants *Ω*_{2} is *Ω*_{1} (respectively, *Ω*_{3}), on the other hand, we have *Γ*_{1} (respectively, *Γ*_{2}) and decays exponentially fast towards the interior of *Ω*_{1} (respectively, *Ω*_{3}).

We are now in a position to derive the desired integral representation for the total fields *v*_{j} in (A 1). Our derivations consider at first the bounded domains
*R*>0 is large enough that *B*_{R} contains all of the surface defects, as illustrated in figure 9. Our bounded-domain calculations use the curves *n*, as depicted in figure 9. To facilitate repeated use of Green’s third identity in our derivations, we follow [22] and letting *C*,
** r**−

**′|=|**

*r***′|−**

*r***′⋅**

*r***/|**

*r***′|+**

*r**O*(|

**′|**

*r*^{−1}).

With reference to figure 9, and in view of Green’s third identity applied to *v*_{j} in (A 1). Since, in view of equations (A 3)–(A 5), these functions can be expressed as linear combinations of *N*-layer versions (A 15) then follow directly from the limiting expressions thus found.

Case *j*=2 for *v*_{2} remains bounded throughout *Ω*_{2} (as it follows from equation (A 7)), we conclude that the term *I*_{2} involving the integral over

Case *j*=1,3 for *I*_{j} that involve integrals over the semicircular curves *j*=1 and *j*=3, we have |** r**′|=

*R*+

*O*(1) and

*θ*

_{j}are as shown in figure 9. Since

*v*

^{rad}

_{j}in (A 5) for

*j*=1,3 satisfies the Sommerfeld radiation condition (A6), utilizing standard arguments [16], it can be shown that

Case *j*=1,3 for *m*=1,…,*M*_{j}), where *θ*_{1}∈[0,*π*], we obtain
*θ*∈[0,*π*], and we thus conclude that *j*=1,3, as

Case *j*=1 for *α*∈(−*π*,0), we have 0<*θ*_{1}−*α*<2*π*. Thus, there is only one point of stationary phase within the domain of integration at *θ*_{1}=*α*+*π*. A straightforward application of the method of stationary phase [17] then yields
*R*^{−1}.)

Case *j*=3 for *α*∈(−*π*,−*π*/2] or *α*∈(−*π*/2,0)) and (c) *α*∈(−*π*/2,0). Under this assumption, *O*(*R*^{−1/2}). The stationary point at *θ*=−*π* in the second integral, on the other hand, leads to
*α*∈(−*π*,*π*/2], it can be shown that *α*′∈(−*π*,0) is determined by the Snell’s law

Therefore, taking the limit as *N*=3, and where *δ*=0 if *δ*=1 if

### Remark A.1.

The total field representation (A 13) can easily be extended to problems of scattering by defects in the presence of layer media composed by *N*>3 layers; the result is

- Received March 2, 2017.
- Accepted May 11, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.