6 Tidal torques in an eccentric disc

We have used the above formalism to estimate eccentricity excitation/damping rates for protoplanets in eccentric discs which can be described using the normal modes calculated above. Before giving details we give a brief summary of our results. Although we consider eccentricities which may substantially exceed H/r, we shall still suppose them sufficiently small $\le \sim 0.2$that we may consider there to be an equilibrium eccentricity as a function of disc radius. Then a torque calculation is characterized by a disc eccentricity, e, and an equilibrium eccentricity which we may also consider to be the protoplanet eccentricity $e_{\rm p}.$ We shall also for the most part restrict consideration to the case when the protoplanet and disc orbits are aligned in equilibrium $(\varpi =0)$. In general when $e_{\rm p} \ll e$ we find excitation of the protoplanet orbit eccentricity, $t_{\rm e} >0,$while for $e_{\rm p} \gg e,$ the eccentricity damps as expected, $t_{\rm e} <0.$The transition between these regimes occurs when e and $e_{\rm p}$ are approximately equal. For $e_{\rm p}$ significantly larger than e, provided $\vert(\Omega_{\rm p} - \omega_{\rm pg}) t_{\rm e}\vert$ is significantly greater than unity, there is an equilibrium solution with apsidal line slightly rotated from zero. The orbit may then suffer significantly reduced or even reversed torques for $e_{\rm p}$ sufficiently large.

Figure 12: In this figure $t_{\rm J}$ in yr is plotted as a function of protoplanet equilibrium eccentricity (negative values correspond to torque reversal ) for $N_{\rm t} =1.$The short dashed curves correspond to an assumed disc eccentricity e= 0.07 and the short dashed long dashed curves to e = 0.05. The disc aspect ratio H/r was taken as 0.07.The crosses were obtained for e=0.05 and H/r=0.05,while the squares correspond to e=0.07 and H/r =0.05.In the above cases disc and protoplanet orbit apsidal lines were taken to be aligned. The diamonds give $t_{\rm J}$ for assumed antialignment and e = H/r = 0.07.Note the weakening of the disc interaction at the higher protoplanet eccentricities due to larger relative velocities especially in the antialigned case, to the disc.

When $e_{\rm p}$ is significantly less than a sufficiently large disc eccentricity e, the protoplanet orbital eccentricity can grow until the orbit ceases to be aligned with the disc and precesses through a full $(2\pi)$ cycle at which point it then damps. Again inward orbital migration may be reduced or reversed. Many of these features can be traced to the fact that when the equilibrium eccentricity is such that $e_{\rm p} =e$ the situation in many ways replaces the equilibrium eccentricity solution $e_{\rm p}=0$ in an axisymmetric disc. This we show below

Figure 13: In this figure $t_{\rm e}$ in yr is plotted as a function of protoplanet equilibrium eccentricity (negative values correspond to circularization) for $N_{\rm t} =1$. The short dashed curves correspond to an assumed disc eccentricity e= 0.07 and the short dashed long dashed curves to e = 0.05. The disc aspect ratio H/r was taken as 0.07.The crosses were obtained for e=0.05 and H/r=0.05,while the squares correspond to e=0.07 and H/r =0.05.In the above cases disc and protoplanet orbit apsidal lines were taken to be aligned. The diamonds give $t_{\rm e}$ for an orbit with apsidal line antialigned with that of the disc and e = H/r = 0.07.Note the eccentricity grows for protoplanet eccentricity less than e in the aligned case while there is always decay in the antialigned case.

We use the expressions for $r, \varphi, R,$ and $\nu,$given by Eqs. (37), (39) in the forcing potential (38). In the case when $e = e_{\rm p},$ and the protoplanet orbit is almost aligned with the disc with very small $\vert\varpi\vert,$the strongest interaction occurs when r=R, and $\varphi =\nu.$ This corresponds to $a=a_{\rm p},$ and $M= \omega t + \varpi.$ Performing a first order Taylor expansion about the point of maximum interaction the forcing potential becomes

$$\Phi_{\rm p}' = -{Gm_{\rm p}\over \sqrt{\vert(a - a_{\rm p})^2+4aa_{\rm p}\sin^2(M -\omega t -\varpi) +b^2\vert}} \cdot$$ (50)

Here $\vert a-a_{\rm p}\vert$ and $a\vert M -\omega t\vert$ are considered as small compared to a and terms of order e times these small quantities have been neglected. Equation (50) is to be considered in comparison to the similar expression appropriate to an axisymmetric disc

$$\Phi_{\rm p}' = -{Gm_{\rm p}\over \sqrt{\vert(r - R)^2+4rR\sin^2(\varphi -\nu) +b^2\vert}} \cdot$$ (51)

Given that we have already shown that a,M behave just like cylindrical coordinates $(r,\varphi)$ the tidal torque calculation in an eccentric disc with $e_{\rm p} =e$ should give the same results as an axisymmetric disc with $e_{\rm p} =0.$Note too that only one term should remain in the sum (42) such that k₁ =m, and k₂ = -m.Then there is zero rate of change of eccentricity. Thus the situation of aligned orbits such that $e_{\rm p} =e$behaves much like the case with $e_{\rm p}=0$ in an axisymmetric disc in that it is one of steady eccentricity. The orbit then migrates at the same rate as a circular orbit in an axisymmetric disc. This is essentially what we have found on application of the torque formulae.

6.1 Numerical results

We have calculated $t_{\rm J}/N_{\rm t}$ and $t_{\rm e}/N_{\rm t}$ by summing the contributions from appropriate resonances. The normalizing factor

$$N_{\rm t} =(r/5~ {\rm AU})^{-1/2}(m_{\rm p}/(3~M_{\oplus}))^{-1}(\Sigma/56)^{-1}.$$

Here the distance is measured in units of 5 AU, the protoplanet mass in units of 3 earth masses and the disc surface density in units of $56~\rm gm~ cm^{-2}.$ We have also taken the central mass to be one solar mass. The softening parameter was taken to be $b=H/\sqrt{2}.$ Our results are consistent with those of Ward (1997) in the limit where both e and $e_{\rm p}$ tend to zero. We plot $t_{\rm J}/N_{\rm t}$ as a function of equilibrium protoplanet orbit eccentricity $e_{\rm p}$ for e=0.05and e=0.07, when disc and protoplanet apsidal lines are aligned in Fig. 12. Results for aspect ratios 0.05 and 0.07 are presented.

We plot $t_{\rm e}$ for $N_{\rm t} =1$ in Fig. 13. The trends in all cases are similar and are that for small $e_{\rm p}$ the protoplanet eccentricity grows in the aligned case while inward migration occurs with $t_{\rm J} \sim 2\times 10^6$ yr. The eccentricity growth reverses for $e_{\rm p} \sim e$ indicating an equilibrium in accordance with the discussion above. We comment that for very small $e_{\rm p}$ and finite e, $de_{\rm p}/{\rm d}t \propto e.$For larger $e_{\rm p},$ $t_{\rm J}$ and $t_{\rm e}$ increase and $t_{\rm J}$eventually changes sign for $e_{\rm p}$ exceeding $\sim$0.1. At these eccentricities $\vert t_{\rm e}\vert \sim 2\times 10^4$ yr. We make the comment that the same dimensionless units can be used for $t_{\rm J}$ and $t_{\rm e}$ as for the disc models introduced in Sect. 2 and thus the same scaling to make results applicable to different radii may be used.

In the antialigned case the disc protoplanet interaction is much weaker (see Figs. 12 and 13). Note too that the interaction with the disc weakens in general for larger $e_{\rm p}$ because of the larger relative velocity of the protoplanet with respect to the disc. This results in larger values of $\vert t_{\rm e}\vert$ and $\vert t_{\rm J}\vert.$Additional calculations have shown that, as expected, these values become independent of the orientation of the apsidal line when $e_{\rm p} \gg e.$

Thus the indications are that for modest eccentricities exceeding a few H/r for both disc and protoplanet the tidal interaction may differ significantly from the circular disc and small protoplanet eccentricity case. We now consider applications to the normal modes calculated in this paper.

The equilibrium protoplanet eccentricities associated with the two highest frequency modes calculated in the case of disc model A with two interior protoplanets are shown in Fig. 2 while those corresponding to disc model B are illustrated in Fig. 3. These modes are associated with significant internal protoplanet eccentricities and can be thought of as giving the disc response to external forcing. Equilibrium eccentricities corresponding to the normal modes are also plotted in Figs. 2 and 3. These are generally larger than the disc eccentricities in the inner parts of the disc.

Equilibrium protoplanet eccentricities associated with the modes calculated in the case of disc model A with one interior protoplanet are shown in Fig. 8 while those corresponding to disc model B are illustrated in Fig. 9. Equilibrium eccentricities in the case of isolated disc model A are plotted in Fig. 10 while those corresponding to disc model B are given in Fig. 11. In all of these cases the form of the equilibrium eccentricity curve tracks that of the corresponding normal mode according to increasing number of nodes. Thus the mode with the largest number of nodes has associated equilibrium eccentricity with the largest number of nodes.

By comparing the equilibrium eccentricities with their corresponding normal modes one sees that the modes with the smallest frequencies or pattern speeds in absolute magnitude tend to have high equilibrium eccentricities several times larger than the disc eccentricity. These correspond to the longest wavelength curves in Figs. 8 and 10 with corresponding modes plotted in Figs. 4 and 6. From our discussion above these are expected to facilitate high embedded protoplanet eccentricities. The reason the protoplanet eccentricity is significantly larger than the disc eccentricity for these modes is that, for the disc mode the effects of the nonaxisymmetric forces due to self-gravity and pressure which drive the eccentricity (see Eq. (18)) tend to cancel. However, the protoplanet is subject only to self-gravity with no cancelling effects from pressure forces. Therefore the equilibrium eccentricity is larger. Note that these low frequency modes are essentially disc modes and are associated with low interior protoplanet eccentricities when the latter are present. Such disc modes may also be associated with high embedded protoplanet eccentricities through secular resonances (see Sect. 4 above). In our models such resonances occur when the pattern speed is slightly prograde. An example is shown in Fig. 9. This occurs for the one protoplanet model with disc B at about 20 times the radius of the disc inner edge. Such a resonance also occurs when disc model A is used but for a higher order mode at smaller radii (see Fig. 8).

To give numerical examples we first consider the two protoplanet model with disc A as an approximation to the Upsilon Andromedae system. Taking the situation represented in Fig. 2 with inner boundary disc eccentricity 0.1, for a semi-major axis of 2.57 AU for the outer protoplanet in the inner cavity, being 0.6 of the inner boundary radius, r=1.5 corresponds to 6.4 AU. For the mode with $\Omega_{\rm p} = 6.99\times 10^{-4}$ and disc model A, $\Sigma=136~\rm gm~ cm^{-2}$ at r = 1.5.At this radius e=0.07, and the equilibrium eccentricity $e_{\rm p}=0.1.$ From Fig. 12 the orbital migration rate is very small and possibly outwards. We also find $\vert t_{\rm e}\vert = 10^4 (3M_{\oplus}/m_{\rm p})$ yr and $\vert(\omega_{\rm pg} -\Omega_{\rm p}) t_{\rm e}\vert \sim 30.0 (M_{\oplus}/m_{\rm p}).$Thus in this case the eccentric disc has significant effects on the tidal torques acting on an embedded protoplanet.

The eccentricity of the outermost protoplanet orbit in the inner cavity would be $\sim$0.1. The currently observed value which is three or so times larger would be attained if the disc inner boundary was scaled to be at twice the outer protoplanet semi-major axis. Very similar results for the migration and circularization rates to those obtained above would still apply.

As an example to illustrate a case involving a very low frequency mode we consider the one protoplanet model with disc A. The modes are plotted in Fig. 4 and the equilibrium eccentricities in Fig. 8.

For the mode with $\Omega_{\rm p} = 4.02\times 10^{-6},$we find that at r=10 $e_{\rm p} =2.5 e.$Supposing that r=10 corresponds to 5 AU, $\Sigma=426~\rm gm~ cm^{-2}.$ From Fig. 12 we find that for e=0.05, there is torque reversal. We also again find that $\vert t_{\rm e}\vert \sim 3\times 10^3 (3M_{\oplus}/m_{\rm p})$ yr and $\vert(\omega_{\rm pg} -\Omega_{\rm p}) t_{\rm e}\vert \sim 30.0 (M_{\oplus}/m_{\rm p}).$

These examples indicate that protoplanets embedded in eccentric discs will be in eccentric orbits and that based on resonant torque calculations orbital migration slows down and may even reverse to become outward when the protoplanet eccentricity is sufficiently large. However, we should emphasize that these calculations are approximate and somewhat uncertain due to the cancellation of torques arising at inner and outer Lindblad resonances. To further examine the issue of protoplanet disc tidal interaction at a large eccentricity ( compared to H/r and e) we present below a simpler calculation based on local dynamical friction which should apply in the appropriate limit and which is in essential qualitative agreement with the resonant torque calculations.

6.2 Dynamical friction calculation

We here consider the case when $e_{\rm p}$ significantly exceeds H/rbut is still also significantly less than unity. This can be realized in sufficiently thin discs. In this situation the motion of the protoplanet through the disc is supersonic so we neglect pressure forces. In addition the scale of the response in both space and time becomes local (even though the protoplanet may move globally through the disc). For example, for a length scale H, the response time scale would be ${\sim} H/(e_{\rm p}r\Omega) \ll \Omega^{-1}.$

Accordingly we work in a reference frame moving instantaneously with the protoplanet in which the disc material appears to move with velocity ${\vec v}.$ Adopting local Cartesian coordinates, we suppose that the perturbing potential due to the protoplanet may be written as Fourier integral

$$\Phi_{\rm p}'({\vec r}) =\int \Phi_{\rm p}'({\vec k})\exp(i{\vec k}{\vec r}) {\rm d}^N{\vec k}.$$ (52)

Here we shall consider both the two dimensional case (N=2)and the three dimensional case (N=3). For the two dimensional case with $\Phi_{\rm p}'({\vec r}) =-Gm_{\rm p}/\sqrt{r^2+b^2},$ $\Phi_{\rm p}'({\vec k}) = -Gm_{\rm p}\exp(-bk)/(2\pi k),$ with $k=\vert{\vec k}\vert.$In the three dimensional case with $\Phi_{\rm p}'({\vec r}) =-Gm_{\rm p}/r,$ $\Phi_{\rm p}'({\vec k}) = -Gm_{\rm p}/(2\pi^2 k^2).$

In each case, assuming a local steady state, the velocity ${\vec v}'$ induced by the protoplanet is found from

$${\vec v}\nabla {\vec v}' = -\nabla \Phi_{\rm p}' .$$ (53)

We also have

$$\nabla \cdot(D {\vec v}')= -\nabla(\cdot D'{\vec v}),$$ (54)

where D denotes the density (N=3) or the surface density (N=2).In this case where only inertial terms are retained the analysis is the same as would be performed for collisionless particles (e.g. Tremaine & Weinberg 1984). In this context we note that this type of calculation is also applicable to the frictional interaction with material in planetesimal form as well as with the gas disc provided the surface density and disc thickness are appropriately specified.

The rate of change of disc momentum ${\vec{\dot P}},$which gives rise to a frictional force on the protoplanet acting in the direction of it's relative velocity, may then be calculated from

$${\vec {\dot P}} = -\int D'\nabla \Phi_{\rm p}' {\rm d}^N{\vec r}.$$ (55)

Performing an integration by parts and working in terms of the Fourier transforms of perturbations, one obtains

$${\vec {\dot P}} \cdot {\vec v}= - (2\pi)^N \int D \Phi_{\rm ... ...\vec k}) i{\vec k}\cdot {\vec v}'({\vec k}){\rm d}^N{\vec k} .$$ (56)

From (53) one obtains

$${\vec v}'({\vec k})\cdot{\vec k} = {-k^2\Phi_{\rm p}'({\vec k})\over {\vec v}\cdot{\vec k}}\cdot$$ (57)

In order to perform the integral (56) one has to apply a Landau prescription by adding an infinitessimally small negative imaginary part to the denominator ${\vec v}\cdot{\vec k}.$One then finds

$${\vec {\dot P}} \cdot {\vec v}= -{\pi D (Gm_{\rm p})^2\over v}{\cal Q},$$ (58)

where for N=3, ${\cal Q}=4\ln(k_{\max}/k_{\min}),$and for N=2, ${\cal Q}= 1/(2b).$ Here $(k_{\max},k_{\min})$ are the usual upper and lower wavenumber cut offs (e.g. Tremaine & Weinberg 1984). Here reasonable values are $k_{\max} = v^2/(Gm_{\rm p}),$ and $k_{\min} = 1/H.$The two and three dimensional cases are thus of the same form. The logarithmic factor is generally of order unity. Comparison of the two and three dimensional cases suggests that $\rho \ln(k_{\max}/k_{\min}) =\Sigma/(8b).$Thus the adopted softening parameter b should be somewhat smaller than H.

Using (58) we may evaluate, remembering that ${\vec v}$ is the relative velocity between disc and protoplanet, the average rate of change of angular momentum of the protoplanet from

$$\left \langle {{\rm d}J\over {\rm d}t}\right \rangle = {\omeg... ...imes{\vec v})\cdot{\hat{\vec k}}\over v^3} {\cal Q} {\rm d}t ,$$ (59)

with the integral being taken round the orbit and ${\hat{\vec k}}$ being the unit vector in the direction normal to the disc.

In the two dimensional case for $D \equiv \Sigma \propto r^{-3/2}$one obtains for small $e_{\rm p}$ and e =0.

$${1\over J}\left\langle{{\rm d}J\over {\rm d}t}\right\rangle \... ... \Sigma(a_{\rm p}) a_{\rm p}^3 \over 2 b e_{\rm p} M_*^2} \cdot$$ (60)

As this is positive it corresponds to outward migration as long as the eccentricity can be maintained. This occurs because the disc flow tends to speed up the protoplanet at apocentre where most time is spent. Numerically for $b=H/\sqrt{2}$

$${t_{\rm J}\over N_{\rm t}} = - 2.6\times 10^6\left({ H\over 0... ...right)^2 \left({ e_{\rm p} a_{\rm p}\over H}\right) {\rm yr}.$$ (61)

While this agrees in form and is of comparable magnitude to what is obtained from resonant torque calculations (see Fig. 12) it is impractical to perform the latter at higher eccentricities where better agreement might be attained because of the small scale of the interaction (even at the eccentricities plotted, over 10⁶ resonances were included).

We may also use the above formalism to calculate the mean rate of change of orbital energy for the protoplanet to be

 

$$\left \langle {{\rm d}E\over {\rm d}t}\right \rangle =$$

$$~~~~ -{\omega\over 2\pi} \int{\pi D (Gm_{\rm p})^2 {\vec v}\cdot ... ...GM_* r^{-1/2} {\hat {\mbox{\boldmath$\varphi$ }}})\over v^3} {\cal Q} {\rm d}t,$$ (62)

with $E = - {GM_* m_{\rm p}\over 2 a_{\rm p}}$ and ${\hat {\mbox{\boldmath$\varphi$ }}}$ being the unit vector in the azimuthal direction.

Performing the integration we find for small $e_{\rm p}$ and assuming $\Sigma \propto r^{-n}$ that

$${1\over E}\left\langle{{\rm d}E\over {\rm d}t}\right\rangle ... ...p}^3 \over 2 b e_{\rm p} M_*^2} \left(-3.01 - 7.17n \right) .$$ (63)

From this we see that, for this circular disc case, the mean rate of change of orbital energy increases for n < -0.42.This simplified calculation indicates outward migration for density profiles that do not increase too rapidly inwards in line with the idea that the effect is caused by the interaction at apocentre.

We further comment that the more sensitive dependence on the softening parameter means that two dimensional calculations of the type carried out here, require a precise specification of this parameter that correctly represents three dimensional effects in order for them to be very accurate. Thus two dimensional torque calculations and use of torque formulae such as (43) suffer from a number of uncertainties which can be of comparable importance.