Free Access
Issue
A&A
Volume 541, May 2012
Article Number A123
Number of page(s) 13
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201118737
Published online 11 May 2012

© ESO, 2012

1. Introduction

Numerical simulations of accretion disks are often performed in the two-dimensional (2D) thin-disk approximation, because a full 3D treatment with high resolution is still a computationally demanding endeavor requiring a lot of patience. In numerical disk models with embedded planets and/or self-gravity the gravitational potential (force) needs to be smoothed because it diverges for very short mutual distances. For that purpose the potential is softened to avoid the singularities that arise, for example, from pointlike objects, such as planets embedded in disks. In full 3D simulations this smoothing may be required solely for stability purposes and can be chosen to be as small as the given numerical resolution allows. In contrast, 2D disk simulations are typically based on a vertical averaging procedure that leads to a physically required smoothing. Ideally, this smoothing should be calculated in a way that the 2D simulations mimic the 3D case as closely as possible. The most often used potential smoothing has a Plummer form with Ψ ∝ −1/(r2 + ϵ2)1/2, where r is the distance to the gravitating object and ϵ is a suitably chosen smoothing or softening length. If a vertically stratified disk has a thickness H, we would expect that somehow ϵ should depend suitably on H. In a sense, the potential is “diluted” in this case due to the disk’s finite thickness.

For the planet-disk problem, Miyoshi et al. (1999) have shown with local shearing sheet simulations that owing to this dilution effect, the total torque exerted on a planet in a three dimensional disk is only about 43% of the torque obtained in the thin, unsmoothed 2D case. These authors also showed that the strength of the one-sided torque depends on the value of the smoothing, where larger ϵ lead to smaller torques. Later, Masset (2002) has studied the smoothing problem in greater detail. He has shown that good agreement of the total 2D and 3D Lindblad torques can be obtained for smoothing lengths of ϵ = 0.75H, where H is the vertical scale height of the disk, see Eq. (11) below. Masset found that for the Lindblad torques this optimum ϵ/H value is independent of the planet mass and the thickness of the disk. On the contrary, for the corotation torques that are generated by material moving on horseshoe orbits, he found that the required smoothing depends on the ratio RH/H, where RH is the Hill radius of the planet. Masset concluded that there is no “magic” value of ϵ that generates overall agreement of 2D with 3D results.

Based on these studies, in 2D planet-disk simulations the parameter ϵ is typically chosen such that the total torque acting on the planet, which determines the important migration speed, is approximately equal to that obtained through 3D (linear) analysis, for example by Tanaka et al. (2002). This argument has led to the choice of ϵ ≈ 0.3−0.6H, a value very often used in these simulations (Masset 2002; de Val-Borro et al. 2006; Paardekooper & Papaloizou 2009). However, recently a very small smoothing has been advocated for 2D planet-disk simulation (Dong et al. 2011). A smoothing based on a vertical integration using Gaussian density profiles has been used by Li et al. (2005, 2009), but they provided no details on the methodology and accuracy.

For self-gravitating disks, the conditions for fragmentation have recently attracted much attention in the context of planet formation via gravitational instability. Typically, studies are performed in full 3D (see e.g. Meru & Bate (2011b) and references therein). However, to save computational effort, 2D simulations present an interesting alternative in this context (Paardekooper et al. 2011). Here, the accurate treatment of self-gravity is very important.

The incorporation of self-gravity in 2D thin-disk calculations can be achieved by using fast Fourier transforms. Baruteau & Masset (2008) presented a method where the force is calculated directly using a smoothing length that scales linearly with radius, and they used the same smoothing, ϵ = 0.3H, for the planet and self-gravity. Li et al. (2005) used this method to calculate the potential in the disk’s midplane. The simultaneous treatment of an embedded planet and disk-self-gravity can be important because the latter may influence the migration properties of the planet (Pierens & Huré 2005; Baruteau & Masset 2008). For global self-gravitating disks the treatment of the potential has been analyzed in more detail by Huré & Pierens (2009), who calculated the required smoothing by comparing 3D and 2D disk models with specified density stratifications. For that purpose they compared the midplane value of the 3D potential to the 2D case and estimated from this the required smoothing length. Additionally, they considered the whole extent of the disk for their integration. They found that ϵ needs to be reduced for close separations  ≈H, while for long distances it approaches a finite value. Huré & Pierens (2009) give an extended list of smoothing prescriptions used in the literature and we refer the reader to their paper.

Here, we will reanalyze the required smoothing in 2D disk simulations. We show that the force to be used in 2D has to be obtained by performing suitable vertical averages of the force. For the planet-disk case we extend the work by Masset (2002) and compare in detail the torque density to full 3D simulations. We present a method to approximately include the change in vertical stratification induced by the presence of the planet. With respect to self-gravitating disks we follow the work by Huré & Pierens (2009) and calculate the optimum smoothing length ϵ by performing a vertical averaging procedure for the force between two disk elements. We show that this will be important for fragmentation of gravitational unstable disks.

In the next two sections we present the vertical averaging procedure and our unperturbed isothermal disk model. In Sect. 4 we analyze the potential of an embedded planet followed by the self-gravitating case. In Sect. 6 we summarize and conclude.

2. The vertical averaging procedure

Throughout the paper we assume that the disk lies in the z = 0 plane and work in a cylindrical coordinate system (r,ϕ,z). Starting from the full 3D hydrodynamic equations the vertically averaged equations, describing the disk evolution in the r − ϕ plane, are obtained by integrating over the vertical direction. The continuity equation then reads (1)where v and ρ are the 3D velocity and density, respectively. This is typically written as (2)where (3)is the surface density of the disk and (4)the vertically averaged 2D velocity in the disk’s plane. The integrals extend over all z ranging from −∞ to  +∞. The integrated momentum equation, considering only pressure and gravitational forces, then reads (5)with the vertically integrated pressure P = ∫p   dz, where p is the 3D pressure. From Eq. (5) it is obvious that the change in the velocities is determined by the specific force (or acceleration) that is given by the ratio of force density (the integral in Eq. (5)) and surface density Σ, i.e. (6)We here deal with the correct treatment of the last term in Eq. (5), which describes the averaging of the potential that can be given, for example, by a central star, an embedded planet, or self-gravity (7)From the derivation of the momentum Eq. (5) it is clear that not the potential in the midplane matters, but a suitable averaging of the gravitational force. In the following we will analyze this averaging using density stratifications in the vertical isothermal case. Even though accretions disks are known not to be isothermal, we nevertheless prefer at this stage to use this assumption because it is still a commonly used approach that avoids solving a more complex energy equation. It leads to a Gaussian density profile that is also in more general cases a reasonable first-order approximation. Additionally, the irradiation of a central light source tends to make disks more isothermal in the vertical direction. More realistic structures will be treated in future work.

3. The locally isothermal accretion disk model

The vertical structure of accretion disks is determined through the hydrostatic balance in the direction perpendicular to the disk’s midplane. In this section we present first an idealized situation, where the disk is non-self-gravitating and its structure is not influenced by an embedded planet. This allows for analytic evaluation of integrals. Then we will consider more general cases. Using the thin-disk approximation and a gravitational potential generated by a point mass M in the center, i.e. Ψ = −GM/(r2 + z2)1/2, the vertical hydrostatic equation for long distances r from the star can be written as (8)To obtain the vertical stratification for the density, ρ, from this equation requires an equation of state p = p(ρ,T) and a detailed thermal balance, by considering, for example, internal (viscous) heating, stellar irradiation and radiative transport through the disk. To simplify matters, often the so-called locally isothermal approximation is applied for the (numerical) study of embedded planets in disks or for self-gravitating disks. Using this approach, a costly solution of the energy equation is avoided by specifying a priori a fixed radial temperature distribution T(r). At each radius, however, the temperature is assumed to be isothermal in the vertical direction. Using (9)where the (isothermal) sound speed cs is now independent of z, we obtain for the vertical disk structure a Gauss profile (10)where the vertical height H is defined as (11)and ρ0(r,ϕ) is the density in the midplane of the disk. For the vertically integrated surface density Σ one obtains in this case (12)In 2D, r − ϕ, simulations of disks this surface density Σ(r,ϕ) is one of the evolved quantities.

4. The potential of an embedded planet

To illustrate the necessity of potential- (or force-)smoothing we consider in this section embedded planets in locally isothermal disks that are non-self-gravitating. First we will look at an unperturbed disk whose structure is determined by the distance from the host star and is given by Eq. (10). Then we include the effects of an embedded planet on the disk structure.

4.1. Unperturbed disk

We consider the potential of a point like planet with mass Mp embedded in a 3D locally isothermal disk with a fixed Gaussian density profile, ρG. The planetary potential is given by (13)where r is the vector pointing to the location in the disk and rp is the vector pointing to the planet. The vector s pointing from the disk element to the projected position of the planet, i.e. in the z = 0 plane, is denoted by (14)where ez is the unit vector in z direction and  ⟨ ·,· ⟩  denotes the scalar product.

thumbnail Fig. 1

Geometry of protoplanetary disk with embedded planet. We calculated the gravitational force from the planet (blue) on a vertical slice of the disk (red). For that purpose the gradient of the potential has to be vertically integrated along the dashed line that goes through the cell center. To obtain the total force exerted on the shaded disk segment, this value has to be multiplied by its area in the r − ϕ plane.

Open with DEXTER

The force acting on each disk element is then calculated from the gradient of the potential. Since it acts along the vector s, we can write for the modulus of the vertically averaged force density (15)The force density (with units force per area) in Eq. (15) has to be evaluated at the centers of all individual grid cells, as illustrated in Fig. 1. From there, either the total force of a disk element can be calculated by multiplying Fp(s) with the disk element’s area (rΔrΔϕ), or one computes the specific force fp = Fp/Σ, which can be used directly to update the velocities, see Eq. (6). For the density we write (16)where we have assumed that the vertical dependence of ρ is a function of the quantity z2/H2, as stated in Eq. (10). Substituting y = z/s in Eq. (15) and using a vertical stratification according to Eq. (16), we find (17)The dimensionless function Ip(s) is defined through an integral over half the disk (18)where the normalized vertical distance y and the quantity c are given by (19)For the standard Gaussian vertical profile, i.e. , the integral Ip(s) can be expressed as (20)where Kn(x) are the modified Bessel functions of the second kind. For illustrating purposes we present evaluations of the force correction function Ip for simpler polynomial density stratifications in Appendix B.

In 2D numerical simulations of disks the above averaging procedure is typically not performed. Instead, an equivalent 2D potential is used in the momentum equation such that (21)We point out that a 2D potential with the property of equation Eq. (21) cannot be the result of an averaging procedure in general as in Eq. (5) because for realistic densities (22)For a simple smoothed version is often used, which reads (23)where ϵp is the smoothing length to the otherwise point mass potential, introduced to avoid numerical problems at the location of the planet. We will refer to this functional form of Ψ as the ϵ-potential, although it is sometimes named Plummer-potential as well. The force acting on each disk element is then calculated from the gradient of the potential.

A finite ϵp regularises the potential and guarantees stable numerical evolution. Additionally, it serves to account for the vertical stratification of the disk. Comparing torques acting on a planet in 2D and 3D simulations it has been suggested for the Lindblad torques that ϵp should be on the order of the vertical disk height, specifically 0.7H, see Masset (2002). He pointed out, however, that for the corotation torques a lower value may be appropriate. Hence, often a value of ϵ = 0.6H is chosen. Here, we calculated the correct smoothing by a vertical average assuming an isothermal, vertically stratified disk. This will lead to a distance-dependent smoothing.

thumbnail Fig. 2

Force correction function Ip(s) (solid line) resulting from an integration over the vertical structure of the disk, see Eqs. ((17), (20)). Additional curves indicate the corresponding function for the 2D potential (equivalent to Eq. (24)) using different values of the smoothing parameter ϵp.

Open with DEXTER

In Fig. 2 we compare the force correction obtained from the vertically averaging procedure and the 2D smoothed ϵ-potential. Specifically, we plot the function Ip(s) for the Gaussian density profile together with the corresponding function for the 2D potential, which reads (24)Because we are mainly interested in distances s up to a few H, we assume, that the disk height H does not change with radius. The unsmoothed ϵp = 0 potential diverges as 1/s for short distances from the planet, leading to a 1/s2 force as expected for a point mass. Since the value of Ip(s) remains finite for s → 0, we see that a vertically extended disk reduces the divergence of the force to 1/s. This is a clear indication that in 2D simulations the potential has to be smoothed for physical reasons alone, and that the assumption of a point mass potential will greatly overestimate the forces. As expected, the ϵ-potential strongly reduces the force for s → 0 and always yields regular conditions at the location of the planet. Additionally, it appears that the value ϵp = 0.6H overestimates the potential slightly for s ≳ H (green curve in Fig. 2).

thumbnail Fig. 3

Specific gravitational force Fp/Σ as a function of distance from the planet. Different approximations are shown: a) an ideal point mass (solid-red) that falls off as 1/r2; b) the exact averaged force according to Eqs. (25) and (20); c) the force according to the smoothed potential of Eq. (23) using here ϵp = 0.7H; and d) a numerically integrated averaged force using 10 grid points in the vertical direction and a maximum zmax = 5H. The force is normalized such that GMp = 1.

Open with DEXTER

Using the Gaussian isothermal density distribution, ρG, and the surface density Σ from Eq. (12), the force density of the planet acting on the disk can be written as (25)This expression could in principle be used directly in numerical simulations of planet-disk interaction. However, even though Ip(s) can be solved in terms of Bessel-functions, in computational hydrodynamics this evaluation is not very efficient, because it has to be calculated once per timestep at each grid point. A possibility is to solve the integral Ip(s) in Eq. (18) directly numerically using a limited number of vertical grid cells. In this case the integral is converted into a sum, where the exponential can be pre-calculated and tabulated at the corresponding nodes. For our purposes we found that only 10 vertical grid points give an adequate solution (see below). This is shown in Fig. 3 where we plot the specific force (acceleration) exerted by a planet on a disk element that is a distance s away. The force for a point mass falls off as 1/s2 while for small radii the exact vertically averaged force shows a 1/s behavior. Two approximations are displayed as well: the curve for the ϵ-potential given in Eq. (23), where we used a constant ϵp = 0.7H. The numerically averaged curve refers to a numerical integration of Ip(s) using 10 grid points in the z-range  [0,5H] . Note that because of the finite discretization no additional smoothing is required. Increasing the number of grid points increases the agreement with the exact averaged force even more for shorter distances s. The scaling of the distance with 1/H and the normalization of the force make the plot independent of the used vertical thickness of the disk. While the ϵ-potential agrees well for s ≳ H with the exact averaged force, our numerical approximation agrees to much shorter distances. In Appendix B we show that simplified density distributions can lead to an equally good agreement. However, the advantage of the described vertical numerical integration procedure lies in its speed and in the fact that it leads directly to a regularized force for very short s. Also, it is easily generalizable, as will be shown in the next section.

Because we can calculate Ip(s) for the Gaussian profile numerically to any required accuracy, we can estimate the optimum ϵp(s) value for the smoothed potential in Eq. (23) for each point to obtain agreement with the exact averaged force. In Fig. 4 we display the correct ϵp(s) value against the distance. The range of ϵp in which the smoothed potential produces an error of less than 10% in the force compared to the correct value is illustrated by the shaded region. Obviously for very short distances it is important to use the correct ϵp and for long distances the exact value of ϵp does not play an important role. Through a Taylor expansion for s → ∞ of the denominator in Eq. (18) of Ip(s) it can be shown that (26)which can be expected because the disk has a vertical extent on the order H.

thumbnail Fig. 4

Optimum value of ϵp in the smoothing potential of Eq. (23) as a function distance s from the planet. The shaded area illustrates which values of ϵp result in a force error of less than 10% with respect to the exact averaged value.

Open with DEXTER

4.2. Taking the planet into account

thumbnail Fig. 5

Exact (see Eq. (29)) and approximate (see Eq. (30)) vertical density profiles for a disk with aspect ratio h = 0.05 for short distance s to a planet with mass ratio q = 6 × 10-5. The exact profiles are shown by the solid lines and the corresponding approximate ones in the same color but with dashed lines. For comparison the unperturbed Gaussian profile (see Eq. (10)) is shown in black for this H/r.

Open with DEXTER

In the previous section we have analyzed the forces acting on the planet considering an unperturbed disk with a given scale height as determined by the star, see Eq. (11). However, an embedded planet changes the disk structure, which will lead, for example, to a reduced thickness in the vicinity of the planet. This plays an important role for the torques acting on the planet. To estimate this effect we start from the vertical hydrostatic Eq. (8) taking an additional planet into account, (27)where s is again the distance from the planet. For a vertically isothermal disk this can be integrated (28)where we assumed for the stellar contribution z ≪ r, as before. Using the previous vertical thickness H (Eq. (11)) and the mass ratio q = Mp/M, this can be written as (29)Because s is on the same order as z in the neighborhood of the planet, this equation cannot be simplified further. To still obtain an estimate of the expected effects, we approximate the vertical density stratification by (30)where we define the reduced scale height near the planet (31)This approximation for the vertical behavior of ρp leads to the correct limits for the integrated surface density Σ in the limits for very short and long s and is a reasonable approximation in between. From the definition of Hp in Eq. (31) we find that the distance st where the two scale heights are equal, i.e. Hp = H, is given by (32)where h is the relative scale height, h = H/r, of the disk. The location st denotes approximately a transitional distance from the planet. For s ≲ st the standard approximation of a Gaussian density distribution with the scale height H is no longer valid, and the influence of the planet dominates.

Figure 5 shows the exact (Eq. (29)) and the approximate (30) vertical density profiles for different distances to the planet and the comparison with the unperturbed Gaussian profile. For very short distances to the planet the effective height of the disk is much lower than in the unperturbed case. Obviously, is only an approximation to ρp, but it captures the change in thickness of the disk, as induced by the planet, very well. We will show below that using in the force calculations yields a very good approximation to the exact case.

thumbnail Fig. 6

Specific force Fp/Σ for different approximations. The first three curves are identical to Fig. 3 and are shown for comparison. The “exact” potential using ρp is shown in yellow and the approximates solution, using the simplified expression for the density and only 10 vertical grid points, is shown in pink. A tapering function near the center was used for the latter to regularize the force.

Open with DEXTER

To test how this new density-scaling influences the previous isothermal results for the planet’s gravity, we calculated a corrected vertically averaged force using the modified density ρp according to (33)This integral has to be calculated numerically using an approximate zmax. Depending on the model to be calculated (either “exact” or approximate force), we substituted either ρp from Eq. (29) or the approximate from Eq. (30). Since the presence of the planet alters the vertical height of the disk, zmax depends on the distance s from the planet. To estimate zmax we first define an effective vertical scale height (34)which is an interpolation of the value at short and long distances from the planet. We choose to take zmax = 6Heff in the “exact” numerical evaluation where we use the density ρp and 1000 grid cells. In contrast, we apply zmax = 3Heff in the approximate numerical evaluation using and only 10 vertical grid cells, see also Appendix A.

In Fig. 6 we display the results for the vertically integrated specific force Fp/Σ in various approximations. The first three curves (red, blue, green) correspond to those shown in Fig. 3 as an illustration. The “exact” potential using ρp is shown in yellow and the approximate solution using in pink. The presence of the planet reduces the scale height of the disk, which leads to an enhancement of the force above the Gaussian, i.e., it diverges as s-2 for short distances. To regularize the approximate force, we added a tapering function that reduces it to zero in the vicinity of the planet (pink dashed-dotted curve) such that it can be used in hydrodynamic simulations, see Appendix A. Clearly the range of applicability for short s is much improved over the simple ϵ-potential. In contrast to the Gaussian approximation that uses the fixed thickness H, the force correction (with respect to the pure point mass potential) depends now on the mass of the planet as well, which enters through Hp.

4.3. Numerical simulations of planet-disk interactions

To test the formulation of the force, we performed numerical simulations of embedded planets in two and three dimensions where we use an isothermal equation of state. For this purpose we solved the 2D and 3D hydrodynamic equations for a viscous gas. We used a setup very similar to that of Kley & Crida (2008) and Kley et al. (2009). The planet with a mass ratios Mp/M = 6 × 10-5 is embedded at r = 1 in the disk with radial extent of 0.4−2.5. The disk is locally isothermal such that the aspect ratio H/r is constant. This implies a radial temperature gradient of T(r) ∝ r-1, while in the vertical direction T is constant. As a consequence, the unperturbed vertical density structure (in the 3D simulations) is Gaussian along the z-axis. The surface density falls off with radius as Σ ∝ r−1/2. We used a constant kinematic viscosity coefficient of ν = 10-5 in dimensionless units. This setup is such that without the planet the disk is exactly in equilibrium and would not evolve with time. At the radial boundaries damping boundary conditions were applied (de Val-Borro et al. 2006). In our simulations we used H/r = 0.05 and evolved the disk after the planet’s insertion for about 100 orbits. To test our improved treatment of the force, we varied the planet mass and the scale height of the disk in comparison models. For the standard parameter, h = 0.05 and q = 6 × 10-5, we found for the transition distance st ≈ 0.35H.

The embedded planet disturbs the disk and torques are exerted on it by the disk through gravitational back-reaction. These might lead to a change in the planet’s orbital parameter. The strength of these torques will depend on the applied smoothing of the gravitational force. To illustrate the effect, it is convenient to study the radial torque distribution per unit disk mass dΓ(r)/dm, which we define here, following D’Angelo & Lubow (2010), such that the total torque Γtot is given as (35)In other words, dΓ(r) is the torque exerted by a disk annulus of width dr located at the radius r and having the mass dm. As dΓ(r)/dm scales with the mass ratio squared and as (H/r)-4, we rescale our results accordingly in units of (36)where the index p denotes that the quantities are evaluated at the planet’s position, with the semi-major axis ap.

thumbnail Fig. 7

Specific radial torque density in units of (dΓ/dm)0 (see Eq. (36)) for embedded planet models using the ϵ potential with various ϵp. The simulations use a given H = 0.05, and the x-axis refers to the radial coordinate where ap is the semi-major axis of the planet. The planet to star mass-ratio is q = 6 × 10-5. Other parameters of the simulations are stated in Sect. 4.3.

Open with DEXTER

thumbnail Fig. 8

Specific radial torque density in units of (dΓ/dm)0 for embedded planet models for 2D models using ϵ-potential (red curve) and a numerically integrated force (see Eq. (25)) that takes into account a Gaussian vertical density distribution (blue curve). The green curve refers to a full 3D model using the same physical model parameter as the 2D models. The 3D result is adapted from Kley et al. (2009) where the setup and numerics is described in more detail.

Open with DEXTER

In Fig. 7 we plot dΓ(r)/dm obtained from 2D simulations using the ϵ-potential for the gravity of the planet and the standard fixed scale height H for the disk. The torque in this and the similar following plots are scaled to (dΓ(r)/dm)0 as given in Eq. (36). Results for five values of the smoothing length are presented. Obviously, the value of ϵp has great impact on the amplitude of the torque density, and making the correct choice is important. We point out that the differences in the torque density also influence the total torques that determine the important migration rate. Because Γtot consists of positive and negative contributions of similar magnitude, even small errors in dΓ(r)/dm can lead to large errors in Γtot. For the range of ϵ displayed in Fig. 7 we find a variation of Γtot larger than about a factor of 4. The best agreement in this case of H/r = 0.05 with 3D results is obtained for ϵ in between 0.6−0.7H.

In Fig. 8 we compare two different force treatments for 2D simulations to a full 3D run. Results for the ϵp = 0.7H potential are given by the red curve, and the blue curve corresponds to the vertically averaged force according to Eq. (25), which assumes a Gaussian vertical density profile. In magnitude the ϵp = 0.7H potential represents the vertically averaged results reasonably well, but it behaves differently close to the planet. The additional green curve corresponds to a full 3D (locally) isothermal simulation as presented in Kley et al. (2009) (their Fig. 10, purple curve). The 3D runs use an identical physical setup and a more realistic cubic-potential with a smoothing of rsm = 0.5. Both 2D results are on the same order as the 3D run, but show small deviations that can lead to larger variations in the total torque, end hence the migration rate, because the positive and negative contributions to the radial torque density are of similar magnitude.

thumbnail Fig. 9

Specific radial torque density in units of (dΓ/dm)0 for 2D and 3D embedded planet models using different mass ratio and disk temperature. The first three curves represent models with the same disk temperature H/r = 0.05. The red curve corresponds to our standard model, the blue curve has a reduced planet mass q = 3 × 10-5, and the third (green curve) corresponds to a full 3D run. The next two curves (yellow and purple) refer to models with H/r = 0.037. The 2D runs used our approximate density distribution in evaluating the gravitational force. The 3D results were adapted from Kley et al. (2009) where the setup and numerics is described in more detail.

Open with DEXTER

To test the applicability of our new procedure for treating the gravity in 2D simulations, we performed runs with different mass ratio and temperature in the disk. The results of 2D and 3D simulations are shown in Fig. 9, where the torque density dΓ(r)/dm is plotted. The first set of simulations refers to H/r = 0.05, where we compared the standard model in 2D and 3D to a run with half the planet mass. The next two curves show results of 2D and 3D simulations for H/r = 0.037. Despite the indicated differences in the parameter, all models used the same physical setup as described above. The 2D runs used our approximate density distribution in evaluating the gravitational force. The 3D results are adapted from Kley et al. (2009) where the setup and numerics is described in more detail. Firstly, all five curves show very similar overall behavior, confirming the scaling with (dΓ(r)/dm)0. The reduced amplitude of the H/r = 0.037 models is due to the onset of gap formation. Secondly, the agreement of the 2D and 3D runs is very good indeed. For example, upon varying the scale height, the change in shape of the curves is identical in 2D and 3D runs (yellow and purple) curve. We point out that the value of ϵ to obtain the best agreement of the total torque Γtot in 2D and 3D simulations may depend on the value of H, because of the influence of the planet. Hence, it is always advisable to perform the simulations using the vertical integration of the force.

To additionally validate our simulations we compare in Fig. 10 our 2D and 3D results obtained for the standard model to an analytic fit by D’Angelo & Lubow (2010) for the same disk parameter. Although the fit has been developed for a smaller planet mass, the agreement of the two 3D results is excellent. This is interesting because our planet mass of 20 MEarth is already in the range where non-linear effects should set in. The 2D torque shows the same amplitude, but small differences are visible just inside the planet. In isothermal disks the flow close to the planet is not in exact hydrostatic equilibrium anymore. However, full 3D high-resolution isothermal simulations (D’Angelo et al. 2003) have shown that vertical velocities are only large inside the Roche region of the planet. This explains the good agreement of the 2D approximation with the full 3D case.

thumbnail Fig. 10

Specific radial torque density in units of (dΓ/dm)0 for 2D and 3D embedded planet models using different mass ratio and disk temperature. The first two curves are identical to those in the previous figure, and the third is a fit presented by D’Angelo & Lubow (2010) corresponding to the model with T(r) ∝ r-1 and Σ(r) ∝ r−1/2.

Open with DEXTER

5. The potential of a self-gravitating disk

Now we turn to full self-gravitating disks where a smoothing of the gravitational potential is required as well to account for the finite thickness of the disk. The potential at a point r generated by the whole self-gravitating disk is given by (37)The smoothing required to obtain the potential in the midplane has been analyzed for this situation by Huré & Pierens (2009). However, if one is interested in problems of fragmentation, it is more convenient to analyze the force between to individual elements (segments) of the disk. Let us consider the force between two such disk segments that are a separated by the distance s, see Fig. 11. The potential at the location r generated by a disk element located at r′ which is a projected distance s away is given by (38)where dA′ is the surface element of the disk in the r − ϕ-plane at the point r′. The vector s is defined in analogy to Eq. (14) and illustrated in Fig. 11.

The force at the location r generated by this vertically extended disk element is calculated from the gradient of the potential. The vertically averaged force density can then be written in analogy to Eq. (15) as (39) The evaluation of this integral depends on the vertical stratification of the density at both locations r and r′. As before, we consider locally isothermal disks. For weakly self-gravitating disks the density structure is then given by the Gaussian form in Eq. (10). However, similar to an embedded planet the disk’s self-gravity will modify the vertical profile. Following our previous treatment of the embedded planet, we first analyzed the smoothing required for a disk that has an unperturbed Gauss profile and will then allow for modifications.

thumbnail Fig. 11

Geometry of a protoplanetary disk for calculations with self-gravity. We calculate the gravitational force exerted by a vertical slice of the disk (blue) on another vertical slice of the disk (red), separated by a distance s. Two vertical integrations have to be performed along the dashed lines that go through the cell centers. To obtain the total force between the two segments, this value has to be multiplied by the corresponding areas, see also Eq. (39).

Open with DEXTER

5.1. Unperturbed disk

For a locally isothermal disk and with Eq. (16) we obtain (40)where we defined the function Isg(s) by (41)where the normalized vertical distance y and the quantity c are given by Eq. (19) again assuming a constant H. This integral cannot be calculated directly for the standard Gaussian profile, and we will evaluate it numerically.

In 2D numerical simulations usually a simple smoothed potential is used instead of calculating the correct averaging. This ϵsg-potential reads as (42)where ϵsg is the smoothing length. The force acting on each disk element is then calculated from the gradient of .

thumbnail Fig. 12

Force correction function Isg(s) multiplied by s resulting from an integration over the vertical structure of the disk, see Eqs. ((40), (41)). Additional curves indicate the corresponding function for the 2D potential (Eq. (24)) using different values of the smoothing parameter ϵsg.

Open with DEXTER

In Fig. 12 we compare the force correction obtained from the vertically averaging procedure and the 2D smoothed potential for a disk height that does not change with radius. Because Isg(s) diverges for short distances, we multiply it by s. Because we are interested in local effects with distances s up to a few H, we assumed for the plot that the disk height H is constant. Comparing this result to the corresponding force correction function for the planet in Fig. 2, it is obvious that for the self-gravitating case a larger smoothing is required for the ϵ-potential than in the planet case. Here, a value of ϵ = 1.2 seems to be a good choice. Lower values, already ϵ = 0.6, considerably overestimate the force. We attribute this larger required smoothing with the double vertical averaging that has to be performed in this case.

thumbnail Fig. 13

Optimum value of ϵsg in the smoothing potential of Eq. (42) as a function of distance s using a Gaussian vertical stratification. The colored area illustrates which values of ϵsg result in an error of the force correction value of less than 10%.

Open with DEXTER

From the numerically calculated Isg(s) for the Gaussian profile we can calculate the best ϵsg(s) value for the smoothing potential in Eq. (42) at each distance s to obtain the right force correction value. In Fig. 13 we display the optimum ϵsg(s) value versus distance. The range of ϵsg over which the smoothing potential produces an error of less than 10   % is shaded. Cleary for short distances it is crucial to use the correct value of ϵsg, whereas for long distances the influence of ϵsg becomes negligible. Because we now have to account twice for the vertical extent of the disk here compared to the planet case, we obtain a higher limiting value of (43)This value was obtained through numerical calculation up to the sixth significant digit for s = 1000H, and we verified it through a Taylor expansion of the denominator in Eq. (42). As before, it is interesting that the required optimum smoothing remains finite, even for very long distances. This result agrees with Huré & Pierens (2009).

5.2. Taking the disk into account

Now we consider the correction for a disk where self-gravity has modified the vertical structure. For that purpose, we considered the case of a pure self-gravitating disk, neglecting the gravitational potential of the central mass. Assuming a large disk with a slowly varying surface density Σ(r), then the derivate with respect to z of the disk potential given in Eq. (37) simplifies (Mestel 1963) to (44)and consequently the hydrostatic Eq. (8) changes to (45)For a vertically isothermal disk this can be integrated (Spitzer 1942) to (46)where the vertical scale height Hsg is defined by (47)with the Toomre parameter Q (Toomre 1964). Figure 14 shows the changed vertical density profiles in the self-gravitating case compared to the unperturbed Gaussian for H = Hsg or Q = 1, respectively. The self-gravitating profile is steeper and therefore, for equal surface density, the mass is consequently located nearer to the midplane of the disk. This is expected because the vertical component of the disk’s gravitational potential is, with 2πGΣ(r) ~ GM/r2, by a factor of about r/H larger than in Eq. (8), and so the mass is more concentrated toward the midplane.

thumbnail Fig. 14

Vertical density profiles of the disk with for a self-gravitating disk (see Eq. (46)) with Q = 1. For comparison the unperturbed Gaussian profile (see Eq. (10)) is shown in red.

Open with DEXTER

thumbnail Fig. 15

Optimum value of ϵsg in the smoothing potential of Eq. (42) versus distance s using a vertical stratification caused by self-gravity with Q = 1. The colored area illustrates which values of ϵsg result in an error of the force correction value of less than 10%.

Open with DEXTER

Now we can calculate the force correction function Isg(s) of Eq. (41) with the self-gravitating vertical density profile ρsg. In analogy to the unperturbed case, this can be used to calculate an optimum ϵsg for the ϵ-potential. In Fig. 15 we display the correct ϵsg(s) value against distance. The optimum ϵsg for this value of Q is always lower than our previous unperturbed Gaussian case. For the limit s → ∞ we find a value about 10% lower. This is in consistency with the lower effective vertical scale height, because more mass is located near the midplane.

In most self-gravitating disks the vertical structure will be affect by both the central mass object and the self-gravity of the disk. Then it is to be expected that the correct ϵsg is a value between the two extreme cases. The combined situation where self-gravity and the central mass both contribute has been considered by Lodato (2007). For clarity, we treat the two cases separately here.

5.3. Numerical simulations of self-gravitating disks

Table 1

Parameters of the disk model.

To demonstrate the influence of ϵsg in numerical simulations, we analyzed the effects of different values of ϵsg on the fragmentation conclusions of a self-gravitating disk. For that purpose, we adopted a disk model from Baruteau et al. (2011). Table 1 shows all important disk parameters. We simulated the disk twice using the ADSG version of the FARGO code (Baruteau & Masset 2008; Masset 2000) for 5000 years with values of 0.6H and 0.006H for ϵsg. Both models include the self-gravity of the disk and a simple β-cooling model, which is defined by (48)where E is the internal energy, Ω the angular velocity and β a constant. The disk must cool fast enough to be able to fragment, which is otherwise prevented by compressional heating. Gammie (2001) showed that this is the case for β ≲ 3 and Rice et al. (2005) found later a dependency from the equation of state for β and suggested a β ≲ 7 for γ = 5/3.

thumbnail Fig. 16

Radial dependency of surface density (top), midplane temperature (middle) and Toomre parameter (bottom) of the two disk models with ϵsg = 0.6H and ϵsg = 0.006H for different timestamps. At t = 0 a both simulations match because they have identical start conditions. In the first 500 years they behave very similar but then their further evolution diverges drastically.

Open with DEXTER

thumbnail Fig. 17

Surface density of two disks after 5000 years which started from the same initial condition but with different values for the smoothing parameter ϵsg of the gravitational potential. The gray circle shows the computational domain.

Open with DEXTER

Meru & Bate (2011a) pointed out that the previous results had not converged with increasing resolution and that the critical value, βcrit, may be higher than previously thought. They measured an βcrit of  ~ 18 for their highest resolution. Because we use β = 20 in our models, we do not expect them to fragment in any case.

Another stability criterion can be described by the Toomre stability parameter Q, which is defined by (49)where cs is the sound speed in the disk and κ the epicyclic frequency. Toomre (1964) showed that for values of Q ≥ 1 an axisymmetric disk should be stable. As shown in Fig. 16, the modeled disks have an initial Q value of 1.85−2.2 depending on the radial distance to the star, and thus should not be fragmenting.

In the first 500 years both disks cool down in the central parts of the disk from about 120 K to about 30 K, resulting in Q values of    1 at the inner edge of the disk, which means that they are not Toomre-stable anymore, but still should not fragment because the cooling constant β is higher than the critical value required for fragmentation. After 500 years both simulations start to differ. The ϵsg = 0.006H model fragments within about 500 years in the inner region of the disk whereas, the ϵsg = 0.6H model needs about 1000 years to start developing small spiral arms in the inner region. After 5000 years two fragments have survived in the ϵsg = 0.006H disk with fragment masses of 16.5   MJup and 8.3   MJup. The ϵsg = 0.6H model only shows spiral arms and no signs of fragmentation. Even after 50000 years we did not observe any fragments. Figure 17 shows the surface density distribution of both models after 5000 years.

In Sect. 5.1 we suggested an ϵsg on the order of unity to obtain results that can compare to 3D simulations. Because the ϵsg = 0.6H model did not fragment as predicted by the stability criteria, this seems to support the validity of our estimate for ϵsg. The very short ϵsg in the ϵsg = 0.006H model overestimates the gravitational forces on short distances (see Fig. 12) and therefore excites disk fragmentation.

6. Summary and conclusions

We analyzed the smoothing of gravity in thin 2D disk simulations for the embedded planet and self-gravitating case. Starting from the vertically averaged hydrodynamic equations, we first showed that the gravitational force has to be calculated using a density-weighted average of the 3D force, see Eq. (5). Because this depends on the density distribution, there cannot be a general equivalent 2D version of the potential (Eq. (22)). To be able to explicitly calculate the averaged force, we first used a locally isothermal disk structure in which the vertical density stratification is Gaussian.

For the embedded planet case the resulting force can be calculated analytically. In full 2D hydrodynamic simulations we compared the resulting torque density acting on the planet for the ϵ-potential (23) and the “exact” averaged force. We found that the overall magnitude of the torque is best modeled using a smoothing of ϵ = 0.7H, while there remain significant differences to the full 3D case, also in the total torque. Taking the modification of the density stratification induced by the planet into account leads to a much reduced vertical thickness of the disk in the vicinity of the planet. We presented a simplified analytical form for the modified vertical density stratification with a planet, the details of the numerical implementation are given in the appendix. Using this in 2D simulations leads to very good agreement of the torque density with full 3D calculations of embedded planets on circular orbits. By varying the disk height and the mass of the planet, we showed that the torque density scales as expected with (dΓ(r)/dm)0, see Eq. (36).

Because the modified density approximation that includes the planet is based on vertical hydrostatic equilibrium, it is not clear, however, whether it is valid for planets on non-circular or inclined orbits. Here, the variations occur on the orbital timescale and to establish hydrostatic equilibrium, the thermal timescale must be on the same order. In the situation of multiple planets that may interact strongly, the same restrictions may apply. Despite these restrictions, we believe that using this prescription will enhance the accuracy of 2D simulations considerably. We expect that our procedure can be generalized to the radiative case using suitable vertical averages, but this needs to be developed.

For the self-gravitating case we showed that the required smoothing, ϵ ≈ H, is even larger than in the planetary case. We attribute this to the necessity of a double averaging over the vertical height of the disk. Owing to the complex integration, the integrals cannot be solved analytically in this case. Taking into account self-gravity lowers the required smoothing because the vertical scale height is reduced due to the additional gravity. In more strongly self-gravitating systems, which have a Toomre parameter Q ≈ 1, non-axisymmetric features may occur. Because the standard self-gravity solvers require a smoothing that scales with radius, one has to take an approximate average in non-axisymmetric situations. The same applies for disks that are close to the fragmentation limit. As shown by our last example, the choice of smoothing may affect the conclusions on whether the disk will fragment or not. Through detailed comparisons with full 3D simulations a suitable smoothing may be found.

Acknowledgments

We thank Clément Baruteau and Sijme-Jan Paardekooper for useful discussions. Tobias Müller received financial support from the Carl-Zeiss-Stiftung. Farzana Meru and Wilhelm Kley acknowledge the support of the German Research Foundation (DFG) through grant KL 650/8-2 within the Collaborative Research Group FOR 759: The formation of Planets: The Critical First Growth Phase, and Farzana Meru is supported by the ETH Zurich Postdoctoral Fellowship Programm as well as by the Marie Curie Actions for People COFUND program. Most of the simulations were performed on the bwGRiD cluster in Tübingen, which is funded by the Ministry for Education and Research of Germany and the Ministry for Science, Research and Arts of the state Baden-Württemberg, and the cluster of the Forschergruppe FOR 759 “The Formation of Planets: The Critical First Growth Phase” funded by the Deutsche Forschungsgemeinschaft. Finally, we thank the referee for the very constructive and and helpful comments.

References

Appendix A: Integration of the force density

Here, we briefly outline a numerically fast and convenient method to vertically integrate the force for the embedded planet case. Specifically, we plan to evaluate the force density Fp in Eq. (33), which reads (A.1)As pointed out in the text, we used for the (approximate) numerical integration a maximum z of zmax = 3Heff, with Heff given by Eq. (34). The interval  [0,zmax]  is divided into Nz equal intervals with the size Δz = zmax/Nz. The integral in Eq. (A.1) is replaced with the following sum (A.2)where the Nz nodes are located at zk = (k − 1/2)Δz. We introduce a small smoothing, rs, here to keep the sum regular at short distances s. In the simulation presented in the paper we used rs = 0.1RHill throughout, which is much shorter than the unperturbed vertical height H. In the 2D hydrodynamic simulations we used Nz = 10, which makes the method feasible numerically.

The specific force, or acceleration (Fp/Σ), is then obtained by dividing through the integrated surface density (A.3)using the same nodes. In 2D simulations Σ is one of the evolved quantities and this relation can be used to calculate the otherwise unknown midplane density ρ0. Because the vertical number of grid points is very small (Nz = 10), the method is reasonably fast and can be used in numerical planet-disk simulations. Additionally, it only needs to be evaluated in the vicinity of the planet and could be omitted farther out. Nevertheless, despite the coarseness of the integration, the agreement with the “exact” force and torque is excellent. For numerical stability we smooth the resulting force using the following tapering function (A.4)with the tapering cutoff-length rt, for which we used in our simulations rt = 0.2RHill. The purple curve in Fig. 6 exactly corresponds to the procedure described here, using the stated parameter.

Finally, we point out that for the stellar contribution to the potential (A.5)a very similar vertical averaging needs to be performed. However, because the vertical profile is always Gaussian, the nodes can be precomputed and the exponential function has only to be evaluated once for the Nz nodes per grid point and timestep.

Appendix B: Approximate vertical density profiles

To simplify some estimates and obtain an idea of the functional behavior of the forces, it is useful to study simpler vertical density stratifications. Here, we present results for a parabolic and quartic behavior. The corresponding density stratifications read (B.1)for the parabolic form, and (B.2)for the quartic form. The vertical heights H(2) and H(4) of the models are specified such that the corresponding surface and midplane densities match the isothermal case, see Eq. (12). We obtain and . Figure B.1 shows all three vertical density profiles in comparison.

thumbnail Fig. B.1

Gaussian vertical profile ρz (solid line) compared with the parabolic vertical profile and fourth-order vertical profile . They all have the same area below the curves and yield the identical surface density.

Open with DEXTER

thumbnail Fig. B.2

Force correction functions for different vertical density profiles. Ip(s) is the numerically calculated force correction function for the Gaussian vertical profile ρz. and are the analytically calculated force correction functions for the parabolic vertical profile  and the fourth-order vertical profile .

Open with DEXTER

These density stratifications are then used to calculate the vertically integrated force and to obtain the corresponding force correction functions. Here, the vertical integrations extend to that zmax value where the density vanishes. For the two distributions we find and . For the second-order integral we find (B.3)and for the fourth-order integral (B.4)Here, we defined c2 = s/H(2) and c4 = s/H(4). As Fig. B.2 illustrates, for the unperturbed disk (without a planet) these functions agree reasonably well with the Gaussian value also for smaller separations s. For long distances all studied force correction functions ( and ) approach each other. These simpler profiles may be also useful in the study of self-gravitating disks.

All Tables

Table 1

Parameters of the disk model.

All Figures

thumbnail Fig. 1

Geometry of protoplanetary disk with embedded planet. We calculated the gravitational force from the planet (blue) on a vertical slice of the disk (red). For that purpose the gradient of the potential has to be vertically integrated along the dashed line that goes through the cell center. To obtain the total force exerted on the shaded disk segment, this value has to be multiplied by its area in the r − ϕ plane.

Open with DEXTER
In the text
thumbnail Fig. 2

Force correction function Ip(s) (solid line) resulting from an integration over the vertical structure of the disk, see Eqs. ((17), (20)). Additional curves indicate the corresponding function for the 2D potential (equivalent to Eq. (24)) using different values of the smoothing parameter ϵp.

Open with DEXTER
In the text
thumbnail Fig. 3

Specific gravitational force Fp/Σ as a function of distance from the planet. Different approximations are shown: a) an ideal point mass (solid-red) that falls off as 1/r2; b) the exact averaged force according to Eqs. (25) and (20); c) the force according to the smoothed potential of Eq. (23) using here ϵp = 0.7H; and d) a numerically integrated averaged force using 10 grid points in the vertical direction and a maximum zmax = 5H. The force is normalized such that GMp = 1.

Open with DEXTER
In the text
thumbnail Fig. 4

Optimum value of ϵp in the smoothing potential of Eq. (23) as a function distance s from the planet. The shaded area illustrates which values of ϵp result in a force error of less than 10% with respect to the exact averaged value.

Open with DEXTER
In the text
thumbnail Fig. 5

Exact (see Eq. (29)) and approximate (see Eq. (30)) vertical density profiles for a disk with aspect ratio h = 0.05 for short distance s to a planet with mass ratio q = 6 × 10-5. The exact profiles are shown by the solid lines and the corresponding approximate ones in the same color but with dashed lines. For comparison the unperturbed Gaussian profile (see Eq. (10)) is shown in black for this H/r.

Open with DEXTER
In the text
thumbnail Fig. 6

Specific force Fp/Σ for different approximations. The first three curves are identical to Fig. 3 and are shown for comparison. The “exact” potential using ρp is shown in yellow and the approximates solution, using the simplified expression for the density and only 10 vertical grid points, is shown in pink. A tapering function near the center was used for the latter to regularize the force.

Open with DEXTER
In the text
thumbnail Fig. 7

Specific radial torque density in units of (dΓ/dm)0 (see Eq. (36)) for embedded planet models using the ϵ potential with various ϵp. The simulations use a given H = 0.05, and the x-axis refers to the radial coordinate where ap is the semi-major axis of the planet. The planet to star mass-ratio is q = 6 × 10-5. Other parameters of the simulations are stated in Sect. 4.3.

Open with DEXTER
In the text
thumbnail Fig. 8

Specific radial torque density in units of (dΓ/dm)0 for embedded planet models for 2D models using ϵ-potential (red curve) and a numerically integrated force (see Eq. (25)) that takes into account a Gaussian vertical density distribution (blue curve). The green curve refers to a full 3D model using the same physical model parameter as the 2D models. The 3D result is adapted from Kley et al. (2009) where the setup and numerics is described in more detail.

Open with DEXTER
In the text
thumbnail Fig. 9

Specific radial torque density in units of (dΓ/dm)0 for 2D and 3D embedded planet models using different mass ratio and disk temperature. The first three curves represent models with the same disk temperature H/r = 0.05. The red curve corresponds to our standard model, the blue curve has a reduced planet mass q = 3 × 10-5, and the third (green curve) corresponds to a full 3D run. The next two curves (yellow and purple) refer to models with H/r = 0.037. The 2D runs used our approximate density distribution in evaluating the gravitational force. The 3D results were adapted from Kley et al. (2009) where the setup and numerics is described in more detail.

Open with DEXTER
In the text
thumbnail Fig. 10

Specific radial torque density in units of (dΓ/dm)0 for 2D and 3D embedded planet models using different mass ratio and disk temperature. The first two curves are identical to those in the previous figure, and the third is a fit presented by D’Angelo & Lubow (2010) corresponding to the model with T(r) ∝ r-1 and Σ(r) ∝ r−1/2.

Open with DEXTER
In the text
thumbnail Fig. 11

Geometry of a protoplanetary disk for calculations with self-gravity. We calculate the gravitational force exerted by a vertical slice of the disk (blue) on another vertical slice of the disk (red), separated by a distance s. Two vertical integrations have to be performed along the dashed lines that go through the cell centers. To obtain the total force between the two segments, this value has to be multiplied by the corresponding areas, see also Eq. (39).

Open with DEXTER
In the text
thumbnail Fig. 12

Force correction function Isg(s) multiplied by s resulting from an integration over the vertical structure of the disk, see Eqs. ((40), (41)). Additional curves indicate the corresponding function for the 2D potential (Eq. (24)) using different values of the smoothing parameter ϵsg.

Open with DEXTER
In the text
thumbnail Fig. 13

Optimum value of ϵsg in the smoothing potential of Eq. (42) as a function of distance s using a Gaussian vertical stratification. The colored area illustrates which values of ϵsg result in an error of the force correction value of less than 10%.

Open with DEXTER
In the text
thumbnail Fig. 14

Vertical density profiles of the disk with for a self-gravitating disk (see Eq. (46)) with Q = 1. For comparison the unperturbed Gaussian profile (see Eq. (10)) is shown in red.

Open with DEXTER
In the text
thumbnail Fig. 15

Optimum value of ϵsg in the smoothing potential of Eq. (42) versus distance s using a vertical stratification caused by self-gravity with Q = 1. The colored area illustrates which values of ϵsg result in an error of the force correction value of less than 10%.

Open with DEXTER
In the text
thumbnail Fig. 16

Radial dependency of surface density (top), midplane temperature (middle) and Toomre parameter (bottom) of the two disk models with ϵsg = 0.6H and ϵsg = 0.006H for different timestamps. At t = 0 a both simulations match because they have identical start conditions. In the first 500 years they behave very similar but then their further evolution diverges drastically.

Open with DEXTER
In the text
thumbnail Fig. 17

Surface density of two disks after 5000 years which started from the same initial condition but with different values for the smoothing parameter ϵsg of the gravitational potential. The gray circle shows the computational domain.

Open with DEXTER
In the text
thumbnail Fig. B.1

Gaussian vertical profile ρz (solid line) compared with the parabolic vertical profile and fourth-order vertical profile . They all have the same area below the curves and yield the identical surface density.

Open with DEXTER
In the text
thumbnail Fig. B.2

Force correction functions for different vertical density profiles. Ip(s) is the numerically calculated force correction function for the Gaussian vertical profile ρz. and are the analytically calculated force correction functions for the parabolic vertical profile  and the fourth-order vertical profile .

Open with DEXTER
In the text

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.