A&A 434, 1-15 (2005)
DOI: 10.1051/0004-6361:20034194

Solutions of the axi-symmetric Poisson equation from elliptic integrals

I. Numerical splitting methods

J.-M. Huré^1,2

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

Abstract
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 $\rho \propto z^n$. Improvements and extensions (such as the non-axi-symmetric case) are possible and discussed.

Key words: hydrodynamics - methods: numerical

   
1 Introduction

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

$$\Delta \Psi = 4 \pi G \rho,$$ (1)

where $\Psi$ is the self-induced potential and $\rho$ is the mass density. Cases where this equation has an analytical, tractable solution are rare or correspond to very simplified situations (e.g. Binney & Tremaine 1987), because astrophysical fluids have i) complex geometries/shapes; ii) very small/very large aspect ratios; iii) at least bi-dimensional density gradients; and iv) open boundaries. Numerical solutions, difficult to compute in general, are very powerful as they give a mean to investigate real systems. Methods of inversion of Eq. (1), also called "Poisson solvers'', are customarily classified into two categories (Müller & Steinmetz 1995): difference methods and integral methods. In difference methods, Eq. (1) is solved for $\Psi$ on a grid, requiring boundary conditions. Dirichlet/Neumann boundary conditions are exactly known at infinity (the potential and derivatives vanish), but dealing with an infinite domain is not convenient, unless a space mapping is considered (e.g. Clément 1974). Boundary conditions can easily be computed with great accuracy just outside the source. Internal solutions found with difference methods have generally low or limited precision (even with exact boundary conditions), but algorithms are extremely fast once initialized. One may then prefer the integral method where $\Psi$ is obtained from the Poisson integral

$$\Psi( \vec{r}) = - G \iiint_{\cal V}{\frac{\rho(\vec{r'}) {\rm d}^3 r' }{\vert \vec{r'} - \vec{r}\vert}},$$ (2)

with the advantage that the summation over the gas volume ${\cal V}$ naturally stops at the system boundaries. This approach not only produces missing Dirichlet boundary conditions (Bodo & Curir 1992; Störzer 1993; Cohl & Tohline 1999), but can be used in the entire domain. Integral methods are often rejected because of point mass singularities (i.e. when $\vec{r} \rightarrow \vec{r'}$ in Eq. (2)) which make the above integral improper everywhere inside the fluid. Often, people circumvent the difficulty by softened gravity (e.g. Adams et al. 1989). The second reason - probably the main one - is their possible prohibitive execution time. But when singularities are properly treated, potential values can be determined with high accuracy, orders of magnitudes larger than finite-difference based methods allow. With the era of fast computers and especially of parallel programming (e.g. Borges & Daripa 2001), integral methods should be more and more accessible, and it seems justified to develop efficient tools accordingly. Spectral methods - a third category - are probably the most efficient techniques (Greengard & Lee 1996; Christopher et al. 1997; Chen et al. 2000; Lai et al. 2002), at least for compact systems. This approach, possibly combined with mapping and space compactification, is commonly employed to model rotating stars (Clément 1974; Bonazzola et al. 1998), but it remains exceptional for toroidal configurations (Ansorg et al. 2003a,b). Expanded Green functions are not convenient in practice because these are alternate and infinite series which do not converge inside sources as often outlined (e.g. Clément 1974; Störzer 1993; Cohl & Tohline 1999; see also Cohl 1999). A systematic cutoff in the expansion (whatever the truncation order, generally a low order) always creates non-uniform and unverifiable errors.

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 $(R,\Phi)-$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

$$\vec{g}=-\nabla \Psi$$ (3)

is more useful than the potential itself. An exception is for instance the theory of Maclaurin spheroids and Dyson-Wong rings which requires only the potential (e.g. Hachisu 1986). Because numerical derivatives have poor precision, the field components obtained from $\Psi$ through Eq. (3) would be inaccurate in most cases, with roughly two or three correct digits at best. As we shall show here, the relative precision of the potential and of field components when considered separately can be very high (provided special care is taken in the treatment of point mass singularities), and can even attain computer precision. Indeed, the computer precision can generally not be reached with finite differentiating. Thus, depending on the actual problem, it can be interesting to have two different methods, one for $\Psi$ and one for $\vec{g}$; the latter is obtained from

$$\vec{g}( \vec{r}) = G \iiint_{\cal V}{\frac{\rho(\vec{r'}) (\... ... \vec{r}) {\rm d}^3r' }{\vert \vec{r'} - \vec{r}\vert^3}}\cdot$$ (4)

To our knowledge, the concepts and techniques presented here have never yet been published, at least in the astrophysical context, and bring new material in many respects: an integral approach, accuracy, range of application, simplicity, analytical treatment. We are not aware of any method aiming to estimate Eqs. (2) or (4) rigorously for toroidal configurations (see Bonazzola & Schneider 1974). Our goal is to provide the community with efficient tools to investigate the structure, stability and dynamical evolution of self-gravitating tori, discs and rings $^{\ref{note:tdr}}$, which are ubiquitous objects in the Universe. Obviously, this applies not only to Gravitation but also to Electrostatics.

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 $\vec{r} \rightarrow \vec{r'}$ 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

   
2 The field: Theoretical grounds

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 $\vec{r}(R,\Phi,Z)$ for field points, and $\vec{r'}(a,\phi,z)$ for source points, even though field points and source points can be confused.

2.1 Self-gravitating field components from elliptic integrals: Formal expressions

An infinitely thin (i.e. topologically one dimensional) homogeneous, circular loop with radius a, linear mass $\lambda$ and located at altitude z as pictured in Fig. 1, generates at any field point a gravity field

$$\vec{g}^{\rm loop}=\lambda \vec{\kappa} = \lambda \left(\begi... ...ay}{c} \kappa_R\\ \kappa_\Phi\\ \kappa_Z \end{array}\right),$$ (5)

where $\kappa_R$, $\kappa_\Phi$ and $\kappa_Z$ are components of the Poisson kernel $\vec{\kappa}$. It can be shown from Eq. (2) that these quantities are respectively given by (see also Durand 1964)

$$\kappa_R = \frac{G}{R} \sqrt{\frac{a}{R}} k \left[ E(k)- K(k) + \frac{(a-R)k^2 E(k)}{2 a{k'}^2} \right],$$ (6)

$$\kappa_\Phi=0$$ (7)

due to the axial symmetry, and

$$\kappa_Z = \frac{ G \left( z - Z\right)}{2 R \sqrt{aR}} \frac{k^3 E(k)}{{k'}^2},$$ (8)

where K and E are the complete elliptic integrals of the first and second kinds respectively (see a plot in Fig. 2), k is their modulus defined by

$$k^2 = \frac{4 aR}{(a+R)^2+(z-Z)^2}, \qquad 0 \le k \le 1,$$ (9)

and $k'= \sqrt{1 - k^2}$ is the complementary modulus. On the loop, k=1. This occurs when simultaneously R=a and Z=z, with the result that $\kappa_R \rightarrow \infty$ and $\kappa_Z \rightarrow \infty$. This is the "loop singularity'' (i.e. the axial symmetry analogue of the point mass singularity). From a mathematical point of view, the singularity is due to the complete elliptic integral of the first kind as $k \rightarrow 1$ (a logarithmic singularity), as well as to the hyperbolic form 1/k'. As soon as k < 1, Poisson kernels are regular. Finally, $k \rightarrow 0$ far from the loop, and the components of the kernel steadily vanish.

	Figure 2: Complete elliptic integrals K and E, and the "regular'' complete integral of the first kind $K^{\rm reg}$ 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 density $\rho(a,z)$ is given by

$$\vec{g}^{\rm }= \iint_{\cal S}{\rho(a,z)\vec{\kappa} {\rm d}a{\rm d}z}$$ (10)

where the double integral extends all over the cross-section ${\cal S}$ of the system (i.e. cross-section in the (a,z)-plane). The loop singularity is present everywhere inside the source, making the numerical determination of the field $\vec{g}$ from Eq. (10) tricky. Indeed, no reliable field values can be derived from integral methods without caution.

	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

   
2.2 Notations and assumptions

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: $a_{\rm in} > 0$ denotes the inner edge of the system, $a_{\rm out} > a_{\rm in}$ is the outer edge, and $\vert z\vert=H(a) \ge 0$ is the equation of the surface. Inside the fluid, the mass density can be given either in the form of an analytical expression $\rho(a,z)$, or as a grid by a two entry table $\rho_{n,m} = \rho(a_n,z_m)$. 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 (a_n,z_m) 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. $\rho(a,H)=0$ with $\partial_z \rho(a,H) \ne \infty$, 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 $^{\ref{note:tdr}}$.

   
3 Vertical integration of the kernel associated with the field: "Density splitting method''

3.1 Coincident cylinder and coincident point

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

$$\vec{g}=\int_{a_{\rm in}}^{a_{\rm out}}{\vec{\breve{\kappa}}{\rm d}a},$$ (11)

where

$$\vec{\breve{\kappa}}=\int_{-H}^{H}{\rho(a,z)\vec{\kappa} {\rm d}z}.$$ (12)

This comes down to replacing the continuous medium by a series of N cylinders (see Fig. 3), each with radius a_n where

$$a_1 \equiv a_{\rm in} \le a_n \le a_N \equiv a_{\rm out}, \qquad n \in [1,N]$$ (13)

and local semi-thickness H(a_n). In the following, $\vec{\breve{\kappa}}$ is called the "secondary kernel''. It corresponds to the field induced by an infinitely thin cylinder. Among these N cylinders, we call the "`coincident cylinder" the cylinder that contains the actual field point. Its radius is a=R. There, the semi-thickness is $H(R) \equiv h$. Further, we call the "coincident point'' the point belonging to the coincident cylinder and which coincides with the actual field point. It has coordinates a=R and z=Z. Thus, at the coincident point, k=1 by construction and the kernel is singular (see again Fig. 3).

   
3.2 Field due to an infinitely thin, homogeneous cylinder

The exact expression for Eq. (12) is known for a vertically homogeneous distribution, that is for $\rho(a,z)=\rho_0(a)$. Actually, we have for any field point (Durand 1964)

 

$$\rho_0 \int_{-H}^{H}{\vec{\kappa}{\rm d}z} -2 G \rho_0 \sqrt{\frac{a}{R}} \left(\begin{array}{c} \displaysty... ...(Z+H)}{2R} k^- \Lambda(m,k^-)\\ [4mm] k^+ K(k^+) - k^-K(k^-)\end{array}\right)$$ = (14)

$$\equiv \vec{\breve{\kappa}}^{(0)} ,$$

where $\Lambda(m,k)$ is defined by

$$\Lambda (m,k) = K(k)-\varpi \Pi(m,k),$$ (15)

$\Pi$ is the complete elliptic integral of the third kind $^{\ref{note:ekpi}}$, k^- and k⁺ correspond respectively to the bottom and top edges of the cylinder (see Fig. 3), namely

$${k^\pm} = \sqrt{\frac{4aR}{(a+R)^2+(\pm H-Z)^2}},$$ (16)

m is defined by

$$m=\frac{2\sqrt{a R}}{a+R},$$ (17)

and

$$\varpi=\frac{a-R}{a+R}\cdot$$ (18)

We have $0 \le {k^\pm} \le m \le 1$, $\varpi^2=1-m^2$ and $\vert\varpi\vert \le 1$. Note that the expression for $\vec{\breve{\kappa}}^{(0)}$ simplifies noticeably for the coincident cylinder since for m=1, we have

$$\lim_{m \rightarrow 1, \; a \ge R} \varpi \Pi(m,k) = + \frac{\pi}{2k'},$$ (19)

and so, Eq. (14) becomes

$$\vec{\breve{\kappa}}^{(0)} = - 2 G \rho_0 \left(\begin{array}... ...{k'}^-}\right]\\ \\ k^+ K(k^+) - k^-K(k^-)\end{array}\right).$$ (20)

   
3.3 Singularity removal: The "density splitting'' method

At the coincident point, the kernel $\vec{\kappa}$ diverges. So does the integrand $\rho\vec{\kappa}$ in Eq. (12). In this case, we re-write the density in the form

$$\rho(a,z) = \rho_0 + \delta \rho(a,z)$$ (21)

where $\rho_0 = \rho(a,Z)$ is the density at the field point, and $\delta \rho(a,z)$ is the "residual'' density profile along the coincident cylinder. By virtue of the superposition principle, we have from Eqs. (12) and (21)

$$\vec{\breve{\kappa}} = \vec{\breve{\kappa}}^{(0)} + \vec{\breve{\kappa}}^{\rm res.},$$ (22)

where $\vec{\breve{\kappa}}^{(0)}$ corresponds to a vertically homogeneous cylinder (it is given by Eq. (14)), and $\vec{\breve{\kappa}}^{\rm res.}$ is the "residual'' secondary kernel

$$\vec{\breve{\kappa}}^{\rm res.} = \int_{-H}^{H}{ (\delta \rho) \vec{\kappa} {\rm d}z},$$ (23)

with

 

$$(\delta \rho) \vec{\kappa} = G (\delta \rho) \left(\begin{array}{... ...style\frac{2}{(z-Z)\sqrt{1+ \left(\frac{z-Z}{2R}\right)^2}} \end{array}\right).$$ (24)

It can be shown (see the proof in the Appendix A) that $(\delta \rho) \vec{\kappa}$ is finite at the coincident point, namely

$$\left\{\begin{array}{l} (\delta \rho) \kappa_R = 0 \\ \\ \d... ...artial \rho }{\partial z}\right\vert _{z=Z} \end{array}\right.$$ (25)

for k=1. This means that the integrands in Eq. (23) are regular over the whole range [-h,h]. It follows that $\vec{\breve{\kappa}}^{\rm res.}$ can be estimated with standard integration schemes, for any a, any R, any Z and for any mass density profile. Since $\vec{\breve{\kappa}}^{(0)}$ is known, the secondary kernel $\vec{\breve{\kappa}}$ can be easily determined. Note that vertical partial derivatives $\partial_z \rho$ are in principle needed in this procedure (unless a dynamial mesh is implemented; see Sect. 3.5).

	Figure 4: Illustration of the density splitting method. The cylinder a) has radius a=1, semi-thickness H=1, and the mass density $\rho$ follows a parabolic law b). The field point has arbitrarily coordinates R=a and $Z=\frac{H}{2}=\frac{h}{2}$. There, both integrands $\rho \kappa _R$ and $\rho \kappa _Z$ diverge c).
Open with DEXTER

3.4 Simple illustration of the method

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 $Z=\frac{H}{2}=\frac{h}{2}$. There, both integrands $\rho \kappa _R$and  $\rho \kappa _Z$ diverge because of the loop singularity. These two quantities are displayed in Fig. 4c. At the field point, the density takes the value $\rho(R,Z)=\rho_0$. The residual density $\delta \rho = \rho -\rho _0$ is plotted in Fig. 5a. The "field'' due to a cylinder with constant density $\rho _0$ (that is $\vec{\breve{\kappa}}^{(0)}$) is known in a closed form. So, the contribution due to $\delta \rho$ remains to be computed numerically. Components of the residual integrand $(\delta \rho) \vec{\kappa}$ 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 $\vec{\breve{\kappa}}^{\rm res.}$, 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 $\delta \rho = \rho -\rho _0$ vanishes at the field point, where the mass density at the field point $\rho _0$. 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

   
3.5 Vertical sampling: A dynamical mesh

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

$$w=\left\{\begin{array}{ll} \displaystyle+ \left(\frac{z-Z}{2h... ... \\ & \mbox{(below the coincident point)} \end{array}\right.$$ (26)

where s is a parameter. We shall discuss later what this choice implies in terms of the physical resolution (i.e. the smallest size that can be resolved vertically). Provided 0 < s < 1, a regular spacing in the w-variable produces the required enhancement of the number of source points around z=Z. In other words, for M sources points located in the domain [-h,h], we have for the coincident cylinder

$$\left\{\begin{array}{l} z_m = Z \pm 2h \left( w_m \right)^{1/... ... \;\; w_M=\left(\frac{Z-h}{2h}\right)^s,\\ \end{array}\right.$$ (27)

for $m=2,\dots,M-1$, where one of the z_m equals Z. For most situations, we have noticed that the choice $s=\frac{1}{2}$ appears the best choice for both $\breve{\kappa}_R$ and $\breve{\kappa}_Z$, and we work with this default value in the following. Note that the integral $\int \dots {\rm d}z$ in Eq. (23) can be converted into the form $\int \dots {\rm d}w$. Since $\frac{{\rm d}z}{{\rm d}w}=0$ at w=0, we have $(\delta \rho) \vec{\breve{\kappa}}= \vec{0}$ at the coincident point. As a consequence, the vertical derivatives $\partial_z \rho$ are no longuer needed (see Eq. (25), and see Appendix A, Eq. (A.4)).

Table 1: Two mass density/shape models considered in this study. For model 9 1, an exact expression for both $\breve{\kappa}_Z$ and g_Z can be derived analytically from Eqs. (8) and (9) provided $H=\rm const.$ (see the Appendix B, Eq. (B.2); the kernel and field are finite for $k^\pm \ne 1$). 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 $\rho _1 \propto 1/a^2$ for $R \ll a$ whereas $\rho _1 \sim \rm const.$ for $R \gg a$.

Figure 6: Top panel: error index $\epsilon(\breve{\kappa}_Z)$ versus the number of vertical source points M, obtained for a cylinder with radius a=1 and density profile $\rho _1$ (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 $Z=\frac{H}{2}=\frac{h}{2}$. 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

   
3.6 Method checking

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. $\vec{\breve{\kappa}}^{(0)} =\vec{0}$). We shall restrict ourselves to the vertical component because of the simple form of the kernel $\kappa_Z$ 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 $\rho _1$, denoted $\breve{\kappa}^{\rm ref.}_Z$, to values obtained from the density splitting method. In the following, comparisons are made with the "error index'' $\epsilon(x)$ defined as

$$\epsilon \equiv \left\{\begin{array}{ll} \epsilon_{\rm CP}, ... ...{\rm ref.}} \right\vert, & \mbox{otherwise} \end{array}\right.$$ (28)

where $\epsilon_{\rm CP}$ corresponds to the typical computer precision (generally $\la$-14 for double precision), x is some quantity and $x^{\rm ref.}$ denotes the reference value (supposed to be better than x). Figure 6a displays the error index $\epsilon(\breve{\kappa}_Z)$ versus $\log_{10} M$. As above, we have considered a cylinder with radius a=1, but four different aspect ratios $\frac{H}{a}=10^{-4}$, 10^-2, 0.1 and 1. We see that a relative precision ranging from $\sim$10^-3 to 10^-12 (depending on the aspect ratio $\frac{H}{a}$) is reached for only 10 source points, which is remarkable. This is essential regarding the computing time. Unsurprisingly, the smaller the aspect ratio, the larger the precision for a given M. We also see that the precision is limited by effects due to the computer maximal precision for about a hundred points in the case with $\frac{H}{a} = 0.01$. Due to the quadrature scheme, the method response is roughly given by the formula

$$\epsilon(\breve{\kappa}_Z) \approx 2 -6 \log_{10} M + 2 \log_{10} \frac{h}{a},$$ (29)

which means almost two more digits fixed for two times more points. Note that 2nd order centered schemes would give $\frac{{\rm d}\epsilon}{{\rm d} \ln M} \sim - 1$ for the field components from Eqs. (1) and (3), and whatever the vertical extent.

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 $\breve{\kappa}^{\rm ref.}_R$ and $\breve{\kappa}^{\rm ref.}_Z$, 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 $R \ne a$ (and $H \ne h$), gives much better results (i.e. $\epsilon_{\rm CP}$ is reached for fewer source points), unsurprisingly.

Figure 7: Typical shape of secondary kernels $\breve{\kappa}_R$ ( left panels) and $\breve{\kappa}_Z$ ( right panels) versus a at two different scales for a flat disc, with inner edge $a_{\rm in}=1$, outer edge $a_{\rm out}=10^5$ 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

   
4 Field values

   
4.1 Radial sampling

The radial integration of the secondary kernel $\vec{\breve{\kappa}}$ 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 $[a_{\rm in}, a_{\rm out}]$ in an adequate manner (see Fig. 3), and using the same quadrature rule. Let us call v(a) this function. For systems with $H \sim a_{\rm out}$ (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:

• a "local'' contribution, due to the matter located in the neighborhood of the coincident cylinder. The corresponding secondary kernel is very narrow but has large amplitude;
• a "long-range'' contribution, due to the matter located far away. The secondary kernel is very wide, but has a tiny amplitude.

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 $m \approx 1$) but also far from it. This is achieved with the following piece-wise sampling-function

$$v=\left\{\begin{array}{ll} \displaystyle -\ln \left( \frac{R}... ...m out} \\ & \mbox{(right-side, long-range)} \end{array}\right.$$ (31)

where $\ell$ and q are free parameters. In particular, $\ell$ must typically represent a few vertical local thicknesses or a fraction of it (depending on the system aspect ratio and radial, density gradients), that is $\ell \sim h$. With $0 < q < \frac{1}{2}$, the number of cylinders is increased around the coincident cylinder with a regular v-mesh, namely

$$\left\{\begin{array}{ll} a_n = R \left(1+e^{\pm v_n}\right)^{... ... \right]^q&\mbox{(local domain)}, \\ [4mm] \end{array}\right.$$ (32)

for $n=2,\dots,N-1$, where one of the a_n equals R.

In practice, the value $q=\frac{1}{2}$ is well suited for general use. We impose only the number $N_{\rm local}$ 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 $N_{\rm left}$ and $N_{\rm right}$ 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.

4.2 Method checking

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 $\epsilon (g_Z)$ (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 $g^{\rm ref.}_Z$ 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 $\sim$10^-11 on average (this would be $\sim$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 $N \times M$. Again, the remarkable issue is that, with a few points only, the relative precision is already good. For instance for system 1, a $16 \times 16$ grid yields $\epsilon(g_Z) \la -5$ (this is achieved with a $\sim$ $200 \times 16$ 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 $\epsilon (g_Z)$ 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 $Z=\frac{H}{2}$. 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 g_R and g_Z 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

$$\epsilon \approx 2 - 6 \log_{10} M \pm1,$$ (33)

for each field component. Obviously, density/shape models considered here as tests are smooth functions, and so these performances have some limits. In other hand, any method that would not behave well in the presence of smooth distributions would apply in special conditions.

	Figure 9: Same legend and same conditions as for Fig. 8 (except $q=\frac{1}{4}$) but for a geometrically thin and flared disc (system 3) with density/shape model 2 (see Table 1).
Open with DEXTER

   
5 Numerical treatment of the self-gravitating potential: The "double splitting'' method

   
5.1 Kernels and potentials

Let us now consider the potential. Following the formalism employed in Sect. 2, the potential $\Psi$ due to an infinitely thin, homogeneous circular loop (see Fig. 1) with linear mass $\lambda$ and radius a is (e.g. Durand 1964)

$$\Psi^{\rm loop} = \lambda(a,z) \kappa_\Psi$$ (34)

where the kernel $\kappa_\Psi$ is given by

$$\kappa_\Psi = - 2 G \sqrt{\frac{a}{R}}k K(k).$$ (35)

Then, for any tri-dimensional, axially symmetric system with mass density $\rho(a,z)$, the potential can be written in the form

$$\Psi=\int_{a_{\rm in}}^{a_{\rm out}}{\breve{\kappa}_\Psi {\rm d}a}$$ (36)

where

$$\breve{\kappa}_\Psi=\int_{-H(a)}^{H(a)}{\rho\kappa_\Psi {\rm d}z}$$ (37)

is the secondary kernel.

	Figure 10: Same legend and same conditions as for Figs. 8 and 9 but for the potential $\Psi$ (see text for more details).
Open with DEXTER

5.2 Kernel splitting

We see from Eq. (35) that the straightforward determination of the secondary kernel $\breve{\kappa}_\Psi$ 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 $\rho_0(a)$ as considered in Sect. 3.2 does not help here. Since the asymptotic form of the K-function as $k \rightarrow 1$ is known, we first set

$$K^{\rm reg}(k) = K(k) + \ln k'$$ (38)

where $K^{\rm reg}(k)$ denotes the "regular'' complete elliptic integral of the first kind. This function is plotted versus k in Fig. 2. It is a smooth, regular function in the range [0,1]. Given this "kernel splitting'', Eq. (37) is

$$\breve{\kappa}_\Psi = \breve{\kappa}_\Psi^{\rm reg.} + \breve{\kappa}_\Psi^{\rm sin.}$$ (39)

where

$$\breve{\kappa}_\Psi^{\rm reg.} = - 2 G \sqrt{\frac{a}{R}} \int_{-H}^{H}{\rho k K^{\rm reg}(k) {\rm d}z}$$ (40)

is the regular part, and

$$\breve{\kappa}_\Psi^{\rm sin.} = 2 G \sqrt{\frac{a}{R}} \int_{-H}^{H}{\rho k \ln k' {\rm d}z}$$ (41)

is the remaining, singular part. We see that $\breve{\kappa}_\Psi^{\rm reg.}$ can be determined numerically without problems since the integrand in Eq. (40) is regular in any case. The singularity remains but it is now in Eq. (41), in an explicit form.

5.3 Density splitting

We can now use the analogue of the "density splitting'' method to compute $\breve{\kappa}_\Psi^{\rm sin.}$. We set $\tilde{\rho}=\rho k$ and split this quantity into two terms according to

$$\tilde{\rho}(a,z) = \left. \tilde{\rho} \right\vert _{k=m} + \delta \tilde{\rho}(a,z).$$ (42)

For the coincident cylinder (i.e. m=1), this relation becomes

$$\tilde{\rho}(a,z) = \rho_0 + \delta \tilde{\rho}(R,z)$$ (43)

where $\rho_0 = \rho(a,Z)$ is the mass density just at the field/coincident point, and $\delta \tilde{\rho}$ is the residual profile, i.e. deviation to the local value. Under these conditions, Eq. (41) is

 

$$\breve{\kappa}_\Psi^{\rm sin.} = 2 G \sqrt{\frac{a}{R}} \left. \t... ... + 2 G \sqrt{\frac{a}{R}} \int_{-H}^{H}{(\delta \tilde{\rho}) \ln k' {\rm d}z}.$$ (44)

The first term in the right-hand-side of this equation is known (see Appendix D). The second term (the "residual kernel''  $\breve{\kappa_\Psi}^{\rm res.}$) can easily be determined numerically using standard quadrature rules since $\delta ( \rho k) \ln k' =0$ at the field point despite $\ln k' \rightarrow \infty$ (the proof that can be made from a Taylor expansion about z=Z; see for instance Appendix A).

   
5.4 Method checking

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 $\breve{\kappa}_\Psi$ (relative precisions are similar to Sect. 3.6). Figure 10 is the analogue of Figs. 8 and 9: it displays the error index $\epsilon (\Psi )$ 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 $\breve{\kappa}^{\rm ref.}_\Psi$, 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 $\rho$ ( top) and potential $\Psi$ ( middle) and error index  $\epsilon (\Psi )$ ( bottom) for an inviscid, self-gravitating polytrope at equilibrium with index $\gamma =2$ and axis ratio 0.5.
Open with DEXTER

   
5.5 An example of astrophysical application: Equilibrium of a compressible Dyson-ring

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)

$$\left\{\begin{array}{ll} \displaystyle \frac{\gamma}{\gamma-1... ...cal E}\\ \\ \Delta \Psi = 4 \pi G \rho\\ \end{array}\right.$$ (45)

where $0 \le \gamma \ne \infty$ is the polytropic index, $\Omega$ is the rotation frequency at any point of the fluid, and ${\cal E}$ and $c_{\rm s}$ are constants. Once $\Omega$, $c_{\rm s}$, $\gamma$ and ${\cal E}$ are prescribed, the solution of the above system can be determined using the Self-Consistent Field (SCF) method (Ostriker & Mark 1968). Figure 11 shows the equilibrium mass density, potential and associated error index $\epsilon (\Psi )$ for a polytropic index $\gamma =2$ and axis ratio $\frac{a_{\rm in}}{a_{\rm out}}=0.5$, a case considered in Hachisu (1986). In this example, $\Psi$ is deliberately computed with a small size, dynamical grid containing about $N \times M \sim 32 \times 32$ density points. The error index is compatible with results shown in Fig. 10: the relative precision on the potential inside the source cannot be better than $\sim$10^-6 for the actual grid size. Our results are however in good agreement with Hachisu's work within a few percent typically, despite a different technique and in particular much fewer source points. Higher precisions are obtained with finer grids.

   
5.6 Potential and field on the z-axis

On the z-axis, $\kappa_\Psi$ and $\kappa_Z$ 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

$$\breve{\kappa}_Z = 2 \pi G \int_{-H}^H{\frac{\rho a(z-Z){\rm d}z}{\left[a^2+ (z-Z)^2\right]^{3/2}}}$$ (46)

and

$$\breve{\kappa}_\Psi=-2 \pi G \int_{-H}^{H}{\frac{\rho a {\rm d}z}{\sqrt{a^2+ (z-Z)^2}}}$$ (47)

respectively. It follows that the homogeneous contributions are simply given by

$$\left\{\begin{array}{ll} \displaystyle \breve{\kappa}_Z = - 2... ...rm arcsinh} \left( \frac{H-Z}{H+Z} \right)& \end{array}\right.$$ (48)

and the corresponding residual secondary kernels are

$$\left\{\begin{array}{ll} \displaystyle \breve{\kappa}_Z = 2 \... ...delta \rho) a dz}{\sqrt{a^2+ (z-Z)^2}}}&\\ \end{array}\right.$$ (49)

where $\delta \rho = \rho(a,z) - \rho_0(a,Z)$. Figure 12 gives a concrete illustration of this case.

	Figure 12: Accuracy of the vertical field gZ on the z-axis for a flat, homogneous disc with inner edge $a_{\rm in}=0$, outer edge $a_{\rm out}=1$, and three different semi-thickness H. The altitude of the field point is Z=H/2.
Open with DEXTER

   
6 Resolutions and efficiency

6.1 Vertical resolution: Thick discs vs. thin discs

There is a relationship between the number of source points M on a given cylinder, the s-parameter and the physical resolution  ${\cal R}_Z$ in the vertical direction, namely

$${\cal R}_Z = \frac{\Delta z}{2h}$$ (50)

for the coincident cylinder, where $\Delta z=z_{m+1}-z_m$ is the distance between two neighboring source points ( $m \in [1,M-1]$). From Eq. (26), we have

$$s\Delta z = \pm 2h (\pm w)^{\frac{1-s}{s}} \Delta w$$ (51)

where

$$\Delta w = \frac{1}{(2h)^s} \left[ (h-Z)^s+(h+Z)^s \right].$$ (52)

Two extreme cases can be considered: the resolution at the field point (where $w \sim \Delta w$), and the resolution far from it (i.e. $w \rightarrow w^\pm \ne \Delta w$). Table 3 gives the expected expression for ${\cal R}_Z$ depending on the local disc aspect ratio $\frac{h}{R}$, and gives a few typical examples. We see that there is always a large resolution at the field point, much larger than allowed in finite-difference methods for the same M. This always ensures that the secondary kernel is determined with a certain accuracy.

Table 3: Theoretical formulae for the physical resolution ${\cal R}_Z$ 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 ${\cal R}_R$ in the radial direction ($\ell =h$ for all examples). $N_{\rm left}$, $N_{\rm local}$ and $N_{\rm right}$ denote the number of cylinders populating the left-side, long-range, the local and right-side, long-range contributions respectively (see the Appendix C).

6.2 Radial resolution: Extended discs/rings vs. tori

A similar analysis can be led in the radial direction where the physical resolution is defined by

$${\cal R}_R=\frac{\Delta a}{L},$$ (53)

where $L=a_{\rm out}-a_{\rm in}$ is the radial extent of the system, and $\Delta a=a_{n+1}-a_n$ is the distance between two adjacent cylinders ( $n \in [1,N-1]$). From Eq. (31), we find

$$\displaystyle\Delta a=\Delta v \times \left\{\begin{array}{ll... ...\mbox{for} \; R+\ell \le a \le a_{\rm out} \end{array}\right.$$ (54)

where $\Delta v$ may depend on the domain. Table 4 gives the theoretical formula for ${\cal R}_R$ depending on the location inside the system. Five key-positions are considered: near the inner edge, at the field point, near the outer edge, and at the connection between the long-range and local contributions. A few typical values of given. Again, due to the dynamical mesh, the resolution is very large at the field point. This is specially pronounced for very extended and flat systems since ${\cal R}_R \ll h/L$.

   
6.3 Efficiency: Computing time/precision ratio

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 $^{\ref{note:arg}}$. On this basis, best methods must i) be "self-starting''; ii) have low execution times $\tau$; iii) be highly accurate (i.e. low $\epsilon$ 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 $\eta$ of a method can be defined as

$$\frac{1}{\eta} = \tau \times 10^\epsilon$$ (55)

and the larger $\eta$, the more efficient the method. Let us estimate $\eta$ in both cases, assuming for simplicity a compact system (i.e. $H \sim L$) and a square grid (i.e. M=N). A 5-point finite difference method is characterized by a computing time $\tau^{\rm FDM} \propto M^2 \ln M$ (e.g. Swarztrauber & Sweet 1996), once initialized. The preliminary determination of boundary conditions (BCs) from an integral method implies an additional time $\tau^{\rm BC} \propto M^3$ generally much larger than for the inversion of the Laplacian. Their relative error corresponds to $\epsilon^{\rm FDM} \sim M^{-2}$ for the potential (and only M^-1 for the field components). For integral methods, $\tau^{\rm IM} \propto M^4$ in general. With the actual quadrature scheme, we have $\epsilon^{\rm IM, ~ this ~ work} \propto M^{-6}$ (see previous sections). So, regarding the potential, efficiencies are

$$\eta^{\rm FDM} \propto \left\{\begin{array}{ll} \displaystyle... ...c{1}{M} & \mbox{inversion including BCs}\\ \end{array}\right.$$ (56)

and

$$\eta^{\rm IM, ~ this ~ work} \propto M^2$$ (57)

respectively. As a result,

$$\eta^{\rm IM, ~ this ~ work} \gg \eta^{\rm FDM}$$ (58)

making the present method more efficient in the limit of large M (and this is true for any integral method with $\frac{{\rm d} \epsilon}{{\rm d}\log_{10} M} < -4$). We have verified this issue experimentally for different configurations, different shapes and density profiles. Figure 13 gives the efficiency for both methods $^{\ref{note:cpu}}$ observed for system 1 and density/shape model 2 with $\rho_0(a)=\rm const.$ as already considered (see Tables 1 and 2, and Figs. 8-10), forcing M=N. These results are typical for compact systems, like tori and thick discs. In this example, we see that the present method is more efficient as soon as $M=N \gtrsim 16$ if we account for the fact that FDMs are not self-starting methods which is always the case for systems of interest here. Regarding field values, a similar analysis shows that $\eta^{\rm IM, ~ this ~ work}$ is unchanged whereas $\epsilon^{\rm FDM} \sim M^{-2}$, with the consequence that Ineq. (58) still holds and even is better satisfied.

	Figure 13: Experimental efficiencies $\eta$ to compute the potential in a square, $M \times M$ 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 $\eta^{\rm IM, ~ this work}/\eta^{\rm FDM}$ gets larger i) as the system aspect ratio H/a decreases, and ii) as the system extension $L=a_{\rm out}-a_{\rm in}$ 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 $\epsilon^{\rm FDM} \approx 2 \log_{10} \left[ {\rm Sup} \{ \frac{H}{M}, \frac{L}{N} \} \right]$ for the potential. This means a considerable increase of the computing time. For instance, for system 2 where $L/a_{\rm in}=10^5$ and $H/a_{\rm in}=1$ (see Table 2 and Fig. 10), $N \simeq 10^6$ is required to obtain everywhere a relative precision on $\Psi$ of the order of $\sim$$1 \%$ only. As a result, $\tau^{\rm FDM}$ is inevitably multiplied by a huge factor ($\sim$10¹⁰ typically). Such a runaway is totally avoided in the present method, thanks to the dynamical mesh (note that $1 \%$ is attained with less than a hundred points in the radial direction). As shown in previous sections, N never greatly exceeds M, so that $\tau^{\rm IM, ~ this ~ work}$ is not very sensitive to the system size.

7 Concluding remarks

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

$$\rho(a,\phi,z) = \rho_0 + \delta \rho(a,\phi,z)$$ (59)

where $\rho_0 \equiv \rho (R,\Phi,Z)$ is the mass density at the field point $(R,\Phi,Z)$, and $\delta \rho$ is the non-axially-symmetric, residual distribution which satisfies $\delta \rho(R,\Phi,Z)=0$ by construction. It then "suffices'' to integrate this residual distribution according to Eq. (4) (which has a regular integrand when $\vec{r'}=\vec{r}$), augmented with the contribution due to the axially symmetric cylinder with density $\rho _0$, from Eq. (14).

Acknowledgements

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.

   
Appendix A: Values of $\vec{(\delta \rho)} \vec{\kappa}$ at the coincident point

For a correct numerical treatment of Eq. (23), it is necessary to know exactly $(\delta \rho) \vec{\kappa}$ at the coincident point. This can be done by considering Eq. (24) in the limit where $Z \rightarrow z$. Since $\delta \rho=0$ by construction at the coincident point, a Taylor expansion of $\delta \rho$ around z=Z at all orders gives

$$\delta \rho (z-Z) \left. \frac{\partial \delta \rho}{\partial z} \right\vert ... ...Z)^2\left.\frac{\partial^2 \delta \rho}{\partial z^2}\right\vert _{z=Z} + \dots$$ =

$$(z-Z) \left[ \left.\frac{\partial \delta \rho}{\partial z} \right... ...z-Z)\left.\frac{\partial^2 \delta \rho}{\partial z^2}\right\vert _{z=Z} \right.$$ =

$$\left. + (z-Z)\sum_3^{\infty}{\frac{(z-Z)^{i-2}}{i!}\left.\frac{\partial^i \delta \rho}{\partial z^i}\right\vert _{z=Z}} \right]$$ (A.1)

where $\left.\frac{\partial^i \delta \rho}{\partial z^i}\right\vert _{z=Z} = \left.\frac{\partial^i \rho}{\partial z^i}\right\vert _{z=Z}$ due to Eq. (21). For m=1, we have

$$\lim_{z \rightarrow Z} K(k) = \lim_{k \rightarrow 1} K(k) = \ln \frac{4}{k'}$$ (A.2)

and $\vert z-Z\vert = k' \sqrt{4R^2+(z-Z)^2}$. Since E is bounded, we have

$$\left\{\begin{array}{l} \lim_{z\rightarrow Z} \; (\delta \rho... ...artial \rho }{\partial z}\right\vert _{z=Z} \end{array}\right.$$ (A.3)

with the consequence that $(\delta \rho) \vec{\kappa}$ is finite at the coincident point, with

 

$$\delta \rho \left(\begin{array}{c} \kappa_R\\ \kappa_Z \end{arra... ...G\left. \frac{\partial \rho }{\partial z}\right\vert _{z=Z} \end{array}\right).$$ (A.4)

For $z \ne Z$ along the coincident cylinder, $(\delta \rho) \vec{\kappa}$ is given by Eq. (24).

   
Appendix B: Density/secondary kernels test-pairs

For density model 1 given in Table 1, the secondary kernel associated with the vertical component is

$$\breve{\kappa}_Z= -2G\rho_0\int_{-H}^{H}{\frac{k{\rm d}k}{1-k^2}} = G\rho_0\ln \left( \frac{1-{k^+}^2}{1-{k^-}^2} \right)\cdot$$ (B.1)

This secondary kernel can then be integrated analytically in the radial direction, depending on the shape H(a). For $H=\rm const.$ (i.e. a flat system), we have

 

$$\rho_0\int_{a_{\rm in}}^{a_{\rm out}}{\breve{\kappa}_Z {\rm d}a}$$                                     g_Z =

$$G \rho_0\left[(a-R) \ln \frac{(a-R)^2+(H-Z)^2}{(a-R)^2+(H+Z)^2} \right.$$ = (B.2)

$$- (a+R) \ln \frac{(a+R)^2+(H-Z)^2}{(a+R)^2+(H+Z)^2}$$

$$+ 2(H-Z) \left( {\rm atan ~}\frac{a-R}{H-Z} - {\rm atan ~}\frac{a+R}{H-Z}\right)$$

$$\left. -2(H+Z) \left( {\rm atan ~}\frac{a-R}{H+Z} - {\rm atan ~}\frac{a+R}{H+Z}\right) \right]_{a_{\rm in}}^{a_{\rm out}}\cdot$$

   
Appendix C: Input parameters for the radial sampling

The total number of cylinder is defined as follows. First, we decide the number of cylinders standing in the domain $[R-\ell, R+\ell]$, say $N_{\rm local}$ cylinders. The regular spacing of the v in this range then yields the location a_n of the cylinders that form the local contribution. In the local contribution, Eq. (54) reads

$$\Delta a = \frac{(a-R)}{2q}\frac{\Delta v}{v}\frac{1}{1-(\pm v)^{1/q}}$$ (C.1)

depending on the sign of v. For a regular sampling in v, $\Delta v$ is given by

$$\Delta v = \frac{2}{N_{\rm local}-1} \left(\frac{n^2}{1+n^2} \right)^q$$ (C.2)

where $N_{\rm local}$ is the total number of cylinders in the local contribution. In particular, $\Delta a$ can be computed at $a = R \pm \ell$.

Second, we determine the number of cylinders populating the long-range domains ( $N_{\rm left}$ and $N_{\rm right}$ for the left-side and right-side respectively), such that the radial sampling $\{a_n\}$ has continuous derivatives at $a = R \pm \ell$. We find

$$\left\{\begin{array}{ll} \displaystyle N_{\rm left}-1 = \frac... ...ta a} \ln \frac{a_{\rm out}-R}{\ell}\cdot & \end{array}\right.$$ (C.3)

The radius a_n of these remaining cylinders is then obtained from regular spacing of the v-variable (see Eq. (31)). The total number of cylinders is then

$$N=N_{\rm left}+N_{\rm right}+N_{\rm local}-2.$$ (C.4)

   
Appendix D: Analytical formula

We have

$$\int{\ln k' {\rm d}z} = \frac{1}{2} \int{\ln \frac{(a-R)^2+(z-Z)^2}{(a+R)^2+(z-Z)^2} \; {\rm d}z}.$$ (D.1)

Since

$$\int{\ln \left( x^2+b^2 \right) {\rm d}x} = x \ln \left( x^2+b^2\right) - 2 \left( x - b ~ {\rm atan} \frac{x}{b} \right),$$ (D.2)

we have

$$2\int_{-H}^H{\ln k' {\rm d}z} (H-Z) \ln {k'}^+ - (H+Z) \ln {k'}^-$$ =

$$+ 2 \left[ (a-R) \left( {\rm atan ~}\frac{H-Z}{a-R} - {\rm atan ~}\frac{H+Z}{a-R} \right)\right.$$

$$- \left. (a+R) \left( {\rm atan ~}\frac{H-Z}{a+R} - {\rm atan ~}\frac{H+Z}{a+R} \right) \right]$$ (D.3)

where ${k'}^\pm = \sqrt{1-{k^\pm}^2} = \sqrt{\frac{(a-R)^2+(H \mp Z)^2}{(a+R)^2+(H \mp Z)^2}}\cdot$

References

• Abramowitz, M., & Stegun, I. A. 1964, in Handbook of mathematical functions with examples, graphs and mathematical tables, NBS, Sect. 17
• Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959 [NASA ADS] [CrossRef] (In the text)
• Adams, J., Swarztrauber, P., & Sweet, R. 1980, fishpack/NETLIB/NCAR
• Ansorg, M., Kleinwächter, A., & Meinel, R. 2003a, MNRAS, 339, 515 [NASA ADS] [CrossRef] (In the text)
• Ansorg, M., Kleinwächter, A., & Meinel, R. 2003b, ApJ, 582, L87 [NASA ADS] [CrossRef] (In the text)
• Bodo, G., & Curir, A. 1992, A&A, 253, 318 [NASA ADS] (In the text)
• Bonazzola, S., & Schneider, J. 1974, ApJ, 191, 273 [NASA ADS] [CrossRef] (In the text)
• Bonazzola, S., Gourgoulhon, E, & Marck, J.-A. 1998, Phys. Rev. D, 58, 104020 [NASA ADS] [CrossRef] (In the text)
• Borges, L., & Daripa, P. 2001, JCP, 169, 151 [NASA ADS] (In the text)
• Binney, J., & Tremaine, S. 1987, in Galactic dynamics (Princeton, NJ: Princeton University Press) (In the text)
• Chen, H., Su, Y., & Shizgal, B. D. 2000, JCP, 160, 453 [NASA ADS] (In the text)
• Christopher, I., Knorr, G., Shoucri, M., & Bertrand, P. 1997, JCP, 131, 323 (In the text)
• Clément, M. J. 1974, ApJ, 194, 709 [NASA ADS] [CrossRef] (In the text)
• Cohl, H. S., & Tohline, J. E. 1999, ApJ, 527, 86 [NASA ADS] [CrossRef] (In the text)
• Cohl, H. S. 1999, Ph.D. Thesis
• Durand, E. 1964, in Electrostatique. I. Les distributions (Masson) (In the text)
• Gill, P. E., & Miller, G. F. 1972, Comput. J., 80 (In the text)
• Greengard, L., & Lee, J. Y. 1996, JCP, 125, 415 (In the text)
• Hachisu, I. 1986, ApJS, 62, 46 [NASA ADS] (In the text)
• Lai, M., Lin W., & Wang W. 2002, JNA, 22, 537 (In the text)
• Müller, E., & Steinmetz, M. 1995, Comp. Phys. Commun., 89, 45 [NASA ADS] (In the text)
• Pierens, A., & Huré, J.-M. 2005, A&A, 434, 17 [EDP Sciences] [CrossRef] (Paper II) (In the text)
• Störzer, H. 1993, ApJ, 271, 25 [NASA ADS] (In the text)
• Swarztrauber, P., & Sweet, R. 1996, in Handbook of fluid dynamics and Fluid Machinery, ed. J. A. Schetz, & A. E. Fuhs (New-York: J. Wiley & Sons) (In the text)