Issue 
A&A
Volume 554, June 2013



Article Number  A45  
Number of page(s)  9  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201321146  
Published online  03 June 2013 
Exact, singularityfree recasting of the Newtonian potential in continuous media
^{1} Univ. Bordeaux, LAB, UMR 5804, 33270 Floirac, France
^{2} CNRS, LAB, UMR 5804, 33270 Floirac, France
email: jeanmarc.hure@obs.ubordeaux1.fr
Received: 22 January 2013
Accepted: 17 April 2013
The gravitational potential is a key function involved in many astrophysical problems. Its evaluation inside continuous media from Newton’s law is known to be challenging because of the diverging kernel 1/ r − r′. This difficulty is generally treated with avoidance techniques (e.g. multipole expansions, softening length), which are themselves not without drawbacks. In this article, we present a new path that basically fixes the pointmass singularity problem in systems with at least two dimensions. It consists of recasting the gravitational potential ψ in an equivalent integrodifferential form,
where (q_{1},q_{2}) is a pair of independent spatial variables (linear and/or angular), f is a known function, and ℋ is an auxiliary scalar function. In contrast with ψ, this “hyperpotential” ℋ is the convolution of the mass density with a finite amplitude kernel κ. We show that closedform expressions for κ can be directly deduced from the potential of homogeneous sheets. We then give a few formulae appropriate to the Cartesian, cylindrical and spherical coordinate systems, including axial symmetry. The method is essentially not limited, either on the geometry of the source or on the distribution, and its implementation is straightforward. Several tests based upon simple quadrature/differentiation schemes are presented (the homogeneous rectangular sheet, cuboid and disk, the Maclaurin disk and a truncated LaneEmden solution). Compared with a direct summation, the extra computational cost is low and the gain is real: no truncated series, no free parameter, and a relative accuracy better than 1% for typically 16 nodes per spatial direction using the most basic numerical schemes.
Key words: gravitation / methods: analytical / methods: numerical
© ESO, 2013
1. Introduction
Gravitation plays an essential role in the evolution of most astrophysical systems, from aggregates and dusty planetary rings to rotating stars, supermassive black holes in active nuclei and galactic clusters (Hachisu 1986; Kozhanov 2004; Colwell et al. 2006; King 2010; Comito et al. 2011). In the investigation of various dynamical problems and equilibrium configurations from first integrals and energy equations, the potential appears as a fundamental scalar function. In continuous media, it is naturally accessible through an integral, namely (1)where the kernel sweeps aways the point mass singularity – a direct consequence of Newton’s inverse square law. Indeed, the potential integral is convergent for most density distributions of physical interest (Kellogg 1929; Durand 1953; Binney & Tremaine 1987). The singularity problem is inherent in the discretizationcounting technique usually adopted. By dividing the system into small massive elements and summing over all individual contributions, it is difficult to estimate precisely the influence of any small element upon itself (i.e., “selfgravity”), which is not ameliorate a lot by increasing the resolution.
The point mass singularity can be avoided in various ways that are more or less faithful to Newton’s law. The multipole expansion of the kernel is one of the most valuable theoretical tools in potential theory (Kellogg 1929; Durand 1953; Cohl et al. 2001; Aksenov 1999). It is extremely efficient outside the material domain and a few terms often suffice to reach computer precision. Inside and even in close neighborhood, however, the convergence of the series is known to be poor because r/r′ ≈ 1. Because the series is an alternate series, convergence is much delayed and truncations are critical (Clement 1974). Low convergence is a common property of multipole expansions and is observed in various contexts other than gravitation (Wuensche 1975; Kosov & Popelier 2000; Gramada & Bourne 2010). Users of multipole expansions generally need to incorporate a large number of terms – tens to hundred typically – before accuracy becomes acceptable (Hachisu 1986; Stone & Norman 1992; Mach & Malec 2012). Because the number of integrals to estimate is equal to the number of terms, the computational time increases linearly. Another option to derive ψ is the Poisson equation, which is rapidly solved with specific algorithms (Stone & Norman 1992; Storzer 1993; Spotz 1995; Briggs et al. 2000; Matsumoto & Hanawa 2003; Jusélius & Sundholm 2007; Guillet & Teyssier 2011). Nevertheless, Poissonsolvers are not always “selfstarting”, meaning they require precise boundary conditions only the integral approach can furnish. Another drawback is the shape of astrophysical bodies, which are often complex and not systematically match the numerical meshes (for techniques based on mapping, see Grandclément et al. 2001; Reese 2006). In particular, the Poisson equation is threedimensional by nature, and not wellsuited to problems in one and two dimensions.
In this paper, we describe a novel path for determining the Newtonian potential of a continuous system by recasting Eq. (1). The new form does not involve any singular kernel, series, or “softening length”, but just the crossderivative of the mass density convolved with a finite amplitude kernel. The recasting is exact and general in the sense that i) it preserves the Newtonian character of the interactions at all scales; and ii) it applies to any density distribution and morphology (shape and number of dimensions larger than one). This is therefore a new tool for both numerical applications and theoretical investigations in various domains of astrophyscis (e.g., simulations, generation of approximations, determination of potential/density pairs) and physics as well. This paper goes beyond the analysis presented in Huré & Dieckmann (2012), which was restricted to axial symmetry. We consider here i) a generic treatment of the regularization step, regardless of the system of coordinates; ii) a full threedimensional approach; iii) a simple recipe to determine the finite amplitude kernel; and iv) a direct application to the Cartesian, cylindrical, and spherical coordinates.
The paper is organized as follows. In Sect. 2, we recall the integral expression for the Newtonian potential of a continuous distribution. We formally describe the recasting of the potential integral based upon the properties of Newton’s law (symmetry and independant spaces). The application to the Cartesian, cylindrical, and spherical coordinate systems is the aim of Sect. 3. In particular, we illustrate the method by considering a few testcases, mostly of astrophysical interest, by using deliberately loworder numerical schemes (our goal is not to perform a critical study of the most efficient techniques for quadratures and differentiations). A conclusion summarizes the results and mentions possible issues to consider next. A few appendices contain formulae and demonstrations.
2. Recasting of the potential integral
The Newtonian potential at a point P(r) in space of a body is given by Eq. (1). The integral extends over the material domain Ω′ (including its boundary), i.e., ψ(r) ≡ ψ(r;Ω′), dm′(r′) is the elementary mass at P′(r′) ∈ Ω′,  r − r′  = PP′, and is the constant of gravity. The configuration is illustrated by Fig. 1. As is well known, PP′ vanishes everywhere inside Ω′, making the kernel singular, while the potential is, most of the time, a finite function of space (see e.g. Kellogg 1929).
Fig. 1 Field point P(r), source point P′(r′) belonging to the material domain Ω′, and elementary mass dm′. Inside Ω′, the separation PP′ vanish. 

Open with DEXTER 
2.1. Idea behind recasting and strategy
As Eq. (1) shows, there are two spaces in potential theory: i) the space of field points where the potential is requested (hereafter, the Pspace); and ii) the space of source points that describes the source (hereafter, the P′space). These spaces are superimposed in practice – this is the physical space –, but are decoupled mathematically. Indeed, when estimating the potential from Eq. (1), r is held fixed while the integration is performed in the P′space. The pointmass singularity, of hyperbolictype, can be regularized using two successive integrations in orthogonal directions (as a proof, note that the potential of flat or curved homogeneous sheets is generally a finite function of space and source parameters). The idea is then to integrate the Newton kernel in the Pspace until the singularity is finally suppressed. This operation is necessarily possible since the Newton kernel is symmetrical with respect to r and r′. Concretely, if P has coordinates (q_{1},q_{2},q_{3}), then the regularization step becomes^{1}(2)where f is introduced for convenience (see below) and it is a function of P only. At this stage, the coordinates (q_{1}′,q_{2}′,q_{3}′) of P′ are regarded as parameters, and κ^{q1q2} must be a function of r. Since the regularization is performed in the Pspace, it is made regardless of the mass distribution, which is especially attractive. The new kernel κ^{q1q2} (hereafter the “hyperkernel”) has, by construction, a finite amplitude and can be convolved with the mass density. This is the convolution step: (3)where the factor is introduced for convenience (see below). This integral produces an auxiliary scalar function, ℋ (hereafter, the “hyperpotential”). The Newtonian potential is then recovered by reversing the regularizationstep. This is the recovering step: (4)The advantage of this approach is twofold: the singularity is circumvented, and at the same time, it is accounted for exactly. In practice, the absence of diverging kernel renders step 2 easier than with the Newton kernel. Because convolutions produce smooth functions, step 3 is also expected to be uncomplicated. Step 1 is by far the most critical, but it is made once only provided the hyperkernel is analytical (there is no interest in the recasting if κ^{q1q2} is to be determined by numerical means). The gravitational potential is finally found from steps 2 and 3. The extracost is therefore low: there is only an additional differentiation compared to the classical approach, but the singularity is correctly managed.
2.2. Note on Chandrasekhar and Lebovitz superpotentials
Our approach may evoke some aspects of the theory developped in Chandrasekhar & Lebovitz (1962) and subsequent papers (see also Chandrasekhar 1987). These authors have shown that the Newtonian potential is the trace of a symmetric tensor potential, but is also the source term of a Poisson equation, namely (5)where the solution is (6)It is clear that ℋ and χ (called “superpotential”) share two common properties: i) they are the convolution of the density field by a finite amplitude kernel; and ii) they exactly reproduce the gravitational potential by partial differencing. However, the present recasting differs from the theory of superpotentials on the following points:

i)
ψ is determined here through a single secondorder partial derivative (not three);

ii)
the recasting is not limited to 3Dproblems, but works for 2Dproblems as well;

iii)
it is not really specific to a particular system of coordinates, while Cartesian coordinates are mostly used in Chandrasekhar & Lebovitz (1962).
In some sense, the hyperpotential ℋ is some kind of optimized version of the χfunction, especially designed for numerical applications which is, initially, our main motivation.
3. Results and examples
3.1. Derivation of hyperkernels, and the link with the potential of homogeneous sheets
According to Eq. (2), the recasting depends on the capability to determine analytically an expression for κ^{q1q2} associated with a given pair (q_{1},q_{2}) of coordinates, preferably a closedform. There are many possibilities, in particular because of the presence of the function f(q_{1},q_{2},q_{3}), which adds a degree of freedom. If we define f such that (7)is an area element, then (8)where is a surface q_{3} = const. in the Pspace. We see that Eq. (8) is nothing but, up to a factor , the formula for the gravitational potential of a homogeneous surface with unit surface density dm/d^{2}A, except that the role of the Pspace and P′space is exchanged. To obtain a formula for κ^{q1q2}, it is sufficient to extract from the list of known potential/density pairs those that correspond to a bidimensional distribution (i.e., a sheet) and constant surface density. As Eq. (8) suggests, there is no special constraint on the shape and size of the sheet (i.e., flat, curved, rectangular, circular, etc.). Nevertheless, from a practical point of view, it seems preferable that and Ω′ be geometrically “compatible” enough to facilitate the convolutionstep. Since there is certain freedom in selecting the integral bounds in Eq. (8), it is also better to consider a finite size (and finite mass) sheet.
3.2. An easy implementation
For each point P where ψ is requested, the sequence of operations is the following:

1.
computing the hyperkernel κ^{q1q2} from an appropriate formula (this depends on the coordinate system; see next section);

2.
estimating the hyperpotential ℋ from Eq. (3) for the actual density distribution (a volume density ρ or a surface density Σ). A quadrature scheme is needed. This is the first place where numerical errors are generated;

3.
estimating the crossderivative of ℋ. A differentiation scheme is needed. This requires determining hyperpotentialvalues in the vicinity of the actual point P. This is the second place where numerical errors are generated.
Clearly, all techniques curently used to compute ψ from Eq. (1) can be employed for ℋ. Because ℋ is a convolution (see below), fast specific algorithms coupled with highperfomance differentiation schemes can probably be envisaged. Tests presented in the following will employ the most basic schemes, which have produced good results.
Pairs (q_{1},q_{2}), function f, and associated area element d^{2}A for the three most popular coordinate systems.
Fig. 2 Notations for the Cartesian, cylindrical, and spherical coordinates: Pspace (left) and P′space (right). 

Open with DEXTER 
Fig. 3 The hyperkernel κ^{xy} in the (x,y)plane in the vicinity of point P′ (magnified by a factor 5), for z = z′ and three different values of y − y′ labelled on the curves. The Newton kernel which diverges for  r − r′  = 0 is shown for comparison. 

Open with DEXTER 
Fig. 4 Error index ϵ on potential values in the plane of the homogeneous square sheet with vertices at (boundary in white). The mean error index is − 2.52 (dashed line). 

Open with DEXTER 
3.3. Results: the case of Cartesian coordinates
Table 1 lists pairs (q_{1},q_{2}), the associated area element d^{2}A, and the function f appropriate for the Cartesian, cylindrical and spherical coordinates that are most often used. Other coordinate systems and geometries can obviously be considered as well. The notations are summarized in Fig. 2. In Cartesian coordinates, the recasting that corresponds to the pair (q_{1},q_{2}) = (x,y) is (9)and other pairs can be considered by permutation. As argued above, we can determine κ^{xy} by considering a surface q_{3} ≡ z = const., i.e., a flat horizontal sheet, and the most natural choice is the rectangular shape. The formula for the potential is known in that case (e.g. Durand 1953). It is reproduced in Appendix A. The hyperkernel is deduced by exchanging P and P′, which leads to Eq. (A.4). Figure 3 displays κ^{xy} in the (x,y)plane in the vicinity of the point P′(x′,y′), for three different values of y − y′ and for z − z′ = 0. The Newton kernel 1/  r − r′  is also shown for comparison. We can see that the hyperkernel is a smooth function with finite amplitude. In particular, it is zero at zero relative separation. Any type of mass density profile can be injected in Eq. (9), bi or threedimensional, uniform or not.
We illustrate the potential recasting with two simple examples. As a first test, we consider a square sheet in the (x,y)plane with constant surface density Σ_{0} (or dm′ = Σ_{0}dx′dy′), length unity, and centered on the origin, with vertices at . It is discretized on a N′ × M′ grid with regular spacing in each direction. The Pgrid is made of N × M points with uniform spacing as well. The quadratures and partial derivatives are determined by secondorder schemes that are among the most basic ones (e.g. Press et al. 1992). Figure 4 shows the error index (10)between ψ determined from Eq. (9) and the reference potential ψ_{e} (2 × 10^{16} is for doubleprecision computations). This case is for N′ = M′ = 32 and N = M = N′ + 2, which leaves just one point outside the sheet, left and right, bottom and top. We see that the relative error is rather uniform inside the material domain, of about 0.1%. Close to the edges of the sheet, the error rises slightly. Outside Ω′, the accuracy is still uniform, but typically better by an order of magnitude.
As a second test, we consider a cube with uniform density ρ_{0} (i.e., dm′ = ρ_{0}dx′dy′dz), length unity, and vertices at . The numerical setup is the same as for the sheet, and z = 0. The reference potential is also known for this 3D body (MacMillan 1930; Waldvogel 1976). Figure 5 displays the error index versus x and y in the cube’s midplane. Again, we notice that the relative error is uniform, with 0.2% typically inside the body, and a factor 10 better outside. Edge effects are less marked than in 2D. The integration of the kernel in the third direction smoothes the errors, and the potential is now derivable when crossing the lateral faces of the cube. The Fortran 90 program used in these two examples is available upon request.
Fig. 5 Same conditions and same color code as for Fig. 4, but for the homogeneous cuboid with vertices at . The mean value is − 3.00 (dashed line). 

Open with DEXTER 
3.4. The hyperpotential is a convolution product
As shown in Appendix A, when we set X = x − x′, Y = y − y′ and Z = z − z′, the hyperkernel κ^{xy} becomes where , and so the hyperpotential writes in the 3D case (11)Since ρ = 0 outside Ω′, the integral bounds can be safely changed for ± ∞. We then conclude that ℋ is a convolution product. This is expected because the potential itself is a convolution product (Binney & Tremaine 1987; Hackbusch et al. 2010). We have (12)This result is independent of the coordinate system, namely: (13)
3.5. Results in cylindrical and spherical coordinates
In curved geometries, there are apparently fewer options. The reason is that a few formula for the potential of canonical surfaces are missing yet, and it is hard to find closedform expressions for the hyperkernel by direct integration of Eq. (2). One can probably use a series representation instead, but any truncation is expected to produce an approximate potential. Surfaces q_{3} = const. of particular interest are (see Fig. 2 and Table 1)

1.
in cylindrical coordinates:

(a)
a piece of a hollow cylinder (surfaceR = const.);

(b)
a meridional sheet (surface θ = const.);

(c)
a polar sector (surface z = 0);

(a)

2.
in spherical coordinates:

(a)
a piece of spherical shell (surface r = const.);

(b)
a meridional sector (surface θ = const.);

(c)
a piece of cone (surface φ = const.).

(a)
3.6. Axial symmetry
If the source is axially symmetrical (i.e., ∂_{θ′}ρ = 0), the integration of the Newton kernel in the P′space over the polar angle θ′ ∈ [0,2π] leads to (14)where K is the complete elliptic integral of the first kind, δ^{2} = (R′ + R)^{2} + (z − z′)^{2}, and . Because this function is hyperbolically singular when k → 1, we can determine an hyperpotential by an integration in the Pspace. Again, there are not many options because we lack formulae for the potential of the hollow cylinder, for the cone, and for the piece of spherical shell (see Table 1). Fortunately, there is the formula for the potential of the circular disk, i.e., a closedform for (Durand 1953; Krough et al. 1982; Lass & Blitzer 1983; Huré 2012). We can therefore deduce an axially symmetrical hyperkernel κ^{R} by exchanging the role of P and P′ (see also Appendix C). This leads to Eq. (C.6). The potential is also axially symmetrical, and it finally writes (15)where (16)We now present three last examples. Figure 6 is for the circular disk with radius unity. The Pgrid and the P′grid are made of N′ = 32 and N = N′ + 1 points equally spaced in R^{2} and R^{′2} respectively – this radial scale is natural in this type of problem, both for the convolution and for the derivative. As above, the two grids coincide inside Ω′ (there is just one point outside the disk) and the numerical schemes are secondorder. The error index is shown at four different altitudes, including for the disk midplane (i.e. z = z′). The top panel shows the homogeneous disk. We see that the relative error is, on average, of about 0.1%, with again a slight degradation near the edge. The bottom panel shows the Maclaurin disk where . For this inhomogeneous case, the reference solution is taken from Schulz (2009). We see that the relative error is about 5 × 10^{3}, which is not as good as in the homogeneous case. The error is almost insensitive to the altitude from the disk plane. This is due to the actual surface density profile (quadrature schemes are generally not very efficient to manage infinite derivatives).
Fig. 6 Error index for the homogeneous disk (top panel), and for the (inhomogeneous) Maclaurin disk (bottom panel). The setup is the same in both cases (edge at R′ = 1), and N′ = 32. The computational grid has N = 33 points, but the same spacing. 

Open with DEXTER 
Fig. 7 Same legend and same color code as for Fig. 5, but a radially inhomogeneous sphere (boundary in white; see text). The mean index is −2.83 (dashed line). The Pgrid (R,z) is a square, larger than the sphere’s radius by 10%. 

Open with DEXTER 
The last example is an inhomogeneous sphere with radius unity. The potential/density pair is the following (17)which corresponds to the solution of the LaneEmden equation with polytropic exponent γ = 2 (or index n = 1), truncated at the first zero. The hyperpotential in spherical coordinates writes (18)where κ^{R} is the cylindrical hyperkernel κ^{R}. This is a typical example where the surface (a disk) and the domain Ω′ (concentric shells) are somewhat disconnected. Here, the sphere is discretized into N′ × M′ points equally spaced in the (r′,φ′)plane. Potential values are determined from Eq. (15) in the (R^{2},z)plane at N × M points equally spaced; the computational box is larger than the sphere’s radius by 10%. Figure 7 shows the error index for N′ = M′ = 32 and N = M = 32. We see that the deviation is remarkably homogeneous inside and outside the sphere. The relative error is about 0.2% inside and outside the material domain. This is comparable to the case of the cube.
Centrally symmetrical configurations can also be treated by using an hyperkernel, but there is nothing really new here (see Appendix E).
4. Summary and concluding remarks
We demonstrated that the Newtonian potential of continuous bodies can be determined from the partial crossdifferentiation of the mass density convolved with a finite amplitude kernel (a hyperkernel), regardless of any coordinate system. The recasting of Newton’s integral is free of singularity and exact, and it applies to any type of two and threedimensional systems. Provided the hyperkernel is analytical, the extracost with respect to direct estimations is weak or negligible: it is only N operations vs. N^{2} in a grid with N points. It is much lower if the method is only used to generate boundary conditions (and coupled with Poissonsolvers based on FFTs). The gain in accuracy is huge since i) direct estimates cannot avoid errors; ii) there is no free parameter; and iii) there is no trunctaed series. We have given a few examples in Cartesian, cylindrical, and spherical coordinates that prove the efficiency even with loworder quadrature and differentiation schemes. The recasting is therefore very attractive for any numerical applications. It is also a new tool for investigating various theoretical problems and derive potential/density pairs or approximations.
This work can therefore be continued and improved in several ways. Knowing analytical expressions for the hyperkernel associated with a given pair of orthogonal coordinates is the critical point of the method. We have shown that hyperkernels can be directly generated by considering the potential of homogeneous sheets, while there are doutblessly other possibilities. For Cartesian coordinates, all hyperkernels are known. In cylindrical and spherical geometries, a few closedform expressions are apparently lacking yet (hollow cylinder and cone for instance). It would therefore be interesting to investigate this kind of question. The formula for the polar sector should be helpful for most astrophysical applications however, such as for modelling rotating fluids. Other coordinate systems and geometries can be envisaged. For instance, the potential of inhomogeneous elliptic bodies can be determined, as done under central symmetry (see the end of Sect. 3.6), from the theory of thin homeoids (e.g. MacMillan 1930). This show the importance of seeking new potential/density pairs associated to 2Dsystems. It would also be interesting to analyze in more detail the numerical implementation of the method. There is obviously a wide panel of techniques at our disposal to perform quadratures/convolutions and differentiations, finite differences (as considered here), spectral methods, etc.
Finally, applications exceed the astrophysical context of gravitation. The approach is obviously suited to electrostatics and to electromagnetism since the potential vector is (19)where u is the velocity of electric charges. It is also transposable to incompressible hydrodynamics where the pressure p obeys a Poisson equation too, (20)where u is the fluid velocity.
Under invariance, the hyperbolic singularity can be converted into a logarithmic singularity that is subsequently regularized by a single integration, i.e., , but this is a special case (see Sect. 3).
Acknowledgments
It is a pleasure to thank A. Dieckmann, M. Gazeau, F. Hersant, D. Pfenniger, and A. Pierens. I sincerely thank the referee for valuable scientific comments and advice to improve the organization of the paper.
References
 Aksenov, A. G. 1999, Astron. Lett., 25, 185 [NASA ADS] [Google Scholar]
 Binney, J., & Tremaine, S. 1987, Galac. Dyn. (Princeton, NJ: Princeton University Press), 747 [Google Scholar]
 Briggs, W. L., Henson, V. E., & McCormick, S. F. 2000, A multigrid tutorial, 2nd edn. (Philadelphia, PA, USA: Society for Industrial and Applied Mathematics) [Google Scholar]
 Chandrasekhar, S. 1987, Ellipsoidal figures of equilibrium (New York: Dover) [Google Scholar]
 Chandrasekhar, S., & Lebovitz, N. R. 1962, ApJ, 135, 238 [NASA ADS] [CrossRef] [Google Scholar]
 Clement, M. J. 1974, ApJ, 194, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Cohl, H. S., Rau, A. R. P., Tohline, J. E., et al. 2001, Phys. Rev. A, 64, 052509 [NASA ADS] [CrossRef] [Google Scholar]
 Colwell, J. E., Esposito, L. W., & Sremčević, M. 2006, Geophys. Res. Lett., 33, 7201 [NASA ADS] [CrossRef] [Google Scholar]
 Comito, C., Tanga, P., Paolicchi, P., et al. 2011, Mem. Soc. Astron. It. Suppl., 16, 84 [Google Scholar]
 Durand, E. 1953, Electrostatique. Les distributions (Ed. Masson), I [Google Scholar]
 Gramada, A., & Bourne, P. 2010, in APS California Section Meeting Abstracts, D4001 [Google Scholar]
 Grandclément, P., Bonazzola, S., Gourgoulhon, E., & Marck, J.A. 2001, J. Comp. Phys., 170, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Guillet, T., & Teyssier, R. 2011, J. Comput. Phys., 230, 4756 [NASA ADS] [CrossRef] [Google Scholar]
 Hachisu, I. 1986, ApJS, 61, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Hackbusch, W., Naraparaju, K. K., & Schneider, J. 2010, J. Numer. Math., 18, 257 [CrossRef] [Google Scholar]
 Huré, J.M. 2012, Celest. Mech. Dyn. Astron., 114, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Huré, J.M., & Dieckmann, A. 2012, A&A, 541, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jusélius, J., & Sundholm, D. 2007, J. Chem. Phys., 126, 4101 [Google Scholar]
 Kellogg, O. D. 1929, Foundations of Potential Theory (NewYork: Frederick Ungar Publishing Company) [Google Scholar]
 King, A. R. 2010, MNRAS, 408, L95 [NASA ADS] [CrossRef] [Google Scholar]
 Kosov, D. S., & Popelier, P. L. A. 2000, J. Chem. Phys., 113, 3969 [NASA ADS] [CrossRef] [Google Scholar]
 Kozhanov, T. S. 2004, in Order and Chaos in Stellar and Planetary Systems, eds. G. G. Byrd, K. V. Kholshevnikov, A. A. Myllri, I. I. Nikiforov, & V. V. Orlov, ASP Conf. Ser., 316, 284 [Google Scholar]
 Krough, F. T., Ng, E. W., & Snyder, W. V. 1982, Celest. Mech., 26, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Lass, H., & Blitzer, L. 1983, Celest. Mech., 30, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Mach, P., & Malec, E. 2012, A&A, 541, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MacMillan, W. 1930, The theory of the potential, Theoretical mechanics (McGrawHill Book Company, inc.), 2 [Google Scholar]
 Matsumoto, T., & Hanawa, T. 2003, ApJ, 583, 296 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN. The art of scientific computing [Google Scholar]
 Reese, D. 2006, Thèse, Université Paul Sabatier – Toulouse III [Google Scholar]
 Schulz, E. 2009, ApJ, 693, 1310 [NASA ADS] [CrossRef] [Google Scholar]
 Spotz, W. F. 1995, HighOrder Compact Finite Difference Schemes for Computational Mechanics, P.H.D. Thesis, University of Texas at Austin, Austin, TX [Google Scholar]
 Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [NASA ADS] [CrossRef] [Google Scholar]
 Storzer, H. 1993, A&A, 271, 25 [NASA ADS] [Google Scholar]
 Waldvogel, J. 1976, ZAMP, 27, 867 [NASA ADS] [CrossRef] [Google Scholar]
 Wuensche, A. 1975, ZAMP, 55, 301 [Google Scholar]
Appendix A: Hyperkernel for the rectangular sheet (Cartesian coordinates)
The notations are those of Fig. 2. Up to a factor , the potential of a homogeneous rectangular sheet with unity surface density is found from the integral (A.1)where the bounds represents the coordinates of the four corners of the rectangular sheet. A closedform expression is found for instance in Durand (1953). If we set X = x − x′, Y = y − y′ and Z = Z − z′, we have and dx′dy′ = dXdY. The indefinite integral is (A.2)For the present purpose, we only need to generate the Newton kernel by a mixed partial derivative ∂^{2}/∂x∂y, and so, we have a certain liberty in choosing the most convenient integral bounds. Here, we take x and x′ for the integral over x′, and y and y′ for the integral over y. from this, we obtain (A.3)\arraycolsep1.75ptConsequently, the hyperkernel is obtained from this expression by exchanging the variables x′ and x, and y′ and y. We find (A.4)To generate κ^{yz} associated with a potential expressed as ∂^{2}ℋ/∂y∂z, we perform the permutations (x,x′) ↔ (z,z′) in the formulae above. To generate κ^{xz} associated with a potential expressed as ∂^{2}ℋ/∂x∂z, we perform the permutations (y,y′) ↔ (z,z′).
Appendix B: Hyperkernel for the meridional sheet (cylindrical coordinates)
The notations are those of Fig. 2. First, we use Eq. (A.4) and perform the permutations (y,y′) ↔ (z,z′). We find (B.1)Then, we apply a counterclockwise rotation (i.e., positive trigonometric sense) around the zaxis by an angle θ′. We obtain the meridional sheet. The hyperkernel then writes (B.2)and the potential is given by (B.3)
Appendix C: Hyperkernel for the polar sector (cylindrical coordinates)
The notations are those of Fig. 2. Again up to a factor , the potential of a circular sector is found from the formula (C.1)where (C.2)This double integral has been calculated in Huré (2012). The indefinite form is (C.3)where F(φ,k), E(φ,k) and Π(φ,m^{2},k) are the incomplete elliptic integral of the first, second, and third kinds, respectively, δ = (R′ + R)^{2} + ζ^{2}, ζ = z − z′, , 2β = π − (θ − θ′). To generate κ^{Rθ}, it is sufficient to consider the following integral bounds: 0 and R′ for the radial integration, and θ − π and θ′ for the angular part. Next, the source point and the field point are exchanged (note that δ, k, ζ^{2} and m are not impacted). We finally obtain (C.4)where 2β′ = π − (θ′ − θ) = 2π − 2β. The potential is then given by (C.5)Under axial symmetry, we have (C.6)where E, K and Π are the complete elliptic integrals of the first, second, and third kinds, respectively.
Appendix D: Hyperkernel for the meridional sector (spherical coordinates)
The notations are those of Fig. 2. Up to a factor , the potential of a meridional sector defined by θ′ = const., is found from the double integral (D.1)with convenient bounds. We can calculate this expression directly from the formula for the polar sector in cylindrical coordinates (see Appendix C). First, we write the relative separation  r − r′  in the following form (D.2)where (D.3)Then, we notice that the similarity between Eqs. (C.1) and (D.1) is perfect if we make the following substitutions (D.4)and then (D.5)From Eq. (C.3), we see that Eq. (D.1) becomes (D.6)We then derive from Eq. (C.4) by making the same substitutions, and the potential is given by (D.7)
Appendix E: Central symmetry
When ∂_{r′}ρ = 0, we can derive an hyperkernel by considering the potential of a spherical shell with unity surface density. From the Gauss theorem, we easily find (still up to ) where H is the Heaviside function. The hyperkernel κ^{shell} is then obtained by exchanging r and r′ in this expression. In this case, we have simply ψ = ℋ with (E.1)where the integral bounds are the inner radius and outer radius of the sphere. This result is well known (e.g. Binney & Tremaine 1987).
All Tables
Pairs (q_{1},q_{2}), function f, and associated area element d^{2}A for the three most popular coordinate systems.
All Figures
Fig. 1 Field point P(r), source point P′(r′) belonging to the material domain Ω′, and elementary mass dm′. Inside Ω′, the separation PP′ vanish. 

Open with DEXTER  
In the text 
Fig. 2 Notations for the Cartesian, cylindrical, and spherical coordinates: Pspace (left) and P′space (right). 

Open with DEXTER  
In the text 
Fig. 3 The hyperkernel κ^{xy} in the (x,y)plane in the vicinity of point P′ (magnified by a factor 5), for z = z′ and three different values of y − y′ labelled on the curves. The Newton kernel which diverges for  r − r′  = 0 is shown for comparison. 

Open with DEXTER  
In the text 
Fig. 4 Error index ϵ on potential values in the plane of the homogeneous square sheet with vertices at (boundary in white). The mean error index is − 2.52 (dashed line). 

Open with DEXTER  
In the text 
Fig. 5 Same conditions and same color code as for Fig. 4, but for the homogeneous cuboid with vertices at . The mean value is − 3.00 (dashed line). 

Open with DEXTER  
In the text 
Fig. 6 Error index for the homogeneous disk (top panel), and for the (inhomogeneous) Maclaurin disk (bottom panel). The setup is the same in both cases (edge at R′ = 1), and N′ = 32. The computational grid has N = 33 points, but the same spacing. 

Open with DEXTER  
In the text 
Fig. 7 Same legend and same color code as for Fig. 5, but a radially inhomogeneous sphere (boundary in white; see text). The mean index is −2.83 (dashed line). The Pgrid (R,z) is a square, larger than the sphere’s radius by 10%. 

Open with DEXTER  
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.