Free Access
 Issue A&A Volume 500, Number 2, June III 2009 617 - 620 Astrophysical processes https://doi.org/10.1051/0004-6361/200911806 08 April 2009

## Self-gravity at the scale of the polar cell

J.-M. Huré1,2 - A. Pierens3 - F. Hersant1,2

1 - Université de Bordeaux, Observatoire Aquitain des Sciences de l'Univers, 2 rue de l'Observatoire, BP 89, 33271 Floirac cedex, France
2 - CNRS/INSU/UMR 5804/LAB, BP 89, 33271 Floirac Cedex, France
3 - LAL-IMCCE/USTL, 1 impasse de l'Observatoire, 59000 Lille, France

Received 6 February 2009 / Accepted 5 March 2009

Abstract
We present the exact calculus of the gravitational potential and acceleration along the symmetry axis of a plane, homogeneous, polar cell as a function of mean radius , radial extension , and opening angle . Accurate approximations are derived in the limit of high numerical resolution at the geometrical mean of the inner and outer radii (a key-position in current FFT-based Poisson solvers). Our results are the full extension of the approximate formula given in the textbook of Binney & Tremaine to all resolutions. We also clarify definitely the question about the existence (or not) of self-forces in polar cells. We find that there is always a self-force at radius  except if the shape factor , asymptotically. Such cells are therefore well suited to build a polar mesh for high resolution simulations of self-gravitating media in two dimensions. A by-product of this study is a newly discovered indefinite integral involving complete elliptic integral of the first kind over modulus.

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

## 1 Introduction

The structure and dynamical evolution of astrophysical discs is governed mainly by gravity. In this context, the numerical computation of the gravitational potential and forces of discs is of fundamental importance and the improvement in accuracy, resolution and computing time remains an interesting challenge (e.g., Li et al. 2008; Londrillo 2004; Matsumoto & Hanawa 2003; Huré 2005; Mueller & Steinmetz 1995; Jusélius & Sundholm 2007). In many (if not all) hydrodynamical simulations of flat, self-gravitating discs (e.g., Baruteau & Masset 2008; Huber & Pfenniger 2001; Zhang et al. 2008), the potential is derived everywhere in the disc by means of Fast Fourier Transforms (e.g., Binney & Tremaine 1987). This is an efficient technique, which can be extended to deduce directly accelerations (Baruteau & Masset 2008). A weak point however is in the determination of the self-gravitating component (i.e., the effect of a cell on itself), which involves an improper integral. An approximation of this self-potential was proposed in Binney & Tremaine (1987), but for a cell in which the surface density varies with the polar radius R as R-3/2. In the same conditions, Baruteau & Masset (2008) argued that the radial acceleration is expected to vanish at the center of each computational cell. This assertion is certainly correct asymptotically for homogeneous cells as their size becomes smaller and smaller. We stress that because Newton's law goes like R-2, the self-gravitating component generally provides a significant contribution to the total potential or force, and must therefore be treated as properly as possible, especially at high resolution. Any error, even relatively small, may introduce artifacts or biases in models.

In this paper, we report the exact expressions for both the potential and radial acceleration due to a homogeneous polar cell inside the cell itself and study the high-resolution limit. We recall in Sect. 2 the integral expression for the potential along the axis of a polar cell which can be decomposed into two parts: a potential created by an homogeneous disc and the other a potential due to horseshoe-like surface. Sections 3 and 4 are devoted to the determinatons of these two contibutions. The results for the polar cell (potential and acceleration) are obtained in Sect. 5. We discuss the case of high resolutions in Sect. 6 and show that there is always a self-acceleration except if the cell has a special shape. These results are important since they contribute to the reliability and accuracy of current simulations of self-gravitating media on a polar mesh.

## 2 The polar mesh

We consider a planar, homogeneous, polar cell with inner polar radius a1, outer radius , lower azimuth , and upper azimuth as shown in Fig. 1. At any position in the plane of this cell, the gravitational potential is given by Newton's generalized formula:

 (1)

where is the surface density. This cell has a single axis of symmetry defined by . If we introduce the quantity:

 (2)

as well as the new variable such that , the potential along its axis of symmetry is given by:

 (3)

where . The integral over can be decomposed into two integrals with 0 as lower bound. We then have:

 (4)

where

 (5)

is the incomplete elliptic integral of the first kind, m is the modulus, is the amplitude, and is the complete elliptic integral of the first kind. The interpretation of this relation is presented in Fig. 2 and follows from the superposition principle: the first part of the integral of Eq. (4) is the potential due to a flat disc, namely:

 (6)

and the second part is that of a flat, horseshoe-like surface:

 (7)

 Figure 1: A polar cell with inner radius a1, outer radius a2, radial extension , and opening angle . Open with DEXTER
 Figure 2: A polar cell is the superposition of a disc and a horseshoe with negative surface density. Open with DEXTER
We can easily derived a as a function of m from Eq. (2), as well as the corresponding derivative . We have, respectively:

 (8)

and

 (9)

where is the complementary modulus, +m' is for (i.e., matter is located right to the position R) and -m' is for (i.e., matter is located to the left). From Eqs. (2) and (9), we deduce that the self-gravitating potential in the cell along its symmetry axis is:

 (10)

We note that this expression is rarely seen in the literature.

## 3 The homogeneous disc

### 3.1 Classical derivation

To calculate , one classically changes the modulus m of by setting where . The new modulus x is either or depending on R. So, the potential due to an homogeneous disc is given by:

 (12)

where region (I) is for , region (II) is for , and region (III) is for (see Fig. 2). Since the indefinite integrals in Eqs. (12) are known, a close-form expression for  can finally be deduced. It is:

 (14)

### 3.2 A new indefinite integral?

From a numerical point of view, it is preferable to use u and v as modulus of the complete elliptic integrals rather than m (see Sect. 6). However, an interesting issue is raised if we perform a change of modulus in Eqs. (14), to restore m as the modulus of and . We find:

 (15)

If we now compare these expressions with the first part of Eq. (10), we conclude that:

 (16)

for any modulus k and complementary modulus k'. This indefinite integral is probably new (Jeffrey 2009).

## 4 The homogeneous horseshoe

To calculate , we have attempted to proceed as for the disc, by applying a change of modulus. We failed, mainly because, as mentioned, the new amplitude  of F resulting from Eq. (11) now becomes a function of the modulus. This introduces a severe mathematical difficulty that we were unable to circumvent. The answer is however obtained by considering Eq. (16) with incomplete elliptic integrals. Using the partial derivative of F and E with respect to their modulus (keeping the amplitude constant), we find that:

 (18)

where

 (19)

and . In fact, can be expressed exactly in terms of basic functions. Actually, after some algebra, we have:

 (20)

It follows from Eqs. (10), (18), (19), and (20) that the potential along the axis of the homogeneous horseshoe is given by:

 (21)

As for Eq. (14), the expression for region (II) can be simplified.

 Figure 3: The potential along the symmetry axis of a homogeneous polar cell with geometrical parameters a1=0.95, and ( and ). The minimum lies at . Also shown are the potential of the associated disc and horseshoe. Open with DEXTER

## 5 Potential and acceleration in the polar cell

The potential inside a polar mesh and along its symmetry axis is found from Eqs. (15b) and (21b). It is:

 (22)

where with . This expression is exact. Figure 3 displays , and their difference  versus R for a typical cell. By deriving this formula with respect to R, we obtain the relation:

 (23)

which is the opposite of the gravitational acceleration gR along the cell axis. If we set and where

 (24)

is the potential due to the cell at the origin R=0, and

 (25)

then Eq. (23) takes the form:

 (26)

This is the generalization of the Ordinary Differential Equation (ODE) derived by Huré & Hersant (2007) to a polar cell with opening angle (which becomes a disc for ). This ODE can be easily solved numerically using standard schemes since  is known both at the center (see Eq. (24)) and at infinity where . The full generalization to polar cells where the surface density  is a power-law of the radius, i.e., , seems straightforward.

## 6 Approximations at high numerical resolutions and the magic'' shape factor

At high resolution, polar cells tends to become segments, squares, or arcs depending on the shape factor  of the cell defined by:

 (27)

where is the center of physical cells. When both (i.e. ) and , we can expand , , , and over m' (e.g., Gradshteyn & Ryzhik 1965; Carlson & Gustafson 1985) to obtain an approximation for and gR as a function of R and the cell geometrical parameters . A radius plays a key role in current FFT-based Poisson solvers (see references in Sect. 1): this is the geometrical mean of the inner and outer radii, namely . This radius is the center of computational cells when a radial logarithmic mapping is applied to the entire physical mesh. For , the two modulus m1 and m2 are equal (and so are the complementary), and we have to first order. After some algebra, we find:

 (28)

and

 (29)

Figure 4 compares exact values of and given by Eqs. (22) and (23) with approximate values obtained from Eqs. (28) and (29) respectively, for three factors . We note that, for , relative errors increase. This is due to the low dynamic'' of the variable m as it approaches unity (and ). In this case, it is advisable to use u and v as the moduli of elliptic integrals instead of m to restore the accuracy in Eq. (22).

 Figure 4: Top panels: exact potential and acceleration in a polar cell at radius for three shape factors ( circles) as well as their approximations ( lines) by Eqs. (28) and (29). Bottom panels: decimal logarithm of the relative error. Also displayed are the results obtained with the magic'' value . Open with DEXTER

The potential always has its minimum inside the cell at a certain radius reagardless of the opening angle . This radius is found by equating the right-hand side of Eq. (23) to zero. It depends on the parameters , , and . We expect for large opening angles (or low azimuthal resolutions) and for small ones (high azimuthal resolutions). When , as the polar cell becomes a rectangle. Since , there is in general a self-force at at high azimuthal resolution. However, the gravitational acceleration exceptionally cancels out if the shape factor has the remarkable value . This is the single root of Eq. (29). For , we have . Thus, if the computational mesh contains cells all with shape factor , then there is no need to compute the gravitational force at , since it is zero in the high resolution limit.

More generally, since the two moduli m1 and m2 appearing in Eq. (23) are functions of the ratio a2/a1, it is always possible to find a shape factor for which the acceleration vanishes at a fixed value of R/a1. It can therefore be advantageous to build the polar mesh accordingly. For instance, let NR and be the numbers of cells in the radial and azimuthal directions, respectively, and be the radius of the disc inner edge. Once the shape factor has been numerically determined such that gR=0 (this is  in the high resolution limit), the radius of the outer edge is found from the relation:

 (30)

With such a mesh, there is no gravitational acceleration at the center of computational cells.

## 7 Concluding remarks

We have reported exact expressions for the self-gravitating potential and acceleration in homogeneous polar cells whatever their shape. The high resolution limit has been analyzed and reliable approximations have been derived. We confirm the validity of Binney & Tremaine (1987)'s prescription for the singular term K(0,0) to first order. We have shown that there is always a self-acceleration at the center of computational cells. However, this acceleration can cancel out by an appropriate choice of cell shape. This work contributes to the treatment of self-gravity in hydrodynamical simulations, and especially for discs.

Acknowledgements
It is a pleasure to thank P. Barge, C. Baruteau and C. Surville for valuable discussions, as well as A. Jeffrey, co-editor of Tables of integrals, series and products'' by Gradshteyn & Ryzhik (1965).

## Footnotes

... changes
The transformation is (e.g. Gradshteyn & Ryzhik 1965):

 (11)

... known
In particular, we have for any (Gradshteyn & Ryzhik 1965):

 (13)

where is the complete elliptic integral of the second kind and is the complementary modulus.
... derivative
In particular, we have (e.g. Gradshteyn & Ryzhik 1965):

 (17)

where is the incomplete elliptic integral of the second kind.
... work
A Fortran 90 package called PolarCELL is available at the following address: http://www.obs.u-bordeaux1.fr/radio/JMHure/intro2polarcell.html

## All Figures

 Figure 1: A polar cell with inner radius a1, outer radius a2, radial extension , and opening angle . Open with DEXTER In the text

 Figure 2: A polar cell is the superposition of a disc and a horseshoe with negative surface density. Open with DEXTER In the text

 Figure 3: The potential along the symmetry axis of a homogeneous polar cell with geometrical parameters a1=0.95, and ( and ). The minimum lies at . Also shown are the potential of the associated disc and horseshoe. Open with DEXTER In the text

 Figure 4: Top panels: exact potential and acceleration in a polar cell at radius for three shape factors ( circles) as well as their approximations ( lines) by Eqs. (28) and (29). Bottom panels: decimal logarithm of the relative error. Also displayed are the results obtained with the magic'' value . Open with DEXTER In the text