A&A 485, 877-895 (2008)
DOI: 10.1051/0004-6361:20078702

On disc protoplanet interactions in a non-barotropic disc with thermal diffusion

S.-J. Paardekooper - J. C. B. Papaloizou

DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK

Received 18 September 2007 / Accepted 25 April 2008

Abstract
Aims. We study the disc planet interactions of low-mass protoplanets embedded in a circumstellar disc. We extend the standard theory of planet migration from the usual locally isothermal assumption to include non-barotropic effects, focusing on the validity of linear theory.
Methods. We compared solutions of the linear equations with results from non-linear hydrodynamic simulations, where in both cases we adopted a background entropy gradient and solved the energy equation.
Results. We show that the migration behaviour of embedded planets depends critically on the background radial entropy gradient in the disc. The presence of such a gradient not only changes the corotation torque on the planet, but also always guarantees a departure from linear behaviour, which gives a singular density response at corotation, in the absence of thermal or viscous diffusion. A negative entropy gradient tends to give rise to positive, non-linear corotation torques apparently produced as material executes horseshoe turns at approximately constant entropy. These torques have no counterpart in linear theory, but can be strong enough to push the planet outwards until saturation starts to occur after a horseshoe libration period. Increased thermal diffusion acts to reduce these non-linear torques, but, at the same time, it can help to prevent their saturation. In combination with a small kinematic viscosity that is able to maintain a smooth density profile the positive torque could be sustained.

Key words: hydrodynamics - planets and satellites: formation

1 Introduction

During their formation process inside circumstellar discs, planets can change their orbital parameters by gravitational interaction with the gaseous disc (Goldreich & Tremaine 1980). Torques generated at various resonances can promote or damp eccentricity growth (Goldreich & Sari 2003), and change the semi-major axis of the forming planet (Ward 1997). One can distinguish torques due to waves, generated at Lindblad resonances, which propagate away from the planet, and corotation torques, generated near the orbit of the planet, where material on average corotates with the planet (see Goldreich & Tremaine 1979). The usual result of these torques is that the planet migrates inwards, towards the central star, at a rate that is determined by the planet mass and the disc parameters (see Ward 1997).

In terms of planet and disc mass, three types of migration can be distinguished. Low-mass planets, whose gravitational influence is not strong enough to overcome pressure effects, generate a linear disc response that gives rise to an inward migration rate that is proportional to the planets mass. This is called type I migration (Ward 1997). For standard Solar nebula parameters, type I migration applies to planets up to several Earth masses ( ${{M}_\oplus }$).

Higher-mass planets excite non-linear waves, and are able to tidally truncate the disc, opening an annular gap around their orbit (Lin & Papaloizou 1979,1986). Migration occurring when such a gap is maintained is referred to as type II migration (Ward 1997), and it takes the planet inward on a viscous time scale. The minimum mass for opening a gap, again for standard disc parameters, is approximately one Jupiter mass.

Both types of migration have been studied extensively through numerical hydrodynamical simulations, in a two-dimensional set-up (Nelson et al. 2000; D'Angelo et al. 2002) as well as in three dimensions (Bate et al. 2003; D'Angelo et al. 2003b), giving good agreement in the respective mass regimes. In the intermediate mass regime, depending on the disc viscosity and surface density distribution, the corotation torque may get a non-linear boost to the extent that it determines the sign of the total torque (Masset et al. 2006).

Corotation torques are also responsible for a third type of migration, sometimes called runaway migration (Masset & Papaloizou 2003). This type III migration (Artymowicz 2004) can be very fast, and applies to planets that open up a partial gap and are embedded in a massive disc (Masset & Papaloizou 2003).

In all studies referred to above, a simplified equation of state, for which the pressure depends on the density and the local radius only was used. The dependence on radius comes about from adopting a fixed radial temperature profile. As at any location, the pressure depends only on the density, we describe such an equation of state as being locally barotropic. The use of such an equation of state removes the need to solve the energy equation and thus makes the equations more tractable to tackle both analytically and numerically. The question whether this approach is valid has only been addressed fairly recently through numerical simulations that do include the energy equation (Klahr & Kley 2006; D'Angelo et al. 2003a), focusing on high-mass planets only. Morohoshi & Tanaka (2003) studied the effect of optically thin cooling on disc-planet interactions using a local approach, while Jang-Condell & Sasselov (2005) calculated the torque on low-mass planets analytically in discs with a realistic temperature profile.

Paardekooper & Mellema (2006a) showed that relaxing the barotropic assumption can change the migration behaviour of low-mass planets dramatically. There it was shown, using three-dimensional radiation-hydrodynamical simulations, that, for the disc parameters adopted, low-mass planets migrate outwards through the action of a strong positive corotation torque. Due to the excessive computational overheads associated with such simulations, it was not feasible to do a parameter survey for different planet masses and disc models.

In this paper, we aim to clarify the origin of the positive corotation torque, its dependence on planet mass and disc parameters together with the relationship of this type of disc planet interaction to the linear perturbation theory that has been successfully applied to the type I migration regime.

To ease the computational demands we focus on non-barotropic two-dimensional discs (commonly termed cylindrical discs, see Hawley 2000) which incorporate thermal diffusion rather than the more complex radiative transfer.

The plan of the paper is as follows. In Sect. 2 we review the governing equations and the disc models used. Section 3 is devoted to the linear theory of the interaction of a low-mass planet with a gaseous disc, and we solve the resulting equations in Sect. 4. In order to remove singularities we adopted the well known Landau prescription together with a small thermal diffusivity applied in some cases. Although migration torques could be calculated in the limit of zero dissipative effects, the density response became singular implying that non-linear effects may always be important in the corotation region.

After describing the numerical set-up in Sect. 5, we perform fully non-linear hydrodynamical simulations in Sect. 6. Indeed we find that a non linear small scale coorbital density structure is set up at early times. We go on to investigate the effect of thermal diffusivity, finding that a non zero value is required to ensure the coorbital region is adequately resolved. The associated torque is then indeed found to be positive, for favourable density gradients, and sufficiently strong and negative entropy gradients, for about one horseshoe libration period after which saturation sets in. In Sect. 7 we present a simple non-linear model of the corotation region based on material executing horseshoe turns with constant entropy, which can result in a density structure that leads to a temporary positive torque on the planet when there is a background negative entropy gradient and compare it with simulations.

Whether the positive torque can be sustained for longer times will depend on the evolution of the planetary orbit and disc and whether these circumstances provide an appropriate corotational flow that can resupply appropriate low entropy material. To investigate this aspect, we study the long-term evolution of the torque in Sect. 8, and in the case of a particular example, show that for an appropriate rate of thermal diffusion and a small viscously driven mass flow through corotation, the entropy-related corotation torque can be sustained for several libration periods.

We go on to give a discussion of our findings in Sect. 9 and conclude in Sect. 10.

   
2 Basic equations and disc models

The basic equations are those of the conservation of mass, momentum and energy for a two dimensional cylindrical disc in a frame rotating with angular velocity  $\Omega_{\rm p}$ in the form

$$% {\partial \rho \over \partial t}= -\nabla\cdot(\rho {\vec v})$$ (1)

$$% {{\rm D} {\vec v} \over {\rm D} t} +2\Omega_{\rm p}{\hat {\vec k}} \times {\vec v} = -{1\over \rho}\nabla P- \nabla\Phi$$ (2)

$$% \rho{ {\rm D} E \over {\rm D}t} - {P\over \rho} { {\rm D} \... ...equiv \rho T {{\rm D} S \over {\rm D}t}= -\nabla\cdot{\vec F}.$$ (3)

Here, $\rho$ denotes the density, ${\vec v} = (v_r, v_{\varphi})$ the velocity, ${\hat {\vec k}}$ denotes the unit vector in the vertical direction, P the pressure and $\Phi =\Phi_{\rm G} +\Omega_{\rm p}^2r^2/2$ where $\Phi_{\rm G}$ is the gravitational potential. The convective derivative is defined by

$$% {{\rm D}\over {\rm D}t} \equiv {\partial \over \partial t}+ {\vec v}\cdot \nabla,$$ (4)

and ${\vec F} = -K \nabla T$ is the heat flux, with K being the thermal conductivity and T being the temperature (of course radiation transport may be incorporated in this formalism). In our numerical simulations, we choose K=K(r) such that the initial temperature profile gives $\nabla\cdot{\vec F} =0$. We adopt an ideal gas equation of state such that

$$% P=R_{\rm g}\rho T,$$ (5)

with $R_{\rm g}$ being the gas constant. The internal energy per unit mass is given by

$$% E =\frac{P}{(\gamma-1)\rho},$$ (6)

where $\gamma$ is the constant adiabatic exponent and S = $R_{\rm g}\log (P/\rho^{\gamma})/(\gamma-1)$ is the entropy per unit mass.

We adopt a cylindrical polar coordinate system $(r,\varphi)$ with origin (r=0) located at the central mass. The self-gravity of the disc is neglected. Thus the gravitational potential is assumed to be due to the central mass and perturbing planet such that $\Phi_{\rm G} = \Phi_{\rm G0}+ \Phi_{\rm Gp},$ where

$$% \Phi_{\rm G0} = -{GM_*\over r},$$ (7)

and

 

$$% \Phi_{\rm Gp} = {-GM_{\rm p}\over\sqrt{r^2+r_{\rm p}^2-2rr_{\rm... ..._{\rm p}^2}} +{GM_{\rm p}r \cos(\varphi-\varphi_{\rm p})\over r_{\rm p}^2}\cdot$$ (8)

In the above the last term is the indirect term, M_* denotes the mass of the central object, with $M_{\rm p},$ $r_{\rm p}$, and $\varphi_{\rm p}$ denoting the mass, radial and angular coordinate of the protoplanet respectively. We also allow for a gravitational softening parameter b.

We remark that there are several limiting cases in the above description. When heat transport is neglected, (K=0), the system is adiabatic. When the energy equation is dropped and the temperature is taken to be a fixed function of r we have the usual locally isothermal limit. When $\Omega_{\rm p}=0$ the reference frame is non-rotating but non-inertial as the origin accelerates together with the central mass. This is accounted for by the indirect term in the potential. Numerical calculations are most conveniently performed in a frame corotating with the protoplanet. Then $\Omega_{\rm p}$ becomes the circular Keplerian angular velocity at radius $r_{\rm p}$. At a general radius, r, this is given by $\Omega_{\rm K}= (GM_*/r^3)^{1/2}$.

The formulation given above applies to a cylindrical disc model where there is no vertical stratification or dependence on the vertical coordinate. The equations apply to a cylindrical system with a constant vertical semi-thickness, H₀ which may be specified arbitrarily as there is no explicit dependence on it. The density, $\rho$, is taken to have the same spatial variation as the surface density to which it is related through $\Sigma =2\rho H_0.$For a thin disc of the type we wish to consider, at any location there is a putative local vertical semi-thickness $H = c_{\rm s}/\Omega_{\rm K}$, with $c_{\rm s}=\sqrt{P/\rho}$ being the local isothermal sound speed that is a function of position and associated with the neglected vertical stratification. When considering the physical state variables associated with any such location we regard H₀ as having been adjusted to coincide with H.

2.1 Equilibrium models

In this paper we work with equilibrium model discs into which a perturbing protoplanet is subsequently inserted. These are such that density has a power law dependence of the form $\rho \propto r^{-\alpha}$, with $\alpha$ being constant. We also adopt $T \propto r^{-1}$, giving a putative semi-thickness $H \propto r$, then we have $P \propto r^{-\alpha-1}$. The local gas angular velocity required for radial hydrostatic equilibrium is given by $\Omega = \Omega_{\rm K}(1-(\alpha +1)(H/r)^2)$. Thus for these models $\Omega \propto \Omega_{\rm K}$ and the epicyclic frequency $\kappa =\Omega$. Here we recall in general that $\kappa^2 = (2\Omega/r) {\rm d}(r^2\Omega)/{\rm d}r$. We also adopt the representative value $\gamma =1.4$ unless stated otherwise.

   
3 Linear theory

We here consider the response of the disc to the forcing of a protoplanet in a circular orbit assuming that the induced perturbations are small so that the basic equations may be linearised. We perform a Fourier decomposition of the perturbing potential such that in the non-rotating frame

$$% \Phi_{\rm Gp} = {\cal R}_{\rm e} \left(\sum_{m=0}^{\infty} ... ...m Gp}m}\exp{({\rm i}m\varphi-{\rm i}m\Omega_{\rm K}t)}\right),$$ (9)

where ${\cal R}_{\rm e}$ denotes that the real part is to be taken. Here, as well as in the rest of this section, $\Omega_{\rm K}$ is short for $\Omega_{\rm K}(r_{\rm p})$.

Similarly, the perturbation response of state variables that are non-zero in the equilibrium, indicated with a prime, for each value of m, has an exponential factor $\exp{({\rm i}m\varphi-{\rm i}m\Omega_{\rm K}t)}$ with an amplitude depending only on r. Equations for these amplitudes are obtained by linearizing the basic equations. We begin by considering the adiabatic limit in which the thermal diffusivity is set to zero. Then we linearize the adiabatic condition

$$% {{\rm D}S \over {\rm D}t} = {\partial S \over \partial t}+{\vec v}\cdot \nabla S =0,$$ (10)

to obtain

$$% S' = -{v_r\over {\rm i}{\overline \sigma}} {{\rm d} S \over {\rm d}r},$$ (11)

where ${\overline \sigma} = m(\Omega-\Omega_{\rm K})$. Expressed in terms of the density and pressure perturbations Eq. (11) takes the equivalent form

$$% {P'\over \gamma P} -{ \rho'\over \rho}= -{v_r{\cal A}\over {\rm i}{\overline \sigma}g_r},$$ (12)

where $g_r = -(1/\rho)({\rm d}P/{\rm d}r)$ and the square of the Brunt-Väisälä frequency is given by

$$% {\cal A} = -{1\over \rho} {{\rm d} P \over {\rm d}r} \left(... ... {\rm d}r} -{1\over \rho}{{\rm d} \rho \over {\rm d}r}\right),$$ (13)

which can be negative in regions of formal local convective instability. Equation (12) can be used to express the density perturbation in terms of the pressure and velocity perturbations.

Linearization of the two components of the equation of motion together with the continuity equation then yield, after elimination of the azimuthal velocity perturbation, a pair of first order ordinary differential equations for the quantities $Q= P'/P^{1/\gamma}$ and $U =r P^{1/\gamma} v_r$ which may be written in the form

$$% {{\rm d}Q\over {\rm d}r} = C_1v_r +C_2Q + S_1,$$ (14)

$$% {{\rm d}U \over {\rm d}r} = D_1Q + D_2v_r + S_2,$$ (15)

where the coefficients are given by

$$C_1 = -{{\rm i}\rho P^{-(1/\gamma)}} ({\overline \sigma} -(\kappa^2 + {\cal A})/{\overline \sigma}),$$

$$C_2 = -(2 m \Omega)/(r {\overline \sigma}),$$

$$D_1 = -{\rm i} r({\overline \sigma}^2 - \gamma P m^2/(r^2\rho))/ ({\overline \sigma}\gamma P^{1-2/\gamma}),$$

$$D_2 = (P^{1/\gamma} m \kappa^2)/(2 \Omega {\overline \sigma}),$$

$$S_1 = -{\rho P^{(-1/\gamma)}}[({\rm d} \Phi_{{\rm Gp}m}/{\rm d}r) + (2 m\Omega)\Phi_{{\rm Gp}m}/(r {\overline \sigma})],$$

$$S_2 = {\rm i}( P^{1/\gamma} m^2)/(r {\overline \sigma}) \Phi_{{\rm Gp}m}.$$

3.1 The corotation singularity

The perturbation of the protoplanet causes the excitation of outgoing density waves that are associated with a conserved wave action or angular momentum flux. This causes a torque to act on the protoplanet. However, angular momentum exchange between protoplanet and disc may also occur directly at corotation where ${\overline \sigma}=0$. In linear theory, this type of exchange is localised at corotation through the operation of a corotation singularity. In order to study the domain near corotation where ${\overline \sigma}=0$, we neglect ${\overline \sigma}^2$ in the first set of brackets in the expression for D₁ above. Then we may derive a single equation for U from Eqs. (14) and (15) in the form

 

$$% { P^{2/\gamma}\over \rho r} {{\rm d}\over {\rm d}r}\left( {\rho... ...m i}\Phi_{{\rm Gp}m}m^2 {\cal A} P^{1/\gamma}\over {\overline \sigma}rg_r}\cdot$$ (16)

We see that Eq. (16) is in general singular at corotation with a second order pole at ${\overline \sigma}=0$. This will result in angular momentum exchange with the perturber. In order to be singularity free we require both that the entropy gradient be zero or equivalently ${\cal A}=0$, and that ${\rm d}(\rho\kappa^2/ (2\Omega P^{2/\gamma})/{\rm d}r =0$. The latter condition is the generalisation of the condition that the gradient of specific vorticity, or vortensity, given by

$$% \frac{{\rm d}}{{\rm d}r}\left(\frac{\kappa^2}{2\Omega\rho}\right),$$ (17)

should be zero that applies in the strictly barotropic case.

3.2 The Landau prescription

In order to solve the forced response problem, we have to deal with the corotation singularity. Physically this is resolved by including the dissipative effects of heat transport and viscosity. However, by increasing the order of the system of equations, this significantly increases the complexity of the problem. For many purposes this can be avoided by using the Landau prescription. According to this the forcing frequency is given a small positive imaginary part such that ${\overline \sigma} \rightarrow m(\Omega-\Omega_{\rm K}-{\rm i}\epsilon\Omega_{\rm K}),$ where $\epsilon$ is very small. In practice, as is the case with the problem considered here, the torques between disc and protoplanet converge as $\epsilon \rightarrow 0$ even though the solution becomes singular.

The latter feature means that even though torques may be convergent the response becomes non-linear near corotation making the linear theory invalid. The time $1/(\Omega_{\rm K}\epsilon)$ may be interpreted as the time scale over which the perturbing potential is turned on. Thus a non-linear breakdown for $\epsilon$ below a certain value, $\epsilon_0$, means that, viewed from the point of view of an initial value problem, the linear solution should be considered as being valid only for a finite time period, which may be estimated as to be of order  $1/(\Omega_{\rm K}\epsilon_0)$.

For the problem on hand, non-linear breakdown occurs mainly through terms associated with the entropy gradient or ${\cal A}$, these being associated with a second order pole. This type of singularity may be removed by introducing heat diffusion so that we would expect that for a sufficiently high thermal diffusivity, the validity of linear theory should be restored. This we explore below.

  
3.3 The effect of a small thermal diffusivity

In order to investigate the effect of thermal diffusion on the corotation singularity, we now consider the modification of the condition given by Eq. (11) that occurs when a small thermal diffusivity is present.

In this case we must linearize the full energy equation, Eq. (3), which leads to

$$% {\rm i}{\overline \sigma}S' + v_r {{\rm d} S \over {\rm d}r}= {K\over \rho T}\nabla^2 T'.$$ (18)

Here, because of rapid variation near corotation, when considering terms involving K, we neglect everything other than the second derivative term. Assuming that the variation of the pressure perturbation can be neglected, (this approach can be validated by inspection of the form of the solutions), we then have $S' = (\partial S/\partial T)_{\rm p} T'$, and accordingly

$$% {\rm i}{\overline \sigma}S' + v_r {{\rm d} S \over {\rm d}r}= {K\over C_{\rm p} \rho }\nabla^2 S',$$ (19)

where $C_{\rm p} =\gamma R_{\rm g}/(\gamma-1)$ is the specific heat at constant pressure. As we are interested in a local region around corotation, we set $r=r_{\rm p}+x$, where $x \ll r_{\rm p}$, and ${\overline \sigma}\rightarrow -3m \Omega_{\rm K}x/(2r_{\rm p}).$ In addition we replace $\nabla^2 \rightarrow ({\rm d}^2/{\rm d}x^2 -m^2/r_{\rm p}^2)$.

Then Eq. (19) can be written in the form

$$% {{\rm d}^2 S'\over {\rm d}z^2} -{\rm i}\beta z S' = F,$$ (20)

where $z= x+2Dm{\rm i}/(3r_{\rm p}\Omega_{\rm K})$, $\beta = -3m\Omega_{\rm K}/(2r_{\rm p}D)$, $D=K/(\rho C_{\rm p})$, and $F = v_r ({\rm d} S / {\rm d}r)/D$, the last three quantities being assumed constant. Equation (20) can be solved in terms of in-homogeneous Airy functions (e.g. Abramowitz & Stegun 1970) in the form

$$% S'=-\left({3^{1/3}\pi F\over (-\beta)^{2/3}}\right) Hi(\zeta)$$ (21)

with

$$% \zeta = -{2({\rm i}{\overline \sigma}+\epsilon_1) (9m/2D_{\rm e})^{1/3}\over 3m \Omega_{\rm K}}$$ (22)

and

$$% \epsilon_1 = m^{2} D_{\rm e}\Omega_{\rm K},$$ (23)

with $D_{\rm e}$ being the dimensionless diffusivity $D_{\rm e}= D/(r_{\rm p}^{2} \Omega_{\rm K}).$ Here all quantities are evaluated at corotation, $r=r_{\rm p}$, apart from ${\overline \sigma} \rightarrow -3m \Omega_{\rm K}x/(2r_{\rm p})$.

The in-homogeneous Airy function $Hi(\zeta)$ (see Abramowitz & Stegun 1970) is defined by

$$% Hi(\zeta) ={1\over \pi} \int_0^{\infty} \exp{(-k^3 +k\zeta)} {\rm d}k.$$ (24)

Using this solution, S' is no longer singular at corotation. By comparing Eq. (11) for S' with the solution determined by Eq. (21) we can see how to remove the corotation singularity.

It is apparent that, where it appears multiplied by the entropy gradient, the singular quantity

$$% {\Omega\over \overline \sigma} \rightarrow {2\pi {\rm i}\ov... ...n_1)(9m/2D_{\rm e}) ^{1/3}\over 3m \Omega_{\rm K}}\right)\cdot$$ (25)

In this way the divergence at ${\overline \sigma}=0$ associated with terms that involve the entropy gradient is removed. When we implemented cases with thermal diffusivity as described below, wherever $\overline \sigma$ otherwise occurred, we applied the Landau prescription with $\epsilon$ = 10^-5.

From Eq. (20), the length scale associated with the diffusively controlled corotation region is $\vert\beta\vert^{-1/3} \sim r_{\rm p}(D_{\rm e}/m)^{1/3}$. Thus for sufficiently large $D_{\rm e}$, we expect that the effects of corotation associated with a radial entropy gradient should be reduced such that the validity of linear theory is restored as far as these are concerned. Additional restrictions may result from a radial vortensity gradient.

   
4 Linear response calculations

We have solved Eqs. (14) and (15) for a variety of disc models, softening parameters and values of m. We have used both the Landau prescription and the combination of that, together with a finite thermal diffusivity used to deal with terms associated with the entropy gradient, to deal with the corotation singularity. The equations were integrated using a fifth order Runge-Kutta scheme and outgoing wave conditions were applied. In practice these were applied at radii that became closer to the protoplanet as m increased.

For convenience throughout, we adopt a system of units for which $r_{\rm p}=1$, $\Omega_{\rm K}(r_{\rm p}) =1$, and the density of the unperturbed disc at the protoplanet orbital radius $\rho(r_{\rm p})=1$. The orbital period of the protoplanet is then $2\pi$. In all cases the aspect ratio  $c_{\rm s}/(r\Omega)$ was taken to be 0.05.

Figure 1: The functions U and Q divided by the protoplanet to central star mass ratio are plotted for $\rho _0 \propto r^{-1/2},$ b=0.03 and m=2 in dimensionless units. Real parts are plotted with dashed lines, imaginary parts with full lines. Curves for the Landau parameters $\epsilon = 10^{-5}$ $\epsilon = 3.\times 10^{-2}$ are shown. The response is very similar for these cases but the latter case is slightly more damped which is manifest by somewhat smaller amplitude swings.

Open with DEXTER

Figure 2: The density response divided by the protoplanet to central star mass ratio, for $\rho _0 \propto r^{1/2},$ b=0.03 and m=2 in dimensionless units. Real part dashed line, imaginary part full line. For the left panel the Landau parameter $\epsilon = 10^{-5}$ and for the right panel $\epsilon = 6$ $\times$ 10-2. In the former case the amplitudes of the real and imaginary parts reach extreme values of order 105.

Open with DEXTER

Figure 3: The cumulative torque in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$ acting on the protoplanet plotted as a function of dimensionless radius for $\rho _0 \propto r^{-1/2}$, b=0.03 and m=2. For the left panel plots are given for the Landau parameters $\epsilon = 10^{-5}$ (dashed line) $\epsilon = 10^{-4}$ (dotted line) $\epsilon = 3$ $\times$ 10-2 (full line). For the right panel plots are given for the dimensionless thermal diffusivities $D_{\rm e}= 10^{-6}$ (dotted line) and $D_{\rm e}= 10^{-5}$ (full line). The cumulative torque is defined to be zero at the outer boundary and equal the total torque at the inner boundary. In all of these cases the total torque is positive.

Open with DEXTER

4.1 Results

The functions U and Q are plotted for the disc with $\rho _0 \propto r^{-1/2}$, b=0.03, being comparable to H/r, and m=2 in Fig. 1. In these cases a Landau prescription was used to deal with the corotation singularity. We show curves for the Landau parameters $\epsilon = 10^{-5}$ and $\epsilon = 3$ $\times$ 10^-2. As is indicated in Fig. 1, the response functions are very similar in these cases for which the Landau parameter varies by a factor of 3000. The $\epsilon = 3.\times 10^{-2}$ case not unexpectedly shows slightly more damping.

However, the density response exhibits different behaviour. This is plotted for the Landau parameters $\epsilon = 10^{-5}$ and $\epsilon = 6$ $\times$ 10^-2 in Fig. 2, for a disc with $\alpha =-1/2$ (very similar behaviour occurs for $\alpha =1/2$). This response shows a strong singularity as $\epsilon$ is decreased to zero. The extreme values reached are $\propto$ $1/\epsilon$. As the singularity is approached the radial width $\propto$$\epsilon$. The effect of this is that although the pressure and velocity response as well as the migration torque converge as $\epsilon \rightarrow 0$, the density response becomes increasingly singular. This means that linear theory always breaks down in that limit. This breakdown is in fact associated with terms $\propto$${\cal A},$ the entropy gradient and it can be removed by considering thermal diffusion. To deal with this, one can use the relative smoothness of the pressure and velocity perturbations to obtain a simple description governed by a second order ordinary differential equation (see above). In order to remain in a regime for which linear theory is valid as far as these are concerned it is essential to either incorporate a non zero thermal diffusivity or restrict consideration to values of $\epsilon$ that are high enough for linear theory to be valid. The latter viewpoint can be regarded as restricting consideration of evolution times for an initial value problem to values smaller than $\sim$ $1/(\epsilon\Omega_{\rm K}(r_{\rm p}))$.

4.2 Migration torque

The torque acting on the protoplanet is calculated directly by summing the torque integrals evaluated for each Fourier component in the form

$$% {\cal T}(r_{\rm in}) = 2H{\cal I}_m\left[\int^{r_{\rm out}}_{r_{\rm in}}\pi m r \rho' \Phi_{{\rm Gp}m}{\rm d}r\right],$$ (26)

where the integral is taken over the unperturbed disc domain $(r_{\rm in},r_{\rm out})$ and ${\cal I}_m$ denotes that the imaginary part is to be taken. The factor 2H takes account of the vertical direction by assuming the cylindrical disc model to apply over a vertical extent 2H. We use the above to define a cumulative torque as a function of r for each m. Thus

$$% {\cal T}(r) = \int^{r_{\rm out}}_{r} T_{\rm p}(r){\rm d}r,$$ (27)

where $T_{\rm p}(r)$ is the torque density. Then the total torque acting on the protoplanet is ${\cal T}(r_{\rm in})$.

The cumulative torque acting on the protoplanet is plotted as a function of dimensionless radius for $\rho _0 \propto r^{-1/2}$, b=0.03 and m=2 in Fig. 3. Plots are given for cases using the Landau prescription with parameters $\epsilon = 10^{-5}$, $\epsilon = 10^{-4}$ and $\epsilon = 3$ $\times$ 10^-2. For comparison plots are given for the dimensionless thermal diffusivities $D_{\rm e}= 10^{-6}$ (dotted line) and $D_{\rm e}= 10^{-5}$ implemented as described above. In all of these cases the net torque is positive. We remark that the curves for $\epsilon = 10^{-5}$ and $\epsilon = 10^{-4}$ are very close as are those for $D_{\rm e}= 10^{-6}$ and $D_{\rm e}= 10^{-5}$. This indicates good convergence of the torques as dissipative effects are reduced.

To emphasise this feature we show the cumulative torque acting on the protoplanet as a function of m for $\rho _0 \propto r^{-1/2}$ for the two softening parameters b=0.03 and b=0.01 in Fig. 4. The softening parameter that should be used is uncertain but it should be of the order of the aspect ratio to simulate 3D effects. Plots are given for the Landau parameters $\epsilon = 10^{-5}$, $\epsilon = 10^{-4}$, $\epsilon = 10^{-3}$ and $\epsilon = 3$ $\times$ 10^-2. The total torques are always negative corresponding to inward migration. Again the convergence for small $\epsilon$ is good but note that torques of larger magnitude are obtained for the smaller softening parameter and larger m need to be considered.

Similar results are obtained for different disc models and when thermal diffusion is employed as illustrated in Fig. 5. There we plot cumulative torques as a function of m for $\rho _0 \propto r^{-3/2}$ and b=0.03. We show also the corresponding plots for $\rho _0 \propto r^{-1/2}$, b=0.03 but with the implementation of thermal diffusion.

Figure 4: The cumulative torque in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$ acting on the protoplanet as a function of m for $\rho _0 \propto r^{-1/2},$ b=0.03 ( left panel) and b=0.01 ( right panel). For both panels plots are given for the Landau parameters $\epsilon = 10^{-5}$ (full line), $\epsilon = 10^{-4}$ (dotted line), $\epsilon = 10^{-3}$ (dashed line), $\epsilon = 3$ $\times$ 10-2 (dot dashed line).

Open with DEXTER

Figure 5: The cumulative torque in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$ acting on the protoplanet as a function of m for $\rho _0 \propto r^{-3/2}$ and b=0.03 ( left panel) and for $\rho _0 \propto r^{-1/2}$, b=0.03 but with thermal diffusion ( right panel). For the left panel plots are given for the Landau parameters $\epsilon = 10^{-5}$ (full line), $\epsilon = 10^{-4}$ (dotted line), $\epsilon = 10^{-3}$ (dashed line), $\epsilon = 3$ $\times$ 10-2 (dot dashed line). In the right panel plots are given for the dimensionless thermal diffusivities $D_{\rm e}= 10^{-6}$ (dotted line) and $D_{\rm e}= 10^{-5}$ (full line). The curves almost coincide for these cases.

Open with DEXTER

As indicated above, linear theory always breaks down for sufficiently small dissipation parameters, and it is important to estimate parameter regimes where it could be valid. We have done this both when a pure Landau prescription is used and also when thermal diffusion is applied. To do this we calculate the full response by summing over m and determine the condition that the perturbation to the entropy gradient be less than or equal to the same magnitude as the unperturbed value. Clearly non-linear effects may set in for weaker perturbations than those that result in this condition being marginal. Thus we estimate that non-linearity sets in for smaller diffusion coefficients than those obtained when this condition is marginally satisfied. We show a contour plot of the full linear density response for $\rho _0 \propto r^{-1/2}$ and thermal diffusivity $D_{\rm e}= 10^{-6}$ in Fig. 6. From this and other similar responses we obtain the condition

$$% \left({0.03\over b}\right)^{1/3} {(1.55\times (q/10^{-5})) \over {\left(D_{\rm e}/10^{-6}\right)^{2/3}}} < 1.$$ (28)

Similarly for a pure Landau prescription, we obtain

$$% \left({0.03\over b}\right)^{1/4} {2.5\times 10^{7}q \over \left(\epsilon/10^{-3}\right)^{2}} < 1.$$ (29)

Interestingly, from Eq. (28) we estimate that for protoplanets in the earth mass range, if $D_{\rm e}< 10^{-6}$, there should be departures from linear theory. Similarly from Eq. (29), we estimate protoplanets in this mass range to be in the non-linear regime if $\epsilon <$ $\sim$10^-2. This indicates that linear theory may only be relevant for short evolutionary times in cases with very low dissipation and then from Fig. 4 we would expect that the full linear torque is never established. Although the above estimates are uncertain, we comment that the scaling implied by Eq. (29), that for the same degree of non-linearity $\epsilon \propto q^{1/2}$, implies that the time required to attain the same degree of departure from linear theory should scale as q^-1/2.

   
5 Numerical hydrodynamical simulations

In this section, we describe the set-up for our numerical hydrodynamical simulations of two-dimensional gas discs with embedded planets, using a non-barotropic equation of state. The initial disc models and system of units are the same as those employed for the linear calculations. The softening parameter is b =0.03. This approximately corresponds to the putative disc semi-thickness and the use of such a softening parameter approximates the result of appropriate vertical averaging of the gravitational potential. There is no explicit viscosity in the model unless otherwise stated. In Sect. 6.3 we include explicit heat diffusion. For locally isothermal simulations, we drop the energy equation and instead use a fixed temperature profile that gives rise to a constant aspect ratio in the initial state. For models including the energy equation we have used different temperature profiles to vary the entropy gradient while keeping the density gradient constant. Unless stated otherwise, we slowly build up the mass of the planet during the first three orbits to avoid to avoid transients due to the sudden introduction of the planet into the disc. For low-mass planets in isothermal simulations this is usually unnecessary, but the non-linearities in the flow associated with an entropy gradient can introduce some artifacts (see below).

Figure 6: Contour plot of the full linear density response for $\rho _0 \propto r^{-1/2}$ and thermal diffusivity $D_{\rm e}= 10^{-6}$ in dimensionless units. The linear theory is only likely to be valid for planet masses significantly less than an earth mass in this case (see Eq. (28)). Note the narrow high and low density regions on the corotation circle near the protoplanet. These, also seen in the non-linear simulations, respectively lead and trail the protolanet giving rise to a positive corotation torque. However, the effects of this torque are not strong enough to counteract the negative Lindblad torques contributed by the high density wakes.

Open with DEXTER

	Figure 7: Total torque in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$ as a function of time for two different planet masses and two density profiles for the disc. The disc was kept locally isothermal, and a viscosity $\nu =10^{-5}$ $r_{\rm p}^2 \Omega _{\rm p}$. The dotted lines indicate the results of our linear calculations.
Open with DEXTER

5.1 Code description

We use the RODEO method (Paardekooper & Mellema 2006b) to evolve the two-dimensional Euler equations in cylindrical coordinates  $(r,\varphi)$. The grid consisted of 512 radial cells, equally spaced between r=0.4 and r=1.8, and 2048 azimuthal cells, equally spaced over $2\pi$. We used non-reflective boundary conditions (for details see Paardekooper & Mellema 2006b). The energy equation is a straightforward extension of the method (see Eulderink & Mellema 1995), which was also used in Paardekooper & Mellema (2006a) to perform three-dimensional, radiation-hydrodynamical simulations of embedded planets. In this paper, we have not used the radiation diffusion solver. Heat diffusion is incorporated the same way as Navier-Stokes viscosity in Paardekooper & Mellema (2006b).

In calculating the torque on the planet, we include all disc material, including that within the Hill sphere of the planet. For low-mass planets considered here, this should not influence the results.

5.2 Comparison to previous results

Before embarking on a study of non isentropic discs in Sect. 6, we compare our results for locally isothermal discs with previous work and expectation from linear theory in this section.

Different kinds of comparisons of torque calculations with linear theory have been reported in the literature. Masset et al. (2006) compare locally isothermal three dimensional calculations with the linear results of Tanaka et al. (2002). This comparison is also reported in Papaloizou et al. (2007). In this comparison the expression for the linear torque given by

$$% T_L = - (1.35+ 0.55\alpha)q^2(H/r)^{-2}\Sigma(r_{\rm p})\Omega_{\rm K}^2r_{\rm p}^4$$ (30)

is used. On the other hand the results of two dimensional simulations using a softened planet gravitational potential, as is done here, have also been compared to those given by Eq. (30) and it has been found that an approximate match may be achieved with an appropriate choice of softening parameter (Nelson & Papaloizou 2004). We shall show that our results are fully consistent with those reported in Nelson & Papaloizou (2004) and Eq. (30) as well as our own customised linear calculations below. We note that the comparisons referred to above are done with constant kinematic viscosity coefficients, $\nu = 10^{-5}r_{\rm p}^2\Omega_{\rm p}$, and so we adopt this value here to make our comparisons. We comment that for this value of $\nu$, corotation torques are expected to be largely unsaturated.

	Figure 8: Overview of $\rho r^{1/2}$ for a q=1.26 $\times$ 10-5 planet (4  ${{M}_\oplus }$ around a Solar mass star) after 20 orbits. Left panel: locally isothermal equation of state. Right panel: adiabatic equation of state.
Open with DEXTER

We have performed locally isothermal simulations of planets with different mass ratio q for different values of $\alpha$, H/r=0.05 and the stated value of $\nu$.

In Fig. 7 we show the total torque, divided by the mass ratio squared, as a function of time for planets with q=6.34 $\times$ 10^-6 and q=1.26 $\times$ 10^-5 for two different values of $\alpha$. In these cases the torques attain almost steady values after about ten orbits with any remaining transients due to the initial conditions producing only small effects. If the response were fully linear, there would be no differences between planets of different masses. It is clear that the disc with $\alpha=3/2$, which has no vortensity gradient, is consistent with the interaction being linear. The torques agree to within 10$\%$ with the value obtained from our linear calculations, as well as with that obtained from Eq. (30). For the case with $\alpha =1/2$, the differences between the two planet masses are more pronounced We find the magnitude of the torque divided by q² to be about $10\%$ larger for the smaller mass indicating the presence of weak non linear effects. For that mass the torques again agree to within 10$\%$ with the value obtained from our linear calculations, as well as that found using Eq. (30).

In agreement with Nelson & Papaloizou (2004), torques can be fitted by Eq. (30) to within $10{-}15\%$ for the mass range considered for an appropriate softening parameter. In our case this would be close to b=0.6H/r the value we used for almost all cases. Nelson & Papaloizou (2004) adopted b=0.7H/r, which is not a significant difference in this context.

   
6 The development of positive corotation torques resulting from initial entropy gradients

In this section we focus on the development of positive torques arising from initial entropy gradients within the first horseshoe libration period, dealing with the possibility of sustaining such torques over longer time scales in later sections. We will consider various disc models, ranging from one with $\alpha =-1/2$, which has strong gradients in entropy and vortensity and a surface density that increases outwards, to discs with constant entropy and vortensity and outwardly decreasing surface density.

6.1 An illustrative case

We begin by considering a planet with q=1.26 $\times$ 10^-5 (corresponding to a 4  ${{M}_\oplus }$ planet around a Solar mass star). This closely resembles the case studied in Paardekooper & Mellema (2006a). The initial density and temperature structure is characterised by $\alpha =1/2$ and constant H/r. The power law index for the entropy is $\xi= {\rm d}(\log P/\rho^{\gamma})/{\rm d}\log r=-0.8$, which corresponds to $\gamma =1.4$. In Fig. 8 we give contour plots of the density (multiplied by r^1/2) for an isothermal equation of state and an adiabatic simulation after 20 orbits. The adiabatic case can be compared with the plot in Fig. 6 obtained from the linear theory with some heat diffusion added to regularise the density perturbation.

Figure 9: Total torque in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$ as a function of time for a q=1.26 $\times$ 10-5 planet (4  ${{M}_\oplus }$ around a Solar mass star), using a locally isothermal equation of state (solid line) and solving the full energy equation, for two different radial entropy gradients $\xi={\rm d}\log(P/\rho^{\gamma})/{\rm d}\log r$ (dotted line: $\xi =0$, dashed line: $\xi =-1.2$). All models have $\alpha =-1/2$.

Open with DEXTER

The wave structure differs slightly between the simulated adiabatic and isothermal cases. Note that the colour scale has been chosen to highlight corotation features, but still the differences between the spiral waves can be seen. The pitch angle differs such that the waves are more open and weaker in the adiabatic case. This is due to the fact that linear waves are isentropic, leading to a higher sound speed and less compression for equal forcing in the adiabatic case. This results in a reduction in the strength of the Lindblad or wave torques in adiabatic simulations compared to isothermal ones.

The two simulations in Fig. 8 also differ around corotation (near $r=r_{\rm p}$). The adiabatic simulation shows a strong, narrow density feature which is associated with a positive perturbation (light shade) for $\varphi > \varphi_{\rm p}$ and associated with a negative perturbation (dark shade) for $\varphi < \varphi_{\rm p}$. As the planets are moving upwards in Fig. 8, it is clear that this feature acts as to accelerate the planet in its orbit, causing outward migration. It is also found in the linear calculations but at reduced strength (see Fig. 6). In Sect. 7 we will see that this is because linear theory fails to take into account the horseshoe turns. We will see below that when non-linear simulations are undertaken, this corotation torque can easily come to dominate the Lindblad torques, and that the sign depends on the direction of the entropy gradient in the disc. The isothermal simulation also shows a corotation feature, which is due to the radial vortensity gradient in the disc. Although the sign of this corotation torque is positive, it is not strong enough to overcome the negative wave torque. For planets of higher mass the situation may be different (Masset et al. 2006).

In Fig. 9 we compare the time evolution of the total torque for an isothermal and an adiabatic simulation with $\alpha =-1/2$, together with a run with the same density profile, but with a temperature profile such that the entropy is constant. In the case with a negative entropy gradient, we find a positive torque on the planet after approximately 10 orbits. After that time, the torque continues to rise steadily in all cases. This is when the corotation torque is set up, see Fig. 10 which shows the torque associated with the m=2 component of the density perturbation for the case with negative entropy gradient. The torque difference between 20 and 10 orbits arises entirely at corotation. In all cases shown in Figs. 9 and 10, the corotation torque is expected to be positive due to the positive vortensity gradient and negative (or zero) entropy gradient.

However, in linear theory the magnitude of the corotation torque has been found to be not enough to produce a positive torque and thus smaller than the wave torque. We note that one finds from Tanaka et al. (2002) that the positive linear corotation torque for $\alpha =-1/2$ in the isothermal case is expected to be $54\%$ of the wave torque in magnitude making the total torque $46\%$ of the latter. Also, the linear torques are expected to reach their final value on a dynamical time scale, which is independent of the planet mass ratio, rather than over 20 orbits as is seen in Fig. 9. The precise time to reach the maximum depends on the time scale on which the planet is introduced, but it is always significantly more than a dynamical time scale. Below, we will see that this time depends on the planet mass ratio and so this behavior must be associated with non-linearities in the flow and we interpret it as arising from torques generated after horseshoe turns as described in Sect. 7 below.

After approximately 35 orbits, the torque starts to drop again in all cases. This is also an indication that corotation torques play a major role, because it is expected that they start to saturate after half a libration period, which is close to 35 orbits for this planet. We note that the torque then starts to oscillate on the libration time scale and, ultimately after the action of non-linear effects and diffusion, the disc settles down to a state with zero entropy and vortensity gradients around corotation in an average sense. Therefore, the corotation torque vanishes unless some form of diffusive or evolutionary process can restore the entropy gradient within one libration period. We will come back to this issue in Sect. 8.

	Figure 10: The cumulative torque in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$ acting on the protoplanet plotted as a function of dimensionless radius for $\rho _0 \propto r^{1/2},$ $\xi =-1.2$ and m=2 at different times. The cumulative torque is defined to be zero at the outer boundary and equal the total torque at the inner boundary.
Open with DEXTER

6.2 Non-linear evolution

Even at early times the evolution of the torque is dominated by non-linear effects. In Fig. 11 we show the time evolution of the torque for three different planet masses. Torques have been scaled by q^-2, so that if the response was purely linear all curves should lie on top of each other. The time scale for the rise towards positive values clearly depends on the mass ratio q, which would not be the case if the effect was linear. The time to reach the same degree of non-linearity is consistent with being proportional with q^-1/2, as implied by Eq. (29). Similarly, the time scale on which saturation sets in, resulting in the torque dropping below zero again, depends on q. Therefore, although there is a linear effect associated with the entropy gradient, the phenomenon that results in the torque switching sign is clearly non-linear.

This is further illustrated in Fig. 12, where we show the m=2 component of the surface density perturbation due to a q=3.17 $\times$ 10^-6 planet (corresponding to a 1  ${{M}_\oplus }$ planet around a Solar mass star) after 10 and 20 orbits. The imaginary parts (solid lines) again are directly related to the torque. At the boundaries of the domain we see clear signatures of outgoing waves. The difference between t=20 and t=10 arises near corotation only (see also Fig. 10). At t=10, the perturbation is reasonably smooth and resembles the fully developed linear case away from corotation (see Fig. 2). As time progresses, however, the corotation feature becomes narrower and deeper, until it reaches a final width. This finite final width is again indicative of non-linear behaviour, and from Fig. 12 it is clear that this non-linear structure is responsible for the cumulative torque on the planet reaching positive values. The ripples on the edges of the corotation region in Fig. 12 are of numerical origin and have to do with the strong non-linearity in the density gradient. They can be reduced by slowly ramping up the mass of the planet (which was done for the runs in Fig. 12) but in the absence of heat diffusion they never go away completely.

	Figure 11: Total torque in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$ as a function of time for 3 different planet masses and an adiabatic equation of state with $\xi =-1.2$ and $\alpha =-1/2$. The torque is divided by q2 to filter out purely linear effects.
Open with DEXTER

Figure 12: Real (dashed lines) and imaginary (solid lines) parts of the m=2Fourier component of the density perturbation due to a q=3.17 $\times$ 10-6 planet after 20 (black) and 10 (red) orbits in an adiabatic disc with $\xi =-1.2$ and $\alpha =-1/2$. The inset shows a close-up on the corotation region. A comparison with results obtained from linear theory indicates the form of a linear response apart from in the corotation region.

Open with DEXTER

	Figure 13: Imaginary parts of the m=2 Fourier component of the density perturbation due to a q=1.25 $\times$ 10-5 planet for different radial entropy profiles with $\alpha =-1/2$. For the negative entropy slope, we show three different planet masses.
Open with DEXTER

In Fig. 13 we show the imaginary part of the m=2 Fourier component of the density for the three planet masses (black lines). These allow us to make an easy estimate of the width of the horseshoe region $x_{\rm s}$. Interestingly, if we define the horseshoe region to extend to the distance where the torque is half the minimum, we find that $x_{\rm s}=r_{\rm p}\sqrt{2q/3b}$ to within a few percent. Also, it can be deduced from Fig. 11 that the time scales for which saturation sets are consistent with libration period  $8\pi/(3x_{\rm s})$ set by this value, which one may obtain from considering the gravitational and centrifugal equipotentials alone. In this context we comment that from a consideration of horseshoe streamlines, Masset et al. (2006) estimated $x_{\rm s}\sim r_{\rm p}\sqrt{8r_{\rm p}\vert\Phi_{{\rm Gps}}\vert/(3GM_*)},$ where $\Phi_{{\rm Gps}}$ is the gravitational potential of the planet evaluated at the point on the horseshoe separatrix closest to the planet (as we use this expression only to indicate a scaling, we neglect an additional enthalpy contribution to the potential). According to this we expect $x_{\rm s}=r_{\rm p}\sqrt{2qr_{\rm p}/(3L)}$ for an appropriate length scale L. As Masset et al. (2006) point out, in practical cases, this should be determined by both the softening length and the scale height, H. But note that the softening length is generally chosen so that these are closely related. For example these authors provide the estimate

$$% x_{\rm s}=1.16 r_{\rm p}\sqrt{qr_{\rm p}/H}.$$ (31)

This corresponds to L=0.495H which therefore agrees with the estimate above which used L=0.6H=b to within ten percent. They indicate that their estimate fits locally isothermal simulation results at smaller softening lengths and that this should be used to give $x_{\rm s}$ when $b\rightarrow 0$.

Incorporating the above ideas, following Masset et al. (2006) from now on we define the horseshoe width through

$$% x_{\rm s}=\lambda r_{\rm p}\sqrt{\frac{2q r_{\rm p}\vert\Phi_{{\rm Gps}}\vert}{3GM_{\rm p}}},$$ (32)

where $\lambda$ is expected to be of order unity and can be determined from simulations.

	Figure 14: Imaginary parts of the m=2 Fourier component of the density perturbation due to a q=1.25 $\times$ 10-5 planet for different resolutions with $\xi =-1.2$ and $\alpha =-1/2$, together with a run with a steeper radial density profile ( $\alpha =1/2$, with the same value for $\gamma$).
Open with DEXTER

The total torque, divided by q², depends on the mass of the planet (see Fig. 11). The reason for this is at least partly through the form of the density profile around corotation. For the torque to be proportional to q², we would expect that the mass associated with the feature, or any of its Fourier components, should scale as q. From Fig. 13, both the width and the depth of the feature are seen to depend on mass but in such a way that the area increases more weakly than q². But note that the small scale structure apparent in the density perturbation around corotation may well depend strongly on numerics, and therefore also the total torque. We will show below that this problem can be overcome by adding a small amount of thermal diffusion.

In Fig. 14 we show the corotation feature for three different resolutions. It appears that the corotation feature gets stronger with increasing resolution, and it is not clear if the torque has converged for our highest resolution. A common feature for all runs is the positive total torque. The total torque is also remarkably robust to variations in the flux limiter, showing differences of less that 5% between the minmod limiter and the soft superbee limiter (see Paardekooper & Mellema 2006b, for their definition). From Fig. 14 it is clear that there are strong features in the density that have a width of one grid cell for all resolutions. At the two highest resolutions, we resolve all relevant scales: the pressure scale height, the Hill radius of the planet and the width of the horseshoe region. The appearance of these small scale features even at these high resolutions suggests that the torque is limited by numerical diffusion. We will return to this issue below, where we discuss models with explicit heat diffusion.

	Figure 15: Total torque in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$ as a function of time for $\rho _0 \propto r^{-1/2}$ and q=1.26 $\times$ 10-5. The solid line shows an isothermal run, the dotted line an adiabatic run.
Open with DEXTER

Masset et al. (2006) noted non-linear behaviour in the barotropic case above a certain mass. Our analysis suggests that in the non-barotropic case, for sufficiently low thermal diffusivity, and for sufficiently long integration times, non-linear effects are always dominant (see Eqs. (28) and (29)). In addition, the discussion of the linear problem indicates that, as the state variables (but not their gradients) converge as $\epsilon \rightarrow 0$, non-linear effects should be much milder in the barotropic case. This is confirmed by Fig. 9, where the non-linear evolution takes the torque on the planet to positive values only in the non barotropic case.

In Fig. 13, we show the results for different entropy gradients, with the planet mass ratio being fixed at q=1.26 $\times$ 10^-5. The run with no entropy gradient resembles the Lorentz profile of the linear case (see the right panel of Fig. 2), confirming that it is the entropy gradient that drives the non-linearity. A positive entropy gradient changes the sign of the corotation feature, making the corotation torque on the planet negative. The run with no entropy gradient is intermediate between the positive and negative entropy gradient runs and indicates a linear dependence of this torque on $\xi$. However, because we measure total torques arising from a combination of physical processes, all of which vary when the disc structure is modified by altering the entropy gradient, this is not conclusive, because entropy-related torque contributions may not always be reliably identified.

The run for which we show the corotation feature in Fig. 14 behaves differently when $\alpha$ is changed from -1/2 considered so far to $\alpha =1/2$. In Fig. 15 we show the evolution of the total torque, comparing an isothermal run with an adiabatic run with $\xi=-0.8$. For this disc, the total torque never becomes positive, a indicating a strong dependence of the migration rate on the density profile of the unperturbed disc. In fact, there are three effects conspiring to make the torque negative in the $\alpha =1/2$ case. First of all, the wave or Lindblad torques are expected to be almost the same being only slightly weaker (see Tanaka et al. 2002). Second, the positive contribution from the linear barotropic corotation torque, that would apply to the locally isothermal case if the temperature variation across the corotation region may be neglected, is expected to be weaker by a factor of two because the radial vortensity gradient is less strong (Goldreich & Tremaine 1979). Not only will the linear torque be more negative, there is also no expected strong non-linear boost for the barotropic corotation torque in these circumstances. This is clear from the solid line in Fig. 15, which does not show the strong rise that was apparent in Figs. 9 and 11. For this example the final torque is also in reasonable agreement with that obtained from the linear analysis which is also indicative of less important non-linear effects at corotation. Finally the radial entropy gradient is less negative than in the case with $\alpha =-1/2$, which reduces the positive contribution of the non-barotropic corotation torque. The result is a negative torque on the planet, and therefore inward migration, albeit at a slower rate due to the positive contribution of the entropy gradient to the torque.

	Figure 16: Total torque in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$ as a function of time on a q=3.17 $\times$ 10-6 planet for different values of the heat diffusion coefficient $D_{\rm e}$, in a disc with $\alpha =-1/2$ and $\xi =-1.2$.
Open with DEXTER

   
6.3 Effect of heat diffusion

In this section we add a small thermal diffusivity into the model, characterised by the value of $D_{\rm e}$. We expect that for a high enough value of $D_{\rm e}$, in the vicinity of the corotation radius, all non-linear behaviour due a radial entropy gradient disappears and the model should behave as if it was locally isothermal. In Fig. 16 we show the time dependence of the total torque on a q=3.17 $\times$ 10^-6 planet for different values of $D_{\rm e}$ with $\xi =-1.2$ and $\alpha =-1/2$. It is clear that heat diffusion decreases the strength of the non-linear part of the torque. However, for times less than 5 orbits, the behaviour is independent of $D_{\rm e}$. This is consistent with the solution being in the linear regime during this phase (see Sect. 3 above).

For $D_{\rm e}=7$ $\times$ 10^-7 the torque during the later non-linear phase is still positive, but less so than in the adiabatic case ( $D_{\rm e}=0$). Thus this value of $D_{\rm e}$ is enough to affect the non-linear solution, which is consistent with the estimates made in Sect. 3 (see Eq. (28)). Increasing $D_{\rm e}$ by a factor of 10 gives an even weaker positive torque on the planet. The fact that after 5 orbits there is a constant shift in torque between this case and the isothermal case indicates that all non-linear effects due to the entropy gradient have disappeared. The difference between these cases is simply due to different Lindblad or wave torques, which is supported by the fact that after 5 orbits the torque ratio is about equal to $\gamma$. Together with the mildly non-linear barotropic part of the corotation torque (which is not affected by heat diffusion, but which may still involve significant density gradients in the coorbital region, see Sect. 8) this is enough to switch the sign of the torque to positive values.

	Figure 17: Real (dashed lines) and imaginary (solid lines) parts of the m=2 Fourier component of the density perturbation due to a q=3.17 $\times$ 10-6 planet after 40 orbits, for different values of the heat diffusion coefficient $D_{\rm e}$, in a disc with $\alpha =-1/2$ and $\xi =-1.2$. The inset shows a close-up on the corotation region.
Open with DEXTER

To illustrate the behaviour of the entropy-related torque, we show the m=2 Fourier component of the density perturbation in Fig. 17 for the cases $D_{\rm e}=0$, $D_{\rm e}=7$ $\times$ 10^-7 and $D_{\rm e}=7$ $\times$ 10^-6. It is clear that the effects of heat diffusion are localised around corotation, where the temperature gradients are strongest. The close-up around corotation shows a smooth Lorentz-profile for the most diffusive run, resembling the linear case in a qualitative way (see the right panel of Fig. 2). Note that Eq. (28) predicts that non-linear effects should be present for $D_{\rm e}<3.4$ $\times$ 10^-7, while from Fig. 16 we see that for a diffusion coefficient that is twice as high non-linear effects associated with the entropy gradient are still present. Bearing in mind that Eq. (28) gives only a rough estimate, and that non-linear effects may indeed set in for a somewhat lower values of $D_{\rm e}$, the agreement is reasonable.

Even in the less diffusive case the sharp features at corotation that can be seen in Fig. 17 for the $D_{\rm e}=0$ case are smoothed out. In fact, for $D_{\rm e}=7$ $\times$ 10^-7 we find the same density response for a simulation with half the resolution, and also the total torques agree to within 5$\%$. Therefore, the torque is numerically well-defined, unlike the case with $D_{\rm e}=0$. If we now compare different masses with a minimum amount of diffusion, enough to smooth out the sharp features but at the same time keep the non-linear behaviour, we find that the total torque scales almost as q². Note that in order to have the same amount of non-linearity in the problem, we need to use different values of $D_{\rm e}$ for different masses (see Eq. (28)). The results are shown in Fig. 18 for $\alpha =1/2$ and $\xi=-0.8$, where we have rescaled the time for the low-mass planet by a factor of 2. This removes the q^-1/2 scaling for the time to set up the torque. Interestingly, the torque plateau in these cases is reduced by a factor of 2 compared to the case with $D_{\rm e}=0$, being negative and about 1/2 of the linear value. This analysis supports the idea that the departure from the scaling of the torque with q² is due to non-linear effects in the corotation region.

   
7 A simple non-linear model

Because linear theory always breaks down for sufficiently small $D_{\rm e}$, it is necessary to develop a non-linear picture to understand the results of the numerical simulations presented above in Sect. 6. In this section, we discuss a very simple non-linear model that can give at least qualitative insight into the torque behaviour. It is a non-linear model because a finite width of the corotation region, in which the gas travels on horseshoe orbits, is adopted. Such a region is not present in linear theory. The model only accounts for the breakdown of linear theory and the build-up of a non-linear contribution to the corotation torque in the early stages after a planet is inserted into the disc. The later expected saturation of the torque is not considered.

	Figure 18: Total torque in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$ as a function of time for $\rho _0 \propto r^{-1/2}$ and two different planet masses. For the q=3.17 $\times$ 10-6 planet, we have used $D_{\rm e}=7$ $\times$ 10-7, for the q=1.26 $\times$ 10-5 planet we have used $D_{\rm e}=5.6$ $\times$ 10-6.
Open with DEXTER

A description of the way the model works is as follows. Viewed in the frame rotating with the planet, material exterior to the planet approaches the azimuth of the planet along approximately linear trajectories. It executes a horseshoe turn, subsequently moving away from the planet when it will eventually will have received a radial displacement equal to twice the initial distance from the orbital radius of the planet (see the streamlines plotted in Fig. 19). If the planet were absent material would shear past without a turn. Thus the effect of the horseshoe turns is to produce physical changes in the material that correspond to a radial displacement equal to the width of the horseshoe region after it has encountered the planet. The horseshoe region is separated from the remainder of the disc by a separatrix.

When the material conserves entropy as it moves, such a displacement produces a density change. For example, if the disc entropy per unit mass decreases outwards, when material is displaced inwards, the density will increase. Here we shall assume the displacement occurs such that fluid near separatrices in the horseshoe region maintains pressure balance with its surroundings and by extension that the Eulerian change in the pressure may be neglected. In this context we comment that it is found both from linear theory and non-linear simulations that the relative density changes always greatly exceed the relative pressure changes (see above and below). We note that the dynamics we have described are clearly seen in the simulations. In Fig. 19 we show a close up of the corotation region near the planet for a 4  ${{M}_\oplus }$ planet. The density structure indicates an increase inside the horseshoe region that is bounded by the inner separatrix. This is also indicated by the form of the streamlines superposed on the plot. Note that these are slightly distorted by the grid resolution. We comment that the change in density that occurs after a horseshoe turn is expected to result in a slight asymmetry in the positions of the horseshoe legs on either side of the corotation radius away from the planet. But for the relative density changes considered here (typically of order one percent) this is a very minor effect.

The horseshoe turns commence immediately after a planet is inserted into a disc and produce a torque that builds up over time in addition to that expected from linear perturbation theory. Departures from linear theory are thus produced early on. As regions close to the planet are the most important in establishing the torque, a plateau is reached fairly quickly when these regions attain a quasi-steady state. This lasts for the order of a libration period after which material recirculates to the location of the planet and the torque is expected to saturate. The model described here accounts for the attainment of the torque plateau but not the final saturation.

To describe the model quantitatively, we use coordinates  $(x,\theta)$ that are centred on the planet: $x=(r-r_{\rm p})/r_{\rm p}$, $\theta=\varphi-\varphi_{\rm p}$. As a fluid element turns to move from one side of the planet to the other, once it has moved there is ultimately a jump from x to -x; after this the streamlines are along lines of constant x. Here we neglect minor effects that may arise from the slight offset of the centre of the horseshoe region from the planet's orbital radius due to the background pressure gradient.

For simplicity, we concentrate on the part of the corotation torque that is due to the entropy gradient only. In general one also expects an additional contribution if the gradient of specific vorticity is non-zero. Because most of the torque arises close to the planet we expect the major part of it to be established before a full libration period and that this could be determined locally. Below we estimate the maximum torque arising from the initial horseshoe turns as described above.

	Figure 19: Density, multiplied by r3/2, after 10 orbits of a q=1.26 $\times$ 10-5 planet embedded in a disc with $\rho _0 \propto r^{-3/2}$, $T \propto r^{-1}$ and $\xi =-0.85$. Approximate locations of the streamlines are indicated by the black lines.
Open with DEXTER

Because entropy is conserved along streamlines and we assume the change Eulerian in pressure is zero, we can write for the perturbed density (we focus on a jump from the horseshoe leg exterior to the planet that has $\theta>0$ to the corresponding leg interior to the planet, located in the upper segment of Fig. 19):

$$% \rho'= -2x\frac{\xi}{\gamma}\rho_0,$$ (33)

where $\xi={\rm d}\log(P/\rho^{\gamma})/{\rm d}\log r$. The torque on the disc can be found by integrating over the disc domain that has participated in a horseshoe turn:

$$% \mathcal T = -2Hr_{\rm p}^2 \int \rho' \frac{\partial \Phi_{{\rm Gp}}}{\partial \theta}{\rm d}\theta {\rm d}x.$$ (34)

The integral over $\theta$ is made simple by adopting as integration limits $\theta = 0$ and $\theta$ large enough such that the force due to the planet is small, which leaves us with:

$$% \mathcal T=4H\rho_0r_{\rm p}\frac{\xi}{\gamma}{\rm G}M_{\rm p}\int_0^{x_{\rm s}}\frac{1}{\sqrt{x^2+b^2}}x{\rm d}x,$$ (35)

where we have approximated the planet potential by:

$$% \Phi_{{\rm Gp}}= -\frac{{\rm G}M_{\rm p}}{r_{\rm p}\sqrt{x^2+\theta^2+b^2}}\cdot$$ (36)

Thus we arrive at:

$$% \mathcal T=4{\rm G}M_{\rm p}H \rho_0 \frac{\xi}{\gamma}\frac{x_{\rm s}^2}{r_{\rm p}b+\sqrt{(r_{\rm p}b)^2+x_{\rm s}^2}}\cdot$$ (37)

The torque on the planet has the opposite sign and may be written in the form

$$% \mathcal{T}_{\rm p}= -\frac{4\xi x_{\rm s}^2}{\gamma q \lef... ...}^2b^2}\right)} q^2\Sigma_{\rm p}r_{\rm p}^3 \Omega_{\rm p}^2,$$ (38)

where we have introduced a factor of 2 to account the additional contribution arising from the horseshoe legs on the other side of the planet. Thus we assume a symmetry such that these behave in a corresponding manner but with a density deficit rather than increase produced near the separatrix after a turn, so also producing a positive torque.

We comment here that the above torque depends on the softening parameter, b through evaluation of the planet potential at the point on the separatrix closest to the planet, here taken to be the actual planet location. We may accordingly replace it with the magnitude of this potential denoted by $\vert\Phi_{{\rm Gps}}\vert= GM_{\rm p}/(r_{\rm p}b)$. Then the torque may be written in the form

$$% \mathcal{T}_{\rm p}= \frac{-4\xi x_{\rm s}^2 \vert\Phi_{{\r... ...)^2}+1}\right)} q^2\Sigma_{\rm p}r_{\rm p}^4 \Omega_{\rm p}^2.$$ (39)

In this form the expression for the torque may be used in the situation where the point on the separatrix closest to the planet is separated or offset from the planet position by an angular extent  $\theta_{\rm s}$. From Eq. (36) we see that this is equivalent to replacing bby $\sqrt{b^2+\theta_{\rm s}^2}$. Now we use Eq. (32) to substitute for $x_{\rm s}$. The expression for the torque then becomes

$$% \mathcal{T}_{\rm p}=\frac{-\frac{8}{3}\frac{\xi}{\gamma}\le... ...}\right)^3+1}} q^2\Sigma_{\rm p}r_{\rm p}^4 \Omega_{\rm p}^2,$$ (40)

where $\lambda$ is the dimensionless parameter defining the horseshoe width that is expected to be of order unity that was introduced above (see Eq. (32)).

Figure 20: Left panel: total torque, in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$, on a q=1.26 $\times$ 10-5 planet embedded in a disc with $\rho _0 \propto r^{-3/2}$, $T \propto r^{-1}$ and $\xi =-0.85$. A small thermal diffusion coefficient $D_{\rm e}=10^{-7}$ was used to regularise the solution near corotation. Results for two values of the softening parameter b are shown in black and red, together with a simulation with b=0.03 and $D_{\rm e}=0$ (green line). Right panel: azimuthal variation of the density at r=0.99 at three different times for b=0.03 and $D_{\rm e}=10^{-7}$.

Open with DEXTER

Equation (38) predicts a torque that is proportional to the entropy gradient of the unperturbed disc, and is positive for negative entropy gradients. The torque that is reached is comparable in magnitude to the isothermal wave, or Lindblad, torque for the expected value of $b \approx H/r$ (Tanaka et al. 2002). The torque, as given by Eq. (38), does not diverge for $b\rightarrow 0$. With any separation or offset of the type described above being neglected, setting b=0 in Eq. (38) gives

$$% \mathcal{T}_{\rm p}= -\frac{4\xi x_{\rm s}}{\gamma } q\Sigma_{\rm p}r_{\rm p}^3 \Omega_{\rm p}^2.$$ (41)

Note that in this limit the only stream parameter that enters this expression is the horseshoe width $x_{\rm s}$. The torque is then increased by a factor $x_{\rm s}/( r_{\rm p}b+\sqrt{x_{\rm s}^2+r_{\rm p}^2b^2})$ compared to the case $b \ne 0$. Using Eq. (31) to estimate $x_{\rm s}$, this amounts to a factor of 4when b=0.6H/r for q=1.26 $\times$ 10^-5.

However, its use in this limit is inappropriate. Following Masset et al. (2006), we note that when vertical stratification is considered it would be expected to produce an effective softening parameter $b=z/r_{\rm p}$ for material moving at height z. Adopting this and performing a vertical average of Eq. (38) for the above example would indicate $b \sim 0.5H/r$, being comparable to values adopted here.

7.1 Time development of the non-linear torque

We studied the time dependence of the disc response as a result of the insertion of a planet of 4  ${{M}_\oplus }$into a disc with $\rho _0 \propto r^{-3/2}$ and $T \propto r^{-1}$. The model was adiabatic with $D_{\rm e}=0$ and $\xi =-0.85$, in this case corresponding to $\gamma =1.1$. Note that, as mentioned in Paardekooper & Mellema (2008), there is no smooth transition to the locally isothermal torque for $\gamma \rightarrow 1$, because in the latter case entropy is not conserved along streamlines (if there is a radial temperature gradient). In fact, for a negative density slope, a lower value of $\gamma$ tends to make the entropy gradient more negative, and therefore the associated corotation torque more positive.

Unlike for the other simulations considered in this paper, in order to study the time development of the non-linear torque in a way corresponding to the discussion of the simple model, the planet was initially immersed in the disc with its full mass. Thus there was no period of slow build up for these cases. From the streamlines illustrated in Fig. 19 and discussed above in Sect. 7, the inner and outer separatrices for a model with b =0.03, away from the planet, lie on circles centred on the central object with dimensionless radii 0.982 and 1.012 respectively.

In order to study the development of the high-density region leading the planet which is responsible for the non-linear torque, we plot the density as a function of azimuthal angle, $\varphi$, as viewed along the radius r=0.99 after 5, 10 and 20 orbits in the right panel of Fig. 20. The planet is located at $\varphi = \pi$ and for higher values of $\varphi$, this radial location lies in the required high density region. The left panel of the same figure shows the total torque as a function of time (black line). A small thermal diffusivity, $D_{\rm e}=10^{-7}$ was used to regularise the solution around corotation. For comparison, we also show the total torque for a simulation with $D_{\rm e}=0$ (green line). The torque in this case rises to a slightly higher value at t=10, followed by a slight decay. This effect is due to the sudden introduction of the planet into the disc, and it can be reduced by introducing a small thermal diffusivity, as shown in Fig. 20, or by slowly ramping up the planet's mass.

We see that as time progresses an approximately uniform density increase by a factor of about 1+5 $\times$ 10^-2 is established. This is entirely consistent with the factor $\delta \rho/\rho = 1.55 x/r$, with x being the distance to the coorbital radius, expected from the discussion above. The uniform density increase is also consistent with a uniform pressure equal to the unperturbed value being attained at the fixed radius.

Let us now check that the time development is consistent with fluid elements moving from one side of the horseshoe region to the other. Assuming this to be the case, taking account only of Keplerian shear so neglecting any time required to cross close to the planet, it would take 5 orbits to set up the high density region out to $\varphi =3.4$, corresponding to 5 disc scale heights from the planet, 10 orbits to set it up out to $\varphi =3.64$ and 20 orbits out to $\varphi=4.14$. However, Fig. 20 shows that the density profile is not fully established after 5 orbits, but that after 10 and 20 orbits the profile is indeed established out to the expected locations. This indicates that the evolution is relatively slow only in the regions close to the planet. This could be related to passage close to the separatrices where low velocities are expected.

Note too that half of the total contribution to the corotation torque we have considered here is attained after only 5 orbits which is much less than the libration period of about 75 orbits which is again consistent with most of the contribution coming from close to the planet. The magnitude of the non-linear contribution to the torque predicted by Eq. (38) is approximately 1000 in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$, which agrees perfectly with the difference between the minimum (or also approximately the linear torque) and the maximum torque in the left panel of Fig. 20. This suggests that there is no additional torque cut-off close to the planet, since one was not included in the non-linear model presented here. From Fig. 19 can be seen that indeed the corotation feature extends to about a smoothing length from the planet in azimuth. Note too that although the associated relative density change is only a few percent, on account of the small radial scale the relative change to the density gradient is at least of order unity. Additional non-linearity not represented in the simple model may be present and possibly be responsible for the variations of the ultimate ratio of torque to mass ratio squared seen in Fig. 11. The fact that these ratios are affected by small diffusion coefficients in the manner illustrated in Fig. 18 indicates localised effects in the coorbital zone are responsible. A possibility is that because density jumps across the horseshoe separatrix are higher for larger masses, because the horseshoe region is wider, non-linear effects are also relatively stronger.

Also shown in the left panel of Fig. 20 is the total torque for a run with a smaller value of the smoothing length b=0.02. From Eq. (40), with $\lambda=1$, we expect the non-linear torque to be a factor 2 stronger in this case, and this is entirely consistent with Fig. 20. We point out that from this trend the total torque would be expected to be always negative for $b \ge H/r$.

	Figure 21: Long-term evolution of the total torque, in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$, on a q=1.26 $\times$ 10-5 planet in a disc with $\rho _0 \propto r^{-3/2}$, $T \propto r^{-1}$ and $\xi =-0.85$, for three different thermal diffusivities.
Open with DEXTER

   
8 Long-term evolution

We now turn to the question of the behaviour of the torque over longer time scales, focusing on a disc with no vortensity gradient. We initialise the temperature $T \propto r^{-1}$, and again take $\xi =-0.85$. As we need to follow the evolution of the torque for several libration cycles, which, for a q=1.26 $\times$ 10^-5, amounts to several hundred orbital periods, we lowered the resolution by a factor 2 in all directions to keep the problem computationally tractable. As mentioned in Sect. 6.3, for a finite value for the thermal diffusion coefficient the torques agreed very well for the first 30 orbits. We focus on a q=1.26 $\times$ 10^-5 planet throughout, and, to keep matters as simple as possible, do not ramp up the mass of the planet over the first 3 dynamical time scales as was done in the previous sections. Again, for a finite thermal diffusion coefficient, this does not influence the torques at all, while for $D_{\rm e}=0$ differences remain well within 10$\%$.

In order for a diffusive process to prevent the corotation torque from saturating, it is necessary that the diffusion time scale across the horseshoe region is smaller than the libration time scale, or:

$$% \frac{x_{\rm s}^2}{r_{\rm p}^2 \Omega_{\rm p}D_{\rm e}} < \... ... = \frac{8}{3} \frac{r_{\rm p}}{x_{\rm s}}\Omega_{\rm p}^{-1}.$$ (42)

Using our estimate Eq. (32) for $x_{\rm s}$, which has been shown to be a good approximation, we find for q=1.26 $\times$ 10^-5 and b=0.03 that $D_{\rm e}> 3/8(x_{\rm s}/ r_{\rm p})^3 \approx 3$ $\times$ 10^-6 to allow thermal diffusion the possibility of restoring the entropy gradient within one libration cycle. On the other hand, according to Fig. 18, for $D_{\rm e}=5.6$ $\times$ 10^-6 the disc response is still sufficiently non-linear that the torque is changed considerably. Therefore, there exists a range of $D_{\rm e}$ for which there may be a possibility of a sustained positive torque contribution from corotation effects related to the entropy gradient. However, thermal diffusion only affects the temperature gradient directly, and it may be the case that a consequent adjustment of the density in the horseshoe region allows it to settle towards a state of zero torque contribution.

In Fig. 21, we show the long-term evolution of the total torque for three values of $D_{\rm e}$. The case of zero diffusion shows the clearest libration cycles, diminishing in amplitude until the system is expected to settle into a state of vanishing entropy gradient in the mean and therefore zero corotation torque. This behaviour is very similar to that described by Ward (2007) for the barotropic case with a radial vortensity gradient. Interestingly, both diffusive runs approach the same value of the torque, which indicates that indeed the density is redistributed in such a way that the corotation torque vanishes.

The final negative torque of magnitude $\sim$1200 compares with the linear value (based on the original density distribution) of $\sim$700. For comparison the result of Tanaka et al. (2002) obtained by use of Eq. (30) for the original density distribution gives a magnitude of $\sim$900. As indicated above this was in much better agreement with the torque obtained for the locally isothermal case with the same original density distribution. The most likely reason for the greater discrepancy in this case is that the original density gradient in the coorbital region and its neighbourhood has been modified as indicated in Fig. 22 in an asymmetric manner.

In the final state, the corotation region contains an asymmetric density structure. This can be seen in the two leftmost panels of Fig. 22, where we show the density after 160 orbits for the cases with low and high diffusion. In the leftmost panel, it is clear that the low density ridge created in the horseshoe turn below the planet returns on the other side. For the higher value of $D_{\rm e}$, a similar situation arises, but since diffusion for this case is actually strong enough to affect the density jump directly, the effect is harder to spot. From Fig. 21, however, we see that in this case approximately the same total torque is attained.

All of this indicates that in order to sustain the positive torque, we need to restore the initial density profile as well as the initial temperature profile in the coorbital region. This can be done by including a small kinematic viscosity. In Fig. 23, we show the torque evolution for the same models as in Fig. 21, but now including a viscosity $\nu =10^{-6}$  $r_{\rm p}^2 \Omega _{\rm p}$. The mass flow induced by this small value of $\nu$ is not strong enough to affect the results; in essence the viscosity only acts on the sharp density features near corotation. We see that for a small viscosity in combination with a small thermal diffusivity the positive torque can be sustained. The magnitude strongly depends on the magnitude of the diffusion coefficients, as is illustrated by the dashed line in Fig. 23. A simulation with no heat diffusion shows similar behaviour as the runs depicted in Fig. 21, but the final value of the torque agrees better with linear theory, presumably because of the smoother density profile. This indicates that indeed we need to restore both the original temperature gradient (through heat diffusion) and smooth the density gradient (though a small viscosity) in the coorbital region in order to sustain the positive torque.

These results indicate that the entropy-related corotation torque alone is indeed able to reverse the sign of the total torque over longer time scales than the libration period. Although the torque in Fig. 23 is only slightly positive at later times, note that this model represents a case for which the entropy-torque alone has to overcome the full Lindblad torque. A higher value of $\gamma$, which reduces the Lindblad torque, and a flatter density profile, which gives rise to a positive corotation contribution due to the vortensity gradient, may help to push the torque to higher positive values. On the other hand increasing $\gamma$ may act to reduce the entropy gradient. The investigation of how the vortensity and entropy corotation torques operate together will be the subject of a future study.

	Figure 22: Density, multiplied by r3/2, for three models for q=1.26 $\times$ 10-5, $\rho _0 \propto r^{-3/2}$, $T \propto r^{-1}$ and $\xi =-0.85$ after 160 orbits, focusing on the horseshoe turn leading the planet. Left panel: $D_{\rm e}= 10^{-6}$, $\nu =0$. Middle panel: $D_{\rm e}= 10^{-5}$, $\nu =0$. Right panel: $D_{\rm e}=\nu =10^{-6}$.
Open with DEXTER

	Figure 23: Long-term evolution of the total torque, in units of $q^2 \Sigma (r_{\rm p})\Omega _{\rm K}^2r_{\rm p}^4$, on a q=1.26 $\times$ 10-5 planet in a disc with $\rho _0 \propto r^{-3/2}$, $T \propto r^{-1}$ and $\gamma =1.1$, for three different thermal diffusivities. All models have a small kinematic viscosity of $\nu =10^{-6}$  $r_{\rm p}^2 \Omega _{\rm p}$.
Open with DEXTER

   
9 Discussion

In this paper, we have considered the torques induced by disc-planet interaction when low-mass planets are introduced into protoplanetary discs with an initial radial entropy gradient and either solved the energy equation or adopted an adiabatic or locally isothermal equation of state. We first considered the linear theory of the interaction. As the response is singular it needs to be regularised. This was done by adopting the well-known Landau prescription, together with a small thermal diffusivity in some cases. This allows migration torques on the planet to be calculated, but the density response diverges as dissipative effects are reduced to zero implying that non-linear effects will ultimately be important. We estimated that non-linear effects would be significant for dimensionless thermal diffusivities $D_{\rm e} < 10^{-(5-6)}$ and that in such cases the linear theory would only be valid for some finite time after insertion of the planet. This was confirmed by the non-linear simulations. It was found from these that a non-linear density structure developed around corotation at early times, that in some cases with a negative entropy gradient, led to a reversal of the sign of the torque leading to outward migration. This positive torque survives for a horseshoe libration period after which saturation occurs and the torque again becomes negative.

We presented a simple non-linear model to account for the non-linear coorbital density structures as being produced after the horseshoe turns undergone by material while conserving its entropy in Sect. 7. Numerically we found that in cases with no diffusivity, the coorbital density structure showed variations on the grid scale and therefore in order to get numerical convergence at limited resolution, a small non-zero thermal diffusivity needs to be included. Increasing the thermal diffusivity brings the migration torque closer to the linear prediction although some non-linearity associated with the torque connected with the specific vorticity gradient may remain.

The non-linear and temporary nature of the physical effect that causes the torque to become positive makes it inappropriate to attempt to deduce an applicable linear torque formula, as was done in the barotropic case (Goldreich & Tremaine 1979). All that can be said is that this non-linear part of the torque depends on the entropy gradient. We have found that, for a disc with H/r=0.05, a density gradient less steep than $\rho _0 \propto r^{-3/2}$ results in a positive torque on the planet. However, this depends also on the adopted value of the smoothing length, the value of $\gamma$ and the amount of thermal diffusion. A smaller value for b will enhance such non-linear effects, while thermal diffusion acts in the opposite direction.

Our simple non-linear model of Sect. 7 does a good job in qualitatively describing the results of the numerical simulations. Because the three-dimensional simulations of Paardekooper & Mellema (2006a) show the same slow rise of the torque, which is indicative of non-linear effects, it is likely that a similar model applies in 3D.

The strong dependence of low-mass planet migration on the radial entropy gradient in the gas disc raises two important questions. First of all: do realistic, self-consistent disc models with detailed radiative transfer allow for negative entropy gradients? Second: is the thermal diffusivity low enough to support non-linearity at corotation? The second question was already answered by the simulations of Paardekooper & Mellema (2006a), which showed that, at least in the inner parts of protoplanetary discs, radiative transport does not provide enough diffusion to restore linearity. As opacity in general decreases with radius, the situation will be different in the outer parts of the disc.

The first question is more difficult to answer. Detailed models of D'Alessio et al. (1998,1999,2001) show a rich variety in temperature and density profiles, with indeed several cases with strongly negative entropy gradients, especially near opacity jumps. These may be preferential radii for halting planet migration, a conclusion also reached by Menou & Goodman (2004) for different reasons. Although it is beyond the scope of this paper to investigate the occurrence of negative entropy gradients in detail, it is important to realise that they do appear in detailed models.

In the inviscid limit, the corotation region is a closed system, and can therefore not continue to exert a torque on the planet indefinitely. The first signs of saturation of the torque are already visible in, for example, Fig. 9 and the complete process is illustrated in Fig. 21. In order to sustain this torque for longer times, the physical conditions in the coorbital region need to be regenerated. Heat diffusion can act to restore the original entropy gradient, but we have found that the form of the surface density gradient also needs to be smoothed locally, which can be achieved by including a small kinematic viscosity.

At this point we stress the sensitivity of torque behaviour derived from effects at corotation to small values of thermal diffusivity and viscosity. This has already been indicated by Masset (2002) who found that small values of viscosity acting in conjunction with effects at corotation could have significant effects on the migration of planets but in a larger mass range than that considered here.

   
10 Summary and conclusion

In this paper, we have studied the effect of including the energy equation in the theory of disc planet interactions undergone by low-mass planets. We have shown that the wave torque on the planet is reduced by a factor $\gamma$ compared to locally isothermal discs. It is the corotation torque, however, that is the most strikingly different from the isothermal case.

We have found that migration behaviour depends critically on the radial entropy gradient in the disc. The linear estimate of the corotation torque depends on this, but the entropy gradient is also responsible for a departure from linearity, giving rise to a strong corotation torque that is positive for a negative entropy gradient. For a strong enough entropy gradient these corotation torques can dominate over the wave torques. This can be understood in a qualitative way by considering simplified horseshoe dynamics. Thermal diffusion reduces non-linearity, and for high enough values of the thermal diffusion coefficient $D_{\rm e}$ we recover modified linear wave torques and a mildly non-linear barotropic corotation torque.

Although for a planet on a fixed orbit and with no viscosity, the positive torques saturate after a libration period, we have shown that, for a small thermal diffusivity and viscosity the entropy-related corotation torque offers the possibility of a sustained reversal of the direction or stalling of migration of low-mass protoplanets. Of course extensive future investigations that take into account three dimensional effects and allow the planet to migrate are required in order to evaluate it.

References

 
• Abramowitz, M., & Stegun, I. A. 1970, Handbook of Mathematical Functions (New York: Dover) In the text  
• Artymowicz, P. 2004, in Debris Disks and the Formation of Planets, ed. L. Caroff, L. J. Moon, D. Backman, & E. Praton, ASP Conf. Ser., 324, 39 In the text  
• Bate, M. R., Lubow, S. H., Ogilvie, G. I., & Miller, K. A. 2003, MNRAS, 341, 213 In the text  
• D'Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411 In the text  
• D'Alessio, P., Calvet, N., Hartmann, L., Lizano, S., & Cantó, J. 1999, ApJ, 527, 893 In the text  
• D'Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321 In the text  
• D'Angelo, G., Henning, T., & Kley, W. 2002, A&A, 385, 647 In the text  
• D'Angelo, G., Henning, T., & Kley, W. 2003a, ApJ, 599, 548 In the text  
• D'Angelo, G., Kley, W., & Henning, T. 2003b, ApJ, 586, 540 In the text  
• Eulderink, F., & Mellema, G. 1995, A&AS, 110, 587 In the text  
• Goldreich, P., & Sari, R. 2003, ApJ, 585, 1024 In the text  
• Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 In the text  
• Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 In the text  
• Hawley, J. F. 2000, ApJ, 528, 462 In the text  
• Jang-Condell, H., & Sasselov, D. D. 2005, ApJ, 619, 1123 In the text  
• Klahr, H., & Kley, W. 2006, A&A, 445, 747 In the text  
• Lin, D. N. C., & Papaloizou, J. 1979, MNRAS, 186, 799 In the text  
• Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846 In the text  
• Masset, F. S. 2002, A&A, 387, 605 In the text  
• Masset, F. S., & Papaloizou, J. C. B. 2003, ApJ, 588, 494 In the text  
• Masset, F. S., D'Angelo, G., & Kley, W. 2006, ApJ, 652, 730 In the text  
• Menou, K., & Goodman, J. 2004, ApJ, 606, 520 In the text  
• Morohoshi, K., & Tanaka, H. 2003, MNRAS, 346, 915 In the text  
• Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849 In the text  
• Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18 In the text  
• Paardekooper, S.-J., & Mellema, G. 2006a, A&A, 459, L17 In the text  
• Paardekooper, S.-J., & Mellema, G. 2006b, A&A, 450, 1203 In the text  
• Paardekooper, S.-J., & Mellema, G. 2008, A&A, 478, 245 In the text  
• Papaloizou, J. C. B., Nelson, R. P., Kley, W., Masset, F. S., & Artymowicz, P. 2007, in Protostars and Planets V, 655 In the text  
• Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 In the text  
• Ward, W. R. 1997, Icarus, 126, 261 In the text  
• Ward, W. R. 2007, in Lunar and Planetary Institute Conference Abstracts, 38, 2289 In the text