5 Response of an eccentric disc to a low mass protoplanet

To calculate the tidal response of an eccentric disc to a perturbing protoplanet it is convenient to introduce a new orthogonal coordinate system $(a,\lambda),$ with a being the local semi-major axis and $\lambda$being the orthogonal angular coordinate. Then taking pericentre in the disc to lie along $\varphi =0,$ $a=r(1+e\cos\varphi)/(1-e^2),$with the eccentricity being a function of a.We assume the eccentricity is small enough that we may work to first order in it and also that it has significant changes only on a global length scale comparable to r. Then we may make the replacement $e(a) \equiv e(r).$In this approximation a can be expressed in terms of cylindrical coordinates in the form $a=r(1+e\cos\varphi)$ and the orthogonal angular coordinate, $\lambda$ is given by

$$\lambda= \varphi +\int\left({e\over r}{\rm d}r\right)\sin\varphi,$$ (33)

where we may leave the integral as indefinite. It is also convenient to work in a frame in which the disc flow appears stationary. As the pattern speed as measured in an inertial frame is very small compared to the disc rotation speed, we shall neglect centrifugal and coriolis forces.

In a general orthogonal coordinate system in which the disc appears to be in a steady state, the equations of motion for the velocity $( v_{\rm a}, v_{\lambda})$may then be written (see Appendix 2)

  

$${\partial v_{\rm a} \over\partial t}+ v_{\rm a}\vert\nabla a\vert... ...partial \left({v_{\rm a}\over \vert\nabla a\vert} \right)\over\partial \lambda}$$

$$~~~~-{v_{\lambda}^2 \over 2} \vert\nabla a\vert\vert\nabla \lambd... ...\nabla a\vert^3 {\partial \left(\vert\nabla a\vert^{-2} \right)\over\partial a}$$

$$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~= -{\vert\nabla a\vert\over \Si... ...artial \Pi \over\partial a}- \vert\nabla a\vert{\partial \Phi \over\partial a},$$ (34)

$${\partial v_{\lambda} \over\partial t}+ v_{\lambda}\vert\nabla \l... ...tial \left( {v_{\lambda}\over \vert\nabla \lambda\vert} \right)\over\partial a}$$

$$~~~~-{v_{\rm a}^2 \over 2} \vert\nabla a\vert^2\vert\nabla \lambd... ...t^3 {\partial \left(\vert\nabla \lambda\vert^{-2}\right) \over\partial \lambda}$$

$$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~= -{\vert\nabla \lambda\vert\ov... ...al \lambda}- \vert\nabla \lambda\vert{\partial \Phi \over\partial \lambda}\cdot$$ (35)

We perform a response calculation by linearizing about a steady state in which the disc gas is taken to be in elliptical Keplerian orbits with $v_{\rm a} =0,$and to first order in e, $v_{\lambda}=(GM_*/a)^{1/2}(1+e\cos\lambda).$In a thin disc forces due pressure and self-gravity provide only a small correction of order (H/r)². We consider the situation when the disc eccentricity, e(r), protoplanet eccentricity, $e_{\rm p},$and H/rcan be considered small and comparable in magnitude.

When considering the response to an embedded protoplanet, we are interested in responses with a scale H and azimuthal mode number $m \sim r/H.$ Here for the time being $\lambda$ is regarded as the azimuthal angle. We perform a response calculation by linearizing Eqs. (34) and (35) about the steady state described above. We may expand the linearized equations in powers of the disc eccentricity.

Here we shall work only to lowest order and neglect terms of order e times smaller than the dominant ones which, linearizing (34) and (35) directly and denoting perturbation quantities by a prime, are of order $\Omega v_{\rm a}'.$In this scheme we may replace $\vert\nabla a\vert$ and $\vert\nabla \lambda \vert$ by unity except where the latter occurs in the combination $\vert\nabla \lambda \vert v_{\lambda}{\partial \over \partial \lambda},$ with $\vert v_{\lambda}\vert$ here denoting the unperturbed velocity in the disc. This is because the contribution of terms of first order in the eccentricity to this operator leads to quantities of order $m e \vert v_{\lambda}\vert v_{\rm a}'/r$ in the linearized equations. For $m \sim r/H, \sim 1/e$ these are comparable to the dominant terms. All other eccentricity contributions are smaller by a factor at most ${\sim} e$.

We have to first order in eccentricity

$$v_{\lambda}\vert\nabla \lambda \vert = (GM_*/a^3)^{1/2}\left(... ... + \int \left({e\over r}{\rm d}r\right)\cos\varphi \right)\cdot$$ (36)

Thus if we define a new angle M through $M= \lambda - (2e+\int\left({e\over r}{\rm d}r\right)\sin\lambda,$the operator $\vert\nabla \lambda \vert v_{\lambda}{\partial \over \partial \lambda},$ becomes $=(GM_*/a^3)^{1/2}{\partial \over \partial M}$. The angle M is just the mean anomaly as can be seen from the fact that to first order in e,

$$r=a(1-e\cos M), \varphi= M+2e\sin M.$$ (37)

For the motion of a particular disc particle we have $M=\Omega(t-t_0),$ where t₀ denotes the time of periastron passage.

To lowest order in e, the coordinates (a,M) behave just like the cylindrical coordinates $(r,\varphi)$ in the sense that one may make the replacements $r \rightarrow a, \varphi \rightarrow \lambda$ in the standard linearized equations expressed in cylindrical coordinates. However, recall that in the forcing potential $\Phi_{\rm p}',$ we must make the replacement (37) which can be thought of evaluation on an eccentric disc orbit.

5.1 Protoplanet forcing potential

Neglecting the indirect term, which does not contribute to the Fourier components with large $m \sim r/H$of interest, the protoplanet forcing potential is

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

Here, the coordinates of the protoplanet are $(R,\nu)$ and b is the softening parameter which is introduced to take account of the finite vertical thickness of the disc. In general $b \sim H$ (e.g. Artymowicz 1993).

For an eccentric protoplanet orbit with semi-major axis $a_{\rm p}$ and eccentricity $e_{\rm p},$ we have

$$R=a_{\rm p}(1-e_{\rm p}\cos (\omega t )), \nu= \omega t +\varpi +2e_{\rm p}\sin(\omega t)$$ (39)

where the mean motion is $\omega ,$ we have taken pericentre passage to be at t=0, and $\varpi$ is the longtitude of pericentre which is not necessarily zero as we allow for the apsidal line of the orbit not to be lined up with that of the disc.

Using (37) and (39), we perform the Fourier decomposition

  $$\begin{eqnarray*}\Phi_{\rm p}' =\sum_{m,n,n_1} C_{m,n,n_1} \exp i\left((n-m)\omega t \right.\\ \left.+ m(M-\varpi)+n_1M\right). \end{eqnarray*}$$

Here, the sum is over positive values of mand both positive and negative values of n and n₁.The convention is that the real part is to be taken to give the physical potential. The coefficients C_m,n,n1are real and given by

$${4 \pi^3 C_{m,n,n_1} \over Gm_{\rm p} \omega} =$$

$$~~~-\int {\exp-i\left((n-m)\omega t - m\varpi+(n_1+m)M\right) \ov... ...vert r^2+R^2-2rR\cos(\varphi -\nu) +b^2\vert}} {\rm d}M {\rm d}\varpi {\rm d}t.$$ (40)

The integral is over a $2\pi$ cycle in each of the angles  $\varpi, M, \omega t.$Note that using (37) and (39) to express the disc and protoplanet coordinates makes C_m,n,n1 a function of $a,e,a_{\rm p},$ and $e_{\rm p}.$

An important aspect of the linearized problem expressed in terms of the coordinates a,M is that the separable coordinates are t and M.Accordingly a separable response occurs to a combination of terms in the Fourier decomposition with fixed (n₁+m) =k₁ and (n-m)=k₂.

Thus for such a particular forcing term with fixed (k₁,k₂)

$$\Phi_{\rm p}' = f_{k_1,k_2} \exp i\left(k_2\omega t +k_1 M\right),$$ (41)

where

$$f_{k_1,k_2} = \sum_{m,n,n_1} \delta_{n_1+m,k_1} \delta_{n-m,k_2} C_{m,n,n1} \exp(-im\varpi) .$$ (42)

Here $\delta$ denotes the Kronnecker delta.

A novel feature is that, unlike in the case of a forced axisymmetric disc, a separated forcing potential component in general depends on $\varpi.$This is because the angle between the apsidal lines of the disc and protoplanet orbits can have physical significance. But note that when both protoplanet and disc eccentricity are zero only one term can survive in (42) when m = k₁ = - k₂ and then $\varpi$ appears only as a redundant complex phase.

For the general forcing term (41) the pattern speed $\Omega_{\rm p} = -k_2\omega /k_1.$ We shall consider the situation when the disc surface density $\Sigma \propto r^{-3/2}$ as is the case in our models away from the boundaries. Then corotation resonances may be ignored and the main interaction is through wave excitation at the Lindblad resonances (e.g. Goldreich & Tremaine 1978). These occur when $k_1(\Omega -\Omega_{\rm p}) = \pm \Omega \sqrt{1+\xi^2}$ with $\xi = k_1c/(r\Omega)$ (e.g. Artymowicz 1993). For k₁ >0, the positive sign applies to the outer Lindblad resonance (OLR) and the negative sign to the inner Lindblad resonance (ILR). In both cases a wave is excited that propagates away from the protoplanet. The waves are associated with an outward energy and angular momentum flux. In the case of the ILR, the background rotates faster than the wave so that as it dissipates, energy and angular momentum are transferred to the protoplanet orbit. In the case of the OLR, the background rotates more slowly so that energy and angular momentum are removed from the protoplanet orbit.

Because the linearized response problem is formally identical to the one obtained for an axisymmetric disc $(a \rightarrow r,\; M \rightarrow \varphi)$, consistent with he other approximations we have made, we can evaluate the outward energy flow rate associated with outward propagating waves using the approximate expressions developed by Artymowicz (1993) and Ward (1997) which depend only on the disc state variables evaluated at the respective resonances. We comment that there are uncertainties associated with the use of these expressions in evaluating orbital evolution rates especially when they are derived by summing torque contributions with varying signs. However, cancellation effects do not appear to be significant (Artymowicz 1993; Ward 1997) and as they are supported by general considerations, we believe that the main features of the results derived to be correct. Following this procedure, the outward energy flow rate associated with outward propagating waves is given by

$${{\rm d} E_{\rm d}\over {\rm d}t} ={\pi^2 \Sigma r^2 \over 3\... ...i^2}(1+4\xi^2) }\left\vert {\cal{E}}_{k_1,k_2} \right\vert^2 ,$$ (43)

where

$${\cal{E}}_{k_1,k_2} = {{\rm d} f_{k_1,k_2}\over {\rm d}r}+ {2m^2(\Omega-\Omega_{\rm p})f_{k_1,k_2} \over \Omega r}\cdot$$

However, unlike in an axisymmetric disc, the outward angular momentum flow in the disc associated with the protoplanet, ${\rm d}J_{\rm d}/{\rm d}t \ne (1/\Omega_{\rm p}){\rm d}E_{\rm d}/{\rm d}t .$

One can find the outward angular momentum flow rate produced directly by the protoplanet associated with the linear response using

$${{\rm d} J_{\rm d}\over {\rm d}t} = \int \Sigma' {\partial \Phi_{\rm p}'\over \partial \varpi} r{\rm d}r{\rm d}\varphi,$$ (44)

the integral being taken over the disc. Recalculation of the torque formula then gives

$${{\rm d} J_{\rm d}\over {\rm d}t}=Re\left({\pi^2 \Sigma r^2 \... ...(1+4\xi^2) } {\cal{E}}_{k_1,k_2}{\cal{T}}^*_{k_1,k_2}\right) ,$$ (45)

where

$${\cal{T}}_{k_1,k_2} = {{\rm d} g_{k_1,k_2}\over {\rm d}r}+ {2m^2(\Omega-\Omega_{\rm p})g_{k_1,k_2} \over \Omega r} ,$$

with

$$g_{k_1,k_2}=\sum_{m,n,n_1} m\delta_{n_1+m,k_1} \delta_{n-m,k_2} C_{m,n,n1} \exp(-im\varpi) .$$ (46)

The rate of change of protoplanet eccentricity, $e_{\rm p}$, induced by each term can be found by using the orbit equation

$${e_{\rm p} \omega J\over (1-e_{\rm p}^2)^{3/2}}{ {\rm d}e_{\r... ...t} -\omega { {\rm d}J\over {\rm d}t} (1-e_{\rm p}^2)^{-1/2} ,$$ (47)

where E, and J are the orbital energy and angular momentum of the protoplanet and ${\rm d}J/{\rm d}t = \pm {\rm d}J_{\rm d}/{\rm d}t, {\rm d}E/{\rm d}t = \pm {\rm d}E_{\rm d}/{\rm d}t,$with the positive sign applying for an ILR and the negative sign for an OLR. The total rate of change is found by summing contributions from each resonance occurring for all values of  (k₁,k₂). The eccentricity evolution time-scale, which might be associated either with excitation, when it is positive, or damping, when it is negative, is then

$$t_{\rm e} = \frac{e_{\rm p}}{{\rm d}e_{\rm p}/{\rm d}t}\cdot$$ (48)

We define the time for the angular momentum to decay by a factor of e as

$$t_{\rm J} = - \frac{J}{{\rm d}J/{\rm d}t},$$ (49)

such that a positive value means inward migration.