Planetary migration with steep temperature gradients

Migration of protoplanets with surfaces through discs with steep temperature gradients

Ben A. Ayliffe and Matthew R. Bate
School of Physics, University of Exeter, Stocker Road, Exeter EX4 4QL
Centre for Stellar and Planetary Astrophysics & School of Mathematical Sciences, Monash University, Melbourne Vic. 3800, Australia
July 12, 2019

We perform three-dimensional self-gravitating radiative transfer simulations of protoplanet migration in circumstellar discs to explore the impact upon migration of the radial temperature profiles in these discs. We model protoplanets with masses ranging between 10-100 , in discs with surface density profiles of , and temperature profiles of the form , where ranges . We find that steep () temperature profiles lead to outward migration of low mass protoplanets in interstellar grain opacity discs, but in more optically thin discs the migration is always inwards. The trend in migration rates with changing obtained from our models shows good agreement with those obtained using recent analytic descriptions which include consideration of the co-orbital torques and their saturation. We find that switching between two models of the protoplanet, one in which accretion acts by evacuating gas and one in which gas piles up on a surface to form an atmosphere, leads to a small shift in the migration rates. If comparing these models in discs with conditions which lead to a marginally inward migration, the small shift can lead to outward migration. However, the direction and speed of migration is dominated by disc conditions rather than by the specific prescription used to model the flow near the protoplanet.

planets and satellites: formation – radiative transfer – methods: numerical – hydrodynamics – planetary systems: formation

1 Introduction

Investigations into the interactions of circumstellar discs and embedded protoplanets were first conducted by Goldreich & Tremaine (1980). They found that the exchange of angular momentum between the two should lead to migration of the protoplanet (latterly known as Type I migration), but their work did not suggest in which direction this migration would be. The analysis of Ward (1986) for a non-self gravitating, two-dimensional disc, showed that if the disc had a negative temperature gradient (i.e. temperature reduces with increasing heliocentric radius) then an embedded protoplanet should migrate inwards. When Hot Jupiters were discovered (Mayor & Queloz, 1995) and the difficulties with their in situ formation became apparent, this inward migration seemed to be an ideal mechanism by which to reconcile existing formation scenarios with these planets’ small final orbital radii. However, the timescales of growth and migration did not at first compare favourably. Migration timescales were so short that protoplanets should plummet into their stars before they were able to grow to any considerable mass. Only by reaching masses sufficient to open disc gaps can protoplanets move from fast Type I migration to slower Type II migration, which proceeds at the disc’s viscous evolution timescale. Once slowed, the planets only have to survive until the disc gas dissipates, thought to occur at a stellar age of less than 6 Myrs (Haisch, Lada, & Lada, 2001), at which point migration due to planet-gas interactions necessarily ceases.

Since Ward (1986) further analytical descriptions (i.e. Ward 1997; Tanaka, Takeuchi, & Ward 2002), and numerical modelling (Korycansky & Pollack 1993; Nelson et al. 2000; Masset 2002; D’Angelo, Henning, & Kley 2002; Bate et al. 2003; D’Angelo, Kley, & Henning 2003; Alibert, Mordasini, & Benz 2004; D’Angelo, Bate, & Lubow 2005; Klahr & Kley 2006; D’Angelo & Lubow 2008; Li et al. 2009) of planet migration have continued to find fast inward rates in the Type I regime ( 100 ). However, these works have generally considered locally-isothermal conditions. A large number of works, both analytical and numerical, have now been devoted to exploring the impact of more complex thermodynamics upon planet migration. The first of such works was published by Morohoshi & Tanaka (2003), who performed local calculations focussing on a low-mass planet’s interaction with its parent disc. They found that radiative cooling led to a non-axisymmetric mass distribution about the planet. This resulted in an additional torque beyond the commonly considered differential Lindblad torques, and altered the migration rate relative to a purely isothermal calculation. Menou & Goodman (2004) found that a planet’s migration rate was sensitive to the temperature and opacity structure of the disc through which it travelled, and under certain conditions the migration rate could be very much slower than that achieved in an isothermal model. This was closely followed by Jang-Condell & Sasselov (2005) who also identified a slow down in Type I migration rates upon the introduction of radiative transfer.

It was Paardekooper & Mellema (2006) who first found cases of outward migration due to a coorbital torque in their three-dimensional radiative transfer models. They found that this coorbital torque could be prevented from saturating if the radiation diffusion timescale was shorter than the libration timescale of gas in the horseshoe region. This enabled temperature asymmetries to be maintained beyond a single libration period. They also showed that reducing the opacity (further shortening the radiation diffusion timescale) of their disc could lead to results in agreement with previous isothermal calculations. This was further described by Baruteau & Masset (2008), who make clear that the radiation diffusion timescale must be greater than the period of a single horseshoe orbit to avoid returning to an isothermal like migration, and that outward migration can be related to the disc’s radial entropy gradient. Many numerical models have since also found evidence of outward migration in the Type I regime (Paardekooper & Mellema 2008; Paardekooper & Papaloizou 2008; Kley & Crida 2008; Kley, Bitsch, & Klahr 2009; Ayliffe & Bate 2010; Yamada & Inaba 2011). Periods of outward migration can help to increase the overall migration timescale of a forming planet embedded in a disc, as is required by synthesis models to explain the population of exoplanets that has been observed (Ida & Lin, 2008; Mordasini, Alibert, & Benz, 2009).

More recently there has been an extension of analytical descriptions of the non-linear coorbital torque, called the horseshoe drag, to describe it in both its unsaturated and saturated states (Paardekooper et al. 2010; Paardekooper, Baruteau, & Kley 2011; Masset & Casoli 2010). These works give expressions for the total torque acting on a planet due to both Lindblad torques and the coorbital component. As such they can describe planet migration for a large range of scenarios, taking into account the disc’s density and temperature profiles, as well as its thermal diffusivity. Outward migration is fastest in discs with steep radial temperature profiles, becoming more marginal at typical gradients such as .

This paper is intended to explore protoplanet migration in a series of discs with radial temperature gradients from to . We conduct three dimensional global disc models, using smoothed particle hydrodynamics (SPH), that include self-gravity and which use a planetary surface to allow modelling of gas flow to well within the Hill sphere and the self-consistent formation of an atmosphere. We conduct a few isothermal models to investigate whether or not the temperature profile alone has an impact on the migration rate, but otherwise our calculations all include radiative transfer using a flux-limited diffusion approximation. In Section 2, we describe our computational method, in Section 3 we explain how we obtained our results which are then presented in Section 4. Section 5 discusses these results, whilst a summary and our conclusions are given in Section 6.

2 Computational Method

The calculations described herein have been performed using a three-dimensional SPH code. This SPH code has its origins in a version first developed by Benz (1990; Benz et al. 1990) but it has undergone substantial modification in subsequent years. Energy and entropy are conserved to timestepping accuracy by use of the variable smoothing length formalism of Springel & Hernquist (2002) and Monaghan (2002) with our specific implementation being described in Price & Bate (2007). Gravitational forces are calculated and neighbouring particles are found using a binary tree. Radiative transfer is modelled in the two temperature (gas, , and radiation, ) flux-limited diffusion approximation using the method developed by Whitehouse, Bate, & Monaghan (2005) and Whitehouse & Bate (2006). Integration of the SPH equations is achieved using a second-order Runge-Kutta-Fehlberg integrator with particles having individual timesteps (Bate, 1995). The code has been parallelised by M. Bate using OpenMP.

2.1 Equation of state and radiative transfer

We present a few calculations performed using a locally-isothermal equation of state, as well as many more calculations which include radiative transfer. In locally-isothermal models the temperature of the gas in the disc remains as a function of radius throughout the calculations. For the radiation hydrodynamical calculations we use the ideal gas equation of state, where is the gas constant, is the density, is the gas temperature, and is the mean molecular mass. The equation of state takes into account the translational, rotational, and vibrational degrees of freedom of molecular hydrogen (assuming a 3:1 mix of ortho- and para-hydrogen; see Boley et al. 2007). It also includes the dissociation of molecular hydrogen, and the ionisations of hydrogen and helium. The hydrogen and helium mass fractions are and , respectively, whilst the contribution of metals to the equation of state is neglected. More details on the implementation of the equation of state can be found in Whitehouse & Bate (2006).

The radiative transfer in these calculations is performed using the flux-limited diffusion approximation, as implemented by Whitehouse et al. (2005) and Whitehouse & Bate (2006), in which work and artificial viscosity (including both bulk and shear components) increase the thermal energy of the gas. Work done on the radiation field increases the radiative energy which can be transported via flux-limited diffusion. The energy transfer between the gas and radiation fields is dependent upon their relative temperatures, the gas density, and the opacity, . Energy is lost from the system by the radiation field into the vacuum surrounding the disc; this is performed numerically through the use of a radiation boundary. This boundary is positioned at a height from the disc midplane at which the radiation path through the gas from this height outwards has an optical depth of ; i.e. the radiation can be expected to reach the vacuum. As a result, two layers of particles, one above and one below the disc, are made to comprise the radiation boundary. These particles are compelled to follow the initial temperature profile of the disc, causing them to function as energy sinks which impose a minimum temperature to which the disc can cool (see Ayliffe & Bate, 2010, for a fuller description).

The use of steep radial temperature profiles in some of our models in this work leads to very high temperatures towards the inner boundary of these circumstellar discs. This material forms an unrealistically hot annulus at the inner boundary, from which radiation diffuses out through the disc, changing its structure out to the environs of the embedded protoplanet. To remedy this situation, the disc out to 2 au is compelled to maintain its initial temperature profile throughout the simulation, which quashes any diffusion in this region. This rule was applied in all the calculations for consistency, including the shallower temperature profile calculations where it has no discernible effect.

The opacities used here are those of Pollack et al. (1985) and Alexander (1975) (the IVa King model), with the former providing the grain opacities and the latter the gas opacities at temperatures beyond the grain sublimation point. In some calculations the grain opacity is reduced by a factor of 100 to emulate possible modification of the population due to agglomeration processes (see Ayliffe & Bate, 2009, for further details).

2.2 Disc models

We model a protoplanetary disc with radial bounds of 0.1 - 3  (0.52 - 15.6 au), where  is the initial orbital radius of our embedded protoplanets, taking a value of 5.2 au. The disc is represented by 2 million SPH particles, a number found to deliver satisfactory resolution (see Ayliffe & Bate, 2010). This leads to smoothing lengths at the disc midplane at  of (away from the planet), where is the disc scaleheight and equals 0.05 at . Of course, considerably better resolution is obtained directly surrounding the protoplanet where the densities become much higher. The surface density of the disc goes as , where the value at  is 75  to allow comparison with previous works (e.g. Lubow, Seibert, & Artymowicz 1999; Bate et al. 2003). At the centre of the disc is a fixed potential representing a 1  star. We implement an inner boundary for the disc that prevents gas flow onto the star, which due to the lack of a self-limiting mechanism would otherwise artificially rapidly drain material from the disc. The outer edge of the disc is bordered by ghost particles which represent the disc beyond 15.6 au, preventing the disc from shear spreading into a vacuum (for more details see Ayliffe & Bate, 2010). The initial discs were evolved in the absence of a planet until any transience resulting from settling had dissipated, which required just over 4 orbits of the disc’s outer edge. Self-gravity is included in all the calculations described in this paper. The impact of including self-gravity was discussed in Ayliffe & Bate (2010), where for these relatively low mass discs it was found to make only a marginal difference to migration rates, slightly slowing inward migration. However, self-gravity is essential for building a self-consistent atmosphere within the Hill radius.

Figure 1: Internal energy (, top panels) and temperature (, bottom panels) profiles along the midplane of unperturbed discs. The dashed lines show the profiles that are given when it is that is established with the given radial dependence. Conversely, the solid lines show the profiles that are produced when it is which is compelled to follow the stated profile. For the shallow profile given in the left panels, the temperature profile that results from setting is similar to that obtained by forcing . However, where a steeper profile is used, as in the right panels, the temperatures reached are sufficient to bring about substantial changes in the specific heat capacity of the gas (i.e. due to dissociation of molecular hydrogen, excitation of rotational and vibrational modes of molecular hydrogen) such that the proportionality of energy and temperature breaks down; we use heat capacities calculated using Boley et al. (2007). As such in this work it is the temperature profile that is set directly (), for which the corresponding energy is then found.

We perform calculations with a number of radial temperature profiles described by where 0, 0.5, 1, 1.5, and 2, and 73 K. Due to the use of steeper temperature profiles in this work, which results in higher temperatures at smaller heliocentric radii, the method used to set the disc temperatures is different to that used in Ayliffe & Bate (2010). In this previous work the initial temperature distribution was set assuming that the heat capacity at constant volume was independent of temperature, such that , where is the gas’s specific internal energy. This is a good approximation so long as the rotational and vibrational modes of are not excited (i.e. ), which is true for K. It was then the distribution of that was actually prescribed to follow an profile, such that ; this case is shown by the dashed lines in the left panels of Fig. 1. The resulting temperature profile is similar to the desired profile as can be seen by comparing the solid line in the lower left panel, for which is set explicitly, with the dashed line that results from setting . The higher temperatures in some of the disc models considered in this work fall beyond the range in which this assumption can readily be made, as shown by the right hand panels of Fig. 1. As such, for these models we set the initial disc temperatures and boundary temperatures to follow the desired profile (the solid lines in lower panels of Fig. 1) and then the appropriate internal energy for the gas is found (solid lines in top panels of Fig. 1) using tabulated heat capacities based upon Boley et al. (2007) which are used in our radiative transfer implementation. We note that at the planet’s orbital radius the disc’s temperature is such that . Therefore we adopt this value for when making comparisons with analytic models.

The viscosity within the disc is introduced in a parameterised form, developed for SPH by Monaghan & Gingold (1983), which conserves linear and angular momentum, and was modified to deal with high Mach number shocks by Monaghan (1992). It is a function of two parameters, the -term establishes a viscosity to damp subsonic velocity oscillations that may be produced in the wake of a shock front, whilst the -term prevents particle interpenetration in supersonic shocks. This viscosity is only applied when the gas is under compression, and should ideally only act near shock fronts. We use the switch developed by Morris & Monaghan (1997) to try to reduce the action of artificial viscosity where the cause is not a shock. This attributes an -value to each particle, which is modified based on the local pressure gradient, such that it increases to a maximum of 1 at a shock, and falls away rapidly to 0.1 as the particle moves away; this significantly reduces unwanted viscous dissipation in the disc.

3 Calculations

3.1 Planet model and migration measurements

Into the protoplanetary discs we embed a protoplanetary core, or partially formed protoplanet. In all the calculations discussed here the protoplanet is represented by a point mass which is free to move, and in all but a few cases (discussed below) the point mass is enwrapped by a surface of radius 0.03 times the Hill radius () (Ayliffe & Bate, 2009, 2010). Gas is able to pile up on this surface to build an atmosphere, and so represent gas accretion in a reasonably realistic manner. As in our previous work, a smooth start to the calculations is achieved by embedding a protoplanet of radius which then shrinks exponentially to the desired radius during the first orbit of the protoplanet.

The point mass is free to move through the disc, enabling us to measure its rate of migration. This measurement is made by using a linear fit to the last 25 orbits of orbital radius evolution, in what is typically a 50 orbit total evolution. The first 25 orbits are given up to ensure that any transient disruption to the disc caused by the sudden introduction of a planet have subsided, and to allow gas in the horseshoe region to reach a quasi-equilibrium state. For a 10  planet the libration time is orbits, for a 333  planet it is orbits, whilst for the a 33  planet it is 29 orbits. As such, our lowest mass protoplanet mass models have completed a single libration period at the end of the calculations. In Ayliffe & Bate (2010) we presented our results using migration timescales, , which are calculated as the time a planet of fixed mass would take to migrate from its starting orbital radius (5.2 AU in all of the models discussed) in to the central star; . This has no real meaning if the protoplanet is migrating outwards, so in this paper we present results using simply the migration rate, , in units of per orbital period at (which for  au is 11.86 years).

We further refine the presentation of our results by quantifying the uncertainties that arise due to our method of measuring the migration rates. To achieve this, in addition to fitting over the entirety of the last 25 orbits, we also divide this interval into 3 smaller overlapping time domains to which individual linear fits are made. For a migration trend that is well represented by a linear fit, these 3 migration rates should all be similar to both one another and to the rate measured over the entire 25 orbits. For less linear trends the measured rates give an indication of the possible variation that could be obtained by measuring at different points in the orbital evolution. We use the maximum and minimum migration rates obtained from these shorter time domain fits to provide error bars for each migration rate that is presented in this work.

In Ayliffe & Bate (2010) we noted a difference between the migration rates measured when modelling a protoplanet with the surface described above and an accreting point mass. To explore this further we conduct a small number of calculations in which the protoplanet is represented as a point mass which accretes gas that comes within , where in each case. The particles representing this accreted gas are removed from the calculation and their properties are used to modify those of the point mass, ensuring conservation of mass, linear momentum, and angular momentum (see Bate, Bonnell, & Price, 1995).

3.2 Opacity to diffusivity

Masset & Casoli (2010) developed an analytic expression to calculate the torques acting on protoplanets embedded in discs with temperature profiles described by the relation , such as the discs we have modelled in this work. They also include terms that describe the impact on the torques of the disc’s thermal diffusivity, which in our radiative transfer calculations is effectively varied through changes in the opacity. In order that we might compare our numerical models with Masset & Casoli’s analytic expression, we calculate the thermal diffusivity in our disc at  in the absence of a protoplanet. In the flux-limited diffusion approximation of radiative transfer used in this work the flux in optically thick regions is given by


where is the speed of light, is the Rosseland mean opacity, is the gas density, and in which is the radiation density constant, . This allows us to restate equation 1 as


Masset & Casoli give an expression for the heat flux in terms of a form of thermal diffusivity, , as


Comparing equations 2 and 3 we obtain an expression for , which is related to a true thermal diffusivity, , by


which is similar to equation 33 in Masset & Casoli (2010), where is the adiabatic index. The difference between this equation and that of Masset & Casoli lies in their use of in the denominator of equation 4 as they are describing a two-dimensional disc, whilst for three-dimensions we use . Making use of the approximation that gives our final equation


where has a value of 0.05 at  regardless of the chosen disc temperature profile. With this equation, taking (which is appropriate at ), and values from our unperturbed initial disc we are able to calculate the thermal diffusivities in models of differing opacity, . For the 100% interstellar grain opacity disc this gives , and the 1% opacity disc has a value of . These diffusivities correspond with characteristic disc cooling timescales of tens of orbital periods and less than an orbital period, respectively, where this timescale is linearly dependent upon the chosen diffusivity. These timescales indicate that in the higher diffusivity (lower opacity) calculations the imposed boundary will dominate the temperature structure of the disc, but this is not the case in the lower diffusivity (higher opacity) models.

It may also be of note that in these three-dimensional models, the dominant direction for radiative diffusion within the horseshoe region is vertically. For the high opacity calculations, vertical diffusion is found to be more than three times as efficient as radial diffusion.

4 Results

4.1 Is the surface treatment important?

First we address a comment made in Ayliffe & Bate (2010) in which we suggested that the modelling of migrating protoplanets using a surface treatment rather than an accreting point mass might significantly alter their migration rates, even changing the direction. This was inspired by a single result for a 10  protoplanet in a disc with a temperature profile exponent of , in which the migration was seen to change from outward to inward when a crude evacuating sink particle was replaced with a protoplanet with a surface in a high opacity disc (shown in the top panel of Fig. 2). To explore this further we have conducted a series of calculations spanning   in a disc with a steeper temperature profile, . It was envisaged that such a steep profile would lead to outward migration for a low mass protoplanet modelled with either of our treatments by increasing the significance of the coorbital torques, making their effects more apparent. The results of these calculations are shown in the lower panel of Fig. 2. It can be seen that using 1% interstellar grain opacities (plus symbols), even with this steep temperature profile, the protoplanet migration is predominantly inwards, with the 10  case being the only exception in exhibiting an almost stalled outwards migration. However, at 100% interstellar grain opacities those protoplanets with masses   are all migrating outwards, regardless of the chosen protoplanet treatment. Comparing the migration of these low mass protoplanets with those seen in the top panel of Fig. 2 illustrates the increased importance of the coorbital torques in the disc with a steeper temperature gradient.

Figure 2: The upper panel shows migration rates, which were presented as migration timescales in Ayliffe & Bate (2010), for protoplanets modelled using surfaces (symbols encircled, red) embedded in discs with . These discs were of either 100% (asterisks) or 1% (plus symbols) interstellar grain opacity. For a 333  protoplanet changing the opacity has such a small impact upon migration that the two rates appear indiscernibly overlaid above. The lower panel shows migration rates for protoplanets embedded in a circumstellar disc for which , some modelled with evacuating sinks (un-encircled symbols, black) and others using sinks with surfaces (once again, encircled symbols, red). All subsequent figures include only protoplanets modelled using sinks with surfaces.

In our new calculations (bottom panel) the 10  and 20  protoplanet migration rates obtained with an evacuating sink (un-encircled, black) are more rapidly outwards or less rapidly inwards than their counterparts modelled using sinks with surfaces (encircled, red); the reason for this is unclear, though it is in agreement with the sole case shown in the top panel for our previous calculations. However, at masses   this pattern reverses regardless of the disc opacity. Bate et al. (2003) showed that it was for such masses that a protoplanet might begin to form a non-negligible gap in the disc. The evacuation of a gap reduces the mass in the coorbital region and as a result reduces the magnitude of coorbital torques. In Ayliffe & Bate (2010) we found that between the two treatments of the protoplanet used here, it was the evacuating sink particles that developed clearer gaps as a result of drawing in material artificially fast relative to protoplanets modelled with surfaces. This suggests that protoplanets of   modelled with evacuating sinks migrate faster inwards (slower outwards) than there equivalents modelled with surfaces because they reduce the coorbital torques, which act outwards, to a greater extent. For a 100  protoplanet the coorbital region is sufficiently depleted with both treatments to lead to inward migration in a high opacity disc.

As to our findings in Ayliffe & Bate (2010) that changing the protoplanet treatment can reverse the direction of migration for a low mass protoplanet, we must clarify that this appears to be a result of the marginal nature of the migration of a 10  protoplanet in a high opacity disc with . The change in the forces acting when comparing a protoplanet modelled using an evacuating sink and one modelled using a surface, happened to be sufficient in this case to change the direction of migration, though more typically it is likely to give a different rate with the same direction (as seen in Fig. 2). The essential point is that as a protoplanet’s migration rate is found to be at least somewhat dependent upon the method by which its accretion is modelled, it is best that the adopted treatment be as realistic as possible. We therefore remain convinced that using a surface to allow gas-pile-up as a model of accretion is the best method for these purposes.

Figure 3: Migration rates of four different mass protoplanets; across and down the panels 10, 20, 33, and 50  (as marked). Each protoplanet was modelled in a range of circumstellar discs with different temperature profiles, , where the temperature profile exponent is marked on the x-axis. Models in which interstellar grain opacities are used (asterisks) lead to slower or reversed migration when compared with the 1% interstellar grain opacity models (plus symbols). Masset & Casoli (2010) give an expression for the total torque acting upon a protoplanet which includes the dependency upon the disc’s temperature profile and its thermal diffusivity, as well as the planet’s mass. The resulting migration rates are plotted in each panel, one for a diffusivity akin to our 100% opacity calculation (solid line), and one for the 1% opacity case (dot-dash line). Migration rates obtained from Tanaka et al. (2002) are independent of the temperature profile, and so appear in each panel as a horizontal dashed line, the rate of which is determined by the disc surface density and the protoplanet mass.

4.2 Impacts of the temperature profile

Our primary goal in this work was to explore the impact of the temperature profile in a circumstellar disc upon the migration of protoplanets therein. In particular we wished to test the predictions of recent analytic models, which by their nature are unable to include details of the gas flow near the protoplanet; Fig. 3 contains the results of our simulations to this end. Shown are the migration rates for protoplanets with masses 10, 22, 33, and 50 , embedded in discs with both high and low opacities, and a range of values of . A horizontal dashed line in each panel marks the migration rate obtained using Tanaka et al. (2002)’s analytic form, which depends on the disc surface density and protoplanet mass, but is independent of the value of . It can be seen that the low opacity models (1% interstellar grain opacity, marked with plus symbols) display some scatter, but generally centre around this analytic description. All of the low opacity calculations lead to inward migration, as would be expected for locally-isothermal calculations.

Also marked in Fig. 3 are analytic migration rates based on recent work by Masset & Casoli (2010). In their description there is a dependence on both and on the thermal diffusivity of the disc. As such we plot lines to correspond with both opacities used in calculating our results, where the diffusivities used are those calculated in Section 3.2; note that high opacities correspond to low thermal diffusivities and vice versa. In the case of the high opacity calculations (100% interstellar grain opacity, marked with asterisks) it can be seen that as the steepness of the temperature profile is increased, the inward migration rates slow towards zero before the direction changes to outwards and the rates begin to increase. In each panel of Fig. 3, excepting the 10  case for which the uncertainties are greatest, it is apparent that the rate of this change is in good agreement with the gradient given by the low thermal diffusivity (high opacity) result of Masset & Casoli’s analytic description (solid line). The analytic form also gives reasonable agreement (given the coarse spacing of in our results) as to the temperature profile at which a given mass protoplanet in a high opacity disc will transition from inwards to outwards migration. The formula of Masset & Casoli is also dependent upon the disc’s viscosity, which for plotting here we have used a Shakura & Sunyaev (1973) -value of 0.004, which is equivalent to the magnitude of the SPH viscosity in the vicinity of the protoplanet in our models. Note also that this analytic description is valid only up to masses for which the half-width of the horseshoe region, , is less than the disc thickness at . Adopting the adiabatic approximation of Baruteau & Masset (2008) for gives a maximum mass of 40 , meaning that the 50  cases shown in Fig. 3 are slightly beyond the range of validity.

Using Masset & Casoli’s analytic form with values that correspond to our low opacity models indicates a reversal of the migration trend with increasing , resulting in migration accelerating inwards for steeper profiles. However, the scatter in migration rates measured from our low opacity calculations is large, and prohibits us from concluding anything other that that the magnitudes show reasonable agreement. True locally-isothermal models were also conducted to establish whether they revealed a trend with the disc temperature profile. The results of these models are shown in Fig. 4, where once again there appears to be no definite trend in the change of migration rates with increasing , and the variation in the rates is small, varying by less than a factor of 2 from the minimum to maximum rate. As such we are unable to say whether or not the migration trend at low opacities given by Masset & Casoli (2010) is real.

Figure 4: Migration rates for a 33  protoplanet in true locally-isothermal discs with varying radial temperature profiles. There is a small range of migration rates, varying by less than a factor of two from the minimum to the maximum rate. No significant trend is evident.

Considering the immutable properties of our disc model, we find that only protoplanets   may migrate outwards for high values of . This can be seen in Fig. 5 where migration rates for a 100  protoplanet are plotted, and it is clear that no combination of opacity and is able to cause outward migration. Included in the figure are the same analytic descriptions used in Fig. 3, but here the Masset & Casoli (2010) line is well beyond the mass regime for which it is valid. For such a high mass protoplanet, the evacuation of the coorbital region should effectively end the significant influence of the coorbital torques. As such these effects are overestimated by the analytic model, which is consistent with the fact that our high opacity results all lie well below the solid line. Noticeably the analytic model, even with the overestimation of the coorbital torques which should favour outward migration, does not predict such migration for this high mass protoplanet (for ) , suggesting it is not possible in discs such as those modelled here.

Figure 5: Migration rates of a 100  protoplanet in a range of circumstellar discs with different temperature profiles, , where the temperature profile exponent is marked on the x-axis. Models in which interstellar grain opacities are used (asterisks) generally lead to very slightly slower inward migration than their equivalent 1% interstellar grain opacity (plus symbols) models. However, the trend with , which was quite pronounced for the high opacity models of lower mass protoplanets, is more subdued in this case. It also fails to follow the trend predicted using Masset & Casoli (2010) for a high opacity disc (solid line), though at 100  we are applying the analytic model well beyond its range of validity; the analytic trend for a low opacity disc is marked with a dot-dashed line. The change in disc opacity leads to a smaller change in the protoplanet migration rate for this relatively high mass protoplanet when compared with the lower mass protoplanets shown in Fig. 3. Rates from Tanaka et al. (2002) are shown as a horizontal dashed line.

4.3 Impact of absolute temperature

In this paper, we have studied the impact of the disc’s temperature profile on migration rates. In all these calculations, the absolute disc temperature at the planet’s orbital radius is the same. As discussed above we find reasonable agreement with the analytic description of Masset & Casoli. However, using Masset & Casoli’s analytic description it is possible to examine the migration rate’s dependence upon the absolute temperature at ; this is presented in Fig. 6. The migration rates, calculated assuming a disc with a diffusivity of (equivalent to our 100% opacity cases), are shown in the left hand panels. For each of the four protoplanet masses plotted ( 10, 20, 33, 50 ) there is a temperature (scaleheight) at which the outward migration rate achieves a maximum, with this temperature scaling as (within the mass range for which the model is valid,  ). The right hand panels of Fig. 6 show the same cases for a higher diffusivity of (equivalent to our 1% opacity cases), for which conditions there is no outward migration, regardless of the chosen absolute temperature of the disc, and there is no turnover in the rate of migration, though these rates still vary significantly with temperature.

These analytic results imply that as well as the temperature profile, the chosen initial disc temperature is an important factor in determining the migration behaviour of the various low mass planets. Temperatures in different discs around various types of star may be expected to exhibit temperatures throughout the range shown here of 10-230K at 5.2 au, suggesting further modelling should be carried out for a range of disc conditions to explore the utility of analytic descriptions like Masset & Casoli (2010), which if valid, can be applied more generally than the results of numerical modelling.

Figure 6: Migration rates of protoplanets in discs of differing absolute temperature, calculated using the analytic description of Masset & Casoli (2010). Rates calculated for conditions equivalent to our 100% opacity models are shown on the left hand side, and 1% opacity on the right. The temperature profiles are and in the top and bottom rows respectively. The four lines in each panel from left to right correspond to protoplanets with masses 10, 20, 33, and 50 , respectively. The vertical dashed line marks the temperature at used in our numerical models.

5 Discussion

First we address the low opacity calculations that we have performed. The resulting migration rates from these calculations appear scattered, showing no evident trend with changes in the disc temperature profile. This can be seen in Fig. 7 which gathers all the low opacity results shown in Fig. 3 into a single panel, as well as including the analytic rates for an isothermal disc calculated using Tanaka et al. (2002) (dashed lines) and rates based upon the recent work of D’Angelo & Lubow (2010) (solid lines) for a locally-isothermal disc. At these low opacities we expect the migration results from our models to be similar to those predicted with isothermal analytic expressions as we found in Ayliffe & Bate (2010). Indeed, despite the scatter, these migration results do reveal a trend with increasing mass that is matched by Tanaka et al. (2002)’s values. However, D’Angelo & Lubow have found that in the locally-isothermal (rather than globally isothermal) regime protoplanet migration rates depend upon the parent disc’s temperature profile, with steeper negative temperature profiles (higher values of ) leading to more rapid inward migration. This trend is not seen in our low opacity calculation results within which, as mentioned previously, it is not possible to establish a trend that is distinguishable from the numerical uncertainties; this is also the case in the locally-isothermal calculations shown in Fig. 4. However, the variation in migration rates obtained from D’Angelo & Lubow (2010) across the range of considered in this work is relatively small, less than a factor of 2 for each of the masses shown in Fig. 7. Our rather large uncertainties prevent us from definitively ruling out the model.

Figure 7: Migration rates of four different mass protoplanets through low opacity discs of differing temperature profiles. The masses are 10  (plus symbols), 20  (asterisks), 33  (diamonds), and 50  (triangles). Also included are analytic predictions. Those of Tanaka et al. (2002) without a dependence, for a three-dimensional isothermal disc, appear as horizontal dashed lines. The solid lines show the rates predicted using the dependent expression devised by D’Angelo & Lubow (2010) for a three-dimensional locally-isothermal disc. The masses that these lines correspond to are marked on the plot above the corresponding pairs.
Figure 8: Migration rates based on the total torques calculated using analytic expressions for a 33  protoplanet in a disc with properties equivalent to our models. The dashed line indicates the expected locally-isothermal result for such a protoplanet in a three-dimensional disc, obtained using D’Angelo & Lubow (2010). The solid lines are for rates obtained using Masset & Casoli (2010), where the colours progressing through the rainbow correspond to differing diffusivities, with the highest shown in red (labelled C). Three labelled thicker lines identify the trend with diffusivity, with the corresponding values given in the key. The locally-isothermal and highest diffusivity cases share the same gradient, with magnitudes differing by .

Returning to D’Angelo & Lubow (2010)’s analytic form for an isothermal disc, it is interesting to note that the acceleration of inward migration with increasing is in qualitative agreement with the description of Masset & Casoli (2010) discs with high thermal diffusivity. As a brief aside we compare these two analytic models with which we draw comparisons. Fig. 8 illustrates the migration rates calculated using D’Angelo & Lubow’s isothermal formula for a 33  protoplanet (dashed line), with the rates obtained using Masset & Casoli’s formula for a whole range of thermal diffusivities also shown. It can be seen that at the highest diffusivity (thick red line labelled C) the gradient of the change of migration rate with increasing for the two descriptions is extremely similar, and they are displaced in magnitude by a factor of , equivalent to , as was found in Paardekooper & Papaloizou (2008). This similar behaviour is to be expected, with a high thermal diffusivity leading to a disc that resembles a locally-isothermal model, and responds to different temperature profiles in the same manner. As seen earlier when considering our low opacity calculations, the manner in which the migration rate of a protoplanet changes with the discs thermal diffusivity is dependent upon the temperature profile of the disc. For discs with , as the thermal diffusivity is increased from its lowest value (thick purple line labelled A) the migration rate of a protoplanet can be expected to initially increase in the outward direction, or slow down if the migration is inwards. Beyond a certain diffusivity, shown by the thick blue line labelled B in Fig. 8, the trend is reversed with increasing diffusivity slowing outward migration, and accelerating inward migration. This behaviour is a result of the edge term of the Horsehoe drag in Masset & Casoli’s model. In future it would be desirable to explore a broader opacity space with greater refinement using our numerical models to determine whether we can capture the effects of the edge terms of the Horseshoe drag.

In our high opacity (low thermal diffusivity) models we find that increasing from 0 leads first to slower inward migration and then accelerating outward migration for protoplanets of  . This is in agreement with the analytic work discussed above. Similar outward migration has been found in the numerical models of Paardekooper & Papaloizou (2008), Baruteau & Masset (2008), Kley et al. (2009), and Paardekooper & Papaloizou (2009), and has been associated with a density enhancement ahead of the planet and a deficit behind it in the horsehoe region. In each of these models the feature has been clearly identified in two-dimensional modelling, but such a feature is likely to be much less apparent in three-dimensional models where gas is not constrained to move in a single plane. The lefthand panels of Fig. 9 present plots of averaged over 25 orbits for a disc of containing a 33  protoplanet. The leftmost panel is a 100% opacity calculation in which the protoplanet exhibits outward migration, whilst the neighbouring panel is for a 1% opacity calculation that leads to inward migration. Relative to the low opacity case there is a more substantial deficit of material trailing the protoplanet (negative ) in the horseshoe region of the high opacity case. The rightmost panel of Fig. 9 shows a ratio of the two surface density perturbation plots, and makes it clear that there is a greater density contrast between the leading and trailing positions in the horseshoe region of the high opacity outwardly migrating case. This density structure is reminiscent of that found in the two-dimensional models, though much less clearly defined.

Figure 9: The two left panels show the surface density normalised by the density gradient of the initial unperturbed disc. Respectively the panels show the 100% and 1% opacity models with and a protoplanet mass of 33 . The right hand panel shows the ratio of the two left panels, revealing how the structures differ between the optically thick and thin models. This makes it clear that relative to the low opacity (inward migrating) model the high opacity (outward migrating) model has a greater contrast between leading (positive ) and trailing (negative ) densities, thus explaining the differing directions of migration.

6 Summary

We have conducted a series of calculations to explore the impact of disc radial temperature profiles on the rates of protoplanet migration. The surface density profile remains unchanged in each calculation, thus each temperature profile gives rise to a different entropy profile through the disc. In agreement with analytic predictions, we find that steeper temperature profiles, , (and so entropy gradients) are much more disposed to bring about outwards migration for low mass ( ) protoplanets. There remains a dependence on the disc’s opacity, which determines how the disc and the embedded protoplanet are able to cool, which can change the discs structure and the resulting torques. In interstellar grain opacity (low thermal diffusivity) discs, low mass protoplanets are found to migrate outwards, whilst in discs with opacities of 1% this level we find no cases of outward migration. The dependence of migration rates upon the disc’s temperature profile in high opacity discs is in agreement with the predictions of recent analytic work by Masset & Casoli (2010).

For the low opacity calculations we find migration rates broadly inline with the rates predicted by Tanaka et al. (2002). These conditions can be represented using the formulae of Masset & Casoli (2010) with a suitably high thermal diffusivity. Such an approach suggests that a given protoplanet should migrate more rapidly inwards in a disc with a more steeply negative temperature profile. We found this to qualitatively agree with recent work by D’Angelo & Lubow (2010) who find a similar trend for a locally-isothermal disc. However, our efforts to numerically test these predictions were hampered by large numerical uncertainties in the results of our low opacity calculations, leaving us unable to determine the veracity of the analytic predictions.

We also expanded our investigation into the importance of the model used to represent the protoplanet. Using a surface that enables accreted gas to pile up and form an atmosphere results in somewhat different migration rates to those obtained when using a more crude accreting sink particle model. For the lowest-mass protoplanets ( ) the change in migration rate with protoplanet treatment was found to be small, though for borderline cases between inward and outward migration it could be sufficient to change the direction. Those protoplanets of   using the less realistic accreting sink model are found to migrate more rapidly inwards as they more effectively clear the coorbital region of material than equivalent surface models, more rapidly diminishing the positive coorbital torques. In general the effect of changing the protoplanet treatment is small compared to the differences brought about by altering the disc properties.


We thank the anonymous referee for prompting some further useful investigation. The calculations reported here were performed using the University of Exeter’s SGI Altix ICE 8200 supercomputer. Much of the analysis was conducted making use of SPLASH (Price, 2007), a visualisation tool for SPH that is publicly available at MRB is grateful for the support of a Philip Leverhulme Prize and a EURYI Award which also supported BAA. This work, conducted as part of the award “The formation of stars and planets: Radiation hydrodynamical and magnetohydrodynamical simulations” made under the European Heads of Research Councils and European Science Foundation EURYI (European Young Investigator) Awards scheme, was supported by funds from the Participating Organisations of EURYI and the EC Sixth Framework Programme.


  • Alexander (1975) Alexander D. R., 1975, ApJS, 29, 363
  • Alibert et al. (2004) Alibert Y., Mordasini C., Benz W., 2004, A&A, 417, L25
  • Ayliffe & Bate (2009) Ayliffe B. A., Bate M. R., 2009, MNRAS, 393, 49
  • Ayliffe & Bate (2010) Ayliffe B. A., Bate M. R., 2010, MNRAS, 408, 876
  • Baruteau & Masset (2008) Baruteau C., Masset F., 2008, ApJ, 672, 1054
  • Bate (1995) Bate M., 1995, PhD thesis, PhD thesis, Univ. Cambridge, (1995)
  • Bate et al. (1995) Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362
  • Bate et al. (2003) Bate M. R., Lubow S. H., Ogilvie G. I., Miller K. A., 2003, MNRAS, 341, 213
  • Benz (1990) Benz W., 1990, in Numerical Modelling of Nonlinear Stellar Pulsations Problems and Prospects, J. R. Buchler, ed., pp. 269–+
  • Benz et al. (1990) Benz W., Cameron A. G. W., Press W. H., Bowers R. L., 1990, ApJ, 348, 647
  • Boley et al. (2007) Boley A. C., Hartquist T. W., Durisen R. H., Michael S., 2007, ApJ, 656, L89
  • D’Angelo et al. (2005) D’Angelo G., Bate M. R., Lubow S. H., 2005, MNRAS, 358, 316
  • D’Angelo et al. (2002) D’Angelo G., Henning T., Kley W., 2002, A&A, 385, 647
  • D’Angelo et al. (2003) D’Angelo G., Kley W., Henning T., 2003, ApJ, 586, 540
  • D’Angelo & Lubow (2008) D’Angelo G., Lubow S. H., 2008, ApJ, 685, 560
  • D’Angelo & Lubow (2010) D’Angelo G., Lubow S. H., 2010, ApJ, 724, 730
  • Goldreich & Tremaine (1980) Goldreich P., Tremaine S., 1980, ApJ, 241, 425
  • Haisch et al. (2001) Haisch Jr. K. E., Lada E. A., Lada C. J., 2001, ApJ, 553, L153
  • Ida & Lin (2008) Ida S., Lin D. N. C., 2008, ApJ, 673, 487
  • Jang-Condell & Sasselov (2005) Jang-Condell H., Sasselov D. D., 2005, ApJ, 619, 1123
  • Klahr & Kley (2006) Klahr H., Kley W., 2006, A&A, 445, 747
  • Kley et al. (2009) Kley W., Bitsch B., Klahr H., 2009, A&A, 506, 971
  • Kley & Crida (2008) Kley W., Crida A., 2008, A&A, 487, L9
  • Korycansky & Pollack (1993) Korycansky D. G., Pollack J. B., 1993, Icarus, 102, 150
  • Li et al. (2009) Li H., Lubow S. H., Li S., Lin D. N. C., 2009, ApJ, 690, L52
  • Lubow et al. (1999) Lubow S. H., Seibert M., Artymowicz P., 1999, ApJ, 526, 1001
  • Masset (2002) Masset F. S., 2002, A&A, 387, 605
  • Masset & Casoli (2010) Masset F. S., Casoli J., 2010, ApJ, 723, 1393
  • Mayor & Queloz (1995) Mayor M., Queloz D., 1995, Nature, 378, 355
  • Menou & Goodman (2004) Menou K., Goodman J., 2004, ApJ, 606, 520
  • Monaghan (1992) Monaghan J. J., 1992, ARA&A, 30, 543
  • Monaghan (2002) Monaghan J. J., 2002, MNRAS, 335, 843
  • Monaghan & Gingold (1983) Monaghan J. J., Gingold R. A., 1983, Journal of Computational Physics, 52, 374
  • Mordasini et al. (2009) Mordasini C., Alibert Y., Benz W., 2009, A&A, 501, 1139
  • Morohoshi & Tanaka (2003) Morohoshi K., Tanaka H., 2003, MNRAS, 346, 915
  • Morris & Monaghan (1997) Morris J. P., Monaghan J. J., 1997, Journal of Computational Physics, 136, 41
  • Nelson et al. (2000) Nelson R. P., Papaloizou J. C. B., Masset F., Kley W., 2000, MNRAS, 318, 18
  • Paardekooper et al. (2010) Paardekooper S., Baruteau C., Crida A., Kley W., 2010, MNRAS, 401, 1950
  • Paardekooper et al. (2011) Paardekooper S., Baruteau C., Kley W., 2011, MNRAS, 410, 293
  • Paardekooper & Mellema (2006) Paardekooper S., Mellema G., 2006, A&A, 459, L17
  • Paardekooper & Mellema (2008) Paardekooper S., Mellema G., 2008, A&A, 478, 245
  • Paardekooper & Papaloizou (2008) Paardekooper S., Papaloizou J. C. B., 2008, A&A, 485, 877
  • Paardekooper & Papaloizou (2009) Paardekooper S., Papaloizou J. C. B., 2009, MNRAS, 394, 2283
  • Pollack et al. (1985) Pollack J. B., McKay C. P., Christofferson B. M., 1985, Icarus, 64, 471
  • Price (2007) Price D. J., 2007, Publications of the Astronomical Society of Australia, 24, 159
  • Price & Bate (2007) Price D. J., Bate M. R., 2007, MNRAS, 377, 77
  • Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  • Springel & Hernquist (2002) Springel V., Hernquist L., 2002, MNRAS, 333, 649
  • Tanaka et al. (2002) Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257
  • Ward (1986) Ward W. R., 1986, Icarus, 67, 164
  • Ward (1997) Ward W. R., 1997, Icarus, 126, 261
  • Whitehouse & Bate (2006) Whitehouse S. C., Bate M. R., 2006, MNRAS, 367, 32
  • Whitehouse et al. (2005) Whitehouse S. C., Bate M. R., Monaghan J. J., 2005, MNRAS, 364, 1367
  • Yamada & Inaba (2011) Yamada K., Inaba S., 2011, MNRAS, 411, 184
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description