Issue 
A&A
Volume 541, May 2012



Article Number  A130  
Number of page(s)  6  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201118443  
Published online  14 May 2012 
A substitute for the singular Green kernel in the Newtonian potential of celestial bodies
^{1} Univ. Bordeaux, LAB, UMR 5804, 33270 Floirac, France
^{2} CNRS, LAB, UMR 5804, 33270 Floirac, France
email: jeanmarc.hure@obs.ubordeaux1.fr
^{3} ELSA, Physikalisches Institut der Universität Bonn, Nussallee 12, 53115 Bonn, Germany
email: dieckman@physik.unibonn.de
Received: 13 November 2011
Accepted: 1 March 2012
The “point mass singularity” inherent in Newton’s law for gravitation represents a major difficulty in accurately determining the potential and forces inside continuous bodies. Here we report a simple and efficient analytical method to bypass the singular Green kernel 1/r − r′ inside the source without altering the nature of the interaction. We build an equivalent kernel made up of a “cool kernel”, which is fully regular (and contains the longrange − GM/r asymptotic behavior), and the gradient of a “hyperkernel”, which is also regular. Compared to the initial kernel, these two components are easily integrated over the source volume using standard numerical techniques. The demonstration is presented for threedimensional distributions in cylindrical coordinates, which are wellsuited to describing rotating bodies (stars, discs, asteroids, etc.) as commonly found in the Universe. An example of implementation is given. The case of axial symmetry is treated in detail, and the accuracy is checked by considering an exact potential/surface density pair corresponding to a flat circular disc. This framework provides new tools to keep or even improve the physical realism of models and simulations of selfgravitating systems, and represents, for some of them, a conclusive alternative to softened gravity.
Key words: gravitation / methods: analytical / methods: numerical
© ESO, 2012
1. Introduction
As a direct consequence of Newton’s law for gravitation (Newton 1760; Kellogg 1929), the potential of any continuous distribution of matter inside a volume at a point P(r) of space is given by (1)where ρ(r′) is the mass density at , and d^{3}τ is the elementary volume^{1}. In general, this is a threedimensional (3D) converging integral. The presence of the Green kernel 1/r − r′ is known to represent a difficulty in calculating ψ everywhere inside and very close to since this function diverges as r → r′. This singularity is classically avoided by converting the Green function into an infinite series (e.g. Kellogg 1929; Cohl & Tohline 1999). Although exact, series expansions suffer from a low convergence rate since these are alternating series; besides, the number of integrals increases linearly with the number of terms considered up to the truncation order. These problems collectively constitute a real practical difficulty (Clement 1974; Stone & Norman 1992). The proper treatment of the singularity is the subject of a longstanding challenge. Shifting the Pgrid and the P′grid relative to each other or raising the numerical resolution around the singularity are the most natural techniques (Stemwedel et al. 1990), but these are of limited efficiency. When possible, the separate treatment of the asymptotic form of the singularity gives very good results (e.g. Ansorg et al. 2003; Huré 2005), although this approach renders the global treatment somewhat complex. The difficulty can also be tackled by introducing a “softening length” (Hockney & Eastwood 1988; Adams et al. 1989). This recipe – widespread in disc simulations – must however be seen as nothing but a crude approximation that cannot be used for accurate modeling under a certain scale. Whatever the prescription for the softening length, which is generally linked to the resolution or smallest physical length scale (e.g. Huré & Pierens 2009), the inferred force field is globally weaker than in the Newtonian case, and the evolution and stability of gaseous systems is inevitably impacted in a nontrivial manner (Romeo 1998; SommerLarsen et al. 1998; Adams et al. 1989; ElZant 1998). The Poisson equation of course provides another other way to derive ψ numerically (Kellogg 1929; Durand 1953). This approach requires accurate boundary or interior/matching conditions only accessible through Eq. (1), and complex geometries are not always easy to manage (Grandclément et al. 2001).
In this paper, we present a new means to evaluate Eq. (1) that avoids the singularity, and, at the same time, properly accounts for it. This is achieved by replacing the singular kernel by an equivalent and regular, twoterm form, one term being the gradient of a new scalar potential^{2}. This is the aim of Sects. 2 and 3. Kernel equivalence is fundamental to preserving the Newtonian character of the interaction on all scales, and for any separation (in particular at longrange). This reformulation is designed to be efficient within sources, and is not expected to surpass usual methods outside sources. From a practical point of view, the potential easily becomes accessible as diverging kernels have disappeared from the volume integrals.
Although not specific to a given system of coordinates, the calculus is developed in cylindrical coordinates for which the equivalent kernel takes a nominal form, in particular under axial symmetry (this equivalent kernel is probably not unique; see Sect. 4). It can be applied as it is to all rotating gaseous/solid bodies (stars, discs, planets, asteroids, etc.) in either steady state or not, and for various applications (Dermott 1979; Hachisu 1986b; Baruteau & Masset 2008). A basic, sixstep algorithm is reported in Sect. 5, together with a numerical experiment using the most simple quadrature and differentiation rules. There is no special assumption about the distribution of matter in space (density field and geometry or shape), making the method general, and transposable to domains of physics other than gravitation. Some interesting perspectives are listed in the last section.
Fig. 1 Typical configuration for the gravitating, celestial body, and associated notations. See note 3 for the definition of points P′′ and Q. 
2. Splitting of the Green kernel
Fig. 2 Dimensionless Green kernel (left) and cool kernel (right) versus a/R for ζ = 0 and labeled on the curves. The Green kernel diverges hyperbolically when r → r′ (which corresponds to a → R and here), in contrast to the cool kernel, which remains bounded (and even vanishes at the singularity). 
We consider a volume of space continuously filled with matter, as depicted in Fig. 1. Using cylindrical coordinates, with P′(a,θ′,z) referring to source points, and P(R,θ,Z) to space points, the above integral for the Newtonian potential becomes (2)where d^{3}τ = adadθ′dz is the elementary volume, is the relative altitude, and (5)Although Δ → 0 inside , the potential is generally a finite quantity (i.e. integration is a regularizing process; e.g. Kellogg 1929; Durand 1953). We now set (6)This quantity is finite everywhere, and nonzero except on the polar axis (see below). Assuming R > 0, we have (7)and so (8)where we have defined (9)The Green kernel is then split into two terms. The term is always regular^{3}; we call this term the cool kernel in the following. Figure 2 displays R/Δ and R/Δ_{⋆} versus a/R around the singularity. We clearly see that the amplitude of Δ_{⋆} is bounded, in contrast to Δ. The second term in Eq. (8) is still singular when Δ → 0, but its integration over the material volume is expected to produce a regular field. The idea is to generate this singular kernel from the gradient of a regular function κ (hereafter called hyperkernel), and to integrate it over the volume . As the two spaces (R,θ,Z) and (a,θ′,z) are independent, the derivative may be drawn before the integral. This reasoning can be summarized as where the gradient ∇ is to be taken with respect to one of the three variables R, θ, or Z. There are then three possible hyperkernels.
The existence of the hyperkernel is not guaranteed a priori, and it is of interest only if it is available in a closed form. The investigation indeed shows that the nominal form is obtained by considering the vertical gradient (i.e. ∇ ≡ ∂_{Z}). This may be due to the special role that the Zaxis plays a in cylindrical coordinates. We therefore do not discuss in detail any of the other two options, although these might be useful in certain circumstances.
3. The singular term as the vertical gradient of a hyperkernel
To get the hyperkernel, we consider the integration of the singular term in Eq. (8) with respect to Z. By using the intermediate variable t = ζ/Δ, we find after some algebra (11)where (12)with 0 ≤ m ≤ 1. We finally get (13)We could obviously add to κ any function of a, R, θ and θ′, but this is not necessary here as we take its Zgradient. The Green kernel is then given by the equivalent form (14)It is necessary to verify that κ is regular. This is straightforward since Δ > ζ as soon as R > 0. In other words, if m → 1 and , then ζ/Δ → ± 1.
Fig. 3 Main steps in the computation of the potential from Eqs. (25) and (27) in a typical case: a) the integral of the cool kernel, b) the integral of the hyperkernel, c) its vertical gradient, and d) the potential as the sum of maps a) and c). Here, the body is an axially symmetric torus with square cross section (boundary indicated with a white line). As clearly visible in graphs a)–c), the treatment differs slightly on the polar axis (first column of pixels). See the text for the numerical setup. 
If we now multiply Eq. (8) by ρ(a,θ′,z)d^{3}τ – which does not depend on Z – and integrate over the material volume , we find that (15)where the partial derivative now operates on the integral. Up to a factor − G, this expression is precisely the potential defined by Eq. (1), and it is both exact and general. It depends on neither the body’s shape nor on the distribution of its mass density. It applies not only to volume distributions, but also to surface distributions (see Sect. 6) and linear distributions. The first integral in Eq. (15) is then the cool potential associated with the cool kernel, and the second term is the vertical gradient of a hyperpotential.
On the polar axis (i.e. R = 0) Δ = δ, and so Eq. (8) does not help us to treat the singularity when ζ = 0 and a = 0. In this case, we have (16)provided that a > 0. The potential can then be written (17)where here there is no cool kernel. Besides, we see that lim_{a → 0} aκ = 0.
There is no continuity between the two different expressions for the cool kernel, the one valid at R = 0 and the other valid at R → 0 (this is also true for the hyperkernel). This is no problem as long as we do not have to consider the radial gradient of κ (see below the numerical experiment).
Finally, we note that, as we work with an equivalent form of the Green kernel, the potential found from Eq. (15) automatically has the right asymptotic property, and varies like M/r sufficiently far away from the body. At large relative distance (i.e. R ≫ a and Z ≫ z), we have . At the lowest order, one finds that (18)which means that the longrange behavior is exclusively contained in the cool kernel.
4. The case of axiallysymmetric bodies
Axiallysymmetric bodies constitute an important class of astrophysical objects (Chandrasekhar 1973; Hachisu 1986a). Interestingly enough, in problems where ∂_{θ′}ρ = 0, we can rewrite the above expressions in a more compact form in terms of elliptic integrals. The first integral in the righthandside of Eq. (15) becomes^{4}(19)where (20)is the complete elliptic integral of the second kind, is the modulus (with k ∈ [0,1] ), and the double integration runs over the meridional cross section of the body. The integration over the polar angle θ′ of the hyperkernel gives: (21)where (22)is the incomplete elliptic integral of the first kind, and (23)is the incomplete elliptic integral of the third kind (m is the parameter and ). Over the whole circle (i.e. ), this yields the axially symmetric potential^{5}: (25)where H is defined for convenience by (26) with and . The presence of the Hfunction is actually expected here since ∂_{z}Δ and ∂_{Z}Δ are linked (Trova et al. 2012).
On the polar axis, we get (27) where κ is, in this case, given by Eq. (16). We note that, in this axially symmetric case, the potential could be determined through Eq. (19) (i.e. by using the cool kernel only, and no hyperkernel), but the integrand still contains a hyperbolic divergence as a → 0 and ζ → 0 which is not easy to manage. This is why it seems much better to consider Eq. (27), as the logarithmic divergence of the hyperkernel (i.e. the asinh term) is cancelled out by the elementary volume when a is close to 0 (i.e. lim_{a → 0} aκ = 0).
5. An example of implementation
We now present a first numerical experiment to briefly describe the main steps of the method, and demonstrate its simple implementation. For a full 3D body, the main steps^{6} of any numerical estimate can be summarized as

1.
discretize the material source on a 3Dgrid, in the form of quadruplets ;

2.
select a point P(R,θ,Z) in space;

3.
for each point of the source, compute the cool kernel, and the hyperkernel κ from Eqs. (9) and (13) or Eq. (16);
 4.

5.
determine the vertical gradient of the hyperpotential;

6.
add this derivative to the cool potential, and multiply by − G, and reiterate steps 2 to 6 to generate a potential map. The components of the associated gravitational force are deduced as usual from the three gradients of ψ.
We have considered a homogeneous, axially symmetric torus with a square meridional cross section with . The mass density (ρ = 1 inside and 0 outside) is defined on a regular mesh, with N nodes in each direction. The computational grid also consists of a square box (R,Z) ∈ [0,2] × [−1,1] with L nodes per direction and regular spacing, therefore encompassing the body. The double integrals in Eqs. (25) and (27) are computed at each node of the computational grid through the twopoint trapezoidal rule. The partial derivative of the hyperpotential is estimated through secondorder finite differences. These basic schemes are very easy to implement and the above sixstep procedure contains no pitfall. We take N = L = 31 here. Figures 3a–d display, respectively, the integral of the cool kernel, the integral of the hyperkernel κ, its vertical gradient, and the total potential, obtained by adding the first and third maps. The boundary of the toroidal body is superimposed on these maps. We note that the vertical gradient of the hyperpotential makes the geometrical cross section rise above the background, and filters the curvature effects (which are enhanced by the integral of the cool kernel).
The computing time is typical of integral methods. Under the conditions of the present example, the integration of the cool kernel and hyperkernel (step 4 above) requires 2N^{2} elementary operations per point of space. To get the potential at the nodes of a L × L square grid, we then find 2N^{2}L^{2} (the time needed to determine the vertical gradients is negligible in comparison). This is therefore much smaller than that is usually obtained from an expanded Green function by a factor equal to the number of terms up to the truncation order (which can be as high as a few hundred; see e.g. Hachisu 1986a; Stone & Norman 1992).
Fig. 4 Same legend as for Fig. 3 but for the flat, axially symmetric disc with surface density Σ ∝ (1 − a^{2})^{3/2}, inner edge a_{in} = 0 and outer edge a_{out} = 1. The disc is indicated with a white line. The exact potential ψ_{e} (left) is derived from the formula of Schulz (2009). The associated error index (right) is determined once ψ is computed from Eqs. (29) and (30). The mesh size is corresponding to N = 31 radial points. 
6. Checking the accuracy
The numerical accuracy is sensitive to various ingredients, such as the quadrature and differentiation schemes. To a lesser extent, it also depends on the mass density distribution ρ and equation of the boundary , which may generate additional difficulties in the calculations (interpolation of data points, infinite derivatives at edges, etc.). We present a second numerical test illustrating the accuracy of the method by considering an exact potential/density pair. Not many configurations correspond to finite mass and finite size systems (e.g. Binney & Tremaine 1987). When matter is gathered in a plane (i.e. a flat disc), Eq. (15) is directly transposable by setting ρ(a,θ′,z) = Σ(a,θ′)δ(z) and integrating over z. Under axial symmetry, we respectively get^{7} from Eqs. (25) and (27) (29)for R > 0, and (30)where a_{in} and a_{out} denote the discs inner and outer edges, and (31)By setting a_{in} = 0 and for a ∈ [0,a_{out}] , one gets one of the three cases analyzed by Schulz (2009). For such a distribution, the associated potential, ψ_{e}, is known exactly in a closedform for any point of space. Figure 4 gives ψ_{e} as well as the error index ϵ = log _{10}1 − ψ/ψ_{e} where ψ is determined from Eqs. (29) and (30) in the same conditions as above (we set a_{out} = 1). We see that the potential outside and especially inside the disc is wellreproduced. The relative error, on the order of 10^{3}, agrees with the secondorder of the schemes at the actual mesh size of . The accuracy can be tuned by changing the quadrature and differentation schemes.
7. Concluding remarks
We have reformulated the Green kernel appearing in potential problems to circumvent the singularity and, at the same time, properly account for it. As a consequence, the gravitational potential of any celestial body, regardless of its shape and matter density distribution, becomes directly accessible through two “classical” volume integrals, followed by a partial derivative. The method is applicable to 3D, fully inhomogeneous systems, as well as to surface and line distributions. It is especially efficient inside distributions where most approaches exhibit a real practical complexity, converge very slowly, or produce spurious errors. The presence of regular kernels ensures that the method is stable and easy to implement. This should encourage modellers to abandon various integration techniques that do not “faithfully” reproduce the Newtonian character of the potential and forces. In the context of discs for instance, this method appears to be a real alternative to softened gravity, which remains a free parameter, nonNewtonian theory. As stressed, it is probably possible to determine other cool kernel/hyperkernel pairs (for instance, by considering the hyperkernel as a radial/angular gradient), but the one presented in the body of this paper seems the simplest one. It is in particular wellsuited to axially symmetric configurations.
This study needs to be continued in several respects, including the analysis of the mathematical properties of the cool kernel and hyperkernel and their physical meanings, as well as the derivation of the equivalent kernel in other systems of coordinates (e.g. Cartesian and spherical coordinate systems). In addition, it would be interesting to expand the two kernels in Eq. (15) in series, and compare their properties with the expansion of the Green function in Legendre polynomials, inside as well as outside the body. The cool kernel/hyperkernel pair is also interesting as a new starting point to generating various kinds of approximations. Apart from the astrophysical context where there are so many applications about gravitation, this technique is also transposable to other kinds of problems involving improper integrals. This can be, for instance in electromagnetism, the determination of the potential vector A and associated magnetic field induced by current densities (Jackson 1998; Cohl & Tohline 1999). These points will be touched on in forthcoming papers.
Some aspects of the theory of tensor potentials are developped by Chandrasekhar (1973) in the context of rotating, selfgravitating fluids.
It is the inverse of a distance, and its value can be interpreted geometrically by noting that (see Fig. 1) (10) where Q(R,θ′ + π,Z) is a point of space, diametrically opposite to P′, and P′′(a,θ + π,z) is a point, diametrically opposite to P, that belongs to a fictitious source.
If we perform the Zderivative and rearrange terms, we recover the wellknown expression (Durand 1953), namely (24) whose kernel is logarithmically singular.
These two formulae can be compared with the classical expression (Durand 1953; Binney & Tremaine 1987): (28) which is singular inside the source.
Acknowledgments
It is a pleasure to thank J. Braine, M. Gazeau, F. Hersant, and A. Trova as well as the second referee for comments and suggestions about infinitely flat systems. J.M. Huré is grateful to the CNU, Sect. 34, for supporting a sixmonths fulltime research project through CRCT2011 funding delivered by the MESR.
References
 Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959 [NASA ADS] [CrossRef] [Google Scholar]
 Ansorg, M., Kleinwächter, A., & Meinel, R. 2003, MNRAS, 339, 515 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., & Masset, F. 2008, ApJ, 678, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton, NJ: Princeton University Press, 747) [Google Scholar]
 Chandrasekhar, S. 1973, Ellipsoidal figures of equilibrium [Google Scholar]
 Clement, M. J. 1974, ApJ, 194, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Cohl, H. S., & Tohline, J. E. 1999, ApJ, 527, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Dermott, S. F. 1979, Icarus, 37, 575 [NASA ADS] [CrossRef] [Google Scholar]
 Durand, E. 1953, Electrostatique, Vol. I, Les distributions (Masson) [Google Scholar]
 ElZant, A. A. 1998, A&A, 331, 782 [NASA ADS] [Google Scholar]
 Grandclément, P., Bonazzola, S., Gourgoulhon, E., & Marck, J.A. 2001, J. Comp. Phys., 170, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Hachisu, I. 1986a, ApJS, 61, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Hachisu, I. 1986b, ApJS, 62, 461 [NASA ADS] [CrossRef] [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation Using Particles, ed. R. W. Hockney, & J. W. Eastwood [Google Scholar]
 Huré, J.M. 2005, A&A, 434, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Huré, J., & Pierens, A. 2009, A&A, 507, 573 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jackson, J. D. 1998, Classical Electrodynamics, 3rd edn. [Google Scholar]
 Kellogg, O. D. 1929, Foundations of Potential Theory (NewYork: Frederick Ungar Publishing Company) [Google Scholar]
 Newton, I. 1760 [Google Scholar]
 Romeo, A. B. 1998, A&A, 335, 922 [NASA ADS] [Google Scholar]
 Schulz, E. 2009, ApJ, 693, 1310 [NASA ADS] [CrossRef] [Google Scholar]
 SommerLarsen, J., Vedel, H., & Hellsten, U. 1998, MNRAS, 294, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Stemwedel, S. W., Yuan, C., & Cassen, P. 1990, ApJ, 351, 206 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [Google Scholar]
 Trova, A., Huré, J. M., & Hersant, F. 2012, MNRAS, submitted [Google Scholar]
All Figures
Fig. 1 Typical configuration for the gravitating, celestial body, and associated notations. See note 3 for the definition of points P′′ and Q. 

In the text 
Fig. 2 Dimensionless Green kernel (left) and cool kernel (right) versus a/R for ζ = 0 and labeled on the curves. The Green kernel diverges hyperbolically when r → r′ (which corresponds to a → R and here), in contrast to the cool kernel, which remains bounded (and even vanishes at the singularity). 

In the text 
Fig. 3 Main steps in the computation of the potential from Eqs. (25) and (27) in a typical case: a) the integral of the cool kernel, b) the integral of the hyperkernel, c) its vertical gradient, and d) the potential as the sum of maps a) and c). Here, the body is an axially symmetric torus with square cross section (boundary indicated with a white line). As clearly visible in graphs a)–c), the treatment differs slightly on the polar axis (first column of pixels). See the text for the numerical setup. 

In the text 
Fig. 4 Same legend as for Fig. 3 but for the flat, axially symmetric disc with surface density Σ ∝ (1 − a^{2})^{3/2}, inner edge a_{in} = 0 and outer edge a_{out} = 1. The disc is indicated with a white line. The exact potential ψ_{e} (left) is derived from the formula of Schulz (2009). The associated error index (right) is determined once ψ is computed from Eqs. (29) and (30). The mesh size is corresponding to N = 31 radial points. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.