EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 507 / No 1 (November III 2009) A&A, 507 1 (2009) 573-579 Full HTML
Free Access
Issue
A&A
Volume 507, Number 1, November III 2009
Page(s) 573 - 579
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/200912348
Published online 03 September 2009

A&A 507, 573-579 (2009)

A local prescription for the softening length in self-gravitating gaseous discs

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

1 - Université de Bordeaux, OASU, 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 17 April 2009 / Accepted 2 July 2009

Abstract
In 2D-simulations of self-gravitating gaseous discs, the potential is often computed in the framework of ``softened gravity'' initially designed for N-body codes. In this special context, the role of the softening length $\lambda $ is twofold: i) to avoid numerical singularities in the integral representation of the potential (i.e., arising when the separation $\vert\vec{r} -\vec{r}'\vert \rightarrow 0$); and ii) to account for stratification of matter in the direction perpendicular to the disc mid-plane. So far, most studies have considered $\lambda $ as a free parameter and various values or formulae have been proposed without much mathematical justification. In this paper, we demonstrate by means of a rigorous calculus that it is possible to define $\lambda $ such that the gravitational potential of a flat disc coincides at order zero with that of a geometrically thin disc of the same surface density. Our prescription for $\lambda $, valid in the local, axisymmetric limit, has the required properties i) and ii). It is mainly an analytical function of the radius and disc thickness, and is sensitive to the vertical stratification. For mass density profiles considered (namely, profiles expandable over even powers of the altitude), we find that $\lambda $: i) is independant of the numerical mesh, ii) is always a fraction of the local thickness H; iii) goes through a minimum at the singularity (i.e., at null separation); and iv) is such that $0.13 \la \lambda/H \la 0.29$ typically (depending on the separation and on density profile). These results should help us to improve the quality of 2D- and 3D-simulations of gaseous discs in several respects (physical realism, accuracy, and computing time).

Key words: accretion, accretion discs - gravitation - methods: numerical

1 Introduction

The computation of gravitational potentials and forces is a critical step in many astrophysical problems where gravity plays an important role. Because Newton's law of gravitation diverges as the relative separation $\vert\vec{r}_i-\vec{r}_j\vert$ between particles vanishes, one classically adds a positive constant $\lambda $ referred as the ``softening length''. This technique of singularity avoidance, initially developped to prevent binary collisions and particle evaporation in N-body simulations, has led to the concept of ``softened gravity'' (Sommer-Larsen et al. 1998; Nelson 2006; Hockney & Eastwood 1988; Anthony & Carlberg 1988). Soon, people realized that the use of a softening length introduces a noticeable bias in models, especially for the stability properties of stellar systems, and there have been several attempts to search for the most appropriate $\lambda $-value (e.g., Sommer-Larsen et al. 1998; Dehnen 2001; Romeo 1994,1997).

Softened gravity has also been employed in the computation of the gravitational potential of gaseous discs (e.g., Sterzik et al. 1995; Adams et al. 1989; Laughlin et al. 1997; Morishima & Saio 1994; Li et al. 2009; Laughlin et al. 1998; Shu et al. 1990; Tremaine 2001; Papaloizou & Lin 1989; Baruteau & Masset 2008; Caunt & Tagger 2001; Saio & Yoshii 1990) with the additional justification that the softening length takes into account the vertical stratification of matter. Various prescriptions, generally in the form of a function of the cylindrical radius and/or of the disc parameters, have been adopted (see references above) without convergence towards a ``universal formula''. As for systems of particles, it is recognized that the softening length can dramatically affect the stability of gaseous discs (again, see references hereabove). It is therefore fundamental to ask i) whether or not softened gravity can definitively help to determine or mimic the Newtonian potential of a gaseous disc with a certain level of accuracy; and ii) eventually which formula - if one exists - is the most appropriate. This is the aim of this paper. By equating the softened potential of a flat disc to that of a geometrically thin disc, we find that we can define a softening length $\lambda $ which has the required properties. We therefore report the first reliable formula for the softening length on the basis of approximate but rigorous calculus. In the axisymmetric limit, $\lambda $ is found to be a sharp function of the relative separation $\vert\vec{r} -\vec{r}'\vert$; it is a fraction (of the order of $\frac{1}{2e}$ at the singularity) of the disc local thickness (see also Baruteau & Masset 2008; Huré & Pierens 2006), and it does not depend on the numerical resolution (see Li et al. 2009). The formula, easy to implement into hydrodynamical codes of self-gravitating 2D- and 3D-discs, will enable to increase the degree of reaslim of simulations, both by preserving the Newtonian character of the potential and force field, and by accounting for the vertical structure.

The paper is organized as follows. In Sect. 2, we recall the expression for the midplane gravitational potential of a geometrically thin disc as well as that of a flat (i.e., zero thickness) disc. Because of the hypothesis of axial symmetry, the integral kernel contains a complete elliptic integral of the first kind that can be appropriately expanded at the singularity $\vert\vec{r} -\vec{r}'\vert=0$ and its neighborhood. In Sect. 3, we estimate the effect of vertical stratification by performing the analytical integration along the z-direction. Various mass density profiles corresponding to finite size discs are considered, namely the homogeneous case and mixtures of even powers of the altitude. In Sect. 4, we show that this calculus naturally leads to a local prescription for the softening length $\lambda $. We present a table showing the great diversity (and incoherence) of prescriptions used so far and derive the general relation for $\lambda $. We discuss the sensitivity of the softening length to the separation and disc aspect ratio. In Sect. 5, we summarize the results and suggest a few possible extensions and generalizations of this work.

2 Mid-plane gravitational potential in both flat discs and geometrically thin discs: local treatment

2.1 Notation and background

We consider two axially symmetric discs (see Fig. 1): i) a flat (i.e., zero thickness) disc with inner edge $a_{\rm in}$, outer edge  $a_{\rm out}$, and total surface density  $\Sigma _{\rm tot}$; and ii) a geometrically thin disc with the same edges and same total surface density, but local thickness H = 2h>0. For the geometrically thin disc, the condition

\begin{displaymath}\left(\frac{H}{a}\right)^2 \ll 1 \end{displaymath} (1)

is assumed at any radius $a \in [a_{\rm in},a_{\rm out}]$ (Pringle 1981). Moreover, we have

\begin{displaymath}\int_{z_-}^{z+}{\rho(z) {\rm d}z} = \Sigma_{\rm tot}, \end{displaymath} (2)

where $\rho$ is the mass density at the altitude z from the mid-plane, z- is the altitude of the bottom of the disc, and z+ =z-+2h is for the top. In general, $\rho$, $\Sigma _{\rm tot}$, z-, and h are expected to depend on the radius a. In the following, we shall also consider symmetry with respect to the mid-plane, which is the case in most disc models (not true for warped discs for instance). This means that z-+z+=0, and z+=h. For $z+ \rightarrow 0$ for all a, the two discs therefore become equivalent.
\begin{figure} \par\includegraphics[width=8.9cm]{12348fig1.eps} \end{figure} Figure 1:

A flat disc and a geometrically thin disc, axially symmetric and symmetric with respect to the mid-plane. The disc edges are $a_{\rm in}$ and $a_{\rm out}$, h is the semi-thickness, $\Sigma _{\rm tot}$ is the total surface density, and a and R are cylindrical radii.

Open with DEXTER

For the thin disc, the mid-plane gravitational potential at radius R is given by the double integral (e.g. Durand 1953)

\begin{displaymath}\Psi^{\rm thin}(R) = -2G\int_{a_{\rm in}}^{a_{\rm out}}{\sqrt... ...{a}{R}} {\rm d}a\; \int_{-h}^{h}{\rho k \vec{K}(k) {\rm d}z}}, \end{displaymath} (3)

where

\begin{displaymath}k = \frac{2\sqrt{aR}}{\sqrt{(a+R)^2+z^2}} \end{displaymath} (4)

is the modulus and $\vec{K}$ is the complete elliptic integral of the first kind (Legendre form). For the flat disc, we have $\rho(z) = \Sigma_{\rm tot}\delta(z)$ everywhere. Consequently, the potential at radius R is

\begin{displaymath}\Psi^{\rm flat}(R)=-2G \int_{a_{\rm in}}^{a_{\rm out}}{\sqrt{\frac{a}{R}}\; \Sigma_{\rm tot}m \vec{K}(m) {\rm d}a}, \end{displaymath} (5)

where

\begin{displaymath}m=\frac{2\sqrt{aR}}{a+R} \end{displaymath} (6)

is the associated modulus (note that k=m for z=0). Both Eqs. (3) and (5) involve a logarithmic singularity when the modulus of the elliptic integral approaches unity. In spite of this divergence, the potential is generally finite at any radius. It is interesting to see that Eq. (3) can also be written as

\begin{displaymath}\Psi^{\rm thin}(R) = -2G\int_{a_{\rm in}}^{a_{\rm out}}{\sqrt{\frac{a}{R}} \Sigma_{\rm tot}\chi {\rm d}a}, \end{displaymath} (7)

where

\begin{displaymath}\Sigma_{\rm tot}\chi = \int_{-h}^h{\rho k \vec{K}(k) {\rm d}z}. \end{displaymath} (8)

Clearly, $\chi $ is a function of a, R, and h and contains the effects of vertical stratification. So, if $\chi $ is known preliminarily at all radii and for a given disc structure (edges, thickness, mass density profile), $\Psi^{\rm thin}$ is found from a single integral over the radius (as for  $\Psi^{\rm flat}$). It is therefore advantageous to have a one-dimensional formula that describes a bi-dimensional distribution of matter (in the axially symmetric case), especially in terms of computing time. By comparing Eqs. (5) and (7), it is tempting to set $\chi = m_{\rm s}\vec{K}(m_{\rm s})$, where $m_{\rm s}$ is a certain modulus to be determined (see Sect. 4).

2.2 Expansion around the singularity

The presence of the function $\vec{K}$ does not allow us to derive $\chi $ analytically for any mass density distribution. The expansion of $\vec{K}(k)$ appropriate for treating the logarithmic singularity is (e.g. Gradshteyn & Ryzhik 1965)

\begin{displaymath}\vec{K}(k)= \sum_{n=0}^{\infty}{\beta_n P_n(k') {k'}^{2n}}, \end{displaymath} (9)

where $k' = \sqrt{1-k^2}$ is the complementary modulus, and

\begin{displaymath}\left\{ \begin{array}{l} P_0= \ln \frac{4}{k'},\\ P_{n+1}(k'... ...}\right]^2, \qquad{\rm for }\qquad n \ge 0. \end{array}\right. \end{displaymath} (10)

By construction, this expansion is efficient when $k' \rightarrow 0$, which corresponds to $k \rightarrow 1$. It can be shown that the condition ${k'}^2 \ll 1$ implies that

\begin{displaymath}\left(\frac{h}{a+R}\right)^2 \ll 1, \end{displaymath} (11)

and

\begin{displaymath}\frac{\vert R-a\vert}{h} \ll \frac{a}{h}\cdot \end{displaymath} (12)

The first inequality is automatically fulfilled within the geometrically thin disc approximation (see Eq. (1)). The second inequality means that the present calculus is valid only at the singularity R=a and for a few disc thicknesses in radius. This is precisely the radial domain where the numerical determination of  $\Psi^{\rm thin}$ is tricky (e.g., Stemwedel et al. 1990), because of the divergence of $\vec{K}(k)$. We note that there is no special constraint on the local shape h(a) of the disc provided that this remains geometrically thin.

To the lowest order, $\vec{K}(k) \approx \ln\frac{4}{k'}$ and $k \approx 1$, and Eq. (8) becomes:

\begin{displaymath}\Sigma_{\rm tot}\chi \approx 2 \int_0^{h}{\rho \ln \frac{4}{k'}{\rm d}z}, \end{displaymath} (13)

with an error of ${\cal O}({k'}^2 \ln k')$.

3 Effect of vertical stratification: an estimate of the $\chi $-function

3.1 Vertically homogeneous discs

We can estimate $\chi $ from Eq. (13) in a few particular cases, for instance when the disc is vertically homogeneous ($\chi $ is denoted $\chi _0$ in this case). If $\rho$ does not depend on z, then $\Sigma_{\rm tot}=2 \rho h \equiv \Sigma_0$ and we find that

\begin{displaymath}\chi_0 = \ln 4 - \ln k'_\pm-\ln f_0, \end{displaymath} (14)

where

\begin{displaymath}\left\{ \begin{array}{l} {k'_\pm}^2 = \frac{(a-R)^2+h^2}{(a+R... ..._- - \frac{1}{\eta_+} {\rm arctan ~}\eta_+. \end{array}\right. \end{displaymath} (15)

This is a complicated function of R/a, where the local aspect ratio 2h/a is the unique parameter. If we define the variable

\begin{displaymath}x=\frac{R-a}{h} \end{displaymath} (16)

which measures the radial separation from the singularity R=a in units of the disc semi-thickness h, and

\begin{displaymath}\epsilon =\frac{h}{2a}, \end{displaymath} (17)

which is the quarter of the aspect ratio, we find that
                             $\displaystyle \chi_0$ = $\displaystyle \ln 4 - \ln \sqrt{\frac{(1+x^2)(1+\epsilon x)^2}{\epsilon^2+(1+\epsilon x)^2}}\nonumber$  
    $\displaystyle -x ~ {\rm arctan ~}\frac{1}{x} + \frac{1+ \epsilon x}{\epsilon} ~ {\rm arctan ~}\frac{\epsilon}{1+ \epsilon x}\cdot$ (18)

The function $\chi(x)$ is plotted in Fig. 2 for h/a=0.1 (that is $\epsilon=0.05$). To order zero around x=0, Eq. (14) reads

\begin{displaymath}\chi_0 \approx 1 + \ln \frac{4}{\epsilon} + \epsilon x - \frac{\pi}{2} \vert x\vert - \ln \sqrt{1+x^2} + x ~ {\rm arctan ~}x. \end{displaymath} (19)

We see that $\chi _0$ reaches a maximum value of $1 + \ln \frac{4}{\epsilon} \gg 1$ at R=a, and is slightly asymmetric with respect to x=0, since we have $\chi_0(x=+1) \approx \chi_0(x=-1)+ 2\epsilon$.

\begin{figure} \par\includegraphics[width=8.9cm, bb=44 51 706 529, clip=]{12348fig2.eps} \end{figure} Figure 2:

Variation of $\chi _0$ with x according to Eq. (14).

Open with DEXTER

3.2 Effect of vertical stratification

We can probe the sensitivity of $\chi $ to vertical stratification by considering the following profiles:

\begin{displaymath}\rho(z)=\left\{ \begin{array}{ll} \rho_0 \left[1-\left(\frac... ...\vert z\vert \le h,\\ 0 & {\rm otherwise}, \end{array}\right. \end{displaymath} (20)

where $\rho_0$ is the density at the disc mid-plane (possibly a function of a) and $q \ge 1$ is an integer. These correspond to finite size discs. When q=1, the profile is close to a Gaussian distribution, a case found for instance in vertically isothermal disc at hydrostatic equilibrium (e.g., Chiang & Goldreich 1997; Hirose et al. 2006; Edgar & Quillen 2008). As q increases, the vertical profile becomes flatter and flatter. For $q \rightarrow \infty $, we recover the homogeneous case considered above.

We can calculate $\chi $ again from Eq. (13), but using Eq. (20). From Eq. (2), we get $\Sigma_{\rm tot}= \frac{2q}{2q+1} \Sigma_0$, and then

\begin{displaymath}\chi % = \chi_0 + \frac{1}{2q} \left( \chi_0 - \chi_q \right), \end{displaymath} (21)

where we have defined

\begin{displaymath}\chi_q = (2q+1) \int_0^{1}{\left(\frac{z}{h}\right)^{2q} k \vec{K}(k) {\rm d} \frac{z}{h}}\cdot \end{displaymath} (22)

We can perform the integration using the expansion of $\vec{K}(k)$ around k'=0, as above. To the lowest-order, we have

\begin{displaymath}\chi_q = \ln 4 - \ln k'_\pm - \ln f_q, \end{displaymath} (23)

where
                      $\displaystyle \ln f_q$ = $\displaystyle (-1)^q \left[ \frac{{\rm arctan ~}\eta_-}{\eta_-^{2q+1}} - \frac{{\rm arctan ~}\eta_+}{\eta_+^{2q+1}} \right.$  
    $\displaystyle \left. -\sum_{k=0}^q{(-1)^k \frac{\eta_-^{2(k-q)}-\eta_+^{2(k-q)}}{2k+1}} \right]\cdot$ (24)

This expression can be rewritten in a different form since each sum represent the first q+1 terms of the expansion of the ${\rm arctan ~}$ function. In particular, a form appropriate to numerical computations around the singularity is
                                $\displaystyle \ln f_q$ = $\displaystyle (-1)^q \left[ \eta_-^{-2q} \left( \frac{\pi}{2\vert\eta_-\vert} - \frac{1}{\eta_-} {\rm arctan ~}\frac{1}{\eta_-} \right) \right.$  
    $\displaystyle \left. -\sum_{k=0}^q{\frac{(-1)^k}{2k+1} \eta_-^{2(k-q)}} \right] + \eta_+^2\sum_{k=0}^\infty{\frac{(-1)^k \eta_+^{2k}}{2(k+q)+3} }\cdot$ (25)

We see that $\chi _q$ is a function of R/a, and h/a is the parameter. Using again the variable x and the $\epsilon$-parameter, we find
                                    $\displaystyle \chi_q$ = $\displaystyle \ln 4 - \ln \sqrt{\frac{(1+x^2)(1+\epsilon x)^2}{\epsilon^2+(1+\epsilon x)^2}}$  
    $\displaystyle - (-1)^q \left[ x^{2q} \left( \frac{\pi}{2}\vert x\vert - x ~ {\r... ...ctan ~}x \right) -\sum_{k=0}^q{\frac{(-1)^k}{2k+1} x^{2(q-k)}} \right]\nonumber$  
    $\displaystyle - \left(\frac{\epsilon}{1+\epsilon x}\right)^2\sum_{k=0}^\infty{\frac{(-1)^k }{2(k+q)+3} \left(\frac{\epsilon}{1+\epsilon x}\right)^{2k}}\cdot$ (26)

The function $\chi _q$ is displayed versus x in Fig. 3 for h/a=0.1 and different values of q. Around x=0, Eq. (26) becomes
$\displaystyle \chi_q \approx \frac{1}{2q+1} + \ln \frac{4}{\epsilon} + \epsilon... ... - (-1)^q x^{2q} \left(\frac{\pi}{2}\vert x\vert - x ~ {\rm arctan ~}x \right),$     (27)

and we note that

\begin{displaymath}\chi_0 - \chi_q = \frac{2q}{2q+1} + \left[(-1)^q x^{2q} -1 \r... ... \left( \frac{\pi}{2}\vert x\vert - x~{\rm arctan ~}x \right). \end{displaymath} (28)

The maximum of $\chi _q$ still occurs at R=a (see Fig. 3):

\begin{displaymath}\chi_{q, \rm max} \approx \chi_{0, \rm max} - \frac{2q}{2q+1}\cdot \end{displaymath} (29)

Again, we note that the difference $\chi_q(R=a+h) - \chi_q(R=a-h) \approx 2 \epsilon$ is insensitive to q at the actual order of precision.

Finally, for the vertical profile defined by Eqs. (20), $\chi $ follows from Eqs. (18), (21), and (26). We see that $\chi $ reaches a maximum at R=a, and

\begin{displaymath}\chi_{\rm max} \approx 1 + \ln \frac{4}{\epsilon} + \frac{1}{2q+1}, \end{displaymath} (30)

for $q \ge 1$. For large values of the q-parameter, $\chi \rightarrow \chi_0$. Figure 4 displays $\chi $ versus x for h/a=0.1 and different values of q. It is remarkable that $\chi $ is very weakly sensitive to the vertical profile and that $\chi $ remains very close to $\chi _0$ (the homogeneous case), especially for $\vert x\vert \ga 1$.
\begin{figure} \par\includegraphics[width=8.9cm, bb=44 51 706 529, clip=]{12348fig3.eps} \end{figure} Figure 3:

Variation of $\chi _q$ with x according to Eq. (23) for h/a= 0.1 and different values of q, namely q=0 (the homogeneous case), q=1 (the parabolic case), q=2, 3, 5 and $q \rightarrow \infty $.

Open with DEXTER

3.3 Full generalization

The generalization of the above result is relatively straightforward for profiles of the form

\begin{displaymath}\rho(z)= \rho_0 \sum_{n=0}^\infty{c_n \left(\frac{z}{h}\right)^{2n}}, \end{displaymath} (31)

where cn are real coefficients. In this case, we find that

\begin{displaymath}\chi = \frac{1}{\Sigma_{\rm tot}}\sum_{n=0}^\infty{\frac{c_n ... ...\Sigma_{\rm tot}= \Sigma_0\sum_{n=0}^\infty{\frac{c_n}{2n+1}}, \end{displaymath} (32)

where $\chi_n$ is found from Eq. (26). Again, $\chi $ is a function of x, and $\epsilon$ is the parameter. It is especially convenient to write $\chi $ as

\begin{displaymath}\chi = \chi_0 + \delta \chi_0, \end{displaymath} (33)

where $\delta \chi_0$ is the deviation relative to the homogeneous case. We note that $\delta \chi_0$ is generally small compared to $\chi _0$ as long as the bulk of the mass is localized in the disc mid-plane. We conclude that $\chi $ is determined mainly by the homogeneous contribution $\chi _0$, provided that the matter is gathered around the mid-plane. For $x \rightarrow 0$, we get from Eqs. (19) and (27)
$\displaystyle \delta \chi_0 = - \frac{\Sigma_0}{\Sigma_{\rm tot}} \sum_{n=0}^\i... ... \right] \left( \frac{\pi}{2}\vert x\vert - x~{\rm arctan ~}x \right) \right\}.$     (34)

We see that $\chi $ reaches a maximum value at x=0. This value is
$\displaystyle \chi_{\rm max} = \chi_{0, \rm max} - \displaystyle\frac{\Sigma_0}{\Sigma_{\rm tot}} \sum_{n=1}^\infty{\frac{2n c_n}{(2n+1)^2}}\cdot$     (35)

\begin{figure} \par\includegraphics[width=8.9cm, bb=14 51 706 529, clip=]{12348fig4.eps}\end{figure} Figure 4:

$\chi $ versus x for the mass density profile defined by Eq. (20). The cosine profile discussed in Sect 4.2 is also shown.

Open with DEXTER

4 Application to softened gravity: a prescription for the softening length

4.1 The avoidance of point mass singularities

As mentioned, the complete elliptic integral $\vec{K}$ exhibits a logarithmic divergence as its modulus  $k \rightarrow 1$. This is the direct consequence of the point mass singularity (i.e., $\vert\vec{r} -\vec{r}'\vert \rightarrow 0$). The determination of accurate potentials from Eqs. (3) and (5) is therefore not straightforward and requires a careful treatment of improper integrals (e.g., Stemwedel et al. 1990; Huré 2005; Huré & Pierens 2005). This technical difficulty is usually circumvented by changing the relative separation according to

\begin{displaymath}\vert\vec{r}-\vec{r'}\vert \quad \leftarrow \quad \sqrt{\vert\vec{r}-\vec{r'}\vert^2 + \lambda^2}, \end{displaymath} (36)

where $\lambda \ne 0$is a constant known as the ``softening length''. The main drawback is that softened gravity modifies Newton's law for gravitation both on short and long ranges. It lowers the magnitude of forces, enhances stability, and introduces a bias in models that is not easy to measure and interpret (for stellar and gas discs see, e.g., Papaloizou & Lin 1989; Sommer-Larsen et al. 1998; Dehnen 2001; Romeo 1994; Saio & Yoshii 1990).

In the case of gaseous discs of interest here, the modification of the relative distances according to Eq. (36) changes the expressions for  $\Psi^{\rm flat}$ and  $\Psi^{\rm thin}$. It is however easy to show that the associated ``softened'' potentials denoted  $\Psi^{\rm flat}_{\rm s}$ and $\Psi^{\rm thin}_{\rm s}$, respectively, can still be written in terms of Eqs. (3) and (5), respectively, provided that the modulus k is replaced by

\begin{displaymath}k_{\rm s}= \frac{2\sqrt{aR}}{\sqrt{(a+R)^2+z^2 + \lambda^2}}, \end{displaymath} (37)

and m is replaced by

\begin{displaymath}m_{\rm s}= \frac{2\sqrt{aR}}{\sqrt{(a+R)^2 + \lambda^2}}\cdot \end{displaymath} (38)

In disc models and simulations, one never (or rarely) computes $\Psi^{\rm thin}$, or its fully asymmetric/tri-dimensional version (see however Li et al. 2009). Instead, one computes $\Psi^{\rm flat}$ in the framework of softened gravity, that is $\Psi^{\rm flat}_{\rm s}$. The softening length $\lambda $ appearing in Eq. (38) must therefore be prescribed. We note that solving Eq. (38) for $\lambda $ leads to

\begin{displaymath}\frac{\lambda}{h} % = \sqrt{\frac{{m_{\rm s}'}^2}{1- {m_{\rm s}'}^2} \frac{4aR}{h^2} - \frac{(a-R)^2}{h^2}}, \end{displaymath} (39)

where $m_{\rm s}' = \sqrt{1-m_{\rm s}^2}$ is the complementary modulus.

Table 1:  Various prescriptions for the softening length adopted in some simulations of self-gravitating gaseous discs (see also Fig. 1).

4.2 A prescription for the softening length $\lambda $

Many prescriptions for $\lambda $ have been proposed. In general, this is not a constant but a certain function of the radius and/or disc parameters. Table 1 gathers a few formulae for $\lambda $used by different authors over twenty years. Although not exhaustive, this list clearly shows that there is no trend in magnitude and variation in space (and possible dependency on the disc parameters). The results obtained in Sect. 3 can help substantially to define the appropriate prescription for $\lambda $. We can see from Eqs. (5) and (7) that the gravitational potential in the mid-plane of a thin disc is equivalent to the softened potential caused by a flat disc provided that $\lambda $ is the root of the equation

\begin{displaymath}\Sigma_{\rm tot}m_{\rm s}\vec{K}(m_{\rm s}) -\int_{-h}^{+h}{\rho k \vec{K}(k) {\rm d}z} =0, \end{displaymath} (40)

for all R. Only a numerical approach can yield the exact value of $\lambda $, if it exists. However, a good approximation to this root can be obtained by considering the expansion of the complete elliptic integral of the first kind over its complementary modulus, as considered in Sect. 3. To the lowest order, Eq. (40) becomes

\begin{displaymath}\ln \frac{4}{m_{\rm s}'} - \chi =0, \end{displaymath} (41)

where $\chi $ is given by Eq. (32). In other words, this is

\begin{displaymath}m_{\rm s}'=4 {\rm e}^{-\chi}. \end{displaymath} (42)

We therefore conclude that the appropriate prescription for $\lambda $, valid over a few disc thicknesses around the singularity R=a, is given by Eq. (39), where  $m_{\rm s}'$ is found from Eq. (42).
\begin{figure} \par\includegraphics[width=8.9cm, bb=32 51 706 522, clip=]{12348fig5.eps} \end{figure} Figure 5:

Softening length normalized to the local semi-thickness h versus x as computed from Eq. (39) for h/a=0.1 in the homogeneous case.

Open with DEXTER

4.3 Results for various stratifications

It follows that $\lambda $ is proportional to h, and depends both on x and $\epsilon$. It is also sensitive to the vertical stratification through the function $\chi $. Figure 5 displays $\lambda/h$ versus x for h/a=0.1 in the homogenous case (i.e., for $\chi \equiv \chi_0$, or $\delta \chi_0=0$). We see that, in the range of validity, $\lambda $ goes through a minimum for x=0. There, we have $k'_\pm \approx h/2a$ and $f_0 \approx 1/e$ and then the minimum is

\begin{displaymath}\frac{\lambda}{h} \approx \frac{1}{e} \times \frac{1}{\sqrt{1... ...ft[ 1 + \left(\frac{\epsilon}{\sqrt{2}e}\right)^2 \right]\cdot \end{displaymath} (43)

This value is almost insensitive to the disc aspect ratio, provided that the disc is geometrically thin. This is shown in Fig. 6, where we plot $\lambda/h$ at the minimum versus $\epsilon$. We note however that the main variations of $\lambda $ occur over a radial range of the order of 2h (i.e., the total disc thickness). The thinner the disc, the sharper the variation around the singularity.

Figure 7 displays $\lambda/h$ for the mass density profile given by Eq. (20) and for $q=\{1,2,3,5, \infty \}$. At x=0, we have the general formula

\begin{displaymath}\frac{\lambda}{h} \approx \frac{{\rm e}^{-(1 + \delta \chi_0)... ...ft[\frac{h}{2 a} {\rm e}^{- (1+\delta \chi_0)}\right]^2}}\cdot \end{displaymath} (44)

Since $\delta \chi_0=\frac{1}{2q+1}$ in this case, we also have

\begin{displaymath}\frac{\lambda}{h} \approx {\rm e}^{-\frac{2q+2}{2q+1}} = \frac{1}{e} \times {\rm e}^{-\frac{1}{2q+1}}. \end{displaymath} (45)

\begin{figure} \par\includegraphics[width=8.9cm, bb=-2 44 706 532, clip=]{12348fig6.eps} \end{figure} Figure 6:

The softening length at the minimum x=0 for the vertically homogeneous profile.

Open with DEXTER

\begin{figure} \par\includegraphics[width=8.9cm, bb=0 51 706 522, clip=]{12348fig7.eps} \end{figure} Figure 7:

Softening length normalized to the local semi-thickness h versus x as computed from Eq. (39) for h/a=0.1. The density profile is according to Eq. (20) with $q=\{1,2,3,5, \infty \}$. The cosine profile is also shown.

Open with DEXTER

4.4 An example of vertical mass density profile expandable in infinite series

To illustrate the generality and power of the result, we consider a cosine profile of the form $\rho(z)=\cos \frac{\pi z}{2 h}$which is a typical example where matter distribution is expandable in infinite series of the altitude. Actually, this profile can also be expanded by means of Eq. (31) with the following coefficients:

\begin{displaymath}c_n=\frac{(-1)^n}{(2n)!} \left(\frac{\pi}{2}\right)^{2n}\cdot \end{displaymath} (46)

The corresponding function $\chi $ is then deduced from Eq. (32) and $\lambda $ is computed from Eqs. (39) and (42). Results are displayed in Figs. 4 and 7 for $\chi $ and $\lambda/h$ respectively. In particular, at x=0, we have

\begin{displaymath}\delta \chi_0= - \pi \sum_1^\infty{\frac{(-1)^q q}{(2q+1)^2 (2q)!} \left(\frac{\pi}{2}\right)^{2q}} = 0.371... \end{displaymath} (47)

which results in $\frac{\lambda}{h} \approx {\rm e}^{-1} \times 0.690 \approx 0.2539$ (this value would be obtained from Eq. (45) for $q \approx 0.849$, which is close to the parabolic case).

5 Concluding remarks

In this paper, we have reported the first reliable prescription for the softening length $\lambda $to be used in the numerical determination of the gravitational potential of geometrically thin discs within the framework of softened gravity. This expression has been found by rigorously comparing the Newtonian potential of a geometrically thin disc (of finite thickness) with the softened potential of a flat disc. This is a function of the radius and disc local aspect ratio, and also depends on vertical stratification. It is accurate at the singularity and around (typically a few disc thicknesses in radius). Although this formula is valid only locally, and obtained in the axisymmetric limit, it should help to improve the quality (realism, accuracy, and computing time) of 2D- and 3D-simulations[*]. If necessary, the accuracy of the prescription can be improved by considering higher order terms in the expanded complete elliptic integrals of the first kind.

Finally, it would then be interesting to generalize our results i) to mass density profiles that extend to infinity, such as the Gaussian profile (which corresponds to vertically isothermal discs); and ii) to the entire disc since this prescription is expected to be inaccurate far from the singularity.

Acknowledgements
It is a pleasure to thank C. Baruteau, D. Bernard, A. Collioud, F. Hersant and A. Romeo for valuable comments.

References

Footnotes

...3D-simulations[*]
A Fortran 90 package called SingLe is available at: http://www.obs.u-bordeaux1.fr/radio/JMHure/intro2single.php

All Tables

Table 1:   Various prescriptions for the softening length adopted in some simulations of self-gravitating gaseous discs (see also Fig. 1).

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.