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/00046361/201118737  
Published online  11 May 2012 
Treating gravity in thindisk simulations
^{1}
Institut für Astronomie & Astrophysik, Universität
Tübingen, Auf der Morgenstelle
10, 72076
Tübingen, Germany
email: Tobias_Mueller@twam.info
^{2}
Institute for Astronomy, ETH Zürich, WolfgangPauliStrasse 27,
8093
Zürich,
Switzerland
Received:
23
December
2011
Accepted:
6
March
2012
Context. In 2Dsimulations of thin gaseous disks with embedded planets or selfgravity the gravitational potential needs to be smoothed to avoid singularities in the numerical evaluation of the gravitational potential or force. The softening prescription used in 2D needs to be adjusted properly to correctly resemble the realistic case of vertically extended 3D disks.
Aims. We analyze the embedded planet and the selfgravity case and provide a method to evaluate the required smoothing in 2D simulations of thin disks.
Methods. Starting from the averaged hydrodynamic equations and using a vertically isothermal disk model, we calculated the force to be used in 2D simulations. We compared our results to the often used Plummer form of the potential, which runs as ∝ 1/(r^{2} + ϵ^{2})^{1/2}. For that purpose we computed the required smoothing length ϵ as a function of distance r to the planet or to a disk element within a selfgravitating disk.
Results. We find that for longer distances ϵ is determined solely by the vertical disk thickness H. For the planet case we find that outside r ≈ H a value of ϵ = 0.7H describes the averaged force very well, while in the selfgravitating disk the value needs to be higher, ϵ = 1.2H. For shorter distances the smoothing needs to be reduced significantly. Comparing torque densities of 3D and 2D simulations we show that the modification to the vertical density stratification as induced by an embedded planet needs to be taken into account to obtain agreeing results.
Conclusions. It is very important to use the correct value of ϵ in 2D simulations to obtain a realistic outcome. In disk fragmentation simulations the choice of ϵ can determine whether a disk will fragment or not. Because a wrong smoothing length can change even the direction of migration, it is very important to include the effect of the planet on the local scale height in 2D planetdisk simulations. We provide an approximate and fast method for this purpose that agrees very well with full 3D simulations.
Key words: accretion, accretion disks / planets and satellites: formation / hydrodynamics / methods: numerical / protoplanetary disks
© ESO, 2012
1. Introduction
Numerical simulations of accretion disks are often performed in the twodimensional (2D) thindisk 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 selfgravity 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/(r^{2} + ϵ^{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 planetdisk 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 onesided 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 R_{H}/H, where R_{H} 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 planetdisk 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 ValBorro et al. 2006; Paardekooper & Papaloizou 2009). However, recently a very small smoothing has been advocated for 2D planetdisk 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 selfgravitating 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 selfgravity is very important.
The incorporation of selfgravity in 2D thindisk 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 selfgravity. Li et al. (2005) used this method to calculate the potential in the disk’s midplane. The simultaneous treatment of an embedded planet and diskselfgravity can be important because the latter may influence the migration properties of the planet (Pierens & Huré 2005; Baruteau & Masset 2008). For global selfgravitating 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 planetdisk 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 selfgravitating 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 selfgravitating 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 selfgravity (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 firstorder 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 nonselfgravitating 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 thindisk approximation and a gravitational potential generated by a point mass M_{∗} in the center, i.e. Ψ_{∗} = −GM_{∗}/(r^{2} + z^{2})^{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 socalled locally isothermal approximation is applied for the (numerical) study of embedded planets in disks or for selfgravitating 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 c_{s} 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 nonselfgravitating. 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 M_{p} 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 r_{p} 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 e_{z} is the unit vector in z direction and ⟨ ·,· ⟩ denotes the scalar product.
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. 
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 F_{p}(s) with the disk element’s area (rΔrΔϕ), or one computes the specific force f_{p} = F_{p}/Σ, 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 z^{2}/H^{2}, 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 I_{p}(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 I_{p}(s) can be expressed as (20)where K_{n}(x) are the modified Bessel functions of the second kind. For illustrating purposes we present evaluations of the force correction function I_{p} 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 Plummerpotential 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 distancedependent smoothing.
Fig. 2 Force correction function I_{p}(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}. 
In Fig. 2 we compare the force correction obtained from the vertically averaging procedure and the 2D smoothed ϵpotential. Specifically, we plot the function I_{p}(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/s^{2} force as expected for a point mass. Since the value of I_{p}(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).
Fig. 3 Specific gravitational force F_{p}/Σ as a function of distance from the planet. Different approximations are shown: a) an ideal point mass (solidred) that falls off as 1/r^{2}; 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 z_{max} = 5H. The force is normalized such that GM_{p} = 1. 
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 planetdisk interaction. However, even though I_{p}(s) can be solved in terms of Besselfunctions, 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 I_{p}(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 precalculated 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/s^{2} 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 I_{p}(s) using 10 grid points in the zrange [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 I_{p}(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 I_{p}(s) it can be shown that (26)which can be expected because the disk has a vertical extent on the order H.
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. 
4.2. Taking the planet into account
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. 
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 = M_{p}/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 H_{p} in Eq. (31) we find that the distance s_{t} where the two scale heights are equal, i.e. H_{p} = H, is given by (32)where h is the relative scale height, h = H/r, of the disk. The location s_{t} denotes approximately a transitional distance from the planet. For s ≲ s_{t} 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.
Fig. 6 Specific force F_{p}/Σ 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. 
To test how this new densityscaling 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 z_{max}. 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, z_{max} depends on the distance s from the planet. To estimate z_{max} 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 z_{max} = 6H_{eff} in the “exact” numerical evaluation where we use the density ρ_{p} and 1000 grid cells. In contrast, we apply z_{max} = 3H_{eff} 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 F_{p}/Σ 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 dasheddotted 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 H_{p}.
4.3. Numerical simulations of planetdisk 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 M_{p}/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 zaxis. 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 ValBorro 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 s_{t} ≈ 0.35H.
The embedded planet disturbs the disk and torques are exerted on it by the disk through gravitational backreaction. 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 semimajor axis a_{p}.
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 xaxis refers to the radial coordinate where a_{p} is the semimajor axis of the planet. The planet to star massratio is q = 6 × 10^{5}. Other parameters of the simulations are stated in Sect. 4.3. 
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. 
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 cubicpotential with a smoothing of r_{sm} = 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.
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. 
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 M_{Earth} is already in the range where nonlinear 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 highresolution 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.
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}. 
5. The potential of a selfgravitating disk
Now we turn to full selfgravitating 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 selfgravitating 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 selfgravitating disks the density structure is then given by the Gaussian form in Eq. (10). However, similar to an embedded planet the disk’s selfgravity 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.
Fig. 11 Geometry of a protoplanetary disk for calculations with selfgravity. 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). 
5.1. Unperturbed disk
For a locally isothermal disk and with Eq. (16) we obtain (40)where we defined the function I_{sg}(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 .
Fig. 12 Force correction function I_{sg}(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}. 
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 I_{sg}(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 selfgravitating 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.
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%. 
From the numerically calculated I_{sg}(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 selfgravity has modified the vertical structure. For that purpose, we considered the case of a pure selfgravitating 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 H_{sg} is defined by (47)with the Toomre parameter Q (Toomre 1964). Figure 14 shows the changed vertical density profiles in the selfgravitating case compared to the unperturbed Gaussian for H = H_{sg} or Q = 1, respectively. The selfgravitating 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/r^{2}, by a factor of about r/H larger than in Eq. (8), and so the mass is more concentrated toward the midplane.
Fig. 14 Vertical density profiles of the disk with for a selfgravitating disk (see Eq. (46)) with Q = 1. For comparison the unperturbed Gaussian profile (see Eq. (10)) is shown in red. 
Fig. 15 Optimum value of ϵ_{sg} in the smoothing potential of Eq. (42) versus distance s using a vertical stratification caused by selfgravity with Q = 1. The colored area illustrates which values of ϵ_{sg} result in an error of the force correction value of less than 10%. 
Now we can calculate the force correction function I_{sg}(s) of Eq. (41) with the selfgravitating 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 selfgravitating disks the vertical structure will be affect by both the central mass object and the selfgravity 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 selfgravity 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 selfgravitating disks
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 selfgravitating 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 selfgravity 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.
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. 
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. 
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 c_{s} 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 Toomrestable 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 M_{Jup} and 8.3 M_{Jup}. 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 selfgravitating case. Starting from the vertically averaged hydrodynamic equations, we first showed that the gravitational force has to be calculated using a densityweighted 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 noncircular 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 selfgravitating 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 selfgravity lowers the required smoothing because the vertical scale height is reduced due to the additional gravity. In more strongly selfgravitating systems, which have a Toomre parameter Q ≈ 1, nonaxisymmetric features may occur. Because the standard selfgravity solvers require a smoothing that scales with radius, one has to take an approximate average in nonaxisymmetric 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 SijmeJan Paardekooper for useful discussions. Tobias Müller received financial support from the CarlZeissStiftung. Farzana Meru and Wilhelm Kley acknowledge the support of the German Research Foundation (DFG) through grant KL 650/82 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 BadenWü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
 Baruteau, C., & Masset, F. 2008, ApJ, 678, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., Meru, F., & Paardekooper, S.J. 2011, MNRAS, 416, 1971 [NASA ADS] [CrossRef] [Google Scholar]
 D’Angelo, G., & Lubow, S. H. 2010, ApJ, 724, 730 [NASA ADS] [CrossRef] [Google Scholar]
 D’Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540 [NASA ADS] [CrossRef] [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Dong, R., Rafikov, R. R., Stone, J. M., & Petrovich, C. 2011, ApJ, 741, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F. 2001, ApJ, 553, 174 [NASA ADS] [CrossRef] [Google Scholar]
 Huré, J.M., & Pierens, A. 2009, A&A, 507, 573 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., & Crida, A. 2008, A&A, 487, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, H., Li, S., Koller, J., et al. 2005, ApJ, 624, 1003 [NASA ADS] [CrossRef] [Google Scholar]
 Li, H., Lubow, S. H., Li, S., & Lin, D. N. C. 2009, ApJ, 690, L52 [NASA ADS] [CrossRef] [Google Scholar]
 Lodato, G. 2007, Nuovo Cimento Riv. Ser., 30, 293 [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masset, F. S. 2002, A&A, 387, 605 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meru, F., & Bate, M. R. 2011a, MNRAS, 411, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Meru, F., & Bate, M. R. 2011b, MNRAS, 410, 559 [NASA ADS] [CrossRef] [Google Scholar]
 Mestel, L. 1963, MNRAS, 126, 553 [NASA ADS] [CrossRef] [Google Scholar]
 Miyoshi, K., Takeuchi, T., Tanaka, H., & Ida, S. 1999, ApJ, 516, 451 [NASA ADS] [CrossRef] [Google Scholar]
 Paardekooper, S.J., & Papaloizou, J. C. B. 2009, MNRAS, 394, 2283 [NASA ADS] [CrossRef] [Google Scholar]
 Paardekooper, S.J., Baruteau, C., & Meru, F. 2011, MNRAS, 416, L65 [NASA ADS] [CrossRef] [Google Scholar]
 Pierens, A., & Huré, J.M. 2005, A&A, 433, L37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rice, W. K. M., Lodato, G., & Armitage, P. J. 2005, MNRAS, 364, L56 [NASA ADS] [CrossRef] [Google Scholar]
 Spitzer, Jr., L. 1942, ApJ, 95, 329 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1964, ApJ, 139, 1217 [NASA ADS] [CrossRef] [Google Scholar]
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 F_{p} in Eq. (33), which reads (A.1)As pointed out in the text, we used for the (approximate) numerical integration a maximum z of z_{max} = 3H_{eff}, with H_{eff} given by Eq. (34). The interval [0,z_{max}] is divided into N_{z} equal intervals with the size Δz = z_{max}/N_{z}. The integral in Eq. (A.1) is replaced with the following sum (A.2)where the N_{z} nodes are located at z_{k} = (k − 1/2)Δz. We introduce a small smoothing, r_{s}, here to keep the sum regular at short distances s. In the simulation presented in the paper we used r_{s} = 0.1R_{Hill} throughout, which is much shorter than the unperturbed vertical height H. In the 2D hydrodynamic simulations we used N_{z} = 10, which makes the method feasible numerically.
The specific force, or acceleration (F_{p}/Σ), 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 (N_{z} = 10), the method is reasonably fast and can be used in numerical planetdisk 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 cutofflength r_{t}, for which we used in our simulations r_{t} = 0.2R_{Hill}. 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 N_{z} 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.
Fig. B.1 Gaussian vertical profile ρ_{z} (solid line) compared with the parabolic vertical profile and fourthorder vertical profile . They all have the same area below the curves and yield the identical surface density. 
Fig. B.2 Force correction functions for different vertical density profiles. I_{p}(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 fourthorder vertical profile . 
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 z_{max} value where the density vanishes. For the two distributions we find and . For the secondorder integral we find (B.3)and for the fourthorder integral (B.4)Here, we defined c_{2} = s/H^{(2)} and c_{4} = 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 selfgravitating disks.
All Tables
All Figures
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. 

In the text 
Fig. 2 Force correction function I_{p}(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}. 

In the text 
Fig. 3 Specific gravitational force F_{p}/Σ as a function of distance from the planet. Different approximations are shown: a) an ideal point mass (solidred) that falls off as 1/r^{2}; 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 z_{max} = 5H. The force is normalized such that GM_{p} = 1. 

In the text 
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. 

In the text 
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. 

In the text 
Fig. 6 Specific force F_{p}/Σ 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. 

In the text 
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 xaxis refers to the radial coordinate where a_{p} is the semimajor axis of the planet. The planet to star massratio is q = 6 × 10^{5}. Other parameters of the simulations are stated in Sect. 4.3. 

In the text 
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. 

In the text 
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. 

In the text 
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}. 

In the text 
Fig. 11 Geometry of a protoplanetary disk for calculations with selfgravity. 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). 

In the text 
Fig. 12 Force correction function I_{sg}(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}. 

In the text 
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%. 

In the text 
Fig. 14 Vertical density profiles of the disk with for a selfgravitating disk (see Eq. (46)) with Q = 1. For comparison the unperturbed Gaussian profile (see Eq. (10)) is shown in red. 

In the text 
Fig. 15 Optimum value of ϵ_{sg} in the smoothing potential of Eq. (42) versus distance s using a vertical stratification caused by selfgravity with Q = 1. The colored area illustrates which values of ϵ_{sg} result in an error of the force correction value of less than 10%. 

In the text 
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. 

In the text 
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. 

In the text 
Fig. B.1 Gaussian vertical profile ρ_{z} (solid line) compared with the parabolic vertical profile and fourthorder vertical profile . They all have the same area below the curves and yield the identical surface density. 

In the text 
Fig. B.2 Force correction functions for different vertical density profiles. I_{p}(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 fourthorder vertical profile . 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.