5 A simple analytic model

In order to help understand some of the simulation results described above, we consider a simplified model system of primary, orbiting companion and disc. We suppose the companion orbits inside a disc with truncated inner edge. We make the additional simplification that the central star remains fixed at the origin of the coordinate system. The disc itself is taken to be inviscid with a simple relation between vertically integrated pressure, Pand surface density, $\Sigma.$Much of the complex nonlinear dynamics involving the balance between tides and viscosity determining the structure of the inner disc edge is thus omitted. However, enough content remains in order to display an instability leading to eccentricity growth in both disc and companion orbit, as well as the excitation of an outward propagating density wave with half the pattern speed. The existence of this has been described above.

5.1 Slowly varying modes with m = 1

We are concerned here with the situation when the companion orbit becomes eccentric. In general it will precess with a period very much longer than the orbital period. In an inertial frame the time averaged orbit then appears as a perturbation of the unperturbed circular orbit with azimuthal mode number m=1.

5.2 Secular perturbation theory

Because we are interested in modes that appear to change slowly in the inertial frame, we adopt the approach of secular perturbation theory. We thus perform a time average over an interval long compared to the orbital period but short compared to the orbital precession period. The companion then acts gravitationally as if it mass were continuously distributed along its orbit with line density inversely proportional to its speed. In an average sense it behaves like a continuous mass distribution and can be regarded as an infinitesimally thin annular extension of the disc. In this respect we suppose that we have a disc with well defined inner edge, interior to which there is a gap or hole in which the companion orbits. In the secular theory, we shall regard the companion as equivalent to an annular extension of the disc orbiting inside the gap and consider an eccentric companion orbit as a non axisymmetric perturbation of the initially axisymmetric combined disc-companion system. In this way an eccentric disc is produced also as comprising the normal mode. Good coupling results in eccentricities of comparable magnitudes in both components, while weak coupling enables them to exhibit eccentricities separately. The condition for good coupling is that the precession period induced in disc orbits by the companion and its own pressure and self-gravity be comparable to the precession period induced in the companion orbit by the disc. This is met when the disc mass, $M_{\rm d},$ contained in a radial scale comparable to the gap radius (or 1:3 resonance radius as taken below) is comparable to the mass of the companion. When these masses are very different, the situation will favour either an eccentric companion orbit or an eccentric disc.

This occurs because the slowly precessing nonaxisymmetric pattern in the disc makes it behave similarly to a planet of mass $M_{\rm d}.$From standard secular perturbation theory (Brouwer & Clemence 1961) eccentricity is freely exchanged between this planet and the companion when their orbital angular momenta are comparable, a situation that occurs here when the masses are comparable.

5.3 Negative angular momentum modes

The slowly varying m=1 modes described above are found to be negative angular momentum modes. This means that if they induce an angular momentum transfer to outer disc material from the region where they are located, they will grow unstably. Within the framework of secular perturbation theory of an initially axisymmetric disc this does not occur. However, if the high frequency components of the companion potential are now allowed to act, their coupling with the m=1 secular mode induces an angular momentum loss and hence growth. This occurs through the coupling of the m=1 low frequency mode with an m=1 potential component with pattern speed equal to the orbital frequency. This produces a forcing potential with m=2 and pattern speed one half the orbital frequency. As we shall see this can produce a density wave resonantly excited at the 1:3 resonance located at $r=2.08 r_{\rm p},$$r_{\rm p}$ being the companion semi-major axis, provided that resonance is contained within the disc, and hence mode growth. A favourable situation occurs when the gap or inner hole is large enough that possibly stabilizing corotation resonances are absent (Goldreich & Tremaine 1980; Artymowicz 1993; Lin & Papaloizou 1993a). When growth occurs we find two limiting cases. When the companion mass dominates that of the disc, an eccentric disc is produced. This situation is illustrated by our run in which the companion was maintained on a fixed circular orbit. In the opposite case, an eccentric companion orbit is expected as has been previously discussed by e.g. Artymowicz (1992) and Lin & Papaloizou (1993a). In the intermediate case eccentricities in both components are produced.

5.4 Slowly varying modes involving disc and companion

We consider linear perturbations of a two dimensional flat axisymmetric disc and time averaged planet. The dependence on time, t, and azimuthal angle, $\varphi,$ is taken to be through a factor $\exp{ i (m\varphi +\sigma t)}$ henceforth taken as read. Here the azimuthal mode number m=1 and the mode frequency $\sigma$ is such that $\vert\sigma\vert \ll \Omega,$ where $\Omega$ is the angular velocity in the unperturbed disc where the unperturbed surface density is $\Sigma.$ The Lagrangian displacement ${\rm\mbox{${\bf\xi}$ }} = (\xi_{\rm r}, \xi_{\varphi})$obeys the perturbed equations of motion

$$-(\sigma +\Omega)^2\xi_{\rm r} -2i\Omega(\sigma +\Omega)\xi_{... ...ver {\rm d}r}\xi_{\rm r} = -{\partial W \over \partial r}\cdot$$ (1)

$$-(\sigma +\Omega)^2\xi_{\varphi} +2i\Omega(\sigma +\Omega)\xi_{\rm r} = -{ i W \over r}\cdot$$ (2)

Equations (1) and (2) can be combined to give

$$\left(-(\sigma +\Omega)^2 + \kappa^2)\right)\xi_{\rm r} =-{\p... ...W \over \partial r}-{ 2\Omega W \over r (\sigma +\Omega)}\cdot$$ (3)

Here $W= \Sigma' c^2/\Sigma +\Psi',$ where $\Sigma'$ and $\Psi'$ are the surface density and gravitational potential perturbations respectively. The vertically integrated pressure P is assumed to be a function of $\Sigma,$ $c^2 = {\rm d}P/{\rm d}\Sigma$ is the square of the sound speed, and $\kappa^2= 2\Omega r^{-1} {\rm d}(r^2 \Omega) /{\rm d}r$ is the square of the epicyclic frequency. The perturbation to the gravitational potential may be written (temporarily restoring the $\varphi$ dependence)

$$\Psi' = G\int{\nabla\cdot\left(\Sigma({\bf r'}) {\rm\mbox{${\... ...\vert{\bf r} -{\bf r'}\right\vert} r'{\rm d}r'{\rm d}\varphi',$$ (4)

where this and similar integrals are taken over the active mass distribution. We may apply (4) to both the distributed disc mass and the separated, time averaged, but initially radially localized companion mass for which

$$\Sigma ={ M_{\rm p} \delta( r -r_{\rm p}) \over 2\pi r_{\rm p}}\cdot$$ (5)

Here $\delta$ denotes a delta function and $M_{\rm p}$ the companion mass. Before so doing, we use the fact that $\vert\sigma\vert/\Omega$ is small to make some simplifications. In this limit and for a thin disc for which forces due to internal pressure and self-gravity are small, Eq. (2) becomes

$$\xi_{\varphi} =2i\xi_{\rm r},$$ (6)

which is the limiting form for purely epicyclic motion. We use this to evaluate the pressure and self-gravity terms on the right hand side of (3). The error in doing this, for a gravitationally stable disc, is measured by the greater of $\vert\sigma\vert/\Omega$ and H/r, where H is the semi-thickness of the disc. Then we have

$$\Sigma' = -r{{\rm d} (r^{-1} \Sigma \xi_{\rm r})\over {\rm d}r}\cdot$$ (7)

Equation (4) can then be written

$$\Psi' = G\int\left(r{{\rm d} (r^{-1} \Sigma \xi_{\rm r})\over {\rm d}r}\right)_{r=r'} K(r,r') r'{\rm d}r',$$ (8)

with

$$K_m(r,r') = \int^{2\pi}_0 {\cos(m\varphi)\over \sqrt{(r^2+r'^2 -2rr'\cos(\varphi))}}{\rm d}\varphi,$$ (9)

and noting that for m=1 we shall drop the subscript from K_m(r,r'),thus $K_1(r,r') \equiv K(r,r').$ We finally obtain a normal mode equation for $\xi_{\rm r}$ by using (8), and (7) in (3) in the form

 

$$2\left(\sigma + \omega_{\rm p} \right) \Omega r \xi_{\rm r} =-{1\... ...{r^3 c^2\over \Sigma}{{\rm d} (r^{-1} \Sigma \xi_{\rm r})\over {\rm d}r}\right)$$

$$+ {1\over r}{{\rm d}\left( r^2\Psi'\right) \over {\rm d} r}\cdot$$ (10)

Here we have used the fact that $\vert\sigma\vert/\Omega \ll 1,$and the precession frequency $\omega_{\rm p} = \Omega -\kappa.$

5.5 Disc and companion system

Although Eq. (10) applies to the joint disc and companion system, we shall separate the contributions to the potential from the companion and disc, writing

$$\Psi'= \Psi'_{\rm d} + \Psi'_{\rm p},$$ (11)

where $\Psi'_{\rm d} ,$ and $\Psi'_{\rm p}$ are the contributions from the disc and planet obtained from the appropriate surface density distributions respectively. Using (5), we find

$$\Psi'_{\rm p}(r) = - {G M_{\rm p} \xi_{\rm r}(r_{\rm p}) \ove... ... \partial r_{\rm p}}\left( K(r,r_{\rm p}) r_{\rm p}^2 \right).$$ (12)

Here the eccentricity of the companion orbit is related to the displacement by $e = \xi_{\rm r}(r_{\rm p})/r_{\rm p}.$

5.6 Disc modes

If we replace $\Psi'$ in (10) by the disc contribution $\Psi'_{\rm d}$ only, we obtain an equation for the normal modes of the disc only. In fact under reasonable boundary conditions that $\xi_{\rm r}$ is well behaved at the inner boundary and vanishes at large distances, (10) gives a self-adjoint eigenvalue problem for the eigenvalue $\sigma.$ The right hand side gives a self-adjoint operator with weight $\Sigma.$ Denoting, the eigenvalues by $\sigma_j,$and corresponding eigenfunctions, $\xi_{\rm r}$ by u_j, j= 1,2,...,we have the orthogonality condition

$$\int \Sigma \Omega r u_k^* u_j {\rm d}r = N_j \delta_{kj},$$ (13)

with this and later integrals being taken over the disc. We comment that the eigenvalues are all real and the local dispersion relation associated with (10) is the low frequency limit of the well known relation for spiral density waves (Lin & Shu 1964)

$$2\left(\sigma + \omega_{\rm p} \right) \Omega = -2\pi G \Sigma \vert k\vert + c^2k^2,$$ (14)

where k is the radial wavenumber and of course m=1.Each mode propagates with a prograde pattern speed $-\sigma_j,$which is also the precession frequency.

5.7 The effect of the companion

The effect of the companion may be incorporated by adding $\Psi'_{\rm p}$into $\Psi'$ in (10). Regarding the companion potential so added as giving an external forcing term, Eq. (10) may be solved using the standard method of eigenfunction expansion to give interior to the disc

$$\xi_{\rm r} = - \sum^{\infty}_{j=1} {u_j\over 2 N_j (\sigma -... ...\rm p} {{\rm d} (r^{-1} \Sigma u^*_j)\over {\rm d}r} {\rm d}r.$$ (15)

5.8 Combined modes

Omitting the pressure term, Eq. (10) gives for the equation of motion of the companion

$$2\left[\sigma + \omega_{\rm p}(r_{\rm p}) \right] \Omega(r_{\... ...si'_{\rm d}(r) \right) \over {\rm d} r} \right]_{r=r_{\rm p}}$$ (16)

with

$$\Psi_{\rm d}'(r) = -G\int \frac{\Sigma (r') \xi_{\rm r} (r')}... ...ial \over \partial r'}\left(K(r,r') (r')^{2}\right) {\rm d}r'.$$ (17)

Combining (16), (15) and (17) gives an equation for the eigenvalue $\sigma$ for the joint normal mode in the form

$$\left[\sigma + \omega_{\rm p}(r_{\rm p}) \right] \Omega(r_{\r... ...t^2 M_{\rm p} \over 8\pi r_{\rm p}^2 N_j (\sigma - \sigma_j)},$$ (18)

with the coupling coefficient given by

 

$$C_j = G\int (r' r_{\rm p} )^{-1} \Sigma (r') u_j (r') {\partial^2 \over \partial r' r_{\rm p} }$$

$$\times\left(K(r_{\rm p},r') (r' r_{\rm p} )^{2}\right) {\rm d}r'.$$ (19)

The schematic behaviour of the eigenvalues is that, when the coupling coefficient C_j is negligible, the solutions are $\sigma= \sigma_j, j=1,2. ..$and $\sigma = -\omega_{\rm p}(r_{\rm p}).$ In this case the companion and disc are decoupled and precess independently. In contrast when the coupling is non negligible the precession of companion orbit and disc are linked. One induces eccentricity in the other. To estimate when this occurs, one needs that $(\vert C_j\vert^2 M_{\rm p} )/( 8\pi r_{\rm p}^2 N_j \Omega)$ be of the same order as the square of the difference in disc and companion orbit precession frequencies $(\omega_{\rm p} + \sigma_j)^2.$If the companion and disc mutually induce eccentricity and precession, simple estimates, assuming there is a single radial scale length, indicate this requires that $M_{\rm d} \sim M_{\rm p}$ for good coupling.

We comment that although the gravity of the disc acting on both itself and the companion has been included in the above formalism, the self-interaction only affects the precise determination of the $\sigma_j$ and u_j.The above discussion still applies if that is neglected, as long as the interaction with the companion is retained. This is the case for the simulations presented here.

5.9 Mode angular momentum

The effective angular momentum content in the joint mode can be evaluated using standard methods (e.g. Goodman & Ryu 1992; Lin & Papaloizou 1993b). For the low frequency modes considered here, we obtain after a straightforward calculation

$${\cal J} = -{1\over 2}\int \Sigma \vert\xi_{\rm r}\vert^2\Omega r{\rm d}r{\rm d}\varphi ,$$ (20)

where the integral is over the entire system. We comment that, using (5), the contribution from the companion to ${\cal J}$ is ${\cal J}_{\rm p} = - 1/2 M_{\rm p} e^2\sqrt{(GM_*)r_{\rm p}}.$This is negative corresponding to a negative angular momentum mode and it is also minus the radial action. If angular momentum is drained from the mode, the companion orbit eccentricity and radial action grow. In general the total angular momentum content is

$${\cal J} = -{1\over 2} M_{\rm p} e^2\sqrt{(GM_*)r_{\rm p}} -... ... \Sigma \vert\xi_{\rm r}\vert^2\Omega r{\rm d}r{\rm d}\varphi ,$$ (21)

where the second contribution now comes from the disc only, involving the radial component of the Lagrangian displacement $\xi_{\rm r} \equiv er$ there.

  
5.10 Angular momentum loss due to resonant torques

The modes discussed above, in the time averaging approximation, do not lose angular momentum. However, only the time averaged potential due to the companion has been included. When the effects of the full potential are included, additional perturbation of the disc occurs. This can lead to angular momentum loss through resonant torques (Goldreich & Tremaine 1978). In that case, the mode eccentricities grow and instability occurs. The full potential due to the companion located at $(r_{\rm p}, \varphi_{\rm p})$ can be written

$$\Psi_{\rm p}' = -{GM_{\rm p}\over \pi} \sum_{m=0}^{\infty}{K_... ...m p}) \cos m(\varphi -\varphi_{\rm p})\over ( 1+ \delta_{m0})}.$$ (22)

When the disc and companion orbit are eccentric, with small eccentricities, to evaluate the companion potential at the location of a perturbed fluid element in the disc, we set $r \rightarrow r +\vert\xi_{\rm r}\vert\cos(\varphi), \varphi \rightarrow \varp... ...{\rm p}\cos(\omega t), \varphi_{\rm p} \rightarrow \omega t -2e\sin(\omega t) .$Here, arbitrary phases are chosen such that all ellipses have apsidal lines aligned at $\varphi=t=0,$consistently with the joint mode. The orbital frequency is $\omega$and we have neglected the small precession frequency. The simulations suggest an alignment such the eccentricities of the disc and orbit are of opposite sign. In this case the disc pericentre is closest to the orbital apocentre. The contribution to the companion potential that is first order in the eccentricities $\Psi_{\rm pe}'$ is given by

 

$${\Psi_{\rm pe}'\over GM_{\rm p}} = \sum_{m=0}^{\infty}{\cos m\beta \over \pi (1+\delta_{m0})}$$ -

$$\times \left( \vert\xi_{\rm r}\vert\cos\varphi {\partial \over \partial ... ..._{\rm p}\cos\omega t{\partial \over \partial r_{\rm p}}\right) K_m(r,r_{\rm p})$$

$$\sum_{m=1}^{\infty} mK_m(r,r_{\rm p}) {\sin m\beta \over \pi}$$ +

$$\times \left[-\vert\xi_{\varphi}\vert\sin\varphi/r +2e\sin(\omega t) \right],$$ (23)

where $\beta = \varphi -\omega t.$ We comment that for zero eccentricities, all time dependent components of the companion potential have pattern speed $\omega.$ In the outer disc, there are possible outer Lindblad resonances where $\Omega =m\omega/(m+1), m= 1,2...$ However, as the inner edge approaches the 1:2 resonance, only the resonance with m=1 remains eventually. When eccentricities are included, contributions occur in (23) with pattern speeds $(m -1) \omega/m, m = 2, 3, ...$ The outer Lindblad resonances associated with these are at $\Omega = (m - 1) \omega/(m+1).$ For m=2, this gives the smallest value of $\Omega = \omega/3$ which occurs at the 1:3 resonance. The corresponding pattern speed of the resonant forcing is $\omega/2.$ As the material at the outer Lindblad resonance rotates more slowly than the pattern speed, resonant torques cause angular momentum loss leading to eccentricity growth of the joint disc-companion system. For discs with inner edges approaching the 1:2 resonance such that the disc surface density becomes depressed there, potentially damping corotation torques (Goldreich & Tremaine 1980) and coorbital Lindblad torques (Artymowicz 1993) are weakened. Here we shall assume that only the 1:3 resonance needs to be considered. The resonant forcing occurs through terms produced by a nonlinear coupling between the secular m=1 mode and disturbances with m=1moving with the orbital pattern speed $\omega.$ We begin by calculating the resonant forcing using (23) as a forcing potential on the unperturbed axisymmetric disc. However, it is important to note that not all the effective forcing terms are included in this way. They arise also by coupling between the secular m=1 mode and the tidal distortion of the disc produced by the circular orbit tide propagating with a pattern speed $\omega.$Evaluation of these is lengthy and depends on the detailed disc model which has here been simplified. Instead of such evaluation, we replace the disc eccentricity in (23) by ${\overline{e_{\rm d}}}$to indicate additional forcing. Although such forcing is not necessarily through a potential, it can still be taken into account by defining a suitable ${\overline{e_{\rm d}}}$ which is proportional to the disc eccentricity amplitude (see Eq. [24]). This procedure will enable eccentricity growth rates to be qualitatively discussed and estimated. The component of (23) with pattern speed $\omega /2,$ is then given by $\Psi_{\rm pe2}' (r,r_{\rm p})\cos(2\varphi -\omega t ),$where

 

$${ \pi \Psi_{\rm pe2}'\over GM_{\rm p}}= - {\overline{e_{\rm d}}}r... ...} {\partial K(r,r_{\rm p}) \over \partial r } -{K(r,r_{\rm p}) \over r }\right)$$

$$- er_{\rm p}\left({1\over 2}{\partial K_2(r,r_{\rm p}) \over \partial r_{\rm p}} +{2 K_2 (r,r_{\rm p}) \over r_{\rm p} } \right).$$ (24)

Here we have set $\vert\xi_{\rm r}\vert= \vert\xi_{\varphi}\vert/2 = {\overline{e_{\rm d}}} r.$ We calculate the rate of loss of angular momentum as a result of resonant torques using the torque formula of Goldreich & Tremaine (1978, 1980) to obtain

$${{\rm d} {\cal J} \over {\rm d}t} = -{ 2\pi^2 r^2\Sigma\over ... ...2}' \over \partial r} - {4 \Psi_{\rm pe2}'\over r} \right )^2,$$ (25)

where r and the disc quantities are evaluated at the 1:3 resonance.

The angular momentum is transported outwards via a wave with pattern speed $\omega/2.$ It has to be extracted from the m=1 joint mode and the companion orbit which have combined together to give the resonant forcing.

Extraction from an eccentric companion orbit can occur directly through the excited m=2 wave. Extraction from the m=1 mode can occur through a recoupling of the m=2 wave with the tidal potential with pattern speed $\omega,$to produce a time independent forcing with m=1. Similarly extraction from the orbit can occur through recoupling of the m=2wave to the m=1 joint mode to produce an m=1 disturbance with the circular orbit pattern speed $\omega.$

If the resonantly excited m=2 wave carries an angular momentum $\Delta J$, the associated energy is $\omega\Delta J/2.$If this energy is supplied by the circular orbit tide, then the angular momentum supplied along with it will be $\Delta J/2,$leaving $\Delta J/2$ to be extracted from the joint m=1 mode.

The above discussion suggests the growth rate of this mode can be obtained by assuming half the resonantly induced angular momentum loss is extracted from it.

The growth rate of the combined mode is then estimated as

$$\gamma = {1\over 4{\cal J}}{{\rm d} {\cal J} \over {\rm d}t},$$ (26)

with

$${\cal J}= -{1\over 2} M_{\rm p} e^2\sqrt{(GM_*)r_{\rm p}} -{... ...r 2}\int \Sigma e_{\rm d}^2 r^3\Omega {\rm d}r{\rm d}\varphi .$$ (27)

Evaluating (25) we obtain

 

$${\gamma\over \omega} {M_{\rm d} M_{\rm p}\over M_{*}^2}\left({r_{\rm p}\over r}\right)^8$$ =

$$\times {9\pi \left[(4{\overline{e_{\rm d}}} - {2r\over3} {{\rm d}{\overl... ...\pi \Sigma e_{\rm d}^2 r^3\Omega {\rm d}r/(M_{\rm p}\omega r_{\rm p}^2)\right]}$$ (28)

where $M_{\rm d} =\pi\Sigma r^2.$

Note that here we neglect any change in the companion semi-major axis due to the time dependent terms in the disc potential. To set the magnitudes of typical growth rates, we first set the terms involving the disc eccentricity $e_{\rm d}, {\overline{e_{\rm d}}}$to zero. We find from (28) that for the 1:3 resonance with $r=2.08 r_{\rm p},$ $\gamma^{-1} = M_*^2/(19 M_{\rm d}M_{\rm p})P_{0},$ with P₀ being the orbital period. Thus for $M_{\rm p}= 2M_{\rm d} = 0.03M_*,$ here we allow a factor of two surface density enhancement above the initial value at the 1:3 resonance, we find a growth time $\sim$120 companion orbital periods. This is comparable, but a factor of $\sim$3 smaller, than the eccentricity growth rate obtained for calculation N5, and is consistent with the idea that the orbital eccentricity driving arises through effects due to both the orbital eccentricity itself and the eccentricity induced in the disc, the latter actually giving stronger effects. We comment that for these parameters the disc and companion orbit eccentricities are likely to be coupled. This is because typically $M_{\rm d} \sim M_{\rm p}.$ Then the precession frequency induced by the companion in externally orbiting disc matter is given by $\omega_{\rm p} = 0.75 \Omega(r) (r_{\rm p}/r)^2 (M_{\rm p}/M_*).$Similarly the precession frequency induced by the disc (assumed to have constant surface density) in the companion orbit is $\omega_{\rm p} = \Omega(r_{\rm p}) (r_{\rm p}/r)^3 (M_{\rm d}/M_*).$These are comparable leading to the likely setting up of a joint mode. We comment that this condition is the same as requiring that the companion orbital angular momentum and the orbital angular momentum contained in $M_{\rm d}$ assumed in circular orbit be comparable.

Under the conditions of the simulations considered here, $M_{\rm d} \sim M_{\rm p}.$Then it is also likely that the radial action in the disc is comparable to that of the companion and so modifies the growth rate. However, as long as they are comparable our estimate is probably reasonable. Note that $\gamma^{-1} \propto M_*^2/(M_{\rm d}M_{\rm p})P_{0}$as long as $M_{\rm p}$ is comparable to or less than $M_{\rm d}.$If the companion mass dominates, so fixing it in a circular orbit, then e=0, in (28) and $\gamma^{-1} \propto M_*^2/(M_{\rm p}^2)P_{0 },$attaining a limiting value independent of the disc mass. Our simulation with fixed companion orbit indicates that the growth rate does not change very much as $M_{\rm d}$ is reduced significantly below $M_{\rm p}.$