Fast Multipole Method For 3-D Helmholtz Equation In Layered Media

TL;DR

This paper introduces a fast multipole method tailored for efficiently computing wave interactions in 3-D layered media, combining classic FMM with new expansions for reaction fields, achieving near-linear complexity.

Contribution

It develops a novel FMM that integrates reaction field components for layered media using new multipole and local expansions, improving efficiency over traditional free space FMM.

Findings

Achieves $O(N)$ complexity for low-frequency wave interactions.

Demonstrates fast convergence of multipole expansions.

Validates efficiency through numerical experiments.

Abstract

In this paper, a fast multipole method (FMM) is proposed to compute long-range interactions of wave sources embedded in 3-D layered media. The layered media Green's function for the Helmholtz equation, which satisfies the transmission conditions at material interfaces, is decomposed into a free space component and four types of reaction field components arising from wave reflections and transmissions through the layered media. The proposed algorithm is a combination of the classic FMM for the free space component and FMMs specifically designed for the four types reaction components, made possible by new multipole expansions (MEs) and local expansions (LEs) as well as the multipole-to-local translation (M2L) operators for the reaction field components. { Moreover, equivalent polarization source can be defined for each reaction component based on the convergence analysis of its ME. The FMMs for the reaction components, implemented with the target particles and equivalent polarization sources, are found to be much more efficient than the classic FMM for the free space component due to the fact that the equivalent polarization sources and the target particles are always separated by a material interface.} As a result, the FMM algorithm developed for layered media has a similar computational cost as that for the free space. Numerical results validate the fast convergence of the MEs and the $O(N)$ complexity of the FMM for interactions of low-frequency wave sources in 3-D layered media.

Tables

Table 1. Table 1: Comparison of CPU time with multiple cores ( p = 4 𝑝 4 p=4 ).

cores	$N$	time for all ${\left\{{\Phi }_{\ell }^{free}\right\}}_{\ell =0}^2$	time for $\{{\Phi }_{00}^{\uparrow \uparrow },{\Phi }_{11}^{\uparrow \uparrow },{\Phi }_{22}^{\downarrow \downarrow }\}$
1	618256	46.01	4.60
	1128556	120.9	12.13
	1862568	282.4	13.88
	2861288	306.9	15.22
6	618256	8.75	2.77
	1128556	21.95	4.00
	1862568	52.19	4.57
	2861288	57.04	5.15
36	618256	2.29	2.74
	1128556	5.00	3.21
	1862568	11.45	3.59
	2861288	12.60	3.90

Equations

$$Y_n^m(\theta ,\phi )=(-1{)}^m\sqrt{\frac{2n+1}{4\pi }\frac{(n-m)!}{(n+m)!}}{P}_n^m(\cos \theta )e^{\mathrm{i}m\phi }={\hat{P}}_n^m(\cos \theta )e^{\mathrm{i}m\phi },$$

$$z_{\nu }(\omega )=j_{\nu }(\omega ),y_{\nu }(\omega ),h_{\nu }^{(1)}(\omega ),or{h}_{\nu }^{(2)}(\omega ).$$

$$z_0(k{r}_2)=\{\begin{array}{l}4\pi \sum _{n=0}^{\infty }\sum _{m=-n}^n(-1{)}^n{z}_n(kb)\overline{Y_n^m(\alpha ,\beta )}{j}_n(k{r}_1)Y_n^m({\theta }_1,{\phi }_1),\mathrm{if}{r}_1<b,\\ 4\pi \sum _{n=0}^{\infty }\sum _{m=-n}^n(-1{)}^n{j}_n(kb)\overline{Y_n^m(\alpha ,\beta )}{z}_n(k{r}_1)Y_n^m({\theta }_1,{\phi }_1),\mathrm{if}{r}_1>b.\end{array}$$

$$j_n(k{r}_2)Y_n^m({\theta }_2,{\phi }_2)=\sum _{\nu =0}^{\infty }\sum _{\mu =-\nu }^{\nu }{\hat{S}}_{n\nu }^{m\mu }(\mathbf{b})j_{\nu }(k{r}_1)Y_{\nu }^{\mu }({\theta }_1,{\phi }_1),$$

$${\hat{S}}_{n\nu }^{m\mu }(\mathbf{b})=4\pi (-1{)}^m{\mathrm{i}}^{\nu -n}\sum _{q=0}^{\infty }{\mathrm{i}}^q{j}_q(kb)\overline{Y_q^{\mu -m}(\alpha ,\beta )}\mathcal{G}(n,m;\nu ,-\mu ;q),$$

$$h_n^{(1)}(k{r}_2)Y_n^m({\theta }_2,{\phi }_2)=\sum _{\nu =0}^{\infty }\sum _{\mu =-\nu }^{\nu }{S}_{n\nu }^{m\mu }(\mathbf{b})j_{\nu }(k{r}_1)Y_{\nu }^{\mu }({\theta }_1,{\phi }_1),$$

$$h_n^{(1)}(k{r}_2)Y_n^m({\theta }_2,{\phi }_2)=\sum _{\nu =0}^{\infty }\sum _{\mu =-\nu }^{\nu }{\hat{S}}_{n\nu }^{m\mu }(\mathbf{b})h_{\nu }^{(1)}(k{r}_1)Y_{\nu }^{\mu }({\theta }_1,{\phi }_1),$$

$$\begin{array}{rl}{S}_{n\nu }^{m\mu }(\mathbf{b})& =4\pi (-1{)}^m{\mathrm{i}}^{\nu -n}\sum _{q=0}^{\infty }{\mathrm{i}}^q{h}_q^{(1)}(kb)\overline{Y_q^{\mu -m}(\alpha ,\beta )}\mathcal{G}(n,m;\nu ,-\mu ;q),\end{array}$$

$$\mathcal{G}(n,m;\nu ,\mu ;q)=(-1{)}^{m+\mu }\mathcal{S}\left(\begin{array}{ccc}n& \nu & q\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}n& \nu & q\\ m& \mu & -m-\mu \end{array}\right),$$

$$h_0^{(1)}(k|\mathbf{r}-{\mathbf{r}}^{\prime }|)=\sum _{|n|=0}^{\infty }\sum _{m=-n}^n{M}_{nm}{h}_n^{(1)}(k{r}_s)Y_n^m({\theta }_s,{\phi }_s),$$

$$h_0^{(1)}(k|\mathbf{r}-{\mathbf{r}}^{\prime }|)=\sum _{|n|=0}^{\infty }\sum _{m=-n}^n{L}_{nm}{j}_n(k{r}_t)Y_n^m({\theta }_t,{\phi }_t),$$

$$M_{nm}=4\pi j_n(k{r}_s^{\prime })\overline{Y_n^m({\theta }_s^{\prime },{\phi }_s^{\prime })},L_{nm}=4\pi h_n^{(1)}(k{r}_t^{\prime })\overline{Y_n^m({\theta }_t^{\prime },{\phi }_t^{\prime })},$$

$$L_{nm}=\sum _{|\nu |=0}^{\infty }\sum _{\mu =-\nu }^{\nu }{S}_{n\nu }^{m\mu }({\mathbf{r}}_c^t-{\mathbf{r}}_c^s)M_{\nu \mu }.$$

$$\begin{array}{rl}{\tilde{M}}_{nm}=& \sum _{\nu =0}^{\infty }\sum _{\mu =-\nu }^{\nu }\overline{{\hat{S}}_{n\nu }^{m\mu }({\mathbf{r}}_c^s-{\tilde{\mathbf{r}}}_c^s)}{M}_{\nu \mu },{\tilde{L}}_{nm}=\sum _{n=0}^{\infty }\sum _{\mu =-\nu }^{\nu }{\hat{S}}_{\nu n}^{\mu m}({\tilde{\mathbf{r}}}_c^t-{\mathbf{r}}_c^t)L_{\nu \mu },\end{array}$$

$${\tilde{M}}_{nm}=4\pi j_n(k{\tilde{r}}_s^{\prime })\overline{Y_n^m({\tilde{\theta }}_s^{\prime },{\tilde{\phi }}_s^{\prime })},{\tilde{L}}_{nm}=4\pi h_n^{(1)}(k{\tilde{r}}_t^{\prime })\overline{Y_n^m({\tilde{\theta }}_t^{\prime },{\tilde{\phi }}_t^{\prime })}$$

$$h_0^{(1)}(k|\mathbf{r}|)=\frac{1}{2k\pi }{\int }_0^{\infty }{\int }_0^{2\pi }{k}_{\rho }{e}^{\mathrm{i}{k}_{\rho }(x\cos \alpha +y\sin \alpha )}\frac{e^{\mathrm{i}{k}_z|z|}}{k_z}d\alpha d{k}_{\rho },$$

$$\begin{array}{r}{h}_0^{(1)}(k|\mathbf{r}-{\mathbf{r}}^{\prime }|)=\frac{1}{2k\pi }{\int }_0^{\infty }{\int }_0^{2\pi }{k}_{\rho }\frac{e^{\mathrm{i}\mathbf{k}\cdot (\mathbf{r}-{\mathbf{r}}_c^s)}{e}^{-\mathrm{i}\mathbf{k}\cdot ({\mathbf{r}}^{\prime }-{\mathbf{r}}_c^s)}}{k_z}d\alpha d{k}_{\rho },\\ h_0^{(1)}(k|\mathbf{r}-{\mathbf{r}}^{\prime }|)=\frac{1}{2k\pi }{\int }_0^{\infty }{\int }_0^{2\pi }{k}_{\rho }\frac{e^{\mathrm{i}\mathbf{k}\cdot (\mathbf{r}-{\mathbf{r}}_c^t)}{e}^{-\mathrm{i}\mathbf{k}\cdot ({\mathbf{r}}^{\prime }-{\mathbf{r}}_c^t)}}{k_z}d\alpha d{k}_{\rho }\end{array}$$

$$\begin{split}e^{{\rm i}\bm{k}\cdot{\bm{r}}}=&\sum\limits_{n=0}^{\infty}\sum\limits_{m=-n}^{n}4\pi{\rm i}^{n}j_{n}(kr)\overline{Y_{n}^{m}(\theta,\varphi)}\widehat{P}_{n}^{m}\Big{(}\frac{k_{z}}{k}\Big{)}e^{{\rm i}m\alpha}\\ =&\sum\limits_{n=0}^{\infty}\sum\limits_{m=-n}^{n}4\pi{\rm i}^{n}j_{n}(kr)Y_{n}^{m}(\theta,\varphi)\widehat{P}_{n}^{m}\Big{(}\frac{k_{z}}{k}\Big{)}e^{-{\rm i}m\alpha},\end{split}$$

$$e^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}}=e^{\mathrm{i}(\sqrt{k^2-k_z^2}(x\cos \alpha +y\sin \alpha )+k_{z}z)},$$

$${\hat{P}}_n^m(x)=\frac{1}{2^{n}n!}\sqrt{\frac{2n+1}{4\pi }\frac{(n-m)!}{(n+m)!}}(1-x^2{)}^{\frac{m}{2}}\frac{d^{n+m}}{d{x}^{n+m}}(x^2-1{)}^n.$$

$${\hat{P}}_n^m(z)=\frac{1}{2^{n}n!}\sqrt{\frac{2n+1}{4\pi }\frac{(n-m)!}{(n+m)!}}(1-z^2{)}^{\frac{m}{2}}\frac{d^{n+m}}{d{z}^{n+m}}(z^2-1{)}^n,z\in \mathbb{C},$$

$${\hat{P}}_n^m(x)=\sqrt{1-x^2}{Q}_{n-1}(x),$$

$${\hat{P}}_n^m(z)=\sqrt{1-z^2}{Q}_{n-1}(z),\forall z\in \mathbb{C},$$

$$z+1=r_1{e}^{\mathrm{i}{\theta }_1},-\pi <{\theta }_1\le \pi ,z-1=r_2{e}^{\mathrm{i}{\theta }_2},-\pi <{\theta }_2\le \pi ,$$

$$\sqrt{1-z^2}=-\mathrm{i}\sqrt{r_1{r}_2}{e}^{\mathrm{i}\frac{{\theta }_1+{\theta }_2}{2}}.$$

$$\underset{ϵ\to 0^{+}}{\lim }\sqrt{1-(x+\mathrm{i}ϵ{)}^2}=\sqrt{1-x^2},\underset{ϵ\to 0^{-}}{\lim }\sqrt{1-(x+\mathrm{i}ϵ{)}^2}=-\sqrt{1-x^2}.$$

$$\underset{ϵ\to 0^{+}}{\lim }{\hat{P}}_n^m(x+\mathrm{i}ϵ)={\hat{P}}_n^m(x),\underset{ϵ\to 0^{-}}{\lim }{\hat{P}}_n^m(x+\mathrm{i}ϵ)=-{\hat{P}}_n^m(x),\forall x\in [-1,1].$$

$${\hat{P}}_n^m(-z)=(-1{)}^{n+m}{\hat{P}}_n^m(z),{\hat{P}}_n^{-m}(z)=(-1{)}^m{\hat{P}}_n^m(z),$$

$$\begin{array}{r}{\hat{P}}_0^0(z)=\frac{1}{4\pi },{\hat{P}}_n^n(z)=-\sqrt{\frac{2n+1}{2n}}\mathrm{i}\sqrt{r_1{r}_2}{e}^{\mathrm{i}\frac{{\theta }_1+{\theta }_2}{2}}{\hat{P}}_{n-1}^{n-1}(z),n=1,2,\cdots ,\\ {\hat{P}}_{m+k}^m(z)=a_{m+k}^{m}z{\hat{P}}_{m+k-1}^m(z)-b_{m+k}^m{\hat{P}}_{m+k-2}^m(z),m=0,1,\cdots ,k=1,2,\cdots ,\end{array}$$

$$a_n^m=\sqrt{\frac{(2n-1)(2n+1)}{(n-m)(n+m)}},b_n^m=\sqrt{\frac{(2n+1)(n-m-1)(n+m-1)}{(2n-3)(n+m)(n-m)}}.$$

Full text

{AMS}

15A15, 15A09, 15A23

1 Introduction

The fast multipole method (FMM) has been a revolutionary development in modern computational algorithms for treating many-body interactions. In fact, it was considered one of the top 10 algorithms in the 20th century [10]. The FMM can reduce the $O(N^{2})$ cost of computing long-range interactions (Columbic electrostatics, wave scattering) among $N$ particles (or sources) to $O(N)$ or $O(N\log N)$. Such a capability of the FMM has had a tremendous impact on modern computational biology, astronomy, and computational acoustics and electromagnetics, among many other applications in sciences and engineering. The original FMMs developed by Greengard and Rokhlin [16, 17] for particles in the free space are based on multipole expansions (MEs) for a Green’s function $G(\bm{r},\bm{r}^{\prime})$=$G(\bm{r}-\bm{r}^{\prime})$ to achieve a low-rank representation for the long-range interactions between sources. The MEs for the wave interactions were made possible by the Graf’s addition theorem for Bessel functions. As a simple way to view the ME, we split the argument of the Green’s function, the zeroth order Bessel function for wave interactions, as $G(\bm{r}-\bm{r}^{\prime})=G((\bm{r}-\bm{r}_{c})+(\bm{r}_{c}-\bm{r}^{\prime}))$ where $\bm{r}_{c}$ is selected as the center of many sources $\bm{r}^{\prime}$ while $\bm{r}$ is any field location far away (well-separated) from all sources. The Graf’s addition theorem gives an expansion of $G$ in terms of separable terms involving ($\bm{r}-\bm{r}_{c}$) and ($\bm{r}_{c}-\bm{r}^{\prime}$), respectively. Terms involving ($\bm{r}-\bm{r}_{c}$) will appear as the higher order multipoles, i.e., higher order Bessel functions, while terms with ($\bm{r}_{c}-\bm{r}^{\prime}$) for each source $\bm{r}^{\prime}$ contributes to the ME expansion coefficients. For a far field location ${\bm{r},}$ such an expansion exhibits exponential convergence, thus only a small number of $p$ terms, therefore $p$ ME-coefficients, are needed, resulting in a $p$-term low-rank approximation for the far field of many sources.

As the original FMM was developed based on far field approximations by MEs obtained by the Graf’s addition theorem applied to free space Green’s functions, we can understand the difficulties of extending this approach to sources embedded in layered media, which are ubiquitous in computer engineering, geophysical, and medical image applications. For those applications, Green’s functions for layered media (layered Green’s functions) are preferred to describe the interactions to avoid introducing artificial unknowns on the infinite material interfaces. Unfortunately, in those cases, no theory like Graf’s addition theorem is available for layered Green’s functions. For this reason, the ME-based FMM of Greengard and Rokhlin has not been extended to layered Green’s functions due to the lack of corresponding multipole expansions of the far field of wave sources in general multi-layered media. In our recent work [9], MEs for the case of an impedance 2-D half space were obtained by using an analytical image representation of the layered Green’s function, which however are not available for general multi-layered media. Then, in our more recent work [38], MEs and multipole to local (M2L) translation operators were developed for 2-D Helmholtz equations in general layered media. In a different approach, an inhomogeneous plane wave fast method was developed [19] by approximating the Green’s function with plane waves sampled along a steepest descent path in the complex wave number space. So far, to handle the wave interaction of sources embedded in general layered media using layered Green’s functions, some other fast methods have been proposed such as FMM using Taylor expansion-based low-rank representation of Green’s function [32][31], and cylindrical wave decomposition of the Green’s function together with 2-D FMMs for cylindrical waves [8]. On the other hand, kernel independent compression techniques [37] [33] could also be considered for the layered Green’s functions.

In this paper, we will develop MEs and local expansions (LEs) for general layered Green’s function for 3-D Helmholtz equations as well as relevant M2L operators (thus extending our previous results for 2-D Helmholtz equations [38]), providing the key ingredients in the hierarchical design of the 3-D FMM. The layered Green’s function is decomposed into a free space part and four reaction fields arising from wave reflections and transmissions through the layered media. Our approach relies on two technical identities, the first one expresses the layered Green’s functions in terms of Sommerfeld integral involving plane waves, and the second one is the Funk-Hecke identity, which expresses plane waves in terms of cylindrical waves (Bessel functions). With the separable property of plane waves as well as the Funk-Hecke identity, we are able to derive the MEs and LEs for the layered Green’s functions in 3-D as well as the M2L operators. As shown in [38], the convergence of the far field for the reaction field component of the layered Green’s function in fact depends on the distance between the target and the location of some kind of polarization source, which can be defined for each reaction component (see (61) and Fig. 5). Therefore, in the implementation of the FMM for the reaction field components, the original and the polarization sources will be combined and embedded into a rectangular box, upon which the oct-tree structure will be built and the FMM can be implemented. As a result, the FMM for the free space can be extended straightforwardly to layered Green’s functions. Numerical results will show the fast convergence of the MEs and LEs for the layered Green’s function as well as the $O(N)$ efficiency for 3-D low-frequency acoustic wave interactions.

The rest of the paper is organized as follows. In section 2, we will re-derive the ME, LE and relevant M2L operators for the free space Green’s function using the new approach discussed above. The same technique will be then applied to layered Green’s functions. Section 3 gives the derivation of ME, LE and M2L operators for layered Green’s functions, given as Sommerfeld integrations of plane waves. A FMM based on the new ME, LE and M2L operator for reaction components is then presented. Section 4 gives numerical results for sources in three layers media. Various efficiency comparisons are given to show the performance of the proposed FMM. Finally, conclusion and discussions are given in Section 5.

2 A new derivation for $h$- and $j$-expansions of the free space Green’s function of 3-D Helmholtz equation and translations

In this section, we first review the $h$- and $j$-expansions, namely the multipole and local expansions, for the free space Green’s function of the Helmholtz equation and the translation from a $h$-expansion to a $j$-expansion. They are the key formulas in the free space FMM and can be derived by using the well-known addition theorems for regular and singular wave functions. Then, we will introduce a new derivation for the $h$- and $j$-expansions by using an integral representation of Hankel functions. This new technique will be applied to derive multipole and local expansions for reaction components of layered Green’s functions in the next section.

2.1 The $h$- and $j$-expansions of free space Green’s function

Let us first recall the addition theorems for regular and singular wave functions. Suppose $\bm{r}_{j}=(r_{j},\theta_{j},\varphi_{j})$ is the position vector of a point $P$ in spherical coordinates with respect to a given center $O_{j}$ for $j=1,2$, and $\bm{b}=(b,\alpha,\beta)$ is the position vector of $O_{1}$ with respect to $O_{2}$, such that $\bm{r}_{2}=\bm{r}_{1}+\bm{b}$. Meanwhile, the spherical harmonics are defined as

[TABLE]

where $0\leq|m|\leq n$ and $P_{n}^{m}(\cos\theta)$ is the associated Legendre function and $\widehat{P}_{n}^{m}(\cos\theta)$ is its normalized version. We will use the following addition theorems [24].

Theorem 2.1*.*

Let $z_{\nu}(\omega)$ be any spherical Bessel function of order $\nu$, that is

[TABLE]

Let $\bm{r}_{2}=\bm{r}_{1}+\bm{b}$. Then

[TABLE]

The special case of the above theorem when $z_{0}=j_{0}$ was proved by Clebsch in 1863, see [34, p.363]. According to Watson, another special case with $z_{0}=y_{0}$ was due to Gegenbauer.

Theorem 2.2*.*

Let $\bm{r}_{2}=\bm{r}_{1}+\bm{b}$. Then

[TABLE]

where

[TABLE]

with $\mathcal{G}(n,m;\nu,-\mu;q)$ being the Gaunt coefficient.

Theorem 2.3*.*

Let $\bm{r}_{2}=\bm{r}_{1}+\bm{b}$. Then

[TABLE]

for $r_{1}<b$, and

[TABLE]

for $r_{1}>b$, where $\widehat{S}_{n\nu}^{m\mu}(\bm{b})$ is given by (4) and

[TABLE]

and $\mathcal{G}(n,m;\nu,-\mu;q)$ is a Gaunt coefficient.

The Gaunt coefficient $\mathcal{G}(n,m;\nu,-\mu;q)$ is defined using the Wigner $3-j$ symbol as follows

[TABLE]

where $\mathcal{S}=\sqrt{(2n+1)(2\nu+1)(2q+1)/(4\pi)}$. Although we have explicit formulas (4) and (7) for the separation matrices $\widehat{S}_{n\nu}^{m\mu}(\bm{b})$ and ${S}_{n\nu}^{m\mu}(\bm{b})$, they are too complicated to be used directly for practical computations. Fortunately, recurrence formulas are available for their computations (cf. [6, 18]).

With these addition theorems, we can present the $h$- and $j$-expansions, which are also named as multipole and local expansions in the FMM. Consider the free space Green’s function of the Helmholtz equation with a source and a target at $\bm{r}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})$ and $\bm{r}=(x,y,z)$, respectively. By using the addition Theorem 2.1, we have an $h$-expansion (multipole expansion) with respect to a (source) center $\bm{r}_{c}^{s}$:

[TABLE]

and a $j$-expansion (local expansion) with respect to a (target) center $\bm{r}_{c}^{t}$:

[TABLE]

where

[TABLE]

$\bm{r}_{c}^{s}=(x_{c}^{s},y_{c}^{s},z_{c}^{s})$ is the source center close to $\bm{r}^{\prime}$, $\bm{r}_{c}^{t}=(x_{c}^{t},y_{c}^{t},z_{c}^{t})$ is the target center close to $\bm{r}$, $(r_{s},\theta_{s},\varphi_{s})$, $(r_{t},\theta_{t},\varphi_{t})$ are the spherical coordinates of $\bm{r}-\bm{r}_{c}^{s}$ and $\bm{r}-\bm{r}_{c}^{t}$, respectively, and $(r_{s}^{\prime},\theta_{s}^{\prime},\varphi_{s}^{\prime})$, $(r_{t}^{\prime},\theta_{t}^{\prime},\varphi_{t}^{\prime})$ are the spherical coordinates of $\bm{r}^{\prime}-\bm{r}_{c}^{s}$ and $\bm{r}^{\prime}-\bm{r}_{c}^{t}$, respectively (Fig. 1).

Applying addition Theorem 2.3 to $h_{n}^{(1)}(kr_{s})Y_{n}^{m}(\theta_{s},\varphi_{s})$ in (9), the translation from the $h$-expansion (9) to the $j$-expansion (10) is given by

[TABLE]

Similarly, we can shift the centers of $h$- and $j$-expansions via the following translations,

[TABLE]

where

[TABLE]

are the $h$- and $j$-expansion coefficients at new centers $\tilde{\bm{r}}_{c}^{s}$ and $\tilde{\bm{r}}_{c}^{t}$, respectively.

An important feature in (9)-(10) is that the source and target coordinates are separated, which will be important in deriving low rank approximations for far fields in the FMM (cf. [16, 17]). However, this target-source separation can also be achieved in the Fourier spectral domain. A new derivation via Fourier domain for (9) and (10) will be given by using an integral representation of $h_{0}^{(1)}(k|\bm{r}|)$. Moreover, this approach can also be used to derive multipole and local expansions for the reaction components of Green’s functions in layered media later.

2.2 A new derivation of $h$- and $j$-expansions

For a spherical wave, we have the well-known Sommerfeld identity

[TABLE]

where $k_{z}=\sqrt{k^{2}-k_{\rho}^{2}}$. It is necessary to point out that the identity (15) is for a lossless medium (i.e., $k$ is real) and needs to be considered as a limiting case of a lossy medium. The contour for the integral of $k_{\rho}$ is usually set by deforming the real axis to the fourth quadrant. In order to satisfy the outgoing wave radiation condition in the integrand, we have to ensure that $\mathfrak{Re}k_{z}>0$ and $\mathfrak{Im}k_{z}>0$. With the Sommerfeld identity (15), we can make a source-target separation in the spectral domain as follows

[TABLE]

for $z\geq z^{\prime}$ where $\bm{k}=(k_{\rho}\cos\alpha,k_{\rho}\sin\alpha,k_{z}).$ Without loss of generality, here we only consider the case $z\geq z^{\prime}$ for an illustration.

Next, we evoke the well-known Funk-Hecke formula [34, 24]

[TABLE]

which gives spherical harmonic expansions for plane waves and will be used for the plane waves inside (16). This classic Funk-Hecke formula only works for propagating plane wave, i.e., $0<k_{\rho}\leq k$ where all components of the propagation vector $\bm{k}$ are real. However, the Fourier spectral representations (16) involve not only propagating but also evanescent plane waves, i.e., $k_{\rho}>k$ where $k_{z}$ is purely imaginary. Therefore, for the spherical harmonic expansion of the evanescent plane wave, we need to extend the range of Funk-Hecke formula from $k_{z}\in\{\omega|0\leq\omega\leq k\}$ to $k_{z}\in\{\omega|0\leq\omega\leq k\}\cup\{\omega|\omega={\rm i}y,\,y\geq 0\}$. Fortunately, we can show that formula (17) holds for all $k_{z}\in\mathbb{C}$ by choosing an appropriate branch for the square root function. To see this, the left hand side of (17) can be written as

[TABLE]

while the right hand side only has $k_{z}$ as the argument of the normalized associated Legendre function. For the normalized associated Legendre function $\widehat{P}_{n}^{m}(x)$, we have $\widehat{P}_{n}^{-m}(x)=(-1)^{m}\widehat{P}_{n}^{m}(x),m>0$ and also the Rodrigues’ formula

[TABLE]

If $m$ is even, the analytic extension of $\widehat{P}_{n}^{m}(x)$ into the complex plane defined by

[TABLE]

is a polynomial of $z$ and hence is an entire function. For odd $m$, $\widehat{P}_{n}^{m}(x)$ can be written as

[TABLE]

where $Q_{n-1}(x)$ is a polynomial of at most $(n-1)$-th degree. Therefore, a complex extension defined by

[TABLE]

involves the multivalued function $\sqrt{1-z^{2}}$ similarly as in (18). In order to have an analytic continuation for the Funk-Hecke formula (17) for complex $k_{z}$, we need to choose an appropriate branch in (18) and (22). Let us consider polar forms

[TABLE]

where $\theta_{1}$ and $\theta_{2}$ are the principal values of the arguments of complex numbers $z+1$ and $z-1$. Then, we cut the complex plane from $-1$ to $+1$ along the real axis and choose the branch

[TABLE]

For any $x\in[-1,1]$, using this branch leads to

[TABLE]

The extensions $e^{{\rm i}\bm{k}\cdot\bm{r}}$ and $\widehat{P}_{n}^{m}(z)$ defined by using (23) enjoy the same property, e.g.,

[TABLE]

We note that the extension of associated Legendre function also satisfies

[TABLE]

and the recurrence formulas

[TABLE]

where

[TABLE]

The above recurrence formulas will be used in the implementation of the FMM.

Next, for the proof of the Funk-Hecke formula in the complex plane, we need the following lemmas.

Lemma 2.4*.*

For any real number $a\geq 0$, there holds

[TABLE]

where $j_{n}(a)=\sqrt{\frac{\pi}{2a}}J_{n+1/2}(a)$ is the spherical Bessel function of the first kind, $P_{n}(z)$ is the Legendre polynomial extended to the complex plane.

Proof 2.5*.*

Recall the series (cf. [27, 10.60.7])

[TABLE]

we can see that (28) holds for all $z\in[-1,1]$. Next, we consider its extension to the whole complex plane. Apparently, $e^{{\rm i}az}$ is an entire function of $z$. Meanwhile, the spherical Bessel function $j_{n}(a)$ has the following upper bound (cf. [1, 9.1.62])

[TABLE]

The Legendre polynomial $P_{n}(z)$ is a polynomial of degree $n$ with $n$ distinct roots $\{z_{j}\}_{j=1}^{n}$ in the interval $[-1,1]$. Therefore,

[TABLE]

here the estimate $a_{n}=\frac{(2n)!}{2^{n}(n!)^{2}}\leq 2^{n}$ for the coefficient of $z^{n}$ is used. These upper bounds for $j_{n}(a)$ and $P_{n}(z)$ give an estimate

[TABLE]

It is easy to show that the series on the righthand side of (28) converges uniformly in any compact set $D\subset\mathbb{C}$ and hence converges to an entire function of $z$. By the analytic extension theory, we have the proof.

By using the branch defined in (23) for the square roots, we have the extension of the well-known Legendre addition theorem [35, p.395].

Lemma 2.6*.*

Let $\bm{w}=(\sqrt{1-w^{2}}\cos\alpha,\sqrt{1-w^{2}}\sin\alpha,w)$ be a vector with complex entries, $\theta,\varphi$ be the azimuthal angle and polar angles of a unit vector $\hat{\bm{r}}$. Define

[TABLE]

then

[TABLE]

for all $w\in\mathbb{C}$.

Proof 2.7*.*

If $w\in[-1,1]$ is real, we have $\beta(w)=\cos\langle\bm{w},\hat{\bm{r}}\rangle$ and (34) is the well-known Legendre addition theorem. Define

[TABLE]

By using the branch defined in (23) for the square roots in $\beta(w)$ and $\widehat{P}_{n}^{m}(w)$, $P_{n}(\beta(w))$ and $g_{n}(w)$ are extended to the complex plane. Moreover, they are analytic in $\mathbb{C}\setminus[-1,1]$ and satisfy limit properties

[TABLE]

for all $x\in[-1,1]$. That implies $P_{n}(\beta(w))-g_{n}(w)$ takes angular boundary values zero on the set $[-1,1]$ (taking limit from the upper complex plane). By using Luzin-Privalov theorem [11], we conclude that $P_{n}(\beta(w))-g_{n}(w)=0$ in $\mathbb{C}\setminus[-1,1]$.

Proposition 2.8*.*

Given $\bm{r}=(x,y,z)\in\mathbb{R}^{3}$, $k>0$, $\alpha\in[0,2\pi)$ and denoted by $(r,\theta,\varphi)$ the spherical coordinates of $\bm{r}$, $\bm{k}=(\sqrt{k^{2}-k_{z}^{2}}\cos\alpha,\sqrt{k^{2}-k_{z}^{2}}\sin\alpha,k_{z})$ is a vector of complex entries. Choosing the branch (23) for $\sqrt{k^{2}-k_{z}^{2}}$ in $e^{{\rm i}\bm{k}\cdot{\bm{r}}}$ and $\widehat{P}_{n}^{m}(\frac{k_{z}}{k})$, then

[TABLE]

holds for all $k_{z}\in\mathbb{C}$, where

[TABLE]

Proof 2.9*.*

Define

[TABLE]

then, $kr\beta=\bm{k}\cdot\bm{r}$. Let $a=kr$, $z=\beta$ in (28), we have

[TABLE]

Finally, the extension of Funk-Hecke formula (36) follows by applying Lemma 2.6.

To show the convergence rate of the expansion (36), we consider three different cases with $\bm{r}=(0.3,0.4,z)$, $\alpha=\frac{\pi}{4}$. The convergence rates against $p$ are depicted in Fig. 2 for the relative $\ell_{2}$ error. Spectral convergence rates against $p$ are observed while the expansion converges more slowly if one of $|\bm{k}|$, $|\bm{r}|$ and $|k_{z}|$ gets larger.

Applying the spherical harmonic expansion (36) to exponential functions $e^{-{\rm i}\bm{k}\cdot(\bm{r}^{\prime}-\bm{r}_{c}^{s})}$ and $e^{{\rm i}\bm{k}\cdot(\bm{r}-\bm{r}_{c}^{t})}$ in (16) gives

[TABLE]

and

[TABLE]

for $z\geq z^{\prime}$, where $M_{nm}$ is defined in (11) and

[TABLE]

For the convergence of the Sommerfeld-type integrals in the above expansions, we only consider centers such that their $z$-coordinates satisfy $z_{c}^{s}<z$ and $z_{c}^{t}>z^{\prime}$. Recall the identity

[TABLE]

for $z\geq 0$, we see that (38) and (39) are exactly the $h$-expansion (9) and $j$-expansion (10) for the case of $z\geq z^{\prime}$.

To derive the translation from the $h$-expansion (38) to the $j$-expansion (39), we perform a further splitting in (38)

[TABLE]

and apply expansion (36) to obtain the following translation

[TABLE]

which implies an integral representation of $S_{n\nu}^{m\mu}(\bm{r}_{c}^{t}-\bm{r}_{c}^{s})$ in (12). In order to ensure the convergence of the Sommerfeld-type integral in the translation operator, the $z$-coordinates of the centers are also assumed to satisfy $z_{c}^{t}>z_{c}^{s}$.

3 FMM for 3-D Helmholtz equation in layered media

In this section, the multipole and local expansions for the reaction components of layered media Green’s function of 3-D Helmholtz equation will be derived by using the techniques introduced in the last section. Based on these expansions and relevant translation operators, FMM for 3-D Helmholtz Green’s function in layered media can be developed.

3.1 Green’s function of Helmholtz equation in layered media

Let us first review the integral representation of the layered Green’s function derived in [32]. Consider a layered medium consisting of $L$ interfaces located at $z=d_{\ell},\ell=0,1,\cdots,L-1$ as shown in Fig. 3.

A point source at $\bm{r}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})$ is located in the $\ell^{\prime}$-th layer ($d_{\ell^{\prime}}<z^{\prime}<d_{\ell^{\prime}-1}$). The layered Green’s function for the Helmholtz equation satisfies

[TABLE]

at a field point $\bm{r}=(x,y,z)$ in the $\ell$-th layer ($d_{\ell}<z<d_{\ell}-1$) where $\delta(\bm{r},\bm{r}^{\prime})$ is the Dirac delta function and $k_{\ell}$ is the wave number in the $\ell$-th layer. The system can be solved analytically in the Fourier $(k_{x},k_{y})$-domain for each layer in $z$ by imposing transmission conditions at the interface between $\ell$-th and $(\ell-1)$-th layer ($z=d_{\ell-1})$ as well as the decay conditions in the top and bottom-most layers as $z\rightarrow\pm\infty$ [8]. Inside a given layer, say, the $\ell$-th layer, the solution has the form

[TABLE]

where

[TABLE]

is called the reaction field in the $\ell$-th layer due to wave reflections and transmissions by the layered media. We can see that the reaction field $u_{\ell\ell^{\prime}}^{\rm r}(\bm{r},\bm{r}^{\prime})$ has up-going and down-going components inside intermediate layers $(0<\ell<L)$. Only up-going or down-going component is required in the top- and bottom-most layers, respectively. The up- and down-going components have Sommerfeld integral representations

[TABLE]

where $\bm{k}_{\alpha}=(k_{\rho}\cos\alpha,k_{\rho}\sin\alpha)=(k_{x},k_{y})$, $\bm{\rho}=(x,y)$, $\bm{\rho}^{\prime}=(x^{\prime},y^{\prime})$, $k_{\ell z}=\sqrt{k_{\ell}^{2}-k_{\rho}^{2}}$,

[TABLE]

The reaction densities $\sigma_{\ell\ell^{\prime}}^{\uparrow\uparrow}(k_{\rho}),\sigma_{\ell\ell^{\prime}}^{\uparrow\downarrow}(k_{\rho}),\sigma_{\ell\ell^{\prime}}^{\downarrow\uparrow}(k_{\rho}),\sigma_{\ell\ell^{\prime}}^{\downarrow\downarrow}(k_{\rho})$ only depend on the layer structure and wave numbers. Equation (46)-(47) are general formulas which are applicable to multi-layered media. Here, we give explicit formulas for the reaction densities of the three layers media with source in the middle layer, i.e.,

[TABLE]

where

[TABLE]

Detailed derivation of the general formulas (45)-(47) and corresponding reaction densities for Green’s function in layered media can be found in [32]. Although the density functions depend on the configuration of the layered media, they can be calculated at run time at any given point by solving one $2\times 2$ linear system and multiplying the results by some $2\times 2$ matrices.

3.2 Multipole and local expansions for general reaction component

Consider the reaction field in the middle layers, i.e.,

[TABLE]

Define

[TABLE]

By expressions in (46) and (47) we have a further decomposition

[TABLE]

and each component has a Sommerfeld-type integral representation:

[TABLE]

Here and in the rest of this paper, $*$ stands for any one of arrows $\uparrow$, $\downarrow$, e.g., $u_{\ell\ell^{\prime}}^{**}$ can be any of the four reaction components $u_{\ell\ell^{\prime}}^{\uparrow\uparrow},u_{\ell\ell^{\prime}}^{\uparrow\downarrow},u_{\ell\ell^{\prime}}^{\downarrow\uparrow},u_{\ell\ell^{\prime}}^{\downarrow\downarrow}$. Note that the source and target coordinates are only involved in the exponential functions $\mathcal{E}_{\ell\ell^{\prime}}^{**}(\bm{r},\bm{r}^{\prime})$. It is easy to make the following source-target separations

[TABLE]

by inserting the source center $\bm{r}_{c}^{s}$, and

[TABLE]

by inserting the target center $\bm{r}_{c}^{t}$. Here $\bm{\rho}_{c}^{s}=(x_{c}^{s},y_{c}^{s})$, $\bm{\rho}_{c}^{t}=(x_{c}^{t},y_{c}^{t})$ are the $x,y$-coordinates of the centers $\bm{r}_{c}^{s}$ and $\bm{r}_{c}^{t}$. Moreover, Funk-Hecke formula (36) gives spherical harmonic expansions for the following plane waves:

[TABLE]

where the fact that $Y_{n}^{m}(\pi-\theta,\varphi)=(-1)^{n+m}Y_{n}^{m}(\theta,\varphi)$ is used. Applying these expansions and source-target separation (52) and (53), we obtain

[TABLE]

and

[TABLE]

Substituting the above identities into (51), we obtain the following multipole expansion:

[TABLE]

with the expansion functions given by Sommerfeld-type integrals

[TABLE]

It is worthy to point out that $\mathcal{F}_{nm}^{**}(\bm{r},\bm{r}_{c}^{s})$ is only defined for $\bm{r}$ in the $\ell$-th layer. Therefore, it could be seen as singular function for $\bm{r}$ outside the $\ell$-th layer as reaction field produced by polarization charges there (more discussion on this issue will be given below). That is why we keep using the “multipole expansion” for expansion (54). Similarly, we obtain local expansion

[TABLE]

with coefficients given by

[TABLE]

According to the definition of $\mathcal{E}_{\ell\ell^{\prime}}^{**}(\bm{r},\bm{r}_{c}^{s})$ and $\mathcal{E}_{\ell\ell^{\prime}}^{**}(\bm{r}_{c}^{t},\bm{r}^{\prime})$ in (49), the centers $\bm{r}_{c}^{s}$ and $\bm{r}_{c}^{t}$ have to satisfy

[TABLE]

to ensure the exponential decay in $\mathcal{E}_{\ell\ell^{\prime}}^{**}(\bm{r},\bm{r}_{c}^{s})$ and $\mathcal{E}_{\ell\ell^{\prime}}^{**}(\bm{r}_{c}^{t},\bm{r}^{\prime})$ as $k_{\rho}\rightarrow\infty$ and hence the convergence of the corresponding Sommerfeld-type integrals in (55) and (57). These restriction can be met in practice, since we are considering targets in the $\ell$-th layer and sources in the $\ell^{\prime}$-th layer.

3.3 Reaction components, associated equivalent polarization sources, and multipole and local expansions

It is well known that the $h$-expansion (9) and the $j$-expansion (10) for the free space Green’s function have convergence rates of order \mathcal{O}\Big{(}\Big{(}\frac{|\bm{r}^{\prime}-\bm{r}_{c}^{s}|}{|\bm{r}-\bm{r}_{c}^{s}|}\Big{)}^{p}\Big{)} and \mathcal{O}\Big{(}\Big{(}\frac{|\bm{r}-\bm{r}_{c}^{t}|}{|\bm{r}_{c}^{t}-\bm{r}^{\prime}|}\Big{)}^{p}\Big{)}, respectively. That is both multipole and local expansion converge exponentially as the target moving away from the source. These convergence results are the key for the success of the hierarchical tree structure design in FMM. However, multipole and local expansions (54)-(57) for reaction components have a different convergence behavior. According to the convergence analysis for multipole and local expansions of 2-D Green’s function of Helmholtz equation in layered media (cf. [38]), we expect convergence estimates for the ME at $\bm{r}_{c}^{s}$ and the LE at $\bm{r}_{c}^{t}$ for the reaction components

[TABLE]

respectively, where

[TABLE]

Suppose $\bm{r}_{c}^{s}$ is a center close to source $\bm{r}^{\prime}$ with a fixed distance $|\bm{r}^{\prime}-\bm{r}_{c}^{s}|$, (59) indicates an important fact that the error of the truncated multipole expansion is not determined by the Euclidean distance between source center $\bm{r}_{c}^{s}$ and target $\bm{r}$ as in the free space case. Actually, the distances along $z$-direction have been replaced by summations of the distances between $\bm{r}$, $\bm{r}_{c}^{s}$ and corresponding nearest interfaces of the layered media. Similar conclusion can also be obtained for local expansion (56).

Remark 1*.*

There are two special cases, i.e., $d^{\uparrow\downarrow}(\bm{r},\bm{r}_{c}^{s})=|\bm{r}-\bm{r}_{c}^{s}|$ if $\bm{r}$ and $\bm{r}_{c}^{s}$ are in the $\ell$-th and $(\ell+1)$-th layer and $d^{\downarrow\uparrow}(\bm{r},\bm{r}_{c}^{s})=|\bm{r}-\bm{r}_{c}^{s}|$ if $\bm{r}$ and $\bm{r}_{c}^{s}$ are in the $\ell$-th and $(\ell-1)$-th layer. Therefore, multipole and local expansions of $u_{\ell\ell+1}^{\uparrow\downarrow}(\bm{r},\bm{r}_{c}^{s})$ and $u_{\ell\ell-1}^{\downarrow\uparrow}(\bm{r},\bm{r}_{c}^{s})$ have the same convergence behavior as that of free space Green’s function.

We will give some numerical examples to show the convergence behavior of the multipole expansions in (59). Consider the multipole expansions of $u_{11}^{\uparrow\uparrow}$ and $u_{11}^{\uparrow\downarrow}$ in a three-layer media with $k_{0}=0.8$, $k_{1}=1.5$, $k_{2}=2.0$, $d_{0}=0,d_{1}=-2.0$. In all the following examples, we fix $x=x^{\prime}=0.625$, $y=y^{\prime}=0.5$, $x_{c}^{s}=1.0,\;y_{c}^{s}=0.5$ for the coordinates of the target $\bm{r}$, source $\bm{r}^{\prime}$ and source center $\bm{r}_{c}^{s}$, respectively. Moreover, we will have the distance $|\bm{r}^{\prime}-\bm{r}_{c}^{s}|=0.375$ fixed by keeping $z^{\prime}=z_{c}^{s}$. For both components, we shall test three groups of $z$-coordinates given as follows

[TABLE]

The relative errors against truncation number $p$ are depicted in Fig. 4. We also use a linear polynomial to fit the value $a$ in the expected exponential convergence $\mathcal{O}(a^{p})$ via

[TABLE]

The fitted values of $a$ for all cases in Fig. 4 are given in the captions of the sub-figures. The results clearly show that multipole expansion of $u_{11}^{**}$ converges faster as the distance $d^{**}(\bm{r},\bm{r}^{\prime})$ rather than Euclidean distance $|\bm{r}-\bm{r}^{\prime}|$ increases. We also see that multipole expansion of $u_{11}^{\uparrow\downarrow}$ converges much faster than that of $u_{11}^{\uparrow\uparrow}$ as $d^{\uparrow\downarrow}(\bm{r},\bm{r}_{c}^{s})=d^{\uparrow\uparrow}(\bm{r},\bm{r}_{c}^{s})$. This can be explained by the extra exponential decay term inside the density $\sigma_{11}^{\uparrow\downarrow}(k_{\rho})$, see Fig. 7.

Polarization source for reaction field. In extending the FMM to the reaction field component of the layered Green’s function, the hierarchical tree structure design as in free space FMM relies on using the Euclidean distance between source and target to determine either direct computation or multipole and local expansions are used for the computation of source-target interactions. Therefore, multipole and local expansions as given in (54)-(57) are generally not compatible with the hierarchical design of FMM. Considering the convergence behavior of ME in (59) and (60), we introduce equivalent polarization sources for the four types of reaction fields (see. Fig. 5)

[TABLE]

which satisfy

[TABLE]

Also, we define

[TABLE]

and

[TABLE]

Recall the expressions (49), we can verify that

[TABLE]

Therefore, the reaction components can be re-expressed using equivalent polarization sources as follows

[TABLE]

Note that $\tilde{u}_{\ell\ell^{\prime}}^{\uparrow*}$ and $\tilde{u}_{\ell\ell^{\prime}}^{\downarrow*}$ have the same form as the two special cases $u_{\ell\ell+1}^{\uparrow\downarrow}$ and $u_{\ell\ell-1}^{\downarrow\uparrow}$, respectively, except for wave number $k_{\ell^{\prime}z}$ and densities $\sigma_{\ell\ell^{\prime}}^{**}$ in the $\ell^{\prime}$-th layer. Following the same derivation before, we obtain the following multipole expansion

[TABLE]

at equivalent polarization source centers $\bm{r}_{c}^{**}$ and local expansion

[TABLE]

at target center $\bm{r}_{c}^{t}$. Here, $(r_{s}^{**},\theta_{s}^{**},\varphi_{s}^{**})$ are spherical coordinates of $\bm{r}^{\prime}_{**}-\bm{r}_{c}^{**}$, $\widetilde{\mathcal{F}}_{nm}^{**}(\bm{r},\bm{r}_{c}^{**})$ are represented by Sommerfeld-type integrals

[TABLE]

and the local expansion coefficients are given by

[TABLE]

Similar to the restrictions in (58), the centers $r_{c}^{**}$ and $\bm{r}_{c}^{t}$ are required to satisfy

[TABLE]

Recall convergence results (59) and the fact (62), we conclude that ME (67) now satisfies

[TABLE]

As a result, the Euclidean distance $|\bm{r}-\bm{r}^{\prime}_{**}|$ can be used to determine if multipole expansions (67) are good approximations to the far reaction field. Similarly, the local expansion (68) for $\tilde{u}_{\ell\ell^{\prime}}^{**}(\bm{r},\bm{r}^{\prime}_{**})$ have a convergence of order $\mathcal{O}\left(\left(\frac{|\bm{r}-\bm{r}_{c}^{t}|}{|\bm{r}_{c}^{t}-\bm{r}^{\prime}_{**}|}\right)^{p}\right)$. These convergence results ensure that the hierarchical design can be applied in FMM with kernels $\tilde{u}_{\ell\ell^{\prime}}^{**}(\bm{r},\bm{r}^{\prime}_{**})$, ME (67) and LE (68).

Next, we discuss the center shifting and translation for ME (67) and LE (68). A desirable feature of the expansions of reaction components discussed above is that the formula (67) for the ME coefficients and the formula (68) for the LE have exactly the same form as the formulas of $h$-expansion coefficients and $j$-expansion for free space Green’s function. Therefore, we can see that center shifting for multipole and local expansions are exactly the same as free space case given in (13).

We only need to derive the translation operator from ME (67) to LE (68). Recall the definition of exponential functions in (64), $\widetilde{\mathcal{E}}_{\ell\ell^{\prime}}^{\uparrow}(\bm{r},\bm{r}_{c}^{\uparrow*})$ and $\widetilde{\mathcal{E}}_{\ell\ell^{\prime}}^{\downarrow}(\bm{r},\bm{r}_{c}^{\downarrow*})$ have the following splitting

[TABLE]

Applying Funk-Hecke formula (36) again, we obtain

[TABLE]

where $\mu=(n^{\prime}+m^{\prime})(\bmod\;2)$. Substituting into (67), the multipole expansion is translated to local expansion (68) via

[TABLE]

and the multipole-to-local translation operators are given in integral forms as follows

[TABLE]

where

[TABLE]

We note that the convergence of the Sommerfeld-type integrals in (74) is ensured by the conditions in (71).

Remark 2*.*

Interpretation of Polarization Sources: In special situations such as the half space (a two layer medium) with an impedance boundary condition on the interface, the reaction fields can be expressed in terms of those from point and line image charges located on the opposite side of the interface from the targets [9] [3]. However, in multilayered cases considered here, the reaction fields, given in terms of complicated integral expressions, in general can not be expressed in terms of explicit image charges. However, we can still introduce the polarization sources with locations given in (61) and Fig. 5, which will represent the effect of the reaction fields and are given locations based on the convergence behaviors of the MEs and LEs in (59) for the corresponding reaction field components. For the practical FMM implementation purpose, the locations of the polarization sources so defined enable us to decide whether the corresponding MEs can be used for the far field of the reaction components based on the distance between the far field location and the polarization sources, and thus, make the extension of FMMs to source interactions in layered media straightforward.

3.4 A FMM algorithm for sources in layered media

Let $\mathscr{P}_{\ell}=\{(Q_{\ell j},\bm{r}_{\ell j}),$ $j=1,2,\cdots,N_{\ell}\}$ be a group of source particles distributed in the $\ell$-th layer of a multi-layered medium with $L+1$ layers (see Fig. 3). The interactions between all $N:=N_{0}+N_{1}+\cdots+N_{L}$ particles are given by the summations

[TABLE]

for $i=1,2,\cdots,N_{\ell}$, $\ell=0,1,\cdots,L$, where

[TABLE]

Here, $u_{\ell\ell^{\prime}}^{**}(\bm{r},\bm{r}^{\prime})$ are the reaction components of the layered Green’s function in the $\ell$-th layer due to a point source in the $\ell^{\prime}$-th layer. We omit the factor $\frac{{\rm i}k_{\ell}}{4\pi}$ in $u_{\ell\ell^{\prime}}^{**}(\bm{r},\bm{r}^{\prime})$ to maintain consistency with the spherical Hankel function used for the free space component. General formulas for $u_{\ell\ell^{\prime}}^{**}(\bm{r},\bm{r}^{\prime})$ are given in (46)-(47).

Since the reaction components of Green’s function in multi-layer media have different expressions (46) for source and target particles in different layers, it is necessary to perform calculation individually for interactions between any two groups of particles among the $L+1$ groups $\{\mathscr{P}_{\ell}\}_{\ell=0}^{L}$. Here, we only consider the contribution from the reaction components and that from free space components will be calculated using free space FMM. Without loss of generality, let us focus on the computation of the $\uparrow\uparrow$ component in $\ell$-th layer due to sources in $\ell^{\prime}$-th layer, i.e., $\Phi_{\ell\ell^{\prime}}^{\uparrow\uparrow}(\bm{r}_{\ell i})$ for all $i=1,2,\cdots,N_{\ell}$ where we have source particles $\mathscr{P}_{\ell^{\prime}}$ and target particles $\mathscr{P}_{\ell}$. According to the discussion in the last section, we use the equivalent polarization sources (see. Fig. 6)

[TABLE]

and re-express $\Phi_{\ell\ell^{\prime}}^{\uparrow\uparrow}(\bm{r}_{\ell i})$ as

[TABLE]

A very important feature of the equivalent polarization sources $\{\bm{r}_{\ell^{\prime}j}^{\uparrow\uparrow}\}_{j=1}^{N_{\ell^{\prime}}}$ is that they are all separated from the corresponding targets $\{\bm{r}_{\ell i}\}_{i=1}^{N_{\ell}}$ by the interface $z=d_{\ell}$. As a matter of fact, they are located on different side of the interface $z=d_{\ell}$, see. Fig. 6 (right). Equivalent polarization sources (61) defined for other reaction components also have this property.

The framework of the free space FMM with ME (67), LE (68), M2L translation (73)-(74) and free space ME and LE center shifting allow us to use FMM for the computation of reaction component $\Phi_{\ell\ell^{\prime}}^{\uparrow\uparrow}(\bm{r}_{\ell i})$. FMM for other reaction components can be treated in the same manner. In the FMM for reaction components, for each reaction component, a corresponding large box is defined to include all equivalent polarization sources and target particles where the adaptive tree structure will be built by a bisection procedure, see. Fig. 6 (right). Note that the validity of the ME (67), LE (68) and M2L translation (73) used in the algorithm imposes restrictions (71) on the centers, accordingly. This can be ensured by setting the largest box for the specific reaction component to be equally divided by the interface between equivalent polarization sources and targets, see. Fig. 6 (right). Thus, the largest box for the FMM implementation will be different for different reaction component. With this setting, all source and target boxes of level higher than zeroth level in the adaptive tree structure will have centers below or above the interface $z=d_{\ell}$, accordingly. The fast multipole algorithm for the computation of the reaction component $\Phi_{\ell\ell^{\prime}}^{\uparrow\uparrow}(\bm{r}_{\ell i})$ is summarized in Algorithm 1.

All the interaction given by (75) will be obtained by calculating all components and summing them up.

3.5 Implementation details for efficiency

In this section, we will discuss some implementation details regarding the computation of double integrals involved in the multipole and local expansions and multipole-to-local translations. Especially, they can be simplified by using the following identity

[TABLE]

In particular, multipole expansion functions in (69) can be simplified as

[TABLE]

and the expression (70) for local expansion coefficients can be simplified as

[TABLE]

where $({\rho}_{s}^{**},\varphi_{s}^{**})$ and $({\rho}_{t}^{**},\varphi_{t}^{**})$ are polar coordinates of $\bm{r}^{\prime}_{**}-\bm{r}_{c}^{**}$ and $\bm{r}^{\prime}_{**}-\bm{r}_{c}^{t}$ projected on $xy$ plane. Moreover, the multipole to local translation (74) can be simplified as

[TABLE]

where $C_{n^{\prime}mm^{\prime}}^{**}=(-1)^{n^{\prime}}{\rm i}^{m^{\prime}-m}e^{{\rm i}(m^{\prime}-m)\varphi_{ts}^{**}}$ and $(\rho_{ts}^{**},{\varphi}_{ts}^{**})$ is the polar coordinates of $\bm{r}_{c}^{t}-\bm{r}_{c}^{**}$ projected on $xy$ plane. Denoting

[TABLE]

then we have

[TABLE]

The FMM demands efficient computation of Sommerfeld-type integrals $I_{nm,n^{\prime}m^{\prime}}^{**}$ defined in (80). It is clear that they have oscillatory integrands with pole singularities in $\sigma_{\ell\ell^{\prime}}^{**}(k_{\rho})$ due to the existence of surface waves. For a long time, much effort has been made on the computation of this type of integrals, including using ideas from high-frequency asymptotics, rational approximation, contour deformation (cf. [3, 4, 8, 26, 29]), complex images (cf.[15, 29, 25, 2]), and methods based on special functions (cf. [20]) or physical images (cf. [22, 23, 28, 21]). These integrals are convergent when the target and source particles are not exactly on the interfaces of a layered medium. Contour deformation with high order quadrature rules could be used for direct numerical computation at runtime. However, this becomes prohibitively expensive due to a large number of integrals needed in the FMM. In fact, $O(p^{4})$ integrals will be required for each source box to target box translation. Moreover, the involved integrand decays more slowly as the order of the involved associated Legendre function increases. The length of contour needs to be very long to obtain a required accuracy.

Note that $I_{nm,n^{\prime}m^{\prime}}^{**}(\rho,z,z^{\prime})$ is a smooth function with respect to $(\rho,z,z^{\prime})$. It is feasible to make a pre-computed table on a fine grid and then use interpolation to obtain approximations for the integrals. If we make pre-computed tables for all $|m|\leq n,|m^{\prime}|\leq n^{\prime}$, $n,n^{\prime}=0,1,\cdots,p$, there will be $(p+1)^{4}$ 3-D tables to be stored. However, the number of pre-computed tables can be reduced to $4(p+1)(2p+1)$.

Let

[TABLE]

then

[TABLE]

Recall the Rodrigues’ formula

[TABLE]

for $0\leq m\leq n$, and use the fact $k_{\ell z}=\sqrt{k_{\ell}^{2}-k_{\rho}^{2}}$ and (83), we obtain

[TABLE]

where

[TABLE]

Furthermore

[TABLE]

where $\mu=(n+m)(\bmod\;2)$, $\nu=(n^{\prime}+m^{\prime})(\bmod\;2)$,

[TABLE]

For $-n\leq m<0$ or $-n^{\prime}\leq m^{\prime}<0$, similar formulas can be obtained by using the fact $\widehat{P}_{n}^{m}(z)=(-1)^{m}\widehat{P}_{n}^{-m}(z)$ for $m<0$. In summary, we have

[TABLE]

for all $n,n^{\prime}=0,1,\cdots,$ and $-n\leq m\leq n$, $-n^{\prime}\leq m^{\prime}\leq n^{\prime}$ where $\mu=(n+|m|)(\bmod\;2)$, $\nu=(n^{\prime}+|m^{\prime}|)(\bmod\;2)$, r=\big{\lfloor}(n+|m|)/2\big{\rfloor},\quad r^{\prime}=\big{\lfloor}(n^{\prime}+|m^{\prime}|)/2\big{\rfloor}, and

[TABLE]

Here, we also use the notation $A_{nn^{\prime}mm^{\prime}}^{s}$ since it is equal to the coefficients defined in (88) for $m,m^{\prime}\geq 0$. Define integrals

[TABLE]

Then

[TABLE]

where $\mu=(n+|m|)(\bmod\;2)$, $\nu=(n^{\prime}+|m^{\prime}|)(\bmod\;2)$. We pre-compute integrals $\mathcal{S}_{nm,\mu\nu}^{**}(\rho,z,z^{\prime})$ on a 3-D grid $\{\rho_{i},z_{j},z_{k}^{\prime}\}$ in the domain of interest for all $n=0,1,\cdots,p$; $m=0,1,\cdots,2n+1$; $\mu,\nu=0,1$. Then, a polynomial interpolation is performed for the computation of integrals in the FMM.

To compute Sommerfeld-type integrals $\mathcal{S}_{nm,\mu\nu}^{**}(\rho,z,z^{\prime})$, it is typical to deform the integration contour by pushing it away from the real line into the fourth quadrant of the complex $k_{\rho}$-plane to avoid branch points and poles in the integrand, see Fig. 7. Here, we use a piecewise smooth contour which consists of two segments:

[TABLE]

We truncate $\Gamma_{2}$ at a point $t_{max}>0$, where the integrand has decayed to a user specified tolerance. As an example, we plot the density $\sigma_{11}^{**}(k_{\rho})$ along $\Gamma_{2}$ (see, Fig. 7 (c)). The three layers case with $k_{0}=1.5,k_{1}=0.8,k_{2}=2.0,d_{0}=0,d_{1}=-2.0$ and density given in Section 3 is used. We can see that the density $\sigma_{11}^{**}(k_{\rho})$ is bounded.

4 Numerical results

In this section, we present numerical results to demonstrate the performance of the proposed FMM for time harmonic wave scattering in layered media. This algorithm is implemented based on an open-source adaptive FMM package DASHMM [14] on a workstation with two Xeon E5-2699 v4 2.2 GHz processors (each has 22 cores) and 500GB RAM using the gcc compiler version 6.3.

We test the problem in three layers media with interfaces placed at $z_{0}=0$, $z_{1}=-1.2$. Particles are set to be uniformly distributed in irregular domains which are obtained by shifting the domain determined by $r=0.5-a+\frac{a}{8}(35\cos^{4}\theta-30\cos^{2}\theta+3)$ with $a=0.1,0.15,0.05$ to new centers $(0,0,0.6)$, $(0,0,-0.6)$ and $(0,0,-1.8)$, respectively (see Fig. 8 (a) for the cross section of the domains). All particles are generated by keeping the uniform distributed particles in a larger cube within corresponding irregular domains. In the layered media, the wave numbers are $k_{0}=1.2$, $k_{1}=1.5$, $k_{2}=1.8$. Let $\widetilde{\Phi}_{\ell}(\bm{r}_{\ell i})$ be the approximated values of $\Phi_{\ell}(\bm{r}_{\ell i})$ calculated by FMM. For accuracy test, we put $N=912+640+1296$ particles in the irregular domains in three layers see Fig. 8 (a). Convergence rates for the relative $\ell^{2}$ error ($Err_{2}$) and relative maximum error ($Err_{max}$) against $p$ are depicted in Fig. 8 (b). The CPU time for the computation of all three free space components $\{\Phi^{free}_{\ell}(\bm{r}_{\ell i})\}_{\ell=0}^{2}$, three selected reaction components $\{\Phi^{\uparrow\uparrow}_{00},\Phi^{\uparrow\uparrow}_{11},\Phi^{\downarrow\downarrow}_{22}\}$ and all sixteen reaction components $\Phi^{**}_{\ell\ell^{\prime}}(\bm{r}_{\ell i})$ with truncation $p=4$ are compared in Fig. 8 (c) for up to 3 millions particles. It shows that all of them have an $O(N)$ complexity while the CPU time for the computation of reaction components has a much smaller linear scaling constant due to the fact that most of the equivalent polarization sources are well-separated with the targets. CPU time with multiple cores is given in Table 1 and it shows that, due to the small amount of CPU time in computing the reaction components, the speedup of the parallel computing is mainly decided by the computation of the free space components. Here, we only use parallel implementation within the computation of each component, and all reaction components are computed independently. Therefore, it is straightforward to implement a version of the code which computes all components in parallel.

Pre-computed translation operator tables: Given the truncation order $p$, there will be $4(p+1)(2p+1)$ 3-D tables to be pre-computed, whose size depends on the required accuracy. In our numerical tests, we use tables with $30\times 30\times 30$ mesh with an interpolation polynomial of order $3$. Then, each table have $91^{3}$ integrals to be calculated. For each integral, the infinite interval is truncated at $t_{max}=100$ and a mesh with 51 sub-intervals and 30 Gaussian quadrature points in each sub-interval is used. It takes about 8.5 seconds on our workstation to compute one 3-D table as all the integrals can be calculated in parallel.

5 Conclusion and discussion

In this paper, we presented an $O(N)$ fast multipole method for the calculation of the discretized integral operator for the low-frequency Helmholtz equation in 3-D layered media. The layered media Green’s function is decomposed into a free space and four types of reaction field component(s). Using the spectral form of the layered Green’s functions and the Funk-Hecke identity, we developed new multipole expansion of $O(p^{2})$ terms for the far field of the reaction components, which can be associated with polarization sources at specific locations for the reaction field components. Multipole to local translation operators are also developed for the reaction fields. As a result, the traditional ME-based FMM can be applied to both the free space part and the reaction fields once the polarization sources are combined with the original sources. Due to the separation of the polarization and the original source charges by a material interface, the computational cost from the reaction field parts is only a fraction of that of the FMM for the free space part. For a given layered structure, besides the one time pre-computations of interpolation tables for the translation operator, computing the wave interactions of many sources in layered media costs the same as that for the wave interactions in the free space.

The complexity of the FMM on $L$ for a $L-$layer medium will scale as $O(L^{2})$ as the reaction component expression is different for different combination of source and target layers, the source-target interactions will need to be computed, separately. However, the FMMs for different group of reaction components can be done in parallel, and in practice $L$ is usually on the order of 10. Also, to reduce the pre-computation time of the tables for the M2L translation matrices, which involve the Sommerfeld integrals (90) with smooth integrands along the deformed contour, faster algorithms can be developed by using better quadrature rule (e.g., Clenshaw-Curtis-Filon-type method) and the recurrence formula of Bessel functions (see details in [39]).

The FMM developed here only considered low-frequency Helmholtz equations. Similar to the free space case, the multipole-expansion based FMM in this paper can not handle high frequency wave source interactions well as the number of terms in the MEs will increase in proportion to $k*D$ [12] where $k$ is the wave number and $D$ is the size of the scatterer. To address this difficulty, plane waves expansions could be introduced to yield diagonal forms for the translation operators in the FMM to produce an $O(N\log N)$ fast algorithm [30][7]. Though, the plane wave expansion may suffer a low-frequency breakdown when interaction between sources within sub-wavelength distance is computed for large objects with small fine structures. Therefore, it is important to develop a stable FMM applicable to a broad range of frequencies from low to high wave numbers without the low-frequency break down while still applicable to high frequency scattering. Various methods have been proposed toward this goal again for the free space case, including hybrid approaches combining multipole and plane wave expansions [5][36] as well as FMMs using stable plane wave expansions [13] and inhomogeneous plane waves [19].

As a future work, as the Green’s functions for the layered media are given in terms of plane waves via Sommerfeld integration, we will study broadband FMMs for the layered media. For an immediate task, we shall carry out error analysis for the new MEs and M2L operators for the reaction components for the Helmholtz equations in 3-D layered media, extending our results for the 2-D Helmholtz equations [38].

Acknowledgments

The authors thank Min Hyung Cho for helpful discussions.