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

TL;DR

This paper introduces a fast, accurate, and flexible Windowed Green Function method for solving acoustic and electromagnetic scattering problems in layered media, avoiding expensive integral computations.

Contribution

The paper develops a novel WGF approach that efficiently evaluates oscillatory integrals in layered media without Sommerfeld integrals, improving speed and accuracy.

Findings

Numerical errors decrease faster than any negative power of window size.

Method is up to thousands of times faster than traditional Sommerfeld integral methods.

High accuracy achieved with less computational expense.

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.

Tables

Table 1. Table 1: Computing times required by the WGF and LGF methods to construct the system matrices for the numerical solution of the problem of scattering of a plane-wave by a semi-circular cavity or radius a = 1 𝑎 1 a=1 on a three-layer medium with wavenumbers k 1 = κ subscript 𝑘 1 𝜅 k_{1}=\kappa , k 2 = 2 ​ κ subscript 𝑘 2 2 𝜅 k_{2}=2\kappa and k 3 = 3 ​ κ subscript 𝑘 3 3 𝜅 k_{3}=3\kappa , with κ = 2 j 𝜅 superscript 2 𝑗 \kappa=2^{j} , j = 1 , … ​ 5 𝑗 1 … 5 j=1,\ldots 5 . All the two-dimensional runs reported in this paper were performed using a Matlab implementation of our algorithms in a MacBook Air laptop (early 2014 model).

	WGF method					LGF method
$\kappa$	$2$	$4$	$8$	16	32	$2$	$4$	$8$	16	32
Number of	1232	1272	1348	1496	1800	68	148	300	596	1204
unknowns
Matrix	3.44	3.53	3.98	5.78	7.29	6.49	22.15	82.86	319.46	$1900$
construction (s)

Equations

$$\Delta u_j+k_j^2{u}_j=0\text{mbox}in{\Omega }_j,j=1,\dots ,N,$$

$$u_j=u_{j+1}\text{mbox}and\frac{\partial u_j}{\partial n}={\nu }_j\frac{\partial u_{j+1}}{\partial n}\text{mbox}on{\Gamma }_j,j=1,\dots ,N-1,$$

$$u_j^p(x,y)=A_j{e}^{i{k}_{1x}x}\left\{e^{-i{k}_{jy}y}+{\tilde{R}}_{j,j+1}{e}^{i{k}_{jy}(y+2d_j)}\right\}\text{mbox}in{D}_j,1\le j\le N,$$

$$A_j$$

$${\tilde{R}}_{j,j+1}$$

$$R_{j,j+1}=\frac{k_{jy}-{\nu }_j{k}_{j+1,y}}{k_{jy}+{\nu }_j{k}_{j+1,y}}\text{mbox}and{T}_{j,j+1}=\frac{2k_{jy}}{k_{jy}+{\nu }_j{k}_{j+1,y}},$$

$$\begin{array}{rl}{S}_j^t[\varphi ](\mathbf{r})={\int }_{{\Gamma }_{j-1}}{G}_{k_j}(\mathbf{r},{\mathbf{r}}^{\prime })\varphi ({\mathbf{r}}^{\prime })\mathrm{d}{s}_{{\mathbf{r}}^{\prime }},& D_j^t[\varphi ](\mathbf{r})={\int }_{{\Gamma }_{j-1}}\frac{\partial G_{k_j}}{\partial n_{{\mathbf{r}}^{\prime }}}(\mathbf{r},{\mathbf{r}}^{\prime })\varphi ({\mathbf{r}}^{\prime })\mathrm{d}{s}_{{\mathbf{r}}^{\prime }},\\ S_j^b[\varphi ](\mathbf{r})={\int }_{{\Gamma }_j}{G}_{k_j}(\mathbf{r},{\mathbf{r}}^{\prime })\varphi ({\mathbf{r}}^{\prime })\mathrm{d}{s}_{{\mathbf{r}}^{\prime }},& D_j^b[\varphi ](\mathbf{r})={\int }_{{\Gamma }_j}\frac{\partial G_{k_j}}{\partial n_{{\mathbf{r}}^{\prime }}}(\mathbf{r},{\mathbf{r}}^{\prime })\varphi ({\mathbf{r}}^{\prime })\mathrm{d}{s}_{{\mathbf{r}}^{\prime }},\end{array}$$

$$\begin{array}{rl}{K}_j^t[\varphi ](\mathbf{r})={\int }_{{\Gamma }_{j-1}}\frac{\partial G_{k_j}}{\partial n_{\mathbf{r}}}(\mathbf{r},{\mathbf{r}}^{\prime })\varphi ({\mathbf{r}}^{\prime })\mathrm{d}{s}_{{\mathbf{r}}^{\prime }}& N_j^t[\varphi ](\mathbf{r})={\int }_{{\Gamma }_{j-1}}\frac{{\partial }^2{G}_{k_j}}{\partial n_{\mathbf{r}}\partial n_{{\mathbf{r}}^{\prime }}}(\mathbf{r},{\mathbf{r}}^{\prime })\varphi ({\mathbf{r}}^{\prime })\mathrm{d}{s}_{{\mathbf{r}}^{\prime }},\\ K_j^b[\varphi ](\mathbf{r})={\int }_{{\Gamma }_j}\frac{\partial G_{k_j}}{\partial n_{\mathbf{r}}}(\mathbf{r},{\mathbf{r}}^{\prime })\varphi ({\mathbf{r}}^{\prime })\mathrm{d}{s}_{{\mathbf{r}}^{\prime }}& N_j^b[\varphi ](\mathbf{r})={\int }_{{\Gamma }_j}\frac{{\partial }^2{G}_{k_j}}{\partial n_{\mathbf{r}}\partial n_{{\mathbf{r}}^{\prime }}}(\mathbf{r},{\mathbf{r}}^{\prime })\varphi ({\mathbf{r}}^{\prime })\mathrm{d}{s}_{{\mathbf{r}}^{\prime }},\end{array}$$

$${\phi }_j=u_{j+1}\text{mbox}and{\psi }_j=\frac{\partial u_{j+1}}{\partial n}\text{mbox}on{\Gamma }_j.$$

$$\boldsymbol{\phi}_{j}=\left[\varphi_{j},\psi_{j}\right]^{T},\quad\boldsymbol{\phi}^{\mathrm{inc}}=\left[u^{\mathrm{inc}}\big{|}_{\Gamma_{1}},\frac{\partial u^{\mathrm{inc}}}{\partial n}\Big{|}_{\Gamma_{1}}\right]^{T}\quad\mbox{and}\quad\boldsymbol{\phi}^{\parallel}=\left[u_{N}^{\parallel}\big{|}_{\Gamma_{N-1}},\frac{\partial u^{\parallel}_{N}}{\partial n}\Big{|}_{\Gamma_{N-1}}\right]^{T}$$

$$\begin{split}&\operatorname{E}_{j}=\left[\begin{array}[]{cc}1&0\\ 0&\frac{1+\nu_{j}}{2}\end{array}\right],\;\;\;\;\operatorname{T}_{j}=\left[\begin{array}[]{cc}D_{j+1}^{t}-D_{j}^{b}&-S_{j+1}^{t}+\nu_{j}S_{j}^{b}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ N_{j+1}^{t}-N_{j}^{b}&-K_{j+1}^{t}+\nu_{j}K_{j}^{b}\end{array}\right],\ \operatorname{L}_{j}=\left[\begin{array}[]{cc}D_{j}^{t}&-S_{j}^{t}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ N_{j}^{t}&-K_{j}^{t}\end{array}\right]\\ &\mbox{and }\operatorname{R}_{j}=\left[\begin{array}[]{cc}-D_{j+1}^{b}&\nu_{j+1}S_{j+1}^{b}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ -N_{j+1}^{b}&\nu_{j+1}K_{j+1}^{b}\end{array}\right]\mbox{ (all operators evaluated at observation points $\boldsymbol{r}$ on $\Gamma_{j}$).}\end{split}$$

$$\begin{array}[]{lclll}u_{1}(\boldsymbol{r})&=&D_{1}^{b}[\varphi_{1}](\boldsymbol{r})-\nu_{1}S_{1}^{b}[\psi_{1}](\boldsymbol{r})+u^{\mathrm{inc}}(\boldsymbol{r}),&\vskip 3.0pt plus 1.0pt minus 1.0pt\\ u_{j}(\boldsymbol{r})&=&D_{j}^{b}[\varphi_{j}](\boldsymbol{r})-\nu_{j}S_{j}^{b}[\psi_{j}](\boldsymbol{r})\vskip 3.0pt plus 1.0pt minus 1.0pt\\ &&-D_{j}^{t}[\varphi_{j-1}](\boldsymbol{r})+S_{j}^{t}[\psi_{j-1}](\boldsymbol{r}),&j=2,\ldots,N-1,\vskip 3.0pt plus 1.0pt minus 1.0pt\\ u_{N}(\boldsymbol{r})&=&-D_{N}^{t}[\varphi_{N-1}](\boldsymbol{r})+S_{N}^{t}[\psi_{N-1}](\boldsymbol{r})+u^{\parallel}_{N}(\boldsymbol{r}),&\end{array}$$

$$E_1{\mathbf{\varphi }}_1+T_1\left[{\mathbf{\varphi }}_1\right]+R_1\left[{\mathbf{\varphi }}_2\right]={\mathbf{\varphi }}^{\mathrm{inc}}\text{mbox}on{\Gamma }_1.$$

$$E_j{\mathbf{\varphi }}_j+L_j\left[{\mathbf{\varphi }}_{j-1}\right]+T_j\left[{\mathbf{\varphi }}_j\right]+R_j\left[{\mathbf{\varphi }}_{j+1}\right]=0\text{mbox}on{\Gamma }_j,j=2,\dots ,N-2$$

$$E_{N-1}{\mathbf{\varphi }}_{N-1}+L_{N-1}\left[{\mathbf{\varphi }}_{N-2}\right]+T_{N-1}\left[{\mathbf{\varphi }}_{N-1}\right]={\mathbf{\varphi }}^{\parallel }\text{mbox}on{\Gamma }_{N-1}.$$

$$\mathbf{\varphi }=[{\varphi }_1,{\psi }_1,{\phi }_2,{\psi }_2,\cdots ,{\phi }_{N-1},{\psi }_{N-1}{]}^T:\Gamma \to {\mathbb{C}}^{2(N-1)}.$$

$$\boldsymbol{\phi}^{\mathrm{inc}}=\left[u^{\mathrm{inc}}|_{\Gamma_{1}},\frac{\partial u^{\mathrm{inc}}}{\partial n}\Big{|}_{\Gamma_{1}},0,0,\cdots,0,0,u^{\parallel}_{N}|_{\Gamma_{N-1}},\frac{\partial u^{\parallel}_{N}}{\partial n}\Big{|}_{\Gamma_{N-1}}\right]^{T}:\Gamma\to\mathbb{C}^{2(N-1)}.$$

$$\mathcal{E}=\left[\begin{array}[]{cccc}\operatorname{E}_{1}&&&\\ &\operatorname{E}_{2}&&\\ &&\ddots&\\ &&&\operatorname{E}_{N-1}\end{array}\right]\quad\mbox{and}\quad\mathcal{T}_{\Gamma}=\left[\begin{array}[]{ccccc}\operatorname{T}_{1}&\operatorname{R}_{1}&&&\\ \operatorname{L}_{2}&\operatorname{T}_{2}&\operatorname{R}_{2}&&\\ &\operatorname{L}_{3}&\ddots&\ddots\\ &&\ddots&\ddots&\operatorname{R}_{N-2}\\ &&&\operatorname{L}_{N-1}&\operatorname{T}_{N-1}\end{array}\right]$$

$$\mathcal{E}\mathbf{\varphi }+{\mathcal{T}}_{\Gamma }\left[\mathbf{\varphi }\right]={\mathbf{\varphi }}^{\mathrm{inc}}\text{mbox}on\Gamma .$$

$$\mathcal{W}_{A}(\boldsymbol{r}_{1},\boldsymbol{r}_{2},\ldots,\boldsymbol{r}_{N-1})=\left[\begin{array}[]{cccc}w_{A}(x_{1})\operatorname{I}&&&\\ &w_{A}(x_{2})\operatorname{I}&&\\ &&\ddots&\\ &&&w_{A}(x_{N-1})\operatorname{I}\end{array}\right],\quad\boldsymbol{r}_{j}=(x_{j},y_{j})\in\Gamma_{j},$$

$$w_A(x)=\eta (x/A;c,1),$$

$$\eta(t;t_{0},t_{1})=\left\{\begin{array}[]{cll}1,&|t|\leq t_{0},\vskip 6.0pt plus 2.0pt minus 2.0pt\\ \!\!\exp\left(\displaystyle\frac{2\operatorname{e}^{-1/u}}{u-1}\right),&\displaystyle t_{0}<|t|<t_{1},u=\frac{|t|-t_{0}}{t_{1}-t_{0}},\vskip 6.0pt plus 2.0pt minus 2.0pt\\ 0,&|t|>t_{1}.\end{array}\right.$$

$$\mathcal{E}\mathbf{\varphi }+{\mathcal{T}}_{\Gamma }\left[{\mathcal{W}}_A\mathbf{\varphi }\right]={\mathbf{\varphi }}^{\mathrm{inc}}-{\mathcal{T}}_{\Gamma }\left[(\mathcal{I}-{\mathcal{W}}_A)\mathbf{\varphi }\right]\text{mbox}on{\Gamma }_A$$

$${\Gamma }_A=\prod _{j=1}^{N-1}{\Gamma }_{j,A}.$$

$${\tilde{\Gamma }}_A={\Gamma }_A\cap \prod _{j=1}^{N-1}\{(x_j,y_j)\in {\Gamma }_j:w_A(x_j)=1\}$$

$$\mathcal{E}{\mathbf{\varphi }}^w+{\mathcal{T}}_{\Gamma }\left[{\mathcal{W}}_A{\mathbf{\varphi }}^w\right]={\mathbf{\varphi }}^{\mathrm{inc}}-{\mathcal{T}}_P[(\mathcal{I}-{\mathcal{W}}_A){\mathbf{\varphi }}^p]\text{mbox}on{\Gamma }_A$$

$$\mathcal{E}{\mathbf{\varphi }}^w+{\mathcal{T}}_{\Gamma }\left[{\mathcal{W}}_A{\mathbf{\varphi }}^w\right]={\mathbf{\varphi }}^{\mathrm{inc}}+{\mathcal{T}}_P\left[{\mathcal{W}}_A{\mathbf{\varphi }}^p\right]-{\mathcal{T}}_P\left[{\mathbf{\varphi }}^p\right]\text{mbox}on{\Gamma }_A.$$

$$\boldsymbol{\mu}_{j}=\delta_{1,j}\boldsymbol{\phi}^{\mathrm{inc}}+\delta_{N-1,j}\boldsymbol{\phi}^{\parallel}-\left\{\begin{array}[]{clc}\operatorname{E}_{j}\boldsymbol{\phi}^{p}_{j}&\mbox{on}\quad\Gamma_{j}\cap P_{j},\vskip 6.0pt plus 2.0pt minus 2.0pt\\ \left[\begin{array}[]{cc}u^{p}\\ \nabla u^{p}\cdot n\end{array}\right]&\mbox{on}\quad\Gamma_{j}\cap(D_{j}\cup D_{j+1}),\end{array}\right.$$

$$u_j^d=u_j-{\tilde{u}}_j^p\text{mbox}in{\Omega }_j(j=1,\dots ,N-1),$$

$${\tilde{u}}_j^p(x,y)=A_j{e}^{i{k}_{1x}x}\left\{e^{-i{k}_{jy}y}+{\tilde{R}}_{j,j+1}{e}^{i{k}_{jy}(y+2d_j)}\right\}\text{mbox}in{\Omega }_j,1\le j\le N.$$

$$\begin{array}[]{lcll}\tilde{u}^{p}_{1}(\boldsymbol{r})&=&D_{1}^{b}\left[\tilde{\varphi}^{p}_{1}+f_{1}\right](\boldsymbol{r})-\nu_{1}S_{1}^{b}\left[\tilde{\psi}^{p}_{1}+g_{1}\right](\boldsymbol{r})+u^{\mathrm{inc}}(\boldsymbol{r}),&\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \tilde{u}^{p}_{j}(\boldsymbol{r})&=&D_{j}^{b}\left[\tilde{\varphi}^{p}_{j}+f_{j}\right](\boldsymbol{r})-\nu_{j}S_{j}^{b}\left[\tilde{\psi}^{p}_{j}+g_{j}\right](\boldsymbol{r})\vskip 3.0pt plus 1.0pt minus 1.0pt\\ &&-D_{j}^{t}\left[\tilde{\varphi}^{p}_{j-1}\right](\boldsymbol{r})+S_{j}^{t}\left[\tilde{\psi}^{p}_{j-1}\right](\boldsymbol{r}),&j=2,\ldots,N-1,\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \tilde{u}^{p}_{N}(\boldsymbol{r})&=&-D_{N}^{t}\left[\tilde{\varphi}^{p}_{N-1}\right](\boldsymbol{r})+S_{N}^{t}\left[\tilde{\psi}^{p}_{N-1}\right](\boldsymbol{r})+u^{\parallel}_{N}(\boldsymbol{r}),&\end{array}$$

Full text

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 1b—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 $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 Section 7. The unperturbed configuration, which is shown in Figure 1a for reference, consists of $N$ planar layers given by $D_{j}=\mathbb{R}\times(-d_{j},-d_{j-1})$ for $j=1,\ldots,N$. The planar boundary at the interface between the layers $D_{j}$ and $D_{j+1}$ is denoted by $P_{j}=\mathbb{R}\times\{-d_{j}\}$ ($j=1,\ldots,N-1$). The corresponding perturbed layers and their boundaries will be denoted by $\Omega_{j}$, $j=1,\ldots N$ and $\Gamma_{j}$, $j=1,\ldots,N-1,$ respectively; naturally, it is assumed that $\Gamma_{j}\cap\Gamma_{i}=\emptyset$.

Letting $\omega>0$ and $c_{j}>0$ denote the angular frequency and the speed of sound in the layer $\Omega_{j}$, the wavenumber $k_{j}$ in that layer is given by $k_{j}=\omega/c_{j}$. Assuming e.g. an incident plane wave of the form $u^{\mathrm{inc}}(\boldsymbol{r})=\operatorname{e}^{ik_{1}(x\cos\alpha+y\sin\alpha)}$ (where $\boldsymbol{r}=(x,y)$ and where $-\pi<\alpha<0$ denotes the angle of incidence measured with respect to the horizontal) and letting $u$ denote the acoustic pressure, the restrictions $u_{j}=u|_{\Omega_{j}}$ of the total field $u$ to the domains $\Omega_{j}$ ($j=1,\dots,N$) satisfy the homogeneous Helmholtz equation

[TABLE]

together with the transmission conditions

[TABLE]

where $\nu_{j}=\varrho_{j}/\varrho_{j+1}$ where $\varrho_{j}$ denotes the fluid density in $\Omega_{j}$. For definiteness, here and throughout this paper the unit normal $n=n(\boldsymbol{r})$ for $\boldsymbol{r}\in\Gamma_{j}$ is assumed to point into $\Omega_{j}$.

As is well known, a closed form expression exists [8, 26] for the total field $u^{p}$ throughout space ($u^{p}=u^{p}_{j}$ in $D_{j}$, $j=1,\ldots,N$), that results as a plane wave $u^{\mathrm{inc}}$ impinges on the planar layer medium $D=\bigcup_{j=1}^{N}D_{j}$. In detail, letting $k_{1x}=k_{1}\cos\alpha$ and $k_{jy}=\sqrt{k_{j}^{2}-k_{1x}^{2}}$, $j=1,\ldots,N$ (where the complex square-root is defined in such a way that $\mathrm{Im}\,k_{jy}\geq 0$, which, noting that $\mathrm{Im}\,k_{j}^{2}\geq 0$, requires $\mathrm{Re}\,k_{jy}\geq 0$ as well), the planar-medium solution $u_{j}^{p}$ in $D_{j}$ is given by

[TABLE]

in terms of certain generalized reflection coefficients $\widetilde{R}_{j,j+1}$ and amplitudes $A_{j}$. The amplitudes and the generalized reflection coefficients can be obtained recursively by means of the relations

[TABLE]

and

[TABLE]

in terms of the reflection and transmission coefficients

[TABLE]

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 $\Gamma_{j}$, $j=1,\ldots,N-1$. As in the contribution [4] we utilize the single- and double-layer potentials

[TABLE]

which are defined for $\boldsymbol{r}\in\mathbb{R}^{2}$ and are expressed in terms of improper integrals whose convergence is conditioned upon the oscillatory behavior of the integrand. Here we have called $G_{k_{j}}(\boldsymbol{r},\boldsymbol{r}^{\prime})=\frac{i}{4}H_{0}^{(1)}(k_{j}|\boldsymbol{r}-\boldsymbol{r}^{\prime}|)$ the free-space Green function for the Helmholtz equation with wavenumber $k_{j}$. Additionally, we define the integral operators

[TABLE]

where the evaluation point $\boldsymbol{r}$ belongs to either $\Gamma_{j}$ or $\Gamma_{j-1}$.

In order to formulate an integral equation for the unknown interface values we define the unknown density functions $\varphi_{j}:\Gamma_{j}\to\mathbb{C}$ and $\psi_{j}:\Gamma_{j}\to\mathbb{C}$ ($j=1,\ldots,N-1$) by

[TABLE]

Additionally we define the vector density functions

[TABLE]

where $u_{N}^{\parallel}$ is defined in (56), and the matrix operators

[TABLE]

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

[TABLE]

in terms of the interface values (12). Therefore, evaluating $u_{1}+u_{2}$ and $\partial u_{1}/\partial n+\partial u_{2}/\partial n$ on $\Gamma_{1}$ from the boundary values on $\Gamma_{1}$ of the expressions in (15) and their normal derivatives, and using the notations (13) and (14), we obtain the $j=1$ interface equation

[TABLE]

(Note that, of course, the calculations leading to equations (16c) 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 (16c) are expressed in terms of a single column vector function $\boldsymbol{\phi}$ (defined on the Cartesian product $\Gamma=\prod_{j=1}^{N-1}\Gamma_{j}$ of the curves $\Gamma_{j}$) whose $j$-entry equals the density pair $\boldsymbol{\phi}_{j}=[\varphi_{j},\psi_{j}]^{T}:\Gamma_{j}\to\mathbb{C}^{2}$ for $j=1,\dots,N-1$. We may thus write

[TABLE]

Similarly we define

[TABLE]

With a slight notational abuse we will write $\boldsymbol{\phi}=[\boldsymbol{\phi}_{1},\boldsymbol{\phi}_{2},\ldots,\boldsymbol{\phi}_{N-1}]^{T}=[\phi_{1},\psi_{1},\varphi_{2},\psi_{2},\cdots,\varphi_{N-1},\psi_{N-1}]^{T}$. More generally, given arbitrary vectors $\boldsymbol{\mu}_{j}=[\alpha_{j},\beta_{j}]^{T}:\Gamma_{j}\to\mathbb{C}^{2}$ for $j=1,\dots,N-1$ we will use the “block-vector” notation $\boldsymbol{\mu}=[\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2},\ldots,\boldsymbol{\mu}_{N-1}]^{T}=[\alpha_{1},\beta_{1},\alpha_{2},\beta_{2},\cdots,\alpha_{N-1},\beta_{N-1}]^{T}:\Gamma\to\mathbb{C}^{2(N-1)}$.

Using the operators

[TABLE]

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

[TABLE]

4 Windowed integral equations

Following [4], in this section we introduce rapidly-convergent windowed versions of the integral formulation (18). In order to do so we utilize the $(N-1)\times(N-1)$ block-diagonal matrix-valued window function $\mathcal{W}_{A}:\prod_{j=1}^{N-1}\Gamma_{j}\mapsto\mathbb{R}^{2(N-1)\times 2(N-1)}$ given by

[TABLE]

in terms of the two-by-two identity matrix $\operatorname{I}$ and the smooth window function

[TABLE]

where $0<c<1$ and where

[TABLE]

Clearly $\eta$ 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 $\boldsymbol{r}=(x,y)\in\mathbb{R}^{2}$ equals the set $[-A,A]\times\mathbb{R}$. Note that the parameter $c$, which controls the steepness of the rise of the window function $w_{A}$, is not displayed as part of the notation $w_{A}$.

(While different values $A_{j}$ of the window-size $A$ , $j=1,\dots,N-1$, could in principle be used for the various layer interfaces and corresponding block entries in (19)—possibly utilizing smaller (resp. larger) values $A_{j}$ 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 $A$ is used for all the interfaces.)

In order to produce a windowed version of equation (18) we proceed in two stages. At first the integrand is multiplied by the window matrix $\mathcal{W}_{A}$ and the equation is restricted to the windowed region $\Gamma_{A}=\{(\boldsymbol{r}_{1},\dots,\boldsymbol{r}_{N-1})\in\Gamma:\boldsymbol{r}_{j}=(x_{j},y_{j})\mbox{ and }w_{A}(x_{j})\neq 0\mbox{ for all }j\}\subset\mathbb{R}^{2(N-1)}$—so that, moving the remainder of the windowed integral operator to the right-hand side and letting $\mathcal{I}$ denote the identity matrix of dimension $2(N-1)\times 2(N-1)$, the exact relation

[TABLE]

results. Note that defining $\Gamma_{j,A}=\Gamma_{j}\cap\{w_{A}\neq 0\}=\Gamma_{j}\cap\left\{[-A,A]\times\mathbb{R}\right\}$ we have

[TABLE]

A successful implementation of the WGF idea requires use of an accurate substitute for the quantity $\mathcal{T}_{\Gamma}\left[(\mathcal{I}-\mathcal{W}_{A})\boldsymbol{\phi}\right]$ throughout $\Gamma_{A}$, which does not depend on knowledge of the unknown density $\boldsymbol{\phi}$ (cf. [4]). In order to obtain such an approximation we introduce an operator $\mathcal{T}_{P}$ which is defined just like $\mathcal{T}_{\Gamma}$ in (17) but in terms of potentials (10) and operators (11) given by integrals on the flat interfaces $P_{j}$ depicted in Figure 1a. Since $(\mathcal{I}-\mathcal{W}_{A})$ vanishes wherever $\Gamma$ differs from $P=\prod_{j=1}^{N-1}P_{j}$, we clearly have $\mathcal{T}_{\Gamma}\left[(\mathcal{I}-\mathcal{W}_{A})\boldsymbol{\phi}\right]=\mathcal{T}_{P}\left[(\mathcal{I}-\mathcal{W}_{A})\boldsymbol{\phi}\right]$. Additionally, we consider the aforementioned scalar densities $\varphi^{p}_{j}=u^{p}_{j+1}|_{\Gamma_{j}}$ and $\psi^{p}_{j}=\partial u^{p}_{j+1}/\partial n|_{\Gamma_{j}}$ on $P_{j}$ ($j=1,\ldots,N-1$) that are associated with the planarly layered medium $P$. As shown in [4], letting $[\boldsymbol{\phi}^{p}]_{j}=[\varphi^{p}_{j},\psi^{p}_{j}]^{T}$ ($j=1,\ldots,N-1$), substitution of $\boldsymbol{\phi}$ by $\boldsymbol{\phi}^{p}$ on the right-hand-side of (25) results in errors that decay super-algebraically fast as $A\to\infty$ within the subset

[TABLE]

of $\Gamma_{A}$ wherein the window function $w_{A}$ equals one. Indeed, even though $\boldsymbol{\phi}$ may differ significantly from $\boldsymbol{\phi}^{p}$, 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

[TABLE]

which we re-express in the form

[TABLE]

As shown in what follows, the right-hand term $\mathcal{T}_{P}\left[\boldsymbol{\phi}^{p}\right]$ in (26) can be expressed in closed form, and thus, using numerical integration over the bounded domain $\Gamma_{A}$ to produce the term $\mathcal{T}_{P}\left[\mathcal{W}_{A}\boldsymbol{\phi}^{p}\right]$, the complete right-hand side can be efficiently evaluated for any given $\boldsymbol{x}\in\Gamma_{A}$.

A closed-form expression for $\boldsymbol{\mu}=\mathcal{T}_{P}\left[\boldsymbol{\phi}^{p}\right]$ (cf. Remark 3.1) can be obtained via an application of Green’s formula: using (52) with $C=P_{j}$ ($j=1,\dots,N-1$), equations (57) yield the desired relations:

[TABLE]

for $j=1,\ldots,N-1,$ where $\delta_{i,j}$ denotes the Kronecker delta symbol.

As demonstrated in Section 7 through a variety of numerical examples, the vector density function $\boldsymbol{\phi}^{w}$, which is the solution of the windowed integral equation (26), converges super-algebraically fast to the exact solution $\boldsymbol{\phi}$ of (18) within $\Gamma^{1}_{A}$ as the window size $A>0$ increases. This observation can be justified via arguments analogous to those presented in [4].

Remark 4.1**.**

The difference $N^{t}_{j+1}-N^{b}_{j}$ of hypersingular operators that appears in the definition of the diagonal blocks $\operatorname{T}_{j}$ of $\mathcal{T}_{\Gamma}$ 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 $u^{w}$ of the solution $u$ of (1)–(2) near the localized defects. In order to obtain this approximation we consider the “defect” field

[TABLE]

given by the difference between the total field $u_{j}$ and the planar-structure total field

[TABLE]

Note that $\tilde{u}_{j}^{p}$ is given in $\Omega_{j}$ by the expressions on the right-hand side of equation (3).

Subtracting the integral representation

[TABLE]

—which follows as equation (57) is applied to $\tilde{u}^{p}_{j}$—from the integral representation (15) we obtain the exact integral relations

[TABLE]

for the defect fields. These relations can be used to evaluate the defect fields $u^{d}_{j}$ in terms of the solution $\boldsymbol{\boldsymbol{\phi}}_{j}=[\varphi_{j},\psi_{j}]^{T}$ of the integral equation (18) together with the planar-structure total fields

[TABLE]

and the jumps

[TABLE]

Note that, importantly, for each $j$, $1\leq j\leq N$, the functions $f_{j}$ and $g_{j}$ vanish outside the $j$-th portion $\Gamma_{j}\setminus\Pi_{j}$ of the boundary of the localized defects.

Relying on the WGF solutions $\boldsymbol{\phi}^{w}$ of equation (26) and applying a windowing procedure similar to the one used in the previous section, a highly-accurate approximation to the defect near-fields (31) results. In detail, substitution of $\varphi_{j}$ by $w_{A}\varphi_{j}^{w}+(1-w_{A})\tilde{\varphi}^{p}_{j}$ and $\psi_{j}$ by $w_{A}\psi_{j}^{w}+(1-w_{A})\tilde{\psi}^{p}_{j}$ in (31) yields the approximate expressions

[TABLE]

for the defect field $u^{d}_{j}$. The desired approximation $u^{w}_{j}$ for the total field $u_{j}$ then follows from (28):

[TABLE]

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

[TABLE]

containing the localized defects—with uniformly small errors, as $A\to\infty$, within every bounded subset of $\widetilde{\Omega}_{A}$. 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 (34)–(35) only provide uniformly accurate approximations within bounded subsets of $\widetilde{\Omega}_{A}$. But, once accurate defect fields $u^{d,w}_{j}$ ($j=1,\ldots,N$) have been obtained within $\widetilde{\Omega}_{A}$, correspondingly accurate far-field values for the solution $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 $\widetilde{\Omega}_{A}$. In detail, defining the defect field $u^{d}=u^{d}(\boldsymbol{r})$ to equal $u^{d}_{j}(\boldsymbol{r})$ for $\boldsymbol{r}\in\Omega_{j}$ ($j=1,\ldots,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 [24, Lemma 4.2.6]

[TABLE]

which is valid for $\boldsymbol{r}$ everywhere outside $S$. Note that the necessary values of $u^{d}_{j}$ and their normal derivatives on $S$ can be computed by means of (35)—since, by construction, $S$ lies inside the region where (35) provides an accurate approximation of the field $u^{d}_{j}$.

The far-field approximation $u^{f}$ of the defect field $u^{d}$ as $\boldsymbol{r}\to\infty$ in any direction is then obtained by replacing the layer Green function $H$ and its normal derivative $\partial H/\partial n_{\boldsymbol{r}^{\prime}}$ in (37) by the respective first-order $|\boldsymbol{r}|\to\infty$ asymptotic expansions $H^{f}$ and $\partial H^{f}/\partial n_{\boldsymbol{r}^{\prime}}$—which can be obtained for the $N$-layer case (as illustrated in [2, 8, 24, 4, 3] 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 $\partial H/\partial n_{\boldsymbol{r}^{\prime}}$ coincides with $\partial H^{f}/\partial n_{\boldsymbol{r}^{\prime}}$ can be verified by direct inspection of these two quantities.) The far field $u^{f}$ is thus given by

[TABLE]

It is important to note that, unlike the layer Green function $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. In order to evaluate $H^{f}$ we first note that, for a source point $\boldsymbol{r}^{\prime}=(x^{\prime},y^{\prime})\in D_{j}$ and a target point $\boldsymbol{r}=(r\cos\theta,r\sin\theta)\in D_{1}$ ($\theta\in[0,\pi]$), the layer Green function $H$ is given by the contour integral [2, 8, 24, 3]

[TABLE]

Here, letting $\gamma_{j}(\xi)=\sqrt{\xi^{2}-k_{j}^{2}}$, $j=1,2,3$, we have set

[TABLE]

The determination of physically admissible branches of the functions $\gamma_{j}(\xi)=\sqrt{\xi^{2}-k_{j}^{2}}=\sqrt{\xi-k_{j}}\sqrt{\xi+k_{j}}$ requires adequate selection of branch cuts. Relevant branches, which must be selected to insure the Green function satisfies outgoing radiation condition for the layered structure, are given by $-3\pi/2\leq\arg(\xi-k_{j})<\pi/2$ for $\sqrt{\xi-k_{j}}$ and $-\pi/2\leq\arg(\xi+k_{j})<3\pi/2$ for $\sqrt{\xi+k_{j}}$. The branch cut stemming from the point $\xi=k_{1}=k_{3}$ is in fact the only branch cut in the domain of definition of the functions $p_{j}(\xi,\boldsymbol{r}^{\prime})/f(\xi)$ ($j=1,2,3$), as it can be shown that these are even functions of $\gamma_{2}$ [8]. The branch cuts and Sommerfeld contour ${\mathit{SC}}$ utilized in (39) 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 [2]. Analysis of the phase function $\phi$ (40a) readily shows that there is only one saddle point on the real axis at $\xi_{0}=k_{1}\cos\theta$ and that the path of steepest descent $\mathit{SD}$ that passes through that point, which is given by the expression $\mathrm{Im}\,\phi(\xi)=k_{1}$, also intersects the real axis at $\xi=k_{1}/\cos\theta$. Furthermore, from the definition of the function $\gamma_{1}$ it can shown that

[TABLE]

for $\xi$ on $\mathit{SD}$. This information suffices to sketch the paths of steepest descent that are displayed in Figure 3.

In order to produce asymptotic expansions of the integrals (39) we then proceed to deform the Sommerfeld contour $\mathit{SC}$ to the steepest descent contour $\mathit{SD}$ (Figure 3). Considering the saddle point at $\xi_{0}$ and taking into account the poles of the integrand $p_{j}(\xi,\boldsymbol{r}^{\prime})/q(\xi)$ at the points $\xi_{p}$ which are enclosed by the curves $\mathit{SD}$ and $\mathit{SC}$, we obtain the following expression for the far-field form of the Green function for the two-dimensional slab:

[TABLE]

as $|\boldsymbol{r}|\to\infty$. Note that for $\cos\theta>0$ (resp. $\cos\theta<0$) only the real poles contained in the set $I=(0,k_{1}\cos\theta)\cup(k_{1}/\cos\theta,\infty)$ (resp. $I=(-\infty,k_{1}/\cos\theta)\cup(k_{1}\cos\theta,0)$) produce contributions which do not decay exponentially.

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 expressions.

7 Numerical examples

This section presents a set of two- and three-dimensional numerical examples that demonstrate the character of the proposed multi-layer WGF methodology. For the sake of definiteness a window function $w_{A}$ (20) 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 [23, 24] (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 (26) by means of the Nyström method described in [10, Section 3.5]. The three-dimensional solutions, in turn, were obtained by means of the algorithm presented in [5].

Our first example concerns the structure depicted in the left portion of Figure 4, in which semi-circular defects of radii $a=1$ are placed at the planar interfaces $P_{1}=\mathbb{R}\times\{0\}$ and $P_{2}=\mathbb{R}\times\{-3/2\}$ of a three-layer medium with wavenumbers $k_{1}=10$, $k_{2}=20$ and $k_{3}=30$. The right portion of Figure 4 displays the maximum relative errors (in log-log scale) in the total field produced by the WGF method on the surface of the semi-circular defects (the curves marked in red in Figure 4 left) for various windows sizes $A>0$ and incidences $\alpha$. 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\lambda$ 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.

In order to compare the computational cost of the LGF method [23] 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 semi-circular 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^{\mathrm{inc}}$ with $\alpha=-\pi/6$ illuminates the structure. Five sets of wavenumbers given by $k_{1}=\kappa,$ $k_{2}=2\kappa$ and $k_{3}=3\kappa$ with $\kappa=2^{j}$, $j=1,\ldots,5$ are considered. The resulting problems of scattering are then solved by employing a Nyström discretization of the WGF equations (26), and a numerical version of (35) 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 [23] (see also [24]). 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. In order 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 $\kappa$. The maximum of the absolute value of the difference between the LGF and WGF solutions (using $A=8\lambda$) on the surface of the defect is no larger than $10^{-4}$ in all the examples considered. Remarkably, in the $\kappa=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 (1st column) and by the LGF method (2nd column) respectively for $\kappa=32$. The fields are evaluated in the rectangular region $[-3,3]\times[-7/2,2]$ at an uniform grid of $280\times 200$ points. Note that, as it follows from consideration of the figure captions, the WGF near field evaluation procedure is up to 1200 times faster than the corresponding LGF near field evaluation procedure—in spite of the fact that a (larger) window size $A=16\lambda$ had to be used to produce accurate near fields throughout the plotted region.

Figure 6, in turn, compares the far-field patterns

[TABLE]

produced by the WGF and LGF algorithms for a semi-circular cavity in a three-layer medium with wavenumbers $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 (42), where $u^{f}$ is given by (38) with WGF defect fields $u^{d}=u^{d,w}_{j}$ in $\Omega_{j}$, $j=1,\ldots,N$ (equation (35)). 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 [23]. We have verified that, as expected, the accuracy of the WGF far-field patterns is comparable to the accuracy of the corresponding defect fields $u^{d,w}_{j}$ within the region $\tilde{\Omega}_{A}$.

Figure 7 displays near fields resulting from the WGF method, with window size $A=12\lambda$, for a structure consisting of nested circular surface defects in a nine-layer medium with planar interfaces $P_{j}=\mathbb{R}\times\{(j-1)/5\}$, $j=1,\ldots,8$. The corresponding wavenumbers are $k_{2j-1}=15$ for $j=1,\ldots,5$ and $k_{2j}=30$ for $j=1,\ldots 4$. The structure is illuminated by plane-waves with two different incidence angles. A 112-second 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 right.

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

Acknowledgements

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 Appendix: 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(\boldsymbol{r})=v_{j}(\boldsymbol{r})$ for $\boldsymbol{r}\in\Omega_{j}$, $j=1,\ldots,N$, where, letting $u^{d}_{j}$ and $\tilde{u}^{p}_{j}$ be defined in (28) and (29), respectively, we have either

[TABLE]

The presentation is restricted to two-dimensional configurations. A related (modified) representation, which can similarly be utilized for all purposes necessary in this paper, can be obtained analogously—albeit with certain additional considerations, as detailed in [15] in the three-dimensional sound-hard case; cf. also [14]. For simplicity, the presentation is further restricted to three-layer structures, but the extension to $N$-layer structures is straightforward.

Our derivations utilize three local polar-coordinate systems, each one of which is associated with one of the layers $\Omega_{j}$. These coordinate systems are centered at $(0,-d_{1})$, $(0,-(d_{1}+d_{2})/2)$ and $(0,-d_{2})$ and, thus, the radial variables are given by

[TABLE]

in terms of the global Cartesian coordinates $x$ and $y$, as illustrated in Figure 9.

Additionally, some of the subsequent derivations utilize the decomposition

[TABLE]

where letting $k_{1x}=k_{1}\cos\alpha$ and $k_{jy}=\sqrt{k_{j}^{2}-k^{2}_{1x}}$, the up-going and down-going plane-waves $u_{j}^{\uparrow}$ and $u_{j}^{\downarrow}$ are given by

[TABLE]

respectively. Here the constants $p_{j}=\operatorname{e}^{2ik_{jy}d_{j}}A_{j}\widetilde{R}_{j,j+1}$ and $q_{j}=A_{j}$ are expressed in terms of the amplitudes $A_{j}$ and the generalized reflection coefficients $\widetilde{R}_{j,j+1}$ defined in (6) and (9). Note that $u_{1}^{\downarrow}=u^{\mathrm{inc}}$ and $u_{3}^{\uparrow}=0$. The defect field $v^{d}_{j}$, on the other hand, is given by [16, 20]

[TABLE]

in terms of radiative and guided wave fields $v_{j}^{\mathrm{rad}}$ and $v_{j}^{\mathrm{gui}}$ which, letting $\beta_{2}=1$ and $\beta_{j}=2/3$ for $j=1,3$ (see (44)), verify

[TABLE]

and

[TABLE]

Here $v^{m}_{j}$ denote the guided modes

[TABLE]

which are expressed in terms of the so-called propagation constants $\xi_{j}^{m}>0$, and $\gamma_{j}^{m}=\sqrt{(\xi_{j}^{m})^{2}-k_{j}^{2}}$, $m=1,\ldots,M_{j}$. The propagation constants $\xi_{j}^{m}$ equal the real poles (sometimes called surface wave poles [9, 7]) of the corresponding three-layer Green function in spectral form. The condition for the existence of the propagative modes in the inner layer $\Omega_{2}$ is $k_{1}<\xi_{2}^{m}<k_{2}$. For the outer layer $\Omega_{1}$ (resp. $\Omega_{3}$), on the other hand, we have $\xi_{1}^{m}=\xi_{2}^{m}$ (resp. $\xi_{3}^{m}=\xi_{2}^{m}$) and the guided-mode condition is $\xi_{1}^{m}>k_{1}$ (resp. $\xi_{3}^{m}>k_{3}$). Thus $v_{1}^{m}$ (resp. $v_{3}^{m}$) corresponds to a surface wave that travels along the interface $\Gamma_{1}$ (resp. $\Gamma_{2}$) and decays exponentially fast towards the interior of $\Omega_{1}$ (resp. $\Omega_{3}$).

We are now in a position to derive the desired integral representation for the total fields $v_{j}$ in (43). Our derivations consider at first the bounded domains

[TABLE]

where $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 $\Gamma^{R}_{2,l}$, $\Gamma^{R}_{2,r}$, $S^{R}_{1}$ and $S^{R}_{3}$ and corresponding normals $n$, as depicted in Figure 9. In order to facilitate repeated use of Green’s third identity in our derivations we follow [14] and, letting $G_{k_{j}}(\boldsymbol{r},\boldsymbol{r}^{\prime})=\frac{i}{4}H_{0}^{(1)}(k_{j}|\boldsymbol{r}-\boldsymbol{r}^{\prime}|)$ we define, for a given curve $C$,

[TABLE]

In what follows, finally, we make frequent use of the $|\boldsymbol{r}^{\prime}|\rightarrow\infty$ asymptotic relations

[TABLE]

that follow directly from the corresponding asymptotic expressions for the Hankel function [17] together with easily verified identity $|\boldsymbol{r}-\boldsymbol{r}^{\prime}|=|\boldsymbol{r}^{\prime}|-\boldsymbol{r}^{\prime}\cdot\boldsymbol{r}/|\boldsymbol{r}^{\prime}|+O\left(|\boldsymbol{r}^{\prime}|^{-1}\right)$.

With reference to Figure 9, and in view of Green’s third identity applied to $\Omega_{1}^{R}$ and its boundary $\partial\Omega_{1}^{R}=\Gamma_{1}^{R}\cup S_{1}^{R}$ we obtain the bounded-domain integral representation

[TABLE]

result. In order to complete our calculations it suffices to evaluate the limiting values as $R\to\infty$ for the various integral quantities in (54c) and for each one of the functions $v_{j}$ in (43). Since, in view of equations (45), (46) and (47), these functions can be expressed as linear combinations of $u_{j}^{\uparrow}$, $u_{j}^{\downarrow}$, $v_{j}^{\mathrm{gui}}$ and $v_{j}^{\mathrm{rad}}$ in what follows we obtain the corresponding limiting values for each one of these functions. The desired representation formulae (55) as well as their $N$-layer versions (57) then follow directly from the limiting expressions thus found.

Case $j=2$ for $v_{2}=u_{2}^{\uparrow}$, $v_{2}=u_{2}^{\downarrow}$, $v_{2}=v_{2}^{\mathrm{gui}}$ and $v_{2}=v_{2}^{\mathrm{rad}}$. In view of the decay of the integral kernels (53) and the fact that the total field $v_{2}$ remains bounded throughout $\Omega_{2}$ (as it follows from equation (49)), we conclude that the term $I_{2}$ involving the integral over $\Gamma^{R}_{2,r}\cup\Gamma_{2,\ell}^{R}$ tends to zero as $R\to\infty$.

Case $j=1,3$ for $v_{j}=v_{j}^{\mathrm{rad}}$. In order to estimate the terms $I_{j}$ that involve integrals over the semi-circular curves $S_{1}^{R}$ and $S_{3}^{R}$, in turn, we note that for $\boldsymbol{r}^{\prime}\in S_{j}^{R}$ with $j=1$ and $j=3$ we have $|\boldsymbol{r}^{\prime}|=R+O(1)$ and $\widehat{\boldsymbol{r}^{\prime}}=(\cos\theta_{j},\sin\theta_{j})+O(R^{-1})$ as $R\to\infty$—where the the angles $\theta_{j}$ are as shown in Figure 9. Since $v^{\mathrm{rad}}_{j}$ in (47) for $j=1,3$ satisfies the Sommerfeld radiation condition (48), utilizing standard arguments [11] it can be shown that

[TABLE]

Case $j=1,3$ for $v_{j}=v_{j}^{\mathrm{gui}}$. We now consider the quantity $I_{1}\left[v^{\mathrm{gui}}_{1};S_{1}^{R}\right]$, which, according to (49), is given by a linear combination of terms of the form $I_{1}\left[v_{1}^{m};S_{1}^{R}\right]$ ($m=1,\ldots,M_{j}$), where $v_{1}^{m}(\boldsymbol{r})=\operatorname{e}^{-\gamma_{1}^{m}|y|+i\xi_{1}^{m}|x|}$ with $\gamma_{1}^{m}>0,$. From (53) and the fact that $v_{1}^{m}(\boldsymbol{r}^{\prime})=\operatorname{e}^{-\gamma_{1}^{m}R\sin\theta_{1}+i\xi_{1}^{m}R|\cos\theta_{1}|}$ for $\boldsymbol{r}^{\prime}=R(\cos\theta_{1},\sin\theta_{1})\in S_{1}^{R}$, $\theta_{1}\in[0,\pi]$, we obtain

[TABLE]

[TABLE]

as $R\to\infty$. Therefore

[TABLE]

The integral in the expression on the right-hand-side can be bounded utilizing the inequality $\sin\theta\geq 2\theta/\pi$, $\theta\in[0,\pi]$, and we thus conclude that $I_{1}\left[v_{1}^{m};S_{1}^{R}\right]=O\left(R^{-1/2}\right)$ as $R\to\infty$. Similarly, it can be shown that $I_{3}\left[v_{3}^{m};S_{3}^{R}\right]=O\left(R^{-1/2}\right)$, and consequently, in view of (49), we conclude that $I_{j}\left[v_{j}^{\mathrm{gui}};S_{j}^{R}\right]=O\left(R^{-1/3}\right)$, $j=1,3,$ as $R\to\infty$.

Case $j=1$ for $v_{1}=u_{1}^{\uparrow}$ and $v_{1}=u_{1}^{\downarrow}$. We consider the term $I_{1}\left[u_{1}^{\uparrow};S_{1}^{R}\right]$ first, where $u^{\uparrow}_{1}(\boldsymbol{r})=\operatorname{e}^{ik_{1x}x+ik_{1y}y}$, or, in polar coordinates, $u^{\uparrow}_{1}(\boldsymbol{r}^{\prime})=\operatorname{e}^{ik_{1}R\cos(\theta_{1}+\alpha)}$ for $\boldsymbol{r}^{\prime}\in S_{1}^{R}$. Then, integration by parts yields

[TABLE]

as $R\to\infty$. In order to deal with the term $I_{1}\left[u_{1}^{\downarrow};S^{R}_{1}\right]$ with $u^{\downarrow}_{1}(\boldsymbol{r})=\operatorname{e}^{ik_{1x}x-ik_{1y}y}$ we proceed similarly. Using the polar form $u^{\downarrow}_{1}(\boldsymbol{r}^{\prime})=\operatorname{e}^{ik_{1}R\cos(\theta_{1}-\alpha)}$ for $\boldsymbol{r}^{\prime}\in S_{1}^{R}$ of the down-going wave we obtain

[TABLE]

Note that since $\alpha\in(-\pi,0)$ we have $0<\theta_{1}-\alpha<2\pi$. Thus, there is only one point of stationary phase within the domain of integration at $\theta_{1}=\alpha+\pi$. A straightforward application of the method of stationary phase [2] then yields

[TABLE]

(Notice that integrating by parts yields that the limit points of the integral give rise to contributions that decay as $R^{-1}$.)

Case $j=3$ for $v_{j}=u_{j}^{\downarrow}$. This case concerns the term $I_{3}\left[u^{\downarrow}_{3};S_{3}^{R}\right]$ with $u^{\downarrow}_{3}(\boldsymbol{r}^{\prime})=\operatorname{e}^{ik_{3x}x-ik_{3y}y}$, where $k_{3x}=k_{1}\cos\alpha$ and $k_{3y}=\sqrt{k_{3}^{2}-k_{3x}^{2}}$. We distinguish three possible cases, namely: (a) $k_{3}<k_{1}|\cos\alpha|$, (b) $k_{3}=k_{1}|\cos\alpha|$ ($k_{3}=-k_{1}\cos\alpha$ for $\alpha\in(-\pi,-\pi/2]$ or $k_{3}=k_{1}\cos\alpha$ for $\alpha\in(-\pi/2,0)$), and (c) $k_{3}>k_{1}|\cos\alpha|$. Since in case (a) we have $k_{3y}=i\sqrt{k_{1}^{2}\cos^{2}\alpha-k_{3}^{2}}$, a calculation completely analogous to the one carried in the estimation of the term $I_{1}\left[u_{1}^{m};S_{1}^{R}\right]$ shows that $I_{3}\left[u^{\downarrow}_{3};S_{3}^{R}\right]=O\left(R^{-1/2}\right)$. In case (b), in turn, we first consider $\alpha\in(-\pi/2,0)$. Under this assumption $u_{3}^{\downarrow}(\boldsymbol{r}^{\prime})=\operatorname{e}^{ik_{3}R\cos\theta_{3}}$ for $\boldsymbol{r}^{\prime}\in S^{R}_{3}$, and consequently

[TABLE]

Splitting the integration domain and using the identity $\cos\theta-1=-\sin^{2}\theta/(1+\cos\theta)$ we obtain

[TABLE]

Integration by parts yields that the first integral above amounts to a quantity of order $O\left(R^{-1/2}\right)$. The stationary point at $\theta=-\pi$ in the second integral, on the other hand, leads to

[TABLE]

Similarly, in the case $k_{3}=-k_{1}\cos\alpha$ for $\alpha\in(-\pi,\pi/2]$ it can be shown that $I_{3}\left(u_{3}^{\downarrow};S_{3}^{R}\right)=-\operatorname{e}^{-ik_{3}x}/2+O\left(R^{-1/2}\right)$. In the case (c), finally, $u_{3}^{\downarrow}(\boldsymbol{r}^{\prime})=a\operatorname{e}^{ik_{3}R\cos(\theta_{3}-\alpha^{\prime})}$, $\boldsymbol{r}^{\prime}\in S_{3}^{R}$, where $a=\operatorname{e}^{-ik_{3}d_{2}\sin\alpha^{\prime}}$ and where the angle $\alpha^{\prime}\in(-\pi,0)$ is determined by the Snell’s law $k_{3}\cos\alpha^{\prime}=k_{1}\cos\alpha$. Thus, once again, integration by parts yields

[TABLE]

Therefore, taking the limit as $R\to\infty$ in (54c) we obtain

[TABLE]

where $u^{\parallel}_{3}$ in (55i) is given by

[TABLE]

with $N=3$, and where $\delta=0$ if $v_{j}=v_{j}^{d}$, and $\delta=1$ if $v_{j}=\tilde{u}_{j}^{p}+v_{j}^{d}$ or $v_{j}=\tilde{u}_{j}^{p}$.

Remark A.1**.**

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

[TABLE]