The Newtonian potential of thin disks
^{1}
Université de Bordeaux, OASU, 2 rue de l’Observatoire, BP 89, 33271 Floirac Cedex, France
^{2}
CNRS, UMR 5804, LAB, 2 rue de l’Observatoire, BP 89, 33271 Floirac Cedex, France
email: jeanmarc.hure@obs.ubordeaux1.fr
Received: 1 October 2010
Accepted: 30 March 2011
The onedimensional, ordinary differential equation (ODE) that satisfies the midplane gravitational potential of truncated, flat powerlaw disks is extended to the whole physical space. It is shown that thickness effects (i.e. nonflatness) can be easily accounted for by implementing an appropriate “softening length” λ. The solution of this “softened ODE” has the following properties: i) it is regular at the edges (finite radial accelerations); ii) it possesses the correct longrange properties; iii) it matches the Newtonian potential of a geometrically thin disk very well; and iv) it tends continuously to the flat disk solution in the limit λ → 0. As illustrated by many examples, the ODE, subject to exact Dirichlet conditions, can be solved numerically with efficiency for any given colatitude at secondorder from center to infinity using radial mapping. This approach is therefore particularly wellsuited to generating grids of gravitational forces in order to study particles moving under the field of a gravitating disk as found in various contexts (active nuclei, stellar systems, young stellar objects). Extension to nonpowerlaw surface density profiles is straightforward through superposition. Grids can be produced upon request.
Key words: accretion, accretion disks / gravitation / methods: analytical / methods: numerical
© ESO, 2011
1. Introduction
Flattened astrophysical objects like disks are known to genrally produce gravitational fields weaker than spherical bodies of comparable mass. However, the gravity of lowmass disks evolving on long time scales may play a significant role in their own dynamics and environment (Goldreich & Tremaine 1978; Šubr & Karas 2005). Gravitational forces from disks are not easily accessible by numerical computation, and this domain still represents an interesting challenge. For various reasons (misknowledge of boundary conditions, sensitivity and inaccuracy of solutions, relevant physical scales, kernel singularities, computational cost, etc.), neither the Poisson equation nor Newton’s integral law offers a simple and straightforward tool, and each must be handled with some caution. Truncated expansions of solutions generally suffer from inaccuracy and instability (Clement 1974; Hachisu 1986). Softened gravity for continuous systems may help in some circumstances, but the influence of the softening length – a freeparameter, classically – is spurious, and it fundamentally destroys the Newtonian character of the gravitational interaction (Hockney & Eastwood 1988; Adams et al. 1989). Each disk configuration (symmetry, edges, mass profile, shape, etc.) must therefore be investigated individually for a given application.
Geometrically thin disks probably constitute the main class of astrophysical disks. These exhibit various shapes, density profiles and sizes. For those orbiting a massive central object (star or black hole), a selfsimilar behavior may develop secularly, leading to a mass density profile varying close to a power law of the radius. Such a profile is widely supported by theory and it is a typical initial ingredient of numerical simulations (e.g. Pringle 1981; Edgar 2007). Even in quasiKeplerian rotation, thin disks can be influenced by their own gravity. Huré & Hersant (2007, hereafter Paper I) have shown that the midplane gravitational potential of flat, powerlaw disks obey an ordinary differential equation (ODE) accounting for edges which are usually ignored (BisnovatyiKogan 1975; Goodman & Evans 1999). Analytical solutions in the form of very rapidly converging series have been reported in Huré et al. (2008). At the same time, Huré & Pierens (2009) have shown that the model of “softened gravity” offers a good framework for determining the Newtonian potential of thin disks (whatever the mass density profile), provided the “softening length” takes a very specific form, locally. In the present paper, we show that the ODE for the gravitational potential described by Huré & Hersant (2007) can be extended to the whole physical space, and we combine this result with an appropriate softening length to describe the potential of disks of non vanishing thickness.
The paper is organized as follows. We recall the basic configuration and useful formulae of potential in flat, powerlaw disks in Sect. 2. The derivation of the generalized ODE, nondimensionalization and asymptotic behavior of solutions, are found in Sect. 3. The numerical solutions are given in Sect. 4. Thickness effects, including the introduction of the softening length, are discussed in Sect 5. The last section is devoted to a conclusion.
2. Theoretical grounds and notation
Following Huré & Hersant (2007, hereafter Paper I), we consider a flat axisymmetrical disk with inner edge a_{in} ≥ 0, outer edge a_{out} > a_{in} (see Fig. 1), and a powerlaw surface density of the form (1)where a is the cylindrical radius, a_{0} some reference radius, and Σ_{0} ≡ Σ(a_{0}) the corresponding surface density. This profile can serve as a basis for defining more complex mass distributions, by mixing power laws with different indices (positive and negative). For such a disk, the gravitational potential in space is given exactly by the expression (Durand 1953) (2)where (3)is the complete elliptic integral of the first kind, and (4)is the modulus (0 ≤ k ≤ 1), R and Z are the cylindrical coordinates, and r is the spherical radius (i.e. r^{2} = R^{2} + Z^{2}). Known properties about this configuration are the followings. The integral in Eq. (2) has a diverging kernel as soon as the modulus k reaches unity. This occurs everywhere inside the disk. Standard quadrature schemes fail to give accurate potential values unless a specific treatment is considered (Huré & Pierens 2005). The potential is not a power law of the radius, because of edges. A closed form for Eq. (2) exists only in the case of infinitely extended disks (BisnovatyiKogan 1975; Goodman & Evans 1999) and for finite size disks with constant surface density (Durand 1953; Lass & Blitzer 1983). The potential is finite everywhere, except when a_{in} = 0 and s + 1 < 0, but has an infinite radial gradient in the midplane at edge crossing, as a consequence of flatness (Durand 1953; Mestel 1963). Finally, the disk mass is finite except when a_{in} = 0 and s + 2 < 0, or when a_{out} → ∞ and s + 2 > 0.
It has been established in Paper I that the potential ψ given by Eq. (2) obeys, in the disk midplane (i.e. for Z = 0), a differential equation of the form (5)where Λ is fully analytical. Solutions in the form of rapidly converging series were reported in Huré et al. (2008).
3. Generalized ODE
The starting point of the generalization of Eq. (5) is the equation of the line containing the origin and making an angle θ ∈ [0,π] (the colatitude) with the zaxis (see Fig. 1). This equation is simply (6)By setting u = a/R in Eq. (4), and using Eq. (6), we get the relation (7)We then see that, once θ is specified, u is a function of the modulus k only. This property was the condition for the existence of Eq. (5). It still holds here. We can therefore proceed as in Paper I: the partial derivative of the potential with respect to the radius can be expressed as a function of the potential itself. Here, we consider the derivative with respect to the spherical radius r instead of the cylindrical radius R (although this choice is not really important). The full demonstration is given in Appendix A. Strictly speaking, we obtain a partial differential equation (PDE) that becomes, for a given colatitude θ, an ordinary differential equation (ODE). It writes as (8)where Λ is defined by (9)with u_{out} = a_{out}/R, and u_{in} = a_{in}/R. This generalizes to the whole space the differential equation reported in Paper I, which is valid only in the disk plane (i.e. the case ). The second member Λ is an analytical function of a single spatial coordinate r. It now depends on the colatitude θ, and on the four parameters s, a_{0}, a_{in}, and a_{out}. The main difference comes from the Zdependent modulus k.
Fig. 1 Configuration for the flat, finite size disk with radius at the edges a_{in} and a_{out}; r and θ are spherical coordinates, R and Z are cylindrical coordinates. 

Open with DEXTER 
We can make the ODE scalefree by setting^{1}, and where ψ_{0} is a constant. Then, we get the nondimensional ODE: (10)where (11)Potential values are known at the two edges of the disk in the form of rapidly converging series (Huré et al. 2008), and can therefore be used for adimensioning. A more convenient value for ψ_{0} is probably the potential at the origin which is exact. From Eq. (2), we have (12)where we have set (definition different than in Paper I) (13)and Δ = a_{in}/a_{out} is the axis ratio of the disk. Thus, with ψ_{0} ≡ ψ_{c}, the inhomogeneous term in the nondimensional ODE is (14)When a_{in} = 0 (no inner edge; the disk extends down to the center), the last term vanishes. When a_{out} → ∞, this is the first term. The asymptotic properties of the solution are discussed in Appendix B.
4. Numerical solution in space with mapping
There are many methods for solving a firstorder ODE numerically. This needs one boundary condition for an initial value problem (IVP), but two conditions for a two boundary value problem (TBVP). The potential is known at four places: at the origin, at the two edges, and at infinity. The origin is particularly well suited to finding the potential in space in a given direction r, as the corresponding boundary condition is the same regardless of the colatitude. Values at the disk edges are useful mainly for determining the potential inside the disk or in its neighborhood. For some applications, potential values are required at large distances from the disk (e.g. Šubr et al. 2004; Šubr & Karas 2005). In this case, it can be interesting to use the trivial boundary condition at infinity (ψ → 0 as r → ∞). This is easily performed by mapping the whole space. For instance, the transformation^{2}(15)where A is a positive constant, maps the range into the compact domain . The “mapped ODE” reads (16)where we have set (17)and (18)It is then easy to discretize Eq. (16) on a grid, regular in ζ, and to solve the associated linear system by standard techniques. Actually, a centered secondorder space discretization based on N + 1 mesh points yields the system: with ζ_{0} = 0, , and ζ_{n} = nδζ. For s = −1, it can be advantageous to rewrite the ODE in a slightly different form, for instance by using the auxiliary function , which satisfies the following ODE: (19)with the boundary conditions and .
To illustrate this numerical part, we show in Fig. 2 the solution obtained at a few colatitudes with N = 1000 and (see note 2) and for the following parameters: Δ = 0.1, a_{out}/a_{0} = 1, and s = −1.5 (due to adimensioning, it is not necessary to specify Σ_{0}). Figures 3 to 5 correspond to power indices s = {−1, −0.5,0 } .
Fig. 2 Potential versus ζ due to a flat disk with Δ = 0.1, a_{out}/a_{0} = 1 and power index s = −1.5, for several colatitudes θ. Values are normalized to the central value ψ_{c}. The grid is such that A = π/4 and N = 1000. Infinity stands at ζ = 2, and the edges are located at ζ ≈ 0.1269 and 1 (for ). Although not visible on the figure, the potential is not derivable at both edges (case with ; see Sect. 2). 

Open with DEXTER 
Fig. 3 Same legend as for Fig. 2 but for s = −1. In this case, we used the auxiliary function and Eq. (19). 

Open with DEXTER 
Fig. 4 Same legend as for Fig. 2 but for s = −0.5. 

Open with DEXTER 
Fig. 5 Same legend as for Fig. 2 but for s = 0 corresponding to an homogeneous disk. 

Open with DEXTER 
The time required to solve the above system is linear with N. A full grid is then generated in a time that is proportional to N × M, where M is the number of mesh points in the θdirection. With a classical laptop, this takes about 15s to generate a grid with N^{2} = 10^{6} mesh points (including the update of the inhomogeneous term Q). This is much faster than what can be obtained from the multipole expansion, by a factor proportional to the number of Legendre polynomials to be considered, which is generally very large – typically, several hundreds inside the source (Clement 1974).
Fig. 6 Logarithm of the error relative to exact values for the homogeneous case (same conditions as for Fig. 5). The grid has N = 1000 and M = 100 (signed ζvariable). 

Open with DEXTER 
Accuracy can be properly checked in the homogeneous case since the potential is known in a closed form (Durand 1953; Lass & Blitzer 1983). Figure 6 displays the error relative to the exact value. The radial grid and the disk parameters are the same as for Fig. 5 (we used M = 100 for the sampling in θ). We see that the agreement is very good in the whole space, less than 10^{6} on average for this radial resolution, and it is rather uniform in the whole space (however, with a certain deterioration in the midplane due to edges; see Sect. 2). The accuracy only depends on the discretization of the ODE, that is, only on N (not on M). Increasing N obviously improves potential values.
5. Thickness effects and softening length
The flat disk hypothesis is wellsuited to many studies, both theoretical and numerical ones. Its physical realism is, however, limited. For instance, gas or particles present at each edge would feel an infinite acceleration there (see again Sect. 2). Such singularity would disappear by considering an extra dimension (or with a vanishing density profile at the disk edges). Thickness effects (matter present off the midplane) can therefore be important, not only for thermodynamic reasons, but also from a dynamical point of view. We have not yet investigated the existence of an ODE corresponding to geometrically thin (i.e. nonflat) disks. However, we can reproduce the Newtonian potential due to a thin disk with a one dimensional differential equation of the form of Eq. (8) by a “softened” potential^{3}. In this context, Huré & Pierens (2009) showed that the midplane potential of a vertically stratified disk can be approximated by an equation resembling Eq. (2) provided the modulus k is changed for the “softened modulus”: (20)where λ is called the “softening length” (generally considered as a free parameter). They found that the best expression for λ that preserves the Newtonian character of the disk potential takes the form (21)where h is the local semithickness of the disk, and g is a bounded, slowing varying function of the radius. This result is limited to disks with small aspect ratios, locally, i.e. to geometrically thin disks in the classical sense. This result has confirmed the common idea that λ is a certain fraction of the disk thickness. In the absence of mass density gradients in the direction perpendicular to the disk equatorial plane, g goes through a sharp minimum of e^{1} ≈ 0.368 at R = a (i.e. for u = 1). The sensitivity with stratification appears weak since for a vertical parabolic profile, this minimum is e^{−4/3} ≈ 0.264 and the sharpness of the peak is unchanged. Appendix C reproduces the result preliminarily derived in the limit of g far from kernel singularity (i.e. for u → 0 or ∞). In particular, we find that g tends asymptotically to 3^{ − 1/2} ≈ 0.577 in the homogeneous case, while it is 5^{ − 1/2} ≈ 0.447 in the parabolic case. We conclude that, for a given vertical profile, the function g does not vary much, except in the vicinity of R = a where it decreases by a factor of ~37%. As λ does not depend on the precise shape h(a) provided the disk has a small aspect ratio (see Huré & Pierens 2009), this conclusion seems quite robust. This point is confirmed by the numerical experiments. In the following, we therefore neglect any radial variation of g, and consider that it is a constant whose nominal value will be selected a posteriori (and this assumption will be justified). Now, reversing Eq. (20) as done before and using Eq. (21), we find at a given colatitude (22)where we have set ϵ = h/a. It means that, provided the aspect ratio of the disk ϵ is a constant, the ODE can still be formed. After some algebra and with the same nondimensionalization as considered above, we actually find (23)where: (24)where k_{λ,in} ≡ k_{λ} at a = a_{in}, and k_{λ,out} ≡ k_{λ} at a = a_{out}. We see that S_{λ} only differs from Eq. (14) by the presence of softened modulus k_{λ} (instead of k). Since we have (25)it follows that S_{λ} is no longer singular at the disk edges. As a result, the radial acceleration at the edges becomes finite. The role of λ is to mimic an extra dimension (the vertical one).
The softened potential entering Eq. (23) can then be solved numerically in the whole space, as done for the flat disk. In this process, the inner boundary value is slightly lower and must be changed to (26)By looking at Eqs. (8) and (23), we see that there is a perfect continuity between these two ODEs (with and without softening length) as λ → 0. It means that as λ → 0 (including the value at the origin). Also, it is worth noting that the explicit dependence with λ must be regarded as a dependence with the disk thickness.
Fig. 7 Comparison of different potentials: Newtonian potential of a geometrically disk with ϵ = 0.1 (red circles), Newtonian potential of a flat disk (dashed line), and softened potential with (thin lines). Values are normalized to their central value. In all cases, Δ = 0.1, a_{out}/a_{0} = 1, and s = −1.5 (and the same disk mass). 

Open with DEXTER 
To check the reliability of this approach, we first compared the softened potential determined numerically from Eq. (23) with the Newtonian potential of a geometrically thin disk with the same edges, same surface density (and mass), and semithickness h ∝ a (precisely ϵ = 0.1) (hereafter thin disk configuration A). This reference is computed from the splitting method (Huré 2005), which is very accurate. A typical result, limited to the vicinity of the inner edge, is shown in Fig. 7. We see that the agreement between the two curves is very good. There is a very weak sensitivity to the gparameter. The best agreement is obtained for g ≈ 0.486, which is very close to the average of bounds. Interestingly enough, thickness effects mainly shift the potential curve, leaving the radial gradients almost unchanged (except at the edges). Figure 8 shows the logarithm of the relative error in the entire midplane. The error is low, less than 0.1% inside the disk. With , the agreement is even better at a large distance from the disk (but not inside it). Figure 9 displays the error for all radii and all colatitudes. We conclude that the global solution of the softened ODE is very close to the Newtonian potential of a thin disk. The largest deviations are observed in the neighborhood of the disk surface (i.e. Z ≈ h). This is expected since the softened modulus has been approximated here (but this point can be revisited if necessary).
Fig. 8 Relative deviation between the softened potential determined in the midplane from Eq. (23) with and the Newtonian potential of a geometrically thin disk with same edges and same mass (see also Fig. 7). 

Open with DEXTER 
Fig. 9 Error map. Same as for Fig. 8 but at all colatitudes. 

Open with DEXTER 
Up to now, only disks with a constant aspect ratio and a homogeneous vertical mass density profile have been considered. The generality of the method can be studied easily by considering other disk configurations, in particular, configurations for which the conditions for the derivation of the ODE are violated. We illustrate this point by considering two other thin disk configurations:

configuration B: same as configuration A but the disk is stronglyflared, with h/a = ϵ(a/a_{out})^{0.25} (this is an extreme case if we refer to most diskmodels). In this case, the expression for the softening length is stillvalid, but the conditions for deriving the softened ODE areviolated since λ/h is not a constant;

configuration C: same as configuration A but the mass density profile is parabolic, vertically. In this case, the ODE is fully valid, but g must be set to the appropriate value (in between e^{−4/2} and 5^{−1/2}; see above).
In all cases, the disk is geometrically thin at the outer edge (with ϵ = 0.1), which is required to keep the softening length valid). We performed the same comparisons as for configuration A and the results are summarized in Fig. 10. We find that the relative deviation between the Newtonian potential of the thin disk and the softened potential (found from the softened ODE) typically varies as shown in Fig. 8, and is much less than 1% inside the disk. This is also the case outside the disk, except for configuration B for which the relative error gradually rises from R = a_{out} and reaches 10% at infinity. We therefore conclude that using the softened ODE to mimic the Newtonial potential of a thin disk with a precision of a few percent is fully justified as soon as the disk is geometrically thin at all radii. The best values are reproduced when the flaring angle of the thin disk is close to zero.
Fig. 10 Same legend as for Fig. 8 but for configuration A, B, and C. For clarity, potential values are shifted upward for configuration B, and downward for configuration C. 

Open with DEXTER 
6. Concluding remarks
In this article, we have shown that the gravitational potential of a flat, powerlaw disk can be numerically determined in the whole physical space by solving, for any colatitude, a linear system resulting from the secondorder discretization of a firstorder ordinary differential equation (ODE). The computing time is especially low because it is proportional to the number of mesh points. It is also demonstrated that the softened potential based on the prescription for the softening length by Huré & Pierens (2009) obeys a similar ODE whose solution agrees very well with the potential of a geometrically thin disk with same size, mass, and edges and a constant aspect ratio. We briefly tackle the limits of the method, and conclude that the Newtonian potential is reproduced in the full space, provided i) the disk is geometrically thin at all radii, which justifies the use of λ proposed by Huré & Pierens (2009), and ii) the local flaring remains moderate (i.e. h/a does not vary much with the radius), hence supporting the derivation of the ODE. The implementation of a softening length not only mimics a vertical stratification, but also totally removes edge singularities typical of flat disks with sharp edges. If the surface density does not obey a power law but can be written in the form of a series of power laws (e.g. Taylor expansion), the potential is then obtained easily by superposition (see Paper I). The computing time is then proportional to the number of terms in the series. This approach is especially efficient at generating grids of forces. This offers reliable tools for investigating the dynamics of particles in a system containing a perturbing, gravitating thin disk (e.g. Šubr & Karas 2005) and possibly for linking back various families of trajectories to the mass distribution in the disk. The two components of the gravitational force are found from potential values by finite differences. The radial component can, however, be found directly from the ODE. For a unit mass, this component is simply (27)It would be interesting to investigate the existence of a firstorder differential equation associated with the colatitudinal gradient of the potential, complementary to the radial ODE reported here. This is the subject of ongoing work.
For the concept of softened gravity in point mass systems, see Hockney & Eastwood (1988), and see Adams et al. (1989) for the case of gas disks.
Acknowledgments
It is a pleasure to thank Masters students D. Bernard and J. Lambert. The referee is acknowledged.
References
 Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959 [NASA ADS] [CrossRef] (In the text)
 BisnovatyiKogan, G. S. 1975, SvA, 1, 177 (In the text)
 Clement, M. J. 1974, ApJ, 194, 709 [NASA ADS] [CrossRef] (In the text)
 Durand, E. 1953, Electrostatique, Vol. I, Les distributions (Ed. Masson) (In the text)
 Edgar, R. G. 2007, ApJ, 663, 1325 [NASA ADS] [CrossRef] (In the text)
 Goldreich, P., & Tremaine, S. 1978, ApJ, 222, 850 [NASA ADS] [CrossRef] (In the text)
 Goodman, J., & Evans, N. W. 1999, MNRAS, 309, 599 [NASA ADS] [CrossRef] (In the text)
 Gradshteyn, I. S., & Ryzhik, I. M. 1965, Table of integrals, series and products 4th ed., ed. Yu. V. Geronimus (New York: Academic Press) (In the text)
 Hachisu, I. 1986, ApJS, 61, 479 [NASA ADS] [CrossRef] (In the text)
 Hockney, R. W., & Eastwood, J. W. 1988, Computer simulation using particles, ed. R. W. Hockney, & J. W. Eastwood (In the text)
 Huré, J., & Hersant, F. 2007, A&A, 467, 907 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Huré, J., & Pierens, A. 2009, A&A, 507, 573 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Huré, J., Hersant, F., Carreau, C., & Busset, J. 2008, A&A, 490, 477 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Huré, J.M. 2005, A&A, 434, 1 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Huré, J.M., & Pierens, A. 2005, ApJ, 624, 289 [NASA ADS] [CrossRef] (In the text)
 Lass, H., & Blitzer, L. 1983, Cel. Mech., 30, 225 [NASA ADS] [CrossRef] (In the text)
 Mestel, L. 1963, MNRAS, 126, 553 [NASA ADS] [CrossRef] (In the text)
 Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] (In the text)
 Šubr, L., & Karas, V. 2005, in RAGtime 6/7: Workshops on black holes and neutron stars, ed. S. Hledík, & Z. Stuchlík, 281 (In the text)
 Šubr, L., Karas, V., & Huré, J. 2004, MNRAS, 354, 1177 [NASA ADS] [CrossRef] (In the text)
Appendix A: Derivation of the generalized ODE
Starting from Eq. (7), and assuming θ constant, the derivative of u with respect to k is (A.1)This derivative is therefore a function of only k. The potential in Eq. (2) can then be written as the product of a function of the spherical radius r by an integral over k, namely (A.2)where (A.3)The new integral bounds, k_{in} and k_{out}, correspond to the disk edges a_{in} and a_{out}, and they are found from Eq. (4). In particular, we have (A.4)and (A.5)respectively. For constant colatitude θ, the derivative of ψ with respect to the spherical radius is then given by (A.6)This quantity is just the opposite of the radial acceleration due to the disk. Using an elementary derivation rule (see Paper I), we can calculate the right handside of this equation. We find (A.7)where (A.8)and similarly for k_{in}. When a is held constant, we have (A.9)and so (A.10)It follows that (A.11)where Λ is defined by: (A.12)u_{out} = a_{out}/Rand u_{in} = a_{in}/R. This ordinary differential equation (ODE) is the generalization to the whole space of the ODE reported in Huré & Hersant (2007), which was valid only in the disk plane (i.e. the case ). This expression differs mainly by the presence of the spherical radius r (instead of R) and by the presence of the Zdependent modulus k (instead of m). As in Paper I, the second member Λ is an analytical function of a single spatial coordinate r. It now depends on the colatitude θ and on the four parameters s, a_{0}, a_{in}, and a_{out}.
Appendix B: Asymptotic properties
The differential equation possesses the right properties both at small and at large distances. At a large distance from the disk, the modulus in the complete elliptic integrals tends to zero. We have (B.1)and then (B.2)Since the total mass of the disk is given by the expression (B.3)the above ODE can then be rearranged into (B.4)As this equation must be satisfied for any s, we must have (B.5)which is the expected behavior (the disk is no longer distinguishable from a point mass). This also implies that (B.6)At a short distance around the origin (i.e. r ≪ a_{in}), we can perform a second order expansion of the Sterm by expanding the elliptic integral accordingly (Gradshteyn & Ryzhik 1965). We find (B.7)and then the ODE becomes (B.8)whose solution is (B.9)At second order, the potential in the inner domain (r ≪ a_{in}) is quadratic with the cylindrical radius R (while the gravitational acceleration is linear).
Appendix C: The softening length at large relative separations
As in Huré & Pierens (2009), we consider a vertical stratification of the form (C.1)for z ≤ h (and 0 elsewhere), where ρ_{0} is the density at the disk midplane, h the local semithickness (both a function of the radius a in general), and q ≥ 1 an integer. Homogeneous profiles correspond to q → ∞, whereas the parabolic profile is obtained for q = 1. We can compute the contribution of vertical stratification to the potential from the integral (C.2)where (C.3)The approximations k ≈ 1 and K(k) ≈ ln4/k′, assumed in Huré & Pierens (2009), is only valid for a ≈ R. To derive the asymptotic behavior at large relative separation, it is sufficient to consider the case k → 0, that is, at secondorder (C.4)as well as . We then find (C.5)and so (C.6) The softening length is then deduced by equating , where m_{s} is the softened modulus: (C.7)We then have (C.8)
All Figures
Fig. 1 Configuration for the flat, finite size disk with radius at the edges a_{in} and a_{out}; r and θ are spherical coordinates, R and Z are cylindrical coordinates. 

Open with DEXTER  
In the text 
Fig. 2 Potential versus ζ due to a flat disk with Δ = 0.1, a_{out}/a_{0} = 1 and power index s = −1.5, for several colatitudes θ. Values are normalized to the central value ψ_{c}. The grid is such that A = π/4 and N = 1000. Infinity stands at ζ = 2, and the edges are located at ζ ≈ 0.1269 and 1 (for ). Although not visible on the figure, the potential is not derivable at both edges (case with ; see Sect. 2). 

Open with DEXTER  
In the text 
Fig. 3 Same legend as for Fig. 2 but for s = −1. In this case, we used the auxiliary function and Eq. (19). 

Open with DEXTER  
In the text 
Fig. 4 Same legend as for Fig. 2 but for s = −0.5. 

Open with DEXTER  
In the text 
Fig. 5 Same legend as for Fig. 2 but for s = 0 corresponding to an homogeneous disk. 

Open with DEXTER  
In the text 
Fig. 6 Logarithm of the error relative to exact values for the homogeneous case (same conditions as for Fig. 5). The grid has N = 1000 and M = 100 (signed ζvariable). 

Open with DEXTER  
In the text 
Fig. 7 Comparison of different potentials: Newtonian potential of a geometrically disk with ϵ = 0.1 (red circles), Newtonian potential of a flat disk (dashed line), and softened potential with (thin lines). Values are normalized to their central value. In all cases, Δ = 0.1, a_{out}/a_{0} = 1, and s = −1.5 (and the same disk mass). 

Open with DEXTER  
In the text 
Fig. 8 Relative deviation between the softened potential determined in the midplane from Eq. (23) with and the Newtonian potential of a geometrically thin disk with same edges and same mass (see also Fig. 7). 

Open with DEXTER  
In the text 
Fig. 9 Error map. Same as for Fig. 8 but at all colatitudes. 

Open with DEXTER  
In the text 
Fig. 10 Same legend as for Fig. 8 but for configuration A, B, and C. For clarity, potential values are shifted upward for configuration B, and downward for configuration C. 

Open with DEXTER  
In the text 