EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 541 (May 2012) A&A, 541 (2012) A130 Full HTML
Free Access
Issue
A&A
Volume 541, May 2012
Article Number A130
Number of page(s) 6
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201118443
Published online 14 May 2012

© ESO, 2012

1. Introduction

As a direct consequence of Newton’s law for gravitation (Newton 1760; Kellogg 1929), the potential of any continuous distribution of matter inside a volume at a point P(r) of space is given by (1)where ρ(r′) is the mass density at , and d3τ is the elementary volume1. In general, this is a three-dimensional (3D) converging integral. The presence of the Green kernel 1/|r − r′| is known to represent a difficulty in calculating ψ everywhere inside and very close to since this function diverges as r → r′. This singularity is classically avoided by converting the Green function into an infinite series (e.g. Kellogg 1929; Cohl & Tohline 1999). Although exact, series expansions suffer from a low convergence rate since these are alternating series; besides, the number of integrals increases linearly with the number of terms considered up to the truncation order. These problems collectively constitute a real practical difficulty (Clement 1974; Stone & Norman 1992). The proper treatment of the singularity is the subject of a longstanding challenge. Shifting the P-grid and the P′-grid relative to each other or raising the numerical resolution around the singularity are the most natural techniques (Stemwedel et al. 1990), but these are of limited efficiency. When possible, the separate treatment of the asymptotic form of the singularity gives very good results (e.g. Ansorg et al. 2003; Huré 2005), although this approach renders the global treatment somewhat complex. The difficulty can also be tackled by introducing a “softening length” (Hockney & Eastwood 1988; Adams et al. 1989). This recipe – widespread in disc simulations – must however be seen as nothing but a crude approximation that cannot be used for accurate modeling under a certain scale. Whatever the prescription for the softening length, which is generally linked to the resolution or smallest physical length scale (e.g. Huré & Pierens 2009), the inferred force field is globally weaker than in the Newtonian case, and the evolution and stability of gaseous systems is inevitably impacted in a non-trivial manner (Romeo 1998; Sommer-Larsen et al. 1998; Adams et al. 1989; El-Zant 1998). The Poisson equation of course provides another other way to derive ψ numerically (Kellogg 1929; Durand 1953). This approach requires accurate boundary or interior/matching conditions only accessible through Eq. (1), and complex geometries are not always easy to manage (Grandclément et al. 2001).

In this paper, we present a new means to evaluate Eq. (1) that avoids the singularity, and, at the same time, properly accounts for it. This is achieved by replacing the singular kernel by an equivalent and regular, two-term form, one term being the gradient of a new scalar potential2. This is the aim of Sects. 2 and 3. Kernel equivalence is fundamental to preserving the Newtonian character of the interaction on all scales, and for any separation (in particular at long-range). This reformulation is designed to be efficient within sources, and is not expected to surpass usual methods outside sources. From a practical point of view, the potential easily becomes accessible as diverging kernels have disappeared from the volume integrals.

Although not specific to a given system of coordinates, the calculus is developed in cylindrical coordinates for which the equivalent kernel takes a nominal form, in particular under axial symmetry (this equivalent kernel is probably not unique; see Sect. 4). It can be applied as it is to all rotating gaseous/solid bodies (stars, discs, planets, asteroids, etc.) in either steady state or not, and for various applications (Dermott 1979; Hachisu 1986b; Baruteau & Masset 2008). A basic, six-step algorithm is reported in Sect. 5, together with a numerical experiment using the most simple quadrature and differentiation rules. There is no special assumption about the distribution of matter in space (density field and geometry or shape), making the method general, and transposable to domains of physics other than gravitation. Some interesting perspectives are listed in the last section.

thumbnail Fig. 1

Typical configuration for the gravitating, celestial body, and associated notations. See note 3 for the definition of points P′′ and Q.

2. Splitting of the Green kernel

thumbnail Fig. 2

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

We consider a volume of space continuously filled with matter, as depicted in Fig. 1. Using cylindrical coordinates, with P′(a,θ′,z) referring to source points, and P(R,θ,Z) to space points, the above integral for the Newtonian potential becomes (2)where d3τ = adadθ′dz is the elementary volume, is the relative altitude, and (5)Although Δ → 0 inside , the potential is generally a finite quantity (i.e. integration is a regularizing process; e.g. Kellogg 1929; Durand 1953). We now set (6)This quantity is finite everywhere, and non-zero except on the polar axis (see below). Assuming R > 0, we have (7)and so (8)where we have defined (9)The Green kernel is then split into two terms. The term is always regular3; we call this term the cool kernel in the following. Figure 2 displays R/Δ and R/Δ⋆ versus a/R around the singularity. We clearly see that the amplitude of Δ⋆ is bounded, in contrast to Δ. The second term in Eq. (8) is still singular when Δ → 0, but its integration over the material volume is expected to produce a regular field. The idea is to generate this singular kernel from the gradient of a regular function κ (hereafter called hyperkernel), and to integrate it over the volume . As the two spaces (R,θ,Z) and (a,θ′,z) are independent, the derivative may be drawn before the integral. This reasoning can be summarized as where the gradient ∇ is to be taken with respect to one of the three variables R, θ, or Z. There are then three possible hyperkernels.

The existence of the hyperkernel is not guaranteed a priori, and it is of interest only if it is available in a closed form. The investigation indeed shows that the nominal form is obtained by considering the vertical gradient (i.e. ∇ ≡ ∂Z). This may be due to the special role that the Z-axis plays a in cylindrical coordinates. We therefore do not discuss in detail any of the other two options, although these might be useful in certain circumstances.

3. The singular term as the vertical gradient of a hyperkernel

To get the hyperkernel, we consider the integration of the singular term in Eq. (8) with respect to Z. By using the intermediate variable t = ζ/Δ, we find after some algebra (11)where (12)with 0 ≤ m ≤ 1. We finally get (13)We could obviously add to κ any function of a, R, θ and θ′, but this is not necessary here as we take its Z-gradient. The Green kernel is then given by the equivalent form (14)It is necessary to verify that κ is regular. This is straightforward since Δ > |ζ| as soon as R > 0. In other words, if m → 1 and , then ζ/Δ →  ± 1.

thumbnail Fig. 3

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

If we now multiply Eq. (8) by ρ(a,θ′,z)d3τ – which does not depend on Z – and integrate over the material volume , we find that (15)where the partial derivative now operates on the integral. Up to a factor  − G, this expression is precisely the potential defined by Eq. (1), and it is both exact and general. It depends on neither the body’s shape nor on the distribution of its mass density. It applies not only to volume distributions, but also to surface distributions (see Sect. 6) and linear distributions. The first integral in Eq. (15) is then the cool potential associated with the cool kernel, and the second term is the vertical gradient of a hyperpotential.

On the polar axis (i.e. R = 0) Δ = δ, and so Eq. (8) does not help us to treat the singularity when ζ = 0 and a = 0. In this case, we have (16)provided that a > 0. The potential can then be written (17)where here there is no cool kernel. Besides, we see that lima → 0   aκ    =    0.

There is no continuity between the two different expressions for the cool kernel, the one valid at R = 0 and the other valid at R → 0 (this is also true for the hyperkernel). This is no problem as long as we do not have to consider the radial gradient of κ (see below the numerical experiment).

Finally, we note that, as we work with an equivalent form of the Green kernel, the potential found from Eq. (15) automatically has the right asymptotic property, and varies like M/r sufficiently far away from the body. At large relative distance (i.e. R ≫ a and Z ≫ z), we have . At the lowest order, one finds that (18)which means that the long-range behavior is exclusively contained in the cool kernel.

4. The case of axially-symmetric bodies

Axially-symmetric bodies constitute an important class of astrophysical objects (Chandrasekhar 1973; Hachisu 1986a). Interestingly enough, in problems where ∂θ′ρ = 0, we can rewrite the above expressions in a more compact form in terms of elliptic integrals. The first integral in the right-hand-side of Eq. (15) becomes4(19)where (20)is the complete elliptic integral of the second kind, is the modulus (with k ∈ [0,1] ), and the double integration runs over the meridional cross section of the body. The integration over the polar angle θ′ of the hyperkernel gives: (21)where (22)is the incomplete elliptic integral of the first kind, and (23)is the incomplete elliptic integral of the third kind (m is the parameter and ). Over the whole circle (i.e. ), this yields the axially symmetric potential5: (25)where H is defined for convenience by (26) with and . The presence of the H-function is actually expected here since ∂zΔ and ∂ZΔ are linked (Trova et al. 2012).

On the polar axis, we get (27) where κ is, in this case, given by Eq. (16). We note that, in this axially symmetric case, the potential could be determined through Eq. (19) (i.e. by using the cool kernel only, and no hyperkernel), but the integrand still contains a hyperbolic divergence as a → 0 and ζ → 0 which is not easy to manage. This is why it seems much better to consider Eq. (27), as the logarithmic divergence of the hyperkernel (i.e. the    asinh term) is cancelled out by the elementary volume when a is close to 0 (i.e. lima → 0   aκ    =    0).

5. An example of implementation

We now present a first numerical experiment to briefly describe the main steps of the method, and demonstrate its simple implementation. For a full 3D body, the main steps6 of any numerical estimate can be summarized as

  • 1.

    discretize the material source on a 3D-grid, in the form of quadruplets ;

  • 2.

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

  • 3.

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

  • 4.

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

  • 5.

    determine the vertical gradient of the hyperpotential;

  • 6.

    add this derivative to the cool potential, and multiply by  − G, and reiterate steps 2 to 6 to generate a potential map. The components of the associated gravitational force are deduced as usual from the three gradients of ψ.

We have considered a homogeneous, axially symmetric torus with a square meridional cross section with . The mass density (ρ = 1 inside and 0 outside) is defined on a regular mesh, with N nodes in each direction. The computational grid also consists of a square box (R,Z) ∈ [0,2]  × [−1,1]  with L nodes per direction and regular spacing, therefore encompassing the body. The double integrals in Eqs. (25) and (27) are computed at each node of the computational grid through the two-point trapezoidal rule. The partial derivative of the hyperpotential is estimated through second-order finite differences. These basic schemes are very easy to implement and the above six-step procedure contains no pitfall. We take N = L = 31 here. Figures 3a–d display, respectively, the integral of the cool kernel, the integral of the hyperkernel κ, its vertical gradient, and the total potential, obtained by adding the first and third maps. The boundary of the toroidal body is superimposed on these maps. We note that the vertical gradient of the hyperpotential makes the geometrical cross section rise above the background, and filters the curvature effects (which are enhanced by the integral of the cool kernel).

The computing time is typical of integral methods. Under the conditions of the present example, the integration of the cool kernel and hyperkernel (step 4 above) requires 2N2 elementary operations per point of space. To get the potential at the nodes of a L × L square grid, we then find 2N2L2 (the time needed to determine the vertical gradients is negligible in comparison). This is therefore much smaller than that is usually obtained from an expanded Green function by a factor equal to the number of terms up to the truncation order (which can be as high as a few hundred; see e.g. Hachisu 1986a; Stone & Norman 1992).

thumbnail Fig. 4

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

6. Checking the accuracy

The numerical accuracy is sensitive to various ingredients, such as the quadrature and differentiation schemes. To a lesser extent, it also depends on the mass density distribution ρ and equation of the boundary , which may generate additional difficulties in the calculations (interpolation of data points, infinite derivatives at edges, etc.). We present a second numerical test illustrating the accuracy of the method by considering an exact potential/density pair. Not many configurations correspond to finite mass and finite size systems (e.g. Binney & Tremaine 1987). When matter is gathered in a plane (i.e. a flat disc), Eq. (15) is directly transposable by setting ρ(a,θ′,z) = Σ(a,θ′)δ(z) and integrating over z. Under axial symmetry, we respectively get7 from Eqs. (25) and (27) (29)for R > 0, and (30)where ain and aout denote the discs inner and outer edges, and (31)By setting ain = 0 and for a ∈ [0,aout] , one gets one of the three cases analyzed by Schulz (2009). For such a distribution, the associated potential, ψe, is known exactly in a closed-form for any point of space. Figure 4 gives ψe as well as the error index ϵ = log 10|1 − ψ/ψe| where ψ is determined from Eqs. (29) and (30) in the same conditions as above (we set aout = 1). We see that the potential outside and especially inside the disc is well-reproduced. The relative error, on the order of 10-3, agrees with the second-order of the schemes at the actual mesh size of . The accuracy can be tuned by changing the quadrature and differentation schemes.

7. Concluding remarks

We have reformulated the Green kernel appearing in potential problems to circumvent the singularity and, at the same time, properly account for it. As a consequence, the gravitational potential of any celestial body, regardless of its shape and matter density distribution, becomes directly accessible through two “classical” volume integrals, followed by a partial derivative. The method is applicable to 3D, fully inhomogeneous systems, as well as to surface and line distributions. It is especially efficient inside distributions where most approaches exhibit a real practical complexity, converge very slowly, or produce spurious errors. The presence of regular kernels ensures that the method is stable and easy to implement. This should encourage modellers to abandon various integration techniques that do not “faithfully” reproduce the Newtonian character of the potential and forces. In the context of discs for instance, this method appears to be a real alternative to softened gravity, which remains a free parameter, non-Newtonian theory. As stressed, it is probably possible to determine other cool kernel/hyperkernel pairs (for instance, by considering the hyperkernel as a radial/angular gradient), but the one presented in the body of this paper seems the simplest one. It is in particular well-suited to axially symmetric configurations.

This study needs to be continued in several respects, including the analysis of the mathematical properties of the cool kernel and hyperkernel and their physical meanings, as well as the derivation of the equivalent kernel in other systems of coordinates (e.g. Cartesian and spherical coordinate systems). In addition, it would be interesting to expand the two kernels in Eq. (15) in series, and compare their properties with the expansion of the Green function in Legendre polynomials, inside as well as outside the body. The cool kernel/hyperkernel pair is also interesting as a new starting point to generating various kinds of approximations. Apart from the astrophysical context where there are so many applications about gravitation, this technique is also transposable to other kinds of problems involving improper integrals. This can be, for instance in electromagnetism, the determination of the potential vector A and associated magnetic field induced by current densities (Jackson 1998; Cohl & Tohline 1999). These points will be touched on in forthcoming papers.

1

This form holds in electrostatics, where ρ is the density of electric charges, and the constant is (instead of  − G).

2

Some aspects of the theory of tensor potentials are developped by Chandrasekhar (1973) in the context of rotating, self-gravitating fluids.

3

It is the inverse of a distance, and its value can be interpreted geometrically by noting that (see Fig. 1) (10) where Q(R,θ′ + π,Z) is a point of space, diametrically opposite to P′, and P′′(a,θ + π,z) is a point, diametrically opposite to P, that belongs to a fictitious source.

4

A factor of two is due to dθ′/dφ, and another factor of two contained in the modulus k comes from symmetry consideration (i.e. matter located at θ′ ∈ [θ,θ + π]  provides the same contribution as matter located at θ′ ∈ [θ,θ − π]).

5

If we perform the Z-derivative and rearrange terms, we recover the well-known expression (Durand 1953), namely (24) whose kernel is logarithmically singular.

6

Step 1 should be executed once and for all, except for systems that evolve with time.

7

These two formulae can be compared with the classical expression (Durand 1953; Binney & Tremaine 1987): (28) which is singular inside the source.

Acknowledgments

It is a pleasure to thank J. Braine, M. Gazeau, F. Hersant, and A. Trova as well as the second referee for comments and suggestions about infinitely flat systems. J.-M. Huré is grateful to the CNU, Sect. 34, for supporting a six-months full-time research project through CRCT-2011 funding delivered by the MESR.

References

  1. Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959 [NASA ADS] [CrossRef] [Google Scholar]
  2. Ansorg, M., Kleinwächter, A., & Meinel, R. 2003, MNRAS, 339, 515 [NASA ADS] [CrossRef] [Google Scholar]
  3. Baruteau, C., & Masset, F. 2008, ApJ, 678, 483 [NASA ADS] [CrossRef] [Google Scholar]
  4. Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton, NJ: Princeton University Press, 747) [Google Scholar]
  5. Chandrasekhar, S. 1973, Ellipsoidal figures of equilibrium [Google Scholar]
  6. Clement, M. J. 1974, ApJ, 194, 709 [NASA ADS] [CrossRef] [Google Scholar]
  7. Cohl, H. S., & Tohline, J. E. 1999, ApJ, 527, 86 [NASA ADS] [CrossRef] [Google Scholar]
  8. Dermott, S. F. 1979, Icarus, 37, 575 [NASA ADS] [CrossRef] [Google Scholar]
  9. Durand, E. 1953, Electrostatique, Vol. I, Les distributions (Masson) [Google Scholar]
  10. El-Zant, A. A. 1998, A&A, 331, 782 [NASA ADS] [Google Scholar]
  11. Grandclément, P., Bonazzola, S., Gourgoulhon, E., & Marck, J.-A. 2001, J. Comp. Phys., 170, 231 [NASA ADS] [CrossRef] [Google Scholar]
  12. Hachisu, I. 1986a, ApJS, 61, 479 [NASA ADS] [CrossRef] [Google Scholar]
  13. Hachisu, I. 1986b, ApJS, 62, 461 [NASA ADS] [CrossRef] [Google Scholar]
  14. Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation Using Particles, ed. R. W. Hockney, & J. W. Eastwood [Google Scholar]
  15. Huré, J.-M. 2005, A&A, 434, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Huré, J., & Pierens, A. 2009, A&A, 507, 573 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Jackson, J. D. 1998, Classical Electrodynamics, 3rd edn. [Google Scholar]
  18. Kellogg, O. D. 1929, Foundations of Potential Theory (New-York: Frederick Ungar Publishing Company) [Google Scholar]
  19. Newton, I. 1760 [Google Scholar]
  20. Romeo, A. B. 1998, A&A, 335, 922 [NASA ADS] [Google Scholar]
  21. Schulz, E. 2009, ApJ, 693, 1310 [NASA ADS] [CrossRef] [Google Scholar]
  22. Sommer-Larsen, J., Vedel, H., & Hellsten, U. 1998, MNRAS, 294, 485 [NASA ADS] [CrossRef] [Google Scholar]
  23. Stemwedel, S. W., Yuan, C., & Cassen, P. 1990, ApJ, 351, 206 [NASA ADS] [CrossRef] [Google Scholar]
  24. Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [Google Scholar]
  25. Trova, A., Huré, J. M., & Hersant, F. 2012, MNRAS, submitted [Google Scholar]

All Figures

thumbnail Fig. 1

Typical configuration for the gravitating, celestial body, and associated notations. See note 3 for the definition of points P′′ and Q.
In the text
thumbnail Fig. 2

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

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

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

Current usage metrics show cumulative count of Article Views (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.

Initial download of the metrics may take a while.