Volume 507, Number 1, November III 2009
|Page(s)||573 - 579|
|Section||Numerical methods and codes|
|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
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 is twofold: i) to avoid numerical singularities in the integral representation of the potential (i.e., arising when the separation ); and ii) to account for stratification of matter in the direction perpendicular to the disc mid-plane. So far, most studies have considered 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 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 , 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 : 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 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
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 between particles vanishes, one classically adds a positive constant 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 -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 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, is found to be a sharp function of the relative separation ; it is a fraction (of the order of 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 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 . We present a table showing the great diversity (and incoherence) of prescriptions used so far and derive the general relation for . 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
We consider two axially symmetric discs (see Fig. 1):
i) a flat (i.e., zero thickness) disc with inner edge ,
and total surface density
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
is assumed at any radius (Pringle 1981). Moreover, we have
where 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, , , 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 for all a, the two discs therefore become equivalent.
A flat disc and a geometrically thin disc, axially symmetric and symmetric with respect to the mid-plane. The disc edges are and , h is the semi-thickness, 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)
is the modulus and is the complete elliptic integral of the first kind (Legendre form). For the flat disc, we have everywhere. Consequently, the potential at radius R is
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
Clearly, is a function of a, R, and h and contains the effects of vertical stratification. So, if is known preliminarily at all radii and for a given disc structure (edges, thickness, mass density profile), is found from a single integral over the radius (as for ). 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 , where is a certain modulus to be determined (see Sect. 4).
The presence of the function
does not allow us to derive
analytically for any mass density distribution. The expansion of
appropriate for treating the logarithmic singularity is (e.g. Gradshteyn &
where is the complementary modulus, and
By construction, this expansion is efficient when , which corresponds to . It can be shown that the condition implies that
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 is tricky (e.g., Stemwedel et al. 1990), because of the divergence of . 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,
and Eq. (8)
with an error of .
We can estimate
from Eq. (13)
in a few particular cases, for instance when the disc is vertically
in this case). If
does not depend on z, then
and we find that
This is a complicated function of R/a, where the local aspect ratio 2h/a is the unique parameter. If we define the variable
which measures the radial separation from the singularity R=a in units of the disc semi-thickness h, and
which is the quarter of the aspect ratio, we find that
The function is plotted in Fig. 2 for h/a=0.1 (that is ). To order zero around x=0, Eq. (14) reads
We see that reaches a maximum value of at R=a, and is slightly asymmetric with respect to x=0, since we have .
We can probe the sensitivity of
to vertical stratification by considering the following profiles:
where is the density at the disc mid-plane (possibly a function of a) and 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 , we recover the homogeneous case considered above.
where we have defined
We can perform the integration using the expansion of around k'=0, as above. To the lowest-order, we have
This expression can be rewritten in a different form since each sum represent the first q+1 terms of the expansion of the function. In particular, a form appropriate to numerical computations around the singularity is
We see that is a function of R/a, and h/a is the parameter. Using again the variable x and the -parameter, we find
The function is displayed versus x in Fig. 3 for h/a=0.1 and different values of q. Around x=0, Eq. (26) becomes
and we note that
The maximum of still occurs at R=a (see Fig. 3):
Again, we note that the difference is insensitive to q at the actual order of precision.
for . For large values of the q-parameter, . Figure 4 displays versus x for h/a=0.1 and different values of q. It is remarkable that is very weakly sensitive to the vertical profile and that remains very close to (the homogeneous case), especially for .
Variation of 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 .
|Open with DEXTER|
The generalization of the above result is relatively straightforward
for profiles of the form
where cn are real coefficients. In this case, we find that
where is found from Eq. (26). Again, is a function of x, and is the parameter. It is especially convenient to write as
where is the deviation relative to the homogeneous case. We note that is generally small compared to as long as the bulk of the mass is localized in the disc mid-plane. We conclude that is determined mainly by the homogeneous contribution , provided that the matter is gathered around the mid-plane. For , we get from Eqs. (19) and (27)
We see that reaches a maximum value at x=0. This value is
|Open with DEXTER|
As mentioned, the complete elliptic integral
exhibits a logarithmic divergence as its modulus
This is the direct consequence of the point mass singularity (i.e., ). The
determination of accurate potentials from Eqs. (3)
is therefore not straightforward and requires a careful treatment of
improper integrals (e.g., Stemwedel et al.
& Pierens 2005). This technical difficulty is usually
circumvented by changing the relative separation according to
where 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
It is however easy to show that the associated ``softened'' potentials
respectively, can still be written in terms of Eqs. (3)
respectively, provided that the modulus k
is replaced by
and m is replaced by
In disc models and simulations, one never (or rarely) computes , or its fully asymmetric/tri-dimensional version (see however Li et al. 2009). Instead, one computes in the framework of softened gravity, that is . The softening length appearing in Eq. (38) must therefore be prescribed. We note that solving Eq. (38) for leads to
where is the complementary modulus.
Many prescriptions for
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 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
We can see from Eqs. (5)
that the gravitational potential in the mid-plane of a thin
disc is equivalent to the softened potential caused by a flat disc
is the root of the equation
for all R. Only a numerical approach can yield the exact value of , 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
where is given by Eq. (32). In other words, this is
We therefore conclude that the appropriate prescription for , valid over a few disc thicknesses around the singularity R=a, is given by Eq. (39), where is found from Eq. (42).
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|
It follows that
is proportional to h, and depends both on x
It is also sensitive to the vertical stratification through the
versus x for h/a=0.1
in the homogenous case (i.e., for ,
or ). We see
that, in the range of validity,
goes through a minimum for x=0. There, we have
and then the minimum is
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 at the minimum versus . We note however that the main variations of 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.
Since in this case, we also have
The softening length at the minimum x=0 for the vertically homogeneous profile.
|Open with DEXTER|
|Open with DEXTER|
To illustrate the generality and power of the result, we consider a
cosine profile of the form 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:
The corresponding function is then deduced from Eq. (32) and is computed from Eqs. (39) and (42). Results are displayed in Figs. 4 and 7 for and respectively. In particular, at x=0, we have
which results in (this value would be obtained from Eq. (45) for , which is close to the parabolic case).
In this paper, we have reported the first reliable prescription for the softening length 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.
- Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959 [CrossRef] [NASA ADS]
- Anthony, D. M., & Carlberg, R. G. 1988, ApJ, 332, 637 [CrossRef] [NASA ADS]
- Baruteau, C., & Masset, F. 2008, ApJ, 678, 483 [CrossRef] [NASA ADS]
- Caunt, S. E., & Tagger, M. 2001, A&A, 367, 1095 [EDP Sciences] [CrossRef] [NASA ADS]
- Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [CrossRef] [NASA ADS]
- Dehnen, W. 2001, MNRAS, 324, 273 [CrossRef] [NASA ADS]
- Durand, E. 1953, Electrostatique, Vol. I. Les distributions (Masson)
- Edgar, R. G., & Quillen, A. C. 2008, MNRAS, 387, 387 [CrossRef] [NASA ADS]
- Gradshteyn, I. S., & Ryzhik, I. M. 1965, Table of integrals, series and products (New York: Academic Press), 4th edn., ed. Yu.V. Geronimus, & M. Yu. Tseytlin
- Hirose, S., Krolik, J. H., & Stone, J. M. 2006, ApJ, 640, 901 [CrossRef] [NASA ADS]
- Hockney, R. W., & Eastwood, J. W. 1988, Computer simulation using particles, ed. R. W. Hockney, & J. W. Eastwood
- Huré, J.-M. 2005, A&A, 434, 1 [EDP Sciences] [CrossRef] [NASA ADS]
- Huré, J.-M., & Pierens, A. 2005, ApJ, 624, 289 [CrossRef] [NASA ADS]
- Huré, J.-M., & Pierens, A. 2006, in SF2A-2006: Semaine de l'Astrophysique Francaise, ed. D. Barret, F. Casoli, G. Lagache, A. Lecavelier, & L. Pagani, 105
- Laughlin, G., Korchagin, V., & Adams, F. C. 1997, ApJ, 477, 410 [CrossRef] [NASA ADS]
- Laughlin, G., Korchagin, V., & Adams, F. C. 1998, ApJ, 504, 945 [CrossRef] [NASA ADS]
- Li, S., Buoni, M. J., & Li, H. 2009, ApJS, 181, 244 [CrossRef] [NASA ADS]
- Morishima, T., & Saio, H. 1994, MNRAS, 267, 766 [NASA ADS]
- Nelson, A. F. 2006, MNRAS, 373, 1039 [CrossRef] [NASA ADS]
- Papaloizou, J. C. B., & Lin, D. N. C. 1989, ApJ, 344, 645 [CrossRef] [NASA ADS]
- Pringle, J. E. 1981, ARA&A, 19, 137 [CrossRef] [NASA ADS]
- Romeo, A. B. 1994, A&A, 286, 799 [NASA ADS]
- Romeo, A. B. 1997, A&A, 324, 523 [NASA ADS]
- Saio, H., & Yoshii, Y. 1990, ApJ, 363, 40 [CrossRef] [NASA ADS]
- Shu, F. H., Tremaine, S., Adams, F. C., & Ruden, S. P. 1990, ApJ, 358, 495 [CrossRef] [NASA ADS]
- Sommer-Larsen, J., Vedel, H., & Hellsten, U. 1998, MNRAS, 294, 485 [CrossRef] [NASA ADS]
- Stemwedel, S. W., Yuan, C., & Cassen, P. 1990, ApJ, 351, 206 [CrossRef] [NASA ADS]
- Sterzik, M. F., Herold, H., Ruder, H., & Willerding, E. 1995, Planet. Space Sci., 43, 259 [CrossRef] [NASA ADS]
- Tremaine, S. 2001, AJ, 121, 1776 [CrossRef] [NASA ADS]
- A Fortran 90 package called SingLe is available at: http://www.obs.u-bordeaux1.fr/radio/JMHure/intro2single.php
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.