Many problems in various physical applications are reduced to expressions of the form φ R = ∭ V ρ r | r − R | d 3 r . (1)

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–3] and geophysics [4–7] to the physics of micro- and nanostructures [8–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 B α R ∝ ∂ 2 ∂ R α ∂ R β ∭ V M β r | r − R | d 3 r , (2) 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 E α R ∝ ∂ 2 ∂ R α ∂ R β ∭ V P β r | r − R | d 3 r . (3)

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, u R ∝ ∂ ∂ R ∭ V ε * r | r − R | d 3 r , (4) 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.

With regard to the charged polyhedra, until recently the analytical solutions in closed forms were found for constant [1, 13], linearly [3, 5, 14–16], quadratically [17–20], or cubically [21, 22] coordinate-dependent charge density. Similar solutions for the elastic strain were obtained in works [8, 23–26]. In a special case of a prismatic body, analytic formulas in closed forms are available for polynomials up to fourth degree [27]. For an arbitrary polynomial charge distribution within a polyhedron, a method of analytical determination of the potential was outlined in Ref. [25] and represented in an algorithmic (recursive) form in Refs. [8, 28]. Similar results in two dimensions can be found in Ref. [29].

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 G R , k = ∭ V exp k ⋅ r − R | r − R | d 3 r (5) 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: φ R = G R , ∂ ∂ R ρ R , (6) where differential operator ∂/∂ R 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 G R , k = G 0 R + G α 1 R k α + G α β 2 R k α k β + ⋯ , (7) where G ⁽⁰⁾(R), etc. are expansion coefficients, one obtains from Eq. 6 that φ R = G 0 R ρ R + G α 1 R ∂ ρ R ∂ R α + G α β 2 R ∂ 2 ρ R ∂ R α ∂ R β + ⋯ (8)

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 ⁽⁰⁾(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 Δ φ R = − 4 π ρ R χ R , (9) 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 A ( φ ) , B ( φ ) , C ( φ ) multiplied by simple but non-analytic functions Ω _f (R), 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.

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 Ω f R = ∬ face  f n f ⋅ R − r | R − r | 3 d S , (10) 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 L e R = ∫ edge  e 1 | R − r | d l , (11) 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: G R , k = ∑ f A r f − R , k , n f Ω f R + ∑ f , e B r e − R , k , n f , b f e L e R + ∑ f , e , v C r v − R , k , n f , b f e , l e v | r v − R | , (12) where summations are over polyhedron faces f; over pairs (f, e) of adjacent face f and edge e; and over triples (f, e, v) of adjacent face f, edge e and vertex v. A , B and C are universal functions that do not depend on the choice of the polyhedron.

Expressions for functions A , B and C are (see details of their derivation in Supplementary Section S3): A r , k , n = z 2 ∑ s = 0 ∞ ∑ u = 0 ∞ k z z s − k ‖ 2 z 2 u s ! 2 u + 1 ! s + 2 u + 2 , (13) where z = r ⋅n, k _z = k ⋅n, k ‖ 2 = k 2 − k z 2 , B r , k , n , b = k y z 3 ∑ s = 0 ∞ ∑ u = 0 ∞ k z z s − k ‖ 2 z 2 u s ! 2 u + 2 ! s + 2 u + 3 + y z ∑ q = 0 ∞ ∑ s = 0 ∞ ∑ t = 0 ∞ k y y q k z z s − k x 2 r ⊥ 2 t 2 t − 1 ! ! q ! s ! 2 t ! q + 2 t + 1 q + s + 2 t + 2 2 t ! ! − k x 2 z 3 ∑ q = 0 ∞ ∑ s = 0 ∞ ∑ t = 0 ∞ ∑ u = 0 ∞ k y y q k z z s − k x 2 r ⊥ 2 t − k ‖ 2 z 2 u q + 2 t ! 2 t − 1 ! ! q ! s ! 2 t ! q + 2 t + 2 u + 3 ! 2 t ! ! y q + s + 2 t + 2 u + 4 + k y r ⊥ 2 q + 2 t + 1 q + 2 t + 2 u + 4 q + s + 2 t + 2 u + 5 2 t + 2 , (14) where y = r ⋅b, z = r ⋅n, r ⊥ 2 = y 2 + z 2 , k _y = k ⋅b, k _z = k ⋅n, k x 2 = k 2 − k y 2 − k z 2 , k ‖ 2 = k x 2 + k y 2 , C r , k , n , b , l = k x y z ∑ p = 0 ∞ ∑ q = 0 ∞ ∑ s = 0 ∞ ∑ t = 0 ∞ k x x p k y y q k z z s − k x 2 r ⊥ 2 t p + 2 t + 1 ! q ! s ! p + q + 2 t + 2 p + q + s + 2 t + 3 × p + 2 t ! ! p − 1 ! ! p ! ! p + 2 t + 1 ! ! − k x z 3 ∑ p = 0 ∞ ∑ q = 0 ∞ ∑ s = 0 ∞ ∑ t = 0 ∞ ∑ u = 0 ∞ k x x p k y y q k z z s − k x 2 r ⊥ 2 t − k ‖ 2 z 2 u p + q + 2 t ! p + 2 t ! q ! s ! p + q + 2 t + 2 u + 3 ! p + 2 t ! ! p − 1 ! ! p ! ! p + 2 t + 1 ! ! × [ k x 2 y p + q + 2 t + 1 p + 2 t + 1 p + q + 2 t + 2 u + 4 p + q + s + 2 t + 2 u + 5 − k y p + q + s + 2 t + 2 u + 4 ] , (15) 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 x 2 + k y 2 , 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 non-analytical 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] L e R = ln | R − r e 1 | + | R − r e 2 | + | r e 1 − r e 2 | | R − r e 1 | + | R − r e 2 | − | r e 1 − r e 2 | , (16) 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]: tan Ω f R 2 = a × b ⋅ c | a ‖ b ‖ c | + a ⋅ b | c | + a ⋅ c | b | + b ⋅ c | a | , (17) where a = r _f1 − R, b = r _f2 − R, c = r _f3 − R, and r _f1, r _f2, r _f3 are radius-vectors 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 Ω f R = 2 a t a n 2 P , Q , (18) 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 infinite-series 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.

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 ρ R = exp k ⋅ R . (19)

Substituting it into Eq. 1 and comparing with Eq. 5, one can find that φ R = exp k ⋅ R G R , k . (20)

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 Δ φ x , y , z = δ − 2 φ x + δ , y , z + φ x − δ , y , z + φ x , y + δ , z + φ x , y − δ , z + φ x , y , z + δ + φ x , y , z − δ − 6 φ x , y , z (21) 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–34], for example.

Expression (20) holds also for complex values of k. Hence, the potential of a polyhedral body with density distribution ρ R = cos k ⋅ R (23) is equal to φ R = Re exp i k ⋅ R G R , i k . (24) Similarly, a polyhedral body with density distribution ρ R = sin k ⋅ R (25) produces the potential φ R = Im exp i k ⋅ R G R , i k . (26)

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 ⁰, k ¹ and k ² in Eqs 12–15 for functions A , B and C : A = z 2 1 2 − k ‖ 2 z 2 24 + k z z 3 + k z 2 z 2 8 + O k 3 , (27) B = k y z 3 1 6 + k z z 8 + y z 1 2 − k x 2 r ⊥ 2 48 + k z z 3 + k z 2 z 2 8 + k y y 6 + k y y k z z 12 + k y 2 y 2 24 ) − k x 2 z 3 y 24 + O k 3 , (28) C = k x y z 1 6 + k z z 8 + k y y 12 + k x x 48 + k x z 3 k y 24 + O k 3 . (29)

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 ρ (∂ ³ ρ/∂R _α ∂R _β ∂R _γ , 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 ∂/∂ R which acts on density distribution ρ(R). The components k _x , k _y and k _z are thus replaced with l _α ∂/∂R _α , b _α ∂/∂R _α and n _α ∂/∂R _α , correspondingly. In Eq. 27, factor k ‖ 2 = k 2 − k z 2 is replaced with (δ _αβ − n _α n _β ) ∂ ²/∂R _α ∂R _β , where δ _αβ is the Kronecker delta. In Eq. 28, factor k x 2 = k 2 − k y 2 − k z 2 is replaced with (δ _αβ − b _α b _β − n _α n _β ) ∂ ²/∂R _α ∂R _β .

As a result, we obtain an expression for potential φ(R) in the form of a sum over polyhedron faces, edges and vertices: φ R = ∑ f A φ R , z f , n f Ω f R + ∑ f , e B φ R , y f e , z f , n f , b f e L e R + ∑ f , e , v C φ R , x e v , y f e , z f , n f , b f e , l e v | r v − R | , (30) where x _ev = l _ev ⋅ (r _v − R), y _fe = b _fe ⋅ (r _e − R), z _f = n _f ⋅ (r _f − R), and functions A ( φ ) , B ( φ ) , C ( φ ) are defined as A φ R , z , n = z 2 2 ρ R + z 3 n α 3 ∂ ρ R ∂ R α + z 4 n α n β 6 − δ α β 24 ∂ 2 ρ R ∂ R α ∂ R β , (31) B φ R , y , z , n , b = y z 2 ρ R + z y 2 + z 2 6 b α + y z 3 n α ∂ ρ R ∂ R α + z y 2 + z 2 16 y b α b β + y 2 + z 2 8 z b α n β + y 2 + 9 z 2 48 y n α n β − y 2 + 3 z 2 48 y δ α β ∂ 2 ρ R ∂ R α ∂ R β , (32) C φ R , x , y , z , n , b , l = y z l α 6 ∂ ρ R ∂ R α + z l α x y 48 l β + 2 y 2 + z 2 24 b β + y z 8 n β ∂ 2 ρ R ∂ R α ∂ R β . (33)

Equations 30–33 provide an exact expression for potential φ(R) in a closed form. This expression is valid both within the polyhedral body and outside it, if the density distribution ρ(R) in the body is a quadratic function of coordinates. In the simplified cases of constant or linear function ρ(R), this result is reduced to previously published formulas [1, 8].

Just as for the exponential density distribution, we have performed a numerical test of potential φ(R) obtained by Eqs 30–33. The polyhedron geometry is chosen according to Figure 3A with h = 1 and l = 2. We choose a quadratic polynomial ρ, all coefficients of which are non-zero: ρ(X, Y, Z) = 1 + 2X + 3Y + 4Z + 5X ² + 6Y ² + 7Z ² + 8XY + 9XZ + 10YZ. Potential φ is shown in Figure 3D, and its Laplacian Δφ calculated by difference scheme (21) with δ = 10⁻⁴—in Figure 3E. We have checked that, with accuracy ±10^–4, Δφ = 0 outside the pyramid, and Δφ = − 4πρ within it.

The same recipe can be applied to any polynomial density distribution. For example, if density distribution ρ(X, Y, Z) contains cubic terms (X ³, X ² Y and so on), then it is necessary to take into account also terms proportional to k ³ in series expansions (12)–(15) for functions A , B and C . This, after the replacement k → ∂/∂ R, results in new terms A 3 ( φ ) , B 3 ( φ ) , C 3 ( φ ) that one should add to functions A ( φ ) , B ( φ ) , C ( φ ) , correspondingly: A 3 φ = z 5 n α n β n γ 15 − δ β γ 30 ∂ 3 ρ R ∂ R α ∂ R β ∂ R γ , (34) B 3 φ = z y 2 + z 2 5 y 2 + z 2 240 b α b β b γ + y z y 2 + z 2 20 b α b β n γ + y 2 + z 2 y 2 + 5 z 2 80 b α n β n γ + y z y 2 + 5 z 2 60 n α n β n γ − y 2 + z 2 2 80 δ α β b γ − y z y 2 + 3 z 2 60 δ α β n γ ∂ 3 ρ R ∂ R α ∂ R β ∂ R γ , (35) C 3 φ = z l α y x 2 − 2 y 2 − 5 z 2 360 l β l γ + y 3 y 2 + z 2 120 b β b γ + y z 2 20 n β n γ + x 3 y 2 + z 2 240 l β b γ + x y z 60 l β n γ + z 2 y 2 + z 2 30 b β n γ ∂ 3 ρ R ∂ R α ∂ R β ∂ R γ . (36)

Equations 30–33 with additions 34–36 to functions A ( φ ) , B ( φ ) , C ( φ ) constitute an exact expression for potential φ in a closed form, provided that density distribution ρ(R) within a polyhedron is a cubic function of coordinates.

This method, based on the results of Section 2 and on Eq. 6, allows us to construct explicit expressions for potentials of arbitrary polyhedra, in which the density distribution functions are polynomials of any degree. One can see that complexity of these expressions grows rapidly (though not exponentially) with the degree of polynomial ρ(R). Our method overcomes this difficulty by “packing” all the complexity into expressions (13)–(15) for functions A , B and C .

We propose a method of constructing analytical, closed-form expressions for Newtonian potentials of massive polyhedral bodies, in which the density distributions ρ(X, Y, Z) are polynomials of coordinates X, Y and Z. The method is based on the power-series representation for generating function G(R, k) presented in Section 2. Explicit expressions for the potential in the case of quadratic or cubic polynomial ρ(X, Y, Z) are derived in Section 3. For a quadratic coordinate dependence of density ρ, the potential can be found by Eqs 30–33. For a cubic coordinate dependence of the density, the same equations with additions Eqs 34–36 to the right-hand sides of Eqs 31–33 are valid. It is possible to extend this scheme to polynomials of any degree, and use it for any polyhedral shape of the body. Potential applications of the proposed method are spread from astronomy to nanotechnology.