# A torque formula for non-isothermal Type I planetary migration - II. Effects of diffusion

S.-J. Paardekooper, C. Baruteau, W. Kley

► To cite this version:

S.-J. Paardekooper, C. Baruteau, W. Kley. A torque formula for non-isothermal Type I planetary migration - II. Effects of diffusion. Monthly Notices of the Royal Astronomical Society, 2011, 410 (1), pp.293-303. 10.1111/j.1365-2966.2010.17442.x. hal-02369776

HAL Id: hal-02369776

https://hal.science/hal-02369776

Submitted on 30 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A torque formula for non-isothermal Type I planetary migration – II. Effects of diffusion

S.-J. Paardekooper,<fnref n="1" />* C. Baruteau<fnref n="2" /> and W. Kley<fnref n="3" />

Accepted 2010 July 28. Received 2010 July 28; in original form 2010 May 12

## ABSTRACT

We study the effects of diffusion on the non-linear corotation torque, or horseshoe drag, in the two-dimensional limit, focusing on low-mass planets for which the width of the horseshoe region is much smaller than the scaleheight of the disc. In the absence of diffusion, the non-linear corotation torque saturates, leaving only the Lindblad torque. Diffusion of heat and momentum can act to sustain the corotation torque. In the limit of very strong diffusion, the linear corotation torque is recovered. For the case of thermal diffusion, this limit corresponds to having a locally isothermal equation of state. We present some simple models that are able to capture the dependence of the torque on diffusive processes to within 20 per cent of the numerical simulations.

**Key words:** planets and satellites: formation – planet–disc interactions.

## 1 INTRODUCTION

Low-mass planets, embedded in circumstellar gas discs, are subject to orbital migration (Goldreich & Tremaine 1979, 1980) through disc tides. Until they can acquire a massive gas envelope, planets comparable in mass to the Earth can therefore significantly move away from their place of birth. It was long thought that this Type I migration was always directed inwards, and alarmingly fast (Ward 1997; Tanaka, Takeuchi & Ward 2002), which suggested that the survival of long-period planets could be problematic within this theory. Recent work has indicated that the situation may be different if the disc's thermodynamics is taken into account properly (Paardekooper & Mellema 2006b; Kley & Crida 2008; Kley, Bitsch & Klahr 2009), in which case outward migration is possible. It was subsequently realized that this migration behaviour is due to an entropy-related corotation torque (Paardekooper & Mellema 2008; Baruteau & Masset 2008a) and is non-linear in nature (Paardekooper & Papaloizou 2008).

In Paardekooper et al. (2010, hereafter Paper I), we studied the non-linear corotation torque, or horseshoe drag (Ward 1991), in the two-dimensional adiabatic limit. It was shown that this horseshoe drag is dominated by a term proportional to a radial entropy gradient in the unperturbed disc and can lead to outward migration if the entropy decreases outwards. In the absence of any diffusion, however, the corotation torque is a transient phenomenon and prone to saturation (Quinn & Goodman 1986; Ward 2007). The necessary gradients [vortensity in the (locally) isothermal limit, entropy and

vortensity for the non-isothermal case] within the horseshoe region need to be restored on a libration time-scale in order for the corotation torque to be sustained (Masset 2001; Ogilvie & Lubow 2003). It was shown in Paardekooper & Papaloizou (2008) that for non-isothermal discs, both thermal diffusion and viscosity are needed to sustain the non-linear corotation torque.

Apart from possibly keeping the corotation torque unsaturated, diffusion of heat and momentum can also, in the limit of high diffusion, act to reduce the non-linearity of the torque. It was already noted in Masset (2002) that there exists a cut-off for the isothermal horseshoe drag at high viscosity. Paardekooper & Papaloizou (2009a) showed that in this limit, the corotation torque returns to its linear value. For non-barotropic discs, the limit of high thermal diffusion effectively corresponds to the locally isothermal case.

Depending on radial profiles of density, temperature, viscosity and thermal diffusion, the latter mainly determined by temperature and opacity, a wealth of Type I migration behaviour is possible therefore. While the sign of the unsaturated torque is independent of planet mass, saturation does depend on the mass of the planet, which can lead to different migration directions for different planet masses. Moreover, since the level of saturation effectively determines whether the non-linear corotation torque can dominate Type I migration, it is important to study the effects of diffusion on low-mass planet migration.

In this paper, we study Type I migration in non-barotropic discs in the presence of thermal and viscous diffusion. We start in Section 2 by reviewing our disc models and numerical methods. In Section 3, we briefly review the results of Paper I on the unsaturated values of the torque. We then move on to the effects of viscous diffusion in Section 4 by considering the simple isothermal case, first looking at

<footnote n="1">¹DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA</footnote>

<footnote n="2">²UCO/Lick Observatory, UC Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="3">³Institut für Astronomie & Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany</footnote>

<footnote>*E-mail: s.paardekooper@damtp.cam.ac.uk</footnote>
saturation and then moving to the limit of high diffusion in which the linear torque is recovered. We extend these results to the non-barotropic case in Section 5. We discuss our results in Section 6 and conclude in Section 7.

## 2 NUMERICAL METHODS

For the basic equations and details of the equilibrium disc models, we refer the reader to Paper I. As in Paper I, we study two-dimensional power-law discs in surface density $\Sigma$ (with index $-\alpha$) and temperature $T$ (with index $-\beta$). We either study isothermal discs, for which the pressure $P = c_s^2 \Sigma$, with $c_s$ being the isothermal sound speed, or solve the energy equation and take a perfect gas equation of state, $P = (\gamma - 1)\epsilon$, where $\epsilon$ is the internal energy density and $\gamma$ the adiabatic exponent. The scaleheight of the disc $H = c_s/\Omega$ is taken to be small compared to the orbital radius. Typically, $h = H/r = 0.05$ at the location of the planet.

We have used the same numerical methods as in Paper I. Fast Advection in Rotating Gaseous Objects (FARGO; Masset 2000a,b) is based on the van Leer upwind algorithm on a staggered mesh. The time-step calculation takes out the average angular velocity, which is dominated by Keplerian rotation, so that the time-step is limited by the residual angular velocity rather than the much larger Keplerian velocity. An energy equation solver was introduced in Baruteau & Masset (2008a). Diffusion in the energy equation is obtained by adding a term proportional to the Laplacian of the entropy:

$$ \frac{\partial \epsilon}{\partial t} + \nabla \cdot \epsilon \mathbf{v} = -(\gamma - 1)\epsilon \nabla \cdot \mathbf{v} + \chi \epsilon \nabla^2 \log \frac{P}{\Sigma^\gamma}, \quad (1) $$

where $\mathbf{v}$ is the velocity and $\chi$ a thermal diffusion coefficient. This equation is strictly equivalent to an advection-diffusion equation for the gas entropy $p/\Sigma^\gamma$. Using this form of the energy equation is a common approach in studies of turbulent energy transport (see Shakura, Sunyaev & Zilitinkevich 1978). We take $\chi$ to be a constant throughout the disc.

Roe solver for Disc Embedded Objects (RODEO; Paardekooper & Mellema 2006a) is a finite volume method, using an approximate Riemann solver (Roe 1981) to calculate fluxes between computational cells. Both codes were shown to yield similar values for the unsaturated horseshoe drag in the adiabatic case (see Paper I), with RODEO usually producing slightly larger values (up to 10 per cent). In RODEO, explicit heat diffusion is implemented by solving

$$ \frac{\partial T}{\partial t} = -\nabla \cdot (\chi \nabla T). \quad (2) $$

For RODEO, we choose $\chi$ to be a function of radius such that the initial temperature profile is a stationary solution. We have found that this only becomes important when the diffusion time-scale over the horseshoe width becomes comparable to the dynamical time-scale. Equation (2) is formally only correct for constant density, but therefore has the advantage of affecting the temperature only, while viscosity affects only density structures. This way, thermal and viscous diffusion can be treated separately.

The two different approaches for obtaining diffusion were found to yield similar results, with differences less than 10 per cent. In both cases, the diffusion equations are solved explicitly, which leads to a restriction of the time-step. Note that we have not included viscous dissipation in the equations: we assume that the background temperature profile is set by an equilibrium between viscous heating and radiative cooling (Kley & Crida 2008) and that viscous and thermal diffusion will act to restore the original profile.

We take the viscosity $\nu$ to be a power law in radius such that the initial surface density profile is a stationary solution. It is easy to verify that this requires $\nu \propto r^{\alpha-1/2}$.

## 3 UNSATURATED TORQUES

Throughout this paper, we will work with a fixed gravitational smoothing length for the planet's potential, $b/h = 0.4$. We discuss possible effects of different smoothing lengths in Section 5.7. In Paper I, we found the following unsaturated torques for this smoothing length, which we repeat here for convenience: the Lindblad torque:

$$ \gamma \Gamma_L / \Gamma_0 = -2.5 - 1.7\beta + 0.1\alpha, \quad (3) $$

the barotropic part of the horseshoe drag

$$ \gamma \Gamma_{hs,baro} / \Gamma_0 = 1.1 \left( \frac{3}{2} - \alpha \right) \quad (4) $$

and the entropy-related part of the horseshoe drag

$$ \gamma \Gamma_{hs,ent} / \Gamma_0 = 7.9 \frac{\xi}{\gamma}, \quad (5) $$

where $\xi = \beta - (\gamma - 1)\alpha$ is the negative of the power-law index of the entropy. The barotropic part of the linear corotation torque reads as

$$ \gamma \Gamma_{c,lin,baro} / \Gamma_0 = 0.7 \left( \frac{3}{2} - \alpha \right), \quad (6) $$

and the entropy-related part of the linear corotation torque is given by

$$ \gamma \Gamma_{c,lin,ent} / \Gamma_0 = \left( 2.2 - \frac{1.4}{\gamma} \right) \xi. \quad (7) $$

As in Paper I, all torques are normalized to $\Gamma_0 = (q/h)^2 \Sigma_p r_p^4 \Omega_p^2$, where subscripts $X_p$ indicate quantity $X$ evaluated at the orbital radius of the planet $r = r_p$.

## 4 ISOTHERMAL DISCS

### 4.1 A simple saturation model

Masset (2001) discussed a model for corotation torque saturation. There are a few drawbacks to this approach, however, the major one being that the dependence of the corotation torque on background gradients has to be put in by hand. This actually leads to three possible saturation laws. Another difficulty is that since Masset (2001) considers a constant viscosity $\nu$ and a constant background surface density, the unperturbed model has a radial velocity profile (i.e. it is not in equilibrium). At large viscosities, the co-orbital region for low-mass planets can be deformed, resulting in a Type III flow structure. Below, we present a version of the Masset (2001) model that is tailored to the problem at hand.

Assuming Keplerian rotation everywhere, the evolution of the surface density is given by

$$ \frac{\partial \Sigma}{\partial t} = \frac{3}{r} \frac{\partial}{\partial r} \left( \frac{1}{r\Omega} \frac{\partial}{\partial r} (r^2 \nu \Sigma \Omega) \right) + \frac{|\Omega - \Omega_p|}{2\pi} \left( \Sigma' \frac{\Omega}{\Omega'} - \Sigma \right) \Pi \left( \frac{r - r_p}{2r_p x_s} \right), \quad (8) $$

where the second term on the right-hand side accounts for changes due to the horseshoe turns (see Masset 2001). Primed quantities indicate that they should be evaluated at $r' = 2r_p - r$, i.e. on the
other side of the horseshoe turn, and $\Pi$ denotes the rectangular function:

$$ \Pi(t) = \begin{cases} 1 & |t| < \frac{1}{2}, \\ 0 & |t| \ge \frac{1}{2}. \end{cases} \qquad (9) $$

The use of the rectangular function ensures that changes due to horseshoe turns only occur for $|r - r_p|/r_p < x_s$, where $x_s$ is the non-dimensional half-width of the horseshoe region. Changing variables to $y = \Sigma/\Omega$ and $x = (r - r_p)/r_p$, we have

$$ \frac{\partial y(x, t)}{\partial t} = \frac{3v_p}{r_p^2} (1+x)^{1/2} \frac{\partial}{\partial x} \left( (1+x)^{1/2} \frac{\partial}{\partial x} \left( (1+x)^{\alpha-3/2} y(x, t) \right) \right) + \frac{|\Omega - \Omega_p|}{2\pi} (y(-x, t) - y(x, t)) \Pi \left( \frac{x}{2x_s} \right). \quad (10) $$

Note that $y$ is proportional to the inverse of the specific vorticity, or vortensity, under the approximation of $\Omega$ being strictly the Keplerian angular velocity everywhere. Equation (10) then shows that the vortensity changes with time due to diffusion (first term on the right-hand side) and due to advection in the horseshoe region (second term on the right-hand side).

From now on, we will look for stationary solutions. Changing variables for a second time, with $z = 2(1+x)^{1/2} - 2$ and $f = (1+z/2)^{2\alpha-3}y(z)$, we find

$$ \frac{3v_p}{r_p^2} \frac{\partial^2 f(z)}{\partial z^2} + \frac{3\Omega_p}{4\pi} |z| \left[ \left( 1 - \left(\frac{3}{2} - \alpha\right)z \right) f(-z) - \left( 1 + \left(\frac{3}{2} - \alpha\right)z \right) f(z) \right] \Pi \left( \frac{z}{2z_s} \right) = 0, \quad (11) $$

where we have approximated the second term of equation (10) for $x \ll 1$ (or, equivalently, $z \ll 1$). Here, $z_s = z(x_s)$. Following Masset (2001), we split $f$ into an even part, $\tilde{f}$, and an odd part, $\tilde{f}$, given by

$$ \tilde{f}(z) = \frac{1}{2} (f(z) + f(-z)) \quad (12) $$

and

$$ \tilde{f}(z) = \frac{1}{2} (f(z) - f(-z)). \quad (13) $$

We then finally obtain

$$ \frac{\partial^2 f(z)}{\partial z^2} - \frac{r_p^2 \Omega_p}{2\pi v_p} |z| \left( \tilde{f}(z) + \left(\frac{3}{2} - \alpha\right) z \tilde{f}(z) \right) \Pi \left( \frac{z}{2z_s} \right) = 0. \quad (14) $$

The even part of equation (14) gives $\tilde{f} = f_p$, if we demand that $\tilde{f}$ is bounded for all $z$. Note that $f = f_p$ exactly corresponds to the unperturbed profile of $y(x)$. Defining

$$ k = \frac{r_p^2 \Omega_p}{2\pi v_p}, \quad (15) $$

the odd part of equation (14) gives

$$ \frac{\partial^2 \tilde{f}(z)}{\partial z^2} - kz \left( \tilde{f}(z) + \left(\frac{3}{2} - \alpha\right) z \tilde{f}(z) \right) = 0, \quad (16) $$

for $0 \le z \le z_s$, subject to boundary conditions at $z = 0$ and $z = z_s$.

Since $\tilde{f}$ is odd, we must have $\tilde{f}(0) = 0$. The most obvious outer boundary condition, to ensure continuity of $\tilde{f}$, is to take $\tilde{f}(z_s) = 0$.

This way, vortensity perturbations are localized to the horseshoe region only. However, for $x_s \ll h$, we expect pressure effects to spread vortensity perturbations over a length-scale $H$ (Casoli & Masset 2009). In other words, assuming Keplerian rotation everywhere means that vortensity perturbations will not be in pressure equilibrium, while for $x_s \ll h$ we expect pressure equilibrium to be maintained at all times.

We can incorporate this in the simple model by demanding that

$$ \Omega^2 = \Omega_K^2 + \frac{c_s^2}{r\Sigma} \frac{\partial \Sigma}{\partial r}. \quad (17) $$

However, using equation (17) in equation (8) makes the problem very difficult to solve. We can still make progress if we assume that Keplerian rotation still holds for $x < x_s$, while we use equation (17) in equation (8) for $x > x_s$. In order to match the two solutions at $x = x_s$, we solve equation (16) with boundary conditions $f(0) = 0$, $f(z_s) = f_s$, which give us the inner solution

$$ \tilde{f}_{\text{in}}(z) = \left(\frac{3}{2} - \alpha\right) f_p \left( \sqrt{z z_s} \frac{I_{1/3}\left(\frac{2}{3}\sqrt{k z^3}\right)}{I_{1/3}\left(\frac{2}{3}\sqrt{k z_s^3}\right)} - z \right) + f_s \sqrt{\frac{z}{z_s}} \frac{I_{1/3}\left(\frac{2}{3}\sqrt{k z^3}\right)}{I_{1/3}\left(\frac{2}{3}\sqrt{k z_s^3}\right)}, \quad (18) $$

where $I$ denotes a modified Bessel function. We fix $f_s$ by demanding that the inner solution matches in a smooth way to the outer solution. From the above equation, it is clear that saturation is governed by a single parameter $p$:

$$ p = 2\sqrt{kx_s^3/3}. \quad (19) $$

The outer solution $f_{\text{out}}$ will decay on a scale $H$ towards the equilibrium solution (Casoli & Masset 2009). There is no need to calculate it exactly, as long as we have $x_s \ll h$. We know that the odd part of the derivative of $f_{\text{out}}$ will approach zero on a length-scale $H$, so for $x_s \ll h$ we can set it to zero. We can therefore fix $f_s$ by demanding that the derivative of the inner solution is zero at $x_s$. Note that this is the same boundary condition as used in Masset (2001). This leads to

$$ f_s = -\left(\frac{3}{2} - \alpha\right) z_s f_p \frac{3p I_{4/3}(p)}{2I_{1/3}(p) + 3p I_{4/3}(p)}. \quad (20) $$

The resulting solution for $y(x)$ is shown in Fig. 1, for three different values of $p$, where we have chosen the outer profile for $x > x_s$ to approach the equilibrium profile exponentially on a scale $H$. The exact form does not matter since it does not affect the torque. Note that all profiles are indeed smooth at $x = x_s$.

The horseshoe drag in an isothermal disc is given by

$$ \Gamma_{\text{hs,baro}} = 3 \Sigma_p r_p^4 \Omega_p^2 \int_0^{\lambda_s} \frac{\bar{y}(x)}{y_p} x^2 dx, \quad (21) $$

where $\bar{y}(x)$ is the odd part of $y(x)$.

The torque on the planet can be obtained by again integrating over all material that has undergone a horseshoe turn. This is why the exact form of the outer solution is of no importance: since that material never executes a horseshoe turn, it does not contribute to the torque. We then get a multiplier for the total horseshoe drag:

$$ F(p) = \frac{8 I_{4/3}(p)}{3p I_{1/3}(p) + \frac{9}{2} p^2 I_{4/3}(p)}. \quad (22) $$

Equation (22) is shown in Fig. 2, together with numerical results obtained with FARGO. The model matches the numerical results to within 20 per cent for $p > 1$. For smaller values of $p$, or larger
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The function $y(x) = \Sigma(x)/\Omega(x)$, for three different saturation parameters $p$. The inset shows a larger spatial domain, where the exponential decay of the perturbation is apparent.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The function $F(p)$. The solid curve represents equation (22). Symbols indicate the results of numerical simulations using FARGO of a $q = 1.26 \times 10^{-5}$ planet embedded in a $h = 0.05$ disc. The standard model has $\alpha = 1/2$, $\beta = 1$ and $b/h = 0.4$.</figcaption></figure>

values of $v_p$, viscosity starts to affect the horseshoe turn itself, a regime we will deal with in Section 4.2. It can be shown that equation (22) is identical to equation (20) in Masset (2002). Note, however, that in our case we automatically get the correct density profile dependence, which removes the ambiguity in the choice of $F$ in Masset (2001). We comment that $F(p)$ may be approximated to within 5 per cent by a simpler form

$$F(p) = \frac{1}{1 + (p/1.3)^2} \qquad (23)$$

## 4.2 The limit of high $v_p$

In the limit $p \ll 1$, the assumption that the horseshoe turn over time is much smaller than the viscous diffusion time-scale across the horseshoe region breaks down. In this regime, the horseshoe turn itself is modified by viscous effects, and specific vorticity is no longer conserved during the turn. It was shown in Paardekooper & Papaloizou (2009a) that the non-linear horseshoe drag is then replaced by the linear corotation torque, which is in general much weaker. Because of the strong viscosity required, saturation does not play a role in this regime, making the torque fully linear. It is important to understand this regime, since in the non-barotropic case it corresponds to the transition between the optically thick, mostly adiabatic, inner disc and the optically thin, locally isothermal, outer disc, with viscosity replaced by thermal diffusion.

The regime in which the viscosity is strong enough to affect the horseshoe turns, but not strong enough to make the torque fully linear, is very difficult to model. All we can say is that in the limit $p \to 0$, there are basically no turns as far as the flow is concerned, since the viscosity is strong enough for the flow to forget where it came from, and the corotation torque will be fully linear. In the opposite limit, the horseshoe turns will prevent the linear torque from being set up, and the corotation torque will be given by the (possibly saturated) horseshoe drag. As the viscosity increases, the horseshoe drag decreases and the linear corotation torque increases, and it is not a priori clear that these processes occur at the same rate. In fact, one may expect the rates to be different, since the time-scale to set up the horseshoe drag, a fraction of the libration time-scale, is longer than the time-scale associated with the linear corotation torque, which is a dynamical time-scale.

We therefore make the assumption that the corotation torque is given by

$$\Gamma_{c,baro} = G\Gamma_{hs,baro} + (1-K)\Gamma_{c,lin,baro}, \qquad (24)$$

with $0 \le G, K \le 1$. In the limit of low viscosity (but still high enough to keep the corotation torque unsaturated), we expect $G, K \to 1$, and the corotation torque is given by the horseshoe drag. In the limit of high viscosity, $G, K \to 0$, and the torque will be fully linear. As mentioned above, we allow the horseshoe drag to decrease at a different rate than the linear torque increases towards high viscosities by having $K$ different from $G$. Below, we will derive a general form for $G$, from which we will extract $K$ using a time-scale argument.

Consider the contribution of a single streamline, making a horseshoe turn from $-x$ to $x$, to the horseshoe drag. It was shown in Paardekooper & Papaloizou (2009a) that the time-scale to set up the horseshoe drag is approximately $\tau_{turn} = 3\tau_{lib}/20$. When this is comparable to the diffusion time-scale across $x$, viscosity is able to affect the horseshoe drag. We take this into account by taking the contribution of the streamline to the horseshoe drag to be

$$d\Gamma_{hs,baro} = 3 \left(\frac{3}{2} - \alpha\right) \Sigma_p r_p^4 \Omega_p^2 x^3 \mathcal{F} \left(\tau \frac{x^3}{x_s^3}\right) dx, \qquad (25)$$

with $\tau = \tau_{diff}/\tau_{turn} = 45\pi p^2/4$ [using $\tau_{turn} = 2/(5x_s\Omega_p)$ and $\tau_{diff} = x_s^2 r_p^2/\nu_p$] and $\mathcal{F}(t)$ such that $\mathcal{F}(0) = 0$ and $\mathcal{F}(t) \to 1$ for $t \to \infty$. For $\tau_{turn} \ll \tau_{diff}$, this gives us the ordinary horseshoe drag, while for $\tau_{turn} \gg \tau_{diff}$, the contribution of the streamline goes to zero.

To keep all integrals tractable, we choose a particularly simple form for $\mathcal{F}$:

$$\mathcal{F}(t) = \begin{cases} (t/\tau_0)^l & t < \tau_0 \\ 1 & t \ge \tau_0, \end{cases} \qquad (26)$$
for some $l, \tau_0 > 0$. We have obtained good agreement with numerical experiments with $l = 3/4$ and $\tau_0 = 2$. It is straightforward to show that in general

$$G(\tau) = 4 \int_{0}^{1} w^3 \mathcal{F}(\tau w^3) dw = \frac{4}{3} \tau^{-4/3} \int_{0}^{\tau} t^{1/3} \mathcal{F}(t) dt, \quad (27)$$

which for our particular form for $\mathcal{F}$ yields

$$G(\tau) = \begin{cases} \frac{4}{4+3l} \left(\frac{\tau}{\tau_0}\right)^l & \tau < \tau_0 \\ 1 - \frac{3l}{4+3l} \left(\frac{\tau}{\tau_0}\right)^{-4/3} & \tau \ge \tau_0. \end{cases} \quad (28)$$

This can be written in terms of $p$ rather than $\tau$:

$$G(p) = \begin{cases} \frac{4}{4+3l} \left( \frac{45\pi}{4\tau_0} \right)^l p^{2l} & p < \sqrt{\frac{4\tau_0}{45\pi}} \\ 1 - \frac{3l}{4+3l} \left( \frac{4\tau_0}{45\pi} \right)^{4/3} p^{-8/3} & p \ge \sqrt{\frac{4\tau_0}{45\pi}}. \end{cases} \quad (29)$$

For $l = 3/4, \tau_0 = 2$, we finally have

$$G(p) = \begin{cases} \frac{16}{25} \left( \frac{45\pi}{8} \right)^{3/4} p^{3/2} & p < \sqrt{\frac{8}{45\pi}} \\ 1 - \frac{9}{25} \left( \frac{28}{45\pi} \right)^{4/3} p^{-8/3} & p \ge \sqrt{\frac{8}{45\pi}}. \end{cases} \quad (30)$$

We now assume that $K$ is similar to $G$, but with a different value of $\tau_0$. Good results were obtained with $\tau_0 = 7$, so that the linear torque decreases slower towards lower viscosities than the horseshoe drag increases. This leads to a maximum corotation torque that is in fact larger than the unsaturated horseshoe drag, a feature that is already apparent in the numerical results shown in Fig. 2 (see also Baruteau & Lin 2010). For this different value of $\tau_0$, we find

$$K(p) = \begin{cases} \frac{16}{25} \left( \frac{45\pi}{28} \right)^{3/4} p^{3/2} & p < \sqrt{\frac{28}{45\pi}} \\ 1 - \frac{9}{25} \left( \frac{28}{45\pi} \right)^{4/3} p^{-8/3} & p \ge \sqrt{\frac{28}{45\pi}}. \end{cases} \quad (31)$$

We can now construct the complete model for the corotation torque, including both saturation and the cut-off at high viscosity:

$$\Gamma_{c,baro} = G(p)F(p)\Gamma_{hs,baro} + (1 - K(p))\Gamma_{c,lin,baro}, \quad (32)$$

which uses equations (22), (30) and (31) and the torques from equations (4) and (6). We compare the complete model to numerical results in Fig. 3, where we have used the expressions for the horseshoe drag and linear corotation torque from Paper I. As in the isothermal case, for the numerical results, $x_s$ was determined using a dichotomic search of the separatrices, with $x_s$ being calculated as the geometric average of the horseshoe half-widths at $\pm 1$ rad (Casoli & Masset 2009). Equation (32) captures the main features of the numerical results quite well. We have checked that this is true for different density and temperature profiles.

# 5 NON-ISOTHERMAL DISCS

## 5.1 Disc models

We now consider the behaviour of non-barotropic effects under the influence of diffusion. The disc models are similar to those in the isothermal case, but in addition we solve the energy equation. We stress again that we do not include viscous energy dissipation in the equations; we assume that dissipation, together with radiative cooling, sets the background temperature profile.

## 5.2 Thermal diffusion

Since we now have an extra physical quantity that is evolved, namely the temperature, we consider thermal diffusion in addition to viscous

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Corotation torque as a function of viscosity (through the $p$ parameter). Symbols indicate numerical results, obtained with FARGO, for a $q = 1.26 \times 10^{-5}$ planet in a locally isothermal disc with $\alpha = 1/2$, $\beta = 1$, $h = 0.05$ and $b/h = 0.4$. The solid curve denotes the model of equation (32).</figcaption></figure>

diffusion. Since we want the disc to be in thermal equilibrium, i.e. viscous heating should compensate for radiative losses, we can constrain the range of values we need to consider for the thermal diffusion coefficient.

Radiative diffusion gives rise to a thermal conductivity

$$K_R = \frac{acT^3}{3\kappa\rho} = \frac{4\sigma T^3}{3\kappa\rho}, \quad (33)$$

where $a$ is the radiation constant, $c$ the speed of light, $\sigma$ the Stefan-Boltzmann constant and $\kappa$ the opacity. The corresponding thermal diffusion coefficient $\chi$ is

$$\chi = \frac{K_R}{c_v \rho} = \frac{4(\gamma - 1)\sigma T^4}{3\kappa\rho P} = \frac{4\gamma(\gamma - 1)\sigma T^4}{3\kappa\rho^2 H^2 \Omega^2}, \quad (34)$$

where we have used the internal energy $\epsilon = c_v \rho T = P/(\gamma - 1)$, with $c_v$ being the specific heat at constant volume.

Thermal equilibrium requires that (Kley & Crida 2008)

$$\sigma T_{\text{eff}}^4 = \frac{9}{4} \nu \Sigma \Omega^2, \quad (35)$$

where $T_{\text{eff}}$ is the effective temperature (Hubeny 1990):

$$T_{\text{eff}}^4 = T^4 / \tau_{\text{eff}}, \quad (36)$$

with an effective optical depth

$$\tau_{\text{eff}} = \frac{3\tau}{8} + \frac{\sqrt{3}}{4} + \frac{1}{4\tau}, \quad (37)$$

with $\tau = \kappa \rho H$ being the vertical optical depth.

Using equation (35) in equation (34) leads to

$$\frac{\chi}{\nu} = Pr^{-1} = \frac{9}{4}\gamma(\gamma - 1)\left(1 + \frac{2\sqrt{3}}{3\tau} + \frac{2}{3\tau^2}\right), \quad (38)$$

where $Pr$ is the Prandtl number, and we have taken $\Sigma = 2\rho H$. We therefore have that in the limit of $\tau \to \infty$, $Pr$ will be of order unity, with in general $Pr < 1$. In other words, thermal and viscous diffusion will be of equal magnitude in the optically thick inner regions of protoplanetary discs. In the outer regions, thermal diffusion will
come to dominate, leading eventually to isothermal behaviour for $\tau \ll 1$ (Paardekooper & Mellema 2006b). Therefore, we only have to consider cases for which $\chi \ge \nu$ or, equivalently, $Pr \le 1$. Note that this also holds in the presence of other heating sources such as external irradiation: in thermal equilibrium, any non-viscous heating source will require an increase in $\chi$ to radiate this energy away. Additional sources of heating therefore necessarily decrease $Pr$ even further.

We first consider a barotropic disc, including thermal diffusion, and show that thermal diffusion does not affect the barotropic part of the corotation torque. Next, we study the case $v_p = \chi_p$, with $\chi_p$ being the thermal diffusion coefficient at the location of the planet, and show how this case is different to the barotropic case. We then adjust the model to incorporate the case $v_p \le \chi_p$.

## 5.3 Constant entropy disc

A disc with $P = P(\Sigma) \propto \Sigma^\gamma$ (or constant entropy) initially behaves in a similar way to an isothermal disc. Since the entropy is constant, there is no entropy-related horseshoe drag, and since vortensity is conserved along streamlines in barotropic flows, the vortensity-related horseshoe drag is exactly the same as in the isothermal case. The only difference between the isothermal and constant entropy cases is the value of $\gamma$, which affects not only the total Lindblad and corotation torques (see Paper I), but also the saturation properties through its effect on $x_s$. Since $x_s$ decreases for increasing $\gamma$, for a given value of $v_p$, the horseshoe drag in a constant entropy disc with $\gamma > 1$ will appear less saturated compared to the isothermal case.

We show in Fig. 4 that the time evolution of the torque in a constant entropy disc is not affected by thermal diffusion. Since we do not apply viscous heating or radiative cooling, entropy is conserved for $\chi_p = 0$, keeping the disc, and the associated torque on the planet, fully barotropic. When applying thermal diffusion, entropy is no longer strictly conserved, but since both thermal and viscous diffusion act to restore the original density and temperature profiles, the disc essentially remains barotropic. The torque on the

<figure><img src="image_6.png" /><figcaption>Figure 4. Total torque on a $q = 1.26 \times 10^{-5}$ planet, $b/h = 0.4$, embedded in a $h = 0.05$ disc with $\alpha = -1/2$, $\beta = -1/5$ and $\gamma = 7/5$, making $\xi = 0$. Two values for thermal diffusion are considered: $\chi_p = 0$ (solid line) and $\chi_p = 10^{-6} r_p^2 \Omega_p$ (dotted line). Results were obtained with RODEO.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Corotation torque as a function of viscosity (through the $p$ parameter). Symbols indicate numerical results for a $q = 1.26 \times 10^{-5}$ planet in a disc with $\alpha = 1/2$, $\beta = 1$, $h = 0.05$, $\gamma = 7/5$ and $b/h = 0.4$. The Prandtl number equals unity in all runs. The solid curve denotes the model of equation (32), with the appropriate values of the barotropic and non-barotropic linear and non-linear unsaturated torques, and the dashed curve denotes the modified model discussed in Section 5.4 (see also equation 51).</figcaption></figure>

planet is therefore very similar in both cases. This is not true for very large values of $\chi_p$: when thermal diffusion is strong enough to affect the horseshoe turn itself, we approach the isothermal regime. This limit is considered in Section 5.5.

## 5.4 The case $v_p = \chi_p$

We now consider the case $Pr = 1$. We have used the same disc model as for Fig. 3, but now solving the energy equation with $\gamma = 7/5$. Note that this disc model has negative gradients in both entropy and vortensity, making the total horseshoe drag positive, and dominated by the entropy-related part (see Paper I). The results are displayed in Fig. 5. The solid curve denotes the model of equation (32), with the horseshoe drag and linear corotation torque from Paper I. It is clear that the torque saturates faster than predicted by equation (32). This has to do with the fact that the entropy-related part of the horseshoe drag essentially originates from a single streamline close to the separatrix. At this location, advection is at its maximum efficiency, since the velocity with respect to the planet is largest, and diffusion has a hard time to keep up. At the same time, the entropy-related part needs both thermal and viscous diffusion to remain unsaturated. Thermal diffusion alone cannot provide fresh angular momentum to the horseshoe region. Even if thermal diffusion (or entropy diffusion) can restore the entropy gradient sufficiently, if viscosity is not strong enough to diffuse the resulting vortensity production at the separatrices, the torque will still saturate (Paardekooper & Palaoizou 2008). In this sense, the entropy torque can be doubly saturated.

Rather than trying to incorporate these effects into the simple isothermal model presented in Section 4.1, we try to adjust the model, guided by the numerical results to allow for the effects described above. Both can be incorporated by making the substitution

$$F(p) \rightarrow F(p_v)F(p_\chi), \qquad (39)$$
where $p_\chi \propto p_v \sqrt{v_p/\chi_p}$ is the saturation parameter related to thermal diffusion and $p_v$ is the saturation parameter associated with viscosity, in equation (32). Since advection is most efficient in the region where the entropy-related torque is generated, we expect $p_\chi > p_v$ for $v_p = \chi_p$. Good results were obtained with $p_\chi = 3p_v\sqrt{v_p/\chi_p}/2$, or

$$p_{\chi} = \sqrt{\frac{r_{\mathrm{p}}^{2} \Omega_{\mathrm{p}} \chi_{\mathrm{s}}^{3}}{2\pi \chi_{\mathrm{p}}}} \quad (40)$$

We can see indeed that the dashed curve in Fig. 5, obtained with equation (39), gives a good fit to the data at large $p_v$. Note that substitution (39) needs to be made for the entropy-related part of the horseshoe drag only.

The cut-off at high $v$ and $\chi$ should also be modified. Note that in this regime, advection being most efficient actually helps to keep the torque non-linear at higher levels of diffusion. Rather than defining new functions $G$ and $K$, we have found that the numerical results can be fitted by making the substitutions

$$G(p) \rightarrow \sqrt{G(p_v)G(p_\chi)} \quad (41)$$

and

$$1 - K(p) \rightarrow \sqrt{(1 - K(p_v))(1 - K(p_\chi))} \quad (42)$$

Here again, these substitutions only concern the entropy-related part of the horseshoe drag. The resulting model is shown by the dashed curve in Fig. 5. Note that at higher levels of diffusion, the RODEO results differ by 20 per cent from results obtained with FARGO. We have verified that half of this is due to the different diffusion implementation: entropy diffusion implemented in RODEO gave lower values for the torque by approximately 10 per cent at high diffusion. The remaining 10 per cent can be attributed to a higher value of the unsaturated torque using RODEO (see Paper I). This dependence on the nature of the diffusion should be kept in mind when interpreting these results in terms of radiative discs, where $\chi$ is due to radiative diffusion. The error in the model compared to the simulations can thus be up to 20 per cent.

## 5.5 The general case $v_p \le \chi_p$

We now consider the general case, still with the restriction $Pr \le 1$. We have run a grid of models with different $v_p$ and $\chi_p$, and the results are displayed in Fig. 6. All numerical results were obtained with RODEO, using thermal diffusion.

For all values of $p_v$, $Pr = 1$ gives the lowest value for the corotation torque. Going up in thermal diffusion, at a fixed value of $p_v$, the torque increases towards the original model, denoted by the solid curve. This is consistent with the effect discussed in the previous section, that a higher thermal diffusion coefficient is needed to fight advection of entropy along the separatrix. For all values of $\chi_p$, the improved model gives a reasonable fit to the numerical results, with all deviations being less than approximately 20 per cent.

For $\chi_p > 10^{-5}$, the corotation torque starts to decrease again. This is due to thermal diffusion starting to affect the horseshoe turns. For $\chi_p$ large enough, we expect to enter the locally isothermal regime, where the corotation torque is given by the barotropic, non-linear horseshoe drag plus the linear, entropy-related corotation torque. Since the linear corotation torque is in general smaller than its non-linear counterpart, the corotation torque will decrease towards higher values of $\chi_p$. This is illustrated in Fig. 7, where we show the time evolution of the corotation torque at a fixed value of $v_p =$

$10^{-5}$ for different values of $\chi_p$. The torque decreases with $\chi_p$, but the barotropic part of the corotation torque remains non-linear. The barotropic horseshoe drag plus the linear entropy-related corotation torque amount to $1.4\Gamma_0$. However, at these high values of $\chi_p$, we also need to consider the effect of thermal diffusion on the effective value of $\gamma$, which will approach unity for large values of $\chi_p$. This we discuss next.

## 5.6 Lindblad torque

The Lindblad torque is due to a superposition of linear waves and does not strongly depend on viscosity, for reasonable values of the viscosity $v$. Thermal diffusion, on the other hand, can play an important role and is necessary to recover the isothermal result for large values of the thermal diffusivity $\chi_p$. The difference between the isothermal and adiabatic Lindblad torque basically results from a difference in sound speeds. Since $\Gamma_L \propto c_s^{-2}$, $\Gamma_{L,iso}/\Gamma_{L,adi} = \gamma$.

Consider the dispersion relation for plane waves $\propto \exp[i(\omega t - kx)]$ (Mihalas & Mihalas 1984):

$$\omega^2 = c_s^2 k^2 \frac{1 - i\chi k^2 / \omega}{1 - i\gamma \chi k^2 / \omega}, \quad (43)$$

with $c_s$ being the adiabatic sound speed. From this, one can derive a phase speed $v_p$:

$$\frac{v_p^2}{c_s^2} = \frac{2Q}{\gamma Q + \frac{1}{2}\sqrt{2\sqrt{(\gamma^2 Q^2 - 1)^2 + (4Q - 2\gamma Q)^2 + 2\gamma^2 Q^2 - 2}}} \quad (44)$$

with $Q = \chi\omega/c_s^2$. In the limit of no thermal diffusion ($Q = 0$), $v_p \rightarrow c_s$: waves propagate at the adiabatic sound speed. For $Q \rightarrow \infty$, $v_p \rightarrow c_s/\gamma$: waves propagate at the isothermal sound speed.

The torque for any Fourier component will now scale as $v_p^{-2}$, which makes the isothermal torque a factor of $\gamma$ stronger than the adiabatic torque. In principle, $\omega = m\Omega_p$ for the problem at hand, and one should therefore consider the effects of diffusion on each Fourier component separately. For simplicity, we consider a single, representative value of $m$ to describe the effects of diffusion. Since most of the torque comes from $m \sim 2/3h$, we set

$$Q = \frac{2\chi_p \Omega_p}{3hc_s^2} = \frac{2\chi_p}{3h^3 r_p^2 \Omega_p} \quad (45)$$

It is convenient to define an effective $\gamma$:

$$\gamma_{\text{eff}} = \frac{2Q\gamma}{\gamma Q + \frac{1}{2}\sqrt{2\sqrt{(\gamma^2 Q^2 + 1)^2 - 16Q^2(\gamma - 1) + 2\gamma^2 Q^2 - 2}}} \quad (46)$$

so that

$$\gamma_{\text{eff}}\Gamma_L / \Gamma_0 = -2.5 - 1.7\beta + 0.1\alpha, \quad (47)$$

which then replaces equation (3). Note that the horseshoe drag also scales with $\gamma_{\text{eff}}$, since the horseshoe width depends on $\gamma$ (Paardekooper & Papalozou 2009b). Equation (46) is plotted in Fig. 8 as a function of $\alpha_\chi = \chi_p/(c_sH)$, for three different scale-heights $h$. It is clear that a substantial thermal diffusion coefficient is required in order to modify the Lindblad torque (typically $\chi_p = 10^{-3}r_p^2\Omega_p$). We have overplotted numerical results obtained with RODEO, for a disc with $h = 0.05$, $\alpha = 3/2$, $\beta = 1$ and $\gamma = 5/3$. Since this disc has no entropy or vortensity gradients, the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Corotation torque as a function of viscosity (through the *p* parameter). Symbols indicate numerical results for a *q* = 1.26 × 10⁻⁵ planet in a disc with *α* = 1/2, *β* = 1, *h* = 0.05, *γ* = 7/5 and *b/h* = 0.4, for different values of thermal diffusion *χ*<sub>p</sub>. The solid curve denotes the model of equation (32) (saturation of the non-barotropic horseshoe drag is similar to the barotropic part), with the appropriate values for the barotropic and non-barotropic parts of the linear and non-linear unsaturated torques. The remaining curves indicate the modified model discussed in Section 5.4 for different values of *χ*<sub>p</sub> (see also equation 51). Different panels highlight different values of *χ*<sub>p</sub>, with the corresponding model prediction. Numerical results were obtained with RODEO.</figcaption></figure>

corotation torque vanishes. We can therefore measure $\gamma_{\text{eff}}$ from the
total torque, which should be given by equation (3) with $\gamma = \gamma_{\text{eff}}$.
Since for equation (46) we have chosen a single value of $m$ to signal
the transition from isothermal to adiabatic, we expect the numeri-
cal results to show a smoother transition between $\gamma_{\text{eff}} = 5/3$ and
$\gamma_{\text{eff}} = 1$. This is indeed what is observed in Fig. 8. Note, however,
that the model does a good job in predicting where the transition
takes place. We have checked that this also holds for different values
of $h$. Around this transition, the maximum error in $\gamma_{\text{eff}}$ observed is
20 per cent. Note that, since all unsaturated torques are proportional
to $\gamma_{\text{eff}}^{-1}$ (see below), the error in $\gamma_{\text{eff}}$ only affects the magnitude of
the total torque, not the sign. The dependence of $\gamma_{\text{eff}}$ on thermal
diffusion will therefore not affect planet migration in a qualitative
way, but it is necessary to connect the adiabatic to the isothermal
regime.

5.7 A modified torque formula

Before we combine the results of the previous sections into a single
torque formula, we restate the assumptions made in deriving these
results. We have worked in the two-dimensional approximation,
considered laminar discs only and kept the planet on a fixed circular
orbit. We have ignored the self-gravity of the disc, assumed thermal
equilibrium and considered a simplified energy equation. It is clear
that future studies should be aimed at clarifying the role of magnetic
turbulence and radiative energy transport in keeping the corotation
torque unsaturated.

We have mainly focused on $q = 1.26 \times 10^{-5}$ using a softening
parameter $b/h = 0.4$, which seems to be an appropriate value to
mimic 3D effects (Paardekooper & Papaloizou 2009a). In the model,
the only place where the mass and the softening parameter enter is
<figure><img src="image_9.png" /><figcaption>Figure 7. Corotation torque on a $q = 1.26 \times 10^{-5}$ planet in a disc with $\alpha = 1/2$, $\beta = 1$, $h = 0.05$, $\gamma = 7/5$ and $b/h = 0.4$, $v_p = 10^{-5}$ for different values of thermal diffusion $\chi_p$. Results were obtained with RODEO.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. The values of $\gamma_{eff}$ from equation (46), as a function of thermal diffusion, here represented by $\alpha_{\chi} = \chi_p/(c_s\Omega_p)$. Curves denote equation (46) for different values of $h$; open circles indicate numerical results obtained with RODEO for $h = 0.05$.</figcaption></figure>

in the values of the unsaturated torques, which have a $q^2$ dependence and depend on $b/h$ (Paardekooper & Papaloizou 2009a), and in the value of $x_s$, which can be written as (Paardekooper & Papaloizou 2009b)

$$x_s = C(b/h) \sqrt{\frac{q}{h}}, \qquad (48)$$

where C can be approximated by a power law around $b/h = 0.4$:

$$C = \frac{1.1}{\gamma_{\text{eff}}^{1/4}} \left( \frac{0.4}{b/h} \right)^{-1/4} \qquad (49)$$

(see Paper I). The dependence on softening of the linear torques was studied in Paardekooper & Papaloizou (2009a) and of $x_s$ and the non-barotropic horseshoe drag in Paper I. Given a different value of the softening, it is straightforward to adjust the formulae below by using equations (48) and (49). The value of $x_s$ then features in the saturation parameters $p_v$ (given by equation 19) and $p_{\chi}$ (equation 40). Note, in particular, that the functional form of F, G and K will not change for different values of $b/h$.

The total torque on a low-mass planet is given by the Lindblad torque plus the corotation torque:

$$\Gamma = \Gamma_L + \Gamma_c, \qquad (50)$$

with $\Gamma_L$ given by equation (3), but with $\gamma \to \gamma_{eff}$ (see equation 47).

The corotation torque is split into a barotropic part and an entropy-related part:

$$\Gamma_c = \Gamma_{c,\text{baro}} + \Gamma_{c,\text{ent}}. \qquad (51)$$

The barotropic part of the horseshoe drag is not affected by thermal diffusion; therefore,

$$\Gamma_{c,\text{baro}} = \Gamma_{\text{hs,baro}} F(p_v)G(p_v) + (1 - K(p_v))\Gamma_{c,\text{lin,baro}}, \qquad (52)$$

where $\Gamma_{hs,baro}$ is the fully unsaturated horseshoe drag (equation 4, but with $\gamma \to \gamma_{eff}$), $\Gamma_{c,lin,baro}$ is the linear barotropic corotation torque (equation 6, but with $\gamma \to \gamma_{eff}$), $F(p)$ governs saturation (see equation 22) and $G(p)$ and $K(p)$ govern the cut-off at high viscosity (see equations 30 and 31).

For the non-barotropic, entropy-related corotation torque, we need the modifications discussed in Section 5.4:

$$\Gamma_{c,\text{ent}} = \Gamma_{\text{hs,ent}} F(p_v) F(p_{\chi}) \sqrt{G(p_v)G(p_{\chi})} + \sqrt{(1 - K(p_v))(1 - K(p_{\chi}))} \Gamma_{c,\text{lin,ent}}, \qquad (53)$$

where $\Gamma_{hs,ent}$ is the fully unsaturated, entropy-related part of the horseshoe drag (equation 5, but with $\gamma \to \gamma_{eff}$), $\Gamma_{c,lin,ent}$ is the linear, entropy-related part of the corotation torque (equation 7, but with $\gamma \to \gamma_{eff}$) and $p_{\chi}$ is the saturation parameter associated with thermal diffusion (see equation 40).

From Fig. 6, we estimate that the maximum error in the corotation torque is approximately 20 per cent compared to the simulations. From Paper I, we know that the error in the Lindblad torque is smaller in general, so for the total torque we expect differences of 20 per cent between model and simulations as well, except close to the transition between isothermal and adiabatic behaviour, where the error in the Lindblad torque can be 20 per cent as well (see Fig. 8).

# 6 DISCUSSION

We have presented numerical results and some simple models that capture the behaviour of the torque due to Lindblad resonances and horseshoe drag on low-mass planets embedded in gaseous discs in the presence of viscous and thermal diffusion. Several simplifications of the real physical system were made to keep the problem tractable.

First, we have worked with vertical integrated quantities only, which make the problem two-dimensional. The entropy-related horseshoe drag seems to play an important role in three-dimensional discs as well (Kley et al. 2009), and it has been shown that, in an isothermal set-up, the effect of the horseshoe drag is stronger in three dimensions (Masset, D’Angelo & Kley 2006b). Clearly, a three-dimensional model of the horseshoe region is desirable. This would also get rid of the dependence of the torque on the softening parameter that is necessarily part of a two-dimensional model (see
Paper I). The results of Kley et al. (2009) indicate, however, that the main features of three-dimensional simulations can be reproduced nicely within the two-dimensional approximation.

We have only considered laminar, viscous discs. In reality, discs are expected to be turbulent in regions where the gas couples to the magnetic field of the disc (Balbus & Hawley 1991). The interaction of turbulence with horseshoe turns is an area that is still largely unexplored. For the isothermal case, with a simplified turbulence model based on stochastic forcing, Baruteau & Lin (2010) showed that turbulent models show similar behaviour to laminar models. It remains to be seen how full magnetic turbulence, in a non-isothermal setting, affects the horseshoe drag.

We have modelled radiative cooling as thermal diffusion, with a constant thermal diffusion coefficient. In radiative models, $\chi$ will depend strongly on temperature and opacity. Moreover, radiative diffusion will probably be most efficient in the vertical direction, while in our two-dimensional approximation thermal diffusion is restricted to the plane of the orbit of the planet. It was noted in Kley & Crida (2008) that this is in fact an unfavourable case for sustaining the corotation torque compared to vertical cooling in radiative models. Note that we have chosen $\chi$ such that the thermal diffusion time-scale is identical to the vertical cooling time-scale. Since both act to restore the original temperature profile, it is to be expected that including heating and cooling (without thermal diffusion) will lead to similar effects to thermal diffusion alone. When acting together, they may reduce the effective Prandtl number since the flow then has two ways of restoring the original temperature profile.

We have taken the background gradients of temperature and surface density to be given. In reality, they will be determined by a balance of radiative heating and cooling and viscous dissipation. It was recently shown in Lyra, Paardekooper & Mac Low (2010) that discs where these profiles are calculated self-consistently allow for outward Type I migration to occur for a large fraction of the lifetime of the disc. Of course, since the mass of the disc decreases with time, radiative cooling becomes more and more efficient, leading in the end to a situation where the entropy-related torques are fully linear.

The strong dependence of the horseshoe drag on the local temperature gradient makes it relatively easy to trap planets at special locations in the disc. While in the isothermal case, a surface density jump (Masset, Morbidelli, Crida & Ferreira 2006a) or at least a positive surface density gradient (Paardekooper & Papaloizou 2009a) is required, in the non-isothermal case a trap is readily established if a local negative temperature gradient exists. Even if the disc as a whole does not permit outward Type I migration, low-mass planets may stop their inward journey when reaching such a temperature trap.

If the disc parameters allow for outward migration, a natural planet trap occurs at the transition between the optically thick and optically thin regions of the disc. Outside this radius, we have isothermal, and therefore inward, Type I migration, while inside this radius, we have adiabatic, and therefore outward, Type I migration. The evolution of this equilibrium radius was explored in more detail in Lyra et al. (2010). Because of the dependence of the corotation torque on the saturation parameter $p$, which depends on the planet mass through $x_s$, different planet masses will have different equilibrium radii. The impact of this migration behaviour on the gravitational interactions between planets is still largely unexplored.

We have focused on low-mass planets, comparable to the Earth. While planets of approximately 1 Jupiter mass start to open up deep annular gaps in the disc (Lin & Papaloizou 1986), there exists an interesting intermediate regime where planets comparable to Neptune experience an enhanced horseshoe drag (Masset et al. 2006b). This

was interpreted in Paardekooper & Papaloizou (2009b) as being due to the ineffectiveness of the Lindblad wake to push against the horseshoe region, leading to a larger value of the horseshoe width at these masses. This means that these planets are more readily subject to outward migration, leading to mass segregation even without considering saturation.

The planet was kept on a fixed circular orbit throughout this paper, and we have ignored any radial viscous flow in the gas. If the relative radial speed of the planet with respect to the gas is large enough, the streamline topology will be modified, leading to a flow structure similar to that encountered in Type III migration (Masset & Papaloizou 2003; Pepliński, Artymowicz & Mellema 2008). Although in this case the discussion on the corotation torque will need to be modified, we comment that a very massive disc, with a Toomre $Q$ of order unity, is needed for a low-mass planet to migrate fast enough to achieve this situation.

Finally, we have ignored the self-gravity of the disc. The impact of self-gravity on Type I migration was discussed in Pierens & Huré (2005) and Baruteau & Masset (2008b), finding an acceleration of inward migration due to a shift in the Lindblad resonances. This should be taken into account especially in massive discs. A quantitative estimate of this effect can be inferred from fig. 8 of Baruteau & Masset (2008b).

## 7 CONCLUSIONS

We have presented a simple model that can capture the complicated behaviour of the corotation torque with respect to viscous and thermal diffusion. The formulae presented in Section 5.7 show good agreement with two-dimensional hydrodynamic simulations to within 20 per cent. The strong dependence of the total torque on temperature gradients and the level of viscous and thermal diffusion allow for a wealth of possible outcomes of Type I migration, depending on global disc parameters.

## ACKNOWLEDGMENTS

We thank the referee for an insightful report that helped us improve the clarity of the paper. It is a pleasure to thank A. Crida and K. Kretke for a careful reading of the manuscript. SJP acknowledges support from STFC in the form of a post-doctoral fellowship. Simulations using RODEO were performed using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk), provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England. FARGO simulations were performed using the Pleiades cluster at UC Santa Cruz.

## REFERENCES

* Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214

* Baruteau C., Lin D. N. C., 2010, ApJ, 709, 759

* Baruteau C., Masset F., 2008a, ApJ, 672, 1054

* Baruteau C., Masset F., 2008b, ApJ, 678, 483

* Casoli J., Masset F. S., 2009, ApJ, 703, 845

* Goldreich P., Tremaine S., 1979, ApJ, 233, 857

* Goldreich P., Tremaine S., 1980, ApJ, 241, 425

* Hubeny I., 1990, ApJ, 351, 632

* Kley W., Crida A., 2008, A&A, 487, L9

* Kley W., Bitsch B., Klahr H., 2009, A&A, 506, 971

* Lin D. N. C., Papaloizou J., 1986, ApJ, 307, 395

* Lyra W., Paardekooper S., Mac Low M., 2010, ApJ, 715, L68

* Masset F., 2000a, A&AS, 141, 165
Masset F. S., 2000b, in Garzón G., Eiroa C., de Winter D., Mahoney T. J., eds, ASP Conf. Ser. Vol. 219, Disks, Planetesimals, and Planets, FARGO: A Fast Eulerian Transport Algorithm for Differentially Rotating Disks. Astron. Soc. Pac., San Francisco, p. 75

Masset F. S., 2001, ApJ, 558, 453

Masset F. S., 2002, A&A, 387, 605

Masset F. S., Papaloizou J. C. B., 2003, ApJ, 588, 494

Masset F. S., Morbidelli A., Crida A., Ferreira J., 2006a, ApJ, 642, 478

Masset F. S., D'Angelo G., Kley W., 2006b, ApJ, 652, 730

Mihalas D., Mihalas B. W., 1984, Foundations of Radiation Hydrodynamics. Oxford Univ. Press, New York

Ogilvie G. I., Lubow S. H., 2003, ApJ, 587, 398

Paardekooper S.-J., Mellema G., 2006a, A&A, 450, 1203

Paardekooper S.-J., Mellema G., 2006b, A&A, 459, L17

Paardekooper S.-J., Mellema G., 2008, A&A, 478, 245

Paardekooper S.-J., Papaloizou J. C. B., 2008, A&A, 485, 877

Paardekooper S.-J., Papaloizou J. C. B., 2009a, MNRAS, 394, 2283

Paardekooper S.-J., Papaloizou J. C. B., 2009b, MNRAS, 394, 2297

Paardekooper S.-J., Baruteau C., Crida A., Kley W., 2010, MNRAS, 401, 1950 (Paper I)

Pepliński A., Artymowicz P., Mellema G., 2008, MNRAS, 386, 179

Pierens A., Huré J., 2005, A&A, 433, L37

Quinn P. J., Goodman J., 1986, ApJ, 309, 472

Roe P. L., 1981, J. Comput. Phys, 43, 357

Shakura N. I., Sunyaev R. A., Zilitinkevich S. S., 1978, A&A, 62, 179

Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257

Ward W. R., 1991, Lunar and Planetary Institute Conf. Abstracts, Horsehoe Orbit Drag. p. 1463

Ward W. R., 1997, Icarus, 126, 261

Ward W. R., 2007, Lunar and Planetary Institute Conf. Abstracts Vol. 38, A Streamline Model of Horseshoe Torque Saturation. p. 2289

This paper has been typeset from a *TEX/LATEX* file prepared by the author.