Adapting free-space fast multipole method for layered media Green's function: algorithm and analysis

TL;DR

This paper develops an adapted fast multipole method for layered media Green's functions, enabling efficient convolution computations by leveraging translated matrix forms and asymptotic analysis, with validated convergence and quadrature rules.

Contribution

It introduces a novel algorithm that adapts the free-space fast multipole method to layered media, using translated matrix forms and alternative integral representations.

Findings

Efficient convolution computation for layered media Green's functions.

Asymptotic decay estimates for multipole and local expansions.

Validated convergence and quadrature rules for the proposed method.

Abstract

In this paper, we present a numerical algorithm for the accurate and efficient computation of the convolution of the frequency domain layered media Green's function with a given density function. Instead of compressing the convolution matrix directly as in the classical fast multipole method, fast direct solvers, and fast H-matrix algorithms, the new algorithm considers a translated form of the original matrix so that many existing building blocks from the highly optimized free-space fast multipole method can be easily adapted to the Sommerfeld integral representations of the layered media Green's function. An asymptotic analysis is performed on the Sommerfeld integrals for large orders to provide an estimate of the decay rate in the new "multipole" and "local" expansions. In order to avoid the highly oscillatory integrand in the original Sommerfeld integral representations when the source and target are close to each other, or when they are both close to the interface in the scattered field, mathematically equivalent alternative direction integral representations are introduced. The convergence of the multipole and local expansions and formulas and quadrature rules for the original and alternative direction integral representations are numerically validated.

Equations

$$g_2^t({\mathbf{x}}_1,{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{\sqrt{{\lambda }^2-k_2^2}{y}_1}{e}^{i\lambda x_1}{e}^{-\sqrt{{\lambda }^2-k_1^2}{y}_0}{e}^{-i\lambda x_0}\frac{1}{4\pi \sqrt{{\lambda }^2-k_2^2}}{\sigma }_2^t(\lambda )d\lambda$$

$${\sigma }_2^t(\lambda )=\frac{e^{d\sqrt{{\lambda }^2-k_2^2}}({\lambda }^2+\sqrt{{\lambda }^2-k_2^2}\sqrt{{\lambda }^2-k_3^2}-k_2^2)}{\sinh (d\sqrt{{\lambda }^2-k_2^2})({\lambda }^2+\sqrt{{\lambda }^2-k_1^2}\sqrt{{\lambda }^2-k_3^2}-k_2^2)+\sqrt{{\lambda }^2-k_2^2}(\sqrt{{\lambda }^2-k_1^2}+\sqrt{{\lambda }^2-k_3^2})\cosh (d\sqrt{{\lambda }^2-k_2^2})}$$

$$g^s(\mathbf{x},{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k^2}(y+d)}{e}^{i\lambda x}{e}^{\pm \sqrt{{\lambda }^2-k_0^2}{y}_0}{e}^{-i\lambda x_0}\frac{\sigma (\lambda )}{4\pi \sqrt{{\lambda }^2-k^2}}d\lambda$$

$$\varphi (\mathbf{x})$$

$$={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k^2}(y+d)}{e}^{i\lambda x}\left(\sum _{j=1}^N{q}_j{e}^{\sqrt{{\lambda }^2-k_0^2}{y}_j}{e}^{-i\lambda x_j}\right)\frac{\sigma (\lambda )}{4\pi \sqrt{{\lambda }^2-k^2}}d\lambda .$$

$$\varphi (\mathbf{x})=\sum _{p=-\infty }^{\infty }{M}_p{\Phi }_p(x,y)$$

$$M_p=\sum _{j=1}^N{q}_j{J}_p(k_0{r}_j)e^{-ip{\theta }_j},$$

$${\Phi }_p(x,y)={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k^2}(y-y_c^s+d)}{e}^{i\lambda (x-x_c^s)}{\left(\frac{\lambda -\sqrt{{\lambda }^2-k_0^2}}{k_0}\right)}^p\frac{\tilde{\sigma }(\lambda )}{4\pi \sqrt{{\lambda }^2-k^2}}d\lambda ,$$

$$e^{ikr\cos \theta }=\sum _{m=-\infty }^{\infty }{i}^m{J}_m(kr)e^{im\theta }$$

$$\sum _{j=1}^N{q}_j{e}^{\sqrt{{\lambda }^2-k_0^2}(y_j-y_c^s)}{e}^{-i\lambda (x_j-x_c^s)}$$

$${\tilde{M}}_p=\sum _{q=-\infty }^{\infty }{M}_{p-q}{J}_q(k_0{r}_{12})e^{iq{\theta }_{12}}$$

$$\varphi (\mathbf{x})=\sum _{p=-\infty }^{\infty }{L}_p{J}_p(kr)e^{ip\theta }$$

$$e^{i(\lambda (x-x_c^t)+\sqrt{k^2-{\lambda }^2}(y-y_c^t))}=\sum _{m=-\infty }^{\infty }{i}^m{J}_m(kr)e^{im\theta }{\left(\frac{\lambda -i\sqrt{k^2-{\lambda }^2}}{k}\right)}^m$$

$$L_p=\sum _{q=-\infty }^{\infty }{A}_{p,q}{M}_q$$

$$A_{p,q}={\int }_{-\infty }^{\infty }{i}^p$$

$$\times {\left(\frac{\lambda +\sqrt{{\lambda }^2-k^2}}{k}\right)}^p{\left(\frac{\lambda -\sqrt{{\lambda }^2-k_0^2}}{k_0}\right)}^q\frac{\tilde{\sigma }(\lambda )}{4\pi \sqrt{{\lambda }^2-k^2}}d\lambda .$$

$${\tilde{L}}_p=\sum _{q=-\infty }^{\infty }{L}_{p-q}{J}_q(k{r}_{12})e^{-iq({\theta }_{12}-\pi )}$$

$$g(\mathbf{x},{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k^2}y}{e}^{i\lambda x}\frac{1}{4\pi \sqrt{{\lambda }^2-k^2}}d\lambda$$

$$g(\mathbf{x},{\mathbf{x}}_0)=\frac{1}{4\pi }{\int }_0^{\pi }{e}^{ik(y\sin \theta -x\cos \theta )}d\theta +\frac{1}{2\pi i}{\int }_0^{\infty }{e}^{-ty}\cos (\sqrt{t^2+k^2}x)dt.$$

$$g(\mathbf{x},{\mathbf{x}}_0)=\frac{1}{4\pi }{\int }_0^{\pi }{e}^{ik(-y\cos \theta +x\sin \theta )}d\theta +\frac{1}{2\pi i}{\int }_0^{\infty }{e}^{-tx}\cos (\sqrt{t^2+k^2}y)dt,$$

$$G(\mathbf{x},{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k^2}(y+d)}{e}^{i\lambda x}{e}^{\pm \sqrt{{\lambda }^2-k_0^2}{y}_0}{e}^{-i\lambda x_0}\frac{\sigma (\lambda )}{4\pi \sqrt{{\lambda }^2-k^2}}d\lambda$$

$$g(\mathbf{x},{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k^2}y}{e}^{i\lambda x}{e}^{\sqrt{{\lambda }^2-k^2}{y}_0}{e}^{-i\lambda x_0}\frac{1}{4\pi \sqrt{{\lambda }^2-k^2}}d\lambda .$$

$$g^s(\mathbf{x},{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k^2}y}{e}^{i\lambda x}{e}^{-\sqrt{{\lambda }^2-k^2}{y}_0}{e}^{-i\lambda x_0}\frac{-1}{4\pi \sqrt{{\lambda }^2-k^2}}d\lambda .$$

$$\frac{\partial u}{\partial y}-i\alpha u=0$$

$$g^s(\mathbf{x},{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k^2}y}{e}^{i\lambda x}{e}^{-\sqrt{{\lambda }^2-k^2}{y}_0}{e}^{-i\lambda x_0}\frac{1}{4\pi \sqrt{{\lambda }^2-k^2}}\frac{\sqrt{{\lambda }^2-k^2}+ik\alpha }{\sqrt{{\lambda }^2-k^2}-ik\alpha }d\lambda .$$

$$g_1^s(\mathbf{x},{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k_1^2}y}{e}^{i\lambda x}{e}^{-\sqrt{{\lambda }^2-k_1^2}{y}_0}{e}^{-i\lambda x_0}\frac{1}{4\pi \sqrt{{\lambda }^2-k_1^2}}{\sigma }_1(\lambda )d\lambda ,$$

$$g_2^t(\mathbf{x},{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{\sqrt{{\lambda }^2-k_2^2}y}{e}^{i\lambda x}{e}^{-\sqrt{{\lambda }^2-k_1^2}{y}_0}{e}^{-i\lambda x_0}\frac{1}{4\pi \sqrt{{\lambda }^2-k_2^2}}{\sigma }_2^t(\lambda )d\lambda ,$$

$$g_2^b(\mathbf{x},{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k_2^2}(y+2d)}{e}^{i\lambda x}{e}^{-\sqrt{{\lambda }^2-k_1^2}{y}_0}{e}^{-i\lambda x_0}\frac{1}{4\pi \sqrt{{\lambda }^2-k_2^2}}{\sigma }_2^b(\lambda )d\lambda ,$$

$$g_3^s(\mathbf{x},{\mathbf{x}}_0)={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k_3^2}(y+2d)}{e}^{i\lambda x}{e}^{-\sqrt{{\lambda }^2-k_1^2}{y}_0}{e}^{-i\lambda x_0}\frac{1}{4\pi \sqrt{{\lambda }^2-k_3^2}}{\sigma }_3(\lambda )d\lambda$$

$$\left(\begin{array}[]{rrrr}-1&\frac{\sqrt{\lambda^{2}-k_{1}^{2}}}{\sqrt{\lambda^{2}-k_{2}^{2}}}&\frac{e^{-d\sqrt{\lambda^{2}-k_{2}^{2}}}\sqrt{\lambda^{2}-k_{1}^{2}}}{\sqrt{\lambda^{2}-k_{2}^{2}}}&0\\ 0&e^{-d\sqrt{\lambda^{2}-k_{2}^{2}}}&1&-\frac{\sqrt{\lambda^{2}-k_{2}^{2}}}{\sqrt{\lambda^{2}-k_{3}^{2}}}\\ 1&1&-e^{-d\sqrt{\lambda^{2}-k_{2}^{2}}}&0\\ 0&e^{-d\sqrt{\lambda^{2}-k_{2}^{2}}}&-1&-1\\ \end{array}\right)\left(\begin{array}[]{c}\sigma_{1}(\lambda)\\ \sigma_{2}^{t}(\lambda)\\ \sigma_{2}^{b}(\lambda)\\ \sigma_{3}(\lambda)\end{array}\right)=\left(\begin{array}[]{c}1\\ 0\\ 1\\ 0\end{array}\right),$$

$$\left(\begin{array}[]{l}\frac{\sinh\left(d\sqrt{\lambda^{2}-k_{2}^{2}}\right)\left(-\lambda^{2}+\sqrt{\lambda^{2}-k_{1}^{2}}\sqrt{\lambda^{2}-k_{3}^{2}}+k_{2}^{2}\right)+\sqrt{\lambda^{2}-k_{2}^{2}}\left(\sqrt{\lambda^{2}-k_{1}^{2}}-\sqrt{\lambda^{2}-k_{3}^{2}}\right)\cosh\left(d\sqrt{\lambda^{2}-k_{2}^{2}}\right)}{\sinh\left(d\sqrt{\lambda^{2}-k_{2}^{2}}\right)\left(\lambda^{2}+\sqrt{\lambda^{2}-k_{1}^{2}}\sqrt{\lambda^{2}-k_{3}^{2}}-k_{2}^{2}\right)+\sqrt{\lambda^{2}-k_{2}^{2}}\left(\sqrt{\lambda^{2}-k_{1}^{2}}+\sqrt{\lambda^{2}-k_{3}^{2}}\right)\cosh\left(d\sqrt{\lambda^{2}-k_{2}^{2}}\right)}\\ \frac{e^{d\sqrt{\lambda^{2}-k_{2}^{2}}}\left(\lambda^{2}+\sqrt{\lambda^{2}-k_{2}^{2}}\sqrt{\lambda^{2}-k_{3}^{2}}-k_{2}^{2}\right)}{\sinh\left(d\sqrt{\lambda^{2}-k_{2}^{2}}\right)\left(\lambda^{2}+\sqrt{\lambda^{2}-k_{1}^{2}}\sqrt{\lambda^{2}-k_{3}^{2}}-k_{2}^{2}\right)+\sqrt{\lambda^{2}-k_{2}^{2}}\left(\sqrt{\lambda^{2}-k_{1}^{2}}+\sqrt{\lambda^{2}-k_{3}^{2}}\right)\cosh\left(d\sqrt{\lambda^{2}-k_{2}^{2}}\right)}\\ -\frac{e^{d\sqrt{\lambda^{2}-k_{2}^{2}}}\left(-\lambda^{2}+\sqrt{\lambda^{2}-k_{2}^{2}}\sqrt{\lambda^{2}-k_{3}^{2}}+k_{2}^{2}\right)}{\sinh\left(d\sqrt{\lambda^{2}-k_{2}^{2}}\right)\left(\lambda^{2}+\sqrt{\lambda^{2}-k_{1}^{2}}\sqrt{\lambda^{2}-k_{3}^{2}}-k_{2}^{2}\right)+\sqrt{\lambda^{2}-k_{2}^{2}}\left(\sqrt{\lambda^{2}-k_{1}^{2}}+\sqrt{\lambda^{2}-k_{3}^{2}}\right)\cosh\left(d\sqrt{\lambda^{2}-k_{2}^{2}}\right)}\\ \frac{2\sqrt{\lambda^{2}-k_{2}^{2}}\sqrt{\lambda^{2}-k_{3}^{2}}e^{d\sqrt{\lambda^{2}-k_{3}^{2}}}}{\sinh\left(d\sqrt{\lambda^{2}-k_{2}^{2}}\right)\left(\lambda^{2}+\sqrt{\lambda^{2}-k_{1}^{2}}\sqrt{\lambda^{2}-k_{3}^{2}}-k_{2}^{2}\right)+\sqrt{\lambda^{2}-k_{2}^{2}}\left(\sqrt{\lambda^{2}-k_{1}^{2}}+\sqrt{\lambda^{2}-k_{3}^{2}}\right)\cosh\left(d\sqrt{\lambda^{2}-k_{2}^{2}}\right)}\end{array}\right).$$

$$\varphi (\mathbf{x})=\sum _{p=-\infty }^{\infty }{M}_p{\Phi }_p(x,y),M_p=J_p(k_{0}r)e^{-ip\theta }$$

$${\Phi }_p(\mathbf{x})={\int }_{-\infty }^{\infty }{e}^{-\sqrt{{\lambda }^2-k^2}(y-y_c^s+d)}{e}^{i\lambda (x-x_c^s)}{\left(\frac{\lambda -\sqrt{{\lambda }^2-k_0^2}}{k_0}\right)}^p\frac{\tilde{\sigma }(\lambda )}{4\pi \sqrt{{\lambda }^2-k^2}}d\lambda$$

Full text

1 Introduction

The efficient numerical simulation of waves in layered media in the frequency domain is a challenging task in scientific computing. One of the key numerical difficulties arises from the Sommerfeld integral representation for a general layered media Green’s function [1, 2]. For example, consider a simple three-layered medium with layer interfaces at $y=0$ and $y=-d$, and assume $k_{1},k_{2}$, and $k_{3}$ are the wave numbers in each layer, respectively. When the source is located in the top layer $\mathbf{x}_{0}=(x_{0},y_{0})$ with $y_{0}>0$, the scattered field at $\mathbf{x}_{1}=(x_{1},y_{1})$ in the middle layer ($-d<y_{1}<0$) has the reflected field from the top and bottom interfaces. The reflected field from the top interface $g_{2}^{t}$ (formula for the bottom interface is similar) is represented by

[TABLE]

where

[TABLE]

is referred to as the density on the top interface found by matching the continuity conditions at the layer interfaces [3]. A discretization of the integral equation description of the wave field in layered media leads to a linear system and each entry in the coefficient matrix requires the evaluation of one or more of these Sommerfeld integrals [4, 5], which is very expensive especially when the source and target are close to the same interface, e.g., when both $y_{1}$ and $y_{0}$ are close to the upper ($y=0$) interface and $|x_{1}-x_{0}|\gg|y_{1}-y_{0}|$ in the function $g_{2}^{t}$. Note that even the state-of-the-art fast direct solvers or $\cal{H}$-matrix algorithms [6, 7, 8, 9] require a sufficient number of samples of the coefficient matrix entries before the matrix can be effectively compressed for further numerical linear algebra operations.

This paper is aimed at completely resolving two of the many existing challenges in the simulation of waves in layered media: (a) an optimal fast algorithm for the application of the layered media Green’s function to a given density function, and (b) an effective numerical scheme to compute the layered media Green’s function when the source and target are close to each other, or close to the media interface. Firstly, we present a general algorithm framework for an optimal fast solver. The new algorithm compresses a transformed version of the original matrix, and both the expansions and translation operators are derived using mathematical analysis techniques. The transformed matrix approach made possible for the new algorithm to use many well-optimized numerical tools from existing free-space fast multipole method (FMM) with minor or no changes. This approach is different from existing FMM, fast direct solver, and fast $\mathcal{H}$-matrix algorithms which compress the coefficient matrix directly. The error analysis of the expansion is based on the asymptotic analysis of integrals which can be generalized to layered media Green’s functions in three dimensions. Secondly, by introducing different integration contours for the Sommerfeld integral representation of the layered media Green’s function, the new method can also effectively and accurately evaluate the interactions when either the source, or target, or both are close to the interface between different layers. The resulting algorithm complexity is asymptotically optimal $\mathcal{O}(N)$ in the low frequency regime, with a prefactor close to that of the well-developed free-space FMM algorithms.

To present the ideas, we restrict our attention to the two-dimensional (2-D) layered media in this paper and organize the paper as follows. In Section 2, we first summarize all the necessary building blocks for the new algorithm and present the pseudocode and complexity analysis. Section 3 provides several examples of these Green’s functions. In Section 4, we present a detailed analysis of the numerical algorithm, including the asymptotic expansion based truncation error analysis of the new “multipole” (far-field) and “local” expansions, and the mathematically equivalent alternative direction Sommerfeld integral representations of the original layered media Green’s function. Section 5 presents numerical results designed to validate the algorithm analysis. We summarize our results and discuss future work in Section 6.

2 Adapting Free-space FMM for 2-D Layered Media Green’s Functions

We first present the algorithm framework. Detailed analysis of each algorithm component will be covered in Section 4. This section is written for readers with sufficient knowledge of the classical FMM [10], which is becoming a standard topic in scientific computing.

2.1 Summary of Algorithm Building Blocks

We use notations and terminologies commonly adopted by the FMM community, and present the building blocks in the same order as they appear in the algorithm.

Layered Media Green’s Functions.

The layered media Green’s function consists of the free-space and scattered field parts. As the free-space Green’s function is well-studied, we focus on a general form of the scattered field Green’s function. For a target point $\mathbf{x}=(x,y)$ in the layer with wave number $k$ and a source point $\mathbf{x}_{0}=(x_{0},y_{0})$ in the layer with wave number $k_{0}$, the scattered field Green’s function is given by

[TABLE]

where $d$ is a constant and $\sigma(\lambda)$ is independent of $\mathbf{x}$ and $\mathbf{x}_{0}$ and converges to a constant when $\lambda\to\pm\infty$. When the “+” sign is used in $e^{\pm\sqrt{\lambda^{2}-k_{0}^{2}}y_{0}}$, $y+d-y_{0}$ is assumed to be positive in order to guarantee integrability. Sample Sommerfeld integral representations of layered media Green’s functions for acoustic waves and time harmonic Maxwell’s equations (vector Helmholtz equations) can be found in [3, 11, 12, 13] and in Section 3. In the following, the mathematical formulas for the “$+$” sign are presented. The formulas for the “$-$” sign can be derived in the same way and are thus omitted.

Multipole Expansion and Source-to-Multipole (S2M) Translation Operator. Assume there are $N$ sources located in a box centered at $\mathbf{x}_{c}^{s}=(x^{s}_{c},y^{s}_{c})$, each carries a charge $q_{j}$ located at $\mathbf{x}_{j}=(x_{j},y_{j})$, their contributions to the far-field location $\mathbf{x}=(x,y)$ are given by

[TABLE]

Definition 1** (Multipole Expansion)**

The far-field contributions at $\mathbf{x}$ due to charges $q_{j}$ located at $\mathbf{x}_{j}$, $j=1,2,\cdots N$ in a box centered at $\mathbf{x}_{c}^{s}$ in Eq. (2) have the multipole expansion:

[TABLE]

where the multipole coefficient is given by

[TABLE]

$(r_{j},\theta_{j})$* are the polar coordinates of $\mathbf{x}_{j}$ with respect to $\mathbf{x}_{c}^{s}$, $J_{p}$ is the Bessel function of order $p$,*

[TABLE]

and ${\tilde{\sigma}(\lambda)=\sigma(\lambda)e^{-\sqrt{\lambda^{2}-k^{2}}y^{s}_{c}+\sqrt{\lambda^{2}-k_{0}^{2}}y^{s}_{c}}}$.

Eq. (4) is derived from the Jacobi-Anger formula [14]

[TABLE]

and the multipole coefficient $M_{p}$ is the same as the multipole coefficient for the free-space Helmholtz kernel [15]. We refer to forming the multipole expansion as the Source-to-Multipole (S2M) operator. This definition simply states that by changing the basis (matrix transformation) to $\Phi_{p}$, existing free-space S2M operator for the Helmholtz equation with the same wave number $k_{0}$ can be used directly to derive the coefficients $M_{p}$ of the compressed representation in Eq. (3). Note that the asymptotic properties of $\tilde{\sigma}(\lambda)=\sigma(\lambda)e^{-\sqrt{\lambda^{2}-k^{2}}y^{s}_{c}+\sqrt{\lambda^{2}-k_{0}^{2}}y^{s}_{c}}$ remain the same as the original $\sigma(\lambda)$ as $\lambda\to\pm\infty$.

Remark: By introducing the change of variable $z=\frac{\lambda-\sqrt{\lambda^{2}-k_{0}^{2}}}{k_{0}}$ in Eq. (5), the expansion $\sum_{p=-\infty}^{\infty}M_{p}z^{p}$ is the Laurent expansion of the function

[TABLE]

and $M_{p}$ is the expansion coefficient independent of $\lambda$ (or $z$).

At the low-frequency regime when both $k_{0}$ and $k$ are small, the number of terms required in the truncated layered media multipole expansion for a prescribed error tolerance is approximately the same as that in the free-space Laplace FMM. We leave the truncation error analysis of the multipole (and local) expansions to Section 4.

Multipole-to-Multipole (M2M) Translation Operator. In FMM, multipole expansion of the parent is constructed by translating that of the children. Translating the center of a multipole expansion from child box to its parent box is referred to as the Multipole-to-Multipole (M2M) translation operator. As the multipole coefficients for the layered media in Eq. (4) are the same as those for the free-space Green’s function for the Helmholtz equation with wave number $k_{0}$ for both the parent and child boxes, we therefore have the following lemma.

Lemma 1** (M2M)**

The M2M translation operator for the layered media Green’s function in Eq. (4) is the same as the M2M operator for the free-space Green’s function of the Helmholtz equation with wave number $k_{0}$. The parent’s multipole coefficients $\widetilde{M}_{p}$ are given by

[TABLE]

where $M_{p-q}$ are the child’s multipole coefficients and $(r_{12},\theta_{12})$ are the polar coordinates of the child’s center with respect to the parent’s center.

Therefore, the free-space M2M translation operator can be used without any change to obtain the multipole expansions for all the boxes on the hierarchical tree structure for the layered media Green’s function.

Local Expansion and Multipole-to-Local (M2L) Translation Operator. Notice that the potential field $\phi(\mathbf{x})$ satisfies the Helmholtz equation with wave number $k$, we therefore use the same Bessel function based expansion as that for the free-space Green’s function with wave number $k$ to compress the received far-field contributions into a local expansion of the target box centered at $\mathbf{x}_{c}^{t}=(x_{c}^{t},y_{c}^{t})$.

Definition 2** (Local Expansion)**

The potential function $\phi$ due to the far-field source contributions can be compressed into a local expansion

[TABLE]

where $(r,\theta)$ are the polar coordinates of $\mathbf{x}$ with respect to the target box center $\mathbf{x}_{c}^{t}$, $L_{p}$ is called the local expansion coefficient. Evaluating the local expansion at a target point is referred to as the Local-to-Target (L2T) translation operator.

Similar to the free-space FMM, the compressed far-field multipole expansion of the source box centered at $\mathbf{x}_{c}^{s}$ given by Eq. (3) can be translated into a local expansion in Eq. (7) of the target box centered at $\mathbf{x}_{c}^{t}$, by plugging the plane wave formula

[TABLE]

in the basis function $\Phi_{p}$ in Eq. (5), where $(r,\theta)$ are the polar coordinates of $\mathbf{x}$ with respect to the target box center at $\mathbf{x}_{c}^{t}$. Translating the multipole expansion to a local expansion is referred to as the Multipole-to-Local (M2L) translation operator and we have the following lemma.

Lemma 2** (M2L)**

The local expansion coefficients of a target box centered at $\mathbf{x}_{c}^{t}$, due to the contributions from particles in a source box centered at $\mathbf{x}_{c}^{s}$ described by its multipole expansion in Eq. (3), can be computed using the M2L translation matrix $\mathrm{A}=\{A_{p,q}\}$ using

[TABLE]

where

[TABLE]

Note that when $|y_{c}^{t}-y^{s}_{c}+d|\ll|x_{c}^{t}-x^{s}_{c}|$, the evaluation of $A_{p,q}$ can be very costly due to the highly oscillatory term $e^{i\lambda(x_{c}^{t}-x^{s}_{c})}$. This difficulty can be resolved by using the alternative direction integral representation of the translation operator that will be discussed in Sec. 4. The translation matrix $\mathrm{A}$ can be precomputed for optimal efficiency when the geometries of the layered media and embedded objects are fixed, or computed on the fly using the Gauss and Laguerre quadrature rules for the alternative direction integral representations.

Local-to-Local (L2L) Translation Operator. The local expansion of the parent can be translated to its children using the Local-to-Local (L2L) translation operator given by the following lemma [15, 16].

Lemma 3** (L2L)**

As the basis for the local expansion of the layered media Green’s function is the same as that for the free-space Green’s function of the Helmholtz equation with wave number $k$, the L2L translation operator for the layered media Green’s function in Eq. (7) is therefore the same as the L2L operator of the free-space FMM for the Helmholtz equation with wave number $k$. The child’s local coefficients $\widetilde{L}_{p}$ are given by

[TABLE]

where $L_{p-q}$ are the parent’s local expansion coefficients and $(r_{12},\theta_{12})$ are the polar coordinates of the parent’s box center with respect to the child’s box center.

This lemma states that the free-space L2L translation operator can be used without change to derive the local expansion of the child box from its parent’s for the layered media case.

Direct Source-to-Target (S2T) Interactions. For neighboring boxes, the interaction of the source and target can be handled in two ways: (1) by evaluating the Sommerfeld integral directly, or (2) for the scattered field contribution from sources in a source box, when its multipole expansion is also valid in its neighboring target box, the source box’s multipole expansion can be translated and merged into the target box’s local expansion using the same M2L translation operator as we discussed previously for the interaction list boxes, which will be evaluated later using the very efficient L2T operator for the leaf boxes. Note that when this happens, the FMM tree structure should be modified accordingly to further accelerate the computation.

One numerical difficulty of evaluating the layered media Green’s function for direct S2T interaction or entries of the translation matrix $\mathrm{A}$ in M2L comes from the oscillatory term $e^{i\lambda x}$ when the other exponential terms in the integrand decay slowly. For example, in the free-space Green’s function,

[TABLE]

where the source point is located at the origin $(x_{0},y_{0})=(0,0)$ and the target point is located at the first quadrant ($x>0$ and $y>0$). In the new version of FMM [16, 17, 18], it is referred to as the “north” plane wave expansion because this formula is valid when $y>0$. However, when $y$ is very close to the line $y=0$ and $x\gg 0$, a huge number of quadrature points has to be used to resolve the oscillatory term $e^{i\lambda x}$ due to the very slow decay of $e^{-\sqrt{\lambda^{2}-{k}^{2}}y}$. Similar problems arise when evaluating the direct interactions of the source and target points that are close to the interface of the layered media and the M2L translation matrix $\mathrm{A}$ in Eq. (9).

To understand the origin of this problem, we divide the Sommerfeld integral representation into the propagating ($|\lambda|<k$) and evanescent ($|\lambda|<k$) parts, and after the change of variables as in [17], Eq. (12) can be rewritten as

[TABLE]

The first and second integrals are called the propagating part and evanescent part, respectively. For the propagating part, the number of required Gauss quadrature points to evaluate the integral only depends on $k$ and $r=\sqrt{x^{2}+y^{2}}$, which is normal. However, when $y\ll x$, the integrand in the evanescent part $e^{-ty}\cos(\sqrt{t^{2}+k^{2}}x)$ requires many Laguerre quadrature points to resolve the oscillatory term $\cos(\sqrt{t^{2}+k^{2}}x)$ because $e^{-ty}$ decays slowly. This problem can be resolved by using the equivalent “east” plane wave representation for $x>0$

[TABLE]

which can be derived from the “north” plane wave representation using contour integration as will be discussed in Section 4. Since $y\ll x$, the evanescent part of the alternative direction integral can be evaluated using a small number of Laguerre quadrature points. Similar representations can also be derived for the “south” and “west” directions for general layered media Green’s functions, either using the method of images, contour integration, or applying the integral transforms directly in the alternative direction.

2.2 Algorithm Pseudocode

We present the algorithm pseudocode in Algorithm 1. Compared with the classical FMM for free-space kernels, the adapted FMM for layered media Green’s function only differs in the M2L and S2T subroutines, and in the number of expansion terms when the wave numbers are different. All other subroutines and functions from free-space FMM can be adopted with minor or no changes by the layered media FMM.

2.3 Algorithm Complexity

We compare the algorithm complexity of the layered media FMM (LM-FMM) with that of the free-space FMM (FS-FMM).

In the upward pass, only the free-space S2M and M2M translation operators are used. Thus, if the number of expansion terms is the same as that of the free-space case, the LM-FMM has the same number of operations as the FS-FMM in the upward pass. In the downward pass, the L2L operator of the LM-FMM has the same complexity as that of the FS-FMM. However, the M2L operator requires more operations when the M2L translation matrix is computed on the fly, as at least one integral has to be evaluated to find each entry of the translation matrix. Note that these integrals can be evaluated efficiently when the alternative direction integral representations are used (further discussed in Sec. 4). Moreover, the translation matrix can often be re-used by many boxes in the same level in the hierarchical tree structure. On the other hand, when the translation matrix $\mathrm{A}$ is precomputed, the algorithm complexity in the downward pass is about the same as that in the classical FS-FMM. In Step 4, the evaluation of the local expansion in the LM-FMM has the same complexity as that in the FS-FMM. For the direct source to target (S2T) interactions, more operations are needed than FS-FMM because the alternative direction Sommerfeld integrals have to be evaluated. In summary, when the local direct interaction operations are not counted (which can be very efficient on parallel computers), the translation matrix $\mathrm{A}$ is precomputed, and the number of expansion terms is the same as that of the FS-FMM, the total number of operations of the LM-FMM is about the same as that of the FS-FMM. In Sec. 4, we show that for many settings, the number of expansion terms of the LM-FMM is in the same order as that of the FS-FMM, for example, when the wave numbers of different layers are all in the low-frequency regime.

3 Examples of 2-D Layered Media Green’s Functions

We consider the potential at a target point $\mathbf{x}=(x,y)$ due to a source charge with density $q_{0}$ located at $\mathbf{x}_{0}=(x_{0},y_{0})$. The source and target may locate in the same or different layers with the $j^{th}$ interface located at $y=y_{j}$. The wave numbers are $k_{0}$ for the source layer and $k$ for the target layer, respectively.

To derive the layered media Green’s function, a Fourier transform is usually performed in the $x$-direction, reducing the 2-D Helmholtz equation to an ODE system, which can be solved analytically with some unknown density functions in the Sommerfeld integral representations. The density functions are obtained by solving a linear system of algebraic equations to match the interface conditions. In this paper, we consider the Green’s function in a rather general form

[TABLE]

where $d$ is a constant and $\sigma(\lambda)$ converges to a constant when $\lambda\to\pm\infty$. Many 2-D layered media Green’s functions are in this particular form. In the following discussions, we refer to the first two exponential terms $e^{-\sqrt{\lambda^{2}-k^{2}}(y+d)}e^{i\lambda x}$ as the target term, the third and fourth exponential terms $e^{\pm\sqrt{\lambda^{2}-k_{0}^{2}}y_{0}}e^{-i\lambda x_{0}}$ as the source term, and $\sigma(\lambda)$ as the image term (for reasons which will be explained later).

In this section, we present a few examples of the layered media Green’s functions, all are in the form of Eq. (13).

Example 1: Free-space Green’s Function. The first example is the free-space Green’s function for the Helmholtz equation with wave number $k$. For a source point $\mathbf{x}_{0}=(x_{0},y_{0})$ and a target point $\mathbf{x}=(x,y)$ with $y-y_{0}>0$, the free-space Green’s function is given by the Sommerfeld integral of the form

[TABLE]

Example 2: Two-layered Medium with Zero Dirichlet Interface Condition. The second example is the half-space problem with Dirichlet condition (total field is [math]) at the layer interface located at $y=0$. We assume the source is located at $\mathbf{x}_{0}=(x_{0},y_{0})$ on the upper half plane ($y_{0}>0$). Using the method of images, the scattered field for a target point located at $\mathbf{x}=(x,y)$ on the upper half plane ($y>0$) can be represented as

[TABLE]

To satisfy the Dirichlet boundary condition, an image charge is added at the location $(x_{0},-y_{0})$ and $-1$ is the charge of the image source. Interested readers are referred to [19] for more details.

Example 3: Two-layered Medium with Impedance Interface Condition. The third example is the half-space Green’s function with the impedance boundary condition

[TABLE]

at the layer interface $y=0$. The scattered field at a target point $\mathbf{x}=(x,y)$ due to the source at $\mathbf{x}_{0}=(x_{0},y_{0})$, where both points are located at the upper half plane with wave number $k$, is given by

[TABLE]

In [13, 20], it was shown that the term $\frac{\sqrt{\lambda^{2}-k^{2}}+ik\alpha}{\sqrt{\lambda^{2}-k^{2}}-ik\alpha}$ can be derived using the method of complex images, we therefore refer to $\sigma(\lambda)=\frac{\sqrt{\lambda^{2}-k^{2}}+ik\alpha}{\sqrt{\lambda^{2}-k^{2}}-ik\alpha}$ term as the image term. This Green’s function was also discussed in [19]. Note that when $\lambda\to\pm\infty$, $\sigma\to 1$.

Example 4: Three-layered Medium with Transmission Condition. The last example we consider in this paper is the Green’s function for a three-layered medium with the transmission condition in [3], where the layer interfaces are located at $y=0$ and $y=-d$. Let a source point be located in the top layer at $\mathbf{x}_{0}=(x_{0},y_{0})$ and a target point $\mathbf{x}=(x,y)$. In the first layer ($y>0$), the scattered field $g_{1}^{s}$ is the field reflected from bottom layers. In the second layer ($-d<y<0$), the scattered field consists of the contribution from top and bottom interfaces $g_{2}^{s}=g_{2}^{t}+g_{2}^{b}$. In the third layer ($y<-d$), the scattered field $g_{3}^{s}$ is the transmitted field from the source in the first layer. By matching the continuity of the field, the scattered field in each layer can be represented as

[TABLE]

where $(\sigma_{1}(\lambda),\sigma_{2}^{t}(\lambda),\sigma_{2}^{b}(\lambda),\sigma_{3}(\lambda))^{T}$ are the solutions of the linear system

[TABLE]

and the solutions $(\sigma_{1}(\lambda),\sigma_{2}^{t}(\lambda),\sigma_{2}^{b}(\lambda),\sigma_{3}(\lambda))^{T}$ are explicitly given by

[TABLE]

Note that when $\lambda\to\pm\infty$, all the $\sigma$ functions converge to constants.

The form of the Green’s function in Eq. (13) is not surprising. Clearly the first target term satisfies the Helmholtz equation at the target layer with wave number $k$ and the source term satisfies the Helmholtz equation at the source layer with wave number $k_{0}$. The third term is independent of the variables $\mathbf{x}_{0}$ and $\mathbf{x}$, and we collect all the exponential growth or decay terms in the constant $d$ in the target term so the image term converges to a constant when $\lambda\to\pm\infty$.

4 Analysis of Layered Media Fast Multipole Method

We present a detailed analysis of the algorithm for layered media Green’s function in this section. We focus on the following two topics: (a) the number of terms in the multipole and local expansions and truncation errors, and (b) evaluation of the local direct interactions and M2L translation operators using the mathematically equivalent alternative direction Sommerfeld integral representation.

4.1 Truncating the Multipole and Local Expansions

We study the decay rates of the terms in the multipole and local expansions by considering the setting of a single source with a unit charge located at $\mathbf{x}_{0}=(x_{0},y_{0})$. The multipole expansion describes the potential at far-field locations as a function of $\mathbf{x}=(x,y)$ due to a charge in the source box, namely,

[TABLE]

where $(r,\theta)$ are the polar coordinates of the point $(x_{0},y_{0})$ with respect to the source box center $\mathbf{x}_{c}^{s}=(x_{c}^{s},y_{c}^{s})$, the basis function $\Phi_{p}(\mathbf{x})$ is given by

[TABLE]

and we define the “modified” distance between $\mathbf{x}$ and $\mathbf{x}_{c}^{s}$ of the multipole expansion as $\rho=\sqrt{(y-y^{s}_{c}+d)^{2}+(x-x^{s}_{c})^{2}}$. The multipole expansion can be considered as a compressed representation of the source box’s contribution to be sent to the far-field locations. Similarly, the local expansion associated with each target box compresses the received far-field contributions and describes the potential as a function of $(x-x_{c}^{t},y-y_{c}^{t})$ given by

[TABLE]

where $(\tilde{r},\tilde{\theta})$ are the polar coordinates of the point $(x,y)$ with respect to the target box center $\mathbf{x}_{c}^{t}=(x_{c}^{t},y_{c}^{t})$ and the basis function $\Psi_{p}(\mathbf{x}_{0})$ is given by

[TABLE]

where ${\accentset{\approx}{\sigma}(\lambda)=\sigma(\lambda)e^{-\sqrt{\lambda^{2}-k^{2}}y^{t}_{c}+\sqrt{\lambda^{2}-k_{0}^{2}}y^{t}_{c}}.}$ We define the “modified” distance between $\mathbf{x}_{0}$ and $\mathbf{x}_{c}^{t}$ of the local expansion as $\tilde{\rho}=\sqrt{(y^{t}_{c}-y_{0}+d)^{2}+(x^{t}_{c}-x_{0})^{2}}$.

Similar to the FS-FMM where the truncation errors of the multipole and local expansions are determined by the decay rate of $|J_{p}(kr)H_{p}(k\rho)|$, we study the decay rate of the term $|J_{p}(k_{0}r)\Phi_{p}(x,y)|$ for the multipole expansion, and $|J_{p}(kr)\Psi_{p}(x_{0},y_{0})|$ for the local expansion for different physical parameter settings. When all other variables are fixed, $J_{p}\to 0$ and $H_{p},\Phi_{p},\Psi_{p}\to\infty$ as $p\to\infty$. Therefore, it is necessary to understand the asymptotic behavior of these functions for large $p$ values.

Asymptotic Forms of Bessel Functions for Large Order. The asymptotic expansion for large order Bessel functions is a well-studied topic. We cite the following well-known results from [14], which are valid for fixed $z$ when $\nu\to\infty$.

[TABLE]

[TABLE]

Therefore for the free-space Green’s function,

[TABLE]

and the truncation errors of both the free-space multipole and local expansions decay exponentially as $p\to\infty$ [15].

Asymptotic Approximation of Integrals. The asymptotic expansion for the Bessel functions can be derived from the integral representations of these special functions. Note that these integral representations are similar to the layered media Green’s functions, therefore the same asymptotic analysis techniques can be applied, and the results can be used to derive more precise error bounds when truncating the layered media multipole and local expansions. We demonstrate the ideas using the following examples.

We start with the integral representation of the Bessel function [14],

[TABLE]

Note that when $k<n$, the integral $\int_{-\pi}^{\pi}(\cos\theta)^{k}\cos(n\theta)d\theta=0$. By applying the Stirling’s formula $n!\sim\sqrt{2\pi n}(\frac{n}{e})^{n}$ to the leading order of the expansion, one recovers the asymptotic form of the Bessel function for large orders in Eq. (21). This asymptotic expansion is valid for all $z$ values and provides more accurate estimate when studying the truncation error in the multipole and local expansions, especially in the low frequency regime when $k_{0}r$ and $kr$ are small.

For the layered media Green’s function, both integrals in the basis functions $\Phi_{p}$ and $\Psi_{p}$ can be formulated as

[TABLE]

where $k_{1}$ is the wave number in the target layer and $k_{2}$ is the wave number in the source layer. The integral can be divided into three parts, when (a) $|\lambda|<\min\{k_{1},k_{2}\}$, (b) $|\lambda|>\max\{k_{1},k_{2}\}$, and (c) $\min\{k_{1},k_{2}\}\leq|\lambda|\leq\max\{k_{1},k_{2}\}$. The asymptotic expansion for each part can be derived using existing asymptotic analysis techniques for integrals [21]. When $k_{1}=k_{2}$, the integral in (a) is often referred to as the “propagating” part and the integral in (c) becomes the “evanescent” part. The asymptotic properties of the basis functions $\Phi_{p}$ and $\Psi_{p}$ are determined by the “propagating” part for large $||\mathbf{x}||$ values with fixed $p$, and by the “evanescent” part for large $p$ values with fixed $||\mathbf{x}||$. To understand the truncation errors in the multipole and local expansions of the layered media Green’s function, we therefore focus on the evanescent part and demonstrate the asymptotic analysis for the simplified integral

[TABLE]

where we assume $k_{1}>k_{2}$. The case when $k_{1}<k_{2}$ can be analyzed in a similar way. Instead of deriving the asymptotic expansion of Eq. (24) directly, we adopt the following steps to further simplify the integral to a more standard form that is commonly used when analyzing the asymptotic behavior of the Bessel functions. First, using the polar coordinates $(r,\theta)$ of $(x,y)$, we rewrite the $\{x,y\}$-related exponential part as

[TABLE]

Second, we define a new variable $u=\frac{t\sin(\theta)}{k_{1}}-\frac{i\cos(\theta)\sqrt{k_{1}^{2}+t^{2}}}{k_{1}}$. Clearly, the new integral for the $u$-variable is on a complex contour, not on the real axis. Third, using contour integration and residue theorem, we can rewrite the integral on the complex contour back to the sum of an integral on the real axis and an easy-to-analyze integral on a line segment in the complex plane (which will be explained in detail in Section 4.2). Finally, we focus on the asymptotic expansion of the dominating integral

[TABLE]

where $\beta(u)=k_{1}\left(\sqrt{u^{2}+1}\sin(\theta)-iu\cos(\theta)\right)$ and $\accentset{\approx}{\sigma}(u)=\tilde{\sigma}(t(u))$. For any fixed $z=k_{1}r$, one approach to find the asymptotic behavior of this integral is to first perform another change of variable $\lambda=k_{1}ru=zu$ (or $u=\frac{\lambda}{k_{1}r}=\frac{\lambda}{z}$) and consider the new integral

[TABLE]

where $\tilde{\beta}(z)=\sqrt{\lambda^{2}+z^{2}}\sin(\theta)-i\lambda\cos(\theta).$ Under proper conditions, the function $\left(\tilde{\beta}(z)+\sqrt{\tilde{\beta}^{2}(z)-(\frac{k_{2}}{k_{1}}z)^{2}}\right)^{p}\tilde{g}(\frac{\lambda}{z})$ is an analytic function for $z$ values on the right half complex plane away from the origin, we can consider its Taylor expansion

[TABLE]

as a function of $z$. Then, the integral representation of evanescent part can be derived as

[TABLE]

The leading order term of the expansion for very large $p$ values is approximately

[TABLE]

for some constant $\theta_{0}$, which has the same asymptotic properties as the Hankel function $H_{p}(k_{2}r)e^{ip\theta_{0}}.$ Note that the expansion can be used to study the properties of both the cases when $p\to\pm\infty$ or when $z$ is small (low-frequency regime). Without presenting the details, we summarize our results in the following theorem.

Theorem 1

(1)* The multipole and local expansions satisfy the following convergence estimates for large $p$ values,*

[TABLE]

*where $r$ is the distance from the particle to its box center and $\rho$ (or $\tilde{\rho}$) is the “modified” distance between the box center and far-field point. Both $d$ and $\pm$ sign are considered to correct the $\rho$ (or $\tilde{\rho}$) value. For example, $\rho=\sqrt{(y-y_{c}^{s}+d)^{2}+(x-x_{c}^{s})^{2}}$ for $\Phi(x,y)$ in Eq. (5) and $\tilde{\rho}=\sqrt{(y_{c}^{t}+d-y_{0})^{2}+(x_{c}^{t}-x_{0})^{2}}$ for $\Psi(x_{0},y_{0})$ in Eq. (20).

(2) When the wave numbers $k_{1}$ and $k_{2}$ are small, the required number of terms in the multipole and local expansions for a prescribed accuracy in the layered media FMM is approximately the same as that in the free-space FMM.*

This theorem simply states that both the multipole and local expansions in the layered media FMM are exponentially convergent and presents the asymptotic convergence rates. In the low-frequency regime, it provides more precise estimates of the number of terms required.

Remark: Note that using the “modified” distance, the truncated multipole expansion of the scattered field Green’s function may become valid for a neighboring (including self) target box, and therefore can be translated to its neighbor’s local expansion using the M2L operator instead of the more expensive S2T operator.

Laplace Transform and Complex Images. For many layered media Green’s functions, one can justify that the inverse Laplace transform of $\accentset{\approx}{\sigma}(z)$ can be derived using the Fourier-Mellin integral formula

[TABLE]

where the integration is done along the vertical line $\text{Re}(z)=\gamma$ in the complex plane and $\gamma$ is greater than the real part of all singularities of $\accentset{\approx}{\sigma}(z)$, so that the image part in the layered media Green’s function can be represented as

[TABLE]

for some $\gamma_{0}$ which may or may not be the same as $\gamma$. Plugging this representation in the original layered media Green’s function in Eq. (13), we have

[TABLE]

In this new representation, $f(t)$ can be considered as the complex image located at $(x_{0},\pm y_{0}-t)$. For many layered media, it is sufficient to analyze the convergence of the multipole and local expansions of the layered media Green’s functions for each $t$-mode where the new source box center is at $(x_{c}^{s},\pm y_{c}^{s}-(t-\gamma_{0}))$. In [20], this complex image approach is applied to the 2-D half-space layered medium Green’s function with impedance boundary conditions (see Example 3 of Sec. 3), and it becomes straightforward to verify using the complex images that for the same prescribed accuracy requirement, the number of expansion terms for the layered medium case is no more than that for the free-space case. Another application of the complex image approach is for the direct source-target (S2T) interactions. For a neighboring source box, when all the complex images of the scattered field are well-separated from the target box, the source box’s multipole expansion becomes applicable, and it is more efficient to translate the multipole expansion to a local expansion of the target box. This technique is used in [20] to compute the scattered field part of the source-target interactions.

Precomputed Tables for Number of Expansion Terms. In practice, both the asymptotic expansion and complex image approach only give a very rough estimate of the number of terms required in different expansions of the layered media Green’s function. A more practical approach is to precompute a table (or table of tables) for different layered media settings. This approach is problem dependent. We are constructing such tables for several real world applications, and results will be reported in the future.

4.2 Alternative Direction Sommerfeld Integral Representations

Another difficulty in simulating waves in layered media is when the source and target are very close to each other or when the source and target are close to the interface of different layers for the scattered field. In these cases, the computation of layered media Green’s function and M2L translation operator becomes extremely expensive. For example, when $y-y_{0}>0$ is close to zero and $|x-x_{0}|$ is relatively a big number in the Sommerfeld integral representation of free-space Green’s function in Eq. (14) or the half-space layered medium Green’s function with impedance boundary conditions in Eq. (17), the exponential term $e^{-\sqrt{\lambda^{2}-k^{2}}(y-y_{0})}$ decays slowly and the highly oscillatory term $e^{i\lambda(x-x_{0})}$ must be sufficiently sampled. A similar problem occurs in the M2L operator in Eq. (9). When the ratio of $(y_{c}^{t}-y_{c}^{s}+d)/|x_{c}^{t}-x_{c}^{s}|$ is small, a wide range of $\lambda$ values have to be sampled before the integrand decays to zero sufficiently.

The numerical difficulty is not from any inherent properties of the original physical problem, it is the result of using the inefficient integral representation in the numerical computation. For the same Hankel function $H_{0}(\beta r)$, which is the Green’s function of the free-space Helmholtz equation, a common practice is to divide ${\bf R}^{2}$ into four overlapping regions–North, South, East, West–corresponding to points $(x,y)\in{\bf R}^{2}$ with $y>0,y<0,x>0,x<0$, respectively. In each region, the plane wave representation of $H_{0}(\beta r)$ takes the following forms:

[TABLE]

where $\rho_{\beta}(t)=\sqrt{t^{2}+\beta^{2}}$ and $(r,\theta)$ are the polar coordinates of the point $(x,y)$. For higher order Hankel functions $H_{l}(\beta r)e^{il\theta}$, we have the following integral representations.

[TABLE]

where the first and the second integrals in each formula are called the propagating part and evanescent part, respectively. These directional representations have been applied in the new version of FMM in both two and three dimensions [16, 18, 22, 23]. They are also effective tools to compute the lattice sums of the free-space Green’s functions [24, 25]. Clearly in these formulas, the oscillation of the integrand in the propagating part is controlled by $\beta r$ and hence no numerical quadrature issues arise. For the evanescent part at the overlapping regions, e.g., when both $x,y>0$, both the north and east integral representations are valid and can be applied, however their numerical properties are very different. Clearly when $y/x\gg 1$, the north formula is preferred, and when $y/x\ll 1$, the east formula can be computed more efficiently using existing quadrature techniques.

The alternative direction integral can be derived using contour integration and residue theorem in two dimensions. To demonstrate how this technique works, we consider the contour integral for $x,y>0$ in fig. 1

[TABLE]

Because there are no singularities or branch cut points inside the contour $\Gamma$, we have

[TABLE]

When $\epsilon\to 0$ and $R\to\infty$,

[TABLE]

and

[TABLE]

Therefore, the evanescent part in $\int_{II}$ can be computed using $-\left(\int_{I}+\int_{III}\right)$. In fact, the east plane wave representation can be re-derived from the north one using this approach for $x,y>0$.

Finally, we point out that the proper choice of the integration contour is problem dependent. We present more details in Section 5 for other layered media Green’s functions.

5 Preliminary Numerical Experiments

In this section, we present numerical experiments to validate our theoretical analysis of the general numerical framework. Matlab and Mathematica codes are developed to numerically validate the analysis presented in Section 4.

Alternative Direction Sommerfeld Integral Representation. We have studied and validated the alternative direction Sommerfeld integral formulas for several layered media Green’s functions using Mathematica. Different direction plane wave representations of the free-space Green’s function presented in Section 4.2 can be readily found from existing literature [24, 25], therefore we focus on the results for the half-space layered medium with impedance boundary condition. The three layered medium Green’s functions can be handled in a very similar way, we skip the detailed formulas and interested readers are referred to the Mathematica files for these formulas as well as their validations.

For the half-space layered medium Green’s function, we focus on the evanescent part of the Green’s function given by

[TABLE]

To avoid the pole on the imaginary axis, we have numerically tested the following two contours.

In the first contour (left of fig. 2), a positive $c$ value is chosen so that

[TABLE]

As the first term integrates from [math] to $c$ on a finite (and reasonably small by proper choice of $c$) line segment (labeled $IV$), it can therefore be efficiently evaluated using standard Gauss quadrature rules. For the second integral, as the sum of the contour integrals on $I+II+III$ is [math] and the integral on $II$ approaches [math] when $R\to\infty$, the alternative direction representation of the integral

[TABLE]

is given by

[TABLE]

This representation is numerically validated using Mathematica’s NIntegrate with options AccuracyGoal$\to$20, PrecisionGoal$\to$20, WorkingPrecision$\to$60, MaxRecursion$\to$100, Method$\to$DoubleExponential for different $x,y>0$ values.

In the second contour (right of fig. 2), we assume $(r,\theta)$ are the polar coordinates of $(x,y)$, and perform the change of variable $u=\frac{t}{k}\sin\theta-i\frac{\sqrt{t^{2}+k^{2}}}{k}\cos\theta$, then the evanescent part becomes

[TABLE]

where the contour $I$ is a curve defined by $z(t)=\frac{t}{k}\sin\theta-i\frac{\sqrt{t^{2}+k^{2}}}{k}$, $t\in[0,\infty)$. As the integral on $II$ approaches [math] when $R\to\infty$, the evanescent part becomes the (negative) sum of the integrals on $III$ and $IV$ given by

[TABLE]

and

[TABLE]

where $\phi(u)=ik\cos\theta\sqrt{1+u^{2}}+uk\sin\theta$, $\psi(u)=(c+i\cos\theta)u-i\cos\theta$, and $c$ is a constant to be optimized so that both the Gauss quadrature applied on $IV$ and Laguerre quadrature on $III$ converge rapidly. These formulas are also validated using Mathematica for different $(x,y)$ and $c$ values. Note that in the second contour, the oscillatory term $e^{i\lambda x}$ or $e^{i\lambda y}$ is completely removed from the integrand, at the cost of a new rational function on a contour closer to the singularities.

We have tested numerically the convergence of the quadrature rules for the original and mathematically equivalent alternative direction representations of the half-space layered medium Green’s function. In fig. 3, we set $x=1$, $y=0.1$, $k=1$, and $c=2$, and present the numerical errors when different numbers of nodes are used in the quadrature rules for different integrals. The reference solutions are computed using Mathematica requesting more than $20$ correct digits. On the left of fig. 3, we present the accuracy when different numbers of Laguerre quadrature nodes are used for the original evanescent part in Eq. (42), the integral on $III$ of the first contour in Eq. (43), and the integral on $III$ of the second contours in Eq. (44). As $y\ll x$, the original integral converges slowly due to the oscillatory term $e^{i\sqrt{t^{2}+k^{2}}x}$. The Laguerre quadrature for the alternative direction integrals, on the other hand, converges much faster. For a fair comparison, we also present the convergence of the quadrature rules for the integrals on the finite line segments. When Legendre polynomial based Gauss quadrature is applied to the integral $IV$ on both the first and second contours, for double precision requirement, under the current setting of $c$, approximately $17$ nodes are required for the integral $IV$ on the first integral, while around $200$ Gauss nodes are required for the integral $IV$ on the second contour. In our numerical implementation, as the formulas for the first contour is easy to derive and manipulate, we therefore adopt the first contour for the evanescent part of a general layered media Green’s function.

Convergence of Multipole and Local Expansions. As we have discussed in Sec. 4.1, the convergence of the multipole and local expansions is determined by the ratio of $r/\rho$, where $r$ is the distance between the source (target) and center of the source (target) box in the multipole (local) expansion, and $\rho$ is the “modified” distance between the far-field target (source) and center of the source (target) box (constant $d$ and $\pm$ sign are considered). We have numerically validated the analysis. In fig. 4, we present the results for the half-space layered medium with impedance boundary conditions using settings $x_{target}-x_{c}^{s}=2$, $y_{target}+y_{c}^{s}=3$, $r=1.5$, $k=0.1$, $\alpha=1$, and the modified distance $\rho=\sqrt{(x_{target}-x_{c}^{s})^{2}+(y_{target}+y_{c}^{s})^{2}}$.

On the left of the figure, we compare the propagating and evanescent parts of the basis $\Phi_{p}$ for different $p$ values. Clearly, for large order $p$, the propagating part can be neglected when analyzing the truncation errors. On the right of the figure, we plot the ratio $(J_{p+1}\Phi_{p+1})/(J_{p}\Phi_{p})$ when $p\to\pm\infty$. When $|p|$ increases, the ratio approaches the constant $r/\rho$ (dashed green line).

A similar analysis is also performed for the three-layered medium Green’s function with wave numbers $k_{1}$, $k_{2}$, and $k_{3}$ in each layer. In the numerical experiment, we set $k_{1}=1$, $k_{2}=3$, $k_{3}=1$ and consider the contribution from upper layer ($y=0$)

[TABLE]

where

$\sigma_{2}^{t}(\lambda)=\frac{e^{d\sqrt{\lambda^{2}-k_{2}^{2}}}(\lambda^{2}+\sqrt{\lambda^{2}-k_{2}^{2}}\sqrt{\lambda^{2}-k_{3}^{2}}-k_{2}^{2})}{\sinh(d\sqrt{\lambda^{2}-k_{2}^{2}})(\lambda^{2}+\sqrt{\lambda^{2}-k_{1}^{2}}\sqrt{\lambda^{2}-k_{3}^{2}}-k_{2}^{2})+\sqrt{\lambda^{2}-k_{2}^{2}}(\sqrt{\lambda^{2}-k_{1}^{2}}+\sqrt{\lambda^{2}-k_{3}^{2}})\cosh(d\sqrt{\lambda^{2}-k_{2}^{2}})}.$

We study the convergence of the local expansion where the basis function $\Psi_{p}(x_{0},y_{0})$ is given by

[TABLE]

We set $x_{0}-x_{c}^{t}=2$, $y_{0}-y_{c}^{t}=3$, therefore the modified distance $\rho=\sqrt{2^{2}+3^{3}}$. The distance between the target and center of the target box is $r=1.5$. We neglect the propagating part in the layered medium Green’s function that is very small compared with the evanescent part for large orders, and only consider the integral from $k_{2}$ to $\infty$ in the evanescent part of the basis function $\Psi_{p}$. In fig. 5, we show the ratio $(J_{p+1}\Psi_{p+1})/(J_{p}\Psi_{p})$, which clearly converges to $r/\rho\approx 0.416$ when $p\to\infty$.

6 Conclusion and Generalization

In this paper, we present a general numerical framework for the efficient application of the layered media Green’s function to a given density function. Instead of constructing and compressing the matrix directly, which involves the expensive evaluations of one or more Sommerfeld type integrals for each matrix entry, the new algorithm considers a transformed matrix, so existing fast algorithms for the free-space Green’s function can be readily adapted for better algorithm efficiency. Theoretical analysis on the convergence of the new expansions and alternative direction Sommerfeld integral representations to accelerate the convergence of the numerical quadrature rules are provided and numerically validated. Similar to deriving the layered media Green’s functions, the detailed translations, alternative direction Sommerfeld integral representations, and number of terms in the expansions all depend on the geometric settings and physical parameters, especially in three dimensions. We have studied a few examples of such layered media Green’s functions in this paper, and we are working on both the analysis and implementation details for other important settings from domain applications. In particular, we are studying the optimal alternative direction Sommerfeld integral representations and more accurate estimate of the number of expansion terms in different scenarios. Results along these directions will be presented in the future.

Acknowledgement

J. Huang was supported by the NSF grant DMS1821093, and the work was finished while he was visiting professors at the King Abdullah University of Science and Technology, National Center for Theoretical Sciences (NCTS) in Taiwan, Mathematical Center for Interdisciplinary Research of Soochow University, and Institute for Mathematical Sciences of the National University of Singapore. M.H. Cho was supported by a grant from the Simons Foundation (No. 404499).