ON THE HORSESHOE DRAG OF A LOW-MASS PLANET. II. MIGRATION IN ADIABATIC DISKS

We consider a planet of mass M_p orbiting a star of mass M_* on a fixed circular orbit of radius a and of angular frequency Ω_p. The planet is immersed in a gaseous protoplanetary disk, such that the planetary orbit is coplanar with the disk and prograde. We use Σ to denote the disk's surface density, P to denote the vertically integrated pressure, and T to denote the vertically averaged temperature. We assume that the disk's gas follows the ideal gas law, which reads

$$\begin{equation} P=\frac{{\cal R}\Sigma T}{\mu }, \end{equation} \tag{ 1 }$$

in which ${\cal R}$ is the ideal gas constant and μ is the mean molecular weight. Therefore, the disk isothermal sound speed is

$$\begin{equation} c_s^{\rm iso}=\sqrt{\frac{{\cal R}T}{\mu }}, \end{equation} \tag{ 2 }$$

while adiabatic sound waves propagate at the speed

$$\begin{equation} c_s=\sqrt{\frac{\gamma {\cal R}T}{\mu }}, \end{equation} \tag{ 3 }$$

where γ is the usual adiabatic index. We shall consider the gas entropy in two flavors. The first of these, which we denote with a lower s, is oftentimes used in the astrophysical literature owing to its simplicity; it reads

$$\begin{equation} s = \frac{P}{\Sigma ^\gamma }. \end{equation} \tag{ 4 }$$

Nevertheless, we shall also need an expression for the entropy that is compliant with the original definition of entropy (the entropy exchanged is the heat exchanged divided by the temperature). We denote the latter with an upper S:

$$\begin{equation} S = \frac{{\cal R}}{(\gamma -1)\mu }\log \left(\frac{P}{\Sigma ^\gamma }\right) = \frac{{\cal R}}{(\gamma -1)\mu }\log s. \end{equation} \tag{ 5 }$$

Another thermodynamics variable that we shall need is the enthalpy η, which reads for an ideal gas

$$\begin{equation} \eta = \frac{\cal R}{\mu }\cdot \frac{\gamma }{\gamma -1}T. \end{equation} \tag{ 6 }$$

In what follows we shall assume, for the sake of simplicity, that ${\cal R}/\mu = 1$, which amounts to changing the units of temperature, entropy, and enthalpy, without loss of generality.

We denote by T₀ the temperature of the unperturbed disk, which we assume to be uniform.

We identify a location in the disk by its distance r to the star and its azimuth ϕ with respect to the planet. The gas has a radial velocity v_r and an azimuthal velocity v_ϕ in the frame corotating with the planet (hence its angular frequency in an inertial frame is Ω = v_ϕ/r + Ω_p). The radius of corotation r_c is the location in the unperturbed disk where the material has the same angular frequency as the planet: Ω(r_c) = Ω_p. We will make use of the distance x to corotation: x = r − r_c. We also consider the disk pressure scaleheight H = c_s/Ω and the aspect ratio h = H/r.

The gravitational potential exerted on the disk can be decomposed as the sum of the stellar potential Φ_* = −GM_*/r, of the planetary potential Φ_p, and of the indirect potential Φ_i. The expressions of these last two terms are, respectively,

$$\begin{equation} \Phi _p = -\frac{GM_p}{(r^2-2ra\cos \phi +a^2+\epsilon ^2)^{1/2}}, \end{equation} \tag{ 7 }$$

where is the softening length of the planetary potential, and

$$\begin{equation} \Phi _i = \frac{GM_p}{a^2}r\cos \phi =qa\Omega _p^2r\cos \phi, \end{equation} \tag{ 8 }$$

where q = M_p/M_*.

Using the notation of Baruteau & Masset (2008), we define the gradient of vortensity across corotation as

$$\begin{equation} {\cal V}=\frac{d\log \left(\Sigma /\omega \right)}{d\log r}=\frac{3}{2}+\frac{\Sigma _c}{r_c} \left.\frac{d\Sigma }{dr}\right|_{r_c}, \end{equation} \tag{ 9 }$$

where ω = (1/r)∂_r(r²Ω) is the vertical component of the flow vorticity and Σ_c is the unperturbed surface density at the corotation radius. Note that instead of the vortensity, which is ω/Σ, we shall oftentimes consider its inverse, Σ/ω. One reason for this choice is that if we consider the flow to be that of the unperturbed disk, the perturbations of the inverse of the vortensity can directly be converted into perturbations of density, whose impact on the torque is straightforward. Finally, we define the gradient of entropy across corotation as

$$\begin{equation} {\cal S}=\frac{1}{\gamma }\frac{d\log s_0}{d\log r}=\frac{\gamma -1}{\gamma }r\partial _rS_0. \end{equation} \tag{ 10 }$$

The governing equations of the flow are the equation of continuity, the Euler equation, and the energy equation, together with the closure relationship provided by the equation of state. The equation of continuity, in the frame corotating with the planet, reads

$$\begin{equation} \partial _t\Sigma +\frac{1}{r}\partial _r(\Sigma rv_r)+\frac{1}{r}\partial _\phi (\Sigma v_\phi)=0. \end{equation} \tag{ 11 }$$

The Euler equations in r and ϕ, respectively, read as

$$\begin{equation} \partial _tv_r+v_r\partial _rv_r+\frac{v_\phi }{r}\partial _\phi v_r-r\Omega _p^2-2\Omega _pv_\phi -\frac{v_\phi ^2}{r} =-\frac{\partial _rP}{\Sigma }-\partial _r\Phi, \end{equation} \tag{ 12 }$$

and

$$\begin{equation} D_tj=-\frac{\partial _\phi p}{\Sigma }-\partial _\phi \Phi, \end{equation} \tag{ 13 }$$

where $D_t\equiv \partial _t+v_r\partial _r+\frac{v_\phi }{r}\partial _\phi$ and j = r²Ω is the specific angular momentum. The equation of state that we adopt is that of ideal gases, which reads P = ΣT. In that case, the internal energy density is e = p/(γ − 1), and the energy equation reads

$$\begin{equation} \Sigma D_t\left(\frac{e}{\Sigma }\right)=-p\vec{\nabla }.\vec{v}. \end{equation} \tag{ 14 }$$

This equation does not include source or sink terms of energy, as we assume the flow to be adiabatic.

Since different authors represent the horseshoe region in different ways, we find it useful to refer to properties of this region in a manner independent of the representation. Orienting the azimuth rotation-wise and assuming that, in the corotating frame, the planet is located at azimuth ϕ = 0, we say that something occurs in front of the planet if it occurs at ϕ>0 and behind or on the rear side of the planet if it occurs at ϕ < 0. Similarly, we refer to the upstream part of the horseshoe streamlines as the set of fluid elements that have not yet performed a horseshoe U-turn, whereas the downstream part corresponds to the set of fluid elements that have already performed their U-turn.

Our main assumptions are as follows.

• 1.  
  We assume that there is only one stagnation point in the vicinity of the planet. This is usually not the case in barotropic situations, in which one generally has two X-stagnation points on the corotation, on each side of the planet (see, e.g., Masset et al. 2006, or Paper I). In the adiabatic case, nevertheless, we usually have only one stagnation point, unless either the entropy gradient is very small or the potential's softening length is very small. We shall therefore assume in what follows that there is no ambiguity when we refer to the stagnation point. We will make some additional remarks on this assumption in Section 8.1.
• 2.  
  We assume that the initial temperature field is uniform. This assumption is very similar to the assumption of global isothermality that is required in the isothermal case to carry out a rigorous horseshoe drag calculation (see Paper I). This assumption is required in order for the quantity that we shall define in Section 2.5 to be conserved along a fluid element's path. We will mention in Section 8.3 how our results can be generalized to the case of an arbitrary temperature profile.
• 3.  
  Finally, we assume that the flow is in a steady state in the frame corotating with the planet and that the flow upstream of the horseshoe U-turns is unperturbed. We note that these two assumptions are contradictory: reaching a steady state implies that the flow has executed many horseshoe librations; hence, the upstream flow cannot be that of the unperturbed disk. Nevertheless, we restrict ourselves, in this analysis, to a study of the unsaturated torque value. This implies that we consider the flow on a timescale (1) longer than what it takes to execute a horseshoe U-turn, so that the flow has reached a steady state in some region of interest enclosing the planet and (2) shorter than the horseshoe libration time, so that in the region of interest any fluid element has executed at most one horseshoe U-turn.

We consider the Jacobi-like Bernoulli constant whose expression, in a steady state, is

$$\begin{equation} B_J=\frac{r^2(\Omega -\Omega _p)^2+v_r^2}{2}+\Phi _*+\Phi _p+\Phi _i -\frac{1}{2}r^2\Omega _p^2+\eta. \end{equation} \tag{ 15 }$$

Since B_J is conserved along the path of a fluid element, the quantity defined as

$$\begin{equation} G=B_J-T_0S \end{equation} \tag{ 16 }$$

is also conserved in an isentropic flow. Far from the planet, assuming a purely azimuthal motion, we can write, from Equations (15) and (16),

$$\begin{equation} \partial _rG=(\Omega -\Omega _p)r\omega +(T-T_0)\partial _rS, \end{equation} \tag{ 17 }$$

where we have used the rotational equilibrium, which can be deduced from Equation (12) by letting v_r ≡ 0, and which reads

$$\begin{equation} \partial _r\phi _*+\partial _r\eta =r\Omega ^2+T\partial _rS. \end{equation} \tag{ 18 }$$

As a consequence, in the unperturbed flow (T ≡ T₀), G is maximal at corotation, as is B_J in an isothermal disk, and it admits the following expansion to second order in x:

$$\begin{equation} G=G_c+2AB x^2, \end{equation} \tag{ 19 }$$

where G_c is the value of G at corotation, A = (1/2)rdΩ/dr is the first Oort's constant, and B = (1/2r)d(r²Ω)/dr is the second Oort's constants, which are to be evaluated at corotation.

We consider a domain D of the disk consisting of an angular sector centered on the star that encloses the planet, as depicted in Figure 1. The azimuth of the boundaries (I) and (III) is chosen sufficiently large so that the material that crosses them can be considered to have a purely circular motion, while the boundaries (II) and (IV) enclose the horseshoe region, and are chosen much further away from corotation than the horseshoe separatrices, so that they can be considered as circular. We note that in a real situation, this last requirement is hindered by the wake, which corresponds to epicyclic motion excited at Lindblad resonances. Nevertheless we discard this behavior in the present work, which amounts to considering an isolated corotation region. We want to avoid the flux of material into the domain D through the boundaries (II) and (IV), so we adopt streamlines to define them. These two boundaries are therefore associated with some value of the G invariant. For convenience, we use the same value for G on both sides, and we denote it as G_∞.⁴

Zoom In Zoom Out Reset image size

The torque exerted on the planet by the material included in the domain D reads

$$\begin{equation} \Gamma = \int \!\int _D \Sigma (\partial _\phi \Phi)\,d\phi \,rdr\,, \end{equation} \tag{ 20 }$$

We can transform the integrand of Equation (20) as follows:

$$\begin{equation} \Sigma (\partial _\phi \Phi)r=-r\partial _\phi P-\partial _\phi (v_\phi \Sigma j)-\partial _r(r\Sigma v_r j), \end{equation} \tag{ 21 }$$

where we have made use of the assumption of a steady state in the corotating frame and used Equations (11) and (13).

Since the boundaries of the domain D are sufficiently far from the planet that we can neglect radial motions, Equation (20) can be recast as

$$\begin{equation} \Gamma = \left[\int _{r^-_{\infty }}^{r^+_{\infty }} [rP+r(\Omega -\Omega _p)\Sigma j]dr\right]_F^R, \end{equation} \tag{ 22 }$$

where r⁺_∞ (r⁻_∞) is the radius of the outer (inner) streamline that has G = G_∞. The R and F notation at the right bracket of Equation (22) means that the integration over r has to be performed, respectively, at the rear or front side of the domain D (boundaries III and I in Figure 1). In Equation (22), the first term of the integrand represents the pressure torque exerted on the material enclosed within the domain D, while the second term represents the budget of angular momentum brought to this region by advection. Since the flow is steady in the corotating frame, the angular momentum of this domain is constant in time and the torque is therefore integrally transmitted to the planet. We recognize in the second term of the integrand of Equation (22) the classical horseshoe drag expression (Ward 1991, 1992; Masset 2001; Masset & Papaloizou 2003; Masset et al. 2006; Paardekooper & Papaloizou 2009a), but this equation also shows the pressure contribution, which has been overlooked in previous analysis.

We assume that in the unperturbed flow, far from the planet, lines of constant entropy and of constant G coincide. This is not exactly true: in the unperturbed disk the lines of constant entropy are circles centered on the primary while, owing to the presence of the potential's indirect term given by Equation (8), the iso-G lines have a variable distance to corotation in the horseshoe region (an extreme example being provided by the tadpole separatrix). Our assumption is therefore more appropriate for a situation in which we discard the indirect term of the potential or for a shearing sheet situation. Nevertheless, we shall make hereafter this assumption, anticipating that our torque expression does not depend sensitively on it.

This assumption allows us to define unambiguously the derivative ∂_GS, which can be obtained from the measure of S and G on two neighboring streamlines, since these variables are intrinsic to the line. The quantity ∂_GS is therefore also intrinsic to the streamline, and hence conserved along it.

In the unperturbed disk, this derivative has the following expression:

$$\begin{equation} \partial _GS=\frac{\partial _rS_0}{4ABx}. \end{equation} \tag{ 23 }$$

It therefore diverges at corotation.

Mass conservation yields a relationship which will be useful to evaluate the production of vortensity during horseshoe U-turns. We consider a horseshoe stream tube, of width δG. The mass flux across the tube is

$$\begin{equation} \dot{M}=\frac{\delta G}{\partial _rG}r(\Omega -\Omega _p)\Sigma. \end{equation} \tag{ 24 }$$

We can transform Equation (17) to write

$$\begin{equation} (1-\delta T\cdot \partial _GS)\partial _rG=(\Omega -\Omega _p)r\omega, \end{equation} \tag{ 25 }$$

where δT = T − T₀. This eventually gives the mass flux

$$\begin{equation} \dot{M}=\delta G\frac{\Sigma }{\omega }(1-\delta T\partial _GS). \end{equation} \tag{ 26 }$$

We can therefore write the following relationship:

$$\begin{equation} \left.\frac{\Sigma }{\omega }(1-\delta T\partial _GS)\right|_d=\left.\frac{\Sigma }{\omega }\right|_u, \end{equation} \tag{ 27 }$$

where the subscript d (u), respectively, denotes the quantities downstream (upstream) of a horseshoe U-turn.

We now perform a change of variable in the integrals of Equation (22). We choose as the new variable of integration the value of the G invariant. We therefore have to split the intervals of integration into intervals over which the dependence of G on r is continuous and monotonic. Anticipating the results exposed in the following sections, we note that the value of the G invariant is not necessarily the same for the front and rear separatrices, as shown in Figure 2. We denote their respective values with G₊ and G₋.

Zoom In Zoom Out Reset image size

Using Equation (25), Equation (22) is therefore recast as

$$\begin{eqnarray} \Gamma &=& \left[ \int _{r_{\infty }^-}^{r_{\infty }^+}rpdr +\int _{G_{\infty }}^{G_c} (1-\delta T\partial _GS)\frac{\Sigma }{\omega }jdG\right.\\ \nonumber &&\left.+\int _{G_c}^{G_-}\! (1-\delta T\partial _GS)\frac{\Sigma }{\omega }jdG +\!\int _{G_+}^{G_{\infty }}\! (1-\delta T\partial _GS)\frac{\Sigma }{\omega }jdG \right]_R\\ \nonumber &&-\left[ \int _{G_{\infty }}^{G_-} (1-\delta T\partial _GS)\frac{\Sigma }{\omega }jdG +\int _{G_+}^{G_c} (1-\delta T\partial _GS)\frac{\Sigma }{\omega }jdG\right.\\ \nonumber &&\left.+\int _{G_c}^{G_{\infty }} (1-\delta T\partial _GS)\frac{\Sigma }{\omega }jdG + \int _{r_{\infty }^-}^{r_{\infty }^+}rpdr\right]_F.\\ &&\nonumber \end{eqnarray} \tag{ 28 }$$

We split the integration on the downstream side of the horseshoe U-turns in order to account for the possible discontinuity of G at the separatrices. More precisely, the first integral over G in Equation (28) corresponds to an integral over intervals (1) and (2) (see Figure 2) and the two following integrals (from G_c to G₋ and G₊ to G_∞) correspond, respectively, to intervals (3) and (4). In a similar manner, on the front side, the intervals of integration over G are, in order of appearance, (5), (6) and eventually (7) and (8) together. We also note that in Equation (28) we have substituted the total pressure P by the perturbed pressure p, since the unperturbed pressure field is axisymmetric and does not contribute to the pressure torque.

Other comments are in order regarding the transformation of Equation (22) into Equation (28). The integrand in Equation (22) is regular, and although it may display discontinuities (we shall see that Ω is discontinuous at the downstream separatrices, as well as Σ) it does not contain any singularity. In contrast, the integrands of Equation (28) feature the vorticity, which is singular wherever Ω is discontinuous. Since the dependence of G upon r on the different intervals of integration is regular, the integrals on G must exclude any singularity (simply because the contribution of the integrand of Equation (22) over any vanishingly small interval in r is vanishingly small, and the same must be true of the equivalent integral in G). As a consequence, any integral of Equation (28) which has a boundary corresponding to a downstream separatrix must be understood as excluding the possible singularities at this boundary. For instance, the second integral over G in Equation (28) must be understood as

$$\begin{equation} \int _{G_c}^{G_-}(1-\delta T\partial _GS)\frac{\Sigma }{\omega } jdG\equiv \lim _{\delta G\rightarrow 0^+} \int _{G_c}^{G_-+\delta G}(1-\delta T\partial _GS)\frac{\Sigma }{\omega } jdG. \end{equation} \tag{ 29 }$$

We also note that including such singularities would be mathematically ill-defined anyway, because they lie on the edge of the integration domain.

Equation (28) can be further simplified as explained below. The specific angular momentum j of a fluid element that has some value of the G invariant can be written as

$$\begin{equation} j=j_0+\delta j, \end{equation} \tag{ 30 }$$

where j₀ is the specific angular momentum of fluid elements in the axisymmetric, unperturbed flow which have the same value of G and δj is the difference, which we wish to evaluate. Azimuthally far from the planet, Equation (16) can be recast as

$$\begin{equation} G=\frac{r^2(\Omega -\Omega _p)^2}{2}+\Phi _*-\frac{1}{2}r^2\Omega _p^2+\eta -T_0S, \end{equation} \tag{ 31 }$$

where we have neglected the indirect term (as discussed in Section 3.2). Differentiating Equation (31), we can write, to first order in the perturbation,

$$\begin{equation} 0 =r\delta r(\Omega ^2-2\Omega \Omega _p)+r^2\delta \Omega (\Omega -\Omega _p) +\partial _r\Phi _*\delta r+\frac{p}{\Sigma _0}. \end{equation} \tag{ 32 }$$

The left-hand side of Equation (32) cancels out because, as stated above, we consider a fluid element that has the same value of G in the unperturbed flow and in the perturbed one. Using the rotational equilibrium of the unperturbed flow, and writing δP = ∂_rP₀δr + p, we can recast Equation (32) as

$$\begin{equation} 0 = (\Omega -\Omega _p)\delta j+\frac{p}{\Sigma _0}. \end{equation} \tag{ 33 }$$

We therefore have an expression for δj that involves the perturbation of pressure p. Since we consider the variation of angular momentum of a fluid element for a constant value of G, it is natural to change to G the variable of the integrals of the pressure torque. The generic integrand of Equation (28) therefore becomes

$$\begin{eqnarray} \frac{\Sigma }{\omega }(1-\delta T\partial _GS)j+\frac{rp}{\partial _rG}&=& \frac{1}{\partial _rG}[\Sigma r(\Omega -\Omega _p)j_0(G)\nonumber\\ && + \Sigma r(\Omega -\Omega _p)\delta j+rp]\nonumber\\ &=&\frac{\Sigma }{\omega }(1-\delta T\partial _GS)j_0(G)\nonumber\\ && +\frac{1}{\partial _rG}\frac{\Sigma _0-\Sigma }{\Sigma _0}rp, \end{eqnarray} \tag{ 34 }$$

where we have used Equations (25) and (33).

The last term on the right-hand side of Equation (34) cancels out to first order in the perturbation. Furthermore, the parity in ϕ of the perturbed density and of the perturbed pressure (which will be studied in Section 5.3) is such that the rear and front contributions of this term to Equation (28) cancel out. Therefore we can discard the pressure torque in Equation (28), provided we use the unperturbed value of the specific angular momentum j₀. Since G has a maximum at corotation, the notation j₀(G) contains an ambiguity, as one has to specify if j₀ has to be evaluated inside or outside of corotation. In what follows, we use the notation j^±₀ to remove this ambiguity and drop the G dependence for the sake of brevity.

The fact that one can neglect the pressure torque and use the value of the unperturbed angular momentum has been demonstrated in Paper I for isothermal disks. It is linked (1) to the fact that the evanescent pressure waves excited by the perturbation of the flow on the downstream side of the horseshoe U-turns do not alter the total horseshoe drag, even though they redistribute the perturbation of density and therefore change the local torque density and (2) to the fact that the impact on the torque of the rear-front asymmetry due to the feedback of the evanescent waves cancels out.

Using the simplification described above, we transform Equation (28) as follows:

$$\begin{eqnarray} \Gamma &=& \bigg[ \int _{G_{\infty }}^{G_c} \left.\frac{\Sigma }{\omega }\right|_{0^-}j_0^-dG +\int _{G_c}^{G_-} \left.\frac{\Sigma }{\omega }\right|_{0^-}j_0^+dG\nonumber\\ && +\int _{G_+}^{G_{\infty }} \left.\frac{\Sigma }{\omega }\right|_{0^+}j_0^+dG \bigg]_R -\bigg[ \int _{G_{\infty }}^{G_-} \left.\frac{\Sigma }{\omega }\right|_{0^-}j_0^-dG \nonumber\\ && +\int _{G_+}^{G_c} \left.\frac{\Sigma }{\omega }\right|_{0^+}j_0^-dG +\int _{G_c}^{G_{\infty }} \left.\frac{\Sigma }{\omega }\right|_{0^+}j_0^+dG \bigg]_F.\\ &&\nonumber \end{eqnarray} \tag{ 35 }$$

Apart from discarding the pressure contribution and changing j into j₀ in Equation (35), we have introduced an additional simplification by making use of Equation (27) and by referring to the unperturbed, upstream value of Σ/ω(1 − δT∂_GS). We explicitly specify by the use of the index 0^± whether the unperturbed vortensity must be considered outside or inside corotation.

We note (see Figure 2) that the integrals over the intervals (1) and (5), on the one hand, and over the intervals (4) and (8) on the other hand, exactly cancel each other; hence, the torque expression can further be reduced to

$$\begin{equation} \Gamma = -\int _{G_-}^{G_c}\left.\frac{\Sigma }{\omega }\right|_{0^-}\Delta j_0dG +\int _{G_+}^{G_c}\left.\frac{\Sigma }{\omega }\right|_{0^+}\Delta j_0dG, \end{equation} \tag{ 36 }$$

where Δj₀ = j⁺₀ − j⁻₀. The corotation torque has therefore an expression very similar to the usual (barotropic) case, when integrated over a Bernoulli-like variable (Masset & Papaloizou 2003), except that the domain of integration does not have the same limits on the front side and on the rear side. Denoting by G_s the arithmetic mean of the neighboring values G₋ and G₊, we can transform Equation (36) into

$$\begin{equation} \Gamma \approx \int _{G_s}^{G_c} \left(\left.\frac{\Sigma }{\omega }\right|_{0^+} - \left.\frac{\Sigma }{\omega }\right|_{0^-} \right)\Delta j_0dG -\frac{\Sigma _0}{2B}\Delta j_0^s\Delta G, \end{equation} \tag{ 37 }$$

where Δj^s₀ is the value of Δj₀ for G = G_s and ΔG = G₊ − G₋. The left term on the right-hand side of Equation (37) is the standard horseshoe drag, which scales with the vortensity gradient, and the right term, which we hereafter denote with Γ₁, is an additional contribution that arises by relaxing the barotropic hypothesis. Evaluating this contribution therefore amounts to evaluating ΔG. This is the purpose of the following section.

The discontinuity in G that we wish to evaluate arises from the presence, at the downstream separatrices, of a contact discontinuity (see Figure 2). Formally, an inspection of Equations (15) and (16) reveals that this discontinuity can come, far from the planet, from a discontinuity in the enthalpy, in the entropy, and in the azimuthal velocity (the other terms being continuous across the separatrices). The situation becomes simpler at the stagnation point, where the velocity vanishes and the discontinuity in G is entirely accounted for by a discontinuity of entropy and enthalpy. We therefore have

$$\begin{equation} \Delta G = \eta _+-\eta _--T_0(S_+-S_-), \end{equation} \tag{ 38 }$$

where the + (−) sign refers to quantities outside, or in front of (inside, or behind) the stagnation point. Equation (38) can be transformed into

$$\begin{equation} \Delta G = \big(s_+^{1/\gamma }-s_-^{1/\gamma }\big)\frac{\gamma }{\gamma -1}P_s^{(\gamma -1)/\gamma } -\frac{T_0}{\gamma -1}(\log s_+-\log s_-), \end{equation} \tag{ 39 }$$

where P_s is the value of the pressure at the stagnation point (this value is well defined at this point in contrast to the other thermodynamics variables since the pressure field is continuous). Using the first order expansions: s^1/γ₊ − s^1/γ₋ ≈ (s₊ − s₋)s^1/γ−1₀/γ and log s₊ − log s₋ ≈ (s₊ − s₋)/s₀, we are led to

$$\begin{equation} \Delta G = \frac{1}{\gamma -1} s_0^{1/\gamma -1} \big(P_s^{(\gamma -1)/\gamma }-P_0^{(\gamma -1)/\gamma }\big) (s_+-s_-), \end{equation} \tag{ 40 }$$

where P₀ and s₀ are, respectively, the pressure and the entropy at the location of the stagnation point in the unperturbed disk.

In the particular case of a small perturbation, |P_s − P₀| ≪ P₀, the discontinuity in G takes the form

$$\begin{equation} \Delta G=\frac{1}{\gamma }\frac{p_s}{\Sigma _0}\frac{s_+-s_-}{s_0}, \end{equation} \tag{ 41 }$$

where p_s = P_s − P₀ is the pressure perturbation at the stagnation point.

In order to proceed in the evaluation of ΔG, we use the conservation of the G invariant along the upstream separatrices. At this stage of derivation, we note that numerical simulations reveal that the horseshoe region in the adiabatic regime is rather asymmetric, so that the distance of the separatrices depends on the quadrant under consideration (rear versus front and upstream versus downstream). Denoting with x^R_u the distance of the upstream separatrix from corotation on the rear side, we have, by evaluating G, respectively, on that separatrix far from the planet, by means of Equation (19), and at the stagnation point

$$\begin{equation} G_-=G_c+2AB\big(x_u^R\big)^2 \end{equation} \tag{ 42 }$$

and

$$\begin{equation} G_-=G_c-\eta _c+\eta _-+\Phi _p^s+T_0(S_c-S_-), \end{equation} \tag{ 43 }$$

where Φ^s_p is the planetary potential at the stagnation point, η_c and S_c are, respectively, the enthalpy and entropy at corotation in the unperturbed disk, and where we have used the relationship G_c = Φ_*(r_c) − (1/2)r²_cΩ²_p + η_c − T₀S_c to derive Equation (43). From Equations (42) and (43), we deduce

$$\begin{equation} 2AB\big(x_u^R\big)^2=\eta _--\eta _c+\Phi _p^s+T_0(S_c-S_-). \end{equation} \tag{ 44 }$$

In a similar manner, denoting with x^F_u the distance of the upstream separatrix from corotation on the front side, we have

$$\begin{equation} 2AB\big(x_u^F\big)^2=\eta _+-\eta _c+\Phi _p^s+T_0(S_c-S_+). \end{equation} \tag{ 45 }$$

From Equation (40), we have

$$\begin{equation} \Delta G=\frac{1}{\gamma }\left[\left(\frac{s_0}{s_\pm }\right)^{1/\gamma }\eta _\pm -\eta _c\right] \frac{s_+-s_-}{s_0}; \end{equation} \tag{ 46 }$$

hence, we have

$$\begin{equation} 2\Delta G=\frac{1}{\gamma }(\eta _++\eta _--2\eta _c)\frac{s_+-s_-}{s_0}, \end{equation} \tag{ 47 }$$

where we have discarded the additional term that scales with ${\cal S}^2$. Similarly, adding Equations (44) and (45) yields

$$\begin{equation} 2AB\big[\big(x_u^R\big)^2+\big(x_u^F\big)^2\big]=\eta _++\eta _--2\eta _c+2\Phi _p, \end{equation} \tag{ 48 }$$

where we have omitted the term T₀(2S_c − S₋ − S₊), which scales as h²Ω²x², and is hence negligible compared to the left-hand side. Using Equations (47) and (48), we obtain

$$\begin{equation} \Delta G=\frac{1}{\gamma }\left\lbrace AB\big[\big(x_u^R\big)^2+\big(x_u^F\big)^2\big]-\Phi _p^s\right\rbrace \frac{s_+-s_-}{s_0}. \end{equation} \tag{ 49 }$$

Since the entropy is conserved during the advection along the upstream separatrices, we have

$$\begin{equation} \frac{s_+-s_-}{s_0}=\gamma {\cal S}\frac{x_u^F+x_u^R}{a}. \end{equation} \tag{ 50 }$$

The discontinuity of G finally reads

$$\begin{equation} \Delta G = \left\lbrace AB\big[\big(x_u^R\big)^2+\big(x_u^F\big)^2\big]-\Phi _p^s\right\rbrace {\cal S}\frac{x_u^F+x_u^R}{a} \end{equation} \tag{ 51 }$$

and the associated torque excess is, using Equation (37),

$$\begin{equation} \Gamma _{1}=-\big(x_u^F+x_u^R\big)^2\Sigma _0{\cal S}\left\lbrace AB\big[\big(x_u^R\big)^2+\big(x_u^F\big)^2\big]-\Phi _p^s\right\rbrace, \end{equation} \tag{ 52 }$$

where we have used the equality Δj^s₀ = 2Ba(x^R_u + x^F_u). Since P_s > P₀, ΔG has the same sign as ${\cal S}$, and the torque excess scales with the negative of the entropy gradient, as already observed in previous works (Baruteau & Masset 2008; Paardekooper & Papaloizou 2008). We can render Equation (52) a bit more compact by defining the averaged upstream half-width:

$$\begin{equation} x_s=\sqrt{\frac{\big(x_u^R\big)^2+\big(x_u^F\big)^2}{2}}. \end{equation} \tag{ 53 }$$

Considering that the relative asymmetry between the rear and front part of the horseshoe, albeit larger than typically observed in barotropic situations (see Paper I), remains moderate, we have

$$\begin{equation} x_s\approx \frac{1}{2}\big(x_u^R+x_u^F\big), \end{equation} \tag{ 54 }$$

and therefore

$$\begin{equation} \Gamma _{1}\approx -4x_s^2\Sigma _0{\cal S}\left(2ABx_s^2-\Phi _p^s\right). \end{equation} \tag{ 55 }$$

This expression displays explicitly the dependence of the torque excess on the entropy gradient. As we shall see hereafter, it also brings in a dependence on the potential's softening length, as both x_s and Φ^s_p depend on this quantity.

We used the code FARGO⁵ (Masset 2000a, 2000b; Baruteau & Masset 2008) in its adiabatic and isothermal versions. Unless otherwise stated, the polar mesh extends from r_min = 0.4 to r_max = 1.8, and has a resolution of 1200 × 1200, the zone boundaries being evenly spaced. The planet is located at r = 1, and it is held on a fixed circular orbit. Its mass ratio to the primary is q = 3 × 10⁻⁶ (i.e., the planet has one Earth mass, if the central object has a solar mass). The disk aspect ratio at r = 1 is h = 0.05. The adiabatic index of the gas is γ = 1.4. Damping boundary conditions (as described in de Val-Borro et al. 2006) are used at the inner and outer edges of the mesh. The frame corotates with the planet, and the Coriolis force is conservatively implemented (Kley 1998). The surface density is a power law of the radius: Σ(r) ∝ r^−σ.

In our first set of calculations, the slope of the surface density is set to σ = +3/2 (hence the vortensity-related part of the corotation torque cancels out), and the potential softening length ranges from $\hat{\epsilon }=0.1$ to 10, where $\hat{\epsilon }=\epsilon /H$ is a dimensionless value of the softening length. The logarithm of $\hat{\epsilon }$ is evenly spaced, and 21 calculations are performed, so that $\hat{\epsilon }=10^{i/10-1}$, with 0 ⩽ i ⩽ 20. These calculations were performed both for an adiabatic flow (for brevity, we denote them with EA_i,0⩽i⩽20) and for an isothermal flow (we denote them with EI_i,0⩽i⩽20). Since the initial temperature field is flat, the index of the entropy gradient is ${\cal S}=\sigma (\gamma -1)/\gamma =+3/7$. We therefore expect a negative torque excess in the adiabatic case.

In our second set of calculations, the slope of the surface density ranges from σ = −4 to σ = 4, the softening length being kept fixed at $\hat{\epsilon }=0.3$. We have divided the slope interval into 41 calculations such that σ_i = −4 + 0.2 × i(0 ⩽ i ⩽ 40). These calculations were performed both for an adiabatic flow (we denote them with SA_i,0⩽i⩽40) and for an isothermal flow (we denote them with SI_i,0⩽i⩽40). Since the initial temperature field is flat, the index of the entropy gradient is

$$\begin{equation} {\cal S}=\frac{\gamma -1}{\gamma }\sigma. \end{equation} \tag{ 56 }$$

It therefore ranges from −8/7 to +8/7.

As underlined in Section 3.4, the horseshoe region in adiabatic calculations is more asymmetric than in the isothermal case. It also has a different dependence, more complex, on the disk parameters. This can be seen in Figure 3, in which we plot the different values of x_s inferred from a streamline analysis on the runs of the series SA_i. In particular, the naive expectation that

$$\begin{equation} x_s^{\rm adi}= x_s^{\rm iso}/\gamma ^{1/4} \end{equation} \tag{ 57 }$$

is true only in the barotropic case $({\cal S}=0)$. For nonvanishing entropy gradients, the averaged upstream half-width exhibits a V-shaped dependence on ${\cal S}$, whereas the isothermal case is insensitive to this parameter.

Zoom In Zoom Out Reset image size

This more complex dependence of the width of the horseshoe region has a consequence on the evaluation of the torque excess. A simple estimate of the latter is

$$\begin{equation} \Gamma ^r_1=\Gamma ^A-\frac{\Gamma ^I}{\gamma }, \end{equation} \tag{ 58 }$$

but this estimate assumes x_s,adi = x_s,iso/γ^1/4, so that one gets rid of the vortensity-related part of the horseshoe drag, retaining only the excess. This is not exactly true, since the adiabatic case has a horseshoe region in general larger than expected. In order to fully eliminate the vortensity part of the estimate, one should therefore use the expression

$$\begin{equation} \Gamma ^c_1=\Gamma ^A-\frac{\Gamma ^I}{\gamma }-\frac{3}{4}{\cal V}\Omega _p^2\Sigma _0\left(x_{s, \rm adi}^4-\frac{x_{s, \rm iso}^4}{\gamma }\right). \end{equation} \tag{ 59 }$$

We shall hereafter refer to the estimate given by Equation (58) as the rough or global estimate (hence the r superscript) and to the estimate of Equation (59) as the corrected estimate (hence the c superscript).

We compare the torque excess, when $\hat{\epsilon }$ varies, as provided by Equation (55) and as provided by a direct comparison of the torque value in adiabatic and isothermal calculations. The results are presented in Figure 4. We see that there is a good overall agreement between the direct measure of Γ^c₁ and the result of Equation (55), as the relative difference is at most ∼10% over the range $0.1\le \hat{\epsilon }\le 1$. The value of x_s, as measured on the output at t = 40 orbits and ϕ = ±1 rad, is used in the evaluation of Equation (55). Similarly, the location of the stagnation point is determined, and used to evaluate the planetary potential that features in the bracket of Equation (55). We note that in all the runs of this series, there is a unique stagnation point in the planet's vicinity; hence, there is no ambiguity in its determination.

Zoom In Zoom Out Reset image size

There is an approximate $\hat{\epsilon }^{-1}$ scaling of the torque excess in the range $0.1\le \hat{\epsilon }\le 1$ (Baruteau 2008). Figure 4 shows that this trend results partly from the increase of x_s as $\hat{\epsilon }$ decreases (dotted curve) and partly from the increase (in absolute value) of the planetary potential as the softening length decreases (dot-dashed line).

We now evaluate how well Equation (55) accounts for the torque excess in the second series of runs (SA and SI; see Section 4.1), in which we vary the entropy gradient.

The results are presented in Figure 5. They show a satisfactory agreement between the corrected estimate of the torque excess (given by Equation (59), and represented by the thick dashed line) and the theoretical prediction of Equation (55), represented by the thick solid line. This figure also shows that the azimuth of the stagnation point oscillates as the entropy gradient varies. The oscillations are significantly large, and can be unambiguously identified as a finite resolution issue, for the following reasons.

• 1.  
  When the entropy gradient varies, the corotation radius is slightly shifted. If one displays the azimuth of the stagnation point as a function of the corotation radius, the period of the oscillations is equal to the radial resolution.
• 2.  
  We have carried out higher resolution calculations (2000 × 2000, not shown here, as they were performed over a small number of orbits and a steady state was not reached in the vicinity of the planet), which were ran for 10 orbits only (the location of the stagnation point hardly varies afterward). They also yield oscillations of the azimuth of the stagnation point with a shorter period and smaller amplitude.

Zoom In Zoom Out Reset image size

Given the large amplitude of the oscillations of the stagnation point for the resolution that we were able to use for this series of runs, we need to smooth the $\phi _s({\cal S})$ function. This is done simply by performing a linear regression fit, independently for the rear and front stagnation points, as shown in the inset plot of Figure 5. The fitted value is then used in Equation (55), which yields the solid thick curve.

The value of ${\cal S}=+3/7$ used in the series of runs EA and EI corresponds to a location where the error due to the oscillations of the azimuth of the stagnation point is small. We therefore understand from this plot why there is a good agreement between the measured and predicted adiabatic excess in Figure 4. The fact that there is a good agreement for all values of $\hat{\epsilon }$ in Figure 4 suggests that the error on the location of the stagnation point is essentially due to the location of the corotation with respect to the mesh (which does not depend on the softening length). We recover the fact that, for $\hat{\epsilon }=0.3$, |Γ₁| ∼ 1.5 × 10⁻¹¹.

It is noteworthy that, on the left part of the plot of Figure 5, there is a sizable difference between the rough and corrected estimates of the excess (thick dashed curve and dotted curve). The reason for this is straightforward, since for these values of the entropy gradient, the surface density increases steeply outward; hence, the vortensity-related corotation torque is large and positive. For large values of $|{\cal S}|$, the adiabatic horseshoe half-width is significantly larger than x^iso_s/γ^1/4 (see Section 4.2); therefore, a fair fraction of the adiabatic excess is actually due to the boost of the vortensity-related torque.

The discontinuity of the G invariant that we evaluated in Section 3.4 at the stagnation point also exists along the downstream separatrices, as can be seen in Figure 2. Far from the planet, one can write the expression of the G discontinuity across a downstream separatrix as

$$\begin{equation} \Delta G=\frac{1}{\gamma }\frac{p^{\prime }_s}{\Sigma _0}\frac{s_+-s_-}{s_0}+ \frac{r^2}{2}[(\Omega _+-\Omega _p)^2-(\Omega _--\Omega _p)^2], \end{equation} \tag{ 60 }$$

where we have performed the same series of transformations as in Equations (38)–(41), where p'_s is the perturbed pressure at the separatrix and Ω₊ (Ω₋) is the angular velocity on the outside (inside) of the separatrix. If we assume in a first approximation that the perturbation of pressure at the separatrix, far from the planet, is much smaller than at the stagnation point, most of the G discontinuity is ensured by the jump of azimuthal velocity, that is, by a vorticity sheet. Denoting with Δv = r(Ω₊ − Ω₋) the magnitude of the vorticity sheet (i.e., the singular part of the flow's vorticity reads ω = Δvδ(x − x_s), where δ is Dirac's delta function) and assuming that the jump is sufficiently small so that we can write

$$\begin{equation} \Omega _++\Omega _--2\Omega _p=\frac{4Ax}{a}, \end{equation} \tag{ 61 }$$

we are led to

$$\begin{equation} \Delta v=\frac{\Delta G}{2Ax}= \left\lbrace \begin{array}{ll}\displaystyle \frac{\Gamma _1}{4a|A|x_s^2\Sigma _0}&\mbox{ at the rear separatrix}\\ \displaystyle -\frac{\Gamma _1}{4a|A|x_s^2\Sigma _0}&\mbox{ at the front separatrix} \end{array}\right. \end{equation} \tag{ 62 }$$

Therefore, if there is a positive torque excess (hence ${\cal S}<0$), then there is a negative (positive) vorticity sheet in front of the planet (at the rear of the planet) and vice versa. In order to illustrate this, we show in Figure 6 the vorticity field in a run with ${\cal S}=1/3$ and similar characteristics as those of Section 4.1, except that the resolution was set to 2000 zones in azimuth and 8000 zones in radius. We measure an adiabatic torque excess Γ₁ = −1.14 × 10⁻¹¹ M_*Ω²_pr²_p, by performing a run with similar characteristics and an isothermal equation of state. This allows us to predict the jump of angular velocity at the downstream separatrices by the use of Equation (62). The results are depicted in Figure 7, which shows the profile of perturbed azimuthal velocity in the rear and front of the planet as well as the prediction of Equation (62). We note that the jumps are not symmetric, although there is a good agreement between the theoretical prediction of the velocity jump, on the one hand, and the average value of the rear and front velocity jumps, on the other hand. One reason for this is that the perturbed pressure at the separatrices is not negligible compared to the perturbed pressure at the stagnation point. Another one is the asymmetry of the horseshoe region: both separatrices do not have the same value of |(Ω₊ + Ω₋)/2 − Ω_p|.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that for the run presented here, for which ${\cal S}>0$, we have a negative torque excess, hence a negative perturbation of surface density in front of the planet, and a positive one behind the planet (Baruteau & Masset 2008). The fact that we have a positive perturbation of azimuthal velocity in front of the planet (dashed curve of Figure 7) is compatible with an expansion of the material and consequently a decrease of the surface density, because this perturbed velocity acts along with the Keplerian shear. Reciprocally, we also have a positive perturbation of velocity behind the planet, therefore opposed to the Keplerian shear, and compatible with compression of material, hence an increase of the surface density.

Finally, we mention that the vortensity is ill-defined at the separatrices. The vorticity is definable: the integral of its singular part yields the azimuthal velocity jump. However, the vortensity would appear, on the separatrix, as the ratio of a singularity (that of the vorticity) and of a noncontinuous function (owing to the contact discontinuity, the surface density undergoes a jump at the separatrix). Such a ratio is mathematically undefined. In contrast, in the locally isothermal case (see Paper I), the surface density field is continuous and the vortensity is defined everywhere.

In contrast to a barotropic situation, for which the vortensity is necessarily conserved along the streamlines (Lovelace et al. 1999), some thermodynamical driving of the vortensity may occur during the horseshoe U-turns. We also note that far from these U-turns, the vortensity is again conserved, since all quantities depend only on r; hence, the pressure and density gradients are both radial, and hence aligned.

Using Equation (23) and (27), we are left with:

$$\begin{equation} \left.\frac{\Sigma }{\omega }\right|_d = \left.\frac{\Sigma }{\omega }\right|_u \left(1-\delta T\frac{\partial _rS_0}{4ABx}\right). \end{equation} \tag{ 63 }$$

We comment that we have used the relationship ∂_SG = −∂_rS₀/(4ABx) to obtain Equation (63). The bracket of this equation indeed corresponds to material downstream of a horseshoe U-turn. During a U-turn, ∂_GS is conserved, while x, to lowest order, reverses its sign, and Equation (23) has been written for the unperturbed disk, that is to say upstream of horseshoe U-turns. The perturbation of temperature, within the horseshoe region, can be easily worked out assuming the conservation of entropy and an unperturbed pressure. This yields (Baruteau & Masset 2008)

$$\begin{equation} \delta T=-2T_0{\cal S}\frac{x}{r_p}. \end{equation} \tag{ 64 }$$

Using Equations (63) and (64) we get the production of vortensity during a U-turn:

$$\begin{equation} \delta \left(\frac{\Sigma }{\omega }\right)=\left.\frac{\Sigma }{\omega }\right|_d -\left.\frac{\Sigma }{\omega }\right|_u\approx \frac{\gamma T_0{\cal S}^2}{2(\gamma -1)ABr_p^2}\cdot \frac{\Sigma _0}{2B}. \end{equation} \tag{ 65 }$$

We illustrate this in Figure 8, which shows that there is a satisfactory agreement between the vortensity on the downstream side (here, outside, at the rear of the planet) and the vortensity predicted by Equation (65). We also note in this figure how the vortensity near corotation differs from the expected value, both on the inside and on the outside. On the outside, this is expected because the material very close to corotation did not have sufficient time to perform a horseshoe U-turn. On the inside, for 0.996 ⩽ r ⩽ 0.999, libration on tadpole streamlines brings vortensity from the outside, thereby yielding a perturbation of vortensity in the upstream region. More generally, the fact that streamlines are not exactly circular alters the initial profile of vortensity in this region. We also note the spike at the outer separatrix, resulting from the vorticity sheet presented in Section 5.1, where we also pointed out that, formally, the vortensity is ill-defined on the downstream separatrices.

Zoom In Zoom Out Reset image size

We can make the following observations regarding Equation (65).

• 1.  
  There is always a positive production of vortensity during the U-turns (hence a decrease of Σ/ω), regardless of the sign of ${\cal S}$ and regardless of the location (at the rear or in front of the planet).
• 2.  
  The production of vortensity is independent of the planet mass, and it occurs uniformly all over the downstream sides of the horseshoe region. While wider horseshoe U-turns bring the material closer to the planet, where the driving of vortensity is more efficient because the gradients of pressure and density are larger, they are performed faster than the more narrow U-turns, which occur further away from the planet. Equation (65) shows that both effects cancel out, resulting in a flat profile of vortensity perturbation.
• 3.  
  If we translate the production of vortensity into a perturbation of density, assuming an unperturbed velocity field, then this perturbation has the same sign on both sides of the planet. This is compatible with the fact that the bulk term of the horseshoe drag, that is, the first term of Equation (37), is the same as in a barotropic situation. The weak thermodynamical driving of vortensity observed within the horseshoe region has therefore no impact on the torque.
• 4.  
  The increase of vortensity after the first horseshoe U-turn is followed by a decrease of the same magnitude at the following U-turn (in the idealized, inviscid situation considered here). This can be seen from Equation (27), in which we revert the d and u indices. The vortensity of a given fluid element therefore oscillates between its initial upstream value and its subsequent downstream value, and it does not undergo any drift on the long term.

In this section, we wish to determine the shape and amplitude of the pressure disturbances which arise in the co-orbital region as the result of the horseshoe dynamics. In Paper I, we have seen that the perturbation of vortensity was entirely accounted for by evanescent waves, so that the density or pressure response in the co-orbital region was smoothed over a length-scale H. Here, we expect that the perturbation will be split into a localized part, which ends at the contact discontinuity at the separatrix, and a smooth part, corresponding to the launch of evanescent waves.

Far from the planet, we assume purely circular motion; hence, we can describe the disk's perturbed state with the perturbed azimuthal velocity δv_ϕ(x), the perturbed surface density δΣ(x), and the perturbed pressure δP(x). We search the relationship between this set of perturbed functions and the perturbations of vortensity δw(x) and entropy δs(x). We note that x has an arbitrary sign for the rest of this section. If we consider a radial profile in front of the planet, where the downstream separatrix lies inside of corotation, then x_s is meant to be a negative quantity.

Using Equation (12), the disk's rotational equilibrium yields

$$\begin{equation} -2\Omega\delta v_\phi +\frac{\partial _x\delta P}{\Sigma _0} +\frac{\partial _xp_0}{\Sigma _0}\frac{\delta \Sigma }{\Sigma _0}=0. \end{equation} \tag{ 66 }$$

Anticipating what follows, we note that while the radial derivative of an unperturbed quantity ξ₀ is of order ξ₀/r, the radial derivative of a perturbed quantity δξ is of order δξ/H. The second term of Equation (66) is therefore r/H times larger than the third one. This equation can therefore be simplified as

$$\begin{equation} -2\Omega _0\delta v_\phi +\frac{\partial _x\delta P}{\Sigma _0}=0 \end{equation} \tag{ 67 }$$

Denoting by l = Σ/ω the inverse of the vortensity, we can write

$$\begin{equation} \delta l=\frac{\delta \Sigma }{\omega _0}-\frac{\Sigma _0}{\omega _0^2}\partial _x\delta v_\phi, \end{equation} \tag{ 68 }$$

while the perturbation of entropy reads

$$\begin{equation} \delta s=\frac{\delta P}{\Sigma _0^\gamma }-\frac{\gamma s_0}{\Sigma _0}\delta \Sigma. \end{equation} \tag{ 69 }$$

Deriving Equation (67) with respect to x by keeping only the derivatives of rapidly varying terms, that is, the perturbations, and using Equations (68) and (69) so as to keep only the variable δP, we obtain the following differential equation:

$$\begin{equation} \delta P-\frac{c_s^2}{\kappa ^2}\partial _{x^2}^2\delta P = c_s^2\Sigma _0 \frac{\delta u}{u_0}, \end{equation} \tag{ 70 }$$

where κ = (2Ω₀ω₀)^1/2 is the epicyclic frequency, u = s^1/γl, and where we recall that c_s is the adiabatic speed of sound (see Section 2.1). The general solution of Equation (70) is the convolution product of its right-hand side by the Green's kernel K(x), which is the solution of

$$\begin{equation} \delta P-\frac{c_s^2}{\kappa ^2}\partial _{x^2}^2\delta P = \delta (x), \end{equation} \tag{ 71 }$$

and whose expression is

$$\begin{equation} K(x) = (2\gamma ^{1/2}H)^{-1}e^{-|x|/(\gamma ^{1/2}H)}, \end{equation} \tag{ 72 }$$

where we have specialized to the Keplerian case, for which H ≡ c^iso_s/Ω = c_s/(γ^1/2κ). This kernel represents the pressure response to a singular perturbation of u at x = 0, of weight ∫u(x)dx = u₀/(c²_sΣ₀). We note that K(x) has a unitary weight:

$$\begin{equation} \int _{-\infty }^{+\infty }K(x)dx = 1. \end{equation} \tag{ 73 }$$

Once the pressure response is known, one can infer the density response from Equation (69):

$$\begin{equation} \delta \Sigma =\frac{\Sigma _0}{\gamma }\left(\frac{\delta P}{P_0}-\frac{\delta s}{s_0}\right). \end{equation} \tag{ 74 }$$

Denoting δΣ_d = δP/c²_s and δΣ_l = −Σ₀δs/(γs₀), we have

$$\begin{equation} \delta \Sigma =\delta \Sigma _d+\delta \Sigma _l. \end{equation} \tag{ 75 }$$

In Equation (75), the subscript d stands for diffuse, since this component of the surface density scales with the perturbed pressure, which is smoothed over a length-scale H by the action of evanescent waves, as it involves the convolution of an arbitrary function by the kernel K defined in Equation (72). Similarly, the subscript l stands for localized, since this component of the surface density scales with the perturbation of entropy, which is advected by the flow, and remains localized within the horseshoe region.

Since δP is the convolution product of K by c²_sΣ₀δu/u, we have

$$\begin{equation} \delta \Sigma _d=K\ast \Sigma _0 \frac{\delta u}{u_0}, \end{equation} \tag{ 76 }$$

where * denotes the convolution product. Similarly, we have

$$\begin{equation} \delta \Sigma _l=-\frac{\Sigma _0}{\gamma }\frac{\delta s}{s_0}. \end{equation} \tag{ 77 }$$

The expression of the entropy perturbation is straightforward, as the latter is simply advected by the horseshoe dynamics. Since $s\propto r^{\gamma {\cal S}}$, we have $\delta s/s_0=-2\gamma {\cal S}x/r_p$; hence (Baruteau & Masset 2008)

$$\begin{equation} \delta \Sigma _l=2\Sigma _0{\cal S}x/r_p. \end{equation} \tag{ 78 }$$

The expression of the perturbation of u involves three components: the singular contribution of the vorticity sheet at the separatrix, the advection of the vortensity and entropy by the horseshoe dynamics, and the production of vortensity given by Equation (65). We have mentioned in Section 5.1 that the vortensity is ill-defined on the downstream separatrices, since the surface density is discontinuous there. Nevertheless, the jumps of surface density are of order x_s/a; hence to lowest order, the vortensity can be estimated to be Δvδ(x − x_s)/Σ₀ at the separatrices. We can therefore estimate the perturbation in u:

$$\begin{equation} \frac{\delta u}{u_0}=-\frac{\Delta v\delta (x-x_s)}{\omega _0}-\frac{2x}{r_p}({\cal S}+{\cal V}) +\frac{\delta l_{\rm prod}}{l_0}, \end{equation} \tag{ 79 }$$

where we have used the fact that $u\propto r^{{\cal S}+{\cal V}}$ and where, using Equation (65):

$$\begin{equation} \frac{\delta l_{\rm prod}}{l_0}=\frac{\gamma T_0{\cal S}^2}{2(\gamma -1)ABr_p^2}. \end{equation} \tag{ 80 }$$

For the sake of brevity, in what follows, we will keep the notation δl_prod, remembering that it corresponds to a uniform production of vortensity that scales with ${\cal S}^2$, and which has the same sign on both sides of the planet. The perturbation of density finally reads

$$\begin{eqnarray} \delta \Sigma &=&K\ast \left[-l_0\Delta v\delta (x-x_s)-\frac{2x}{r_p}\Sigma _0({\cal S}+{\cal V}) +\frac{\Sigma _0\delta l_{\rm prod}}{l_0}\right]\nonumber\\[-3pt] &&+\frac{2x}{r_p}\Sigma _0{\cal S}. \end{eqnarray} \tag{ 81 }$$

Equation (81) features two kinds of terms as follows.

• 1.  
  Bulk terms (second to fourth term), which are defined for 0 ⩽ x ⩽ x_s, prior to a possible convolution by K.
• 2.  
  An edge term (the first one), which is defined at x = x_s, prior to the convolution by K.

We better understand why the bulk of the horseshoe drag has the same expression as in the barotropic case. The convolution by K, which describes the spread of disturbances by evanescent waves, does not change the linear mass (∫_xδΣ(x)dx) of the perturbation, since K has a unit weight. Therefore, the torque exerted on the planet by the diffuse and localized stripes of perturbed density is the same as if all the disturbances remained localized, that is, as if we omitted the convolution product in Equation (81). If we discard the singularity, which corresponds to the edge term (and, as we shall see, to the adiabatic torque excess), and the production of vortensity, which has the same sign on both sides of the planet and does not affect the torque, we are left with the middle term of the bracket of Equation (81), and the trailing, localized term, which simplify as

$$\begin{equation} \delta \Sigma =-\frac{2x}{r_p}\Sigma _0{\cal V}, \end{equation} \tag{ 82 }$$

exactly as in the barotropic case, prior to the convolution by Green's kernel (see Paper I). These properties are schematically depicted in Figure 9.

Zoom In Zoom Out Reset image size

This first term of Equation (81) corresponds to a single evanescent wave excited at the separatrix. It is the perturbed density associated with this wave that exerts the adiabatic torque excess, since it is the only component that scales with the entropy gradient, owing to the simplification mentioned above. This is in agreement with Equation (37), which shows that the torque excess comes from an edge effect rather than a bulk effect, since it relies on a minute difference between the limit of the integration domain between the rear and front sides. We can check from Equation (62) that this edge term has the correct sign. Since the perturbed density associated with this wave has a sign opposite of Δv, the perturbed density has the same sign as Γ₁ in front of the planet and an opposite sign behind the planet, as expected.

A consequence of these properties is that, in contrast to early expectations, it is not the localized component of the surface density that exerts the torque excess. The latter does exert a torque on the planet, which can be estimated either by means of partial horseshoe drag calculations (Baruteau & Masset 2008) or by direct summation (Paardekooper & Papaloizou 2008). It is found to have the same sign as the excess and the same order of magnitude in usual setups. Nevertheless, this partial torque has to be added to the torque exerted by the diffuse component that scales as ${\cal S}+{\cal V}$, and the result is a corotation torque that scales only with the vortensity gradient. This can be expected on general grounds. The splitting of the perturbation into a localized component (an entropy wave) and pressure-supported waves is very similar to the procedure used by numericists who use characteristic tracing to predict Riemann states, except that here the pressure-supported waves are evanescent instead of propagative. The entropy perturbation regulates how much of the perturbation goes into the entropy wave, but it is the vortensity perturbation that dictates the total amount of mass perturbation, and it is therefore logical that the bulk of the corotation torque scales with the vortensity gradient.

Owing to the simplicity of the functions involved in the convolution product, it is possible to write explicitly the expression of the pressure perturbation. This is done in the Appendix. We show in Figure 10 that the resulting expression, given by Equation (A5), satisfactorily reproduces the perturbations of pressure observed in numerical simulations. Namely, we transform the profile of perturbed pressure from a numerical simulation as follows, in order to filter out the term arising from the production of vortensity, which is poorly reproduced by a rectangular function, as has been shown in Figure 8:

$$\begin{equation} \delta P_{\rm sym}(x) = \frac{1}{2}[\delta P(\phi =2\pi -1,x)-\delta P(\phi =1,-x)]. \end{equation} \tag{ 83 }$$

We then evaluate the theoretical expression of the perturbed pressure given by Equation (A5), in which we use for Γ₁ the estimate provided by Equation (59) (corresponding to the thick dashed curve of Figure 5) and for x_s the arithmetic mean of x^d_F and x^d_R, determined by a streamline analysis. Owing to the symmetrization performed in Equation (83), we discard the term in δl_prod in Equation (A5). The correct agreement between the observed and predicted pressure profiles indirectly confirms that it is indeed the single evanescent wave excited at the separatrix that is responsible for the torque excess. We also recover a feature noticeable in Figure 5. While the plots a and f of Figure 10 correspond to the same absolute value of the entropy gradient (namely 8/7), the single evanescent wave has a much larger amplitude in the case ${\cal S}>0$ (i.e., for the plot f), almost of a factor of 2, as can be seen from the distance between the dotted and dashed curves. This is compatible with the fact that in Figure 5, we found that the corrected torque excess, in absolute value, is almost a factor of 2 larger at ${\cal S}=8/7$ than at ${\cal S}=-8/7$.

Zoom In Zoom Out Reset image size

We first establish the scaling of the torque excess as a function of the disk's and planet's parameters. Equation (55) shows that the torque excess involves the half-width of the horseshoe region and the potential at the stagnation point. In order to remain in the framework of a separable expression, we neglect the dependence of x_s on the entropy gradient, mentioned in Section 4.2, and we use the barotropic approximation (i.e., for ${\cal S}=0$) to infer the half-width. Using Equation (57), we obtain

$$\begin{equation} x_s^{\rm adi}=\frac{C(\hat{\epsilon })}{\gamma ^{1/4}}a\sqrt{\frac{q}{h}}, \end{equation} \tag{ 91 }$$

where $C(\hat{\epsilon })$ is a dimensionless constant that depends on the softening length (Paardekooper & Papaloizou 2009b). Examination of Figure 3 shows that this approximation leads to errors of at most ∼30%. We check in what follows that the potential at the stagnation point can also be considered as essentially independent of the entropy gradient.

7.1.1. Potential at the Stagnation Point

Evaluating the planetary potential at the stagnation point requires to know the location of this point. Writing the Euler equations in a steady state, in a system of units in which the unit of length is H = c_s/Ω, one can realize that, as long as the planet is deeply embedded (i.e., aq^1/3 ≪ H), the azimuth of the stagnation point (which lies almost on corotation) has to scale as $hf(\hat{\epsilon })$, where f is a function that we shall determine. We therefore have

$$\begin{equation} \Phi _p ^s=-\frac{GM_p}{H[f(\hat{\epsilon })^2+\hat{\epsilon }^2]^{1/2}}. \end{equation} \tag{ 92 }$$

Figure 13 shows the location of the stagnation point as a function of $\hat{\epsilon }$ in an adiabatic flow with γ = 1.4, obtained from the series of numerical simulations EA_0⩽i⩽20 (see Section 4.1) and a similar series with a different entropy gradient. Under $\hat{\epsilon }\simeq 1.0$, the distance of the stagnation point to the planet is larger than the softening length, and the denominator of Equation (92) is dominated by the function $f(\hat{\epsilon })$, which is found to hardly depend on ${\cal S}$, as both series essentially coincide for $\hat{\epsilon }< 1.0$. Reciprocally, for $\hat{\epsilon }>1.0$, the stagnation point lies within a softening length from the planet, and the potential at the stagnation point approximately scales as $\hat{\epsilon }^{-1}$, and is therefore also independent of the entropy gradient.

Zoom In Zoom Out Reset image size

7.1.2. General Expression of the Torque Excess

Using Equations (55), (91) and (92), we obtain

$$\begin{equation} \Gamma _1=-A(\hat{\epsilon }){\cal S}\Gamma _0, \end{equation} \tag{ 93 }$$

where Γ₀ = Σ₀Ω²_pq²a⁴/h² and

$$\begin{equation} A(\hat{\epsilon }) = \frac{4C^2(\hat{\epsilon })}{\sqrt{1.4}}\left[\frac{1}{\sqrt{f(\hat{\epsilon })^2+\hat{\epsilon }^2}}-\frac{3}{8}\frac{C^2(\hat{\epsilon })}{\sqrt{1.4}}\right], \end{equation} \tag{ 94 }$$

where we have specialized to the Keplerian case (A = −3Ω/4, B = Ω/4). Note that we explicitly restrict ourselves to the case γ = 1.4. Although Equation (91) applies for arbitrary values of γ, the azimuth of the stagnation point, considered in Section 7.1.1, may depend on the adiabatic index. Keeping an explicit dependence in γ in Equation (94) could therefore be misleading. The adiabatic torque excess, given by Equation (93), has a scaling similar to the barotropic corotation torque or the differential Lindblad torque (Tanaka et al. 2002). The only differences with the barotropic corotation torque is that the excess scales with ${\cal S}$ rather than ${\cal V}$ and that it involves a different dimensionless function of $\hat{\epsilon }$.

In order to determine this function, we use the series SA_i, as explained above. We restrict ourselves to the range $0.1\le \hat{\epsilon }\le 1$. Values of $\hat{\epsilon }$ larger than 1 are not relevant, whereas values smaller than 0.1 are difficult to study owing to the very high resolution required. As mentioned in Section 4.3, the torque excess exhibits a dependence on $\hat{\epsilon }^{-1}$ with a good approximation over the range of $\hat{\epsilon }$ considered (a dependence partly due to the increase of x_s at small softening length and partly due to the increase of the planetary potential at the stagnation point at small softening length). Fitting the torque excess over the range of $\hat{\epsilon }$ considered, we obtain

$$\begin{equation} A(\hat{\epsilon }) \approx \frac{1.3}{\hat{\epsilon }}. \end{equation} \tag{ 95 }$$

Therefore, the expression of the torque excess is

$$\begin{equation} \Gamma _1=-\frac{1.3{\cal S}}{\hat{\epsilon }}\Sigma _0\Omega _p^2q^2a^4h^{-2}. \end{equation} \tag{ 96 }$$

One can see in Figures 4 and 5 that this expression (which is represented by the gray dashed line) is in good agreement with the global adiabatic excess. It is also in good agreement with the results of Baruteau & Masset (2008).

If one considers that a two-dimensional case with softening length represents the layer at altitude z = ± of a three-dimensional case (see, e.g., Masset 2002), one can perform an integration of the torque over in order to get an approximate value of the three-dimensional torque. A layer comprised between altitude z and z + dz has a surface density $(\Sigma _0/\sqrt{2\pi }H)\exp [-z^2/(2H^2)]dz$; hence, the three-dimensional adiabatic torque excess should be approximately given by

$$\begin{eqnarray} \Gamma _{3D}&\approx & -2\int _0^{+\infty }\frac{\Sigma _0}{\sqrt{2\pi }H}\exp [-z^2/(2H^2)] A\left(\frac{z}{H}\right)\nonumber \\ && \times{\cal S}\Omega _p^2\left(\frac{q}{h}\right)^2r_p^4dz\nonumber \\ &\approx &-\sqrt{\frac{2}{\pi }}{\cal S}\Gamma _0\int _0^{+\infty }e^{-\hat{\epsilon }^2/2}A(\hat{\epsilon })d\hat{\epsilon }. \end{eqnarray} \tag{ 97 }$$

While this kind of calculation can be performed relatively easily for the vortensity-related torque (for which the dimensionless function of $\hat{\epsilon }$ converges when $\hat{\epsilon }\rightarrow 0$) or for the differential Lindblad torque, it is a risky exercise in the case of the adiabatic torque excess, since the approximate expression of Equation (95) diverges for $\hat{\epsilon }\rightarrow 0$. For the setup that we consider, we have tried to perform calculations at higher resolution and smaller smoothing lengths. The smallest smoothing length for which we have a steady situation is = 0.07H (corresponding to the left square in Figure 13). Under this value, we end up with very messy, unsteady situations, in which the flow is strongly affected by the presence of numerous vortices that drift along the separatrices. We can therefore only provide a conservative estimate of Equation (97), by truncating the integral at $\hat{\epsilon }=0.07$. We obtain

$$\begin{equation} \Gamma _1^{3D}\approx -2.8{\cal S}\Gamma _0. \end{equation} \tag{ 98 }$$

In order to estimate the total torque in the three-dimensional case, we need to know the differential Lindblad torque and the vortensity-related part of the horseshoe drag. The first of those is given by Tanaka et al. (2002), and is, with our notation

$$\begin{equation} \Gamma _{LR}=-(2.34-0.1\alpha)\Gamma _0, \end{equation} \tag{ 99 }$$

for a disk without a temperature gradient. We note that we do not use the total torque value given by Tanaka et al. (2002), but only the differential Lindblad torque. The total torque value indeed contains the linear corotation torque, which is not the corotation torque exerted in a steady state (Paardekooper & Papaloizou 2009a). Instead, we consider separately the vortensity-related horseshoe drag and apply to it a treatment similar to that of Equation (97), which is straightforward in this case, because of the absence of divergence when $\hat{\epsilon }\rightarrow 0$. The two-dimensional component of the horseshoe drag that scales with the vortensity gradient reads (Ward 1991; Masset 2001)

$$\begin{equation} \Gamma _V=\frac{3}{4}{\cal V}C^4(\hat{\epsilon })\Gamma _0; \end{equation} \tag{ 100 }$$

hence, its three-dimensional estimate is

$$\begin{equation} \Gamma _V^{3D}\approx \frac{3}{\sqrt{8\pi }}{\cal V}\Gamma _0\int _0^{+\infty }e^{-\hat{\epsilon }^2/2} C^4(\hat{\epsilon })d\hat{\epsilon }\approx 0.93{\cal V}\Gamma _0, \end{equation} \tag{ 101 }$$

where we have used the series SI_i to tabulate $C(\hat{\epsilon })$. Therefore, the total torque acting on a low-mass planet in a three-dimensional disk reads

$$\begin{equation} \Gamma =[-(2.34-0.1\alpha)+0.93{\cal V}-2.8{\cal S}]\Gamma _0, \end{equation} \tag{ 102 }$$

where we recall that the last term is a conservative estimate, and is likely to actually have a numerical coefficient larger than 2.8. We see that for $|{\cal S}|\sim |{\cal V}|\sim 1$, the term associated with the entropy gradient dominates the horseshoe drag; hence, the adiabatic torque excess should efficiently halt or reverse migration in disks with moderate, negative entropy gradients, in agreement with the result of Paardekooper & Mellema (2006), who found that planetary migration can indeed be reversed in a three-dimensional disk.

As stated in Section 2.4, we have assumed in this work that there is only one X-point in the vicinity of the planet. We find that this is indeed the case in most situations. We can nevertheless get two X-points in the vicinity of the planet, if one of the following conditions is fulfilled.

• 1.  
  The entropy gradient vanishes, or is very small.
• 2.  
  The softening length of the potential is very small.
• 3.  
  The flow has not reached a steady state and is observed before a U-turn timescale.

While the first and last conditions are not relevant to the present work, which assumes a steady state and a sizable entropy gradient, the second case has some importance, as it gives indications of what can happen near the equatorial plane in a three-dimensional case. In the cases in which we have two X-points that subsist during the whole simulation, we get very time-dependent torque values. This can be interpreted as entropy trapped in the small libration island defined by the two X-points. The perturbation of density associated with the entropy distribution, which librates in this island, gives rise to large torque variation, as it corresponds to material very close to the planet. We note that this is quite in contrast with what we have claimed in Section 5.3, where we have argued that it is the vortensity perturbation that yields the bulk torque. Whereas this is true for the downstream horseshoe stripes which lie at some distance of the planet, here the planet lies within the libration island, and the fact that the density response is spread radially or not definitely matters. We also note that the present discussion offers a number of similarities with the discussion in Paper I about the flow topology as a function of the temperature gradient and of the softening length, except that here it is the entropy gradient that plays a role rather than the temperature gradient.

When the flow is observed at an early stage, it may display two X-stagnation points, which are generally not on the same separatrix. We observe that once the flow has settled in a steady state, it is the outermost stagnation point (the one that is on the separatrix that lies further from corotation) that subsists. It turns out that this point lies behind the planet for ${\cal S}<0$ and is in front of the planet for ${\cal S}>0$. While this fact is unimportant for a planet in a fixed circular orbit, where only the distance of the stagnation point to the planet matters for evaluating the torque excess, it is crucial when the planet migrates, as the stagnation point is shifted from the position that it has on a fixed circular orbit. The consequences of this are discussed in the following section.

Since the adiabatic torque excess depends on the flow properties at the stagnation point of the horseshoe region, and since, as we shall see below, the location of the latter depends on the drift rate, we expect that a feedback loop can be established between the total torque and the migration rate, much as in type III migration (Masset & Papaloizou 2003). In contrast to type III migration, however, the feedback should generally be negative, therefore acting to decrease the absolute value of the drift rate. It should also be of relatively minor importance, except in very massive disks, with Toomre's Q parameters close to unity. Consider the following numerical experiment: we repeat the runs SA₁₀ and SI₁₀, except that we impose an inward disk drift $\dot{r}_d=-1.9\times 10^{-5}$. We note that, if the planet was allowed to freely migrate in run SA₁₀, it would do so at the rate $\dot{a}=1.9\times 10^{-5}$. Rather than releasing the planet, we impose an inward drift of the disk, of the same magnitude. This technique has already been used for type III migration (see Masset & Papaloizou 2003, Section 5.6). We see in Figure 14 that the torque in the isothermal case does not depend on the drift rate, while there is a noticeable difference in the adiabatic case, and that the torque tends toward a somehow smaller value. The reason for this can be understood from the profile of radial velocity at corotation, depicted in Figure 15. If the planet migrates outward, the stagnation point (in a frame that moves radially with the planet) is located where the horseshoe dynamics endows fluid elements with a radial velocity equal to that of the planet (similarly, if the planet is kept fixed, the new stagnation point is expected to lie where the radial velocity in the initial case is $-\dot{r}_d$.) Since fluid elements move outward behind the stagnation point, the azimuth of the stagnation point decreases; hence the latter recedes from the planet, since it is located behind the planet in the case considered. A consequence of this recession is a lower perturbed pressure at the new position of the stagnation point, hence a lower value of the adiabatic torque excess. A similar decrease (in absolute value) is also expected in the case ${\cal S}>0$. The stagnation point is then in the front of the planet (see Section 8.1), and the total torque acting on the latter is negative. The stagnation point in the migrating case is therefore expected where the horseshoe U-turns are performed inward, that is to say at the front of the stagnation point of the fixed case, and therefore further again from the planet. We have checked this with additional simulations, not reproduced here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A number of comments are in order as follows.

• 1.  
  Figure 15 shows a poor agreement between the radial velocity for the disk's drift case and the original radial velocity offset by $-\dot{r}_d$, in the region of the stagnation point (−0.07 ⩽ ϕ ⩽ −0.03). This is presumably due, at least in part, to the resolution of the grid. The direction of the shift of the stagnation point and its order of magnitude are nevertheless compatible with the new stagnation point being determined by the location where, in the initial flow, $v_r=\dot{a}$.
• 2.  
  Since both v_r and $\dot{a}$ scale with the planet mass (for the range of planetary masses considered here), the shift of the stagnation point should not depend on the planetary mass.
• 3.  
  Since $\dot{a}$ scales with the disk mass, the shift of the stagnation point of a freely migrating planet should increase with the disk mass. The magnitude of this effect, for the disk considered in the numerical experiment described above, is at about 1/20th of the total torque. The disk has a Toomre parameter Q ≈ 8. This shows that the effect we describe is of marginal importance, except in the most massive disks, for which Q ≳ 1.
• 4.  
  The drift rate in the steady state is given by

  $$\begin{equation} 2Ba\dot{a}=\gamma +\frac{\partial \gamma }{\partial \dot{a}}\dot{a}, \end{equation} \tag{ 103 }$$

  where we assume that the torque excess, to lowest order, has an affine dependence on $\dot{a}$ and γ is the specific torque acting on the planet. We therefore obtain

  $$\begin{equation} \dot{a}=\frac{\gamma }{2Ba-\partial \gamma /\partial \dot{a}}. \end{equation} \tag{ 104 }$$

  In the particular case described above, we have $\partial \gamma /\partial \dot{a} \approx -0.05 \cdot 2Ba$; hence, the steady-state drift is almost equal to that dictated by the torque value measured in the disk drift case (i.e., further iterations with a drift rate of the disk given by the new torque value would hardly change the result).

The analysis presented in this work suffers from a restriction very similar to the restriction of the isothermal case, which is that the temperature profile is assumed to be flat. We have nevertheless performed many additional calculations in which we relaxed this constraint, and found that our main results essentially hold in disks with arbitrary temperature profiles.

• 1.  
  The global torque excess, defined by Equation (57), scales with the entropy gradient. The isothermal runs which we use to apply Equation (57) in the general case have the same profiles as their adiabatic counterparts; hence they are locally isothermal runs, which have singular evanescent waves excited at the separatrices, in contrast to the globally isothermal runs (see Paper I). These components, which scale with the temperature gradient, should therefore be also present in the general, adiabatic case.
• 2.  
  The half-width of the horseshoe region exhibits the same characteristic V-shaped dependence on ${\cal S}$, as shown in Figure 3.

The adiabatic torque excess, as a function of the entropy gradient, displays little scatter, and is therefore essentially a one-to-one correspondence. The results presented here can therefore be used to predict the torque in a general situation. They indicate that, in the general case, the horseshoe drag consists of three terms: the classical vortensity-related horseshoe drag, the adiabatic torque excess presented in this work, and an additional term which scales with the temperature gradient (see Paper I). The latter corresponds to the creation of vortensity in the vicinity of the planet, near the stagnation point, where the flow is slow and the driving of vortensity is efficient. We comment that the adiabatic torque excess can also be seen as a singular creation of vortensity at the downstream separatrices, with the word of caution that vortensity is mathematically ill-defined at these separatrices, because of the presence of a contact discontinuity (so that one has to consider the low order estimate ω/Σ₀ rather than ω/Σ).

The analysis presented above relies upon the existence of a steady state and on that of a unique stagnation point. It therefore does not apply to the case of a turbulent disk. The adiabatic excess has been interpreted, however, as due to a constitutive asymmetry of the horseshoe region, acquired by the upstream flow under the action of the pressure gradient on overdense or underdense regions. To some extent, this effect should therefore persist in a turbulent disk, since we expect to find, in a turbulent disk, regions that are, in average, underdense or overdense with respect to the barotropic case, as a consequence of the advection of entropy. This issue deserves significant further work, since most of the regions of the disks where planetary migration is supposed to take place should be turbulent. In particular, it would be of interest to know whether the time-averaged torque excess is equal to the steady-state estimate.

In the absence of any process that can transfer angular momentum between the horseshoe region and the rest of the disk, the horseshoe drag is bound to saturate. A nonvanishing time-averaged value of this torque would indeed imply a sustained transfer of angular momentum from the planet to the horseshoe region, where it would accumulate. Viscous torques, exerted at the separatrices between the horseshoe region and the outer or inner disk, can ensure this transfer of angular momentum (Masset 2001). In the case of an isothermal disk, viscous diffusion tends to restore the large-scale gradient of vortensity across the horseshoe region. If it can do so in less than a libration time, the vortensity upstream of the horseshoe U-turn is that of the unperturbed disk, and the horseshoe drag is permanently equal to its initial value. In the case of an adiabatic flow, the vortensity-related part of the horseshoe drag follows the same picture, and a sufficiently large viscosity is required to avoid the saturation of this component of the drag. The adiabatic torque excess, however, results from the advection of entropy and from the existence, at the stagnation point, of a discontinuity of the entropy brought from the inside and from the outside. Sustaining the torque excess therefore implies that the discontinuity of entropy at the stagnation point is maintained. Paardekooper & Papaloizou (2008) considered thermal diffusion as the dissipative process that fulfills this function, while Kley & Crida (2008) considered realistic radiative effects. In both cases there can be a sustained, positive torque acting on the planet. We however note the following facts.

• 1.  
  In any case, a finite amount of viscosity is required. Thermal diffusion or radiative effects, which do not feature in Euler's equation, cannot ensure the transfer of angular momentum out of the horseshoe region.
• 2.  
  There is no theoretical expression of the steady-state horseshoe drag that takes into account the balance between saturation and the dissipative processes that prevent it. Such an expression is definitely required for planet population synthesis based on one-dimensional models of migration.

We also comment that while saturation can be interpreted as the phase mixing of the contribution of different horseshoe streamlines in the barotropic case (Balmforth & Korycansky 2001), the case of the adiabatic torque excess is slightly different because there is essentially one streamline involved in this excess: the separatrix. It is therefore possible that the torque excess oscillates over a very long time or indefinitely. Numerical simulations show that the oscillations of the torque damp after a few libration times (Baruteau 2008). It is not clear, however, whether this is a physical effect or the consequence of numerical diffusivity. This again appeals for the need of simulations at very high resolution.

We finally mention the systematic appearance of a vortex at one of the downstream separatrices (the one that has a negative singularity of vorticity). It would be of interest to study the impact of this vortex on the saturation of the adiabatic excess. On the one hand, it provides a mixing of the entropy in the vicinity of the separatrix and, on the other hand, it may assist the horseshoe region in exchanging angular momentum with the rest of the disk.