Influence of viscosity and the adiabatic index on planetary migration
Key Words.:accretion disks – planet formation – hydrodynamics – radiative transport – planet disk interactions – adiabatic index
Context: The strength and direction of migration of low mass embedded planets depends on the disk’s thermodynamic state. It has been shown that in active disks, where the internal dissipation is balanced by radiative transport, migration can be directed outwards, a process which extends the lifetime of growing embryos. Very important parameters determining the structure of disks, and hence the direction of migration, are the viscosity and the adiabatic index.
Aims: In this paper we investigate the influence of different viscosity prescriptions (-type and constant) and adiabatic indices on disk structures. We then determine how this affects the migration rate of planets embedded in such disks.
Methods: We perform three-dimensional numerical simulations of accretion disks with embedded planets. We use the explicit/implicit hydrodynamical code NIRVANA that includes full tensor viscosity and radiation transport in the flux-limited diffusion approximation, as well as a proper equation of state for molecular hydrogen. The migration of embedded 20 planets is studied.
Results: Low-viscosity disks have cooler temperatures and the migration rates of embedded planets tend toward the isothermal limit. Hence, in these disks, planets migrate inwards even in the fully radiative case. The effect of outward migration can only be sustained if the viscosity in the disk is large. Overall, the differences between the treatments for the equation of state seem to play a more important role in disks with higher viscosity. A change in the adiabatic index and in the viscosity changes the zero-torque radius that separates inward from outward migration.
Conclusions: For larger viscosities, temperatures in the disk become higher and the zero-torque radius moves to larger radii, allowing outward migration of a 20- planet to persist over an extended radial range. In combination with large disk masses, this may allow for an extended period of the outward migration of growing protoplanetary cores.
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; Baruteau & Masset 2008; Paardekooper & Papaloizou 2008; Paardekooper & Mellema 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; Masset & Casoli 2010; Paardekooper et al. 2011), 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 and . The simulation of an embedded planet in a global magnetorotationally unstable accretion disk through direct magnetohydrodynamical 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 towards 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 K (Decampli et al. 1978; Boley et al. 2007). As a result, the adiabatic index of the gas will transition from to 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 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 Section 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 Section 3. These structural changes influence the migration rate of an embedded planet, which is discussed in Section 4. In Section 5 the influence of viscosity and of the different gas mixtures on the zero-torque radius are investigated. We then summarize and conclude in Section 6.
2 Numerics and 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, 2011).
The three-dimensional () computational domain (with active cells) consists of a complete annulus of the protoplanetary disk centered on the star, extending from to in units of . For simulations with planets at larger distances from the central star, the computational grid is set to have an inner computational boundary at and an outer boundary at . In the vertical direction the annulus extends above the disk’s midplane, meaning . As we assume that the disk is symmetric for the upper and lower parts, we apply symmetric boundary conditions at . Here denotes the polar angle of the used spherical polar coordinate system measured from the polar axis. The central star has one solar mass , and the total disk mass inside  is . 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 - directions.
The planet is then embedded at . For the planet’s gravity, we use the cubic potential (Klahr & Kley 2006; Kley et al. 2009) with . 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 plane) reduces the required computational effort substantially. The evolution from the initial isothermal state towards 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., ) 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 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 to 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 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 ) the temperature at is about 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.
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 , where 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:
where is typically (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 and . As will be discussed later in the manuscript, we define the half-disk thickness as
where is the isothermal sound speed, and 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 2011). 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 () is needed for the calculation, as the lower half is directly symmetric to the upper half for these conditions.
The torques acting on a 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:
Here is the planetary mass, denotes the distance of the disk element to the planet, and is the smoothing length of the potential, usually a fraction of the Hill radius, . The construction of the planetary potential is in such a way that, for distances larger than , the potential matches the correct potential. For , the potential is smoothed by a cubic polynomial. The parameter is equal to 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 .
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.
3.1 Influence of a constant adiabatic index
For the following simulations, we limit ourselves to constant values of in the range of , which is larger than the range of the adiabatic index in Fig. 1. We include values smaller than 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 cm/s. 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 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 . 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 .
Radial density and temperature gradients result in vortensity and entropy () 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 ), 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 , the aspect ratio is smaller at 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 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, ). Such a feature in the disk’s structure is the opposite of what one would expect. When taking stellar irradiation into account an increasing profile for larger radii may be possible, e.g. Dullemond & Dominik (2004). However, the disk’s opacity can greatly influence the profile.
With increasing distance from the central star, the temperature is decreasing. For temperatures less than 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 , the opacity decreases. This is exactly the point in the disk () where the aspect ratio begins to decrease. Please note that the height of the disk is still increasing, but at a lower rate, when K compared to K. Hence, the structure of the disk is a result of the opacity law that is used in these calculations.
The profiles for the disks with 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.
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 and for a constant viscosity of cm/s, a value used in our previous simulations (Kley et al. 2009; Bitsch & Kley 2010, 2011).
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 cm/s, the radial slope of temperature was at , while it now depends on the chosen -value. The slope of temperature , where , is displayed in Table 1 for .
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 therefore influences the torque acting on an embedded planet quite severely (Ayliffe & Bate 2011).
The aspect ratio 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 profiles of our constant viscosity simulation match quite well with the 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 , 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 there are differences in the density and temperature profiles between disks that have different ortho-para ratios. For , 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 K for ortho-para ratios O:P(equi), this behavior is consistent with the results of the constant simulations.
In Fig. 6 the aspect ratios 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 case, the aspect ratio profile is smaller and approximately constant. For the more viscous disks, the aspect ratio increases for small , but then decreases with increasing 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 compared to the O:P(3:1). The aspect ratio differs by about for the high viscosity cases at . This difference at this location can influence the torque acting on embedded planets, as seen in isothermal disks.
A higher disk viscosity leads to a larger temperature in the midplane. Recently, Bitsch & Kley (2011) 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 (see Fig. 4).
We display in Fig. 7 the velocity in the -direction for the disks with two different ortho-para gas mixtures, for the and O:P(equi) disk, and for the fixed disk. In black we indicate and in red the 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 (2011) 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 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 ) 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 discs 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 , 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 (), the fluctuations in the upper layers mostly vanish, except in the low-density region at large , and are totally diminished for (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 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 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:
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 , convection can be present in the disk, where is defined as
Here, 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 disk (bottom in Fig. 7) and the disk with O:P(equi) (3rd from top in Fig. 7).
In the case of , , 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 is about a factor of . For the disk, only the inner regions feature . However, fluctuations in the velocity have not been observed (Fig. 7). In addition, the difference between and is at maximum a factor of , which is much lower than in the case. If convection is present in the disk, then it must be much less vigorous than seen in the 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 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 ), the torques have been averaged over planetary orbits.
4.1 Constant adiabatic index
In Fig. 9 we show the total torque acting on planets that are on circular orbits and are embedded in disks that have different adiabatic indices. The viscosity is held constant at cm/s. The torque acting on the planet is positive when the adiabatic index is larger than , while it is negative, indicating inward migration, for . The maximum of the torque is at , and the torque decreases for lower and higher ’s. The most negative torque is at .
For , 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 to compared to the maximum torque at . 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 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.
In Fig. 10, the radial torque density , which is the torque exerted by an annulus of width and related to the total torque by
is displayed for planets embedded in fully radiative disks with a constant adiabatic index. All displayed torque densities (with ) feature a spike just below 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 , 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. 11 shows surface density images of our fully radiative disk simulations with an embedded . The adiabatic indices for these images are (top), (middle) and (bottom). We note that the case is discussed in great detail in Kley et al. (2009). For the 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 disk, on the other hand, shows no sign of convection. In fact it is very similar to the (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 , 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).
4.2 Influence of viscosity
In Section 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 planet in disks with different viscosities. Planets in disks with (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 , the torque acting on the planet becomes positive, and for even larger viscosities (), the magnitude of the torque begins to converge. The torque acting on the planet with constant viscosity of cm/s is similar to the one with . This is consistent with the similarities between their temperature and density profiles (Fig. 4).
In Fig. 13 we display the radial torque density acting on a planet in disks with different viscosity. In the 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 and torque density patterns are consistent with the that of the constant viscosity of cm/s. The spike in the distribution near is nearly identical for the constant viscosity and simulations, while the spike is a little bit higher. However, there are some differences. For smaller and larger distances from the central star ( and ) the torque density is larger for the constant viscosity case compared to , although all simulations follow the same trend. The underlying Lindblad torque distribution is also smaller for the disk, due to a difference in the overall disk structure (see Fig. 4). The reduction of 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 , 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).
In Fig. 14 we display the surface densities of several disks with an embedded planet. Different viscosities are shown to gain more insight into the origin of the torque density profiles. For very low viscosities (), 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.
The 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 ( and ) 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 ( and ) is not so clear in the simulation compared to the cm/s 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 simulation at that distance from the central star, thus explaining the higher total torque.
For , 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 ) 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 planets in disks is displayed for different viscosity and gas treatments. For low viscosities () 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 , 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 Section 4.1.
Fig. 16 shows the torque density . As all total torques are positive (Fig. 15), each simulation shows a spike on top of the Lindblad torque in the distribution just below . 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 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 disk with the O:P(equi) gas.
The torque density in the O:P(equi) disk with a constant viscosity is slightly smaller around 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 . 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).
In Fig. 17 the surface density for the and simulations with the two different ortho-para gas ratios are displayed.
The surface density is a little bit higher in front of the planet ( and ) and a little bit lower behind the planet ( and ) 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 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 (), the density in the spiral waves exerted by the planet is higher compared to the high viscosity case (), 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 disks compared with (for the same gas ortho-para ratios).
When , 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 . 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 planets decreases with increasing distance from the central star (Bitsch & Kley 2011). At some point in the disk, the torque becomes zero (zero-torque radius ) and the planet ceases to migrate due to planet-disk interactions. The zero-torque radius is strongly dependent on planetary mass and on gradients in the disk. In addition the location of depends on the mass of the disk (Bitsch & Kley 2011), and based on the results presented in Section 4, we also expect to depend on the viscosity and the equation of state of the gas.
In Figure 18 we show the torque acting on 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 . The data of the simulations with constant and constant viscosity are taken from Bitsch & Kley (2011). Because a high viscosity is needed to keep torques that are acting on a planet from saturating, high-viscosity disks should have their at larger orbital distances than low-viscosity disks. This is indeed the case in Figure 18, where for the simulation is closer to the star than in either the or cm/s simulation. For these viscosities, 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 . Table 2 summarizes our findings.
To understand these results in more detail, consider Figures 1 and 10. The torque acting on an embedded 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 and then increases until at . As the temperature decreases with increasing distance in the disk, the result that O:P(equi) simulations have a larger 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 than the O:P(3:1) simulations, and usually have larger 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 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 also increase for increasing viscosity. The simulation with a constant viscosity of cm/s seems to lie between the and simulations, which has allowed us to make comparisons with the earlier simulations of Kley et al. (2009) and Bitsch & Kley (2010).
The aspect ratio, , of the disk is approximately independent of the index . Only for very small values, , 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 at small disk radii, followed by decreasing ratio at large disk radii. A decreasing 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 profile is broadly consistent with results from isothermal disks, where embedded planets migrate inward faster for smaller (Tanaka et al. 2002). The slope of the temperature also influences the torque acting on an embedded planet. In many theoretical formulae (Masset & Casoli 2010; Paardekooper et al. 2011) the parameter 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 , we find that the torque acting on a 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 (), the torque becomes positive and seems to stall for . Although the density, temperature, and aspect ratio profiles for the simulation match quite well with the constant viscosity simulation of an unperturbed disk, there is a difference of about in the total torque acting on an embedded planet, emphasizing further that the treatment for viscosity affects planet migration. Because the -viscosity is dependent on the temperature in the disk (), 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 than for a fixed . 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 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 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 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.
Acknowledgements.B. 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. A. C. 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.
- offprints: B. Bitsch,
- Ayliffe, B. A. & Bate, M. R. 2011, ArXiv e-prints
- Ayliffe, B. A. & Bate, R. 2010, MNRAS, 408, 876
- Baruteau, C. & Masset, F. 2008, ApJ, 672, 1054
- Bitsch, B. & Kley, W. 2010, A&A, 523, A30
- Bitsch, B. & Kley, W. 2011, A&A, 530, A41
- Bitsch, B. & Kley, W. 2011, A&A, 536, A77
- Boley, A. C., Hartquist, T. W., Durisen, R. H., & Michael, S. 2007, ApJ, 656, L89
- Boley, A. C., Payne, M. J., Corder, S., et al. 2012, ApJ, 750, L21
- Crida, A., Baruteau, C., Kley, W., & Masset, F. 2009, A&A, 502, 679
- Crida, A., Sándor, Z., & Kley, W. 2008, A&A, 483, 325
- Decampli, W. M., Cameron, A. G. W., Bodenheimer, P., & Black, D. C. 1978, ApJ, 223, 854
- Dullemond, C. P. & Dominik, C. 2004, A&A, 417
- Flaig, M., Ruoff, P., Kley, W., & Kissmann, R. 2012, MNRAS
- Gullbring, E., Hartmann, L., Briceno, C., & Calvet, N. 1998, ApJ, 492, 323
- Hartigan, P., Edwards, S., & Ghandour, L. 1995, ApJ, 452, 736
- Kalas, P., Graham, J. R., Chiang, E., et al. 2008, Science, 322, 1345
- Klahr, H. & Kley, W. 2006, A&A, 445, 747
- Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971
- Kley, W. & Crida, A. 2008, A&A, 487, L9
- Kley, W., Müller, T. A. W., Kolb, S. M., Benítez-Llambay, P., & Masset, F. 2012, A&A, 546, A99
- Levermore, C. D. & Pomraning, G. C. 1981, ApJ, 248, 321
- Lin, D. N. C. & Papaloizou, J. C. B. 1985, in Protostars and Planets II, 981
- Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348
- Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010, Nature, 468, 1080
- Masset, F. 2000, A&AS, 141, 165
- Masset, F. & Casoli, J. 2010, ApJ, 723
- Masset, F. S. 2001, ApJ, 558, 453
- Masset, F. S. 2002, A&A, 387, 605
- Morbidelli, A., Crida, A., Masset, F., & Nelson, R. 2008, A&A, 478, 929
- Nelson, R., Gressel, O., & Umurhan, O. 2012, ArXiv e-prints
- Paardekooper, S., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950+
- Paardekooper, S. J., Baruteau, C., & Kley, W. 2011, MNRAS, 293
- Paardekooper, S.-J. & Mellema, G. 2006, A&A, 459, L17
- Paardekooper, S.-J. & Mellema, G. 2008, A&A, 478, 245
- Paardekooper, S.-J. & Papaloizou, J. C. B. 2008, A&A, 485, 877
- Ruden, S. P. & Pollack, J. B. 1991, ApJ, 375, 740
- Sándor, Z., Lyra, W., & Dullemond, C. P. 2011, ApJ, L9
- Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
- Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257
- Ziegler, U. & Yorke, H. 1997, Computer Physics Communications, 101, 54