A Heterogeneous FMM for 2-D Layered Media Helmholtz Equation I: Two & Three Layers Cases

TL;DR

This paper introduces a heterogeneous fast multipole method (H-FMM) for efficiently solving 2-D Helmholtz equations in layered media, focusing on two and three layers, with rigorous error analysis and promising numerical results.

Contribution

The paper develops a new H-FMM that leverages free-space multipole expansions for layered media, enabling easy adaptation and improved efficiency for multi-layered Helmholtz problems.

Findings

Demonstrates efficiency and accuracy of the H-FMM in numerical experiments.

Reduces computational cost for layered media Green's function evaluations.

Shows potential for extension to multi-layered media with minor modifications.

Abstract

In this paper, we will introduce a new heterogeneous fast multipole method (H-FMM) for 2-D Helmholtz equation in layered media. To illustrate the main algorithm ideas, we focus on the case of two and three layers in this work. The key compression step in the H-FMM is based on a fact that the multipole expansion for the sources of the free-space Green's function can be used also to compress the far field of the sources of the layered-media or domain Green's function, and a similar result exists for the translation operators for the multipole and local expansions. The mathematical error analysis is shown rigorously by an image representation of the Sommerfeld spectral form of the domain Green's function. As a result, in the H-FMM algorithm, both the "multipole-to-multipole" and "local-to-local" translation operators are the same as those in the free-space case, allowing easy adaptation of existing free-space FMM. All the spatially variant information of the domain Green's function are collected into the "multipole-to-local" translations and therefore the FMM becomes "heterogeneous". The compressed representation further reduces the cost of evaluating the domain Green's function when computing the local direct interactions. Preliminary numerical experiments are presented to demonstrate the efficiency and accuracy of the algorithm with much improved performance over some existing methods for inhomogeneous media. Furthermore, we also show that, due to the equivalence between the complex line image representation and Sommerfeld integral representation of layered media Green's function, the new algorithm can be generalized to multi-layered media with minor modification where details for compression formulas, translation operators, and bookkeeping strategies will be addressed in a subsequent paper.

Figures

Tables

Table 1. Table 1: Accuracy results with different expansion terms for k = 0.1 𝑘 0.1 {k}=0.1 and k = 1 𝑘 1 {k}=1 . Reference solution is computed with p = 39 𝑝 39 p=39 .

$p$	Error for $k=0.1$	Error for $k=1$
$E_5$	$1.23\times {10}^{-4}$	$1.43\times {10}^{-4}$
$E_{10}$	$2.73\times {10}^{-6}$	$3.81\times {10}^{-6}$
$E_{20}$	$2.06\times {10}^{-9}$	$2.85\times {10}^{-9}$
$E_{30}$	$1.19\times {10}^{-11}$	$1.65\times {10}^{-11}$

Table 2. Table 2: CPU time (seconds) for different N 𝑁 N using p = 39 𝑝 39 p=39 and k = 0.1 𝑘 0.1 {k}=0.1

$N$	100	6400	10000	90000	360000	640000	810000	1000000
CPU time	0.01	0.67	1.19	10.92	46.58	100.85	116.03	135.05

Equations

$$(\Delta +k^2)u(\mathbf{x})=0$$

$$\underset{r\to \infty }{\lim }\sqrt{r}\left(\frac{\partial }{\partial r}u(\mathbf{x})-iku(\mathbf{x})\right)=0,$$

$$g(\mathbf{x},{\mathbf{x}}_0)=\frac{i}{4}{H}_0^{(1)}(k||\mathbf{x}-{\mathbf{x}}_0||)$$

$$-(\Delta +k^2)g(\mathbf{x},{\mathbf{x}}_0)=\delta (\mathbf{x}-{\mathbf{x}}_0)$$

$$\underset{r\to \infty }{\lim }\sqrt{r}\left(\frac{\partial }{\partial r}g(\mathbf{x},{\mathbf{x}}_0)-ikg(\mathbf{x},{\mathbf{x}}_0)\right)=0,$$

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

$$g(\mathbf{x},{\mathbf{x}}_0)=$$

$$=$$

$$+\frac{1}{4\pi }{\int }_0^{\infty }\frac{e^{-t|y-y_0|}}{\sqrt{t^2+k^2}}\left(e^{i\sqrt{t^2+k^2}\left(x-x_0\right)}+e^{-i\sqrt{t^2+k^2}\left(x-x_0\right)}\right)dt$$

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

$$u_{{\mathbf{x}}_0}^s(\mathbf{x})={\int }_0^{\infty }g(\mathbf{x},{\mathbf{x}}_0^{im}-s\hat{y})\tau (s)ds,$$

$$\tau (s)=\delta (s)+\mu (s),s>0,$$

$$\mu (s)=2i\alpha e^{i\alpha \cdot s}.$$

$$u_{{\mathbf{x}}_0}^s(\mathbf{x})=g(\mathbf{x},{\mathbf{x}}_0^{im})+{\int }_0^{\infty }g(\mathbf{x},{\mathbf{x}}_0^{im}-s\hat{y})\mu (s)ds,$$

$$u_{{\mathbf{x}}_0}(\mathbf{x})$$

$$\equiv g(\mathbf{x},{\mathbf{x}}_0)+\left(g(\mathbf{x},{\mathbf{x}}_0^{im})+{\int }_0^{\infty }g(\mathbf{x},{\mathbf{x}}_0^{im}-s\hat{y})\mu (s)ds\right).$$

$${\int }_0^{\infty }g(\mathbf{x},{\mathbf{x}}_0^{im}-\eta \hat{y})e^{i\alpha \cdot s}ds$$

$$={\int }_0^{\infty }\left[\frac{1}{4\pi }{\int }_{-\infty }^{\infty }\frac{e^{-\sqrt{{\lambda }^2-k^2}|y+y_0+s|}}{\sqrt{{\lambda }^2-k^2}}{e}^{i\lambda (x-x_0)}d\lambda \right]e^{i\alpha \cdot s}ds$$

$$=\frac{1}{4\pi }{\int }_{-\infty }^{\infty }\frac{e^{-\sqrt{{\lambda }^2-k^2}\left(y+y_0\right)}{e}^{i\lambda (x-x_0)}}{\sqrt{{\lambda }^2-k^2}}\left[{\int }_0^{\infty }{e}^{-\sqrt{{\lambda }^2-k^2}s}{e}^{i\alpha \cdot s}ds\right]d\lambda$$

$$=\frac{1}{4\pi }{\int }_{-\infty }^{\infty }\frac{e^{-\sqrt{{\lambda }^2-k^2}\left(y+y_0\right)}{e}^{i\lambda (x-x_0)}}{\sqrt{{\lambda }^2-k^2}}\frac{1}{\sqrt{{\lambda }^2-k^2}-i\alpha }d\lambda .$$

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

$$\hat{\sigma }(\lambda )=\frac{\sqrt{{\lambda }^2-k^2}+i\alpha }{\sqrt{{\lambda }^2-k^2}-i\alpha },$$

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

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

$$u_2^t(\mathbf{x})$$

$$u_2^b(\mathbf{x})$$

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

$$\left(\begin{array}[c]{cccc}\frac{1}{\sqrt{\lambda^{2}-{k^{2}_{1}}}}&\frac{-1}{\sqrt{\lambda^{2}-{k^{2}_{2}}}}&-\frac{e^{-\sqrt{\lambda^{2}-{k^{2}_{2}}}d}}{\sqrt{\lambda^{2}-{k_{2}}^{2}}}&0\\ 0&\frac{e^{-\sqrt{\lambda^{2}-{k^{2}_{2}}}d}}{\sqrt{\lambda^{2}-{k^{2}_{2}}}}&\frac{1}{\sqrt{\lambda^{2}-{k^{2}_{2}}}}&\frac{-1}{\sqrt{\lambda^{2}-{k^{2}_{3}}}}\\ 1&1&-e^{-\sqrt{\lambda^{2}-{k^{2}_{2}}}d}&0\\ 0&e^{-\sqrt{\lambda^{2}-{k^{2}_{2}}}d}&-1&-1\end{array}\right)\left(\begin{array}[c]{c}\sigma_{1}(\lambda)\\ \sigma_{2}^{+}(\lambda)\\ \sigma_{2}^{-}(\lambda)\\ \sigma_{3}(\lambda)\end{array}\right)=\left(\begin{array}[c]{c}\frac{-1}{\sqrt{\lambda^{2}-{k^{2}_{1}}}}\\ 0\\ 1\\ 0\end{array}\right).$$

$$\varphi (\mathbf{x})=\int u_{{\mathbf{x}}_0}({\mathbf{x})\rho (\mathbf{x}}_0)d{\mathbf{x}}_0,$$

$$u(\mathbf{x})\approx \frac{i}{4}\sum _{p=-P}^P{\alpha }_p{H}_p(k|\mathbf{x}-{\mathbf{x}}_c|)e^{ip{\theta }_c},$$

$${\alpha }_p=\sum _{j=1}^N{q}_j{e}^{-ip{\theta }_j}{J}_p(k{\rho }_j),$$

$$u(\mathbf{x})=\sum _{j=1}^N{q}_j{u}_{{\mathbf{x}}_j}(\mathbf{x})=\sum _{j=1}^N{q}_{j}g(\mathbf{x},{\mathbf{x}}_j)+\sum _{j=1}^N{q}_j{u}_{{\mathbf{x}}_j}^s(\mathbf{x}).$$

$$u^s(\mathbf{x})=\sum _{j=1}^N{q}_j{u}_{{\mathbf{x}}_j}^s(\mathbf{x})$$

$$=\frac{i}{4}\sum _{j=1}^N{q}_j\left(H_0^{(1)}(k|\mathbf{x}-{\mathbf{x}}_j^{im}|)+{\int }_0^{\infty }{H}_0^{(1)}(k|\mathbf{x}-({\mathbf{x}}_j^{im}-s\hat{y})|)\mu (s)ds\right)$$

Full text

{AMS}

65R20, 65Z05, 78M25

1 Introduction

To compute the interactions of the electromagnetic or acoustic waves with objects of complex geometry embedded in the multi-layered media, an attractive numerical method in the engineering community is to reformulate the frequency domain Helmholtz equation as a boundary integral equation (BIE) using the layered-media Green’s function where the unknowns are only defined on the surface of the objects. Unlike the translation-invariant free-space Green’s function for wave scattering in homogeneous media, the layered-media Green’s function incorporates the interface and far-field boundary conditions and becomes a spatially variant function. In this work, we will refer to the layered-media Green’s function as the domain Green’s function or domain kernel function. Subsequently, the boundary integral equations are discretized using proper numerical integration techniques, for instance, the trapezoidal rule with end point corrections [3, 27] in two dimensions or the Quadrature by Expansion (QBX) technique in higher dimensions [29], resulting in a dense linear system where the matrix describes how the discretized source and target particles interact through the domain Green’s function. In the numerical solver step, an important building block is the efficient application of this matrix to a given vector representing the source contributions, to derive the potential field due to the domain Green’s function interactions at the target particle locations.

There exist several strategies to compute integral operators of the spatially variant domain Green’s function efficiently. For simple geometries, for instance, the half-space or spheres, one technique is to represent the domain Green’s function contribution as the sum of the free-space Green’s function contributions from both the original source and some image points (the spatial variant properties are incorporated into the locations of the images). This approximation allows the direct application of existing free-space fast matrix vector multiplication algorithms specially designed for the free-space translation-invariant kernels, including the well-developed fast multipole method (FMM) packages in [10, 17, 19, 37]. Representing works along this direction include the classical Kelvin image for the half-space problem (or spheres) for the perfect conducting media [32], so the spatially variant domain Green’s function simply consists of two free-space Coulomb potentials, one from the source charge and one from its image. In the case of dielectric inhomogeneity, such as a spherical cavity embedded in a dielectric medium, the reaction field from the media can also be approximated by a small number of image charges [7][6]. Unfortunately, for more complex geometry, the image approximations are extremely hard to derive or non-existent, and for a few special cases including the multi-layered media, the domain Green’s functions are customarily derived as Sommerfeld integral formulas using integral transformations. Even when the image approximation formulas are available (e.g., the two-layered media Helmholtz equation), a large number of images is usually required. In [34], to approximate the interaction of $6,400$ particles described by the domain Green’s function of the 2D Helmholtz equation with half-space impedance boundary condition, a total of $1,122,960$ additional images were introduced in a hybrid approach, which combines the image and Sommerfeld integral representations. Other efforts to speed up the computation of integral operator for layered media Green’s functions include the inhomogeneuous plane wave method [26] and the cylindrical wave decomposition of the Green’s function in 3-D and the 2D-FMM [11].

Another approach is to compress the matrix describing the domain Green’s function interactions directly using the fast direct solvers (FDS) [20, 25] or closely related $\mathcal{H}$-matrix theory [22, 24], where the low-rank structures of the sub-matrices are derived and processed recursively on a hierarchical tree structure using purely numerical linear algebra techniques. However, evaluating the domain Green’s functions (entries in the matrix) involves very expensive computation of the Sommerfeld type integrals, and the compression stage of the FDS is expensive and memory intensive. It is worth mentioning that the FMM, FDS, and $\mathcal{H}$-matrix are all hierarchical algorithms that recursively compress the information in a system to low-rank or low-dimensional forms and transmit the compressed information non-locally on a hierarchical tree structure. In this paper, we apply this hierarchical algorithm design philosophy to multi-layered media domain Green’s functions and present a new hierarchical algorithm for evaluating the spatially variant domain Green’ s function interactions. Our algorithm shares many common features with FMM and FDS algorithms, especially in the information transmission patterns on the tree structure: the compressed representations are transmitted through an upward pass from leaf to parent nodes on a hierarchical tree structure, collected by interacting nodes and stored as “local expansions”, and then transmitted to the children nodes in a downward pass.

We start from the 2-D half-space (two-layered) problem where the domain Green’s function can be explicitly represented with the Sommerfeld integrals or complex line images. The “complex line image” representation intuitively reveals how the compression of the interaction matrix can be performed analytically on a “transformed” matrix which only involves the free-space Green’s function, and provides rigorous error analysis using available analytical results from the classical free-space FMM. As the compressed representation separates the spatially variant components and spatially invariant free-space kernels in the domain Green’s function, both the “multipole-to-multipole” and “local-to-local” translations from the existing free-space Helmholtz FMM algorithm can be easily adapted. Unlike the classical spatially invariant FMM algorithm, all the spatially variant information are collected in the “multipole-to-local” translations and the new algorithm becomes spatially variant. We refer to this new algorithm as the Heterogeneous FMM due to the heterogeneous nature of the “multipole-to-local” translations and the use of the free-space Green’s function and similar translations on the hierarchical tree structure from the classical FMM. We present the algorithm structure and demonstrate its accuracy and efficiency by comparing with the hybrid method in Ref. [34] for handling inhomogeneous media. More interestingly, by relating the results from the complex line image representation to the equivalent Sommerfeld integral representation, we demonstrate how the algorithm can be generalized to a particular three-layered media setting where the complex image representation becomes too complicate to derive and only the Sommerfeld integral representation via integral transformations [33] is available.

This paper is organized as follows. In Sec. 2, we present both the complex line image and Sommerfeld integral representations of the free-space and domain Green’s functions for the 2-D Helmholtz equation with half-space impedance boundary condition, and present the Sommerfeld integral representation of the domain Green’s function for a particular three-layered media setting where the analytic image formula is not available. In Sec. 3, we focus on the half-space impedance boundary condition problem, and present the hierarchical algorithm for the efficient evaluation of the spatially variant domain Green’s function interactions. We will discuss the hierarchical tree structure, compression of the complex image representations as “multipole expansions”, compression of the local interactions to allow more efficient evaluations of the integrals, adaptation of the spatially invariant “multipole-to-multipole” and “local-to-local” translations from existing free-space FMM, analytical formulas for the heterogeneous “multipole-to-local” translations and their efficient evaluations, and present the algorithm structure and some implementation details. Preliminary numerical results are presented in Sec. 4 to demonstrate the algorithm accuracy and efficiency. Although it is easier to describe the two-layered media algorithm and perform the error analysis purely using the complex line image representation in Sec. 3, we also present the mathematically equivalent results for the Sommerfeld integral representation. Comparing these two mathematically equivalent representations, we introduce a numerical framework solely based on the Sommerfeld integral representation in Sec. 5 for a particular three-layered media setting, which allows direct extension to multi-layered media. Finally, we summarize our results in Sec. 6.

2 2-D Helmholtz Equation in Multi-layered Media

We present both the complex line image and Sommerfeld integral representations of the domain Green’s functions for the 2-D Helmholtz equation in half-space with impedance boundary condition, and the Sommerfeld integral representation for a particular three-layered media setting where an analytical expression for the image representation may not be available.

2.1 Free-space Green’s function

Consider the 2-D Helmholtz equation in free-space

[TABLE]

with the Sommerfeld radiation condition at $\infty$

[TABLE]

where $\mathbf{x}=(x,y)$ and $r=||\mathbf{x}||$. Its Green’s function is given by the $0^{th}$ order Hankel function of the first kind as

[TABLE]

which solves the equation

[TABLE]

with the Sommerfeld radiation condition

[TABLE]

where $\delta$ is the 2-D Dirac delta function, $\mathbf{x}_{0}=(x_{0},y_{0})$, $r=||\mathbf{x}-\mathbf{x}_{0}||$, ${k}$ is the wave number, and $i=\sqrt{-1}$.

The free-space Green’s function can be found in the frequency (spectral) domain by taking the Fourier transform of Eq. (2) in the $x$-direction and solving the resulting ODE in the $y$-direction to give its spectral representation

[TABLE]

This representation is often referred to as the Sommerfeld identity, which can be separated into the propagating and evanescent modes for wave number variable $|\lambda|<k$ (propagating modes) and $|\lambda|>k$ (evanescent modes as $|y|\rightarrow\infty$), respectively, to arrive at the following form after some changes of variables

[TABLE]

for $|y-y_{0}|>0$.

The free-space Green’s function is commonly used in the potential theory, where solutions of the Helmholtz equation are represented as combinations of volume and/or layer potentials defined as the convolution of the Green’s function or its derivatives with certain density functions either over the volume or surface area of a given object. Theoretical properties of the free-space Green’s function and corresponding potentials are also well-established in existing literature [14, 30], and their efficient evaluations can be carried out using fast algorithms such as the well-developed wide-band fast multipole method [8, 10].

2.2 Domain Green’s function for two-layered media

In layered media, it is usually possible to derive the spatially variant domain Green’s function analytically either using the method of images (complex image representation) or applying the integral transforms (e.g., Laplace and Fourier transforms) to derive the spectral domain representation (Sommerfeld integral representation). In this subsection, we focus on the 2-D half-space Helmholtz equation with the impedance boundary condition

[TABLE]

which is imposed on the interface defined by $y=0$, and present the complex image and Sommerfeld integral representations from existing literature (e.g., see [28]).

Complex Image Representation. We first present the image representation of the domain Green’s function. The domain Green’s function for the interaction of a point source located at $\mathbf{x}_{0}$ with a target point at $\mathbf{x}$ is usually decomposed as the sum of the free-space interaction of the source and target points and contribution from a scattered field $u_{\mathbf{x}_{0}}^{s}(\mathbf{x})$ which can be explicitly represented in the two-layered media as complex image contributions of the free-space kernel as

[TABLE]

where $\mathbf{x}_{0}^{im}=(x_{0},-y_{0})$, $\hat{y}=(0,1)$, and $\tau(s)$ is the complex image charge density distribution. By applying the impedance boundary condition, the image function $\tau(s)$ can be explicitly found by (see [34])

[TABLE]

where a point image is indicated by the Dirac delta distribution $\delta(s)$ and a line image $\mu(s)$ is complex and

[TABLE]

As a result, we have

[TABLE]

where the first term on the right hand side represents the contribution from the point-image source, and the second term represents the contributions from the line-images. Therefore, the domain Green’s function $u_{\mathbf{x}_{0}}(\mathbf{x})$ for the half-space Helmholtz equation with impedance boundary condition can be represented in terms of the free-space Green’s function $g(\mathbf{x},\mathbf{x}_{0})$ as

[TABLE]

Sommerfeld Integral Representation. The scattered field in (9) involves an integration of an oscillatory line image density $\mu(s)=2i\alpha e^{i\alpha\cdot s}$, which cannot be handled efficiently with numerical quadratures directly as in the case for the Laplace equation in [7]. However, using the Sommerfeld identity for $g(\mathbf{x},\mathbf{x}_{0})$ in Eq. (3), we can resolve this difficulty with an analytic integration of the $s$ variable as follows:

[TABLE]

Plugging Eqs. (11) and (3) into Eq. (10), we obtain the following spectral domain representation for the scattered field (assume $y>0$)

[TABLE]

or by defining

[TABLE]

we have

[TABLE]

where $\hat{\sigma}(\lambda)$ is independent of $\mathbf{x}$ and $\mathbf{x}_{0}$.

2.3 Domain Green’s function for multi-layered media

For multi-layered media, the explicit forms of the complex image representations are in general unavailable and the domain Green’s functions are customarily expressed in terms of the Sommerfeld integrals [9, 13, 16, 33]. We leave the detailed derivation of the Sommerfeld integral type domain Green’s functions for different multi-layered media settings to a subsequent paper. In the following, we cite the results from Ref. [31] for a particular three-layered media setting, and compare its Sommerfeld integral representation of the domain Green’s function with that of the two-layered media in Eq. (14).

Assume a point source is located at $\mathbf{x}_{0}=(x_{0},y_{0})$ in the top layer with a wave number $k_{1}$ and the two interfaces of the three-layered media are located at $y=0$ and $y=-d$, respectively. The Sommerfeld integral representation for the scattered field in the top layer ($y>0$) can be represented as

[TABLE]

where the unknown function $\sigma_{1}(\lambda)$ will be determined later. In the middle layer with a wave number $k_{2}$, the scattered field $u_{2}^{s}$ can be written as the sum of the contributions $u_{2}^{t}$ (from upper interface) and $u_{2}^{b}$ (from lower interface) as

[TABLE]

and in the bottom layer with wave number $k_{3}$, we have

[TABLE]

where $\sigma_{2}^{+}(\lambda)$, $\sigma_{2}^{-}(\lambda)$, and $\sigma_{3}(\lambda)$ are unknown quantities which are associated with layer reflection and transmission coefficients of waves in spectral domain. When the interface conditions are given by $[u]=0$ and $[\frac{\partial u}{\partial n}]=0$, these quantities can be determined by solving the linear system

[TABLE]

When all the source and target points are located in the top layer, we see that the domain Green’s function $u_{1}^{s}$ is similar to Eq. (14) but with a different $\sigma$ function, the nature of which will be revealed later.

One additional complexity of the multi-layered media computation is that the source and target points may be located in different layers. We leave the technical details for different domain Green’s functions and their compressions, translations, and bookkeeping strategies for such cases in a subsequent paper, and only focus on the Green’s function in the form of Eq. (15) in this paper. This simplifies the notations while presenting the main idea of the new heterogeneous FMM algorithm, and still covers important applications for scattering of multiple objects over conducting surface (the two-layered media with impedance boundary condition) and layered media (meta-surfaces over layered substrate in material sciences).

2.4 Domain Green’s function in integral equation methods

In most integral equation formulations of the Helmholtz equation, unlike the translation invariant free-space Green’s function that only depends on the distance of $\mathbf{x}$ and $\mathbf{x}_{0}$, the domain Green’s function for general complex geometry brings complications for being a two variable function and its values are simply no longer translation invariant but spatially variant. As a result, the computation and evaluation of the domain Green’s function are more expensive than finding the solution of the original differential equation. There are a few exceptions, including the simulation of the layered-media Helmholtz equation where the interface or boundary is infinite but regular. For such cases, if using the free-space Green’s function, the resulting potentials will involve the evaluations of integrals (potentials) over infinite interfaces. However, as the geometry is regular, one can analytically derive the spatially variant domain Green’s function, in the form of a Sommerfeld integral representation using integral transformations such as the Laplace and Fourier transforms.

There are many advantages by using the domain Green’s function in forming the integral equation method (IEM) for the multi-layered media problem, for instance, the interface conditions are naturally enforced by the domain Green’s function and no unknowns are necessary on the layer interfaces. However, the numerical solution of the integral equation poses many challenges and is still an active research topic. In addition to problems common to all integral equation approaches such as the design of high order quadrature and derivation of well-conditioned systems, the IEM for layered media using the domain Green’s function has its specific challenges. In particular, the evaluation of the domain Green’s function interactions with large number of source and target points is expensive for either the complex image representation, or the Sommerfeld integral, or even the optimized hybrid representations. This implies that explicitly constructing the discretized interaction matrix is also expensive, and therefore matrix compression using the FDS will be costly where purely numerical linear algebra techniques are applied. This paper focuses on the fast application of the domain Green’s function to a given density function $\rho(\mathbf{x}_{0})$ as in

[TABLE]

where the integral either represents a volume potential or a surface layer potential. After discretization, the resulting linear algebra question becomes how to efficiently calculate the matrix-vector multiplication of $A\mathbf{v}$ where entries in the matrix $A$ are given by the domain Green’s function $A_{i,j}=u_{\mathbf{x}_{j}}(\mathbf{x}_{i})$. The main results of this paper include (a) the analysis-based low-rank compression of the matrix $A$, which is not directly performed on the matrix itself as in the FDS methods, but on a closely related matrix after certain transformations; (b) how the compressed representations can be transmitted through the hierarchical tree structure using analysis-based translation operators, and whenever possible, utilizing existing translation operators for the free-space kernels; and (c) the selected compression and translation strategies allow the implementation of a “heterogeneous” FMM algorithm for the layered media by an easy adaptation of existing free-space FMM codes.

3 Algorithm for Two-layered Media

In this section, we present the technical details of a fast hierarchical algorithm for the two-layered media domain Green’s function. The algorithm is similar in structure to that of FMM and the compression stage of FDS, and is developed by considering the design philosophy of the hierarchical modeling technique. This technique identifies any low-rank, or low-dimensional, or other compact features in a given system, recursively collects the compressed representations from children to parents, and transmits the information between different nodes on a hierarchical tree structure using properly compressed translation operators. It is worth mentioning that the resulting hierarchical models are often re-expressed as recursive algorithms, which can be easily interfaced with existing dynamical schedulers from High-Performance Computing (HPC) community for optimal parallel efficiency [1, 4, 18].

In addition to FMM and FDS, different aspects of the hierarchical modeling technique have been known and addressed by many researchers previously. Examples include the classical fast Fourier transform (FFT) [15] where the Halving Lemma shows how data can be compressed and the odd-even term splitting of the polynomials creates a hierarchical tree to allow recursively processing the compressed information efficiently; the multigrid method (MG) [5, 23] where the hierarchical tree structure is formed via adaptively refining the computational domain, and data compression and transmission are performed using the relaxation (smoother) and projection (restriction) operators by analyzing the frequency domain behaviors of the error functions between different levels of the (adaptive) tree to effectively reduce the high frequency errors. When there are $n$ terms (FFT) in the polynomial or $n$ approximately uniformly distributed particles (MG or FMM), the depth of the hierarchical tree is normally $\mathcal{O}(\log n)$ and the number of tree nodes is approximately $\mathcal{O}(n)$. Therefore, if each level only requires $\mathcal{O}(n)$ operations (e.g., FFT), the algorithm complexity will be $\mathcal{O}(n\log n)$. If each tree node only requires a constant amount of operations (e.g. MG or FMM), the algorithm complexity will be asymptotically optimal $\mathcal{O}(n)$. In this section, we describe our algorithm following the design guidelines of the hierarchical modeling technique.

3.1 Adaptive hierarchical tree structure

Consider the Helmholtz equation in 2-D with half-space impedance boundary condition for the scattering of a finite-sized object with a complex geometry in the upper half plane $y>0$. A surface integral equation can be derived to give the scattering solution as layer potentials through a convolution of the domain Green’s function or its derivatives with some unknown density functions over the object’s surface. We assume the surface is “discretized” into a number of particles via proper numerical integration techniques. In the hierarchical modeling technique, a spatial adaptive hierarchical tree is first generated. In our algorithm, the tree structure is identical to that in FMM or FDS and is generated by a recursive partition to divide the particle-occupant region into nested cubical boxes, where the root box is the smallest bounding box that contains the entire particle set. Without loss of generality, the root box is normalized to size $1$ along each side. The root box is partitioned equally along each dimension. The partition continues recursively on the resulting box until the box contains no more than $s$ beads, at which point it becomes a leaf node. Empty boxes encountered during partition are pruned off. In our implementation, the value $s$ is chosen depending on the size of the particle set and other factors to allow an optimal performance.

3.2 Low rank compression

The low-rank structure for well-separated source and target points has been extensively studied for the free-space FMM and FDS algorithms. Consider $N$ sources with strength $q_{j}$ placed at $\mathbf{x}_{j}=(x_{j},y_{j})$ in a circle centered at $\mathbf{x}_{c}=(x_{c},y_{c})$ with radius $R$ on top of the half-space, and suppose we are interested in the field at $\mathbf{x}$ due to all the source points given by $u(\mathbf{x})=\sum_{j=1}^{N}q_{j}g(\mathbf{x},\mathbf{x}_{j})$ where $g(\mathbf{x},\mathbf{x}_{j})$ is the free-space Green’s function contribution. We say $\mathbf{x}$ is well-separated from the sources if the distance between $\mathbf{x}$ and the source circle center $\mathbf{x}_{c}$ is at least $3R$, see Fig. 1.

Free-space Green’s function compression. Using Graf’s addition theorem [2], the free-space Green’s function interaction of well-separated sources $\mathbf{x}_{j}$ and target $\mathbf{x}$ can be compressed as a “multipole expansion” given by

[TABLE]

where

[TABLE]

$\theta_{c}$ is the polar angle of $\mathbf{x}-\mathbf{x}_{c}$, $(\rho_{j},\theta_{j})$ are the polar coordinates of the complex number $\mathbf{x}_{j}-\mathbf{x}_{c}$, and the number of terms $P$ is a constant independent of the number of the sources $N$.

Domain Green’s function compression: complex image representation. For the half-space problem with an impedance boundary condition, the field at $\mathbf{x}$ due to all the source points is

[TABLE]

As the first term on the right hand side represents the free-space Green’s function interaction and is already compressed in Eq. (20), we focus on the compression of the second term representing the scattered field. When using the complex image representation, it can be compressed simply as follows:

[TABLE]

where $\mathbf{x}_{j}^{im}=(x_{j},-y_{j})$ are the coordinates of the point-image charge, $\bar{\alpha}_{p}$ is the complex conjugate of the free-space multipole coefficient $\alpha_{p}$ in Eq. (21) (see [34]), and

[TABLE]

Therefore, the “multipole expansion” for both the original and image sources is

[TABLE]

which is the key formula behind the proposed new heterogeneous FMM for layered media.

We emphasize that the number of terms $P$ for the scattered field expansion is the same as the one in the free-space expansion for the same accuracy requirement. This can be rigorously justified by the observation that when the original sources $\mathbf{x}_{j}$ (in the green circle centered at $\mathbf{x}_{c}$) in Fig. 1 are well-separated from the target point $\mathbf{x}$ (in the blue circle centered at $\mathbf{x}_{c}^{l}$), all the corresponding point-images $\mathbf{x}_{j}^{im}$ (in the circle centered at $\mathbf{x}_{c}^{im}$) are also well-separated from $\mathbf{x}$, and this is true also for the set of line-images on the rays starting from $\mathbf{x}_{c}^{im}$ with the same $s$ value.

Remark 1: Eq. (26) suggests that for the domain Green’s function interactions, when the source and target clusters are well-separated, it is possible to only compress the translation invariant free-space Green’s function using a $P$-term multipole expansion with coefficients $\alpha_{p}$ as in Eq. (20) for a prescribed accuracy requirement, and all other related information in Eq. (26) can be recovered from $\alpha_{p}$ to the same accuracy. Also, unlike the FDS, the compression of the domain Green’s function is not performed directly on the matrix entries, but on another matrix after some spatially variant transformations implicitly described in Eq. (26), and these transformations involve the complex conjugation operator. Finally, as the compression is only on the free-space Green’s function, deriving parent’s compressed representation and corresponding error analysis are exactly the same as those in the classical FMM algorithms, where the “multipole-to-multipole” translation can be used without any modification. We therefore skip the details of this translation operator in this paper.

Domain Green’s function compression: Sommerfeld integral representation. For multi-layered media, deriving the image representation of the domain Green’s function becomes very complicated, and most existing techniques apply the integral transformations to obtain the Sommerfeld representation directly. Instead of deriving the compression formula using the Sommerfeld integral domain Green’s function in Eq. (12), we reformulate its compressed image representation in Eq. (26) as an equivalent Sommerfeld integral representation for the sources of $N$ particles. Then, it is compared directly with the uncompressed Sommerfeld representations. This comparison provides insight on a direct compression of the domain Green’s function in the spectral domain.

We start from the Sommerfeld representation of $H_{n}({{k}}\rho)e^{in\theta}$. Applying the relation

[TABLE]

where $(\rho,\theta)$ are the polar coordinates of the complex number $x+iy$, and using the Sommerfeld representation of $H_{0}({k}\rho)$ given in Eq. (3), we have for $y>0$,

[TABLE]

Plugging this representation into Eq. (26), and integrating the $s$ variable analytically, we have the compressed Sommerfeld representation directly as

[TABLE]

Comparing with the uncompressed Sommerfeld representation by adding up the Sommerfeld representation of the scattered field in Eq. (12) for each source in $u^{s}(\mathbf{x})=\sum_{j=1}^{N}q_{j}u_{\mathbf{x}_{j}}^{s}(\mathbf{x})$, we have

[TABLE]

We further notice that $\bar{\alpha}_{p}$ is independent of $\lambda$, and the compressed term in Eq. (28) is the Laurent expansion in $z=\frac{\lambda-\sqrt{\lambda^{2}-{k}^{2}}}{{k}}$ of the uncompressed term in Eq. (29).

Remark 2: Comparing Eq. (28) with Eq. (29), we can identify the roles of different terms in the Sommerfeld representation of the domain Green’s function. In particular, we see that the conjugates of the free-space multipole expansion coefficients are the same as the Laurent expansion ones in the compressed representation. This observation reveals how the domain Green’s function can be compressed directly when the line image $\mu(s)$ in Eq. (8) is complicated or unavailable. For instance, for the top layer Sommerfeld domain Green’s function in Eq. (15) for the three layered media setting, as the terms “free-space info” and “uncompressed” have the same structure as in Eq. (29), and the term “image info” is independent of $\mathbf{x}$ and $\mathbf{x}_{0}$, we can therefore simply compute the free-space multipole expansion either directly from the sources, or through the free-space “multipole-to-multipole” translations, and the results will directly give a compressed Sommerfeld representation similar to Eq. (28). In Sec. 5, we apply this observation and present an algorithm framework purely based on the Sommerfeld representation.

Remark 3: The error analysis of the direct compression of the Sommerfeld integral representation described in Eqs. (28) and (29) is not an easy task. Luckily for the 2-D half-space problem, the error analysis becomes trivial when performed in the physical domain using the image representation. Note that the error analysis only requires that the target box and all the image sources are well-separated, which can be easily carried out for the multi-layered case in Eq. (15) using repeated image reflections, without knowing the exact density and location of the image.

Local expansions for received information. We have so far discussed how the domain Green’s function can be compressed in different ways so the compressed representations can be transmitted to the “receiving” nodes on the hierarchical tree structure. In the hierarchical algorithms, the received information is also stored in some compressed compact form. For the Helmholtz equation, the most convenient and commonly used form is the Bessel functions based “local expansion” as in the classical FMM algorithms. We adopt this representation in this paper to store the received information, so the “local-to-local” translations in the free-space FMM algorithms can be used without any modification. Other compact forms are also possible, for example, the “equivalent source” representation in the kernel independent FMM algorithms [35, 36].

3.3 Translations on the hierarchical tree structure

We discuss how the compressed representations can be transmitted on the hierarchical tree structure in this section. As our selected “multipole” and “local” representations of the compressed domain Green’s function are the same as those for the translation invariant free-space Green’s function, existing “multipole-to-multipole” and “local-to-local” translations in the free-space FMM can be applied without any modification. We therefore focus on the “multipole-to-local” (M2L) translation operator, and study how the multipole expansion of the compressed domain Green’s function can be converted to local expansions.

We start from the following well-known M2L translation operator for the free-space Green’s function. Consider the same source points $\mathbf{x}_{j}$, $j=1,\cdots,N$ described in Fig. 1 and the compressed representation of the free-space kernel in Eq. (20). Then the potential $u(\mathbf{x})$ can be translated to a local expansion using Graf’s addition theorem as

[TABLE]

where the coefficients are

[TABLE]

$\theta$ is the polar angle of $\mathbf{x}-\mathbf{x}_{c}^{l}$, and $(\rho_{l},\theta_{l})$ are the polar coordinates of $\mathbf{x}_{c}^{l}-\mathbf{x}_{c}$. Because the complex image representation of the domain Green’s function is given in terms of the free-space Green’s function, we can therefore plug the free-space M2L translation formula in the compressed image representation of the scattered field

[TABLE]

to derive its local expansion given by

[TABLE]

where the local expansion coefficients are given by

[TABLE]

$\theta$ is the polar angle of $\mathbf{x}-\mathbf{x}_{c}^{l}$, and

[TABLE]

The local expansion for $u(\mathbf{x})$ is simply the sum of the free-space Green’s function and scattered field local expansions. As the translation operator from the compressed “multipole coefficients” $\{\alpha_{p}\}$ to the local coefficients $\{\beta_{p}^{s}\}$ involves the complex conjugate operator, for notation reasons, instead of combining the free-space with the complex image contributions in one single translation, we only construct the mapping matrix $A$ for the scattered field,

[TABLE]

where

[TABLE]

Notice that the integrand in Eq. (38) is highly oscillatory for large $s$ and its numerical computation usually requires special treatment, for instance, by choosing different integration contours. In the following, we use the Sommerfeld representations of $H_{n}({k}\rho)e^{in\theta}$ in Eq. (27), but separate it to the propagating and evanescent terms as in Eq. (4), and reformulate $H_{m-p}({{k}}\tilde{\rho}_{im})e^{i(m-p)\tilde{\theta}_{im}}$ using the plane wave representation

[TABLE]

where $(x,y)$ are the Cartesian coordinates of $(\tilde{\rho}_{im},\tilde{\theta}_{im})$ and

[TABLE]

for $y>0$. This plane wave representation was also used to diagonalize the M2L translation operator in the new version of the low frequency FMM for the free-space Green’s function in Ref. [21]. We skip the similar formula for $H_{m-p}({k}\hat{\tilde{\rho}}_{im})e^{i(m-p)\hat{\tilde{\theta}}_{im}}$, and present the translation matrix explicitly as

[TABLE]

after integrating the $s$ variable analytically. In the numerical evaluation, as the integral of the propagating term is over a finite interval, high order Gauss quadrature can be applied and for the evanescent term, the generalized Laguerre quadrature with weight function $t^{n}e^{-t}$ is used.

The translation matrix $A$ in Eq. (40) has several special features. Unlike in the classical FMM algorithms, it depends on the $x$ and $y$ and is therefore spatially variant. However, for the fixed $x$ and $y$, it is not a two variable function of $m$ and $p$ and only depends on $m-p$. In the numerical implementation, the matrix can be either computed on-the-fly using high order Gauss and Laguerre quadratures, or precomputed and stored. We can estimate the required storage in the algorithm as follows: the translation operators $A_{p,m}$ are needed for all levels of the tree. For a fixed box, translation matrix consist of $4p$ complex values (as it is only a function of $m-p$) and there are a total of no more than $7\cdot 7=49$ surrounding boxes representing the well-separated “receiving” boxes. For the two-layered media case, the matrix also depends on the $y$-coordinate of the center of the box as the translation operator takes different values as their images change. Thus, we can conclude that at tree level $l$, there are $2^{l}$ different values of $y$-coordinates, and for each $y$-coordinate $49$ possible well-separated boxes that requires $4p$ complex values. Therefore, the total required storage for a system with $L$-levels is approximately $(2^{L+1}\cdot 49\cdot 4p)\cdot 16$ bytes, which is very small compared with the required storage for different expansions.

3.4 Accelerated evaluation of local direct interactions

We consider the submatrix representing the “local direct” interactions in this section. For a source box with $N_{b}$ particles located at $\{\mathbf{x}_{j}\}_{j=1}^{N_{b}}$, its domain Green’s function contribution to a target point $\mathbf{x}$ in a neighboring box is defined as (using the complex image representation)

[TABLE]

This formula shows the entries in one row of the submatrix. Further compression of the submatrix is usually impossible as it is not low-rank. However, it is still possible to take advantage of the compressed scattered field representations of the domain Green’s function given in Eq. (23) (complex image representation) or Eq. (28) (Sommerfeld integral representation), so the entries in the submatrices can be evaluated more efficiently.

In Fig. 2, we show a source box sitting next to a target box. In this figure, notice that most of the line-images are well-separated from the target box. We can therefore choose an appropriate constant $C$ and cut the line-images into two parts: those that are well-separated from the target box and those that are not. The evaluation of Eq. (41) can be divided as $u^{d}(\mathbf{x})=I+II$, where

[TABLE]

The first summation $I$ is computed directly using high order quadrature for the finite size integral. For the second summation $II$, because $\mathbf{x}_{j}^{im}-s\hat{y}$ is far away from the target point, the computation can be accelerated using the available source box multipole expansion as follows.

[TABLE]

where

[TABLE]

and the translation matrix is given by

[TABLE]

which can be evaluated efficiently using the corresponding Sommerfeld integral representation.

In the hierarchical tree structure, most boxes are well-separated from the interface $y=0$. This implies that $C=0$ for most direct interactions of the source and target boxes, and the separation can be simplified as

[TABLE]

where both the point- and line-image contributions belong to $II$. In this case, the corresponding translation operator becomes the same as in Eq. (40). In the numerical simulation, all the translation matrices can be either precomputed or computed on-the-fly using high order quadrature for the Sommerfeld integral representation.

The selected compression schemes and translations allow easy adaptation of existing fast multipole algorithms for computing the domain Green’s function interactions of the 2-D half-space Helmholtz equation with impedance boundary condition. Our solver is based on the wide-band FMM [10, 12]. We present the pseudo-code of our algorithm in Algorithm 1.

In our current implementation, all the tables are precomputed using Mathematica requesting more than 20 digits accuracy. Compared with the original free-space FMM algorithm, the domain Green’s function FMM only requires a small portion of additional cost, as demonstrated in the next section.

4 Numerical Results

We present some preliminary numerical results in this section to demonstrate the performance of the new heterogeneous FMM algorithm for the two-layered media with the interface placed at $y=0$, and set $\alpha=1$ in the impedance boundary condition. We assume the source and target points are the same set of $N$ particles located in a unit box centered at $(0,1.5)$ as shown in Fig. 3. The numerical simulations are performed on a desktop with 3.7 GHz Xeon E5 processor and 32GB RAM using the gcc compiler version 4.9.3. All the required translation tables are precomputed using Mathematica.

As the analytical solution is not available for this problem, we first check the algorithm accuracy by studying how the errors change as a function of the number of expansion terms $p$. We consider the example with $N=100\times 100$ particles uniformly distributed in the box. A reference solution is computed using $p=39$ (which should provide results with approximately $12$-digit accuracy) and by setting $L=3$ in the hierarchical tree structure. The new heterogeneous FMM algorithm took about $1.19$ seconds to derive the reference solutions. The accuracy results are presented in Table 1 for $p=5,10,20,30$ and wave numbers ${k}=0.1$ and ${k}=1$, respectively, where the error $E_{p}$ using $p$ terms in the expansion is defined as

[TABLE]

In Fig. 5(a), we plot how the error $E_{p}$ decays as a function of $p$ for ${k}=0.1$. We see that the error dependency on the number of terms $p$ to compress the domain Green’s function is similar to that in existing FMM analysis for the free-space kernels.

We demonstrate the algorithm efficiency by presenting the CPU times in Table 2 for different numbers of source/target points $N$ from $100$ to $1,000,000$ for ${k}=0.1$. A $\log$-$\log$ plot of the CPU time is also presented in Fig. 5(b), which clearly shows the linear scaling of the new heterogeneous FMM algorithm. For comparisons, estimated results of direct computations (using CPU times for $N=100$ and $6,400$) and the ideal linear scaling curve are also presented. Similar experiments are performed for ${k}=1$ and results are almost identical to that when ${k}=0.1$ and are therefore omitted in this paper.

5 Heterogeneous FMM for Three Layers

In this section, we study a three-layered media setting where the source and target points are all located in the top layer and their interactions are described by Eq. (15), and present the framework of a hierarchical algorithm based on the Sommerfeld integral representation. This framework allows direct generalization to multiple layers with analysis and algorithm design for more general settings. All the technical details of the compressions, translation operators, and bookkeeping strategies for different numbers of layers, boundary conditions, and locations of the source and target points will be presented in a subsequent paper.

Compressed Representation: We start from the Sommerfeld integral representation of the domain Green’s function for sources and target points located in the top layer described in Eq. (15). Similar to Eqs. (28) and (29) in Sec. 3.2 for the two-layered case, we first derive the compressed “multipole” expansion in the form

[TABLE]

where $(x_{c},y_{c})$ represents the center of the source box, and the coefficients $\bar{\alpha}_{p}$ are either computed directly using the free-space “particle-to-multipole” translation described in Eq. (21) for a childless box, or computed using the free-space “multipole-to-multipole” translation operator for a parent box.

Remark 4: We expect that the error analysis for the three layers, in fact for multiple layers for this matter, should follow in a manner similar to the case of two layers because the point and line images in the domain Green’s function used to compute the far field will stay far away from the far field box as illustrated in Fig. 1. We will present the construction of the images for multiple layers and the mathematical error analysis of compression and translation operators in the follow-up paper.

Moreover, it should be noted that the error analysis can also be performed without using the image representation, by observing that the decay rate of the series expansion is determined by the ratio of $||(x,y)-(x_{c},-y_{c})||$ and the size of the source box. As the source and target boxes are well-separated and the distance $||(x,y)-(x_{c},-y_{c})||>||(x,y)-(x_{c},y_{c})||$, we conclude that the number of terms required in the expansion is no more than that in the free-space.

Multipole-to-local Translation: Similar to the two-layered media, we use the Hankel function based free-space “local” expansions to collect the far-field compressed representation. There are several ways to derive the transformation matrix from the multipole coefficients $\bar{\alpha}_{p}$ in Eq. (48) to the local coefficients. For instance, one can use the identity (see [2])

[TABLE]

to derive the local expansion coefficients after shifting the multipole expansion in Eq. (48) to the new center of the target box. This approach was also used in the new version of FMM algorithm for the Helmholtz equation [21]. A much simpler approach is to compare with the two-layered media formula and modify Eq. (40) to directly derive the translation matrix as

[TABLE]

The translation matrix can be either precomputed, or computed on-the-fly using high order Gauss and Laguerre quadratures. Once the local expansion is available, the “local-to-local” translations are the same as the free-space ones.

Accelerating the evaluation of local direct interactions. Similar to the two-layered case, the local direct interactions can be separated into the free-space interactions and scattered field due to the images. However, finding the corresponding explicit image representation from the Sommerfeld integral usually requires a numerical integral transformation and how to reduce its computational complexity is still a research topic. We therefore use the following strategies to accelerate the evaluation of local direct interactions: (a) when the source box is located at least one box size away from the interface $y=0$, the compressed multipole expansion in Eq. (48) for the scattered field is then valid at the target box, we therefore perform a “multipole-to-local” translation using a precomputed table, which will add the scattered field contribution to the local expansion of the box; (b) otherwise, we compute the direct interactions using numerical integration schemes for the domain Green’s function.

Remark 5: Mathematically, the cut-off idea is the same as separating the term $\sigma_{1}$ in Eq. (15) as the sum of two terms, so that one term contains a factor $e^{-\sqrt{\lambda^{2}-k_{1}^{2}}C}$ which makes the truncated “multipole expansion” in Eq. (48) valid for a prescribed accuracy requirement, while the other term when transformed into the physical domain, only involved a finite integral where the “images” are located in a small region in the neighborhood of the target box. The numerical separation of $\sigma_{1}$ may be more difficult than computing the Sommerfeld representation directly.

Algorithm Summary. The hierarchical algorithm for the special three-layered media setting can be summarized as follows. It uses the same tree structure, free-space compressed multipole coefficient, multipole-to-multipole, and local-to-local translation operators as the two-layered algorithm. Also, when the source box is well separated from the interface $y=0$, its local direct interactions can be accelerated similarly as in the two-layered media case. There are only two necessary revisions of the two-layered media algorithm: (a) a different set of tables for the “multipole-to-local” translation operators are used; and (b) when the “images” in the scattered field are no longer well-separated from the target boxes when evaluate the local direct interactions, instead of using a cut-off distance, the domain Green’s function is evaluated directly using some numerical integration schemes. As the algorithm structure is almost identical to the two-layered case, we neglect the details in this paper.

Remark 6: In the numerical experiment, when the configuration of the source and target points in the top layer is the same as that in the upper layer of the half-space problem in Fig. 3, the three-layered media algorithm only differs from two-layered one in the precomputed tables, and both the accuracy and efficiency are therefore the same and we neglect the details.

6 Conclusion and Future Work

We present a heterogeneous FMM for the efficient calculation of the discretized integral operator for the Helmholtz equation in two and three layers media with details on the case of two-layered media with impedance boundary conditions. The two-layered media setting allows the use of both the complex line image and Sommerfeld integral representations to compress the domain Green’s function and to derive the translation operators analytically. Instead of compressing the interaction matrix directly, the complex line image representation intuitively reveals how a “transformed” matrix can be compressed through a procedure that only involves the free-space Green’s function, and provides rigorous error bounds by using existing free-space FMM results.

Unlike the fast direct solvers, the compression is performed analytically on a “transformed” matrix which allows the easy adaptation of existing free-space FMM packages. Also disimilar to the classical FMM, the “multipole-to-local” translation operators are spatially variant, thus the translation operators in the FMM becomes heterogeneous. Numerical experiments show that the new hierarchical algorithm provides significant improvement over existing hybrid methods [34] in two-layered media settings. By connecting the complex line image and Sommerfeld representations in the two-layered media algorithm, we demonstrate how the compressions and translations can be performed directly on the Sommerfeld representation for a particular setting in the three-layered media for the case when all the source and targets points are located in the top layer.

This paper focuses on the intuitions through the two-layered and three-layered media setting. In a subsequent paper, we will present the H-FMM for multiple layered media including constructions and error analysis of compressions, translations as well as bookkeeping strategies, for both the scalar Helmholtz equations in acoustic studies and the multi-layered media dyadic Green’s function for the Maxwell’s equations [13].

Finally, it is interesting to compare the analysis based compressions with those using purely numerical linear algebra techniques as in the fast direct solvers, to understand how the efficiencies of both compressions can be further improved. Research along these directions will be explored in the future.