Laplace Green’s Functions for Infinite Ground Planes with Local Roughness

Nail A. Gumerov and Ramani Duraiswami
University of Maryland Institute for Advanced Computer Studies, College Park, MD 20742

Abstract

The Green’s functions for the Laplace equation respectively satisfying the Dirichlet and Neumann boundary conditions on the upper side of an infinite plane with a circular hole are introduced and constructed. These functions enables solution of the boundary value problems in domains where the hole is closed by any surface. This approach enables accounting for arbitrary positive and negative ground elevations inside the domain of interest, which, generally, is not possible to achieve using the regular method of images. Such problems appear in electrostatics, however, the methods developed apply to other domains where the Laplace or Poisson equations govern. Integral and series representations of the Green’s functions are provided. An efficient computational technique based on the boundary element method with fast multipole acceleration is developed. A numerical study of some benchmark problems is presented.

keywords:

Green’s Function, Laplace’s Equation , Electrostatics , infinite ground

MSC:

[2010] 00-01, 99-00

1 Introduction

Many problems of physical interest involve fields satisfying Laplace’s equation with appropriate boundary conditions. Often, the boundaries of the domain might include surfaces of infinite extent, which for analytical convenience, are often considered to be infinite planes, so that they can be treated via the method of images. Examples include electrostatics, magnetostatics, conduction (or diffusion) from a surface at a fixed temperature (or concentration), potential flow of an ideal fluid over a body. Motivating the work reported here are the problems of electric and magnetic fields near an infinite ground, where the surface elevation departs from the plane in both the positive and negative directions, e.g., a ground with bumps and troughs.

Computation of electrostatic fields in scales where the effects of ground should be taken into account can be an important task for a number of applications, for example, for modeling fields in urban environments [1], fields generated by the transmission power lines [2], [3], and fields generated by lightning over rough ocean surface [4]. Scattering of higher frequency electromagnetic waves from rough surfaces is also of interest [5]. Note that if the wavelength is much larger than the scale of the roughness, electrostatic approximation for characterization of scattering properties of the surface is also valid.

Most of numerical methods used for such computations either neglect the presence of objects outside the computational domain or impose some boundary conditions on the domain surface (e.g., see monograph [6] and review [7]). The boundary elements can deliver solutions which provide decay of the electric potential at the infinity, while solving a boundary integral equation for the charge distribution only for the objects located inside the computational domain. It is relatively easy to account the effects of the infinite flat ground in such computations using the method of images [8]. Formally, this can be done by replacing the free space Green’s function in the boundary integral equation by the Green’s function for the semispace. This function takes zero at any point on the ground plane and so any convolution with such a Green’s function produces a solution satisfying the boundary condition on the ground. Some techniques for construction of Green’s functions for the Laplace equation in different 2D domains are developed (e.g., [9]), but they formulated in terms of complex variables and conformal mappings, which makes difficult to extend them for 3D.

In practice, the situation frequently appears to be more complicated as the ground has elevations and dips. The term “ground” can also include buildings or other construction, and can be modeled using surface boundary element meshes. In any case, such modeling is limited to the area where the solution is needed and detailed representation of the environment can be provided only inside the computational domain. The question then is how to take into account the effect of the environment outside the computational domain? One possibility is to assume periodic boundary conditions. In this case some “periodization” techniques for solution of the Laplace equation can be used (e.g., see [10]; quasi-periodic Green’s functions can be found in [11]). Another approach can be to model the environment outside some ball $\mathbf{B}$ of radius $R$ as a flat ground, while provide a detailed meshing of all objects and the ground inside $\mathbf{B}$ (sometimes called “locally rough surface”). In the present paper we take the latter approach.

A simple solution of the problem can be found if all points on the ground inside $\mathbf{B}$ are located above or on the level of the flat ground at infinity (which we can denote as zero level). In this case the method of images can be applied to any objects in the domain. However, if there exist some dips, or points below the zero level the method of images is not applicable, since it requires placing the image source on the side of the ground plane opposite to the side where the actual source is located, which creates a physically unacceptable (singular) solution. In a number of papers treatment of such situation using boundary integrals for the Helmholtz equation was considered (e.g., [12], [13]). In this paper, we take a different approach and derive an analogue of the Green’s function, which accounts for the effect of infinite flat ground plane outside the computational domain, and which has no singularities besides the source location inside $\mathbf{B}$ . This function is suitable for solution of problems with arbitrary location of boundaries inside $\mathbf{B}$ . We implemented the method and provide some numerical examples. Note that this method can also be applied to the Helmholtz equation.

The problem described above formally reduces to solution of the Laplace equation with the Dirichlet boundary conditions on the ground. In practice, there appear also situations when the zero Neumann boundary conditions should be imposed on an infinite surface. For example, if we wish to determine underwater electric fields generated by some sources, this can be formulated as the problem with zero Neumann boundary conditions on the sea level (since the electric conductivity of the sea water is much larger than that of the air). The ocean surface can be rough, so that can be taken into account for points close to the surface, while the ocean surface can be considered as flat in the far field (also, see [14], [15], [16]). Another example of the Neumann problem appears in aerodynamics, where the effects of ground should be taken into account to determine correctly the lift force for objects flying close to the ground [17]. For flat surfaces the Neumann problem can be also solved using the method of images, which causes the same problems as for the Dirichlet problem described above. The Green’s functions for the Dirichlet and Neumann boundary conditions are closely related. In this paper we address both problems.

2 Statement of the problem

Refer to caption — Figure 1: The problem geometry and notation. A solution needs to be obtained in domain $\Omega_{0}$ bounded by the object surface $S$ , ground surface $S_{g0}$ and a hemisphere of radius $R_{0}$ . The ground surface is flat outside $\Omega_{0}$ . An extended computational domain $\Omega_{e}$ , which is a part of $\Omega$ located inside a sphere of radius $R_{e}$ can be introduced for efficient computations.

First, we consider the Dirichlet problem. Consider an infinite ground surface $S_{g}$ on which the electric potential $\phi$ is zero and let there be some objects on it. Let the surface be $S$ the value $V$ of the potential is given (here and below we use dimensionless equations). In the electrostatics approximation, we have the following Dirichlet problem for the Laplace equation,

	$\displaystyle\nabla^{2}\phi\left(\mathbf{y}\right)$	$\displaystyle=$	$\displaystyle 0,\quad\mathbf{y}\in\Omega,$		(1)
	$\displaystyle\left.\phi\right\|_{\mathbf{y\in}S_{g}}$	$\displaystyle=$	$\displaystyle 0,\quad\left.\phi\right\|_{\mathbf{y\in}S}=V,\quad\left.\phi\right\|_{\left\|\mathbf{y}\right\|\mathbf{=}\infty}=0,$

where $\Omega$ is an infinite domain in $\mathbb{R}^{3}$ bounded by $S$ and $S_{g}$ . The solution of this problem can be obtained using indirect boundary integral method, which represents the potential in the form of the single layer potential

\phi\left(\mathbf{y}\right)=\int_{S\cup S_{g}}\sigma\left(\mathbf{x}\right)G\left(\mathbf{y},\mathbf{x}\right)dS\left(\mathbf{x}\right),\quad G\left(\mathbf{y},\mathbf{x}\right)=\frac{1}{4\pi\left|\mathbf{y}-\mathbf{x}\right|},

(2)

where $G\left(\mathbf{y},\mathbf{x}\right)$ is the free-space Green’s function for the Laplace equation, and $\sigma\left(\mathbf{x}\right)$ is the surface charge density, which can be obtained from the solution of the boundary integral equation,

\int_{S\cup S_{g}}\sigma\left(\mathbf{x}\right)G\left(\mathbf{y},\mathbf{x}\right)dS\left(\mathbf{x}\right)=\left\{\begin{array}[]{c}V\left(\mathbf{y}\right),\quad\mathbf{y}\in S,\\ 0,\quad\mathbf{y}\in S_{g}.\end{array}\right.

(3)

Assume further that $S_{g}$ can be modelled by a flat horizontal surface (zero elevation, $z=0$ ), except for a finite part near the source, and in the region in which we wish to compute the solution in details. The ground plane in this region may have arbitrary positive and negative elevations. We assume that this region is contained inside a ball $\mathbf{B}_{0}$ of radius $R_{0}$ centered at the origin of the reference frame. It is also assumed that surface $S$ is located inside the ball. The ball partitions $\Omega$ into domains $\Omega_{0}$ and $\Omega_{\infty}$ and $S_{g}$ into $S_{g0}$ and $S_{g\infty}$ (see Fig. 1),

	$\displaystyle\Omega_{0}$	$\displaystyle=$	$\displaystyle\left\{\mathbf{y}\in\Omega,\quad\left\|\mathbf{y}\right\|<R_{0}\right\},\quad\Omega_{\infty}=\left\{\mathbf{y}\in\Omega,\quad\left\|\mathbf{y}\right\|\geqslant R_{0}\right\}=\left\{\mathbf{y},\quad\left\|\mathbf{y}\right\|\geqslant R_{0},\text{ }z>0\right\},$		(4)
	$\displaystyle S_{g0}$	$\displaystyle=$	$\displaystyle\left\{\mathbf{y}\in S_{g},\quad\left\|\mathbf{y}\right\|<R_{0}\right\},\quad S_{g\infty}=\left\{\mathbf{y}\in S_{g},\quad\left\|\mathbf{y}\right\|\geqslant R_{0}\right\}=\left\{\mathbf{y},\quad\left\|\mathbf{y}\right\|\geqslant R_{0},\text{ }z=0\right\}.\mathbf{\quad}$

Note now that $S_{g\infty}$ is an open surface, for which we have an “upper” side, $S_{g\infty}^{+}$ , faced towards $z>0$ , and a “lower” side, $S_{g\infty}^{-}$ faced towards $z<0$ . The boundary conditions need be satisfied on the upper side only. The Green’s function $G^{(D)}\left(\mathbf{y},\mathbf{x}\right)$ for the zero Dirichlet boundary conditions on $S_{g\infty}^{+}$ can be introduced as the solution of the problem

	$\displaystyle\nabla_{\mathbf{y}}^{2}G^{(D)}\left(\mathbf{y},\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle-\delta\left(\mathbf{y-x}\right),\quad\quad\mathbf{y,x}\in\mathbb{R}^{3}\mathbf{,}$		(5)
	$\displaystyle\left.G^{(D)}\left(\mathbf{y},\mathbf{x}\right)\right\|_{\mathbf{y\in}S_{g\infty}^{+},\mathbf{y}\neq\mathbf{x}}$	$\displaystyle=$	$\displaystyle 0,\text{\quad}\left.G^{(D)}\left(\mathbf{y},\mathbf{x}\right)\right\|_{\left\|\mathbf{y}\right\|\mathbf{=}\infty,\mathbf{y}\neq\mathbf{x}}=0.$

The solution of the Dirichlet problem (1) then can be sought in the form

\phi\left(\mathbf{y}\right)=\int_{S\cup S_{g0}}\sigma^{(D)}\left(\mathbf{x}\right)G^{(D)}\left(\mathbf{y},\mathbf{x}\right)dS\left(\mathbf{x}\right),

(6)

where $\sigma^{(D)}\left(\mathbf{x}\right)$ satisfies the boundary integral equation

\int_{S\cup S_{g0}}\sigma^{(D)}\left(\mathbf{x}\right)G^{(D)}\left(\mathbf{y},\mathbf{x}\right)dS\left(\mathbf{x}\right)=\left\{\begin{array}[]{c}V\left(\mathbf{y}\right),\quad\mathbf{y}\in S,\\ 0,\quad\mathbf{y}\in S_{g0}.\end{array}\right.

(7)

We seek to derive analytical expressions for, and methods for computation of, the function $G^{(D)}\left(\mathbf{y},\mathbf{x}\right)$ .

Similarly, we can also get the Green’s function $G^{(N)}\left(\mathbf{y},\mathbf{x}\right)$ for the zero Neumann boundary conditions on $S_{g\infty}$ , which can be found by solving

	$\displaystyle\nabla_{\mathbf{y}}^{2}G^{(N)}\left(\mathbf{y},\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle-\delta\left(\mathbf{y-x}\right),\quad\quad\mathbf{y,x}\in\mathbb{R}^{3}\mathbf{,}\quad$		(8)
	$\displaystyle\left.\frac{\partial}{\partial n\left(\mathbf{y}\right)}G^{(N)}\left(\mathbf{y},\mathbf{x}\right)\right\|_{\mathbf{y\in}S_{g\infty}^{+},\mathbf{y}\neq\mathbf{x}}$	$\displaystyle=$	$\displaystyle 0,\text{\quad}\left.G^{(N)}\left(\mathbf{y},\mathbf{x}\right)\right\|_{\left\|\mathbf{y}\right\|\mathbf{=}\infty,\mathbf{y}\neq\mathbf{x}}=0.$

Construction of this function allows one to solve the problem with mixed boundary conditions

	$\displaystyle\nabla^{2}\phi\left(\mathbf{y}\right)$	$\displaystyle=$	$\displaystyle 0,\quad\mathbf{y}\in\Omega,$		(9)
	$\displaystyle\left.\frac{\partial}{\partial n\left(\mathbf{y}\right)}\phi\right\|_{\mathbf{y\in}S_{g}}$	$\displaystyle=$	$\displaystyle 0,\quad\left.\phi\right\|_{\mathbf{y\in}S}=V,\quad\left.\phi\right\|_{\left\|\mathbf{y}\right\|\mathbf{=}\infty}=0,$

in form (6), where superscript $\left(D\right)$ should be changed to $\left(N\right)$ and the charge density can be determined from the boundary integral equation

	$\displaystyle\int_{S\cup S_{g0}}\sigma^{(N)}\left(\mathbf{x}\right)G^{(N)}\left(\mathbf{y},\mathbf{x}\right)dS\left(\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle V\left(\mathbf{y}\right),\quad\mathbf{y}\in S,$		(10)
	$\displaystyle\frac{\partial}{\partial n\left(\mathbf{y}\right)}\int_{S\cup S_{g0}}\sigma^{(N)}\left(\mathbf{x}\right)G^{(N)}\left(\mathbf{y},\mathbf{x}\right)dS\left(\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle 0,\quad\mathbf{y}\in S_{g0}.$

Note that functions $G^{(D)}\left(\mathbf{y},\mathbf{x}\right)$ and $G^{(N)}\left(\mathbf{y},\mathbf{x}\right)$ defined above are not unique, since the boundary condition is specified only on $S_{g\infty}^{+}$ , while an arbitrary boundary condition can be imposed on $S_{g\infty}^{-}$ . Depending on the type of condition on $S_{g\infty}^{-}$ the functions may or may not satisfy the symmetry typical for Green’s functions for Dirichlet problems for closed or two-sided open surfaces with homogeneous boundary conditions [18]. Despite this non-uniqueness any solution of problems (5) and (8) delivers the Green’s function suitable for representation of solutions in form (6 ), where the charge density is determined from Eqs (7) and (10). Indeed, these relations guarantee that the solution satisfies the Laplace equation and the boundary conditions, and, so, as $G^{(D)}$ or $G^{(N)}$ are selected the solution is unique.

We also remark that, despite $G^{(D)}\left(\mathbf{y},\mathbf{x}\right)$ and $G^{(N)}\left(\mathbf{y},\mathbf{x}\right)$ can be computed for any location of the source ( $\mathbf{x}$ ) and evaluation ( $\mathbf{y}$ ) points, for boundary integral representation of the solution we only consider the case when the source $\mathbf{x}$ is located inside the ball $\mathbf{B}$ . Moreover, $\mathbf{y}$ should be also located in $\mathbf{B}$ to solve the boundary integral equations (7) and (10). As soon as the charge density is determined from these equations the values of Green’s functions at arbitrary $\mathbf{y}$ can be used.

3 Solution

3.1 Single and double layer potentials

Solutions of the Dirichlet and Neumann problems (5) and (8) can be found using single and double layer potentials,

	$\displaystyle G^{(D)}\left(\mathbf{y},\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle G\left(\mathbf{y},\mathbf{x}\right)+K^{(D)}\left(\mathbf{y},\mathbf{x}\right),\quad K^{(D)}\left(\mathbf{y,x}\right)=\int_{S_{g\infty}}\mu_{D}\left(\mathbf{x}^{\prime}\mathbf{,x}\right)\frac{\partial G\left(\mathbf{y},\mathbf{x}^{\prime}\right)}{\partial n\left(\mathbf{x}^{\prime}\right)}dS\left(\mathbf{x}^{\prime}\right),$		(11)
	$\displaystyle G^{(N)}\left(\mathbf{y},\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle G\left(\mathbf{y},\mathbf{x}\right)+K^{(N)}\left(\mathbf{y},\mathbf{x}\right),\quad K^{(N)}\left(\mathbf{y,x}\right)=\int_{S_{g\infty}}\sigma_{N}\left(\mathbf{x}^{\prime}\mathbf{,x}\right)G\left(\mathbf{y},\mathbf{x}^{\prime}\right)dS\left(\mathbf{x}^{\prime}\right).\quad$

Here $\sigma_{N}\left(\mathbf{x}^{\prime}\mathbf{,x}\right)$ and $\mu_{D}\left(\mathbf{x}^{\prime}\mathbf{,x}\right)$ are the single and double layer densities, which should be found from the boundary conditions. These conditions show that

	$\displaystyle\left.K^{(D)}\left(\mathbf{y},\mathbf{x}\right)\right\|_{\mathbf{y\in}S_{g\infty}^{+}}$	$\displaystyle=$	$\displaystyle-\left.G\left(\mathbf{y},\mathbf{x}\right)\right\|_{\mathbf{y\in}S_{g\infty}},$		(12)
	$\displaystyle\left.\frac{\partial K^{(N)}\left(\mathbf{y},\mathbf{x}\right)}{\partial n\left(\mathbf{y}\right)}\right\|_{\mathbf{y\in}S_{g\infty}^{+}}$	$\displaystyle=$	$\displaystyle-\left.\frac{\partial G\left(\mathbf{y},\mathbf{x}\right)}{\partial n\left(\mathbf{y}\right)}\right\|_{\mathbf{y\in}S_{g\infty}}.$

The single and double layer potentials (11) obey the following symmetry property

\left.K^{(D)}\left(\mathbf{y,x}\right)\right|_{\mathbf{y\in}S_{g\infty}^{-}}=-\left.K^{(D)}\left(\mathbf{y,x}\right)\right|_{\mathbf{y\in}S_{g\infty}^{+}},\quad\left.K^{(N)}\left(\mathbf{y,x}\right)\right|_{\mathbf{y\in}S_{g\infty}^{-}}=\left.K^{(N)}\left(\mathbf{y,x}\right)\right|_{\mathbf{y\in}S_{g\infty}^{+}}.

(13)

Furthermore, we note that the single and double layer densities are related to the jumps of the potentials and their normal derivatives on the surface as

	$\displaystyle\mu_{D}\left(\mathbf{y,x}\right)$	$\displaystyle=$	$\displaystyle\left.K^{(D)}\left(\mathbf{y,x}\right)\right\|_{\mathbf{y\in}S_{g\infty}^{+}}-\left.K^{(D)}\left(\mathbf{y,x}\right)\right\|_{\mathbf{y\in}S_{g\infty}^{-}},$		(14)
	$\displaystyle\sigma_{N}\left(\mathbf{y},\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle\left.\frac{\partial K^{(N)}\left(\mathbf{y},\mathbf{x}\right)}{\partial n\left(\mathbf{y}\right)}\right\|_{\mathbf{y\in}S_{g\infty}^{-}}-\left.\frac{\partial K^{(N)}\left(\mathbf{y},\mathbf{x}\right)}{\partial n\left(\mathbf{y}\right)}\right\|_{\mathbf{y\in}S_{g\infty}^{+}}.$

Symmetry (13) then provides

	$\displaystyle\mu_{D}\left(\mathbf{y,x}\right)$	$\displaystyle=$	$\displaystyle 2\left.K^{(D)}\left(\mathbf{y},\mathbf{x}\right)\right\|_{\mathbf{y\in}S_{g\infty}^{+}},$		(15)
	$\displaystyle\sigma_{N}\left(\mathbf{y},\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle-2\left.\frac{\partial K^{(N)}\left(\mathbf{y},\mathbf{x}\right)}{\partial n\left(\mathbf{y}\right)}\right\|_{\mathbf{y\in}S_{g\infty}^{+}}.$

Boundary conditions (12) complete the solution, which can be written in the integral form

	$\displaystyle K^{(D)}\left(\mathbf{y,x}\right)$	$\displaystyle=$	$\displaystyle-2\int_{S_{g\infty}}G\left(\mathbf{x}^{\prime}\mathbf{,x}\right)\frac{\partial G\left(\mathbf{y},\mathbf{x}^{\prime}\right)}{\partial n\left(\mathbf{x}^{\prime}\right)}dS\left(\mathbf{x}^{\prime}\right),$		(16)
	$\displaystyle K^{(N)}\left(\mathbf{y,x}\right)$	$\displaystyle=$	$\displaystyle 2\int_{S_{g\infty}}G\left(\mathbf{y},\mathbf{x}^{\prime}\right)\frac{\partial G\left(\mathbf{x}^{\prime},\mathbf{x}\right)}{\partial n\left(\mathbf{x}^{\prime}\right)}dS\left(\mathbf{x}^{\prime}\right).\quad$

These equations show that $\mathbf{x}$ and $\mathbf{y}$ enter the expression in a non-symmetric way. However, due to the symmetry of the free-space Green’s function, we have

K^{(N)}\left(\mathbf{y,x}\right)=-K^{(D)}\left(\mathbf{x,y}\right),

(17)

Therefore, we consider only computation of the function $K^{(D)}\left(\mathbf{y,x}\right)$ .

3.2 Integral representations

Equations (16) already provide an integral representation of the solution. Since the domain of integration $S_{\infty}$ is simple, the integrals can be written in a more explicit form. Let us consider the Cartesian, cylindrical, and spherical coordinates of the points,

\mathbf{r}=\left(x,y,z\right)=\left(\rho\cos\varphi,\rho\sin\varphi,z\right)=r\left(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta\right),

(18)

and mark the coordinates of $\mathbf{x,y,}$ and $\mathbf{x}^{\prime}$ with subscripts $x$ , $y$ , and the prime, respectively. Note then that for points $\mathbf{x}^{\prime}$ on the surface $z^{\prime}=0$ , so

	$\displaystyle\left.\frac{\partial}{\partial z^{\prime}}\frac{1}{\left\|\mathbf{x}^{\prime}-\mathbf{y}\right\|}\right\|_{z^{\prime}=0}$	$\displaystyle=$	$\displaystyle\left.\frac{z_{y}}{\left\|\mathbf{x}^{\prime}-\mathbf{y}\right\|^{3}}\right\|_{z^{\prime}=0},$		(19)
	$\displaystyle\left.\left\|\mathbf{x}-\mathbf{x}^{\prime}\right\|\right\|_{z^{\prime}=0}$	$\displaystyle=$	$\displaystyle\left.\left(\left\|\mathbf{x}^{\prime}\right\|^{2}-2\mathbf{x}^{\prime}\cdot\mathbf{x}+\left\|\mathbf{x}\right\|^{2}\right)^{1/2}\right\|_{z^{\prime}=0}=\left(\rho^{\prime 2}-2\rho_{x}\rho^{\prime}\cos\left(\varphi^{\prime}-\varphi_{x}\right)+r_{x}^{2}\right)^{1/2}.$

Using the expression for the free-space Green’s function (2), we obtain from Eq. (16)

	$\displaystyle K^{(D)}\left(\mathbf{y},\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle K^{(D)}\left(\rho_{y},\varphi_{y},z_{y};\rho_{x},\varphi_{x},z_{x}\right)=-\frac{z_{y}}{8\pi^{2}}\int_{0}^{2\pi}\int_{R_{0}}^{\infty}fd\rho^{\prime}d\varphi^{\prime},$		(20)
	$\displaystyle f$	$\displaystyle=$	$\displaystyle\frac{\rho^{\prime}}{\left(\rho^{\prime 2}-2\rho_{y}\rho^{\prime}\cos\left(\varphi^{\prime}-\varphi_{y}\right)+r_{y}^{2}\right)^{3/2}\left(\rho^{\prime 2}-2\rho_{x}\rho^{\prime}\cos\left(\varphi^{\prime}-\varphi_{x}\right)+r_{x}^{2}\right)^{1/2}}.$

The integrand has the following asymptotic behavior for large $\rho^{\prime}$

f=\frac{1}{\rho^{\prime 3}}+\frac{1}{\rho^{\prime 4}}\left(\rho_{x}\cos\left(\varphi^{\prime}-\varphi_{x}\right)+3\rho_{y}\cos\left(\varphi^{\prime}-\varphi_{y}\right)\right)+O\left(\frac{r_{x}^{2}+r_{y}^{2}}{\rho^{\prime 5}}\right).

(21)

This shows that the integral converges. The infinite domain can be truncated and we have

K^{(D)}=K_{t}^{(D)}-\frac{z_{y}}{8\pi R_{\infty}^{2}}\left(1+O\left(\frac{r_{x}^{2}+r_{y}^{2}}{R_{\infty}^{2}}\right)\right),\quad K_{t}^{(D)}=-\frac{z_{y}}{8\pi^{2}}\int_{0}^{2\pi}\int_{R_{0}}^{R_{\infty}}fd\rho^{\prime}d\varphi^{\prime},

(22)

where $R_{\infty}$ is some cutoff radius. The second term of the asymptotic expansion in Eq. (21) does not contribute to the integral, since being integrated over the angle it produces zero. A closed form for numerical integration can be obtain by substituting the integration variable $\rho^{\prime}=R_{0}/\eta$ . We have then from Eq. (20)

	$\displaystyle K^{(D)}\left(\mathbf{y},\mathbf{x}\right)=-\frac{R_{0}^{2}z_{y}}{8\pi^{2}}\int_{0}^{2\pi}\int_{0}^{1}\frac{\eta d\eta d\varphi^{\prime}}{Q\left(\eta,\varphi^{\prime}-\varphi_{y},\rho_{y},r_{y},R_{0}\right)^{3}Q\left(\eta,\varphi^{\prime}-\varphi_{x},\rho_{x},r_{x},R_{0}\right)},$		(23)
	$\displaystyle Q\left(\eta,\varphi,\rho,r,R\right)=\left(r^{2}\eta^{2}-2R\rho\eta\cos\varphi+R^{2}\right)^{1/2}.$

3.3 Spherical harmonic expansions

While the integral representations of $K^{(D)}\left(\mathbf{y},\mathbf{x}\right)$ enable computation of this function for arbitrary locations of the source and evaluation points, numerical integration can be relatively expensive and more efficient way to compute the integral can be of practical interest. Equation (7) shows that to solve the boundary integral equations $K^{(D)}\left(\mathbf{y},\mathbf{x}\right)$ should be evaluated only for the case $\mathbf{x,}$ $\mathbf{y}\in\Omega_{0}$ . In this case the integral can be computed using spherical harmonic expansion. For a single source located at $\mathbf{x}\in\Omega_{0}$ and $\mathbf{x}^{\prime}\in\Omega_{\infty}$ the following expansion of the free space Green’s function holds [19], [20]:

\!G\left(\mathbf{x,x}^{\prime}\right)\!=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\frac{1}{2n+1}R_{n}^{-m}\left(\mathbf{x}\right)S_{n}^{m}\left(\mathbf{x}^{\prime}\right),\quad\!\left|\mathbf{x}^{\prime}\right|\!>R_{0}>\left|\mathbf{x}\right|.\!

(24)

Here $R_{n}^{m}\left(\mathbf{r}\right)$ and $S_{n}^{m}\left(\mathbf{r}\right)$ are the regular and singular (at $r=0)$ spherical basis functions

R_{n}^{m}\left(\mathbf{r}\right)=r^{n}Y_{n}^{m}\left(\theta,\varphi\right),\quad S_{n}^{m}\left(\mathbf{r}\right)=r^{-n-1}Y_{n}^{m}\left(\theta,\varphi\right),

(25)

where $Y_{n}^{m}\left(\theta,\varphi\right)$ are the orthonormal spherical harmonics,

	$\displaystyle Y_{n}^{m}\left(\theta,\varphi\right)$	$\displaystyle=$	$\displaystyle N_{n}^{m}P_{n}^{\left\|m\right\|}(\xi)e^{im\varphi},\quad\xi=\cos\theta,\quad$		(26)
	$\displaystyle N_{n}^{m}$	$\displaystyle=$	$\displaystyle(-1)^{m}\sqrt{\frac{2n+1}{4\pi}\frac{(n-\left\|m\right\|)!}{(n+\left\|m\right\|)!}},\quad n=0,1,2,...,\quad m=-n,...,n.$

where $P_{n}^{m}\left(\xi\right)$ are the associated Legendre functions defined by the Rodrigue’s formula (see [21]),

P_{n}^{m}\left(\xi\right)=\frac{\left(-1\right)^{m}\left(1-\xi^{2}\right)^{m/2}}{2^{n}n!}\frac{d^{n+m}}{d\xi^{n+m}}\left(1-\xi^{2}\right)^{n},\quad n\geqslant 0,\quad m\geqslant 0.

(27)

The derivatives of the basis functions can be expressed via the basis functions (e.g., [20])

\frac{\partial}{\partial z}S_{n}^{m}\left(\mathbf{r}\right)=-\left(2n+1\right)a_{n}^{m}S_{n+1}^{m}\left(\mathbf{r}\right),\quad

(28)

where

	$\displaystyle a_{n}^{m}$	$\displaystyle=$	$\displaystyle a_{n}^{-m}=\sqrt{\frac{(n+1+m)(n+1-m)}{\left(2n+1\right)\left(2n+3\right)}},\quad\text{for }n\geqslant\left\|m\right\|.$		(29)
	$\displaystyle a_{n}^{m}$	$\displaystyle=$	$\displaystyle 0,\quad\text{for }n<\left\|m\right\|.$

Thus, we have

\!\frac{\partial G\left(\mathbf{y,x}^{\prime}\right)}{\partial z^{\prime}}=-\sum_{n=0}^{\infty}\sum_{m=-n}^{n}a_{n}^{m}R_{n}^{-m}\left(\mathbf{y}\right)S_{n+1}^{m}\left(\mathbf{x}^{\prime}\right),\quad\!\left|\mathbf{x}^{\prime}\right|\!>R_{0}>\left|\mathbf{y}\right|,\quad z^{\prime}=z(\mathbf{x}^{\prime}).

(30)

For the computation of the integrals in (16) only the values for $z^{\prime}=0$ , or $\theta^{\prime}=\pi/2$ are needed. Using the value for $P_{n}^{|m|}(0)$ [21], we have from (26)

$\displaystyle Y_{n}^{m}\left(\frac{\pi}{2},\varphi\right)$	$\displaystyle=$	$\displaystyle L_{n}^{m}e^{im\varphi},$	(31)
$\displaystyle L_{n}^{m}$	$\displaystyle=$	$\displaystyle N_{n}^{m}P_{n}^{\|m\|}(0)=l_{n}^{m}\frac{\left(-1\right)^{\left(n+\left\|m\right\|\right)/2}\left(n+\left\|m\right\|\right)!}{2^{n}\left(\frac{n-\left\|m\right\|}{2}\right)!\left(\frac{n+\left\|m\right\|}{2}\right)!}N_{n}^{m},$
$\displaystyle l_{n}^{m}$	$\displaystyle=$	$\displaystyle\left\{\begin{array}[]{c}0,\quad n+m=2l+1,\quad l=0,1,...\\ 1,\quad n+m=2l,\quad l=0,1,...\end{array}\right..$	(34)

Hence,

\left.S_{n}^{m}\left(\mathbf{r}\right)\right|_{\mathbf{r\in}S_{\infty}}=L_{n}^{m}\rho^{-n-1}e^{im\varphi},\quad\left(\left.\rho\right|_{\mathbf{r\in}S_{\infty}}=\left.r\sin\theta\right|_{\mathbf{r\in}S_{\infty}}=\left.r\right|_{\mathbf{r\in}S_{\infty}}\right)

(35)

This shows that

	$\displaystyle\int_{S_{\infty}}S_{n^{\prime}+1}^{m^{\prime}}\left(\mathbf{r}\right)S_{n}^{m}\left(\mathbf{r}\right)dS\left(\mathbf{r}\right)$	$\displaystyle=$	$\displaystyle L_{n^{\prime}+1}^{m^{\prime}}L_{n}^{m}\int_{R}^{\infty}\rho^{-n-n^{\prime}-3}\rho d\rho\int_{0}^{2\pi}e^{i\left(m+m^{\prime}\right)\varphi}d\varphi$
		$\displaystyle=$	$\displaystyle 2\pi\delta_{m^{\prime},-m}\frac{L_{n^{\prime}+1}^{m^{\prime}}L_{n}^{m}}{n+n^{\prime}+1}R^{-n-n^{\prime}-1}$

where $\delta_{m^{\prime}m}$ is the Kronecker delta.

Using Eqs (24), (30), (35), and (3.3) we obtain from Eq. (16)

	$\displaystyle K^{(D)}\left(\mathbf{y,x}\right)$	$\displaystyle=$	$\displaystyle\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\sum_{n^{\prime}=0}^{\infty}\sum_{m^{\prime}=-n^{\prime}}^{n^{\prime}}\frac{2a_{n^{\prime}}^{m^{\prime}}}{2n+1}R_{n}^{-m}\left(\mathbf{x}\right)R_{n^{\prime}}^{-m^{\prime}}\left(\mathbf{y}\right)\int_{S_{\infty}}S_{n^{\prime}+1}^{m^{\prime}}\left(\mathbf{x}^{\prime}\right)S_{n}^{m}\left(\mathbf{x}^{\prime}\right)dS\left(\mathbf{x}^{\prime}\right)$
		$\displaystyle=$	$\displaystyle\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\sum_{n^{\prime}=0}^{\infty}I_{n^{\prime}n}^{m}R_{n}^{-m}\left(\mathbf{x}\right)R_{n^{\prime}}^{m}\left(\mathbf{y}\right)=\sum_{m=-\infty}^{\infty}\sum_{n=\left\|m\right\|+1}^{\infty}\sum_{n^{\prime}=\left\|m\right\|}^{\infty}I_{nn^{\prime}}^{m}R_{n^{\prime}}^{-m}\left(\mathbf{x}\right)R_{n}^{m}\left(\mathbf{y}\right)$

where

I_{nn^{\prime}}^{m}=J_{nn^{\prime}}^{m}R^{-n-n^{\prime}-1},\quad J_{nn^{\prime}}^{m}=\frac{4\pi a_{n}^{m}L_{n+1}^{m}L_{n^{\prime}}^{m}}{\left(2n^{\prime}+1\right)\left(n+n^{\prime}+1\right)}.

(38)

Here we used symmetries $a_{n}^{m}=a_{n}^{-m}$ and $L_{n}^{m}=L_{n}^{-m}$ . Summation in the last sum of Eq (3.3) can start at $n=\left|m\right|$ , but at this value $L_{n+1}^{m}=0$ (see Eq. (31)). Note that $I_{nn^{\prime}}^{m}=I_{nn^{\prime}}^{-m}$ , but the coefficients are non-symmetric with respect to the lower indices, $I_{nn^{\prime}}^{m}\neq I_{n^{\prime}n}^{m}$ . Indeed, we have

I_{n^{\prime}n}^{m}I_{nn^{\prime}}^{m}=\frac{\left(4\pi\right)^{2}a_{n}^{m}a_{n^{\prime}}^{m}R^{-2n-2n^{\prime}-2}}{\left(2n+1\right)\left(2n^{\prime}+1\right)\left(n+n^{\prime}+1\right)^{2}}L_{n}^{m}L_{n+1}^{m}L_{n^{\prime}}^{m}L_{n^{\prime}+1}^{m}=0.

(39)

This is due to the fact $L_{n}^{m}L_{n+1}^{m}=0$ (see definition (31)). Hence, if at some $n^{\prime}$ and $n$ we have $I_{n^{\prime}n}^{m}\neq 0$ (and this is true as not all $I_{n^{\prime}n}^{m}$ are zeros), then due to Eq. (39) $I_{nn^{\prime}}^{m}=0$ , means $I_{n^{\prime}n}^{m}\neq$ $I_{nn^{\prime}}^{m}$ . Non-symmetry of the coefficients also shows that $K^{(D)}\left(\mathbf{y,x}\right)$ is non-symmetric, $K^{(D)}\left(\mathbf{y,x}\right)\neq K^{(D)}\left(\mathbf{x,y}\right)$ since the spherical harmonics form a complete orthogonal basis on a sphere (and $R_{n}^{m}\left(\mathbf{y}\right)$ are proportional to them).

4 Numerical methods

4.1 Dimensionless forms

Since the Laplace equation does not have intrinsic length scales, we have for arbitrary $R$ (see Eqs (23) and (3.3))

$\displaystyle K^{(D)}\left(\mathbf{y},\mathbf{x;}R\right)$	$\displaystyle=$	$\displaystyle\frac{1}{R}\widetilde{K}^{(D)}\left(\frac{\mathbf{y}}{R},\frac{\mathbf{x}}{R}\right),$	(40)
$\displaystyle\widetilde{K}^{(D)}\left(\mathbf{y},\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle-\frac{z_{y}}{8\pi^{2}}\int_{0}^{2\pi}\int_{0}^{1}\frac{\eta d\eta d\varphi^{\prime}}{Q\left(\eta,\varphi^{\prime}-\varphi_{y},\rho_{y},r_{y},1\right)^{3}Q\left(\eta,\varphi^{\prime}-\varphi_{x},\rho_{x},r_{x},1\right)},$
$\displaystyle\widetilde{K}^{(D)}\left(\mathbf{y},\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle\sum_{n=0}^{\infty}\sum_{m=-n}^{n}U_{n}^{m}\left(\mathbf{x}\right)R_{n}^{m}\left(\mathbf{y}\right),\quad U_{n}^{m}\left(\mathbf{x}\right)=\sum_{n^{\prime}=\left\|m\right\|}^{\infty}J_{nn^{\prime}}^{m}R_{n^{\prime}}^{-m}\left(\mathbf{x}\right),\quad\left(\left\|\mathbf{y}\right\|<1,\left\|\mathbf{x}\right\|<1\right).$

where $\widetilde{K}^{(D)}$ is the dimensionless function, which does not depend on parameter $R$ , which is seen from its integral representation. We put an extra argument $R$ to the dimensional function $K^{(D)}$ to show its parametric dependence on $R$ . The series for $\widetilde{K}^{(D)}\left(\mathbf{y},\mathbf{x}\right)$ and $U_{n}^{m}\left(\mathbf{x}\right)$ converge for $\mathbf{x}$ and $\mathbf{y}$ located inside the unit ball.

4.2 Improvement of convergence

Our numerical experiments show that computation of the kernel function via double integral representation (40) provides accurate results, but such computations are somehow slow as they require evaluation of integrals for every pair $\mathbf{x}$ and $\mathbf{y}$ . On the other hand, recursive computations of the spherical basis functions provide a fast procedure.

For practical use all infinite series should be truncated with some truncation number $p$ ( $n<p,$ $n^{\prime}<p$ ), which should be selected based on the error bounds. The radius of convergence of series (40) is $1$ both, with respect to $\mathbf{x}$ and $\mathbf{y}$ . This means that for points in $\Omega_{0}$ located closely to the boundary with domain $\Omega_{\infty}$ the truncation number can be large, which may reduce the efficiency of computations.

Several tricks can be proposed to improve the convergence. First, one can simply implement integral and series representations of the kernel and switch between algorithms based on the values $\left|\mathbf{y}\right|$ and $\left|\mathbf{x}\right|$ . Second, we can effectively reduce the size of the domain, where $K^{(D)}\left(\mathbf{y},\mathbf{x}\right)$ should be computed. This is based on the observation that for the evaluation points located on plane $z=0$ we have from Eqs (25), (26), (31), and (38)

J_{nn^{\prime}}^{m}\left.R_{n}^{m}\left(\mathbf{y}\right)\right|_{z_{y}=0}=\frac{4\pi a_{n}^{m}L_{n^{\prime}}^{m}}{\left(2n^{\prime}+1\right)\left(n+n^{\prime}+1\right)}\left(L_{n+1}^{m}L_{n}^{m}\right)\rho_{y}^{n}e^{-im\varphi_{y}}=0,

(41)

due to $L_{n+1}^{m}L_{n}^{m}=0$ . Hence,

K^{(D)}\left(\mathbf{y},\mathbf{x;}R\right)=0,\quad z_{y}=0.

(42)

Also, relation (42) follows from the integral representation (20).

Consider now an extended computational domain, which can be denoted as $\Omega_{e}$ , and which includes $\Omega_{0}$ and the part of $\Omega_{\infty}$ located between the concentric spheres of radii $R_{0}$ and $R_{e}>R_{0}.$ Respectively, $S_{ge}$ denotes the part of ground surface $S_{g}$ in domain $\Omega_{e}$ (see Fig. 1). All formulae derived hold in case if we replace $\Omega_{0}$ with $\Omega_{e}$ , $R_{0}$ with $R_{e}$ and $S_{g0}$ with $S_{ge}$ . So, instead of Eq. (6) we should have

\phi\left(\mathbf{y}\right)=\int_{S\cup S_{ge}}\sigma_{e}^{(D)}\left(\mathbf{x}\right)\left[G\left(\mathbf{y},\mathbf{x}\right)+K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{e}\right)\right]dS\left(\mathbf{x}\right).

(43)

The charge density $\sigma_{e}^{(D)}\left(\mathbf{x}\right)$ can be determined from the boundary conditions

	$\displaystyle\int_{S\cup S_{ge}}\sigma_{e}^{(D)}\left(\mathbf{x}\right)\left[G\left(\mathbf{y},\mathbf{x}\right)+K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{e}\right)\right]dS\left(\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle\left\{\begin{array}[]{c}V\left(\mathbf{y}\right),\quad\mathbf{y}\in S,\\ 0,\quad\mathbf{y}\in S_{g0},\end{array}\right.,$		(46)
	$\displaystyle\int_{S\cup S_{ge}}\sigma_{e}^{(D)}\left(\mathbf{x}\right)G\left(\mathbf{y},\mathbf{x}\right)dS\left(\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle 0,\quad\mathbf{y}\in S_{ge}\backslash S_{g0}.$

The latter holds since $K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{e}\right)=0,$ for $\mathbf{y}\in S_{ge}\backslash S_{g0}$ . Hence, both for the field computations (43) and determination of the charge density only values of $K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{e}\right)$ for $\mathbf{y}\in\Omega_{0}$ are needed. On the other hand, $K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{e}\right)$ should be computed for $\mathbf{x}\in\Omega_{e}$ with a notice that for $\mathbf{x}\in S_{ge}\backslash S_{g0}$ we have $z_{x}=0.$ We developed an efficient recursive procedure for such computations described in the next subsection.

4.3 Computations for the ground points on the extended domain

So, the problem for use of the extended domain is to compute function $U_{n}^{m}\left(\mathbf{x}\right)$ represented in Eq. (40) by infinite series at $z_{x}$ $=0$ , at $\xi\in\left[R_{0}/R_{e},1\right),$ where $\xi=\left|\mathbf{x}\right|$ . Integral representation of $U_{n}^{m}\left(\mathbf{x}\right)$ can be obtained by substituting of Eq. (30) into integral (16). In the dimensionless form we have

\widetilde{K}^{(D)}\left(\mathbf{y,x}\right)=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}U_{n}^{m}\left(\mathbf{x}\right)R_{n}^{m}\left(\mathbf{y}\right)=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}2a_{n}^{m}R_{n}^{m}\left(\mathbf{y}\right)\int_{1}^{\infty}\int_{0}^{2\pi}G\left(\mathbf{x},\mathbf{x}^{\prime}\right)S_{n+1}^{-m}\left(\mathbf{x}^{\prime}\right)dS\left(\mathbf{x}^{\prime}\right).\quad

(47)

So,

\left.U_{n}^{m}\left(\mathbf{x}\right)\right|_{z_{x}=0}=\frac{1}{2\pi}a_{n}^{m}L_{n+1}^{m}e^{-im\varphi_{x}}u_{n}^{m}\left(\xi\right),\quad\xi=\rho_{x},

(48)

	$\displaystyle u_{n}^{m}\left(\xi\right)$	$\displaystyle=$	$\displaystyle\int_{1}^{\infty}\frac{d\eta}{\eta^{n+1}}\int_{0}^{2\pi}\frac{e^{im\varphi}d\varphi}{\left(\eta^{2}-2\xi\eta\cos\varphi+\xi^{2}\right)^{1/2}}=\frac{1}{\xi^{n+1}}\int_{0}^{\xi}\zeta^{n}w_{m}\left(\zeta\right)d\zeta,$		(49)
	$\displaystyle w_{m}\left(\xi\right)$	$\displaystyle=$	$\displaystyle\int_{0}^{2\pi}\frac{e^{im\varphi}d\varphi}{\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}},\quad u_{n}^{-m}\left(\xi\right)=u_{n}^{m}\left(\xi\right),\quad w_{-m}\left(\xi\right)=w_{m}\left(\xi\right).$

These functions can be expressed via the complete elliptic integrals of the first and second kinds, $K$ and $E$ (see[21]), and recursive relations (to shorten notation we drop the argument $\xi$ for functions $u_{n}^{m}\left(\xi\right)$ and $w_{m}\left(\xi\right)$ )

$\displaystyle u_{n}^{0}$	$\displaystyle=$	$\displaystyle\frac{1}{n^{2}\xi^{2}}\left[4E\left(\xi^{2}\right)-4n\left(1-\xi^{2}\right)K\left(\xi^{2}\right)+\left(n-1\right)^{2}u_{n-2}^{0}\right],\quad n=1,3,5,...$	(50)
$\displaystyle u_{n+1}^{1}$	$\displaystyle=$	$\displaystyle\frac{1}{\left(2n+3\right)\xi}\left[\left(n+1\right)u_{n}^{0}+\left(n+2\right)\xi^{2}u_{n+2}^{0}+4\left(1-\xi^{2}\right)K\left(\xi^{2}\right)-8E\left(\xi^{2}\right)\right],\quad n=1,3,5,...$	(51)
$\displaystyle u_{n+1}^{m+1}$	$\displaystyle=$	$\displaystyle\frac{2}{\left(2n+3\right)\xi}\left[\left(n+1\right)u_{n}^{m}+\left(n+2\right)\xi^{2}u_{n+2}^{m}-v_{m}\right]-u_{n+1}^{m-1},\text{ }m=1,2,...,\,n=m+1,m+3,...$	(52)
$\displaystyle v_{m}$	$\displaystyle=$	$\displaystyle\left(1+\xi^{2}\right)w_{m}-\xi\left(w_{m+1}+w_{m-1}\right),\quad m=1,2,...$	(53)
$\displaystyle w_{m}$	$\displaystyle=$	$\displaystyle\frac{1}{\xi}\left(1+\xi^{2}\right)\frac{2m-2}{2m-1}w_{m-1}-\frac{2m-3}{2m-1}w_{m-2},\quad m=2,3,...,$	(54)
$\displaystyle w_{0}$	$\displaystyle=$	$\displaystyle 4K\left(\xi^{2}\right),\quad w_{1}=\frac{4}{\xi}\left[K\left(\xi^{2}\right)-E\left(\xi^{2}\right)\right],$	(55)
$\displaystyle K\left(\mu\right)$	$\displaystyle=$	$\displaystyle\int_{0}^{\pi/2}\left(1-\mu\sin^{2}\varphi\right)^{-1/2}d\varphi,\quad E\left(\mu\right)=\int_{0}^{\pi/2}\left(1-\mu\sin^{2}\varphi\right)^{1/2}d\varphi.$	(56)

Appendix A describes how these relations are obtained. We remark that because $L_{n+1}^{m}=0$ when $n+m$ is an even number, to obtain $U_{n}^{m}$ in Eq. (48) we need to compute functions $u_{n}^{m}$ only at odd values of $n+m$ , which is reflected in the above recurrences. The recurrence for $u_{n}^{m}$ works in layers in subscript $n$ . So, first, we obtain $u_{n}^{0}$ for all $n$ needed, then $u_{n}^{1}$ , and so on until $u_{n}^{m}$ at maximum $m$ needed. For truncation number $p$ , we have max $\left(n\right)=p-1$ , max $\left(m\right)=p-2$ . As recurrence (52) show that the maximum computed $n$ at $m=p-3$ should be $p$ , while at $m=0$ it should be $2p-3$ . It is also remarkable that recurrence (50) can be used first at $n=1$ , which provides the value of $u_{1}^{0}$ , and then applied for other $u_{n}^{0}$ in the sequence.

4.4 Boundary elements

Equation (46) can be solved using the method of moments (MoM) or collocation boundary element methods (BEM). In the present study we use the BEM. The surface is discretized by triangular elements, the integrals involving Green’s function are computed analytically assuming that the charge density is constant over the elements (constant panel approximation). The integrals involving kernel $K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{c}\right)$ were computed using center panel quadrature,

\int_{\Delta_{j}}\sigma\left(\mathbf{x}\right)K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{c}\right)dS\left(\mathbf{x}\right)\approx\sigma_{j}w_{j}K^{(D)}\left(\mathbf{y},\mathbf{x}_{j}^{c}\mathbf{;}R_{c}\right),

(57)

where $\sigma_{j}$ is the mean value of the charge density, while $w_{j}$ and $\mathbf{x}_{j}^{c}$ are the area and the centroid of triangle $\Delta_{j}$ . In the standard BEM realization the integrals over the kernels with Green’s functions were computed analytically.

\int_{\Delta_{j}}\sigma\left(\mathbf{x}\right)G\left(\mathbf{y},\mathbf{x}\right)dS\left(\mathbf{x}\right)\approx\sigma_{j}L_{j}\left(\mathbf{y}\right),

(58)

$\displaystyle L_{j}\left(\mathbf{y}\right)$	$\displaystyle=$	$\displaystyle\int_{\Delta_{j}}G\left(\mathbf{y},\mathbf{x}\right)dS\left(\mathbf{x}\right)=\sum_{q=1}^{3}\left[H\left(l_{jq}-x_{jq},h_{j},z_{jq}\right)-H\left(-x_{jq},h_{j},z_{jq}\right)\right],$	(59)
$\displaystyle x_{jq}$	$\displaystyle=$	$\displaystyle\left(\mathbf{y-x}_{jq}\right)\cdot\mathbf{i}_{jq},\quad h_{j}=\left\|\left(\mathbf{y-x}_{j1}\right)\cdot\mathbf{n}_{j}\right\|,\quad z_{jq}=\left(\mathbf{y-x}_{jq}\right)\cdot\mathbf{n}_{jq},$
$\displaystyle\mathbf{i}_{jq}$	$\displaystyle=$	$\displaystyle\frac{1}{l_{jq}}\left(\mathbf{x}_{j,q(\mathop{\mathrm{m}od}3)+1}-\mathbf{x}_{jq}\right),\quad l_{jq}=\left\|\mathbf{x}_{j,q(\mathop{\mathrm{m}od}3)+1}-\mathbf{x}_{jq}\right\|,\quad\mathbf{n}_{jq}=\mathbf{i}_{jq}\times\mathbf{n}_{j},$
$\displaystyle\quad H\left(x,y,z\right)$	$\displaystyle=$	$\displaystyle y\left(\arctan\frac{x}{z}-\arctan\frac{yx}{zr}\right)-z\ln\left\|r+x\right\|,\quad r=\sqrt{x^{2}+y^{2}+z^{2}},$

where $\mathbf{n}_{j}$ is the normal to triangle $\Delta_{j}$ directed into the computational domain.

4.5 Fast multipole method

For large problems the BEM can be accelerated using the fast multipole method (FMM) (e.g., [22]). As shown, the kernel $K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{c}\right)$ is either zero or can be factored for any pair of points inside the computational domain of radius $R_{c}$ . This shows that boundary integral equation (46) can be written as

	$\displaystyle\int_{S\cup S_{ge}}\!\!\!\!\!\!\!\!\!\!\sigma_{e}^{(D)}\left(\mathbf{x}\right)G\left(\mathbf{y},\mathbf{x}\right)dS\left(\mathbf{x}\right)+$		(60)
	$\displaystyle\frac{1}{R_{e}}\sum_{n=0}^{p-1}\!\sum_{m=-n}^{n}R_{n}^{m}\left(\frac{\mathbf{y}}{R_{e}}\right)\int_{S\cup S_{ge}}\!\!\!\!\!\!\!\!\!\!\sigma_{e}^{(D)}\left(\mathbf{x}\right)U_{n}^{m}\left(\frac{\mathbf{x}}{R_{e}}\right)dS\left(\mathbf{x}\right)=\left\{\begin{array}[]{c}V\left(\mathbf{y}\right),\quad\mathbf{y}\in S,\\ 0,\quad\mathbf{y}\in S_{ge},\end{array}\right.$		(63)

where the infinite sum over $n$ is truncated. The first integral in the left hand side can be computed using the FMM with the free space Green’s function. In our implementation, we used some nearfield radius $r_{nf}$ . So, for pairs of points $\left(\mathbf{y},\mathbf{x}\right)$ located closer than $r_{nf}$ the integrals over triangles were computed analytically using Eqs (58) and (59), while for larger distances the integrals were computed using center point approximation, similarly to (57).

The most important feature of the global kernel factorization is that the number of operations to compute the matrix-vector product for the kernel is $O\left(p^{2}N\right)$ , not $O\left(N^{2}\right)$ , since all integrals over $\mathbf{x}$ can be computed independently from the evaluation points $\mathbf{y}$ ( $N$ is the number of boundary elements). So, computation of this sum has the asymptotic complexity consistent with the FMM, and the method can be applied for solution of large problems, where the overall linear system is solved iteratively.

4.6 Real basis functions

Kernel $\widetilde{K}^{(D)}$ is real, but is presented in Eq. (40) via sums of complex valued functions $R_{n}^{m}$ . A fast procedure can be developed based on recursive computations of real basis functions (see [23]),

\underline{R_{n}^{m}}\left(\mathbf{r}\right)=\frac{\left(-1\right)^{n+m}}{(n+\left|m\right|)!}r^{n}P_{n}^{\left|m\right|}(\cos\theta)\left\{\begin{array}[]{c}\cos m\varphi,\quad m\geqslant 0,\\ \sin m\varphi,\quad m<0,\end{array}\right.

(64)

which can be computed recursively as

$\displaystyle\underline{R_{0}^{0}}$	$\displaystyle=$	$\displaystyle 1,\quad\underline{R_{1}^{1}}=-\frac{1}{2}x,\quad\underline{R_{1}^{-1}}=\frac{1}{2}y,$	(65)
$\displaystyle\underline{R_{m}^{m}}$	$\displaystyle=$	$\displaystyle-\frac{1}{2m}\left(x\underline{R_{m-1}^{m-1}}+y\underline{R_{m-1}^{-\left(m-1\right)}}\right),\quad m=2,3,...,$
$\displaystyle\quad\underline{R_{m}^{-m}}$	$\displaystyle=$	$\displaystyle\frac{1}{2m}\left(y\underline{R_{m-1}^{m-1}}-x\underline{R_{m-1}^{-\left(m-1\right)}}\right),\quad m=2,3,...,$
$\displaystyle\underline{R_{m+1}^{\pm m}}$	$\displaystyle=$	$\displaystyle-z\underline{R_{m}^{\pm m}},\quad m=0,1,...,$
$\displaystyle\underline{R_{n+1}^{\pm m}}$	$\displaystyle=$	$\displaystyle-\frac{\left(2n+1\right)z\underline{R_{n}^{\pm m}}+r^{2}\underline{R_{n-1}^{\pm m}}}{\left(n+1\right)^{2}-m^{2}},\quad m=0,1,...,\text{ }n=m+1,m+2,....\quad$

Equations (40), (25), (26), (29), (31),and (38) show then that

	$\displaystyle\widetilde{K}^{(D)}\left(\mathbf{y},\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\underline{U_{n}^{m}}\left(\mathbf{x}\right)\underline{R_{n}^{m}}\left(\mathbf{y}\right),\quad\underline{U_{n}^{m}}\left(\mathbf{x}\right)=-\frac{2-\delta_{m,0}}{4\pi}\nu_{n+1}^{m}\sum_{n^{\prime}=\left\|m\right\|}^{\infty}\frac{\nu_{n^{\prime}}^{m}}{n^{\prime}+n+1}\underline{R_{n^{\prime}}^{m}}\left(\mathbf{x}\right),$		(66)
	$\displaystyle\nu_{n}^{m}$	$\displaystyle=$	$\displaystyle l_{n}^{m}\left(-1\right)^{\left(n+\left\|m\right\|\right)/2}(n-\left\|m\right\|-1)!!(n+\left\|m\right\|-1)!!.$

Also, using Eq. (48) we have for points $z_{x}=0$

\underline{U_{n}^{m}}\left(\mathbf{x}\right)=-\frac{2-\delta_{m,0}}{8\pi^{2}}\nu_{n+1}^{m}u_{n}^{m}\left(\xi\right)\left\{\begin{array}[]{c}\cos\left(m\varphi_{x}\right),\quad m\geqslant 0,\\ \sin\left(m\varphi_{x}\right),\quad m<0.\end{array}\right.

(67)

5 Complexity and optimization

5.1 Kernel computations

Assume there are $N$ sources and $M$ receivers in domain $\Omega_{0}$ . Truncation of the series representation of $K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{e}\right)$ for $\mathbf{y\in}$ $\Omega_{0}$ with truncation number $p$ produces errors $\epsilon\sim(R_{0}/R_{e})^{p}$ . If there are $N_{e}$ sources in the extended domain $\Omega_{e}$ then the computational cost for factored $K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{e}\right)$ can be estimated as

C_{fact}=O\left(Mp^{2}\right)+O\left(Np^{3}\right)+O\left((N_{e}-N)p^{2}\right),

(68)

where the first term is the cost to compute $p^{2}$ functions $R_{n}^{m}$ $\left(\mathbf{y/}R_{e}\right)$ , the second term is the cost to compute $p^{2}$ functions $U_{n}^{m}$ $\left(\mathbf{x/}R_{e}\right)$ for $\mathbf{x}\in\Omega_{0}$ using series (66), and the third term is the cost to compute $p^{2}$ functions $U_{n}^{m}$ $\left(\mathbf{x/}R_{e}\right)$ for $\mathbf{x}\in\Omega_{e}$ $\backslash$ $\Omega_{0}$ $\left(z_{x}=0\right)$ using Eq. (67) and recursions (50)-(56) for $u_{n}^{m}\left(\xi\right)$ . Assume now that the sources are uniformly distributed in domain $\Omega_{e}$ $\backslash$ $\Omega_{0}$ , with the density $\rho$ points per surface area. In this case $N_{e}-N=\pi\rho\left(R_{e}^{2}-R_{0}^{2}\right)$ . For prescribed accuracy $\epsilon$ we have $R_{e}=R_{0}\epsilon^{-1/p}$ , which brings the following expression for the computational cost

C_{fact}=AMp^{2}+BNp^{3}+C\pi\rho R_{0}^{2}\left(\epsilon^{-2/p}-1\right)p^{2},

(69)

where $A,B$ , and $C$ are some asymptotic constants.

Considering $C_{fact}$ as a function of $p$ , we can see that

\frac{dC_{fact}}{dp}=2AMp+3BNp^{2}+2C\pi\rho R_{0}^{2}p\left[\left(\epsilon^{-2/p}-1\right)-\epsilon^{-2/p}\frac{1}{p}\ln\frac{1}{\epsilon}\right].

(70)

The term in the square brackets has a single zero at $p=p_{c}\left(\epsilon\right)$ and is positive at

p>p_{c}\left(\epsilon\right)=\frac{1}{\alpha}\ln\frac{1}{\epsilon},\quad\alpha=0.79681213...\approx 0.8.

(71)

So, function $C_{fact}\left(p\right)$ grows monotonically at $p>p_{c}\left(\epsilon\right)$ for any values of $M,N,$ and $\pi\rho R_{0}^{2}.$ The function can grow or decay at $p<p_{c}\left(\epsilon\right)$ depending on $M,N,$ and $\pi\rho R_{0}^{2}$ . Note now, that $p=p_{c}\left(\epsilon\right)$ corresponds to

\left(\frac{R_{e}}{R_{0}}\right)_{c}=\epsilon^{-1/p_{c}}=e^{\alpha}=\beta\approx 2.2.

(72)

This means, that $C_{fact}$ being considered as a function of $R_{e}$ decays in the region $R_{0}<R_{e}\leqslant\beta R_{0}$ and may continue to decay or start to grow at $R_{e}>\beta R_{0}$ . It is remarkable, that $\beta$ does not depend on $\epsilon$ , so for computations using the factored kernel the optimal radius of the extended domain should be always larger than $\beta R_{0}$ . Also, it can be noted that the values of $p_{c}\left(\epsilon\right)$ are rather small. For example, at $\epsilon=10^{-4}$ we have $p_{c}\left(\epsilon\right)\approx$ $11.6$ , which means that $p=12$ should be sufficient to achieve the required accuracy at $R_{e}=\beta R_{0}$ . This also shows that at large enough $N$ and $M$ computational costs in Eq.(69) are much smaller than that for kernel computations using the integral representation, which formally is $O\left(NM\right)$ with the asymptotic constant that can be substantially large as the cost of numerical quadrature also depends on the required accuracy.

To validate these findings we conducted some numerical experiments. First, we did error estimates of the factored solutions using comparisons with the error-controlled integral computations (Matlab integral2 function). For this purpose due to the symmetry of the problem, all receivers can be placed on the intersection of a unit hemisphere (“bump”) with plane $y=0$ . On the other hand, the sources can be distributed over the hemisphere surface and in the extension of the domain (see Fig. 2). Figure 3 shows the relative $L_{2}$ -norm errors, $\epsilon_{2}$ ,

\epsilon_{2}=\frac{\left\|f-f_{ref}\right\|_{2}}{\left\|f_{ref}\right\|_{2}},

(73)

where $f$ is the computed and $f_{ref}$ is the reference solution for all computed source-receiver pairs as the filled contours in plane $\left[R_{e}/R_{0},p\right]$ . It is seen that the boundaries of the regions agree well with the error approximation $R_{e}=R_{0}\epsilon_{2}^{-1/p}$ . Figure 4 is computed by finding for each $R_{e}/R_{0}$ the smallest $p$ at which the computed $\epsilon_{2}$ is smaller than the prescribed accuracy and plotting the wall-clock time required for kernel computation via factorization. It is seen that dependences computed at different prescribed $\epsilon_{2}$ have global minima in the range $2\leqslant R_{e}/R_{0}\leqslant 3.5$ at $10^{-2}\geqslant\epsilon_{2}\geqslant 10^{-8}$ . This agrees well with the theoretical estimates provided above. Note that non-monotonous character of the curves shown on this plot is due to discrete change of parameters $R_{e}/R_{0}$ and $p$ and thresholding (the curves show the prescribed error curves, not the actual errors, which can be close to or several times smaller than the prescribed errors). Also, this figure shows that optimization of the extended domain size can speed up computations several times.

5.2 Optimization of the Boundary Element Solution

The extension of the domain can substantially impact the complexity and memory requirements of the BEM and it is preferable not to extend the domain too much. Assume that $M=N$ and $R_{e}=R_{0}\left(1+\delta\right)$ , $\delta\ll 1$ . In this case, $p\sim\delta^{-1}\ln\left(1/\epsilon\right)$ , the first term in Eq. (69) can be neglected, and we obtain

C_{fact}\sim B\frac{N}{\delta^{3}}\ln^{3}\frac{1}{\epsilon}\left(1+\frac{\pi R_{0}^{2}}{A_{0}}\frac{2C}{B}\frac{\delta^{2}}{\ln\frac{1}{\epsilon}}\right)\sim B\frac{N}{\delta^{3}}\ln^{3}\frac{1}{\epsilon}.

(74)

Here we assumed that the density of sampling of the domain extension is the same as the density of sampling of surface $S\cup S_{g0}$ of area $A_{0}$ , means $\rho=N/A_{0}$ . The last term in the parentheses then can be neglected at $\delta\rightarrow 0$ , due to $\pi R_{0}^{2}\leqslant A_{0}$ , $\ln\frac{1}{\epsilon}>1,$ and assumed $C\sim B.$

There are different versions of the BEM including those using iterative methods and the FMM, which have different computational costs, e.g., $O\left(N^{\alpha}\right)$ , $\alpha\geqslant 1$ , for problems of size $N.$ So, the cost of the BEM with the extended domain kernel computing at $\delta\ll 1$ can be estimated as

	$\displaystyle C_{BEM}$	$\displaystyle=$	$\displaystyle C_{fact}+DN_{e}^{\alpha}=C_{fact}+DN^{\alpha}\left(1+\frac{N_{e}-N}{N}\right)^{\alpha}$
		$\displaystyle\sim$	$\displaystyle C_{fact}+DN^{\alpha}\left(1+\alpha\frac{\pi R_{0}^{2}}{A_{0}}\left(\frac{R_{e}^{2}}{R_{0}^{2}}-1\right)\right)\sim DN^{\alpha}+B\frac{N}{\delta^{3}}\ln^{3}\frac{1}{\epsilon}+2DN^{\alpha}\alpha\frac{\pi R_{0}^{2}}{A_{0}}\delta.$

Here $D$ is some asymptotic constant related to the BEM cost. While it appears that the cost of the last term is much smaller than the cost of the first term, we retained it, as the first term is constant, while the last term reflects the main mechanism of increase of the cost at the increasing $\delta$ . Indeed, taking derivative of this expression we can see that $C_{BEM}$ has a global minimum,

	$\displaystyle\frac{dC_{BEM}}{d\delta}$	$\displaystyle=$	$\displaystyle-3B\frac{N}{\delta^{4}}\ln^{3}\frac{1}{\epsilon}+2\alpha D\frac{\pi R_{0}^{2}}{A_{0}}N^{\alpha},\quad\frac{dC_{BEM}}{d\delta}\left(\delta_{opt}\right)=0,\quad$		(76)
	$\displaystyle\delta_{opt}$	$\displaystyle=$	$\displaystyle\frac{3B}{2\alpha D}\frac{A_{0}}{\pi R_{0}^{2}}\left(\frac{1}{N^{\alpha-1}}\ln^{3}\frac{1}{\epsilon}\right)^{1/4}=O\left(\frac{1}{N^{\left(\alpha-1\right)/4}}\ln^{3/4}\frac{1}{\epsilon}\right).$

This shows that, indeed, at $\alpha>1$ and $N\rightarrow\infty$ and fixed other parameters the optimal size of the extension, $\delta_{opt}$ , turns to zero. For the regular BEM with exact solvers we have $\alpha=3$ . However, when $N$ is not very high the cost of the FMM can be limited by the cost of boundary integral computing, in which case $\alpha=2$ . Also, this can be achieved with iterative solvers. In the case of the FMM acceleration $\alpha=1$ or close to 1. In this case, $\delta_{opt}$ should not depend on $N$ , and the only reason why $\delta_{opt}\ll 1$ can be related to large values of asymptotic constant $D$ . Usually, this is true for the FMM.

6 BEM computing

For illustrations of the kernel can be used with the boundary elements we consider the following benchmark problems for validation. We assumed that the ground elevation $z$ is the following function of $x$ and $y$ ,

S_{g}:\left\{(x,y,z),\quad z=\left\{\begin{array}[]{c}\pm\sqrt{1-\rho^{2}},\quad\rho<1,\\ 0,\quad\rho\geqslant 1.\end{array}\right.,\quad\rho=\sqrt{x^{2}+y^{2}}\right\}.

(77)

Here the sign “+” corresponds to a “bump”, while the sign “-” corresponds to a “dip”. Assume then that the electric field is generated by a monopole source (charge) of unit intensity located at $\mathbf{x}_{s}=\left(0,0,h\right)$ , and the domain $\Omega_{0}$ is bounded by a sphere of radius $R_{0}=h$ ( $h>1$ ) centered at the origin for the “bump” and , $R_{0}=1$ for the “dip” ( $\left|h\right|<1$ ). Figure 5 shows the meshes used in computations.

For validation purposes, the solution of the problem for the bump can be found analytically using the method of images. In this case, the image of the bump is a hemisphere, which together with the “image” bump forms a sphere of radius 1 on which the potential is zero. The image source of negative intensity (sink) is located at $\mathbf{x}_{s}^{\ast}=\left(0,0,-h\right)$ , and we have the following expression for the potential

\phi\left(\mathbf{y}\right)=G\left(\mathbf{y},\mathbf{x}_{s}\right)-G\left(\mathbf{y},\mathbf{x}_{s}^{\ast}\right)-\frac{1}{h}G\left(\mathbf{y},\mathbf{x}_{r}^{+}\right)+\frac{1}{h}G\left(\mathbf{y},\mathbf{x}_{r}^{-}\right),\quad\mathbf{x}_{r}^{\pm}=\left(0,0,\pm\frac{1}{h}\right),

(78)

where $\mathbf{x}_{r}^{+}$ and $\mathbf{x}_{r}^{-}$ are the images of $\mathbf{x}_{s}$ and $\mathbf{x}_{s}^{\ast}$ with respect to the sphere. We are not aware about analytical solutions for the case of “dip” for which the method of images cannot be applied.

The boundary integral solution of the problem can be represented in the form

\phi\left(\mathbf{y}\right)=G\left(\mathbf{y},\mathbf{x}_{s}\right)+K^{(D)}\left(\mathbf{y},\mathbf{x}_{s}\mathbf{;}R_{e}\right)+\int_{S_{ge}}\sigma^{(D)}\left(\mathbf{x}\right)\left[G\left(\mathbf{y}-\mathbf{x}\right)+K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{e}\right)\right]dS\left(\mathbf{x}\right),

(79)

where the charge density can be found from

\int_{S_{ge}}\sigma^{(D)}\left(\mathbf{x}\right)\left[G\left(\mathbf{y},\mathbf{x}\right)+K^{(D)}\left(\mathbf{y},\mathbf{x;}R_{e}\right)\right]dS\left(\mathbf{x}\right)=-G\left(\mathbf{y},\mathbf{x}_{s}\right)-K^{(D)}\left(\mathbf{y},\mathbf{x}_{s}\mathbf{;}R_{e}\right),\quad\mathbf{y}\in S_{ge}.

(80)

In the case of “bump” we computed and compared four solutions: 1) analytical; 2) regular BEM using method of images (two symmetric sources and a sphere) (“BEMimage”); 3) regular BEM, where the domain ground plane was simply truncated and the effect of the infinite plane ignored (“BEM”); 4) BEM with the kernel accounting for the infinite plane as described by Eqs (79) and (80) (“BEMinf”). In case of “dip” only solutions 3) and 4) were computed. In all cases we varied $R_{e}$ . The truncation number $p$ for “BEMinf” method was selected based on the theoretical estimate, $p=\log\epsilon/\log(R_{0}/R_{e})$ , where the prescribed error $\epsilon$ varied in range $10^{-3}-10^{-6}$ .

Figure 6 illustrates the effect of the domain extension on the performance for different prescribed accuracies and the BEM complexities for the “bump” case ( $h=2$ ). Note that here the FMM accelerated BEM was implemented in the way that the boundary integrals were computed within some cutoff radius of the evaluation point using Matlab, while the FMM procedure was executed in the form of a high performance library, so the cost of the FMM BEM is mainly determined by the Matlab integral computations, which is consistent with the kernel evaluation also implemented in Matlab. Also, in the regular BEM the solution time was much smaller than the integral computation time. So, the methods shown can be characterized as $O\left(N^{2}\right)$ for the BEM and $O\left(N\right)$ for the BEM/FMM. It is seen that the curves form minima (sometimes a few local minima), at relatively low $\delta=R_{e}/R_{0}-1$ . According to the theory (see Eq. (76)), increasing the prescribed accuracy results in an increase of the optimal $R_{e}$ , which is seen on the figure. For $\delta$ substantially larger than $\delta_{opt}$ the cost of kernel computation becomes negligible in comparison to the BEM costs, and the curves computed for different $\epsilon$ become close to each other.

Figure 7 shows the analytical solution (78) (the induced part, $\phi_{ind}\left(\mathbf{y}\right)=\phi\left(\mathbf{y}\right)-G\left(\mathbf{y},\mathbf{x}_{s}\right)$ ), and the absolute errors in the computational domain obtained using different methods of solution for the “bump” case, $h=2$ , $\delta=0.0935$ , $N_{0}=6401$ , and $N_{e}=7661$ , where $N_{0}$ and $N_{e}$ are the number of boundary elements in $\Omega_{0}$ and $\Omega_{e}$ , respectively. The numerical solution using image sources is the most accurate and its relative $L_{2}$ -norm error is $\epsilon_{2}=2.3\cdot 10^{-3}$ . The BEMinf shows error $\epsilon_{2}=4.5\cdot 10^{-3}$ (at the prescribed accuracy $\epsilon=10^{-4}$ ) while the BEM shows $\epsilon_{2}=3.7\cdot 10^{-2}$ , which is one order of magnitude larger. This figure shows also different error distributions for different methods inside the computational domain.

Figure 8 shows the reference solution $\left(\phi_{ind}\left(\mathbf{y}\right)\right)$ for the “dip” case and the absolute errors in the computational domain obtained using BEMinf and the regular BEM, $h=0.5$ , $\delta=0.124$ , $N_{0}=1592$ , and $N_{e}=2017$ . The reference solution was computed using BEMinf with parameters $\delta=0.5$ , $N_{0}=6401$ , and $N_{e}=14038$ and prescribed $\epsilon=10^{-6}$ . The BEMinf shows error $\epsilon_{2}=4.7\cdot 10^{-4}$ (at the prescribed accuracy $\epsilon=10^{-4}$ ) while the BEM shows $\epsilon_{2}=5.6\cdot 10^{-2}$ , which is two orders of magnitude larger. Similarly to the case of “bump” one can observe different error distributions for the BEMinf and the BEM. Similarly to the case of “bump” the largest BEMinf errors are observed near the edge of domain $\Omega_{0}$ coinciding with the ground plane.

Figure 9 plots $\epsilon_{2}$ for the “dip” case as a function of parameter $R_{e}/R_{0}$ . It is seen that the error is practically constant for BEMinf, while it decays inversely proportionally to the volume of the computational domain for the regular BEM. It is clear then that the results of BEMinf can be obtained by simple extension of the computational domain, but for substantially higher computational costs. For example, to achieve two orders of magnitude in error reduction one should have $R_{e}\sim 5R_{0}$ , which may require one order of magnitude increase in the number of boundary elements and increase the computational costs 10-1000 times, based on the solution method used (scaled as $O\left(N^{\alpha}\right)$ , $1\leqslant\alpha\leqslant 3$ ).

7 Conclusion

We introduced analogues of Green’s function for an infinite plane with a circular hole. We also proposed and tested efficient methods for computing of these functions, which may have broad applications to model infinite ground or ocean surfaces in the problems requiring Laplace equation solvers. While such functions can be expressed in the forms of finite integrals, their computation can be costly. However, a computationally efficient factorization of these functions can be done, which drastically reduces the computational costs, and, in fact, creates rather small overheads to the problems, which simply ignore the existence of ground outside the computational domain. Also, we conducted some optimization study and provided actual boundary element computations for some benchmark problems. The study shows that the accuracy of the solution using the analogues of Green’s function improves substantially compared to a simple domain truncation, which also can become computationally costly for improved accuracy computations.

8 Acknowledgements

This work is supported by Cooperative Research Agreement (W911NF1420118) between the University of Maryland and the Army Research Laboratory, with David Hull, Ross Adelman and Steven Vinci as Technical monitors.

9 References

References

[1] R. Adelman, Extremely large, wide-area power-line models, Supercomputing 2017, https://sc17.supercomputing.org/SC17%20Archive/tech_poster/poster_files/post251s2-file3.pdf.
[2] M. D’Amore, M. S. Sarto, Simulation models of a dissipative transmission line above a lossy ground for a wide-frequency range. I. Single conductor configuration, IEEE Trans. Electromag. Compat. 38(2) (1996) 127-138. https://doi.org/10.1109/15.494615.
[3] M. Trlep, A. Hamler, M. Jesenik, B. Stumberger, Electric field distribution under transmission lines dependent on ground surface, IEEE Trans. Mag. 45(3) (2009) 1748-1751. https://doi.org/10.1109/TMAG.2009.2012806.
[4] Q. Zhang, J. Yang, D. Li, Z. Wang, Propagation effects of a fractal rough ocean surface on the vertical electric field generated by lightning return strokes, J. Electrostat. 70(1)(2012) 54-59. https://doi.org/10.1016/j.elstat.2011.10.003.
[5] K. S. Chen, T.-D. Wu, L. Tsang, Q. Li, J. Shi, A.K. Fung, Emission of rough surfaces calculated by the integral equation method with comparison to three-dimensional moment method simulations, IEEE Trans. Geo. Rem. Sensing 41(1) (2003) 90-101. https://doi.org/10.1109/TGRS.2002.807587.
[6] D. Givoli, Numerical Methods for Problems in Infinite Domains, Elsevier, Amsterdam, 1992.
[7] S.V. Tsynkov, Numerical solution of problems on unbounded domains. A review, Appl. Num. Math. 27(4) (1998) 465-532. https://doi.org/10.1016/S0168-9274(98)00025-7.
[8] J.D. Jackson, Classical Electrodynamics, third ed., John Wiley and Sons, New York, 1998.
[9] D. Crowdy, J. Marshall, Green’s functions for Laplace’s equation in multiply connected domains, IMA J. Appl. Math. 72 (2007) 278-301. https://doi.org/10.1093/imamat/hxm007.
[10] N.A. Gumerov and R. Duraiswami, A method to compute periodic sums, J. Comput. Phys. 272 (2014) 307-326. https://doi.org/10.1016/j.jcp.2014.04.039.
[11] A. Moroz, Quasi-periodic Green’s functions of the Helmholtz and Laplace equations, J. Phys. A: Math. Gen. 39(36) (2006) 11247. https://doi.org/10.1088%2F0305-4470%2F39%2F36%2F009.
[12] G. Bao, J. Lin, Near-field imaging of the surface displacement on an infinite ground plane, Inverse Problems & Imaging, 7(2) (2013), 377-396. https:/doi.org/10.3934/ipi.2013.7.377.
[13] S. N. Chandler-Wilde, E. Heinemeyer, R. Potthast, Acoustic scattering by mildly rough unbounded surfaces in three dimensions, SIAM J. Appl. Math. 66 (2006) 1002–1026. https://doi.org/10.1137/050635262.
[14] J.J. Holmes, Past, present, and future of underwater sensor arrays to measure the electromagnetic field signatures of naval vessels, Marine Tech. Soc. J. 49(6) (2015) 123-133. https://doi.org/10.4031/MTSJ.49.6.1.
[15] R. Yue, P. Hu, J. Zhang, The influence of the seawater and seabed interface on the underwater low frequency electromagnetic field signatures, 2016 IEEE/OES China Ocean Acoustics (COA), Harbin (2016) 1-7. https://doi.org/10.1109/COA.2016.7535694.
[16] X. Wang, Q. Xu, J. Zhang, Simulating underwater electric field signal of ship using the boundary element method, Progress In Electromag. Res. M 76 (2018) 43-54. https:/doi.org/10.2528/PIERM18092706.
[17] K.V. Rozhdestvensky, Aerodynamics of a Lifting System in Extreme Ground Effect, Springer, Berlin-Heidelberg, 2000.
[18] D. G. Duffy, Green’s Functions with Applications, CRC Press, Taylor & Francis Group, FL, 2015.
[19] P.M. Morse, H. Feshbach, Methods of Theoretical Physics, McGraw-Hill, NY, 1953.
[20] N.A. Gumerov, R. Duraiswamin, Comparison of the efficiency of translation operators used in the fast multipole method for the 3D Laplace equation, UMIACS TR 2005-09, Also issued as Computer Science Technical Report CS-TR-# 4701, University of Maryland, College Park, 2005. https://drum.lib.umd.edu/bitstream/handle/1903/3023/LaplaceTranslation_CSTR_4701.pdf.
[21] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington D.C., 1972.
[22] R. Adelman, N.A Gumerov, and R. Duraiswami, FMM/GPU-accelerated boundary element method for computational magnetics and electrostatics, IEEE Trans. Mag. 53(12) (2017), 7002311. https://doi.org/10.1109/TMAG.2017.2725951.
[23] N.A. Gumerov and R. Duraiswami, Fast multipole methods on graphics processors, J. Comput. Phys. 227 (2008), 8290-8313. https://doi.org/10.1016/j.jcp.2008.05.023.

Appendix A Derivation of Recurrence Relations

Recurrence (52) can be derived by comparing two representations of integral

	$\displaystyle I_{n}^{m}$	$\displaystyle=$	$\displaystyle\frac{1}{\xi^{n+1}}\int_{0}^{2\pi}e^{im\varphi}d\varphi\int_{0}^{\xi}\left(1-2\zeta\cos\varphi+\zeta^{2}\right)^{1/2}\zeta^{n}d\zeta$
		$\displaystyle=$	$\displaystyle u_{n}^{m}+\xi^{2}u_{n+2}^{m}-\xi\left(u_{n+1}^{m-1}+u_{n+1}^{m+1}\right),$

which can be checked using definition of $u_{n}^{m}\left(\xi\right)$ (49), and integration by parts:

	$\displaystyle I_{n}^{m}$	$\displaystyle=$	$\displaystyle\frac{1}{\left(n+1\right)\xi^{n+1}}\int_{0}^{2\pi}e^{im\varphi}d\varphi\int_{0}^{\xi}\left(1-2\zeta\cos\varphi+\zeta^{2}\right)^{1/2}d\zeta^{n+1}$
		$\displaystyle=$	$\displaystyle\frac{1}{n+1}\left[v_{m}-\xi^{2}u_{n+2}^{m}+\frac{1}{2}\xi\left(u_{n+1}^{m-1}+u_{n+1}^{m+1}\right)\right],\quad v_{m}=\int_{0}^{2\pi}e^{im\varphi}\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}d\varphi.$

Equation (53) can be checked directly using integral representation of $v_{m}$ and definition of $w_{m}\left(\xi\right).$

To obtain recurrence (54) we note the following transformation

	$\displaystyle i\int_{0}^{2\pi}\frac{e^{i\left(m-1\right)\varphi}\sin\varphi d\varphi}{\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}}=\frac{i}{\xi}\int_{0}^{2\pi}e^{i\left(m-1\right)\varphi}d\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}$		(83)
	$\displaystyle=-\frac{i}{\xi}\int_{0}^{2\pi}\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}de^{i\left(m-1\right)\varphi}=\left(m-1\right)\frac{1}{\xi}\int_{0}^{2\pi}\frac{e^{i\left(m-1\right)\varphi}\left(1-2\xi\cos\varphi+\xi^{2}\right)}{\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}}d\varphi.$

Using this expression, we obtain

$\displaystyle w_{m}$	$\displaystyle=$	$\displaystyle\int_{0}^{2\pi}\frac{e^{i\left(m-1\right)\varphi}e^{i\varphi}d\varphi}{\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}}=\int_{0}^{2\pi}\frac{e^{i\left(m-1\right)\varphi}\left(\cos\varphi+i\sin\varphi\right)d\varphi}{\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}}$
	$\displaystyle=$	$\displaystyle\left(m-1\right)\frac{1+\xi^{2}}{\xi}\int_{0}^{2\pi}\frac{e^{i\left(m-1\right)\varphi}d\varphi}{\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}}-\frac{2m-3}{2}\int_{0}^{2\pi}\frac{\left(e^{i\varphi}+e^{-i\varphi}\right)e^{i\left(m-1\right)\varphi}}{\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}}d\varphi$
	$\displaystyle=$	$\displaystyle\left(m-1\right)\frac{1+\xi^{2}}{\xi}w_{m-1}-\frac{2m-3}{2}\left(w_{m}+w_{m-2}\right).$

It is not difficult to see that Eq. (54) follows from this expression.

Functions $w_{0}$ and $w_{1}$ can be expressed immediately via complete elliptic integrals. Indeed, we have

	$\displaystyle w_{m}\left(\xi\right)$	$\displaystyle=$	$\displaystyle\int_{\pi}^{3\pi}\frac{e^{im\varphi}d\varphi}{\left(1-2\xi\cos\varphi+\xi^{2}\right)^{1/2}}=\left(-1\right)^{m}\int_{0}^{2\pi}\frac{e^{im\varphi}d\varphi}{\left(1+2\xi\cos\varphi+\xi^{2}\right)^{1/2}}$
		$\displaystyle=$	$\displaystyle 2\left(-1\right)^{m}\int_{0}^{\pi}\frac{\cos m\varphi d\varphi}{\left(1+2\xi\cos\varphi+\xi^{2}\right)^{1/2}}=\frac{4\left(-1\right)^{m}}{1+\xi}\int_{0}^{\pi/2}\frac{\cos\left(2m\psi\right)d\psi}{\left(1-\mu\sin^{2}\psi\right)^{1/2}},\quad\mu=\frac{4\xi}{\left(1+\xi\right)^{2}}.$

So,

w_{0}\left(\xi\right)=\frac{4}{1+\xi}K\left(\mu\right),\quad w_{1}\left(\xi\right)=\frac{4}{1+\xi}\left[\left(\frac{2}{\mu}-1\right)K\left(\mu\right)-\frac{2}{\mu}E\left(\mu\right)\right],

(86)

where $K$ and $E$ are defined by Eq. (56). The Landen transformation of elliptic functions [21] shows that

	$\displaystyle K\left(\mu\right)$	$\displaystyle=$	$\displaystyle\frac{2}{1+\mu_{1}^{1/2}}K\left(\mu_{2}\right),\quad\mu_{1}=1-\mu,\quad\mu_{2}=\left(\frac{1-\mu_{1}^{1/2}}{1+\mu_{1}^{1/2}}\right)^{2}$		(87)
	$\displaystyle\quad E\left(\mu\right)$	$\displaystyle=$	$\displaystyle\left(1+\mu_{1}^{1/2}\right)E\left(\mu_{2}\right)-\frac{2\mu_{1}^{1/2}}{1+\mu_{1}^{1/2}}K\left(\mu_{2}\right).$

Substituting this into Eq. (86) and using the value of $\mu$ from Eq. (A) we obtain Eq. (55).

Note now that relation (52) is also valid in case $m=0$ . In this case we can use symmetry $u_{n}^{-1}=u_{n}^{1}$ . Since Eq. (53) also holds in this case, we obtain using symmetry $w_{-1}=w_{1}$ and Eq. (55),

v_{0}=\left(1+\xi^{2}\right)w_{0}-2\xi w_{1}=8E\left(\xi^{2}\right)-4\left(1-\xi^{2}\right)K\left(\xi^{2}\right),

(88)

and confirms recurrence (51).

So, the last relation to prove is Eq. (50). This relation can be obtained using differentiation relations for the complete elliptic integral,

K\left(\mu\right)=\frac{d}{d\mu}\left[\mu\left(1-\mu\right)\frac{dK\left(\mu\right)}{d\mu}\right],\quad\frac{dK\left(\mu\right)}{d\mu}=\frac{E\left(\mu\right)}{2\mu\left(1-\mu\right)}-\frac{K\left(\mu\right)}{2\mu},

(89)

and integration by parts,

$\displaystyle u_{n}^{0}$	$\displaystyle=$	$\displaystyle\frac{4}{\xi^{n+1}}\int_{0}^{\xi}\zeta^{n}K\left(\zeta^{2}\right)d\zeta=\frac{8}{\xi^{n+1}}\int_{0}^{\xi^{2}}\mu^{\left(n-1\right)/2}\frac{d}{d\mu}\left(\mu\left(1-\mu\right)\frac{dK\left(\mu\right)}{d\mu}\right)d\mu$
	$\displaystyle=$	$\displaystyle\frac{8}{\xi^{n+1}}\left[\xi^{n+1}\left(1-\xi^{2}\right)\frac{dK\left(\xi^{2}\right)}{d\left(\xi^{2}\right)}-\frac{n-1}{2}\int_{0}^{\xi^{2}}\mu^{\left(n-1\right)/2}\left(1-\mu\right)dK\left(\mu\right)\right]$
	$\displaystyle=$	$\displaystyle\frac{4}{\xi^{2}}\left[E\left(\xi^{2}\right)-n\left(1-\xi^{2}\right)K\left(\xi^{2}\right)+\frac{\left(n-1\right)^{2}}{4}u_{n-2}^{0}\right]-\left(n^{2}-1\right)u_{n}^{0}.$

Equation (50) is just a rearrangement of terms in this relation.

	$\displaystyle\left.K^{(D)}\left(\mathbf{y},\mathbf{x}\right)\right\|_{\mathbf{y\in}S_{g\infty}^{+}}$	$\displaystyle=$	$\displaystyle-\left.G\left(\mathbf{y},\mathbf{x}\right)\right\|_{\mathbf{y\in}S_{g\infty}},$		(12)
	$\displaystyle\left.\frac{\partial K^{(N)}\left(\mathbf{y},\mathbf{x}\right)}{\partial n\left(\mathbf{y}\right)}\right\|_{\mathbf{y\in}S_{g\infty}^{+}}$	$\displaystyle=$	$\displaystyle-\left.\frac{\partial G\left(\mathbf{y},\mathbf{x}\right)}{\partial n\left(\mathbf{y}\right)}\right\|_{\mathbf{y\in}S_{g\infty}}.$

	$\displaystyle\left.\frac{\partial}{\partial z^{\prime}}\frac{1}{\left\|\mathbf{x}^{\prime}-\mathbf{y}\right\|}\right\|_{z^{\prime}=0}$	$\displaystyle=$	$\displaystyle\left.\frac{z_{y}}{\left\|\mathbf{x}^{\prime}-\mathbf{y}\right\|^{3}}\right\|_{z^{\prime}=0},$		(19)
	$\displaystyle\left.\left\|\mathbf{x}-\mathbf{x}^{\prime}\right\|\right\|_{z^{\prime}=0}$	$\displaystyle=$	$\displaystyle\left.\left(\left\|\mathbf{x}^{\prime}\right\|^{2}-2\mathbf{x}^{\prime}\cdot\mathbf{x}+\left\|\mathbf{x}\right\|^{2}\right)^{1/2}\right\|_{z^{\prime}=0}=\left(\rho^{\prime 2}-2\rho_{x}\rho^{\prime}\cos\left(\varphi^{\prime}-\varphi_{x}\right)+r_{x}^{2}\right)^{1/2}.$