A&A 434, 1-15 (2005)
1 - LUTh (CNRS UMR 8102), Observatoire de Paris-Meudon, Place Jules Janssen, 92195 Meudon Cedex, France
2 - Université Paris 7 Denis Diderot, 2 place Jussieu, 75251 Paris Cedex 05, France
Received 14 August 2003 / Accepted 10 December 2004
In a series of two papers, we present numerical, integral-based methods to compute accurately the self-gravitating field and potential induced by tri-dimensional, axially symmetric fluids, with a special regard for tori, discs and rings. This first article is concerned with a fully numerical approach. Complex shapes, small/large aspect ratios, important density gradients and compact/extended systems can be accounted for. Loop singularities in the Poisson integrals are carefully treated from kernel splitting and/or density splitting. Field components are obtained from density splitting: the local density field is separated into a vertically homogeneous contribution for which the integrable singularity is known in a closed form, plus a "residual'' contribution which yields a regular integrand. This technique is exact in the vertically homogeneous limit. The potential is computed from double splitting: kernel splitting (to isolate explicitly the singular function), followed by density splitting. In each two directions, numerical quadratures are performed using a dynamical mesh (very efficient for flat and extended systems), combined with a high-order, irregular-spacing scheme. Global performances are demonstrated through a few test-configurations. The accuracy is potentially very high, and can reach the computer precision (for simple geometries and smooth density profiles). In terms of computing time to precision ratio, present methods are asymptotically more efficient than classical finite-difference methods. These can be employed either as a "Poisson-solver'' or only to compute a sub-set of Dirichlet/Neumann-type boundary conditions in order to initialize spectral/finite-difference/finite element methods. As an example of an astrophysical application, we determine the structure of an uniformly rotating polytrope belonging to the compressible, Dyson-Wong sequence. In a second paper (Paper II), we show that singularities are integrable analytically for geometrically thin systems where the density is of the form . Improvements and extensions (such as the non-axi-symmetric case) are possible and discussed.
Key words: hydrodynamics - methods: numerical
Reliable models and simulations of self-gravitating systems like the cosmological fluid, galaxies, stars, interstellar clouds, and tori, discs and rings need accurate solutions of the Poisson equation
There is obviously a wide astrophysical literature related to this topic, but most publications i) are of theoretical nature; ii) concern bi-dimensional distributions where matter is located in the plane (i.e. zero-thickness or so-called "razor'' thin objects); and iii) are of limited application or/and limited precision. In comparison, there is indeed little numerical work (Hachisu 1986; Bodo & Curir 1992; Störzer 1993; Cohl 1999; Ansorg et al. 2003b). In a series of two papers, we present powerful methods to determine numerically the field and the potential induced by axially symmetric, tri-dimensional mass distributions inside sources, including source boundaries (with a straightforward extension to non-axial symmetry). So, their use is twofold: i) as a Poisson-solver; and ii) as a boundary condition-solver. These methods are made for continuous media, not for particles. No disc model will be produced: we are interested only on potential and field values given a mass density in space. Note that for most applications, the gravitating field
In Sect. 2, we recall the exact expressions for the field induced by an isolated, homogeneous loop of matter, and by a tri-dimensional system as well. Additional notation, definitions and assumptions are introduced. In Sect. 3, we explain the density splitting method in order to estimate accurately integrable singularities occurring when in Eq. (4). We then discuss the dynamical sampling of source points along the vertical direction, and demonstrate the reliability of this technique through a few simple examples. Section 4 deals with the integration of kernels along the radial direction. Sampling of source points is again discussed, and the accuracy of the global method is checked. The treatment for the potential is explained in Sect. 5. Tests are reported. A classical example of astrophysical application of this Poisson-solver is given. The expected radial and vertical physical resolutions are estimated in Sect. 6. The efficiency of the method in terms of computing time to precision ratio is discussed and compared to that of finite-difference methods. Concluding remarks are found in the last section.
In a second paper (Pierens & Huré 2005; hereafter Paper II), we show that integrable singularities can be estimated analytically in the case of systems with small aspect ratio, namely slim discs, thin discs and rings which form an important sub-class of toroidal configurations.
|Figure 1: The method is based on the exact expressions for the field and potential induced by a homogeneous loop of matter (see Eqs. (5)-(8)). Poisson kernels are regular everywhere in space except on the loop where they diverge. This singularity still remains in two and three dimensions, but it is integrable.|
|Open with DEXTER|
We use polar cylindrical coordinates throughout. We call "field point'' a point where the field is requested. We call "source point'' a point which belongs to the matter distribution. To avoid any confusion, we employ two different coordinate systems as usual in potential theory, namely for field points, and for source points, even though field points and source points can be confused.
An infinitely thin (i.e. topologically one dimensional) homogeneous, circular loop with radius a, linear mass
and located at altitude z as pictured in Fig. 1, generates at any field point a gravity field
|Figure 2: Complete elliptic integrals K and E, and the "regular'' complete integral of the first kind defined by Eq. (38), versus the modulus k.|
|Open with DEXTER|
Under these conditions, the self-gravitating field due to any axi-symmetric, tri-dimensional system with mass
is given by
|Figure 3: The axially symmetric system ( left) is divided into N cylinders. One of these cylinders, named the "coincident cylinder'', contains the actual field point. On this coincident cylinder ( right), one of the source points, named the "coincident point'', coincides with the field point. For the coincident point, k=1 by construction and the kernels are singular.|
|Open with DEXTER|
Before going into the details of the computational method, we introduce a few additional notations and assumptions. A typical configuration is shown in Fig. 3: denotes the inner edge of the system, is the outer edge, and is the equation of the surface. Inside the fluid, the mass density can be given either in the form of an analytical expression , or as a grid by a two entry table . If, as will be proposed here, a dynamical mesh is implemented (see Sect. 3.5 and Sect. 4.1), interpolations in this density grid will be necessary. This means that the selection of source points (an,zm) is not known in advance but varies from one field point to another.
The shape (and subsequently the aspect ratio H/a) of the system is not imposed, that is H(a) can take any value. The mass density at the boundary is not constrained either, except that it must vanish continuously at the surface, i.e. with , as expected in any real system. All these assumptions are quite natural for astrophysical media, making the method of general use. In particular, it applies to thick discs, tori, thin discs and rings .
Double integrations in Eq. (10) can be performed in two different orders. For toroidal configurations, integrating first following the z-direction, and second following the a-direction, seems the most appropriate. In other words, given the notations introduced above, Eq. (10) reads
The exact expression for Eq. (12) is known for a vertically homogeneous distribution, that is for
Actually, we have for any field point (Durand 1964)
At the coincident point, the kernel
diverges. So does the integrand
in Eq. (12). In this case, we re-write the density in the form
|Figure 4: Illustration of the density splitting method. The cylinder a) has radius a=1, semi-thickness H=1, and the mass density follows a parabolic law b). The field point has arbitrarily coordinates R=a and . There, both integrands and diverge c).|
|Open with DEXTER|
Let us now illustrate the method step by step, graphically. Figure 4a shows a cylinder with radius a=1 and total height 2H=2, carrying a parabolic density profile that satisfies the assumptions listed in Sect. 2.2. This profile is plotted in Fig. 4b. As field point, we choose arbitrarily the point with coordinates R=a and . There, both integrands and diverge because of the loop singularity. These two quantities are displayed in Fig. 4c. At the field point, the density takes the value . The residual density is plotted in Fig. 5a. The "field'' due to a cylinder with constant density (that is ) is known in a closed form. So, the contribution due to remains to be computed numerically. Components of the residual integrand are shown in Fig. 5b. We see that these are regular functions, in particular at the field point. Subsequently, standard quadrature schemes can furnish the residual secondary kernel , and the total secondary kernel is found from Eq. (22).
|Figure 5: Illustration of the density splitting method (continued). By construction a), the "residual'' profile vanishes at the field point, where the mass density at the field point . The "residual'' integrands b) are finite, and can be integrated numerically with standard schemes. In the homogeneous limit, the method is exact.|
|Open with DEXTER|
We now have a problem of numerical quadrature. As usual, efficient rules that minimize the amount of source points and maximize the precision should be preferred. Several techniques can be considered, such as local quadrature rules or collocation methods. Here, we use a sampling function, denoted w(z), in order to enhance the number of source points around the coincident point in combination with the high-order, irregular step quadrature rule of Gill & Miller (1972).
Many choices are possible for w(z). Preliminary tests have shown that an appropriate function is
Table 1: Two mass density/shape models considered in this study. For model 9 1, an exact expression for both and gZ can be derived analytically from Eqs. (8) and (9) provided (see the Appendix B, Eq. (B.2); the kernel and field are finite for ). Model 2 is a vertically parabolic profile. In this case, there is no solution in a closed form, whatever the shape of the system. Note that for whereas for .
|Figure 6: Top panel: error index versus the number of vertical source points M, obtained for a cylinder with radius a=1 and density profile (see Table 1). Four semi-thicknesses are considered: H=10-4, 0.01, 0.1 and 1. The field point has coordinates R=a, and . Bottom panel: error index obtained in the same conditions as above for h=1 ( plane line) compared with density/shape model 2 ( symbols).|
|Open with DEXTER|
To check the density splitting method at this intermediate level (i.e. before performing the radial integration) as well as its accuracy, we need vertical test profiles for which the corresponding secondary kernel is exactly known. For a proper and fair check, it is necessary that test-densities do not vanish just at the coincident point. This would imply that there is no kernel singularity, and thus no homogeneous contribution (i.e.
). We shall restrict ourselves to the vertical component because of the simple form of the kernel
that enables us to easily find "density/secondary kernel'' pairs. Also, to shorten the discussion, we report a single test, although we have successfully checked many. We give in Table 1 two models for the density profile and shape. Let us first consider model 1. We can thus compare the reference value corresponding to ,
to values obtained from the density splitting method. In the following, comparisons are made with the "error index''
Figure 6b shows the results obtained in the same conditions as above but for the density/shape model 2 reported in Table 1 which corresponds to a parabolic vertical profile, an astrophysically more realistic case. The error index for both the radial and vertical components is plotted. Since exact values of the secondary kernels are not known in this case, we take as references and , values computed for very large M such that all the available digits are fixed and stationary (because of the order of the present quadrature scheme, doubling the resolution would suffice). We see that there is almost no difference between models 1 and 2. We have checked many density profiles and have noticed that the performances of the splitting method indeed hold. Applying the splitting method outside the cylinder, that is for (and ), gives much better results (i.e. is reached for fewer source points), unsurprisingly.
|Figure 7: Typical shape of secondary kernels ( left panels) and ( right panels) versus a at two different scales for a flat disc, with inner edge , outer edge and thickness 2H=2. The area under each curve ( shaded zones) gives the field components. The location of the N cylinders as specified by the v-function is marked by circles. The radial extension of the peaks is a few local thicknesses. Units are arbitrary.|
|Open with DEXTER|
The radial integration of the secondary kernel according to Eq. (11) is the second and last step. This operation is less tricky since it has regular components. We shall proceed as in the previous section, through a sampling function in order to set the position of the N cylinders in the range in an adequate manner (see Fig. 3), and using the same quadrature rule. Let us call v(a) this function. For systems with (i.e. for compact systems like tori), we could employ a formula similar to w(z). However, such a choice is not appropriate for general purposes, especially for geometrically thin and extended systems. We stress that the self-gravitating field is determined by two contributions:
This is illustrated in Fig. 7. Both the local and long-range contributions can have comparable magnitude. As a consequence, the v-function must be such that the sampling of cylinders is good enough not only in the neighborhood of the coincident cylinder (i.e. for
but also far from it. This is achieved with the following piece-wise sampling-function
In practice, the value is well suited for general use. We impose only the number of cylinders populating the local domain, and then determine the number of cylinders belonging to the left-side and right-side long-range domains (denoted and respectively) in such a way that v has continuous derivatives at the long-range/local connection. Thus, N is an output quantity. The detail of these settings is found in Appendix C.
We have checked the method for various tri-dimensional configurations. Here are a few typical examples involving the three systems listed in Table 2. Systems 1 and 2 are perfectly flat objects, whereas system 3 is a flared, extended disc. Figure 8 shows the error index (see Eq. (28)) and the number N of cylinders for systems 1 and 2 with density/shape model 1 (see Table 1). In this case, the reference field is known. Computations have been performed for M=8, 16, 32, 64 and 128 vertical source points (for all cylinders). We see that the accuracy is fully tunable, and can be extremely high. For system 1, a hundred points in each direction yield a relative precision of 10-11 on average (this would be 10-2 using a 3-point, finite difference scheme). Due to the present quadrature rule, the relative precision is improved by two orders of magnitude when both N and M are doubled. The more compact the system, the more efficient the method for a given grid size . Again, the remarkable issue is that, with a few points only, the relative precision is already good. For instance for system 1, a grid yields (this is achieved with a size grid for system 2).
Table 2: Three systems to check the full method.
|Figure 8: Number N of cylinders ( top panels) and error index on the vertical field ( middle panels) for the density/shape model 1 (see Table 1) and systems 1 ( left) and 2 ( right) for which characteristics are listed in Table 2. The altitude of the field point is . The vertical field is also given ( bottom panels).|
|Open with DEXTER|
Figure 9 displays the results obtained in the same conditions but for system 3 and for model 2 for the density and shape. This case is typical of a certain class of astrophysical discs. As in Sect. 3.6, unknown reference values for gR and gZ have been estimated using a larger number of source points. We see that the global trends are identical, suggesting the "universal'' character of the method performance. In general, we find that the accuracy is mainly sensitive to M (and weakly to N, thanks to the dynamical mesh), with roughly
|Figure 9: Same legend and same conditions as for Fig. 8 (except ) but for a geometrically thin and flared disc (system 3) with density/shape model 2 (see Table 1).|
|Open with DEXTER|
Let us now consider the potential. Following the formalism employed in Sect. 2, the potential
due to an infinitely thin, homogeneous circular loop (see Fig. 1) with linear mass
and radius a is (e.g. Durand 1964)
|Figure 10: Same legend and same conditions as for Figs. 8 and 9 but for the potential (see text for more details).|
|Open with DEXTER|
We see from Eq. (35) that the straightforward determination of the secondary kernel
is problematic because of the loop singularity.
In order to compute this quantity everywhere and especially at the coincident point, we shall use a slightly different approach than for the field. Actually, we failed to determine an analytical expression for a density field, non-vanishing at the coincident point and such that Eq. (37) is known in a closed form. The choice of a z-independent mass density
as considered in Sect. 3.2 does not help here. Since the asymptotic form of the K-function as
is known, we first set
We can now use the analogue of the "density splitting'' method to compute
and split this quantity into two terms according to
We give here a few examples to prove the efficiency of the double splitting method for the potential. To abridge the discussion, we shall not report tests made on the secondary kernel (relative precisions are similar to Sect. 3.6). Figure 10 is the analogue of Figs. 8 and 9: it displays the error index obtained for systems 1, 2 and 3 as already considered before (see Table 2). For these three systems, the "exact'' potential being not known, we take as reference potentials , values determined by enhancing the numerical resolutions. Graphs show that, in all cases, the relative precision on the potential is excellent with a few source points. Issues raised for the field components hold here. In particular, Eq. (33) applies, with the same weak dependence on N.
|Figure 11: Mass density ( top) and potential ( middle) and error index ( bottom) for an inviscid, self-gravitating polytrope at equilibrium with index and axis ratio 0.5.|
|Open with DEXTER|
To illustrate the method in a concrete case, we consider the classical problem of uniformly rotating, inviscid polytropes. For such a system, the equilibrium structure is governed by the following two equations (e.g. Hachisu 1986)
On the z-axis,
are no more singular, for any a, with the consequence that the secondary kernels can easily be determined, in principle without the need for density splitting (see however note 10). From Eqs. (12) and (37), we have
|Figure 12: Accuracy of the vertical field gZ on the z-axis for a flat, homogneous disc with inner edge , outer edge , and three different semi-thickness H. The altitude of the field point is Z=H/2.|
|Open with DEXTER|
There is a relationship between the number of source points M on a given cylinder, the s-parameter and the physical resolution
in the vertical direction, namely
Table 3: Theoretical formulae for the physical resolution in the vertical direction as functions of the number M of source points,s-parameter and disc semi-thickness h. The last column corresponds to a 3-point, finite-difference method.
Table 4: Same legend as for Table 3 but for the physical resolution in the radial direction ( for all examples). , and denote the number of cylinders populating the left-side, long-range, the local and right-side, long-range contributions respectively (see the Appendix C).
A similar analysis can be led in the radial direction where the physical resolution is defined by
Intrinsically, finite difference methods (FDMs) are known to be extremely fast compared to integral methods (IMs). But comparing methods together has sense only if their accuracy as well as their ability for finding solutions ab nihilo for various kinds of systems are considered too
On this basis, best methods must i) be "self-starting''; ii) have low execution times ;
iii) be highly accurate (i.e. low
values); and iv) be able to treat any shape. For astrophysical discs, tori and rings, conditions i) and iv) are not satisfied by FDMs (because FDMs cannot start without boundary conditions, and because grid meshes generaly do not exactly fit the system surfaces). Regarding points ii) and iii), the efficiency
of a method can be defined as
|Figure 13: Experimental efficiencies to compute the potential in a square, computational grid with a 5-point finite difference method (FDM) and with the present (double splitting) method for system 1 and density/shape model 2 (see text).|
|Open with DEXTER|
The ratio gets larger i) as the system aspect ratio H/a decreases, and ii) as the system extension increases. In particular, Eq. (57) does not account for the dynamical mesh. For large and flat systems, FDMs need a drastic reduction of the mesh size in the radial direction in order to reach a given precision, since for the potential. This means a considerable increase of the computing time. For instance, for system 2 where and (see Table 2 and Fig. 10), is required to obtain everywhere a relative precision on of the order of only. As a result, is inevitably multiplied by a huge factor (1010 typically). Such a runaway is totally avoided in the present method, thanks to the dynamical mesh (note that is attained with less than a hundred points in the radial direction). As shown in previous sections, N never greatly exceeds M, so that is not very sensitive to the system size.
We have presented efficient numerical (splitting) methods to compute the field and the potential inside axially symmetric sources. The implementation of a dynamical mesh together with a high-order quadrature scheme enables a very accurate determination of Poisson integrals for a wide variety of system shapes and sizes. Performances and convincing tests have been reported. We have derived all necessary formulae to make the implementation of the method into various codes easier. Improvement, modifications and extensions of this work are obviously possible. For instance, better results can probably be obtained by considering different quadrature rules and/or other sampling functions. In particular, it would be interesting to enhance the concentration of source points not only around the field point, but in all regions of the system where the density field could exhibit important gradient variations. We have always assumed (for simplicity) that all cylinders have vertically the same number M of source points. Since the field and the potential decrease quite rapidly with the separation, it is conceivable to use a small number of source points for cylinders located (very) far away from the coincident cylinder (depending on the density distribution). This can be advantageous in terms of computing time. The challenge is then to find the most appropriate distribution M(a) without losing the precision.
The major assumption made here is axial symmetry which, for some problems, is too restrictive. The splitting method can easily be extended to non-axisymmetric configurations. Actually, in a fully tri-dimensional case, Eq. (21) can be generalized as follows
I am grateful to S. Bonazzola for enlightening discussions both on mathematical and physical aspects of the problem. I thank my colleagues R. Grappin, V. Karas, J. Léorat, M. Mouchet, D. Pelat, A. Pierens and L. Subr for many comments and suggestions. I thank the referee Howard Cohl for valuable comments and suggestions.
For a correct numerical treatment of Eq. (23), it is necessary to know exactly
at the coincident point. This can be done by considering Eq. (24) in the limit where
by construction at the coincident point, a Taylor expansion of
around z=Z at all orders gives
For density model 1 given in Table 1, the secondary kernel associated with the vertical component is
The total number of cylinder is defined as follows. First, we decide the number of cylinders standing in the domain
cylinders. The regular spacing of the v in this range then yields the location an of the cylinders that form the local contribution.
In the local contribution, Eq. (54) reads
Second, we determine the number of cylinders populating the long-range domains (
for the left-side and right-side respectively), such that the radial sampling
has continuous derivatives at