# Image Theory for Neumann Functions in the Prolate Spheroidal Geometry

Authors: Changfeng Xue; Robert Edmiston; Shaozhong Deng
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

Interior and exterior Neumann functions for the Laplace operator are derived in terms of prolate spheroidal harmonics with the homogeneous, constant, and nonconstant inhomogeneous boundary conditions. For the interior Neumann functions, an image system is developed to consist of a point image, a line image extending from the point image to infinity along the radial coordinate curve, and a symmetric surface image on the confocal prolate spheroid that passes through the point image. On the other hand, for the exterior Neumann functions, an image system is developed to consist of a point image, a focal line image of uniform density, another line image extending from the point image to the focal line along the radial coordinate curve, and also a symmetric surface image on the confocal prolate spheroid that passes through the point image.


---

## Body

## 1. Introduction

LetS be a smooth closed surface in R3, with Ω being its interior and Ωc its exterior, respectively. An interior Neumann function for the Laplace operator is the solution of the following boundary value problem for the potential Nir,rs:(1a)ΔNir,rs=δr-rs,r∈Ω,(1b)∂∂nNir,rs=gSr,r∈S,where Δ is the Laplace operator, rs is a fixed point in the open domain Ω, δr is the Dirac delta function, ∂/∂n is the outward normal derivative on the surface S, and gSr is a given function specified on the surface S satisfying the constraint(2)∮SgSrds=1.This constraint follows by applying the divergence theorem to (1a). If we further demand the normal derivative to be constant on S, then we arrive at probably the most common boundary condition used in developing a Neumann-Green’s function [1]:(3)∂∂nNir,rs≡1S,r∈S,where S stands for the total surface area of S.The solutions of a Neumann problem(4a)Δur=fr,r∈Ω,(4b)∂∂nur=ψr,r∈S,when they exist, have the following integral representation (see [2], p.286-287):(5)ur=u+∫ΩNir,r′fr′dr′-∮SNir,r′ψr′ds′,where u stands for theweighted mean value of the solution u on S; namely, u=∮SgSrurds, which can be chosen to be zero to simplify the equation. All other solutions to (4a) and (4b) can be obtained by adding an arbitrary constant to this solution.Likely, an exterior Neumann function for the Laplace operator is the solution of the following boundary value problem for the potentialNer,rs: (6a)ΔNer,rs=δr-rs,r∈Ωc,(6b)∂∂nNer,rs=-gSr,r∈S,(6c)Ner,rs=O1/r2,r⟶+∞,where rs now is a given point in the open domain Ωc.Neumann functions are analogous to Green’s functions for Dirichlet problems, so they are often also called Green’s function for the Neumann problem or Green’s function of the second kind. While Dirichlet-Green’s functions are generally used for electrostatic problems where the potential is specified on bounding surfaces, Neumann-Green’s functions are often useful for finding temperature distributions where the bounding surfaces have specified currents. Similar to Dirichlet-Green’s functions, a Neumann-Green function can also be decomposed into a singular and a regular part as [2](7)Nr,rs=-14πr-rs+Rr,rs,where R(r,rs) is a harmonic function that secures the satisfaction of the boundary conditions. In what follows, Rr,rs is also called the reflected part of the Neumann function.An image system for a Neumann function is a system of fictitious sources inside the complement of the fundamental domain that produces a potential equal to the reflected part of the Neumann function. These fictitious sources are commonly called images because they are located not in the real domain of interest for the problem but in its complement. In general, such an image system is not unique. There are all types of images, from isolated point images to continuous distributions of images on lines, curves, and surfaces to combinations of these images. While image theories for Dirichlet-Green’s functions have been studied quite extensively in the literature [3–11], much less work has been done with image theories for Neumann-Green’s functions.In their interesting article [12], Dassios and Sten studied the Neumann function with the constant boundary condition and the corresponding image system in spherical and ellipsoidal geometry. For example, they have found that, in the spherical case, an image system for the exterior spherical Neumann function may consist of (i) a point image at the origin with strength -1, (ii) a point image at the conventional Kelvin image point rk=(a/rs)2rs with strength a/rs, where a is the radius of the sphere, and (iii) a uniform line image with strength -1/a on the line segment between the origin and the Kelvin image point rk. On the other hand, an image system for the interior spherical Neumann function may consist of (i) a point image at rk with strength a/rs and (ii) a line image that extends from rk radially to infinity with field-point dependent charge density.In the present work, we shall study Neumann functions and their image systems for the Laplace operator in the prolate spheroidal geometry. Although in theory the case to be considered here is a special case of that studied in [12], the authors still believe it is much beneficial to study it separately because the Neumann functions and their image systems in the ellipsoidal geometry have to be constructed using ellipsoidal harmonics, while those in the prolate spheroidal geometry can be constructed using prolate spheroidal harmonics, but the ellipsoidal harmonics are much more complicated to handle than the spheroidal ones. In addition, it seems that [12] contains several mistakes of omission and commission [13]. The present paper is organized as follows. Section 2 provides a brief introduction to the prolate spheroidal coordinates and to the prolate spheroidal harmonics. Then, interior and exterior Neumann functions and their image systems for the Laplace operator in the prolate spheroidal geometry are developed in Sections 3 and 4, respectively. Finally, some concluding remarks are given in Section 5.


## 2. Elements of Prolate Spheroidal Harmonics

Prolate spheroidS is defined by(8)x2+y2a2+z2b2=1,where b>a>0. Here, the focal symmetry axis of the spheroid is aligned with the z-axis, and the interfocal distance is 2c with c=b2-a2; see Figure 1.Figure 1
A prolate spheroid is centered at the origin, its focal axis is aligned with thez-axis, and its interfocal distance is 2c. In terms of the prolate spheroidal coordinates ξ,η,ϕ defined in the main text, the prolate spheroid is determined by ξ=ξb. A unit source is located at rs=ξs,ηs,0 inside the spheroid.In this paper, the prolate spheroidal coordinatesξ,η,ϕ are defined through (9a)x=cξ2-11-η2cos⁡ϕ,(9b)y=cξ2-11-η2sin⁡ϕ,(9c)z=cξη,where ξ∈1,+∞ is the radial variable, η∈-1,1 is the angular variable, and ϕ∈0,2π is the azimuthal variable, respectively. The prolate spheroid (8) is called the reference prolate spheroid of this coordinate system. The coordinate surface of constant ξ is a prolate spheroid confocal to the prolate spheroid (8), and such a confocal prolate spheroid, denoted by Sξ, can be written in the Cartesian coordinates as(10)x2+y2c2ξ2-1+z2c2ξ2=1.In particular, ξ=ξb with ξb=b/c is the prolate spheroid (8), and ξ=1 corresponds to the focal line connecting the two focal points, respectively. The coordinate surfaces of constant η are the two sheets of a hyperboloid of revolution with focal points z=±c, and in particular η=1, η=-1, and η=0 correspond to the positive z-axis beyond the focal point 0,0,c, the negative z-axis beyond the focal point 0,0,-c, and the xy-plane, respectively. The coordinate surface of constant ϕ is a plane through the z-axis at an angle ϕ to the xz-plane. The metric coefficients for the prolate spheroidal coordinates are given by(11a)hξ⁡ξ,η=cξ2-η2ξ2-1,(11b)hηξ,η=cξ2-η21-η2,(11c)hϕξ,η=cξ2-11-η2.The interior prolate spheroidal harmonics arePnmξPnmηcos⁡mϕ and PnmξPnmηsin⁡mϕ, and the exterior prolate spheroidal harmonics are QnmξPnmηcos⁡mϕ and QnmξPnmηsin⁡mϕ, for n=0,1,… and m=0,1,…,n, where Pnmx and Qnmx are the associated Legendre functions of the first and second kinds, respectively. In particular, Pn0x≡Pnx and Qn0x≡Qnx are the Legendre functions of the first and second kinds, respectively. Accordingly, the surface prolate spheroidal harmonics are Cnmη,ϕ=Pnmηcos⁡mϕ and Snmη,ϕ=Pnmηsin⁡mϕ, for n=0,1,… and m=0,1,…,n, which are orthogonal over a confocal prolate spheroid Sξ with respect to the following geometrical weighting function:(12)wξη=1c2ξ2-η2ξ2-1.In fact, we have the following orthogonality relation: (13a)∮SξSnmη,ϕSNMη,ϕwξηdsξη,ϕ=γmnδnNδmM,M>0,m>0,(13b)∮SξSnmη,ϕCNMη,ϕwξηdsξη,ϕ=0,(13c)∮SξCnmη,ϕCNMη,ϕwξηdsξη,ϕ=γmnδnNδmM,where δij is the Kronecker delta, the differential surface element on Sξ is(14)dsξη,ϕ=hηhϕdηdϕ=1wξηdηdϕ,and the normalization constants γmn are(15)γmn=2n+m!2n+1n-m!1+δm0π.The surface prolate spheroidal harmonics form a complete set of eigenfunctions over a prolate spheroidSξ. Therefore, any smooth function f defined over the prolate spheroid Sξ can be expanded in terms of the surface prolate spheroidal harmonics. Furthermore, if f is even with respect to ϕ, then it can be expanded in terms of only the even surface prolate spheroidal harmonics; namely,(16)fη,ϕ=∑n=0∞∑m=0ncmnCnmη,ϕ,with the coefficients cmn given by(17)cmn=1γmn∮Sξwξηfη,ϕCnmη,ϕdsξη,ϕ.Finally, the expansion of1/r-rs in the prolate spheroidal coordinates is [1, 14, 15](18)1r-rs=1c∑n=0∞∑m=0nHmnPnmξ<Qnmξ>PnmηsPnmηcos⁡mϕ-ϕs,where ξ<=min⁡ξ,ξs, ξ>=max⁡(ξ,ξs), and(19)Hmn=2n+12-δm0-1mn-m!n+m!2.


## 3. Interior Neumann Functions and Their Image Systems

Given a unit source located at pointrs=xs,ys,zs=ξs,ηs,ϕs inside the prolate spheroid S, that is, 1≤ξs<ξb, we assume, without loss of any generality, that the source is in the xz-plane; that is, ys=ϕs=0. We consider interior Neumann functions with several different Neumann boundary conditions.


### 3.1. Nonconstant Boundary Condition

The interior Neumann function considered in this case, denoted byNair,rs, is the solution of the following boundary value problem: (20a)ΔNir,rs=δr-rs,1≤ξ<ξb,(20b)∂∂nNir,rs=14πwξbη,ξ=ξb,where the normal derivative on S is set to be a constant multiple of the geometrical weighting function wξη on S. Compared to the boundary condition (3) which assumes a constant normal derivative on S, the boundary condition (20b) not only sounds more physical but also leads to a simpler mathematical formulation for the Neumann function. In fact, 1/S is the average of 1/4πwξbη on the prolate spheroid S since(21)14π∮Swξbηds=1.This identity also guarantees that such a Neumann function exists. In this case, the weighted mean value of a function u on the surface S is u=∮Swξbηurds/(4π).Due to the azimuthal symmetry of the system, the interior Neumann functionNair,rs can be expressed in terms of the even interior prolate spheroidal harmonics as(22)Nair,rs=-14πr-rs+14πc∑n=0∞∑m=0nAmnHmnPnmξsPnmηsPnmξCnmη,ϕ,where Amn are constants that have to be determined in such a way so as to satisfy the boundary condition (20b). Using (18), we can rewrite Nair,rs in the neighborhood of the boundary S, where ξs≤ξ≤ξb, as(23)Nair,rs=-14πc∑n=0∞∑m=0nHmnPnmξsPnmηsQnmξ-AmnPnmξCnmη,ϕ.Since the surface normal derivative at any pointξb,η,ϕ on S is(24)∂∂n=1hξξb,η∂∂ξ=a2cwξbη∂∂ξ,the boundary condition (20b) demands that the following equation holds at every point on S:(25)-a2c2wξbη∑n=0∞∑m=0nHmnPnmξsPnmηsQnm′ξb-AmnPnm′ξbCnmη,ϕ=wξbη.Obviously,A00 is the arbitrary constant of the solution Nair,rs since P0′x≡0. Indeed, integrating (25) over S and using the orthogonality relation (13a), (13b), and (13c), we obtain(26)Q0′ξb=-c2a2,which is an identity that can be verified using Q0x=ln⁡(x+1/(x-1))/2.Next, for eachi≥1 and j=0,…,i, multiplying both sides of (25) by Cijη,ϕ, integrating it over S, and using the orthogonality relation (13a), (13b), and (13c), we obtain(27)-a2c2γjiHjiPijξsPijηsQij′ξb-AjiPij′ξb=∮SwξbηCijη,ϕdsξbη,ϕ=0,from which we get(28)Amn=Qnm′ξbPnm′ξb,n≥1.Hence, the interior Neumann function Nair,rs is(29)Nair,rs=-14πr-rs+A004πc+14πc∑n=1∞∑m=0nHmnQnm′ξbPnm′ξbPnmξsPnmηsPnmξCnmη,ϕ.The arbitrary constantA00 may be fixed by demanding one additional condition, say Nai0,rs=0. In this case, noting that the origin in the prolate spheroidal coordinates corresponds to the point 1,0,0 and that(30)Pnm±1=±1n,ifm=0,0,ifm>0,Pn0=-1n/2n-1!!n!!,ifnis  even,0,ifnis  odd,we get(31)A00=crs-∑k=1∞4k+1Q2k′ξbP2k′ξbP2kξsP2kηs-1k2k-1!!2k!!.


### 3.2. Constant Nonzero Boundary Condition

The interior Neumann function with a constant nonzero normal derivative onS, denoted by Nbir,rs, is defined as the solution of the differential equation (20a) with the boundary condition(32)∂∂nNir,rs=1S,ξ=ξb,where S stands for the total surface area of S. In this case, the term u in (5) is the average value of the solution u on S; namely, 〈u〉=∮Surds/S.The interior Neumann functionNbir,rs is also given by (22), but now the boundary condition (32) demands that the following equation holds at every point on S:(33)-a24πc2wξbη∑n=0∞∑m=0nHmnPnmξsPnmηsQnm′ξb-AmnPnm′ξbCnmη,ϕ=1S.Clearly,A00 still serves as the arbitrary constant of the interior Neumann function. Indeed, integrating (33) over S and using (13a), (13b), and (13c), we still have identity (26). Next, for each i≥1 and j=0,…,i, multiplying both sides of (33) by Cijη,ϕ, integrating it over S, and using the orthogonality relation (13a), (13b), and (13c), we obtain(34)-a24πc2γjiHjiPijξsPijηsQij′ξb-AjiPij′ξb=1S∮SCijη,ϕdsξbη,ϕ.The right-hand side of (34), however, is no longer zero for all i≥1 and j=0,…,i. Instead, noting that(35)∮SCijη,ϕdsξbη,ϕ=ac∫-11Pijηξb2-η2dη∫02πcos⁡jϕdϕ,we have(36)∮SCijη,ϕdsξbη,ϕ=2πac∫-11Piηξb2-η2dη,ifj=0,iis  even,0,otherwise.Therefore, when m>0 or when m=0 and n≥1 is odd, we still have (28), but when m=0 and n≥1 is even, we instead have(37)A0n=Qn′ξbPn′ξb+2πc3∫-11Pnηξb2-η2dηaSPnξsPnηsPn′ξb.


### 3.3. Homogeneous Boundary Condition

If we demand the homogeneous boundary condition onS, then the interior Neumann function, denoted by Ncir,rs, is defined as the solution of the following boundary value problem [16]: (38a)ΔNir,rs=δr-rs-1/Ω,1≤ξ<ξb,(38b)∂∂nNir,rs=0,ξ=ξb,where Ω=(4π/3)a2b stands for the volume of Ω. The extra term 1/Ω in the differential equation (38a) is needed to guarantee the existence of such a Neumann function, and in this case the solution to the Neumann problem (4a) and (4b) is still given by (5) but with the term u representing the average value of the solution u in the volume Ω; namely, u=∫Ωu(r)dr/Ω.Letu0r be a function in Ω defined as(39)u0r=r26Ω=x2+y2+z26Ω=c2ξ2+η2-18πa2b.Then, it is easy to see that Δu0r=1/Ω and that(40)∂∂nu0rξ=ξb=18πa2ba2cwξbη∂∂ξc2ξ2+η2-1ξ=ξb=14πwξbη.Therefore, the Neumann function Ncir,rs can be obtained by subtracting u0r from the Neumann function Nair,rs; namely, we have(41)Ncir,rs=Nair,rs-r26Ω.


### 3.4. Image Systems for Interior Neumann Functions

Let us now turn to the construction of image systems for the interior Neumann functions. We will focus ourselves on the interior Neumann functionNair,rs given by (29), but the basic idea and the key conclusions apply to other interior Neumann functions as well.More specifically, we wish to build an image system that represents the following part of the interior Neumann functionNair,rs:(42)Nai,refr,rs=14πc∑n=1∞∑m=0nHmnQnm′ξbPnm′ξbPnmξsPnmηsPnmξCnmη,ϕ.As pointed out earlier, such an image system is not unique. By reasoning as Dassios and Sten did in the case of the ellipsoidal geometry [12], here we consider an image system consisting of a point image with strength Q at some exterior point rk=xk,0,zk=ξk,ηk,0 (since the source is in the xz-plane, naturally the only point image should be in the xz-plane as well) as a continuous one-dimensional line image extending from the point rk to infinity along the radial coordinate curve C:η,ϕ=ηk,0 with density(43)ρt=ρξ,ηk,0=αqξhξξ,ηk,where α is some constant and qξ is some continuous function in ξk,+∞, and as a continuous two-dimensional surface image on the confocal prolate spheroid Sξk:ξ=ξk with density(44)ds=dξk,η,ϕ=wξkη∑n=2∞∑m=0ndmnCnmη,ϕ,where dmn are constants to be determined later. By reasoning as we did in the case of Dirichlet-Green’s function [11], the vanishing of the n=0 term in ds implies that the total strength (or in electrostatics, the total surface “charge”) on Sξk is zero, and the vanishing of the n=1 terms in ds implies that the distribution on Sξk is symmetric with its centroid at the origin. Such an image system is graphically illustrated in Figure 2.Figure 2
Illustration of the image system for the interior Neumann functionNair,rs.As pointed out by Dassios [10], the existence of the continuous one-dimensional distribution of images in the proposed image system is characteristic of the Neumann boundary condition, which in fact was shown 70 years ago by Weiss who studied image systems through applications of Kelvin’s transformation in electricity, magnetism, and hydrodynamics [17, 18]. Dassios further provided an intuitive explanation of the one-dimensional distribution of images reflecting the physics of the underlying problem, that is, the Neumann boundary condition. Since the derivative is the limit of the difference of the solution between two points and the image interpretation of such difference demands a sequence of point images with gaps that are proportional to the distance between the two points, it then becomes clear that the set of point images evolves into a continuous curve as the two points approach each other in order to define the derivative. On the other hand, the existence of the continuous two-dimensional distribution of images over a closed surface in the image system is a direct consequence of the three-dimensional character of the spheroidal or ellipsoidal geometry as it compares with the essentially one-dimensional character of the spherical geometry in which, as mentioned in Section 1, point and line images are enough to represent interior Neumann functions.The potential in the interior of the spheroidS generated by this image system is(45)Nai,imr=-Q4πr-rk-14π∫Cρt′r-t′dlt′-14π∮Sξkds′r-s′dsξkη′,ϕ′,where t′=ξ′,ηk,0, s′=ξk,η′,ϕ′, and dlt′ is the differential length element on C; that is, dlt′=hξξ′,ηkdξ′. Therefore, we have(46)Nai,imr=-Q4πr-rk-α4π∫ξk+∞qξ′r-t′dξ′-14π∮Sξkdξk,η′,ϕ′r-s′dsξkη′,ϕ′.Using (18), we have(47)Q4πr-rk=Q4πc∑n=0∞∑m=0nHmnQnmξkPnmηkPnmξCnmη,ϕ,α4π∫ξk+∞qξ′r-t′dξ′=14πc∑n=0∞∑m=0nHmnα∫ξk+∞qξ′Qnmξ′dξ′PnmηkPnmξCnmη,ϕ.Furthermore, using (18), (44), and the orthogonality relation (13a), (13b), and (13c), we obtain(48)14π∮Sξkdξk,η′,ϕ′r-s′dsξkη′,ϕ′=14πc∑n=2∞∑m=0nHmnγmndmnQnmξkPnmξCnmη,ϕ.So, demanding Nai,refr,rs=Nai,imr, we have(49)14πc∑n=0∞∑m=0nHmnQQnmξk+α∫ξk+∞qξ′Qnmξ′dξ′PnmηkPnmξCnmη,ϕ+14πc∑n=2∞∑m=0nHmnγmndmnQnmξkPnmξCnmη,ϕ=-14πc∑n=1∞∑m=0nHmnQnm′ξbPnm′ξbPnmξsPnmηsPnmξCnmη,ϕ.First, comparing the monopole(n=0) term in both sides of (49), we have(50)QQ0ξk+α∫ξk+∞qξ′Q0ξ′dξ′=0,from which we get(51)Q=-αQ0ξk∫ξk+∞qξ′Q0ξ′dξ′.The functionqξ in ξk,+∞ is something one can choose as long as the corresponding improper integrals in (49) all exist. For instance, we can choose(52)qξ=1-1ξQ0ξ.In this case, the fact that the improper integral in (51) exists can be verified straightforwardly and the fact that the improper integrals in (49) for n≥1 all exist is due to Qnmξ=O1/ξn+1 as ξ→+∞. In fact, we have(53)∫ξk+∞qξ′Q0ξ′dξ′=ξk2lnξk+1ξk-1+12lnξk2-1ξk2-1.Next, comparing the dipole(n=1) terms in both sides of (49), we obtain for m=0,1(54)α∫ξk+∞qξ′Q1mξ′-Q1mξkQ0ξkQ0ξ′dξ′P1mηk=-Q1m′ξbP1m′ξbP1mξsP1mηs.Using (9a), (9b), and (9c) and noting that ϕk=ϕs=0, we obtain(55a)xk=-P11ξkαg11ξkQ11′ξbP11′ξbxs,(55b)zk=-P1ξkαg1ξkQ1′ξbP1′ξbzs,where g1ξk≡g10ξk and(56)g1mξk=∫ξk+∞qξ′Q1mξ′-Q1mξkQ0ξkQ0ξ′dξ′,m=0,1.Since the point imageQ lies on the confocal prolate spheroid Sξk, we obtain(57)xs2c2ξk2-1P11ξk2α2g11ξk2Q11′ξb2P11′ξb2+zs2c2ξk2P1ξk2α2g1ξk2Q1′ξb2P1′ξb2=1,which can be rewritten as(58)xs2a2P11ξb2α2g11ξk2Q11′ξb2P11′ξb2+zs2b2P1ξb2α2g1ξk2Q1′ξb2P1′ξb2=1.Equation (58) provides us with a nonlinear algebraic equation for the radial coordinate ξk of the point image Q. Unfortunately, this equation does not necessarily have a solution in ξb,+∞. However, if we choose the constant α as one satisfying(59)α2=β2xs2a2P11ξb2g11ξb2Q11′ξb2P11′ξb2+zs2b2P1ξb2g1ξb2Q1′ξb2P1′ξb2,where β is a constant such that β2>1, then (58) has a solution in ξb,+∞. Indeed, define an auxiliary function fξ:ξb,+∞→R as(60)fξ=xs2a2P11ξb2α2g11ξ2Q11′ξb2P11′ξb2+zs2b2P1ξb2α2g1ξ2Q1′ξb2P1′ξb2-1.First, it can be demonstrated that f(ξ) is continuous on ξb,+∞. Then, note that, for m=0,1, as ξ→+∞, we have g1mξ→0 and thus fξ→+∞. Next, it is easy to see that fξb=1/β2-1<0. Therefore, (58) must have a solution in ξb,+∞. In addition, we note that the constant β should be chosen in such a way so as to guarantee that the point image is close to the boundary if the source is close to the boundary. For example, we can choose β=ξb/ξsnorβ=Q0ξs/Q0ξbn for some constant n>0.We were unable to prove the uniqueness of such a solution, but numerical investigations have suggested that the solution of (58) in ξb,+∞ should be unique. Figure 3 shows the graph of fξ for the case of a=1, b=2 (so ξb=2/3), and rs=0.5,0,1 with β=Q0(ξs)/Q0(ξb). In this case, the unique solution is about 1.197158. In general, Newton’s method can be used to solve (58) for ξk. Once ξk is found, we can use (55a) and (55b) to calculate xk and zk and consequently ηk, and we can use (51) to calculate the strength Q of the point image, respectively.Figure 3
The graph offξ for the case of a=1, b=2, and rs=0.5,0,1.Finally, comparing then≥2 terms in both sides of (49), we obtain(61)γmndmnQnmξk=-Qnm′ξbPnm′ξbPnmξsPnmηs-QQnmξk+α∫ξk+∞qξ′Qnmξ′dξ′Pnmηk.Therefore, the coefficients of the surface image density (44) are given by(62)dmn=-1γmnQnmξkQnm′ξbPnm′ξbPnmξsPnmηs+QQnmξk+α∫ξk+∞qξ′Qnmξ′dξ′Pnmηk.


## 3.1. Nonconstant Boundary Condition

The interior Neumann function considered in this case, denoted byNair,rs, is the solution of the following boundary value problem: (20a)ΔNir,rs=δr-rs,1≤ξ<ξb,(20b)∂∂nNir,rs=14πwξbη,ξ=ξb,where the normal derivative on S is set to be a constant multiple of the geometrical weighting function wξη on S. Compared to the boundary condition (3) which assumes a constant normal derivative on S, the boundary condition (20b) not only sounds more physical but also leads to a simpler mathematical formulation for the Neumann function. In fact, 1/S is the average of 1/4πwξbη on the prolate spheroid S since(21)14π∮Swξbηds=1.This identity also guarantees that such a Neumann function exists. In this case, the weighted mean value of a function u on the surface S is u=∮Swξbηurds/(4π).Due to the azimuthal symmetry of the system, the interior Neumann functionNair,rs can be expressed in terms of the even interior prolate spheroidal harmonics as(22)Nair,rs=-14πr-rs+14πc∑n=0∞∑m=0nAmnHmnPnmξsPnmηsPnmξCnmη,ϕ,where Amn are constants that have to be determined in such a way so as to satisfy the boundary condition (20b). Using (18), we can rewrite Nair,rs in the neighborhood of the boundary S, where ξs≤ξ≤ξb, as(23)Nair,rs=-14πc∑n=0∞∑m=0nHmnPnmξsPnmηsQnmξ-AmnPnmξCnmη,ϕ.Since the surface normal derivative at any pointξb,η,ϕ on S is(24)∂∂n=1hξξb,η∂∂ξ=a2cwξbη∂∂ξ,the boundary condition (20b) demands that the following equation holds at every point on S:(25)-a2c2wξbη∑n=0∞∑m=0nHmnPnmξsPnmηsQnm′ξb-AmnPnm′ξbCnmη,ϕ=wξbη.Obviously,A00 is the arbitrary constant of the solution Nair,rs since P0′x≡0. Indeed, integrating (25) over S and using the orthogonality relation (13a), (13b), and (13c), we obtain(26)Q0′ξb=-c2a2,which is an identity that can be verified using Q0x=ln⁡(x+1/(x-1))/2.Next, for eachi≥1 and j=0,…,i, multiplying both sides of (25) by Cijη,ϕ, integrating it over S, and using the orthogonality relation (13a), (13b), and (13c), we obtain(27)-a2c2γjiHjiPijξsPijηsQij′ξb-AjiPij′ξb=∮SwξbηCijη,ϕdsξbη,ϕ=0,from which we get(28)Amn=Qnm′ξbPnm′ξb,n≥1.Hence, the interior Neumann function Nair,rs is(29)Nair,rs=-14πr-rs+A004πc+14πc∑n=1∞∑m=0nHmnQnm′ξbPnm′ξbPnmξsPnmηsPnmξCnmη,ϕ.The arbitrary constantA00 may be fixed by demanding one additional condition, say Nai0,rs=0. In this case, noting that the origin in the prolate spheroidal coordinates corresponds to the point 1,0,0 and that(30)Pnm±1=±1n,ifm=0,0,ifm>0,Pn0=-1n/2n-1!!n!!,ifnis  even,0,ifnis  odd,we get(31)A00=crs-∑k=1∞4k+1Q2k′ξbP2k′ξbP2kξsP2kηs-1k2k-1!!2k!!.


## 3.2. Constant Nonzero Boundary Condition

The interior Neumann function with a constant nonzero normal derivative onS, denoted by Nbir,rs, is defined as the solution of the differential equation (20a) with the boundary condition(32)∂∂nNir,rs=1S,ξ=ξb,where S stands for the total surface area of S. In this case, the term u in (5) is the average value of the solution u on S; namely, 〈u〉=∮Surds/S.The interior Neumann functionNbir,rs is also given by (22), but now the boundary condition (32) demands that the following equation holds at every point on S:(33)-a24πc2wξbη∑n=0∞∑m=0nHmnPnmξsPnmηsQnm′ξb-AmnPnm′ξbCnmη,ϕ=1S.Clearly,A00 still serves as the arbitrary constant of the interior Neumann function. Indeed, integrating (33) over S and using (13a), (13b), and (13c), we still have identity (26). Next, for each i≥1 and j=0,…,i, multiplying both sides of (33) by Cijη,ϕ, integrating it over S, and using the orthogonality relation (13a), (13b), and (13c), we obtain(34)-a24πc2γjiHjiPijξsPijηsQij′ξb-AjiPij′ξb=1S∮SCijη,ϕdsξbη,ϕ.The right-hand side of (34), however, is no longer zero for all i≥1 and j=0,…,i. Instead, noting that(35)∮SCijη,ϕdsξbη,ϕ=ac∫-11Pijηξb2-η2dη∫02πcos⁡jϕdϕ,we have(36)∮SCijη,ϕdsξbη,ϕ=2πac∫-11Piηξb2-η2dη,ifj=0,iis  even,0,otherwise.Therefore, when m>0 or when m=0 and n≥1 is odd, we still have (28), but when m=0 and n≥1 is even, we instead have(37)A0n=Qn′ξbPn′ξb+2πc3∫-11Pnηξb2-η2dηaSPnξsPnηsPn′ξb.


## 3.3. Homogeneous Boundary Condition

If we demand the homogeneous boundary condition onS, then the interior Neumann function, denoted by Ncir,rs, is defined as the solution of the following boundary value problem [16]: (38a)ΔNir,rs=δr-rs-1/Ω,1≤ξ<ξb,(38b)∂∂nNir,rs=0,ξ=ξb,where Ω=(4π/3)a2b stands for the volume of Ω. The extra term 1/Ω in the differential equation (38a) is needed to guarantee the existence of such a Neumann function, and in this case the solution to the Neumann problem (4a) and (4b) is still given by (5) but with the term u representing the average value of the solution u in the volume Ω; namely, u=∫Ωu(r)dr/Ω.Letu0r be a function in Ω defined as(39)u0r=r26Ω=x2+y2+z26Ω=c2ξ2+η2-18πa2b.Then, it is easy to see that Δu0r=1/Ω and that(40)∂∂nu0rξ=ξb=18πa2ba2cwξbη∂∂ξc2ξ2+η2-1ξ=ξb=14πwξbη.Therefore, the Neumann function Ncir,rs can be obtained by subtracting u0r from the Neumann function Nair,rs; namely, we have(41)Ncir,rs=Nair,rs-r26Ω.


## 3.4. Image Systems for Interior Neumann Functions

Let us now turn to the construction of image systems for the interior Neumann functions. We will focus ourselves on the interior Neumann functionNair,rs given by (29), but the basic idea and the key conclusions apply to other interior Neumann functions as well.More specifically, we wish to build an image system that represents the following part of the interior Neumann functionNair,rs:(42)Nai,refr,rs=14πc∑n=1∞∑m=0nHmnQnm′ξbPnm′ξbPnmξsPnmηsPnmξCnmη,ϕ.As pointed out earlier, such an image system is not unique. By reasoning as Dassios and Sten did in the case of the ellipsoidal geometry [12], here we consider an image system consisting of a point image with strength Q at some exterior point rk=xk,0,zk=ξk,ηk,0 (since the source is in the xz-plane, naturally the only point image should be in the xz-plane as well) as a continuous one-dimensional line image extending from the point rk to infinity along the radial coordinate curve C:η,ϕ=ηk,0 with density(43)ρt=ρξ,ηk,0=αqξhξξ,ηk,where α is some constant and qξ is some continuous function in ξk,+∞, and as a continuous two-dimensional surface image on the confocal prolate spheroid Sξk:ξ=ξk with density(44)ds=dξk,η,ϕ=wξkη∑n=2∞∑m=0ndmnCnmη,ϕ,where dmn are constants to be determined later. By reasoning as we did in the case of Dirichlet-Green’s function [11], the vanishing of the n=0 term in ds implies that the total strength (or in electrostatics, the total surface “charge”) on Sξk is zero, and the vanishing of the n=1 terms in ds implies that the distribution on Sξk is symmetric with its centroid at the origin. Such an image system is graphically illustrated in Figure 2.Figure 2
Illustration of the image system for the interior Neumann functionNair,rs.As pointed out by Dassios [10], the existence of the continuous one-dimensional distribution of images in the proposed image system is characteristic of the Neumann boundary condition, which in fact was shown 70 years ago by Weiss who studied image systems through applications of Kelvin’s transformation in electricity, magnetism, and hydrodynamics [17, 18]. Dassios further provided an intuitive explanation of the one-dimensional distribution of images reflecting the physics of the underlying problem, that is, the Neumann boundary condition. Since the derivative is the limit of the difference of the solution between two points and the image interpretation of such difference demands a sequence of point images with gaps that are proportional to the distance between the two points, it then becomes clear that the set of point images evolves into a continuous curve as the two points approach each other in order to define the derivative. On the other hand, the existence of the continuous two-dimensional distribution of images over a closed surface in the image system is a direct consequence of the three-dimensional character of the spheroidal or ellipsoidal geometry as it compares with the essentially one-dimensional character of the spherical geometry in which, as mentioned in Section 1, point and line images are enough to represent interior Neumann functions.The potential in the interior of the spheroidS generated by this image system is(45)Nai,imr=-Q4πr-rk-14π∫Cρt′r-t′dlt′-14π∮Sξkds′r-s′dsξkη′,ϕ′,where t′=ξ′,ηk,0, s′=ξk,η′,ϕ′, and dlt′ is the differential length element on C; that is, dlt′=hξξ′,ηkdξ′. Therefore, we have(46)Nai,imr=-Q4πr-rk-α4π∫ξk+∞qξ′r-t′dξ′-14π∮Sξkdξk,η′,ϕ′r-s′dsξkη′,ϕ′.Using (18), we have(47)Q4πr-rk=Q4πc∑n=0∞∑m=0nHmnQnmξkPnmηkPnmξCnmη,ϕ,α4π∫ξk+∞qξ′r-t′dξ′=14πc∑n=0∞∑m=0nHmnα∫ξk+∞qξ′Qnmξ′dξ′PnmηkPnmξCnmη,ϕ.Furthermore, using (18), (44), and the orthogonality relation (13a), (13b), and (13c), we obtain(48)14π∮Sξkdξk,η′,ϕ′r-s′dsξkη′,ϕ′=14πc∑n=2∞∑m=0nHmnγmndmnQnmξkPnmξCnmη,ϕ.So, demanding Nai,refr,rs=Nai,imr, we have(49)14πc∑n=0∞∑m=0nHmnQQnmξk+α∫ξk+∞qξ′Qnmξ′dξ′PnmηkPnmξCnmη,ϕ+14πc∑n=2∞∑m=0nHmnγmndmnQnmξkPnmξCnmη,ϕ=-14πc∑n=1∞∑m=0nHmnQnm′ξbPnm′ξbPnmξsPnmηsPnmξCnmη,ϕ.First, comparing the monopole(n=0) term in both sides of (49), we have(50)QQ0ξk+α∫ξk+∞qξ′Q0ξ′dξ′=0,from which we get(51)Q=-αQ0ξk∫ξk+∞qξ′Q0ξ′dξ′.The functionqξ in ξk,+∞ is something one can choose as long as the corresponding improper integrals in (49) all exist. For instance, we can choose(52)qξ=1-1ξQ0ξ.In this case, the fact that the improper integral in (51) exists can be verified straightforwardly and the fact that the improper integrals in (49) for n≥1 all exist is due to Qnmξ=O1/ξn+1 as ξ→+∞. In fact, we have(53)∫ξk+∞qξ′Q0ξ′dξ′=ξk2lnξk+1ξk-1+12lnξk2-1ξk2-1.Next, comparing the dipole(n=1) terms in both sides of (49), we obtain for m=0,1(54)α∫ξk+∞qξ′Q1mξ′-Q1mξkQ0ξkQ0ξ′dξ′P1mηk=-Q1m′ξbP1m′ξbP1mξsP1mηs.Using (9a), (9b), and (9c) and noting that ϕk=ϕs=0, we obtain(55a)xk=-P11ξkαg11ξkQ11′ξbP11′ξbxs,(55b)zk=-P1ξkαg1ξkQ1′ξbP1′ξbzs,where g1ξk≡g10ξk and(56)g1mξk=∫ξk+∞qξ′Q1mξ′-Q1mξkQ0ξkQ0ξ′dξ′,m=0,1.Since the point imageQ lies on the confocal prolate spheroid Sξk, we obtain(57)xs2c2ξk2-1P11ξk2α2g11ξk2Q11′ξb2P11′ξb2+zs2c2ξk2P1ξk2α2g1ξk2Q1′ξb2P1′ξb2=1,which can be rewritten as(58)xs2a2P11ξb2α2g11ξk2Q11′ξb2P11′ξb2+zs2b2P1ξb2α2g1ξk2Q1′ξb2P1′ξb2=1.Equation (58) provides us with a nonlinear algebraic equation for the radial coordinate ξk of the point image Q. Unfortunately, this equation does not necessarily have a solution in ξb,+∞. However, if we choose the constant α as one satisfying(59)α2=β2xs2a2P11ξb2g11ξb2Q11′ξb2P11′ξb2+zs2b2P1ξb2g1ξb2Q1′ξb2P1′ξb2,where β is a constant such that β2>1, then (58) has a solution in ξb,+∞. Indeed, define an auxiliary function fξ:ξb,+∞→R as(60)fξ=xs2a2P11ξb2α2g11ξ2Q11′ξb2P11′ξb2+zs2b2P1ξb2α2g1ξ2Q1′ξb2P1′ξb2-1.First, it can be demonstrated that f(ξ) is continuous on ξb,+∞. Then, note that, for m=0,1, as ξ→+∞, we have g1mξ→0 and thus fξ→+∞. Next, it is easy to see that fξb=1/β2-1<0. Therefore, (58) must have a solution in ξb,+∞. In addition, we note that the constant β should be chosen in such a way so as to guarantee that the point image is close to the boundary if the source is close to the boundary. For example, we can choose β=ξb/ξsnorβ=Q0ξs/Q0ξbn for some constant n>0.We were unable to prove the uniqueness of such a solution, but numerical investigations have suggested that the solution of (58) in ξb,+∞ should be unique. Figure 3 shows the graph of fξ for the case of a=1, b=2 (so ξb=2/3), and rs=0.5,0,1 with β=Q0(ξs)/Q0(ξb). In this case, the unique solution is about 1.197158. In general, Newton’s method can be used to solve (58) for ξk. Once ξk is found, we can use (55a) and (55b) to calculate xk and zk and consequently ηk, and we can use (51) to calculate the strength Q of the point image, respectively.Figure 3
The graph offξ for the case of a=1, b=2, and rs=0.5,0,1.Finally, comparing then≥2 terms in both sides of (49), we obtain(61)γmndmnQnmξk=-Qnm′ξbPnm′ξbPnmξsPnmηs-QQnmξk+α∫ξk+∞qξ′Qnmξ′dξ′Pnmηk.Therefore, the coefficients of the surface image density (44) are given by(62)dmn=-1γmnQnmξkQnm′ξbPnm′ξbPnmξsPnmηs+QQnmξk+α∫ξk+∞qξ′Qnmξ′dξ′Pnmηk.


## 4. Exterior Neumann Functions and Their Image Systems

Given a unit source located at pointrs=xs,ys,zs=ξs,ηs,ϕs outside the prolate spheroid S, so ξs>ξb>1. Again, we assume the source is in the xz-plane so ys=ϕs=0. We shall consider exterior Neumann functions with several different boundary conditions as well.


### 4.1. Nonconstant Boundary Condition

The exterior Neumann function considered in this case, denoted byNaer,rs, is the solution of the following boundary value problem: (63a)ΔNer,rs=δr-rs,ξ>ξb,(63b)∂∂nNer,rs=-14πwξbη,ξ=ξb,(63c)Ner,rs=O1/r2,r⟶+∞.Due to the azimuthal symmetry of the system, the exterior Neumann functionNaer,rs can be expressed in terms of the even exterior prolate spheroidal harmonics as(64)Naer,rs=-14πr-rs+14πc∑n=0∞∑m=0nBmnHmnQnmξsPnmηsQnmξCnmη,ϕ,where Bmn are constants that have to be determined by the boundary condition (63b). Using (18), we can rewrite Naer,rs in the neighborhood of the boundary S, where ξb≤ξ≤ξs, as(65)Naer,rs=-14πc∑n=0∞∑m=0nHmnQnmξsPnmηsPnmξ-BmnQnmξCnmη,ϕ.The boundary condition (63b) demands that the following equation holds at every point on S:(66)a2c2wξbη∑n=0∞∑m=0nHmnQnmξsPnmηsPnm′ξb-BmnQnm′ξbCnmη,ϕ=wξbη.Integrating (66) over S and using the orthogonality relation (13a), (13b), and (13c), we obtain(67)-a2c2γ00H00Q0ξsP0ηsQ0′ξbB00=∮Swξbηdsξbη,ϕ=γ00.Therefore, we have(68)B00=-c2a2Q0ξsP0ηsQ0′ξb=1Q0ξs.Next, for eachi≥1 and j=0,…,i, multiplying both sides of (66) by Cijη,ϕ, integrating it over S, and using (13a), (13b), and (13c), we obtain(69)a2c2γjiHjiQijξsPijηsPij′ξb-BjiQij′ξb=∮SwξbηCijη,ϕdsξbη,ϕ=0,from which we get(70)Bmn=Pnm′ξbQnm′ξb,n≥1.Hence, the exterior Neumann function Naer,rs is given by(71)Naer,rs=-14πr-rs+Q0ξ4πc+14πc∑n=1∞∑m=0nHmnPnm′ξbQnm′ξbQnmξsPnmηsQnmξCnmη,ϕ.


### 4.2. Constant Nonzero Boundary Condition

The exterior Neumann function with a constant nonzero normal derivative onS, denoted by Nber,rs, is defined as the solution of the differential equation (63a) with the boundary conditions (63c) and(72)∂∂nNer,rs=-1S,ξ=ξb.In this case, the exterior Neumann function Nber,rs is still given by (64), but now the boundary condition (72) demands that the following equation holds at every point on S:(73)a24πc2wξbη∑n=0∞∑m=0nHmnQnmξsPnmηsPnm′ξb-BmnQnm′ξbCnmη,ϕ=1S.Integrating (73) over S and using the orthogonality relation (13a), (13b), and (13c), we obtain(74)-a24πc2γ00H00Q0ξsP0ηsQ0′ξbB00=1S∮Sdsξbη,ϕ=1.Therefore, we still have B00=1/Q0ξs. Next, for each i≥1 and j=0,…,i, multiplying both sides of (73) by Cijη,ϕ, integrating it over S, and using (13a), (13b), and (13c), we obtain(75)a24πc2γjiHjiQijξsPijηsPij′ξb-BjiQij′ξb=1S∮SCijη,ϕdsξbη,ϕ.As shown by (36), the right-hand side of (75) is no longer zero for all i≥1 and j=0,…,i. When m>0 or when m=0 and n≥1 is odd, we still have (70); however, when m=0 and n≥1 is even, we instead have(76)B0n=Pn′ξbQn′ξb-2πc3∫-11Pnηξb2-η2dηaSQnξsPnηsQn′ξb.Note that (76) applies to the case of n=m=0, yielding B00=1/Q0ξs.


### 4.3. Homogeneous Boundary Condition

The exterior Neumann function with the homogeneous boundary condition onS, denoted by Ncer,rs, is defined as the solution of the differential equation (63a) with the boundary conditions (63c) and(77)∂∂nNer,rs=0,ξ=ξb.In this case, the exterior Neumann function Ncer,rs is still given by (64), but now the boundary condition (77) demands that the following equation holds at every point on S:(78)a24πc2wξbη∑n=0∞∑m=0nHmnQnmξsPnmηsPnm′ξb-BmnQnm′ξbCnmη,ϕ=0,from which we get B00=0 and (70). Therefore, the exterior Neumann function Ncer,rs is given by(79)Ncer,rs=-14πr-rs+14πc∑n=1∞∑m=0nHmnPnm′ξbQnm′ξbQnmξsPnmηsQnmξCnmη,ϕ.


### 4.4. Image Systems for Exterior Neumann Functions

Again, we will focus ourselves on an image system for the exterior Neumann functionNaer,rs given by (71). In particular, let us build an image system that can represent the reflected part of the Neumann function Naer,rs; namely,(80)Nae,refr,rs=Q0ξ4πc+14πc∑n=1∞∑m=0nHmnPnm′ξbQnm′ξbQnmξsPnmηsQnmξCnmη,ϕ.First of all, the potential Q0ξ/4πc is generated by a focal line image of uniform density:(81)μz=-12c,-c≤z≤c.Indeed, the potential generated by this focal line image is(82)-14π∫-ccμzr-0,z,0dz=18π∫-111r-1,η′,0dη′=18πc∑n=0∞2n+1∫-11Pnη′dη′Pn1QnξPnη=Q0ξ4πc.The strength (or in electrostatics, the total “charge”) of this focal line image is(83)∫-ccμzdz=-∫-cc12cdz=-1,indicating that this focal line image corresponds to the point image of strength -1 that has to be put at the origin in the case of the sphere, as mentioned in Section 1.Next, working as in the interior case, we conjecture that an image system for the remaining expansion inNae,refr,rs may consist of a point image with strength Q at some interior point rk=xk,0,zk=ξk,ηk,0, a line image extending from the point rk to the focal line along the radial coordinate curve C:η,ϕ=ηk,0 with density(84)ρt=ρξ,ηk,0=αqξhξξ,ηk,where α is some constant and qξ is some continuous function specified in 1,ξk, and a surface image on the confocal prolate spheroid Sξk:ξ=ξk with density(85)ds=dξk,η,ϕ=wξkη∑n=2∞∑m=0ndmnCnmη,ϕ,where dmn are constants to be determined later. Again, the surface image has a total strength of zero and is symmetric with its centroid at the origin. The proposed image system for the exterior Neumann function Naer,rs is graphically illustrated in Figure 4.Figure 4
Illustration of the image system for the exterior Neumann functionNaer,rs.The potential in the exterior of the spheroid generated by this image system (including the focal line image of uniform densityμz) is(86)Nae,imr=Q0ξ4πc-Q4πr-rk-14π∫Cρt′r-t′dlt′-14π∮Sξkds′r-s′dsξkη′,ϕ′,where t′=ξ′,ηk,0 and s′=ξk,η′,ϕ′; that is,(87)Nae,imr=Q0ξ4πc-Q4πr-rk-α4π∫1ξkqξ′r-t′dξ′-14π∮Sξkdξk,η′,ϕ′r-s′dsξkη′,ϕ′.Using (18), we have(88)Q4πr-rk=Q4πc∑n=0∞∑m=0nHmnPnmξkPnmηkQnmξCnmη,ϕ,and (89)α4π∫1ξkqξ′r-t′dξ′=14πc∑n=0∞∑m=0nHmnα∫1ξkqξ′Pnmξ′dξ′PnmηkQnmξCnmη,ϕ.Furthermore, using (18), (85), and the orthogonality relation (13a), (13b), and (13c), we can obtain(90)14π∮Sξkdξk,η′,ϕ′r-s′dsξkη′,ϕ′=14πc∑n=2∞∑m=0nHmnγmndmnPnmξkQnmξCnmη,ϕ.Demanding Nae,imr=Nae,refr,rs, we have(91)14πc∑n=0∞∑m=0nHmnQPnmξk+α∫1ξkqξ′Pnmξ′dξ′PnmηkQnmξCnmη,ϕ+14πc∑n=2∞∑m=0nHmnγmndmnPnmξkQnmξCnmη,ϕ=-14πc∑n=1∞∑m=0nHmnPnm′ξbQnm′ξbQnmξsPnmηsQnmξCnmη,ϕ.First, comparing the monopole(n=0) term in both sides of (91), we obtain(92)Q=-α∫1ξkqξ′dξ′.Next, comparing the dipole (n=1) terms in both sides of (91), we get for m=0,1(93)αg1mξkP1mηk=-P1m′ξbQ1m′ξbQ1mξsP1mηs,where(94)g1mξk=∫1ξkqξ′P1mξ′-P1mξkdξ′,m=0,1.Finally, comparing the n≥2 terms in both sides of (91), we obtain(95)γmndmnPnmξk=-Pnm′ξbQnm′ξbQnmξsPnmηs-QPnmξk+α∫1ξkqξ′Pnmξ′dξ′Pnmηk.Therefore, the coefficients of the surface image density (85) are given by(96)dmn=-1γmnPnmξkPnm′ξbQnm′ξbQnmξsPnmηs+QPnmξk+α∫1ξkqξ′Pnmξ′dξ′Pnmηk.To finalize the image system, we have essentially now three unknown values plus one unknown function to decide:Q, ξk, ηk, and qξ as a continuous function in 1,ξk. They are intercorrelated through the underdetermined system of (92) and (93), so there could be infinitely many different solutions, each giving us a different image system. Below we consider an approach that is again similar to that proposed by Dassios and Sten in [12]; namely, we choose the unknown function qξ and then determine the other three unknown values. Since, in this case, only proper integrals on 1,ξk are involved, any continuous function qξ on 1,ξk can guarantee the existence of all integrals in (91). For instance, we can simply choose(97)qξ≡1,ξ∈1,ξb.In this case, by (92), we immediately have(98)Q=α1-ξk.On the other hand, (93) can be rewritten as (99a)xk=-P11ξkαg11ξkQ11ξsP11ξsP11′ξbQ11′ξbxs,(99b)zk=-P1ξkαg1ξkQ1ξsP1ξsP1′ξbQ1′ξbzs,where(100)g1ξk=ξk-ξk22-12,g11ξk=ξk2-1-ξkξk2-12-ln⁡ξk2-1+ξk2.Since the point image Q lies on the confocal prolate spheroid Sξk, we obtain(101)xs2c2ξk2-1P11ξk2α2g11ξk2Q11ξs2P11ξs2P11′ξb2Q11′ξb2+zs2c2ξk2P1ξk2α2g1ξk2Q1ξs2P1ξs2P1′ξb2Q1′ξb2=1,which can be rewritten as(102)xs2a2P11ξb2α2g11ξk2Q11ξs2P11ξs2P11′ξb2Q11′ξb2+zs2b2P1ξb2α2g1ξk2Q1ξs2P1ξs2P1′ξb2Q1′ξb2=1.Equation (102) is a nonlinear algebraic equation for the radial coordinate ξk of the point image Q. Unfortunately, this equation does not necessarily have a solution in (1,ξb). However, if we choose the constant α as one satisfying(103)α2=β2xs2a2P11ξb2g11ξb2Q11ξs2P11ξs2P11′ξb2Q11′ξb2+zs2b2P1ξb2g1ξb2Q1ξs2P1ξs2P1′ξb2Q1′ξb2,where β is a constant such that β2>1, then (102) has a solution in 1,ξb. Indeed, define an auxiliary function gξ:1,ξb→R as (104)gξ=xs2a2P11ξb2α2g11ξ2Q11ξs2P11ξs2P11′ξb2Q11′ξb2+zs2b2P1ξb2α2g1ξ2Q1ξs2P1ξs2P1′ξb2Q1′ξb2-1.First of all, g(ξ) is continuous on (1,ξb). Indeed, when q(ξ)≡1, it is easy to see from (94) that both g1(ξ) and g11(ξ) are continuous on 1<ξ<ξb and that g1(ξ)<0 and g11(ξ)<0 for 1<ξ<ξb since Pnm(x) is strictly increasing on (1,+∞) for all n≥1. Then, note that, for m=0,1, as ξ→1+, we have g1m(ξ)→0 and thus g(ξ)→+∞. Next, it is easy to see that g(ξb)=1/β2-1<0. Therefore, (102) must have a solution in (1,ξb). In addition, we note again that the constant β should be chosen in such a way so as to guarantee that the point image is close to the boundary if the source is close to the boundary. For example, we can choose β=ξs/ξbn or β=Q0(ξb)/Q0(ξs)n for some constant n>0.We were unable to prove the uniqueness of such a solution, but numerical investigations have suggested that (102) has a unique solution in (1,ξb). Figure 5 shows the graph of g(ξ) for the case of a=1, b=2 (so ξb=2/3), and rs=(1,0,2) with β=Q0(ξb)/Q0(ξs). In this case, the unique solution is about 1.12575182. Once ξk is found, we can use (99a) and (99b) to calculate xk and zk and consequently ηk, and we can use (92) to calculate the strength Q of the point image, respectively.Figure 5
The graph ofg(ξ) for the case of a=1, b=2, and rs=(1,0,2).


## 4.1. Nonconstant Boundary Condition

The exterior Neumann function considered in this case, denoted byNaer,rs, is the solution of the following boundary value problem: (63a)ΔNer,rs=δr-rs,ξ>ξb,(63b)∂∂nNer,rs=-14πwξbη,ξ=ξb,(63c)Ner,rs=O1/r2,r⟶+∞.Due to the azimuthal symmetry of the system, the exterior Neumann functionNaer,rs can be expressed in terms of the even exterior prolate spheroidal harmonics as(64)Naer,rs=-14πr-rs+14πc∑n=0∞∑m=0nBmnHmnQnmξsPnmηsQnmξCnmη,ϕ,where Bmn are constants that have to be determined by the boundary condition (63b). Using (18), we can rewrite Naer,rs in the neighborhood of the boundary S, where ξb≤ξ≤ξs, as(65)Naer,rs=-14πc∑n=0∞∑m=0nHmnQnmξsPnmηsPnmξ-BmnQnmξCnmη,ϕ.The boundary condition (63b) demands that the following equation holds at every point on S:(66)a2c2wξbη∑n=0∞∑m=0nHmnQnmξsPnmηsPnm′ξb-BmnQnm′ξbCnmη,ϕ=wξbη.Integrating (66) over S and using the orthogonality relation (13a), (13b), and (13c), we obtain(67)-a2c2γ00H00Q0ξsP0ηsQ0′ξbB00=∮Swξbηdsξbη,ϕ=γ00.Therefore, we have(68)B00=-c2a2Q0ξsP0ηsQ0′ξb=1Q0ξs.Next, for eachi≥1 and j=0,…,i, multiplying both sides of (66) by Cijη,ϕ, integrating it over S, and using (13a), (13b), and (13c), we obtain(69)a2c2γjiHjiQijξsPijηsPij′ξb-BjiQij′ξb=∮SwξbηCijη,ϕdsξbη,ϕ=0,from which we get(70)Bmn=Pnm′ξbQnm′ξb,n≥1.Hence, the exterior Neumann function Naer,rs is given by(71)Naer,rs=-14πr-rs+Q0ξ4πc+14πc∑n=1∞∑m=0nHmnPnm′ξbQnm′ξbQnmξsPnmηsQnmξCnmη,ϕ.


## 4.2. Constant Nonzero Boundary Condition

The exterior Neumann function with a constant nonzero normal derivative onS, denoted by Nber,rs, is defined as the solution of the differential equation (63a) with the boundary conditions (63c) and(72)∂∂nNer,rs=-1S,ξ=ξb.In this case, the exterior Neumann function Nber,rs is still given by (64), but now the boundary condition (72) demands that the following equation holds at every point on S:(73)a24πc2wξbη∑n=0∞∑m=0nHmnQnmξsPnmηsPnm′ξb-BmnQnm′ξbCnmη,ϕ=1S.Integrating (73) over S and using the orthogonality relation (13a), (13b), and (13c), we obtain(74)-a24πc2γ00H00Q0ξsP0ηsQ0′ξbB00=1S∮Sdsξbη,ϕ=1.Therefore, we still have B00=1/Q0ξs. Next, for each i≥1 and j=0,…,i, multiplying both sides of (73) by Cijη,ϕ, integrating it over S, and using (13a), (13b), and (13c), we obtain(75)a24πc2γjiHjiQijξsPijηsPij′ξb-BjiQij′ξb=1S∮SCijη,ϕdsξbη,ϕ.As shown by (36), the right-hand side of (75) is no longer zero for all i≥1 and j=0,…,i. When m>0 or when m=0 and n≥1 is odd, we still have (70); however, when m=0 and n≥1 is even, we instead have(76)B0n=Pn′ξbQn′ξb-2πc3∫-11Pnηξb2-η2dηaSQnξsPnηsQn′ξb.Note that (76) applies to the case of n=m=0, yielding B00=1/Q0ξs.


## 4.3. Homogeneous Boundary Condition

The exterior Neumann function with the homogeneous boundary condition onS, denoted by Ncer,rs, is defined as the solution of the differential equation (63a) with the boundary conditions (63c) and(77)∂∂nNer,rs=0,ξ=ξb.In this case, the exterior Neumann function Ncer,rs is still given by (64), but now the boundary condition (77) demands that the following equation holds at every point on S:(78)a24πc2wξbη∑n=0∞∑m=0nHmnQnmξsPnmηsPnm′ξb-BmnQnm′ξbCnmη,ϕ=0,from which we get B00=0 and (70). Therefore, the exterior Neumann function Ncer,rs is given by(79)Ncer,rs=-14πr-rs+14πc∑n=1∞∑m=0nHmnPnm′ξbQnm′ξbQnmξsPnmηsQnmξCnmη,ϕ.


## 4.4. Image Systems for Exterior Neumann Functions

Again, we will focus ourselves on an image system for the exterior Neumann functionNaer,rs given by (71). In particular, let us build an image system that can represent the reflected part of the Neumann function Naer,rs; namely,(80)Nae,refr,rs=Q0ξ4πc+14πc∑n=1∞∑m=0nHmnPnm′ξbQnm′ξbQnmξsPnmηsQnmξCnmη,ϕ.First of all, the potential Q0ξ/4πc is generated by a focal line image of uniform density:(81)μz=-12c,-c≤z≤c.Indeed, the potential generated by this focal line image is(82)-14π∫-ccμzr-0,z,0dz=18π∫-111r-1,η′,0dη′=18πc∑n=0∞2n+1∫-11Pnη′dη′Pn1QnξPnη=Q0ξ4πc.The strength (or in electrostatics, the total “charge”) of this focal line image is(83)∫-ccμzdz=-∫-cc12cdz=-1,indicating that this focal line image corresponds to the point image of strength -1 that has to be put at the origin in the case of the sphere, as mentioned in Section 1.Next, working as in the interior case, we conjecture that an image system for the remaining expansion inNae,refr,rs may consist of a point image with strength Q at some interior point rk=xk,0,zk=ξk,ηk,0, a line image extending from the point rk to the focal line along the radial coordinate curve C:η,ϕ=ηk,0 with density(84)ρt=ρξ,ηk,0=αqξhξξ,ηk,where α is some constant and qξ is some continuous function specified in 1,ξk, and a surface image on the confocal prolate spheroid Sξk:ξ=ξk with density(85)ds=dξk,η,ϕ=wξkη∑n=2∞∑m=0ndmnCnmη,ϕ,where dmn are constants to be determined later. Again, the surface image has a total strength of zero and is symmetric with its centroid at the origin. The proposed image system for the exterior Neumann function Naer,rs is graphically illustrated in Figure 4.Figure 4
Illustration of the image system for the exterior Neumann functionNaer,rs.The potential in the exterior of the spheroid generated by this image system (including the focal line image of uniform densityμz) is(86)Nae,imr=Q0ξ4πc-Q4πr-rk-14π∫Cρt′r-t′dlt′-14π∮Sξkds′r-s′dsξkη′,ϕ′,where t′=ξ′,ηk,0 and s′=ξk,η′,ϕ′; that is,(87)Nae,imr=Q0ξ4πc-Q4πr-rk-α4π∫1ξkqξ′r-t′dξ′-14π∮Sξkdξk,η′,ϕ′r-s′dsξkη′,ϕ′.Using (18), we have(88)Q4πr-rk=Q4πc∑n=0∞∑m=0nHmnPnmξkPnmηkQnmξCnmη,ϕ,and (89)α4π∫1ξkqξ′r-t′dξ′=14πc∑n=0∞∑m=0nHmnα∫1ξkqξ′Pnmξ′dξ′PnmηkQnmξCnmη,ϕ.Furthermore, using (18), (85), and the orthogonality relation (13a), (13b), and (13c), we can obtain(90)14π∮Sξkdξk,η′,ϕ′r-s′dsξkη′,ϕ′=14πc∑n=2∞∑m=0nHmnγmndmnPnmξkQnmξCnmη,ϕ.Demanding Nae,imr=Nae,refr,rs, we have(91)14πc∑n=0∞∑m=0nHmnQPnmξk+α∫1ξkqξ′Pnmξ′dξ′PnmηkQnmξCnmη,ϕ+14πc∑n=2∞∑m=0nHmnγmndmnPnmξkQnmξCnmη,ϕ=-14πc∑n=1∞∑m=0nHmnPnm′ξbQnm′ξbQnmξsPnmηsQnmξCnmη,ϕ.First, comparing the monopole(n=0) term in both sides of (91), we obtain(92)Q=-α∫1ξkqξ′dξ′.Next, comparing the dipole (n=1) terms in both sides of (91), we get for m=0,1(93)αg1mξkP1mηk=-P1m′ξbQ1m′ξbQ1mξsP1mηs,where(94)g1mξk=∫1ξkqξ′P1mξ′-P1mξkdξ′,m=0,1.Finally, comparing the n≥2 terms in both sides of (91), we obtain(95)γmndmnPnmξk=-Pnm′ξbQnm′ξbQnmξsPnmηs-QPnmξk+α∫1ξkqξ′Pnmξ′dξ′Pnmηk.Therefore, the coefficients of the surface image density (85) are given by(96)dmn=-1γmnPnmξkPnm′ξbQnm′ξbQnmξsPnmηs+QPnmξk+α∫1ξkqξ′Pnmξ′dξ′Pnmηk.To finalize the image system, we have essentially now three unknown values plus one unknown function to decide:Q, ξk, ηk, and qξ as a continuous function in 1,ξk. They are intercorrelated through the underdetermined system of (92) and (93), so there could be infinitely many different solutions, each giving us a different image system. Below we consider an approach that is again similar to that proposed by Dassios and Sten in [12]; namely, we choose the unknown function qξ and then determine the other three unknown values. Since, in this case, only proper integrals on 1,ξk are involved, any continuous function qξ on 1,ξk can guarantee the existence of all integrals in (91). For instance, we can simply choose(97)qξ≡1,ξ∈1,ξb.In this case, by (92), we immediately have(98)Q=α1-ξk.On the other hand, (93) can be rewritten as (99a)xk=-P11ξkαg11ξkQ11ξsP11ξsP11′ξbQ11′ξbxs,(99b)zk=-P1ξkαg1ξkQ1ξsP1ξsP1′ξbQ1′ξbzs,where(100)g1ξk=ξk-ξk22-12,g11ξk=ξk2-1-ξkξk2-12-ln⁡ξk2-1+ξk2.Since the point image Q lies on the confocal prolate spheroid Sξk, we obtain(101)xs2c2ξk2-1P11ξk2α2g11ξk2Q11ξs2P11ξs2P11′ξb2Q11′ξb2+zs2c2ξk2P1ξk2α2g1ξk2Q1ξs2P1ξs2P1′ξb2Q1′ξb2=1,which can be rewritten as(102)xs2a2P11ξb2α2g11ξk2Q11ξs2P11ξs2P11′ξb2Q11′ξb2+zs2b2P1ξb2α2g1ξk2Q1ξs2P1ξs2P1′ξb2Q1′ξb2=1.Equation (102) is a nonlinear algebraic equation for the radial coordinate ξk of the point image Q. Unfortunately, this equation does not necessarily have a solution in (1,ξb). However, if we choose the constant α as one satisfying(103)α2=β2xs2a2P11ξb2g11ξb2Q11ξs2P11ξs2P11′ξb2Q11′ξb2+zs2b2P1ξb2g1ξb2Q1ξs2P1ξs2P1′ξb2Q1′ξb2,where β is a constant such that β2>1, then (102) has a solution in 1,ξb. Indeed, define an auxiliary function gξ:1,ξb→R as (104)gξ=xs2a2P11ξb2α2g11ξ2Q11ξs2P11ξs2P11′ξb2Q11′ξb2+zs2b2P1ξb2α2g1ξ2Q1ξs2P1ξs2P1′ξb2Q1′ξb2-1.First of all, g(ξ) is continuous on (1,ξb). Indeed, when q(ξ)≡1, it is easy to see from (94) that both g1(ξ) and g11(ξ) are continuous on 1<ξ<ξb and that g1(ξ)<0 and g11(ξ)<0 for 1<ξ<ξb since Pnm(x) is strictly increasing on (1,+∞) for all n≥1. Then, note that, for m=0,1, as ξ→1+, we have g1m(ξ)→0 and thus g(ξ)→+∞. Next, it is easy to see that g(ξb)=1/β2-1<0. Therefore, (102) must have a solution in (1,ξb). In addition, we note again that the constant β should be chosen in such a way so as to guarantee that the point image is close to the boundary if the source is close to the boundary. For example, we can choose β=ξs/ξbn or β=Q0(ξb)/Q0(ξs)n for some constant n>0.We were unable to prove the uniqueness of such a solution, but numerical investigations have suggested that (102) has a unique solution in (1,ξb). Figure 5 shows the graph of g(ξ) for the case of a=1, b=2 (so ξb=2/3), and rs=(1,0,2) with β=Q0(ξb)/Q0(ξs). In this case, the unique solution is about 1.12575182. Once ξk is found, we can use (99a) and (99b) to calculate xk and zk and consequently ηk, and we can use (92) to calculate the strength Q of the point image, respectively.Figure 5
The graph ofg(ξ) for the case of a=1, b=2, and rs=(1,0,2).


## 5. Conclusions

In the present work, we present Neumann functions for the Laplace operator in the case where the fundamental domain is either the interior or the exterior of a prolate spheroid. Three Neumann boundary conditions are considered, that is, the homogeneous, constant, and nonconstant inhomogeneous boundary conditions. Then, we give image systems for such Neumann functions. For the interior Neumann functions, an image system is developed to consist of a point image, a line image extending from the point image to infinity along the radial coordinate curve, and a symmetric surface image on the confocal prolate spheroid that passes through the point image. On the other hand, for the exterior Neumann functions, an image system is developed to consist of a point image, a focal line image of uniform density, another line image extending from the point image to the focal line along the radial coordinate curve, and also a symmetric surface image on the confocal prolate spheroid that passes through the point image.



---

*Source: 7683929-2018-03-11.xml*