© ESO, 2013

1. Introduction

Migration of a low-mass planet embedded in a fully radiative gaseous disk can be significantly different from migration in an isothermal or purely adiabatic disk (Paardekooper & Mellema 2006, 2008; Baruteau & Masset 2008; Paardekooper & Papaloizou 2008; Kley & Crida 2008; Kley et al. 2009; Ayliffe & Bate 2010). While all authors agree that radiation transport can slow the rate of inward migration, there is still a lack of consensus whether the direction of migration can be outward. Part of this confusion may be due to the sensitivity of the direction and magnitude of migration on global disk parameters (Paardekooper et al. 2010, 2011, Masset & Casoli 2010), including, e.g., the radial disk temperature gradient (Ayliffe & Bate 2011). Different authors also use different viscosities, usually either a constant viscosity or a Shakura & Sunyaev (1973)α-viscosity, with a typical value of α ranging between 10^-4 and 10^-2. The simulation of an embedded planet in a global magnetorotationally unstable accretion disk through direct magnetohydrodynamic calculations including radiative transport is still computationally too expensive. However, the disk structure of magnetorotationally unstable disks including radiative transport has been calculated in the local shearing box approximation (Flaig et al. 2012).

An unperturbed, viscous, fully radiative disk will evolve toward an equilibrium state, where viscous heating is balanced by radiation transport and cooling (as described in, e.g., Kley et al. 2009). This equilibrium state is dependent on the disk mass, the viscosity, and the adiabatic index of the gas. Variations in viscosity change the radial density and temperature profiles, thus changing the profiles of vortensity and entropy. The vortensity and entropy gradients across the horseshoe region, in turn, determine the magnitude and sign of the corotation torque (Baruteau & Masset 2008). If the torques can remain unsaturated due to the action of viscosity, the total torque can be positive, leading to outward migration. This effect is also possible in isothermal disks (Paardekooper & Papaloizou 2008).

There is an additional complication to understanding the migration of low-mass planets in disks. The rotational states of molecular hydrogen are only fully accessible at temperatures  ≳300 K (Decampli et al. 1978; Boley et al. 2007). As a result, the adiabatic index of the gas will transition from γ = 7/5 to 5/3 as the temperature in a disk drops with radial distance. The region over which this variation occurs is precisely where we expect planetary cores and planets to form in the core accretion scenario and begin their initial stages of migration. So far, only a constant adiabatic index of γ = 1.4 has been explored in our previous simulations with radiative transport.

All of these effects on low-mass planet migration are relevant to understanding whether planet traps can exist, i.e., regions in the disk where a protoplanet could experience zero torque. Protoplanets migrating from a smaller radius outwards or from an outer radius inwards would collect at the zero-torque radius, creating areas conducive to planetary mergers, possibly leading to the growth of large cores. A planet trap could also be formed by surface density changes, which could create an enhanced feeding zone for these cores (Morbidelli et al. 2008); however, it is yet unclear how realistic surface density changes are in disks. In contrast, radiation transport might allow traps to exist in even smooth disk structure. As outward migration is dependent on the viscosity, disk mass, and adiabatic index, these parameters will influence the radius and breadth of the zero-torque region. In this work, we present a study on how these effects modify the migration properties of embedded planets.

The paper is organized as follows. In Sect. 2 we give an overview of our numerical methods. We then describe the influence of the adiabatic index, viscosity, and the differences between fixed ortho-to-para and equilibrium gas mixtures on the disk structure in Sect. 3. These structural changes influence the migration rate of an embedded 20 M_Earth planet, which is discussed in Sect. 4. In Sect. 5 the influence of viscosity and of the different gas mixtures on the zero-torque radius are investigated. We then summarize and conclude in Sect. 6.

2. Numerics and setup

2.1. Setup

The protoplanetary disk is modeled as a three-dimensional (3D), non-self-gravitating gas. Fluid motion is described by the Navier-Stokes equations, where the equations are solved numerically using a spatially second-order finite volume method that is based on the code NIRVANA (Ziegler & Yorke 1997). The disk is heated solely by internal viscous dissipation, and is allowed to cool by flux-limited diffusion (FLD, Levermore & Pomraning 1981). The FLD approximation allows internally produced energy to diffuse radiatively through the optically thick regions of the disk and into the optically thin regions, where the energy can be radiated away by free-streaming. The flux-limiter interpolates between the optically thick and thin regimes and ensures that energy loss never exceeds the free-streaming limit. In the code, radiative transport is handled implicitly and it uses the FARGO (Masset 2000) extension as described in Kley et al. (2009). A more detailed description of the modeling and the numerical methodology is provided in our previous papers (Kley et al. 2009; Bitsch & Kley 2010, 2011a).

The three-dimensional (r,θ,φ) computational domain (with 266 × 32 × 768 active cells) consists of a complete annulus of the protoplanetary disk centered on the star, extending from r_min = 0.4 to r_max = 2.5 in units of r₀ = a_Jup = 5.2 AU. For simulations with planets at larger distances from the central star, the computational grid is set to have an inner computational boundary at 0.4r_P and an outer boundary at 2.5r_P. In the vertical direction the annulus extends 7° above the disk’s midplane, meaning 83° < θ < 90°. As we assume that the disk is symmetric for the upper and lower parts, we apply symmetric boundary conditions at θ = 90°. Here θ denotes the polar angle of the used spherical polar coordinate system measured from the polar axis. The central star has one solar mass M_∗ = M_⊙, and the total disk mass inside [r_min,r_max] is M_disk = 0.01 M_⊙. The aspect ratio of the disk is calculated self-consistently from the equilibrium structure, given by viscous internal heating and radiative diffusion. This also determines the surface density gradient in the equilibrium state of the disk. To construct the equilibrium state, we first calculate axisymmetric 2D models in the r − θ directions.

The planet is then embedded at r_P = 1.0 a_Jup. For the planet’s gravity, we use the cubic potential (Klahr & Kley 2006; Kley et al. 2009) with r_sm = 0.5. The planetary potential and its influence on migration is discussed in great detail in Kley et al. (2009). As in our previous simulations we use the opacity law provided by Lin & Papaloizou (1985).

2.2. 2D axisymmetric models

The initialization through an axisymmetric 2D phase (in the r − θ plane) reduces the required computational effort substantially. The evolution from the initial isothermal state toward the disk equilibrium between viscous heating and radiative transport/cooling (henceforth the disk equilibrium state) takes about 100 orbits, if the disk is started in equilibrium for an isothermal gas and constant H/r. The surface density or temperature profiles of the initial (isothermal) state are unimportant as the equilibrium state of the disk solely depends on the disk mass, viscosity, and adiabatic index of the disk.

After reaching the disk equilibrium state, we extend this model to a full 3D simulation by expanding the grid into the φ-direction and by embedding the planet. From this starting configuration it is also possible to investigate the vertical structure of an unperturbed disk in the equilibrium state. This corresponds to a disk with zero total mass flow because the radial boundaries do not allow in- or outflow.

2.3. Adiabatic index

The gas in an accretion disk is primarily molecular hydrogen, which exists as para and ortho-hydrogen for proton spins that are antiparallel and parallel, respectively. If enough protonated ions (e.g., ${\mathit{H}}_{\mathrm{3}}^{\mathrm{+}}$) are available to exchange proton spins on timescales that are shorter than the dynamical time, the para- and ortho-hydrogen should be treated as being in statistical equilibrium (Boley et al. 2007). If the dynamical timescale is shorter than the equilibrium timescale, a fixed ratio should be used, where a 3:1 ortho-to-para mix is common for many astrophysical systems. A different mixture of gas changes the ratio of specific heats, which is the adiabatic index γ. The adiabatic index is dependent of the temperature of the underlying gas (Boley et al. 2007), and astrophysical disks should transition between γ = 5/3 to γ = 7/5 near the frost line. This is represented in Fig. 1, where the adiabatic index’s dependence on temperature is plotted for the equilibrium (blue) and the 3:1 ortho-to-para mix (red) gas. The adiabatic index can span quite a large region in γ for low temperatures. To make things clearer for the reader, we will refer to the equilibrium gas mixture as O:P(equi) and as O:P(3:1) for the 3:1 ortho-to-para mix. This is also indicated in Fig. 1. In our typical disk simulations (with γ = 1.4) the temperature at r = 1.0 a_Jup is about 65 K, which is in a region where the adiabatic index changes rapidly, for the O:P(equi) gas and is close to the temperature where the rotational states first activate for the O:P(3:1) gas. Such variation could affect an embedded planet’s migration rate, as the planet modifies the temperature structure of its surroundings. For the disks and orbital times explored in this study, a fixed ortho-to-para ratio is likely the most physical treatment for H₂, with the ortho-para equilibration timescale becoming comparable to the orbital timescale outside of 40 AU. Nevertheless, we explore both the equilibrium and the 3:1 fixed ratio cases to show limiting behaviors.

Fig. 1 Adiabatic index for the 3:1 ortho to para mix [O:P(3:1)] and equilibrium gas [O:P(equi)] configuration in dependence of the temperature.

2.4. Viscosity

The viscosity in an accretion disk is typically assumed to be either constant or an α-type viscosity (Shakura & Sunyaev 1973), for which we use here $\mathit{\nu }\mathrm{=}\mathit{\alpha }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathrm{\Omega }$, where c_s denotes the local adiabatic sound speed and Ω the local orbital angular velocity. If the accretion disk is in a steady state, the viscosity can be determined by the observed mass accretion rate Ṁ onto the star: $\mathit{Ṁ}\mathrm{\approx }\mathrm{3}\mathit{\pi \nu }\mathrm{\Sigma }\mathrm{=}\mathrm{3}\mathit{\pi \alpha }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathrm{\Omega }\hspace{0.17em}\mathrm{\Sigma }\mathit{,}$(1)where Ṁ is typically 10^-9−10^-7 M_⊙/yr (Hartigan et al. 1995; Gullbring et al. 1998) and Σ is the surface density. The viscosity determines the structure of the disk and thus influences the migration of embedded planets. Based on models of accretion disks, the observed mass accretion rates constrain α to be roughly between 10^-4 and 10^-2. As will be discussed later in the manuscript, we define the half-disk thickness as $\mathit{H}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{iso}}}{\mathrm{\Omega }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Omega }}\hspace{0.17em}\sqrt{\frac{\mathit{P}}{\mathit{\rho }}}\mathrm{·}$(2)where c_iso is the isothermal sound speed, and P and ρ denote the pressure and density in the disk’s midplane.

2.5. Torque calculations

In previous work, we have discussed the calculation of the torque acting on the planet in great detail. Outward migration seems only to be possible when the planet is on a nearly circular orbit (Bitsch & Kley 2010) and close to the midplane of the disk (Bitsch & Kley 2011a). Since we want to investigate the occurrence of outward migration and its implications for planet formation and evolution, we assume that the planets remain in the midplane of the disk and are initialized on circular orbits. Only the upper half of the disk (83° ≤ θ ≤ 90°) is needed for the calculation, as the lower half is directly symmetric to the upper half for these conditions.

The torques acting on a 20 M_Earth planet are calculated to determine the direction of migration. A positive torque indicates outward migration, while a negative torque indicates inward migration. As the planet is simulated as a point mass, the planetary potential needs to be smoothed. We use the cubic potential (Klahr & Kley 2006; Kley et al. 2009) for our calculations: ${\mathrm{\Phi }}_{\mathrm{p}}^{\mathrm{cub}}\mathrm{=}\{\begin{array}{c}\\ \mathrm{-}\frac{{\mathit{m}}_{\mathrm{p}}\mathit{G}}{\mathit{d}}\hspace{0.17em}\left[{\mathrm{(}{\frac{\mathit{d}}{{\mathit{r}}_{\mathrm{sm}}}}^{\mathrm{)}}}^{\mathrm{4}}\mathrm{-}\mathrm{2}{\mathrm{(}{\frac{\mathit{d}}{{\mathit{r}}_{\mathrm{sm}}}}^{\mathrm{)}}}^{\mathrm{3}}\mathrm{+}\mathrm{2}\frac{\mathit{d}}{{\mathit{r}}_{\mathrm{sm}}}\right] & \mathrm{for} \mathit{d}\mathrm{\le }{\mathit{r}}_{\mathrm{sm}}\\ \mathrm{-}\frac{{\mathit{m}}_{\mathrm{p}}\mathit{G}}{\mathit{d}} & \mathrm{for} \mathit{d}\mathit{>}{\mathit{r}}_{\mathrm{sm}}\mathit{.}\end{array}$(3)Here m_P is the planetary mass, d = |r − r_P| denotes the distance of the disk element to the planet, and r_sm is the smoothing length of the potential, usually a fraction of the Hill radius, R_H = (1/3 m_p/M_∗)^1/3. The construction of the planetary potential is in such a way that, for distances larger than r_sm, the potential matches the correct 1/r potential. For d < r_sm, the potential is smoothed by a cubic polynomial. The parameter r_sm is equal to 0.5 R_H in all our simulations.

The gravitational torques acting on the planet are calculated by summing over all grid cells, which are treated in this case as point masses. We also 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 large, probably noisy contributions from the inner parts of the Roche lobe and to disregard material that is possibly gravitationally bound to the planet (Crida et al. 2009). Here we assume (as in our previous papers) a transition radius of 0.8 R_H.

3. The equilibrium disk structure

As the adiabatic index is connected to the sound speed, a change in the adiabatic index can alter the structure of the disk by modifying the disk viscosity, affecting the disk equilibrium state. The effects of these parameters on the disk structure (without an embedded planet) are discussed in this section.

Fig. 2 Densities (top), temperatures (middle), and aspect ratios (bottom) from the 2D axisymmetric disks with different (constant) adiabatic indices, taken after the disks have reached their equilibrium state. Each quantity is evaluated using midplane values of the disk at a radius r = 1 aJup. The kinematic viscosity has been set to ν = 1015 cm2 s-1 for all models.

3.1. Influence of a constant adiabatic index

For the following simulations, we limit ourselves to constant values of γ in the range of 1.05 ≤ γ ≤ 1.8, which is larger than the range of the adiabatic index in Fig. 1. We include γ values smaller than 1.3 to capture the transition to isothermal disks, and higher γ values to capture additional possible stiffness resulting from magnetic fields. In a first set of simulations we study disks with a constant kinematic viscosity of ν = 10¹⁵ cm² s^-1. For this case, we do not expect a change in the disk structure for different values of γ because the disk equilibrium state is determined by a balance between viscous heating and radiative cooling, both of which do not depend of γ. Indeed, as can be seen in Fig. 2, for a wide range of γ, the midplane density and temperature and the resulting aspect ratio at r = 1.0 a_Jup for a 2D axisymmetric disk do not depend on γ. Only for small values of γ is the disk structure different. The reason for this change is the onset of convection for γ ≲ 1.2. As the adiabatic index is lowered, the gas becomes more compressible owing to the larger degrees of freedom that must be excited before the thermodynamic temperature can be raised. This leads to a higher midplane density, a lower midplane temperature, and a lower aspect ratio (see Eq. (2)) in the disk at r = 1.0 a_Jup.

Radial density and temperature gradients result in vortensity and entropy ($\mathit{S}\mathrm{\propto }\frac{\mathit{p}}{{\mathit{\rho }}^{\mathit{\gamma }}}$) gradients. These vortensity and entropy gradients across the horseshoe region determine the magnitude and sign of the corotation torque, which, if positive, can lead to a total positive torque, indicating outward migration in fully radiative disks. As the density decreases and the temperature increases with increasing γ (for γ < 1.2), the strength of the gradient in entropy might change and therefore influence planetary migration. The influence of this effect might be stronger for the O:P(equi) gas disk, as the adiabatic index can become lower compared to the O:P(3:1) gas disk.

The radial variation of the aspect ratio of a fully radiative disk with different adiabatic indices is shown in Fig. 3.

For γ = 1.1, the aspect ratio is smaller at r < 1.5 a_Jup and the fluctuations in this region are a sign of convection. However, such a small adiabatic index is not predicted for accretion disks, as seen by the γ(T) profile in Fig. 1, unless efficient heating and cooling push the disk toward isothermal behavior or the gas is undergoing dissociation or ionization. In contrast, the profiles for a higher adiabatic index show an overall drop in aspect ratio with increasing distance to the central star (after a small increase for smaller distances to the central star, r < 0.7). Such a feature in the disk’s structure is the opposite of what one would expect. When taking stellar irradiation into account an increasing H/r profile for larger radii may be possible, e.g. Dullemond & Dominik (2004). However, the disk’s opacity can greatly influence the H/R profile.

With increasing distance from the central star, the temperature is decreasing. For temperatures less than T ≈ 155 K, the opacity increases with increasing temperature. At higher temperature, the opacity decreases. So, with increasing distance to the central star, the opacity first increases and, when T ≈ 155 K, the opacity decreases. This is exactly the point in the disk (r ≈ 0.7) where the aspect ratio H/r begins to decrease. Please note that the height H of the disk is still increasing, but at a lower rate, when T < 155 K compared to T > 155 K. Hence, the structure of the disk is a result of the opacity law that is used in these calculations.

The H/r profiles for the disks with γ > 1.4 are identical, as the disk equilibrium state is determined by a balance of viscous heating and radiative cooling, neither of which depends on γ, at least for the constant ν assumed for these disks.

Fig. 3 Aspect ratio (H/r) of a fully radiative disk with different adiabatic indices without an embedded planet. The aspect ratio was computed using midplane values of the constant viscosity disks (see Eq. (2)).

3.2. Influence of viscosity

As can be seen in Fig. 4, where we display the midplane density (top) and temperature (bottom) of fully radiative disks in the disk equilibrium state, a change in the viscosity affects the disk’s equilibrium structure, altering the density and temperature profiles of the disk. In this figure we show the densities and temperatures for α-viscosities 0.001 ≤ α ≤ 0.008 and for a constant viscosity of ν = 10¹⁵ cm² s^-1, a value used in our previous simulations (Kley et al. 2009; Bitsch & Kley 2010, 2011a).

Fig. 4 Density (top) and temperature (bottom) profiles of the fully radiative disks with different viscosities, all taken in the disk equilibrium state. Densities and temperatures are measured in the midplanes of the disks.

For increasing α-viscosity, the density in the midplane of the disk becomes smaller, but the differences in density become smaller for higher viscosities.

For the temperature profile, the trend is reversed: a higher viscosity results in a higher temperature in the midplane of the disk. As the viscosity increases, so does the viscous heating. The heating process mainly takes place in the midplane of the disk, as the viscosity is highest there. Not only is the temperature higher in the midplane, but the temperature profiles also change with viscosity. For our constant viscosity of ν = 10¹⁵ cm² s^-1, the radial slope of temperature was  ∝ r^-1.8 at r = 1 a_Jup, while it now depends on the chosen α-value. The slope of temperature β_T, where T ∝ r^−β_T, is displayed in Table 1 for r = 1 a_Jup.

The temperature slope has a direct influence on the migration of embedded planets, e.g., Paardekooper et al. (2010, 2011), as the Lindblad torque and the entropy related corotation torque are very sensitive to it. A change in β_T therefore influences the torque acting on an embedded planet quite severely (Ayliffe & Bate 2011).

The aspect ratio H/r of the disk also increases for increasing viscosity (see Table 1), indicating another factor that influences planetary migration in these disks. For isothermal disks, this has been know for a while (Tanaka et al. 2002).

The density, temperature and H/r profiles of our constant viscosity simulation match quite well with the α = 0.005 simulations, indicating that one could expect a very similar torque acting on a planet embedded in these disks. However, as the α-viscosity is not constant in r, one might also expect some differences in the spiral wave densities and temperatures in these disks. In addition, the area near the planet might be subject to some changes for different viscosities.

3.3. Influence of a varying adiabatic index

In Fig. 5 the temperature (bottom) and density (top) profiles in the midplane of fully radiative disks are displayed. The disks feature different α viscosities and two different ortho-para mixtures of the gas [O:P(equi) and O:P(3:1)]. As expected from simulations with a constant adiabatic index, an increasing α viscosity results in lower densities and higher temperatures in the midplane of the disk.

The gradients in density and temperature are comparable to the gradients for constant adiabatic indices. It also seems that the mixture of the gas does not have a huge influence on the temperature and density profiles for large viscosities. However, for α = 0.001 there are differences in the density and temperature profiles between disks that have different ortho-para ratios. For r < 1.0 a_Jup, the temperature for the disk equilibrium state is lower in the O:P(equi) than in the O:P(3:1). As a result, the density is higher in the O:P(equi) disk. As γ is largely reduced for temperatures around T ≈ 50 K for ortho-para ratios O:P(equi), this behavior is consistent with the results of the constant γ simulations.

Fig. 5 Midplane density (top) and temperature (bottom) profiles for the α viscosity disks with different assumptions for the ortho-para ratio of the gas.

In Fig. 6 the aspect ratios H/r for the same disks are shown. The profiles are very similar to those in disks with a constant adiabatic index and an α viscosity. In the α = 0.001 case, the aspect ratio profile is smaller and approximately constant. For the more viscous disks, the aspect ratio increases for small r, but then decreases with increasing r continuously. A higher viscosity leads to a higher aspect ratio in the disk. The O:P(equi) shows a slightly smaller aspect ratio for viscosities α ≥ 0.004 compared to the O:P(3:1). The aspect ratio differs by about 10% for the high viscosity cases at r = 1.0 a_Jup. This difference at this location can influence the torque acting on embedded planets, as seen in isothermal disks.

Fig. 6 Aspect ratio H/r of the accretion disks for different α viscosities and different assumptions for the ortho-para ratio of the gas.

Fig. 7 Velocities in z-direction for disks with α = 0.001 and O:P(3:1) (top), α = 0.001 and O:P(equi) (2nd from top), α = 0.004 and O:P(equi) (3rd from top), and a fixed γ = 1.1 with constant viscosity (bottom). The black line indicates H(r) in the disk and the red line delineates τ = 1.0 (the optical depth integrated from the top of the grid toward the midplane of the disk). Please note the different color scales and that they only show a subset of the range of values to better highlight the internal structure.

A higher disk viscosity leads to a larger temperature in the midplane. Recently, Bitsch & Kley (2011b) have shown that in fully radiative disks with viscous heating, a convective region arises near the central star, where its radial extent depends on disk mass. A higher disk mass resulted in a larger convective zone in the disk. The density and temperature structures of the disk with constant viscosity seem to correspond to a disk with α ≈ 0.005 (see Fig. 4).

We display in Fig. 7 the velocity in the z-direction for the α = 0.001 disks with two different ortho-para gas mixtures, for the α = 0.004 and O:P(equi) disk, and for the fixed γ = 1.1 disk. In black we indicate H(r) and in red the τ = 1.0 line in the disk (integrated from the top of the grid toward the midplane). We note that when taking both sides of the disk into account, convective eddies cross the midplane. However, Bitsch & Kley (2011b) have shown that the first signs of convection can also be seen in simulations where only one half of the disk is simulated. The extent to which the convective region inside the disk depends on the adiabatic index should be tested in simulations of disk instabilities. As we want to focus here on the influence of the adiabatic index on planet migration, we do not pursue this topic further. We do note that Ruden & Pollack (1991) have shown that a lower adiabatic index can increase the susceptibility of the disk to convection.

For the α = 0.001 disk with either the O:P(3:1) or the O:P(equi) gas, variations in the vertical velocity of the gas in the disk’s upper layers can be seen in Fig. 7. However, these fluctuations are limited to the low-density material, and do not reach all the way down to the high optical depth region. Convection does not appear to be present. Fluctuations such as those in the top two panels of Fig. 7 only occur for the thinner disks (small H/r) that have a very large vertical range in density, which is difficult to resolve numerically (in this case no-density floor is used, possibly increasing the effect). Nelson et al. (2012) stated that vertical instabilities in isothermal disks can arise, if the viscosity is negligible. This vertical shear instability can then drive a non negligible transport of angular momentum through Reynolds stress with a the corresponding viscous alpha value α = 0.001, which is the value of viscosity assumed in the simulations shown in the top two plots in Fig. 7, indicating that the shown fluctuations might be resulting from the vertical shear instability. However, detailed studies of the vertical shear instability for the case of fully radiative disks is beyond the scope of the presented work.

Convection, on the other hand, can only be present in the optically thick regions of the disk, where radiative cooling is inefficient. For larger viscosities (α = 0.004), the fluctuations in the upper layers mostly vanish, except in the low-density region at large r, and are totally diminished for α = 0.008 (not shown). Large viscosities suppress the growth of vertical oscillations, in agreement with the vertical shear instability (Nelson et al. 2012).

On the other hand, the fixed γ = 1.1 disk shows fluctuations reaching all the way down to midplane, where the disk is optically thick. Please also note here that the absolute values of the velocities are larger for the γ = 1.1 disk compared with the other presented simulations. In this case, we have to check if convection is really present in the disk, as we already suspected from the aspect ratio profile of the disk (see Fig. 3).

In order to test whether convection should be present in the disk, we compute the adiabatic gradient as follows: $\mathit{\beta }\mathrm{=}\frac{\mathrm{d}\ln \mathit{T}}{\mathrm{d}\ln \mathit{p}}\mathrm{·}$(4)The adiabatic gradient is evaluated at a two-degree angle from the midplane, which roughly corresponds to the half opening angle of the optical thick parts of the disk.

For β > β_adi, convection can be present in the disk, where β_adi is defined as ${\mathit{\beta }}_{\mathrm{adi}}\mathrm{=}\frac{{\mathit{\gamma }}_{\mathrm{local}}\mathrm{-}\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{local}}}\mathrm{·}$(5)Here, γ_local represents the local adiabatic index, which in our case, is the local adiabatic index along the two-degree angle of the disk. In Fig. 8 we present the the adiabatic gradients of the γ = 1.1 disk (bottom in Fig. 7) and the α = 0.004 disk with O:P(equi) (3rd from top in Fig. 7).

Fig. 8 Adiabatic gradients for the γ = 1.1 disk with constant viscosity (top) and for the α = 0.004 disk with O:P(equi) (bottom). The gradients have been evaluated in the optical thick regions of the disk at an angle of 2° from the midplane.

In the case of γ = 1.1, β > β_adi, indicating that convection should be present in the disk, which corresponds to the vertical fluctuations for this simulation in the optically thick regions of the disk. The strong fluctuations in β indicate that the disk is highly unstable. Please also note that in the inner regions of the disk, the difference between β and β_adi is about a factor of 5. For the α = 0.004 disk, only the inner regions feature β > β_adi. However, fluctuations in the velocity have not been observed (Fig. 7). In addition, the difference between β and β_adi is at maximum a factor of 2, which is much lower than in the γ = 1.1 case. If convection is present in the α = 0.004 disk, then it must be much less vigorous than seen in the γ = 1.1 disk.

4. Influence on planetary migration

Because the above results show that the adiabatic index and the viscosity affect disk profiles, we also expect that they will change the migration rates of embedded planets. Embedding a small mass planet in such a disk will, of course, change the profiles of the disk. However, for most viscosities and adiabatic indices, the changes that the low-mass planet will cause are limited to the vicinity of the planet, so that the overall structure of the disk remains intact. To determine the direction of migration of protoplanets in disks, we embed 20 M_Earth planets in our disks and measure the torque acting on those planets. A positive torque represents outward migration, while a negative torque indicates inward migration. The torque acting on planets on circular, non-inclined orbits is directly proportional to the migration rate. In cases where the planet is embedded in the convective region in the disk (for γ ≲ 1.1), the torques have been averaged over 20 planetary orbits.

4.1. Constant adiabatic index

In Fig. 9 we show the total torque Γ_tot acting on planets that are on circular orbits and are embedded in disks that have different adiabatic indices. The viscosity is held constant at ν = 10¹⁵ cm² s^-1. The torque acting on the planet is positive when the adiabatic index is larger than 1.12 , while it is negative, indicating inward migration, for γ < 1.12. The maximum of the torque is at γ = 1.2, and the torque decreases for lower and higher γ’s. The most negative torque is at γ = 1.075.

For γ ≥ 1.12, changing the adiabatic index influences the overall results of planetary migration only by magnitude, not by direction. At high adiabatic indices, the torque acting on the embedded planet is decreased by a factor of 3 to 4 compared to the maximum torque at γ = 1.2. The tendency of outward migration, however, remains intact. An increase in γ leads to an increase in the adiabatic sound speed, which in turn decreases the amplitude of the total torque. It also reduces the amplitude of the (negative) entropy gradient, which then reduces the amplitude of the (positive) corotation torque. This can be seen by the decrease of the (positive) spike in the torque density distribution (see Fig. 10). For lower adiabatic indices, the direction of migration is reversed as γ approaches the value for isothermal disks. However, the torque acting on planet in a radiative disk with γ = 1.001 does not reflect the torque acting on a planet in an isothermal disk, as the disk structure is different (Paardekooper & Mellema 2008; Kley et al. 2012). On the other hand, the torque has strong fluctuations in time, indicating that the planet has reached the convective zone of the disk.

Fig. 9 Total torque Γtot acting on a 20 MEarth planet embedded at r = 1 in ν = 1015 cm2 s-1 disks for different (constant) adiabatic indices. The torques have been measured in the equilibrium state of the disk.

In Fig. 10, the radial torque density Γ(r), which is the torque exerted by an annulus of width dr and related to the total torque by ${\mathrm{\Gamma }}_{\mathrm{tot}}\mathrm{=}{\mathrm{\int }}_{{\mathit{r}}_{\min }}^{{\mathit{r}}_{\max }}\hspace{0.17em}\mathrm{\Gamma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r,}$(6)is displayed for 20 M_Earth planets embedded in fully radiative disks with a constant adiabatic index. All displayed torque densities Γ(r) (with γ > 1.2) feature a spike just below r ≈ 1.0 a_Jup on top of the underlying Lindblad torque curve, which is the normal pattern for a fully radiative disk with a high adiabatic index (Kley et al. 2009). For increasing γ, the Lindblad torque and the spike in the torque density are reduced, indicating a reduction of the positive corotation torque. This reduction in torque density is also represented by a decrease in the total torque (see Fig. 9). For γ = 1.1, no spike in the torque density is visible, and only the Lindblad torque remains. This also is reflected in the total torque, which is negative (see Fig. 9).

Fig. 10 Profiles of the radial torque density Γ(r) acting on 20 MEarth planets embedded in disks with different adiabatic indices (during the disk equilibrium state). In contrast to the total torque, the radial torque density has not been averaged in time.

Figure 11 shows surface density images of our fully radiative disk simulations with an embedded 20 M_Earth. The adiabatic indices for these images are γ = 1.1 (top), γ = 1.4 (middle) and γ = 1.7 (bottom). We note that the γ = 1.4 case is discussed in great detail in Kley et al. (2009). For the γ = 1.1 simulation, strong variations in the density surrounding the planet are visible. The variations are due to convection inside the disk.

The surface density pattern of the γ = 1.7 disk, on the other hand, shows no sign of convection. In fact it is very similar to the γ = 1.4 (middle in Fig. 11). The only difference is that the density perturbations in front and behind the planet are not as well pronounced. Due to a difference in the adiabatic sound speed, which is  ∝ $\sqrt{\mathit{\gamma }}$, the opening angles of the spiral waves increase with increasing γ. The density in the spiral waves is reduced for increasing γ, resulting in a decrease of the Lindblad torque, as can be seen in the torque density distribution (Fig. 10), and in the total torque (Fig. 9).

Fig. 11 Surface density maps for planets on fixed circular orbits in fully radiative disks with different adiabatic indices (top: γ = 1.1, middle: γ = 1.4, bottom: γ = 1.7) and a constant viscosity of ν = 1015 cm2 s-1. The disruption in the surface density pattern for the γ = 1.1 disk is due to convection inside the disk.

4.2. Influence of viscosity

In Sect. 3.2, we discussed the changes of the structure of planetary disks with different viscosities. Changes in the temperature gradient lead to changes in the torque acting on embedded planets (Paardekooper et al. 2011). This sensitivity is highlighted in Fig. 12, which shows the total torque acting on an embedded 20 M_Earth planet in disks with different viscosities. Planets in disks with α ≤ 0.002 (not displayed) feel a negative torque, indicating inward migration. A decrease of the torque with decreasing viscosity can also be found in isothermal disks (see for example our simulations in Fig. 1 in Bitsch & Kley 2010). The decrease in the total torque is a result of the lower viscosity. If the viscosity is very low, the corotation torques saturate and outward migration can not be supported any longer (for isothermal disks, see Masset 2001). The planet depletes its coorbital region of gas (i.e., begins to open a gap, Fig. 14), which reduces the torque.

For α viscosities with α ≥ 0.003, the torque acting on the planet becomes positive, and for even larger viscosities (α ≥ 0.004), the magnitude of the torque begins to converge. The torque acting on the planet with constant viscosity of ν = 10¹⁵ cm² s^-1 is similar to the one with α ≈ 0.006. This is consistent with the similarities between their temperature and density profiles (Fig. 4).

Fig. 12 Torque acting on 20 MEarth planets embedded in disks with different viscosities (γ = 1.4).

In Fig. 13 we display the radial torque density Γ(r) acting on a 20 M_Earth planet in disks with different viscosity. In the α = 0.001 simulation, the displayed torque density shows a typical Lindblad torque curve, but without the spike at corotation (see below); a spike in the torque curve at corotation is the typical pattern of a fully radiative disk with a high viscosity (Kley et al. 2009). In fact the torque density is very similar to an isothermal disk.

The α = 0.005 and α = 0.006 torque density patterns are consistent with the that of the constant viscosity of ν = 10¹⁵ cm² s^-1. The spike in the distribution near r ≈ 1.0 a_Jup is nearly identical for the constant viscosity and α = 0.006 simulations, while the α = 0.005 spike is a little bit higher. However, there are some differences. For smaller and larger distances from the central star (r < 1.0 a_Jup and r > 1.0 a_Jup) the torque density is larger for the constant viscosity case compared to α = 0.006, although all simulations follow the same trend. The underlying Lindblad torque distribution is also smaller for the α = 0.005 disk, due to a difference in the overall disk structure (see Fig. 4). The reduction of Γ(r) for larger α is partly a consequence of the increased disk temperature. These small differences are the reason why the total torque of the constant viscosity simulations is somewhat larger compared to the α viscosity disks.

When α = 0.008, the Lindblad torque is much smaller compared to the other disks, as the temperature is higher. However, the torque density shows the usual spike, indicating a positive corotation torque, which can lead to outward migration if it over-compensates for the negative Lindblad torque (Kley et al. 2009). One might suspect from the trend of the simulations that even higher viscosities could destroy this effect of outward migration. Moreover, modifications to the horseshoe dynamics for large viscosities could limit or cutoff the torque acting on the planet (Masset 2001, 2002).

Fig. 13 Profiles of the radial torque density Γ(r) acting on 20 MEarth planets embedded in disks with different viscosities (in the disk equilibrium state) and a constant adiabatic index (γ = 1.4).

In Fig. 14 we display the surface densities of several disks with an embedded 20 M_Earth planet. Different viscosities are shown to gain more insight into the origin of the torque density profiles. For very low viscosities (α = 0.001), the planet seems to open a very small or partial gap inside the disk. In addition, the spiral waves of the planet are more dense compared to the constant viscosity simulation (middle panel in Fig. 11), which may be due to the higher density in the midplane of the disk. A low viscosity encourages the formation of gaps in the disk, thus explaining the torque density in Fig. 13.

Fig. 14 Surface density maps for 20 MEarth planets on fixed circular orbits in fully radiative disks (with γ = 1.4). The viscosities in the disks are (from top to bottom): α = 0.001, α = 0.005, and α = 0.008. The snapshots are taken in the equilibrium state.

The α = 0.005 surface density profile seems very similar to the constant viscosity profile we discussed in great detail in Kley et al. (2009). The density increase ahead of the planet (φ > 180° and r < 1.0 a_Jup) is clearly visible in both cases, thus creating a nearly identical “spike” in the torque density distribution (Fig. 13). However, the density decrease behind the planet (φ < 180° and r > 1.0 a_Jup) is not so clear in the α = 0.005 simulation compared to the ν = 10¹⁵ cm² s^-1 simulation (middle panel in Fig. 11). This directly reflects on the torque density, as it is higher for the constant viscosity simulation compared to the α = 0.005 simulation at that distance from the central star, thus explaining the higher total torque.

For α = 0.008, the surface density profile shows the same structure as for the constant viscosity simulation, but the surface density is generally reduced. The spiral waves and the vicinity near the planet show smaller surface densities. The structure ahead and behind the planet is also not as distinctive as in the constant viscosity (or α = 0.005) simulation. This all leads to a smaller curve in the torque density plot (Fig. 13) and to a smaller total torque (Fig. 12).

4.3. Varying adiabatic index

In addition to different viscosities, we now consider models that feature a varying adiabatic index. In Fig. 15 the total torque acting on 20 M_Earth planets in disks is displayed for different viscosity and gas treatments. For low viscosities (α = 0.001) the torque is negative (not displayed), indicating inward migration, and for higher viscosities, the torque is positive indicating outward migration, as expected from simulations with a constant adiabatic index.

For all shown viscosities, the torque acting on the planet is positive, indicating outward migration. In the O:P(equi) disk the torque is higher than the torque for a planet in the O:P(3:1) disk (with the same viscosity). The largest differences between the two ortho-para gas states is observed for the disks with a constant viscosity.

The temperature near the planet is in the region of 80 K ≤ T_Planet ≤ 130 K, which is where the separation of the adiabatic index between O:P(equi) and O:P(3:1) first becomes pronounced for decreasing temperatures (see Fig. 1). This difference in temperature leads to a change in the total torque acting on the planet (see Fig. 9), which can be up to a factor of two in this temperature region. Therefore the total torque in the O:P(equi) disk is higher than in the O:P(3:1) disk. This effect is due to a change in the adiabatic index in the region near the planet. The origin of this effect is also described in Sect. 4.1.

Fig. 15 Torque acting on 20 MEarth planets embedded in disks with different viscosities. The graph features α = 0.004, α = 0.008, and constant viscosity for both gas ortho-para mixtures.

Figure 16 shows the torque density Γ(r). As all total torques are positive (Fig. 15), each simulation shows a spike on top of the Lindblad torque in the Γ(r) distribution just below r ≈ 1 a_P. This spike is largest for O:P(equi) in the constant viscosity simulation, which is also the simulation featuring the highest total torque. The spikes in the Γ(r) distributions are smaller for O:P(3:1) disks compared to O:P(equi) disks, which is also reflected in the total torque (Fig. 15). However, the Lindblad torque seems strongest for the α = 0.004 disk with the O:P(equi) gas.

The torque density in the O:P(equi) disk with a constant viscosity is slightly smaller around r ≈ 1.0 a_Jup compared with the constant viscosity and γ case (see Fig. 13), which leads to a smaller total torque. The torque density spike in the disk with O:P(3:1) is much smaller than in the constant γ case (both for constant viscosity again), which results in a smaller total torque. The same trend of the behavior of the torque density seems to be true for α viscosities with α = 0.008. The spike in the torque density is larger for the O:P(equi) simulations than in the O:P(3:1) simulations. This indicates a larger positive corotation torque, which results in a larger torque acting on the planet (see Fig. 15). This is a consequence of the temperature dependence of the adiabatic index, where γ has smaller values in the O:P(equi) disks, resulting in a larger total torque (see Fig. 9).

Fig. 16 Torque density of 20 MEarth planets embedded in disks with different viscosity. The graph features α = 0.004, α = 0.008, and a constant viscosity for both gas ortho-para mixtures.

In Fig. 17 the surface density for the α = 0.004 and α = 0.008 simulations with the two different ortho-para gas ratios are displayed.

Fig. 17 Surface density maps for 20 MEarth planets on fixed circular orbits in fully radiative disks. The viscosities and gas mixtures in the disks are (from top to bottom): α = 0.004 with O:P(3:1), α = 0.004 with O:P(equi), α = 0.008 with O:P(3:1), and α = 0.008 with O:P(equi). The snapshots are taken in the heating/cooling equilibrium state of the disk.

The surface density is a little bit higher in front of the planet (φ > 180 and r < 1.0 a_Jup) and a little bit lower behind the planet (φ < 180 and r > 1.0 a_Jup) for O:P(equi) compared with O:P(3:1) for both of the viscosities shown, influencing the positive corotation torque. This can also be seen in the Γ(r) distribution (Fig. 16), where the spike becomes larger, resulting in the observed higher total torque for the O:P(equi) models.

For the lower viscosity (α = 0.004), the density in the spiral waves exerted by the planet is higher compared to the high viscosity case (α = 0.008), which could also be observed for disks with a constant γ (Fig. 14). Also, the surface density near the planet is higher for the lower viscosity disks. These two features cause the higher torque in the α = 0.004 disks compared with α = 0.008 (for the same gas ortho-para ratios).

When α = 0.004, the spiral waves in the O:P(equi) disks are much denser compared with O:P(3:1) disks, but it seems to be the other way around for α = 0.008. The difference in the spiral wave density influences the height of the Lindblad torques in Fig. 16.

In total, the differences in disk structure between the two different ortho-para mixtures are larger in a constant viscosity disk than in disks with an α-viscosity. The reason for this may be in the complicated relation that arises when an α viscosity is used, because in this case ν depends explicitly on the temperature.

5. Zero-torque radius

The torque acting on an embedded 20 M_Earth planets decreases with increasing distance from the central star (Bitsch & Kley 2011b). At some point in the disk, the torque becomes zero (zero-torque radius R₀) and the planet ceases to migrate due to planet-disk interactions. The zero-torque radius R₀ is strongly dependent on planetary mass and on gradients in the disk. In addition the location of R₀ depends on the mass of the disk (Bitsch & Kley 2011b), and based on the results presented in Sect. 4, we also expect R₀ to depend on the viscosity and the equation of state of the gas.

In Fig. 18 we show the torque acting on 20 M_Earth planets placed at various radii in fully radiative disks with different viscosities. The figure also features simulations with O:P(equi), O:P(3:1), and a constant adiabatic index γ = 1.4. The data of the simulations with constant γ and constant viscosity are taken from Bitsch & Kley (2011b). Because a high viscosity is needed to keep torques that are acting on a planet from saturating, high-viscosity disks should have their R₀ at larger orbital distances than low-viscosity disks. This is indeed the case in Fig. 18, where R₀ for the α = 0.004 simulation is closer to the star than in either the α = 0.008 or ν = 10¹⁵ cm² s^-1 simulation. For these viscosities, R₀ is usually largest for the constant adiabatic index simulations and is smallest for the O:P(3:1) simulations. The most realistic treatment for the equation of state is likely O:P(3:1), so the other approximations tend to overestimate the location of R₀. Table 2 summarizes our findings.

Fig. 18 Torque acting on 20 MEarth planets in disks with different viscosities. Each plot features the two gas mixtures [O:P(3:1) and O:P(equi)] and a constant γ = 1.4. From top to bottom the viscosity is: α = 0.004, α = 0.008, and ν = 1015 cm2 s-1. The simulations with γ = 1.4 and a constant viscosity are taken from Bitsch & Kley (2011b).

To understand these results in more detail, consider Figs. 1 and 10. The torque acting on an embedded 20 M_Earth planet decreases with increasing adiabatic index. At low temperatures, the adiabatic index for O:P(3:1) is constantly increasing with increasing temperature over the range of temperatures of interest. For O:P(equi), the adiabatic index decreases until γ ≈ 1.3 and then increases until γ = 1.67 at T ≈ 20 K. As the temperature decreases with increasing distance in the disk, the result that O:P(equi) simulations have a larger R₀ than in the O:P(3:1) simulations is consistent with results from the previous sections. Using the same reasoning, the constant adiabatic index simulations should have larger R₀ than the O:P(3:1) simulations, and usually have larger R₀ than in the O:P(equi) case, although the latter is complicated by the break in monotonicity of the O:P(equi) γ temperature dependence.

6. Summary and conclusions

In this paper we have studied the disk structure due to different viscosities (constant α and constant ν, respectively) and changing adiabatic indices. We then investigated how an embedded 20 M_Earth planet migrates in such disks.

The viscosity determines the structure of the disk. We find that the midplane density (temperature) decreases (increases) for increasing α-viscosity. The aspect ratio of the disk and the slope of the temperature profile β_T also increase for increasing viscosity. The simulation with a constant viscosity of ν = 10¹⁵ cm² s^-1 seems to lie between the α = 0.005 and α = 0.006 simulations, which has allowed us to make comparisons with the earlier simulations of Kley et al. (2009) and Bitsch & Kley (2010).

The aspect ratio, H/r, of the disk is approximately independent of the index γ. Only for very small values, γ ≲ 1.2, H/r drops due to the onset of convection in the disk. At very low viscosities, the aspect ratio profile is roughly constant with small fluctuations. However, for higher viscosities, the aspect ratio profile has a complex shape, with a slight increase for increasing r at small disk radii, followed by decreasing ratio at large disk radii. A decreasing H/r profile with increasing distance is the result of the temperature dependence of the opacity.

A change in the aspect ratio influences the migration of embedded planets in radiatively cooling disks. Although we find that the direction of migration can be inward or outward, the effect of the H/r profile is broadly consistent with results from isothermal disks, where embedded planets migrate inward faster for smaller H/r (Tanaka et al. 2002). The slope of the temperature β_T also influences the torque acting on an embedded planet. In many theoretical formulae (Masset & Casoli 2010; Paardekooper et al. 2011) the parameter β_T is used to calculate the Lindblad torque and the entropy-related corotation torque. A change in these parameters therefore changes the torque acting on an embedded planet. We note that in all cases, isothermal or radiative, viscosity is required to prevent corotation torques from saturating (Masset 2001). In fully radiative disks, radiative transport/cooling acts to prevent saturation of the entropy related torque (Kley et al. 2009).

For α < 0.0025, we find that the torque acting on a 20 M_Earth planet is negative (inward migration), consistent with planets embedded in an isothermal disk (Bitsch & Kley 2010). In the low viscosity case, a low-mass planet seems to open a partial gap in the disk, such that corotation torques are weakened. For increasing viscosities (α > 0.003), the torque becomes positive and seems to stall for α ≳ 0.005. Although the density, temperature, and aspect ratio profiles for the α = 0.005 simulation match quite well with the constant viscosity simulation of an unperturbed disk, there is a difference of about 50% in the total torque acting on an embedded 20 M_Earth planet, emphasizing further that the treatment for viscosity affects planet migration. Because the α-viscosity is dependent on the temperature in the disk (${\mathit{c}}_{\mathrm{s}}\mathrm{\propto }\sqrt{\mathit{T}}$), and because the presence of the planet modifies the temperature in its neighborhood, the torques will change. Overall, this results in a smaller torque acting on the embedded planet.

In addition to viscosity, the adiabatic index influences the torque on the planet. In our test simulations with a constant adiabatic index, the torque on a planet is about a factor of three larger for an equation of state with a fixed γ = 7/5 than for a fixed γ = 5/3. For the simulations in which the adiabatic index is allowed to vary with temperature, the torque acting on embedded planets in a disk with O:P(equi) is larger than in simulations with O:P(3:1) for the same viscosity. This picture is complicated though by the connection between the viscosity and the adiabatic index. Higher viscosity results in higher midplane temperatures, which affects the temperature-dependent adiabatic index. The adiabatic index in turn alters the structure of the disk by changing for example the H/r profile. As a result, the location of the zero-torque radius is moved in the disk.

If viscosity and disk temperature is positively correlated with disk mass, we would then expect massive disks to support outward migration of embedded planets. This may be a way of transporting cores to the outer nebula and maintaining low eccentricities, which may be necessary to explain systems such as Fomalhaut (Kalas et al. 2008; Boley et al. 2012) and HR8799 (Marois et al. 2008, 2010). However, the embedded planet would need to be prevented from undergoing rapid gas accretion during the migration. We suggest that this might be an opacity effect that is exacerbated by the high disk temperatures around A stars. The higher temperature will place the zero-torque radius at larger semi-major axes, and the naturally higher mass of A star disks will increase the dust opacity, for constant metallicity. This could prevent cooling and rapid gas capture onto the planetary core until much lower surface densities are reached.

In addition to affecting massive planetary cores, the zero-torque radius, could in principle serve as a location for collecting embryos, perhaps aiding in the formation of a large core (Sándor et al. 2011). However, variations in viscosity as the disk evolves will cause the location of the zero-torque radius to also evolve. This could result in the loss of embryos and changes that result to inward migration, or could open the possibility of providing new feeding zones for young planets. The overall influence of the zero-torque radius on the formation of large cores is further complicated by collisions between planetesimals in the envisaged planet trap, which may influence the opacity. Such a change will affect the H/r profile, which ultimately influences the zero-torque radius.

The O:P(3:1) equation of state is the most realistic for disks explored here. This treatment for the equation of state also leads to the smallest zero-torque radii compared with both the constant γ = 1.4 and the O:P(equi) simulations. The effect is most pronounced in high-viscosity disks, and the inclusion of a proper equation of state is crucial for determining the location of possible planet traps, where planetary embryos or cores can collect.

Acknowledgments

Bitsch has been sponsored through the German D-grid initiative and through the Helmholtz Alliance Planetary Evolution and Life. 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. Boley’s support was provided by a contract with the California Institute of Technology (Caltech) funded by NASA through the Sagan Fellowship Program. We would also like to thank an anonymous referee for the useful comments that helped to improve the manuscript.

References

All Tables

All Figures

	Fig. 1 Adiabatic index for the 3:1 ortho to para mix [O:P(3:1)] and equilibrium gas [O:P(equi)] configuration in dependence of the temperature.
In the text

	Fig. 2 Densities (top), temperatures (middle), and aspect ratios (bottom) from the 2D axisymmetric disks with different (constant) adiabatic indices, taken after the disks have reached their equilibrium state. Each quantity is evaluated using midplane values of the disk at a radius r = 1 aJup. The kinematic viscosity has been set to ν = 1015 cm2 s-1 for all models.
In the text

	Fig. 3 Aspect ratio (H/r) of a fully radiative disk with different adiabatic indices without an embedded planet. The aspect ratio was computed using midplane values of the constant viscosity disks (see Eq. (2)).
In the text

	Fig. 4 Density (top) and temperature (bottom) profiles of the fully radiative disks with different viscosities, all taken in the disk equilibrium state. Densities and temperatures are measured in the midplanes of the disks.
In the text

	Fig. 5 Midplane density (top) and temperature (bottom) profiles for the α viscosity disks with different assumptions for the ortho-para ratio of the gas.
In the text

	Fig. 6 Aspect ratio H/r of the accretion disks for different α viscosities and different assumptions for the ortho-para ratio of the gas.
In the text

Fig. 7 Velocities in z-direction for disks with α = 0.001 and O:P(3:1) (top), α = 0.001 and O:P(equi) (2nd from top), α = 0.004 and O:P(equi) (3rd from top), and a fixed γ = 1.1 with constant viscosity (bottom). The black line indicates H(r) in the disk and the red line delineates τ = 1.0 (the optical depth integrated from the top of the grid toward the midplane of the disk). Please note the different color scales and that they only show a subset of the range of values to better highlight the internal structure.

In the text

	Fig. 8 Adiabatic gradients for the γ = 1.1 disk with constant viscosity (top) and for the α = 0.004 disk with O:P(equi) (bottom). The gradients have been evaluated in the optical thick regions of the disk at an angle of 2° from the midplane.
In the text

	Fig. 9 Total torque Γtot acting on a 20 MEarth planet embedded at r = 1 in ν = 1015 cm2 s-1 disks for different (constant) adiabatic indices. The torques have been measured in the equilibrium state of the disk.
In the text

	Fig. 10 Profiles of the radial torque density Γ(r) acting on 20 MEarth planets embedded in disks with different adiabatic indices (during the disk equilibrium state). In contrast to the total torque, the radial torque density has not been averaged in time.
In the text

	Fig. 11 Surface density maps for planets on fixed circular orbits in fully radiative disks with different adiabatic indices (top: γ = 1.1, middle: γ = 1.4, bottom: γ = 1.7) and a constant viscosity of ν = 1015 cm2 s-1. The disruption in the surface density pattern for the γ = 1.1 disk is due to convection inside the disk.
In the text

	Fig. 12 Torque acting on 20 MEarth planets embedded in disks with different viscosities (γ = 1.4).
In the text

	Fig. 13 Profiles of the radial torque density Γ(r) acting on 20 MEarth planets embedded in disks with different viscosities (in the disk equilibrium state) and a constant adiabatic index (γ = 1.4).
In the text

	Fig. 14 Surface density maps for 20 MEarth planets on fixed circular orbits in fully radiative disks (with γ = 1.4). The viscosities in the disks are (from top to bottom): α = 0.001, α = 0.005, and α = 0.008. The snapshots are taken in the equilibrium state.
In the text

	Fig. 15 Torque acting on 20 MEarth planets embedded in disks with different viscosities. The graph features α = 0.004, α = 0.008, and constant viscosity for both gas ortho-para mixtures.
In the text

	Fig. 16 Torque density of 20 MEarth planets embedded in disks with different viscosity. The graph features α = 0.004, α = 0.008, and a constant viscosity for both gas ortho-para mixtures.
In the text

	Fig. 17 Surface density maps for 20 MEarth planets on fixed circular orbits in fully radiative disks. The viscosities and gas mixtures in the disks are (from top to bottom): α = 0.004 with O:P(3:1), α = 0.004 with O:P(equi), α = 0.008 with O:P(3:1), and α = 0.008 with O:P(equi). The snapshots are taken in the heating/cooling equilibrium state of the disk.
In the text

	Fig. 18 Torque acting on 20 MEarth planets in disks with different viscosities. Each plot features the two gas mixtures [O:P(3:1) and O:P(equi)] and a constant γ = 1.4. From top to bottom the viscosity is: α = 0.004, α = 0.008, and ν = 1015 cm2 s-1. The simulations with γ = 1.4 and a constant viscosity are taken from Bitsch & Kley (2011b).
In the text