© 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/(r² + ϵ²)^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 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 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 $\mathrm{\int }\frac{\mathit{\partial \rho }}{\mathit{\partial t}}\hspace{0.17em}\mathrm{d}\mathit{z}\mathrm{+}\mathrm{\int }\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{z}\mathrm{=}\mathrm{0}\mathit{,}$(1)where v and ρ are the 3D velocity and density, respectively. This is typically written as $\frac{\mathit{\partial }\mathrm{\Sigma }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathrm{\Sigma }{u}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}$(2)where $\mathrm{\Sigma }\mathrm{=}\mathrm{\int }\mathit{\rho }\hspace{0.17em}\mathrm{d}\mathit{z}$(3)is the surface density of the disk and ${u}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Sigma }}\mathrm{\int }\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{z}$(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 $\mathrm{\Sigma }\frac{\mathrm{d}{u}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\mathrm{\nabla }\mathit{P}\mathrm{-}\mathrm{\int }\mathit{\rho }\mathrm{\nabla \Psi }\hspace{0.17em}\mathrm{d}\mathit{z,}$(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. $\mathit{f}\mathrm{=}\mathrm{-}\frac{\mathrm{\int }\mathit{\rho }\mathrm{\nabla \Psi }\hspace{0.17em}\mathrm{d}\mathit{z}}{\mathrm{\Sigma }}\mathrm{·}$(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 $\mathrm{\Psi }\mathrm{=}{\mathrm{\Psi }}_{\mathrm{\ast }}\mathrm{+}{\mathrm{\Psi }}_{\mathrm{p}}\mathrm{+}{\mathrm{\Psi }}_{\mathrm{sg}}\mathit{.}$(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_∗/(r² + z²)^1/2, the vertical hydrostatic equation for long distances r from the star can be written as $\frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathit{\partial p}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{z}}{{\left({\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}}\mathrm{\approx }\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{z}}{{\mathit{r}}^{\mathrm{3}}}\mathrm{\equiv }\mathrm{-}{\mathrm{\Omega }}_{\mathrm{K}}^{\mathrm{2}}\mathit{z}\mathit{.}$(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 $\mathit{p}\mathrm{=}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{,}$(9)where the (isothermal) sound speed c_s is now independent of z, we obtain for the vertical disk structure a Gauss profile ${\mathit{\rho }}_{\mathrm{G}}\mathrm{\left(}\mathit{r,\varphi ,z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{r,\varphi }\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{z}}^{\mathrm{2}}}{{\mathit{H}}^{\mathrm{2}}}}\mathit{,}$(10)where the vertical height H is defined as $\mathit{H}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}}{{\mathrm{\Omega }}_{\mathrm{K}}}\mathit{,}$(11)and ρ₀(r,ϕ) is the density in the midplane of the disk. For the vertically integrated surface density Σ one obtains in this case $\mathrm{\Sigma }\mathrm{=}\mathrm{\int }{\mathit{\rho }}_{\mathrm{0}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{z}}^{\mathrm{2}}}{{\mathit{H}}^{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}\mathit{z}\mathrm{=}\sqrt{\mathrm{2}\mathit{\pi }}\mathit{H}{\mathit{\rho }}_{\mathrm{0}}\mathit{.}$(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 M_p embedded in a 3D locally isothermal disk with a fixed Gaussian density profile, ρ_G. The planetary potential is given by ${\mathrm{\Psi }}_{\mathrm{p}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}}{\left|{r}\mathrm{-}{r}\mathrm{p}\right|}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}}{{\left({\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}}\mathit{,}$(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 ${s}\mathrm{=}\mathrm{(}{r}\mathrm{p}\mathrm{-}{r}\mathrm{)}\mathrm{-}\mathrm{⟨}{r}\mathrm{p}\mathrm{-}{r}\mathit{,}{e}\mathit{z}\mathrm{⟩}\hspace{0.17em}{e}\mathit{z}\mathit{,}$(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 ${\mathit{F}}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{\int }\mathit{\rho }\frac{\mathit{\partial }{\mathrm{\Psi }}_{\mathrm{p}}}{\mathit{\partial s}}\hspace{0.17em}\mathrm{d}\mathit{z}\mathrm{=}\mathrm{-}\mathit{G}{\mathit{M}}_{\mathrm{p}}\mathit{s}\mathrm{\int }\frac{\mathit{\rho }}{\mathrm{(}{\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}\mathit{z}\mathit{.}$(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 ${\mathit{\rho }}_{\mathrm{G}}\mathrm{\left(}\mathit{r,\varphi ,z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{r,\varphi }\mathrm{\right)}\mathrm{·}{\mathit{\rho }}_{\mathit{z}}\mathrm{(}{\mathit{z}}^{\mathrm{2}}\mathit{/}{\mathit{H}}^{\mathrm{2}}\mathrm{)}\mathit{,}$(16)where we have assumed that the vertical dependence of ρ is a function of the quantity z²/H², as stated in Eq. (10). Substituting y = z/s in Eq. (15) and using a vertical stratification according to Eq. (16), we find ${\mathit{F}}_{\mathrm{p}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}{\mathit{\rho }}_{\mathrm{0}}}{\mathit{s}}\mathrm{·}\mathrm{2}{\mathit{I}}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{.}$(17)The dimensionless function I_p(s) is defined through an integral over half the disk ${\mathit{I}}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{{\mathit{\rho }}_{\mathit{z}}\mathrm{\left(}{\mathit{c}}^{\mathrm{2}}{\mathit{y}}^{\mathrm{2}}\mathrm{\right)}}{{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{y}}^{\mathrm{2}}\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}\mathit{y,}$(18)where the normalized vertical distance y and the quantity c are given by $\mathit{y}\mathrm{=}\frac{\mathit{z}}{\mathit{s}} \mathrm{and} \mathit{c}\mathrm{=}\frac{\mathit{s}}{\mathit{H}}\mathrm{·}$(19)For the standard Gaussian vertical profile, i.e. ${\mathit{\rho }}_{\mathit{z}}\mathrm{\left(}{\mathit{c}}^{\mathrm{2}}{\mathit{y}}^{\mathrm{2}}\mathrm{\right)}\mathrm{=}\exp {}^{\mathrm{(}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}{\mathit{y}}^{\mathrm{2}}{}^{\mathrm{)}}$, the integral I_p(s) can be expressed as ${\mathit{I}}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}{\mathit{c}}^{\mathrm{2}}\exp \left(\frac{{\mathit{c}}^{\mathrm{2}}}{\mathrm{4}}\right)\left[{\mathit{K}}_{\mathrm{1}}\left(\frac{{\mathit{c}}^{\mathrm{2}}}{\mathrm{4}}\right)\mathrm{-}{\mathit{K}}_{\mathrm{0}}\left(\frac{{\mathit{c}}^{\mathrm{2}}}{\mathrm{4}}\right)\right]\mathit{,}$(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 $\frac{\mathrm{d}{u}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\mathrm{\nabla }{\mathrm{\Psi }}_{\mathrm{p}}^{\mathrm{2}\mathrm{D}}\mathit{.}$(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 $\mathrm{\nabla }{\mathrm{\Psi }}_{\mathrm{p}}^{\mathrm{2}\mathrm{D}}\hspace{0.17em}\ne \hspace{0.17em}\frac{\mathrm{\int }\mathit{\rho }\mathrm{\nabla \Psi }\hspace{0.17em}\mathrm{d}\mathit{z}}{\mathrm{\Sigma }}\mathrm{·}$(22)For ${\mathrm{\Psi }}_{\mathrm{p}}^{\mathrm{2}\mathrm{D}}$ a simple smoothed version is often used, which reads ${\mathrm{\Psi }}_{\mathrm{p}}^{\mathrm{2}\mathrm{D}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}}{{\mathrm{(}{\mathit{s}}^{\mathrm{2}}\mathrm{+}{{\mathit{ϵ}}_{\mathrm{p}}^{\mathrm{2}}}^{\mathrm{)}}}^{\frac{\mathrm{1}}{\mathrm{2}}}}\mathit{,}$(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.

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.

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 ${\mathit{I}}_{\mathit{ϵ}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\frac{{\mathit{s}}^{\mathrm{2}}}{{\left({\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{ϵ}}^{\mathrm{2}}\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}}\frac{\sqrt{\mathrm{2}\mathit{\pi }}\mathit{H}}{\mathrm{2}}\mathrm{·}$(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² 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 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.

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 ${\mathit{F}}_{\mathrm{p}}^{\mathrm{G}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}\mathrm{\Sigma }}{\mathit{s}}\frac{\mathrm{1}}{\mathit{H}}\hspace{0.17em}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\hspace{0.17em}{\mathit{I}}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{.}$(25)This expression could in principle be used directly in numerical simulations of planet-disk interaction. However, even though I_p(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 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 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/s² 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 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 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 $\underset{\mathit{s}\mathrm{\to }\mathrm{\infty }}{\lim }{\mathit{ϵ}}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\mathit{H,}$(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, $\frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathit{\partial p}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{z}}{{\left({\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}\mathit{z}}{{\left({\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}}\mathit{,}$(27)where s is again the distance from the planet. For a vertically isothermal disk this can be integrated ${\mathit{\rho }}_{\mathrm{p}}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\exp \left\{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{3}}}{\mathit{z}}^{\mathrm{2}}\mathrm{+}\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{s}}\left[\frac{\mathit{s}}{\mathrm{(}{\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{-}\mathrm{1}\right]\right\}\mathit{,}$(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 ${\mathit{\rho }}_{\mathrm{p}}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\exp \left\{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{z}}^{\mathrm{2}}}{{\mathit{H}}^{\mathrm{2}}}\mathrm{+}\mathit{q}\frac{{\mathit{r}}^{\mathrm{3}}}{{\mathit{H}}^{\mathrm{2}}\mathit{s}}\left[\frac{\mathit{s}}{\mathrm{(}{\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{-}\mathrm{1}\right]\right\}\mathit{.}$(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 ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\exp \left\{\mathrm{-}\left(\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{z}}^{\mathrm{2}}}{{\mathit{H}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{\left|}\mathit{z}\mathrm{\right|}}{{\mathit{H}}_{\mathrm{p}}}\right)\right\}\mathit{,}$(30)where we define the reduced scale height near the planet ${\mathit{H}}_{\mathrm{p}}\mathrm{=}\frac{\mathrm{4}{\mathit{s}}^{\mathrm{2}}{\mathit{H}}^{\mathrm{2}}}{\mathit{q}{\mathit{r}}^{\mathrm{3}}}\mathrm{·}$(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 $\frac{{\mathit{s}}_{\mathrm{t}}}{\mathit{r}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\hspace{0.17em}{\left(\frac{\mathit{q}}{\mathit{h}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(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, ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}$ 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 ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}$ in the force calculations yields a very good approximation to the exact case.

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 ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}$ 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 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 ${\mathit{F}}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{G}{\mathit{M}}_{\mathrm{p}}\mathit{s}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{z}}_{\max }}\frac{{\mathit{\rho }}_{\mathrm{p}}}{\mathrm{(}{\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}\mathit{z}\mathit{.}$(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 ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}$ 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 ${\mathit{H}}_{\mathrm{eff}}\mathrm{=}{\left(\frac{\mathrm{1}}{{\mathit{H}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{H}}_{\mathrm{p}}^{\mathrm{2}}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(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 ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}$ 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 ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}$ 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 H_p.

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 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 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 s_t ≈ 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 ${\mathrm{\Gamma }}_{\mathrm{tot}}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{\int }\frac{\mathrm{d\Gamma }}{\mathrm{d}\mathit{m}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathit{r}\mathrm{d}\mathit{r}\mathit{.}$(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 ${\left(\frac{\mathrm{d\Gamma }}{\mathrm{d}\mathit{m}}\right)}_{\mathrm{0}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{p}}^{\mathrm{2}}\mathrm{\left(}{\mathit{a}}_{\mathrm{p}}\mathrm{\right)}\hspace{0.17em}{\mathit{a}}_{\mathrm{p}}^{\mathrm{2}}{\mathit{q}}^{\mathrm{2}}{\left(\frac{\mathit{H}}{{\mathit{a}}_{\mathrm{p}}}\right)}^{-4}\mathit{,}$(36)where the index p denotes that the quantities are evaluated at the planet’s position, with the semi-major 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 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.

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 ${\mathit{F}}_{\mathrm{p}}^{\mathrm{G}}$ (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)₀ 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 ${\mathit{F}}_{\mathrm{p}}^{\mathrm{G}}$ 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 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 ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}$ 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 ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}$ 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)₀. 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 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.

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 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 ${\mathrm{\Psi }}_{\mathrm{sg}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\mathrm{Disk}}\frac{\mathit{G\rho }\mathrm{\left(}{r}\mathrm{\prime }\mathrm{\right)}}{\left|{r}\mathrm{-}{r}\mathrm{\prime }\right|}\hspace{0.17em}\mathrm{d}{r}\mathrm{\prime }\mathit{.}$(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 ${\mathrm{\Psi }}_{\mathrm{sg}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\$\frac{\mathit{G\rho }\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathit{,}{\mathit{\varphi }}^{\mathrm{\prime }}\mathit{,}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\right)}}{{\left({\mathit{s}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{d}{\mathit{A}}^{\mathrm{\prime }}\mathit{,}$(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 $\begin{array}{ccc}{\mathit{F}}_{\mathrm{sg}}\mathrm{\left(}\mathit{s}\mathrm{\right)}& & \mathrm{=}\mathrm{-}\mathrm{\int }\mathit{\rho }\mathrm{\left(}\mathit{r,\varphi ,z}\mathrm{\right)}\frac{\mathit{\partial }{\mathrm{\Psi }}_{\mathrm{sg}}}{\mathit{\partial s}}\hspace{0.17em}\mathrm{d}\mathit{z}\\ & & \mathrm{=}\mathrm{-}\mathit{Gs}\&\frac{\mathit{\rho }\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathit{,}{\mathit{\varphi }}^{\mathrm{\prime }}\mathit{,}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}\mathit{r,\varphi ,z}\mathrm{\right)}}{{\left({\mathit{s}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{d}{\mathit{A}}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{d}\mathit{z}\mathit{.}\end{array}$(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.

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).

5.1. Unperturbed disk

For a locally isothermal disk and with Eq. (16) we obtain ${\mathit{F}}_{\mathrm{sg}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{G}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{r,\varphi }\mathrm{\right)}\hspace{0.17em}”{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathit{,}{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{A}}^{\mathrm{\prime }}\mathrm{·}\mathrm{2}{\mathit{I}}_{\mathrm{sg}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{,}$(40)where we defined the function I_sg(s) by ${\mathit{I}}_{\mathrm{sg}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\frac{{\mathit{\rho }}_{\mathit{z}}\mathrm{\left(}{\mathit{c}}^{\mathrm{2}}{\mathit{y}}^{\mathrm{2}}\mathrm{\right)}\hspace{0.17em}{\mathit{\rho }}_{\mathit{z}}\mathrm{\left(}{\mathit{c}}^{\mathrm{2}}{\mathit{y}}^{\mathrm{\prime }\mathrm{2}}\mathrm{\right)}}{{\left(\mathrm{1}\mathrm{+}\mathrm{(}\mathit{y}\mathrm{-}{\mathit{y}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}{\mathit{y}}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{d}\mathit{y,}$(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 ${\mathrm{\Psi }}_{\mathrm{sg}}^{\mathrm{2}\mathrm{D}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\mathrm{-}”\frac{\mathit{G}\mathrm{\Sigma }\mathrm{\left(}{r}\mathrm{\prime }\mathrm{\right)}}{{\mathrm{(}{\mathit{s}}^{\mathrm{2}}\mathrm{+}{{\mathit{ϵ}}_{\mathrm{sg}}^{\mathrm{2}}}^{\mathrm{)}}}^{\frac{\mathrm{1}}{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}{\mathit{A}}^{\mathrm{\prime }}\mathit{,}$(42)where ϵ_sg is the smoothing length. The force acting on each disk element is then calculated from the gradient of ${\mathrm{\Psi }}_{\mathrm{sg}}^{\mathrm{2}\mathrm{D}}$.

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.

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 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.

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 $\underset{\mathit{s}\mathrm{\to }\mathrm{\infty }}{\lim }{\mathit{ϵ}}_{\mathrm{sg}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{2}}\mathit{H}\mathit{.}$(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 $\frac{\mathit{\partial }{\mathrm{\Psi }}_{\mathrm{sg}}}{\mathit{\partial z}}\mathrm{=}\mathrm{2}\mathit{\pi G}\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}$(44)and consequently the hydrostatic Eq. (8) changes to $\frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathit{\partial p}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{\pi G}\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}$(45)For a vertically isothermal disk this can be integrated (Spitzer 1942) to ${\mathit{\rho }}_{\mathrm{sg}}\mathrm{\left(}\mathit{r,\varphi ,z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{r,\varphi }\mathrm{\right)}\hspace{0.17em}{\cosh }^{-2}\left(\frac{\mathit{z}}{{\mathit{H}}_{\mathrm{sg}}}\right)\mathit{,}$(46)where the vertical scale height H_sg is defined by ${\mathit{H}}_{\mathrm{sg}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{\mathit{\pi G}\mathrm{\Sigma }}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}{\mathrm{\Omega }}_{\mathrm{K}}}{\mathit{\pi G}\mathrm{\Sigma }}\mathit{H}\cong \mathit{QH,}$(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 = H_sg 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/r², 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 self-gravitating 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 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%.

Now we can calculate the force correction function I_sg(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

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 $\frac{\mathit{\partial E}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\frac{\mathit{E}\mathrm{\Omega }}{\mathit{\beta }}\mathit{,}$(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 $\mathit{Q}\mathrm{=}\frac{\mathit{\kappa }{\mathit{c}}_{\mathrm{s}}}{\mathit{\pi G}\mathrm{\Sigma }}\mathit{,}$(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 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   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 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)₀, 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 F_p in Eq. (33), which reads ${\mathit{F}}_{\mathrm{p}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{G}{\mathit{M}}_{\mathrm{p}}\mathit{s}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{z}}_{\max }}\frac{{\mathit{\rho }}_{\mathrm{p}}}{\mathrm{(}{\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}\mathit{z}\mathit{.}$(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 ${\mathrm{\int }}_{\mathrm{0}}^{{\mathit{z}}_{\max }}\frac{{\mathit{\rho }}_{\mathrm{p}}\hspace{0.17em}\mathrm{d}\mathit{z}}{\mathrm{(}{\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}\mathrm{\to }\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{z}}}\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{p}}\mathrm{\left(}{\mathit{z}}_{\mathit{k}}\mathrm{\right)}\hspace{0.17em}\mathrm{\Delta }\mathit{z}}{\mathrm{(}{\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}\mathit{,}$(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 $\mathrm{\Sigma }\mathrm{=}\mathrm{2}\hspace{0.17em}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{z}}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{p}}\mathrm{\left(}{\mathit{z}}_{\mathit{k}}\mathrm{\right)}\hspace{0.17em}\mathrm{\Delta }\mathit{z,}$(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 ρ₀. Because the vertical number of grid points is very small (N_z = 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 ${\mathit{f}}_{\mathrm{taper}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\exp \left[\mathrm{-}\mathrm{(}\mathit{s}\mathrm{-}{\mathit{r}}_{\mathrm{t}}\mathrm{)}\mathit{/}\mathrm{\left(}\mathrm{0.2}{\mathit{r}}_{\mathrm{t}}\mathrm{\right)}\right]\mathrm{+}\mathrm{1}}$(A.4)with the tapering cutoff-length 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 ${\mathrm{\Psi }}_{\mathrm{\ast }}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{\left|{r}\mathrm{-}{r}\mathrm{\ast }\right|}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{{\left({\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}}\mathit{,}$(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 ${\mathit{\rho }}_{\mathit{z}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}\left[\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{z}}{{\mathit{H}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}}\right)}^{\mathrm{2}}\right]$(B.1)for the parabolic form, and ${\mathit{\rho }}_{\mathit{z}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\mathrm{=}\left[\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{z}}{{\mathit{H}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}}\right)}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{16}}{\left(\frac{\mathit{z}}{{\mathit{H}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}}\right)}^{\mathrm{4}}\right]$(B.2)for the quartic form. The vertical heights H⁽²⁾ and H⁽⁴⁾ of the models are specified such that the corresponding surface and midplane densities match the isothermal case, see Eq. (12). We obtain ${\mathit{H}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}\mathrm{3}\mathit{/}\mathrm{4}\sqrt{\mathit{\pi }}\hspace{0.17em}\mathit{H}$ and ${\mathit{H}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\mathrm{=}\mathrm{15}\mathit{/}\mathrm{16}\sqrt{\mathit{\pi }\mathit{/}\mathrm{2}}\hspace{0.17em}\mathit{H}$. 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 ${\mathit{\rho }}_{\mathit{z}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ and fourth-order vertical profile ${\mathit{\rho }}_{\mathit{z}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}$. 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. Ip(s) is the numerically calculated force correction function for the Gaussian vertical profile ρz. ${\mathit{I}}_{\mathrm{p}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathit{s}\mathrm{\right)}$ and ${\mathit{I}}_{\mathrm{p}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\mathrm{\left(}\mathit{s}\mathrm{\right)}$ are the analytically calculated force correction functions for the parabolic vertical profile ${\mathit{\rho }}_{\mathit{z}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ and the fourth-order vertical profile ${\mathit{\rho }}_{\mathit{z}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}$.

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 ${\mathit{z}}_{\max }^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}\sqrt{\mathrm{2}}\hspace{0.17em}{\mathit{H}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ and ${\mathit{z}}_{\max }^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\mathrm{=}\mathrm{2}{\mathit{H}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}$. For the second-order integral we find ${\mathit{I}}_{\mathrm{p}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{c}}_{\mathrm{2}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{c}}_{\mathrm{2}}^{\mathrm{2}}\hspace{0.17em}\mathrm{arcsinh}\left(\frac{\sqrt{\mathrm{2}}}{{\mathit{c}}_{\mathrm{2}}}\right)\mathit{,}$(B.3)and for the fourth-order integral ${\mathit{I}}_{\mathrm{p}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}{\mathit{c}}_{\mathrm{4}}^{\mathrm{2}}\mathrm{+}\mathrm{8}}{\mathrm{16}}\sqrt{{\mathit{c}}_{\mathrm{4}}^{\mathrm{2}}\mathrm{+}\mathrm{4}}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{{\mathit{c}}_{\mathrm{4}}^{\mathrm{2}}}{\mathrm{2}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{3}{\mathit{c}}_{\mathrm{4}}^{\mathrm{2}}}{\mathrm{16}}\right)\hspace{0.17em}\mathrm{arcsinh}\frac{\mathrm{2}}{{\mathit{c}}_{\mathrm{4}}}\mathrm{·}$(B.4)Here, we defined c₂ = s/H⁽²⁾ and c₄ = s/H⁽⁴⁾. 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 (${\mathit{I}}_{\mathrm{p}}\mathit{,}{\mathit{I}}_{\mathit{ϵ}}\mathit{,}{\mathit{I}}_{\mathrm{p}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ and ${\mathit{I}}_{\mathrm{p}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}$) approach each other. These simpler profiles may be also useful in the study of self-gravitating 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 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.
In the text

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.

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 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 ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}$ 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 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.
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 ${\mathit{F}}_{\mathrm{p}}^{\mathrm{G}}$ (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 ${\mathit{\rho }}_{\mathrm{p}}^{\mathrm{a}}$ 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 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).

In the text

	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.
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 self-gravitating 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 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%.
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 ${\mathit{\rho }}_{\mathit{z}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ and fourth-order vertical profile ${\mathit{\rho }}_{\mathit{z}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}$. 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. Ip(s) is the numerically calculated force correction function for the Gaussian vertical profile ρz. ${\mathit{I}}_{\mathrm{p}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathit{s}\mathrm{\right)}$ and ${\mathit{I}}_{\mathrm{p}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\mathrm{\left(}\mathit{s}\mathrm{\right)}$ are the analytically calculated force correction functions for the parabolic vertical profile ${\mathit{\rho }}_{\mathit{z}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ and the fourth-order vertical profile ${\mathit{\rho }}_{\mathit{z}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}$.

In the text