Image Theory for Neumann Functions in the Prolate Spheroidal Geometry

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.


Introduction
Let  be a smooth closed surface in R 3 , 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  i (r, r s ): where Δ is the Laplace operator, r s is a fixed point in the open domain Ω, (r) is the Dirac delta function, / is the outward normal derivative on the surface , and  S (r) is a given function specified on the surface  satisfying the constraint This constraint follows by applying the divergence theorem to (1a).If we further demand the normal derivative to be constant on , then we arrive at probably the most common boundary condition used in developing a Neumann-Green's function [1]: where || stands for the total surface area of .
The solutions of a Neumann problem when they exist, have the following integral representation (see [2], p.286-287): where ⟨⟩ stands for the weighted mean value of the solution  on ; namely, ⟨⟩ = ∮   S (r)(r)d, which can be chosen to 2 Advances in Mathematical Physics 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 potential  e (r, r s ): Δ e (r, r s ) =  (r − r s ) , r ∈ Ω c , (6a)    e (r, r s ) = − S (r) , r ∈ , e (r, r s ) =  (1/ |r| 2 ) , |r| → +∞, where r s 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]  (r, r s ) = − 1 4 r − r s     +  (r, r s ) , where (r, r s ) is a harmonic function that secures the satisfaction of the boundary conditions.In what follows, (r, r s ) 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][4][5][6][7][8][9][10][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 r k = (/ s ) 2 r s with strength / s , where  is the radius of the sphere, and (iii) a uniform line image with strength −1/ on the line segment between the origin and the Kelvin image point r k .On the other hand, an image system for the interior spherical Neumann function may consist of (i) a Figure 1: A prolate spheroid is centered at the origin, its focal axis is aligned with the -axis, and its interfocal distance is 2.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 r s = ( s ,  s , 0) inside the spheroid.
point image at r k with strength / s and (ii) a line image that extends from r k 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.

Elements of Prolate Spheroidal Harmonics
Prolate spheroid  is defined by where  >  > 0. Here, the focal symmetry axis of the spheroid is aligned with the -axis, and the interfocal distance is 2 with  = √  2 −  2 ; see Figure 1.
In this paper, the prolate spheroidal coordinates (, , ) are defined through where  ∈ [ In particular,  =  b with  b = / 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  = ±, and in particular  = 1,  = −1, and  = 0 correspond to the positive -axis beyond the focal point (0, 0, ), the negative -axis beyond the focal point (0, 0, −), and the -plane, respectively.The coordinate surface of constant  is a plane through the -axis at an angle  to the -plane.The metric coefficients for the prolate spheroidal coordinates are given by The In fact, we have the following orthogonality relation: where   is the Kronecker delta, the differential surface element on   is and the normalization constants   are The surface prolate spheroidal harmonics form a complete set of eigenfunctions over a prolate spheroid   .Therefore, any smooth function  defined over the prolate spheroid   can be expanded in terms of the surface prolate spheroidal harmonics.Furthermore, if  is even with respect to , then it can be expanded in terms of only the even surface prolate spheroidal harmonics; namely, with the coefficients   given by Finally, the expansion of 1/|r−r s | in the prolate spheroidal coordinates is [1,14,15] where  < = min(,  s ),  > = max(,  s ), and

Interior Neumann Functions and Their Image Systems
Given a unit source located at point r s = ( s ,  s ,  s ) = ( s ,  s ,  s ) inside the prolate spheroid , that is, 1 ≤  s <  b , we assume, without loss of any generality, that the source is in the -plane; that is,  s =  s = 0. We consider interior Neumann functions with several different Neumann boundary conditions.

Nonconstant Boundary Condition.
The interior Neumann function considered in this case, denoted by  i a (r, r s ), is the solution of the following boundary value problem: where the normal derivative on  is set to be a constant multiple of the geometrical weighting function   () on .
Compared to the boundary condition (3) which assumes a constant normal derivative on , the boundary condition (20b) not only sounds more physical but also leads to a simpler mathematical formulation for the Neumann function.
In fact, 1/|| is the average of (1/4)  b () on the prolate spheroid  since This identity also guarantees that such a Neumann function exists.In this case, the weighted mean value of a function  on the surface Due to the azimuthal symmetry of the system, the interior Neumann function  i a (r, r s ) can be expressed in terms of the even interior prolate spheroidal harmonics as where   are constants that have to be determined in such a way so as to satisfy the boundary condition (20b).Using (18), we can rewrite  i a (r, r s ) in the neighborhood of the boundary , where  s ≤  ≤  b , as Since the surface normal derivative at any point ( b , , ) on  is the boundary condition (20b) demands that the following equation holds at every point on : Obviously,  00 is the arbitrary constant of the solution  i a (r, r s ) since   0 () ≡ 0. Indeed, integrating (25) over  and using the orthogonality relation (13a), (13b), and (13c), we obtain which is an identity that can be verified using  0 () = ln((+ 1)/( − 1))/2.Next, for each  ≥ 1 and  = 0, . . ., , multiplying both sides of (25) by    (, ), integrating it over , and using the orthogonality relation (13a), (13b), and (13c), we obtain from which we get Hence, the interior Neumann function The arbitrary constant  00 may be fixed by demanding one additional condition, say  i a (0, r s ) = 0.In this case, noting that the origin in the prolate spheroidal coordinates corresponds to the point (1, 0, 0) and that we get (31)

Constant Nonzero Boundary Condition. The interior
Neumann function with a constant nonzero normal derivative on , denoted by  i b (r, r s ), is defined as the solution of the differential equation (20a) with the boundary condition where || stands for the total surface area of .In this case, the term ⟨⟩ in ( 5) is the average value of the solution  on ; namely, ⟨⟩ = ∮  (r)d/||.The interior Neumann function  i b (r, r s ) is also given by ( 22), but now the boundary condition (32) demands that the following equation holds at every point on : Clearly,  00 still serves as the arbitrary constant of the interior Neumann function.Indeed, integrating (33) over  and using (13a), (13b), and (13c), we still have identity (26).Next, for each  ≥ 1 and  = 0, . . ., , multiplying both sides of (33) by    (, ), integrating it over , and using the orthogonality relation (13a), (13b), and (13c), we obtain The right-hand side of (34), however, is no longer zero for all  ≥ 1 and  = 0, . . ., .Instead, noting that we have Therefore, when  > 0 or when  = 0 and  ≥ 1 is odd, we still have (28), but when  = 0 and  ≥ 1 is even, we instead have . (37)

Homogeneous Boundary Condition.
If we demand the homogeneous boundary condition on , then the interior Neumann function, denoted by  i c (r, r s ), is defined as the solution of the following boundary value problem [16]: where |Ω| = (4/3) Then, it is easy to see that Δ 0 (r) = 1/|Ω| and that Therefore, the Neumann function  i c (r, r s ) can be obtained by subtracting  0 (r) from the Neumann function  i a (r, r s ); namely, we have More specifically, we wish to build an image system that represents the following part of the interior Neumann function  i a (r, r s ): 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  at some exterior point r k = ( k , 0,  k ) = ( k ,  k , 0) (since the source is in the -plane, naturally the only point image should be in the -plane as well) as a continuous one-dimensional line image extending from the point r k to infinity along the radial coordinate curve  : (, ) = ( k , 0) with density where  is some constant and () is some continuous function in [ k , +∞), and as a continuous two-dimensional surface image on the confocal prolate spheroid   k :  =  k with density

Advances in Mathematical Physics
where   are constants to be determined later.By reasoning as we did in the case of Dirichlet-Green's function [11], the vanishing of the  = 0 term in (s) implies that the total strength (or in electrostatics, the total surface "charge") on   k is zero, and the vanishing of the  = 1 terms in (s) implies that the distribution on   k is symmetric with its centroid at the origin.Such an image system is graphically illustrated in Figure 2.
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 onedimensional 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 spheroid  generated by this image system is where t  = (  ,  k , 0), s  = ( k ,   ,   ), and d(t  ) is the differential length element on ; that is, d(t  ) = ℎ  (  ,  k )d  .Therefore, we have Using ( 18), we have Furthermore, using ( 18), (44), and the orthogonality relation (13a), (13b), and (13c), we obtain First, comparing the monopole ( = 0) term in both sides of (49), we have from which we get The function () in [ k , +∞) is something one can choose as long as the corresponding improper integrals in (49) all exist.For instance, we can choose 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  ≥ 1 all exist is due to    () = (1/ +1 ) as  → +∞.In fact, we have Next, comparing the dipole ( = 1) terms in both sides of (49), we obtain for  = 0, 1 Using (9a), (9b), and (9c) and noting that  k =  s = 0, we obtain where  1 ( k ) ≡  0 1 ( k ) and Since the point image  lies on the confocal prolate spheroid   k , we obtain which can be rewritten as Equation ( 58) provides us with a nonlinear algebraic equation for the radial coordinate  k of the point image .Unfortunately, this equation does not necessarily have a solution in ( b , +∞).However, if we choose the constant  as one satisfying First, it can be demonstrated that () is continuous on ( b , +∞).Then, note that, for  = 0, 1, as  → +∞, we have   1 () → 0 and thus () → +∞.Next, it is easy to see that ( 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 / s )  or  = ( 0 ( s )/ 0 ( b ))  for some constant  > 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 () for the case of  = 1,  = 2 (so  b = 2/ √ 3), and r s = (0.5, 0, 1) with  =  0 ( s )/ 0 ( 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  k and  k and consequently  k , and we can use (51) to calculate the strength  of the point image, respectively.Finally, comparing the  ≥ 2 terms in both sides of (49), we obtain

Advances in Mathematical Physics
Therefore, the coefficients of the surface image density (44) are given by

Exterior Neumann Functions and Their Image Systems
Given a unit source located at point r s = ( s ,  s ,  s ) = ( s ,  s ,  s ) outside the prolate spheroid , so  s >  b > 1.
Again, we assume the source is in the -plane so  s =  s = 0. We shall consider exterior Neumann functions with several different boundary conditions as well.

Nonconstant Boundary Condition.
The exterior Neumann function considered in this case, denoted by  e a (r, r s ), is the solution of the following boundary value problem: Due to the azimuthal symmetry of the system, the exterior Neumann function  e a (r, r s ) can be expressed in terms of the even exterior prolate spheroidal harmonics as where   are constants that have to be determined by the boundary condition (63b).Using (18), we can rewrite  e a (r, r s ) in the neighborhood of the boundary , where  b ≤  ≤  s , as The boundary condition (63b) demands that the following equation holds at every point on : Integrating (66) over  and using the orthogonality relation (13a), (13b), and (13c), we obtain Therefore, we have Next, for each  ≥ 1 and  = 0, . . ., , multiplying both sides of (66) by    (, ), integrating it over , and using (13a), (13b), and (13c), we obtain In this case, the exterior Neumann function  e b (r, r s ) is still given by ( 64), but now the boundary condition (72) demands that the following equation holds at every point on : Integrating (73) over  and using the orthogonality relation (13a), (13b), and (13c), we obtain Therefore, we still have  00 = 1/ 0 ( s ).Next, for each  ≥ 1 and  = 0, . . ., , multiplying both sides of (73) by    (, ), integrating it over , and using (13a), (13b), and (13c), we obtain As shown by (36), the right-hand side of (75) is no longer zero for all  ≥ 1 and  = 0, . . ., .When  > 0 or when  = 0 and  ≥ 1 is odd, we still have (70); however, when  = 0 and  ≥ 1 is even, we instead have Note that (76) applies to the case of  =  = 0, yielding  00 = 1/ 0 ( s ).

Homogeneous Boundary
In this case, the exterior Neumann function  e c (r, r s ) is still given by ( 64), but now the boundary condition (77) demands that the following equation holds at every point on : from which we get  00 = 0 and (70).Therefore, the exterior Neumann function  e c (r, r s ) is given by First of all, the potential  0 ()/(4) is generated by a focal line image of uniform density: Indeed, the potential generated by this focal line image is The strength (or in electrostatics, the total "charge") of this focal line image is 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 in  e, ref a (r, r s ) may consist of a point image with strength  at some interior point r k = ( k , 0,  k ) = ( k ,  k , 0), a line image extending from the point r k to the focal line along the radial coordinate curve  : (, ) = ( k , 0) with density where  is some constant and () is some continuous function specified in where   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  e a (r, r s ) is graphically illustrated in Figure 4.  where  is a constant such that  2 > 1, then (102) has a solution in (1,  b ).Indeed, define an auxiliary function () : (1,  b ) → R as

) 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 function  i a (r, r s ) given by (29), but the basic idea and the key conclusions apply to other interior Neumann functions as well.

Figure 2 :
Figure 2: Illustration of the image system for the interior Neumann function  i a (r, r s ).
Condition.The exterior Neumann function with the homogeneous boundary condition on , denoted by  e c (r, r s ), is defined as the solution of the differential equation (63a) with the boundary conditions (63c) and    e (r, r s ) = 0,  =  b .

Figure 4 :
Figure 4: Illustration of the image system for the exterior Neumann function  e a (r, r s ).