A&A 475, 401-407 (2007)
DOI: 10.1051/0004-6361:20066808

Generation of potential/surface density pairs in flat disks

Power law distributions

J.-M. Huré^1,2 - D. Pelat² - A. Pierens³

1 - Université Bordeaux 1/CNRS/OASU/UMR 5804/LAB, 2 rue de l'Observatoire, BP 89, 33270 Floirac, France
2 - LUTh/Observatoire de Paris-Meudon-Nancay, Place Jules Janssen, 92195 Meudon Cedex, France
3 - Astronomy Unit, Queen Mary, University of London, Mile end Road, London E1 4NS, UK

Received 24 November 2006 / Accepted 5 June 2007

Abstract
Aims. We report a simple method to generate potential/surface density pairs in flat axially symmetric finite size disks.
Methods. Potential/surface density pairs consist of a ``homogeneous'' pair (a closed form expression) corresponding to a uniform disk, and a "residual'' pair. This residual component is converted into an infinite series of integrals over the radial extent of the disk. For a certain class of surface density distributions (like power laws of the radius), this series is fully analytical.
Results. The extraction of the homogeneous pair is equivalent to a convergence acceleration technique, in a matematical sense. In the case of power law distributions (i.e. surface densities of the form $\Sigma(R) \propto R^s$), the convergence rate of the residual series is shown to be cubic inside the source. As a consequence, very accurate potential values are obtained by low order truncation of the series. At zero order, relative errors on potential values do not exceed a few percent typically, and scale with the order N of truncation as 1/N³. This method is superior to the classical multipole expansion whose very slow convergence is often critical for most practical applications.

Key words: gravitation - methods: analytical - accretion, accretion disks - methods: numerical

1 Introduction

The construction of potential/density pairs is a common concern in the context of galactic dynamics (e.g. Binney & Tremaine 1987). More generally, the accurate determination of the gravitational potential corresponding to a given mass distribution is a critical step in astrophysical disk models and numerical simulations. Unfortunately, potential/density pairs are not closed-form expressions, except in a few cases only (de Zeeuw & Pfenniger 1988; Binney & Tremaine 1987; Mestel 1963), and most known pairs involve infinite matter distribution (e.g. Quan 1993; Binney & Tremaine 1987; Cuddeford 1993; Earn 1996), whereas astrophysical disks are finite (in size and mass).

For flat (i.e. zero thickness) and finite disks of interest here, the potential is mainly accessible through the double integration of the surface density weighted by the Green function $1/\vert\vec{r}-\vec{r}'\vert$ (see Binney & Tremaine 1987, for a review). This integral approach is very tricky due to the inevitable occurrence of singularities everywhere inside matter, i.e. when $\vec{r} \rightarrow \vec{r}'$. The difficulty is circumvented when the Green function is converted into series (Kellogg 1929). However, these are infinite and alternate series which do not converge inside sources, although their integration over the distribution does (Durand 1953). The efficiency of the series representation is therefore strongly misleading in practice because series truncations do not lead to a numerically stable description of the potential inside disks (Clement 1974). The errors in potential values can be as large as a few percent on average, even when the order of truncation attain several tens (Stone & Norman 1992). Singularities are no more present when the potential is expressed in terms of integrals of Bessel functions. But the oscillatory behavior of Bessel functions as well as their infinite definition domain pose technical difficulties, even for infinitely extended disks. Potential/density pairs can also be constructed by superposition of homeoids conveniently shrunk along one axis (Binney & Tremaine 1987), or by parametric expansion and ordering of the associated Poisson equation (Ciotti & Bertin 2005).

Huré & Pierens (2005) have reported a numerical method to compute the gravitational potential in flat axially symmetric disks with great accuracy by a correct treatment of the singularity. When the surface density profile is conveniently split into two components (namely, a homogeneous component plus a residual one), the contribution of the singularity to the potential can be accounted for exactly through a closed-form expression. In this paper, we explore this "splitting method'' from a purely theoretical point of view in search for fully analytical solutions (i.e. potential/surface density pairs). In particular, it is shown that the residual component is an infinite series of radial integrals which are known analytically for a certain class of surface density distributions. This is the case of power law profiles which are of great astrophysical interest (e.g. Pringle 1981). The convergence rate of the residual component is then cubic. This is a major issue, when compared with the classical multipole expansion whose converge speed is prohibitively low inside sources.

This paper is organized as follows. In Sect. 2, we briefly recall the density splitting method for one dimensional disks (Huré & Pierens 2005). We expand the residual component as an infinite series of the integrals over the modulus of the complete elliptic integral of the first kind. We then show that the construction of potential/surface density pairs is possible for a certain class of surface density distribution. The impact of the series truncation on potential values is discussed in Sect. 3. The case of surface densities scaling as a power of the radius is investigated in detail in Sect. 4 (a Fortran 90 code is available on request). For this kind of distribution, the convergence rate is cubic. This is established in Sect. 5. The paper ends with a summary and concluding remarks.

	Figure 1: A finite axially symmetric disk.
Open with DEXTER

  
2 Exact potential/surface density pairs in flat axisymmetric disks

2.1 Outline of the density splitting method

The gravitational potential $\Psi$ in the plane of a flat axially symmetric disk as pictured in Fig. 1 is given by the exact expression (e.g. Durand 1953):

$$\Psi(R) = -2G \int_{a_{\rm in}}^{a_{\rm out}}{\Sigma(a)\sqrt{\frac{a}{R}}m {\vec K}(m){\rm d}a},$$ (1)

where

$${\vec K}(m)=\int_0^{\pi/2}{\frac{{\rm d} \phi}{\sqrt{1-m^2 \sin^2 \phi}}}$$ (2)

is the complete elliptic integral of the first kind,

$$m=\frac{2\sqrt{aR}}{a+R}$$ (3)

is the modulus of ${\vec K}$ (with $0 \le m \le 1$), a and R are polar cylindrical radii (a refering to the matter distribution), $a_{\rm in}\ge 0$ is the radius of the inner edge, $a_{\rm out}> a_{\rm in}$ is the radius of the outer edge, $\Sigma(a)$ is the surface density, and G is the Gravitation constant. The function ${\vec K}$ is plotted versus m in Fig. 2; it is characterized by a divergence for $m \rightarrow 1$. This is a logarithmic divergence (see e.g. Mestel 1963; Kellogg 1929), since we have (Gradshteyn & Ryzhik 1965):

$$\lim_{m \rightarrow 1}{\vec K}(m) = \ln \frac{4}{\sqrt{1-m^2}}\left[ 1 + {\cal O}(1-m^2) \right].$$ (4)

Despite this singularity which occurs everywhere inside the disk (i.e. for $a_{\rm in}\le R \le a_{\rm out}$), the potential is generally finite everywhere.

	Figure 2: Complete elliptic integrals of the first and second kinds versus the modulus m.
Open with DEXTER

In order to determine $\Psi$ inside the disk, we set:

$$u=\frac{a}{R} \le 1 \qquad {\rm and} \qquad v = \frac{R}{a} \le 1,$$ (5)

and so, u locates matter of the inner disk $[a_{\rm in},R]$, whereas v is for the outer part $[R,a_{\rm out}]$. Using the transformation (e.g. Gradshteyn & Ryzhik 1965):

$${\vec K}\left(\frac{2\!\sqrt{x}}{1+x}\right) = (1+x) {\vec K}(x),$$ (6)

where $0 \le x \le 1$, Eq. (1) reads:

$$\Psi(R) = -4GR \left[ \int_{u_{\rm in}}^1{\Sigma(u) {\vec K}(... ..._{\rm out}}{\Sigma(v) \frac{{\vec K}(v)}{v^2}{\rm d}v}\right],$$ (7)

where $u_{\rm in}= a_{\rm in}/R \le 1$ and $v_{\rm out}= R/a_{\rm out}\le 1$ correspond to the position of the disk edges. As quoted, $\Sigma$ must now be regarded as a function of u in the first integral, and as a function of v=1/u in the second one. According to the "density splitting'' method (Huré & Pierens 2005), we divide $\Sigma$ into two components, namely:

$$\Sigma(a) = \Sigma_0+\delta \Sigma(a,R),$$ (8)

where $\Sigma_0 \equiv \Sigma(R)$ is the value of the surface density at the radius of the singularity (where the moduli m, u and v are unity), and $\delta \Sigma$ is the "residual'' surface density profile. The potential then takes the form:

$$\Psi(R) = \Psi_0(R) + \delta \Psi(R),$$ (9)

where

$$\Psi_0(R) = -4GR \Sigma_0 \left[ \int_{u_{\rm in}}^1{{\vec K}... ... \int_1^{v_{\rm out}}{\frac{{\vec K}(v)}{v^2}{\rm d}v} \right]$$ (10)

is the potential due to a radially homogeneous disk with constant surface density $\Sigma_0$, and

 

$$\delta \Psi(R) = - 4GR \Bigg[ \int_{u_{\rm in}}^1{\delta \Sigma(u... ... - \int_1^{v_{\rm out}}{\delta \Sigma(v)\frac{{\vec K}(v)}{v^2}{\rm d}v} \Bigg]$$ (11)

is the "residual'' potential associated with $\delta \Sigma$. Note that $\Sigma_0$ varies from place to place in the disk, as well as $\delta \Sigma$. As demonstrated elsewhere (see also Huré 2005), this technique of density splitting can be numerically very accurate, for two reasons. First, $\Psi_0$ is known in a closed form:

$$\Psi_0(R) = -4GR \Sigma_0 \left[ \frac{{\vec E}(v_{\rm out})}... ...E}(u_{\rm in}) + {u_{\rm in}'}^2 {\vec K}(u_{\rm in}) \right],$$ (12)

and is finite everywhere in the disk mid-plane. In Eq. (12), ${u_{\rm in}'}~=~\sqrt{1-u_{\rm in}^2}$ is the complementary modulus and ${\vec E}$ is the complete elliptic integral of the second kind (see Fig. 2):

$${\vec E}(m)=\int_0^{\pi/2}{{\rm d}\phi \sqrt{1-m^2 \sin^2 \phi}}.$$ (13)

Second, each integrand in Eq. (11) now vanishes at the singularity (this is the role of the splitting), making the two integrals easily computable with standard quadrature rules.

2.2 The residual potential from series representation

In order to determine analytical expressions for $\Psi$ on the basis of the splitting method, we consider the expansion of the elliptic integral of the first kind over its modulus (e.g. Gradshteyn & Ryzhik 1965), namely:

$$\left\{ \begin{array}{l} {\vec K}(m)= \frac{\pi}{2} \sum_{n=0... ...t(\frac{2n-1}{2n}\right)^2, \qquad n \ge 1. \end{array}\right.$$ (14)

When this infinite series is inserted into Eq. (11), the residual potential $\delta \Psi$ becomes an infinite series. Let

$$S_N(x_n) = \sum_{n=0}^N{x_n},$$ (15)

be the sum of N+1 terms $x_0, x_1,\dots, x_N$. The residual potential is then given by:

$$\delta \Psi = S_\infty \left(\delta \Psi_n \right),$$ (16)

with for any $n \ge 0$

 

$$\delta \Psi_n = - 2 \pi G R \gamma_n \Bigg[ \int_{u_{\rm in}}^1{\... ...2n+1}{\rm d}u} - \int_1^{v_{\rm out}}{\delta \Sigma(v)v^{2n-2}{\rm d}v} \Bigg].$$ (17)

Note that this expression for $\delta \Psi$ is still exact, provided all the terms $\delta \Psi _n$ of the series are included. It follows that potential/surface density pairs $(\Psi_0 + \delta \Psi,\Sigma_0+\delta \Sigma)$ can easily be constructed for flat disks with finite size (or not, depending on  $u_{\rm in}$ and $v_{\rm out}$). These pairs are analytical if the residual pair ( $\delta \Psi, \delta \Sigma)$ is analytical, that is, if the product of $\Sigma(a)$ by a power law of a is analytically integrable. We immediately see that many simple functions $\Sigma(a)$ are concerned. In particular, if $\Sigma$ is a power law of the radius a - as often met in astrophysical disk models -, then $\delta \Psi _n$ is a power law of the radius too (more precisely, a mixture of four power laws; see Sect. 4). Obviously, pairs determined through Eqs. (16) and (17) should generally remain in the form of an infinite series with no equivalent closed form expression. This is not a big problem, first because working with infinite expansions is common in potential theory, and second because the residual potential  $\delta \Psi$ is a rapidly converging series (see Sect. 5).

  
3 Approximate potential from truncated series

3.1 General remark

The estimate of $\Psi_0$ should generate no significant error since special functions ${\vec E}$ and ${\vec K}$ can be determined within the computer precision from a numerical library. In contrast, the computation of  $\delta \Psi$ is not errorless because, in practice, any series is necessarily truncated at a certain order N. Does the truncation produce large errors in potential values? In order to answer this question, it is important to stress again that the occurrence of the condition  $\delta \Sigma = 0$ enables the neutralization of the logarithmic singularity, making the remaining integrands fully regular. As a consequence, the value of the residual potential $\delta \Psi$ is mainly determined by the shape of the integrands far from the singularity, that is close to the edges. In other words, neither the precise behavior of the complete elliptic integral ${\vec K}$ nor the shape of the residual surface density profile $\delta \Sigma$ around the singularity are critical to estimate $\delta \Psi$ with a certain accuracy. We therefore expect that  $\delta \Psi$ can be accurately determined if ${\vec K}$ is replaced by any approximation which preserves at best the values of this function around the bounds  $u_{\rm in}$ and $v_{\rm out}$ rather than at the radius of the singularity (where m=u=v=1). From this point of view, the series defined by Eq. (14) is fully appropriate as it converges very rapidly for moduli significantly less than unity. Figure 3 shows a typical example where $\Sigma$ is a power law of the radius. We clearly see that the two integrands $\delta \Sigma(u) u {\vec K}(u)$ and $\delta \Sigma(v){\vec K}(v)/v^2$ appearing in Eq. (11) are quite insensitive to the approximation adopted for ${\vec K}$ (one-term, two-term expansions or exact). This conclusion is independent of the $\Sigma$-profile and of R as well.

	Figure 3: Typical shape of the integrands appearing in Eq. (11) computed from different approximations for the function ${\vec K}$. In this example, the disk has inner edge $a_{\rm in}=0.1$ and outer edge $a_{\rm out}=1$, the surface density obeys a power law of the radius with $\Sigma (a) \propto a^{-3/2}$ (that is s=-1.5), and the radius is $R=\frac{1}{2}(a_{\rm in}+a_{\rm out})$.
Open with DEXTER

3.2 The first terms

For a truncation of the series at order N, the potential is given by the approximate value:

$$\Psi_{\rm approx} \equiv \Psi_0+S_N(\delta \Psi_n) \ne \Psi,$$ (18)

and we have asymptotically:

$$\Psi_0 + \lim_{N \rightarrow \infty}{\left[ S_N ( \delta \Psi_n ) \right]} = \Psi.$$ (19)

With only one term in Eq. (14), we have $S_0(\delta \Psi_n) = \delta \Psi_0$ with:

$$\delta \Psi_0 = -2 \pi GR \left[ \int_{u_{\rm in}}^1{\delta ... ...1^{v_{\rm out}}{\delta \Sigma(v)\frac{{\rm d}v}{v^2}} \right],$$ (20)

and so the approximate potential is $\Psi_{\rm approx} = \Psi_0 + \delta \Psi_0$. For a two-term expansion of the complete elliptic integral of the first kind, the residual potential is $S_1(\delta \Psi_n) = \delta \Psi_0+\delta \Psi_1$, where $\delta \Psi _0$ is given by Eq. (20), and:

$$\delta \Psi_1 = -\frac{\pi}{2} G R \left[ \int_{u_{\rm in}}^1... ...d}u} - \int_1^{v_{\rm out}}{\delta \Sigma(v){\rm d}v} \right].$$ (21)

The approximate potential is then $\Psi_{\rm approx.} = \Psi_0+\delta \Psi_0 +\delta \Psi_1$, and so on. Each new term added into the expansion brings a new component $\delta \Psi _n$ defined by Eq. (17). We then expect that the more terms in the expansion, the more accurate the potential. If the series converges rapidly (this is the case of power law distributions; see below), the absolute error made on $\Psi$ is roughly $\vert\delta \Psi_{N+1}\vert$.

  
4 The case of power-law surface density profiles

Let us consider the case where $\Sigma$ is a power-law of the radius, namely:

$$\Sigma = \sigma_0 \left(\frac{a}{a_{\rm out}} \right)^s$$ (22)

where $\sigma_0$ is the surface density at the outer edge of the disk. Power-law surface density solutions are found in many disk models (e.g. Pringle 1981) and match observations (e.g. Guilloteau & Dutrey 1998), with s<0 in general. We then have:

$$\Sigma_0 = \sigma_0 \left(\frac{R}{a_{\rm out}} \right)^s$$ (23)

and so:

$$\left\{ \begin{array}{l} \delta \Sigma = \Sigma_0 \left( u^s ... ...s}}-1\right) \quad {\rm for} \quad v \le 1. \end{array}\right.$$ (24)

	Figure 4: Top: potential ( lines) inside the disk for three different surface density exponents s when computed with a one-term expanded elliptic integral (i.e. $\delta \Psi \equiv \delta \Psi _0$). Reference values $\Psi _{\rm ref}$ ( symbols) are obtained numerically within the computer precision. Bottom: decimal logarithm of the relative error (the average of $\epsilon$ is also indicated).
Open with DEXTER

4.1 Order zero

For a one-term expansion of the complete elliptic integral of the first kind, the residual potential is given by:

$$\delta \Psi_0 = 2\pi G R \Sigma_0 \left[ \int_{u_{\rm in}}^1... ... out}}{\left(v^{-s}-1\right)\frac{{\rm d}v}{v^2}} \right]\cdot$$ (25)

If $s \ne \{-2, -1\}$, this is (see Appendix A for these two cases):

$$\delta \Psi_0 = 2 \pi G R \Sigma_0 \Bigg[\frac{u_{\rm in}^{s+2}-1... ...m in}^2}{2} +\frac{1-v_{\rm out}^{-s-1}}{s+1} + \frac{1}{v_{\rm out}}-1 \Bigg].$$ (26)

Figure 4 displays the approximate potential $\Psi_{\rm approx} = \Psi_0 + \delta \Psi_0$ for disks with inner edge $a_{\rm in}=0.1$, outer edge $a_{\rm out}=1$ and power law exponents s=-1, -1.5 and -2.5. References potential values $\Psi_{\rm ref}(R)$ have been determined within the computer precision from the numerical splitting method (Huré & Pierens 2005) and compared to $\Psi_{\rm approx}$. The decimal logarithm of the relative error $\epsilon \equiv \vert 1- \Psi_{\rm approx}/\Psi_{\rm ref}\vert$ is also shown in Fig. 4. We see that the potential is reproduced with an accuracy better than $\sim$2% (on average) in the entire disk with only one term in the expansion of the ${\vec K}$-function. This is already remarkable. The accuracy is relatively uniform, and it is slightly better at the outer edge. This point is particularly interesting to model young stellar objects or Active Galactic Nuclei. In these systems, the inner disk is dynamically under the control of a massive central object (a proto-star or a black hole), and so there is no need for very accurate potential values. On the other hand, in the outer regions where disks can be self-gravitating (e.g. Huré 2000), a reliable description of the potential is required.

4.2 First order

The accuracy of potential values can be improved by adding a second term. For N=1, we find:

$$\delta \Psi_1 = \frac{\pi}{2} G R \Sigma_0 \left[ \int_{u_{\r... ...+ \int_1^{v_{\rm out}}{\left(v^{-s}-1\right){\rm d}v} \right].$$ (27)

If $s \ne \{-4, 1\}$ (see Appendix B for these two cases), the expression is:

$$\delta \Psi_1 = 2 \pi G R \Sigma_0 \Bigg[ \frac{u_{\rm in}^{s+4}-... ... + \frac{1- v_{\rm out}^{-s+1}}{4(s-1)} + \frac{1 - v_{\rm out}}{4} \Bigg]\cdot$$ (28)

Figure 5 displays the results obtained with a two-term expansion of ${\vec K}$, i.e. $\Psi_{\rm approx} = \Psi_0 + \delta \Psi_0 + \delta \Psi_1$, under the same conditions as for Fig. 4. We see that the relative accuracy of mid-plane potential values (of the order of $0.3\%$ on average) is now better by a factor of $\sim$5.

	Figure 5: Same legend as for Fig. 4 but the residual potential $\delta \Psi$ is determined from a two-term expanion of the complete elliptic integral of the first kind (i.e. $\delta \Psi \equiv \delta \Psi _0+ \delta \Psi _1$).
Open with DEXTER

Such a level of accuracy is probably sufficient for many disk models and applications. However, since the process of adding more and more terms represents no technical difficulty and converges rapidly (see Sect. 5), it can be repeated until a given accuracy is reached. In particular, for $s+2n+2\ne0$ and $s-2n+1 \ne 0$, we have (see the Appendix C for these two cases)

 

$$\delta \Psi_n = 2 \pi G \Sigma_0 R \gamma_n \Bigg[\frac{1-u_{\rm ... ...v_{\rm out}^{2n-s-1}-1}{2n-s-1} - \frac{v_{\rm out}^{2n-1}-1}{2n-1} \Bigg]\cdot$$ (29)

	Figure 6: Decimal logarithm of the relative error on potential values $\Psi =\Psi _0 + \delta \Psi$ versus the order N of the truncation in Eq. (19) for five power law exponents s=-3.5, -3, -2.5, -1.5 and -0.5. The disk has inner edge $a_{\rm in}=0.1$, outer edge  $a_{\rm out}=1$ and the potential is computed at $R=\frac{1}{2}(a_{\rm in}+a_{\rm out})$.
Open with DEXTER

  
5 Convergence rate of the series $\delta \Psi = S_\infty (\delta \Psi _{\textbf {n}})$ in the case of power law distributions

5.1 A cubic convergence ?

Figure 6 shows the rise of the accuracy of potential values $\Psi(R)$ as the order N of the truncation increases for the same disk parameters as already considered and a few exponents s. We find that the relative error $\epsilon \equiv \vert 1- \Psi_{\rm approx}/\Psi_{\rm ref}\vert$ first decreases slowly as N increases. But above a certain rank $n_0 \sim 8$-10, this error approximately scales as 1/N³, before saturation near the computer precision, for $N \gtrsim 10^4$ terms. We have checked that this conclusion holds for a wide variety of configurations and radii R. Convergence appears to be cubic. With this approach, we can easily select the lowest order corresponding to a given level of accuracy. From Fig. 6, we have roughly:

$$\log \epsilon \approx - 3 \log N -2 \pm 0.5,$$ (30)

for $N \gtrsim n_0$, or $N \approx \left(100 \epsilon \right)^{-1/3}$ within about $0.17 \%$.

There are two remarkable cases. For s=0, the potential is fully analytical: $\Psi=\Psi_0$ (i.e. $\delta \Psi = 0$). For s=-3, the relative accuracy rises with N much more drastically than for others exponents and reaches the maximum allowed by the computer for $N \gtrsim n_0$ (see below).

	Figure 7: Left: $\gamma _n$ versus n. Right: $\zeta _n$ and $\tau _n$ (for s=-1.5) versus n.
Open with DEXTER

	Figure 8: Left: series $S_N(\gamma _n)$ versus the truncation order N which diverges as $N \rightarrow \infty$. Right: series $\pi S_N(\zeta _n)/4$ versus the truncation order N which converges towards Catalan's constant (see also Fig. 7).
Open with DEXTER

5.2 Proof

The expansion of ${\vec K}(m)$ through Eq. (14) converges slowly as m approaches unity, and even diverges for m=1, although $\gamma _n$ decreases rapidly as n increases (for large n, $\gamma_n \sim 1/n$). This is shown in Figs. 7 and 8. Note that the divergence of ${\vec K}$ as $m \rightarrow 1$ is the main reason why the determination of the potential inside sources by direct integration is generally not the appropriate way. What about the convergence of $\delta \Psi$? Since $\delta \Psi _n$ comes from the integration of ${\vec K}$ over its modulus, $S_N(\delta \Psi_n)$ is made of terms of the form $\frac{\gamma_n}{2n+1} \equiv \zeta_n$ which decrease faster with n than terms $\gamma _n$ do. The series $\delta \Psi$ is then expected to converge (see the Appendix D). Figure 7 displays the terms $\zeta _n$ versus n, and the series $S_N (\zeta_n)$ is plotted versus N in Fig. 8.

When the series $\delta \Psi$ is truncated at order N, the relative error $\epsilon(N)$ on the potential value is:

$$\epsilon(N) = \left\vert \frac{\Psi_{\rm approx} - \Psi}{\Psi} \right\vert\cdot$$ (31)

From Eqs. (18) and (19), this expression can be re-written in the following form:

$$\epsilon(N) = \left\vert \frac{S_N(\delta \Psi_n) - S_\infty(\delta \Psi_n)}{\Psi}\right\vert\cdot$$ (32)

The quantities $u_{\rm in}$ and $v_{\rm out}$ being always smaller than unity, their contribution in Eq. (29) - an exponential drop - is expected to vanish rapidly as n exceeds a certain rank. This transition occurs at $N~\approx~n_0$; it is clearly visible in Fig. 6 where $n_0 \sim 8$. Thus, for $N \gtrsim n_0$, the leading terms in $\delta \Psi$ form a series $S_N(\tau_n)$, with:

$$\tau_n = \gamma_n \Bigg[ \frac{1}{2(n+1)} - \frac{1}{2(n+1)+s} +\frac{1}{2n-1} - \frac{1}{2n-1-s}\Bigg],$$ (33)

 which can easily be rearranged:

$$\tau_n = - \gamma_n \frac{s(s+3)(4n+1)}{(2n+2)(2n+2+s)(2n-1)(2n-1-s)}\cdot$$ (34)

For large values of n, we have:

$$\tau_n \approx -s(s+3) \frac{\gamma_n}{4 n^3} \left(1 - \frac{3}{4n} + \dots \right).$$ (35)

and so, $\tau _n$ behaves like n^-4, asymptotically. The coefficient $\tau _n$ is displayed versus n for s=-1.5 in Fig. 7. Then, as soon as $N \gg n_0$, the relative error is

$$\epsilon(N) = \frac{2 \pi G \Sigma_0 R}{\Psi} \left(\sum_{n_0}^N{\tau_n} - \sum_{n_0}^\infty{\tau_n} \right).$$ (36)

Using Eq. (35), we have:

$$\epsilon(N) \approx \frac{2 \pi G \Sigma_0 R}{\Psi} \int_\infty^N{\tau_n {\rm d}n}$$

$$\approx - \frac{2 \pi G \Sigma_0 R s(s+3)}{\Psi} \int_\infty^N{\frac{\gamma_n}{4n^3} {\rm d}n}.$$ (37)

Since $\gamma _n$ scales as 1/n for large n (see also Fig. (7)), we have:

$$\epsilon(N) \propto - \frac{\Sigma_0 R s(s+3)}{\Psi} \frac{1}{N^3},$$ (38)

as found "experimentally''. The convergence rate of this series is therefore cubic. This result holds irrespective of the radius R (except at the disk edges) and whatever s, except in the case of the homogeneous disk (i.e. for s=0) and for s=-3 as well. In this latter case, the convergence rate is not limited by the terms $\tau _n$, but by the series asssociated with  $u_{\rm in}$ and $v_{\rm out}$. For s=-3, we have:

$$\delta \Psi_n = 2 \pi G \Sigma_0 R \gamma_n \left[\frac{v_{\r... ...}-\frac{v_{\rm out}^{2n-1} - u_{\rm in}^{2n-1}}{2n-1} \right],$$ (39)

and so the convergence of $\delta \Psi$ is exponential (see Fig. 6), except at the inner and outer edges (see Appendix E).

All numerical experiments have succesfully illustrated these results. This is a major issue, in particular because classical multipole expansions converge very slowly inside sources.

6 Summary and concluding remarks

In this paper, we have reported a method to determine exact and approximate potential/surface density pairs in flat finite size axially symmetrical disks, in the continuity of the numerical splitting method described in Pierens & Huré (2005). At a given radius R, each pair is the sum of two components:

• a "homogeneous'' pair $(\Psi_0,\Sigma_0)$ where $\Sigma_0= \Sigma(R)$ is the local value of the surface density and $\Psi_0$ is given by Eq. (12),
• a "residual'' pair $(\delta \Psi, \delta \Sigma)$, where $\delta \Psi$ is an infinite series (see Eqs. (16) and (17)) and $\delta \Sigma(a) = \Sigma(a) - \Sigma_0$.

In general, the infinite series $\delta \Psi$ cannot be put into a closed form expression. However, in contrast to the multipole expansions commonly met in potential theory, this series has a cubic convergence rate inside the source. This major result is due to the splitting method which is equivalent to a convergence acceleration technique (in the mathematical sense). It means that very good approximations for the potential can be obtained by low order truncations of the series. Pairs can be constructed only for a certain class of the surface density profile $\Sigma(a)$: any function $\Sigma(a)$ whose product by a power law of the radius a is analytically integrable does lead to a potential/surface density pair. This is the case of power laws. It would be interesting to list all functions satisfying this criterion.

We have fully considered power law surface densities. Such profiles are particularly attractive for two main reasons: i) power laws underly many (if not most) astrophysical disk models, and ii) power laws may serve as a basis set to construct more complicated surface density profiles. This new approach along with associated computing tools should enable a much better representation of the gravitational potentials and forces in flat disks.

Acknowledgements

We thank J. Braine for many fruitful comments. We are especially grateful to the referee L. Ciotti whose suggestions and criticisms allowed significant improvement of the paper.

  
Appendix A: Special cases for $\delta \Psi _0$

For s=-1 in Eq. (22), then Eq. (20) reads

$$\delta \Psi_0= 2 \pi G \Sigma_0 R\left\{ \left[\frac{u^2}{2}-... ...+ \left[ \ln v + \frac{1}{v} \right]_1^{v_{\rm out}} \right\}.$$ (A.1)

For s=-2 in Eq. (22), then Eq. (20) reads

$$\delta \Psi_0= 2 \pi G \Sigma_0 R\left\{ \left[\frac{u^2}{2}-... ...}^1 + \left[ v + \frac{1}{v} \right]_1^{v_{\rm out}} \right\}.$$ (A.2)

  
Appendix B: Special cases for $\delta \Psi _1$

For s=-4 in Eq. (22), then Eq. (21) reads

$$\delta \Psi_1= \frac{\pi}{2} G \Sigma_0 R\left\{ \left[\frac{... ...1 + \left[ \frac{v^5}{5} - v \right]_1^{v_{\rm out}} \right\}.$$ (B.1)

For s=+1 in Eq. (22), then Eq. (21) reads

$$\delta \Psi_1= \frac{\pi}{2} G \Sigma_0 R\left\{\left[\frac{u... ...rm in}}^1 + \left[ \ln v - v \right]_1^{v_{\rm out}} \right\}.$$ (B.2)

  
Appendix C: Special cases for $\delta \Psi _n$

For s+2n+2=0 in Eq. (22), then $s-2n+1 \ne 0$ and Eq. (17) reads

$$\delta \Psi_n = 2 \pi G \Sigma_0 R \gamma_n \Bigg\{ \left[\frac{u... ...ac{v^{2n-s-1}}{2n-s-1} - \frac{v^{2n-1}}{2n-1} \right]_1^{v_{\rm out}} \Bigg\}.$$ (C.1)

For $s+2n+2\ne0$ and s-2n+1=0 in Eq. (22), then Eq. (17) reads

$$\delta \Psi_n = 2 \pi G \Sigma_0 R \gamma_n \Bigg\{ \left[\frac{u... ... in}}^1 + \left[ \ln v - \frac{v^{2n-1}}{2n-1} \right]_1^{v_{\rm out}} \Bigg\}.$$ (C.2)

  
Appendix D: Convergence of S $_\infty (\zeta _n)$

From (e.g. Gradshteyn & Ryzhik 1965), we see that

$$S_\infty (\zeta_n) \sum_{n=0}^\infty{\frac{\gamma_n}{2n+1}},$$ =

$$\sum_{n=0}^\infty{\int_0^1{\gamma_n m^{2n} {\rm d}m}},$$ =

$$\int_0^1{\sum_{n=0}^\infty{\gamma_n m^{2n} {\rm d}m}},$$ =

$$\frac{2}{\pi} \int_0^1{{\vec K}(m){\rm d}m},$$ =

$$\equiv \frac{4C}{\pi},$$ (D.1)

where C is Catalan's constant.

  
Appendix E: Convergence rate for s = -3 at the disk edges

At the disk outer edge where $v_{\rm out}=1$, we have:

$$\delta \Psi_n = 2 \pi G \Sigma_0 R \gamma_n \left[\frac{1-\Delta^{2(n+1)}}{2(n+1)}-\frac{1 - \Delta^{2n-1}}{2n-1} \right],$$ (E.1)

where $\Delta=a_{\rm in}/a_{\rm out}$ is the disk axis ratio. Since $\Delta < 1$, the leading terms in $\delta \Psi$ then form a series $S_N(\tau'_n)$, with:

 

$$\tau'_n \gamma_n \left[ \frac{1}{2(n+1)} - \frac{1}{2n-1} \right],$$ =

$$-\gamma_n \frac{3}{2(n+1)(2n-1)}\cdot$$ = (E.2)

For large n, we have $\tau'_n \sim 1/n^3$, and so $\epsilon(N) \propto 1/N^2$. The convergence is therefore quadratic (as at the disk inner edge).

References

• Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton, NJ: Princeton University Press), 747 (In the text)
• Ciotti, L., & Bertin, G. 2005, A&A, 437, 419 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Clement, M. J. 1974, ApJ, 194, 709 [NASA ADS] [CrossRef] (In the text)
• Cuddeford, P. 1993, MNRAS, 262, 1076 [NASA ADS]
• de Zeeuw, T., & Pfenniger, D. 1988, MNRAS, 235, 949 [NASA ADS]
• Durand, E. 1953, Electrostatique, vol. I, Les distributions (Masson) (In the text)
• Earn, D. J. D. 1996, ApJ, 465, 91 [NASA ADS] [CrossRef]
• Gradshteyn, I. S., & Ryzhik, I. M. 1965, Table of integrals, series and products (New York: Academic Press), ed. Yu. V. Geronimus, 4th Ed., M. Yu. Tseytlin (In the text)
• Guilloteau, S., & Dutrey, A. 1998, A&A, 339, 467 [NASA ADS] (In the text)
• Huré, J.-M. 2000, A&A, 358, 378 [NASA ADS] (In the text)
• Huré, J.-M. 2005, A&A, 434, 1 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Huré, J.-M., & Pierens, A. 2005, ApJ, 624, 289 [NASA ADS] [CrossRef] (In the text)
• Kellogg, O. D. 1929, Foundations of Potential Theory (New-York: Frederick Ungar Publishing Company) (In the text)
• Mestel, L. 1963, MNRAS, 126, 553 [NASA ADS]
• Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] (In the text)
• Quan, E. E. 1993, MNRAS, 263, 394 [NASA ADS]
• Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [NASA ADS] [CrossRef] (In the text)