© ESO, 2011

1. Introduction

Planetary migration of embedded low-mass planets in isothermal discs indicates inward migration, so that the planet might be lost in the star before the accretion disc is gone (Tanaka et al. 2002). Recent studies (starting with the work of Paardekooper & Mellema 2006) have shown that the inclusion of radiation transport in planet-disc interaction studies resulted in a slowed down or even reversed migration (Baruteau & Masset 2008a; Paardekooper & Papaloizou 2008; Paardekooper & Mellema 2008; Kley & Crida 2008; Kley et al. 2009; Ayliffe & Bate 2010, 2011).

All authors agree that the inclusion of radiation transport is an important effect that should be considered, however, not all authors find outward migration. The efficiency of outward migration depends on the magnitude of positive corotation torques. These are determined by the local entropy and vortensity gradients. For isothermal and adiabatic discs these gradients cannot be maintained and so-called torque saturation reduces the corotation effects. However, the inclusion of radiative transport and viscosity prevents saturation, and the negative Lindblad torques (caused by the spiral arms) can be overwhelmed by the positive contributions from the corotation region. The magnitude of the effect depends on the planetary mass (Kley & Crida 2008; Kley et al. 2009) and the temperature gradient in the disc. For example, Ayliffe & Bate (2011) found for temperature slopes of β > 1.0 (with T ∝ r^−β) that outward migration is possible even for p lanets with up to 50 M_Earth in 3D simulations. The migration rate increases with an increasing temperature slope. These results are also reflected by the theoretical analysis from Masset & Casoli (2010). Recently, improved theoretical torque formulae for low-mass planets embedded in an adiabatic disc have been presented by Masset & Casoli (2009) and Paardekooper et al. (2010). The most recent improvements consider the necessary inclusions of radiative diffusion and viscosity (Masset & Casoli 2010; Paardekooper et al. 2011). These have been developed for 2D discs where the diffusive effects can operate only in the disc’s plane. Interestingly, radiative cooling alone can lead to similar, even stronger effects (Kley & Crida 2008).

A very important open question is how far the planet will move outwards in a fully radiative disc. When the protoplanets are stopped from migrating outwards at a certain point in the disc, a protoplanet moving inwards from farther out will stop at the same location in the disc. Therefore, a planetary trap inside the disc can be created, which can act as an area of planetary mergers, leading to larger cores. A planetary trap inside the disc created by surface density changes can function as a feeding zone to planetary cores (Morbidelli et al. 2008), but it is unclear how realistic these surface density changes are. The inclusion of radiation transport/cooling in a disc might provide such a trap in a normal disc structure, at least for planets within a given mass range.

In Sándor et al. (2011) the authors show in N-body simulations that because of the inclusion of unsaturated type-I-migration (Paardekooper et al. 2010) large planetary cores (up to 10 M_Earth) can form in protoplanetary discs well before the disc is accreted onto the star.

In population synthesis models that include the standard migration prescription for isothermal discs the migration is too fast, such that the outcome distribution does not match the observed one. Specifically, the type-I-migration rate should be 10 to 1000 times slower than expected from linear analysis (Alibert et al. 2004; Ida & Lin 2008; Mordasini et al. 2009). As pointed out above, this could be provided by the inclusion of radiation transport in hydrodynamical type-I-migration simulations. An updated version of the type-I-migration, by using the formula from Paardekooper et al. (2010), in population synthesis models shows very promising results (Mordasini et al. 2011), but additional simulations are needed to verify these assumptions.

In our previous simulations and studies (Kley et al. 2009; Bitsch & Kley 2010; Bitsch & Kley 2011) we have only considered one standard disc model with a given mass. In reality, of course, protoplanetary discs can have a variety of masses. A fully radiative disc will settle in an equilibrium state between viscous heating and radiative transport/cooling. Given that all other physical parameters are fixed, the resulting disc structure only depends on the disc’s mass and viscosity, because a more massive or more viscous disc is able to heat the disc more effectively. A more massive disc therefore influences the migration rate of embedded planets. Here we focus on low-mass planets that can undergo outward migration in fully radiative discs.

In this paper we extend previous studies and investigate the possible radial range over which outward migration can occur and analyse the influence of the disc’s mass on the migration in detail. An important effect is the onset of convection in the disc, which becomes stronger for more massive discs.

In Sect. 2 we give a short overview of our code and numerical methods, where more details can be found in (Kley et al. 2009). The influences of the distance of the embedded planet to the central star is discussed in Sect. 3. We then analyse the influence of the disc’s mass on the disc’s structure (density, temperature, aspect ratio) and then the influence on migration of embedded low-mass planets in Sect. 4. Convection in the disc is also briefly discussed in Sect. 4. We then summarize and conclude in Sect. 5.

2. Setup of the simulations

The protoplanetary disc is modelled as a three-dimensional (3D), non-self-gravitating gas whose motion is described by the Navier-Stokes equations. We treat the disc as a viscous medium, where the dissipative effects can then be described via the viscous stress-tensor approach. We also assume that the heating of the disc occurs solely through internal viscous dissipation and ignore the influence of additional energy sources (e.g. irradiation form the central star). This internally produced energy is then radiatively diffused through the disc and eventually emitted from its surface. For this process we use the flux-limited diffusion approximation, which allows us to treat the transition from optically thick to thin regions as an approximation. A more detailed description of the modelling and the numerical methodology is provided in our previous papers (Kley et al. 2009; Bitsch & Kley 2010; Bitsch & Kley 2011), and for that purpose we limit ourselves here to present only the most necessary information.

2.1. Physical setup

We solve the Navier-Stokes equations numerically using a spatially second order finite volume method based on the code NIRVANA (Ziegler & Yorke 1997), with implicit radiative transport in the flux-limited diffusion approximation and the FARGO (Masset 2000) extension as described in Kley et al. (2009). We use a spherical polar coordinate system (r,θ,φ), where the computational domain consists of a complete annulus (0 ≤ φ ≤ 2π) of the protoplanetary disc centred on the star, extending from r_min to r_max, which are determined by the location of the planet. For symmetry reasons and because we only use non-inclined planets, we restrict ourselves in the standard setup to the upper half of the disc. Hence, in the vertical direction the annulus extends from the equator up to 7° above the disc’s midplane, meaning 83° < θ < 90°. Here θ denotes the polar angle of our spherical polar coordinate system measured from the polar axis. For the study of convection we use in addition a two-sided disc, see below. The central star has one solar mass M_∗ = M_⊙, and the total disc mass inside [r_min,r_max] is M_disc = 0.01M_⊙, unless stated otherwise in Sect. 4. For our radiative models the temperature and subsequently H/r is calculated self-consistently from the equilibrium structure given by the viscous internal heating and radiative diffusion. We note that for given physics (equation of state, opacity, viscosity) the structure of the disc is solely dependent on its mass, and this is one aspect that we will investigate in this paper.

For the radiative transport we use a one-temperature approach and apply the flux-limited diffusion approximation using analytic Rosseland opacities, for details see Kley et al. (2009). To close the basic system of equations we use an ideal gas equation of state with constant mean molecular weight μ = 2.35 for a standard solar mixture. The adiabatic index is γ = 1.4. For the present study, we use a constant kinematic viscosity coefficient with a value of ν = 10¹⁵ cm²/s, a value that relates to an equivalent α = 0.004 at 5.2 AU for a disc aspect ratio of H/r = 0.05, where ν = αH²Ω_K. In standard dimensionless units with r₀ = a_Jup = 5.2 AU and t₀ = Ω_K(r₀)^-1 we have ν = 10^-5.

2.2. Calculation setup

Our coordinate system rotates at the initial orbital frequency of the planet (at r = r_P). We use an equidistant spherical grid in r,θ,φ with a resolution of (N_r,N_θ,N_φ) = (266,32,768) active cells for our simulations. At r_min and r_max we use damping boundary conditions for all three velocity components to minimize disturbances (wave reflections) at these boundaries. The velocities are relaxed towards their initial state on a time scale of approximately the local orbital period. The angular velocity is relaxed towards the Keplerian values, while the radial velocities at the inner and outer boundaries vanish. Reflecting boundary conditions are applied for the density and temperature in the radial directions. We apply periodic boundary conditions for all variables in the azimuthal direction. In the vertical direction we set outflow boundary conditions for θ_min = 83° (the surface of the computational domain).

In our previous work, we have discussed the calculation of the torque acting on a planet in great detail. Outward migration seems only possible and is strongest when the planet is on nearly circular orbits in the midplane of the disc (Bitsch & Kley 2010; Bitsch & Kley 2011). Additionally, we stated in Bitsch & Kley (2011) that the measured migration rate from planets on fixed orbits is equivalent to the migration rate determined from moving planets in discs. Hence, we consider here primarily planets on fixed circular orbits in the midplane of the disc, and calculate the torque acting on the planet, because the torque represents a direct measurement of migration in this case. For a one-sided disc we use symmetric boundaries at θ_max = 90° (the disc’s midplane). To correctly see the influence of convection in the disc we use two-sided discs for some high-mass discs. For these simulations we used outflow boundaries for both θ_min and θ_max.

The models are initialized with constant temperatures on cylinders with a profile T(s) ∝ s^-1 with s = rsinθ. The initial vertical density stratification is approximately given by a Gaussian where the radial surface density follows a Σ(r) ∝ r^−1/2 profile. In the radial and θ-direction we set the initial velocities to zero, while for the azimuthal component the initial velocity u_φ is given by the equilibrium of gravity, centrifugal acceleration, and the radial pressure gradient. This corresponds to the equilibrium configuration for a purely isothermal disc. Before embedding the planet, we bring the disc into radiative equilibrium by performing first 2D axisymmetric simulations in the r − θ plane. This takes about 100 orbits. We then extend this model to a full 3D simulation by expanding the grid into the φ-direction, and the planet is embedded.

For the gravitational potential of the planet we utilize the cubic potential, where the potential is exact beyond a transition (smoothing) radius r_sm and smoothed by a cubic polynomial inside (Klahr & Kley 2006; Kley et al. 2009). Here we use r_sm = 0.5R_H, where R_H is the Hill radius of the planet.

The torques acting on 20, 25, and 30 M_Earth planets are calculated to determine the direction of migration. The gravitational torques acting on the planet are calculated by integrating over the whole disc, where we apply a tapering function to exclude the inner parts of the Hill sphere of the planet (Crida et al. 2008). This torque-cutoff is necessary to avoid strong, probably noisy contributions from the inner parts of the Roche lobe and to disregard material that is gravitationally bound to the planet (Crida et al. 2009). Here we assume (as in our previous papers) a transition radius of 0.8 Hill radii.

3. Range of outward migration

Previous simulations by several authors (Baruteau & Masset 2008a; Paardekooper & Papaloizou 2008; Paardekooper & Mellema 2008; Kley & Crida 2008; Kley et al. 2009) indicated that outward migration of low-mass planets may be possible during an early evolutionary state of planet formation. However, because the simulations dealt mostly with planets at a given distance from the star, typically 5.2 AU, the radial range over which the migration may be directed outwards has not been addressed in great detail. In Bitsch & Kley (2010) we obtained some results for moving planets but the extent of the outward migration remained unclear.

In order to address this problem, we simulate 20,25, and 30 M_Earth planets on fixed circular orbits embedded in fully radiative discs at various distances from the host star. The planet’s semi-major axis r_p lies in a range of 0.6 ≤ r_P ≤ 5.0r₀, where r₀ = a_Jup = 5.2 AU. With increasing distance from the star, the density and temperature of the disc decrease and at some point in the disc the conditions for outward migration might not be fulfilled anymore.

For this set of simulations with different planet locations we use a disc setup with a density profile such that the total disc mass equals M_Disc = 0.01 M_⊙ for a planet at r_p = 1 and r_min = 0.4 and r_max = 2.5 in units of r₀. The planets are embedded in a way that the distance to the inner edge is always the planets location divided by 2.5, while the distance to the outer edge is the planet’s location to the star multiplied by 2.5. This setup ensures that the radial boundaries are always far enough away so they do not influence our results of embedded planets. Since the overall surface density profile (Σ ∝ r^-0.5) of the different disc models refers to the same physical situation, the total disc mass in the computational domain changes in the same way as the computed domain changes its size. The surface density at a given distance to the central star is constant for all disc models. The rotation frequency of the grid matches the rotation speed of the planet, so that the planet remains at a fixed position inside the computational grid at all times.

In Fig. 1 the specific torque (per planet mass) acting on the three planetary masses at different distance from the central star is displayed. For all planet masses the most extended positive torque (indicating outward migration) is around r ≈ 1.0 a_Jup. At longer distances to the central star, the torque acting on the planets decreases continuously to negative torques, and this transition from positive to negative torques occurs at larger distances for lower planet masses. For the lowest mass planet with 20 M_Earth the transition is at r ≈ 2.4 a_Jup (zero-torque distance to the central star). With even longer distances the torques remain negative but with diminishing strength, indicating inward migration. For higher planetary masses (25 and 30 M_Earth) the zero-torque distance is decreasing (1.9 and 1.4 a_Jup). For shorter distances to the central star, inside the maximum, the torque acting on the planet is reduced again until it reaches about zero for r_P = 0.5 a_Jup for all planetary masses.

In Fig. 2 the torque for the 20 M_Earth planet is again displayed together with semi-analytical results from Paardekooper et al. (2010, 2011). The black boxes refer to the torque formula presented in Paardekooper et al. (2010), which applies for the unsaturated torques in adiabatic discs. It is given by (1)with α and β being the slope of the radial surface density and midplane temperature profile, respectively (2)To calculate the semi-analytic results in the plot we use a constant α = 0.5, and β changes from about 1.7 at r = 1.0 a_Jup to 0.33 at r = 5.0 a_Jup, see also Fig. 7 below. The slope of the entropy profile is then given by ξ = β − (γ − 1)α = 1.5 at r = 1.0 a_Jup. The torque is normalized to

with q the planet/star mass ratio, h the relative disc height (H/r), Σ_P the surface density at the planet’s location and Ω_P the rotation frequency of the planet in the disc. The vertical height H of the disk is defined as H = c_s/Ω where the midplane values are used for the adiabatic sound speed c_s and the angular velocity Ω.

At r = 1.0 a_Jup the formula from Paardekooper et al. (2010) agrees well with our 3D simulations. As can be inferred from Eq. (1), the theoretical torque could never become negative for constant slopes (α,β,ξ). However, as β varies with r, a change from positive to negative torques is possible and indeed occurs at r ≈ 3.5 (black squares in Fig. 2). However, in the range of 1.2 < r < 3.8 the torque formula predicts a much higher torque than our 3D simulations, as is the case for larger distances to the central star, where the formula suggests a larger negative torque compared with our results. It should be noted that Eq. (1) includes Lindblad and corotation torques where the latter include contributions from the vorticity as well as entropy gradient. One should also be aware that Eq. (1) is valid only at the beginning of the evolution when the torques acting on the planet are not saturated. However, at later times the flow settles to an equilibrium state where the torques saturate and Eq. (1) is not valid any more.

	Fig. 1 Torque acting on planets with three different masses embedded in fully radiative discs located at various distances from the star (red for 20 MEarth, blue for 25 MEarth and purple for 30 MEarth). The torques are plotted when planet and disc have reached a quasi-stationary equilibrium where the torque acting on the planet is approximately constant in time.
Open with DEXTER

Fig. 2 Torque acting on a 20 MEarth planet in a fully radiative disc as a function of distance from the star (red plus signs). Additional curves indicate the results from the torque formula published in Paardekooper et al. (2010) (black squares) and Paardekooper et al. (2011) (purple crosses). Additionally we also plot the linear Lindblad torque by the Paardekooper et al. (2011) formula (blue asterisks).

Open with DEXTER

In Fig. 2 we additionally plot results from the improved formula in Paardekooper et al. (2011) (purple crosses). Its derivation is quite complex and we give a brief summary in the appendix. The formula captures the behaviour of the torque caused by Lindblad resonances and horseshoe drag on low-mass planets embedded in gaseous discs in the presence of viscous and thermal diffusion (Paardekooper et al. 2011). The new formula of Paardekooper et al. (2011) including diffusion provides a slightly better fit compared to the adiabatic Paardekooper et al. (2010) formula. The torque is positive for r < 1.8 a_Jup and becomes negative for larger distances to the central star. In general the formula captures the trend of our simulations, but lies consistently below our radiative disks. Also the zero-torque equilibrium radius differs. At the often used reference radius r =  a_Jup the torque formula and our 3D simulations differ by a factor of about three. For r > 2.5 a_Jup our simulations show a negative torque acting on the planet, which continues to grow until about r ≈ 4.0 a_Jup. For these larger radii the difference between our results and the analytic formula increases.

In Fig. 2 the Lindblad torque according to Paardekooper et al. (2011) is also displayed (blue). For larger distances to the star, the disc becomes thinner and the effects of the corotation torque become less important. The total torque should therefore converge towards the Lindblad torque, however the torques from our simulations differ by a factor of three at r = 5.0 a_Jup with the Lindblad torque of Paardekooper et al. (2011). As the Lindblad torque is dependent on the temperature gradient β, Eq. (1), and on a pre-factor, a reduction of the pre-factor as in Masset & Casoli (2010) reduces the Lindblad torque (see also Fig. 3).

	Fig. 3 Torque acting on a 20 MEarth planet in fully radiative discs as a function of distance from the star. The additional curves indicate the estimates from the theoretical formula for viscous, diffusive disks by Masset & Casoli (2010) and for isothermal disks by D’Angelo & Lubow (2010).
Open with DEXTER

In Fig. 3 we also plot the results from Masset & Casoli (2010), who provided an alternative recipe for planetary migration in viscous discs with thermal diffusion. We do not cite the full formulae here, due to their complexity. In paragraph 8 of Masset & Casoli (2010), a summary of the torque formula is given. Additionally, the isothermal results of the 3D simulations of D’Angelo & Lubow (2010) are plotted.

Compared to our 3D results the torque formula of Masset & Casoli (2010) leads to qualitatively different results. The torque predicted by their formula is negative for all investigated distances to the central star. For large distances the torque of Masset & Casoli (2010) matches better with our simulations, but is still about 35% too low. The Lindblad torque of Masset & Casoli (2010) is smaller compared to the Lindblad torque in Paardekooper et al. (2010) due to a different pre-factor in front of the temperature slope β (in Eq. (1)). The 3D isothermal formula by D’Angelo & Lubow (2010) shows a larger Lindblad torque compared to our simulations and to Masset & Casoli (2010). But as the formula is for isothermal discs only, the torque has to be reduced by a factor of γ compared to our radiative simulations or compared to the torque from Masset & Casoli (2010). Taking this into account the torques match quite well for large distances to the central star.

To test our implementation of the torque formula by Masset & Casoli (2010), we applied it to the results of Ayliffe & Bate (2011) using our temperature gradient (see their Fig. 3). However, in contrast to Ayliffe & Bate (2011) we find negative torques when applying that formula. Apparently, when comparing their numerical results to the anytical formula by Masset & Casoli (2010), Ayliffe & Bate applied the latter with a somewhat larger viscosity (approx. by a factor of 2π) than quoted in text. This resulted in the positive torque quoted in their paper (B. Ayliffe, priv. comm.).

There are many reasons for the differences between the simulations and the formulae. The formulae in Paardekooper et al. (2010, 2011) and in Masset & Casoli (2010) were derived for 2D discs only, but the horseshoe drag can be even stronger in three dimensions, as shown for isothermal discs in Masset et al. (2006). The radiative diffusion is also most effective in the vertical direction, meaning that the two-dimensional approximation where diffusion is is restricted to the orbital plane, does not capture the true physical effects (see also Kley & Crida 2008). The formulae were derived for background gradients of temperature and surface density, but as the disc with an embedded planet evolves, the structure of the disc changes and the basic assumptions (gradients in temperature, density, and so on) used to derive the formulae might not be valid any more. This is in particular true for planets in the mass range studied here, because the theory has been developed for lower mass planets, which do not alter their surroundings significantly. A more detailed discussion about the smoothing of the planet in Paardekooper et al. (2011) can be found in Appendix A. Nevertheless, there seems to be a very rough qualitative agreement between the improved torque formula of Paardekooper et al. (2011) and our numerical results.

In Fig. 4 the radial torque density Γ(r) is displayed for 20 M_Earth planets at r_P = 0.6, r_P = 1.0, r_P = 2.5 and r_P = 4.0 a_Jup, where the radius is normalized to the actual planetary distances to allow a direct comparison. In all curves the contributions by the Lindblad torques are clearly visible, positive for r < 1 and negative for r > 1. The contribution responsible for the torque reversal, the “spike” just inside r = 1, is visible only for the r_P = 0.6 and r_P = 1.0 locations, which both show an overall positive total torque. Even though the torque acting on the r_P = 0.6 a_Jup planet is much reduced, the spike in the torque distribution is clearly visible. In Kley et al. (2009) we discussed the origin of the spike for the r_P = 1.0 a_Jup planet. It is an indicator for a density enhancement ahead of the planet just inside of r_P, and thus creating a positive torque.

For the r_P = 2.5 a_Jup case the corotation spike in the fully radiative case in the torque density is not visible any more, only the Lindblad torques are visible. The resulting total torque is about zero, which indicates that the Lindblad torques are just counterbalanced by the corotation effects. For even longer distances to the central star, the torque is identical to the (negative) Lindblad torque, indicating inward migration (see also Fig. 1).

	Fig. 4 Radial torque density Γ(r) acting on 20 MEarth planets at rP = 0.6, rP = 1.0, rP = 2.5 and rP = 4.0 aJup. The distances have been scaled in r-direction to the actual planet location rP. The snapshots of the torque density have been taken in the equilibrium state at 80 Jupiter orbits.
Open with DEXTER

In Kley et al. (2009) we explained in great detail with the help of 2D surface density plots how the torque acting on the planet is created in fully radiative discs. The torque acting on a 20 M_Earth planet on a circular orbit at r_P = a_Jup in a fully radiative disc is positive, indicating outward migration. In Fig. 5 the 2D surface density of the 20 M_Earth planets at r = 0.6,1.0,2.5,4.0 a_Jup (from top to bottom) is displayed. The origin of the structure of the standard r_P = 1.0 a_Jup case (second from top) was described in Kley et al. (2009), and we display this case here again for better reference.

The planet located at r_P = 0.6 a_Jup (top panel) is still prone to outward migration (see positive torque in Fig. 1), but at a slower rate. The surface density distribution shows some significant differences compared to the r_P = 1.0 a_Jup planet. The density increase ahead and inside of the planet (φ > 180° and r < 0.6) is still visible in the r_P = 0.6 a_Jup case, but the density decrease behind the planet (φ < 180° and r > 0.6) is not that clearly visible. Indeed it seems as if the density behind the planet increased in a way that the total torque acting on the planet is reduced dramatically, resulting in reduced positive torque acting on the planet.

For the planet at r_P = 2.5 a_Jup, where the total torque acting on the planet is about zero, the density increase ahead of the planet is no longer visible in the surface density plot, but the decrease behind the planet is clearly visible. It also seems that the planet is able to deplete a larger area around it, possibly owing to the onset of gap formation. For even longer distances to the central star (e.g. r_P = 4.0 a_Jup), the effect becomes ever stronger, and the gap is more pronounced. The density increase in front of the planet is no longer visible at all. The torque acting on this planet is clearly negative, indicating inward migration. With increasing distance from the host star, the opening angle δ of the spiral arms becomes smaller (see Fig. 5), as can be inferred from δ ≈ c_s/v_Kep which scales  ∝ r^-0.35 for our temperature gradient of T(r) ∝ r^-1.7.

	Fig. 5 Surface density maps for a 20 MEarth planet on fixed circular orbit in fully radiative discs at four different locations. The distance from star to planet changes (from top to bottom): rP = 0.6, rP = 1.0, rP = 2.5 and rP = 4.0 (in Jupiter radii). Please note the slightly different colour scale for each plot.
Open with DEXTER

To analyse the extent of outward migration from our disc properties in more detail, it seems useful to compare various time scales: the libration, the radiative and the viscous time scale in the disc. The necessary unsaturated torques needed for sustained outward migration require approximately equal libration and radiative diffusion times (see Baruteau & Masset 2008a; Paardekooper & Papaloizou 2008). For the latter we use in our case the radiative diffusion time scale, τ_rad. We define (Bitsch & Kley 2010) (3)with the diffusion coefficient (4)For the typical diffusion length s we substitute the vertical disc height, i.e. s = H, where we use H = c_s/Ω with the sound speed c_s. The libration time given by (Baruteau & Masset 2008a) (5)where x_s denotes the half-width of the horseshoe-orbit, . Similarly to the radiative diffusion time, the viscous time scale is given by τ_visc = s²/ν, again with s = H, and a constant ν. To compute the time scales all required quantities are evaluated in the midplane of the unperturbed disc at the beginning of the simulations. This applies to the density, temperature, opacity κ(ρ,T), and the sound speed and Ω.

The three time scales are displayed in Fig. 6. For accretion discs that are solely heated internally by viscous dissipation we expect in equilibrium τ_rad ≈ τ_visc. Apparently, this relation is fulfilled well (for r < 1.5 a_Jup). We have plotted the libration time for two different planet masses, 20 and 30 M_Earth. Time scale arguments suggest a most efficient outward migration for equal τ_lib and τ_rad, which is indeed roughly what we find in our 3D simulations. However, the overall range of outward migration is surprisingly broad. Specifically we find that 20 M_Earth planets are prone to outward migration up to about r ≈ 2.4, where the two time scales differ by a factor of 3–4. For 30 M_Earth planets the range of outward migration is substantially smaller and centres directly around equal libration and radiative time scale.

We have checked the above estimates of the torques for a stationary planet with additional simulations of 20 M_Earth planets in discs, starting at r = 2.0 a_Jup and r = 3.0 a_Jup respectively, which were able to move freely inside the disc. The planets gather in this case indeed at the zero torque radius (results are not displayed here).

Fig. 6 Radiative and viscous diffusion time scales that depend on the distance from the central star for our standard disc model. To compute the time scales we used the density and temperature of the midplane at the beginning of the fully radiative simulations (when the disc is in the r − θ equilibrium between viscous heating and radiative transport/cooling). Libration time scales are stated for 20 and 30 MEarth planets.

Open with DEXTER

4. Influence of the disc’s mass

In this section we examine the influence of the disc’s mass on the migration of low-mass planets embedded in these discs. First we compare the relevant physical properties of the discs with different masses, and then we investigate the planetary migration in those discs. We then finally discuss convection inside fully radiative discs.

4.1. Properties of discs with different masses

In our previous work, the disc’s mass was fixed to 0.01 M_⊙. We now extend the range of disc masses from 0.005 M_⊙ to 0.04 M_⊙ (with respect to the standard radial distance, from 0.4−2.5). All models started locally isothermal with H/r = 0.05 and during initial evolution on time this will change to the appropriate equilibrium configurations (between viscous heating and radiative transport/cooling).

In Fig. 7 we display the density, temperature, and the aspect ratio of the equilibrium discs for different disc masses at r = 1.0. Density and temperature, and H/r are evaluated in the disc’s midplane. The results are as expected from our previous simulations. A higher mass of the disc results in a higher density, temperature, and aspect ratio of the disc in the equilibrium state. In the isothermal case, a higher aspect ratio of the disc would result in a slower inward migration of a low-mass planet, see Tanaka et al. (2002) for linear analysis and e.g. Bitsch & Kley (2010) for non-linear simulations. However, low-mass planets in fully radiative discs migrate outwards, so that the linear isothermal approach is not valid any more.

	Fig. 7 Density (top), temperature (middle) and aspect ratio (bottom) of discs with different masses in the initial equilibrium state. All quantities are measured in the disc’s midplane at the reference distance, r = 1.0.
Open with DEXTER

In Fig. 8 the radial distributions of surface density (top) and temperature (bottom) are displayed. The temperature has been measured in the disc’s midplane. By construction, the surface density increases for higher disc masses, while it falls off with increasing distance to the star, on average according to Σ(r) ∝ r^−1/2 as expected for a constant ν. The surface density profiles for the higher mass discs with M_Disc ≥ 0.015 M_⊙ show some fluctuations. With increasing disc mass, these fluctuations become stronger and reach out to a longer distance from the star. While they are quite short for M_Disc ≤ 0.02 M_⊙ and reach only to r ≈ 1.3, they become very strong and reach out to r ≈ 2.3 for M_Disc = 0.04 M_⊙. These fluctuations of the surface density vary in time and are related through convective motions in the disc, see below.

The described fluctuations in the surface density can also be seen in the temperature profiles of discs with different disc masses (bottom panel in Fig. 8). The variabilities of the temperature are not as strong as those for the surface density, nevertheless, they are clearly notable for M_Disc ≥ 0.015 M_⊙ and increase with the disc’s mass. They also, change in time, as does the surface density. A higher disc mass seems to support stronger fluctuations that reach farther out into the disc.

	Fig. 8 Surface density (top) and temperature (bottom) of discs with different masses in the equilibrium state of the disc. The temperature is measured in the disc’s midplane.
Open with DEXTER

Because a higher disc mass results in a higher aspect ratio of the disc, the cut-off of the computed domain (at 7° above and below the midplane) might change the structure of the disc. Additional simulations with larger θ did not change the density and temperature patterns at θ = 83° and θ = 97° for all disc masses (see also Kley et al. 2009). However, a too low boundary substantially changes the distributions and might therefore influence convection in the disc as well.

The changes in the surface density profiles have a direct influence on the migration of an embedded planet. If a planet is embedded in a region in the disc where the fluctuations of the surface density are very strong (i.e. strong convection), the direction of migration might not be clearly determinable, because changes in the surface density profile directly influence the migration. Stationary gradients in surface density profiles can even be used as planetary traps to collect planetary embryos (Morbidelli et al. 2008).

4.2. Influences on the migration of low-mass planets

As mentioned above, a different aspect ratio of the disc will change the rate of planetary migration. In the isothermal case, a higher aspect ratio will result in slower inward migration in theory (Tanaka et al. 2002), which we supported in our previous simulations (Bitsch & Kley 2010; Bitsch & Kley 2011). A higher disc mass results in a higher aspect ratio and temperature of the disc, and we may expect changes in the migration rates. When the planet is farther away from the central star, the temperature and density of the disc are reduced. If the reduction in density and temperature is sufficient, the planets stop their outward migration (see Fig. 1). One might now expect that for discs with very low masses the torque acting on a 20 M_Earth planet at r = 1.0 a_Jup might become negative. Very high temperatures and high densities inside the disc, on the other hand, might influence the outward migration as well. Following the formula of i¸tet2010MNRAS.401.1950P in Eq. (1), one might suspect stronger torques acting on planets in more massive discs for constant α, β and ξ (increase in surface density overcompensates the increase in aspect ratio).

In the top panel of Fig. 9 the total torque acting on 20 M_Earth planets on circular orbits embedded in fully radiative discs with different masses and the theoretical results from Paardekooper et al. (2010) and Paardekooper et al. (2011) are displayed (blue and purple). The torque acting on the planet remains nearly constant within a small interval around our standard disc mass of 0.01 M_⊙. For lower disc masses, the torque drops off very rapidly to even negative values for M_Disc = 0.005 M_⊙, as we expected. For higher disc masses, the torque drops off as well. First at a faster rate (to  ≈0.020 M_⊙), then at a slightly slower rate, until it reaches an about zero torque state for M_Disc = 0.040 M_⊙. This contradicts to our first expectation that planets in more massive discs should experience a higher torque, the reason may be a change in the temperature gradient and the influence of convection.

When looking at the surface density profile displayed in Fig. 8, it is clear that the changes in the surface density may disrupt the very sensitive density pattern near the planet. As the convection cells in the disc change with time, the torque acting on the planet will change as well, giving rise to high fluctuations/oscillations in the total torque acting on the planet. Hence, the torques acting on the planet have been averaged over 20 planetary orbits. After averaging, the torques acting on planets in convective discs show only very low fluctuations.

For disc masses around  ≈0.01 M_⊙ the theoretical formula from Paardekooper et al. (2010, see Eq. (1)) fits our 3D simulations up to a factor of 25%. For very low disc masses (M_disc = 0.005 M_⊙), however, the fit is not as good. This may be because of the reduced disc mass and the consequently changed surface density distribution (which changes the torque acting on the planet), as explained in Sect. 3. For higher disc masses (M_disc ≈ 0.015 M_⊙), the two torque values differ more and more. As the torques of our simulations tend to go to zero, the theoretically predicted adiabatic torques from Paardekooper et al. (2010) become even more extended.

The formula from Paardekooper et al. (2011), which includes the important effects of viscosity and heat diffusion, differs by a factor of three near disc masses of M_disc ≈ 0.01 M_⊙ as one could have expected from the results presented in Fig. 2. For higher disc masses, the torques from Paardekooper et al. (2011) remain nearly unchanged. The torques from our simulations are reduced, but they do not reach negative values. Besides the differences described in Sect. 3, more differences arise from convection in the disc, because convection results in fluctuations in the torque, which is not considered in the analytical formulae. However, convection seems to play a role only for discs with M_disc > 0.02 M_⊙.

4.3. Torque analysis

In the bottom panel of Fig. 9 the radial torque density Γ(r) is displayed for different disc masses. For the lowest disc mass in our simulations, M_Disc = 0.005 M_⊙, the usual spike in the torque density cannot be seen. The spike in the torque density distribution is an indication for a positive torque in a fully radiative disc (Kley et al. 2009). For higher disc masses (up to M_Disc ≈ 0.02 M_⊙), it is is clearly visible. For those disc masses, the total torque is indeed positive (see top panel of Fig. 9), indicating outward migration of the embedded protoplanet.

Fig. 9 Torque acting on a planet located at 5 AU for different disc masses. Top: specific total torque acting on planets (20 MEarth) embedded in discs with different disc masses. The planets embedded in the higher mass discs MDisc > 0.02 M⊙ are in the convective zone in the disc, so that the torque acting on the planet is very noisy and has been averaged over 20 planetary orbits. Additionally, we over-plotted results (blue and purple) from the theoretical torque formulae of Paardekooper et al. (2010, 2011) for 20 MEarth planets. Bottom: radial torque density Γ(r) acting on the planet for different disc masses. For comparison, the vertical line indicates the location of the maximum as found for our standard case.

Open with DEXTER

The torque density for the M_Disc = 0.04 M_⊙ disc seems to indicate a total positive torque acting on the planet, and it is indeed positive at the moment of the snapshot, but as the fluctuations in the surface density change in time, so does the torque acting on an embedded planet. Therefore the torque density for the M_Disc = 0.04 M_⊙ disc at one single moment during the evolution does not necessarily reflect the longterm outcome, if the fluctuations are to strong. This time variation of the total torque and torque density acting on the planet is displayed in Fig. 10. The top panel in Fig. 10 shows the time evolution of the total torque acting on embedded planets for discs with different disc masses. For M_disc < 0.02 M_⊙ the torque acting on the planet is constant in time (after about 50 orbits), while it shows very high fluctuations for M_disc = 0.04 M_⊙. In the bottom panel of Fig. 10 the torque density Γ(r) for a 20 M_Earth planet in a M_disc = 0.04 M_⊙ disc is shown. The torque density is plott ed at different times. Because the total torque was fluctuating very much, it is no surprise to find these fluctuations for the torque density as well. These fluctuations, induced by convection, clearly show that a higher disc mass disturbs the evolution of the torque.

Fig. 10 Top: specific total torque evolution acting on planets (20 MEarth) embedded in discs with different disc masses. The planets embedded in the higher mass discs MDisc > 0.02 M⊙ are in the convective zone in the disc, so that the torque acting on the planet is very noisy compared to the low-mass discs. Bottom: radial torque density Γ(r) acting on the planet embedded in a disc with Mdisc = 0.04 M⊙ at different stages of the evolution.

Open with DEXTER

In Fig. 11 the surface density profiles of fully radiative discs with disc masses of 0.005, 0.015 and 0.04 M_⊙ with embedded 20 M_Earth planets are displayed. The planet in the 0.005 M_⊙ disc generates a very similar surface density structure compared to our standard 0.01 M_⊙ disc (second panel in Fig. 5). The overall density is, of course, reduced (because the disc mass is much lower), but the general pattern of the density increases in front of the planet φ > 180° and r < 1.0 and the decrease behind the planet φ < 180° and r > 1.0 remains intact. However, the density structure relevant for outward migration is not as pronounced as it should be to result in a positive torque acting on the planet.

For planets embedded in discs with higher masses, the picture is quite different. For a disc mass of 0.015 M_⊙ (middle picture in Fig. 11) the density structure in the direction of the star (r < 0.9 a_Jup) is very distorted, but one can still see the density increase ahead and the density decrease behind the planet. The distortion seems to reduce the torque acting on the planet, but the overall torque is still positive, indicating outward migration. For an even higher disc mass (bottom picture in Fig. 11 with M_Disc = 0.040 M_⊙) the distortions in the disc increase more. The density structure, normally seen for low-mass planets in fully radiative discs, is no longer visible at all. The distortions are so strong that the torque acting on the planet becomes about zero, indicating only a low migration rate.

However, in this last case the torque is in a state of constant fluctuations, which complicates realistic predictions about the direction of migration in these massive discs. The fluctuations of the torque have their cause in fluctuations of the density patterns, which indicates that the convective zone inside the disc is enlarged compared to low-massive discs. We observed the phenomenon of convection briefly in our previous work (Kley et al. 2009) for discs with M_Disc = 0.01 M_⊙ as well, but the convective zone did not reach the planet, and thus did not disturb the density pattern around the planet.

	Fig. 11 Surface density maps for planets on fixed circular orbits in fully radiative discs with different disc masses (from top to bottom): MDisc = 0.005 M⊙, MDisc = 0.015 M⊙ and MDisc = 0.040 M⊙. The disruption in the surface density patterns for the higher mass discs are caused by convection inside the disc.
Open with DEXTER

In Fig. 12 we display the radiative diffusion time scale and the libration time scale for discs with different masses. In order to keep the torques acting on the planet unsaturated (to evoke outward migration), the libration and radiative diffusion time scale need to be approximately equal (Baruteau & Masset 2008a; Paardekooper & Papaloizou 2008). However, for convective discs the equilibrium state through cooling through convection is reached only when τ_rad > τ_visc as observed in Fig. 12. The time scales for the 0.005 M_⊙ disc indicate outward migration in a region at r ≈ 2.0 a_Jup, but at r = 1.0 a_Jup the time scales differ by a factor of 4, indicating that the torques are not kept unsaturated in this region of the disc, which in turn indicates inward migration (as presented in in the top figure in Fig. 9).

For the 0.015 M_⊙ disc the time scales are nearly identical at r = 1.0 a_Jup, which indicates outward migration (as can be seen in the top of Fig. 9). However, for longer distances to the central star, the time scales start to differ, which indicates inward migration. In the 0.040 M_⊙ disc, the time scales differ by a factor of 3 at r = 1.0 a_Jup, which indicates inward migration. However, the measured torque acting on the planet is positive, indicating slow outward migration. But because the planet is embedded in a highly convective region in the disc, it is very difficult to predict the motion of the planet correctly by considering only the time scales.

	Fig. 12 Time scales depending on the distance from the central star for 0.005 M⊙, 0.010 M⊙, 0.015 M⊙ and 0.040 M⊙ discs. To compute the time scales we used the density and temperature of the midplane at the beginning of the fully radiative simulations (when the disc is in the r − θ equilibrium between viscous heating and radiative transport/cooling).
Open with DEXTER

4.4. Orbital evolution

In Fig. 13 the evolution of semi-major axis for 20 M_Earth planets in isothermal and fully radiative discs with different disc masses is displayed. The isothermal reference simulations are started from a H/r = 0.037 disc, which represents the H/r value at the planets starting location in the fully radiative regime. In the isothermal disc, no convection is present and therefore the embedded planet migrates as expected. A higher disc mass results in a faster inward migration. However, it seems that for the M_disc = 0.04 M_⊙ disc the type-III-migration regime is hit, because the planet moves inwards very fast.

For a planet on a fixed circular orbit in a fully radiative disc with M_disc = 0.04 M_⊙ we determined a positive torque (see Fig. 9) by averaging in time. However, the total torque acting on the planet undergoes strong fluctuations in time (see Fig. 10). When embedding a 20 M_Earth planet in such a highly convective disc, the evolution pattern should to some extend reflect the strong fluctuations in the torque acting on a planet on a fixed orbit. Indeed, the planet experiences some small kicks in its evolution pattern (see Fig. 13). Interestingly, the planet migrates inwards despite the positive torque acting on a fixed planet. It seems that in the convective region of high-mass discs, the measurement of the torque for planets on fixed orbits is not as reliable as for planets in low-mass discs. Additionally, the time averaging may have to be performed over a longer time span.

	Fig. 13 Evolution of semi-major axis for 20 MEarth planets (starting at r = 1.0 aJup) in isothermal and fully radiative discs with different disc masses.
Open with DEXTER

Fig. 14 Vertical velocities in z-direction for planets on fixed circular orbits in fully radiative discs with different disc masses (from top to bottom) in the disc’s midplane: MDisc = 0.005 M⊙, MDisc = 0.015 M⊙ and MDisc = 0.04 M⊙. The disruption in the velocity patterns for the higher mass discs are caused by convection inside the disc. A positive velocity simulates outward flows, while a negative velocity indicates a flow towards the disc’s midplane.

Open with DEXTER

4.5. Convective zone

In order to investigate whether convection is actually present in the disc, we display the velocities in z-direction (out of midplane of the disc) for different disc masses (M_Disc = 0.005 M_⊙, M_Disc = 0.015 M_⊙ and M_Disc = 0.04 M_⊙) in the disc’s midplane in Fig. 14. These plots represent simulations with a whole disc, meaning 83° ≤ θ ≤ 97°. As the surface density plots indicated, there are no disrupted areas in the velocity patters for the M_Disc = 0.005 M_⊙ simulations. Therefore for this low disc mass we observe no convection near the planet.

However, for the M_Disc = 0.015 M_⊙ case the surface density patterns already indicated that convection is possible in the disc inside of the planet’s distance to the central star. The velocity distributions confirm this result, lines with positive velocities (indicating a flow towards the upper boundary of the disc) alternate with lines with negative velocities (indicating a flow towards the midplane of the disc). These flows are a clear indicator of motion caused by convection in the disc. In the M_Disc = 0.040 M_⊙ case, the fluctuations in the surface density increased dramatically, as did the changes in the velocity pattern. Alternations between positive and negative velocities have increased and indicate a very strong convective region that disturbs the torque acting on the planet (and therefore it’s migration). The flow pattern in the disc is very erratic, making it absolutely necessary to average the torque acting on an embedded planet for many orbits.

For simulations that cover only the upper half of the disc, the convection cells inside the disc end at the disc’s midplane, but in reality these convective cells continue to the lower half of the disc. However, that convection inside the disc changes the behaviour of the planet disc interactions can also be seen for simulations containing only one half of the disc.

In Fig. 15 the velocity in z direction is displayed for a half-size disc (only the upper half of the disc is computed, top figure) and for a full disc. These computations have been performed only for fully radiative axisymmetric 2D discs (in r − θ direction) with a disc mass of M_disc = 0.04 M_⊙ without an embedded planet. The convection in the disc is clearly visible. In both simulations the convection cells in the disc become more symmetric for distances longer than r > 1.25 a_Jup, meaning that the velocity changes from positive to negative only in radial direction and not in the vertical direction as well. For shorter distances to the central star the convection cells are very irregular in both cases.

	Fig. 15 Velocities in z-direction for 0.04 M⊙ discs without planet. In the top panel only the upper half of the disc was simulated, while in the bottom panel both sides of the disc were simulated (with twice the number of grid cells in θ direction). The other simulation parameters are identical.
Open with DEXTER

When putting a 20 M_Earth planet in the 0.04 M_⊙ discs, the structure of the convective cells changes as the disc gets disturbed by the planet. Without the planet, the convective cells are very regular for distances to the central star exceeding r > 1.25 a_Jup. With an embedded planet these regular structures break down and become very irregular, as can be seen in Fig. 16. This effect may caused by to the wake created by the planet inside the disc, which acts as an additional heat source. At shorter distances to the central star, the structure is irregular with or without an embedded planet.

	Fig. 16 Velocities in z-direction for 0.04 M⊙ discs with embedded 20 MEarth planet. In the top panel only the upper half of the disc was simulated, while in the bottom panel both sides of the disc were simulated (with twice the number of grid cells in θ direction). The other simulation parameters are identical.
Open with DEXTER

The velocity patterns for the one sided and two sided discs show a little difference. It seems that the fluctuations in velocity are more centred in midplane for the one sided disc, while they seem to be located near the upper and lower boundaries for the two sided discs. This effect has several reasons. In the one sided disc, material is reflected at the midplane of the disc, which might lead to an increase of fluctuations near the midplane. In the two sided disc, material can flow through the midplane, so that the fluctuations near the midplane are reduced. Furthermore, the one sided disc might be unrealistic if convection is present in the disc.

However, the general structure changes when both halves of the disc are computed, independent of an embedded planet. The convection cells are now moving through the midplane of the disc, which was not possible for simulations of only one half of the disc, see also Kley et al. (1993); Klahr et al. (1999). Therefore, the surface density structure in the two sided disc (with M_disc = 0.04 M_⊙) is slightly different compared to the one sided disc. In the two sided case, the fluctuations in the surface density continue only to  ≈1.5 a_Jup, while they covered the whole disc in the one sided case. The temperature profiles, on the other hand, show no difference at all. We therefore only expect little change in the torque acting on planets embedded in one or two sided discs at r_P = 1.0 a_Jup.

Simulations of embedded 20 M_Earth planets in fully radiative discs that cover both sides of the disc show only very small differences in the velocity pattern in mid-plane of the disc. The fluctuations in time of the torque acting on the planet embedded in M_disc = 0.04 M_⊙ and M_disc = 0.03 M_⊙ discs are comparable (with only small differences in the amplitude of the fluctuations) for both simulations, confirming our previous assumptions. Simulations of planets embedded in lower mass discs show no difference at all (simulations not displayed here), because the convective region does not reach to the planet at all.

Considering that a longer distance from the central star resulted in a turn from outward to inward migration for a 20 M_Earth planet (see Fig. 1) because of a reduction in temperature and density at the given location of the planet, one might argue that the zero-torque distance from the central star might be at longer distances for higher disc masses. Moreover, as the disrupting convective zone in massive discs reaches farther out from the star, outward migration might be possible at larger radii in more massive discs, because the disc’s convective zone stops at longer distance to the central star and the density is still high enough to produce the surface density patter needed for outward migration.

Additional simulations with 20 M_Earth planets in 0.02, 0.025, 0.03 and 0.04 M_⊙ discs with r_p = 2.0, 2.5, 3.0 and 4.0 a_Jup, respectively, confirm our assumption (see Fig. 17). The torques acting on those planets are positive, indicating outward migration, and show no fluctuations in time. The surface density plots also show no sign of convection in the disc at the location of the planet (not displayed here). It seems that outward migration is therefore possible to farther distances from the central star in more massive discs.

	Fig. 17 Specific torque acting on 20 MEarth planets embedded in 0.02, 0.025, 0.03 and 0.04 M⊙ discs at the following distances rp = 2.0, 2.5, 3.0 and 4.0aJup, respectively.
Open with DEXTER

The picture of convection in our disc would change when including stellar irradiation because it would heat the surfaces of the disc in contrast to the applied cooling right now. This would result in less convection in the disc. Because the convective region is a result of the higher surface density (increasing τ_rad) and viscosity in the disc, a reduction in viscosity could prevent convection in the disc. However, a reduction of viscosity also reduces the torque of an embedded planet, so that outward migration might not be possible any more, even for low-mass discs. The influence of viscosity will also be addressed in a next paper in much more detail.

In self-gravitating discs, the torque acting on an embedded planet can differ by a factor of two compared to non self-gravitating discs, as shown by Baruteau & Masset (2008b). These authors also state that self-gravity has no effect on the corotation torque in the linear regime, but our 3D simulations are in the non-linear regime. Therefore the influence of self-gravity on planet migration in fully radiative discs should be investigated in the future.

The Toomre stability criterion can be used to estimate the stability of discs against self-gravity (Toomre 1964). The stability parameter reads (6)where c_s is the sound speed in the disc, κ_ep is the epicyclic frequency, which for Keplerian discs is approximately equal to the angular frequency Ω, Σ is the surface mass density and G is the gravitational constant. In order to achieve stability in discs, the stability parameter must be Q ≫ 1. For all disc masses used in this work, this criterion is fulfilled well, so that the discs are not gravitationally unstable.

Because convection is a 3D effect, 2D simulations (in r-φ direction in the midplane) of fully radiative discs with high discs masses (M_disc > 0.02 M_⊙) cannot capture this effect. Therefore planets embedded in these 2D simulations will not be exposed to these fluctuations and might therefore be inaccurate near the central star because of convection in the disc.

5. Summary and conclusions

We performed full 3D radiation hydrodynamical simulations of low-mass planets embedded in accretion discs at different distances to the central star and for various disc masses.

In the first sequence of our simulations we changed the planetary distance to the central star of embedded planets on circular orbits. With increasing distance to the central star, the torque acting on 20 M_Earth planets embedded in fully radiative discs becomes ever more reduced and it reaches negative torques for longer distances. We find an equilibrium, zero-torque distance, to the central star for 20 M_Earth planets at r ≈ 2.4 a_Jup. This equilibrium distance varies with the planetary mass (for 25 M_Earth planets it is r ≈ 1.9 a_Jup and r ≈ 1.4 a_Jup for 30 M_Earth planets), indicating that a quite extended region in the disc might act as a feeding zone to create even larger planetary cores. The concept of equilibrium radius (zero torque radius) for planetary embryos in fully unsaturated discs has been stated in Lyra et al. (2010) as well and it could easily act as a feeding or collection zone for planetary embryos.

Planets embedded in fully radiative discs migrating outwards create a very sensitive pattern in the surface density distribution. Ahead and inside of the planet a density increase is visible (see the second panel from the top in Fig. 5), which shrinks with increasing distance to the central star. This density enhancement is indeed accountable for the positive torque acting on the planet (also visible as a spike in the radial torque density distribution in Fig. 4), but as the distance to the central star increases, this effect becomes less, so that it cannot overcompensate the negative Lindblad torques any more, which results in inward migration.

We compared our results to the recently developed torque formulae by Paardekooper et al. (2010, 2011) and Masset & Casoli (2010). Paardekooper et al. (2010) includes just the fully unsaturated torques in the inviscid and adiabatic case, shows no torque reversal option for constant β and is as such unphysical when comparing it with the long-term evolution of planets. However, as in our simulations β changes in the disc with distance to the central star (β(r = 1.0 a_Jup) = 1.7 to β(r = 5.0 a_Jup) = 0.33), the torque given by the formula becomes negative for large distances to the central star. But, overall the match of the torque of the formula with our simulations is not good. This formula is only valid in the first orbits after the planet is embedded when the torques are still unsaturated, however, the torques do saturate in time.

However, as expected, the improved version of Paardekooper et al. (2011) that includes viscous and heat diffusion describes our results more accurately. It also features a transition from positive to negative torques, but the negative torques at large distances to the central star are about a factor of two larger than in our simulations. Even though an exact match of the torques has not been achieved, the trend seems to be caputured quite well by the formula.

The torque formulae are derived for 2D discs with thermal diffusion operating only in the disk’s midplane, while our simulations are 3D, and take full account of vertical diffusion as well, which can change the structure of the disc. The formulae were also derived for given gradients in temperature and surface density, but in a real disc the temperature and surface density profiles are disturbed when a planet is embedded in a disc. Because the formulae were derived and checked for a 5 M_Earth planet (with about 20 per cent agreement), the disturbances of such a small planet in a disc are much weaker than for our embedded 20 M_Earth planet, which may give rise to more significant, non-linear disturbances in the temperature and density profiles. All of this may lead to differences between our simulations and the theoretical formulae in the torque acting on the planet.

The formula of Masset & Casoli (2010) does not match our simulations that well. The torques from the formula are always negative, which is in contrast to our simulations, where we find outward migration for r < 2.5 a_Jup. For long distances to the central star, Masset & Casoli (2010) match better than Paardekooper et al. (2011), as the Lindblad torque is reduced. The Lindblad torque of Masset & Casoli (2010) also matches quite well with the isothermal torque of D’Angelo & Lubow (2010), when taking the factor γ into account due to the different sound speeds in isothermal and fully radiative discs. The overall trend to have torques from positive to negative seems best achieved with Paardekooper et al. (2011), but a quantitative agreement is still lacking. In Appendix A we discuss the influence of the smoothing length on the formulae.

For increasing disc masses, the temperature, density (in midplane), and aspect ratio of the disc increases in the equilibrium state where viscous heating and radiative transport/cooling are in balance. The convective zone in the inner discs stretches farther out with increasing disc mass, resulting in high fluctuations of the surface density in our computed domain for discs with a mass higher than M_disc ≈ 0.02 M_⊙.

Starting from a M_disc = 0.01 M_⊙ disc, the torque acting on embedded 20 M_Earth planets decreases for increasing disc masses. As the disc mass increases, the convective zone in the disc stretches farther out from the central star and influences planetary migration. The fluctuations in the disc’s density disrupt the torque acting on the planet on a stationary orbit for high-mass discs in a way that the torque is very irregular and shows high fluctuations as well, making it difficult to determine the correct direction of migration. For lower disc masses, the torque reduces as well, presumably because of the same reasons as the torque is reduced for longer distances to the central star in a discs with M_dics = 0.01 M_⊙.

The formula in Paardekooper et al. (2011) fits within a factor of three with our 3D simulations for planets in discs with M_dics ≈ 0.01 M_⊙. For higher disc masses, the difference between the formula and our 3D simulations becomes smaller. However, the 3D simulations show a decreasing torque for increasing disc mass, while the Paardekooper et al. (2011) formula increases slightly. As the disc mass increases, viscous heating increases and cooling becomes inefficient, which results in a structure similar to an adiabatic disc. In adiabatic discs the corotation torques saturates, resulting in a lower torque acting on the planet, hence the drop of the torque for increasing disc masses. Convection certainly plays a role in more massive discs, but it is unaccounted for in Paardekooper et al. (2011). In more massive discs the convective zone reaches longer distances from the central star, disrupting the density pattern near the embedded planet and thus creating fluctuations in the torque acting on the planet. These disruptions in the density pattern are caused by the convective cells evolving in the disc. These cells also change in time, giving rise to the stronger fluctuations of the torque.

Convection is inefficient for transporting angular momentum (Kley et al. 1993; Lesur & Ogilvie 2010), but the influences of convection on planetary migration are very dramatic, because a planet close to the star in the convective zone of the disc is essentially disrupted. The direction of migration is not clearly determinable any more, but when the planet is farther out in the massive disc and the convection fades away, the direction of migration is easy to specify, indicating outward migration. Therefore the zero-torque radius for migration lies farther out in more massive discs.

Convection is also a 3D effect only and cannot be simulated in 2D discs (in r − φ direction in the midplane). Two-dimensional simulations of planets in massive discs (M_disc > 0.02 M_⊙) might therefore be inaccurate near the central star, because the effects of convection are not considered.

Acknowledgments

B. Bitsch has been sponsored through the German D-grid initiative. W. Kley acknowledges the support through the German Research Foundation (DFG) through grant KL 650/11 within the Collaborative Research Group FOR 759: The formation of Planets: The Critical First Growth Phase. The calculations were performed on systems of the Computer centre of the University of Tübingen (ZDV) and systems operated by the ZDV on behalf of bwGRiD, the grid of the Baden Württemberg state. Very helpful discussions with K.M. Dittkrist about the torque formulae are thankfully acknowledged. Finally, we gratefully acknowledge the helpful and constructive comments of an anonymous referee.

References

Appendix A: Comparison with Paardekooper et al. (2011)

In Sect. 3 we compared our numerical results to the analytical torque formula derived in Paardekooper et al. (2011). As the derivation is rather cumbersome, we present here a brief summary of the relevant contributions to the total torque acting on a planet, so that the reader can follow our calculations. In our notation we closely follow Paardekooper et al. (2011).

The total torque acting on a low-mass planet consists of two main contributions the Lindblad torque, Γ_L, plus the corotation torque, Γ_c(A.1)The Lindblad torque is caused by the action of the induced spiral arms and is given as (Paardekooper & Papaloizou 2008) (A.2)where α denotes the negative slope of the surface density profile Σ ∝ r^−α, β refers to the slope of the temperature profile T ∝ r^−β, and γ is the adiabatic index of the gas.

It is important to notice that all torques listed here are normalized to

with q the planet/star mass ratio, Σ_P the surface density at the planet’s location, Ω_P the rotation frequency of the planet in the disc, and h is in this appendix the isothermal relative disk height.

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

where the first part applies to barotropic flows where the pressure only depends on the density, and it depends on the gradient of the vorticity in the flow; the second part relates to the variations of entropy. Each of them is split again into a linear contribution and a so-called horseshoe drag contribution. This separation is necessary because the two parts are affected differently by the diffusion processes. The barotropic part of the (non-linear) horseshoe drag is given by (A.3)and the entropy-related part of the horseshoe drag is given by (A.4)where ξ = β − (γ − 1.0)α is the negative of the power-law index of the entropy. We note that the total torque formula given by Paardekooper et al. (2010), as summarized in Eq. (1), is exactly the sum Γ_tot = Γ_L + Γ_hs,baro + Γ_hs,ent.

The barotropic part of the linear corotation torque reads as (A.5)and the entropy-related part of the linear corotation torque is given by (A.6)Owing to the difference between the isothermal and adiabatic sound speed, differences in the torque arise. To compensate for this, the adiabatic index γ should be replaced by an effective γ: so that all γ’s in the previous Eqs. (A.2) to (A.6) have to be replaced by γ_eff. The parameter Q is given by

where h = H/r and

with κ beeing the opacity and σ the Stefan-Boltzmann constant. Please note that in Paardekooper et al. (2011) a factor of 4 is missing. The final correction relates to the non-ideal effects of viscosity and heat transfer, which both have to be present to avoid the saturation of the corotation torque. The barotropic part of the horseshoe drag is not affected by thermal diffusion and is only determined by the viscosity. According to Paardekooper et al. (2011) it can be written as

where Γ_hs,baro and Γ_c,lin,baro are given by Eqs. (A.3) and (A.5), but now with γ → γ_eff. F(p) (Eq. (A.7)) governs saturation and G(p) (Eq. (A.8)) and K(p) (Eq. (A.9)) govern the cut-off at high viscosity.

For the non-barotropic, entropy-related corotation torque Paardekooper et al. (2011) find where Γ_hs,ent and Γ_c,lin,ent are given by Eqs. (A.4) and (A.6), again with γ → γ_eff, and p_χ is the saturation parameter associated with thermal diffusion.

The function F(p) is given by (A.7)The function G(p) is given by (A.8)The function K(p) is given by (A.9)The parameters p_ν and p_χ, relate to the strength of viscosity and thermal diffusivity, and are given by where ν_p is the kinematic viscosity and χ_p the thermal conductivity at the planet location. x_s is the half width of the horseshoe, given by

The scaling with ϵ/h breaks down for small softening (ϵ/h < 0.3). All these contributions have to be substituted into Eq. (A.1) to calculate the total torque acting on the planet. In deriving these formulae one has to use a description for the smoothing of the gravitational potential. Here, a standard ϵ potential has been assumed, and a smoothing length of ϵ/h = 0.4 has been made.

To compare with our simulations, we have used the following parameter in the formulae above α = 0.5,β = 1.7,ξ = 1.5. To evaluate the viscosity parameter p_ν we use ν = 10^-5. Please note that ν_p ≈ χ_p for discs in radiative equilibrium.

In Kley et al. (2009) we pointed out that a different planetary smoothing results in different torques. This phenomenon can be up to a factor of two for the ϵ potential with r_sm = 0.5 and the cubic potential with r_sm = 0.5. In Sect. 3 we compared our numerical simulations to the smoothing with ϵ/h = 0.4.

Because our cubic potential with r_sm is deeper in the vicinity of the planet, we used ϵ/h = 0.3, ϵ/h = 0.35 and ϵ/h = 0.5 for the Paardekooper et al. (2011) formula as well. The results are shown in Fig. A.1. In the ϵ/h = 0.5 case, the torque is negative for all planetary distances, which is not a match at all for our simulations, but interestingly agrees better with Masset & Casoli (2010). For ϵ/h = 0.3, the torque is much higher in the positive torque regime and much lower in the negative torque regime. For the ϵ/h = 0.35 case, the formula matches our simulations quite well in the positive torque regime (r < 2.5 a_Jup). However, for larger distances to the central star, the torque of the formula is about a factor of two to four larger than the torque from our simulations. Nevertheless, the match for r < 2.5 a_Jup is quite good. Obviously, the smoothing of the planetary potential seems to have a huge effect on the torque acting on the planet. The effect seems much stronger in the formula compared to different smoothings in our 3D hydrodynamical simulations (Kley et al. 2009) but they are based on 2D simulations, where a larger impact is to be expected.

	Fig. A.1 Specific torque acting on 20 MEarth planets embedded in a 0.01 M⊙ disc. Overplotted are the results of Paardekooper et al. (2011) using various smoothings lengths from ϵ/h = 0.3 and ϵ/h = 0.5.
Open with DEXTER

All Figures

	Fig. 1 Torque acting on planets with three different masses embedded in fully radiative discs located at various distances from the star (red for 20 MEarth, blue for 25 MEarth and purple for 30 MEarth). The torques are plotted when planet and disc have reached a quasi-stationary equilibrium where the torque acting on the planet is approximately constant in time.
Open with DEXTER
In the text

	Fig. 2 Torque acting on a 20 MEarth planet in a fully radiative disc as a function of distance from the star (red plus signs). Additional curves indicate the results from the torque formula published in Paardekooper et al. (2010) (black squares) and Paardekooper et al. (2011) (purple crosses). Additionally we also plot the linear Lindblad torque by the Paardekooper et al. (2011) formula (blue asterisks).
Open with DEXTER
In the text

	Fig. 3 Torque acting on a 20 MEarth planet in fully radiative discs as a function of distance from the star. The additional curves indicate the estimates from the theoretical formula for viscous, diffusive disks by Masset & Casoli (2010) and for isothermal disks by D’Angelo & Lubow (2010).
Open with DEXTER
In the text

	Fig. 4 Radial torque density Γ(r) acting on 20 MEarth planets at rP = 0.6, rP = 1.0, rP = 2.5 and rP = 4.0 aJup. The distances have been scaled in r-direction to the actual planet location rP. The snapshots of the torque density have been taken in the equilibrium state at 80 Jupiter orbits.
Open with DEXTER
In the text

	Fig. 5 Surface density maps for a 20 MEarth planet on fixed circular orbit in fully radiative discs at four different locations. The distance from star to planet changes (from top to bottom): rP = 0.6, rP = 1.0, rP = 2.5 and rP = 4.0 (in Jupiter radii). Please note the slightly different colour scale for each plot.
Open with DEXTER
In the text

	Fig. 6 Radiative and viscous diffusion time scales that depend on the distance from the central star for our standard disc model. To compute the time scales we used the density and temperature of the midplane at the beginning of the fully radiative simulations (when the disc is in the r − θ equilibrium between viscous heating and radiative transport/cooling). Libration time scales are stated for 20 and 30 MEarth planets.
Open with DEXTER
In the text

	Fig. 7 Density (top), temperature (middle) and aspect ratio (bottom) of discs with different masses in the initial equilibrium state. All quantities are measured in the disc’s midplane at the reference distance, r = 1.0.
Open with DEXTER
In the text

	Fig. 8 Surface density (top) and temperature (bottom) of discs with different masses in the equilibrium state of the disc. The temperature is measured in the disc’s midplane.
Open with DEXTER
In the text

	Fig. 9 Torque acting on a planet located at 5 AU for different disc masses. Top: specific total torque acting on planets (20 MEarth) embedded in discs with different disc masses. The planets embedded in the higher mass discs MDisc > 0.02 M⊙ are in the convective zone in the disc, so that the torque acting on the planet is very noisy and has been averaged over 20 planetary orbits. Additionally, we over-plotted results (blue and purple) from the theoretical torque formulae of Paardekooper et al. (2010, 2011) for 20 MEarth planets. Bottom: radial torque density Γ(r) acting on the planet for different disc masses. For comparison, the vertical line indicates the location of the maximum as found for our standard case.
Open with DEXTER
In the text

	Fig. 10 Top: specific total torque evolution acting on planets (20 MEarth) embedded in discs with different disc masses. The planets embedded in the higher mass discs MDisc > 0.02 M⊙ are in the convective zone in the disc, so that the torque acting on the planet is very noisy compared to the low-mass discs. Bottom: radial torque density Γ(r) acting on the planet embedded in a disc with Mdisc = 0.04 M⊙ at different stages of the evolution.
Open with DEXTER
In the text

	Fig. 11 Surface density maps for planets on fixed circular orbits in fully radiative discs with different disc masses (from top to bottom): MDisc = 0.005 M⊙, MDisc = 0.015 M⊙ and MDisc = 0.040 M⊙. The disruption in the surface density patterns for the higher mass discs are caused by convection inside the disc.
Open with DEXTER
In the text

	Fig. 12 Time scales depending on the distance from the central star for 0.005 M⊙, 0.010 M⊙, 0.015 M⊙ and 0.040 M⊙ discs. To compute the time scales we used the density and temperature of the midplane at the beginning of the fully radiative simulations (when the disc is in the r − θ equilibrium between viscous heating and radiative transport/cooling).
Open with DEXTER
In the text

	Fig. 13 Evolution of semi-major axis for 20 MEarth planets (starting at r = 1.0 aJup) in isothermal and fully radiative discs with different disc masses.
Open with DEXTER
In the text

	Fig. 14 Vertical velocities in z-direction for planets on fixed circular orbits in fully radiative discs with different disc masses (from top to bottom) in the disc’s midplane: MDisc = 0.005 M⊙, MDisc = 0.015 M⊙ and MDisc = 0.04 M⊙. The disruption in the velocity patterns for the higher mass discs are caused by convection inside the disc. A positive velocity simulates outward flows, while a negative velocity indicates a flow towards the disc’s midplane.
Open with DEXTER
In the text

	Fig. 15 Velocities in z-direction for 0.04 M⊙ discs without planet. In the top panel only the upper half of the disc was simulated, while in the bottom panel both sides of the disc were simulated (with twice the number of grid cells in θ direction). The other simulation parameters are identical.
Open with DEXTER
In the text

	Fig. 16 Velocities in z-direction for 0.04 M⊙ discs with embedded 20 MEarth planet. In the top panel only the upper half of the disc was simulated, while in the bottom panel both sides of the disc were simulated (with twice the number of grid cells in θ direction). The other simulation parameters are identical.
Open with DEXTER
In the text

	Fig. 17 Specific torque acting on 20 MEarth planets embedded in 0.02, 0.025, 0.03 and 0.04 M⊙ discs at the following distances rp = 2.0, 2.5, 3.0 and 4.0aJup, respectively.
Open with DEXTER
In the text

	Fig. A.1 Specific torque acting on 20 MEarth planets embedded in a 0.01 M⊙ disc. Overplotted are the results of Paardekooper et al. (2011) using various smoothings lengths from ϵ/h = 0.3 and ϵ/h = 0.5.
Open with DEXTER
In the text