Disc–planet eccentricities

Growth of eccentric modes in disc–planet interactions

Jean Teyssandier and Gordon I. Ogilvie
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road,
Cambridge CB3 0WA, United Kingdom
E-mail: jt591@cam.ac.uk
Abstract

We formulate a set of linear equations that describe the behaviour of small eccentricities in a protoplanetary system consisting of a gaseous disc and a planet. Eccentricity propagates through the disc by means of pressure and self-gravity, and is exchanged with the planet via secular interactions. Excitation and damping of eccentricity can occur through Lindblad and corotation resonances, as well as viscosity. We compute normal modes of the coupled disc–planet system in the case of short-period giant planets orbiting inside an inner cavity, possibly carved by the stellar magnetosphere. Three-dimensional effects allow for a mode to be trapped in the inner parts of the disc. This mode can easily grow within the disc’s lifetime. An eccentric mode dominated by the planet can also grow, although less rapidly. We compute the structure and growth rates of these modes and their dependence on the assumed properties of the disc.

keywords:
celestial mechanics – accretion, accretion discs – hydrodynamics – planet-disc interactions – planetary systems: protoplanetary discs

1 Introduction

A significant fraction of extrasolar planets that have so far been discovered are known to have large orbital eccentricities. In particular, radial-velocity surveys have shown that giant planets far enough from their host stars not to suffer strong tidal circularization have eccentricities spanning the whole range from zero to near-unity (see, e.g., Marcy et al., 2005; Butler et al., 2006). Kane et al. (2012) confirmed that this feature was still present in the transit survey obtained with the Kepler mission, and argued that the mean eccentricity increases with planet size, although Plavchan et al. (2014) cautioned that uncertainties in the measurement of stellar radii in the Kepler sample can affect this result. Theories of the formation and evolution of planetary systems have yet to explain this wide distribution of eccentricities.

It is likely that mutual gravitational interactions between planets in a gas-free environment can cause their eccentricities to grow to near-unity values. This can be the result of a violent scattering (Rasio & Ford, 1996; Weidenschilling & Marzari, 1996; Chatterjee et al., 2008; Jurić & Tremaine, 2008) or dynamical relaxation (Papaloizou & Terquem, 2001; Adams & Laughlin, 2003) of a large population of planets. Kozai-Lidov oscillations, caused by a stellar companion (Wu & Murray, 2003; Fabrycky & Tremaine, 2007; Naoz et al., 2012), a planet (Naoz et al., 2011), or a disc (Terquem & Ajmia, 2010; Teyssandier et al., 2013a), can also pump up eccentricities. Finally, chaotic diffusion of angular momentum deficit can increase eccentricities and inclinations over very long timescales, a process known as secular chaos (Wu & Lithwick, 2011). However, the initial conditions required to trigger these mechanisms often require the planets to have a certain amount of eccentricity to start with. It is unclear whether this amount is realistic or not, as it is currently difficult to predict planetary eccentricities at the beginning of the gas-free stage of dynamical evolution.

Another mechanism by which planets could gain eccentricities is their gravitational interaction with the protoplanetary disc in which they formed, in particular through mean-motion resonances between the planet and the disc. Since the resonant interactions involve an exchange of energy and angular momentum between the planet and disc, as well as a dissipative loss of energy from the system, they can affect the orbital eccentricities of both the planet and the disc.

It is now widely accepted that low-mass planets with small eccentricities that are unable to open a gap in a circular disc are expected to experience a relatively rapid eccentricity damping. Ward (1988) showed that material orbiting at the same radius as the planet exerts a strong torque via coorbital resonances. This torque overcomes that of the non-coorbital resonances discussed by Goldreich & Tremaine (1980) and results in eccentricity damping. This result holds as long as the surface density in the coorbital region is not strongly depleted. A refined analytical treatment of this problem was carried out by Artymowicz (1993), and a similar result was obtained semi-analytically by Tanaka & Ward (2004). Numerical simulations of eccentric low-mass planets embedded in 2D and 3D discs (see, e.g., Cresswell et al., 2007; Bitsch & Kley, 2010) have confirmed these findings.

One the other hand, Goldreich & Tremaine (1980) considered a satellite that is massive enough to open a gap in a circular disc (the same applies to gap-opening giant planets in protoplanetary discs), and showed that eccentric Lindblad resonances (ELRs) act to increase the eccentricity of the satellite, while the net effect of eccentric corotation resonances (ECRs) is to decrease the eccentricity. They concluded that, to lowest order in eccentricity, ECRs dominate by a small margin, and prevent eccentricity growth. An important breakthrough was achieved by Goldreich & Sari (2003) and Ogilvie & Lubow (2003) in their study of the saturation of the corotation torque (see also Goldreich & Tremaine, 1981). Goldreich & Sari (2003) noted that, unlike Lindblad resonances, corotation resonances readily undergo a nonlinear saturation. Thus there is a finite-amplitude instability which can occur if the eccentricity exceeds a small critical value. Ogilvie & Lubow (2003) derived an expression for the level of saturation of the corotation torque in a three-dimensional gaseous disc.

The above description makes the significant assumption that the disc remains circular. However, a general Keplerian disc involves nested elliptical orbits, and the secular interaction between the planet and the disc inevitably leads them to exchange eccentricity. Excitation of the eccentricity of a disc has been found to be important in various astrophysical contexts, such as in the rings of Uranus (Goldreich & Tremaine, 1981), or in superhump binary stars (Whitehurst, 1988; Lubow, 1991). A more general description is therefore desirable, in which both the planet and the disc are allowed to have eccentric orbits. Theories of eccentric gaseous discs have been discussed in the past by, e.g., Kato (1983) and Ogilvie (2001). Moreover, Ogilvie & Barker (2014) recently developed a local model of eccentric discs, in which they explored non-linear vertical oscillations and deduced the linear global evolution of eccentricity, and then studied the hydrodynamical instability of such discs in Barker & Ogilvie (2014).

Away from the coorbital region, the disc–planet interaction can be regarded as twofold. First, for the long-term orbit-averaged (secular) interaction, an analogy with planetary dynamics can be made. In celestial mechanics, the secular interaction between planets can be described by an eigenvalue problem, usually known as the Laplace–Lagrange theory (Murray & Dermott, 1999). In this theory, a reversible exchange of eccentricity occurs on long timescales (compared to orbital periods), and the total angular momentum deficit (AMD, a positive-definite measure of the eccentricity of the coupled system) of the system is conserved. (At this level of approximation, the inclinations undergo a similar, although decoupled, evolution.) A similar process can occur in disc–planet interactions, where the disc can be regarded as a continuum of rings of matter interacting with the planet, and with each other. As in the case of a planetary system, AMD is conserved in gravitational secular disc–planet interactions. We remark here that for a planet on an initially eccentric orbit interacting with a circular disc, secular exchange of angular momentum would cause the eccentricity of the planet to decay at first, while that of the disc increases.

The second part of the interaction occurs through mean-motion resonances, leading to an irreversible evolution in which the AMD can either increase or decrease. A theory of mean-motion resonances allowing for the eccentricities and inclinations of both the satellite and the disc to vary was formulated by Ogilvie (2007). Previously, a linear theory of disc–planet interactions was formulated by Lubow & Ogilvie (2001) in the case of inclination, embodying secular and resonant contributions. Unlike previous work (see, e.g., Borderies et al., 1984), the disc was allowed to develop a warped shape, which is likely if the planet and the disc have comparable angular momenta. The linear theory allows for a set of normal modes which can grow or decay according to the balance between resonances (which mostly increase the AMD) and viscous dissipation of the warp (which decreases it). Unlike Borderies et al. (1984), Lubow & Ogilvie (2001) found that, after allowing the disc to become warped, inclination growth was suppressed by viscous damping for typical estimates of viscosity and disc parameters (although Borderies et al. (1984) were concerned with interaction between planetary rings and moons, while Lubow & Ogilvie (2001) were studying protoplanetary discs).

The results of Lubow & Ogilvie (2001) regarding inclination highlight the fact that the secular evolution of the disc and that of the planet(s) are strongly coupled, and should be treated together rather than separately. The same statement can be made for eccentricities: in a system consisting of a disc and one or more planets, having comparable angular momenta, eccentricity is freely exchanged between the various components on secular timescales. To address the origin of eccentricity in protoplanetary systems, it is therefore more appropriate to ask under what conditions the total AMD of the coupled system, rather than the eccentricity of a planet, can grow.

Several numerical simulations of disc–planet interactions (or related systems) have investigated the possibility of eccentricity growth. Artymowicz et al. (1991) found, using smoothed particle hydrodynamics (SPH), that the eccentricity of a binary star can increase through its interaction with a circumbinary disc. In this case the gap is so wide that only the 1:3 ELR is present (and no ECR), so that the damping argument of Goldreich & Tremaine (1980) is not applicable. The mechanism is analogous to that operating in superhump binary stars, where the eccentricity of the circumstellar disc is excited and only the 3:1 resonance is operative. Such excitation was readily observed in grid-based simulations by Kley et al. (2008), with rapid growth rates. In the case of disc–planet interactions, Papaloizou et al. (2001) showed that massive planets (e.g., those of 20 Jupiter masses) could experience eccentricity growth, with the exterior disc also becoming noticeably eccentric. Kley & Dirksen (2006) found that, even if the planet is held on a fixed circular orbit, the disc can become eccentric for planets exceeding about 3 Jupiter masses (for typical disc parameters). Without fixing the planet’s orbit, d’Angelo et al. (2006) found growing eccentricities for planetary masses as small as 1 Jupiter mass. This extension of the work by Artymowicz et al. (1991) and Papaloizou et al. (2001) was made possible through the use of high-resolution two-dimensional numerical simulations over several thousand planetary orbits. In contrast with circumbinary discs, they found that the 1:3 ELR was highly ineffective at pumping the eccentricity, and that the 2:4 and 3:5 ELRs were the most efficient ones, with additional contributions arising from resonances located within the gap. Cresswell et al. (2007) studied the evolution of eccentricity and inclination of non-gap-opening planets, and confirmed that their eccentricities (and inclinations) would be damped by the disc, in good agreement with the semi-analytical work by Tanaka & Ward (2004). This result was later confirmed by Bitsch & Kley (2010) for various planet masses, up to 1 Jupiter mass, in three-dimensional isothermal and fully radiative discs. Perhaps most relevant to the present work is the study by Rice et al. (2008) of a hot Jupiter orbiting inside the magnetospheric cavity of the star. The authors found eccentricity growth for massive planets, which would eventually lead to the destruction of the planet, and possibly explain the paucity of very massive planets on short-period orbits. We will discuss these claims in the present paper. More recently, Dunhill et al. (2013) have claimed, using SPH simulations, that only very massive planets (e.g., those of Jupiter masses) embedded in discs with high surface densities could experience eccentricity growth, although their simulations were limited to a few hundred orbital periods. One the other hand, Duffell & Chiang (2015) recently argued that disc–planet interactions can indeed increase the eccentricity of a Jupiter-mass planet embedded in an isothermal disc. Such growth is possible when the planet opens a deep gap, and only if its initial eccentricity exceeds a threshold value, similarly to what was predicted analytically by Goldreich & Sari (2003). Finally, Tsang et al. (2014) argued in favour of eccentricity growth for gap-opening planets in non-barotropic discs. They showed that when stellar illumination heats the gap to a certain threshold, the corotation torque can be reduced to a point where it is not able to damp eccentricity.

In this paper we present a linear theory of eccentricity evolution in disc–planet systems. Our theory is analogous to the one established for inclination by Lubow & Ogilvie (2001). It is also related to the work by Papaloizou (2002) but differs by including a new treatment of mean-motion resonances between planets and discs (Ogilvie, 2007) and a three-dimensional description of the dynamics of eccentric discs (Ogilvie, 2001, 2008). We aim to compare our results with previous analytical studies (e.g., Goldreich & Sari, 2003) as well as numerical simulations such as that of Rice et al. (2008). The equations describing the behaviour of eccentricity are presented in Section 2. The associated conservation laws and equations for the angular momentum deficit of the system are obtained in Section 3, together with integral relations for the precession rates and the growth rates of normal modes. An analogy with the Schrödinger equation and an application to a numerical solution for a non-self-gravitating disc are presented in Section 5, while a more complete study that includes additional physics is carried out in Section 6. In Section 7 we conduct a numerical exploration of the influence of various physical parameters. In Section 8 we discuss our results, and we conclude in Section 9.

2 Evolutionary equations for an eccentric disc–planet system

In this section we present a set of linear equations that describe the evolution of small eccentricities in a disc–planet system. To a first approximation, fluid elements in an eccentric disc follow elliptical Keplerian orbits. The eccentricity and argument of pericenter can vary smoothly with the orbital semi-major axis . For small gradients in and , the orbits are nested without intersections (Ogilvie, 2001).

Throughout this paper the results are presented in terms of the complex eccentricity . Since and , this variable conveniently describes both the shape and the orientation of the elliptical orbits. With similar notations, the planet has a complex eccentricity . We also denote by , and the semi-major axis, orbital frequency and mass of the planet, with the planet-to-star mass ratio. In addition we denote by the total mass of the disc and .

2.1 Pressure

Goodchild & Ogilvie (2006) derived a linear partial differential equation governing the propagation of a small complex eccentricity in a two-dimensional, non-self-gravitating, adiabatic disc. It has the form of a dispersive wave equation related to the Schrödinger equation. In Appendix A.1 we derive a similar equation for a three-dimensional disc in which the isothermal sound speed is a specified function of radius, as is commonly assumed in numerical simulations of disc–planet systems. A similar method can be used to derive an equation for a 3D adiabatic disc (see also Ogilvie & Barker, 2014). The different models take the following forms:

  • 2D adiabatic model:

    (1)
  • 3D adiabatic model:

    (2)
  • 2D locally isothermal model:

    (3)
  • 3D locally isothermal model:

    (4)

Here is the Keplerian angular velocity around a star of mass , although our equations do take into account the small departure from Keplerian rotation arising from the radial pressure gradient. In addition, is the surface density, is the vertically integrated pressure, and is the adiabatic index. (In the 2D model, and are simply the density and pressure that appear in the two-dimensional fluid-dynamical equations.) Note that in a locally isothermal disc. There is an equivalence between the adiabatic and isothermal models when setting and assuming that is independent of .

Differences between the 2D and 3D models arise from the periodic variation of the vertical gravitational force exerted on a fluid element as it goes along its elliptical orbit. The vertical variation of the radial gravitational force also makes a difference. These effects account for the last term in each of the 3D equations, as well as for the different coefficients in the adiabatic case. In Eq. (2.1) the last term derives from 3D effects, while the preceding ‘non-adiabatic’ term results from the variation of the isothermal sound speed with radius. The very significant role of the 3D term for eccentric discs around Be stars has already been demonstrated by Ogilvie (2008).

2.2 Self-gravity and secular interactions with a companion

The self-gravity of the disc can be in treated in a secular way as a continuum version of the classical Laplace–Lagrange theory (see, e.g., Murray & Dermott, 1999). In Appendix A.2, we show that self-gravity gives rise to the following contribution:

(5)

where

(6)

is a symmetric kernel representing the strength of the gravitational interaction between slightly elliptical rings at radii and . In equation (2.2) the term represents the contribution to the precession of the annulus at arising from the axisymmetric gravitational potential generated by the annulus at , while represents the forcing of the eccentricity of the first annulus by that of the second. We resolve the singularity at by replacing with a softened symmetric kernel:

(7)

where is a dimensionless softening parameter, is a Laplace coefficient,

(8)

and is the smaller solution () of the quadratic equation

(9)

The softening length is then , and, in the limit , . We set to be the scaling factor which appears in the expression for the disc aspect ratio (see the description of the disc model in Section 4.3). This gives an approximation of the averaging of the gravitational potential of the disc over its vertical extent.

Self-gravity is therefore an additional means of propagation of eccentricity, as it gives rise to a global coupling of eccentricity between different parts of the disc. The role of self-gravity in eccentric discs was already discussed by Tremaine (2001).

The secular potential of a planet can be regarded as that of an elliptical ring whose mass is that of the planet. Hence its contribution to the eccentricity evolution is not fundamentally different from that of self-gravity:

(10)

Likewise, the contribution to the eccentricity evolution of the planet is

(11)

We note that in the linear theory of Lubow & Ogilvie (2001), the gravitational couplings in disc–planet systems lead to identical equations for the complex inclination, except for a change of sign and the replacement of with . It is straightforward to generalize the equations to the case of a multiple-planet system.

2.3 Viscosity

The mechanism by which angular momentum is transported in protoplanetary discs remains poorly understood. Even if this is modelled as an effective shear viscosity, its influence on the eccentricity is complicated and can even lead to growth rather than decay in some circumstances (Latter & Ogilvie, 2006). In the present we adopt a simplified model of eccentricity damping. Following Goodchild & Ogilvie (2006), we employ an effective bulk viscosity described by a Shakura-Sunyaev -parametrization, which leads to

(12)

where is a dimensionless parameter such that the effective bulk viscosity is given in terms of the pressure by . This diffusive term reduces the AMD of the system, and is therefore a convenient representation of eccentricity damping. This effective bulk viscosity represents whichever process (thermal or mechanical), apart from resonances, acts to damp the eccentricity, and should not necessarily be regarded as a hydrodynamic viscosity.

2.4 Resonances

When the planet is on an eccentric orbit, its gravitational potential can be expanded in a Fourier series of the form

(13)

where the coefficients can be found, e.g., in Goldreich & Tremaine (1978). Here we simply state that they are proportional to where is the planet’s eccentricity. For a planet on a circular orbit, only terms with appear. At the first order in eccentricity, we keep terms with . Here and in the remaining of Section 2.4 (and also in Section 3), we adopt the convention that the upper sign applies to inner resonances (i.e., a disc interior to the planet) and the lower sign applies to outer resonances (i.e., a disc exterior to the planet).

In a linear theory in which we work to first order in , two types of resonances need to be considered. Eccentric Lindblad resonances (ELRs) correspond to locations in the disc where the perturbing frequency in the rotating frame matches , the epicyclic frequency . We further assume a Keplerian disc, for which . Hence, inner (resp. outer) ELRs are of the form (resp. ) for , where . Eccentric corotation resonances (ECRs) correspond to . Hence they occur at commensurabilities of the form and for , with for an inner ECR and for an outer ECR, respectively. Irreversible growth and decay of eccentricity may occur because of resonant interactions between the planet and the disc.

2.4.1 Eccentric Lindblad resonances

As we show in Appendix B, based on Ogilvie (2007), a single ELR will contribute to the equation for the complex eccentricities of the disc and planet in the following way:

(14)
(15)
3 0.4807 0.8332 1.5397
4 0.6300 2.1863 3.2019
5 0.7114 3.8918 5.1463
6 0.7631 5.8890 7.3456
7 0.7991 8.1406 9.7756
8 0.8255 10.6210 12.4173
9 0.8457 13.3108 15.2555
10 0.8618 16.1950 18.2777
10 1.1604 21.3831 18.9703
9 1.1824 18.2033 15.9081
8 1.2114 15.1989 13.0272
7 1.2515 12.3809 10.3395
6 1.3104 9.7628 7.8592
5 1.4057 7.3617 5.6037
4 1.5874 5.2009 3.5940
3 2.0801 0.6072 1.8490
Table 1: Resonant radii and coefficients for eccentric Lindblad resonances

The resonant radii and the dimensionless coefficients and are displayed in Table 1. These equations show that a single ELR tends to cause an exponential growth of either the eccentricity of the disc or that of the planet, if the other is circular. In general it is a ‘relative’ eccentricity, as measured by the linear combination of , that tends to grow. The radial profile of the resonance is described through the dimensionless function and a resonant width , which represents the broadening of Lindblad resonances by collective effects such as pressure and self-gravity. Meyer-Vernet & Sicardy (1987) showed that the radial distribution of the torque exerted at a Lindblad resonance in a non-self-gravitating disc is given by an Airy function of a scaled radial coordinate. The peak of the Airy function is very close to , which is on the side of the resonance on which a density wave is emitted. We model the torque distribution using a Gaussian function that is off-centred by one resonant width. We find that a Gaussian evaluated at gives a reasonable representation of the Airy function. We recall that, following our convention, the (resp. ) sign is for an inner (resp. outer) Lindblad resonance. As we will discuss further, the off-centred torque distribution has the important effect of pushing more resonances inside the disc.

A simple estimate of the resonant width can found from the dispersion relation for density waves in a Keplerian disc (), where self-gravity is neglected:

(16)

Here is the radial wavenumber, which vanishes at the resonance. Lindblad resonances launch a wave in the disc, which propagates away from the planet. Differentiating equation (16) gives

(17)

Substituting by and assuming that , which means that the wavenumber has increased to at a distance from the resonance, we find

(18)

where is a typical value of the aspect ratio of the disc (see Section 4.3). Indeed, Meyer-Vernet & Sicardy (1987) already showed that the width of Lindblad resonances in a non-self-gravitating disc was proportional to . Lubow (2010) found that the non-zero width of the 3:1 ELR was important for the eccentric mode in superhump binaries. Table 1 shows that resonances with large values of are intrinsically stronger than those with low values of , but accumulate close to the planet, where the surface density is likely to be zero. The dependence of equation (18) on means that resonances with large values of have a decreasing width. In the remainder of this paper, we apply a torque cutoff at which ensures that resonances with higher are ineffective (Goldreich & Tremaine, 1980). Finally, we note that the coefficient is greatly reduced for the 1:3 resonance, because of the contribution of the indirect term in equation (97). Hence we expect this resonance to give only a small contribution to the eccentricity growth in the disc, especially if the planet is fixed on a circular orbit.

2.4.2 Eccentric corotation resonances

Similarly to what we derived for ELRs, a single ECR gives the following contribution:

(19)
(20)

where again, the (resp. ) sign is for an inner (resp. outer) corotation resonance.

2 0.6300 1.0853 0.3906
3 0.7631 2.2367 2.7434
4 0.8255 3.3933 3.9223
5 0.8618 4.5517 5.0930
6 0.8855 5.7109 6.2601
7 0.9023 6.8705 7.4252
8 0.9148 8.0304 8.5890
9 0.9245 9.1905 9.7522
9 1.0817 9.9412 10.5488
8 1.0931 8.7780 9.3887
7 1.1082 7.6141 8.2288
6 1.1292 6.4489 7.0692
5 1.1604 5.2817 5.9099
4 1.2114 4.1107 4.7515
3 1.3104 2.9309 3.5949
2 1.5874 1.7229 0.6200
Table 2: Resonant radii and coefficients for eccentric corotation resonances

Fluid elements in the corotation region are subject to libration around the nominal resonant radius. In the absence of dissipative effects such as viscosity, the corotation region is a closed system with a limited budget of angular momentum to transfer to the planet. Hence, the corotation torque will eventually saturate in a non-dissipative disc. When viscosity is included, the torque will also eventually vanish if the viscous diffusion timescale across the libration region is larger than the libration timescale. Thus, there exists a critical viscosity below which eccentricity growth is possible. In this picture, the resonant width is defined by the size of the libration zone. An estimate of this width was obtained by Goldreich & Tremaine (1981) using a ballistic treatment.

Following Masset & Ogilvie (2004), we adopt the following form:

(21)

where are coefficients which can be found in Tables 1 and 2 of Ogilvie & Lubow (2003), and which we set to unity for simplicity, as most of the uncertainty in our knowledge of the width lies in our choice for (where could be the eccentricity of the planet, that of the disc, or a linear combination of the two). The formula of Masset & Ogilvie (2004) gives the full width of the resonance, and thus should be divided by a factor of to give the correct Gaussian width. In practice, this means that the width of ECRs can be significantly smaller than the width of ELRs.

2.5 Short-range forces

2.5.1 General relativity correction

Relativistic corrections at the first order in post-Newtonian expansion cause an additional, prograde, contribution to the apsidal motion of the planet and the disc. Its effect has long been incorporated and appreciated in the framework of secular planetary dynamics (see, e.g., Adams & Laughlin, 2006). To leading order in eccentricity, it reads

(22)

for the eccentricity in the disc, and

(23)

for the eccentricity of the planet (here is the speed of light in vacuum).

2.5.2 Stellar oblateness

The stellar rotation causes the star to become oblate, generating a force on the orbit that causes it to progress in a prograde manner. This can be easily included in the framework of secular theory (see, e.g., Murray & Dermott, 1999) as

(24)

and

(25)

where is the stellar radius. In addition represents the leading order correction from a spherically symmetric gravitational potential. Values of are poorly known and depend on the stellar spin frequency. Values for the Sun are typically low (about ), but could be higher for fast-rotating T Tauri-type stars.

Although we give the expression for the stellar oblateness, we have not included it in any of the calculations in the remaining of the paper. Here we merely state that for rapidly rotating stars, its influence on the precession of a hot Jupiter could be as important as that of general relativity, but we chose not to include it because of our poor knowledge of . Both effects simply add a prograde contribution to the precession (see equation 39), and the absence of the term does not affect the overall conclusions of the paper.

3 Precession rate and growth rate

3.1 Angular momentum deficit

By analogy with celestial mechanics, it is useful to consider the angular momentum deficit (AMD) of the system. For small eccentricities, the total AMD of the disc is

(26)

In this Section and in Section 3.2, all integrals are carried from to , the radial extent of the disc. With appropriate boundary conditions, the AMD is conserved for adiabatic discs, in the absence of viscous damping or resonances (see Goodchild & Ogilvie, 2006, and details in Appendix C).

A planet with a complex eccentricity has an associated AMD of the form (for small )

(27)

When secular interactions between the disc and the planet are included (equations 2.2 and 2.2), AMD is exchanged between the disc and the planet. For an adiabatic, inviscid disc, the conserved quantity is now the total AMD of the system,

(28)

Including the self-gravity of the disc (equation 2.2) does not affect the AMD conservation. However, resonances will affect the AMD. When resonances are included (equations 2.4.12.4.1 and 2.4.22.4.2), the rate of change of AMD associated with each ELR and ECR is

(29)
(30)

We note that ELRs always act to increase the AMD, with an effect that depends on the linear combination of the complex eccentricities of the planet and the disc near the resonant location (and is approximately proportional to ). Individual ECRs can either increase or decrease the AMD depending on their location and the local vortensity gradient. ECR saturation, albeit a non-linear effect, should be considered here as it can operate at very small eccentricities. Linear equations can be obtained in the two bracketing cases in which ECRs are either completely unsaturated or completely saturated.

Viscosity acts to damp the AMD at a rate

(31)

In order to cause eccentricity growth, the net effect of mean-motion resonances on the AMD should exceed that of viscous damping.

3.2 Normal modes and integral relations

The simplest solutions of our system are normal modes of the form , where is a complex eigenfrequency. Such a mode corresponds to a fixed distribution of elliptical orbits of the disc and planet(s), which precesses at a rate given by the real part of and grows at a rate given by minus the imgainary part of . If the system supports one or more normal mode with a positive growth rate, we may expect the fastest-growing mode to dominate the evolution of the system until it grows into a nonlinear regime.

Using the integral relations we derived in Appendix C, and adding contributions from the planet, self-gravity and short-range forces, we find that the precession rate can be written

(32)

The different terms are defined as follows:

(33)

with for the 2D adiabatic model, for the 3D adiabatic model, and for the locally isothermal model;

(34)

with for the 3D adiabatic model, and otherwise;

(35)

a non-adiabatic contribution for the locally isothermal model only, and for the adiabatic model;

(36)

for 3D adiabatic or locally isothermal models respectively;

(37)

for self-gravity;

(38)

for secular disc–planet interaction, and

(39)

for short-range forces, where

(40)

It can be shown that , and therefore and are positive definite. Gravitational coupling leads to prograde precession, as do the short-range forces. In the 2D case, which was already studied by Goodchild & Ogilvie (2006), we see that the two terms associated with pressure cause a retrograde precession of the mode, provided that the pressure gradient is negative. The 3D term gives rise to a prograde precession rate, and so do the terms coming from the planet and short-range forces. The total precession rate of the mode depends on the net contribution of all these effects.

Similarly, the growth rate reads

(41)

where

(42)
(43)
(44)

and

(45)

for eccentric Lindblad resonances, eccentric corotation resonances, viscosity and non-adiabatic effects (in the locally isothermal model only), respectively. It clearly shows that the contribution from each ELR always acts to increase the growth rate, while viscosity reduces it. The effect of each ECR depends on the local vortensity gradient, but their net effect usually contributes to the damping of eccentricity. It is important to notice that, in expression (41) for the growth rate, the AMD of the mode appears in the denominator. In anticipation of the solutions in terms of normal modes that will be presented in the remainder of the paper, we can expect that, if an eccentric mode is confined in the inner parts of the disc, and decays towards zero in the outer parts, it should present a higher growth rate than a less confined mode, if the source terms due to resonances are comparable.

4 Solution in terms of normal modes

4.1 Numerical methods

In order to solve the eccentricity equations of Section 2, we seek normal modes of the form for the complex eccentricity. We can then solve the resulting eigenvalue problem for the (usually complex) frequencies . This is achieved by dividing the disc into rings, spaced equally in the variable , and associating an eccentricity to each ring. The discrete model that results from this scheme is described in Appendix C. We use LAPACK to solve for the eigenvalues and eigenvectors of the generalized eigenvalue problem. Some spurious solutions can be generated that show oscillations on the grid scale, and which should be discarded. Each mode represents a radial eccentricity distribution, and its precession and growth rates are given by the real and (negative) imaginary part of the eigenfrequency, respectively. As a measure of the mode order, we take the number of zeros of the real part of the eigenvector. While this is strictly valid only for a “classical” Sturm-Liouville problem with real functions of , we also use it as an indicator of the mode order for our more general problem. When necessary, integrals are solved numerically using a Romberg method (Press et al., 1992).

The disc extends from an inner radius to an outer radius . For a disc with free boundaries at both edges, the boundary conditions are

(46)

at and , where is a function of defined in Appendix C. These conditions are automatically satisfied by our prescription for the surface density and pressure, as both fall to zero at each edge (see below).

The grid resolution has to be such as the low-order modes that are of interest here do not show any variation with resolution. We find that dividing the disc into rings gives consistent results. A good internal check on the accuracy can be obtained by comparing the real and imaginary parts of the eigenfrequency with the ones obtained using the integral formulae in Section 3.2.

4.2 Units

We choose to work in physical units (astronomical units, solar masses and years). For most parts, the disc–planet interactions could be described in arbitrary units, but the inclusion of short-range forces such as the general relativistic correction, which depend on the distance to the star, make it more convenient to choose physical units. Unless mentioned otherwise, the precession and growth rates are given in units of the planet’s orbital frequency. We keep this convention even when the planet is not included, in order to help the comparison between the various models. In our fiducial example, the planet’s orbital frequency is calculated at , a distance typical of hot Jupiters. Finally we recall that since we are solving a normal mode problem, the solution can be scaled by any factor and should not be interpreted as the true eccentricity, but as being arbitrarily normalized.

4.3 Disc model

In order to solve the eccentricity equation, one needs a description of the disc structure. Following Bell & Lin (1994), we note that in the inner parts of the disc (from to ), the temperature ranges from 100 to 1000 K. In this region the opacity is mostly proportional to . This allows us to construct a simple disc model in which the surface density is given by

(47)

and the aspect ratio follows

(48)

where and are characteristic values of surface density and radius, and is the limiting value of the aspect ratio at large radius. is defined such that the vertically integrated pressure is

(49)

as follows from vertical hydrostatic equilibrium. With our choice for and , we have

(50)

In the locally isothermal case, we relate the isothermal sound speed to through

(51)

Apart from the factor, this model corresponds to a steady accretion disc with constant alpha viscosity and opacity . In this model, the surface density, aspect ratio and pressure naturally go to zero at the inner edge of the disc. At the outer edge, we also wish to have zero surface density and pressure, which is achieved by using a taper function, in which the drop to zero occurs on a width defined through . Here is taken to be . The value has to be chosen according to the resolution so that quantities smoothly drop off to zero over a few grid points. The exact value has no important effect on the solution. We plot the surface density, aspect ratio and pressure radial profiles in Fig. 1.

In this model, one expects the surface density to be around at (Bell et al., 1997). This implies that is of the order a few times in our units.

Figure 1: Normalized surface density, aspect ratio and pressure as functions of radius (in logarithmic scale).

5 A non-self-gravitating and inviscid disc

In this section we solve the normal mode problem described in Section 4, where the eccentricity is assumed to be of the form . We start with the case of a disc in which viscosity, self-gravity and resonant interactions with the planet are excluded. However, we keep the secular part of the planet–disc interaction. The eccentricity evolution is governed by one of the models described by equations (1) to (2.1), and equation (2.2) with .

We first illustrate the solution for a disc that extends from au to au, with a mass ratio , an aspect ratio , and a viscosity . In the adiabatic case, we take the adiabatic index to be 111Bitsch et al. (2013) pointed out that the adiabatic index can vary from 5/3 to 7/5 with increasing disc temperature. However the exact value of the adiabatic index does not significantly change our results.. The planet is on a fixed circular orbit at au, with a mass ratio of , and we neglect the self-gravity of the disc.

The total mass of the disc is a somewhat arbitrary parameter, since for a given surface density profile it depends on the inner and outer disc radius. Here our mass is chosen such as , in the range of values expected for such discs. It ensures that self-gravity can be neglected. We note here that if the disc–planet interaction is localized in the inner part of the disc, the total disc mass is irrelevant, as most of it is located in the outer part. A more relevant quantity is the local disc mass near the inner edge which we define as . With our disc parameters, the local disc mass is expected to be much lower than the mass of a Jovian planet.

5.1 Bound states in a potential well

In the absence of self-gravity and viscosity, the eccentricity propagates through the disc solely by means of pressure, as a one-armed density wave. The evolutionary equations of Section 2.1 are dispersive wave equations related to the Schrödinger equation of quantum mechanics. A basic understanding of the physics at play can be achieved by rewriting the equation for an eccentric normal mode as a time-independent Schrödinger equation:

(52)

where is a new independent variable related to , is the wavefunction, is an effective potential and is the effective energy eigenvalue. In the following we assume that the sound speed can be written , where is taken to be the characteristic aspect ratio . In Appendix D we present a derivation in which can take a more general form such as the one we use in Equation (48). The simple model that we use here is sufficient to capture the relevant physical mechanism.

In the case of a locally isothermal disc, the transformation required to obtain a Schrödinger equation is

(53)

where . In the 3D case, an additional term is added to .

In the adiabatic case we have

where . In 2D, we have , and , while in 3D we have , , .

In both cases, the shape of the effective potential can help in understanding the trapping of the modes by analogy with bound states of the Schrödinger equation. Typical potentials are plotted in Fig. 2. The results are rather similar in the locally isothermal and adiabatic cases. The profiles show that a deep potential well exists in 3D, which could support a bound state in the form of a confined eccentric mode that decays at . Since as , a bound state has , and is therefore a prograde mode . In 2D, the potential well is less deep, and bound states might not exist. This distinction, which arises from the additional prograde precession due to pressure in a 3D disc, was identified as important by Ogilvie (2008). The potential arising from the planet also contributes to deepening the potential well. It will turn out to be of importance for helping the confinement of 2D growing modes.

Figure 2: Effective potential for the analogous Schrödinger equation, for the 2D and 3D adiabatic (left panel) and locally isothermal (right panel) models. The solid line is the potential without the planet, while the dotted line includes the planet’s contribution. The planet has , , and the disc has .

5.2 Mode analysis for a non-self-gravitating disc

We start by solving the eigenvalue problem for a simple case of a non-self-gravitating, inviscid disc without a planet, short-range forces or self-gravity. We use the parameters given at the beginning of Section 5.

In Fig. 3, we represent the mode shape of the four lowest-order modes. In this simple model, major differences between 2D and 3D adiabatic and isothermal discs appear. In 3D, there exists a mode of order 0 which is confined in the inner parts of the disc, in agreement with what we speculated in Section 5.1. This mode decays on the characteristic length of 1 au, and is insensitive to the outer boundary condition as long as . Higher-order modes do not always decay to zero at large radii, and might therefore be sensitive to the outer radius. However they contain a significantly larger amount of AMD, and are therefore expected to grow less rapidly, as we will confirm in Section 6.1. In Table 3 we show the precession rates of these modes. The disc supports three prograde modes, all the others being retrograde. Although we do not include a planet in this calculation, we still give the precession and growth rates in units of the planet’s orbital frequency at au, since it will help comparing with the next section where we include the planet.

In 2D, there no longer exists a mode which is trapped in the inner parts of the disc (see right panels of Fig. 11). This is especially true in the adiabatic case where all modes show a large eccentricity distribution at large radii. As can be seen in Table 3, 2D discs only support a set of retrograde modes.

Figure 3: The four lowest-order modes of our fiducial example for an inviscid, non-self-gravitating disc without a planet, for the four different models: 3D adiabatic, 2D adiabatic, 3D isothermal and 2D isothermal. The number in the top right corner indicates the order of the corresponding mode. Their precession rates are given in Table 3.
Mode order 3D adiabatic 2D adiabatic 3D isothermal 2D isothermal
0
1
2
3
Table 3: Precession rate of the four lowest-order modes of our fiducial example, for an inviscid, non-self-gravitating disc without planets, in the 2D and 3D abiabatic and isothermal cases. The shape of the corresponding modes is shown in Fig. 3.

5.3 Secular perturbations from a planet

The secular perturbation from the planet gives rise to an additional source of prograde precession in the disc. Previously, we have noted that the planet helps to deepen the effective potential well, and could therefore assist the trapping of a mode in the inner part of the disc (see Section 5.1 and Fig. 2).

Our choice of semi-major axis is such that . The planet is expected to migrate into the cavity via tidal interaction with the disc, until its 1:2 outer Lindblad resonance (OLR) lies just inside the inner edge of the disc, which here would be . The planet should then remain at a constant semi-major axis, a convenient setup in our secular framework. With our choice of surface density profile, motivated by the truncation of the disc by a magnetospheric cavity, the surface density falls to zero at and exerts no torque on the planet within this radius. Instead we expect the planet to mostly interact with the peak of the surface density just outside the inner edge. Therefore the planet should stop slightly closer to the disc than the 1:2 orbital commensurability with the inner edge, and we choose the fiducial value of . We point out that in the study of planets in cavities by Rice et al. (2008), Jupiter-mass planets stalled at orbital distances over the few hundred orbits over which the simulations were performed. We recall that in our secular treatment of disc–planet interactions, it was assumed that the semi-major axis of the planet was constant in time. We discuss the influence of in Section 7.3.

In Fig. 4 we show the mode shape, and the associated precession rates are displayed in Table 4. In 3D, the effect of adding a planet on the mode shape is negligible, as the mode is already trapped in the deep effective potential due to 3D pressure effects and discussed in Section 5.1. The planet enhances the prograde precession of these modes, but does not strongly affect the retrograde modes, which are pure disc modes. In 2D, the planet is observed to reduce the amplitude of the mode in the outer parts of the disc compared to that of the inner parts. Here a significant difference arises between the 2D isothermal and adiabatic cases. In the isothermal case, the additional contribution of the planet to the potential well is enough to confine the lowest-order mode in the inner parts of the disc, while in the adiabatic case this mode ‘leaks’ into the outer disc. When resonant and viscous effects are added in the next section, this will have important consequences for the growth of the mode. In Table 4 we can see that the lowest-order 2D isothermal mode has become prograde owing to the inclusion of the planet, while the 2D adiabatic modes are only weakly affected by the planet (see Table 3 for comparison). In the example that we show here, the perturbation from the planet is not strong enough to cause any mode to precess progradely in the 2D adiabatic calculation. However more massive planets closer to the disc’s inner edge could assist the confinement of a 2D progade mode.

Figure 4: The four lowest-order modes of our fiducial example for an inviscid non-self-gravitating disc secularly perturbed by a planet (but with no resonances or back-reaction on the planet), for the four different models: 3D adiabatic, 2D adiabatic, 3D isothermal and 2D isothermal. The number in the top right corner indicates the order of the corresponding mode. Their precession rates are given in Table 4.
Mode order 3D adiabatic 2D adiabatic 3D isothermal 2D isothermal
0
1
2
3
Table 4: Precession rate of the four lowest-order modes of our fiducial example, for an inviscid, non-self-gravitating disc, with the secular perturbation from a planet (), in the 2D and 3D abiabatic and isothermal cases. The shape of the corresponding modes is shown in Fig. 4.

6 A more realistic example

In this Section, we progressively add more physical effects to our disc model, and show how they affect the shape, precession rate and growth rate of the modes. The physical parameters for the disc and the planet are the same as in Section 5.

6.1 Viscosity and resonant interactions

We first add effects that will contribute to the growth rate: viscosity, and eccentric Lindblad and corotation resonances. The planet is kept on a fixed circular orbit. The shape of the modes is shown in Fig. 5, and the corresponding precession rates and growth rates are displayed in Table 5. The shape of the 3D modes is slightly different from Fig. 4 because the eccentricity now has an imaginary component that shifts the location of local minima away from zero. However, comparing Table 5 with Table 4 shows that the precession rate of these modes remain similar. In 3D, the lowest-order mode, which is confined in the inner parts of the disc, contains a very small amount of AMD. This mode therefore has the highest growth rate. The growth rates are fairly similar in the adiabatic and locally isothermal 3D cases, and the fastest growing mode can grow in orbital periods of the planet, which is short compared to the disc’s lifetime if the planet has an orbital period of 3 to 4 days.

In 2D, the lowest-order modes in the adiabatic case do not present significant differences from what was found in Section 5.3. The precession rates are similar and, because no mode is confined in the inner disc, all modes have relatively low growth rates, much lower than in 3D. However there exists a series of higher-order modes that decay slowly at large radii, but do not show the large-amplitude oscillations of the lowest-order mode. These modes are of order 6 to 10, and all show growth rates of about . We show in Fig. 5 and Table 5 the fastest growing of these modes, which is of order 7.

The 2D isothermal case is more subtle. As can be seen in Table 5, the only prograde mode that was found in Table 4 has now disappeared. The fastest growing mode is of order 2 (see Fig. 5), and seems at odds with other growing modes because its growth rate is faster than its precession rate. These peculiar modes only appear because of the resonant effects, which can impose a very fast growth rate to a mode, which then decays rapidly towards the outer boundary. This decay occurs before the mode has time to propagate to the outer boundary and reflect to form a standing wave. Because this mode is confined in the inner disc it grows faster than the fastest growing mode in the 2D adiabatic case. As we will see in Section 7.1 these modes can also appear in 3D. The analogy with a Schrödinger equation is harder to establish in this case since the potential now contains an imaginary part as well. In particular it is possible for a mode to be trapped in a complex effective potential with a retrograde precession rate and a positive growth rate. The three other modes in the 2D adiabatic regime have an equivalent in Table 4, although the mode order can be different.

Figure 5: The four lowest-order modes of our fiducial example for a viscous, non-self-gravitating disc, including secular and resonant interactions with a planet, for the four different models: 3D adiabatic, 2D adiabatic, 3D isothermal and 2D isothermal. The number in the top right corner indicates the order of the corresponding mode. Their precession rates are given in Table 5.
Model Mode order precession rate growth rate
3D adiabatic 0
1
2
3
2D adiabatic 1
2
2
4
7
3D isothermal 0
1
2
3
2D isothermal 0
1
2
3
Table 5: Precession rates and growth rates for the lowest-order modes of our fiducial example, when viscosity and secular and resonant interactions with the planet are included, in the 2D and 3D abiabatic and isothermal cases. The shape of the corresponding modes is shown in Fig. 5.

Using Equations (3.2)–(45), we can study the contribution of different physical effects to the growth and decay of eccentricity. In Fig. 6 we show, in the case of a 3D adiabatic disc, the normalized contribution from the first few ELRs and ECRs. The contribution of ECRs is actually negative (causing a damping of eccentricity), and we plot its absolute value. The contribution from viscosity is negligible in this case. It clearly appears that the only ECR contributing to the damping is the 1:2. Regarding ELRs, the 1:3 is rather ineffective at exciting the eccentricity. The main contributions come from the 2:4 and 3:5 ELRs, with a smaller contribution from the 4:6. Although the 2:4 ELR shares the same nominal radius as the 1:2 ECR, we recall that the peak of the torque distribution of the ELR is actually slightly shifted away from the planet. In addition, the width of ELRs is significantly larger than the width of ECRs, allowing for a larger part of the disc to contribute to the growth. The consequence of the shift of the centre of the resonance is even more visible when comparing the 4:6 ELR with the 2:3 ECR. The 2:3 ECR lies in the cavity and is inoperative in the disc, while the 4:6 ELR is partly in the disc and contributes to the eccentricity growth.

Figure 6: Contribution to the growth rate of the fastest growing mode in the 3D adiabatic model for various resonances. Eccentric Lindblad resonances (ELRs) are represented as full circles, while eccentric corotation resonances (ECRs) are empty squares. All contributions are normalized so that their sum is one. ECRs give a negative contribution to the growth rate, and we show here their absolute value. The 2:4, 3:5 and 4:6 ELRs contribute to most of the growth rate, much more than the 1:3 ELR. The only ECR to contribute significantly is the 1:2.

6.2 Back-reaction on the planet

We now allow the disc to act back on the eccentricity of the planet (i.e. we do not assume any more in all the equations where appears). Fig. 7 and Table 6 show a few modes of interest for the disc, and the mode that is dominated by the planet (defined as the mode for which the AMD is predominantly in the planet, and labelled “P”). The planet-dominated mode undergoes a small prograde precession, and has a slow growth rate.

We find that the growth rate of the planet-dominated mode is significantly smaller than that of the disc mode, when they both exist. Interestingly, the 2D planet mode grows on a timescale equal to, or larger than that in 3D. This effect will appear consistently for a large range of parameters described in Section 7. In 2D, the planet mode does not propagate very far into the disc, and therefore its AMD is low. On the other hand, the 3D planet mode is more strongly coupled to the disc. It shows some oscillations and propagates further into the disc. Its AMD is therefore larger, reducing its growth rate.

In this calculation, we also include the relativistic precession term. It appears that its effect on the disc is negligible compared to the effect of pressure, but it provides the dominant source of precession on the planet, competing with the precession induced by the gravitational potential of the disc. For instance, in the 3D adiabatic case, the precession rate of the planet-dominated mode is without the GR term, while it is when GR is included.

Figure 7: The four lowest-order modes of our fiducial example for a viscous, non-self-gravitating disc, including secular and resonant interactions with a planet, and back-reaction on the planet, for the four different models: 3D adiabatic, 2D adiabatic, 3D isothermal and 2D isothermal. The relativistic contribution to precession is also included. We show here the distribution of eccentricity in the disc, truncated at the disc’s inner edge . All modes but one have a negligible amount of eccentricity (or AMD) in the planet, and we do not show it here. For the mode dominated by the planet, the eccentricity in the planet is represented by a black dot at the location of the planet (), which links to the corresponding mode in the disc. The number in the top right corner indicates the order of the corresponding mode. Their precession rates are given in Table 6.
Model Mode order precession rate growth rate
3D adiabatic 0
1
3
P
2D adiabatic 0
1
7
P
3D isothermal 0
1
2
P