Subscriber Authentication Point
Free Access
 Issue A&A Volume 554, June 2013 A45 9 Numerical methods and codes https://doi.org/10.1051/0004-6361/201321146 03 June 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 discretization-counting 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., “self-gravity”), 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, Poisson-solvers are not always “self-starting”, 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 three-dimensional by nature, and not well-suited 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 cross-derivative 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 three-dimensional 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 test-cases, mostly of astrophysical interest, by using deliberately low-order 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 P-space); 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 point-mass singularity, of hyperbolic-type, 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 P-space 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 (q1,q2,q3), then the regularization step becomes1(2)where f is introduced for convenience (see below) and it is a function of P only. At this stage, the coordinates (q1′,q2′,q3′) of P′ are regarded as parameters, and κq1q2 must be a function of r. Since the regularization is performed in the P-space, 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 regularization-step. 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 extra-cost 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 second-order partial derivative (not three);

• ii)

the recasting is not limited to 3D-problems, but works for 2D-problems 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 (q1,q2) of coordinates, preferably a closed-form. There are many possibilities, in particular because of the presence of the function f(q1,q2,q3), which adds a degree of freedom. If we define f such that (7)is an area element, then (8)where is a surface q3 = const. in the P-space. 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/d2A, except that the role of the P-space 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 bi-dimensional 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 convolution-step. 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 cross-derivative of ℋ. A differentiation scheme is needed. This requires determining hyperpotential-values 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 high-perfomance differentiation schemes can probably be envisaged. Tests presented in the following will employ the most basic schemes, which have produced good results.

Table 1

Pairs (q1,q2), function f, and associated area element d2A for the three most popular coordinate systems.

 Fig. 2 Notations for the Cartesian, cylindrical, and spherical coordinates: P-space (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 (q1,q2), the associated area element d2A, 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 (q1,q2) = (x,y) is (9)and other pairs can be considered by permutation. As argued above, we can determine κxy by considering a surface q3 ≡ 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 three-dimensional, 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′ = Σ0dx′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 P-grid is made of N × M points with uniform spacing as well. The quadratures and partial derivatives are determined by second-order 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 double-precision 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′ = ρ0dx′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 closed-form 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 q3 = 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);

• 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.).

For cases 1a, 2a, and 2c (with ), the potential is apparently not known in closed-form; this would be helpful. We have no hyperkernel to propose. This question remains open. Case 1b is accessible since the meridional sheet is nothing but a rectangular sheet in the plane (x,z), rotated counter-clockwise by an angle θ′. The formula for κRz can then be deduced from the Cartesian case (see Appendix B); this is Eq. (B.2). Case 1c can also be treated since the potential of a polar sector has been derived in Huré (2012). The formula is reproduced in the Appendix C, and the hyperkernel κ is given by Eq. (C.4). Finally, case 2b is feasible since the meridional sector is a polar sector. We can then use the result for κ established in cylindrical coordinates and apply convenient rotations to derive κ. A more direct calculus is presented in Appendix D.

### 3.6. Axial symmetry

 Fig. 6 Error index for the homogeneous disk (top panel), and for the (inhomogeneous) Maclaurin disk (bottom panel). The set-up 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 P-grid (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 Lane-Emden 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 (R2,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 cross-differentiation 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 three-dimensional systems. Provided the hyperkernel is analytical, the extra-cost with respect to direct estimations is weak or negligible: it is only N operations vs. N2 in a grid with N points. It is much lower if the method is only used to generate boundary conditions (and coupled with Poisson-solvers 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 low-order 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 closed-form 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 2D-systems. 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.

1

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

1. Aksenov, A. G. 1999, Astron. Lett., 25, 185 [NASA ADS] [Google Scholar]
2. Binney, J., & Tremaine, S. 1987, Galac. Dyn. (Princeton, NJ: Princeton University Press), 747 [Google Scholar]
3. 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]
4. Chandrasekhar, S. 1987, Ellipsoidal figures of equilibrium (New York: Dover) [Google Scholar]
5. Chandrasekhar, S., & Lebovitz, N. R. 1962, ApJ, 135, 238 [NASA ADS] [CrossRef] [Google Scholar]
6. Clement, M. J. 1974, ApJ, 194, 709 [NASA ADS] [CrossRef] [Google Scholar]
7. Cohl, H. S., Rau, A. R. P., Tohline, J. E., et al. 2001, Phys. Rev. A, 64, 052509 [NASA ADS] [CrossRef] [Google Scholar]
8. Colwell, J. E., Esposito, L. W., & Sremčević, M. 2006, Geophys. Res. Lett., 33, 7201 [NASA ADS] [CrossRef] [Google Scholar]
9. Comito, C., Tanga, P., Paolicchi, P., et al. 2011, Mem. Soc. Astron. It. Suppl., 16, 84 [Google Scholar]
10. Durand, E. 1953, Electrostatique. Les distributions (Ed. Masson), I [Google Scholar]
11. Gramada, A., & Bourne, P. 2010, in APS California Section Meeting Abstracts, D4001 [Google Scholar]
12. Grandclément, P., Bonazzola, S., Gourgoulhon, E., & Marck, J.-A. 2001, J. Comp. Phys., 170, 231 [NASA ADS] [CrossRef] [Google Scholar]
13. Guillet, T., & Teyssier, R. 2011, J. Comput. Phys., 230, 4756 [NASA ADS] [CrossRef] [Google Scholar]
14. Hachisu, I. 1986, ApJS, 61, 479 [NASA ADS] [CrossRef] [Google Scholar]
15. Hackbusch, W., Naraparaju, K. K., & Schneider, J. 2010, J. Numer. Math., 18, 257 [CrossRef] [Google Scholar]
16. Huré, J.-M. 2012, Celest. Mech. Dyn. Astron., 114, 365 [NASA ADS] [CrossRef] [Google Scholar]
17. Huré, J.-M., & Dieckmann, A. 2012, A&A, 541, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
18. Jusélius, J., & Sundholm, D. 2007, J. Chem. Phys., 126, 4101 [Google Scholar]
19. Kellogg, O. D. 1929, Foundations of Potential Theory (New-York: Frederick Ungar Publishing Company) [Google Scholar]
20. King, A. R. 2010, MNRAS, 408, L95 [NASA ADS] [CrossRef] [Google Scholar]
21. Kosov, D. S., & Popelier, P. L. A. 2000, J. Chem. Phys., 113, 3969 [NASA ADS] [CrossRef] [Google Scholar]
22. 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]
23. Krough, F. T., Ng, E. W., & Snyder, W. V. 1982, Celest. Mech., 26, 395 [NASA ADS] [CrossRef] [Google Scholar]
24. Lass, H., & Blitzer, L. 1983, Celest. Mech., 30, 225 [NASA ADS] [CrossRef] [Google Scholar]
25. Mach, P., & Malec, E. 2012, A&A, 541, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
26. MacMillan, W. 1930, The theory of the potential, Theoretical mechanics (McGraw-Hill Book Company, inc.), 2 [Google Scholar]
27. Matsumoto, T., & Hanawa, T. 2003, ApJ, 583, 296 [NASA ADS] [CrossRef] [Google Scholar]
28. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN. The art of scientific computing [Google Scholar]
29. Reese, D. 2006, Thèse, Université Paul Sabatier – Toulouse III [Google Scholar]
30. Schulz, E. 2009, ApJ, 693, 1310 [NASA ADS] [CrossRef] [Google Scholar]
31. Spotz, W. F. 1995, High-Order Compact Finite Difference Schemes for Computational Mechanics, P.H.D. Thesis, University of Texas at Austin, Austin, TX [Google Scholar]
32. Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [NASA ADS] [CrossRef] [Google Scholar]
34. Waldvogel, J. 1976, ZAMP, 27, 867 [NASA ADS] [CrossRef] [Google Scholar]
35. 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 closed-form 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 counter-clockwise rotation (i.e., positive trigonometric sense) around the z-axis 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 Π(φ,m2,k) are the incomplete elliptic integral of the first, second, and third kinds, respectively, δ = (R′ + R)2 + ζ2, ζ = z − z′, , 2β = π − (θ − θ′). To generate κ, 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

Table 1

Pairs (q1,q2), function f, and associated area element d2A 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: P-space (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 set-up 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 P-grid (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 (full-text 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 48-96 hours after online publication and is updated daily on week days.