ON THE INTERACTION BETWEEN A PROTOPLANETARY DISK AND A PLANET IN AN ECCENTRIC ORBIT: APPLICATION OF DYNAMICAL FRICTION

We consider a homogeneous gaseous slab with surface density Σ₀. We take the Cartesian coordinate system (x, y) in the plane of the slab and consider a point particle with mass M traveling in the y-direction at a constant speed v₀. The basic equations we consider are the (vertically integrated) hydrodynamic equations with a simplified prescription of viscosity,

$$\begin{equation} {\displaystyle \frac{\partial \Sigma }{\partial t}} + \nabla \cdot \left(\Sigma \mathbf {v} \right) = 0, \end{equation} \tag{ 1 }$$

$$\begin{equation} {\displaystyle \frac{\partial \mathbf {v}}{\partial t}} + \mathbf {v}\cdot \nabla \mathbf {v} = -\frac{1}{\Sigma }\nabla {P} + \nu \nabla ^2 \mathbf {v} - \nabla \psi, \end{equation} \tag{ 2 }$$

where Σ is the surface density, $\mathbf {v}$ is the velocity, P is the (vertically integrated) pressure, ν is the viscous coefficient, and ψ is the gravitational potential of the particle. We use a simple form of viscosity in order to keep the problem simple and tractable. We use a simple, isothermal equation of state

$$\begin{equation} P = c^2 \Sigma, \end{equation} \tag{ 3 }$$

where c is the sound speed. We consider an inertial coordinate system in which the particle is always at the origin so the background flow velocity is given by $\mathbf {v}=v_0 \mathbf {e}_y$. For the gravitational potential, we make use of the form that is commonly incorporated in the investigation of the disk–planet interaction in two-dimensional analyses,

$$\begin{equation} \psi = - {\displaystyle \frac{{\it GM}}{\sqrt{x^2 + y^2 + \epsilon ^2}}}, \end{equation} \tag{ 4 }$$

where is the softening length, which is a sizable fraction of the thickness of the slab in order to mimic the three-dimensional effects.

This form of the gravitational potential (4) is the model for the vertically averaged gravitational potential; therefore, the softening length should be of the order of the vertical scale length of the disk. It is to be noted that, with this form of the potential, the gravitational potential close to the point mass (typically, the distance closer than ) is underestimated. If the vertical averaging is taken, the potential close to the mass behaves as ∝log r, where r is the distance from the point mass, in contrast to Equation (4), which behaves as a constant for r ≲ .

However, as we shall show later, the main contribution to the dynamical friction force comes from the region r ∼ , and the contribution from the region r < would affect the results by only some factor. It is also noted that the dependence of the force on the physical parameters can be captured with the form of the potential given by Equation (4). Therefore, we use this form of the potential as a model for the vertically averaged potential in this paper. The outcome of this approximation will be discussed in detail later in this section. We also note that this form of the potential is widely used in the numerical calculations of the disk–planet interaction. One benefit of using the simple form of the gravitational potential given by Equation (4) is that it is possible to obtain an analytic expression for the dynamical friction.

We now perform the linear perturbation analysis to calculate the surface density perturbation induced by the particle's gravity. We use the subscript zero to indicate the background quantities, and we denote all the perturbed quantities by δ, e.g., surface density is given by Σ = Σ₀ + δΣ. We retain the terms up to the first order of the perturbation. Equations (1) and (2) now become

$$\begin{equation} {\displaystyle \frac{\partial }{\partial t}} {\displaystyle \frac{\delta \Sigma }{\Sigma _0}} + v_0 {\displaystyle \frac{\partial }{\partial y}} {\displaystyle \frac{\delta \Sigma }{\Sigma _0}} + {\displaystyle \frac{\partial }{\partial x}} \delta v_x + {\displaystyle \frac{\partial }{\partial y}} \delta v_y = 0, \end{equation} \tag{ 5 }$$

$$\begin{equation} \left[ {\displaystyle \frac{\partial }{\partial t}} + v_0 {\displaystyle \frac{\partial }{\partial y}} - \nu \nabla ^2 \right] \delta v_x = -c^2 {\displaystyle \frac{\partial }{\partial x}} {\displaystyle \frac{\delta \Sigma }{\Sigma _0}} - {\displaystyle \frac{\partial \psi }{\partial x}}, \end{equation} \tag{ 6 }$$

$$\begin{equation} \left[ {\displaystyle \frac{\partial }{\partial t}} + v_0 {\displaystyle \frac{\partial }{\partial y}} - \nu \nabla ^2 \right] \delta v_y = -c^2 {\displaystyle \frac{\partial }{\partial y}} {\displaystyle \frac{\delta \Sigma }{\Sigma _0}} - {\displaystyle \frac{\partial \psi }{\partial y}}. \end{equation} \tag{ 7 }$$

From these equations, we derive a single second-order differential equation for the surface density perturbation,

$$\begin{equation} \left({\displaystyle \frac{\partial }{\partial t}} + v_0 {\displaystyle \frac{\partial }{\partial y}} \right)^2 \alpha - \nu \nabla ^2 \left({\displaystyle \frac{\partial }{\partial t}} + v_0 {\displaystyle \frac{\partial }{\partial y}} \right) \alpha - c^2 \nabla ^2 \alpha = \nabla ^2 \psi, \end{equation} \tag{ 8 }$$

where we have defined α ≡ δΣ/Σ₀, and ∇² = (∂/∂x)² + (∂/∂y)² is the Laplacian.

We now consider a steady state where ∂/∂t = 0. In Ostriker (1999), it is pointed out that it is necessary to perform a time-dependent analysis in order to correctly obtain the dynamical friction force, especially when the particle's velocity is subsonic. This is because the contribution to the force coming from a place very far from the particle is not negligible. However, as we shall show below, the analysis assuming the steady state is adequate in the slab geometry we consider in this section since the perturbation induced at a distant place from the particle does not contribute to the force.

In Appendix A, we explicitly show that the force coming from the time-dependent terms falls as t⁻¹. Here, we briefly show this by the order-of-magnitude estimate. If we consider the place far from the particle, the gravitational potential is given by

$$\begin{equation} \psi \sim {\displaystyle \frac{{\it GM}}{r}}, \end{equation} \tag{ 9 }$$

where r is the distance from the particle. We expect that this gravitational energy is of the same order of magnitude as the perturbed thermal energy of the gas. In the three-dimensional analysis, we therefore expect that

$$\begin{equation} {\displaystyle \frac{\delta \rho }{\rho _0}} \sim {\displaystyle \frac{{\it GM}}{c^2 r}}, \end{equation} \tag{ 10 }$$

where ρ is the gas density. In the slab geometry, we expect that

$$\begin{equation} {\displaystyle \frac{\delta \Sigma }{\Sigma _0}} \sim {\displaystyle \frac{{\it GM}}{c^2 r}}, \end{equation} \tag{ 11 }$$

where Σ is the surface density. In the three-dimensional analysis, the force δF acting on the particle from the gas shell at the distance r with the width δr is

$$\begin{equation} \delta F \sim {\displaystyle \frac{{\it GM} \delta \rho r^2 \delta r}{r^2}} \sim \frac{({\it GM})^2 \rho \delta r}{c^2 r} \propto {\displaystyle \frac{\delta r}{r}}. \end{equation} \tag{ 12 }$$

The force from the shell decays only with r⁻¹. Therefore, when summed over all the shells, the contribution from the distant shell is not negligible.³ In the case of two-dimensional slab geometry, on the other hand, the force from the ring at distance r with width δr is

$$\begin{equation} \delta F \sim {\displaystyle \frac{{\it GM} \delta \Sigma r \delta r}{r^2}} \sim {\displaystyle \frac{({\it GM})^2 \Sigma _0 \delta r}{c^2 r^2}} \propto {\displaystyle \frac{\delta r}{r^2}}. \end{equation} \tag{ 13 }$$

The force from the distant ring decays as r⁻²; therefore, the contribution from the distant ring does not account for the total force when summed over all the rings. We note here that the contribution to the dynamical friction force far from the planet decays as r⁻¹ if we integrate all the rings. This explains why the force decays as t⁻¹ in the subsonic case if we consider the slab geometry. The rings that contribute to the dynamical friction force reside at r ∼ ct, as pointed out by Ostriker (1999).

In this paper, we shall show in detail the derivation of the dynamical friction force exerted on the particle embedded in a gaseous slab using the two-dimensional approximation. Before going on to the details, we show that the dependences of the force on the physical parameters of the gas can be derived through the two-dimensional approximation.

In the realistic case of the gaseous slab with the thickness of L_z(∼ ), the interaction between the gas and the particle can be well approximated by the two-dimensional analysis for the scales larger than L_z (far region). Inside this region (near region), the interaction should be calculated in three dimensions. Integrating Equation (12) between the minimum cutoff scale r_min and γL_z (γ is a factor of the order of unity), the contribution to the dynamical friction force from the near region, F_N, is given by

$$\begin{equation} F_N \sim {\displaystyle \frac{({\it GM})^2 \rho }{c^2}} \log \left({\displaystyle \frac{\gamma L_z}{r_{\rm min}}} \right). \end{equation} \tag{ 14 }$$

The contribution to the force from the far region, F_F, is given by integrating Equation (13) from γL_z to infinity and

$$\begin{equation} F_F \sim {\displaystyle \frac{({\it GM})^2 \Sigma }{c^2 \gamma L_z}} \sim {\displaystyle \frac{({\it GM})^2 \rho }{\gamma c^2}}, \end{equation} \tag{ 15 }$$

where we have used Σ ∼ ρL_z. The integral from the far region does not diverge at r = ∞ since there is not enough mass to contribute to the force in the slab geometry in the far region, in contrast to the case of the homogeneous three-dimensional distribution of the gas. It should be noted that the dependences of the force on physical parameters in F_F and F_N are the same except for the logarithmic factor, which only introduces very weak dependences on L_z and r_min. By choosing the appropriate value of γ, it is possible to have an expression that would reproduce the results of three-dimensional calculations from two-dimensional calculations. Therefore, we expect that the two-dimensional analyses can be used to understand the fundamental aspects of the dynamical friction force. More detailed discussions on the two-dimensional approximation in relation to the application to the disk–planet interaction can be found in Section 2.5.

We now derive the steady-state solution of Equation (8). We make use of the Fourier transform to solve the equations, but our goal is to derive the expression of the force exerted on the particle, which is, of course, the quantity in the real space. For any perturbed quantity f (x, y), we define the Fourier transform by

$$\begin{equation} f(x,y) = {\displaystyle \frac{1}{2\pi }} \int dk_x dk_y \tilde{f}(k_x,k_y) e^{i(k_x x + k_y y)}, \end{equation} \tag{ 16 }$$

and the inverse transform is

$$\begin{equation} \tilde{f}(k_x,k_y) = {\displaystyle \frac{1}{2\pi }} \int dx dy f(x,y) e^{-i(k_x x + k_y y)}. \end{equation} \tag{ 17 }$$

The steady state of surface density perturbation can be derived from Equation (8) in the Fourier space as

$$\begin{equation} \tilde{\alpha } = - {\displaystyle \frac{k^2 \tilde{\psi }}{c^2 k^2 - v_0^2 k_y^2 + i \nu v_0 k_y k^2}}, \end{equation} \tag{ 18 }$$

where k² = k²_x + k_y².

If we know the profile of the perturbed surface density, the dynamical friction force exerted on the particle is calculated by

$$\begin{equation} F_y(x,y) = \int d^2x \delta \Sigma (x,y) {\displaystyle \frac{\partial \psi }{\partial y}}. \end{equation} \tag{ 19 }$$

We note that the x-component of the force is zero because of the symmetry. From the solution in the Fourier space, the dynamical friction force is given by

$$\begin{equation} F_y(x,y) = 2 \int _0^{\infty } dk_y \int _{-\infty }^{\infty } dk_x k_y \tilde{\psi } \mathrm{Im} \left(\tilde{\alpha } \right), \end{equation} \tag{ 20 }$$

where the Fourier transform of the gravitational potential $\tilde{\psi }$ is given by

$$\begin{equation} \tilde{\psi } = - {\displaystyle \frac{{\it GM}}{k}} e^{-k \epsilon }. \end{equation} \tag{ 21 }$$

We comment that in Equation (20), we integrate the contribution to the force coming from all over the space, including the region close to the particle. The contribution to the force coming from the near region is effectively cut off due to the softening. This is actually a reasonable treatment, given the uncertainty of the softening length , and it is possible to derive the dependences of the force on physical parameters. See Section 2.5 for more detailed discussions.

We calculate the dynamical friction force using Equations (18) and (20). By changing the integration variables first from (k_x, k_y) to (k, θ) via

$$\begin{equation} k_x = k \cos \theta \end{equation} \tag{ 22 }$$

and

$$\begin{equation} k_y = k \sin \theta, \end{equation} \tag{ 23 }$$

and then from k to u = νv₀k/c², we can rewrite the integral as

$$\begin{equation} F_y = {\displaystyle \frac{\Sigma _0 ({\it GM})^2}{c^2 \epsilon }} {\displaystyle \frac{\Gamma }{2}} \int _0^{\infty } u I_{\theta } \big(M_0^2, u^2\big) e^{ -\Gamma u } du, \end{equation} \tag{ 24 }$$

where M₀ ≡ v₀/c is the Mach number of the background flow and Γ ≡ 2c/M₀ν. The function I_θ(p, q) that appears in the integral is given by

$$\begin{equation} I_{\theta } (p,q) = \int _{0}^{2\pi } {\displaystyle \frac{\sin ^2\theta }{(1-p\sin ^2\theta)^2 + q \sin ^2\theta }} d\theta. \end{equation} \tag{ 25 }$$

It is actually possible to perform the integration involved in I_θ(p, q) analytically. The explicit form of this is found in Appendix B.

It is now possible to perform the integral (24) numerically to obtain the dynamical friction force in a gaseous slab if we give two dimensionless parameters M₀ and Γ. However, for a moment, we look at some limits to obtain analytic expressions in order to investigate how the dynamical friction force depends on these parameters.

2.3.1. Subsonic Limit

We first look at the subsonic case where M₀ ≪ 1. In this case, we approximate 1 − M²₀sin θ ∼ 1 so the integral (25) becomes

$$\begin{equation} I_{\theta }\big(M_0^2,u^2\big) \sim \int _0^{2\pi } d\theta {\displaystyle \frac{\sin ^2\theta }{1+u^2 \sin \theta }}. \end{equation} \tag{ 26 }$$

Therefore, in the limit of u → 0, I_θ(M²₀, u²) ∼ π while in the limit of u → ∞, I_θ(M²₀, u²) ∼ 2π/u². Connecting these two limits, we approximate the integral by

$$\begin{eqnarray} I_{\theta } \big(M_0^2,u^2\big) = \left\lbrace \begin{array}{@{}cc}\pi & (u<\sqrt{2}) \\ 2\pi /u^2 & (u>\sqrt{2}) \end{array} \right.. \end{eqnarray} \tag{ 27 }$$

The dynamical friction force is obtained from Equation (24). It is possible to perform the integral analytically and we have

$$\begin{equation} F_y = {\displaystyle \frac{\pi }{2}} {\displaystyle \frac{\Sigma _0 ({\it GM})^2}{\epsilon c^2}} \Gamma \left[ {\displaystyle \frac{1-(1+\sqrt{2}\Gamma)e^{-\sqrt{2}\Gamma }}{\Gamma ^2}} - 2 \mathrm{Ei}(-\sqrt{2}\Gamma) \right], \end{equation} \tag{ 28 }$$

where Ei(x) is the exponential integral defined by

$$\begin{equation} \mathrm{Ei}(-x) = -\int _{x}^{\infty } dt {\displaystyle \frac{e^{-t}}{t}}. \end{equation} \tag{ 29 }$$

In the limit of subsonic perturbation, we expect that Γ ≫ 1. We then obtain

$$\begin{equation} F_y \sim {\displaystyle \frac{\pi }{4}} {\displaystyle \frac{\Sigma _0 (GM)^2}{\epsilon c^2}} {\displaystyle \frac{M_0 \nu }{c \epsilon }}. \end{equation} \tag{ 30 }$$

The dynamical friction force is proportional to the speed of the particle and also the viscosity. In the case of an inviscid, steady state with a subsonic particle embedded in the gaseous slab, the dynamical friction force vanishes as indicated by the time-dependent analysis shown in Appendix A.

2.3.2. Supersonic Limit

We now consider the supersonic limit where M₀ ≫ 1. We first consider the behavior of the integral I_θ(M²₀, u²) separately in the two limits, u ≪ M²₀ and u ≫ M²₀.

In the case of u ≪ M²₀, the dominant contribution to the integral comes from θ ∼ θ₀, where θ₀ is given by

$$\begin{equation} \sin ^2{\theta _0}={\displaystyle \frac{1}{M_0^2}}. \end{equation} \tag{ 31 }$$

Approximating sin θ = sin (θ₀ + δθ) to the first order of δθ,

$$\begin{equation} \sin (\theta _0+\delta \theta) \sim \sin \theta _0 + \delta \theta \cos \theta _0, \end{equation} \tag{ 32 }$$

the integral can be approximated by

$$\begin{equation} I_{\theta }\big(M_0^2,u^2\big) \sim {\displaystyle \frac{4}{M_0^2}} \int _{-\infty }^{\infty } d \delta \theta {\displaystyle \frac{1}{4\big(M_0^2-1\big)\delta \theta ^2 + (u/M_0)^2}}, \end{equation} \tag{ 33 }$$

where the factor of four comes from the fact that there are four positions within 0 < θ < 2π where sin ²θ equals 1/M²₀, and the range of integration is extended from −∞ to ∞. The integrand of Equation (33) is the Lorentzian function with width ∼u/M²₀. The condition where this width must be very small compared to the unity (or more exactly 2π, in order to justify the extension of the range of the integration) leads to the condition u ≪ M²₀. Using the formula

$$\begin{equation} \int _{-\infty }^{\infty } dx {\displaystyle \frac{1}{1+x^2}} = \pi, \end{equation} \tag{ 34 }$$

we obtain

$$\begin{equation} I_{\theta }\big(M_0^2,u^2\big) \sim {\displaystyle \frac{2\pi }{uM_0\sqrt{M_0^2-1}}}. \end{equation} \tag{ 35 }$$

In the case of u ≫ M²₀, the width of the Lorentzian is too wide to extend the integration range from (0, 2π) to (− ∞, ∞). In this case, we expect that u²sin ²θ ≫ (1 − M²₀sin ²θ)² is satisfied in the most part of θ. We also approximate 1 − M²₀sin ²θ ∼ −M₀²sin ²θ so we can rewrite the integral as

$$\begin{equation} I_{\theta }\big(M_0^2,u^2\big) \sim \int _0^{2\pi } {\displaystyle \frac{1}{M_0^2 \sin ^2\theta + u^2}}. \end{equation} \tag{ 36 }$$

Using the formula

$$\begin{equation} \int _0^{2\pi } d\theta {\displaystyle \frac{1}{a^2 + \sin ^2 \theta }} = {\displaystyle \frac{2\pi }{a \sqrt{a^2+1}}}, \end{equation} \tag{ 37 }$$

we obtain

$$\begin{equation} I_{\theta }\big(M_0^2,u^2\big) \sim {\displaystyle \frac{2\pi }{u^2}}, \end{equation} \tag{ 38 }$$

where we have used the condition u ≫ M²₀. In summary, in the supersonic limit, the integral (25) can be approximated by

$$\begin{eqnarray} I_{\theta }\big(M_0^2,u^2\big) = \left\lbrace \begin{array}{@{}cc}{\displaystyle \frac{2\pi }{uM_0\sqrt{M_0^2-1}}} & \big(u<M_0\sqrt{M_0^2-1}\big) \\ {\displaystyle \frac{2\pi }{u^2}} & \big(u>M_0\sqrt{M_0^2-1}\big) \end{array} \right.. \end{eqnarray} \tag{ 39 }$$

We can now perform the integral (24) to give the dynamical friction force. The result is

$$\begin{equation} F_y = \pi {\displaystyle \frac{\Sigma _0 ({\it GM})^2}{\epsilon c^2}} \Gamma \left[ {\displaystyle \frac{1-\exp \left(-\Gamma M_0 \sqrt{M_0^2-1}\right)}{\Gamma M_0 \sqrt{M_0^2-1}}} - \mathrm{Ei}\left(-\Gamma M_0 \sqrt{M_0^2-1} \right) \right]. \end{equation} \tag{ 40 }$$

In the supersonic limit, we expect that ΓM²₀ ≫ 1. In this case the expression simplifies to

$$\begin{equation} F_y \sim \pi {\displaystyle \frac{\Sigma _0 ({\it GM})^2}{\epsilon c^2 M_0^2}}. \end{equation} \tag{ 41 }$$

As expected, we obtain that the force is proportional to v⁻²₀.

We note that the dynamical friction force is independent of the viscosity in the supersonic limit. This indicates that the dynamical friction force in the supersonic limit is insensitive to the dissipation process. We conjecture that the dynamical friction force is insensitive to the physical state of the slab (whether the slab is isothermal or radiative) in the case of the supersonic motion. A similar situation happens in the discussion of the Lindblad torque exerted on a planet embedded in a disk (Meyer-Vernet & Sicardy 1987).

We now come to the point where we consider all the values of M₀. We integrate Equation (24) numerically to find the dependence of the dynamical friction force on the physical parameters. As we have seen in the above analytic discussions, there are two dimensionless parameters, Γ and M₀, in the problem at hand. In this section, we use the "Reynolds number," Re, defined by

$$\begin{equation} {\rm Re} = {\displaystyle \frac{c\epsilon }{\nu }} \end{equation} \tag{ 42 }$$

instead of Γ, since the parameter Γ contains both the Mach number and the viscous coefficient.

Figure 1 shows the dependence of the dynamical friction force normalized by Σ₀(GM)²/c² as a function of the Mach number and the Reynolds number. In the subsonic regime, we see the expected behavior from Equation (30) where dynamical friction force decreases as we decrease the Mach number and viscosity. For the supersonic case, we also see the behavior expected from Equation (41) where the force decreases as we increase the velocity but the force only weakly depends on the values of viscosity.

Zoom In Zoom Out Reset image size

Figure 2 shows the dependence of the dynamical friction force on the Mach number in the case of Re = 100, and we also show the formulae in the subsonic and supersonic limits, Equations (30) and (41). These limiting formulae well describe the behaviors of the dynamical friction force. Equation (41) is actually a good approximation of the dynamical friction force even in the case of M₀ ∼ 2.

Zoom In Zoom Out Reset image size

We briefly comment on the divergence at M₀ = 1. In Figures 1 and 2, we see that the dynamical friction force diverges at M₀ = 1. This divergence comes from the matching between the background velocity and the sound speed, and was already seen in the previous linear analyses by Ostriker (1999). It is not possible to avoid this divergence by the effect of viscosity. We consider that nonlinear effects are important here.

We have derived the dynamical friction force exerted on a particle embedded in a homogeneous gas slab using the two-dimensional approximation. We now discuss the validity of this approximation in the subsonic and supersonic cases.

As we have discussed before, the two-dimensional treatment is based on the averaging of the equations in the vertical direction, and the phenomena that occur on the scales larger than the vertical averaging scale can be well approximated by the two-dimensional calculations. In our model, the vertical averaging scale is given by , which is the softening scale of the gravitational potential in Equation (4).

In the subsonic case, the dominant contribution to the integral (24) comes from $\Gamma u \sim \mathcal {O}(1)$, which is $\epsilon k \sim \mathcal {O}(1)$. This indicates that in the subsonic case, the perturbation with the scale comparable to the gravitational softening length becomes important in determining the dynamical friction force. Therefore, it is indicated that the full three-dimensional treatment may be necessary to have more quantitative results.

In the case of the supersonic motion, let us consider the case in which ΓM²₀ = 2v₀/ν ≫ 1 for simplicity. In this case, the integral (24) for the dynamical friction force is approximated by

$$\begin{equation} F_y \propto \int _{0}^{\infty } du {\displaystyle \frac{e^{-\Gamma u}}{M_0^2}} \propto \int _{0}^{\infty } dk e^{-\epsilon k}; \end{equation} \tag{ 43 }$$

therefore, all the scales satisfying Γu ∼ k ≲ 1 contribute to the force. In other words, this includes all the scales from the large scale (small k), where we expect that the two-dimensional approximation is valid, down to the cutoff scale contribute to the force. Since the dependence of the cutoff scale is also present in the case of the supersonic motion, rigorous three-dimensional treatment is necessary to give more quantitative expressions of the dynamical friction force.

However, we expect that it is possible to capture the dependences on physical parameters even in the two-dimensional approximation. In the beginning of Section 2.2, we have estimated the order of the magnitude of the dynamical friction force in Equations (12) and (13), and the contribution from the near region (14) and that from the far region (15) differ by only a logarithmic factor log (L_z/r_min). The dynamical friction force derived in this paper corresponds to the contribution from the far region. The lack of the logarithmic factor is due to the fact that the potential is smoothed in the region r ≲ (see also the discussion after Equation (4)). We note the similarity of the expressions derived by the simple order-of-magnitude estimate (15) and more rigorous calculations (41).

In order to estimate the contribution from the near region, it is necessary to determine the appropriate value of r_min. However, it is not straightforward; there are several publications on this using the nonlinear calculations (e.g., Sánchez-Salcedo & Brandenburg 1999; Kim & Kim 2009). Let us consider, for a moment, the case of the planet embedded in a protoplanetary disk. The naive estimate of r_min is the radius of the planet, and in the case of an Earth-mass planet, it is of the order of 10⁴ km. If we use the typical disk scale height H ∼ 0.05 AU, the value of log (H/r_min) amounts to ∼7. Recent nonlinear calculations (Kim & Kim 2009) suggest that the deviation from the linear value may be smaller due to the shock formation. If we use r_min ∼ GM/v², which is the typical scale of the distance between the particle and the shock in supersonic motion (e.g., Kim & Kim 2009), instead of the planet radius for r_min, the logarithmic factor is ∼3–6 for an Earth-mass body. We shall present further discussions on the nonlinear effects on dynamical friction at the end of Section 3.⁴ In any case, the force derived in this section may be different by some factor due to the two-dimensional approximations, but such a logarithmic factor can be approximately taken into account by choosing the appropriate value for (or γ in Equation (15)).

Despite the uncertainty of the two-dimensional approximation, especially for the value of that should be taken, we still apply this result to the problem of the disk–planet interaction. It is because many of the calculations to date are done in two-dimensional approximations, and they involve the potential of the form given in Equation (4). One goal of this paper is to investigate how well the simple model using the dynamical friction force may be able to describe the numerical results of the disk–planet interaction. More rigorous treatment of the dynamical friction force including the three-dimensional effects and the vertical stratification will be discussed in future publications. It is also necessary to have three-dimensional nonlinear calculations of the disk–planet interaction to compare the model. We note that the discussion on the form of the potential applies not only to the linear perturbation analyses, but also to the nonlinear calculations.

We consider a protoplanetary disk with a central star with mass M_* and surface density profile with r^−p,

$$\begin{equation} \Sigma = \Sigma _0 \left({\displaystyle \frac{r}{r_0}} \right)^{-p}, \end{equation} \tag{ 51 }$$

where Σ₀ is the surface density at r = r₀. We assume for simplicity that the disk is locally isothermal with the temperature profile T∝r^−q. Then, the sound speed c of the disk gas varies as c∝r^−q/2. We assume that the disk is rotating at the Kepler angular frequency and neglect the small difference arising from the pressure gradient. We define the scale height of the disk H by H = c/Ω_K, where c is the sound speed and Ω_K is the Kepler angular frequency. The scale height H therefore varies as H∝r^{(3 − q)/2}. The disk aspect ratio, H/r, is given by

$$\begin{equation} {\displaystyle \frac{H}{r}} = h_0 \left({\displaystyle \frac{r}{r_0}} \right)^{(1-q)/2}, \end{equation} \tag{ 52 }$$

where h₀ is the disk aspect ratio at r = r₀. We consider a planet with mass M_p with semimajor axis a and eccentricity e.

In the later sections, we shall often refer to the "fiducial model," which we define as p = 3/2, q = 1, M_* = M_☉, M_p = M_⊕, h₀ = 0.05, and Σ₀ is chosen in such a way that the mass contained within 5 AU is 2M_J. In this model, the disk aspect ratio is constant throughout the disk. We call this model "fiducial" because it is the model used in the numerical simulations by Cresswell & Nelson (2006); we shall compare our analytic results with their numerical calculations. Note that this is not exactly the same as Minimum Mass Solar Nebula.

The gravitational potential of the planet is given by Equation (4), where (x, y) is now the coordinate centered on the planet. The softening length is given by = ₀H(r) where ₀ is the dimensionless parameter. Cresswell & Nelson (2006) used = 0.5. We note that in our model, the softening parameter varies as the location of the disk if H varies as r.

We first calculate the relative velocity between the planet and the disk. It is readily calculated if we give the planet's semimajor axis and eccentricity, since we assume the disk is in Keplerian rotation.

We assume that the planet mass M_p is negligible compared to the mass of the central star M_* and use the cylindrical coordinate system (r, ϕ, z) with the origin at the central star. We denote the mean motion of the planet by n = (GM_*/a³)^1/2, and the true anomaly by f. The velocity of the planet is then

$$\begin{equation} \mathbf {v}_{\rm p} = {\displaystyle \frac{aen}{\eta }} \sin f \mathbf {e}_r + {\displaystyle \frac{an}{\eta }} \left(e\cos f +1 \right)\mathbf {e}_{\phi }, \end{equation} \tag{ 53 }$$

where $\mathbf {e}_r$ and $\mathbf {e}_{\phi }$ are the unit vectors in r- and ϕ-directions, respectively. The velocity of the gas is

$$\begin{equation} \mathbf {v}_{g} = \sqrt{{\displaystyle \frac{{\it GM}_{\ast }}{r}} } \mathbf {e}_{\phi } = \sqrt{{\displaystyle \frac{n^2 a^3}{r}} } \mathbf {e}_{\phi }. \end{equation} \tag{ 54 }$$

The relative velocity between the gas and the planet is calculated by $\Delta \mathbf {v}= \mathbf {v}_g-\mathbf {v}_{\rm p}$. Now, we define the unit vectors $(\mathbf {e}_{T},\mathbf {e}_N)$, where $\mathbf {e}_T$ is directed toward the velocity of the planet and $\mathbf {e}_N$ is perpendicular to $\mathbf {e}_T$ and in the orbital plane (see Figure 3). Unit vectors $(\mathbf {e}_r,\mathbf {e}_{\phi })$ and $(\mathbf {e}_T,\mathbf {e}_N)$ are related by

$$\begin{equation} \mathbf {e}_r = \mathbf {e}_T \cos \beta - \mathbf {e}_N \sin \beta \end{equation} \tag{ 55 }$$

$$\begin{equation} \mathbf {e}_{\phi } = \mathbf {e}_T \sin \beta + \mathbf {e}_N \cos \beta, \end{equation} \tag{ 56 }$$

where

$$\begin{equation} \cos \beta = {\displaystyle \frac{e \sin f}{\sqrt{1+2e\cos f + e^2}}} \end{equation} \tag{ 57 }$$

and

$$\begin{equation} \sin \beta = {\displaystyle \frac{1 + e \cos f}{\sqrt{1+2e\cos f + e^2}}}. \end{equation} \tag{ 58 }$$

Zoom In Zoom Out Reset image size

In Figure 4, we show the evolution of the relative velocity vector $\Delta \mathbf {v}$ over one orbital period with e = 0.5. The planet feels headwind at the perihelion, and tailwind at the aphelion. Figure 5 shows the amplitude of the relative velocity over one orbital period. It is noted that if the planet's eccentricity exceeds the disk aspect ratio, the relative velocity is always supersonic, regardless of the position in the orbit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We make use of Gauss's equations to calculate the orbital evolution of the planet. For the evolution of the semimajor axis and eccentricity, Gauss's equations read

$$\begin{equation} {\displaystyle \frac{da}{dt}} = {\displaystyle \frac{2 v}{a n^2}} T \end{equation} \tag{ 59 }$$

and

$$\begin{equation} {\displaystyle \frac{de}{dt}} = {\displaystyle \frac{1}{v}} \left[ 2(e + \cos f)T - {\displaystyle \frac{r}{a}} \sin f N \right], \end{equation} \tag{ 60 }$$

where v is the speed of the planet and r is the distance between the planet and the central star. We denote the perturbing force per unit mass by (T, N), where T is the component of the force (per unit mass) in the tangential direction of the orbit (positive direction is the direction of the velocity), and N is the component in the plane of the orbit and perpendicular to T (positive direction is directed inside the orbital ellipse). The evolution equation for the true anomaly is given by

$$\begin{equation} {\displaystyle \frac{df}{dt}} = {\displaystyle \frac{a^2 n \eta }{r^2}} + {\displaystyle \frac{\eta }{nae}} \left[ R \cos f - S \left(1+{\displaystyle \frac{a(1-e^2)}{r}} \sin f \right) \right], \end{equation} \tag{ 61 }$$

where η = (1 − e²)^1/2, R is the component of the force in the radial direction, and S is the component perpendicular to R. Components (T, W) and (R, S) are related by

$$\begin{equation} R = T \cos \beta - N \sin \beta \end{equation} \tag{ 62 }$$

and

$$\begin{equation} S = T \sin \beta + N \cos \beta. \end{equation} \tag{ 63 }$$

For the expressions of the perturbing force (T, W), we use the dynamical friction formula. In this paper, we especially focus on the planet with high eccentricity, and therefore assume that the relative speed between the gas and the planet is always supersonic. Therefore, we use the limiting form of a highly supersonic case, Equation (41). The amplitude of the force is

$$\begin{equation} F = {\displaystyle \frac{\pi \Sigma (r) G^2M_{\rm p}}{\epsilon \Delta v^2}}, \end{equation} \tag{ 64 }$$

where Δv is the relative speed between the gas and the planet. The direction of the force is opposite to the relative velocity vector.

In order to find the long-term evolution of orbital parameters, we average Equations (59)–(61) over one orbital period assuming the planet is in a fixed orbit. Since the variation of a and e is very small over one orbital period, the assumption of a fixed orbit is justified. The averaging is done in such a way that

$$\begin{equation} \overline{{\displaystyle \frac{da}{dt}} } = {\displaystyle \frac{1}{T_{\rm orb}}} \int _0^{T_{\rm orb}} dt {\displaystyle \frac{da}{dt}} = {\displaystyle \frac{1}{T_{\rm orb}}} \int _0^{2\pi } df {\displaystyle \frac{da/dt}{df/dt}}, \end{equation} \tag{ 65 }$$

where T_orb is the time taken over one period. The averaging is done in the same way for eccentricity.

In this section, we compare our results with previous numerical calculations for the fiducial model.

We first show the torque and power exerted on the planet by the disk over one orbit. The power $\mathcal {P}$ exerted on the planet is given by

$$\begin{equation} \mathcal {P} = \mathbf {v}_{\rm p} \cdot \mathbf {F}_{\rm disk} = vT \end{equation} \tag{ 66 }$$

and the torque $\mathcal {T}$ is

$$\begin{equation} \mathcal {T} = \mathbf {e}_z \cdot \left(\mathbf {r}_{\rm p} \times \mathbf {F}_{\rm disk}\right) = rS, \end{equation} \tag{ 67 }$$

where $\mathbf {r}_{\rm p}$, $\mathbf {v}_{\rm p}$, and $\mathbf {F}_{\rm disk}$ are the position vector of the planet from the disk, the velocity of the planet, and the force exerted on the planet, respectively. In Figure 6, we show $\mathcal {P}$ and $\mathcal {T}$ normalized by the planet's energy (the sign inverted),

$$\begin{equation} E_{\rm p}= {\displaystyle \frac{1}{2}} {\displaystyle \frac{{\it GM}_{\ast }}{a}}, \end{equation} \tag{ 68 }$$

and the angular momentum,

$$\begin{equation} L_{\rm p}= \sqrt{{\it GM})_{\ast }a (1-e^2)}, \end{equation} \tag{ 69 }$$

respectively.

Zoom In Zoom Out Reset image size

In the vicinity of the perihelion, the torque and the power are negative because of the headwind toward the planet, and they are positive in the vicinity of the aphelion owing to the tailwind. Similar behavior is observed in the numerical simulations of Cresswell et al. (2007). The behavior of the torque exerted on the planet looks very similar except for the small difference between the location of the perihelion or aphelion and the location of the maximum or minimum of the torque (see Figure 8 of Cresswell et al. 2007). In our model, this small difference in phase does not appear since we assume that the force is proportional to the relative velocity at the location of the planet.

The comparison of our results with Figure 10 of Cresswell et al. (2007) shows that the behavior of the power is very different. This difference cannot be attributed to the sets of parameters they have used different from our fiducial model. This may be because our model for the force exerted on the planet is very simplified. This difference in power leads to the different behavior in da/dt and de/dt over one orbital period, which is shown in Figure 7. We note here that the evolution of eccentricity is qualitatively different from Figure 9 of Cresswell et al. (2007). Fortunately, however, it is possible to obtain quantitatively the same results with previous numerical calculations when we take the average over one orbital period, as shown below.

Zoom In Zoom Out Reset image size

Figure 8 shows the results of migration rate t_m and eccentricity damping rate t_e, which are obtained by averaging over one orbital period. Also plotted in this figure is the formula obtained by Papaloizou & Larwood (2000) given in Equations (48) and (49). For the migration rate, we also show t_m given by Papaloizou & Larwood (2000) times 3/2 for e < 0.5, which describes the results of numerical simulations better as noted by Cresswell & Nelson (2006). For the damping of eccentricity, the formula by Papaloizou & Larwood (2000) fits the results of numerical simulations well. We see that our results for t_e agree very well with the formula by Papaloizou & Larwood (2000) for e ≲ 0.7. For the migration rate, our formulation consistently results in a timescale slower by a factor of two compared to the numerical simulations. Our results deviate from those of Papaloizou & Larwood (2000) for large values of eccentricity. We expect that our treatment is actually better for highly eccentric cases (discussed in later sections), and that this difference is attributed to the breakdown of the expansion of e used by Papaloizou & Larwood (2000).

Zoom In Zoom Out Reset image size

We have seen that our formulation results in a consistent timescale of the evolution of orbital parameters for the fiducial model. We now look at how the timescales behave as we vary the disk model.

In our framework, the force exerted on the planet can be written in the form

$$\begin{equation} T = K r^{-\alpha } {\displaystyle \frac{\Delta v_T}{\Delta v^3}} \end{equation} \tag{ 70 }$$

and

$$\begin{equation} N = K r^{-\alpha } {\displaystyle \frac{\Delta v_N}{\Delta v^3}}, \end{equation} \tag{ 71 }$$

where K and α are constants. For the model of the force we use in this paper, which is Equation (64), if we use the softening parameter proportional to the disk scale height, = ₀H(r), the constant K is given by

$$\begin{equation} K = {\displaystyle \frac{\pi G^2 M_{\rm p}\Sigma _0}{\epsilon _0 H_0 r_0^{-\alpha }}}, \end{equation} \tag{ 72 }$$

and α is given by

$$\begin{equation} \alpha = p + {\displaystyle \frac{3-q}{2}}, \end{equation} \tag{ 73 }$$

where H₀ is the disk scale height at r = r₀. Therefore, models with the same K and α yield the same results. In this section, we show the results of how t_a and t_e depend on the disk profile p, or in more general terms, α. We fix the value of surface density to be Σ₀/(M_*/r²₀) = 10⁻⁴, the planet mass to be M_p/M_* = 10⁻⁶, the disk aspect ratio to be H₀/r₀ = 0.05, and the softening parameter to be ₀ = 0.5. The dependence of these parameters on the timescales t_a and t_e is t∝M⁻¹_pΣ₀⁻¹₀H₀, as expected from the form of K, as long as the force exerted by the disk is sufficiently small compared to that from the central star. We use the disk model with q = 1, which gives α = p + 1, but we note again that the value of α is the key parameter that determines the values of t_a and t_e.

Figure 9 shows t_a and t_e for α = 1, 2, 3, 4 (p = 0, 1, 2, 3, respectively, for q = 1). In Figure 9, we fix a/r₀ = 1 and see how the timescale behaves as we vary the eccentricity. We see that the behaviors of t_a and t_e strongly depend on the imposed disk model. The timescale of the evolution of the semimajor axis increases steadily as we increase eccentricity if α ⩽ 2. If, on the other hand, α ⩾ 2, the values of t_a first increase as we increase e, but decrease at large e. The timescale of the evolution of eccentricity always increases as we increase eccentricity if α ⩽ 3, but there is a maximum of t_e at a certain value of e if α > 3.

Zoom In Zoom Out Reset image size

As we decrease the values of α, another interesting behavior occurs for t_a. In Figure 10, we show the behavior of t_a and t_e for α = −1, 0, 1 (p = −2, −1, 0). We see that the timescale of the eccentricity evolution does not vary much within this parameter range, while t_a is negative for almost all the values of e for α = 0 and −1. The negative sign of t_a indicates that the direction of the semimajor axis evolution is outward.

Zoom In Zoom Out Reset image size

The reason why we obtain the outward semimajor axis evolution can be qualitatively explained as follows. If the index p is negative, the surface density increases as a function of radius. As we have seen before, the planet feels the tailwind at the aphelion and therefore is exerted positive torque by the gas, while it feels headwind at the perihelion and is exerted negative torque. Since the surface density is larger at the aphelion than at the perihelion when p is negative, the planet feels overall positive torque, and therefore the semimajor axis increases.

The disk with negative p seems rather unrealistic. However, the surface density may increase locally as a function of radius. The inner edge of the disk is one example in which such local increase of the surface density is expected. If a planet with a finite eccentricity is placed in such a place, we expect that the semimajor axis can increase as a result of the disk–planet interaction. Recently, Ogihara et al. (2010) have suggested that the planets trapped in mean motion resonances can be stopped at the disk inner edge. This is because at the gap edge, the planet feels the disk only at the place close to the aphelion; therefore, it is exerted the positive torque by the disk. This positive torque is balanced by the negative torque exerted by the other planet in a mean motion resonance.

We have seen how the timescales of the semimajor axis evolution and the eccentricity damping vary as we change the disk parameter. Specifically, we have seen a strong dependence of t_a and t_e on the disk parameter α = p + (3 − q)/2 if the planet eccentricity is high. In this section, we analytically investigate the behaviors of t_a and t_e in the limit of e → 1.

In this section, we approximate df/dt (see Equation (61)) by

$$\begin{equation} {\displaystyle \frac{df}{dt}} = {\displaystyle \frac{a^2 n\eta }{r^2}}, \end{equation} \tag{ 74 }$$

since the contribution from the perturbing force is small compared to the force from the central star. We then obtain

$$\begin{eqnarray} {\displaystyle \frac{da}{df}} &=& {- }K {\displaystyle \frac{2(1-e^2)^{2-\alpha }}{a^{2+\alpha }n^4}} {\displaystyle \frac{e^2 \sin^2\!f + (1+e\cos\!f)^{3/2} (\sqrt{1+e \cos\!f}-1)}{(1+e\cos\!f)^{2-\alpha } \left[ e^2 \sin ^2\!f + (1+\cos \!f)(2+e\cos \!f - 2\sqrt{1+e\cos\!f}) \right]^{3/2}}} \end{eqnarray} \tag{ 75 }$$

and

$$\begin{eqnarray} {\displaystyle \frac{de}{df}} &=& {-}K {\displaystyle \frac{(1-e^2)^{3-\alpha }}{a^{3+\alpha }n^4}} {\displaystyle \frac{1}{(1+e^2+2e\cos \!f)(1+e\cos \!f)^{5/2-\alpha }}} \Bigg \lbrace {\displaystyle \frac{\sin ^2 \!f \left[ 2e^2 (e+\cos \!f)\sqrt{1+e\cos \!f}+e(1-e^2) \right]}{\left[ e^2 \sin ^2\! f + (1+\cos\! f)(2+e\cos \!f - 2\sqrt{1+e\cos \!f}) \right]^{3/2}}} \nonumber \\ && + {\displaystyle \frac{2 (e+\cos \!f)(1+e \cos \!f)^2 (\sqrt{2+e \cos \!f}-1)}{\left[ e^2 \sin ^2 \!f + (1+\cos \!f)(2+e\cos \!f - 2\sqrt{1+e\cos \!f}) \right]^{3/2}}} \Bigg \rbrace. \end{eqnarray} \tag{ 76 }$$

From Figure 7, we see that the force exerted when the planet is near the aphelion and the perihelion is important in determining the orbital evolution of the planet. We therefore look at the places close to the perihelion or aphelion. Around these points, cos f and sin f are approximated by

$$\begin{equation} \cos \!f \sim \pm \left(1 + {\displaystyle \frac{1}{2}} \delta f^2 \right) \end{equation} \tag{ 77 }$$

$$\begin{equation} \sin f\! \sim \pm \delta f, \end{equation} \tag{ 78 }$$

where we write f = δf near the perihelion and f = π + δf near the aphelion, and the upper sign is for the perihelion and the lower sign is for the aphelion. We expand Equations (75) and (76) up to the second order of δf. After the expansion with respect to δf, we take the limit of e → 1. We define

$$\begin{equation} \varepsilon = 1-e \end{equation} \tag{ 79 }$$

and take the lowest order of ε.

The calculations are tedious but straightforward. In the close vicinity of the perihelion, we obtain

$$\begin{equation} {\displaystyle \frac{da}{df}} \sim - K {\displaystyle \frac{\varepsilon ^{2-\alpha }}{a^{2+\alpha }n^4}} {\displaystyle \frac{2(\sqrt{2}-1)}{(3-2\sqrt{2})^{3/2}}} {\displaystyle \frac{1+\delta f^2 (2-\sqrt{2})/8\sqrt{2}}{(1-\delta f^2/4)^{2-\alpha } \lbrace 1+\delta f^2 3(\sqrt{2}-1)/4(3-2\sqrt{2})\rbrace ^{3/2}}} \end{equation} \tag{ 80 }$$

and

$$\begin{equation} {\displaystyle \frac{de}{df}} \sim K {\displaystyle \frac{\varepsilon ^{3-\alpha }}{a^{3+\alpha }n^4}} {\displaystyle \frac{(\sqrt{2}-1)(1-\delta f^2 (5+6\sqrt{2})/4\sqrt{2})}{\sqrt{2}(3-2\sqrt{2})(1-\delta f^2/4)^{7/2-\alpha } (1+\delta f^2 3(\sqrt{2}-1)/4(3-2\sqrt{2})}}. \end{equation} \tag{ 81 }$$

We then integrate over the orbit. We denote the integral by

$$\begin{equation} \overline{{\displaystyle \frac{da}{df}} } = \int _0^{2\pi } {\displaystyle \frac{da}{df}} df. \end{equation} \tag{ 82 }$$

If we integrate over f, the terms in Equations (80) and (81) that contain δf result in a numerical factor of the order of unity. Therefore, close to the perihelion, da/df and de/df averaged over one orbital period are

$$\begin{equation} \overline{{\displaystyle \frac{da}{df}} } \propto a^{4-\alpha } \varepsilon ^{2-\alpha } \end{equation} \tag{ 83 }$$

and

$$\begin{equation} \overline{{\displaystyle \frac{de}{df}} } \propto a^{3-\alpha } \varepsilon ^{3-\alpha }. \end{equation} \tag{ 84 }$$

We now turn our attention to the aphelion. In the close vicinity of the aphelion, we obtain

$$\begin{equation} {\displaystyle \frac{da}{df}} \sim 2^{3-\alpha } K {\displaystyle \frac{1}{a^{2+\alpha }n^4}} {\displaystyle \frac{1-\delta f^2/\varepsilon ^{3/2}}{(1+\delta f^2/2\varepsilon)^{2-\alpha } (1+3\delta f^2/2\varepsilon)^{3/2}}} \end{equation} \tag{ 85 }$$

and

$$\begin{equation} {\displaystyle \frac{de}{df}} \sim -2^{3-\alpha } K {\displaystyle \frac{1}{a^{3+\alpha }n^4 \varepsilon ^{1/2}}} {\displaystyle \frac{1+\delta f^2/\varepsilon ^{3/2}}{(1+\delta f^2/2\varepsilon)^{5/2-\alpha } (1+\delta f^2/\varepsilon ^2)(1+3\delta f^2/2\varepsilon)^{3/2}}}. \end{equation} \tag{ 86 }$$

We now integrate over f near the aphelion. The integration for da/df is rather complicated. We have

$$\begin{equation} \mathcal {I}_a = \int d \delta f {\displaystyle \frac{1-\delta f^2/\varepsilon ^{3/2}}{(1+\delta f^2/2\varepsilon)^{2-\alpha } (1+3\delta f^2/2\varepsilon)^{3/2}}}, \end{equation} \tag{ 87 }$$

and changing the integration variable to X = δf/ε^1/2,

$$\begin{equation} \mathcal {I}_a = \varepsilon ^{1/2} \int dX {\displaystyle \frac{1-X^2/\varepsilon ^{1/2}}{(1+X^2/2)^{2-\alpha } (1+3 X^2 /2)^{3/2}}}. \end{equation} \tag{ 88 }$$

If ε is sufficiently small, the numerator is dominated by the second term, and therefore the resulting $\mathcal {I}_a$ is negative and $\mathcal {I}_a \propto \varepsilon ^0$. If ε is not very small, the numerator is dominated by the first term. Therefore the resulting integral is positive and $\mathcal {I}_a \propto \varepsilon ^{1/2}$. Therefore, the contribution of the aphelion to $\overline{da/df}$ is

$$\begin{eqnarray} \overline{{\displaystyle \frac{da}{df}} } \propto \left\lbrace \begin{array}{@{}cc}a^{4-\alpha } & \mathrm{negative \ for \ small \ \varepsilon } \\ a^{4-\alpha } \varepsilon ^{1/2} & \mathrm{positive \ for \ intermediate \ \varepsilon } \end{array} \right.. \end{eqnarray} \tag{ 89 }$$

The integration for de/df is more straightforward. We have

$$\begin{equation} \mathcal {I}_e = \int d \delta f {\displaystyle \frac{1+\delta f^2/\varepsilon ^{3/2}}{(1+\delta f^2/2\varepsilon)^{5/2-\alpha } (1+\delta f^2/\varepsilon ^2)(1+3\delta f^2/2\varepsilon)^{3/2}}} \end{equation} \tag{ 90 }$$

and changing the variable to X = δf/ε,

$$\begin{equation} \mathcal {I}_e = \varepsilon \int dX {\displaystyle \frac{1+\varepsilon ^{1/2}X^2}{(1+X^2)(1+\varepsilon X^2/2)^{5/2-\alpha } (1+3\varepsilon X^2/2)^{3/2}}}. \end{equation} \tag{ 91 }$$

The terms proportional to ε^1/2 and ε are safely neglected, and therefore we have $\mathcal {I}_e \propto \varepsilon$. Therefore, for the contribution from the aphelion to $\overline{de/df}$, we have

$$\begin{equation} \overline{{\displaystyle \frac{de}{df}} } \propto a^{3-\alpha }\varepsilon ^{1/2}. \end{equation} \tag{ 92 }$$

It still remains to be determined whether contributions from the perihelion or from the aphelion dominate. The relative importance of the contributions from the aphelion and the perihelion depends on the disk model, and we numerically see which contribution is more important as we change α. In Figure 11, we plot the values of $-\overline{da/df}$ and $-\overline{de/df}$ for various disk models. We see that the contribution from the perihelion dominates when α ≳ 2 (corresponding to p ≳ 1 if q = 1), and the expected behaviors of $\overline{da/df}\propto \varepsilon ^{2-\alpha }$ and $\overline{de/df}\propto \varepsilon ^{3-\alpha }$ are observed. If α ≲ 2, the contribution from the aphelion comes into play, and we see $\overline{de/df}$ is roughly proportional to ε^1/2. We note that for the models with α = −1 and α = 0, the direction of semimajor axis evolution is outward for the moderate values of eccentricity (ε > 0.09 for the α = 0 model and ε > 0.05 for α = −1), as we have already seen in the previous section. The inversion of the sign of $\overline{da/df}$ is also expected from the above analytic discussions of the aphelion. We have also checked that the values of $\overline{da/df}$ are always positive if we have a large negative value of p, when the interaction is completely dominated by the contribution from the aphelion.

Zoom In Zoom Out Reset image size

We summarize the dependence of the semimajor axis evolution timescale and eccentricity damping timescale on the disk parameters in the case of a highly eccentric planet. Since $\overline{da/dt}$ and $\overline{da/df}$ are related by

$$\begin{equation} \overline{{\displaystyle \frac{da}{dt}} } = {\displaystyle \frac{1}{T_{\rm orb}}} \overline{{\displaystyle \frac{da}{df}} }, \end{equation} \tag{ 93 }$$

$t_a \propto a^{5/2} \overline{da/df}^{-1}$. In the same way, it is possible to show that $t_e \propto e a^{3/2} \overline{de/df}^{-1} \propto a^{3/2} \overline{de/df}^{-1}$, where we have used e ∼ 1. In the case of α ≳ 2, the contribution from the perihelion dominates and therefore,

$$\begin{equation} t_a \propto a^{\alpha -3/2} \varepsilon ^{\alpha -2} \end{equation} \tag{ 94 }$$

$$\begin{equation} t_e \propto a^{\alpha -3/2} \varepsilon ^{\alpha -3}. \end{equation} \tag{ 95 }$$

The evolution of the semimajor axis is always inward and the eccentricity always damps. If α ≲ 2, we take into account the contribution from the aphelion,

$$\begin{eqnarray} t_a \propto \left\lbrace \begin{array}{@{}cc}a^{\alpha -3/2} & \mathrm{inward \ for \ small \ \varepsilon } \\ a^{\alpha -3/2} \varepsilon ^{-1/2} & \mathrm{outward \ for \ intermediate \ \varepsilon } \end{array} \right. \end{eqnarray} \tag{ 96 }$$

and

$$\begin{equation} t_e \propto a^{\alpha -3/2}\varepsilon ^{-1/2}. \end{equation} \tag{ 97 }$$

It is noted that the outward evolution of the semimajor axis occurs only when α ≲ 0, and eccentricity always damps. The range of ε where the planet experiences the outward semimajor axis evolution increases as the value of α is decreased. In the particular disk model with q = 1, α and p are related by α = p + 1.

We have shown that there is a critical value for the power of the surface density, which determines the behavior of semimajor and eccentricity evolutions for a highly eccentric planet (e ∼ 1). Let us consider the model with q = 1. If p > 2, both eccentricity and semimajor axis evolution timescales become small as we increase the planet's eccentricity. If 1 < p < 2, the semimajor axis evolution timescale becomes shorter for higher eccentricity, while the eccentricity evolution timescale becomes longer. For p < 1, both evolution timescales become longer as we increase the eccentricity of the planet. Further complications arise for the semimajor axis evolution for p ≲ 0, where the contribution from the aphelion comes into play to invert the direction of the evolution.

However, we note that these critical values of p should not be overstated. In the model described above, the critical parameter that determines the behavior of the orbital evolution timescales is α = p + (3 − q)/2. The appearance of this single parameter α is partly due to the prescription of the cutoff of the gravitational force. Note that the dynamical friction force we have used depends on the cutoff scale , and we have used the model in which depends linearly on the disk scale height. The appearance of q in the parameter α depends on this simplified prescription of the cutoff. We also note that the dependence of the timescales on the softening length is different from Papaloizou & Larwood (2000), which is due to the difference in the formulation.

Therefore, we consider that the critical values of p described above only have qualitative meanings. Yet, we expect that the model describes the qualitative behavior of the semimajor axis and the eccentricity evolution of highly eccentric planets. The prediction is readily compared with two-dimensional numerical calculations, although a very large simulation box may be necessary. In order to eliminate the dependence on the softening parameter , it is necessary to perform three-dimensional analyses.

In this section, we briefly discuss the applicability of our model. The most important assumption we have used in the model is the use of the supersonic dynamical friction formula (41), which is derived from linear perturbation analyses. The use of this formula assumes that (1) the relative velocity is supersonic, (2) the flow in the vicinity of the planet is homogeneous, and (3) the steady state is reached instantaneously. We shall now discuss each assumption separately. We will also discuss how nonlinear effects on the dynamical friction force may change our results.

The first assumption is justified if we consider a planet that has high enough eccentricity, as indicated in Figure 5. Typically, when e ≫ H/r, we can safely assume that the flow is supersonic.

For the second assumption, let us consider which scale of the perturbation contributes to the force. In the discussion of Section 2, we have seen that the perturbation of all the scales larger than the cutoff length equally contribute to the dynamical friction force. Let us consider a planet embedded in a disk with Keplerian rotation, located at the perihelion or aphelion. Then, there is a location in the disk where the relative velocity between the gas and the planet becomes zero. We call this radius an "instantaneous corotation radius," r_C, given by, for the perihelion,

$$\begin{equation} {\displaystyle \frac{GM_{\ast }}{r_C}} = {\displaystyle \frac{GM_{\ast }}{a}} {\displaystyle \frac{1+e}{1-e}}. \end{equation} \tag{ 98 }$$

A similar equation can be derived for the aphelion.

The assumption of the homogeneous flow breaks down if the length scale comparable to the distance between the planet's location and the instantaneous corotation radius is important in determining the force acting on the planet. If the planet's eccentricity is large, the distance between the instantaneous corotation radius r_C and the planet's location is of the order of the scale of the disk itself (|r_C − a| ∼ a), which is much larger than the cutoff scale, which is comparable with the disk scale height. Therefore, in this case, we expect that the assumption of the homogeneous flow is justified. If the planet's eccentricity is small, the instantaneous corotation radius comes very close to the planet's orbital radius (r_C ∼ a), and the assumption of a homogeneous medium becomes worse. Typically, we expect that this assumption breaks down if the planet's eccentricity is smaller than the disk aspect ratio, and in such a small eccentricity case, it is necessary to fully take into account the effects of shear around the planet's location.

Regarding the third assumption, the time taken to reach the steady state may be estimated as follows: The length scale that is important in determining the force is of the order of the disk scale height, and the minimum speed of the propagation of the information of the perturbing potential may be the sound speed. This indicates that the time taken to develop the steady state is of the order of the Kepler timescale. However, since the relative velocity is supersonic, the time taken to develop the steady state is less than that, since the background flow will carry the information of the perturber. In any case, this assumption is only marginally satisfied.

As previously discussed, it is known from numerical simulations that the time when the maxima/minima of the torque occurs lags the time of the maxima/minima of the distance between the planet and the central star (e.g., Cresswell et al. 2007). Our model does not give a precise prescription for these effects, and the instantaneous torque or power profiles are different from what we expect from the numerical simulations. However, if we take the average over one orbital period, our model is in good agreement with the previous calculations in the parameter range where both our model and the previous calculations are expected to give reasonable results.

In short, we expect that our formulation is applicable to planets whose eccentricity is larger than the disk aspect ratio, in particular when one takes the average over one orbital period.

Finally, we must make a comment on the use of the dynamical friction formula derived by the linear perturbation analysis. Kim & Kim (2009) investigated the dynamical friction using the nonlinear numerical simulations. They considered the case in which the particle is embedded in a homogeneous, three-dimensional gas flow, and calculated the axisymmetric pattern of the gas flow induced by the particle's gravity. They found that there is a deviation from the linear theory depending on the Mach number M₀ of the flow and the nonlinear parameter

$$\begin{equation} \mathcal {A} = {\displaystyle \frac{GM}{c^2 r_S}}, \end{equation} \tag{ 99 }$$

where r_S is the (effective) radius of the body. The force exerted on the particle is deviated from the linear results by a factor of (η/2)^−0.45, where $\eta =\mathcal {A}/(M_0^2-1)$. In the case of a planet embedded in a protoplanetary disk, $\mathcal {A}$ is of the order of 10–100; therefore, nonlinearity can be important. The correction factor, however, is the quantity of the order of unity, (since the Mach number is of the order of 5–10); therefore, our results may give a reasonable estimate of the timescale of the evolution of the orbital parameter. One possible effect of such nonlinearity on our results is that the dependence of the timescales t_a and t_e on the orbital parameters of the planet can be different, since the force depends on the velocity of the perturbing potential in a different way.

In this section, we derive the fitting formula for t_a and t_e. As seen before, within our framework, t_a, t_e∝M⁻¹_pΣ⁻¹₀₀H₀, and the only disk parameter that controls the timescale is α = p + (3 − q)/2. Therefore, we fit the values of t_a and t_e as a function of the semimajor axis and the eccentricity for a given value of α with 1 ⩽ α ⩽ 4. For the data for fitting, we use the values obtained for 0.3 ⩽ a/r₀ ⩽ 10 and 0.1 ⩽ e ⩽ 0.99.

Motivated by the results of the analytic considerations for e ∼ 1, we use the fitting formula of the form

$$\begin{equation} t_a = C_a(a) e^{\zeta _a(a)} (1-e^{\lambda _a(a)})^{\eta _a(a)} \left({\displaystyle \frac{M_p/M_{\ast }}{10^{-6}}} \right)^{-1} \left({\displaystyle \frac{\Sigma _0/M_{\ast }r_0^{-2}}{10^{-4}}} \right)^{-1} \left({\displaystyle \frac{h_0}{0.05}} \right) \left({\displaystyle \frac{\epsilon _0}{0.5}} \right) \end{equation} \tag{ 100 }$$

Here, C_a, ζ_a, λ_a, and η_a are the function of a of the form

$$\begin{eqnarray} && C_a(a) = C_{a0} \left({\displaystyle \frac{a}{r_0}} \right)^{C_{ap}} \end{eqnarray} \tag{ 101 }$$

$$\begin{eqnarray} && \zeta _a(a) = \zeta _{a0} \left({\displaystyle \frac{a}{r_0}} \right)^{\zeta _{ap}} \end{eqnarray} \tag{ 102 }$$

$$\begin{eqnarray} && \lambda _a(a) = \lambda _{a0} \left({\displaystyle \frac{a}{r_0}} \right)^{\lambda _{ap}} \end{eqnarray} \tag{ 103 }$$

$$\begin{eqnarray} && \eta _a(a) = \eta _{a0} \left({\displaystyle \frac{a}{r_0}} \right)^{\eta _{ap}}. \end{eqnarray} \tag{ 104 }$$

The same form of the fitting function is used for t_e,

$$\begin{equation} t_e = C_e(a) e^{\zeta _e(a)} (1-e^{\lambda _e(a)})^{\eta _e(a)} \left({\displaystyle \frac{M_p/M_{\ast }}{10^{-6}}} \right)^{-1} \left({\displaystyle \frac{\Sigma _0/M_{\ast }r_0^{-2}}{10^{-4}}} \right)^{-1} \left({\displaystyle \frac{h_0}{0.05}} \right) \left({\displaystyle \frac{\epsilon _0}{0.5}} \right), \end{equation} \tag{ 105 }$$

with

$$\begin{eqnarray} && C_e(a) = C_{e0} \left({\displaystyle \frac{a}{r_0}} \right)^{C_{ep}} \end{eqnarray} \tag{ 106 }$$

$$\begin{eqnarray} && \zeta _e(a) = \zeta _{e0} \left({\displaystyle \frac{a}{r_0}} \right)^{\zeta _{ep}} \end{eqnarray} \tag{ 107 }$$

$$\begin{eqnarray} && \lambda _e(a) = \lambda _{e0} \left({\displaystyle \frac{a}{r_0}} \right)^{\lambda _{ep}} \end{eqnarray} \tag{ 108 }$$

$$\begin{eqnarray} && \eta _e(a) = \eta _{e0} \left({\displaystyle \frac{a}{r_0}} \right)^{\eta _{ep}}. \end{eqnarray} \tag{ 109 }$$

We note that in the set of the fitting parameters, C_a0 and C_e0 have the dimension of time, while others are dimensionless.

In Tables 1 and 2, we show the results of fitting for t_a and t_e, respectively. For C_a0 and C_e0, we give the values in the units of year with M_* = M_☉ and r₀ = 1 AU. In the last column of the tables, we show the maximum values of the error of the fitting function, which is calculated as |t_calc − t_fit|/t_calc, where t_calc is the value of t_a or t_e calculated in the model presented in this paper and t_fit is the value of the fitting function. Our fitting formula successfully reproduces the values of t_a and t_e within a 10%–20% error. We also note that a^{α − 3/2} dependence that appears in Equations (94)–(97) can be seen clearly in the fitting parameters C_ap and C_ep.

Here, we note that caution must be paid in using the fitting formulae. In obtaining the fitting parameters, we have used the values with 0.3 ⩽ a/r₀ ⩽ 10 and 0.1 ⩽ e ⩽ 0.99. Our formulation is not applicable to the case of a low-eccentricity planet (typically e ≲ H/r). It is necessary to combine the results of this paper and such formulae given by Tanaka & Ward (2004), which are applicable to the case of low eccentricity.

It is also noted that the fitting formulae are given in the case of 1 ⩽ α ⩽ 4. Our form of fitting functions (100) and (105) fails to describe the change of the sign of t_a that happens in the case of p ≲ −1 (see Figure 11). Although we believe that this range of α covers most of the typical protoplanetary disk parameters, it is necessary to construct the fitting formula if one wants to apply for protoplanetary disk parameters out of the range of this fitting. Even in this case, however, we expect that the model described in this paper, which uses the dynamical friction force to calculate the evolution of orbital parameters, is still applicable.

We have developed a model for the disk–planet interaction that can be used specifically for highly eccentric planets. We have found that in the case of a highly eccentric planet, the disk–planet interaction is essentially described by dynamical friction. We have derived how the evolution timescales of the planet's semimajor axis and eccentricity depend on disk parameters within our models. In this section, we discuss the possible implications for the planet formation theory.

Recently, Machida et al. (2011) suggested that a disk surrounding a newly born protostar is gravitationally unstable, and that it therefore may be possible to form a planet through gravitational instability (see also Inutsuka et al. 2010). Planets born in such a way are expected to have high eccentricity. Also, it is possible to form a planet at a large distance from the central star.

We briefly discuss the evolution timescale of such planets. From the fitting formula of the evolution timescales, it is possible to show that the orbital evolution timescales of the newly born Jupiter-mass planets through the gravitational instability can be very rapid, since the disk is very massive compared to the central star. However, after the accretion phase onto the central star, the stellar mass becomes much larger than the disk mass. If there is a low-mass planet (e.g., a sub-Jupiter mass) that has survived in the early phase and a reasonable amount of eccentricity is maintained, the orbital evolution timescale can be of an order of magnitude comparable with the observed disk lifetime (Haisch et al. 2001). Such planets may be found by the direct imaging observations (e.g., Marois et al. 2008), although small mass planets may be very cold and the observations may be difficult. Nevertheless, a highly eccentric planet at a large orbital separation can be a good observational signature to test the theory of the planet formation and disk–planet interaction.

Another possible mechanism for producing a planet in a highly eccentric orbit is planet–planet scattering (e.g., Chambers et al. 1996; Marzari & Weidenschilling 2002). The scattering between planets can naturally result in highly eccentric planets. From the calculations given in this paper, the evolution timescale of the orbital parameters can be very long even if the planet's orbital plane coincides with the gas disk plane. If the planet's orbital plane is not aligned with the plane of the disk, we expect that the evolution timescale becomes much larger, since the planet interacts with the disk only at a very small part of the orbit. The calculations given in this paper set the lower limit of the timescales of the orbital evolution for such planets. It is therefore possible that the scattered planets remain in a highly eccentric orbit, even if the gas disk remains and the tidal interaction with the central star is not effective.

Another interesting application of our model is the recently suggested "eccentricity trap" mechanism in Ogihara et al. (2010). They have suggested that if there is an inner cavity in a protoplanetary disk, it may be possible to maintain a planet at the inner edge of the disk. This is due to the positive torque acting on the planet at the aphelion, which stops the planet from migrating inward. Although our formulation, which assumes a smooth disk profile, is not directly applicable to the disk with a sharp inner edge, we have already seen that the evolution of the semimajor axis can be outward if there is a steep surface density gradient. It is an interesting extension to apply the method of dynamical friction to a disk with a cavity.