Improved torque formula for low and intermediate mass planetary migration
María Alejandra Jiménez Frédéric S. Masset
astro-ph.EP


firstpage–lastpage
Binary Hermitian Forms and Optimal Embeddings
    Michael Zhao
    December 30, 2023
=============================================





  The migration of planets on nearly circular, non-inclined orbits in
  protoplanetary discs is entirely described by the disc's
  torque. This torque is a complex function of the disc parameters,
  and essentially amounts to the sum of two components: the Lindblad
  torque and the corotation torque. Known torque formulae do not
  reproduce accurately the torque actually experienced in numerical
  simulations by low- and intermediate-mass planets in radiative
  discs. One of the main reasons for this inaccuracy is that these
  formulae have been worked out in two-dimensional analyses. Here we
  revisit the torque formula and update many of its dimensionless
  coefficients by means of tailored, three-dimensional numerical
  simulations. In particular, we derive the dependence of the Lindblad
  torque on the temperature gradient, the dependence of the corotation
  torque on the radial entropy gradient (and work out a suitable
  expression of this gradient in a three-dimensional disc). We also
  work out the dependence of the corotation torque on the radial
  temperature gradient, overlooked so far. Corotation torques are
  known to scale very steeply with the width of the horseshoe
  region. We extend the expression of this width to the domain of
  intermediate mass planets, so that our updated torque formula
  remains valid for planets up to typically several tens of Earth
  masses, provided these relatively massive planets do not
  significantly deplete their coorbital region. Our torque expression
  can be applied to low- and intermediate-mass planets in optically
  thick protoplanetary discs, as well as protomoons embedded in
  circumplanetary discs.





Planetary systems: formation – Planetary systems: protoplanetary discs – Accretion, 
accretion discs – Methods: numerical – Hydrodynamics – Planet-disc interactions 





§ INTRODUCTION
 


An important ingredient for studies of planetary population synthesis
are the so-called migration maps, which are two-dimensional functions
Γ(M_p,r) that provide the torque exerted by the disc on a
planet on a nearly circular orbit, as a function of its mass M_p and
of its distance r to the central star.  A migration map is intrinsic
to a given disc model. Over most of the (M_p,r) domain, the torque
is generally a negative quantity, corresponding to an orbital decay of
the planet. Nevertheless, there may be some regions where the torque
is positive (named islands of outward migration). These regions
broadly occur where the disc's temperature drops faster than r^-1,
and span a mass range from a few Earth masses to potentially several
tens of Earth masses. The outcomes of models of planet population
synthesis depend sensitively on the outline of these islands
<cit.>. Migration maps can be established
in one of two ways: they can either be obtained by three dimensional
(3D) simulations, including all the physical ingredients needed for a
correct description of the disc <cit.>,
or they can be obtained using torque formulae. The first approach is
naturally not suitable to the exploration of a large number of disc
models, owing to its considerable computational cost, and one must
resort to torque formulae <cit.>. Numerical
simulations of low-mass planets embedded in radiative discs allows to
assess the validity of torque formulae. <cit.>
have compared the torque experienced by a 20 M_⊕ planet in a
radiative disc to the predictions of <cit.> and
<cit.>. The former displays a broad agreement
with the outcome of the simulations. The latter, however, was found to
be at odds with the simulations outcome, a discrepancy identified as
too saturated corotation torques, and later resolved when estimates of
the width of the horseshoe region in 3D discs became available
<cit.>. Despite their broad agreement with
simulation outcomes, the torque formulae of <cit.> and
<cit.> lack the accuracy required by models of
planetary population synthesis. One of the main reasons is that these
analyses are based on results of two-dimensional analysis. Even though
they capture most of the mechanisms that contribute to the torque, the
dimensionless coefficient that they feature in many places may be
significantly off.

The purpose of this paper is to update the torque formula of
<cit.> by means of 3D numerical simulations,
successively tailored to determine the value of each of the different
dimensionless coefficients that we want to update. It is organised as
follows. In section <ref> we introduce and
discuss the different torque components that we want to update and we
describe the numerical code and setups that we use for this
purpose. In section <ref> we present our different
results, and provide an updated torque formula in
section <ref>. For the convenience of the
reader mainly interested in the formula, this section is
self-contained.  We compare the results of our updated torque formula
to published numerical simulations in section <ref>,
and draw our conclusions in section <ref>.



§ NOTATION


We consider a planet of mass M_p on a circular orbit of radius r_p
around a central star of mass M_⋆, with angular speed
Ω_p, embedded in a disc of surface density Σ(r) which
follows the power law:

    Σ(r)=Σ_0(r/r_p)^-α.

The disc's midplane temperature T(r) obeys the law:

    T(r)=T_0(r/r_p)^-β.

We denote with ν the kinematic viscosity of the disc, with χ
its thermal diffusivity and with κ its opacity.  The disc's
pressure scale length is H, and is given by:

    H=c_s/Ω_p,

where c_s is the disc's isothermal sound speed:

    c_s=√(%s/%s) RTμ,

R being the constant of ideal gases and μ the mean
molecular weight. We will also make use of the disc's aspect ratio

    h=H/r,

and of the planet-to-star mass ratio

    q=M_p/M_⋆.




§ FOCUS OF OUR TORQUE UPDATE


Both the Lindblad and corotation
torques can be normalised to:

    Γ_0=Σ_0Ω_p^2r_p^4(q/h)^2=Σ_0Ω_p^4r_p^6q^2c_s^-2.

We draw hereafter a list of the different updates that we perform on
the normalised torques, which are dimensionless quantities.


 §.§ Lindblad torque


The normalised Lindblad torque depends on α and β, which
are respectively the slopes of surface density and temperature
<cit.>. Its dependence on α has
already been worked out by <cit.> by means of
semi-analytical calculations in the linear regime for planets in
globally isothermal discs, and we do not revise here this (weak)
dependence. On the contrary, its dependence on the temperature slope
has been worked out by <cit.> by a linear
analysis in two-dimensional discs with a softened potential.

Here we reevaluate the dependence of the Lindblad torque on the
gradient of midplane temperature. We undertake this analysis in
section <ref>. We note that the Lindblad torque
can exhibit a dependence not only on α and β, but also on
αβ <cit.>. This dependence is very
weak, however, and we neglect it. More generally, we seek, and
restrict ourselves to, linear dependencies of the different torque
components on the radial gradients of physical quantities at the
planet's location.

The Lindblad torque scales with the inverse squared of the sound speed
(see Eq. <ref>). It therefore depends on the thermal diffusivity
of the disc: when the latter is large enough that the disturbances
triggered by the planet at Lindblad resonances behave isothermally, it
scales with the inverse squared of the isothermal sound speed. On the
other hand, when the thermal diffusivity is small and the disturbances
behave adiabatically, the Lindblad torque scales with the inverse
squared of the adiabatic sound speed, and is therefore a factor of
γ smaller, where γ is the disc's adiabatic index. This
dependence has been considered by <cit.> and
<cit.>. We do not revise it here, and we will adopt for the
updated expression of the Lindblad torque (see
section <ref>) the dependence given by
<cit.>.



 §.§ Corotation torque


The corotation torque is the sum of three main components: one which
scales with the radial gradient of vortensity, one which scales with
the radial gradient of entropy, and one which scales with the radial
gradient of temperature. In addition, there is a contribution to this
torque arising from the vortensity viscously created at the contact
discontinuities which appear on the downstream separatrices in the
presence of an entropy gradient. These four contributions are
represented in Fig. <ref>. Each of the first three
components is a blend of a linear component and of a horseshoe
drag. Each horseshoe drag is itself the product of the
unsaturated horseshoe drag by a saturation function which is
generally a number between 0 and 1 (corresponding respectively to
a totally saturated and to an unsaturated horseshoe drag).  Each
torque component can in principle have its own saturation function.




We neither update, in this work, the saturation functions nor the
blending coefficients. <cit.> gave an expression
of the former, obtained from first principles by a two-dimensional
analysis of a simplified model of the horseshoe flow. The horseshoe
flow has since been shown to be essentially two-dimensional even in 3D
discs <cit.>. On another hand, the flow arising
from turbulence in three dimensions and that arising from a kinematic
viscosity in a laminar 3D disc are markedly different
<cit.>, so a 3D study of the torque saturation as
a function of kinematic viscosity would probably be pushing the
viscous model of protoplanetary discs beyond its domain of validity.

The blending coefficients were determined through a fit of numerical
experiments by <cit.> and by
<cit.>. <cit.> use the weights
ε and 1-ε respectively for the horseshoe drag
and for the linear corotation torque (where ε depends on
the torque component). <cit.> use the weights 1-K and G, so
that their sum is not necessarily 100 %. This allows them to
reproduce the fact that the horseshoe drag may be larger than its
unsaturated value. In the analysis of <cit.>,
this fact is incorporated at a different place (namely the fact that
the saturation functions can have values slightly larger than
unity). Here we stick to the blending coefficients of
<cit.>.

The vortensity component of the corotation torque does not need any
update. <cit.> have shown that the 3D horseshoe
drag in a barotropic disc has the same expression in a two dimensional
disc, and that the horseshoe region has a vertical separatrix, hence a
constant width on the whole vertical extent of the disc <cit.>. Naturally, an updated width of the
horseshoe region is used to evaluate this torque (more detail on this
is given in section <ref>), but no amendment has
to be done to the horseshoe drag expression. For the sake of
simplicity, we consider here that the gradient of
vortensity[More specifically, the quantity
  dlog(Σ/B)/dlog r, where B, the second Oort's constant, is
  half the flow's vertical vorticity.] is equal to 3/2-α. This
is true in any part of the disc where the surface density and angular
frequency are power laws of the distance to the star.

The temperature component of the corotation torque has been overlooked
in previous torque formulae. Here we obtain the value of this torque
as a by-product of the study of the Lindblad torque of
section <ref> and present its properties in
section <ref>. This torque component arises from
the production of vortensity by the radial temperature gradient, and
appears in locally isothermal calculations. The vortensity thus
produced is concentrated near the (downstream) separatrices. Although
it is not singular as the vortensity arising from non-barotropic
effects, which can formally be represented by a delta-function at the
separatrix <cit.>, it has nearly same spatial
distribution and it is reasonable to expect it to saturate in a
similar manner. We therefore assume that this torque has same
dependence on the viscosity as the entropy torque, and focus here
exclusively on the magnitude of the unsaturated horseshoe drag and
that of the linear torque.

Finally, we reevaluate the dependence of the entropy component on the
entropy gradient, a much needed determination since this component
plays a preponderant role in the appearance of outward migration
islands in migration maps. We also reevaluate how to determine the
entropy gradient in 3D disc. Whereas its expression is straightforward
in 2D discs, one may contemplate different ways of evaluating it in a
3D disc (such as the gradient in the midplane, or of the vertically
averaged quantity, etc.). Here we simply adopt the expression that
yields to a one-to-one relationship (with the least possible
dispersion) of the torque excess in a non-barotropic simulation with
respect to a barotropic simulation. This analysis is presented in
section <ref>.



 §.§ Additional analysis


In addition to the aforementioned reevaluations of key dependencies of
the torque, we have performed two side studies in order to improve the
accuracy and predictive power of the torque formula, which we present
below.



  §.§.§ Width of horseshoe region


The different components of the horseshoe drag scale with the width of
the horseshoe region to the fourth power. It is therefore crucial to
get an accurate estimate of this width. Besides, as soon as the ratio
q/h^3 is larger than ∼ 0.1, the law that gives the width for
small mass planets
<cit.> is no longer
valid and the horseshoe region is actually larger than what is
predicted by the low-mass law, resulting in a boost of the corotation
torque <cit.>. The
threshold of 0.1 is particularly stringent for typical
protoplanetary discs, and translates into 4 M_⊕ only for a
disc with an aspect ratio of h=0.05 (which shows, incidentally, that
the comparisons mentioned in section <ref> were
performed much beyond the domain of validity of the torque formulae
for low-mass planets). In order to relax such a stringent condition on
the planetary mass, we study here the transition between the low-mass
regime, for which the width scales as (q/h)^1/2, and the high-mass
regime, for which, as in the restricted three body problem, the width
scales as q^1/3. This study is presented in
section <ref>.



  §.§.§ Thermal diffusivity


Since the saturation of the corotation torque depends on the value of
the disc's thermal diffusivity, it is important to determine this
value as accurately as possible. We have studied the radial spread of
an initially localised excess of temperature in order to assess the
accuracy of standard estimates which relate the disc's thermal
diffusivity to the disc's temperature, density and opacity. This study
is presented in section <ref>.




 §.§ Numerical details



We conduct all our numerical experiments with the public hydrocode
FARGO3D[<http://fargo.in2p3.fr>]
<cit.> with orbital advection enabled
<cit.>. Here we use a spherical mesh that covers the full
azimuthal range [-π,π]. We denote N_ϕ the number of cells
in azimuth, N_r the number of cells in radius and N_θ the
number of cells in colatitude. We always simulate a half-disc in
colatitude (since the planet is always coplanar with the disc), and
use reflecting boundary conditions at the midplane. All our
calculations involving a planet are performed in the corotating frame.

According to our needs, we either use a locally isothermal equation of
state, or we solve an energy equation. For the first case we use the
closure relationship between the pressure p and density ρ:

    p( r,t)=c_s^2(r)ρ( r,t),

where c_s(r) is the isothermal sound speed, which depends on the
distance to the central star, and which is constant in time. In such
case there is no need to solve an energy equation. In the second case,
we solve the time evolution of the internal equation, and in the
absence of some form of thermal diffusion, the flow behaves
adiabatically. In this case we assume the gas to be ideal, so that the
pressure is given by:

    p( r, t)=(γ-1)e( r,t),

where e is the internal energy density.

Finally, in order to avoid the divergence of the planetary potential Φ_p
in the vicinity of the planet, we soften it over a length scale
ϵ:

    Φ_p(r)=-GM_p/√(|r-r_p|^2+ϵ^2).

In three-dimensional calculations, ϵ has to be chosen small
compared to the pressure scale length, and larger or comparable to the
resolution. In all our numerical experiments we have ϵ=0.1H.



§ RESULTS


We present hereafter the results of the numerical experiments
conducted to address the points listed in
section <ref>.



 §.§ Dependence of the Lindblad torque on the temperature
  gradient


For this experiment, we use a locally isothermal setup, and use the
fact that at larger time the corotation torque saturates, so that the
total torque essentially amounts to the Lindblad torque.

In a locally isothermal setup, the Lindblad torque has
the form:

    Γ_L/Γ_0=-(2.34-0.099α+kβ),

where as discussed in section <ref> the first two
terms of the right hand side come from the work of <cit.>,
and where k is the coefficient that we want to determine. In
principle, it suffices to perform two runs which have different values
of β (all other parameters being the same) to infer k.

Our choice of parameters results from a trade-off between different
effects. We want the planetary mass to be sufficiently large so that
the horseshoe region is well resolved and can be saturated in a
reasonable amount of time. On the other hand, we want the planetary
mass to be sufficiently small so that the planet does not
significantly perturb the disc by carving a gap.

There is another difficulty inherent to this experiment: when the disc
is not globally isothermal (β 0), it is prone to a baroclinic
instability <cit.> which results in considerable
noise in the torque curves, precluding any converged
measurement. However interesting this instability may be, it has for
us a parasitic character, and we aim at getting rid of it, while
preserving an almost complete saturation of the corotation torque. We
prevent the appearance of the instability by using some amount of
viscosity, and determine by dichotomy the minimal amount of viscosity
required to quench the instability. Also, the strength of the
instability scales with |β|, so one has to choose values of
β sufficiently large to allow an accurate measurement of k,
but also sufficiently small to avoid the development of the baroclinic
instability, once the viscosity has been chosen.

With these considerations in mind, we have chosen for this numerical
experiment the following values: q=2.4· 10^-5 (which
corresponds to 8 M_⊕ if the central star has a solar mass) and
ν=7· 10^-8r_p^2Ω_p. The disc's aspect ratio is
h=0.05. Although the kinematic viscosity is quite small, we do not
expect the planet to carve a significant gap in the disc over the time
scale required to achieve the saturation of the corotation torque. The
Hill radius of the planet [r_p(q/3)^1/3] is much smaller than the
pressure length scale, and the expression
3h/[4(q/3)^1/3]+50ν/(r_p^2Ω_p q) is a factor of 2 above
the critical value for gap opening <cit.>. Even under these conditions a gap may eventually be
opened
<cit.>,
but this will occur on time scales much longer than those of our
numerical experiments. As we shall see below, we do find evidence for
secular effects in our torque measurements, which affect our results
at the percent level.

  Our mesh extends from π/2-3h to π/2 in colatitude, over
  N_θ=40 cells, and from 0.6r_p to 1.4r_p in radius, over
  N_r=200 cells. In addition, we have N_ϕ=900 cells. We use
  wave-killing boundary conditions in radius <cit.>, in
  order to avoid reflection of the wake at the radial boundaries. We
  therefore expect our results to be nearly insensitive to the exact
  location of the mesh radial boundaries, as these are located at a
  large number of pressure scale lengths from the orbit. We
  perform two simulations, with β=1/2 and β=-1/2
  respectively. Each one is run over 1000 orbital periods of the
  planet. We choose α=3/2 so that there is no vortensity
  gradient, and the only corotation torque present before saturation
  is the temperature component of the corotation torque.

This choice of parameters implies that the horseshoe region is
resolved over 11 zones, which is sufficient to allow for a nearly
complete saturation of the torque <cit.>. We have
q/h^3=0.19, so that the flow is weakly non-linear. The horseshoe
libration time is ∼ 60 orbits, so that our runs duration allows
for a full saturation of the corotation torque. The ratio of the
viscous time scale across the horseshoe region to the libration time
scale is ∼ 7, which implies that the corotation should reach a
degree of saturation of ∼ 90 %. We show the torque as a function
of time for the two runs in Fig. <ref>. The corotation
torque saturates, leaving essentially the Lindblad torque. We see a
systematic, slow increase of the torque in both cases, which we
attribute to the slow opening of a very shallow dip (we find that the
rate of this increase decays as the viscosity increases).





The value of k can be directly obtained by subtracting the two
normalised torques. We show their difference in
Fig. <ref>. From this figure we infer a value of
k≈ 1.4, comparable to, and marginally smaller than the value
of 1.7 found by <cit.> for 2D discs.



 §.§ Dependence of the corotation torque on the temperature
  gradient


We can further exploit the runs presented at the previous section to
work out the dependence on β of Γ _T^UHD and
Γ _T^lin (see Fig. <ref>).
Over the first half libration time, the torque is the sum of the
Lindblad torque and the unsaturated horseshoe drag. Since there is no
vortensity gradient in these runs, the horseshoe drag is entirely
attributable to the temperature gradient. As we now have an expression
of the Lindblad torque that takes the temperature gradient into
account, we can subtract this torque from the torque measured to get
an estimate of the corotation torque. We can also slightly improve
upon the results of the previous section by taking into account the
residual value of the corotation torque at larger time, as it does not
fully saturate.  We write the torque over the first half libration
time as:

    Γ^(U)_β=Γ_L+Γ _T^UHD,

where the (U) superscript stands for “unsaturated”, and where the
β subscript conveys that in this numerical experiment the only
parameter that is varied is β. In Eq. (<ref>), we do not
have a vortensity component for the reasons exposed above, and we do
not have an entropy component either because the setup is
isothermal. The planetary mass is sufficiently large, and the
viscosity sufficiently small, for the torque to be the (unsaturated)
horseshoe drag <cit.>.

Writing the unsaturated corotation torque as

    Γ _T^UHD =k'βΓ_0,

where k' is the dimensionless coefficient that we
want to determine, we can write:

    Γ^(U)_β/Γ_0=K-kβ+k'β,
 
where K=-2.34+0.099α <cit.>. We can also write the
torque value at larger time as

    Γ^∞_β/Γ_0=K-kβ+ϵ k'β,

where ϵ≪ 1 represents the amount of saturation of the corotation
torque. From these relations we infer:

    k-ϵ k'   =   Γ^∞_-1/2-Γ^∞_1/2/Γ_0≡ΔΓ^∞/Γ_0
    
      k-k'   =   Γ^(U)_-1/2-Γ^(U)_1/2/Γ_0≡ΔΓ^(U)/Γ_0,

hence

    k   =   ΔΓ^∞-ϵΔΓ^(U)/Γ_0(1-ϵ)
    
      k'   =   ΔΓ^∞-ΔΓ^(U)/Γ_0(1-ϵ).

We recover the fact that k=ΔΓ^∞/Γ_0 when
ϵ=0, that is to say when the corotation torque is fully
saturated at larger time. Fig. <ref> shows the torque
behaviour at early time, from which we infer
ΔΓ^(U)=0.3Γ_0. When ϵ=0, this yields 
k'=1.1. Evaluating k and k' when ϵ=0.1 (a reasonable
amount of residual corotation torque, as discussed in
section <ref>), we get k=1.5, and
k'=1.2. This shows that the coefficients that we have obtained for
the temperature dependence of the normalised Lindblad and corotation
torques are accurate to within ∼ 0.1.




We conclude this section by determining in a similar manner the linear
corotation torque. We simply have to amend the parameters so that the
corotation torque remains linear. This is achieved by adopting a much
smaller planet mass and a large viscosity
<cit.>. The results of these new runs are
displayed in Fig. <ref>, in which one can see that indeed
the torques remain nearly constant after a dynamical timescale and do
not show the variation toward a different value on O(10) orbits
typical of the horseshoe drag <cit.>. These
results imply that the temperature component of the linear corotation
torque scales as:

    Γ_T^lin=1.0βΓ_0,

where again the numerical coefficient is determined to within
∼ 0.1.  At this stage, considering that we use the saturation
function and the weights provided by <cit.>, we
have completed our determination of the terms entering in
Γ _T^CR (see Fig. <ref>, third frame).



 §.§ Torque dependence on the radial entropy gradient


We now perform a 3D version of the study presented by
<cit.> and update the terms entering in
Γ_S^CR (see Fig. <ref>,
second frame). We take an adiabatic index γ=1.4, unless
stated otherwise. We consider a large number of disc models (here N=80)
for which each value of the surface density and temperature slopes are
chosen randomly (namely, α is given by a random variable
uniformly distributed over [-1.5,1.5], while β is uniformly
distributed over [-2,2]). For each disc model, two runs are
performed: one with a locally isothermal equation of state, and one
with an adiabatic flow, which provide respectively the torque
Γ_iso and Γ_adi. The parameters and
the duration of the runs are such that the corotation torque is in the
regime of unsaturated horseshoe drag. We then seek an appropriate
linear combination ξ of α and β such that the torque
excess, defined as
ΔΓ=γΓ_adi-Γ_iso, be a
one-to-one map of ξ. We illustrate this process in
Fig. <ref>. We find that an appropriate expression for
ξ is the following:

    ξ=β-0.4α-0.64,

and that in this case we have:

    ΔΓ=4.6ξΓ_0.




By construction, our procedure implies that there is a linear
relationship between ξ and the torque excess, and therefore the
calculation of the adiabatic torque is straightforward once one knows
the isothermal torque. In a 2D situation, the excess has been found to
scale with the radial entropy gradient <cit.>. Here,
its expression resembles the gradient of entropy evaluated in the
midplane[Here this gradient is normalised so that the
  coefficient of β is one.], which is
β-2(γ-1)/(3-γ)α-3(γ-1)/(3-γ)=β-0.5α-0.75,
but differs slightly from it.  By extension, we call ξ the radial
entropy gradient, but it must be kept in mind that this is an abuse of
language, as this quantity has not been obtained by considerations
about the entropy, but by a fit of numerical simulations. As
mentioned in section <ref>, we assume that
this corotation torque component saturates as the entropy-related
torque in 2D discs.

We have not investigated in a systematic manner the
dependency of the different coefficients of Eq. (<ref>) on the
adiabatic index. Nonetheless, for the reader interested in the torque
in a different kind of discs, we mentioned that we have conducted a
similar study for the case γ=5/3, and found in that
case that ξ=β-0.58α-0.85 provides the one-to-one linear
relationship ΔΓ=4.1ξΓ_0.

At this stage we see that the horseshoe drag in the adiabatic case is
the sum of a term that scales with ξ (combination of α and
β) and of the isothermal horseshoe drag, itself sum of a term
that scales with the vortensity gradient (3/2-α in power
law discs) and a term that scales with the temperature gradient
(β) as found in section <ref>. We have
three components and only two slopes (α and β). It may
seem desirable at first glance to simplify the expression for the net
corotation torque into a dependence on α and one on
β. This, as discussed by <cit.>, would
obfuscate the physical meaning of each term. Besides, each component
involves a different distribution of vortensity perturbation within
the horseshoe region, and thus in principle saturates in a different
manner, so the simplification could only be valid for the regime of
unsaturated horseshoe drag, not for the general regime. We shall
therefore keep the three components in the general expression of the
corotation torque.

We have performed a numerical experiment similar to the one described
above, with a lower mass (q=10^-6) and a higher viscosity
(ν=2· 10^-5r_p^2Ω_p), so that the corotation torque
remains linear <cit.>. We find in that case:

    ΔΓ^lin=0.8ξΓ_0.





 §.§ Width of the horseshoe region


As said in section <ref>, the different
components of the horseshoe drag scale with the width of the horseshoe
region to the fourth power. It is therefore important to have an
accurate value of this width, so we study its scaling as a function of
the planet mass in the regimes of low and intermediate masses,
corresponding respectively to q≪ h^3 and q≲ h^3 (we shall
provide slightly more quantitative definitions at the end of this
section). For this purpose we use globally isothermal simulations, and
perform a streamline analysis in the midplane, exploiting the fact
that the horseshoe region has a width which is independent of the
altitude in this case
<cit.>. Our runs have a setup
similar to the setup outlined in section <ref>,
except for the resolution (N_ϕ,N_r,N_θ)=(1320,440,40). In
particular, we still have a null vortensity gradient. <cit.>
have found that the horseshoe region is asymmetric when the
vortensity gradient is finite. Even though we average the width of the
horseshoe region in front and at the rear of the planet in order to
minimise the effect of an asymmetry, we prefer not to introduce a
source of asymmetry of the horseshoe region[The result of
  <cit.> was obtained in two-dimensional discs. No similar
  result has been obtained in three-dimensional discs, so we assume
  that adopting a null vortensity gradient keeps the degree of
  asymmetry to an acceptable level.]. The planet mass is introduced
progressively over two orbital periods, and the width of the horseshoe
region is measured 20 orbits after the introduction of the
planet. We perform a systematic study for three values of the aspect
ratio: h=0.03, h=0.04 and h=0.05. For each of these values, we
perform 30 calculations with a planetary mass varying in a geometric
sequence from 0.05h^3 to 4h^3. The width of the horseshoe region
is determined automatically by dichotomy, averaging the results in
ϕ=1 rad and ϕ=-1 rad. We find that more accurate results are
obtained if the disc is inviscid when q<0.4h^3, and viscous for
larger masses. For lower mass indeed, a viscous drift of the disc can
distort severely the horseshoe region <cit.>, whereas for
high masses, the planet can trigger vortices on the edges of its
horseshoe region if the disc is inviscid, which precludes an accurate
determination of the width of the horseshoe region.

For low masses, we have the scaling
<cit.>:

    x_s=1.05r_p√(%s/%s)qh,

whereas for large masses we recover the scaling found in 2D simulation
with softened potential <cit.>:

    x_s≈ 2.5r_p(q/3)^1/3.

Using the width normalised to the pressure length scale
X_s=x_s/(hr_p) and denoting with Q the planetary mass normalised
to the thermal mass M_th=c_s^3/(GΩ_p):

    Q=M_p/M_th=q/h^3,

we can rewrite Eqs. (<ref>) and (<ref>) as:

    X_s=1.05Q^1/2Q≪ 1

and

    X_s=1.7Q^1/3Q≫ 1.

The lack of an explicit reference to h in these expressions suggests
that the normalised width of the horseshoe region has a universal
dependence on the mass expressed in thermal
masses. Fig. <ref> confirms this expectation. We find
that a function that gives the asymptotic behaviours of
Eq. (<ref>) and (<ref>) and that matches satisfactorily
the width behaviour for intermediate masses is a blend of the
asymptotic values at low- and large-mass with weights respectively
ϵ and 1-ϵ, where ϵ=1/(1+2Q^2). This yields:

    X_s=1.05Q^1/2+3.4Q^7/3/1+2Q^2,

or

    x_s=1.05(q/h)^1/2+3.4q^7/3/h^6/1+2q^2/h^6r_p.


From now on we will use Eq. (<ref>) to evaluate the width of the
horseshoe region.

As a side result, this leads us to a reasonable definition of an
intermediate mass protoplanet as a protoplanet that has a
horseshoe width intermediate between the widths given by the low-mass
and high-mass scaling (Eqs. (<ref>)
and (<ref>)). Examination of Fig. <ref> shows
that this broadly corresponds to 0.2h^3≲ q≲ 2h^3.


We assume that the corotation torque tends toward its linear value
when dissipation (either viscous or thermal) is large, regardless of
the planetary mass. This result has been shown by
<cit.> for low-mass planets, but it may be
questionable for intermediate mass objects. Although we have not
explored this issue in a systematic manner, we have performed
simulations with a planet of half the thermal mass (Q=0.5), and
found that it is subjected, at very large viscosity, to a vortensity
corotation torque comparable to the linear estimate. This result was
obtained considering globally isothermal discs with and without
vortensity gradients[The unsaturated total torque in
  isothermal simulations when Q∼ 0.5 can be off by several
  Γ_0 from its expected value, whereas we found with subsidiary
  calculations of isentropic discs that the total torque is correctly
  reproduced. These different behaviours have also been observed by
  <cit.> and are linked to
  small scale features of the flow. Since the realistic discs in which
  we apply our torque formula are not isothermal, we believe that this
  discrepancy is unimportant, but it should be kept in mind when
  dealing with intermediate mass planets embedded in isothermal
  discs.].  We note, however, that a substantial amount of dissipation
is required for a fall back of the corotation torque of an
intermediate mass object to its linear estimate. The decay toward the
linear regime is obtained for <cit.>:

    hr_pν/Ω_px_s^3≫ 0.1,

and a similar condition holds, obtained by writing χ in stead of
ν, for the decay of the entropy torque to its linear value. We
introduce the disc's alpha parameter α_ss
<cit.>, and write
ν=α_ssr_p^2Ω_ph^2. Using for x_s the
conservative value given by the low-mass estimate, we recast the
condition of Eq. (<ref>) as:

    α_ss≫ 0.1Q^3/2.

For Q=0.2, which corresponds to the mass at which departure from the
small mass regime becomes noticeable, this translates into
α_ss≫ 10^-2. This implies that under most
circumstances the corotation torque exerted on intermediate mass
objects should essentially be the horseshoe drag.




 §.§ Thermal diffusion coefficient


Previous estimates of the thermal diffusion coefficient have been
obtained by using values at the midplane of the disc
<cit.>. We
have undertaken simulations to check whether heat diffusion can indeed
be described by a diffusion coefficient evaluated at the disc's
midplane.

Following <cit.>, we write the thermal
diffusivity, in the absence of scattering, as (their Eq. B.2):


    χ=D/c_vρ=λ4acT^3/ρ^2c_vκ,

where D is the heat diffusion coefficient, λ, the
flux-limiter, is 1/3 in the optically thick regions where the planet
is located, a is the radiation constant and c_v is the specific
heat at constant volume. We can recast Eq. (<ref>) as:

    χ=16(γ-1)σ T^3/3ρ^2( R/μ)κ,

where σ is Stefan's constant. We note that this expression has
oftentimes been transformed into an expression involving the disc's
thickness H and orbital frequency Ω at the denominator. These
expressions, in general, overestimate by a factor γ the thermal
diffusivity of Eq. (<ref>). The equation that describes the
hydrostatic equilibrium of the disc does not involve the adiabatic
index, and it yields the relationship H=c_s/Ω, rather than
H=c_s^adi/Ω=√(γ)c_s/Ω, independently of
the fact that an energy equation is used to describe the flow.
However, we suggest to simply use Eq. (<ref>) rather than the
transformed expression, as the latter, which features all the
variables of Eq. (<ref>) as well as new ones, does not bring any
simplification.

Here, we specifically study the radial diffusion of heat, as it
is mainly this effect which determines (together with viscous
diffusion) the degree of saturation of the horseshoe drag.  Namely, we
study the radial spread of a radially localised temperature excess in
order to determine experimentally the thermal diffusivity of the disc
and compare it to the theoretical estimate.  In this section only we
use a module of radiative transfer in the flux limited diffusion
approximation <cit.> with the flux limiter of
<cit.>, and solve the equation of radiative
energy in addition to that of internal energy (two-temperature
approach) as in <cit.>.  We use a meridional mesh
of size (N_ϕ,N_r,N_θ)=(1,800,43) spanning a radial range of
0.75 au to 3 au, and a range in colatitude
[π/2-0.1,π/2]. The surface density at 1 au is
Σ_0=1700 gcm^-2, the Rosseland opacity is fixed to
1.8 cm^2g^-1, and the kinematic viscosity is
ν=4.46· 10^14 cm^2s^-1.  The slope of surface
density is α=0.3. The temperature profile of the disc at
equilibrium depends on a balance between viscous heating and radiative
losses through the disc's photospheres. We firstly run our initial
setup, which has an arbitrary aspect ratio, in order to relax the disc
toward hydrostatic and thermal equilibria. Once these are reached, we
increase by 1 % the temperature of a three-zone wide ring at
r=1.5 au (over the whole colatitude range), and restart the
simulation (hot ring run). We also perform a “neutral” restart in
which the temperature has not been altered. The radial profile of the
temperature excess obtained by subtracting the temperature profile of
the neutral run from the hot ring run (at the same date), quickly
adopts a nearly Gaussian profile:

    δ T(r,t)∝1/σ(t)exp[-(r-r_0)^2/σ^2(t)].

In linear diffusion theory, σ^2(t)=4χ t. We use this
relationship to infer the value of χ. Fig. <ref> shows
the time behaviour of σ^2, which is nearly linear in time, as
expected. From the slope of the linear regression fit, we can estimate
that the value of χ  at r=r_0 is:

    χ≈ 1.03· 10^15 cm^2 s^-1.

An estimate based on the temperature excess averaged in colatitude,
rather than measured at the midplane, yields nearly exactly the same
result. At r=1.5 au, we have in the neutral run a midplane
temperature T=509 K, a density ρ=4.33·
10^-10 gcm^-3, and R/μ=3.615·
10^7 erg K^-1 g^-1. Eqs. (<ref>) yields a value 27 %
larger than the value found experimentally.



The disc of this numerical experiment is in radiative equilibrium and
its only source of heating is viscous friction. Its temperature
therefore decays from the midplane to low values near the boundaries
in colatitude. We have repeated this numerical experiment with a
nearly isothermal vertical profile, obtained by lowering the viscosity
by a factor of ten and by imposing a temperature at the upper boundary
T_b such that T_b^4=0.9T^4_midplane. The vertical heat
flux is therefore lowered by an order of magnitude with respect to the
previous disc, the temperature is nearly constant with respect to
colatitude, and the resulting disc has same midplane temperature as
the previous one. 

We measured a thermal diffusivity nearly equal to the one of the
previous disc, which shows that this quantity is rather insensitive to
the vertical profile of temperature, and determined by the midplane
temperature.





§ NEW TORQUE FORMULA


We hereafter summarise our results and provide an updated torque
formula for different cases: the general case, the linear regime for
a radiative disc, and the linear regime for an isothermal disc.



 §.§ General case


Usually torque expressions are given in terms of the reference torque
Γ_0 of Eq. (<ref>). For the horseshoe drag, this
expression arises naturally when using the low mass value of the
horseshoe width (Eq. <ref>), and an horseshoe drag expression,
which involves the product Σ_0Ω_p^2x_s^4. Here, however,
we use a more complex law (Eq. <ref>) for the horseshoe width,
so that casting our results in terms of Γ_0, for the different
components of the corotation torque, would lead to cumbersome
expressions. We therefore rather cast the different components of the
horseshoe drag in terms of x_s.

The total torque is a function of the slopes of surface density and
temperature α and β, of the disc's opacity, viscosity,
aspect ratio and surface density (respectively κ, ν, h
and Σ_0) and of the planet orbital frequency, orbital radius
and mass to star mass ratio (respectively Ω_p, r_p and
q). It requires the prior evaluation of the thermal diffusivity
using Eq. (<ref>), and that of the half width of the horseshoe
region x_s, which is given, when the flow is adiabatic or behaves
nearly adiabatically over the horseshoe U-turn timescale, by a
modification of Eq. (<ref>) that uses the adiabatic sound speed,
rather than the isothermal one:

    x_s=1.05(q/h')^1/2+3.4q^7/3/h'^6/1+2q^2/h'^6r_p,

where h'=h√(γ).

The total torque is the sum of the Lindblad torque Γ_L and
corotation torque Γ_C:

    Γ_tot=Γ_L+Γ_C.

The Lindblad torque is given by:

    Γ_L=-(2.34-0.1α+1.5β)Γ_0f(χ/χ_c),

where the first two coefficients in the factor of the right hand side
(2.34 and -0.1) are given by <cit.>, and where an
approximate value of the third coefficient (the factor of β) has
been obtained in section <ref> and subsequently
slightly improved in section <ref>. The function
f(x) and critical diffusivity χ_c are given by
<cit.> and have respectively the expressions:

    f(x)=(x/2)^1/2+1/γ/(x/2)^1/2+1,

and:

    χ_c=r_p^2h^2Ω_p.

The function f(χ/χ_c) can be regarded as the inverse of an
effective adiabatic index γ_eff <cit.>.

We now turn to the evaluation of the four terms that compose the
corotation torque (see Fig. <ref>). The vortensity component
Γ_V^CR of the corotation torque is given by:


    Γ_V^CR=ε_bΓ_V^HD+(1-ε_b)Γ_V^lin.

In this expression, the blending coefficient ε_b, as
advertised in section <ref>, has the expression
worked out by <cit.>:

    ε_b=(1+30hz_ν)^-1,

where z_ν is given by <cit.>:

    z_ν=r_pν/Ω_px_s^3.

The (vortensity components of the) horseshoe drag
Γ^HD_V and linear corotation torque
Γ_V^lin are respectively:

    Γ_V^HD= F_VΓ_V^UHD,

and

    Γ_V^lin=(0.976-0.640α)Σ_0Ω_p^2r_p^4(h')^-2=(0.976-0.640α)Γ_0/γ,

where the numerical coefficients of the right hand side of
Eq. (<ref>) can be obtained from the data given by
<cit.>.  In Eq. (<ref>) the two factors of the right
hand side are respectively the saturation function of the vortensity
component of the horseshoe drag, and the unsaturated horseshoe
drag. The former has the expression <cit.>:

    F_V=8π/3z_ν F(z_ν),

where F(x) is defined as:

    F(x)= {[ 1-x^1/2        ; 4/(27x)         ].

and the latter has the expression:

    Γ^UHD_V=3/4(3/2-α)Σ_0Ω_p^2x_s^4.

This last expression, initially established for two dimensional discs
<cit.>, has been generalised to three
dimensional discs <cit.>. Eqs. (<ref>),
(<ref>), (<ref>), (<ref>) and (<ref>)
correspond to the five cells of the first frame of
Fig. <ref>, from left to right and bottom to top.

Next, we give the expression of the entropy component of the
corotation torque, corresponding to the second frame of
Fig. <ref>. It reads[Note that in
  <cit.>, the thermal diffusivity was denoted by
  κ, whereas here it is denoted by χ, while κ is the
  opacity. As a consequence, ε_κ (z_κ) has been
  changed into ε_χ (z_χ).]:


    Γ_S^CR=ε_νε_χΓ_S^HD+(1-ε_νε_χ)Γ_S^lin.

As for the vortensity torque, we keep for the blending coefficients
the expressions given by <cit.>:

    ε_ν=[1+(6hz_ν)^2]^-1

and

    ε_χ=(1+15hz_χ)^-1,

where z_χ is defined in a similar fashion as z_ν <cit.>:

    z_χ=r_pχ/Ω_px_s^3.

The linear component of the torque, as seen in
section <ref>, is:

    Γ_S^lin=0.8ξΓ_0/γ,

where ξ is given by:

    ξ=β-0.4α-0.64.

The horseshoe drag Γ_S^HD of Eq. (<ref>) is
simply the product of the unsaturated horseshoe drag by the saturation
function of the entropy corotation torque:

    Γ_S^HD= F_SΓ_S^UHD.

In two dimensions, the expression of the unsaturated horseshoe drag
can be written KξΣ_0Ω_p^2x_s^4, where K is a numerical
constant <cit.>. We assume this dependency to hold in three
dimensions, and seek the value of K using the results of
section <ref>. In the numerical exploration that
we performed, we had, using Eq. (<ref>),
x_s=1.086r_p(q/h')^1/2, from which we infer K=3.3. The
unsaturated torque expression is therefore:

    Γ_S^UHD=3.3ξΣ_0Ω_p^2x_s^4.

For the saturation function we keep the dependence given by
<cit.>:

    F_S=1.2×1.4z_χ^1/2×1.8z_ν^1/2,

where x=min(1,x). Eqs. (<ref>), (<ref>),
(<ref>), (<ref>) and (<ref>) correspond to the five
cells of the second frame of Fig. <ref>, from left to right
and bottom to top.

The third frame of Fig. <ref> corresponds to the temperature
component of the corotation torque, that we now evaluate:

    Γ_T^CR=ε_νΓ_T^HD+(1-
    ε_ν)Γ_T^lin.

Since the vortensity distribution associated to this torque component
is roughly similar to that of the entropy torque (it is concentrated
near the separatrices), we assume that it behaves like the entropy
torque as a function of viscosity. However, we discard the dependence
on thermal diffusivity. Indeed, in the limit of a large diffusivity,
the disc behaves isothermally and the temperature torque subsists (as
we saw in section <ref>), on the contrary of the
entropy torque. The linear component that features in
Eq. (<ref>) reads (see section <ref>):

    Γ_T^lin=1.0βΓ_0/γ,

while the horseshoe drag reads:

    Γ_T^HD= F_TΓ_T^UHD.

From the above discussion we adopt:

    F_T=1.2×1.8z_ν^1/2.

For the unsaturated horseshoe drag, we assume the form
K'βΣ_0Ω_p^2x_s^4, i.e. similar to that of the entropy,
but with a different numerical constant and a scaling on the
temperature gradient. The toy model presented by <cit.> suggests
it is reasonable, at least in the low mass limit, as it displays a
scaling in x_s^2Δ P, where Δ P is the perturbation of
pressure at the stagnation point. This quantity, in the low mass
limit, scales as Σ_0Ω_p^2x_s^2 if the stagnation point is
typically at a pressure scale length away from the planet. We seek the
value of K' using the results of
section <ref>. In that section, we had, using
Eq. (<ref>) (rather than Eq. <ref> since the disc was
isothermal), x_s=1.13r_p(q/h)^1/2. Equating the form given above
with that of Eq. (<ref>), we infer K'=k'/1.13^4=0.73 and thus:

    Γ_T^UHD=0.73βΣ_0Ω_p^2x_s^4.

Eqs. (<ref>), (<ref>), (<ref>), (<ref>) and
(<ref>) correspond to the five cells of the third frame of
Fig. <ref>, from left to right and bottom to top.


In a fourth and last stage, we evaluate the additional component,
arising from the viscous creation of vortensity at the abrupt density
jumps that appear at the contact discontinuities at the separatrices
of the horseshoe region. For brevity and definiteness we hereafter
call this term the viscous coupling term (VCT).
<cit.> give an expression of this term (the
second part of the bracket of their Eq. 129). This term corresponds to
an asymptotic value (at larger time) and therefore already embeds its
own saturation function. It naturally scales with the entropy
gradient, as it is this gradient which primarily determines the
magnitude of the density jumps at the separatrices. As emphasised in
section <ref>, the quantity ξ was determined
from a best fit of numerical data, and named entropy gradient on the
grounds of the very close resemblance of the effect studied in that
section with the well-studied two-dimensional process that triggers
the appearance of an additional torque component that scales with the
radial entropy gradient <cit.>. We therefore suggest that
the three-dimensional version of the VCT scales with ξ.  Our
choice of normalisation of ξ is that the coefficient of β is
unity, whereas it is -1/γ with the normalisation of
<cit.>. As this torque behaves as a bulk term (it
implies a smooth distribution of vortensity from corotation to the
separatrices), and since it has no known linear equivalent, we follow
<cit.> and assume it to decay like the vortensity
component of the horseshoe drag, proportionally to ε_b:

    Γ_VCT^CR=4πξ/γΣ_0Ω_p^2x_s^4ϵ_bz_ν[z_ν
          F(z_ν)-z_χ F(z_χ)/z_ν-z_χ]

The full corotation torque is finally given by:

    Γ_C=Γ_V^CR+Γ_S^CR+Γ_T^CR+Γ_VCT^CR.

Eqs. (<ref>) to (<ref>), together with Eq. (<ref>),
provide an expression for the total tidal torque exerted on a low to
intermediate mass planet in circular orbit in an optically thick disc in hydrostatic
and radiative equilibrium.



 §.§ Linear regime in the general case


We give hereafter the simpler expression of the torque in the linear
regime. It reduces to:

    Γ_tot^lin=Γ_L+Γ_V^lin+Γ_S^lin+Γ_T^lin,

which simplifies as:

    Γ_tot^lin=Γ_L+(0.46-0.96α+1.8β)Γ_0/γ,

where Γ_L is given by Eq. (<ref>).



 §.§ Linear regime in an isothermal disc


We can further simplify this expression in the case of an isothermal
disc (in which case we discard the entropy component of the
corotation torque). This yields:

    Γ_tot^lin,iso=-(1.36+0.54α+0.5β)Γ_0.

This result is in good agreement with the result obtained by
<cit.> in three-dimensional
simulations of locally isothermal discs.











§ DISCUSSION


We compare here the torque value provided by our formula of
section <ref> to published studies of the
torque value in radiative discs, either as a function of the planetary
mass or of the orbital radius. Fig. <ref> shows the
comparison to results recently published by
<cit.>. They show a broad agreement between the
torque measured in their simulations and our torque formula. They also
show that the boost of the horseshoe width, described by
Eq. (<ref>), is an essential ingredient of the torque formula
above 10-15 M_⊕ (i.e. for intermediate mass planets),
and that torque formulae based on a low-mass expression of the width
of the horseshoe region hardly yield a migration reversal.

It should be kept in mind that the torque measured in numerical
simulations is subjected to some inaccuracy, as can be seen from the
results for 20 M_⊕ (in the first two plots), for which two
torques values are displayed. This inaccuracy arises from a variety of
numerical effects, such as resolution (in particular the number of
zones over which the radial width of the horseshoe region is
resolved), or the recipe for the torque calculation (in particular
whether it includes all zones or whether it excludes those located
near the planet), etc.  Yet, some points or parts of the torque curves
are outliers with respect to our formula:


  * At low planetary mass (on the left plot), the torque measured is
  somehow below the value predicted by our formula. This behaviour was
  already discussed by <cit.>, and attributed to
  the “cold-finger” effect <cit.>, which is
  observed on low-mass planets when there is thermal diffusion in the
  disc. We note that <cit.> have obtained their
  Eq. (156) — here Eq. (<ref>) — using a fit of numerical
  simulations with different values of the thermal diffusivity, and
  may have unwittingly included a two-dimensional version of the
  “cold-finger” effect. It is therefore unclear to which extent the
  discrepancy found at low planetary mass is due to this effect. We
  also note that for the lowest planetary mass considered here
  (5 M_⊕), <cit.> found that the cold
  finger effect nearly vanishes. This issue requires further work, and
  will be presented elsewhere.

  * The plateau value for the mass range 15-35 M_⊕ on the
  middle plot, or the elbow of the left plot (at 38 M_⊕
  according to our formula, and at ∼ 50 M_⊕ in the results of
  <cit.>) are also discrepant features, which are
  unlikely to be explained by the errors and systematic effects on the
  torque measurement. 

Beside the “cold-finger” effect, which might account for the mismatch
at low mass, there is a number of simplifications in our torque
formula which can explain the residual discrepancies between our
torque formula and the results of numerical simulations.  We neglect
the feed back of the planet's torque on the disc density profile. At
intermediate masses, this feed back tends to create a dip around the
orbit (precursor of the gap that would be carved at larger masses),
which has an effect on the different components of the torque. In
particular the effect can be different on the outer Lindblad torque
and on the inner Lindblad torque, resulting in a non-trivial effect on
the net Lindblad torque. Also, the slight decrease of density can
alter the thermal diffusivity at the planetary orbit, by virtue of
Eq. (<ref>). On the other hand, as we saw in
section <ref>, this equation is accurate only to
within ∼ 25 %, and both the saturation functions and the
blending coefficients between the linear torque and the horseshoe drag
depend sensitively on the value of the thermal diffusivity.  Also, we
use for these functions the analytic dependence worked out in two
dimensions by <cit.>, which can constitute
another reason for the discrepancy. Lastly, we have used for all
components of the corotation torque the actual adiabatic index of the
gas γ, instead of an effective adiabatic index as done for the
Lindblad torque <cit.>. We comment that we
have performed runs with planets of mass 1 M_⊕ to
20 M_⊕ embedded in radiative discs with different opacities
(namely 1.8 cm^2.g^-1 and 5 cm^2.g^-1). We have found
that the width of the horseshoe region at low mass corresponds indeed
to x_s=1.05√(q/h)/γ^1/4 rather than
x_s=1.05√(q/h)/γ^1/4_eff. At larger masses, it
displays a growth above this value, as expected for the transition
from the low-mass to the high-mass branch of the horseshoe
width. Nonetheless, owing to thermal diffusion, the actual law might
be more complex than the one of Eq. (<ref>) and this can have an
impact on the torque value, which is a very sensitive function of
x_s.

We finally compare the torque predicted by our formula to that
measured by <cit.> in Fig. <ref>. For
radii below 2.5 a_Jup, there is a satisfactory agreement
between the formula and the numerical simulations. Beyond that radius,
however, the predicted torque plummets to low values, whereas the
torque measured in the simulations plateaus at an intermediate,
negative value. We note that this happens when the disc's optical
depth drops below value of a few (at r=a_Jup it is
τ_eff∼ 4). Note that the actual optical depth of the
disc is somehow lower than the one displayed, as the opacity law
adopted by the authors was:
κ/(1 cm^2 g^-1)=2.0·
10^-4[T/(1 K)]^2 (Bitsch, private communication). The
heat source of the disc being viscous friction, it is hotter at the
midplane, and therefore its opacity is largest at the midplane. This
example illustrates the fact that Eq. (<ref>), and therefore our
torque formula, are valid when the disc is optically
thick. Fig. <ref> also shows that torque formulae based on the
low-mass scaling of the horseshoe width hardly predict any migration
reversal.

We comment that the dependence of the corotation torque on viscosity
and thermal diffusivity is entirely contained in the saturation
functions and blending coefficients. Since those are not reconsidered
in the present work, our torque has same dependence on viscosity and
thermal diffusivity as the dependence found by
<cit.>. Namely, Figs. 13, 17 and 18 of that work
show the dependence of the vortensity torque in the viscosity, and the
dependence of the entropy torque on the viscosity and thermal
diffusivity. The temperature term of the corotation torque was not
considered in that work, but it exhibits same dependence on viscosity
as the entropy torque.



§ CONCLUSION


We provide an updated formula for the torque experienced by a low mass
(q < 0.2h^3) or intermediate mass planet
(0.2h^3≲ q≲ 2h^3) in circular and non-inclined orbit
in an optically thick disc. Our torque formula agrees reasonably well
with the torque measured in numerical simulations, as shown by
comparison with recently published results. One key ingredient of the
formula is an accurate expression of the width of the horseshoe region
as a function of the planet-to-star mass ratio and of the disc's
aspect ratio.  Using in our torque formula the low-mass scaling
x_s=1.05r_p√(q/h) clearly yields wrong, largely underestimated
results for intermediate mass planets. We have used tailored,
three-dimensional explorations of the parameter space to update
several numerical coefficients that appear in the torque formula. Our
work does not consider the effect found at low planetary mass by
<cit.> in discs with thermal diffusion, and
called by these authors the “cold-finger” effect. A further step
toward accurate torque formulae could be an analytic expression of
this torque component, as well as a relaxation of the different
simplifying assumptions mentioned in section <ref>.




§ ACKNOWLEDGEMENTS

The authors thank Elena Lega and Bertram Bitsch for an extremely  
careful reading of a first draft of this manuscript, and for  
constructive feedback.  
Alejandra Jiménez acknowledges financial support from a UNAM DGAPA
fellowship. The simulations presented in this work were run on a GPU
cluster acquired with CONACyT's grant 178377. The authors acknowledge
support from UNAM's grant PAPIIT IN101616.



mnras

@urlcharsothermakeother $&#_%
@doi@urlcharsother ifnextchar [ @doi@
  @doi@[]
@doi@[#1]#2tempa#1tempaempty http://dx.doi.org/#2 doi:#2http://dx.doi.org/#2 #1
@eprint#1#2@eprint@#1:#2::nil
@eprint@arXiv#1http://arxiv.org/abs/#1 arXiv:#1
@eprint@dblp#1http://dblp.uni-trier.de/rec/bibtex/#1.xml
  dblp:#1
@eprint@#1:#2:#3:#4niltempa #1tempb #2tempc
  #3tempc empty tempc tempb tempb tempa tempb empty tempb arXivifundefined
  mn@eprint@tempbtempb:tempcmn@eprint@tempbtempc

[Baruteau & MassetBaruteau &
  Masset2008]bm08
Baruteau C.,  Masset F.,  2008, @doi [] 10.1086/523667, http://adsabs.harvard.edu/abs/2008ApJ...672.1054B 672, 1054

[Baruteau & MassetBaruteau &
  Masset2013]2013LNP...861..201B
Baruteau C.,  Masset F.,  2013, in Souchay J.,  Mathis S.,   Tokieda
  T.,  eds,  Lecture Notes in Physics, Berlin Springer Verlag Vol. 861, Lecture
  Notes in Physics, Berlin Springer Verlag. p. 201 (@eprint arXiv
  1203.3294), @doi10.1007/978-3-642-32961-6_6

[Benítez-Llambay &
  MassetBenítez-Llambay & Masset2016]2016ApJS..223...11B
Benítez-Llambay P.,  Masset F. S.,  2016, @doi []
  10.3847/0067-0049/223/1/11, http://adsabs.harvard.edu/abs/2016ApJS..223...11B 223, 11

[Bitsch & KleyBitsch &
  Kley2010]2010A   A...523A..30B
Bitsch B.,  Kley W.,  2010, @doi [] 10.1051/0004-6361/201014414,
  http://adsabs.harvard.edu/abs/2010A.26A...523A..30B 523, A30

[Bitsch & KleyBitsch &
  Kley2011]2011A   A...536A..77B
Bitsch B.,  Kley W.,  2011, @doi [] 10.1051/0004-6361/201117202,
  http://adsabs.harvard.edu/abs/2011A.26A...536A..77B 536, A77

[Bitsch, Crida, Morbidelli, Kley  &
  Dobbs-DixonBitsch et al.2013]2013A   A...549A.124B
Bitsch B.,  Crida A.,  Morbidelli A.,  Kley W.,   Dobbs-Dixon I.,
  2013, @doi [] 10.1051/0004-6361/201220159, http://adsabs.harvard.edu/abs/2013A.26A...549A.124B 549, A124

[Bitsch, Morbidelli, Lega  &
  CridaBitsch et al.2014]2014A   A...564A.135B
Bitsch B.,  Morbidelli A.,  Lega E.,   Crida A.,  2014, @doi []
  10.1051/0004-6361/201323007, http://adsabs.harvard.edu/abs/2014A.26A...564A.135B 564, A135

[Casoli & MassetCasoli &
  Masset2009]cm09
Casoli J.,  Masset F. S.,  2009, @doi []
  10.1088/0004-637X/703/1/845, http://adsabs.harvard.edu/abs/2009ApJ...703..845C 703, 845

[Cossou, Raymond, Hersant  &
  PierensCossou et al.2014]2014A   A...569A..56C
Cossou C.,  Raymond S. N.,  Hersant F.,   Pierens A.,  2014, @doi
  [] 10.1051/0004-6361/201424157, http://adsabs.harvard.edu/abs/2014A.26A...569A..56C 569, A56

[Crida, Morbidelli  & MassetCrida
  et al.2006]crida06
Crida A.,  Morbidelli A.,   Masset F.,  2006, @doi [Icarus]
  10.1016/j.icarus.2005.10.007, http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=2006Icar..181..587C   db_key=AST
  181, 587

[D'Angelo & LubowD'Angelo &
  Lubow2010]2010ApJ...724..730D
D'Angelo G.,  Lubow S. H.,  2010, @doi []
  10.1088/0004-637X/724/1/730, http://adsabs.harvard.edu/abs/2010ApJ...724..730D 724, 730

[Dong, Rafikov  & StoneDong
  et al.2011]2011ApJ...741...57D
Dong R.,  Rafikov R. R.,   Stone J. M.,  2011, @doi []
  10.1088/0004-637X/741/1/57, http://adsabs.harvard.edu/abs/2011ApJ...741...57D 741, 57

[DuffellDuffell2015]2015ApJ...806..182D
Duffell P. C.,  2015, @doi [] 10.1088/0004-637X/806/2/182, http://adsabs.harvard.edu/abs/2015ApJ...806..182D 806, 182

[Fromang, Lyra  & MassetFromang
  et al.2011]2011A   A...534A.107F
Fromang S.,  Lyra W.,   Masset F.,  2011, @doi []
  10.1051/0004-6361/201016068, http://adsabs.harvard.edu/abs/2011A

[Fung, Shi  & ChiangFung
  et al.2014]2014ApJ...782...88F
Fung J.,  Shi J.-M.,   Chiang E.,  2014, @doi []
  10.1088/0004-637X/782/2/88, http://adsabs.harvard.edu/abs/2014ApJ...782...88F 782, 88

[Fung, Artymowicz  & WuFung
  et al.2015]2015arXiv150503152F
Fung J.,  Artymowicz P.,   Wu Y.,  2015, @doi []
  10.1088/0004-637X/811/2/101, http://adsabs.harvard.edu/abs/2015ApJ...811..101F 811, 101

[Fung, Masset, Lega  & VelascoFung
  et al.2017]2017AJ....153..124F
Fung J.,  Masset F.,  Lega E.,   Velasco D.,  2017, @doi []
  10.3847/1538-3881/153/3/124, http://adsabs.harvard.edu/abs/2017AJ....153..124F 153, 124

[Kanagawa, Muto, Tanaka, Tanigawa,
  Takeuchi, Tsukagoshi  & MomoseKanagawa
  et al.2015]2015ApJ...806L..15K
Kanagawa K. D.,  Muto T.,  Tanaka H.,  Tanigawa T.,  Takeuchi T.,
  Tsukagoshi T.,   Momose M.,  2015, @doi []
  10.1088/2041-8205/806/1/L15, http://adsabs.harvard.edu/abs/2015ApJ...806L..15K 806, L15

[KlahrKlahr2004]2004ApJ...606.1070K
Klahr H.,  2004, @doi [] 10.1086/383119, http://adsabs.harvard.edu/abs/2004ApJ...606.1070K 606, 1070

[KleyKley1989]1989A   A...208...98K
Kley W.,  1989, , http://adsabs.harvard.edu/abs/1989A

[Kley, Bitsch  & KlahrKley
  et al.2009]2009A   A...506..971K
Kley W.,  Bitsch B.,   Klahr H.,  2009, @doi []
  10.1051/0004-6361/200912072, http://adsabs.harvard.edu/abs/2009A

[Korycansky & PollackKorycansky &
  Pollack1993]1993Icar..102..150K
Korycansky D. G.,  Pollack J. B.,  1993, @doi [Icarus]
  10.1006/icar.1993.1039, http://adsabs.harvard.edu/abs/1993Icar..102..150K 102, 150

[Lega, Crida, Bitsch  &
  MorbidelliLega et al.2014]2014MNRAS.440..683L
Lega E.,  Crida A.,  Bitsch B.,   Morbidelli A.,  2014, @doi
  [] 10.1093/mnras/stu304, http://adsabs.harvard.edu/abs/2014MNRAS.440..683L 440, 683

[Lega, Morbidelli, Bitsch, Crida  &
  SzulágyiLega et al.2015]2015MNRAS.452.1717L
Lega E.,  Morbidelli A.,  Bitsch B.,  Crida A.,   Szulágyi J.,
  2015, @doi [] 10.1093/mnras/stv1385, http://adsabs.harvard.edu/abs/2015MNRAS.452.1717L 452, 1717

[Levermore & PomraningLevermore &
  Pomraning1981]1981ApJ...248..321L
Levermore C. D.,  Pomraning G. C.,  1981, @doi [] 10.1086/159157,
  http://adsabs.harvard.edu/abs/1981ApJ...248..321L 248, 321

[MassetMasset2000]fargo2000
Masset F.,  2000, , http://adsabs.harvard.edu/abs/2000A

[MassetMasset2001]masset01
Masset F. S.,  2001, , http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=2001ApJ...558..453M   db_key=AST
  558, 453

[MassetMasset2002]masset02
Masset F. S.,  2002, , http://cdsads.u-strasbg.fr/cgi-bin/nph-bib_query?bibcode=2002A
  387, 605

[Masset & Benítez-LlambayMasset
  & Benítez-Llambay2016]2016ApJ...817...19M
Masset F. S.,  Benítez-Llambay P.,  2016, @doi []
  10.3847/0004-637X/817/1/19, http://adsabs.harvard.edu/abs/2016ApJ...817...19M 817, 19

[Masset & CasoliMasset &
  Casoli2009]mc09
Masset F. S.,  Casoli J.,  2009, @doi []
  10.1088/0004-637X/703/1/857, http://adsabs.harvard.edu/abs/2009ApJ...703..857M 703, 857

[Masset & CasoliMasset &
  Casoli2010]2010ApJ...723.1393M
Masset F. S.,  Casoli J.,  2010, @doi []
  10.1088/0004-637X/723/2/1393, http://adsabs.harvard.edu/abs/2010ApJ...723.1393M 723, 1393

[Masset & OgilvieMasset &
  Ogilvie2004]mo2004
Masset F. S.,  Ogilvie G. I.,  2004, , http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=2004ApJ...615.1000M   db_key=AST
  615, 1000

[Masset, D'Angelo  & KleyMasset
  et al.2006]mak2006
Masset F. S.,  D'Angelo G.,   Kley W.,  2006, @doi []
  10.1086/507515, http://adsabs.harvard.edu/abs/2006ApJ...652..730M
  652, 730

[Paardekooper & PapaloizouPaardekooper
  & Papaloizou2008]pp08
Paardekooper S.-J.,  Papaloizou J. C. B.,  2008, @doi []
  10.1051/0004-6361:20078702, http://adsabs.harvard.edu/abs/2008A.26A...485..877P 485, 877

[Paardekooper & PapaloizouPaardekooper
  & Papaloizou2009a]2009arXiv0901.2265P
Paardekooper S.-J.,  Papaloizou J. C. B.,  2009a, @doi []
  10.1111/j.1365-2966.2009.14511.x, http://adsabs.harvard.edu/abs/2009MNRAS.394.2283P 394, 2283

[Paardekooper & PapaloizouPaardekooper
  & Papaloizou2009b]2009arXiv0901.2263P
Paardekooper S.-J.,  Papaloizou J. C. B.,  2009b, @doi []
  10.1111/j.1365-2966.2009.14512.x, http://adsabs.harvard.edu/abs/2009MNRAS.394.2297P 394, 2297

[Paardekooper, Baruteau, Crida  &
  KleyPaardekooper et al.2010]pbck10
Paardekooper S.,  Baruteau C.,  Crida A.,   Kley W.,  2010, @doi
  [] 10.1111/j.1365-2966.2009.15782.x, http://adsabs.harvard.edu/abs/2010MNRAS.401.1950P 401, 1950

[Paardekooper, Baruteau  &
  KleyPaardekooper et al.2011]pbk11
Paardekooper S.,  Baruteau C.,   Kley W.,  2011, @doi []
  10.1111/j.1365-2966.2010.17442.x, http://adsabs.harvard.edu/abs/2011MNRAS.410..293P 410, 293

[Shakura & SunyaevShakura &
  Sunyaev1973]ss73
Shakura N. I.,  Sunyaev R. A.,  1973, , http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1973A
  24, 337

[Tanaka, Takeuchi  & WardTanaka
  et al.2002]tanaka2002
Tanaka H.,  Takeuchi T.,   Ward W. R.,  2002, , 565, 1257

[WardWard1986]ww86
Ward W. R.,  1986, @doi [Icarus] 10.1016/0019-1035(86)90182-X, http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1986Icar...67..164W   db_key=AST
  67, 164

[WardWard1991]wlpi91
Ward W. R.,  1991, in Lunar and Planetary Institute Conference Abstracts.
  p. 1463

[de Val-Borro et al.,de Val-Borro
  et al.2006]valborro06
de Val-Borro M.,  et al., 2006, @doi []
  10.1111/j.1365-2966.2006.10488.x, http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=2006MNRAS.tmp..695D   db_key=AST
  pp 695–+









