# A substitute for the singular Green kernel in the Newtonian potential of celestial bodies

J.-M. Huré, Andreas Dieckmann

► To cite this version:

J.-M. Huré, Andreas Dieckmann. A substitute for the singular Green kernel in the Newtonian potential of celestial bodies. Astronomy & Astrophysics - A&A, 2012, 541, pp.A130. 10.1051/0004-6361/201118443. hal-00685411

HAL Id: hal-00685411

https://hal.science/hal-00685411v1

Submitted on 4 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A substitute for the singular Green kernel in the Newtonian potential of celestial bodies

J.-M. Huré<fnref n="1" />,<fnref n="2" /> and A. Dieckmann<fnref n="3" />

Received 13 November 2011 / Accepted 1 March 2012

## ABSTRACT

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 long-range $-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 three-dimensional distributions in cylindrical coordinates, which are well-suited 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 self-gravitating systems, and represents, for some of them, a conclusive alternative to softened gravity.

**Key words.** gravitation – methods: analytical – methods: numerical

## 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 $V$ at a point $P(r)$ of space is given by

$$ \psi(r) = -G \iiint_V \frac{\rho(r')d^3\tau}{|r - r'|}, \quad (1) $$

where $\rho(r')$ is the mass density at $P'(r') \in V$, and $d^3\tau$ is the elementary volume<fnref n="1" />. In general, this is a three-dimensional (3D) converging integral. The presence of the Green kernel $1/|r - r'|$ is known to represent a difficulty in calculating $\psi$ everywhere inside and very close to $V$ since this function diverges as $r \to 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 P-grid 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 non-trivial manner (Romeo 1998; Sommer-Larsen et al. 1998; Adams et al. 1989; El-Zant 1998). The Poisson equation of course provides another way to derive $\psi$ 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, two-term form, one term being the gradient of a new scalar potential<fnref n="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 long-range). This reformulation is designed to

<footnote n="1">¹ Univ. Bordeaux, LAB, UMR 5804, 33270 Floirac, France</footnote>

<footnote n="2">² CNRS, LAB, UMR 5804, 33270 Floirac, France
e-mail: jean-marc.hure@obs.u-bordeaux1.fr</footnote>

<footnote n="3">³ ELSA, Physikalisches Institut der Universität Bonn, Nussallee 12, 53115 Bonn, Germany
e-mail: dieckman@physik.uni-bonn.de</footnote>

<footnote n="1">¹ This form holds in electrostatics, where $\rho$ is the density of electric charges, and the constant is $\frac{1}{4\pi\epsilon_0}$ (instead of $-G$).</footnote>

<footnote n="2">² Some aspects of the theory of tensor potentials are developed by Chandrasekhar (1973) in the context of rotating, self-gravitating fluids.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Typical configuration for the gravitating, celestial body, and associated notations. See note 3 for the definition of points $P''$ and $Q$.</figcaption></figure>

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, six-step 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.

## 2. Splitting of the Green kernel

We consider a volume of space $V$ continuously filled with matter, as depicted in Fig. 1. Using cylindrical coordinates, with $P'(a, \theta', z)$ referring to source points, and $P(R, \theta, Z)$ to space points, the above integral for the Newtonian potential becomes

$$ \psi(R, \theta, Z) = -G \iint\int_V \frac{1}{\Delta} \rho(a, \theta', z) d^3\tau, \quad (2) $$

where $d^3\tau = a da d\theta' dz$ is the elementary volume,

$$ \begin{aligned} \Delta^2 &= |\mathbf{r} - \mathbf{r}'|^2 \\ &= (a + R)^2 + \zeta^2 - 4aR \sin^2 \phi, \end{aligned} \quad (3) $$

$$ \zeta = Z - z, \quad (4) $$

is the relative altitude, and

$$ 2\phi = \pi - (\theta - \theta'). \quad (5) $$

Although $\Delta \to 0$ inside $V$, the potential is generally a finite quantity (i.e. integration is a regularizing process; e.g. Kellogg 1929; Durand 1953). We now set

$$ \delta = \sqrt{(a + R)^2 + \zeta^2} \ge R. \quad (6) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Dimensionless Green kernel (left) and cool kernel (right) versus $a/R$ for $\zeta = 0$ and $\phi = \{0, \frac{\pi}{4}, \frac{2\pi}{3}, \frac{\pi}{2}\}$ labeled on the curves. The Green kernel diverges hyperbolically when $\mathbf{r} \to \mathbf{r}'$ (which corresponds to $a \to R$ and $\phi \to \frac{\pi}{2}$ here), in contrast to the cool kernel, which remains bounded (and even vanishes at the singularity).</figcaption></figure>

This quantity is finite everywhere, and non-zero except on the polar axis (see below). Assuming $R > 0$, we have

$$ \begin{aligned} \frac{\Delta}{\delta^2} &= \frac{1}{\Delta} \frac{\Delta^2}{\delta^2} \\ &= \frac{1}{\Delta} \left( 1 - \frac{4aR \sin^2 \phi}{\delta^2} \right), \end{aligned} \quad (7) $$

and so

$$ \frac{1}{\Delta} = \frac{1}{\Delta_*} + 4aR \sin^2 \phi \times \frac{1}{\Delta \delta^2}, \quad (8) $$

where we have defined

$$ \frac{1}{\Delta_*} = \frac{\Delta}{\delta^2}. \quad (9) $$

The Green kernel is then split into two terms. The term $\frac{1}{\Delta_*}$ is always regular$^3$; we call this term the cool kernel in the following. Figure 2 displays $R/\Delta$ and $R/\Delta_*$ versus $a/R$ around the singularity. We clearly see that the amplitude of $\Delta_*$ is bounded, in contrast to $\Delta$. The second term in Eq. (8) is still singular when $\Delta \to 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 $\kappa$ (hereafter called hyperkernel), and to integrate it over the volume $V$. As the two spaces $(R, \theta, Z)$ and $(a, \theta', z)$ are independent, the derivative may be drawn before the integral. This reasoning can be summarized as

singular kernel $\equiv \nabla$ regular hyperkernel

↓

$$ \iiint_V \text{sing. kernel } d^3\tau = \nabla \iint_V \text{reg. hyperkernel } d^3\tau, $$

$$ \frac{\Delta}{\delta^2} = \frac{|\mathbf{r} - \mathbf{r}'|}{Q P'^2} = \frac{|\mathbf{r} - \mathbf{r}'|}{P P''^2} = \frac{|\mathbf{r} - \mathbf{r}'|}{d^2} < \infty, \quad (10) $$

<footnote>$^3$ It is the inverse of a distance, and its value can be interpreted geometrically by noting that (see Fig. 1)</footnote>

<footnote>where $Q(R, \theta' + \pi, Z)$ is a point of space, diametrically opposite to $P'$, and $P''(a, \theta + \pi, z)$ is a point, diametrically opposite to $P$, that belongs to a fictitious source.</footnote>
where the gradient $\nabla$ is to be taken with respect to one of the three variables $R$, $\theta$, 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. $\nabla \equiv \partial_z$). This may be due to the special role that the $Z$-axis 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 = \zeta/\Delta$, we find after some algebra

$$ \int \frac{dZ}{\delta^2 \Delta} = \frac{1}{(a + R)^2} \int \frac{dt}{1 - m^2 \sin^2 \phi \times t^2}, \quad (11) $$

where

$$ m = \frac{2 \sqrt{aR}}{a + R}, \quad (12) $$

with $0 \le m \le 1$. We finally get

$$ \begin{aligned} \kappa &\equiv 4aR \sin^2 \phi \int \frac{dZ}{\delta^2 \Delta} \\ &= m \sin \phi \operatorname{atanh} \left( \frac{\zeta m \sin \phi}{\Delta} \right). \end{aligned} \quad (13) $$

We could obviously add to $\kappa$ any function of $a$, $R$, $\theta$ and $\theta'$, but this is not necessary here as we take its $Z$-gradient. The Green kernel is then given by the equivalent form

$$ \frac{1}{\Delta} = \frac{1}{\Delta_*} + \partial_Z \kappa. \quad (14) $$

It is necessary to verify that $\kappa$ is regular. This is straightforward since $\Delta > |\zeta|$ as soon as $R > 0$. In other words, if $m \to 1$ and $\phi \to \pi/2$, then $\zeta/\Delta \to \pm 1$.

If we now multiply Eq. (8) by $\rho(a, \theta', z)d^3\tau$ – which does not depend on $Z$ – and integrate over the material volume $V$, we find that

$$ \iiint_V \frac{1}{\Delta} \rho d^3 \tau = \iiint_V \frac{1}{\Delta_*} \rho d^3 \tau + \partial_Z \iiint_V \rho \kappa d^3 \tau, \quad (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$) $\Delta = \delta$, and so Eq. (8) does not help us to treat the singularity when $\zeta = 0$ and $a = 0$. In this case, we have

$$ \begin{aligned} \kappa &\equiv \int \frac{1}{\Delta} dZ \\ &= \operatorname{asinh} \frac{\zeta}{a}, \end{aligned} \quad (16) $$

provided that $a > 0$. The potential can then be written

$$ \psi(0, Z) = -G \partial_Z \iiint_V \rho \kappa d^3 \tau, \quad (17) $$

where here there is no cool kernel. Besides, we see that $\lim_{a \to 0} ak = 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 \to 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 $\kappa$ (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 \gg a$ and $Z \gg z$), we have $\delta \to r = \sqrt{R^2 + Z^2}$. At the lowest order, one finds that

$$ \iiint_V \frac{1}{\Delta_*} \rho d^3 \tau \approx \frac{M}{r}, \quad (18) $$

which means that the long-range behavior is exclusively contained in the cool kernel.

### 4. The case of axially-symmetric bodies

Axially-symmetric bodies constitute an important class of astrophysical objects (Chandrasekhar 1973; Hachisu 1986a). Interestingly enough, in problems where $\partial_{\theta'} \rho = 0$, we can rewrite the above expressions in a more compact form in terms of elliptic integrals. The first integral in the right-hand-side of Eq. (15) becomes<fnref n="4" />

$$ \iiint_V \frac{1}{\Delta_*} \rho d^3 \tau = 2 \iint_S \rho \sqrt{\frac{a}{R}} k \mathbf{E}(k) d\mathbf{a} dz, \quad (19) $$

where

$$ \mathbf{E}(k) = \int_0^{\pi/2} \sqrt{1 - k^2 \sin^2 x} \, dx, \quad (20) $$

is the complete elliptic integral of the second kind, $k = 2\sqrt{aR}/\delta$ is the modulus (with $k \in [0, 1]$), and the double integration runs over the meridional cross section $S$ of the body. The integration over the polar angle $\theta'$ of the hyperkernel gives:

$$ \int_{\phi} \sin \phi \operatorname{atanh} \left( \frac{\zeta m \sin \phi}{\Delta} \right) d\phi = -\cos \phi \operatorname{atanh} \left( \frac{\zeta m \sin \phi}{\Delta} \right) + \frac{\zeta}{m\delta} \left[ F(\phi, k) - m'^2 \Pi(\phi, m, k) \right], \quad (21) $$

where

$$ F(\phi, k) = \int_0^\phi \frac{dx}{\sqrt{1 - k^2 \sin^2 x}}, \quad (22) $$

is the incomplete elliptic integral of the first kind, and

$$ \Pi(\phi, m, k) = \int_0^\phi \frac{dx}{(1 - m^2 \sin^2 x) \sqrt{1 - k^2 \sin^2 x}}, \quad (23) $$

<footnote n="4">⁴ A factor of two is due to $d\theta'/d\phi$, and another factor of two contained in the modulus $k$ comes from symmetry consideration (i.e. matter located at $\theta' \in [\theta, \theta + \pi]$) provides the same contribution as matter located at $\theta' \in [\theta, \theta - \pi]$).</footnote>
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

is the incomplete elliptic integral of the third kind ($m$ is the parameter and $m' = \sqrt{1-m^2}$). Over the whole circle (i.e. $\phi \in [0, \frac{\pi}{2}]$), this yields the axially symmetric potential<fnref n="5" />:

$$ \psi(R, Z) = -2G \iint_S \rho \sqrt{\frac{a}{R}} k \mathbf{E}(k) \mathrm{d}a \mathrm{d}z \\ - 2G \partial_Z \left[ \iint_S \rho \sqrt{\frac{a}{R}} \zeta \mathbf{H}(m, k) \mathrm{d}a \mathrm{d}z \right], \quad (25) $$

where $\mathbf{H}$ is defined for convenience by

$$ \mathbf{H}(m, k) = k \left[ \mathbf{K}(k) - m'^2 \mathbf{\Pi}(m, k) \right], \quad (26) $$

$$ \psi(R, Z) = -2G \iint_S \rho \sqrt{\frac{a}{R}} k \mathbf{K}(k) \mathrm{d}a \mathrm{d}z, \quad (24) $$

whose kernel is logarithmically singular.

with $\mathbf{\Pi}(m, k) = \Pi(\frac{\pi}{2}, m, k)$ and $\mathbf{K}(k) = F(\frac{\pi}{2}, k)$. The presence of the $\mathbf{H}$-function is actually expected here since $\partial_z\Delta$ and $\partial_z\Delta$ are linked (Trova et al. 2012).

On the polar axis, we get

$$ \psi(0, Z) = -2\pi G \partial_Z \iint_S \rho \kappa a \mathrm{d}a \mathrm{d}z, \quad (27) $$

where $\kappa$ 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 \to 0$ and $\zeta \to 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 \to 0} ak = 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

<footnote n="5">⁵ If we perform the Z-derivative and rearrange terms, we recover the well-known expression (Durand 1953), namely</footnote>
implementation. For a full 3D body, the main steps<fnref n="6" /> of any numerical estimate can be summarized as

1. discretize the material source $\mathcal{V}$ on a 3D-grid, in the form of quadruplets $\{P'(a_{i,j,k}, \theta'_{i,j,k}, z_{i,j,k}), \rho_{i,j,k}\}$;

2. select a point $P(R, \theta, Z)$ in space;

3. for each point $P'_{i,j,k}$ of the source, compute the cool kernel, and the hyperkernel $\kappa$ from Eqs. (9) and (13) or Eq. (16);

4. perform the volume integrals in Eq. (15) or (17);

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 $\psi$.

We have considered a homogeneous, axially symmetric torus with a square meridional cross section $S$ with $(a, z) \in [\frac{1}{2}, \frac{3}{2}] \times [-\frac{1}{2}, \frac{1}{2}]$. The mass density ($\rho = 1$ inside $S$ 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) \in [0, 2] \times [-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 two-point trapezoidal rule. The partial derivative of the hyperpotential is estimated through second-order finite differences. These basic schemes are very easy to implement and the above six-step 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 $\kappa$, 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 \times L$ square grid, we then find $2N^2L^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).

## 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 $\rho$ and equation of the boundary $\partial\mathcal{V}$, 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 $\rho(a, \theta', z) = \Sigma(a, \theta')\delta(z)$ and

integrating over $z$. Under axial symmetry, we respectively get<fnref n="7" /> from Eqs. (25) and (27)

$$ \psi(R, Z) = -2G \int_{a_{\text{in}}}^{a_{\text{out}}} \Sigma \sqrt{\frac{a}{R}} k_0 \mathbf{E}(k_0) \mathrm{d}a - 2G \partial_Z \left[ \int_{a_{\text{in}}}^{a_{\text{out}}} \Sigma \sqrt{\frac{a}{R}} Z \mathbf{H}(m, k_0) \mathrm{d}a \right], \quad (29) $$

for $R > 0$, and

$$ \psi(0, Z) = -2\pi G \partial_Z \int_{a_{\text{in}}}^{a_{\text{out}}} \Sigma \sinh\left(\frac{Z}{a}\right) \mathrm{d}a, \quad (30) $$

where $a_{\text{in}}$ and $a_{\text{out}}$ denote the discs inner and outer edges, and

$$ k_0^2 = \frac{4aR}{(a+R)^2 + Z^2}. \quad (31) $$

By setting $a_{\text{in}} = 0$ and $\Sigma \propto (1 - a^2/a_{\text{out}}^2)^{3/2}$ for $a \in [0, a_{\text{out}}]$, one gets one of the three cases analyzed by Schulz (2009). For such a distribution, the associated potential, $\psi_e$, is known exactly in a closed-form for any point of space. Figure 4 gives $\psi_e$ as well as the error index $\epsilon = \log_{10}|1 - \psi/\psi_e|$ where $\psi$ is determined from Eqs. (29) and (30) in the same conditions as above (we set $a_{\text{out}} = 1$). We see that the potential outside and especially inside the disc is well-reproduced. The relative error, on the order of $10^{-3}$, agrees with the second-order of the schemes at the actual mesh size of $\frac{a_{\text{out}}-a_{\text{in}}}{N-1} = \frac{1}{30}$. The accuracy can be tuned by changing the quadrature and differentiation 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, non-Newtonian 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 well-suited 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

$$ \psi(R, Z) = -2G \int_{a_{\text{in}}}^{a_{\text{out}}} \Sigma \sqrt{\frac{a}{R}} k_0 \mathbf{K}(k_0) \mathrm{d}a, \quad (28) $$

which is singular inside the source.

<footnote n="7">⁷ These two formulae can be compared with the classical expression (Durand 1953; Binney & Tremaine 1987):</footnote>

<footnote n="6">⁶ Step 1 should be executed once and for all, except for systems that evolve with time.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Same legend as for Fig. 3 but for the flat, axially symmetric disc with surface density $\Sigma \propto (1 - a^2)^{3/2}$, inner edge $a_{\text{in}} = 0$ and outer edge $a_{\text{out}} = 1$. The disc is indicated with a white line. The exact potential $\psi_c$ (left) is derived from the formula of Schulz (2009). The associated error index (right) is determined once $\psi$ is computed from Eqs. (29) and (30). The mesh size is $1/30$ corresponding to $N = 31$ radial points.</figcaption></figure>

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.

***Acknowledgements.*** 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 six-months full-time research project through CRCT-2011 funding delivered by the MESR.

## References

* Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959

* Ansorg, M., Kleinwächter, A., & Meinel, R. 2003, MNRAS, 339, 515

* Baruteau, C., & Masset, F. 2008, ApJ, 678, 483

* Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton, NJ: Princeton University Press, 747)

* Chandrasekhar, S. 1973, Ellipsoidal figures of equilibrium

* Clement, M. J. 1974, ApJ, 194, 709

* Cohl, H. S., & Tohline, J. E. 1999, ApJ, 527, 86

* Dermott, S. F. 1979, Icarus, 37, 575

* Durand, E. 1953, Electrostatique, Vol. I, Les distributions (Masson)

* El-Zant, A. A. 1998, A&A, 331, 782

* Grandclément, P., Bonazzola, S., Gourgoulhon, E., & Marck, J.-A. 2001, J. Comp. Phys., 170, 231

* Hachisu, I. 1986a, ApJS, 61, 479

* Hachisu, I. 1986b, ApJS, 62, 461

* Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation Using Particles, ed. R. W. Hockney, & J. W. Eastwood

* Huré, J.-M. 2005, A&A, 434, 1

* Huré, J., & Pierens, A. 2009, A&A, 507, 573

* Jackson, J. D. 1998, Classical Electrodynamics, 3rd edn.

* Kellogg, O. D. 1929, Foundations of Potential Theory (New-York: Frederick Ungar Publishing Company)

* Newton, I. 1760

* Romeo, A. B. 1998, A&A, 335, 922

* Schulz, E. 2009, ApJ, 693, 1310

* Sommer-Larsen, J., Vedel, H., & Hellsten, U. 1998, MNRAS, 294, 485

* Stemwedel, S. W., Yuan, C., & Cassen, P. 1990, ApJ, 351, 206

* Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753

* Trova, A., Huré, J. M., & Hersant, F. 2012, MNRAS, submitted