A&A 490, 477-486 (2008)
DOI: 10.1051/0004-6361:200809682

A new equation for the mid-plane potential of power law discs

II. Exact solutions and approximate formulae

J.-M. Huré1,2 - F. Hersant2 - C. Carreau3 - J.-P. Busset1


1 - Université de Bordeaux, LAB, 351 cours de la Libération, Talence 33405, France
2 - CNRS/INSU, UMR 5804/LAB, 2 rue de l'Observatoire, BP 89, 33271 Floirac Cedex, France
3 - La Maurellerie, 37290 Bossay-sur-Claise, France

Received 29 February 2008 / Accepted 31 May 2008

Abstract
Aims. The first-order ordinary differential equation (ODE) that describes the mid-plane gravitational potential in flat finite size discs of surface density $\Sigma(R) \propto R^{s}$ (Huré & Hersant 2007, A&A, 467, 907) is solved exactly in terms of infinite series.
Methods. The formal solution of the ODE is derived and then converted into a series representation by expanding the elliptic integral of the first kind over its modulus before analytical integration.
Results. Inside the disc, the gravitational potential consists of three terms: a power law of radius R with index 1+s, and two infinite series of the variables R and 1/R. The convergence of the series can be accelerated, enabling the construction of reliable approximations. At the lowest-order, the potential inside large astrophysical discs ( $s \sim -1.5 \pm 1$) is described by a very simple formula whose accuracy (a few percent typically) is easily increased by considering successive orders through a recurrence. A basic algorithm is given.
Conclusions. Applications concern all theoretical models and numerical simulations where the influence of disc gravity must be checked and/or reliably taken into account.

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

1 Introduction

Gaseous discs in which the main physical quantities (density, pressure, temperature, thickness, velocity) scale with cylindrical radius as power laws, i.e. ``power-law discs'', represent an important class of theoretical systems. These are used customary to model accretion in evolved binaries (Shakura & Sunyaev 1973; Pringle 1981), circumstellar matter (Edgar 2007; Dubrulle 1992), the environment of massive black holes inside active galactic nuclei (Collin-Souffrin & Dumont 1990; Huré 1998; Semerák 2004) or even the stellar component of some galaxies (Evans 1994; Zhao et al. 1999). In most applications however, power-law discs are truncated either to avoid diverging values at the disc centre (such as density, mass) or in attempting to reproduce the properties of observed discs of finite extension and mass. Although self-similarity is not compatible with the presence of edges, it is generally considered that power laws offer a good description of disc properties in some regions (far from the edges). Note that the presence of sharp edges can be misleading when interpreting observational data (e.g. Hughes et al. 2008). In general, the surface density $\Sigma$ in the outer parts of discs is a decreasing function of the cylindrical radius R. Depending on the models, hypotheses, and objects, we have, for instance, $\Sigma \propto R^{-3/5}$ in binaries (Shakura & Sunyaev 1973), $\Sigma \propto R^{\pm9/20}$ in active galactic nuclei (Collin-Souffrin & Dumont 1990), $\Sigma \propto R^{-1}$ for a Mestel disc (Mestel 1963), or $\Sigma \propto R^{-3/2}$ in circumstellar discs (Piétu et al. 2007). In the context of stationary viscous $\alpha $-discs, a wide range of power-law exponents is allowed since the temperature T, $\Sigma$, and R satisfy the condition (e.g. Pringle 1981):

\begin{displaymath}\Sigma T R^{3/2} ={\rm const}.,
\end{displaymath}

while it is $\Sigma \propto R^{-1/2}$ in $\beta $-discs (Huré et al. 2001).

The calculus of the gravitation potential of finite-size, power-law discs has received little attention yet. Several reasons can be put forward. Solving the Poisson equation or computing the integral of the potential is not a trivial procedure, especially in the presence of edges. It is generally believed that gravity due to low mass discs is unimportant compared with that of a central proto-star or black hole, and cannot be probed (see however Baruteau & Masset 2008). Many studies employ the multi-pole expansion which is known to converge too slowly inside sources to be efficient for the numerical applications (e.g. Stone & Norman 1992; Clement 1974). Huré & Hersant (2007) demonstrated that the mid-plane potential of flat power-law discs obeys an inhomogeneous first-order ordinary differential equation (ODE). In this second paper, we discuss the exact solutions of this ODE for the entire physical range (outside and inside the disc) in terms of infinite series. In particular, it is shown that the mid-plane potential is a combination of a power law for the radius R and two series of the variables R and 1/R. Since these series converge rapidly inside large discs, it is possible to derive reliable approximations by truncating the series at low orders.

This paper is organised as follows. The ODE for the potential is briefly recalled in Sect. 2 and its formal solution is derived in Sect. 3. In Sect. 4, we express the potential at the two disc edges and consider a few special cases. The inside and outside solutions in the form of series are presented in Sect. 5. In Sect. 6, we analyse the potential in the disc inside in detail, and in particular, the power-law contribution. Since all series involved converge rapidly, we are able to derive reliable approximations for the potential; this is done in Sect. 7. We discuss in Sect. 8 the case of discs with no inner and/or outer edge. The paper ends with a few concluding remarks.


  \begin{figure}
\par\includegraphics[width=9cm,clip]{fig1hure09682.eps}
\end{figure} Figure 1: Configuration for a finite-size flat disc.
Open with DEXTER

   
2 The mid-plane potential in power-law discs from an ordinary differential equation (ODE)

Following Huré & Hersant (2007) hereafter Paper I), the mid-plane potential $\psi $ due to a flat power-law disc satisfies the ODE:

 \begin{displaymath}%
\frac{{\rm d} \psi}{{\rm d} \varpi} - (1+s)\frac{\psi}{\varpi} = S(\varpi),
\end{displaymath} (1)

where $\varpi=R/a_{\rm out}$ is the cylindrical radius in units of the radius $a_{\rm out}$ of the disc outer edge, s is the power-law index of the surface density $\Sigma$, namely (generally, s < 0 in astrophysical discs):

\begin{displaymath}%
\Sigma \propto R^s,
\end{displaymath} (2)

and $S(\varpi)$ is the piecewise defined function. Depending on the position in the disc (see Fig. 1), we have:

 \begin{displaymath}%
S(\varpi) = \frac{ 2 \psi_{\rm out}}{\pi \varpi^2}
\left\{ ...
...uad \varpi \ge 1 \qquad {\rm (region~III)},
\end{array}\right.
\end{displaymath} (3)

where $\Delta=a_{\rm in}/a_{\rm out}< 1$ is the axis ratio, $a_{\rm in}$ is the radius of the inner edge, $\psi_{\rm out}=2 \pi G \Sigma_{\rm out}a_{\rm out}$ is a positive constant, $\Sigma_{\rm out}$ is the surface density at the disc outer edge, ${\vec K}$ is the complete elliptic integral of the first kind:

\begin{displaymath}{\vec K}(x)=\int_0^{\pi/2}{\frac{{\rm d} \phi}{\sqrt{1-x^2 \sin^2 \phi}}}, \quad 0 \le x \le 1,
\end{displaymath} (4)

and G is the gravitation constant. The above ODE can in principle be solved in the entire radial domain since boundary conditions (both $\psi $ and the associated acceleration $-{\rm d}_\varpi \psi$) are known precisely at $\varpi=0$ and at $\varpi=\infty$ for power law distributions. Although S is singular at the two disc edges (i.e. for $\varpi \in \{\Delta,1\}$), Eq. (1) is far more tractable in computing $\psi $ than the integral form (e.g. Durand 1953):

 \begin{displaymath}%
\psi(\varpi) = -4 G \Sigma_{\rm out}a_{\rm out}^{-s} \int_{...
...1}}{a+R} {\vec K}\left(\frac{2\sqrt{aR}}{a+R}\right){\rm d}a},
\end{displaymath} (5)

whose integrand is logarithmically singular everywhere inside the disc, or than the Poisson equation, which involves vertical gradients.

   
3 Formal solution of the ODE

A formal solution of Eq. (1) is found by setting (e.g. Rybicki & Lightman 1979):

 \begin{displaymath}%
\bar{\psi}(\varpi) = \varpi^{-(1+s)} \psi(\varpi),
\end{displaymath} (6)

and

 \begin{displaymath}%
\bar{S}(\varpi) = \varpi^{-(1+s)} S(\varpi).
\end{displaymath} (7)

The exact derivative of $\bar{\psi}$ is:

\begin{displaymath}%
\frac{{\rm d} \bar{\psi}}{{\rm d} \varpi} = \varpi^{-(1+s)}...
...\rm d} \psi}{{\rm d} \varpi} - (1+s)\frac{\psi}{\varpi}\right]
\end{displaymath} (8)

where we recognise, inside brackets, the function S. Therefore, we have:

\begin{displaymath}%
\frac{{\rm d} \bar{\psi}}{{\rm d} \varpi} = \bar{S}(\varpi),
\end{displaymath} (9)

whose formal solution for $\bar{\psi}$ is of the form:

\begin{displaymath}%
\bar{\psi}(\varpi) = \bar{\psi}(\varpi_0) + \int_{\varpi_0}^\varpi{\bar{S}(\varpi'){\rm d}\varpi'},
\end{displaymath} (10)

Back-substituting $\bar{\psi}$ and $\bar{S}$ from Eqs. (6) and (7), we find the general expression for the mid-plane potential:

 \begin{displaymath}%
\psi(\varpi) = \varpi^{1+s} \left[ \frac{\psi(\varpi_0)}{\v...
...rpi{\frac{S(\varpi')}{{\varpi'}^{1+s}}{\rm d}\varpi'} \right].
\end{displaymath} (11)

This solution is fully determined in the entire spatial domain or part of it as soon as the potential is known at a given normalised radius $\varpi_0$. We observe that, for $s \ne -1$, $\psi(\varpi)$ is the mixture of a power law of the radius (the first term in the right-hand-side) with exponent (e.g. Bisnovatyi-Kogan 1975; Evans & Read 1998):

 
p = 1+s, (12)

and a complicated function of the radius R (the definite integral).

In the following, we shall analyse Eq. (11) analytically in terms of infinite series by considering in Eq. (3) the expansion of ${\vec K}$ over its modulus (e.g. Gradshteyn & Ryzhik 1965):

 \begin{displaymath}%
\left\{
\begin{array}{l}
{\vec K}(x)= \frac{\pi}{2} \sum_{n...
...t(\frac{2n-1}{2n}\right)^2, \qquad n \ge 1.
\end{array}\right.
\end{displaymath} (13)

   
4 Potential at the disc edges

4.1 Inner edge

To determine $\psi(\varpi)$ from Eq. (11) for $\varpi \in [0,\infty[$, it is sufficient to calculate the potential at the two disc edges[*], that is $\psi(\Delta)$ and $\psi(1)$. For this purpose, we use Eq. (5) and define $v = R/a \le 1$. After some algebra[*] we find:

 \begin{displaymath}%
\psi(\Delta) = 4G \Sigma_{\rm out}\Delta^{1+s} a_{\rm out}\int_1^{\Delta}{\frac{{\vec K}(v)}{v^{2+s}}{\rm d}v}.
\end{displaymath} (15)

From Eq. (13), the potential at the inner edge is:

 \begin{displaymath}%
\psi(\Delta) = -\psi_{\rm out}\Delta^{1+s} \sum_{n=0}^\infty{I_n \left(1-\Delta^{2n-s-1}\right)},
\end{displaymath} (16)

where

\begin{displaymath}%
I_n = \frac{\gamma_n}{2n-s-1},
\end{displaymath} (17)

and assuming $\Delta \ne 0$ (see below). As $\gamma_n \sim$1/n for large n (see Fig. 2), $\psi(\Delta)$ consists of terms that vary asymptotically like $\sim$1/n2. Since $\Delta < 1$, $\psi(\Delta)$ is a converging series.


  \begin{figure}
\par\includegraphics[width=9cm,clip]{fig2hure09682.eps}\end{figure} Figure 2: Coefficient $\gamma _n$ versus n. Terms an and |dn| versus n for s=-1.5.
Open with DEXTER

4.2 Outer edge

At the outer edge, the potential is calculated in a similar manner, but using the variable $u = a/R \le 1$. For $\varpi=1$, Eq. (5) writes1:

 \begin{displaymath}%
\psi(1) = - 4G \Sigma_{\rm out}a_{\rm out}\int_\Delta^1{{\vec K}(u)u^{1+s}{\rm d}u}.
\end{displaymath} (18)

Replacing ${\vec K}(u)$ by its series representation yields:

 \begin{displaymath}%
\psi(1) = - \psi_{\rm out}\sum_{n=0}^\infty{J_n \left(1 - \Delta^{2n+s+2}\right)},
\end{displaymath} (19)

where

\begin{displaymath}%
J_n = \frac{\gamma_n}{2n+s+2},
\end{displaymath} (20)

and by assuming $\Delta \ne 0$ (see below). As for $\psi(\Delta)$ and for the same reasons, $\psi(1)$ is also a converging series.

   
4.3 Special values of s

If the power law exponent s is such that:

 
2n-s-1=0 (21)

at a certain rank $n \equiv n_\Delta$, then $\psi(\Delta)$ must be, in practice, written in a slightly different form. This happens for $s \in {\cal E}_\Delta$ with:

\begin{displaymath}%
{\cal E}_\Delta = \{-1,+1,+3,+5,\dots \}.
\end{displaymath} (22)

Since

 \begin{displaymath}%
\lim_{q \rightarrow 0 }\frac{1 - x^q}{q} = -\ln x,
\end{displaymath} (23)

for x >0 and any q, we have:

\begin{displaymath}%
\lim_{n \rightarrow n_\Delta }I_n(1 - \Delta^{2n-s-1}) = -\gamma_{n_\Delta}\ln \Delta.
\end{displaymath} (24)

Then, if $s \in {\cal E}_\Delta$, the potential at the disc inner edge is given by:

 \begin{displaymath}
%
\psi(\Delta) = -\psi_{\rm out} \Delta^{1+s} \left[
\sum_{...
...Delta^{2n-s-1}\right)} - \gamma_{n_\Delta} \ln \Delta \right].
\end{displaymath} (25)

In a similar way, if s is such that:

 
2n+s+2=0 (26)

at a rank $n \equiv n_1$, then one term in Eq. (19) must be treated separately. This happens for $s \in {\cal E}_1$ where:

\begin{displaymath}%
{\cal E}_1 = \{-2,-4,-6, -8, \dots\}.
\end{displaymath} (27)

From Eq. (23), we have:

 \begin{displaymath}%
\psi(1) = -\psi_{\rm out}\left[ \sum_{\begin{array}{l}
{\sc...
...t(1-\Delta^{2n+s+2}\right) - \gamma_{n_1} \ln \Delta} \right].
\end{displaymath} (28)

We note that s cannot belong simultaneously to set ${\cal E}_\Delta$ and to set  ${\cal E}_1$.

4.4 Cases with $\Delta =0$

The case $\Delta =0$ occurs i) when the disc has no inner hole (i.e. $a_{\rm in}=0$) but finite size, and/or ii) when the disc has an inner edge but is infinitely extended (i.e. $a_{\rm out}\rightarrow \infty$). In the first case, $\psi(\Delta)$ (denoted $\psi_{\rm c}$ in Paper I) becomes the potential at the origin of coordinates and has infinite value as soon as $1+s\le 0$. In the second case, $\psi(1)$ represents the gravitational potential at infinity. It diverges if $s+1 \ge 0$. Figure 3 summarises the ranges of s where the edge surface density, edge potential, and total disc mass are finite. We note that only discs having either i) $a_{\rm in}=0$, $a_{\rm out}\ne \infty$ with s > 0; or ii) $a_{\rm in}>0$, $a_{\rm out}\rightarrow \infty$ with $s \le -2$ are physically meaningful since they are characterised by a finite surface density, a finite total mass and a finite potential. Mestel discs do not belong to these categories (e.g. Hunter et al. 1984; Mestel 1963).


  \begin{figure}
\par\includegraphics[width=9cm,clip]{fig3hure09682.eps}
\end{figure} Figure 3: Ranges of s where the edge surface density, edge potential and total disc mass are infinite.
Open with DEXTER

   
5 Solution of the ODE in the form of series

We see from Eqs. (3) and (13) that the function S can be easily expressed as a series. We have:

 \begin{displaymath}%
S(\varpi) = \frac{\psi_{\rm out}}{\varpi^2} \sum_{n=0}^\inf...
...{s+2+2n} \right) \quad {\rm in~region~III.}
\end{array}\right.
\end{displaymath} (29)

By inserting this general expression into Eq. (11), we find for $s \notin {\cal E}_\Delta \cup {\cal E}_1$ and $\Delta \ne 0$:
 
$\displaystyle %
\psi(\varpi) = \frac{\psi(\varpi_0)}{\varpi_0^{1+s}} \varpi^{1+...
...rac{1}{\varpi^{2n+1}} - \frac{\varpi^{1+s}}{\varpi_0^{2n+2+s}} \right) \right],$     (30)

where the coefficients an and bn are respectively given by:

 \begin{displaymath}%
a_n = I_n \times \left \{
\begin{array}{l}
1-\Delta^{1+s-2n...
...~II},\\ [2mm]
0 \quad {\rm in~region~III},
\end{array}\right.
\end{displaymath} (31)

and

 \begin{displaymath}%
b_n = J_n \times\left \{
\begin{array}{l}
0 \quad {\rm in~r...
...Delta^{2n+s+2}-1 \quad {\rm in~region~III}.
\end{array}\right.
\end{displaymath} (32)

Since $\psi(\Delta)$ and $\psi(1)$ are available (see Sect. 4), we use $\varpi_0 = \Delta$ and $\varpi_0 = 1$ to simplify Eq. (30). Using Eqs. (16) and (19), we thus have:

 \begin{displaymath}%
\frac{\psi(\varpi)}{\psi_{\rm out}} = A \varpi^{1+s} + \sum...
...fty{\left(a_n \varpi^{2n} + \frac{b_n}{\varpi^{2n+1}}\right)},
\end{displaymath} (33)

where:

 \begin{displaymath}%
A =\left \{
\begin{array}{l}
0,\quad {\rm in~region~I},\\ [...
...n~II},\\ [2mm]
0,\quad {\rm in~region~III},
\end{array}\right.
\end{displaymath} (34)

with:

cn = -In -Jn. (35)

We note that A is a function of s.

Figure 4 compares the total potential $\psi(\varpi)$ with the power law contribution (i.e. the term $A \varpi ^{1+s}$) for three typical values of the exponent s in a disc of axis ratio $\Delta =0.1$. We clearly see that, in a finite size disc: i) the gravitational potential is not a power-law function of the radius; ii) a power-law contribution is present inside the disc only; and iii) the power law is not the dominant part of the potential. As expected, spatial self-similarity is broken due to edges. We note that, outside the disc (i.e. in regions I and III), the series coincides with the multi-pole expansions.


  \begin{figure}
\par\includegraphics[width=9cm,clip]{fig4hure09682.eps}
\end{figure} Figure 4: Potential ( plain lines) in a power-law disc with axis ratio $\Delta =0.1$ and $s=\{-2.5,-1.5,-0.5\}$. The contribution of the power law, i.e. the term $A \varpi ^{1+s}$ in Eq. (33), is shown in comparison.
Open with DEXTER

   
6 Potential inside the disc

The determination of the gravitational potential is usually straightforward outside the distribution where different kinds of expansions are efficient in practice (Kellogg 1929). In contrast, it is problematic inside matter where the classical multi-pole approach fails to converge rapidly (e.g. Stone & Norman 1992; Clement 1974).


  \begin{figure}
\par\includegraphics[width=9cm,clip]{fig5hure09682.eps}
\end{figure} Figure 5: Coefficients A ( open circles), $A_\Delta $ ( triangles), and A1 ( stars) versus the power-law exponent s of the surface density. A few remarkable exponents are shown ( open squares): the Mestel (infinite) disc with s=-1 (Mestel 1963), the Hayashi model for the solar nebula with s=-1.5 (Hayashi 1981), and the $\alpha $ and $\beta $-viscosity discs where $s \approx -0.6$ depending on models (Richard & Zahn 1999; Collin-Souffrin & Dumont 1990; Dubrulle 1992; Shakura & Sunyaev 1973).
Open with DEXTER

6.1 Converging series

In region II, we have $b_n = J_n \Delta^{2n+s+2}$. If we set:

\begin{displaymath}%
X = \frac{a_{\rm in}}{R} \equiv \frac{\Delta}{\varpi} \le 1,
\end{displaymath} (36)

then

\begin{displaymath}%
\frac{b_n}{\varpi^{2n+1}} = J_n X^{2n+1} \Delta^{1+s} \le \Delta^{1+s}.
\end{displaymath} (37)

Figure 2 shows the term an versus n for s=-1.5 (in this case, In=Jn). As a consequence, the three series involved in Eq. (33) converge rapidly since i) In and Jn both vary as 1/n2 at large n; ii) all terms are positive for large n; and iii) $\varpi \le 1$ and $X \le 1$. This is interesting for the truncation of series and the construction of reliable approximations (see Sect. 7).

6.2 The coefficient A

The coefficient A is plotted in Fig. 5. It is symmetric with respect to $s=-\frac{3}{2}$ since:

\begin{displaymath}%
\frac{c_n}{\gamma_n} = - \frac{4(4n+1)}{(4n-2s'+1)(4n+2s'+1)},
\end{displaymath} (38)

where $s'=s+\frac{3}{2}$. For certain integer values of |s|, A is strictly zero. This is in particular the case for s=-3 and s=0 (see Appendix A). As Eq. (34) shows, A rises as soon as the exponent s is such that either 2n-s-1 or 2n+s+2 is small (see Sect. 4). Even, if $n = n_\Delta$ (or n=n1), the coefficient A apparently contains a singular term, namely $I_{n_\Delta} \varpi^{1+s}$ (resp. $J_{n_1} \varpi^{1+s}$); however, this singularity exactly cancels with the term $a_{n_\Delta} \varpi^{2 n_\Delta}$ (resp. $b_{n_1} \varpi^{-2n_1 -1}$). In practice, when $s \in {\cal E}_\Delta$, Eq. (33) is no longer valid. Instead, we have:

 \begin{displaymath}%
\frac{\psi(\varpi)}{\psi_{\rm out}} = \left( A_\Delta + \ga...
...n \varpi^{2n}} + \sum_{n=0}^\infty{\frac{b_n}{\varpi^{2n+1}}},
\end{displaymath} (39)

where

 \begin{displaymath}%
A_\Delta =- \sum_{\begin{array}{l}
{\scriptstyle n=0}\\ [-1...
...ne n_\Delta}\end{array}}^\infty{I_n} - \sum_{n=0}^\infty{J_n}.
\end{displaymath} (40)

Similarly, if $s \in {\cal E}_1$, the potential writes

 \begin{displaymath}%
\frac{\psi(\varpi)}{\psi_{\rm out}} = \left( A_1 + \gamma_{...
...tyle n \ne n_1}\end{array}}^\infty{\frac{b_n}{\varpi^{2n+1}}},
\end{displaymath} (41)

where

 \begin{displaymath}%
A_1=- \sum_{n=0}^\infty{I_n} - \sum_{\begin{array}{l}
{\scr...
...}\\ [-1.8mm]
{\scriptstyle n \ne n_1}\end{array}}^\infty{J_n}.
\end{displaymath} (42)

A few values of A, $A_\Delta $, and A1 are listed in Appendix B. We note that $A_\Delta $ (or A1) differs only from A by the term $I_{n_\Delta}$ (resp. Jn1).

   
6.3 Convergence acceleration

Once s is given, the coefficients A, $A_\Delta $, and A1 can be easily determined at the required accuracy. It is also possible to improve the convergence rate of the associated series. This accelerates the computation of the coefficients and makes their dependence with the exponent s more explicit. Convergence acceleration is performed by using the properties of the definite integrals of the complete elliptic integrals of the first and second kinds. The demonstration reported in the Appendix C yields, for $s \notin {\cal E}_\Delta \cup {\cal E}_1$:

 \begin{displaymath}%
A= \frac{1-6C}{\pi} - s(s+2)\sum_{n=0}^\infty{d_n} - (s+1)(s+3)\sum_{n=0}^\infty{e_n},
\end{displaymath} (43)

where

\begin{displaymath}%
\left\{
\begin{array}{l}
d_n = \frac{I_n}{4n^2-1},\\ [2mm]
e_n = \frac{J_n}{4n^2-1}\cdot
\end{array}\right.
\end{displaymath} (44)

and C is the Catalan constant (half the area under the function ${\vec K}$). Numerically, the constant in A is

\begin{displaymath}%
\frac{1-6C}{\pi} \approx -1.4310555380011220.%
\end{displaymath} (45)

Figure 2 shows the coefficient dn for s=-1.5 (in this case, dn=en). It follows that, for large n, dn and en behave like $\sim 1/n^4$ asymptotically, and so A approaches its converged value far more rapidly than by means of Eq. (34).

For $s \in {\cal E}_\Delta$, we then find:

 \begin{displaymath}
%
\begin{array}{l}
A_\Delta = \frac{1-6C}{\pi} + \gamma_{n_\...
... [2.5mm]
\qquad -(s+1)(s+3) \sum_{n=0}^\infty{e_n},
\end{array}\end{displaymath} (46)

instead of Eq. (40), and for $s \in {\cal E}_1$, this is:

 \begin{displaymath}%
\begin{array}{l}
A_1 = \frac{1-6C}{\pi} + \gamma_{n_1} \fra...
...]
{\scriptstyle n \ne n_1}\end{array}}^\infty{e_n},
\end{array}\end{displaymath} (47)

instead of Eq. (42).


  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{fig6hure09682.eps}\end{figure} Figure 6: Relative error in A when computed approximately from Eq. (48).
Open with DEXTER

Depending on the exponent s of interest, a good approximation for A, $A_\Delta $ or A1 can be obtained by considering only the largest terms in the sum, i.e. all terms up to the rank $n \approx \frac{1}{2}(1+s)$ or $n \approx -\frac{1}{2}(s+2)$. For astrophysical discs, s is around -1 meaning that we retain only the first term in Eq. (43). We then find the following approximation:

 \begin{displaymath}%
A \approx \frac{1-6C}{\pi} - \frac{s(s+2)}{s+1} +\frac{(s+1)(s+3)}{s+2},
\end{displaymath} (48)

whose accuracy is better than $3 \%$ for $-2.7 \lesssim s \lesssim -0.3$ as Fig. 6 shows. For $s \in \{-2,-1\}$, we find from Eqs. (46) and (47):

 \begin{displaymath}%
A_\Delta = A_1 \approx \frac{1-6C}{\pi}
\end{displaymath} (49)

which is in good agreement with the converged value given in Table B.


  \begin{figure}
\par\includegraphics[width=9cm,clip]{fig7hure09682.eps}\end{figure} Figure 7: The exact potential for a power law disc with s=-1.5 ( thick line) compared to approximate values ( thin lines) found from Eq. (50) for N=0, 1, 2 and 3 (i.e. see Eqs. (52), (54), (55), etc.). The largest errors are found around edges.
Open with DEXTER

   
7 Approximate formulae inside the disc

Equations (33), (39) and (41) contain three rapidly converging series that can be truncated to derive reliable approximations for the potential. For $s \approx -1$, only a few terms can be considered (see Sect. 6). Although many truncations are possible, we have noticed that the most accurate approximations of $\psi $ for discs[*] are obtained provided the coefficient A (or $A_\Delta $ or A1 depending on s) takes its converged value. Under these circumstances, the N-order approximation for the potential in region II becomes:

 \begin{displaymath}%
\frac{\psi_{\rm app.}^{(N)}(\varpi)}{\psi_{\rm out}} = A \v...
...}^N{\left(a_n \varpi^{2n} + \frac{b_n}{\varpi^{2n+1}}\right)},
\end{displaymath} (50)

which is, in its asymptotic limit:

\begin{displaymath}\lim_{N \rightarrow \infty}\psi_{\rm app.}^{(N)}(\varpi) \equiv \psi(\varpi).
\end{displaymath} (51)

7.1 Zero-order approximation

As argued in Sect. 6.3, a reliable formula for the potential in astrophysical discs (for which $s \approx -1.5 \pm 1$) is obtained by considering only the terms a0 and $b_0/\varpi$, in addition to the power law. At the lowest order, we thus have[*]:

 \begin{displaymath}%
\frac{\psi_{\rm app.}^{(0)}(\varpi)}{\psi_{\rm out}} = \lef...
...elta^{1+s} X}{(2+s)} \quad {\rm otherwise},
\end{array}\right.
\end{displaymath} (52)

where, in this case, $A_1 = A_\Delta \approx -1.386$ (see Appendix B). Figure 7 compares this zero-order approximation with the exact potential for typical disc parameters. It follows that the relative deviation $\vert 1 - \psi_{\rm app.}^{(0)}/\psi\vert$ does not exceed 10%, the deviation being the largest close to the edges.

It is worth noting that the accuracy remains of the same order if coefficient A is determined by Eq. (48). This is convenient when the explicit dependence of $\psi $ on s is required. Under this hypothesis, the potential becomes (see note 4):


 
  $\displaystyle \psi(R) \approx - 2 \pi G \Sigma_{\rm out}a_{\rm out}\times \Bigg\{ \frac{1}{1+s} + \left[ 0.431 - \frac{1}{(1+s)(2+s)} \right]$  
  $\displaystyle \qquad \quad \times \left(\frac{R}{a_{\rm out}}\right)^{1+s} - \f...
...\left(\frac{a_{\rm in}}{a_{\rm out}}\right)^{1+s} \frac{a_{\rm in}}{R} \Bigg\}.$ (53)

Figure 8 shows the accuracy of this formula in the $(\varpi ,s)$-plane. We see that the relative deviation of $\sim$10% observed previously for s=-1.5 holds globally for s roughly in the range[*] [-3,0]. This agrees with the fact that Eq. (48) produces values of A within a few percents for this range of exponents. The deviation can be reduced at the inner and outer edges provided additional terms are included (see below).

7.2 Higher orders

If necessary, more accurate expressions are obtained by accounting gradually for following terms (each acting as a smaller and smaller correction). For N=1, we have:

 \begin{displaymath}%
\frac{\psi_{\rm app.}^{(1)}(\varpi)}{\psi_{\rm out}} = \fra...
... + \frac{\varpi^2}{4(1-s)}
+ \frac{\Delta^{1+s} X^3}{4(4+s)},
\end{displaymath} (54)

and for N=2, this is:

 \begin{displaymath}%
\frac{\psi_{\rm app.}^{(2)}(\varpi)}{\psi_{\rm out}} = \fra...
...rac{9 \varpi^4}{64(3-s)} + \frac{9 \Delta^{1+s} X^5}{64(6+s)},
\end{displaymath} (55)

and so on. We note that Eq. (33) is particularly well suited to numerical computation since $\psi $ can be determined by means of a recurrence procedure. A possible algorithm (not including the treatment of singular cases where $s \in {\cal E}_\Delta \cup {\cal E}_1$) is proposed in Appendix D.


  \begin{figure}
\par\includegraphics[width=9cm,clip]{fig8hure09682.eps}\end{figure} Figure 8: Contour map showing the decimal logarithm of the relative error in the potential in the $(\varpi ,s)$-plane when $\psi $ is approximated by Eq. (53).
Open with DEXTER

   
8 Discs with no inner/outer edges

8.1 Finite disc without inner hole

If the disc has no inner edge but a finite size (i.e. $a_{\rm in}=0$ and $a_{\rm out}\ne \infty$), then the ODE is (see Huré & Hersant 2007):

 \begin{displaymath}%
S(\varpi) = \frac{ 2 \psi_{\rm out}}{\pi \varpi^2}\left \{
...
...{\varpi}\right), \quad {\rm in~region~III}.
\end{array}\right.
\end{displaymath} (56)

It can be verified that the solution is still described by Eqs. (33) and (43), but the coefficients an and bn are:

\begin{displaymath}%
a_n = \left \{
\begin{array}{l}
I_n, \quad {\rm in~region~II}\\ [2mm]
0, \quad {\rm in~region~III}
\end{array}\right.
\end{displaymath} (57)

and

\begin{displaymath}%
b_n = \left \{
\begin{array}{l}
0, \quad {\rm in~region~II}\\ [2mm]
-J_n, \quad {\rm in~region~III}.
\end{array}\right.
\end{displaymath} (58)

8.2 Infinitely extended disc with inner hole

If the disc is infinite but has an inner edge (i.e. $a_{\rm in}>0$ and $a_{\rm out}\rightarrow \infty$), then

\begin{displaymath}%
\frac{{\rm d} \psi}{{\rm d}Y} - (1+s)\frac{\psi}{Y} = S(Y),
\end{displaymath} (59)

where

\begin{displaymath}%
Y = \frac{R}{a_{\rm in}} \equiv \frac{\varpi}{\Delta},
\end{displaymath} (60)

is the new space variable,

 \begin{displaymath}%
S(Y) = - \frac{ 2 \psi_{\rm in}}{\pi Y^2}\left \{
\begin{ar...
...c{1}{Y}\right)\!, \quad {\rm in~region~II},
\end{array}\right.
\end{displaymath} (61)

and $\psi_{\rm in}= 2 \pi G \Sigma_{\rm in}a_{\rm in}$ is a constant (this is not the potential at the inner edge). The analogue of Eq. (33) is:

 \begin{displaymath}%
\frac{\psi(Y)}{\psi_{\rm in}} = A Y^{1+s} + \sum_{n=0}^\infty{\left(a_n Y^{2n} + \frac{b_n}{Y^{2n+1}}\right)},
\end{displaymath} (62)

where A is still given by Eq. (34),

\begin{displaymath}%
a_n =\left \{
\begin{array}{l}
- I_n, \quad {\rm in~region~I}\\ [2mm]
0, \quad {\rm in~region~II}
\end{array}\right.
\end{displaymath} (63)

and

\begin{displaymath}%
b_n =\left \{
\begin{array}{l}
0, \quad {\rm in~region~I}\\ [2mm]
J_n, \quad {\rm in~region~II}.
\end{array}\right.
\end{displaymath} (64)

8.3 Infinitely extended disc

If the disc is infinite, the ODE become homogeneous:

\begin{displaymath}%
\frac{{\rm d} \psi}{{\rm d}R} - (1+s)\frac{\psi}{R} = 0,
\end{displaymath} (65)

and the solution is a power law:

\begin{displaymath}%
\psi(R) = \psi(R_0) \left(\frac{R}{R_0}\right)^{1+s},
\end{displaymath} (66)

where R0 is some reference radius. In this case only, a self-similar surface density can rigorously be associated with a self-similar potential. The presence of edges destroys this property. We note that, for s=-1 (i.e. Mestel's disc), the derivation of the ODE requires a careful treatment. The integral form, i.e. Eq. (5), gives:

\begin{displaymath}%
\psi(R) = - 4 G \Sigma_0 a_0 \left\{ 2C + \frac{\pi}{2} \le...
...{2n} + \ln \frac{a}{R} \right]_{a=a_{\rm out}}^{a=R} \right\}.
\end{displaymath} (67)

We then have

\begin{displaymath}\frac{{\rm d} \psi}{{\rm d} R} = 2 \pi G \Sigma_0 a_0 \frac{1}{R}\cdot
\end{displaymath} (68)

This expression is compatible with Eq. (56a) when $a_{\rm out}\rightarrow \infty$ (in this case, region III no longer exists and we have $S=\psi_{\rm out}/\varpi$).

9 Concluding remarks

In this paper, we have derived an exact expression for the gravitational potential in the plane of flat power-law discs as a solution of the ODE reported in Paper I. This expression is valid over the entire spatial domain and takes into account finite size effects. Inside the disc (the most difficult case to treat in general), it consists of three terms of comparable magnitude: a power law of the cylindrical radius R with index 1+s (where s is the exponent of the surface density) and two series of R and 1/R. In terms of convergence, our expression is by far superior to the multi-pole expansion method. Reliable approximations for the potential can be produced by performing fully controlled truncations. We have shown that the potential can be expressed by means of a simple function of R and s, which is valid to within a few percents in the range of exponents $-3 \lesssim s \lesssim 0$. This formula should be sufficiently accurate for most astrophysical applications. If necessary, more accurate formulae can be developped by including successive terms. These results should help in investigating various phenomena where disc gravity plays a significant role.

An interesting point concerns the case of discs for which the surface density is not a power-law function. As shown in Paper I, it is easy to reproduce numerically the potential when the profile $\Sigma(R)$ is a mixture of power laws. From an analytical point of view however, the construction of a reliable formula for $\psi $ as compact as the one obtained here is not guaranteed at all. For instance, for an expansion of the form:

\begin{displaymath}%
\Sigma(R)=\sum_{s=0}^M{\alpha_s R^s},
\end{displaymath} (69)

where s is an integer, each of the M+1 series $\{a_n,b_n\}$ should be truncated at a rank $N \gtrsim n_\Delta \sim (s+1)/2$ (see Eqs. (33) and (50)), which corresponds to an approximate formula for $\psi $ containing about 2M(M+1) terms. The number of terms to consider can become prohibitively large when several power laws are required.

Acknowledgements
F. Hersant was supported by a CNRS fellowship which is gratefully acknowledged. We thank C. Baruteau. We thank the anonymous referee for valuable comments.

   
Appendix A: Vanishing coefficient A for s $\in $ {-3,0}

For $s\in \{-3,0\}$, the coefficient is:

\begin{displaymath}A = - \sum_0^\infty{\gamma_n \left( \frac{1}{2n+2}+\frac{1}{2n-1} \right)}.
\end{displaymath} (A.1)

We notice that the first sum is the definite integral of u K(u), which is
              $\displaystyle \sum_0^\infty{\frac{\gamma_n}{2n+2}}$ = $\displaystyle \sum_0^\infty{\gamma_n \int_0^1{u^{2n+1}{\rm d}u}},$  
  = $\displaystyle \int_0^1{ {\rm d}u \sum_0^\infty{\gamma_n u^{2n+1}}},$  
  = $\displaystyle \int_0^1{u {\rm d}u \sum_0^\infty{\gamma_n u^{2n}{\rm d}u}},$  
  = $\displaystyle \frac{2}{\pi}\int_0^1{{\vec K}(u) u {\rm d}u},$  
  = $\displaystyle \frac{2}{\pi} \left[ {\vec E}(u) - \left(1 - u^2 \right) {\vec K}(u) \right]_0^1,$  
  = $\displaystyle \frac{2}{\pi}\cdot$ (A.2)

To find the second sum, we compute the definite integral of K(u)/u2 in two ways. First, we have by direct integration (e.g. Gradshteyn & Ryzhik 1965):
         $\displaystyle \frac{2}{\pi}\int_0^1{\frac{{\vec K}(u)}{u^2}{\rm d}u}$ = $\displaystyle \frac{2}{\pi} \left[ - \frac{{\vec E}(u)}{u}\right]_0^1,$  
  = $\displaystyle - \frac{2}{\pi} \left[ {\vec E}(1) - \lim_{u \rightarrow 0} \frac{{\vec E}(u)}{u} \right],$  
  = $\displaystyle - \frac{2}{\pi} + \frac{2}{\pi} \lim_{u \rightarrow 0} \frac{{\vec E}(u)}{u},$  
  = $\displaystyle - \frac{2}{\pi} + \lim_{u \rightarrow 0} \frac{1}{u}\cdot$ (A.3)

Second, by replacing above ${\vec K}$ by its series representation, we also find:
                        $\displaystyle \frac{2}{\pi}\int_0^1{\frac{{\vec K}(u)}{u^2}{\rm d}u}$ = $\displaystyle \int_0^1{{\rm d}u \sum_0^\infty{\gamma_n u^{2n-2}}},$  
  = $\displaystyle \int_0^1{{\rm d}u \left( \frac{\gamma_0}{u^2} + \sum_1^\infty{\gamma_n u^{2n-2}}\right)},$  
  = $\displaystyle - \gamma_0 \left[1 - \lim_{u \rightarrow 0} \frac{1}{u} \right] + \sum_1^\infty{\frac{\gamma_n}{2n-1}},$  
  = $\displaystyle - \gamma_0 + \gamma_0 \lim_{u \rightarrow 0} \frac{1}{u} + \left(...
...ma_0}{-1} + \sum_1^\infty{\frac{\gamma_n}{2n-1}} - \frac{\gamma_0}{-1} \right),$  
  = $\displaystyle \gamma_0 \lim_{u \rightarrow 0} \frac{1}{u} + \sum_0^\infty{\frac{\gamma_n}{2n-1}}\cdot$ (A.4)

Since $\gamma_0=1$, we have

\begin{displaymath}\sum_0^\infty{\frac{\gamma_n}{2n-1}} = - \frac{2}{\pi},
\end{displaymath} (A.5)

and so A=0 for s=-3 and s = 0.

   
Appendix B: Converged values of coefficient A, A1, and A$_\Delta $ for s $\in $ [-4,1]

Table B.1: Values of A (4th column) computed within the computer precision (double precision; about 16-digit) for different power law exponents s or s' (Cols 1 to 3) including those relevant for astrophysical applications (see also Fig. 5). Also given are the relative error on the coefficient (5th column) and the number of iterations (6th column) required.

   
Appendix C: Convergence acceleration for A

The coefficient A is given by a series whose convergence can be accelerated by considering an interesting property of the integral of ${\vec K}$ (e.g. Gradshteyn & Ryzhik 1965), namely:

 \begin{displaymath}\int_0^1{{\vec K}(x){\rm d}x} = \frac{\pi}{2}\sum_0^\infty{\frac{\gamma_n}{2n+1}} = 2C,
\end{displaymath} (C.1)

where C is Catalan's constant. We can then write
                       $\displaystyle \sum_{n=0}^\infty{\frac{\gamma_n}{1+s-2n}}$ = $\displaystyle \sum_{n=0}^\infty{\gamma_n\left(\frac{1}{1+s-2n}+\frac{1}{2n+1}\right)} - \sum_{n=0}^\infty{\frac{\gamma_n}{2n+1}}$  
  = $\displaystyle (s+2) \sum_{n=0}^\infty{\frac{\gamma_n}{(2n+1)(1+s-2n)}} - \frac{4C}{\pi},$ (C.2)

where terms vary like 1/n3 for large n. A second convergence acceleration can be obtained by considering the complete elliptic integral of the second kind:

\begin{displaymath}{\vec E}(x)=\int_0^{\pi/2}{{\rm d} \phi \sqrt{1-x^2 \sin^2 \phi}}, \quad 0 \le x \le 1,
\end{displaymath} (C.3)

whose definite integral over the modulus x is:

 \begin{displaymath}\int_0^1{{\vec E}(x){\rm d}x} = -\frac{\pi}{2}\sum_0^\infty{\frac{\gamma_n}{(2n+1)(2n-1)}} = \frac{1}{2}+C.
\end{displaymath} (C.4)

We then have:


    $\displaystyle \sum_{n=0}^\infty{\frac{\gamma_n}{(2n+1)(1+s-2n)}} = \sum_{n=0}^\infty\gamma_n\left[\frac{1}{(2n+1)(1+s-2n)}\right.$  
    $\displaystyle \qquad \qquad\left. \quad+\frac{1}{(2n+1)(2n-1)}\right] - \sum_{n=0}^\infty{\frac{\gamma_n}{(2n+1)(2n-1)}}$  
    $\displaystyle \quad = s \sum_{n=0}^\infty{\frac{\gamma_n}{(2n+1)(1+s-2n)(2n-1)}} +\frac{2}{\pi}\left(\frac{1}{2}+C \right).$ (C.5)

It follows that:

                          $\displaystyle \sum_{n=0}^\infty{\frac{\gamma_n}{1+s-2n}}$ = $\displaystyle s(s+2) \sum_{n=0}^\infty{\frac{\gamma_n}{(2n+1)(1+s-2n)(2n-1)}}$  
    $\displaystyle + \frac{2+s+2sC}{\pi}\cdot$ (C.6)

The second term in Eq. (34) gives:
                     $\displaystyle \sum_{n=0}^\infty{\frac{\gamma_n}{2n+s+2}}$ = $\displaystyle \sum_{n=0}^\infty\gamma_n\left(\frac{1}{2n+s+2}-\frac{1}{2n+1}\right)$  
    $\displaystyle \quad+ \sum_{n=0}^\infty\frac{\gamma_n}{2n+1}$  
  = $\displaystyle \! -(s+1)\!\! \sum_{n=0}^\infty\frac{\gamma_n}{(2n+1)(2n+s+2)} \!\!+\!\! \frac{4C}{\pi}\cdot$ (C.7)

But, from Eq. (C.4), we have:


    $\displaystyle \sum_{n=0}^\infty{\frac{\gamma_n}{(2n+1)(2n+s+2)}} = \sum_{n=0}^\infty\gamma_n\left[\frac{1}{(2n+1)(2n+s+2)}\right.$  
    $\displaystyle \left.\qquad \qquad-\frac{1}{(2n+1)(2n-1)}\right] + \sum_{n=0}^\infty\frac{\gamma_n}{(2n+1)(2n-1)}$  
    $\displaystyle = \!-\!(s+3)\!\sum_{n=0}^\infty\!{\frac{\gamma_n}{(2n+1)(2n-1)(2n+s+2)}} \!-\! \frac{2}{\pi}\left(\frac{1}{2}+C \right),$ (C.8)

and then:


    $\displaystyle \sum_{n=0}^\infty{\frac{\gamma_n}{2n+s+2}} = (s+1)(s+3)$  
    $\displaystyle \times\sum_{n=0}^\infty{\frac{\gamma_n}{(2n+1)(2n-1)(2n+s+2)}} + \frac{1+s+2(s+3)C}{\pi}\cdot$ (C.9)

Finally, we have:

                                        A = $\displaystyle \sum_{n=0}^\infty{c_n}$  
  = $\displaystyle s(s+2) \sum_{n=0}^\infty{\frac{\gamma_n}{(2n+1)(1+s-2n)(2n-1)}}$  
    $\displaystyle \quad - (s+1)(s+3) \sum_{n=0}^\infty{\frac{\gamma_n}{(2n+1)(2n-1)(2n+s+2)}}$  
    $\displaystyle \quad + \frac{1-6C}{\pi}\cdot$ (C.10)

   
Appendix D: A basic algorithm

A possible algorithm is the following (see the text for cases where $s \in {\cal E}_1 \cup {\cal E}_\Delta$):

1.
disc parameters (see Sect. 2):

- set the inner edge $a_{\rm in}=???$
- set the outer edge $a_{\rm out}=???$
- set the outer surface density $\Sigma_{\rm out}=???$
- set the power law exponent s=???

     $\rightarrow$ compute $\Delta$.
2.
power-law coefficient A:

     $\rightarrow$ compute A (or $A_\Delta $ or A1 depending on s) at the computer accuracy (see Sect. 7); see Eqs. (43), (46) and (47), or use pre-computed converged values (see Appendix B).
3.
radius

- set the radius R = ???

     $\rightarrow$ compute dimensionless variables $\varpi$ and X (depending on the region I, II or III).

4.
initializations

- set $\gamma_0$ to 1; see Eqs. (13),
- set a0 to $\frac{-1}{1+s}$; see Eq. (31),
- set b0 to $\frac{1}{2+s}$; see Eq. (32),
- set $\psi $ to 0.
5.
main loop on n

     $\rightarrow$ compute $\gamma _n$ from $\gamma_{n-1}$,
     $\rightarrow$ compute an from an-1,
     $\rightarrow$ compute bn from bn-1,
     $\rightarrow$ compute $\varpi^{2n}$ and X2n+1,
     $\rightarrow$ update $\psi $ to $\psi+ a_n \varpi^{2n} + b_n \varpi^{-2n-1}$; see Eq. (50).
6.
power-law contribution

     $\rightarrow$ update $\psi $ to $\psi+ A \varpi^{1+s}$; see Eq. (50).
7.
final potential value

     $\rightarrow$ change $\psi $ for $\psi \times 2 \pi G \Sigma_{\rm out}a_{\rm out}$.

The loop ends after N steps; the accuracy of the potential value is then given by the next term (rank N+1).

References

 

Copyright ESO 2008