text:
# 1 Motivation

## 1 Motivation

At the heart of Einstein’s theory of gravity is local diffeomorphism invariance: coordinate systems are unimportant, only diffeomorphism invariants matter. However, a poor choice of local coordinates can sometimes obscure the nature of the global aspects of spacetime, such as horizons or causal boundaries. Indeed, the true nature of the “coordinate singularity” at the Schwarzschild radius eluded Einstein himself, and was only fully illuminated with the discovery of coordinate systems that were regular at the horizon. Coordinate systems that cover larger patches of spacetime are especially useful to have in dealing with physical phenomena that are in some sense nonlocalized.

In this note, we present a simple new coordinate system for de Sitter space, covering the causal future/past of an observer. The new coordinates, which we shall call Painlevé-de Sitter coordinates, are a cross between static coordinates and planar coordinates, and inherit the strengths of each of these. For example, like static coordinates but unlike planar coordinates, the new coordinates have a direction of time that is a Killing vector, making them well-adapted to thermodynamics. On the other hand, like planar coordinates, but unlike static coordinates, Painlevé-de Sitter coordinates continue smoothly through the horizon, and constant time slices are just flat Euclidean space. Painlevé-de Sitter coordinates differ also from Eddington-Finkelstein type coordinates in that the coordinates are all either timelike or spacelike, rather than null.

The combination of a Killing time direction and regularity at the horizon is particularly powerful as it allows one to study across-horizon physics as seen by an observer. A natural application is de Sitter radiation. Heuristically, one envisions de Sitter radiance as arising in much the same way that Hawking radiation does. That is to say, a particle-pair forms just inside the horizon, one member of the pair tunnels across the horizon, and the virtual pair becomes real. To show that this is actually what happens, one would like to compute the amplitude for traversing the horizon; because of their regularity, Painlevé-de Sitter coordinates make the calculation feasible. By directly evaluating the imaginary part of the action, we obtain the emission amplitude associated with a tunneling particle. In the s-wave limit, it is in fact possible to extend the computation to include the effects of self-gravitation. As a result, the de Sitter spectrum turns out to be only approximately thermal. In particular it has a cutoff at high energies. The back-reaction is such that, unlike Schwarzschild black holes, de Sitter space lowers its temperature as it radiates.

## 2 Painlevé-de Sitter Coordinates

Many different coordinate systems are known for de Sitter space (see, e.g., ). One commonly used metric is the static metric, analogous to the familiar Schwarzschild metric for an uncharged black hole. This metric covers a static patch, that part of de Sitter space that an observer at the origin can interact with. The time coordinate, $`t_s`$, corresponds to a timelike Killing vector, which makes it suitable for thermodynamics; thermal equilibrium requires among other things that the spatial metric be at equilibrium. However, the static metric has the limitation that it is only valid upto the horizon; it therefore covers only a very small region of the full space. Another drawback is that in a static background one cannot expect to get radiation, a phenomenon that is manifestly time-reversal asymmetric. Indeed, in early calculations of Hawking radiation from black holes, time-reversal symmetry had to be broken by hand through the introduction of a collapsing surface.

We shall now illustrate the method for obtaining Painlevé-type coordinates, valid across the horizon. Consider then a general static metric of the form

$$ds^2=(1g(r))dt_s^2+\frac{dr^2}{1g(r)}+r^2d\mathrm{\Omega }_{D2}^2.$$
(1)
Here $`g(r)=1`$ corresponds to a horizon; we shall assume for simplicity that there is only one horizon. For de Sitter space, $`g(r)=r^2/l^2`$ and $`r=l`$ marks the horizon. $`r=0`$ could be the worldline of an observer at the origin.

To obtain the new line element, define a new time coordinate, $`t`$, by $`t_s=t+f(r)`$. The function $`f`$ is required to depend only on $`r`$ and not $`t`$, so that the metric remains stationary, i.e. time-translation invariant. Stationarity of the metric automatically implements the desirable property that the time direction be a Killing vector. What other conditions should we impose on $`f`$? Our key requirement is that the metric be regular at the horizon. We can implement this as follows. We know that a radially free-falling observer who falls through the horizon does not detect anything abnormal there; we can therefore choose as our time coordinate the proper time of such an observer. As a corollary, we demand that constant-time slices be flat. We then obtain the condition

$$\frac{1}{1g(r)}(1g(r))(f^{}(r))^2=1.$$
(2)
For de Sitter space, $`g(r)=r^2/l^2`$, and the solution to Eq. (2) yields

$$t_s=t\pm \frac{l}{2}\mathrm{ln}(1r^2/l^2).$$
(3)
As usual, the fact that the transformation is singular at $`r=l`$ merely indicates that the original coordinate $`t_s`$ was ill-defined there. Choosing the minus sign for now, gives us the new coordinates, which we shall call Painlevé-de Sitter coordinates. The desired line element is

$$\overline{)ds^2=(1\frac{r^2}{l^2})dt^22\frac{r}{l}dtdr+dr^2+r^2d\mathrm{\Omega }^2.}$$
(4)
The Painlevé-de Sitter metric has a number of attractive features. First, none of the components of either the metric or the inverse metric diverge at the horizon. Second, by construction constant-time slices are just flat Euclidean space. Third, the generator of $`t`$ is a Killing vector. “Time” becomes spacelike across the horizon, but is nevertheless Killing; this fact can be exploited to compute global charges such as mass  in a natural way. Finally, an observer precisely at the origin does not make any distinction between these coordinates and static coordinates; indeed, the function $`f`$ that distinguishes the two time coordinates vanishes there.

We call this the Painlevé-de Sitter metric because the analogous line element for four-dimensional Schwarzschild black holes is

$$ds^2=\left(1\frac{2M}{r}\right)dt^22\sqrt{\frac{2M}{r}}dtdr+dr^2+r^2d\mathrm{\Omega }_2^2.$$
(5)
This superb though relatively unknown line element was discovered by Painlevé  many years ago; it was rediscovered in a modern context by Kraus and Wilczek . Similar coordinate systems have been found for black holes in anti-de Sitter space .

These coordinates cover the observer’s causal future, or half of the full space. The metric is necessarily nonstatic (that is, not time-reversal invariant) as the causal patch is itself not time-reversal invariant. To obtain the other half of the space, corresponding to the causal past of an observer at the antipode, one chooses the opposite (plus) sign for the off-diagonal component and reverses the sense of time, so that $`t`$ increases to the past. Choosing the other sign in Eq. (3) without reversing the direction of time gives a metric that covers the causal past of the original observer. The relation between Painlevé-de Sitter coordinates and planar coordinates is made clear by the transformation

$$r=\rho e^{t/l}ds^2=dt^2+e^{2t/l}\left(d\rho ^2+\rho ^2d\mathrm{\Omega }_{D2}^2\right).$$
(6)
We see that the new coordinates have the same radial coordinate as static coordinates and the same time coordinate as planar coordinates.

We conclude this section by writing down the radial geodesics in these coordinates. The null radial geodesics obey

$$\frac{dr}{dt}=r/l\pm 1,$$
(7)
where the plus (minus) sign corresponds to rays that go away from (towards) the observer. When the particle is beyond the horizon ($`r>l`$), both ingoing and outgoing trajectories correspond to increasing $`r`$, and the particle cannot (classically) cross the horizon. The general solution is

$$r(t)=l\left(e^{t/l}1\right).$$
(8)
Radially-directed null rays leaving the observer at $`t=0`$ reach the horizon at $`t=l\mathrm{ln}2`$, and reach future null infinity at $`t=\mathrm{}`$. Turning now to massive particles, the radial geodesic equation implies

$$U^r=\frac{dr}{d\tau }=\left(1r^2/l^2\right)m^2p_t,$$
(9)
where $`\tau `$ is the proper time. Here (by stationarity) $`p_t`$ is a constant of the motion, being equal to the energy measured by an observer at $`r=0`$. Note that if we set $`p_t(=mU_t)m`$, then using $`U^2=1`$ we find that

$$U^t=1\tau =t+c,$$
(10)
so the Painlevé-de Sitter time coordinate is nothing more than the proper time along a radial geodesic worldline, such as that of a free-falling observer. Indeed, that is precisely how we arrived at these coordinates in the first place.

## 3 Tunneling Across the de Sitter Horizon

The great utility of having a coordinate system that is well-behaved at the horizon is that one can study across-horizon physics. In this section, we will determine the temperature of de Sitter space by directly computing the rate at which particles tunnel across the horizon. Our computation will parallel the analogous calculation for black holes, in which Hawking radiation is expressed as a tunneling phenomenon .

Now, because of the infinite blueshift near the horizon, the characteristic wavelength of any wavepacket is always arbitrarily small there, so that the geometrical optics limit becomes an especially reliable approximation. This is of course a big plus: the geometrical optics limit allows us to obtain rigorous results directly in the language of particles, rather than having to use the unwieldy and physically less transparent Bogolubov method that is more traditionally used. Moreover, as we shall see, the inclusion of back-reaction effects is also perhaps easier now.

In any event, since we are in the semi-classical limit, we can apply the WKB formula. This relates the tunneling amplitude to the imaginary part of the particle action at stationary phase. (The phase is $`iL𝑑t/\mathrm{}=i(p\dot{x}H)𝑑t/\mathrm{}`$; since energy is real, exponential damping comes from the imaginary part of emission rate $`p𝑑x`$, which in the nonrelativistic limit becomes the usual $`𝑑x\sqrt{2m(V(x)E)}`$. The imaginary part of the action, $`I`$, is thus given by the imaginary part of the momentum integral.) The emission rate, $`\mathrm{\Gamma }`$, is the square of the tunneling amplitude:

$$\mathrm{\Gamma }\mathrm{exp}(2\mathrm{Im}I/\mathrm{})\mathrm{exp}(\beta E).$$
(11)
On the right-hand side, we have equated the emission probability to the Boltzmann factor for a particle of energy $`E`$. To the extent that the exponent depends linearly on the energy, the thermal approximation is justified; we can then identify the inverse temperature as the coefficient $`\beta `$.

We will consider here the s-wave emission of massless particles. Higher partial wave emission is in any case suppressed by $`\mathrm{}`$. In the s-wave, particles are really massless shells. If we imagine a shell to consist of constituent massless particles each of which travels on a radial geodesic, then we see that the motion of the shell itself must follow the radial null geodesic for a particle. That is, it obeys Eq. (7), with the minus sign. We will use these radial geodesics to compute the imaginary part of the action, as follows.

Since the calculation involves a few tricks , we outline it here before putting in the details. First observe that we can formally write the action as

$$\mathrm{Im}I=\mathrm{Im}_{r_i}^{r_f}p_r𝑑r=\mathrm{Im}_{r_i}^{r_f}_0^{p_r}𝑑p_r^{}𝑑r,$$
(12)
where $`p_r`$ is the radial momentum. We expect $`r_i`$ to correspond roughly to the site of pair-creation, which should be slightly outside the horizon. (Note that the second member of the pair contributes nothing to the tunneling rate, since it is always classically allowed and therefore has real action.) We expect $`r_f`$ to be a classical turning point, at which the semi-classical trajectory (i.e. instanton) can join onto a classical-allowed motion. This must be slightly within the horizon, else the particle would not be able to propagate classically from there to the observer. However, the precise limits on the radial integral are unimportant, so long as the range of integration includes the horizon.

We now eliminate the momentum in favor of energy by using Hamilton’s equation

$$\frac{dH}{dp}|_r=\frac{H}{p}=\frac{dr}{dt},$$
(13)
where the Hamiltonian, $`H`$, is the generator of Painlevé time. Hence within the integral over $`r`$, one can trade $`dp`$ for $`dH`$. Without being very careful about signs, the integral over $`H`$ now just gives the particle energy $`E`$. However, substituting Eq. (7) for a particle going radially towards $`r=0`$, the radial integral has a simple pole at the horizon:

$$\mathrm{Im}I=\mathrm{Im}_{r_i}^{r_f}_0^E\frac{dH}{\frac{dr}{dt}}𝑑r=\mathrm{Im}E_{r_i}^{r_f}\frac{ldr}{rl}.$$
(14)
The pole lies along the line of integration, and therefore yields $`\pi i`$ (rather than $`2\pi i`$) times the residue. Again, we postpone consideration of the sign associated with the direction of the contour. We get

$$\mathrm{Im}I=\pi lE.$$
(15)
Consulting Eq. (11), we find that this corresponds to a temperature

$$T_{dS}=\frac{\mathrm{}}{2\pi l},$$
(16)
which is precisely the temperature of de Sitter space. To summarize: Painlevé-de Sitter coordinates have allowed us to compute the radiation rate directly from the particle action, with the action incurring an imaginary part from a pole at the horizon.

Let us now do the calculation more carefully, keeping track of the signs, and including the effects of back-reaction. Perhaps it should be stressed that the reason we are interested in back-reaction is not merely to compute higher order in $`E`$ effects, but because self-gravitation is central to the entire process of across-horizon tunneling. Without self-gravitation the back-of-the-envelope calculation above is puzzling: if this is tunneling, where is the barrier? Put another way, if particles created just inside the horizon have only to tunnel just across – an infinitesimal separation – what characterizes the scale of the tunneling? Recall that in the Schwinger process of electron-positron pair production in an electric field, there is a nonzero separation scale, $`rmc^2/qE`$, between the classically allowed configurations. In the following, we will see that, as with Hawking radiation , self-gravitation resolves these issues. Back-reaction results in a shift of the horizon radius; the finite separation between the initial and final radius is the classically-forbidden region, the barrier.

How does one incorporate back-reaction? In a general situation, this is a notoriously difficult problem, calling for a theory of quantum gravity. Indeed, generically one has to worry about how to consistently quantize the gravitational waves produced by a matter source. However, for the special case of spherical gravity it is possible to integrate out gravity, at least semi-classically. This is because for spherical gravity, Birkhoff’s theorem (more precisely, its generalization to a nonzero cosmological constant) states that the only effect on the geometry that the presence of a spherical shell has, is to provide a junction condition for matching the total mass inside and outside the shell. (In three dimensions, where there are no gravitational waves, it may be possible to compute the emission rates for the higher partial waves as well.)

Since the geometry is different on the two sides of the shell, one can now ask which geometry determines the motion of the shell. (Thus, self-gravitation automatically breaks the principle of equivalence.) The geometry inside the shell is empty de Sitter space, while the geometry outside is that of Schwarzschild-de Sitter space with energy $`E`$. (Technically, empty de Sitter space has a mass too ; what follows is unaffected by this shift.) It is the outside $`E`$-dependent metric that determines the motion of the self-gravitating shell. Consider, for simplicity, $`dS_3`$; generalization to higher dimensions is straightforward. The effective geometry whose radial geodesic determines the motion of the shell has the line element

$$ds_{\mathrm{effective}}^2=(18GEr^2/l^2)dt_s^2+\frac{dr^2}{18GEr^2/l^2}+r^2d\varphi ^2.$$
(17)
Here we have inserted Newton’s constant, and $`E>0`$ is physically the energy of the shell as measured by an observer at $`r=0`$. The corresponding Painlevé metric is now

$$ds_{\mathrm{effective}}^2=(18GEr^2/l^2)dt^22\sqrt{r^2/l^2+8GE}dtdr+dr^2+r^2d\varphi ^2.$$
(18)
The imaginary part of the action is then

$$\mathrm{Im}I=\mathrm{Im}_0^H_{r_i}^{r_f}\frac{drdH^{}}{\sqrt{r^2/l^2+8GE^{}}1},$$
(19)
where we have inserted the radial geodesic derived from the effective metric. Here $`r_i=l`$ is the original radius of the horizon just before pair-creation, while $`r_f`$ is the new radius of the horizon, and is equal to $`l\sqrt{18GE}`$. What matters is that $`r_f<r_i`$. The Feynman prescription for evaluating the sign of the contour is to displace the energy from $`E^{}`$ to $`E^{}i\epsilon `$. Substituting $`u=r^2`$,

$$\mathrm{Im}I=\mathrm{Im}_0^H_{u_i}^{u_f}\frac{du}{2\sqrt{u}}\frac{l\sqrt{u+8Gl^2(E^{}i\epsilon )}+l^2}{u(l^28Gl^2(E^{}i\epsilon ))}𝑑H^{},$$
(20)
we see that the pole lies in the upper-half $`u`$-plane. Doing the $`u`$ integral first we find

$$\mathrm{Im}I=\pi l_0^H\frac{dH^{}}{\sqrt{18GE^{}}}.$$
(21)
Now the total energy of de Sitter space decreases when positive-energy matter is added to it , because of the negative gravitational binding energy. Therefore the Hamiltonian $`H`$ satisfies $`dH=dE`$, giving

$$\mathrm{Im}I=\frac{\pi l}{4G}\left(\sqrt{18GE}1\right),$$
(22)
and the tunneling rate is therefore

$$\overline{)\mathrm{\Gamma }\mathrm{exp}(+\frac{\pi l}{2G\mathrm{}}(\sqrt{18GE}1)).}$$
(23)
When the particle’s energy is small, $`8GE1`$, the square root can be approximated. To linear order in $`GE`$, we recover our previous back-of-the-envelope result, Eq. (16). As a check, we note that the sign has also come out correctly. To this order then, the thermal approximation is a good one. But at higher energies the spectrum cannot be approximated as thermal. Indeed, the spectrum has an ultraviolet cutoff at $`8GE=1`$, beyond which there is no radiation whatsoever. The precise expression, Eq. (23), is related to the change in de Sitter entropy. The entropy of three-dimensional de Sitter space is

$$S=\frac{2\pi r_H}{4G\mathrm{}},$$
(24)
where $`r_H`$ is the horizon radius. Once the shell has been emitted, its energy causes the horizon to shrink to the new radius $`r_H=l\sqrt{18GE}`$. Thus, Eq. (23), can be written as the exponent of the difference in the entropies, $`\mathrm{\Delta }S`$, before and after the particle has been emitted. This also explains the UV cutoff: the horizon cannot shrink past zero. Incidentally, Eq. (23) takes the same form as the corresponding expression for Hawking radiation from a Schwarzschild black hole , once back-reaction effects have been taken into account.

Indeed, one expects Eq. (23) on general grounds. For whatever the ultimate form of the holographic description of de Sitter space, quantum field theory tells us (via Fermi’s Golden Rule) that the rate must be expressible as

$$\mathrm{\Gamma }(if)=|_{fi}|^2\left(\text{phase space factor}\right),$$
(25)
where the first term on the right is the square of the amplitude for the process. The phase space factor is obtained by summing over final states and averaging over initial states. But the number of final states is just the final exponent of the final entropy, while the number of initial states is the exponent of the initial entropy. Hence

$$\mathrm{\Gamma }\frac{e^{S_{\mathrm{final}}}}{e^{S_{\mathrm{initial}}}}=\mathrm{exp}(\mathrm{\Delta }S)$$
(26)
in agreement with our result.

We end this section by noting that, although we have derived de Sitter radiance directly as particles tunneling across the de Sitter horizon, an alternate viewpoint is also possible. In this view, there is no de Sitter horizon and the spacetime is cutoff by a membrane living at the horizon. An observer then interprets the radiation not as tunneling particles, but rather as the spontaneous emissions of the membrane. In  it was shown that the classical equations of motion of the membrane can be derived from an action, and that Euclideanizing this action yields the correct entropy. It would be interesting to see whether the radiation can also be understood in this language.

## 4 On the Thermal Stability of de Sitter Space

The signs in Eq. (23) have the consequence that three-dimensional de Sitter space is thermally stable. For after a particle has been emitted, the new horizon radius is smaller than what it had previously been. The probability for emission of a second particle is now

$$\mathrm{\Gamma }_2\mathrm{exp}\mathrm{\Delta }S=\mathrm{exp}\left(\frac{\pi l}{2G\mathrm{}}\left(\sqrt{18G(E_1+E_2)}\sqrt{18GE_1}\right)\right),$$
(27)
where $`E_1`$ and $`E_2`$ are the energies of the two particles. For small $`E_2`$, this is

$$\mathrm{\Gamma }_2\mathrm{exp}\left(\frac{2\pi l}{\sqrt{18GE_1}}\frac{E_2}{\mathrm{}}\right).$$
(28)
The effective temperature which governs the emission probability of the second particle is therefore

$$T_2=\frac{\mathrm{}\sqrt{18GE_1}}{2\pi l},$$
(29)
which is lower than before. This is to be contrasted with the situation for Schwarzschild black holes, for which there is a runaway explosion as the black hole becomes smaller. In that respect, de Sitter space more closely resembles charged black holes. Note also that the change in the temperature of de Sitter space takes place because of the matter inside it, and not by any change in the cosmological constant, which of course remains constant throughout.

The decrease in the effective temperature holds out the possibility of thermally stabilizing de Sitter space . We started by considering empty de Sitter space. As this radiates and lowers its temperature, the horizon volume fills up with radiation. Eventually, the radiation passes through the origin and leaves the horizon; this causes the horizon radius to increase again, raising the temperature. By detailed balance, a stable thermodynamic equilibrium can be reached.

In higher dimensions ($`D>3`$), the story is somewhat different. The presence of the de Sitter radiation can now lead to the formation of a black hole, which in turn will radiate Hawking radiation outwards. Equilibrium is presumably reached only when the de Sitter and black hole horizons are coincident.

Acknowledgments

It is a pleasure to thank Daniel Kabat, Ivo Savonije, and especially Erik Verlinde for discussions. The author is supported by DOE grant DF-FCO2-94ER40818.


format: text

text:
# Global m=1 modes and migration of protoplanetary cores in eccentric protoplanetary discs

## 1 Introduction

The recent discovery of a number of extrasolar giant planets orbiting around nearby solar–type stars ( Mayor & Queloz, 1995, Marcy & Butler 1998, 2000) reveals that they have masses that are comparable to that of Jupiter, have orbital semi-major axes in the range

$`0.04\mathrm{AU}\stackrel{<}{}a\stackrel{<}{}\mathrm{\hspace{0.33em}\hspace{0.33em}2.5}\mathrm{au}`$, and orbital eccentricities in the range $`0.0\stackrel{<}{}e\stackrel{<}{}\mathrm{\hspace{0.33em}\hspace{0.33em}0.93}.`$

Orbital migration originally considered as a mechanism operating in protoplanetary discs by Goldreich & Tremaine (1980) has been suggested as an explanation for the existence of giant planets close to their central star (eg. Lin, Bodenheimer, & Richardson 1996). Protoplanetary cores are thought to form at several astronomical units and then migrate inwards either before accumulating a gaseous envelope and while in the earth mass range (type I migration) or in the form of a giant planet (type II migration) (see for example Lin & Papaloizou 1993, Ward 1997). In the former case the interaction with the disc is treated in a linearized approximation ( eg. Goldreich & Tremaine 1980) while in the latter nonlinearity is important leading to gap formation ( eg. Bryden et al 1999, Kley 1999). Both types of migration heve been found to occur on a timescale more than one order of magnitude shorter than that required to form giant planets or the expected lifetime of the disc and accordingly the survival of embryo protoplanets is in question ( eg. Ward 1997, Nelson et al 2000 ). Accordingly mechanisms that slow or halt migration such as the entry into a magnetospheric cavity close to the star ( Lin, Bodenheimer, & Richardson 1996) have been suggested. However, the existence of giant planets over a range of semi-major axes suggests that a more general mechanism for halting migration should be sought.

In this paper we shall consider effects produced when the large scale motion of the disc gas deviates from that of pure circular motion about the central mass as could be produced when the disc becomes eccentric as happens when it supports a global $`m=1`$ mode. Such modes could be produced by disc protoplanet interactions (eg. Papaloizou, Nelson & Masset 2001). However, for standard parameters, it was found that an instability leading to an eccentric disc together with an eccentric protoplanet orbit occurred only for large masses exceeding about fifteen Jupiter masses.

In this context note these calculations were done for laminar model discs with anomalously large viscosity coefficient and results may be somewhat different for turbulent discs. These simulations as well as resonant torque calculations appropriate to embedded cores by Papaloizou & Larwood (2000) have indicated possible reversal of orbital migration associated with eccentric orbits with modest eccentricity comparable to the disc aspect ratio.

This leads us to study further the dynamics of protoplanetary cores, massive enough for tidal interactions with the disc to be more important than effects due to gas drag, in an eccentric protoplanetary disc. In this paper we consider large scale $`m=1`$ modes in gaseous protoplanetary discs that correspond to making them eccentric. These modes have a low frequency branch for which the disc gas follows trajectories differing from Keplerian eccentric orbits by small corrections depending on forces due to disc self-gravity and pressure. These modes can be global in that they may vary on a length scale comparable to that of the whole disc even though it might have a large dynamic range. We shall consider a ratio of outer to inner radius of one hundred. These modes are of interest because, even though no definitive excitation mechanism of general applicability has yet been identified, their large scale implies a long life time comparable to the viscous time of the disc making them of potential interest in Astrophysics ( eg. Ogilvie 2001). They have also been considered as a potential source of angular momentum transport by Lee & Goodman (1999) in a tight winding approximation. Furthermore a disc composed of many stars on near Keplerian orbits has been postulated to occur in such objects as the nucleus of M31 (Tremaine 1995).

Here we consider global $`m=1`$ modes for various disc models neglecting viscous processes which are presumed to act over a longer time scale than that appropriate to the phenomena of interest. For global disturbances in protoplanetary discs of the type we consider, inclusion of self-gravity is important. Even though the discs are gravitationally stable, pressure and self-gravity can be equally important on scales comparable to the current radius, $`r,`$ when the Toomre stability parameter $`Qr/H,`$ with $`H`$ being the disc semi-thickness. This condition is satisfied for typical protostellar disc models ( eg. Papaloizou & Terquem 1999). In addition to this, disc self-gravity has to be considered in the description of the motion of embedded protoplanets and it generally causes them to move in eccentric orbits, with eccentricity comparable to or exceeding that in the disc. We find that circularization due to tidal interaction with the disc may play only a minor role if the local test particle precession frequency relative to the disc is large compared to the circularization rate. We also find that for modest disc eccentricities comparable to the aspect ratio $`H/r,`$ disc tidal interactions may differ significantly from those found in axisymmetric discs. This in turn may have important consequences for estimates of orbital migration rates for protoplanetary cores in the earth mass range.

In section 2 we give the Basic equations and the linearized form we use to calculate global $`m=1`$ modes with low pattern speed corresponding to introducing a finite disc eccentricity. We go on to describe the disc models used which may contain protoplanets orbiting in an inner cavity. In section 3) we present the results of normal mode calculations. We go on to discuss the motion of a protoplanet in the earth mass range in an eccentric disc in section 4, determining the equilibrium (non precessing) orbits which maintain apsidal alignment with the disc gas orbits. We then formulate the calculation of the tidal response of an eccentric disc to a low mass protoplanet in section 5 determining the time rate of change of the eccentricity and orbital migration rate. We find that aligned orbits with very similar eccentricity to that of the gas disc may suffer no eccentricity change while undergoing inward migration in general. However, when the non precessing aligned orbit has a significantly higher eccentricity than the disc, as can occur for modes with very small pattern speed, orbital migration may be significantly reduced or reverse from inwards to outwards for the disc models we consider. This finding is supported by a local dynamical friction calculation applicable when the protoplanet eccentricity is much larger than the disc aspect ratio. Finally we go on to discuss our results in section 7.

## 2 Basic equations and their linearized form for large scale global modes with $`m=1.`$

We work in a non rotating cylindrical coordinate system $`(r,\phi )`$ which may initially be considered to be centred on the primary star. The basic equations of motion in a two dimensional approximation which should be appropriate for a large scale description of the disc are taken to be

$$\frac{v_r}{t}+v_r\frac{v_r}{r}+\frac{v_\phi }{r}\frac{v_r}{\phi }\frac{v_\phi ^2}{r}=\frac{1}{\mathrm{\Sigma }}\frac{\mathrm{\Pi }}{r}\frac{\mathrm{\Phi }}{r},$$
(1)
$$\frac{v_\phi }{t}+v_r\frac{v_\phi }{r}+\frac{v_\phi }{r}\frac{v_\phi }{\phi }+\frac{v_\phi v_r}{r}=\frac{1}{r\mathrm{\Sigma }}\frac{\mathrm{\Pi }}{\phi }\frac{1}{r}\frac{\mathrm{\Phi }}{\phi },$$
(2)
In addition we have the two dimensional form of the continuity equation

$$\frac{(r\mathrm{\Sigma })}{t}+\frac{\left(r\mathrm{\Sigma }v_r\right)}{r}+\frac{\left(r\mathrm{\Sigma }v_\phi \right)}{\phi }=0.$$
(3)
Here the velocity is $`𝐯=(v_r,v_\phi ),`$ $`\mathrm{\Pi }=_{\mathrm{}}^{\mathrm{}}P𝑑z`$ represents a vertically integrated pressure which we assume to be a function of the surface density $`\mathrm{\Sigma }`$ defined by

$$\mathrm{\Sigma }=_{\mathrm{}}^{\mathrm{}}\rho 𝑑z.$$
(4)
The sound speed is then given by $`c=\sqrt{d\mathrm{\Pi }/d\mathrm{\Sigma }}.`$ The gravitational potential $`\mathrm{\Phi }=GM_{}/r+\mathrm{\Phi }_D+\mathrm{\Phi }_{ext},`$ with $`G`$ being the gravitational constant, has a point mass contribution arising from the central mass $`M_{},`$ a contribution, $`\mathrm{\Phi }_D,`$ due to the disc and a contribution, $`\mathrm{\Phi }_{ext},`$ due to external protoplanets.

The unperturbed disk is axisymmetric with no radial motion such that the velocity $`𝐯=(0,r\mathrm{\Omega })`$ with $`\mathrm{\Omega }>0.`$ In equilibrium we then have from equation (1)

$$\mathrm{\Omega }^2r=\frac{1}{\mathrm{\Sigma }}\frac{\mathrm{\Pi }}{r}+\frac{\mathrm{\Phi }}{r}$$
(5)
and the gravitational potential due to the disc is given by

$$\mathrm{\Phi }_D=G\mathrm{\Sigma }K(r,r^{})r^{}𝑑r^{},$$
(6)
with

$$K_0(r,r^{})=_0^{2\pi }\frac{1}{\sqrt{(r^2+r^22rr^{}\mathrm{cos}(\phi ))}}𝑑\phi .$$
(7)
We allow for the effects of orbiting external protoplanets which provide the external potential $`\mathrm{\Phi }_{ext}.`$ In this paper we are primarily interested in phenomena which vary on a time scale long compared to a local orbital period. Accordingly we use the time averaged or secular protoplanet perturbing potential which is derived in Appendix 1. For a single protoplanet of mass $`m_p`$ and orbital radius $`r_p,`$ the contribution to the external potential is given by

$$\mathrm{\Phi }_{ext}=G\frac{m_p}{2\pi }K_0(r,r_p).$$
(8)
In a thin disc of the kind considered here, the contributions due to pressure and self-gravity are comparable and small leading to nearly Keplerian rotation which has $`\mathrm{\Omega }r^{3/2}.`$ When $`cr^{1/2},`$ the disc then has a nearly constant putative aspect ratio $`H/r=c/(r\mathrm{\Omega }),`$ $`H`$ being the putative semi thickness.

### 2.1 Linearization

Here we are interested in global $`m=1`$ modes with slowly varying pattern when viewed in the adopted reference frame. To study these, we linearize the basic equations about the equilibrium state denoting perturbations to quantities with a prime.

As usual we assume that the dependence of all perturbations (in cylindrical coordinates) on $`\phi `$ and $`t`$ is through a factor $`\mathrm{exp}i\left(m\left(\phi \mathrm{\Omega }_pt\right)\right).`$ For the slowly varying modes we consider $`\mathrm{\Omega }_p<<\mathrm{\Omega }.`$

The linearized forms of equations (1) ,(2) and (3) are then

$$im(\mathrm{\Omega }\mathrm{\Omega }_p)v_r^{}2\mathrm{\Omega }v_\phi ^{}=\frac{W}{r},$$
(9)
$$im(\mathrm{\Omega }\mathrm{\Omega }_p)v_\phi ^{}+\frac{\kappa ^2}{2\mathrm{\Omega }}v_r^{}=\frac{imW}{r},$$
(10)
and

$$\frac{im(\mathrm{\Omega }\mathrm{\Omega }_p)\mathrm{\Sigma }(W\varphi ^{})}{c^2}=\frac{1}{r}\frac{\left(r\mathrm{\Sigma }v_r^{}\right)}{r}\frac{im\mathrm{\Sigma }v_\phi ^{}}{r}.$$
(11)
Here $`W=\mathrm{\Pi }^{}/\mathrm{\Sigma }+\mathrm{\Phi }^{}=\mathrm{\Sigma }^{}c^2/\mathrm{\Sigma }+\mathrm{\Phi }^{},`$

and $`\kappa ^2=(2\mathrm{\Omega }/r)(d(r^2\mathrm{\Omega })/dr)`$ denotes the square of the epicyclic frequency.

For the gravitational potential perturbation we have $`\mathrm{\Phi }^{}=\mathrm{\Phi }_D^{}+\mathrm{\Phi }_{ext}^{}.`$ Here we shall allow for the contribution of the secular effects of other sources such as protoplanets to the modes through the gravitational potential perturbation $`\mathrm{\Phi }_{ext}^{}.`$

### 2.2  Slowly varying modes with m=1

We consider linear perturbations of an axisymmetric disc that correspond to normal modes with $`m=1`$ that make it eccentric. The modes we consider are such that $`\mathrm{\Omega }_p<<\mathrm{\Omega }`$ and the dominant motion is epicyclic. In this case, to lowest order, we may neglect the pressure and self-gravity term $`W`$ as well as $`\mathrm{\Omega }_p`$ and assume $`\mathrm{\Omega }`$ is keplerian, to obtain from equation (10)

$$iv_\phi ^{}=\frac{1}{2}v_r^{}.$$
(12)
Using the above and introducing the radial Lagrangian displacement $`\xi _r=iv_r^{}/(\mathrm{\Omega }\mathrm{\Omega }_p),`$ equation (11) gives the surface density perturbation in the low frequency limit in the form

$$\mathrm{\Sigma }^{}=r\frac{d(\mathrm{\Sigma }e(r))}{dr},$$
(13)
with $`e(r)\xi _r/r`$ being the disc eccentricity.

The gravitational potential perturbation induced by the disc is

$$\mathrm{\Phi }_D^{}=G\mathrm{\Sigma }^{}K_1(r,r^{})r^{}𝑑r^{},$$
(14)
with

$$K_1(r,r^{})=_0^{2\pi }\frac{\mathrm{cos}(\phi )}{\sqrt{(r^2+r^22rr^{}\mathrm{cos}(\phi ))}}𝑑\phi \frac{\pi r}{r^2},$$
(15)
where the second term corresponds to the indirect term which results from the acceleration of the coordinate system produced by the disc material.

The existence of a mode with $`m=1`$ in the disc causes the orbits of the protoplanets interior to the disc to become eccentric. Noting the different form of indirect term used, the secular perturbing potential derived in Appendix 1 then gives the contribution of an orbiting protoplanet to the external potential perturbation to be

$$\mathrm{\Phi }_p^{}(r)=\frac{Gm_p\xi _r(r_p)}{2\pi r_p^2}\frac{}{r_p}\left(K_1(r,r_p)r_p^2\right).$$
(16)
Here the eccentricity of the companion orbit is related to the displacement by $`e_p=\xi _r(r_p)/r_p.`$ The total external potential perturbation is then found by summing over the perturbing protoplanets:

$$\mathrm{\Phi }_{ext}^{}(r)=\underset{p}{}\mathrm{\Phi }_p^{}(r)$$
(17)
Eliminating $`v_\phi ^{}`$ from equations (9) and (10), using equation (13) for $`\mathrm{\Sigma }^{}`$ and expanding to first order in the small frequencies $`\mathrm{\Omega }_p`$ and $`\omega _p=\mathrm{\Omega }\kappa ,`$ the fluid element orbital precession frequency in the absence of surface density or gravitational potential perturbations, one obtains a normal mode equation relating $`e(r),e_p`$ and $`\mathrm{\Omega }_p`$ in the form (see also Papaloizou, Nelson & Masset, 2001)

$$2\left(\mathrm{\Omega }_p\omega _p\right)\mathrm{\Omega }r^3e(r)=\frac{d}{dr}\left(\frac{r^3c^2}{\mathrm{\Sigma }}\frac{d[\mathrm{\Sigma }e(r)]}{dr}\right)\frac{d\left(r^2\mathrm{\Phi }^{}\right)}{dr},$$
(18)
where $`,c,`$ is the local sound speed.

In addition we have an equation for each protoplanet of the form

$$2\left(\mathrm{\Omega }_p\omega _p\right)\mathrm{\Omega }(r_p)r_p^3e_p=\left[\frac{d\left(r^2\mathrm{\Phi }^{}\right)}{dr}\right]_{r=r_p},$$
(19)
where of course for a particular protoplanet there is no self-interaction term in the sum for $`\mathrm{\Phi }_{ext}^{}.`$ The inclusion of self-gravity in the eigenvalue problem is essential if modes with prograde precession frequency are to be obtained. For typical protoplanetary disc models, self-gravity can be strong enough to induce prograde precession for the long wavelength $`m=1`$ modes considered here.

Normal modes were calculated by discretizing (18) and (19) and formulating a matrix eigenvalue problem on an unequally spaced grid with $`200`$ grid points with intervals increasing in geometric progression (see Terquem & Papaloizou 2000 for consideration of a related problem).

### 2.3  Disc models

We present here results obtained for a disc model with equation of state given by a polytrope of index $`n=1.5.`$ The sound speed is given by

$$c^2=\frac{GM_{}H^2}{R_{in}^3}\left(1\left(R_{in}/r\right)^{10}\right)\left(1\left(r/R_{out}\right)^{10}\right)R_{in}/r.$$
(20)
Here the first two factors provide sharp edges at the disc inner and outer boundaries and $`H/r=0.05`$ is the disc aspect ratio away from these boundaries.

The surface density is given by

$$\mathrm{\Sigma }=\mathrm{\Sigma }_0(c^2)^n.$$
(21)
The boundary regions are chosen to be of order the disc scale height in width and away from these $`\mathrm{\Sigma }r^{3/2}.`$ In the above $`R_{in}`$ is the inner boundary radius which is taken to be the unit of length. The unit of mass is the central mass $`M_{}`$ and the unit of time is $`R_{in}^{3/2}/\sqrt{GM_{}}.`$ The arbitrary scaling factor $`\mathrm{\Sigma }_0`$ was chosen such that in model $`A`$ the total disc mass was $`4.0\times 10^2M_{}.`$ In model $`B`$ the total disc mass was $`4.0\times 10^3.`$ $`R_{out}`$ is the outer boundary radius, here, in order to study large scale modes, taken as $`100R_{in}.`$

The importance of self gravity is measured by the Toomre parameter $`Q=\mathrm{\Omega }c/(\pi G\mathrm{\Sigma }).`$ For model A this has a minimum value of $`5.2`$ while for the lower mass disc model of model B this minimum value is $`52.`$ Another quantity of interest is the local precession frequency of a test particle orbit $`\omega _{pg}`$ in the axisymmetric component of the total gravitational potential ( see equation (29) below).

This is plotted for the two disc models in figure 1. It is small and negative in dimensionless units over much of the discs, corresponding to retrograde precession, scaling with the disc mass as it is determined by the disc self-gravity. As we shall see this behaviour provides scope for secular resonance associated with low frequency modes.

## 3 Normal modes

We have calculated the lowest order global normal modes for two different configurations involving protoplanets orbiting interior to the disc in an inner cavity. The first, approximating the Upsilon Andromedae system involves two protoplanets with mass ratios $`0.00383`$ and $`0.00196`$ orbiting at $`0.6R_{in}`$ and $`0.194R_{in}`$ respectively. However, the disc normal modes do not change much in character if the protoplanet orbital radii are scaled to somewhat smaller fractions of $`R_{in}.`$ But the ratio of protoplanet to disc eccentricity in the joint modes changes more significantly. The second configuration we consider is a single protoplanet with mass ratio $`0.002`$ orbiting at $`0.6R_{in}.`$ Finally we consider a disc with no interior orbiting protoplanets. In all of these cases we consider both model $`A`$ and model $`B`$ discs.

The pattern speeds for the highest frequency normal modes and the protoplanet orbital eccentricities occurring jointly with the normal modes are given in table 1. For the results presented here, the modes are normalized such that the disc eccentricity at the inner boundary is 0.1.

Many of the properties of these modes can be understood with reference to the local dispersion relation for density waves in the low frequency limit (Lin & Shu 1969) in the form

$$2\mathrm{\Omega }(\omega _p\mathrm{\Omega }_p)=2\pi G\mathrm{\Sigma }|k|+c^2k^2,$$
(22)
where $`k`$ is the radial wavenumber. The above indicates that when self-gravity, governed by the first term on the right hand side, is unimportant as occurs for either low mass discs or large $`|k|,`$ $`\mathrm{\Omega }_p`$ is negative. On the other hand positive $`\mathrm{\Omega }_p`$ can occur for low $`|k|`$ or massive discs. Thus in that case the lowest order modes (in terms of numbers of nodes) should be prograde with more of these existing for more massive discs, a trend we find in our results. There can also be very low frequency (in magnitude) modes for which the self-gravity and pressure terms approximately cancel. This occurs when $`|k|=2\pi G\mathrm{\Sigma }/c^22/(QH).`$ For $`Qr/H`$ as in the models here, this is comparable to the radius making a very global mode.

As we calculate normal modes jointly involving protoplanets and disc, the protoplanets have associated eccentricities. In general the $`𝒩`$ highest frequency modes predominantly involve the protoplanets while the others predominantly involve the disc, $`𝒩`$ being the number of protoplanets. The two highest frequency modes for the two protoplanet case are plotted in figure 2 for disc model $`A`$ and in figure 3 for disc model $`B.`$ For the second highest frequency mode calculated for disc model $`A`$, the protoplanet eccentricity ratio is $`1.38,`$ while for model $`B`$ it is $`1.88.`$ The latter result is similar to that given by Chiang, Tabachnik, & Tremaine (2001) who considered the current Upsilon Andromedae system with no disc and concluded it was predominantly in this mode.

The results found here suggest the disc and protoplanet orbits were antialigned. With the outermost protoplanet at $`0.6R_{in}`$ the disc inner boundary and outer protoplanet eccentricities are comparable for model $`A`$ while for model $`B`$ the relative disc eccentricity is $`56`$ percent smaller. However, when the outer protoplanet orbits at $`0.5R_{in}`$ the disc inner boundary eccentricity is only $`1/3`$ that of the protoplanet for model $`A`$ and $`19`$ percent for model $`B.`$ Apart from the eccentricity scaling relative to the interior protoplanets, the spatial form of the eigenfunction in the disc remains almost identical.

Apart fom the two modes with highest frequency, other normal modes are essentially pure disc modes with only small associated protoplanet eccentricities.

The four highest frequency modes for the one protoplanet case with disc model $`A`$ are plotted in figure 4 and in figure 5 for model $`B.`$ Pure disc modes with no protoplanet are plotted in figure 6 for model $`A`$ and figure 7 for model $`B`$. The modes develop increasing numbers of nodes as their frequency decreases but the mode with lowest frequency in absolute magnitude can be very global with only a few nodes out to $`r=R_{out}.`$ We comment that the modes are rather non compressive requiring eccentricities of order unity to provide Lagrangian changes in surface density of order unity. As the motion in the modes in linear theory is assumed epicyclic, this is suggestive that the analysis should be valid as long as the epicyclic approximation is.

Having calculated these modes that may be largely driven by protoplanets in eccentric orbits or possibly self-excited or long lived structures we now go on to consider the migration of low mass protoplanets embedded in such eccentric discs.

## 4  The motion of a low mass protoplanet in an eccentric disc

The evolution of the coplanar orbit of an interior protoplanet in the earth mass range embedded in the disc due to the gravitational interaction with the slowly precessing disc can be found using secular perturbation theory. We regard the embedded protoplanet as a test particle evolving in a prescribed gravitational potential and neglect the influence of the protoplanet on the disc mode. This should be reasonable when as here, the angular momentum of the protoplanet is significantly less than that of the material supporting the normal mode ( eg. Papaloizou et al 2001).

The orbit evolves under the perturbing potential per unit mass,

$$𝒱\mathrm{\Phi }_D+\mathrm{\Phi }_{ext}=\mathrm{\Phi }_0(r)+\mathrm{\Phi }_1(r)\mathrm{cos}(\phi \mathrm{\Omega }_pt).$$
(23)
Here $`\mathrm{\Phi }_0(r)`$ is the axisymmetric component of the potential but excluding that due to the central mass and $`\mathrm{\Phi }_1(r)`$ is the radial amplitude of the $`m=1`$ component arising from the normal mode.

On performing a time average over the protoplanet orbit one obtains the Hamiltonian system:

$$\frac{dh}{dt}=\frac{}{\alpha },\frac{d\alpha }{dt}=\frac{}{h}.$$
(24)
Here $`\alpha =\varpi \mathrm{\Omega }_pt,`$ with $`\varpi `$ being the longtitude of periapse and $`=𝒱\mathrm{\Omega }_ph.`$ The specific angular momentum is $`h=\sqrt{GM_{}a_p(1e_p^2)},`$ $`a_p`$ is the constant semi-major axis, and $`e_p`$ is the eccentricity. To leading powers of $`e_p`$ we have

$$=A(a_p)+B(a_p)e_p^2\mathrm{\Omega }_ph+C(a_p)e_p^4+e_pD(a_p)\mathrm{cos}\alpha .$$
(25)
Hamilton’s equations then give to leading powers of $`e_p`$

$$\frac{d\alpha }{dt}=\frac{2(B+e_p^2(2CB/2))}{\sqrt{GM_{}a_p}}\mathrm{\Omega }_p\frac{D(a_p)\mathrm{cos}\alpha }{\sqrt{GM_{}a_p}e_p}.$$
(26)
$$\frac{de_p}{dt}=\frac{\sqrt{1e_p^2}}{\sqrt{GM_{}a_p}e_p}\frac{}{\alpha }=\frac{D(a_p)}{\sqrt{GM_{}a_p}}\mathrm{sin}\alpha .$$
(27)
Here the precession rate of the apsidal line induced by the axisymmetric potential, to first order in the perturbing potential, expressed as an expansion in powers of $`e_p`$ is given by the first term on the right hand of equation(26) as

$$\omega _{pg}=\frac{2}{\sqrt{GM_{}a_p}}\left(B+e_p^2(2CB/2)\right).$$
(28)
This may also be written in terms of the axisymmetric component of the potential directly as

$$2\omega _{pg}\sqrt{GM_{}a_p^5}=(1+e_p^2)\frac{d^2\mathrm{\Phi }_0}{du^2}+\frac{e_p^2}{8a_p^2}\frac{d^4\mathrm{\Phi }_0}{du^4},$$
(29)
where $`u=1/r,`$ and $`r`$ is taken to be equal to $`a_p.`$

Equilibrium or steady state solutions of (26) and (27) with $`e_p`$ and $`\alpha `$ constant occur when $`\alpha =0`$ or $`\alpha =\pi .`$ Adopting the convention that $`e_p`$ is positive suffices to select one of the latter possibilities as (26) gives an expression from which $`e_p`$ can be determined in the form

$$e_p=\frac{D(a_p)\mathrm{cos}\alpha }{2(B+e_p^2(2CB/2))+\mathrm{\Omega }_p\sqrt{GM_{}a_p}}.$$
(30)
or equivalently

$$e_p=\frac{D(a_p)\mathrm{cos}\alpha }{(\omega _{pg}\mathrm{\Omega }_p)\sqrt{GM_{}a_p}}.$$
(31)
Alternatively one may fix $`\alpha `$ to be $`0`$ or $`\pi `$ and allow $`e_p`$ to change sign as we have done for the normal modes.

An equilibrium solution so determined corresponds to the situation when the protoplanet orbit precesses at the same rate, $`\mathrm{\Omega }_p,`$ as the mass distribution that produces the gravitational potential while maintaining a constant eccentricity. Assuming that $`\omega _{pg}`$ and $`\mathrm{\Omega }_p`$ are of comparable magnitude, in general for modest eccentricity,

$$e_pD(a_p)/(\mathrm{\Omega }_p\sqrt{GM_{}a_p})e_0$$

and is proportional to the magnitude of the nonaxisymmetric potential. An exception occurs when

$$2B+\mathrm{\Omega }_p\sqrt{GM_{}a_p}=0.$$

In this case, the precession frequency of a free orbit with small eccentricity, $`\omega _{pg},`$ matches that of the nonaxisymmetric mass distribution $`\mathrm{\Omega }_p`$ corresponding to a secular resonance. When this occurs equation (30) indicates that $`e_p(e_0)^{1/3}.`$ Accordingly significantly larger equilibrium eccentricities are expected close to a secular resonance. We comment that for fixed $`\alpha `$ $`e_p`$ changes sign as a secular resonance is passed through corresponding to an alignment change through a rotation of the axis of the ellipse through $`\pi .`$

One may also investigate the effect of orbital circularization by adding a term $`e_p/|t_e|`$ to the right hand side of (27), where $`|t_e|`$ is the circularization time. In this case one can still find an equilibrium but with orbital apsidal line rotated. Restricting consideration to the situation away from secular resonance, one finds that the equilibrium eccentricity is reduced by a factor $`\sqrt{1+1/((\omega _{pg}\mathrm{\Omega }_p)^2t_e^2)}`$ and $`|\mathrm{sin}\alpha |=1/\sqrt{1+(\omega _{pg}\mathrm{\Omega }_p)^2t_e^2}.`$ Thus when $`|(\omega _{pg}\mathrm{\Omega }_p)t_e|`$ is large, the effect of the circularization term is to produce a small rotation of the apsidal line of the orbit.

As indicated above equilibrium solutions correspond to the situation where the apsidal precession of the protoplanet orbit is locked to that of the underlying nonaxisymmetric disc.

One can find solutions undergoing small librations in the neighbourhood of equilibrium solutions (eg. Brouwer & Clemence 1961) and when dissipative forces are added these may decay making the attainment of equilibrium solutions natural. As we indicate below, tidal interaction of a protoplanet with the disc may produce such dissipative forces. Accordingly we shall focus on equilibrium solutions in what follows below.

### 4.1 Determination of equilibrium eccentricities corresponding to normal modes

It is a simple matter to determine equilibrium eccentricities for protoplanets moving under the gravitational potential appropriate to a normal mode corresponding to an eccentric disc, with pattern speed $`\mathrm{\Omega }_p,`$ of the type calculated above. We recall the normal mode equation(18) which can be regarded as determining the equilibrium disc eccentricity in the form

$$2\left(\mathrm{\Omega }_p\omega _p\right)\mathrm{\Omega }r^3e(r)=\frac{d}{dr}\left(\frac{r^3c^2}{\mathrm{\Sigma }}\frac{d[\mathrm{\Sigma }e(r)]}{dr}\right)\frac{d\left(r^2\mathrm{\Phi }^{}\right)}{dr}.$$
(32)
If the pressure forces are neglected by dropping the terms involving $`c^2,`$ this is exactly equivalent to equation (31) for determining the protoplanet eccentricity with $`e(r)`$ corresponding to $`e_p`$ and $`\omega _p`$ replaced by $`\omega _{pg}.`$ Thus a local protoplanet equilibrium eccentricity appropriate to a disc mode or response can be determined by using the disc normal mode equation (18) retaining the term resulting from the disc self-gravity but omitting pressure contributions.

This procedure should be applicable provided the eccentricities are not too large. Equilibrium eccentricities are plotted for some of the normal modes we calculated in figures 2 and 3 and also 8 \- 11. In figures 2 and 3 we also plot equilibrium eccentricities but calculated with the nonaxisymmetric contribution to the gravitational potential due to the disc removed which is equivalent to assuming that it remains circular. We see that the assumption that disc remains circular results in a significantly smaller equilibrium eccentricity for disc model A, particularly for the mode with second highest frequency. For the lower mass disc model B, differences arising fom assuming the disc remains circular are less pronounced. We also coment that the occurence of secular resonances in the inner parts of the disc where $`r/R_{in}<10`$ is common.

One expects that dissipative torques produced by

protoplanet-disc tidal interaction may result in the approach to such equilibrium solutions from general initial conditions. This we now discuss.

## 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 $`\phi =0,`$ $`a=r(1+e\mathrm{cos}\phi )/(1e^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)e(r).`$ In this approximation $`a`$ can be expressed in terms of cylindrical coordinates in the form $`a=r(1+e\mathrm{cos}\phi )`$ and the orthogonal angular coordinate, $`\lambda `$ is given by

$$\lambda =\phi +(\frac{e}{r}dr)\mathrm{sin}\phi ,$$
(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_a,v_\lambda )`$ may then be written ( see Appendix 2)

$$\frac{v_a}{t}+v_a|a|^2\frac{}{a}\left(\frac{v_a}{|a|}\right)+v_\lambda |a||\lambda |\frac{\left(\frac{v_a}{|a|}\right)}{\lambda }$$

$$\frac{v_\lambda ^2}{2}|a||\lambda |^2\frac{\left(|\lambda |^2\right)}{a}\frac{v_a^2}{2}|a|^3\frac{\left(|a|^2\right)}{a}$$

$$=\frac{|a|}{\mathrm{\Sigma }}\frac{\mathrm{\Pi }}{a}|a|\frac{\mathrm{\Phi }}{a},$$
(34)
$$\frac{v_\lambda }{t}+v_\lambda |\lambda |^2\frac{\left(\frac{v_\lambda }{|\lambda |}\right)}{\lambda }+v_a|a||\lambda |\frac{\left(\frac{v_\lambda }{|\lambda |}\right)}{a}$$

$$\frac{v_a^2}{2}|a|^2|\lambda |\frac{\left(|a|^2\right)}{\lambda }\frac{v_\lambda ^2}{2}|\lambda |^3\frac{\left(|\lambda |^2\right)}{\lambda }$$

$$=\frac{|\lambda |}{\mathrm{\Sigma }}\frac{\mathrm{\Pi }}{\lambda }|\lambda |\frac{\mathrm{\Phi }}{\lambda },$$
(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_a=0,`$ and to first order in $`e,`$ $`v_\lambda =(GM_{}/a)^{1/2}(1+e\mathrm{cos}\lambda ).`$ In a thin disc forces due pressure and self-gravity provide only a small correction of order $`(H/r)^2.`$ We consider the situation when the disc eccentricity, $`e(r),`$ protoplanet eccentricity, $`e_p,`$ and $`H/r`$ can 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 $`mr/H.`$ Here for the time being $`\lambda `$ is regarded as the azimuthal angle. We perform a response calculation by linearizing equations (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 $`\mathrm{\Omega }v_a^{}.`$ In this scheme we may replace $`|a|`$ and $`|\lambda |`$ by unity except where the latter occurs in the combination $`|\lambda |v_\lambda \frac{}{\lambda },`$ with $`|v_\lambda |`$ 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 $`me|v_\lambda |v_a^{}/r`$ in the linearized equations. For $`mr/H,1/e`$ these are comparable to the dominant terms. All other eccentricity contributions are smaller by a factor at most $`e.`$

We have to first order in eccentricity

$$v_\lambda |\lambda |=(GM_{}/a^3)^{1/2}\left(1+2e\mathrm{cos}\phi +(\frac{e}{r}dr)\mathrm{cos}\phi \right).$$
(36)
Thus if we define a new angle $`M`$ through $`M=\lambda (2e+(\frac{e}{r}dr)\mathrm{sin}\lambda ,`$ the operator $`|\lambda |v_\lambda \frac{}{\lambda },`$

becomes $`=(GM_{}/a^3)^{1/2}\frac{}{M}.`$ The angle $`M`$ is just the mean anomaly as can be seen from the fact that to first order in $`e,`$

$$r=a(1e\mathrm{cos}M),\phi =M+2e\mathrm{sin}M,$$
(37)
For the motion of a particular disc particle we have $`M=\mathrm{\Omega }(tt_0),`$ where $`t_0`$ denotes the time of periastron passage.

To lowest order in $`e,`$ the coordinates $`(a,M)`$ behave just like the cylindrical coordinates $`(r,\phi )`$ in the sense that one may make the replacements $`ra,\phi \lambda `$ in the standard linearized equations expressed in cylindrical coordinates. However, recall that in the forcing potential $`\mathrm{\Phi }_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 $`mr/H`$ of interest, the protoplanet forcing potential is

$$\mathrm{\Phi }_p^{}=\frac{Gm_p}{\sqrt{|r^2+R^22rR\mathrm{cos}(\phi \nu )+b^2|}}.$$
(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 $`bH`$ ( eg. Artymowicz 1993).

For an eccentric protoplanet orbit with semi-major axis $`a_p`$ and eccentricity $`e_p,`$ we have

$$R=a_p(1e_p\mathrm{cos}(\omega t)),\nu =\omega t+\varpi +2e_p\mathrm{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

$$\mathrm{\Phi }_p^{}=\underset{m,n,n_1}{}C_{m,n,n_1}\mathrm{exp}i\left((nm)\omega t+m(M\varpi )+n_1M\right).$$

Here, the sum is over positive values of $`m`$ and both positive and negative values of $`n`$ and $`n_1.`$ The convention is that the real part is to be taken to give the physical potential. The coefficients $`C_{m,n,n_1}`$ are real and given by

$$\frac{4\pi ^3C_{m,n,n_1}}{Gm_p\omega }=$$

$$\frac{\mathrm{exp}i\left((nm)\omega tm\varpi +(n_1+m)M\right)}{\sqrt{|r^2+R^22rR\mathrm{cos}(\phi \nu )+b^2|}}𝑑M𝑑\varpi 𝑑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,n_1}`$ a function of $`a,e,a_p,`$ and $`e_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_1+m)=k_1`$ and $`(nm)=k_2`$.

Thus for such a particular forcing term with fixed $`(k_1,k_2)`$

$$\mathrm{\Phi }_p^{}=f_{k_1,k_2}\mathrm{exp}i\left(k_2\omega t+k_1M\right),$$
(41)
where

$$f_{k_1,k_2}=\underset{m,n,n_1}{}\delta _{n_1+m,k_1}\delta _{nm,k_2}C_{m,n,n1}\mathrm{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_1=k_2`$ and then $`\varpi `$ appears only as a redundant complex phase.

For the general forcing term (41) the pattern speed $`\mathrm{\Omega }_p=k_2\omega /k_1.`$ We shall consider the situation when the disc surface density $`\mathrm{\Sigma }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 (eg. Goldreich & Tremaine 1978). These occur when $`k_1(\mathrm{\Omega }\mathrm{\Omega }_p)=\pm \mathrm{\Omega }\sqrt{1+\xi ^2}`$ with $`\xi =k_1c/(r\mathrm{\Omega })`$ (eg. Artymowicz 1993). For $`k_1>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 $`(ar,M\phi )`$ , 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

$$\frac{dE_d}{dt}=\frac{\pi ^2\mathrm{\Sigma }r^2}{3\mathrm{\Omega }\sqrt{1+\xi ^2}(1+4\xi ^2)}\left|_{k_1,k_2}\right|^2,$$
(43)
where

$$_{k_1,k_2}=\frac{df_{k_1,k_2}}{dr}+\frac{2m^2(\mathrm{\Omega }\mathrm{\Omega }_p)f_{k_1,k_2}}{\mathrm{\Omega }r}.$$

However, unlike in an axisymmetric disc, the outward angular momentum flow in the disc associated with the protoplanet,
$`dJ_d/dt(1/\mathrm{\Omega }_p)dE_d/dt.`$

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

$$\frac{dJ_d}{dt}=\mathrm{\Sigma }^{}\frac{\mathrm{\Phi }_p^{}}{\varpi }r𝑑r𝑑\phi ,$$
(44)
the integral being taken over the disc. Recalculation of the torque formula then gives

$$\frac{dJ_d}{dt}=Re\left(\frac{\pi ^2\mathrm{\Sigma }r^2}{3\mathrm{\Omega }\omega |k_2|\sqrt{1+\xi ^2}(1+4\xi ^2)}_{k_1,k_2}𝒯_{k_1,k_2}^{}\right),$$
(45)
where

$$𝒯_{k_1,k_2}=\frac{dg_{k_1,k_2}}{dr}+\frac{2m^2(\mathrm{\Omega }\mathrm{\Omega }_p)g_{k_1,k_2}}{\mathrm{\Omega }r},$$

with

$$g_{k_1,k_2}=\underset{m,n,n_1}{}m\delta _{n_1+m,k_1}\delta _{nm,k_2}C_{m,n,n1}\mathrm{exp}(im\varpi ).$$
(46)
The rate of change of protoplanet eccentricity, $`e_p`$, induced by each term can be found by using the orbit equation

$$\frac{e_p\omega J}{(1e_p^2)^{3/2}}\frac{de_p}{dt}=\frac{dE}{dt}\omega \frac{dJ}{dt}(1e_p^2)^{1/2},$$
(47)
where $`E,`$ and $`J`$ are the orbital energy and angular momentum of the protoplanet and $`dJ/dt=\pm dJ_d/dt,dE/dt=\pm dE_d/dt,`$ 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_1,k_2).`$

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_e=\frac{e_p}{\mathrm{d}e_p/\mathrm{d}t}.$$
(48)
We define the time for the angular momentum to decay by a factor of $`e`$ as

$$t_J=\frac{J}{\mathrm{d}J/\mathrm{d}t},$$
(49)
such that a positive value means inward migration.

## 6 Tidal torques in an eccentric disc

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

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

We use the expressions for $`r,\phi ,R,`$ and $`\nu ,`$ given by equations (37), (39) in the forcing potential (38). In the case when $`e=e_p,`$ and the protoplanet orbit is almost aligned with the disc with very small $`|\varpi |,`$ the strongest interaction occurs when $`r=R,`$ and $`\phi =\nu .`$ This corresponds to $`a=a_p,`$ and $`M=\omega t+\varpi .`$ Performing a first order Taylor expansion about the point of maximum interaction the forcing potential becomes

$$\mathrm{\Phi }_p^{}=\frac{Gm_p}{\sqrt{|(aa_p)^2+4aa_p\mathrm{sin}^2(M\omega t\varpi )+b^2|}}.$$
(50)
Here $`|aa_p|`$ and $`a|M\omega t|`$ are considered as small compared to $`a`$ and terms of order $`e`$ times these small quantities have been neglected. Equation (50) is to be considered in comparison to the similar expression appropriate to an axisymmetric disc

$$\mathrm{\Phi }_p^{}=\frac{Gm_p}{\sqrt{|(rR)^2+4rR\mathrm{sin}^2(\phi \nu )+b^2|}}.$$
(51)
Given that we have already shown that $`a,M`$ behave just like cylindrical coordinates $`(r,\phi )`$ the tidal torque calculation in an eccentric disc with $`e_p=e`$ should give the same results as an axisymmetric disc with $`e_p=0.`$ Note too that only one term should remain in the sum (42) such that $`k_1=m,`$ and $`k_2=m.`$ Then there is zero rate of change of eccentricity. Thus the situation of aligned orbits such that $`e_p=e`$ behaves much like the case with $`e_p=0`$ in an axisymmetric disc in that it is one of steady eccentricity. The orbit then migrates at the same rate as a circular orbit in an axisymmetric disc. This is essentially what we have found on application of the torque formulae.

### 6.1 Numerical results

We have calculated $`t_J/N_t`$ and $`t_e/N_t`$ by summing the contributions from appropriate resonances. The normalizing factor

$$N_t=(r/5\mathrm{A}\mathrm{U})^{1/2}(m_p/(3M_{}))^1(\mathrm{\Sigma }/56)^1.$$

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

We plot $`t_e`$ for $`N_t=1`$ in figure 13. The trends in all cases are similar and are that for small $`e_p`$ the protoplanet eccentricity grows in the aligned case while inward migration occurs with $`t_J2\times 10^6`$ yr. The eccentricity growth reverses for $`e_pe`$ indicating an equilibrium in accordance with the discussion above. We comment that for very small $`e_p`$ and finite $`e,`$ $`de_p/dte.`$ For larger $`e_p,`$ $`t_J`$ and $`t_e`$ increase and $`t_J`$ eventually changes sign for $`e_p`$ exceeding $`0.1.`$ At these eccentricities $`|t_e|2\times 10^4`$ yr. We make the comment that the same dimensionless units can be used for $`t_J`$ and $`t_e`$ as for the disc models introduced in section 2 and thus the same scaling to make results applicable to different radii may be used.

In the antialigned case the disc protoplanet interaction is much weaker (see figures 12 and 13). Note too that the interaction with the disc weakens in general for larger $`e_p`$ because of the larger relative velocity of the protoplanet with respect to the disc. This results in larger values of $`|t_e|`$ and $`|t_J|.`$ Additional calculations have shown that, as expected, these values become independent of the orientation of the apsidal line when $`e_pe.`$

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

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

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

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

To give numerical examples we first consider the two protoplanet model with disc $`A`$ as an approximation to the Upsilon Andromedae system. Taking the situation represented in figure 2 with inner boundary disc eccentricity $`0.1,`$ for a semi-major axis of $`2.57`$ AU for the outer protoplanet in the inner cavity, being $`0.6`$ of the inner boundary radius, $`r=1.5`$ corresponds to $`6.4`$ AU. For the mode with $`\mathrm{\Omega }_p=6.99\times 10^4`$ and disc model $`A,`$ $`\mathrm{\Sigma }=136gmcm^2`$ at $`r=1.5.`$ At this radius $`e=0.07,`$ and the equilibrium eccentricity $`e_p=0.1.`$ From figure 12 the orbital migration rate is very small and possibly outwards. We also find $`|t_e|=10^4(3M_{}/m_p)`$ yr and $`|(\omega _{pg}\mathrm{\Omega }_p)t_e|30.0(M_{}/m_p).`$ Thus in this case the eccentric disc has significant effects on the tidal torques acting on an embedded protoplanet.

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

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

For the mode with $`\mathrm{\Omega }_p=4.02\times 10^6,`$ we find that at $`r=10`$ $`e_p=2.5e.`$ Supposing that $`r=10`$ corresponds to $`5`$ AU, $`\mathrm{\Sigma }=426gmcm^2.`$ From figure 12 we find that for $`e=0.05,`$ there is torque reversal. We also again find that $`|t_e|3\times 10^3(3M_{}/m_p)`$ yr and $`|(\omega _{pg}\mathrm{\Omega }_p)t_e|30.0(M_{}/m_p).`$

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

### 6.2 Dynamical friction calculation

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

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

$$\mathrm{\Phi }_p^{}(𝐫)=\mathrm{\Phi }_p^{}(𝐤)\mathrm{exp}(i\mathrm{𝐤𝐫})d^N𝐤.$$
(52)
Here we shall consider both the two dimensional case $`(N=2)`$ and the three dimensional case $`(N=3)`$. For the two dimensional case with $`\mathrm{\Phi }_p^{}(𝐫)=Gm_p/\sqrt{r^2+b^2},`$ $`\mathrm{\Phi }_p^{}(𝐤)=Gm_p\mathrm{exp}(bk)/(2\pi k),`$ with $`k=|𝐤|.`$ In the three dimensional case with $`\mathrm{\Phi }_p^{}(𝐫)=Gm_p/r,`$ $`\mathrm{\Phi }_p^{}(𝐤)=Gm_p/(2\pi ^2k^2).`$

In each case, assuming a local steady state, the velocity $`𝐯^{}`$ induced by the protoplanet is found from

$$𝐯𝐯^{}=\mathrm{\Phi }_p^{}.$$
(53)
We also have

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

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

$$\dot{𝐏}=D^{}\mathrm{\Phi }_p^{}d^N𝐫$$
(55)
Performing an integration by parts and working in terms of the Fourier transforms of perturbations, one obtains

$$\dot{𝐏}𝐯=(2\pi )^ND\mathrm{\Phi }_p^{}(𝐤)i𝐤𝐯^{}(𝐤)d^N𝐤.$$
(56)
From (53) one obtains

$$𝐯^{}(𝐤)𝐤=\frac{k^2\mathrm{\Phi }_p^{}(𝐤)}{𝐯𝐤}.$$
(57)
In order to perform the integral (56) one has to apply a Landau prescription by adding an infinitessimally small negative imaginary part to the denominator $`𝐯𝐤.`$ One then finds

$$\dot{𝐏}𝐯=\frac{\pi D(Gm_p)^2}{v}𝒬,$$
(58)
where for $`N=3,`$ $`𝒬=4\mathrm{ln}(k_{max}/k_{min}),`$ and for $`N=2,`$ $`𝒬=1/(2b).`$ Here $`(k_{max},k_{min})`$ are the usual upper and lower wavenumber cut offs (eg. Tremaine & Weinberg 1984). Here reasonable values are $`k_{max}=v^2/(Gm_p),`$ and $`k_{min}=1/H.`$ The two and three dimensional cases are thus of the same form. The logarithmic factor is generally of order unity. Comparison of the two and three dimensional cases suggests that $`\rho \mathrm{ln}(k_{max}/k_{min})=\mathrm{\Sigma }/(8b).`$ Thus the adopted softening parameter $`b`$ should be somewhat smaller than $`H.`$

Using (58) we may evaluate, remembering that $`𝐯`$ is the relative velocity between disc and protoplanet, the average rate of change of angular momentum of the protoplanet from

$$\frac{dJ}{dt}=\frac{\omega }{2\pi }\frac{\pi D(Gm_p)^2(𝐫\times 𝐯)\widehat{𝐤}}{v^3}𝒬𝑑t,$$
(59)
with the integral being taken round the orbit and $`\widehat{𝐤}`$ being the unit vector in the direction normal to the disc.

In the two dimensional case for $`D\mathrm{\Sigma }r^{3/2}`$ one obtains for small $`e_p`$ and $`e=0.`$

$$\frac{1}{J}\frac{dJ}{dt}\frac{\pi m_p\omega \mathrm{\Sigma }(a_p)a_p^3}{2be_pM_{}^2}.$$
(60)
As this is positive it corresponds to outward migration as long as the eccentricity can be maintained. This occurs because the disc flow tends to speed up the protoplanet at apocentre where most time is spent. Numerically for $`b=H/\sqrt{2}`$

$$\frac{t_J}{N_t}=2.6\times 10^6\left(\frac{H}{0.07a_p}\right)^2\left(\frac{e_pa_p}{H}\right)\mathrm{yr}.$$
(61)
While this agrees in form and is of comparable magnitude to what is obtained from resonant torque calculations (see figure 12) it is impractical to perform the latter at higher eccentricities where better agreement might be attained because of the small scale of the interaction (even at the eccentricities plotted, over $`10^6`$ resonances were included).

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

$$\frac{dE}{dt}=\frac{\omega }{2\pi }\frac{\pi D(Gm_p)^2𝐯(𝐯GM_{}r^{1/2}\widehat{𝝋})}{v^3}𝒬𝑑t,$$
(62)
with $`E=\frac{GM_{}m_p}{2a_p}`$ and $`\widehat{𝝋}`$ being the unit vector in the azimuthal direction.

Performing the integration we find for small $`e_p`$ and assuming $`\mathrm{\Sigma }r^n`$ that

$$\frac{1}{E}\frac{dE}{dt}\frac{m_p\omega \mathrm{\Sigma }(a_p)a_p^3}{2be_pM_{}^2}\left(4.167.17n\right).$$
(63)
From this we see that, for this circular disc case, the mean rate of change of orbital energy increases for $`n<0.6.`$ This simplified calculation indicates outward migration for density profiles that do not increase too rapidly inwards in line with the idea that the effect is caused by the interaction at apocentre.

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

## 7 Discussion

In this paper we have calculated global $`m=1`$ modes with low pattern speed corresponding to introducing a finite disc eccentricity. An important aspect was the inclusion of self-gravity which is important for the structure of these modes as well as determining the motion of embedded protoplanets.

We considered disc models that were isolated or contained one or two protoplanets orbiting in an inner cavity. In all cases global modes were found that could be global on scales up to one hundred times the inner cavity radius. The modes could be considered as being of a type that were strongly coupled to the inner protoplanets or essentially free disc modes. In the former case the disc eccentricity could be comparable to that of the protoplanets for up to three times the outer protoplanet orbital semi-major axis with apsidal line antialigned with that of the protoplanet orbits. In the latter case the inner protoplanet orbital eccentricities were small compared to that found in the disc.

We went on to discuss the motion of a protoplanet embedded in an eccentric disc and determined, initially neglecting tidal torques, the equilibrium (non precessing) orbits which maintain apsidal alignment with the disc gas orbits. Equilibrium eccentricities were found to be comparable or possibly even exceed the disc eccentricity. In some cases secular resonance could occur producing particularly large protoplanet eccentricities.

We then formulated the calculation of the response of an eccentric disc to a protoplanet in the earth mass range in order to determine the time rate of change of the eccentricity and orbital migration rate. We found that equilibrium aligned orbits with very similar eccentricity to that of the gas disc may suffer no eccentricity change while undergoing inward migration in general. This was found from the resonant torque calculations but is also expected from direct consideration of the equations governing the tidal response. However, when the non precessing equilibrium aligned orbit has a significantly higher eccentricity than the disc, as can occur generally, but in particular for modes with very small pattern speed, orbital migration may be significantly reduced or reverse from inwards to outwards for the disc models we considered.

Attainment of high eccentricities in this way typically requires the characteristic test particle orbit precession frequency or mode pattern speed to significantly exceed the characteristic orbital circularization rate, a situation more likely for lower mass protoplanets. When tidal circularization dominates, the protoplanet equilibrium eccentricity is reduced while the apsidal line becomes significantly inclined to that of the disc. However, high protoplanet eccentricities could be excited by gravitational interactions between them (Papaloizou & Larwood 2000) and under favourable conditions this effect could act to counter tidal circularization generating significantly higher protoplanet eccentricities than that of the disc. This will be a topic for future investigation.

Although there is some uncertainty in the resonant torque calculations because of the need to sum contributions of different sign, weakening of the tidal interaction is expected on general physical grounds on account of larger protoplanet disc relative velocities at higher eccentricity. This indication of migration reversal at the higher eccentricities was found to be supported by a local dynamical friction calculation applicable in that limit. In this case the interaction near apocentre tends to speed up the protoplanet while the interaction near pericentre tends to slow it down. These effects are of opposite sign but the longer time spent near apocentre results in a net outward migration of the protoplanet for the surface density considered.

Thus the existence of global non circular motions in discs with radial excursions comparable to or exceeding the

semi-thickness may have important consequences for the migration of cores in the earth mass range. While processes of the type considered in this paper are unlikely to lead to the very high eccentricities observed for some giant planets, they may be important in controlling migration during planet formation as well as producing modest eccentricities $`0.2.`$

###### Acknowledgements.

The author thanks the IAP for visitor support and caroline Terquem for valuable and stimulating discussions as well as a carefull reading of a preliminary draft of this paper.

APPENDIX 1

## The time averaged potential due to a perturbing inner planet

The gravitational potential per unit mass due to a planet of mass $`m_p`$ located at $`𝐫_𝐩(r_p,\phi _p)`$ at $`𝐫(r,\phi )`$ is

$$\mathrm{\Phi }=\frac{Gm_p}{\sqrt{r^2+r_p^22rr_p\mathrm{cos}(\theta )}},$$
(64)
with $`\theta =\phi \phi _p.`$

In order to incorporate the effects of protoplanets orbiting interior to the disc we adopt a Jacobi coordinate system. In this system the coordinates of the innermost protoplanet are referred to the central star. The coordinates of the remainder are referred to the centre of gravity of the central mass and all interior protoplanets. The disc is referred to the centre of mass of central star and all inner protoplanets.

This has the following consequences:

For an object interior to the protoplanet with $`r<r_p,`$ the acceleration of the coordinate system due to the protoplanet must be allowed for. This gives rise, correct to first order in $`m_p,`$ to the indirect potential ( eg. Brouwer & Clemence 1961)

$$\mathrm{\Phi }_i=\frac{Gm_pr\mathrm{cos}(\theta )}{r_p^2},$$
(65)
to be added to the potential $`\mathrm{\Phi }.`$

For an object exterior to the protoplanet one must take account of the fact that the coordinate system is now based on the centre of mass of the inner protoplanets and central star. Assuming initially that $`m_p`$ is the only such protoplanet, the central potential is modified to become

$$\mathrm{\Phi }=\frac{GM_{}}{|𝐫+m_p𝐫_𝐩/M_{}|}.$$
(66)
To first order in $`m_p`$ this gives

$$\mathrm{\Phi }=\frac{GM_{}}{r}+\frac{Gm_pr_p\mathrm{cos}\theta }{r^2}.$$
(67)
In this case too the additional potential on the right hand side of the above can be incorporated into the perturbing planet potential and viewed as giving rise to an indirect potential. Thus the form of the perturbing potential

$$\mathrm{\Phi }=\frac{Gm_p}{\sqrt{r^2+r_p^22rr_p\mathrm{cos}(\theta )}}+\frac{Gm_prr_p\mathrm{cos}\theta }{max(r^3,r_p^3)},$$
(68)
incorporates both cases $`r<r_p`$ amd $`r_p<r.`$

In addition, although we included just one inner perturbing protoplanet, because we work to first order in their masses, the principle of linear superposition is valid such that the contributions of many such objects may be linearly superposed.

The perturbing potential due to a single protoplanet may be decomposed as a Fourier expansion in $`\theta .`$

Thus

$$\frac{\mathrm{\Phi }}{Gm_p}=\underset{m=0}{\overset{\mathrm{}}{}}\frac{K_m(r,r_p)}{\pi (1+\delta _{m0})}\mathrm{cos}m\theta ,$$
(69)
with

$$K_m(r,r_p)=\frac{1}{Gm_p}_0^{2\pi }\mathrm{\Phi }\mathrm{cos}m\theta d\theta .$$
(70)
For the problem on hand, namely the study of global $`m=1`$ modes of the disc planet system, we need only consider $`m=0`$ and $`m=1.`$ Then from the above

$$\frac{\pi \mathrm{\Phi }}{Gm_p}=\frac{1}{2}K_0(r,r_p)K_1(r,r_p)\mathrm{cos}\theta ,$$
(71)
with

$$K_0(r,r_p)=_0^{2\pi }\frac{1}{\sqrt{r^2+r_p^22rr_p\mathrm{cos}(\theta )}}𝑑\theta $$
(72)
and

$$K_1(r,r_p)=(_0^{2\pi }\frac{1}{\sqrt{r^2+r_p^22rr_p\mathrm{cos}(\theta )}}\mathrm{cos}\theta d\theta $$

$$\frac{\pi rr_p}{max(r^3,r_p^3)}).$$
(73)
## Time averaged potential for a protoplanet with small eccentricity

We now suppose the protoplanet has a small eccentricity, $`e_p`$ and write for its motion $`r_p=a_p(1+e_p\mathrm{cos}\omega t),\phi _p=\omega t2e_p\mathrm{sin}\omega t,`$ where without loss of generality we assume the apsidal line to be along $`\phi =0`$ where the displacement $`\xi _r(r_p)=e_pa_p`$ at $`t=0.`$

Here the protoplanet semi-major axis is $`a_p`$ and the orbital frequency is $`\omega .`$

Expanding to first order in $`e_p,`$ we find for the single protoplanet perturbing potential

$$\frac{\pi \mathrm{\Phi }}{Gm_p}=\frac{1}{2}K_0(r,a_p)\frac{a_pe_p}{2}\frac{K_0(r,a_p)}{a_p}\mathrm{cos}\omega t$$

$$\left(K_1(r,a_p)+a_pe_p\frac{K_1(r,a_p)}{a_p}\mathrm{cos}\omega t\right)\mathrm{cos}(\phi \omega t+2e_p\mathrm{sin}\omega t).$$

Performing the time average then gives

$$\frac{\mathrm{\Phi }}{Gm_p}=\frac{1}{2\pi }K_0(r,a_p)\frac{e_p}{2a_p\pi }\frac{(a_p^2K_1(r,a_p))}{a_p}\mathrm{cos}\phi .$$
(74)
For small $`e_p`$ we may replace $`a_p`$ by $`r_p`$ and expressing the result in terms of the radial displacement $`\xi _r(r_p),`$ we find

$$\frac{\mathrm{\Phi }}{Gm_p}=\frac{1}{2\pi }K_0(r,r_p)\frac{\xi _r(r_p)}{2r_p^2\pi }\frac{(r_p^2K_1(r,r_p))}{r_p}\mathrm{cos}\phi .$$
(75)
From equation (75) the time averaged potential due to a protoplanet in circular orbit $`(\xi _r(r_p)=0)`$ is

$$\mathrm{\Phi }_{ext}=\frac{Gm_p}{2\pi }K_0(r,r_p).$$
(76)
When $`\xi _r(r_p)0,`$ equation (75) gives the $`m=1`$ component of the perturbing potential as the real part of $`\mathrm{\Phi }_p^{}(r)\mathrm{exp}(i\phi ),`$ with

$$\mathrm{\Phi }_p^{}(r)=\frac{Gm_p\xi _r(r_p)}{2\pi r_p^2}\frac{(r_p^2K_1(r,r_p))}{r_p}.$$
(77)
APPENDIX 2

## The equations of motion in two dimensional orthogonal cordinates

The basic equations of motion (1) written in vector form are

$$\frac{𝐯}{t}+𝐯𝐯=\frac{1}{\mathrm{\Sigma }}\mathrm{\Pi }\mathrm{\Phi }.$$
(78)
We wish to write these in component form in a general two dimensional orthogonal coordinate system $`(a,\lambda )`$ in the Cartesian $`(x,y)`$ plane as introduced in section (5). The unit vectors in these orthogonal coordinate directions are $`𝐢_a=a/|a|`$ and $`𝐢_\lambda =\lambda /|\lambda |`$ respectively. The orthogonal vertical coordinate $`z`$ has an additional associated orthogonal unit vector $`\widehat{𝐤}=𝐢_a\times 𝐢_\lambda .`$ However, there is no dependence on $`z`$ and the velocity component in that direction can be taken to be zero.

Thus we may write $`𝐯(v_a,v_\lambda )`$ or

$$𝐯=v_a𝐢_a+v_\lambda 𝐢_\lambda .$$
(79)
To deal with equation (78), we first use the vector identity

$$𝐯𝐯=\frac{1}{2}(|𝐯|^2)+𝝎\times 𝐯,$$
(80)
where $`𝝎=\times 𝐯`$ is the vorticity (Arfken & Weber 2000).

Because we only wish to consider the equations in a two dimensional limit , only the $`z`$ component of vorticity, $`\omega _z`$ is non zero so that we may write equation (78) as

$$\frac{𝐯}{t}+\omega _z\widehat{𝐤}\times 𝐯+\frac{1}{2}(|𝐯|^2)=\frac{1}{\mathrm{\Sigma }}\mathrm{\Pi }\mathrm{\Phi }.$$
(81)
We may now write equation (81 ) in component form using the identities (see Arfken & Weber 2000)

$$a\frac{}{a}+\lambda \frac{}{\lambda }$$
(82)
and

$$\omega _z=\left(\times 𝐯\right)\widehat{𝐤}=|\lambda ||a|\left(\frac{(\frac{v_\lambda }{|\lambda |})}{a}\frac{(\frac{v_a}{|a|})}{\lambda }\right).$$
(83)
Doing this we obtain

$$\frac{v_a}{t}\omega _zv_\lambda +\frac{1}{2}|a|\frac{(v_a^2+v_\lambda ^2)}{a}=\frac{1}{\mathrm{\Sigma }}|a|\frac{\mathrm{\Pi }}{a}|a|\frac{\mathrm{\Phi }}{a}.$$
(84)
and

$$\frac{v_\lambda }{t}+\omega _zv_a+\frac{1}{2}|\lambda |\frac{(v_a^2+v_\lambda ^2)}{\lambda }=\frac{1}{\mathrm{\Sigma }}|\lambda |\frac{\mathrm{\Pi }}{\lambda }|\lambda |\frac{\mathrm{\Phi }}{\lambda }.$$
(85)
After inserting (83) into (84) and (85) it is a simple matter to obtain equations (34) and (35) as given in section (5).