Generating Function Method for Calculating the Potentials of Inhomogeneous Polyhedra

We propose a method of constructing analytical, closed-form expressions for electrostatic/Newtonian potentials of non-uniform polyhedral bodies, in which the density distributions are polynomials of coordinates. Possible applications of the proposed method are spread from astronomy to nanotechnology. The method is based on the use of the generating function for the potential. Explicit expressions for the potential are derived in the case of quadratic or cubic coordinate dependence of the density within a polyhedral body.


INTRODUCTION
Many problems in various physical applications are reduced to expressions of the form One of them is the determination of the Coulomb (electrostatic) potential from a given distribution of charge. In this respect, ρ(r) is the charge density, and the integration is carried out over the volume V of the charged body. This expression provides also the Newtonian (gravitational) potential of a massive body with a given density distribution ρ(r). Applications of its analytical solutions range from astronomy [1][2][3] and geophysics [4][5][6][7] to the physics of microand nanostructures [8][9][10], as illustrated in Figure 1. Gravitational potentials of asteroids or gravity anomalies have straightforwardly the form of Eq. 1. The magnetic field B(R) of a distributed magnetic moment M(r) (such as a magnetic anomaly or a micromagnet) is proportional to Where α and β are tensor indexes. Similarly, the electric field E(R) of a distributed electric dipole moment P(r) (for example, in a metallic nanoparticle under illumination) is proportional to The elastic displacement vector u(R) due to an inclusion that possesses an eigenstrain ε*(r), such as a lattice misfit in an epitaxial quantum dot, also has a similar form, in the case of an elastically-isotropic medium. Therefore, all the fields considered above can be expressed via potentials having the form of Eq. 1.
There are two generic shapes of the three-dimensional body that enable analytical solutions for potentials φ(R) in closed forms: an ellipsoid [11,12] and a polyhedron [1,12]. Analytical solutions are possible not only for a homogeneously charged ellipsoid or a polyhedron, but also for the case of a polynomial dependence of the charge density ρ(x, y, z) on coordinates [2,8]. We focus here on potentials of polyhedra.
In this work, we present an ultimate expression that unifies all known analytical solutions for Newtonian (or Coulomb) potentials of polyhedron-shaped bodies. For this reason, we introduce a generating function that depends on the geometry of the polyhedron, and on the additional parameter k. One can represent potential φ(R) via the generating function as follows: where differential operator z/zR acts only "to the right," on function ρ(R). The derivation of Eq. 6 is given in Supplementary Section S1. We provide a series expansion of generating function G(R, k) in powers of k. Each term of this expansion is a closed-form expression containing only elementary functions. Writing the series expansion as where G (0) (R), etc. are expansion coefficients, one obtains from Eq. 6 that If ρ(R) is a polynomial, the r.h.s. of Eq. 6 contains only a finite number of terms. Therefore, this equation gives an analytical, closed-form expression for the Newtonian potential φ(R) for any polynomial density distribution ρ(R) within a polyhedral body, provided that coefficients G (0) (R), etc. are expressed in a closed form.
Eq. 8, along with the results of Section 2, provides a method of constructing analytical, closed-form expressions for Newtonian (or electrostatic) potentials of non-uniform polyhedral bodies, in which the density distributions are polynomials of coordinates. The method is universal in the sense that any density distribution can be approximated by a polynomial, as well as any shape of a body can be approximated by a polyhedron. It is worth noting that primitive shapes in the finite-difference method are polyhedra, therefore the results of the present paper can be used in connection with finite-difference schemes.
The results of the present paper provide an insight into analytic properties of solutions of linear inhomogeneous partial differential equations (PDEs) possessing a geometry of a polyhedron. Potential φ(R) obeys Poisson's equation where χ(R) is a characteristic function of the polyhedron: χ(R) 1 for R within the polyhedron, otherwise χ(R) 0. Since function χ(R) is not analytic, the solution φ(R) also must be a non-analytic function of coordinates. However, as we will see later in Eq. 30, the solution can be represented as a sum of analytic functions , L e (R), |r v − R| associated with polyhedron faces, edges and vertices. The latter functions have simple meanings: Ω f (R) is a solid angle subtended by face f at point R, L e (R) is a potential of uniformly charged edge e at R, and |r v − R| is a distance from vertex v to point R. One can consider such a structure of the solution as a useful ansatz for linear PDEs, an inhomogeneous part of which is spread within a polyhedron. The paper is organized as follows. In Section 2, a formula for the generating function in the form of power series is presented. This is the central result of this paper. Then, in Section 3, this formula is tested numerically and applied to calculating the potential. We show how to construct an exact expression in a closed form for potential φ(R) in the case of a polynomial dependence of the charge density on coordinates. Concluding remarks are gathered in Section 4.

SERIES EXPANSION OF THE GENERATING FUNCTION
As shown in Supplementary Section S2, generating function G(R, k) can be represented as a sum of contributions of all faces, edges and vertices of the polyhedron. In this Section, we report on such a representation, and then, in Section 3, we demonstrate how to obtain exact solutions in closed forms on the basis of the expression for the generating function.
First, we introduce some notations for geometrical entities associated with the polyhedron. Let us use symbols f, e and v for enumerating of polyhedron faces, edges and vertices, correspondingly. For each face f, we introduce a unit normal vector n f to the face directed outwards from the polyhedron, see Figure 2. For each pair (f, e) of a face and an edge adjacent to it, we introduce a unit vector b fe that is parallel to face f, perpendicular to edge e, and directed outwards from face f, as shown in Figure 2. For each pair (e, v) of a face and a vertex lying on it, we define unit vector l ev directed from the opposite vertex of edge e towards vertex v ( Figure 2). r f and r e are radius-vectors of arbitrarily chosen points on face f and on edge e, correspondingly. r v is the radius-vector of vertex v. For a given face f, we introduce function where dS is the surface element of face f, and r is the radius-vector of this element. One can easily recognize that Ω f (R) is the signed solid angle subtended by face f at point R: its sign is positive (negative) if the outer (inner) side of the face is seen from point R.
Similarly, for a given edge e, we define function where dl is a line element, and r is its radius-vector. One can interpret L e (R) as a potential at point R of edge e considered as a uniform massive rod with the unit linear mass density. In these notations, generating function G(R, k) acquires the following form derived in Supplementary Section S2: where z r ·n, k z k ·n, where y r ·b, z r ·n, r 2 ⊥ y 2 + z 2 , k y k ·b, k z k ·n, where x r ·l, y r ·b, z r ·n, r 2 ⊥ y 2 + z 2 , k x k ·l, k y k ·b, k z k ·n, k 2 k 2 x + k 2 y , and p!! p(p − 2)(p − 4). . . denotes the double factorial of p [by convention, ( − 1)!! 0!! 1!! 1]. These formulas provide power series expansions of A, B and C as functions of components of vector k. Simultaneously, they are series expansions of A, B and C as functions of r. Series in Eqs 13-15 converge at any values of their parameters, that can be easily proved by the direct comparison test (see Supplementary Section S4). Thus, A, B and C are entire functions of k and r, that is, they have no singularities at any finite k and r.
In contrast, generating function G(R, k) itself is nonanalytical at the surface of the polyhedron. This behavior is related to discontinuity of the Laplacian of potential φ(R) defined by Eq. 1: Δφ(R) − 4πρ(R) within the polyhedron, and Δφ(R) 0 outside. It is evident from Eq. 12 that such a non-analytic behavior of the generating function is due to presence of functions Ω f (R), L e (R) and |r v − R| in the r.h.s. Indeed, Ω f has a discontinuity on face f, L e diverges logarithmically on edge e, and |r v − R| has a cusp-like behavior at vertex v. Hence, representation (12) of generating function G(R, k) allows us to detach its calculational complexity from its non-analytic behavior. All non-analiticity of G(R, k) is "grasped" by functions Ω f , L e and |r v − R| that have simple mathematical expressions in closed forms, see Eqs 16, 17. And all calculational complexity of G(R, k) is contained in fully analytical functions A, B and C.
Functions Ω f (R) and L e (R) can be expressed in closed forms via elementary functions. A convenient representation for L e (R) is [8,13] where r e1 and r e2 are radius-vectors of two ends of edge e. For a triangular face, solid angle Ω f (R) can be represented as follows [30,31]: where a r f1 − R, b r f2 − R, c r f3 − R, and r f1 , r f2 , r f3 are radiusvectors of vertices of face f, counted clockwise as seeing from the outer side of the face. Solid angle Ω f (R) must fall into the range ( − 2π, 2π) and have the same sign as [a ×b] ·c. The appropriate solution of this equation is where P and Q are the numerator and the denominator of the r.h.s. of Eq. 17, and atan2 is a function from the math library of C language [31]. If the face is not triangular, one can either break it into triangles, or use other formulas [13,31] for solid angle Ω f (R). The representation of generating function G(R, k) in the form of Eqs 12-15 constitutes the main result of the present paper. In Section 3, we will demonstrate that the infiniteseries representations for functions A, B and C result in formulas for potential φ(R) containing finite number of terms, when density ρ(R) is a polynomial function of coordinates.

NUMERICAL TESTS AND APPLICATIONS TO EXPONENTIAL, SINUSOIDAL AND POLYNOMIAL DENSITY DISTRIBUTIONS
The expression for generating function G(R, k) presented in Section 2 looks rather complicated. It is necessary therefore to verify it. As a natural numerical test, we have used the generating function to derive potential φ(R) for a certain polyhedron and a certain density distribution ρ(R), and have checked that this potential obeys Poisson's equation.
For the testing purpose, the simplest choice of the density distribution function is Substituting it into Eq. 1 and comparing with Eq. 5, one can find that As an example of a polyhedral body, we choose a pyramid with a square base ( Figure 3A), having height h 1 and lateral size (length of a base edge) l 2. Values of components of vector k are chosen as k x 0.5, k y 0.4 and k z 0.3. When calculating G(R, k), we have kept terms up to ∼ k 13 in power series (13)- (15). Potential φ(R) in the xz-plane passing through the pyramid vertex is shown in Figure 3B. Then, the Laplacian Δφ(R) is calculated numerically by the difference scheme with δ 10 -4 . One can see in Figure 3C that Δφ vanishes outside the pyramid, and depend on coordinate exponentially within it. We have checked that the deviation of the numerical value of Δφ from its exact value Δφ R ( ) −4π exp k · R ( ), within the pyramid 0, outside (22) does not exceed 10 -4 for 1000 randomly chosen points within the cube −2 < x < 2, −2 < y < 2, −2 < z < 2. This deviation is due to errors of the round-off and the difference scheme. The accuracy is sufficient to confirm the correctness of formulas presented in Section 2. The exponential density distribution, Eq. 19, appears in modeling of gravity anomalies [32][33][34], for example. Expression (20) holds also for complex values of k. Hence, the potential of a polyhedral body with density distribution is equal to Similarly, a polyhedral body with density distribution produces the potential Frontiers in Physics | www.frontiersin.org January 2022 | Volume 9 | Article 795693 Now we consider the case of a polynomial density distribution ρ(R) in the body. For definiteness, let ρ(R) be a quadratic function. In order to obtain an expression for potential φ in a closed form, we keep only the terms proportional to k 0 , k 1 and k 2 in Eqs 12-15 for functions A, B and C: B k y z 3 1 6 C k x yz We do not need here the higher terms ( ∼ k 3 , ∼ k 4 and so on), since, according to Eq. 8, they would appear in potential φ as terms containing higher derivatives of density distribution ρ (z 3 ρ/ zR α zR β zR γ , etc.) that vanish in the case of quadratic function ρ(R).
In order to get the potential, as explained in Section 1, each entry of vector k in generating function G(R, k) is replaced with differential operator z/zR which acts on density distribution ρ(R). The components k x , k y and k z are thus replaced with l α z/zR α , b α z/zR α and n α z/zR α , correspondingly. In Eq. 27, factor k 2 k 2 − k 2 z is replaced with (δ αβ − n α n β ) z 2 /zR α zR β , where δ αβ is the Kronecker delta. In Eq. 28, factor k 2 x k 2 − k 2 y − k 2 z is replaced with (δ αβ − b α b β − n α n β ) z 2 /zR α zR β .