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.

1 Introduction

Many problems in various physical applications are reduced to expressions of the form

$\phi \left(\mathbf{R}\right)=\underset{V}{\iiint }\frac{\rho \left(\mathbf{r}\right)}{|\mathbf{r}-\mathbf{R}|}{\mathrm{d}}^3\mathbf{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_{\alpha }\left(\mathbf{R}\right)\propto \frac{{\partial }^2}{\partial R_{\alpha }\partial R_{\beta }}\underset{V}{\iiint }\frac{M_{\beta }\left(\mathbf{r}\right)}{|\mathbf{r}-\mathbf{R}|}{\mathrm{d}}^3\mathbf{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_{\alpha }\left(\mathbf{R}\right)\propto \frac{{\partial }^2}{\partial R_{\alpha }\partial R_{\beta }}\underset{V}{\iiint }\frac{P_{\beta }\left(\mathbf{r}\right)}{|\mathbf{r}-\mathbf{R}|}{\mathrm{d}}^3\mathbf{r}.(3)$

FIGURE 1

FIGURE 1. Possible applications of analytical solutions of Eq. 1: from astronomy (gravitational fields of asteroids) and Earth sciences (gravity and magnetic anomalies) to the physics of micro-and nanostructures (magnetic fields of nanomagnets, plasmonic resonance of metallic nanoparticles, elastic fields of inclusions such as epitaxial quantum dots).

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,

$\mathbf{u}\left(\mathbf{R}\right)\propto \frac{\partial }{\partial \mathbf{R}}\underset{V}{\iiint }\frac{{\epsilon }^{*}\left(\mathbf{r}\right)}{|\mathbf{r}-\mathbf{R}|}{\mathrm{d}}^3\mathbf{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\left(\mathbf{R},\mathbf{k}\right)=\underset{V}{\iiint }\frac{\exp \left[\mathbf{k}\cdot \left(\mathbf{r}-\mathbf{R}\right)\right]}{|\mathbf{r}-\mathbf{R}|}{\mathrm{d}}^3\mathbf{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:

$\phi \left(\mathbf{R}\right)=G\left(\mathbf{R},\frac{\partial }{\partial \mathbf{R}}\right)\rho \left(\mathbf{R}\right),(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\left(\mathbf{R},\mathbf{k}\right)=G^{\left(0\right)}\left(\mathbf{R}\right)+G_{\alpha }^{\left(1\right)}\left(\mathbf{R}\right)k_{\alpha }+G_{\alpha \beta }^{\left(2\right)}\left(\mathbf{R}\right)k_{\alpha }{k}_{\beta }+\cdots ,(7)$

where G⁽⁰⁾(R), etc. are expansion coefficients, one obtains from Eq. 6 that

$\phi \left(\mathbf{R}\right)=G^{\left(0\right)}\left(\mathbf{R}\right)\rho \left(\mathbf{R}\right)+G_{\alpha }^{\left(1\right)}\left(\mathbf{R}\right)\frac{\partial \rho \left(\mathbf{R}\right)}{\partial R_{\alpha }}+G_{\alpha \beta }^{\left(2\right)}\left(\mathbf{R}\right)\frac{{\partial }^2\rho \left(\mathbf{R}\right)}{\partial R_{\alpha }\partial R_{\beta }}+\cdots (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

$\mathrm{\Delta }\phi \left(\mathbf{R}\right)=-4\pi \rho \left(\mathbf{R}\right)\chi \left(\mathbf{R}\right),(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 ${\mathcal{A}}^{(\phi )}$, ${\mathcal{B}}^{(\phi )}$, ${\mathcal{C}}^{(\phi )}$ 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.

2 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

${\mathrm{\Omega }}_f\left(\mathbf{R}\right)=\underset{\text{face\hspace{0.17em}}f}{\iint }\frac{{\mathbf{n}}_f\cdot \left(\mathbf{R}-\mathbf{r}\right)}{|\mathbf{R}-\mathbf{r}{|}^3}\mathrm{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\left(\mathbf{R}\right)=\underset{\text{edge\hspace{0.17em}}e}{\int }\frac{1}{|\mathbf{R}-\mathbf{r}|}\mathrm{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:

$\begin{array}{l}\hfill G\left(\mathbf{R},\mathbf{k}\right)=\sum _f\mathcal{A}\left({\mathbf{r}}_f-\mathbf{R},\mathbf{k},{\mathbf{n}}_f\right){\mathrm{\Omega }}_f\left(\mathbf{R}\right)+\sum _{\left(f,e\right)}\mathcal{B}\left({\mathbf{r}}_e-\mathbf{R},\mathbf{k},{\mathbf{n}}_f,{\mathbf{b}}_{fe}\right)L_e\left(\mathbf{R}\right)\\ \hfill +\sum _{\left(f,e,v\right)}\mathcal{C}\left({\mathbf{r}}_v-\mathbf{R},\mathbf{k},{\mathbf{n}}_f,{\mathbf{b}}_{fe},{\mathbf{l}}_{ev}\right)|{\mathbf{r}}_v-\mathbf{R}|,\end{array}(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. $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ are universal functions that do not depend on the choice of the polyhedron.

FIGURE 2

FIGURE 2. Unit vectors and radius-vectors associated with faces, edges, and vertices of the polyhedron. Vector n_f is the outward normal vector to face f. Vector b_fe lies in the plane of face f and is directed out of this face perpendicular to edge e. Vector l_ev is directed along edge e toward vertex v from the opposite vertex. 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, and R is the radius-vector of the observation point.

Expressions for functions $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ are (see details of their derivation in Supplementary Section S3):

$\mathcal{A}\left(\mathbf{r},\mathbf{k},\mathbf{n}\right)=z^2\sum _{s=0}^{\infty }\sum _{u=0}^{\infty }\frac{{\left(k_{z}z\right)}^s{\left(-k_{\Vert }^2{z}^2\right)}^u}{s!\left(2u+1\right)!\left(s+2u+2\right)},(13)$

where z = r ⋅n, k_z = k ⋅n, $k_{\Vert }^2=k^2-k_z^2$,

$\begin{array}{l}\hfill \mathcal{B}\left(\mathbf{r},\mathbf{k},\mathbf{n},\mathbf{b}\right)=k_y{z}^3\sum _{s=0}^{\infty }\sum _{u=0}^{\infty }\frac{{\left(k_{z}z\right)}^s{\left(-k_{\Vert }^2{z}^2\right)}^u}{s!\left(2u+2\right)!\left(s+2u+3\right)}\\ \hfill +yz\sum _{q=0}^{\infty }\sum _{s=0}^{\infty }\sum _{t=0}^{\infty }\frac{{\left(k_{y}y\right)}^q{\left(k_{z}z\right)}^s{\left(-k_x^2{r}_{\perp }^2\right)}^t\left(2t-1\right)!!}{q!s!\left(2t\right)!\left(q+2t+1\right)\left(q+s+2t+2\right)\left(2t\right)!!}\\ \hfill -k_x^2{z}^3\sum _{q=0}^{\infty }\sum _{s=0}^{\infty }\sum _{t=0}^{\infty }\sum _{u=0}^{\infty }\frac{{\left(k_{y}y\right)}^q{\left(k_{z}z\right)}^s{\left(-k_x^2{r}_{\perp }^2\right)}^t{\left(-k_{\Vert }^2{z}^2\right)}^u\left(q+2t\right)!\left(2t-1\right)!!}{q!s!\left(2t\right)!\left(q+2t+2u+3\right)!\left(2t\right)!!}\left[\frac{y}{q+s+2t+2u+4}\right.\\ \hfill \left.+\frac{k_y{r}_{\perp }^2\left(q+2t+1\right)}{\left(q+2t+2u+4\right)\left(q+s+2t+2u+5\right)\left(2t+2\right)}\right],\end{array}(14)$

where y = r ⋅b, z = r ⋅n, $r_{\perp }^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_{\Vert }^2=k_x^2+k_y^2$,

$\begin{array}{l}\hfill \mathcal{C}\left(\mathbf{r},\mathbf{k},\mathbf{n},\mathbf{b},\mathbf{l}\right)=k_{x}yz\sum _{p=0}^{\infty }\sum _{q=0}^{\infty }\sum _{s=0}^{\infty }\sum _{t=0}^{\infty }\frac{{\left(k_{x}x\right)}^p{\left(k_{y}y\right)}^q{\left(k_{z}z\right)}^s{\left(-k_x^2{r}_{\perp }^2\right)}^t}{\left(p+2t+1\right)!q!s!\left(p+q+2t+2\right)\left(p+q+s+2t+3\right)}\\ \hfill \times \frac{\left(p+2t\right)!!\left(p-1\right)!!}{p!!\left(p+2t+1\right)!!}\\ \hfill -k_x{z}^3\sum _{p=0}^{\infty }\sum _{q=0}^{\infty }\sum _{s=0}^{\infty }\sum _{t=0}^{\infty }\sum _{u=0}^{\infty }\frac{{\left(k_{x}x\right)}^p{\left(k_{y}y\right)}^q{\left(k_{z}z\right)}^s{\left(-k_x^2{r}_{\perp }^2\right)}^t{\left(-k_{\Vert }^2{z}^2\right)}^u\left(p+q+2t\right)!}{\left(p+2t\right)!q!s!\left(p+q+2t+2u+3\right)!}\frac{\left(p+2t\right)!!\left(p-1\right)!!}{p!!\left(p+2t+1\right)!!}\\ \hfill \times \left[\frac{k_x^{2}y\left(p+q+2t+1\right)}{\left(p+2t+1\right)\left(p+q+2t+2u+4\right)\left(p+q+s+2t+2u+5\right)}-\frac{k_y}{p+q+s+2t+2u+4}\right],\end{array}(15)$

where x = r ⋅l, y = r ⋅b, z = r ⋅n, $r_{\perp }^2=y^2+z^2$, k_x = k ⋅l, k_y = k ⋅b, k_z = k ⋅n, $k_{\Vert }^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 $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ as functions of components of vector k. Simultaneously, they are series expansions of $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{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, $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{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 $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{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\left(\mathbf{R}\right)=\ln \frac{|\mathbf{R}-{\mathbf{r}}_{e1}|+|\mathbf{R}-{\mathbf{r}}_{e2}|+|{\mathbf{r}}_{e1}-{\mathbf{r}}_{e2}|}{|\mathbf{R}-{\mathbf{r}}_{e1}|+|\mathbf{R}-{\mathbf{r}}_{e2}|-|{\mathbf{r}}_{e1}-{\mathbf{r}}_{e2}|},(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 \frac{{\mathrm{\Omega }}_f\left(\mathbf{R}\right)}{2}=\frac{\left[\mathbf{a}\times \mathbf{b}\right]\cdot \mathbf{c}}{|\mathbf{a}\Vert \mathbf{b}\Vert \mathbf{c}|+\left(\mathbf{a}\cdot \mathbf{b}\right)|\mathbf{c}|+\left(\mathbf{a}\cdot \mathbf{c}\right)|\mathbf{b}|+\left(\mathbf{b}\cdot \mathbf{c}\right)|\mathbf{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

${\mathrm{\Omega }}_f\left(\mathbf{R}\right)=2\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{2}\left(P,Q\right),(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 $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ result in formulas for potential φ(R) containing finite number of terms, when density ρ(R) is a polynomial function of coordinates.

3 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

$\rho \left(\mathbf{R}\right)=\exp \left(\mathbf{k}\cdot \mathbf{R}\right).(19)$

Substituting it into Eq. 1 and comparing with Eq. 5, one can find that

$\phi \left(\mathbf{R}\right)=\exp \left(\mathbf{k}\cdot \mathbf{R}\right)G\left(\mathbf{R},\mathbf{k}\right).(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 $\sim 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

$\begin{array}{l}\hfill \mathrm{\Delta }\phi \left(x,y,z\right)={\delta }^{-2}\left[\phi \left(x+\delta ,y,z\right)+\phi \left(x-\delta ,y,z\right)+\phi \left(x,y+\delta ,z\right)+\phi \left(x,y-\delta ,z\right)\right.\\ \hfill +\left.\phi \left(x,y,z+\delta \right)+\phi \left(x,y,z-\delta \right)-6\phi \left(x,y,z\right)\right]\end{array}(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

$\mathrm{\Delta }\phi \left(\mathbf{R}\right)=\left\{\begin{array}{cc}-4\pi \exp \left(\mathbf{k}\cdot \mathbf{R}\right),\hfill & \text{within\hspace{0.17em}the\hspace{0.17em}pyramid}\hfill \\ 0,\hfill & \text{outside}\hfill \end{array}\right.(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.

FIGURE 3

FIGURE 3. Examples of potentials of an inhomogeneously charged polyhedron. (A) An example of a polyhedron: a pyramid with height h = 1 and lateral size (length of a base edge) l = 2. (B,C) Potential φ(X, Z) and its Laplacian Δφ(X, Z) in the case of an exponential density distribution ρ(R) =  exp(k ⋅R) with k_x = 0.5, k_y = 0.4 and k_z = 0.3. (D,E) Potential φ(X, Z) and its Laplacian Δφ(X, Z) in the case of a quadratic density distribution ρ(X, Y, Z) = 1 + 2X + 3Y + 4Z + 5X² + 6Y² + 7Z² + 8XY + 9XZ + 10YZ. Plots b, c, d, e are presented in the section by the XZ-plane that passes through the pyramid apex (Y = 0). Green triangles depict the pyramid boundary in the same section.

Expression (20) holds also for complex values of k. Hence, the potential of a polyhedral body with density distribution

$\rho \left(\mathbf{R}\right)=\cos \left(\mathbf{k}\cdot \mathbf{R}\right)(23)$

is equal to

$\phi \left(\mathbf{R}\right)=\text{Re}\left[\exp \left(i\mathbf{k}\cdot \mathbf{R}\right)G\left(\mathbf{R},i\mathbf{k}\right)\right].(24)$

Similarly, a polyhedral body with density distribution

$\rho \left(\mathbf{R}\right)=\sin \left(\mathbf{k}\cdot \mathbf{R}\right)(25)$

produces the potential

$\phi \left(\mathbf{R}\right)=\text{Im}\left[\exp \left(i\mathbf{k}\cdot \mathbf{R}\right)G\left(\mathbf{R},i\mathbf{k}\right)\right].(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 $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$:

$\mathcal{A}=z^2\left(\frac{1}{2}-\frac{k_{\Vert }^2{z}^2}{24}+\frac{k_{z}z}{3}+\frac{k_z^2{z}^2}{8}\right)+O\left(k^3\right),(27)$

$\begin{array}{l}\hfill \mathcal{B}=k_y{z}^3\left(\frac{1}{6}+\frac{k_{z}z}{8}\right)+yz\frac{1}{2}-\frac{k_x^2{r}_{\perp }^2}{48}+\frac{k_{z}z}{3}+\frac{k_z^2{z}^2}{8}+\frac{k_{y}y}{6}+\frac{k_{y}y{k}_{z}z}{12}+\frac{k_y^2{y}^2}{24})-k_x^2{z}^3\frac{y}{24}+O\left(k^3\right),\end{array}(28)$

$\mathcal{C}=k_{x}yz\left(\frac{1}{6}+\frac{k_{z}z}{8}+\frac{k_{y}y}{12}+\frac{k_{x}x}{48}\right)+k_x{z}^3\frac{k_y}{24}+O\left(k^3\right).(29)$

We do not need here the higher terms ($\sim k^3$, $\sim 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_{\Vert }^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:

$\begin{array}{l}\hfill \phi \left(\mathbf{R}\right)=\sum _f{\mathcal{A}}^{\left(\phi \right)}\left(\mathbf{R},z_f,{\mathbf{n}}_f\right){\mathrm{\Omega }}_f\left(\mathbf{R}\right)+\sum _{\left(f,e\right)}{\mathcal{B}}^{\left(\phi \right)}\left(\mathbf{R},y_{fe},z_f,{\mathbf{n}}_f,{\mathbf{b}}_{fe}\right)L_e\left(\mathbf{R}\right)\\ \hfill +\sum _{\left(f,e,v\right)}{\mathcal{C}}^{\left(\phi \right)}\left(\mathbf{R},x_{ev},y_{fe},z_f,{\mathbf{n}}_f,{\mathbf{b}}_{fe},{\mathbf{l}}_{ev}\right)|{\mathbf{r}}_v-\mathbf{R}|,\end{array}(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 ${\mathcal{A}}^{(\phi )}$, ${\mathcal{B}}^{(\phi )}$, ${\mathcal{C}}^{(\phi )}$ are defined as

${\mathcal{A}}^{\left(\phi \right)}\left(\mathbf{R},z,\mathbf{n}\right)=\frac{z^2}{2}\rho \left(\mathbf{R}\right)+\frac{z^3{n}_{\alpha }}{3}\frac{\partial \rho \left(\mathbf{R}\right)}{\partial R_{\alpha }}+z^4\left(\frac{n_{\alpha }{n}_{\beta }}{6}-\frac{{\delta }_{\alpha \beta }}{24}\right)\frac{{\partial }^2\rho \left(\mathbf{R}\right)}{\partial R_{\alpha }\partial R_{\beta }},(31)$

$\begin{array}{l}\hfill {\mathcal{B}}^{\left(\phi \right)}\left(\mathbf{R},y,z,\mathbf{n},\mathbf{b}\right)=\frac{yz}{2}\rho \left(\mathbf{R}\right)+z\left(\frac{y^2+z^2}{6}{b}_{\alpha }+\frac{yz}{3}{n}_{\alpha }\right)\frac{\partial \rho \left(\mathbf{R}\right)}{\partial R_{\alpha }}\\ \hfill +z\frac{y^2+z^2}{16}y{b}_{\alpha }{b}_{\beta }+\frac{y^2+z^2}{8}z{b}_{\alpha }{n}_{\beta }+\frac{y^2+9z^2}{48}y{n}_{\alpha }{n}_{\beta }-\frac{y^2+3z^2}{48}y{\delta }_{\alpha \beta }\frac{{\partial }^2\rho \left(\mathbf{R}\right)}{\partial R_{\alpha }\partial R_{\beta }},\end{array}(32)$

${\mathcal{C}}^{\left(\phi \right)}\left(\mathbf{R},x,y,z,\mathbf{n},\mathbf{b},\mathbf{l}\right)=\frac{yz{l}_{\alpha }}{6}\frac{\partial \rho \left(\mathbf{R}\right)}{\partial R_{\alpha }}+z{l}_{\alpha }\left(\frac{xy}{48}{l}_{\beta }+\frac{2y^2+z^2}{24}{b}_{\beta }+\frac{yz}{8}{n}_{\beta }\right)\frac{{\partial }^2\rho \left(\mathbf{R}\right)}{\partial R_{\alpha }\partial R_{\beta }}.(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 $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$. This, after the replacement k → ∂/∂R, results in new terms ${\mathcal{A}}_3^{(\phi )}$, ${\mathcal{B}}_3^{(\phi )}$, ${\mathcal{C}}_3^{(\phi )}$ that one should add to functions ${\mathcal{A}}^{(\phi )}$, ${\mathcal{B}}^{(\phi )}$, ${\mathcal{C}}^{(\phi )}$, correspondingly:

${\mathcal{A}}_3^{\left(\phi \right)}=z^5{n}_{\alpha }\left(\frac{n_{\beta }{n}_{\gamma }}{15}-\frac{{\delta }_{\beta \gamma }}{30}\right)\frac{{\partial }^3\rho \left(\mathbf{R}\right)}{\partial R_{\alpha }\partial R_{\beta }\partial R_{\gamma }},(34)$

$\begin{array}{l}\hfill {\mathcal{B}}_3^{\left(\phi \right)}=z\left[\frac{\left(y^2+z^2\right)\left(5y^2+z^2\right)}{240}{b}_{\alpha }{b}_{\beta }{b}_{\gamma }+\frac{yz\left(y^2+z^2\right)}{20}{b}_{\alpha }{b}_{\beta }{n}_{\gamma }+\frac{\left(y^2+z^2\right)\left(y^2+5z^2\right)}{80}{b}_{\alpha }{n}_{\beta }{n}_{\gamma }\right.\\ \hfill +\frac{yz\left(y^2+5z^2\right)}{60}{n}_{\alpha }{n}_{\beta }{n}_{\gamma }-\frac{{\left(y^2+z^2\right)}^2}{80}{\delta }_{\alpha \beta }{b}_{\gamma }-\frac{yz\left(y^2+3z^2\right)}{60}{\delta }_{\alpha \beta }{n}_{\gamma }\frac{{\partial }^3\rho \left(\mathbf{R}\right)}{\partial R_{\alpha }\partial R_{\beta }\partial R_{\gamma }},\end{array}(35)$

$\begin{array}{l}\hfill {\mathcal{C}}_3^{\left(\phi \right)}=z{l}_{\alpha }\left[\frac{y\left(x^2-2y^2-5z^2\right)}{360}{l}_{\beta }{l}_{\gamma }+\frac{y\left(3y^2+z^2\right)}{120}{b}_{\beta }{b}_{\gamma }\right.\\ \hfill +\frac{y{z}^2}{20}{n}_{\beta }{n}_{\gamma }+\frac{x\left(3y^2+z^2\right)}{240}{l}_{\beta }{b}_{\gamma }+\frac{xyz}{60}{l}_{\beta }{n}_{\gamma }+\frac{z\left(2y^2+z^2\right)}{30}{b}_{\beta }{n}_{\gamma }\frac{{\partial }^3\rho \left(\mathbf{R}\right)}{\partial R_{\alpha }\partial R_{\beta }\partial R_{\gamma }}.\end{array}(36)$

Equations 30–33 with additions 34–36 to functions ${\mathcal{A}}^{(\phi )}$, ${\mathcal{B}}^{(\phi )}$, ${\mathcal{C}}^{(\phi )}$ 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 $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$.

4 Conclusion

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.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

AN: conceptualization, analytical results, numerical calculations and writing the manuscript.

Funding

This work is funded by the Ministry of Science and Higher Education of the Russian Federation, grant 075-15-2020-797 (13.1902.21.0024).

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2021.795693/full#supplementary-material

References

1. Werner RA. The Gravitational Potential of a Homogeneous Polyhedron or Don't Cut Corners. Celestial Mech Dyn Astr (1994) 59:253–78. doi:10.1007/BF00692875

CrossRef Full Text | Google Scholar

2. Chandrasekhar S. Ellipsoidal Figures of Equilibrium. New York: Dover Books on Mathematics (1987).

Google Scholar

3. Conway JT. Analytical Solution from Vector Potentials for the Gravitational Field of a General Polyhedron. Celest Mech Dyn Astr (2015) 121:17–38. doi:10.1007/s10569-014-9588-x

CrossRef Full Text | Google Scholar

4. Nagy D. The Gravitational Attraction of a Right Rectangular Prism. Geophysics (1966) 31:362–71. doi:10.1190/1.1439779

CrossRef Full Text | Google Scholar

5. Holstein H. Gravimagnetic Anomaly Formulas for Polyhedra of Spatially Linear Media. Geophysics (2003) 68:157–67. doi:10.1190/1.1543203

CrossRef Full Text | Google Scholar

6. Abtahi SM, Pedersen LB, Kamm J, Kalscheuer T. A New Reference Model for 3D Inversion of Airborne Magnetic Data in Hilly Terrain - A Case Study from Northern Sweden. Geophysics (2018) 83:B1–B12. doi:10.1190/geo2016-0331.1

CrossRef Full Text | Google Scholar

7. Nibisha VA, Ramamma B, Sastry SR, Chakravarthi V. Forward Modelling: Magnetic Anomalies of Arbitrarily Magnetized 2D Fault Sources with Analytically Defined Fault Planes. J Earth Syst Sci (2021) 130:130. doi:10.1007/s12040-021-01634-x

CrossRef Full Text | Google Scholar

8. Nenashev AV, Dvurechenskii AV. Elastic Strain Field Due to an Inclusion of a Polyhedral Shape with a Non-uniform Lattice Misfit. J Appl Phys (2017) 121:125102. doi:10.1063/1.4978750

CrossRef Full Text | Google Scholar

9. Eshelby JD. The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems. Proc R Soc Lond A (1957) 241:376–96. doi:10.1098/rspa.1957.0133

CrossRef Full Text | Google Scholar

10. Mura T. Micromechanics of Defects in Solids. Mechanics of Elastic and Inelastic Solids. Dordrecht, Netherlands: Springer (1987).

Google Scholar

11. Landau LD, Lifshitz EM. The Classical Theory of Fields. Course of Theoretical Physics. Amsterdam: Elsevier Science (2013).

Google Scholar

12. MacMillan W. The Theory of the Potential. New York: Dover books on physics and mathematical physics Dover Publications (1958).

Google Scholar

13. Werner R, Scheeres D. Exterior Gravitation of a Polyhedron Derived and Compared with Harmonic and Mascon Gravitation Representations of Asteroid 4769 Castalia. Celestial Mech Dyn Astr (1997) 65:313–44. doi:10.1007/BF00053511

CrossRef Full Text | Google Scholar

14. Pohánka V. Optimum Expression for Computation of the Gravity Field of a Polyhedral Body with Linearly Increasing Density. Geophys Prospecting (1998) 46:391–404. doi:10.1046/j.1365-2478.1998.960335.x

CrossRef Full Text | Google Scholar

15. Hamayun , Prutkin I, Tenzer R. The Optimum Expression for the Gravitational Potential of Polyhedral Bodies Having a Linearly Varying Density Distribution. J Geod (2009) 83:1163–70. doi:10.1007/s00190-009-0334-1 Available at: https://link.springer.com/article/10.1007%2Fs00190-009-0334-1

CrossRef Full Text | Google Scholar

16. D’Urso MG. Gravity Effects of Polyhedral Bodies with Linearly Varying Density. Celest Mech Dyn Astr (2014) 120:349–72. doi:10.1007/s10569-014-9578-z

CrossRef Full Text | Google Scholar

17. Rao DB, Prakash MJ, Babu NR. 3D and 21/2 D Modelling of Gravity Anomalies with Variable Density Contrast. Geophys Prospect (1990) 38:411–22. doi:10.1111/j.1365-2478.1990.tb01854.x

CrossRef Full Text | Google Scholar

18. Gallardo‐Delgado LA, Pérez‐Flores MA, Gómez‐Treviño E. A Versatile Algorithm for Joint 3D Inversion of Gravity and Magnetic Data. Geophysics (2003) 68:949–59. doi:10.1190/1.1581067

CrossRef Full Text | Google Scholar

19. Gokula AP, Sastry RG. Gravitational Attraction of a Vertical Pyramid Model of Flat Top-And-Bottom with Depth-wise Parabolic Density Variation. J Earth Syst Sci (2015) 124:1735–44. doi:10.1007/s12040-015-0633-2

CrossRef Full Text | Google Scholar

20. Ren Z, Chen C, Pan K, Kalscheuer T, Maurer H, Tang J. Gravity Anomalies of Arbitrary 3D Polyhedral Bodies with Horizontal and Vertical Mass Contrasts. Surv Geophys (2017) 38:479–502. doi:10.1007/s10712-016-9395-x

CrossRef Full Text | Google Scholar

21. D’Urso MG, Trotta S. Gravity Anomaly of Polyhedral Bodies Having a Polynomial Density Contrast. Surv Geophys (2017) 38:781–832. doi:10.1007/s10712-017-9411-9

CrossRef Full Text | Google Scholar

22. Ren Z, Zhong Y, Chen C, Tang J, Kalscheuer T, Maurer H, et al. Gravity Gradient Tensor of Arbitrary 3D Polyhedral Bodies with up to Third-Order Polynomial Horizontal and Vertical Mass Contrasts. Surv Geophys (2018) 39:901–35. doi:10.1007/s10712-018-9467-1

CrossRef Full Text | Google Scholar

23. Pearson GS, Faux DA. Analytical Solutions for Strain in Pyramidal Quantum Dots. J Appl Phys (2000) 88:730–6. doi:10.1063/1.373729

CrossRef Full Text | Google Scholar

24. Glas F. Elastic Relaxation of Truncated Pyramidal Quantum Dots and Quantum Wires in a Half Space: An Analytical Calculation. J Appl Phys (2001) 90:3232–41. doi:10.1063/1.1394158

CrossRef Full Text | Google Scholar

25. Kuvshinov BN. Elastic and Piezoelectric fields Due to Polyhedral Inclusions. Int J Sol Structures (2008) 45:1352–84. doi:10.1016/j.ijsolstr.2007.09.024

CrossRef Full Text | Google Scholar

26. Nenashev AV, Dvurechenskii AV. Strain Distribution in Quantum Dot of Arbitrary Polyhedral Shape: Analytical Solution. J Appl Phys (2010) 107:064322. doi:10.1063/1.3357302

CrossRef Full Text | Google Scholar

27. Chen C, Ren Z, Pan K, Tang J, Kalscheuer T, Maurer H, et al. Exact Solutions of the Vertical Gravitational Anomaly for a Polyhedral Prism with Vertical Polynomial Density Contrast of Arbitrary Orders. Geophys J Int (2018) 214:2115–32. doi:10.1093/gji/ggy250

CrossRef Full Text | Google Scholar

28. Ren Z, Chen C, Zhong Y, Chen H, Kalscheuer T, Maurer H, et al. Recursive Analytical Formulae of Gravitational fields and Gradient Tensors for Polyhedral Bodies with Polynomial Density Contrasts of Arbitrary Non-negative Integer Orders. Surv Geophys (2020) 41:695–722. doi:10.1007/s10712-020-09587-4

CrossRef Full Text | Google Scholar

29. D’Urso MG. The Gravity Anomaly of a 2d Polygonal Body Having Density Contrast Given by Polynomial Functions. Surv Geophys (2015) 36:391–425. doi:10.1007/s10712-015-9317-3

CrossRef Full Text | Google Scholar

30. Van Oosterom A, Strackee J. The Solid Angle of a Plane triangle. IEEE Trans Biomed Eng (1983) BME-30:125–6. doi:10.1109/TBME.1983.325207

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Nenashev AV, Dvurechenskii AV. Strain Distribution in Quantum Dot of Arbitrary Polyhedral Shape: Analytical Solution in Closed Form (2007). Available at: https://arxiv.org/pdf/0707.2183.pdf

Google Scholar

32. Cordell L. Gravity Analysis Using an Exponential Density‐depth Function-san Jacinto Graben, California. Geophysics (1973) 38:684–90. doi:10.1190/1.1440367

CrossRef Full Text | Google Scholar

33. Chakravarthi V, Pramod Kumar M, Ramamma B, Rajeswara Sastry S. Automatic Gravity Modeling of Sedimentary Basins by Means of Polygonal Source Geometry and Exponential Density Contrast Variation: Two Space Domain Based Algorithms. J Appl Geophys (2016) 124:54–61. doi:10.1016/j.jappgeo.2015.11.007

CrossRef Full Text | Google Scholar

34. Mallesh K, Chakravarthi V, Ramamma B. 3D Gravity Analysis in the Spatial Domain: Model Simulation by Multiple Polygonal Cross-Sections Coupled with Exponential Density Contrast. Pure Appl Geophys (2019) 176:2497–511. doi:10.1007/s00024-019-02103-9

CrossRef Full Text | Google Scholar