Unified Construction of Green’s functions for Poisson’s equation in inhomogeneous media with diffuse interfaces

Journal of Computational and Applied Mathematics 326 (2017) 296–319

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Unified Construction of Green’s functions for Poisson’s equation in inhomogeneous media with diffuse interfaces Changfeng Xue a , Shaozhong Deng b, * a b

School of Mathematics and Physics, Yancheng Institute of Technology, Yancheng, Jiangsu 224051, PR China Department of Mathematics and Statistics, University of North Carolina at Charlotte, Charlotte, NC 28223, United States

article

info

Article history: Received 27 October 2016 Received in revised form 14 April 2017 Keywords: Green’s function Laplace’s equation Poisson’s equation Diffuse interface

a b s t r a c t Green’s functions for Poisson’s equation in inhomogeneous media with material interfaces have many practical applications. In the present work, we focus on Green’s functions for Poisson’s equation in inhomogeneous media with diffuse material interfaces where a gradual and continuous transition in the material constant is assumed in a small region around the interfaces between different materials. We present a unified general framework for calculating Green’s functions for Poisson’s equation in such inhomogeneous media and the framework can apply to all eleven orthogonal coordinate systems in which the three-dimensional Laplace equation is separable. Within this framework, the idea on how to design the so-called quasi-harmonic diffuse interface is discussed, formulations for building Green’s function for Poisson’s equation in an inhomogeneous medium with such a diffuse interface is elaborated, and a robust numerical method for calculating Green’s functions for Poisson’s equation in inhomogeneous media with general diffuse interfaces is developed. Several practically relevant separable coordinate systems are briefly surveyed at the level of definition and basic facts relevant for implementing the unified framework in these coordinate systems. © 2017 Elsevier B.V. All rights reserved.

1. Introduction In the present paper, we consider the three-dimensional Poisson equation

∇ · c(x)∇ u(x, xs ) = −4πδ (x − xs )

(1)

in the whole space R3 , where δ (·) is the Dirac delta function representing a unit singular source at xs ∈ R3 , and the nonconstant ‘‘material’’ coefficient c(x) > 0 corresponds to an inhomogeneous medium (i.e., different ‘‘materials’’ at different points in space). The solution u(x, xs ) defines Green’s function, often interchangeably denoted as G(x, xs ), for Poisson’s equation in the inhomogeneous medium. In practical applications, such a problem could be regarded e.g. as a problem of calculating static potential field at an arbitrary location due to a point charge in a dielectric material sharing a common boundary, also called an interface, with another dissimilar dielectric material. This problem occurs frequently in many physical, chemical, and biological applications; two representative areas of application include the study of solvation effects on biological macromolecules and the simulation of semiconductor quantum dots/wells; see [1–10] and references therein. The most common case of inhomogeneous media is where c(x) is piecewise constant. e.g. in one region we have one liquid, in another region we have a biological membrane, and in another region we have different liquid, each corresponding to a

*

Corresponding author. Fax: +1 704 687 1392. E-mail address: [email protected] (S. Deng).

http://dx.doi.org/10.1016/j.cam.2017.06.007 0377-0427/© 2017 Elsevier B.V. All rights reserved.

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

297

different material constant [11–13]. While it is relatively simple, this model of the inhomogeneous media is disappointing for at least two reasons. First, the assumption of the sharp jump in the material constant across the interfaces does not sound very physical on a microscopic scale in many practical applications. Second, it leads to several artifacts such as the wellknown divergence of the so-called self-polarization potential in the proximity of the interfaces. Many solutions have been proposed in the context of various applications to overcome these limitations of the sharp interface model the overwhelming majority of which, however, can be classified into one group: the diffuse interface model [14–22]. While they could differ greatly in their complexity and technical details, all diffuse interface models share a common design framework: a gradual and continuous transition in the material constant is assumed in a small region around the interfaces between different materials. In the present paper, we shall focus on Poisson’s equation in inhomogeneous media with diffuse interfaces. More specifically, we shall present a unified framework, a framework that can apply across various geometries, for calculating Green’s functions for Poisson’s equation in inhomogeneous media first with a specific so-called quasi-harmonic diffuse interface and then with a general diffuse interface. The paper is organized as follows. Section 2 provides a short introduction to a general three-dimensional orthogonal coordinate system and to the corresponding harmonic functions at the level of notation. Then Section 3 briefly presents Green’s function for Poisson’s equation in a piecewise-homogeneous medium for completeness of the paper. Next, in Section 4, the quasi-harmonic diffuse interface is introduced and Green’s function for Poisson’s equation in an inhomogeneous medium with such a diffuse interface is constructed. After that, Section 5 presents a robust numerical method for calculating Green’s functions for Poisson’s equation in inhomogeneous media with general diffuse interfaces. Then, several important separable coordinate systems are briefly surveyed in Section 6 at the level of definition and basic facts relevant for implementing the unified framework in these coordinate systems. Finally, a simple case study is presented in Section 7 and a few concluding remarks are given in Section 8. 2. Preliminaries Let (ξ , λ, ψ ) specify a general three-dimensional orthogonal coordinate system. We assume that in this coordinate system, an interface or a boundary surface is one of the coordinate surfaces, say ξ = X , a constant, and the coordinate variable ξ is thus referred to as the radial variable throughout this paper. We restrict ourselves to cases that the interface or the boundary surface is not too ‘‘pathological’’ in shape so that at least one such orthogonal coordinate system can be found. We also assume that Laplace’s equation 1u = 0 is completely separable in the coordinates (ξ , λ, ψ ) although the results presented in this unified framework may also apply to cases that Laplace’s equation is separable only in the radial variable ξ (namely, Laplace’s equation assumes solutions in the form of a product of a function of ξ alone and a function of λ and ψ but not ξ ). In other words, we assume that in this coordinate system, Laplace’s equation 1u = 0 assumes a family of completely separated solutions in the form of a product of three functions, each dependent on only one coordinate variable, namely, Ξ (ξ )Λ(λ)Ψ (ψ ), and that all solutions of the Laplace equation in question can be built up out of linear combinations of the members of this family of separated solutions, say, u=

∑

Ξα (ξ )Λα (λ)Ψα (ψ ),

(2)

α∈J

where J represents a finite or infinite index set. Moreover, we assume that all radial functions Ξα (ξ ) in (2) can be categorized into two groups. The radial functions in the first group are analogous to the radial functions r n , n = 0, 1, . . . , in the spherical harmonic theory in the sense that they numerically stay finite at ‘‘zero’’ but do not go to value 0 at ‘‘infinity’’, while those in the second group are analogous to the radial functions 1/r n+1 , n = 0, 1, . . . , in the spherical harmonic theory in the sense that they numerically go to value 0 at ‘‘infinity’’ but cannot stay finite at ‘‘zero’’. It should be clarified that ‘‘zero’’ in the radial coordinate variable ξ is the smallest value it can take in a particular coordinate system; for examples, ‘‘zero’’ is numerical value 0 in the spherical coordinate system, numerical value 1 in the prolate spheroidal coordinate system, and −∞ in the rectangular coordinate system, respectively, while ‘‘infinity’’ is the largest value the radial coordinate variable ξ can take in the coordinate system, and usually ‘‘infinity’’ is indeed +∞. We denote, respectively, by Eα (ξ ) the radial functions in the first group and by Fα (ξ ) the radial functions in the second group, and we refer them to, respectively, as the radial functions of the first and the second kind throughout the paper. We shall mention that this is indeed true in a great many cases including those to be surveyed in Section 6. On the other hand, we assume that the products Λα (λ)Ψα (ψ ) in (2) are analogous to the spherical harmonics Ynm (θ, φ ) defined on the surface of a sphere so they are referred to as the surface harmonics in the coordinate system (ξ , λ, ψ ) and we let Yα (λ, ψ ) = Λα (λ)Ψα (ψ ). The collection {Yα (λ, ψ ) : α ∈ J } is assumed to be a complete set of orthogonal functions on a surface specified by ξ = X , a constant, and in general the orthogonality relation is with respect to some geometric weighting function defined on the surface. Accordingly, the interior and the exterior harmonics are defined, respectively, as Eα (x) = Eα (ξ )Λα (λ)Ψα (ψ ),

(3a)

298

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

Fig. 1. Illustration of the sharp interface model: the space is divided into an interior region and an exterior region by one material interface Γ : ξ = ξb . The material coefficient in the interior region (ξ ≤ ξb ) is ci , and that in the exterior region (ξ > ξb ) is ce , respectively.

Fα (x) = Fα (ξ )Λα (λ)Ψα (ψ ).

(3b)

Then a solution of the Laplace equation can be expanded in terms of these harmonics. In general, in a region that contains ‘‘zero’’, the solution shall be expanded in terms of only the interior harmonics, and in a region that contains ‘‘infinity’’, the solution shall be expanded in terms of only the exterior harmonics. In other regions, however, the solution shall be expanded in terms of both the interior and the exterior harmonics. Since it is assumed that Laplace’s equation 1u = 0 is completely separable in the coordinates (ξ , λ, ψ ), we write 1 = 1ξ + 1λ + 1ψ where 1ξ , 1λ , and 1ψ stand for the second-order differential operator with respect only to ξ , only to λ, and only to ψ , respectively. We call a function f of ξ that satisfies 1ξ f = 0 a radial harmonic function. Note that f (ξ ) ≡ 1 is always a radial harmonic function. Finally, we assume that in the coordinates (ξ , λ, ψ ), the reciprocal distance is expanded as

⎧∑ ⎪ Kα Eα (xs )Fα (x), ⎪ ⎨ 1 α = ∑ |x − xs | ⎪ ⎪ Kα Fα (xs )Eα (x), ⎩

if ξ ≥ ξs , if ξ ≤ ξs ,

(4)

α

where ξs is the ξ -coordinate of the source and Kα are constant expansion coefficients which in general depend on the normalization constant of the orthogonal surface harmonics Yα (λ, ψ ). To conclude this section, we introduce the following shorthand notations in order to make later formulations concise. We let µα (ξ ) and να (ξ ) be the ratios of the radial functions of the first and the second kind, namely,

µα (ξ ) =

Eα (ξ ) F α (ξ )

and να (ξ ) =

Fα (ξ ) E α (ξ )

,

(5)

and E˜α (ξ ) and F˜α (ξ ) be the logarithmic derivatives of the radial functions, namely, E˜α (ξ ) =

Eα′ (ξ ) E α (ξ )

and

F˜α (ξ ) =

Fα′ (ξ ) F α (ξ )

.

(6)

3. Piecewise-homogeneous media As mentioned earlier, the most common case of inhomogeneous media is where c(x) is piecewise constant. For simplicity, here we consider a piecewise-homogeneous medium with only one material interface Γ specified by the coordinate surface ξ = ξb . Suppose that the space R3 is divided into an interior region and an exterior region by this interface; see Fig. 1, and we assume that c(x) = ci in the interior region and c(x) = ce in the exterior region, respectively, namely,

{ c(x) =

ci , ce ,

if ξ ≤ ξb , if ξ > ξb .

(7)

On the interface Γ : ξ = ξb , both the solution u (scalar potential field) and the normal flux c(x)un are continuous. In addition, the solution u is assumed to be zero at ‘‘infinity’’. So we have the following interface–boundary value problem for Poisson’s equation:

∇ · c(x)∇ u(x, xs ) = −4πδ (x − xs ), u(ξb+ ) = u(ξb− ), ce uξ (ξb+ ) = ci uξ (ξb− ), ( ) 1 u(x, xs ) = O , |x| → ∞. |x|

(8a) (8b) (8c) (8d)

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

299

Green’s function for Poisson’s equation in such a piecewise-homogeneous medium can be built straightforwardly. Let us begin by considering the case when the source is located at point xs = (ξs , λs , ψs ) in the interior region. In this case, the solution must be finite at ‘‘infinity’’ so the solution in the exterior region (ξb ≤ ξ ), denoted as ue , shall be expanded in terms of only the exterior harmonics Fα (x). On the other hand, Green’s function in the interior region (ξ ≤ ξb ), denoted as ui , should be composed of two parts as follows: ui (x, xs ) = GS (x, xs ) + R(x, xs ),

(9)

where the singular part GS (x, xs ) =

1

(10)

ci |x − xs |

carries the behavior at the singular source xs and the regular part R(x, xs ) is a harmonic function that secures the satisfaction of the interface conditions so it shall be expanded in terms of only the interior harmonics Eα (x). In short, the solution ue or ui at point x = (ξ , λ, ψ ) in the exterior or interior region due to a source at point xs in the interior region (so ξb ≥ ξs ) takes the form 1 ∑ (1) Aα Kα Eα (xs )Fα (x), (11a) ue (x, xs ) = √ ci ce α

ui (x, xs ) =

1 ci |x − xs |

+

1 ∑ ci

√

να (ξb )B(1) α Kα Eα (xs )Eα (x).

(11b)

α

Here, factors such as 1/ci , 1/ ci ce , and Kα Eα (xs ) are extracted from expansion coefficients of the harmonics explicitly (1) so that A(1) α and Bα , the remaining parts of the expansion coefficients, depend neither on the source nor on the particular values of ci and ce . Rather, as is given below, they depend only on the interface ξ = ξb and the ratio cr = ci /ce and thus only need to be calculated once even if there are many sources at different locations present. (1) The unknown expansion coefficients A(1) α and Bα can be determined by the interface conditions on the interface ξ = ξb , which, in terms of ue and ui , become ue |ξ =ξ + = ui |ξ =ξ − ,

(12a)

b ⏐ b ⏐ ∂ ue ⏐⏐ ∂ ui ⏐⏐ = c , ce i ∂ξ ⏐ξ =ξ + ∂ξ ⏐ξ =ξ − b

(12b)

b

the supposed expansion of the reciprocal distance (4), and the orthogonality relation of the surface harmonics. Omitting all details, we have

[ √

cr ,

F˜α (ξb ),

−1 √ − cr E˜α (ξb )

](

A(1) α

)

B(1) α

(

1

)

. = √ ˜ cr Fα (ξb )

(13)

Similarly, when the source is located at point xs = (ξs , λs , ψs ) in the exterior region (so ξs ≥ ξb ), the solutions in the two regions take the form ue (x, xs ) =

1 ce |x − xs |

ui (x, xs ) = √

1 ci ce

∑

+

1 ∑ ce

µα (ξb )A(2) α Kα Fα (xs )Fα (x),

(14a)

α

B(2) α Kα Fα (xs )Eα (x),

(14b)

α

(2) where the expansion coefficients A(2) α and Bα are given by

[ √

cr ,

F˜α (ξb ),

−1 − cr E˜α (ξb ) √

](

A(2) α

)

B(2) α

( =

√ ) − cr . −E˜α (ξb )

(15)

4. Quasi-harmonic diffuse interface Now we turn to Green’s function for Poisson’s equation in an inhomogeneous medium with a diffuse interface, a small ‘‘thin’’ region around the boundary between two materials in which a gradual and continuous transition in the material constant c(x) is assumed. Furthermore, we assume the coefficient c(x) depends only on the radial coordinate ξ so c(x) = c(ξ ). More specifically, we assume the space R3 is divided into three regions, referred to as the interior region Ωi , the exterior region Ωe , and a small, thin middle region Ωt in between, respectively. The interface between Ωi and Ωt , denoted by Γi , is specified by the coordinate surface ξ = ξi , and that between Ωt and Ωe , denoted by Γe , is specified by the coordinate surface ξ = ξe , respectively. The middle region Ωt is also referred to as the transition region or diffuse interface throughout the paper. The coefficient c(x) is assumed to be constant in both the interior and the exterior regions, say again c(x) = ci in Ωi

300

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

Fig. 2. Illustration of a diffuse interface model: The material coefficient in the interior region (ξ ≤ ξi ) is ci , that in the exterior region (ξ ≥ ξe ) is ce , while the middle transition region (ξi < ξ < ξe ) assumes a continuous coefficient ct (ξ ) that connects ci and ce at the two interfaces Γi and Γe .

and c(x) ≡ ce in Ωe , respectively, while, as a function of only the radial coordinate variable ξ , it changes continuously (and monotonically) from ci to ce through the transition region; see Fig. 2. Two choices of the coefficient profile in the diffuse interface include the linear profile defined by ci + ce

ct (ξ ) =

2

+

ci − ce

(

ξe + ξi − 2ξ ξe − ξi

)

,

(16)

) ξ − ξi cos π . ξe − ξi

(17)

2

and the cosine-like profile given by ci + ce

ct (ξ ) =

2

+

ci − ce 2

(

On the interfaces Γi and Γe , both the solution u and the normal flux c(x)un are continuous. In addition, the solution u is assumed to be zero at infinity. So we have the following interface–boundary value problem for Poisson’s equation:

∇ · c(x)∇ u(x, xs ) = −4πδ (x − xs ), u(ξi+ ) = u(ξi− ), u(ξe+ ) = u(ξe− ), uξ (ξi+ ) = uξ (ξi− ), uξ (ξe+ ) = uξ (ξe− ), ( ) 1 u(x, xs ) = O , |x| → ∞. |x|

(18a) (18b) (18c) (18d) (18e) (18f)

For convenience, in the sequel we denote by ui , ut , and ue the solution in the interior, the transition, and the exterior region, respectively. Generally speaking, the above interface–boundary value problem cannot be solved analytically. For some special diffuse interface models, however, the series solutions of the interface–boundary value problem (18) in terms of the interior and exterior harmonics could still be built up. In particular, let Υ (ξ ) be one of the non-constant radial harmonic functions for the Laplace operator (Note that 1 is a radial harmonic function.). Then we introduce the following diffuse interface model in which the radially dependent coefficient c(ξ ) takes the form: c(ξ ) =

⎧ ⎨ci ,

c (ξ ) = [α + β Υ (ξ )]2 , ⎩ct , e

if ξ ≤ ξi , if ξi < ξ < ξe , if ξ ≥ ξe ,

(19)

where α and β are constants chosen to ensure the continuity of c(ξ ) at ξi and ξe , namely,

√

α=

s ce − t

√ β=

√

ci

s−t ci −

√

s−t

ce

,

(20a)

,

(20b)

with s = Υ (ξi ),

t = Υ (ξe ).

√

Note that ct (ξ ) defined in (19) is not harmonic by itself, but ct (ξ ), being a linear combination of two radial harmonic eigenfunctions 1 and Υ (ξ ), is harmonic. As such, this diffuse interface is termed the quasi-harmonic diffuse interface. The three diffuse interface models mentioned so far together with the sharp interface model are illustrated in Fig. 3. The interface–boundary value problem (18) with the quasi-harmonic diffuse interface can be solved analytically because now the Poisson equation (18a) in each of the three regions can be transformed into a Poisson equation involving only the

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

301

Fig. 3. Illustration of diffuse interface models together with the sharp interface model, assuming ci > ce .

Laplace operator 1, the key for us to be able to do so being the following theorem [14] whose two-dimensional version can be found in [23]. Theorem 1. If the variable coefficient c(x) satisfies

√ 1 c(x) = 0,

(21)

then the solution u(x) of the quasi-harmonic equation

∇ · c(x)∇ u(x) = 0

(22)

satisfies

1

]

[√

c(x)u(x) = 0,

(23)

and similarly, the solution v (x) of the quasi-elliptic equation

∇ · c(x)∇v (x) = ρ (x)

(24)

satisfies

1

[√

]

ρ (x)

c(x)v (x) = √

c(x)

.

(25)

Proof. Since (22) and (23) can be regarded as a special case of (24) and (25), we only need to prove the latter which can be achieved at mainly by repeatedly using the vector identity:

∇ · (χ A) = A · ∇χ + χ∇ · A,

(26)

where χ and A are a scalar function and a vector function, respectively. First, we have

(√ ) √ √ √ 1 cu = u1 c + 2∇ c · ∇ u + c 1u. √ Since 1 c = 0, we obtain (√ ) √ √ 1 cu = 2∇ c · ∇ u + c 1u.

(27)

(28)

Noting that

√ 1 ∇ c = √ ∇c, 2 c

(29)

we further have

1

(√ )

√ 1 cu = √ ∇ c · ∇ u + c 1u, c

(30)

302

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

from which we get

√

c1

(√ )

cu = ∇ c · ∇ u + c 1u.

(31)

Now, using (24) we obtain

∇ c · ∇ u + c 1u = ∇ u · ∇ c + c ∇ · ∇ u = ∇ · (c ∇ u) = ρ. Finally plugging (32) into (31) and then dividing both sides by

1

[√

(32)

√

c yield

ρ (x) . c(x)u(x) = √

]

c(x)

By the way how ct (ξ ) is constructed in the quasi-harmonic diffuse interface model, we certainly have diffuse interface. Therefore, the Poisson equation (18a) in the transition region can now be rewritten as

1

[√

√

ct (ξ ) = 0 in the

] 4π c(ξ )ut (x, xs ) = − √ δ (x − xs ), c(ξ )

(33)

√

whose solution c(ξ )ut (x, xs ) shall be expanded in terms of both the interior and the exterior harmonics. Consequently, ut (x, xs ) can be expressed as 1 ∑ [Cα Eα (x) + Dα Fα (x)] , ut (x, xs ) = √ c(ξ ) α

(34)

when the source is not located in the transition region, or as 1 1 ∑ [Cα Eα (x) + Dα Fα (x)] , ut (x, xs ) = √ +√ cs c(ξ )|x − xs | c(ξ ) α

(35)

when the source is located in the transition region, where cs = c(ξs ), Cα and Dα are undetermined constant expansion coefficients. 4.1. Case I: The source located inside the interior region When the source is located inside the interior region, namely, when ξs ≤ ξi , by invoking Theorem 1, the Poisson equation (18a) in the three regions becomes 4π 1ui (x, xs ) = − δ (x − xs ), ci ] [√ c(ξ )ut (x, xs ) = 0, 1

1ue (x, xs ) = 0,

if ξ ≤ ξi ,

(36a)

if ξi < ξ < ξe ,

(36b)

if ξ ≥ ξe .

(36c)

At the two boundaries of the diffuse interface, the continuity of the solution and the normal flux require that ui |ξ =ξ − = ut |ξ =ξ + ,

(37a)

⏐ i ⏐i ∂ ui ⏐⏐ ∂ ut ⏐⏐ = , ∂ξ ⏐ξ =ξ − ∂ξ ⏐ξ =ξ + i

(37b)

i

ue |ξ =ξ + = ut |ξ =ξ − , ⏐ e ⏐e ∂ ue ⏐⏐ ∂ ut ⏐⏐ = . ∂ξ ⏐ξ =ξe+ ∂ξ ⏐ξ =ξe−

(37c) (37d)

The solutions in the three regions can be written as (after some constants being factored out from expansion coefficients for the purpose of simplifying mathematical formulations) ue (x, xs ) = √ ui (x, xs ) =

1

∑

ci ce

(38a)

α

1 ci |x − xs |

ut (x, xs ) = √

A(1) α Kα Eα (xs )Fα (x),

1 ci c(ξ )

+

∑ α

1 ∑ ci

B(1) α Kα Eα (xs )Eα (x),

(38b)

α

Kα Eα (xs ) Cα(1) Eα (x) + D(1) α Fα (x) .

[

]

(38c)

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

303

(1) (1) (1) Here, the constant expansion coefficients A(1) α , Bα , Cα , and Dα can be determined by the interface conditions (37), together with the orthogonality relation of the surface harmonics Yα (λ, ψ ), and the expansion of the reciprocal distance (4). Omitting all tedious details, we get

⎛

A(1) α

⎞

⎞ ⎛ −Fα (ξi ) ⎜ B(1) ⎟ ⎜ ⎟ ⎜ 0 ⎟ , M × ⎜ α(1) ⎟ = ⎝ ′ −Fα (ξi )⎠ ⎝Cα ⎠

(39)

0

D(1) α

where the matrix M is 0, ⎢Fα (ξe ),

Eα (ξi ), 0,

⎡

⎢ ⎢ ⎢ M = ⎢ 0, ⎢ ⎢ ⎣ ′ Fα (ξe ),

Eα′ (ξi ),

−Eα (ξi ), −Eα (ξe ), β Υ ′ (ξ i ) −Eα′ (ξi ) + √ Eα (ξi ), ci

−Eα′ (ξe ) +

0,

β Υ ′ (ξe ) Eα (ξe ), √ ce

⎤ −Fα (ξi ) ⎥ −Fα (ξe ) ⎥ ′ ⎥ β Υ (ξi ) −Fα′ (ξi ) + √ F α (ξ i ) ⎥ ⎥. ci ⎥ ⎥ ′ ⎦ β Υ ( ξ ) e −Fα′ (ξe ) + √ Fα (ξe )

(40)

ce

4.2. Case II: the source located inside the exterior region When the source is located inside the exterior region, namely, when ξs ≥ ξe , by invoking Theorem 1, the Poisson equation (18a) in the three regions becomes

1

1ui (x, xs ) = 0, ] c(ξ )ut (x, xs ) = 0,

[√

1ue (x, xs ) = −

4π ce

δ (x − xs ),

if ξ ≤ ξi ,

(41a)

if ξi < ξ < ξe ,

(41b)

if ξ ≥ ξe .

(41c)

Now the solutions in the three regions can be written as ue (x, xs ) =

1

+

ce |x − xs |

ui (x, xs ) = √ ut (x, xs ) = √

1 ci ce

∑

1 ∑ ce

A(2) α Kα Fα (xs )Fα (x),

(42a)

α

B(2) α Kα Fα (xs )Eα (x),

(42b)

α

1

∑

ce c(ξ )

Kα Fα (xs ) Cα(2) Eα (x) + D(2) α Fα (x) ,

[

]

(42c)

α

(2) (2) (2) where the expansion coefficients A(2) α , Bα , Cα , and Dα can be calculated by

⎛

A(2) α

⎞

⎛ ⎞ 0 ⎜ (2) ⎟ ⎜ Bα ⎟ ⎜−Eα (ξe )⎟ ⎟ M×⎜ ⎜C (2) ⎟ = ⎝ 0 ⎠ . ⎝ α ⎠ −Eα′ (ξe ) (2)

(43)

Dα

4.3. Case III: the source located inside the diffuse interface When the source is located inside the diffuse interface, namely, when ξi < ξs < ξe , by invoking Theorem 1, the Poisson equation (18a) in the three regions becomes

1

1ui (x, xs ) = 0, ] 4π c(ξ )ut (x, xs ) = − √ δ (x − xs ), c(ξ ) 1ue (x, xs ) = 0,

[√

if ξ ≤ ξi ,

(44a)

if ξi < ξ < ξe ,

(44b)

if ξ ≥ ξe .

(44c)

Now we can write the solutions in the three regions as ue (x, xs ) = √

1 cs ce

∑ α

(4) Kα A(3) α Fα (xs ) + Aα Eα (xs ) Fα (x),

[

]

(45a)

304

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

ui (x, xs ) = √

1 cs ci

∑

(4) Kα B(3) α Fα (xs ) + Bα Eα (xs ) Eα (x),

]

[

α

1

ut (x, xs ) = √

+√

1

∑

{[

Kα

c c(ξ )|x − xs | cs c(ξ ) α ] } [ s (4) E (x F (x ) + D + D(3) α s ) Fα (x) , α s α α (3)

(3)

(45b)

Cα(3) Fα (xs ) + Cα(4) Eα (xs ) Eα (x)

]

(45c)

(3)

(3)

where the expansion coefficients Aα , Bα , Cα , and Dα are calculated by

⎛

A(3) α

⎞

⎛ ⎞ Eα (ξi ) ⎜ (3) ⎟ ⎜ 0 ⎟ ⎜ Bα ⎟ ⎜ ⎟ ⎟=⎜ ′ β Υ ′ (ξi ) M×⎜ ⎟, ⎜C (3) ⎟ ⎝Eα (ξi ) − √ E α (ξ i ) ⎠ ⎝ α ⎠ ci D(3) α

(46)

0

(4) (4) (4) while A(4) α , Bα , Cα , and Dα are calculated by

⎛

A(4) α

⎞

⎛ ⎜ (4) ⎟ ⎜ ⎜ Bα ⎟ ⎜ ⎟ M×⎜ ⎜C (4) ⎟ = ⎜ ⎝ α ⎠ ⎝ D(4) α

Fα′ (ξe )

−

0 Fα (ξe ) 0 β Υ ′ (ξe )

√

ce

⎞

Eα (ξe )

⎟ ⎟ ⎟. ⎠

(47)

Remark 2. In the above formulations, all constant expansion coefficients depend only on the two boundaries ξ = ξi and ξ = ξe of the diffuse interface as well as the two coefficient constants ci and ce . Therefore, once the diffuse interface is chosen, these coefficients have to be calculated only once even when there are many sources present. 4.4. Improving the formulations Typically, for good numerical accuracy, practical implementation of the series solution requires evaluation of the interior and the exterior harmonic functions of high degrees/orders, which poses an additional numerical challenge. For example, the numerical magnitude of the radial functions of high degrees/orders may easily exceed the capacity of a computer especially when ξ is small for Fα (ξ ) or when ξ is large for Eα (ξ ). Therefore, to avoid potential overflow and to reduce computer cutoff errors, we further carry out more convenient rewriting of the previous formulations so that the matrices of the resulting systems for the expansion coefficients involve only the ratios of the radial functions Eα (ξ ) and Fα (ξ ) and the logarithmic derivatives of these functions. To this end, in addition to the shorthand notations given in (5) and (6), we let

γα = µα (ξi )να (ξe ) =

E α (ξ i ) Eα (ξe )

×

Fα (ξe ) Fα (ξi )

.

(48)

Then, when the source is located inside the interior region, the solutions are ue (x, xs ) = √

1 ci ce 1

∑

A(1) α Kα Eα (xs )Fα (x),

(49a)

α

1 ∑ + να (ξi )B(1) α Kα Eα (xs )Eα (x), ci |x − xs | ci α ∑ [ ] 1 ut (x, xs ) = √ Kα Eα (xs ) να (ξe )Cα(1) Eα (x) + D(1) α Fα (x) , ci c(ξ ) α ui (x, xs ) =

(49b) (49c)

(1) (1) (1) where the expansion coefficients A(1) α , Bα , Cα , and Dα are now calculated by

A(1) α

⎛

⎞

⎛ ⎞ −1 ⎜ (1) ⎟ ⎜ ⎟ B ⎜ 0 ⎟ α ⎟ M(1) × ⎜ ⎜C (1) ⎟ = ⎝−F˜α (ξi )⎠ , ⎝ α ⎠

(50)

0

D(1) α

with the matrix M(1) defined as

⎡ M(1)

0, 1,

⎢ ⎢ ⎢ 0, =⎢ ⎢ ⎢ ⎣ F˜α (ξe ),

1, 0, E˜α (ξi ), 0,

−γα , −1, ] [ β Υ ′ (ξi ) ˜ γα , −Eα (ξi ) + √ ci

−E˜α (ξe ) +

β Υ ′ (ξe ) , √ ce

⎤ −1 ⎥ −1 ⎥ β Υ ′ (ξ i ) ⎥ ⎥. −F˜α (ξi ) + √ ci ⎥ ⎥ β Υ ′ (ξ e ) ⎦ −F˜α (ξe ) + √ ce

(51)

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

305

Similarly, when the source is located inside the exterior region, the solutions are 1

ue (x, xs ) =

ce |x − xs |

ui (x, xs ) = √ ut (x, xs ) = √

1 ce ci

1 ∑

+

∑

ce

µα (ξe )A(2) α Kα Fα (xs )Fα (x),

(52a)

α

B(2) α Kα Fα (xs )Eα (x),

(52b)

α

1

∑

ce c(ξ )

Kα Fα (xs ) Cα(2) Eα (x) + µα (ξi )D(2) α Fα (x) ,

]

[

(52c)

α

(2) (2) (2) where the expansion coefficients A(2) α , Bα , Cα , and Dα are now calculated by

A(2) α

⎛

⎞

⎛ ⎞ 0 ⎜ B(2) ⎟ ⎜ ⎟ ⎜ −1 ⎟ M(2) × ⎜ α ⎟ = ⎝ ⎠, 0 ⎝Cα(2) ⎠ −E˜α (ξe ) (2)

(53)

Dα

with the matrix M(2) defined as

⎡ M(2)

0, 1,

⎢ ⎢ ⎢ 0, =⎢ ⎢ ⎢ ⎣ F˜α (ξe ),

1, 0,

−1, −1, β Υ ′ (ξi ) −E˜α (ξi ) + √ ,

E˜α (ξi ),

ci β Υ (ξe ) ′

0,

−E˜α (ξe ) +

√

ce

,

⎤ −1 ⎥ −γα ⎥ β Υ ′ (ξi ) ⎥ ⎥. −F˜α (ξi ) + √ ⎥ ci ⎥ ] [ ⎦ β Υ ′ (ξe ) ˜ γα −Fα (ξe ) + √

(54)

ce

Finally, when the source is located inside the diffuse interface, the solutions are ue (x, xs ) = √ ui (x, xs ) = √

1 cs ce 1 cs ci

ut (x, xs ) = √

∑

(4) Kα µα (ξi )A(3) α Fα (xs ) + Aα Eα (xs ) Fα (x),

(55a)

(4) Kα B(3) α Fα (xs ) + να (ξe )Bα Eα (xs ) Eα (x),

(55b)

[

]

α

∑

[

]

α

1

+√

1

∑

Kα

{[

] γα Cα(3) Fα (xs ) + να (ξe )Cα(4) Eα (xs ) Eα (x)

cs c(ξ )|x − xs | cs c(ξ ) α [ ] } (3) + µα (ξi )Dα Fα (xs ) + γα D(4) α Eα (xs ) Fα (x) ,

(55c)

(3) (3) (3) where the expansion coefficients A(3) α , Bα , Cα , and Dα are calculated now by

⎛

A(3) α

⎞

⎛

1 0

⎞

⎜ (3) ⎟ ⎜ ⎟ ⎜ Bα ⎟ ⎜ ⎟=⎜ β Υ ′ (ξi ) ⎟ ⎟, ˜ (3) ⎟ ⎝Cα ⎠ ⎝Eα (ξi ) − √ci ⎠

M(1) × ⎜ ⎜

D(3) α

(56)

0

(4) (4) (4) while A(4) α , Bα , Cα , and Dα are calculated by

⎛

A(4) α

⎞

⎛ ⎜ (4) ⎟ ⎜ ⎜ Bα ⎟ ⎜ ⎟ M(2) × ⎜ ⎜C (4) ⎟ = ⎜ ⎝ α ⎠ ⎝˜ D(4) α

0 1 0

⎞

⎟ ⎟ ⎟. ′ β Υ (ξe ) ⎠ Fα (ξe ) − √

(57)

ce

5. Numerical method for general diffuse interfaces In an inhomogeneous medium with a general diffuse interface such as the linear diffuse interface defined by (16) or the cosine-like diffuse interface given by (17), the interface–boundary value problem (18) usually does not admit an easy series solution, if there exists one such solution at all, and thus in general needs to be solved numerically. In doing so, a conventional approach first divides the transition region where a variable coefficient c(x) is assumed into multiple subregions and then, in each one of these subregions, approximates the variable coefficient c(x) by a constant value such as its mean value in this subregion [11]. As the result, the original Poisson equation in an inhomogeneous medium of continuous material constant

306

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

Fig. 4. Approximation of the linear (left) or cosine-like (right) diffuse interface model, assuming ci > ce . Solid line: the original linear or cosine-like model, and dot-dashed line: the approximation.

reduces back to a Poisson equation in a piecewise-homogeneous medium. Although the solution can then be found by using a procedure in analogy to the analysis of transmission lines in the same way as for solving Poisson’s equation in layered spheres [11,24–28], physically this treatment will introduce numerical divergence in the self-polarization potential in the proximity of every sharp interface between subregions, a problem that we are trying to avoid for the first place by utilizing a diffuse interface. Fortunately, empowered by the quasi-harmonic diffuse interface model, we can develop a computationally much more robust (but only slightly more complicated in idea) numerical procedure for solving the interface–boundary value problem (18) with a general diffuse interface. The first step of this numerical procedure is still to divide the transition region specified by ξ ∈ [ξi , ξe ] into multiple, say L − 1, subregions, ξ ∈ [ξl−1 , ξl ], l = 1, 2, . . . , L − 1, with ξ0 = ξi and ξL−1 = ξe . For convenience, also set ξ−1 = ‘‘zero’’ and ξL = ‘‘infinity’’. For each index l = −1, 0, . . . , L, denote by cl the value of the material constant evaluated from c(ξ ) at ξ = ξl , namely, cl = c(ξl ). Note that c−1 = c0 = ci and cL−1 = cL = ce . Then in each subregion specified by ξ ∈ [ξl−1 , ξl ], the coefficient c(ξ ) is approximated by a quasi-harmonic one of the form (19) that connects cl−1 and cl , namely, by cl (ξ ) = [αl + βl Υ (ξ )]2 , where

αl =

sl

√

cl − tl

√

ξl−1 ≤ ξ ≤ ξl , √

cl−1

and βl =

sl − t l

(58)

cl−1 −

√

cl

sl − t l

,

(59)

with sl = Υ (ξl−1 )

and

tl = Υ (ξl ) .

(60)

Fig. 4 shows the approximation of the radially dependent material coefficient c(ξ ) in the linear or cosine-like diffuse interface model when the transition region specified by ξ ∈ [ξi , ξe ] is divided into three subregions. Next, in each subregion with ξ ∈ [ξl−1 , ξl ], the Poisson equation (18a) is approximated as

∇ · cl (ξ )∇ u(x, xs ) = −4πδ (x − xs ), √ which, by noting that cl (ξ ) = 0, and using Theorem 1, can be rewritten as 4π

1 [cl (ξ )u(x, xs )] = − √

cl (ξ )

(61)

δ (x − xs ).

(62)

Consequently, the solution in this subregion, denoted by u(l) , is expressed either as u(l) (x, xs ) = √

1 cl (ξ )

∑

Kα Cα(l) Eα (x) + D(l) α Fα (x)

[

]

(63)

α

if the source is not located inside this subregion, or as u(l) (x, xs ) = √

1 cl (ξs )cl (ξ )|x − xs |

+√

1 cl (ξ )

∑ α

Kα Cα(l) Eα (x) + D(l) α Fα (x)

[

]

(64)

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

307

if otherwise. The coefficients Cα(l) and D(l) α for l = 0, 1, . . . , L are to be determined by the interface conditions for l = 0, 1, . . . , L − 1: u(l−1) |ξ =ξ − = u(l) |ξ =ξ + , l −1

⏐ ∂ u(l−1) ⏐⏐ ∂ξ ⏐

ξ =ξl− −1

(65a)

⏐ l −1 ∂ u(l) ⏐⏐ = . ∂ξ ⏐ξ =ξ +

(65b)

l−1

In fact, the constants Cα(l) and D(l) α can be obtained by a recursive procedure in analogy to the analysis of transmission lines [5,11,25–28]. Here, we pass the straightforward but tedious details involved in matching the interface conditions and reformulating the solution and give only the final formulations for the numerical procedure. In what follows, for each α ∈ J we define

γα,l = µα (ξl−1 )να (ξl ),

l = 0, 1, . . . , L,

and more generally, we let

γα,ij = µα (ξi−1 )να (ξj ),

0 ≤ i ≤ j ≤ L.

In addition, in the following formulations, Rα,l and Tα,l , identified as the interface parameters associated to the interface

ξ = ξl−1 , are given by

(βl − βl−1 ) Υ ′ (ξl−1 ) , 1 ] [ α,l √ cl−1 F˜α (ξl−1 ) − E˜α (ξl−1 ) = , 1α,l [ ] √ = cl−1 F˜α (ξl−1 ) − E˜α (ξl−1 ) − (βl − βl−1 ) Υ ′ (ξl−1 ).

Rα,l =

(66a)

Tα,l

(66b)

1α,l

(l)

(66c)

(l)

On the other hand, Rn and Tn , identified as the interface reflection and the transmission coefficients associated to the interface ξ = ξl−1 , respectively, are given through the following recursive expressions R(lα−1) = Rα,l +

Tα(l−1) = Rα,l +

2 (l−2) γα,l−1 Tα, l Rα (l−2)

1 − Rα,l Rα

γα,l−1

2 (l) Tα, l Tα γα,l

with R(α−1) = 0,

with Tα(L) = 0.

(l)

1 − Rα,l Tα γα,l

(67)

(68)

5.1. Case I: the source located inside Subregion 0: ξs ∈ [‘‘zero’’, ξi ] When the source is located inside Subregion 0, i.e., the interior region of the quasi-harmonic diffuse interface, the solutions in all subregions are 1 ∑ (L) u(L) (x, xs ) = √ Dα Kα Eα (xs )Fα (x), c0 cL

(69a)

∑ [ ] 1 u(l) (x, xs ) = √ Kα Eα (xs ) να (ξl )Cα(l) Eα (x) + Dα(l) Fα (x) , c0 cl (ξ ) α

(69b)

α

u(0) (x, xs ) =

1 c0 |x − xs |

+

1 ∑ c0

να (ξ0 )Cα(0) Kα Eα (xs )Fα (x),

(69c)

α

where l in (69b) takes values between 1 and L − 1, and the coefficients Cα(l) and D(l) α , l = 0, 1, . . . , L, are calculated recursively by D(l) α =

Tα,l

Cα(l) = Tα(l) D(l) α , (0)

(l)

1 − Rα,l Tα γα,l

in which Dα = 1.

D(lα−1) ,

(70a) (70b)

308

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

5.2. Case II: the source located inside Subregion L: ξs ∈ [ξe , ‘‘infinity’’) When the source is located inside Subregion L, i.e., the exterior region of the quasi-harmonic diffuse interface, the solutions in all subregions are u(L) (x, xs ) =

1 cL |x − xs |

+

1 ∑ cL

µα (ξL−1 )D(L) α Kα Fα (xs )Fα (x),

(71a)

α

∑ [ ] 1 u(l) (x, xs ) = √ Kα Fα (xs ) Cα(l) Eα (x) + µα (ξl−1 )Dα(l) Fα (x) , cL cl (ξ ) α

(71b)

1 ∑ (0) Cα Kα Fα (xs )Eα (x), u(0) (x, xs ) = √ cL c0

(71c)

α

where l in (71b) takes values between 1 and L − 1, and the coefficients Cα(l) and Dα(l) , l = 0, 1, . . . , L, are calculated recursively by Cα(l−1) =

Tα,l (l−2)

1 − Rα,l Rα

Cα(l) ,

γα,l−1

(72a)

D(lα−1) = R(lα−2) Cα(l−1) ,

(72b)

in which Cα(L) = 1. 5.3. Case III: the source located inside Subregion m: ξs ∈ [ξm−1 , ξm ], where 0 < m < L When the source is located inside Subregion m with 0 < m < L, the solution in Subregion m is u(m) (x, xs ) = √

1 cs cm (ξ )|x − xs |

+√

(m,1)

+ µα (ξm−1 )Dα

[

1

cs cm (ξ )

∑

Kα

{[ ] γα,m Cα(m,1) Fα (xs ) + να (ξm )Cα(m,2) Eα (xs ) Eα (x)

α

,2) Fα (xs ) + γα,m D(m Eα (xs ) Fα (x) , α

]

}

(73)

in which the ξs -independent expansion coefficients are defined as Cα(m,2) = ,1) D(m = α

Tα(m) (m−1)

1 − Rα

(m)

Tα γα,m

−1) R(m α (m−1)

1 − Rα

(m)

Tα γα,m

,

(74a)

,

(74b)

−1) (m,2) Cα(m,1) = R(m Cα , α

(74c)

,2) ,1) D(m = Tα(m) D(m . α α

(74d)

It should be pointed out that now cs should be interpreted as cs = cl (ξs ). For l = m − 1, . . . , 0, the solution u(l) in Subregion l is

∑ {[ ] 1 u(l) (x, xs ) = √ Kα Cα(l,1) Fα (xs ) + να (ξm )Cα(l,2) Eα (xs ) Eα (x) cs cl (ξ ) α [ ] } + µα (ξl−1 )Dα(l,1) Fα (xs ) + γα,lm Dα(l,2) Eα (xs ) Fα (x) ,

(75)

where Cα(l,1) =

Cˆ α(l) (m−1)

1 − Rα

D(lα,1) = R(lα−1) Cα(l,1) ,

(m)

Tα γα,m

,

Cα(l,2) = Tα(m) Cα(l,1) ,

D(lα,2) = Tα(m) D(lα,1) .

Here, Cˆ α(l) , l = m − 1, . . . , 0, are calculated by the recursive formula (72a) with Cˆ α(m) = 1, namely, m−1

Cˆ α(l) =

∏

Tα,j+1 (j−1)

j=l

1 − Rα,j+1 Rα

γα,j

.

(76)

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

309

On the other hand, for l = m + 1, . . . , L, the solution u(l) is

∑ {[ ] 1 u(l) (x, xs ) = √ Kα γα,ml Cα(l,1) Fα (xs ) + να (ξl )Cα(l,2) Eα (xs ) Eα (x) cs cl (ξ ) α [ ] } + µα (ξm−1 )D(lα,1) Fα (xs ) + Dα(l,2) Eα (xs ) Fα (x) ,

(77)

where

ˆ (l) D α

D(lα,2) =

(m−1)

1 − Rα

Cα(l,2) = Tα(l) D(lα,2) ,

(m)

Tα γα,m

,

−1) (l,2) Dα , D(lα,1) = R(m α −1) (l,2) Cα . Cα(l,1) = R(m α

ˆ (l) ˆ (m) Here, D α , l = m + 1, . . . , L, are calculated by the recursive formula (70a) with Dα = 1, namely, ˆ (l) D α =

l ∏

Tα,j

j=m+1

1 − Rα,j Tα γα,j

(j)

.

(78)

5.4. A unified formulation The numerical solution can be written in a unified form. In short, the solution of the interface–boundary value problem (18) at point x = (ξ , λ, ψ ) in the lth subregion due to a source at point xs = (ξs , λs , ψs ) in the mth subregion, denoted by u(lm) , can be calculated by u(lm) (x, xs ) = √

1 cl (ξ )cm (ξs )

∑

Kα Yα (λs , ψs )Yα (λ, ψ ) ×

α

Wα,lm (ξ , ξs ) 1 − pα,mm qα,mm

,

(79)

where the functions Wα,lm (ξ , ξs ) are given by Wα,lm (ξ , ξs ) = pα,mm Eα (ξ> ) + Fα (ξ> )

[

][

Eα (ξ< ) + qα,lm Fα (ξ< ) Cˆ α(l)

]

(80)

if l < m, Wα,mm (ξ , ξs ) = pα,mm Eα (ξ> ) + Fα (ξ> )

[

][

Eα (ξ< ) + qα,mm Fα (ξ< )

]

(81)

for l = m, and Wα,lm (ξ , ξs ) = pα,lm Eα (ξ> ) + Fα (ξ> )

[

][

ˆ (l) Eα (ξ< ) + qα,mm Fα (ξ< ) D α

]

(82)

when l > m. Here, ξ< = min{ξ , ξs }, ξ> = max{ξ , ξs }, and pα,lm = Tα(m) να (ξl ),

(83a)

−1) qα,lm = R(m µα (ξl−1 ). α

(83b)

6. Survey of separable coordinates In this section, we survey several three-dimensional orthogonal coordinate systems in which Laplace’s equation is separable, at the level of definition and basic facts relevant for implementing the unified framework for calculating Green’s function for Poisson’s equation in an inhomogeneous medium with a diffuse interface in these coordinate systems. These coordinate systems are selected mainly because they have more practical relevance than those separable coordinate systems not surveyed here. 6.1. Spherical coordinates Consider the spherical coordinates (r , θ, φ ) defined through x = r sin θ cos φ,

(84a)

x = r sin θ sin φ,

(84b)

z = r cos θ ,

(84c)

310

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

where r ∈ [0, ∞) is the radial variable, θ ∈ [0, π] is the polar angle variable, and φ ∈ [0, 2π] is the azimuthal variable, respectively. The surface of constant r = X is a sphere of radius X . The Laplace equation in the spherical coordinates is

[

1 ∂

(

r2 ∂r

∂ ∂r

r2

) +

∂ 2 r sin θ ∂θ 1

(

∂ ∂θ

sin θ

) +

] ∂2 Ψ = 0. r 2 sin2 θ ∂φ 2 1

(85)

n m m −(n+1) m The interior (regular) and the exterior (irregular) solid harmonics are Rm Yn (θ, φ ), n (x) = r Yn (θ, φ ) and In (x) = r respectively, for n = 0, 1, . . . and m = 0, 1, . . . , n, where Ynm (θ, φ ) are the (surface) spherical harmonics. The expansion of the reciprocal distance between x = (r , θ, φ ) and xs = (rs , θs , φs ) is [29]

1

=

|x − xs |

∞ ∑ n ∑

n r<

(

κm

)

n+1 r>

n=0 m=0

Pnm (cos θs )Pnm (cos θ ) cos[m(φ − φs )],

(86)

where r< = min{r , rs }, r> = max{r , rs }, and

κm = (2 − δm0 )

(n − m)! (n + m)!

,

(87)

and Pnm (·) is the associated Legendre function of the first kind of degree n and order m. In particular, Pn (x) ≡ Pn0 (x) is the Legendre function of the first kind. Here, the Kronecker delta δij = 1 when i = j and is zero otherwise. In solving (1), without loss of generality, we can assume the source is on the positive polar axis so θs = φs = 0. In this case, due to the azimuthal symmetry, the solution can be expanded using only the interior harmonics r n Pn (cos θ ) and the exterior harmonics r −(n+1) Pn (cos θ ). Accordingly, the expansion of the reciprocal distance (86) becomes 1

=

|x − xs |

∞ ( n ) ∑ r< n+1 r>

n=0

Pn (cos θ ).

(88)

Then the unified framework can be applied with the following substitutions. (ξ , λ, ψ ) : (r , θ , φ ) with r ∈ [0, ∞) α ∈ J : {n : n = 0, 1, . . .} Eα (ξ ) : r n with E0 (ξ ) = 1 Fα (ξ ) : 1/r n+1 with F0 (ξ ) = 1/r Yα (λ, ψ ) : Pn (cos θ ) Kα : Hn = 1 Υ (ξ ) : F0 (ξ ) = 1/r. In particular, the radially dependent coefficient c(x) for the diffuse interface becomes [17]

c(x) =

⎧ ci , ⎪ ⎪ ⎪ ⎨

(

ct (r) =

α+

⎪ ⎪ ⎪ ⎩ce , where

√

s ce − t

α=

√

s−t

ci

β

)2

r

if r ≤ ri ,

,

(89)

if r ≥ re ,

√ ,

if ri < r < re ,

β=

ci −

√

s−t

ce

,

with s=

1 ri

,

t=

1 re

.

6.2. Prolate spheroidal coordinates Given a reference prolate spheroid x2 + y2

+

a2

z2 b2

= 1,

(90)

where b > a > 0. Here, the z-axis is the focal symmetry axis, and the interfocal distance is 2f with f = that we adopt the following definition of the prolate spheroidal coordinates (ξ , η, φ ) defined through x=f

√

y=f

√

(ξ 2 − 1)(1 − η2 ) cos φ, (ξ − 1)(1 − η

z = f ξ η,

2

2 ) sin

φ,

√

b2 − a2 . Suppose (91a) (91b) (91c)

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

311

where ξ ∈ [1, ∞) is the radial variable, η ∈ [−1, 1] is the angular variable, and φ ∈ [0, 2π] is the azimuthal variable, respectively. Note √ that the surface of constant ξ is a prolate spheroid with interfocal distance 2f ; in particular, ξ = ξb with ξb = b/f = b/ b2 − a2 is the prolate spheroid given by (90), and ξ = 1 corresponds to the line between the two focal points, respectively. The Laplace equation in the prolate spheroidal coordinates is

] ∂ 2 ∂ ∂ ξ 2 − η2 ∂2 2 ∂ ( ξ − 1) + (1 − η ) + Ψ = 0. f 2 (ξ 2 − η2 ) ∂ξ ∂ξ ∂η ∂η (ξ 2 − 1)(1 − η2 ) ∂φ 2 [

1

(92)

The interior and the exterior harmonics are Pnm (ξ )Pnm (η)eimφ and Qnm (ξ )Pnm (η)eimφ respectively, for n = 0, 1, . . . and m = 0, 1, . . . , n. Here, Qnm (·) is the associated Legendre function of the second kind of degree n and order m. In particular, Qn (x) ≡ Qn0 (x) is the Legendre function of the second kind. The expansion of the reciprocal distance between x = (ξ , η, φ ) and xs = (ξs , ηs , φs ) is [29] 1

=

|x − xs |

∞ ∑ n ∑

Hmn Pnm (ξ< )Qnm (ξ> )Pnm (η)Pnm (ηs ) cos[m(φ − φs )],

(93)

n=0 m=0

where ξ< = min{ξ , ξs }, ξ> = max{ξ , ξs }, and 1

Hmn =

f

(2n + 1)(2 − δm0 )(−1)m

[

(n − m)!

]2

(n + m)!

.

(94)

In solving (1), without loss of generality, we can assume the source is on the xz-plane so φs = 0. In this case, the solution can be expanded using only the even (with respect to φ ) interior harmonics Pnm (ξ )Pnm (η) cos(mφ ) and the even exterior harmonics Qnm (ξ )Pnm (η) cos(mφ ). Accordingly, the expansion of the reciprocal distance (93) becomes [29–31] 1

=

|x − xs |

∞ ∑ n ∑

Hmn Pnm (ξ< )Qnm (ξ> )Pnm (η)Pnm (ηs ) cos(mφ ).

(95)

n=0 m=0

Then the unified framework can be applied with the following substitutions. (ξ , λ, ψ ) : (ξ , η, φ ) with ξ ∈ [1, ∞) α ∈ J : {(n, m) : n = 0, 1, . . . ; m = 0, 1, . . . , n} Eα (ξ ) : Pnm (ξ ) with E(0,0) (ξ ) = P0 (ξ ) = 1 ( ) Fα (ξ ) :

Qnm (ξ ) with F(0,0) (ξ ) = Q0 (ξ ) =

1 2

ln

ξ +1 ξ −1

Yα (λ, ψ ) : η) cos mφ Kα : Hmn Υ (ξ ) : F(0,0) (ξ ). In particular, the radially dependent coefficient c(x) for the diffuse interface becomes [7] Pnm (

⎧ ci , ⎪ ⎪ ⎨

[ ( )]2 β ξ +1 c(x) = ct (ξ ) = α + ln , ⎪ 2 ξ −1 ⎪ ⎩ ce , where

√

s ce − t

α=

s−t

√

ci

√ ,

ci −

β=

if ξ ≤ ξi , if ξi < ξ < ξe ,

(96)

if ξ ≥ ξe ,

√

ce

s−t

,

with s=

1 2

( ln

) ξi + 1 , ξi − 1

t=

1 2

( ln

) ξe + 1 . ξe − 1

6.3. Oblate spheroidal coordinates Given a reference oblate spheroid x2 + y2 a2

+

z2 b2

= 1,

where 0 < b < a. Here, the z-axis is the focal symmetry axis, and the interfocal distance is 2f with f = that we adopt the following definition of the oblate spheroidal coordinates (ξ , η, φ ) defined through x=f

√

(1 + ξ 2 )(1 − η2 ) cos φ,

(97)

√

a2 − b2 . Suppose (98a)

312

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

√

(1 + ξ 2 )(1 − η2 ) sin φ,

y=f

(98b)

z = f ξ η,

(98c)

where ξ ∈ [0, ∞) is the radial variable, η ∈ [−1, 1] is the angular variable, and φ ∈ [0, 2π ] is the azimuthal variable, respectively. Note √ that the surface of constant ξ is an oblate spheroid with interfocal distance 2f ; in particular, ξ = ξb with ξb = b/f = b/ a2 − b2 is the oblate spheroid given by (97), and ξ = 0 corresponds to the circular disk in the plane z = 0 with radius f , respectively. The Laplace equation in the oblate spheroidal coordinates is

] ∂ 2 ∂ ∂ ξ 2 + η2 ∂2 2 ∂ (ξ + 1) + (1 − η ) + Ψ = 0. f 2 (ξ 2 + η2 ) ∂ξ ∂ξ ∂η ∂η (ξ 2 + 1)(1 − η2 ) ∂φ 2 [

1

(99)

The results for the prolate spheroidal case can be extended to the oblate spheroidal coordinates basically by following √ this rule: replace ξ by iξ where i = −1 and Hmn by Hmn =

i f

(2n + 1)(2 − δm0 )(−1)m

[

(n − m)!

]2

(n + m)!

.

(100)

As an example, the radially dependent coefficient c(x) for the diffuse interface becomes [7]

⎧ ci , ⎪ ⎪ ⎨

if ξ ≤ ξi ,

[ ( )]2 β iξ + 1 c(x) = ct (ξ ) = α + ln , ⎪ 2 iξ − 1 ⎪ ⎩ ce , where

√

s ce − t

α=

√

ci

s−t

√ ,

β=

ci −

if ξi < ξ < ξe ,

(101)

if ξ ≥ ξe ,

√

s−t

ce

,

with s=

1 2

( ln

iξi + 1

)

iξi − 1

,

t=

1 2

( ln

iξ e + 1 iξe − 1

)

.

6.4. Ellipsoidal coordinates Suppose that we adopt Hobson’s formalism for the ellipsoidal coordinates [30]. Then the ellipsoidal coordinates (ξ , µ, ν ) corresponding to point (x, y, z) in the rectangular coordinates, generated by a reference ellipsoid x2

y2

z2

= 1, c2 where a > b > c > 0 are its semi-axes, satisfy a2

x2

+

b2

+

y2

(102)

z2

= 1, (103) λ λ − λ − k2 where λ stands for either √ ξ , µ, or ν . The two √ constants k and h are determined by the semi-focal distances of the reference ellipsoid, namely, k = a2 − c 2 and h = a2 − b2 . The variable ξ in the ellipsoidal coordinate system is called radial and √ 2 2 assumes √ only positive values, namely, ξ ∈ [k, ∞). The surface of ξ = constant is a tri-axial ellipsoid of semi-axes ξ , ξ − h and ξ 2 − k2 ; in particular, ξ = a corresponds to the reference ellipsoid (102). The transformation between the ellipsoidal 2

+

2

+

h2

2

and the rectangular coordinates is [30] x2 = y2 = z2 =

ξ 2 µ2 ν 2 k2 h2

,

(ξ 2 − h2 )(µ2 − h2 )(h2 − ν 2 ) h2 (k2 − h2 ) (ξ 2 − k2 )(k2 − µ2 )(k2 − ν 2 ) k2 (k2 − h2 )

(104a)

,

(104b)

.

(104c)

The Laplace equation is separable in the ellipsoidal coordinates [32]. As a matter of fact, each of the three ellipsoidal variables satisfies the same Lamé differential equation

( 2 )( ) d2 E(λ) ( ) dE(λ) λ − h2 λ2 − k2 + λ 2λ2 − h2 − k2 2 dλ dλ [( ) ] + h2 + k2 q − n(n + 1)λ2 E(λ) = 0, where q is an arbitrary constant to be fixed appropriately [29].

(105)

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

313

p

The solutions of (105), En , are called Lamé functions of the first kind of degree n and order p. Here both indices n and p are positive integers satisfying 2n + 1 ≥ p ≥ 1. For the determination of the Lamé functions, see e.g., [29], for a short list of the Lamé functions of the first kind, see e.g., [33], and for numerical computation of the Lamé functions, see e.g., [34–36]. A general interior ellipsoidal harmonic which is a normal solution of the Laplace equation and is regular at the origin may be written in terms of the Lamé product as Enp (x) = Enp (ξ )Enp (µ)Enp (ν ).

(106)

Similarly, an exterior ellipsoidal harmonic which is regular at infinity is defined as Fnp (x) = Fnp (ξ )Enp (µ)Enp (ν ),

(107)

p Fn (

where ξ ) is the Lamé function of the second kind of degree n and order p, which is related to the corresponding Lamé p function of the first kind En (ξ ) by [30] Fnp (ξ ) = (2n + 1)Enp (ξ )

∞

∫ ξ

dξ ′ √( )( ). [Enp (ξ ′ )]2 ξ ′2 − h2 ξ ′2 − k2

p

(108)

p

Both the interior En (x) and the exterior Fn (x) ellipsoidal harmonics are linearly independent and form a complete set of functions. Moreover, there exists an orthogonality relation of the Lamé functions of the first kind as [33,34] h

∫ 0

p

p′

p

p′

En (µ)En (ν )En′ (µ)En′ (ν ) µ2 − ν 2

k

∫

(

)

√( )( )( )( ) dµdν = γnp δnn′ δpp′ , µ2 − h2 k2 − µ2 h2 − ν 2 k2 − ν 2

h

(109)

p

where the normalization constant γn is h

∫

γnp =

0

p

p

En (µ)En (ν )

]2 (

) µ2 − ν 2 √( )( )( )( ) dµdν. µ2 − h2 k2 − µ2 h2 − ν 2 k2 − ν 2 [

k

∫ h

(110)

The expansion for the reciprocal distance in the ellipsoidal coordinates is [33] 1

|x − xs |

=

∞ 2n +1 ∑ ∑

Knp Enp (ξ< )Fnp (ξ> )Enp (µ)Enp (µs )Enp (ν )Enp (νs ),

(111)

n=0 p=1

where ξ< = min{ξ , ξs }, ξ> = max{ξ , ξs }, and Knp =

π 2(2n + 1)γnp

.

(112)

Then the unified framework can be applied with the following substitutions. (ξ , λ, ψ ) : (ξ , µ, ν ) with ξ ∈ [k, ∞) α ∈ J : {(n, p) : n = 0, 1, . . . ; p = 1, 2, . . . , 2n + 1} p Eα (ξ ) : En (ξ ) with E(0,1) (ξ ) = E01 (ξ ) = 1 p Fα (ξ ) : Fn (ξ ) with F(0,1) (ξ ) = F01 (ξ ) p p Yα (λ, ψ ) : En (µ)En (ν ) Kα : Knp Υ (ξ ) : F(0,1) (ξ ) = F01 (ξ ). In particular, the radially dependent coefficient c(x) for the diffuse interface becomes [4] c(x) = where

⎧ ⎨ci ,

[

√

α=

if ξ ≤ ξi ,

]2

c (ξ ) = α + β F01 (ξ ) , ⎩t ce ,

s ce − t

√

ci

s−t

√ and β =

if ξi < ξ < ξe , if ξ ≥ ξe , ci −

(113)

√

ce

s−t

,

(114)

with s = F01 (ξi )

and

t = F01 (ξe ) .

(115)

Note that in terms of the Jacobi elliptic function sn(·), we have [29] F01 (

ξ) =

dξ ′

∞

∫ ξ

√

(ξ ′2 − k2 )(ξ ′2 − h2 )

=

1 k

sn

−1

(

k h

, ξ k

)

.

(116)

314

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

So 1 ′ F01 (ξ ) = − √ . 2 2 (ξ − k )(ξ 2 − h2 )

(117)

6.5. Circular cylinder coordinates Consider the typical circular cylinder coordinates (ρ, φ, z) in which ρ is the radial coordinate, φ is the polar angle, and z is the height, respectively. Assume that c(x) varies along only the radial direction and that, without loss of any generality, the source is located at xs = (ρs , φs = 0, zs = 0). Then Green’s function G(x, xs ) = G(ρ, φ, z ; ρs ) can be expanded in a Fourier–Bessel form [13] as G(x, xs ) =

∞

∫

dk cos(kz) 0

) ∞ ( ∑ 2 − δn0 2π 2

n=0

ˆ n (k, ρ, ρs ). cos(nφ )G

(118)

ˆ n (k, ρ, ρs ) satisfies the following mono-dimensional Substituting (118) into (1), one can find that the radial Green’s function G differential equation 1 d

[

ρ dρ

ρ c(ρ )

]

d dρ

ˆ n (k, ρ, ρs ) − G

(

k2 +

n2

ρ2

)

ˆ n (k, ρ, ρs ) = − c(ρ )G

4π

ρ

δ (ρ − ρs )

(119)

ˆ n (k, ρ, ρs ) is finite as ρ → 0 and, on the other hand, Gˆ n (k, ρ, ρs ) → 0 as ρ → ∞. with the boundary condition that G ˆ n (k, ρ, ρs ). Likely, in the region where The unified framework can be applied to solve (119) for the radial Green’s function G the source is located, the radial Green’s function is composed of two parts as ˆ n (k, ρ, ρs ) = Gˆ Sn (k, ρ, ρs ) + Rn (k, ρ, ρs ), G

(120)

where 4π

ˆ Sn (k, ρ, ρs ) = G

cs

In (kρ< )Kn (kρ> ),

(121)

with ρ< = min{ρ, ρs } and ρ> = max{ρ, ρs }, carries the behavior at the singular source xs and Rn (k, ρ, ρs ) is a regular part that can be written in terms of In (kρ ) and Kn (kρ ), say, Rn (k, ρ, ρs ) = An In (kρ ) + Bn Kn (kρ ).

(122)

Here, In (ρ ) and Kn (ρ ) are the modified Bessel functions of the first and the second kind which are monotonically increasing ˆ Sn (k, ρ, ρs ) is actually the radial Green’s function in the and decreasing with respect to ρ > 0, respectively. Note that G homogeneous medium of constant coefficient cs = c(ρs ) [13], from which we can obtain the expansion of the reciprocal distance: ∞

∫

1

|x − xs |

dk cos(kz)

= 0

∞ ∑

Kn cos(nφ )In (kρ< )Kn (kρ> ),

(123)

n=0

where Kn =

4 − 2δn0

π

.

(124)

ˆ n (k, ρ, ρs ) For each fixed n = 0, 1, . . . , the unified framework can be applied to solve (119) for the radial Green’s function G with the following substitutions. Radial variable ξ : ρ with ρ ∈ [0, ∞) α ∈ J : {n} Eα (ξ ) : In (kρ ) Fα (ξ ) : Kn (kρ ) Yα (λ, ψ ) : 1 Kα : 4π Υ (ξ ) : ln(ρ ). Note that ln(ρ ) is a harmonic function with radial symmetry. The radially dependent coefficient c(x) for the diffuse interface becomes [14] c(x) =

⎧ ⎨ci ,

c (ρ ) = [α + β ln(ρ )]2 , ⎩ct , e

if ρ ≤ ρi , if ρi < ρ < ρe , if ρ ≥ ρe ,

(125)

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

where

√

α=

s ce − t

√

√

ci

and β =

s−t

ci −

315

√

ce

s−t

,

(126)

with and t = ln (ρe ) .

s = ln (ρi )

(127)

6.6. Rectangular coordinates Consider the rectangular coordinates (x, y, z). Assume that c(x) is a function of only one rectangular coordinate, say the z coordinate. By exploiting the translational symmetry along the x and y directions, one can put the origin √ of the x and y axes at the position of the source (so xs = (0, 0, zs )), introduce the two-dimensional radial coordinate ρ = x2 + y2 , and expand the Green’s function G(x, xs ) = G(ρ, z , zs ) in a Fourier–Bessel form [21]: ∞

∫

G(ρ, z , zs ) =

ˆ , z , zs )k dk, J0 (ρ k)G(k

(128)

0

where J0 is the Bessel function of zero order. Using the naming convention for Green’s functions, we call G(ρ, z , zs ) and ˆ , z , zs ) Green’s functions in the spatial and the spectral domains, or simply ‘‘the spatial Green’s function’’ and ‘‘the spectral G(k Green’s function’’, respectively. Substituting (128) into (1) and using the properties of the Bessel functions, one can find that ˆ , z , zs ) satisfies the following mono-dimensional differential equation the spectral Green’s function G(k d

[ c(z)

dρ

]

d dρ

ˆ , z , zs ) − k2 c(z)G(k ˆ , z , zs ) = −2δ (z − zs ) G(k

(129)

ˆ , z , zs ) → 0 as |z | → ∞. with the boundary condition that G(k ˆ , z , zs ). Likely, in the region The unified framework can be applied to solve (129) for the spectral Green’s function G(k where the source is located, the spectral Green’s function is composed of two parts as ˆ , z , zs ) = Gˆ S (k, z , zs ) + R(k, z , zs ), G(k

(130)

where

ˆ S (k, z , zs ) = G

1 kcs

e−k|z −zs |

(131)

carries the behavior at the singular source xs and R(k, z , zs ) is a regular part that can be written in terms of ekz and e−kz , say, R(k, z , zs ) = Aekz + Be−kz .

(132)

ˆ S (k, z , zs ) is actually the spectral Green’s function in the homogeneous medium of constant coefficient cs = c(zs ), Note that G from which we can obtain the expansion of the reciprocal distance: ∞

∫

1

|x − xs |

J0 (ρ k)e−k|z −zs | dk.

=

(133)

0

ˆ , z , zs ) For each fixed k ∈ (0, ∞), the unified framework can be applied to solve (129) for the spectral Green’s function G(k with the following substitutions. Radial variable ξ : z with z ∈ (−∞, ∞) α ∈ J : {1} Eα (ξ ) : ekz Fα (ξ ) : e−kz Yα (λ, ψ ) : 1 Kα : 1/k Υ (ξ ) : z. Note that z is a harmonic function. The radially dependent coefficient c(x) for the diffuse interface becomes [15] c(x) = where

α=

zi

⎧ ⎨ci ,

c (z) = (α + β z)2 , ⎩ct , e

√

ce − ze zi − ze

√

ci

if z ≤ zi , if zi < z < ze , if ze ≤ z ,

√ ,

β=

ci −

√

zi − ze

ce

.

(134)

316

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

7. A case study on self-polarization energy Semiconductor quantum dots have been a subject of intense theoretical and experimental interest in science and technology [37,38]. For simplicity, in many theoretical studies of quantum dots, macroscopic dielectric constants εi and εe are assigned for the dot and the surrounding matrix, respectively. It is well-known from classical electromagnetism that the presence of a charged particle inside such a quantum dot polarizes the surrounding matrix, inducing charges at the dot surface of the same sign as the source charge if εi > εe (the typical situation) [39]. Consequently, a repulsive (positive) potential energy contribution, commonly called the self-polarization energy, arises due to the mutual interaction between the source and its own induced charges. This quantity is important in deciding the electronic states, energy spectrum, or other characteristics of such semiconductor structures. Here, we present a simple case study of computing the self-polarization energy of charged particles in a spherical quantum dot. Mathematically, the potential energy between two charged particles e and es , located at x and xs respectively, is V (x, xs ) = ees u(x, xs ), where u(x, xs ) is the electrostatic potential that verifies the Poisson equation (1) where c(x) is now the dielectric permittivity function ε (x). The self-polarization energy Vs (x) of a charged particle e located at x can be calculated from V (x, xs ) by taking e = es and x = xs , excluding the direct Coulomb interaction from u(x, xs ), and dividing by 2 as it corresponds to a self-energy. Applying the proposed unified general framework to the spherical coordinates (r , θ, φ ) as described in Section 6.1 yields

⎧ 2 ∞ ( )2(n+1) ∑ e ⎪ (2) rb ⎪ A , ⎪ n ⎪ ⎨ 2εe rb r n=0 Vs (x) = ( )2n ∞ ⎪ e2 ∑ ⎪ r ⎪ ⎪ , B(1) ⎩ n 2εi rb rb

if r ≥ rb , (135) if r ≤ rb ,

n=0

(2)

(1)

where rb is the radius of the dot, and An and Bn are calculated through (13) and (15). Fig. 5 shows the self-polarization energy profile for the case that rb = 10 nm, εi = 12.6 (GaAs), and εe = 1 (vacuum). All illustrative plots below are based on the self-polarization energy of 201 unit charges (in atomic unit) uniformly distributed along the radial direction from r = 0 to r = 20 nm. As can be seen, the self-polarization energy clearly diverges at the dot surface, which is unphysical. This is because, under the sharp interface model, all induced charges are localized at the dot surface of zero width so both the real and the induced charges can coincide at the same location. In fact, as indicated earlier, at the dot surface the sharp transition from εi to εe in the dielectric constant is also unphysical because the experimental growth procedure of a quantum dot cannot guarantee such a sharp permittivity profile without inter-diffusion between the dot and the matrix materials [11]. Similarly, the self-polarization energy of a charged particle e located at x under the quasi-harmonic diffuse interface model (89) can be calculated by

⎧ ∞ ⎪ e2 ∑ (2) ( re )2(n+1) ⎪ ⎪ An , ⎪ ⎪ 2εe re r ⎪ ⎪ n = 0 ⎪ ⎪ ⎪ ( )2n ∞ ⎪ ⎪ e2 ∑ (1) r ⎪ ⎪ B , ⎪ n ⎪ ⎨ 2εi ri ri n=0 [ ∞ Vs (x) = ( )2n+1 ( )2n ∞ 2 ⎪ ⎪ 1 ∑ (3) ri 1 ∑ (4) r ⎪ e ⎪ Cn + Cn ⎪ ⎪ 2ε (r) r re re re ⎪ ⎪ n=0 n=0 ⎪ ⎪ ⎪ ( )2n+1 ] ∞ ∞ ⎪ ⎪ 1 ∑ (3) ( ri )2(n+1) 1 ∑ (4) ri ⎪ ⎪ Dn + Dn , ⎪ ⎩ +r r r r i

(2)

(4)

n=0

n=0

if r ≥ re , if r ≤ ri , (136)

if ri < r < re ,

e

where An to Dn are calculated through (50), (53), (56) and (57). Fig. 6 shows the self-polarization energy profile under the quasi-harmonic model with diffuse interfaces of different thickness. As can be seen, the self-polarization energy still exhibits certain unphysical singularity at both edges of the diffuse interface, precisely where the first derivative of ε (r) is discontinuous, but this singularity is integrable. In addition, as the diffuse interface decreases in size, the self-polarization energy under the quasi-harmonic model reduces to that under the sharp interface model. Fig. 7 shows the self-polarization energy profile under the linear and the cosine-like models with interface thickness of 10 nm, obtained by the numerical method proposed in Section 5. As can be seen, the self-polarization energy under the cosine-like model does not have unphysical singularity at either edge of the diffuse interface, which is attributed to the fact that the first derivative of ε (r) is continuous at both edges of the diffuse interface, suggesting that, in applications where the presence of unphysical singularity in the self-polarization energy may cause issues, diffuse interfaces with smooth radially dependent coefficient profiles should be used.

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

317

Fig. 5. Typical self-polarization energy profile under the sharp interface model.

Fig. 6. Typical self-polarization energy profile under the quasi-harmonic diffuse interface model.

Fig. 7. Typical self-polarization energy profile under the linear and the cosine-like diffuse interface models.

Finally, we plot together in Fig. 8 the self-polarization energy profiles under the cosine-like model with interface thickness of 10 nm, obtained respectively by the proposed numerical method and by the conventional approach without any additional augmentation: first divide the transition region where a variable coefficient c(x) is assumed into multiple subregions and then, in each one of these subregions, approximate the variable coefficient c(x) by a constant value such as its mean value

318

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

Fig. 8. Self-polarization energy profile under the cosine-like diffuse interface model, obtained by the proposed numerical method and the conventional approach, respectively.

in this subregion [11]. As expected, the self-polarization energy profile obtained by the conventional approach exhibits numerical divergence within the transition region, clearly indicating that the numerical method proposed in this paper is much more stable. 8. Conclusions Green’s functions for Poisson’s equation in inhomogeneous media with material interfaces have many practical applications. The most common cases of inhomogeneous media are piecewise-homogeneous media where sharp jumps in the material constant are assumed across the interfaces between materials. While these cases of inhomogeneous media are relatively simple, the assumption of the sharp jumps in the material constant across the interfaces is not physical on a microscopic scale in many practical applications. As such, in the present work, we focused on Green’s functions for Poisson’s equation in inhomogeneous media with diffuse interfaces where a gradual and continuous transition in the material constant is assumed in a small region around the interfaces between different materials instead. We presented a unified framework that can apply to all orthogonal coordinate systems in which the three-dimensional Laplace equation is separable, for calculating Green’s functions for Poisson’s equation in such inhomogeneous media. Within this framework, first the idea on how to design the so-called quasi-harmonic diffuse interface is discussed. Then, formulations for building Green’s function for Poisson’s equation in an inhomogeneous medium with the quasi-harmonic diffuse interface is elaborated. Next, a robust numerical method for calculating Green’s functions for Poisson’s equation in inhomogeneous media with general diffuse interfaces is developed. In addition, we surveyed several practically quite relevant separable coordinate systems to show how the unified framework can be applied to solve Poisson’s equation in an inhomogeneous medium with a diffuse interface in these coordinate systems. Acknowledgment This work was supported by the National Natural Science Foundation of China [Grant Number 11471281]. References [1] R.D. Remigio, K. Mozgawa, H. Cao, V. Weijo, L. Frediani, A polarizable continuum model for molecules at spherical diffuse interfaces, J. Chem. Phys. 144 (2016) 124103. [2] F. Fahrenberger, Z. Xu, C. Holm, Simulation of electric double layers around charged colloids in aqueous solution of variable permittivity, J. Chem. Phys. 141 (2014) 064902. [3] C. Xue, S. Deng, Numerical calculation of electronic states in ellipsoidal finite-potential quantum dots with an off-centered impurity, Physica E 43 (2011) 1118–1126. [4] C. Xue, S. Deng, Three-layer dielectric models for generalized Coulomb potential calculation in ellipsoidal geometry, Phys. Rev. E 83 (2011) 056709. [5] S. Deng, A robust numerical method for self-polarization energy of spherical quantum dots with finite confinement barriers, Comput. Phys. Comm. 181 (2010) 787–799. [6] C. Xue, S. Deng, A robust numerical method for generalized Coulomb and self-polarization potentials of dielectric spheroids, Commun. Comput. Phys. 8 (2010) 374–402. [7] C. Xue, S. Deng, Three-dielectric-layer hybrid solvation model with spheroidal cavities in biomolecular simulations, Phys. Rev. E 81 (2010) 016701. [8] T.A.S. Pereira, J.S. de Sousa, J.A.K. Freire, G.A. Farias, On the interplay between quantum confinement and dielectric mismatch in high-k based quantum wells, J. Appl. Phys. 108 (2010) 054311. [9] C. Curutchet, R. Cammi, B. Mennucci, S. Corni, Self-consistent quantum mechanical model for the description of excitation energy transfers in molecules at interfaces, J. Chem. Phys. 125 (2006) 054710.

C. Xue, S. Deng / Journal of Computational and Applied Mathematics 326 (2017) 296–319

319

[10] L. Frediani, R. Cammi, S. Corni, J. Tomasi, A polarizable continuum model for molecules at diffuse interfaces, J. Chem. Phys. 120 (2004) 3893–3907. [11] P.G. Bolcatto, C.R. Proetto, Partially confined excitons in semiconductor nanocrystals with a finite size dielectric interface, J. Phys.: Condens. Matter. 13 (2001) 319–334. [12] P.G. Bolcatto, C.R. Proetto, Self-polarization energies of semiconductor quantum dots with finite confinement barriers, Phys. Status Solidi 220 (2000) 191–194. [13] J.D. Jackson, Classical Electrodynamics, John Wiley, New York, 1999. [14] C. Xue, Q. Huang, S. Deng, Coulomb Green’s function and image potential near a cylindrical diffuse interface, Comput. Phys. Comm. 197 (2015) 153–168. [15] C. Xue, S. Deng, Coulomb Green’s function and image potential near a planar diffuse interface, revisited, Comput. Phys. Comm. 184 (2013) 51–59. [16] Z. Xu, Three-layer dielectric models with spherical cavities for reaction potential calculations, Procedia Eng. 31 (2012) 1033–1038. [17] P. Qin, Z. Xu, W. Cai, D. Jacobs, Image charge methods for a three-dielectric-layer hybrid solvation model of biomolecules, Commun. Comput. Phys. 6 (2009) 955–977. [18] F. Stern, Image potential near a gradual interface between two dielectrics, Phys. Rev. B 17 (1978) 5009–5015. [19] J.W. Perram, M.N. Barber, Coulomb Green’s functions at diffuse interfaces, Mol. Phys. 28 (1974) 131–150. [20] J.R. Clay, N.S. Goel, F.P. Buff, Electrostatics of diffuse anisotropic interfaces III. Point charge and dipole image potentials for air–water and metal–water interfaces, J. Chem. Phys. 56 (1972) 4245–4255. [21] F.P. Buff, N.G. Goel, Electrostatics of diffuse anisotropic interfaces I. Planar Layer Model, J. Chem. Phys. 51 (1969) 4983–4996. [22] F.P. Buff, N.G. Goel, Electrostatics of diffuse anisotropic interfaces II. Effects of long-range diffuseness, J. Chem. Phys. 51 (1969) 5363–5373. [23] D.L. Clements, Fundamental solutions for second order linear elliptic partial differential equations, Comput. Mech. 22 (1998) 26–31. [24] J.L. Movilla, J. Planelles, Image charges in spherical quantum dots with an off-centered impurity: algorithm and numerical results, Comput. Phys. Comm. 170 (2005) 144–152. [25] A. Sihvola, I.V. Lindell, Polarizability and effective permittivity of layered and continuously inhomogeneous dielectric spheres, J. Electromagn. Waves Appl. 3 (1989) 37–60. [26] I.V. Lindell, M.E. Ermutlu, A.H. Sihvola, Electrostatic image theory for layered dielectric sphere, IEE Proc. H 139 (1992) 186–192. [27] A. Sihvola, I.V. Lindell, Transmission line analogy for calculating the effective permittivity of mixtures with spherical multilayer scatterers, J. Electromagn. Waves Appl. 2 (1988) 741–756. [28] J.L. Movilla, J. Planelles, Off-centering of hydrogenic impurities in quantum dots, Phys. Rev. B 71 (2005) 075319. [29] P.M. Morse, H. Feshbach, Methods of Theoretical Physics, McGraw–Hill, New York, 1953. [30] E.W. Hobson, The Theory of Spherical and Ellipsoidal Harmonics, Cambridge University Press, Cambridge, England, 1931. [31] W.R. Smythe, Static and Dynamic Electricity, third ed., Hemisphere, New York, NY, 1989. [32] G. Dassios, Ellipsoidal Harmonics, Eencyclopedia of Mathematics and Its Aapplications, Vol. 146, Cambridge University Press, 2012. [33] J.C.-E. Sten, Ellipsoidal harmonics and their application in electrostatics, J. Electrostat. 64 (2006) 647–654. [34] R. Garmier, J.-P. Barriot, Ellipsoidal harmonic expansions of the gravitational potential: Theory and application, Celestial Mech. Dynam. Astronom. 79 (2000) 235–275. [35] S. Ritter, On the computation of Lame´ functions, of eigenvalues and eigenfunctions of some potential operators, ZAMM·Z, Z. Angew. Math. Mech. 78 (1998) 66–72. [36] H.-J. Dobner, S. Ritter, Verified computation of Lame´ functions with high accuracy, Computing 60 (1998) 81–89. [37] E. Borovitskaya, M.S. Shur, Quantum Dots, World Scientific, Singapore, 2002. [38] A.W.L. Jacak, P. Hawrylak, Quantum Dots, Springer, Berlin, 1998. [39] J.D. Jackson, Classical Electrodynamics, Wiley, New York, 1962.