© 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 $\begin{array}{ccc}\mathit{\psi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{𝒢}\mathrm{\int }\frac{\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathit{,}& & \end{array}$(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 $\mathrm{|}{r}\mathrm{-}{r}{}^{\mathrm{\prime }}\mathrm{|}^{-1}$ 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).

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 (q₁,q₂,q₃), then the regularization step becomes¹$\begin{array}{ccc}{”}_{{\mathit{q}}_{\mathrm{1}}\mathit{,}{\mathit{q}}_{\mathrm{2}}}\frac{\mathit{f}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{d}{\mathit{q}}_{\mathrm{1}}\mathrm{d}{\mathit{q}}_{\mathrm{2}}}{{\mathrm{PP}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{q}}_{\mathrm{1}}\mathit{,}{\mathit{q}}_{\mathrm{2}}\mathit{,}{\mathit{q}}_{\mathrm{3}}\mathrm{\right)}}\mathrm{\equiv }{\mathit{\kappa }}^{{\mathit{q}}_{\mathrm{1}}{\mathit{q}}_{\mathrm{2}}}\mathrm{\left(}{r}\mathrm{;}{r}\mathrm{\prime }\mathrm{\right)}\mathit{,}& & \end{array}$(2)where f is introduced for convenience (see below) and it is a function of P only. At this stage, the coordinates (q₁′,q₂′,q₃′) of P′ are regarded as parameters, and κ^q₁q₂ 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 κ^q₁q₂ (hereafter the “hyperkernel”) has, by construction, a finite amplitude and can be convolved with the mass density. This is the convolution step: $\begin{array}{ccc}\mathrm{-}\mathrm{𝒢}{\mathrm{\int }}_{{\mathrm{\Omega }}^{\mathrm{\prime }}}{\mathit{\kappa }}^{{\mathit{q}}_{\mathrm{1}}{\mathit{q}}_{\mathrm{2}}}\mathrm{\left(}{r}\mathrm{;}{r}\mathrm{\prime }\mathrm{\right)}\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}\mathrm{\equiv }\mathrm{\mathscr{H}}\mathrm{\left(}{r}\mathrm{;}{\mathrm{\Omega }}^{\mathrm{\prime }}\mathrm{\right)}\mathit{,}& & \end{array}$(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: $\begin{array}{ccc}{\mathit{\partial }}_{{\mathit{q}}_{\mathrm{1}}{\mathit{q}}_{\mathrm{2}}}^{\mathrm{2}}\mathrm{\mathscr{H}}& \mathrm{=}& \mathrm{-}\mathrm{𝒢}{\mathit{\partial }}_{{\mathit{q}}_{\mathrm{1}}{\mathit{q}}_{\mathrm{2}}}^{\mathrm{2}}{\mathrm{\int }}_{{\mathrm{\Omega }}^{\mathrm{\prime }}}{\mathit{\kappa }}^{{\mathit{q}}_{\mathrm{1}}{\mathit{q}}_{\mathrm{2}}}\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}\mathit{,}\\ & \mathrm{=}& \mathrm{-}\mathrm{𝒢}{\mathrm{\int }}_{{\mathrm{\Omega }}^{\mathrm{\prime }}}\mathrm{(}{\mathit{\partial }}_{{\mathit{q}}_{\mathrm{1}}{\mathit{q}}_{\mathrm{2}}}^{\mathrm{2}}{\mathit{\kappa }}^{{\mathit{q}}_{\mathrm{1}}{\mathit{q}}_{\mathrm{2}}}\mathrm{)}\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}\mathit{,}\\ & \mathrm{=}& \mathrm{-}\mathrm{𝒢}\mathit{f}{\mathrm{\int }}_{{\mathrm{\Omega }}^{\mathrm{\prime }}}\frac{\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\\ & \mathrm{=}& \mathit{f\psi }\mathit{.}\end{array}$(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 κ^q₁q₂ 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 $\begin{array}{ccc}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\chi }\mathrm{=}\mathrm{-}\mathrm{2}\mathit{\psi ,}& & \end{array}$(5)where the solution is $\begin{array}{ccc}\mathit{\chi }\mathrm{\equiv }\mathrm{-}\mathrm{𝒢}{\mathrm{\int }}_{{\mathrm{\Omega }}^{\mathrm{\prime }}}\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}\mathit{.}& & \end{array}$(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 κ^q₁q₂ associated with a given pair (q₁,q₂) of coordinates, preferably a closed-form. There are many possibilities, in particular because of the presence of the function f(q₁,q₂,q₃), which adds a degree of freedom. If we define f such that $\begin{array}{ccc}{\mathit{d}}^{\mathrm{2}}\mathit{A}\mathrm{=}\mathit{f}\mathrm{\left(}{\mathit{q}}_{\mathrm{1}}\mathit{,}{\mathit{q}}_{\mathrm{2}}\mathit{,}{\mathit{q}}_{\mathrm{3}}\mathrm{\right)}\mathit{d}{\mathit{q}}_{\mathrm{1}}\mathit{d}{\mathit{q}}_{\mathrm{2}}& & \end{array}$(7)is an area element, then $\begin{array}{ccc}{\mathit{\kappa }}^{{\mathit{q}}_{\mathrm{1}}{\mathit{q}}_{\mathrm{2}}}\mathrm{=}{”}_{\mathrm{𝒮}}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{A}\mathrm{\left(}\mathrm{P}\mathrm{\right)}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathit{,}& & \end{array}$(8)where is a surface q₃ = 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/d²A, except that the role of the P-space and P′-space is exchanged. To obtain a formula for κ^q₁q₂, 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 κ^q₁q₂ 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.

3.3. Results: the case of Cartesian coordinates

Table 1 lists pairs (q₁,q₂), the associated area element d²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₁,q₂) = (x,y) is $\begin{array}{ccc}\mathit{\psi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}{\mathit{\partial }}_{\mathit{xy}}^{\mathrm{2}}{\mathrm{\int }}_{{\mathrm{\Omega }}^{\mathrm{\prime }}}{\mathit{\kappa }}^{\mathit{xy}}\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}\mathit{,}& & \end{array}$(9)and other pairs can be considered by permutation. As argued above, we can determine κ^xy by considering a surface q₃ ≡ 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 Σ₀ (or dm′ = Σ₀dx′dy′), length unity, and centered on the origin, with vertices at $\mathrm{(}\mathrm{\pm }\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\mathrm{\pm }\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\mathrm{0}\mathrm{)}$. 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 $\begin{array}{ccc}\mathit{ϵ}\mathrm{=}\log \left[\begin{array}{}\max \end{array}\left(\mathrm{2}\mathrm{\times }{\mathrm{10}}^{-16}\mathit{,}\left|\mathrm{1}\mathrm{-}\frac{\mathit{\psi }}{{\mathit{\psi }}_{\mathrm{e}}}\right|\right)\right]& & \end{array}$(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 ρ₀ (i.e., dm′ = ρ₀dx′dy′dz), length unity, and vertices at $\mathrm{(}\mathrm{\pm }\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\mathrm{\pm }\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\mathrm{\pm }\frac{\mathrm{1}}{\mathrm{2}}\mathrm{)}$. 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.

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 $\begin{array}{ccc}{\mathit{\kappa }}^{\mathit{xy}}& \mathrm{=}\mathrm{-}\mathit{Z}\hspace{0.17em}\mathrm{atan}\frac{\mathit{XY}}{\mathit{Z}\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathrm{+}\mathit{Y}\ln \frac{\mathit{X}\mathrm{+}\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}{\sqrt{{\mathit{Y}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}}& \\ & \mathrm{+}\mathit{X}\ln \frac{\mathit{Y}\mathrm{+}\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}{\sqrt{{\mathit{X}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}}\mathrm{\equiv }{\mathit{\kappa }}^{\mathit{xy}}\mathrm{\left(}\mathit{X,Y,Z}\mathrm{\right)}\mathit{,}& \end{array}$where $\mathrm{|}{r}\mathrm{-}{r}{}^{\mathrm{\prime }}\mathrm{|}\mathrm{=}\sqrt{{\mathit{X}}^{\mathrm{2}}\mathrm{+}{\mathit{Y}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}$, and so the hyperpotential writes in the 3D case $\begin{array}{ccc}\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}& \mathrm{=}& \${\mathrm{\Omega }}^{\mathrm{\prime }}\mathit{\rho }\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{y}}^{\mathrm{\prime }}\mathit{,}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\right)}\\ & \mathrm{\times }& {\mathit{\kappa }}^{\mathit{xy}}\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}^{\mathrm{\prime }}\mathit{,y}\mathrm{-}{\mathit{y}}^{\mathrm{\prime }}\mathit{,z}\mathrm{-}{\mathit{z}}^{\mathrm{\prime }}\mathrm{)}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathrm{d}{\mathit{y}}^{\mathrm{\prime }}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathit{.}\end{array}$(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 $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{xy}}^{\mathrm{2}}\mathrm{(}\mathit{\rho }\mathrm{\ast }{\mathit{\kappa }}^{\mathit{xy}}\mathrm{)}& \mathrm{=}& \mathit{\rho }\mathrm{\ast }{\mathit{\partial }}_{\mathit{xy}}^{\mathrm{2}}{\mathit{\kappa }}^{\mathit{xy}}\mathrm{=}\mathit{\rho }\mathrm{\ast }\frac{\mathrm{1}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathrm{·}\end{array}$(12)This result is independent of the coordinate system, namely: $\begin{array}{ccc}\mathrm{\mathscr{H}}\mathrm{=}\{\begin{array}{c}\\ \\ \\ \end{array}& & \end{array}$(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 q₃ = 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 ${\mathit{\phi }}^{\mathrm{\prime }}\mathit{<}\frac{\mathit{\pi }}{\mathrm{2}}$), 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 κ^Rθ 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 κ^Rθ established in cylindrical coordinates and apply convenient rotations to derive κ^rφ. A more direct calculus is presented in Appendix D.

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 $\begin{array}{ccc}{\mathrm{\int }}_{\mathrm{2}\mathit{\pi }}\frac{\mathrm{d}{\mathit{\theta }}^{\mathrm{\prime }}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathrm{=}\frac{\mathrm{4}}{\mathit{\delta }}{K}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}& & \end{array}$(14)where K is the complete elliptic integral of the first kind, δ² = (R′ + R)² + (z − z′)², and $\mathit{k\delta }\mathrm{=}\mathrm{2}\sqrt{\mathit{R}{\mathit{R}}^{\mathrm{\prime }}}$. Because this function is hyperbolically singular when k → 1, we can determine an hyperpotential by an integration in the P-space. 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 closed-form for ${}^{\mathrm{\int }}\frac{\mathrm{4}}{\mathit{\delta }}{K}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathit{R}}^{\mathrm{\prime }}\mathrm{d}{\mathit{R}}^{\mathrm{\prime }}$ (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 $\begin{array}{ccc}\mathit{\psi }\mathrm{=}\frac{\mathrm{1}}{\mathit{R}}{\mathit{\partial }}_{\mathit{R}}\mathrm{\mathscr{H}}\mathit{,}& & \end{array}$(15)where $\begin{array}{ccc}\mathrm{\mathscr{H}}\mathrm{=}{\mathrm{\int }}_{{\mathrm{\Omega }}^{\mathrm{\prime }}}{\mathit{\kappa }}^{\mathit{R}}\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}\mathit{.}& & \end{array}$(16)We now present three last examples. Figure 6 is for the circular disk with radius unity. The P-grid and the P′-grid are made of N′ = 32 and N = N′ + 1 points equally spaced in R² 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 second-order. 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 $\mathrm{\Sigma }\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{R}}^{\mathrm{\prime }\mathrm{2}}}$. 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).

The last example is an inhomogeneous sphere with radius unity. The potential/density pair is the following $\begin{array}{ccc}\{\begin{array}{c}\\ & \\ & \\ \\ & \\ & \end{array}& & \end{array}$(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 $\begin{array}{ccc}\mathrm{\mathscr{H}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\sin {\mathit{\phi }}^{\mathrm{\prime }}\mathrm{d}{\mathit{\phi }}^{\mathrm{\prime }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathit{\rho }\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{\kappa }}^{\mathit{R}}{\mathit{r}}^{\mathrm{\prime }\mathrm{2}}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathit{.}& & \end{array}$(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²,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. N² 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 $\begin{array}{ccc}{A}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{\int }\frac{{u}\mathrm{d}{\mathit{q}}^{\mathrm{\prime }}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathit{,}& & \end{array}$(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, $\begin{array}{ccc}\mathit{p}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{\int }\frac{\mathrm{\nabla }\mathrm{·}\left[\mathrm{(}{u}\mathrm{·}{\nabla }\mathrm{)}{u}\right]\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathit{,}& & \end{array}$(20)where u is the fluid velocity.

---

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

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 $\begin{array}{ccc}”\frac{\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathrm{d}{\mathit{y}}^{\mathrm{\prime }}}{\sqrt{\mathrm{(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{-}\mathit{x}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{y}}^{\mathrm{\prime }}\mathrm{-}\mathit{y}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{z}}^{\mathrm{\prime }}\mathrm{-}\mathit{z}{\mathrm{)}}^{\mathrm{2}}}}\mathit{,}& & \end{array}$(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 $\mathrm{|}{r}\mathrm{-}{r}{}^{\mathrm{\prime }}\mathrm{|}\mathrm{=}\sqrt{{\mathit{X}}^{\mathrm{2}}\mathrm{+}{\mathit{Y}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}$ and dx′dy′ = dXdY. The indefinite integral is $\begin{array}{ccc}& & ”\frac{\mathrm{d}\mathit{X}\mathrm{d}\mathit{Y}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathrm{=}\mathit{X}\mathrm{-}\mathit{Y}\ln \left(\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{-}\mathit{X}\right)\mathrm{-}\mathit{X}\ln \left(\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{-}\mathit{Y}\right)\mathrm{-}\mathit{Z}\\ & & \left[\hspace{0.17em}\mathrm{atan}\frac{\mathit{X}}{\mathit{Z}}\mathrm{+}\hspace{0.17em}\mathrm{atan}\frac{\mathit{XY}}{\mathit{Z}\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\right]\mathrm{·}\end{array}$(A.2)For the present purpose, we only need to generate the Newton kernel by a mixed partial derivative ∂²/∂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 $\begin{array}{ccc}& & \int \begin{array}{}{\mathit{x}}^{\mathrm{\prime }}\mathrm{-}\mathit{x}\\ \mathrm{0}\end{array}\mathrm{d}\mathit{X}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{y}}^{\mathrm{\prime }}\mathrm{-}\mathit{y}}\frac{\mathrm{d}\mathit{Y}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathrm{=}\mathrm{-}\mathit{Y}\ln \frac{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{-}\mathit{X}}{\sqrt{{\mathit{Y}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}}\mathrm{-}\mathit{X}\\ & & \ln \frac{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{-}\mathit{Y}}{\sqrt{{\mathit{X}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}}\mathrm{-}\mathit{Z}\hspace{0.17em}\mathrm{atan}\frac{\mathit{XY}}{\mathit{Z}\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathrm{·}\end{array}$(A.3)\arraycolsep1.75ptConsequently, the hyperkernel is obtained from this expression by exchanging the variables x′ and x, and y′ and y. We find $\begin{array}{ccc}{\mathit{\kappa }}^{\mathit{xy}}& \mathrm{=}\mathit{Y}\ln \frac{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{+}\mathit{X}}{\sqrt{{\mathit{Y}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}}\mathrm{+}\mathit{X}\ln \frac{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{+}\mathit{Y}}{\sqrt{{\mathit{X}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}}\mathrm{-}\mathit{Z}\hspace{0.17em}\mathrm{atan}\frac{\mathit{XY}}{\mathit{Z}\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathrm{·}& \end{array}$(A.4)To generate κ^yz associated with a potential expressed as ∂²ℋ/∂y∂z, we perform the permutations (x,x′) ↔ (z,z′) in the formulae above. To generate κ^xz associated with a potential expressed as ∂²ℋ/∂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 $\begin{array}{ccc}{\mathit{\kappa }}^{\mathit{xz}}\mathrm{=}\mathit{Z}\ln \frac{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{+}\mathit{X}}{\sqrt{{\mathit{Y}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}}\mathrm{+}\mathit{X}\ln \frac{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{+}\mathit{Z}}{\sqrt{{\mathit{X}}^{\mathrm{2}}\mathrm{+}{\mathit{Y}}^{\mathrm{2}}}}\mathrm{-}\mathit{Y}\hspace{0.17em}\mathrm{atan}\frac{\mathit{XZ}}{\mathit{Y}\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathrm{·}& & \end{array}$(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 $\begin{array}{ccc}{\mathit{\kappa }}^{\mathit{RZ}}& \mathrm{=}& \mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}^{\mathrm{\prime }}\mathrm{)}\ln \frac{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{-}\mathrm{(}{\mathit{R}}^{\mathrm{\prime }}\mathrm{+}\mathit{R}\cos \mathrm{2}\mathit{\beta }\mathrm{)}}{\sqrt{{\mathit{R}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{2}\mathit{\beta }\mathrm{+}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}}}\\ & & \mathrm{-}\mathrm{(}{\mathit{R}}^{\mathrm{\prime }}\mathrm{+}\mathit{R}\cos \mathrm{2}\mathit{\beta }\mathrm{)}\ln \frac{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}\mathrm{+}\mathit{z}\mathrm{-}{\mathit{z}}^{\mathrm{\prime }}}{\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{R}}^{\mathrm{\prime }\mathrm{2}}\mathrm{-}{\mathrm{2}}^{\mathrm{\prime }}{\mathit{R}}^{\mathrm{\prime }}\mathit{R}\cos \mathrm{2}\mathit{\beta }}}\\ & & \end{array}$(B.2)and the potential is given by $\begin{array}{ccc}\mathit{\psi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{𝒢}{\mathit{\partial }}_{\mathit{Rz}}^{\mathrm{2}}{\mathrm{\int }}_{{\mathrm{\Omega }}^{\mathrm{\prime }}}{\mathit{\kappa }}^{\mathit{Rz}}\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}\mathit{.}& & \end{array}$(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 $\begin{array}{ccc}”\frac{{\mathit{R}}^{\mathrm{\prime }}\mathrm{d}{\mathit{R}}^{\mathrm{\prime }}\mathrm{d}{\mathit{\theta }}^{\mathrm{\prime }}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathit{,}& & \end{array}$(C.1)where $\begin{array}{ccc}\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }{\mathrm{|}}^{\mathrm{2}}\mathrm{=}\mathrm{(}{\mathit{R}}^{\mathrm{\prime }}\mathrm{+}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{4}{\mathit{R}}^{\mathrm{\prime }}\mathit{R}{\cos }^{\mathrm{2}}\left(\frac{\mathit{\theta }\mathrm{-}{\mathit{\theta }}^{\mathrm{\prime }}}{\mathrm{2}}\right)\mathrm{·}& & \end{array}$(C.2)This double integral has been calculated in Huré (2012). The indefinite form is $\begin{array}{ccc}”\frac{{\mathit{R}}^{\mathrm{\prime }}\mathrm{d}{\mathit{R}}^{\mathrm{\prime }}\mathrm{d}{\mathit{\theta }}^{\mathrm{\prime }}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}& \mathrm{=}& \mathit{\delta E}\mathrm{\left(}\mathit{\beta ,k}\mathrm{\right)}\mathrm{+}\frac{{\mathit{R}}^{\mathrm{\prime }\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{2}}}{\mathit{\delta }}\mathit{F}\mathrm{\left(}\mathit{\beta ,k}\mathrm{\right)}\mathrm{+}\frac{{\mathit{R}}^{\mathrm{\prime }}\mathrm{-}\mathit{R}}{{\mathit{R}}^{\mathrm{\prime }}\mathrm{+}\mathit{R}}\frac{{\mathit{\zeta }}^{\mathrm{2}}}{\mathit{\delta }}\\ \mathrm{\Pi }\mathrm{\left(}\mathit{\beta ,}{\mathit{m}}^{\mathrm{2}}\mathit{,k}\mathrm{\right)}& \mathrm{-}& \mathit{R}\sin \mathrm{2}\mathit{\beta }\hspace{0.17em}\mathrm{asinh}\frac{{\mathit{R}}^{\mathrm{\prime }}\mathrm{+}\mathit{R}\cos \mathrm{2}\mathit{\beta }}{\sqrt{{\mathit{\zeta }}^{\mathrm{2}}\mathrm{+}{\mathit{R}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{2}\mathit{\beta }}}\\ & \mathrm{-}& \end{array}$(C.3)where F(φ,k), E(φ,k) and Π(φ,m²,k) are the incomplete elliptic integral of the first, second, and third kinds, respectively, δ = (R′ + R)² + ζ², ζ = z − z′, $\mathit{k\delta }\mathrm{=}\mathrm{2}\sqrt{{\mathit{R}}^{\mathrm{\prime }}\mathit{R}}$, 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, ζ² and m are not impacted). We finally obtain $\begin{array}{ccc}{\mathit{\kappa }}^{\mathit{R\theta }}& \mathrm{=}& \mathit{\delta E}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathit{,k}\mathrm{\right)}\mathrm{+}\frac{{\mathit{R}}^{\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{\prime }\mathrm{2}}}{\mathit{\delta }}\mathit{F}\left({\mathit{\beta }}^{\mathrm{\prime }}\mathit{,k}\right)\mathrm{+}\frac{\mathit{R}\mathrm{-}{\mathit{R}}^{\mathrm{\prime }}}{\mathit{R}\mathrm{+}{\mathit{R}}^{\mathrm{\prime }}}\frac{{\mathit{\zeta }}^{\mathrm{2}}}{\mathit{\delta }}\mathrm{\Pi }\mathrm{(}{\mathit{\beta }}^{\mathrm{\prime }}\mathit{,}{\mathit{m}}^{\mathrm{2}}\mathit{,}{\mathit{k}}^{\mathrm{)}}\\ & \mathrm{+}& {\mathit{R}}^{\mathrm{\prime }}\sin \mathrm{2}\mathit{\beta }(\hspace{0.17em}\mathrm{asinh}\frac{\mathit{R}\mathrm{+}{\mathit{R}}^{\mathrm{\prime }}\cos \mathrm{2}\mathit{\beta }}{\sqrt{{\mathit{\zeta }}^{\mathrm{2}}\mathrm{+}{\mathit{R}}^{\mathrm{\prime }\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{2}\mathit{\beta }}}\\ & \mathrm{-}& \hspace{0.17em}\mathrm{asinh}\frac{{\mathit{R}}^{\mathrm{\prime }}\cos \mathrm{2}\mathit{\beta }}{\sqrt{{\mathit{\zeta }}^{\mathrm{2}}\mathrm{+}{\mathit{R}}^{\mathrm{\prime }\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{2}\mathit{\beta }}})\\ & \mathrm{+}& \mathit{\zeta }\{\hspace{0.17em}\mathrm{atan}\left[\frac{\mathit{\zeta }\mathrm{(}\mathit{R}\mathrm{+}{\mathit{R}}^{\mathrm{\prime }}\cos \mathrm{2}\mathit{\beta }\mathrm{)}}{\mathit{a}\sin \mathrm{2}\mathit{\beta }\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\right]\\ & \mathrm{-}& \hspace{0.17em}\mathrm{atan}\left(\frac{\mathit{\zeta }}{\sqrt{{\mathit{\zeta }}^{\mathrm{2}}\mathrm{+}{\mathit{R}}^{\mathrm{\prime }\mathrm{2}}}}\mathrm{cotan}\mathrm{2}\mathit{\beta }\right)\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{,}\end{array}$(C.4)where 2β′ = π − (θ′ − θ) = 2π − 2β. The potential is then given by $\begin{array}{ccc}\mathit{\psi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{R}}\mathrm{𝒢}{\mathit{\partial }}_{\mathit{R\theta }}^{\mathrm{2}}{\mathrm{\int }}_{{\mathrm{\Omega }}^{\mathrm{\prime }}}{\mathit{\kappa }}^{\mathit{R\theta }}\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}\mathit{.}& & \end{array}$(C.5)Under axial symmetry, we have $\begin{array}{ccc}{\mathit{\kappa }}^{\mathit{R}}\mathrm{=}\mathrm{2}\left[\mathrm{-}\mathit{\pi }\mathrm{\left|}\mathit{\zeta }\mathrm{\right|}{\mathit{ϵ}}^{\mathrm{\prime }}\mathrm{+}\mathit{\delta }{E}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}\frac{{\mathit{R}}^{\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{\prime }\mathrm{2}}}{\mathit{\delta }}{K}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}\frac{{\mathit{\zeta }}^{\mathrm{2}}}{\mathit{\delta }}\frac{\mathit{R}\mathrm{-}{\mathit{R}}^{\mathrm{\prime }}}{\mathit{R}\mathrm{+}{\mathit{R}}^{\mathrm{\prime }}}{\Pi }\mathrm{\left(}{\mathit{m}}^{\mathrm{2}}\mathit{,k}\mathrm{\right)}\right]\mathit{,}& & \end{array}$(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 $\begin{array}{ccc}”\frac{{\mathit{r}}^{\mathrm{\prime }}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathrm{d}{\mathit{\phi }}^{\mathrm{\prime }}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}\mathrm{·}& & \end{array}$(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 $\begin{array}{ccc}\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }{\mathrm{|}}^{\mathrm{2}}\mathrm{=}\mathrm{(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{+}\mathit{r\nu }{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\nu }}^{\mathrm{2}}\mathrm{)}\mathrm{-}\mathrm{4}{\mathit{r}}^{\mathrm{\prime }}\mathit{r\nu }{\sin }^{\mathrm{2}}\mathit{\tau ,}& & \end{array}$(D.2)where $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \\ \\ \end{array}& & \end{array}$(D.3)Then, we notice that the similarity between Eqs. (C.1) and (D.1) is perfect if we make the following substitutions $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \\ \\ \end{array}& & \end{array}$(D.4)and then $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(D.5)From Eq. (C.3), we see that Eq. (D.1) becomes $\begin{array}{ccc}”\frac{{\mathit{r}}^{\mathrm{\prime }}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathrm{d}{\mathit{\phi }}^{\mathrm{\prime }}}{\mathrm{|}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{|}}& \mathrm{=}& \mathit{\delta E}\mathrm{\left(}\mathit{\tau ,k}\mathrm{\right)}\mathrm{+}\frac{{{\mathit{r}}^{\mathrm{\prime }}}^{\mathrm{2}}\mathrm{-}\mathrm{(}\mathit{r\nu }{\mathrm{)}}^{\mathrm{2}}}{\mathit{\delta }}\mathit{F}\mathrm{\left(}\mathit{\tau ,k}\mathrm{\right)}\\ & \mathrm{+}& \frac{{\mathit{r}}^{\mathrm{\prime }}\mathrm{-}\mathit{r\nu }}{{\mathit{r}}^{\mathrm{\prime }}\mathrm{+}\mathit{r\nu }}\frac{{\mathit{r}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\nu }}^{\mathrm{2}}\mathrm{)}}{\mathit{\delta }}\mathrm{\Pi }\mathrm{\left(}\mathit{\tau ,}{\mathit{m}}^{\mathrm{2}}\mathit{,k}\mathrm{\right)}\\ & \mathrm{-}& \mathit{r\nu }\sin \mathrm{2}\mathit{\tau }\hspace{0.17em}\mathrm{asinh}\frac{{\mathit{r}}^{\mathrm{\prime }}\mathrm{+}\mathit{r\nu }\cos \mathrm{2}\mathit{\tau }}{\mathit{r}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\nu }}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathrm{2}\mathit{\tau }}}\\ & \mathrm{-}& \end{array}$(D.6)We then derive ${\mathit{\kappa }}^{\mathit{r\phi }}\mathrm{=}{}^{!}\frac{\mathit{r}\mathrm{d}\mathit{r}\mathrm{d}\mathit{\phi }}{\mathrm{|}{r}\mathrm{-}{r}{}^{\mathrm{\prime }}\mathrm{|}}$ from Eq. (C.4) by making the same substitutions, and the potential is given by $\begin{array}{ccc}\mathit{\psi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{r}}\mathrm{𝒢}{\mathit{\partial }}_{\mathit{r\phi }}^{\mathrm{2}}{\mathrm{\int }}_{{\mathrm{\Omega }}^{\mathrm{\prime }}}{\mathit{\kappa }}^{\mathit{r\phi }}\mathrm{d}{\mathit{m}}^{\mathrm{\prime }}\mathrm{·}& & \end{array}$(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 $\begin{array}{ccc}\mathrm{\mathscr{H}}& \mathrm{=}& \mathrm{-}\mathrm{𝒢}{\mathrm{\int }}_{{\mathit{r}}^{\mathrm{\prime }}}\mathit{\rho }\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{r}}^{\mathrm{\prime }\mathrm{2}}{\mathit{\kappa }}^{\mathrm{shell}}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\\ & \mathrm{=}& \mathrm{-}\mathrm{4}\mathit{\pi }\mathrm{𝒢}\mathrm{\int }\mathit{\rho }\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{r}}^{\mathrm{\prime }}\mathit{H}\mathrm{(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{-}\mathit{r}\mathrm{)}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathrm{-}\frac{\mathrm{4}\mathit{\pi }\mathrm{𝒢}}{\mathit{r}}\mathrm{\int }\mathit{\rho }\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{r}}^{\mathrm{\prime }\mathrm{2}}\mathit{H}\mathrm{(}\mathit{r}\mathrm{-}{\mathit{r}}^{\mathrm{\prime }}\mathrm{)}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathit{,}\end{array}$(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

All Figures