Planet migration: self-gravitating radiation hydrodynamical models of protoplanets with surfaces
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 11 November 2018 (MN L A TEX style file v2.2)
Planet migration: self-gravitating radiation hydrodynamical modelsof protoplanets with surfaces
Ben A. Ayli ff e (cid:63) and Matthew R. Bate † School of Physics, University of Exeter, Stocker Road, Exeter EX4 4QL
11 November 2018
ABSTRACT
We calculate radial migration rates of protoplanets in laminar minimum mass solar neb-ula discs using three-dimensional self-gravitating radiation hydrodynamical (RHD) models.The protoplanets are free to migrate, whereupon their migration rates are measured. For lowmass protoplanets (10-50 M ⊕ ) we find increases in the migration timescales of up to an or-der of magnitude between locally-isothermal and RHD models. In the high-mass regime themigration rates are changed very little.These results are arrived at by calculating migration rates in locally-isothermal models,before sequentially introducing self-gravity, and radiative transfer, allowing us to isolate thee ff ects of the additional physics. We find that using a locally-isothermal equation of state,without self-gravity, we reproduce the migration rates obtained by previous analytic and nu-merical models. We explore the impact of di ff erent protoplanet models, and changes to theirassumed radii, upon migration. The introduction of self-gravity gives a slight reduction ofthe migration rates, whilst the inertial mass problem, which has been proposed for high massprotoplanets with circumplanetary discs, is reproduced. Upon introducing radiative transfer tomodels of low mass protoplanets ( ≈
10 M ⊕ ), modelled as small radius accreting point masses,we find outward migration with a rate of approximately twice the analytic inward rate. How-ever, when modelling such a protoplanet in a more realistic manner, with a surface whichenables the formation of a deep envelope, this outward migration is not seen. Key words: planets and satellites: formation – radiative transfer – methods: numerical –hydrodynamics – planetary systems: formation
The discovery of giant planets in tight orbits around their stars hasintroduced a complication that planet formation models must ad-dress. The formation theories contemporary to the discovery of 51Peg b (Mayor & Queloz 1995), the first such exoplanet discovered,had assumed planet systems with morphologies similar to the solarsystem, and were unable to explain the planet’s creation. The hightemperatures in the inner regions of circumstellar discs, as well asthe dominance of the star’s gravity, prevent formation by gravita-tional instability. Fairing no better, the core accretion model is con-founded by the low solids surface density in such high temperatureregions where no ice can form. However, a mechanism that mighto ff er an explanation for the existence of such Hot Jupiters had al-ready been suggested by Goldreich & Tremaine (1980) who hadstudied the rate of angular momentum and energy transfer betweena disc and an embedded planet. They found that the angular mo- (cid:63) E-mail: ayli ff [email protected] † E-mail: [email protected] mentum exchange was such that a planet’s orbital radius shouldchange significantly on short timescales ( ∼ − years).There are two primary types of planet migration which act indi ff erent mass regimes. Type I acts in the low mass regime ( (cid:46) ⊕ ), and is a result of a low mass planet’s interactions with an-nuli of gas that are in resonance with it, at so called Lindblad res-onances, as well as the corotation resonance. Angular momentumexchange occurs e ffi ciently with Lindblad resonances both insideand outside of the protoplanets orbital radius. The inner Lindbladresonance passes angular momentum to the protoplanet, acting toaccelerate it, whilst the outer Lindblad resonance works in oppos-tion. A linear analysis can be used to describe the e ff ect of Type Imigration, and in doing so it is found that the Lindblad resonancesare not spaced symmetrically inside and outside of the planet’sorbit. Instead each outer Lindblad resonance is somewhat closerthan the corresponding inner resonance, allowing the angular mo-mentum exchange with these outer resonances to dominate (Ward1986). The interaction of a planet with gas in the corotation region,the so called corotation torque, has a magnitude at most of only afew tens of percent of the di ff erential Lindblad torque (Ward 1997;Tanaka et al. 2002) when fully unsaturated. As a result, the net loss c (cid:13) a r X i v : . [ a s t r o - ph . E P ] J un B.A. Ayli ff e & M.R. Bate of angular momentum by the planet to the disc due to the di ff er-ential Lindblad torque causes the planet to migrate inwards (Ward1986).The rates of planet migration derived from analytic models(i.e. Ward 1997; Tanaka et al. 2002) are in agreement with thoseobtained using locally-isothermal numerical models (Korycansky& Pollack 1993; Nelson et al. 2000; Masset 2002; D’Angelo et al.2002; Bate et al. 2003; D’Angelo et al. 2003; Alibert et al. 2004;D’Angelo et al. 2005; Klahr & Kley 2006; D’Angelo & Lubow2008; Paardekooper & Mellema 2008a; Li et al. 2009). These ratesare such that low mass planets very rapidly fall into their stars, onmuch shorter timescales than are required for them to grow to of or-der a Jupiter mass. As such, to explain the observed exoplanet pop-ulation, some modification is required to slow down Type I migra-tion, something beyond what the locally-isothermal models haveconsidered (Alibert et al. 2004, 2005; Ida & Lin 2008; Mordasiniet al. 2009).One such modification is the inclusion of turbulence, dispens-ing with the usual laminar disc assumption. Three-dimensionalmagnetohydrodynamic models performed by Nelson & Papaloizou(2004) focused on the implications for migration of interactions ofprotoplanets within magnetorotationally unstable accretion discs. Itwas found that the disc turbulence led to large scale fluctuations inthe net torque experienced by a protoplanet over the course of itsorbit. As such the change in a protoplanet’s orbital radius becomesa random-walk. The most significant fluctuations were found forprotoplanets of <
30 M ⊕ , for which the time averaged fluctuationswere inconsistent with conventional Type I migration. The increasein a low mass protoplanet’s migration timescale due to the oscil-lation of its orbit suggests turbulence as a route by which to im-prove the survival probability of such bodies (Nelson & Papaloizou2004; Nelson 2005). Another addition to migration models whichcan lead to slower rates is more realistic thermodynamics. Sev-eral numerical models have introduced non-isothermal equations ofstate (Morohoshi & Tanaka 2003; Jang-Condell & Sasselov 2005;Paardekooper & Mellema 2006, 2008b; Baruteau & Masset 2008a;Paardekooper & Papaloizou 2008b; Kley & Crida 2008; Kley et al.2009), and have found that the rate of Type I migration can be re-duced, or its direction reversed.The second type of migration, called Type II, becomes e ff ec-tive once a planet grows to a mass of a few tenths of a Jupiter massand is able to clear its corotation region of material, forming a gapin its parent disc (Ward 1997; Bryden et al. 1999; Kley 1999). Thewidth of the gap is determined by the competing influences of theviscous and pressure forces, which tend to close the gap, and thegravitational torques which work to open it. The planet becomeslocked in the gap it has created, and migrates inwards on the disc’sviscous accretion timescale. Such Type II migration continues onthe viscous timescale whilst the protoplanet mass is less than orcomparable to the mass of the local disc with which it is interacting.A protoplanet that comes to dominate the local mass can no longerbe treated simply as a constituent part of the disc. Rather, its inertiaacts to slow the rate of migration that the disc’s viscous evolutioncan bring about (Syer & Clarke 1995; Ivanov et al. 1999). Ivanovet al. suggested an expression to describe this late stage migrationwhich led to a reduction of the migration rate with increasing pro-toplanet mass and with reduced orbital radius. This gives a slowingmigration that might enable massive protoplanets to avoid fallinginto the central star, and help to form Hot Jupiters. The transitionfrom Type I into Type II migration cannot be treated analyticallydue to the many non-linear processes at work, but is has been mod- elled numerically in both two and three-dimensions (Bryden et al.1999; Kley 1999; D’Angelo et al. 2002; Bate et al. 2003).This paper discusses numerical models conducted to investi-gate planet migration, as has been done previously by numerousauthors in both two and three dimensions (Korycansky & Pollack1993; Nelson et al. 2000; Masset 2002; D’Angelo et al. 2002; Bateet al. 2003; Nelson & Papaloizou 2003; D’Angelo et al. 2003; Alib-ert et al. 2004; Nelson & Papaloizou 2004; Nelson 2005; D’Angeloet al. 2005; Masset et al. 2006; Klahr & Kley 2006; Crida & Mor-bidelli 2007; Paardekooper & Papaloizou 2008b; Paardekooper &Mellema 2008a; D’Angelo & Lubow 2008; Kley & Crida 2008;Pepli´nski et al. 2008a,b; Li et al. 2009; Paardekooper & Papaloizou2009; Crida et al. 2009; Kley et al. 2009; Yu et al. 2010; Baruteau& Lin 2010). Some of the more recent of these works have begun toinclude more physics, such as magnetic fields, and non-isothermalequations of state as discussed above. Other models, conducted byNelson et al. (2000), D’Angelo et al. (2003), Alibert et al. (2004),and Alibert et al. (2005), have included mass and angular momen-tum accretion onto migrating protoplanets, whilst the influence ofself-gravity was explored by Nelson & Benz (2003), Baruteau &Masset (2008b), and analytically by Pierens & Hur´e (2005). Re-gardless of the physics included, it has been common to excludetorques from the immediate vicinity of the planet due to the di ffi -culty of resolving this region in any great detail. This poor resolu-tion also prevents the self-consistent modelling of a protoplanetaryenvelope, which requires both good resolution and the inclusion ofself-gravity and radiative transfer.In this work we conduct three dimensional global disc mod-els, using smoothed particle hydrodynamics (SPH), that include ra-diative transfer, self-gravity, and which use a planetary surface toallow modelling of gas flow to well within the Hill sphere. Withour surface treatment, and the adaptable resolution of SPH, the in-nermost torques can be well resolved and a realistic protoplanetaryenvelope can develop. We conduct models of varying complexity,introducing more physics in stages such that we can determine theimpact of each addition on protoplanet migration. The first modelswe discuss are designed to emulate the treatment used in the gridmodels of Bate et al. (2003), allowing us to test our SPH calcula-tions. We then introduce mass and angular momentum accretion,followed by the addition of self-gravity. Next the surface treatmentis employed to explore the impact of realistic gas flow within theHill sphere. Finally radiative transfer is folded into the mix, and arange of opacities explored.In Section 2, we describe our computational method and itstesting. In Section 3 we present the results of our models, startingwith the simplest locally-isothermal case, and progressing to self-gravitating radiation hydrodynamic cases using our planet surfacetreatment. Section 4 discusses the implications of our results for gi-ant planet formation theory, whilst a summary and our conclusionsare given in Section 5. The calculations described herein have been performed using athree-dimensional SPH code. This SPH code has its origins in a ver-sion first developed by Benz (1990; Benz et al. 1990) but it has un-dergone substantial modification in subsequent years. Energy andentropy are conserved to timestepping accuracy by use of the vari-able smoothing length formalism of Springel & Hernquist (2002)and Monaghan (2002) with our specific implementation being de-scribed in Price & Bate (2007). Gravitational forces are calculated c (cid:13) , 000–000 lanetary migration and growth and neighbouring particles are found using a binary tree. Radiativetransfer is modelled in the two temperature (gas, T g , and radiation, T r ) flux-limited di ff usion approximation using the method devel-oped by Whitehouse et al. (2005) and Whitehouse & Bate (2006).Integration of the SPH equations is achieved using a second-orderRunge-Kutta-Fehlberg integrator with particles having individualtimesteps (Bate 1995). The code has been parallelised by M. Bateusing OpenMP. We model a protoplanetary disc (2 π radians) with radial bounds of0.1 - 3 r p (0.52 - 15.6 AU), where r p is the initial orbital radius ofour embedded protoplanets, taking a value of 5.2 AU. The radialtemperature profile for the disc is initialised as T g ∝ r − , whilstthe surface density of the disc goes as Σ ∝ r − / ; these profilesare equivalent to those of Bate et al. (2003), giving a constant ra-tio of disc scaleheight to radius of H / R ≈ .
05. At r p the initialtemperature is ≈
75K and the surface density is 75 g cm − . Thisleads to a total disc mass of ≈ .
005 M (cid:12) within the boundaries, (or ≈ . (cid:12) within 1.56 - 20.8 AU inline with Bate et al. 2003, or0.01 M (cid:12) within 26 AU inline with D’Angelo et al. 2003). The discbegins with a Keplerian velocity structure, established around a 1M (cid:12) star. It should be noted that when modelled with radiation hy-drodynamics the disc settles to have a slightly steeper temperatureprofile at the midplane, where in our disc it becomes approximately T g ∝ r − . . A similar steepening is found by Kley et al. (2009) intheir radiation hydrodynamical models.The inner edge of the disc is bounded by a potential barrierwhich acts to prevent the disc accreting on to the central star. Theouter edge of the disc is surrounded by a population of ghost parti-cles that act as a pressure barrier to prevent the disc from spreadingdue to shear viscosity. Particles that do move beyond the outer edgeof the computational domain are removed from the calculation. Themigration rates of planets are determined by the structure of the discin their vicinity, particularly the surface density profile. As such,the aim of these boundaries was to establish a broad region of discabout r p throughout which the surface density followed the desiredprofile, with the only modification being due to the planet-disc in-teraction. Maintaining a broad stable region is important to allowany embedded protoplanet to migrate away from its initial positionat r p , and to allow density waves to propagate to large distances,such that their interactions with the protoplanet are included to thepoint of inconsequence.The inner boundary allowed gas to pile up at a small radius,far from the intended initial orbital radius of any introduced pro-toplanets. Once gas has accumulated at the boundary, the concen-tration of material exerts an outward pressure gradient back intothe disc, supporting further gas against inward migration. Fig. 1 il-lustrates the disc profile that develops, and illustrates the successof the boundaries in preserving the desired surface density profile(shown by the long-dashed line) throughout the region of likelymigration, from 0.5 to 2 r p (2.6 - 10.4 AU). The profiles shownare for calculations with 5 × , 2 × and 10 particles. In thehighest resolution calculation, the rarefaction due to the outwardspressure from the boundary is narrowest and shallowest. However,the radial range over which a consistent surface density profile isachieved is not greatly changed between the two higher resolutionmodels. For calculations performed with 5 × particles this is notthe case, with a much greater encroachment of the boundary e ff ectsinto the likely migration region. In the calculations that follow weuse 2 × particles as this was su ffi cient to establish a stable disc Figure 1.
Surface density profiles of a locally-isothermal disc modelledwith 10 (solid line), 2 × (short-dashed line), and 5 × (long-dashedline) particles, in the absence of a planet, after 4 orbits of the disc’s outer-edge. The straight dashed line across the top denotes the desired surfacedensity profile of Σ ( r ) ∝ r − , with a value of 75 g cm − at r = r p . With theimplemented disc boundaries this profile is achieved by the two higher res-olution calculations in the region over which a planet may migrate, whilstin the lowest resolution calculation the e ff ect of the inner boundary reachesto almost r p . in the planet’s vicinity, whilst being considerably faster than using10 particles, for which the benefits were marginal.Upon creating a disc model, there is generally a period of tran-sient structural perturbations as the model settles to a stable state.This settling must be completed before protoplanets are introducedin order that their orbits are not immediately disturbed by perturba-tions unrelated to the interesting protoplanet-disc interactions. Theinitial discs were evolved in the absence of a planet until any tran-sience resulting from settling had dissipated, which required justover 4 orbits of the disc’s outer edge. Two equations of state are used for the calculations presented inthis paper. One is a locally-isothermal equation of state with thetemperature of the gas throughout the disc being a fixed function ofradius (as described in the disc setup). The second is used for ourradiation hydrodynamical calculations, which are conducted usingan ideal gas equation of state p = ρ T g R g /µ where R g is the gas con-stant, ρ is the density, T g is the gas temperature, and µ is the meanmolecular mass. The equation of state takes into account the trans-lational, rotational, and vibrational degrees of freedom of molecu-lar hydrogen (assuming a 3:1 mix of ortho- and para-hydrogen; seeBoley et al. 2007). It also includes the dissociation of molecular hy-drogen, and the ionisations of hydrogen and helium. The hydrogenand helium mass fractions are X = .
70 and Y = .
28, respec-tively, whilst the contribution of metals to the equation of state isneglected. More details on the implementation of the equation ofstate can be found in Whitehouse & Bate (2006). c (cid:13) , 000–000 B.A. Ayli ff e & M.R. Bate
The two temperature (gas and radiation) radiative transfer inthe flux-limited di ff usion approximation employed in this work isimplemented as described by Whitehouse et al. (2005) and White-house & Bate (2006). Briefly, work and artificial viscosity (whichincludes both bulk and shear components) increase the thermal en-ergy of the gas, and work done on the radiation field increases theradiative energy which can be transported via flux-limited di ff u-sion. The energy transfer between the gas and radiation fields isdependent upon their relative temperatures, the gas density, and theopacity, κ . We use interpolation of the opacity tables of Pollack et al. (1985) toprovide the interstellar grain opacities for solar metallicity molec-ular gas, whilst at higher temperatures when the grains have subli-mated we use the tables of Alexander (1975) (the IVa King model)to provide the gas opacities (for further details see Whitehouse& Bate 2006). The grain opacities are based on interstellar abun-dances, and may not hold true in a circumstellar disc where ag-glomeration can strongly modify the grain population, reducingthe opacity. Therefore, we vary the grain opacity in our models tomimic the e ff ects of these processes; the gas opacities of Alexander(1975) that are used at higher temperatures are not scaled. In thispaper we use 100, 10, and 1 per cent interstellar grain opacities (seeAyli ff e & Bate 2009a for further details). Many of the calculations performed here involved radiative trans-fer, and so required a mechanism by which energy could be lostfrom the disc into the encompassing vacuum. This necessity arisesbecause the flux-limited di ff usion scheme transfers energy betweenSPH particles, rendering it unable to radiate into a vacuum wherethere are no particles. By determining the height above and belowthe midplane at which the optical depth into the disc, τ op ≈
1, wecan identify two layers of particles which exist in optically thin re-gions. Compelling particles that fall within these regions to followthe initial temperature profile of the disc allows the energy they re-ceive from the bulk of the disc to be e ff ectively lost, as though intothe vacuum. This initial temperature profile was chosen for pur-poses of comparison with previous work where it is common, andcrudely represents the insolation of the disc.The vertical density structure of the disc approximately fol-lows a Gaussian distribution, integration over which yields the sur-face density, Σ r , in terms of the error function, erf(). The opticaldepth is given by τ op ( a ) = − (cid:82) ∞ a κρ dz . Completing this integrationwith the Gaussian form of ρ , setting τ op =
1, and the upper limit ofintegration set to infinity to represent the vacuum, it is possible torearrange to find the lower limit of integration, depth ‘ a ’, as shownin equation 1. aH = √ − (cid:32) − τ op κ Σ r (cid:33) . (1)This is the depth in the disc from above which radiation can escapeto the vacuum, and so defines the height of our radiation boundary.The boundary height varies with radius, due to the radial depen-dence of Σ r and κ , the values of which are taken from the initialdisc distribution. Once a suitable disc has been established, it is possible to embed aprotoplanet and model its evolution. In all of the models discussedin this paper, the protoplanet is free to move from its initial orbitalradius, and its migration rate is measured, rather than calculatedfrom static torques. How the protoplanet interacts with the disc,and so migrates, depends upon the way in which the protoplanetis modelled, and we employ several methods in this work. In thefirst case gas is not accreted, but instead gas that comes withina specified radius, r acc , of the point mass representing the proto-planet, is removed from the calculation, and its energy, momentum,and mass are all dispensed with. We call the particles used for thistype of protoplanet treatment Killing sink particles, and they mi-grate purely due to torque interactions with the disc, approximatelyrecreating the treatment used in Bate et al. (2003). Typically we set r acc to be some fraction of the Hill radius, R H .The second treatment also models the protoplanet as a sinkparticle which removes gas that comes within r acc , but in this casethe gas properties are added to the sink particle. These sink parti-cles are called Accreting sinks. They are able to absorb the energy,mass, linear momentum, and angular momentum of the accretedgas (see Bate et al. 1995).The third and final protoplanet treatment is one that we haveused in previous work (Ayli ff e & Bate 2009a,b) which uses a sur-face force that allows gas to pile up on the surface of a protoplanet,realistically modelling accretion. This scheme is now applied bywrapping the surface around a gravitating point mass that is freeto move, and which represents the protoplanet’s initial mass. Sucha protoplanet can move in response to changes in the momentumof the point mass and its bound envelope, the extent of which isself-consistently arrived at, rather than being prescribed. The sur-face is applied as a modification to the usual gravitational force inthe protoplanet’s immediate vicinity, taking the form F r = − GM p r − (cid:32) R p − rR p (cid:33) , (2)for r < R p , where R p is the protoplanet radius, r is the radiusfrom the centre of the protoplanet, and M p is the protoplanet’s ini-tial mass. This equation yields zero net force between a gas parti-cle and the gravitating point mass at the surface radius R p . Insideof this radius the force is outwards and increases rapidly with de-creasing radius, such that the inner most gas particles come to restvery close to the surface, with an equilibrium position slightly in-side due to the envelope weight pressing down. As in our previouswork, we provide a smooth start to the calculations by embedding aprotoplanet of radius R p = . r p which then shrinks exponentiallyto the desired radius during the first orbit of the protoplanet. Resolution was considered previously in relation to the disc bound-aries, and establishing the desired density profile in an undisturbeddisc. However, ensuring that the planet-disc interaction is well re-solved in our models, and that the migration rates obtained are in-dependent of resolution, required further testing. To this end, twosimulations were performed which were identical in every respectexcepting their resolution. A 10 M ⊕ protoplanet, modelled with anAccreting sink particle, was used for these calculations. The spiralwaves launched by such a low mass planet are weak, as are the re-sulting torques acting back upon the planet. So whilst a degree ofmigration is expected, insu ffi cient resolution might poorly resolve c (cid:13) , 000–000 lanetary migration and growth Figure 2.
The orbital evolution of a 10 M ⊕ protoplanet modelled usingAccreting sink particles, with accretion radii of ∼ R H (a radius chosento mimic the ZEUS models of Bate et al. 2003). Two resolutions were usedto perform otherwise identical calculations, 10 particles (thick upper line)and 2 × particles (thin lower line). The evolution covers 50 orbits of theprotoplanet using a locally-isothermal equation of state. The model at bothresolutions establishes the same rate of migration after a period of settling,roughly equal to the libration time. these small e ff ects; this makes such a planet an ideal candidate fortesting resolution. Fig. 2 shows the evolution of the planet’s orbitalradius with time for two resolutions.There is an initial period of settling whilst the torques due tothe horseshoe region are established. As was seen in Masset (2002),the outward corotation torques almost balance the inward di ff eren-tial Lindblad torques when the protoplanet is first embedded lead-ing to almost no migration. However, as the model evolves thecorotation torques become partially saturated, reducing their mag-nitude and allowing the Lindblad torques to dominate. As a resultthe protoplanet’s migration is seen to accelerate until the corota-tion torques reach a quasi-equilibrium state, which takes of order alibration time (Kley et al. 2009), given by τ lib = r p x s P p , (3)where x s is the width of the horseshoe region, and P p is the planet’sorbital period (11.8 years in the calculations presented here). Fol-lowing Kley et al., by using an approximation for the width ofthe horseshoe region, taken from the locally-isothermal models ofMasset et al. (2006), it is possible to get an idea of the horseshoesettling time. For a 10 M ⊕ protoplanet τ lib ∼
50 orbits, whilst fora 333 M ⊕ protoplanet τ lib ∼
8. It appears from Fig. 2 that a clearmigration trend is established in around 30 orbits. Once Type I mi-gration is established, the migration rate is consistent between thetwo di ff erent resolution calculations. It is not clear why the initialsettling appears to be resolution dependent, but the agreement ofthe eventual rates gives us confidence in the migration rates ob-tained using 2 × particles. Viscosity is the means by which angular momentum is transportedthrough a laminar circumstellar disc, enabling material to flow intowards the star. The viscosity envisaged is not caused by a molec-ular interaction, which is insu ffi cient to explain the rates of accre-tion inferred from observations. The primary alternative is thoughtto be a magneto-rotational instability (Balbus & Hawley 1991).Within the models presented here the viscosity is not achievedby modelling the responsible processes, but instead by using a pa-rameterised form of viscosity, though somewhat di ff erent to a typ-ical Shakura-Sunyaev α -viscosity (Shakura & Sunyaev 1973). Anartificial viscosity for SPH, designed to conserve linear and angu-lar momentum, was introduced by Monaghan & Gingold (1983),and modified by Monaghan (1992) to deal with high mach num-ber shocks. This artificial viscosity is not designed to represent realviscosity, but instead to allow the modelling of shock phenomena,and to damp numerical noise. The α and β terms set the strengthof the viscosity, with typical values of α = β =
2. The twoterms serve di ff erent functions. The α -term establishes a viscosityto damp subsonic velocity oscillations that may be produced in thewake of a shock front. The β -term prevents particle interpenetrationin supersonic shocks.The viscosity is only applied when the gas is under compres-sion, ideally near shocks, and we use the switch developed by Mor-ris & Monaghan (1997) to try and reduce the action of artificialviscosity where the cause is not a shock. Using this switch meansthat instead of the global α being applied to all the particles, theyeach scale this value (down to a minimum value of α min ) based ontheir proximity to a shock. Setting α min = . ff erential rate of rotation of thedisc leads to a shear viscosity by which angular momentum is ex-changed through the disc. The value of the viscosity at any pointcan be stated as a Shakura-Sunyaev α -viscosity, α SS , using α SS (cid:39) . α hH (4)as given by Lodato & Price (2010), where the coe ffi cient given hereis a factor 2 smaller because the viscous force is only applied be-tween approaching, not receding, particles. Here h is the local aver-age SPH smoothing length, and H is the disc scaleheight. Around r p these quantities have values of approximately 0.016, and 0.05respectively, with ¯ α ≈ .
25 giving α SS ≈ . The migration rates were calculated by taking linear fits of the or-bital radii evolution data (e.g. Fig. 2). Each simulation was allowedto evolve for at least 50 orbital periods, (as was done in Bate et al.2003), of which only the latter 25 were used for fitting purposes.The early orbits are sacrificed to ensure that any transient disrup-tion to the disc caused by the sudden introduction of a planet havesubsided and to allow gas in the horseshoe region to reach a quasi-equilibrium state, as discussed previously. Fig. 3 illustrates the or-bital migration of 4 di ff erent mass protoplanets modelled with Ac-creting sinks and a locally-isothermal equation of state. It can beseen that the final 25 orbits of evolution are settled, and so providea suitable period over which to measure the migration rates. For thehighest mass protoplanets it can be seen that their migration overthe last 25 orbits is linear in time, suggesting that the rate is well c (cid:13) , 000–000 B.A. Ayli ff e & M.R. Bate
Figure 3.
The orbital evolution of 4 protoplanets modelled in locally-isothermal discs, using Accreting sink particles with r acc = . R H . In allinstances the planetary migration is well underway and free from large scalefluctuations after 25 orbits. The migration rates are measured over the last25 orbits. The Jupiter mass protoplanet appears to settle in ≈
10 orbits,inline with the libration time of ≈ established. Deviations from linearity, which might cause a longterm change in the rate are very subtle, and would only be estab-lished on very long timescales, inaccessible to our models due totheir computational expense.There are instances where measuring the migration rate overthe last 25 orbits will lead to faster rates for a stated mass thanshould truly be ascribed to it. In the first 25 orbits, in the Ac-creting sink and point mass with surface calculations, the proto-planet’s mass has increased through accretion. As a result of this,the migration rate tends to accelerate, giving a non-linear orbitalradius evolution. It is therefore the case that linear fits to the latterhalf of this evolution don’t exactly render the gradient expected forthe initial protoplanetary mass. The impact of such mass accretionmost strongly influences the orbital evolution of the lowest masscores, which accrete a higher fraction of their total mass. As a re-sult, we present the migration rates with uncertainties in the mass,portrayed by error bars which stretch across the mass range overwhich the linear fit is taken, with the migration timescale plotted atthe mean mass. In the Accreting sink calculations, the mass of theprotoplanet includes that of all the gas accreted, but for the pointmass with surface calculations the accreted mass is measured asthe mass within the Hill radius. For a majority of the calculations,the mass spread is very small, and thus the mass associated witha given migration rate is well known. The migration timescale, τ ,is then calculated as the time a planet of fixed mass would take tomigrate from its starting orbital radius, (5.2 AU in all of the modelsdiscussed), in to the central star; τ = r p / ˙ r . The rate of migrationis assumed constant in calculating these timescales, allowing readycomparison with previous works. Figure 4.
A reproduction of figure 10 from Bate et al. (2003) showing themigration rates they obtained by considering the torques at work upon pro-toplanets of fixed orbital radius in their models. The dots (with solid errorbars) denote the rates calculated when the torques outside of R H were in-cluded, the lower error bars are the rates when including torques externalto 0.5 R H , and the upper error bars are the rates for just the torques beyond1.5 R H . The analytic Type I model of Tanaka et al. (2002) and the TypeI / II model of Ward (1997) are plotted as short-dashed and long-dashed linesrespectively. The triangles (with dotted error bars) mark the measured mi-gration timescales from the SPH models in which the protoplanet is able tomigrate; the Killing sink particle accretion radii are R H . The lower and up-per error bars for the SPH cases mark the measured rates with Killing sinkparticles of 0.5 and 1.5 R H accretion radii respectively. The SPH modelsshow good agreement with the grid-based ZEUS models. We began by demonstrating that out SPH code can reproduce themigration rates obtained from past studies. These comparisons arebetween the migration rates obtained using our Killing sinks inlocally-isothermal discs, and those obtained in ZEUS finite dif-ference models by Bate et al. (2003), who used the same initialconditions. The Killing sinks are used to approximately imitate theplanet treatment of Bate et al.. We expect our SPH models to re-produce similar migration timescales to these previous grid basedresults, however, Bate et al. calculate their migration timescales byexamining the net torque due to gas outside the planet’s Hill radiusupon the planet, which was on a fixed orbit at 5.2 AU. In makingthis comparison, the implicit assumption is that the torques do notchange significantly when the planet is allowed to migrate ratherthan being held on a fixed orbit. Indeed, over the 50 orbits that aremodelled in the SPH calculations, the planets do not migrate anyconsiderable distance, as can be seen in Fig. 3. A reproduction ofFig. 10 from Bate et al. (2003) is shown in Fig. 4, where the ZEUSresults are plotted with circles and solid error bars, whilst the SPHderived migration timescales are plotted using triangles with dottederror bars. The term error bar is a misnomer, instead the extremesof these bars indicate the migration timescales calculated by in-cluding more (0.5 R H , lower bound) or less (1.5 R H , upper bound) c (cid:13) , 000–000 lanetary migration and growth Figure 5.
A comparison of the protoplanet migration timescales measured in locally-isothermal, non-self-gravitating models using Killing sinks (asterisks),that exclude mass and angular momentum accretion, and otherwise identical calculations using Accreting sink particles (plus symbols) which include suchaccretion. The left panel is for those calculations where r acc = . R H , whilst in the right panel r acc = . R H . The analytic Type I model of Tanaka et al. (2002)is plotted as a dashed line. Note that the Accreting sinks begin with the same four masses as the Killing sinks, but their masses increase substantially duringthe calculations, particularly the low-mass protoplanets. The mass error bars for the Accreting sink particles show the range in mass over which the linear fitto determine the migration rate was applied. Including mass and angular momentum accretion consistently leads to faster rates of migration when measured atequivalent masses. of the gas closest to the protoplanet in the torque calculations. Inthe SPH models these same limits are found by conducting simu-lations with killing radii of r acc = . R H , 1 R H , and 1.5 R H ; notethat this is not exactly equivalent to changing the torque exclusionradius, and the limitations of the comparison will be discussed inthe following paragraphs. For all protoplanet masses the migrationtimescales obtained with the SPH code, using killing radii of R H ,di ff er by less than 35 per cent from the ZEUS based results. TheSPH results also show the expected transition from the linear TypeI migration regime to the non-linear Type II beyond ∼ . J . As aresult, we are confident that our SPH code is capable of modellingmigration proficiently.The spread in the migration timescales we find for the variouskilling radii are larger that those found in the ZEUS calculationswith di ff erent torque exclusion radii. Explanations for these largerspreads can be found in the density and resolutions of the two mod-els. For a given mass protoplanet, there is only one ZEUS modelupon which the three di ff erent torque exclusion radius calculationsare based. This gives a consistent density structure between the dif-ferent cases, but a fixed spatial resolution which is poor near to 0.5 R H . In the SPH calculations, the r acc equivalent to these exclusionradii are each probed in separate models. The evacuated region foran r acc = . R H Killing sink is far larger than for the 0.5 R H case,which changes the circumplanetary density structure between thetwo cases, with lower densities at equivalent radii for the larger r acc .This leads to lower torques acting from the boundary of a Killingsink with r acc = . R H than at an equivalent exclusion radius in theZEUS models; hence the larger migration timescale.For r acc = . R H the SPH resolution is far better than that of the ZEUS calculations, and so the models better discern the gasthat is available to exert torques, leading to faster migration. At allthe radii, regardless of the di ff erent spreads, there is consistently aturn away from Type I to Type II migration at higher masses, andthe absolute migration timescales are within a factor of two or bet-ter at each mass and radius. The figure makes clear that in both theSPH and ZEUS models, the migration timescales reduce when thekilling radius (or torque exclusion radius) is shrunk. This indicatesthat the contributions from torques from within R H are negative, aswas found in Bate et al. (2003). This is an important agreement asthe region within R H is not treated in analytic descriptions of mi-gration, and so is not well understood. It appears that the persistenttorques from within R H are negative, which is in the opposite senseto the corotation torques in a disc with the surface density profilemodelled here. We return to this in Section 3.4. Whilst Killing sinks provided a means by which to test and com-pare our models, they do not allow us to include protoplanetgrowth. To that end we next introduced accretion into the migrationmodels, using our Accreting sinks. Fig. 5 compares the migrationtimescale found for Accreting sink particles in locally-isothermalmodels with those obtained using the Killing sink method pre-viously described. Whilst the Killing sinks maintain their initialmasses throughout, the Accreting sink models lead to increasinglymassive protoplanets, sometimes growing to several times their ini-tial mass. This makes interpreting the relative migration rates some-what di ffi cult, though the broad result is that including accretion c (cid:13)000
A comparison of the protoplanet migration timescales measured in locally-isothermal, non-self-gravitating models using Killing sinks (asterisks),that exclude mass and angular momentum accretion, and otherwise identical calculations using Accreting sink particles (plus symbols) which include suchaccretion. The left panel is for those calculations where r acc = . R H , whilst in the right panel r acc = . R H . The analytic Type I model of Tanaka et al. (2002)is plotted as a dashed line. Note that the Accreting sinks begin with the same four masses as the Killing sinks, but their masses increase substantially duringthe calculations, particularly the low-mass protoplanets. The mass error bars for the Accreting sink particles show the range in mass over which the linear fitto determine the migration rate was applied. Including mass and angular momentum accretion consistently leads to faster rates of migration when measured atequivalent masses. of the gas closest to the protoplanet in the torque calculations. Inthe SPH models these same limits are found by conducting simu-lations with killing radii of r acc = . R H , 1 R H , and 1.5 R H ; notethat this is not exactly equivalent to changing the torque exclusionradius, and the limitations of the comparison will be discussed inthe following paragraphs. For all protoplanet masses the migrationtimescales obtained with the SPH code, using killing radii of R H ,di ff er by less than 35 per cent from the ZEUS based results. TheSPH results also show the expected transition from the linear TypeI migration regime to the non-linear Type II beyond ∼ . J . As aresult, we are confident that our SPH code is capable of modellingmigration proficiently.The spread in the migration timescales we find for the variouskilling radii are larger that those found in the ZEUS calculationswith di ff erent torque exclusion radii. Explanations for these largerspreads can be found in the density and resolutions of the two mod-els. For a given mass protoplanet, there is only one ZEUS modelupon which the three di ff erent torque exclusion radius calculationsare based. This gives a consistent density structure between the dif-ferent cases, but a fixed spatial resolution which is poor near to 0.5 R H . In the SPH calculations, the r acc equivalent to these exclusionradii are each probed in separate models. The evacuated region foran r acc = . R H Killing sink is far larger than for the 0.5 R H case,which changes the circumplanetary density structure between thetwo cases, with lower densities at equivalent radii for the larger r acc .This leads to lower torques acting from the boundary of a Killingsink with r acc = . R H than at an equivalent exclusion radius in theZEUS models; hence the larger migration timescale.For r acc = . R H the SPH resolution is far better than that of the ZEUS calculations, and so the models better discern the gasthat is available to exert torques, leading to faster migration. At allthe radii, regardless of the di ff erent spreads, there is consistently aturn away from Type I to Type II migration at higher masses, andthe absolute migration timescales are within a factor of two or bet-ter at each mass and radius. The figure makes clear that in both theSPH and ZEUS models, the migration timescales reduce when thekilling radius (or torque exclusion radius) is shrunk. This indicatesthat the contributions from torques from within R H are negative, aswas found in Bate et al. (2003). This is an important agreement asthe region within R H is not treated in analytic descriptions of mi-gration, and so is not well understood. It appears that the persistenttorques from within R H are negative, which is in the opposite senseto the corotation torques in a disc with the surface density profilemodelled here. We return to this in Section 3.4. Whilst Killing sinks provided a means by which to test and com-pare our models, they do not allow us to include protoplanetgrowth. To that end we next introduced accretion into the migrationmodels, using our Accreting sinks. Fig. 5 compares the migrationtimescale found for Accreting sink particles in locally-isothermalmodels with those obtained using the Killing sink method pre-viously described. Whilst the Killing sinks maintain their initialmasses throughout, the Accreting sink models lead to increasinglymassive protoplanets, sometimes growing to several times their ini-tial mass. This makes interpreting the relative migration rates some-what di ffi cult, though the broad result is that including accretion c (cid:13)000 , 000–000 B.A. Ayli ff e & M.R. Bate leads to faster migration rates. Comparing an Accreting sink whichhas grown from an initial mass of 10 M ⊕ to beyond 33 M ⊕ with aKilling sink which maintains a mass of 33 M ⊕ throughout its evo-lution must be done with some care as there are several factorscontributing to their di ff ering migration timescales. One factor isthat the Accreting sink has a smaller r acc than the 33 M ⊕ Killingsink model as the former is calculated for its initial mass of 10 M ⊕ ,and so it develops a di ff erent circumplanetary structure. As was dis-cussed in section 3.1 for the di ff erent sized Killing sinks, a smaller r acc leads to faster migration. This factor therefore contributes somefraction of the di ff erence in migration rates seen in Fig. 5.As seen in Fig. 4, and again in comparing the left and rightpanels of Fig. 5, increasing r acc reduces the migration rates ofthe Killing sinks. However, the Accreting sink migration rates arelargely una ff ected by the change in r acc from 0.5 to 1 R H . Thetorques experienced by these sinks are equivalent to those of theKilling sinks. It therefore must be the case that the loss of nega-tive torques from within R H as r acc is increased is countered in thecase of Accreting sinks by greater mass accretion. To test this weperformed a series of calculations in which Accreting sinks wereembedded in discs with which they did not gravitationally inter-act. These sinks accrete only material that they sweep up from theirpath through the disc as they cannot gravitationally draw in gas.This also means that they experience no torques from the disc. Anymigration of these protoplanets is a result solely of the dilution oftheir specific angular momentum due to accretion of pressure sup-ported, and so sub-Keplerian, gas.To test this we numerically integrate from a protoplanet’s ini-tial mass and radius to its final mass, accounting for its accretionof gas with sub-Keplerian angular momentum to compare its finalorbital radius with that obtained from the SPH models. In our un-perturbed disc gas, at the initial orbital radius of the protoplanet,the specific angular momentum of the gas has a value of 0.998 ofthe protoplanet’s specific angular momentum. However, we notefrom the SPH models that the mean radius from which gas is ac-creted is smaller than r p , and so the average specific angular mo-mentum is somewhat smaller than this value. Employing this mea-sured value in our simple integration we obtain final protoplanetorbital radii that are in good agreement with those measured fromthe no-interaction SPH models. Applying this iterative calculationto the Accreting sink models shown in Fig. 5 can account for asurprisingly large fraction of the change in their orbital radii. Fora 10 M ⊕ protoplanet, about 50% of the change can be attributedto the e ff ects of accreting sub-Keplerian material. This implies upto a factor of 2 di ff erence in the migration timescale between anAccreting sink and a Killing sink due to gas accretion.Fig. 6 includes all the di ff erent accretion radii used in our Ac-creting sink calculations. It can be seen that those with radii of0.1 R H (marked as plus symbols) and masses (cid:46) . Jupiter closelyfollow the analytic model of Tanaka et al. (2002). In these casesthe dilution of the angular momentum and the e ff ective change inthe ratio of r acc to R H are minimal due to the small accretion radii,and lower accretion rates. The accretion rates onto these Accretingsinks are still substantial because of the evacuation of gas within r acc , which leads to an inward pressure gradient near their bound-ary where the gas is removed. This e ff ectively pulls material in at afaster rate than would be expected in the absence of the hole. Thisartificial evacuation at the planet boundary not only leads to artifi-cially fast accretion, it also alters the density structure around theprotoplanet out to distances that are large compared to r acc ; this istackled in section 3.4. Figure 6.
Migration timescales for 10, 33, 100, and 333 M ⊕ protoplan-ets undergoing mass and angular momentum accretion. Each protoplanet ismodelled as an Accreting sink particle, with r acc = R H (plus signs), 0.5 R H (diamonds), or 1.0 R H (asterisks) for each mass. The analytic Type Imodel of Tanaka et al. (2002) is plotted as a dashed line. The smallest ac-cretion radii protoplanet models give the migration rates that most closelymatch the analytic model. All of our models discussed to this point have included the gravi-tational interaction of the protoplanet with the disc and vice versa,but have neglected the disc’s self-gravity, that is its influence uponitself. For a low mass disc this is a reasonable approximation asthe central star is so dominant. Indeed the inclusion of self-gravityin low-mass discs has been found to have only a small impact onmigration rates by several authors (Nelson & Benz 2003; Pierens& Hur´e 2005; Baruteau & Masset 2008b; Crida et al. 2009). Com-paring two series of identical calculations, with the inclusion ofself-gravity being the only di ff erentiating characteristic, we foundthere to be a very small, but consistent increase in the migrationtimescales for self-gravity models. Fig. 7 makes evident the smallmagnitude of the change due to self-gravity, in this instance for Ac-creting sinks with r acc = . R H . However, self-gravity remains anessential component in building a self-consistent model.When substantial circumplanetary discs or envelopes form,the inclusion of self-gravity becomes important in delivering ac-curate migration rates. As discussed in Crida et al. (2009), a cir-cumplanetary disc is locked to the protoplanet, and moves with itaround the central star. In the absence of the disc self-gravity thecircumplanetary disc does not experience the torques due to theprotoplanetary disc; the torques which the protoplanet is respond-ing to. As such it does not migrate of its own accord, but must bepulled by the protoplanet, artificially slowing the rate of migration;known as the ‘inertial mass problem’ (Crida et al. 2009). To explorethe impact of this phenomenon we performed two calculations us-ing point masses with surfaces (discussed in the next section) tomodel 333 M ⊕ protoplanets, one model with and one without the c (cid:13) , 000–000 lanetary migration and growth Figure 7.
Comparisons between the migration timescales of protoplanets,modelled in locally-isothermal discs by Accreting sinks with r acc = . R H , in non-self-gravitating discs (plus symbols), and self-gravitating discs(asterisks). As has been found by several authors, the inclusion of disc self-gravity has a very small e ff ect upon the migration timescales. This smallimpact is consistent across all protoplanet masses, slightly increasing themigration timescales. The analytic Type I model of Tanaka et al. (2002) isplotted as a dashed line. inclusion of disc self-gravity. The two resulting migration evolu-tions are shown in Fig. 8. In both cases a protoplanetary disc devel-ops, and in the case without disc self-gravity, this causes the pro-toplanet’s rate of migration to slow by approximately 14 per cent.Having illustrated the e ff ect, we proceed with calculations whichall include self-gravity to avoid what is an unphysical interaction. The final addition to our locally-isothermal models was to use ourplanet surface treatment. This addresses the failing of our Accret-ing sink models by allowing the gas in the protoplanet’s vicinity todevelop a structure free of the undesired e ff ects of evacuating gasnear the protoplanet. Combining this treatment with the inclusionof self-gravity we obtain a model in which the interaction of thedisc and protoplanet is free of any artificially imposed radial limit,such as excluding material within R H . Instead, the radius at whichgas becomes bound to the protoplanet, and stops applying orbit af-fecting torques is achieved self-consistently. Fig. 9 illustrates themidplane density distribution around both an Accreting sink anda point mass with surface (in both cases a 333 M ⊕ protoplanet, ofradius 0.1 R H ). Whilst gas locked in a circumplanetary disc willnot exert migration causing torques, gas beyond such a disc (foundto be (cid:38) R H / ff e & Bate 2009b) that moves inand out of the Hill sphere will. The density distribution in the twocases are not similar until radii of (cid:38) R H , suggesting that the Accret-ing sinks experience smaller torques from within and around R H than the point masses with surfaces. The importance of these localtorques can be seen in the left panel of Fig. 10 which is a map of the Figure 8.
The orbital evolution of a 333 M ⊕ protoplanet modelled using apoint mass with surface of radius 0.03 R H , using a locally-isothermal equa-tion of state. A calculation including disc self-gravity is shown using a nar-row line, whilst the non-self-gravity calculation is shown using a thick line.There is reduction in the migration rate of the protoplanet modelled with-out self-gravity as it develops a circumplanetary disc; a result of the inertialmass problem. The migration rate is slower by ≈
14 per cent.
Figure 9.
Midplane densities around a 333 M ⊕ protoplanet modelled us-ing a locally-isothermal equation of state and self-gravity after 50 orbits.The solid line denotes the densities around a protoplanet modelled using apoint mass with a surface of radius 0.1 R H . The dashed line gives the mid-plane densities around an identical protoplanet modelled by an Accretingsink with r acc = . R H . The evacuation of material at the Accreting sink’sboundary leads to an inward pressure gradient that pushes in gas, givingartificially high accretion rates, and reducing the density in the vicinity ofthe growing planet out to radii of (cid:38) R H .c (cid:13) , 000–000 B.A. Ayli ff e & M.R. Bate
Figure 10.
Torque distributions around a 333 M ⊕ protoplanet modelled by a point mass with a surface of radius 0.03 R H after 50 orbits. The left hand panelis from a locally-isothermal model, whilst the right hand panel is from a radiation hydrodynamic calculation using interstellar opacities; both cases includeself-gravity, and the figures exclude torques from within R H /
3. The Roche lobe is marked on in each case using a white line. All the torques above the whitedashed line are positive, whilst all those below are negative. Introducing radiative transfer smears out the spiral arms, fills in the corotation region, and reducesthe maximum torques which act from within the Roche lobe. The section shown is r p ± R H , which avoids the inner and outer boundary regions. torques acting on the protoplanet in the plane of the disc. The fig-ure excludes torques acting from within R H / R H / c (cid:13) , 000–000 lanetary migration and growth Figure 11.
Sections of the whole disc models plotted as midplane density maps, focussing on the region surrounding the embedded protoplanets, with velocityvectors overplotted. It can be seen that material is flowing into the Hill sphere and back out again, reaching significant depths within this region. This materialis able to exchange angular momentum with the protoplanet before escaping its potential, and so a persistent torque is established from within R H . From leftto right the protoplanet masses are 33, 100, and 333 M ⊕ . From top to bottom the models are locally-isothermal, and then radiation hydrodynamical with 1, 10,and 100 per cent interstellar grain opacities. Fig. 12 includes migration timescales from Accreting sinkmodels and point mass with surface models with accretion / surfaceradii of 0.1 R H . Throughout the low mass regime both types ofmodel lead to results in good agreement with the analytic modelof Type I migration. At and beyond a Saturn mass, where Type IImigration becomes e ff ective, the point mass with surface calcula-tions show slightly faster ( ≈ factor of 2) rates of migration. This is a result of the less evacuated gaps formed in sink with surfacemodels, which results in a faster migration mechanism which isintermediate between Type I and II. The gap is more thoroughlyevacuated when using an Accreting sink because it establishes noenvelope structure to hinder accretion, and a negative pressure gra-dient to encourage it. This strips the corotational region of materialquickly, resulting in a clearer gap at an equivalent time. However, c (cid:13) , 000–000 B.A. Ayli ff e & M.R. Bate in both cases the gap is becoming more evacuated with time, whichwould be expected to slow migration to a Type II rate; this is aresult of the models relatively short evolution of just 50 orbital pe-riods starting from a uniform disc.
The structure of protoplanetary discs, and circumplanetary en-velopes and discs, are the result of competing gravitational andthermal e ff ects. In the locally-isothermal models discussed to thispoint, the thermal e ff ects associated with embedding a protoplanetin a protoplanetary disc have not been accounted for. To addressthis shortcoming we finally introduce radiative transfer and a real-istic equation of state into the migration models. The inclusion ofa better thermodynamic treatment has several potentially importantramifications. The accretion energy released near a growing pro-toplanet warms its surroundings, which in concert with disc self-gravity leads to the development of self-consistent local disc andenvelope structures. For example, a more poorly defined disc gap(opened by high mass planets) develops when radiative transfer isincluded (Fouchet & Mayer 2008; Ayli ff e & Bate 2009a). Com-pressional heating also tends to smear out the spiral density wavesthat are launched by planets embedded in the discs, giving themlower peak densities than in locally-isothermal models.The impact of these changes can be seen in the panels ofFig. 10, which illustrate the torques acting on a 333 M ⊕ protoplanetin both a locally-isothermal model (left panel) and an interstel-lar opacity radiation hydrodynamical (RHD) model (right panel).The peak torques acting near the protoplanet are lower in the RHDcase, though the maxima occur in similar locations (near the Rochelobe), illustrating the impact of the reduced local density. In thelocally-isothermal case, the azimuthal angle at which torques fromthe corotation region fall away is considerably smaller than in theRHD case. This is a result of the more poorly cleared gap when us-ing RHD. Finally, the thermal spreading of the spiral arms can beseen by the broadened regions from which the associated torquesare exerted, though there is no discernible di ff erence in the mag-nitude of the torque acting from wave crest. Densities along a linein the midplane, 10 degrees clockwise of the same 333 M ⊕ proto-planet in both locally-isothermal, and RHD models, are shown inthe right hand panel of Fig. 13. This distribution directly illustratesthe di ff erence in the vacuity of the gap. It is also clear that the den-sity maximum of the (outer) spiral arm is similar in both instances,despite the breadth being altered by the heating. The left panel isan equivalent plot for a 100 M ⊕ protoplanet. In this case the spiralarm density maximum is reduced by about a factor of 2 comparedwith the locally-isothermal case, possibly weakening the di ff eren-tial Lindblad torques. However, the disc gap contains more massthan in the locally-isothermal model, once again leading to a mi-gration process even more similar to Type I, and so faster, than theexpected Type II. The overall impact of using radiation hydrody-namics in this case will depend on the relative impact of these twoopposing e ff ects.The migration rates found in these interstellar grain opacity ra-diation hydrodynamical models are compared with those from oth-erwise equivalent locally-isothermal models in Fig. 14. In the highmass regime, the impact is very small, with migration rates alteredby no more than 30%. The 333 M ⊕ protoplanet, which maintainsthe same spiral arm densities but a less evacuated gap, migratesfaster with RHD. This is consistent with an intermediate migrationprocess, between Types I and II, as a result of the less evacuated gapformed with RHD. Once again, it is important to note that for both Figure 12.
The migration timescales of protoplanets modelled using Ac-creting sinks (asterisks) compared with those obtained using point masseswith surfaces (plus signs), with accretion / surface radii of 0.1 R H . Thesetimescales are from models of locally-isothermal self-gravitating discs. Theanalytic Type I model of Tanaka et al. (2002) is plotted as a dashed line.With small accretion radii, the Accreting sinks give results inline with theanalytic expectation for Type I migration and with the more realistic pointmass with surface models. Beyond a Saturn mass, in the Type II regime, theAccreting sinks show slower migration by factors of ≈ high mass protoplanets the change in migration timescales whenradiative transfer is introduced are very small. In the lower massType I regime, the introduction of radiative transfer leads to sub-stantial increases in the migration timescales. A 50 M ⊕ protoplanettakes almost 4 times longer to migrate into its central star whenmodelled with radiative transfer than it does without, whilst a 33M ⊕ protoplanet takes more than 10 times longer. A 10 M ⊕ pro-toplanet is slowed by more than a factor of 2. In these cases, thedensities around the protoplanets out to several R H are lower whenmodelled with RHD using interstellar grain opacities than in thelocally-isothermal models; about a 33 M ⊕ protoplanet the gas den-sity at the disc midplane at a radius of R H is a factor of 3 lower. Sim-ilar reductions in Type I migration rates have been found in otherrecent work considering non-isothermal discs (e.g. Paardekooper& Mellema 2006). Recent models suggest that the slowing resultsfrom the dependence of the co-orbital torques on the radial entropygradient of the protoplanetary disc, and that under some circum-stances these torques can be su ffi cient to reverse the direction ofmigration (Paardekooper & Papaloizou 2008a; Baruteau & Masset2008b; Kley & Crida 2008; Kley et al. 2009; Paardekooper et al.2010). The reduced migration rates for low mass protoplanets arean important result in trying to increase the survival probability ofproto-giant cores. Specifically it may justify the reduction, or somefraction of the reduction in Type I migration rates assumed in coreaccretion models that seek to ensure planets survive (Alibert et al.2004, 2005; Ida & Lin 2008; Mordasini et al. 2009).In the point mass with surface models the envelope / disc sur-rounding the protoplanets is structured realistically due to the influ- c (cid:13) , 000–000 lanetary migration and growth Figure 13.
Midplane densities along lines 10 degrees clockwise of 100 M ⊕ (left) and 333 M ⊕ (right) protoplanets, modelled with point masses with surfaces of0.03 R H after 50 orbits. The profiles run over r p ± . p . The solid lines mark the profiles for locally-isothermal calculations, whilst the dashed lines are takenfrom radiation hydrodynamics calculations using interstellar grain opacities. In both cases, the introduction of radiative transfer leads to poorer evacuation ofthe corotation region, allowing for greater torques to act from these radii. For the 100 M ⊕ protoplanet the outer spiral arm can be seen to have a lower peakdensity, a result of the heat produced in the shock. Figure 14.
The migration timescales of protoplanets modelled using pointmasses with surfaces of radius 0.03 R H in locally-isothermal (asterisks) andinterstellar grain opacity radiation hydrodynamics (plus symbols) models.The analytic Type I model of Tanaka et al. (2002) is plotted as a dashedline. In the mass regime corresponding to Type I migration, the inclusionof radiative transfer with high opacities leads to slower rates of migration.For a 33 M ⊕ protoplanet the migration timescale is increased by an order ofmagnitude. Such slowed migration may improve the chances of proto-giantcores surviving long enough to grow to giant planet masses. Figure 15.
Scaleheights measured in discs containing a 33 M ⊕ protoplanetmodelled using a point mass with surface of radius 0.03 R H . The resultsof a locally-isothermal model are shown using a solid line, along with theradiation hydrodynamics models using 100 per cent (long-dashed), 10 percent (short-dashed), and 1 per cent (dotted) interstellar grain opacities. Thelocally-isothermal and reduced opacity radiation hydrodynamics models allyield similar scaleheights through the disc, but the full opacity model givesa thicker disc, most notably around the radius of the embedded protoplanet.c (cid:13) , 000–000 B.A. Ayli ff e & M.R. Bate ences of gravity, self-gravity, and thermodynamics. It is thereforepossible to account for the di ff erentiation of the angular momen-tum of accreted material into the planet’s orbital angular momen-tum, and into rotation about the planet. This di ff erentiation has nextto no impact upon the compound protoplanet’s orbital angular mo-mentum. Typically a fraction of ∼ − of the angular momentumof the gas accreted into a circumplanetary disc is diverted into therotation of that disc about the planet.The impact of radiative transfer may go beyond the immedi-ate vicinity of a protoplanet, for example, inflating a protoplanetarydisc’s scaleheight as a result of viscous heating at the midplane.Thicker discs result in smaller torques from the inner and outerLindblad resonances, scaling as ( H / R ) − (Ward 1997), whilst themismatch between the torques scales as H / R . This means that themigration rates scale as ( H / R ) − , and a hotter, thicker disc shouldlead to slower migration (Masset & Kley 2006). Fig. 15 shows themeasured value of H / R against radius for a disc with an embedded33 M ⊕ protoplanet. It can be seen that there is a larger scaleheightat and about the protoplanet’s orbital radius in a radiation hydro-dynamics calculation using interstellar opacities when comparedwith all the other calculations. The scaleheight is a factor of ≈ . Our final results come from varying the protoplanetary disc’s opac-ity to determine its e ff ect upon migration. Fig. 16 shows that thee ff ect of varying the opacity on the migration rates is very smallfor the high mass planets. For the low mass planets there is a trendof increasing migration timescale with increased opacity. Interest-ingly, protoplanets in the low opacity RHD models migrate fasterthan in the locally-isothermal models.The largest change in the migration timescales of low massprotoplanets is seen when employing the interstellar grain opaci-ties, with the change between the 1% and 10% models, and thelocally-isothermal models being somewhat smaller. In all but theinterstellar grain opacity calculations, the densities around a givenprotoplanet are almost identical, and so the local torques shouldbe similar. The spiral arms are also likely to be most broadened,and their peak density most reduced, in the highest opacity cal-culations, with lower opacities leading to structures more similarto those found in the locally-isothermal models. These structures,and maximum local torques can be seen in the torque maps shownin Fig. 17, which are produced for a 33 M ⊕ protoplanet modelledusing a point mass with a surface radius of 0.03 R H . The torquesacting at and around the Roche lobe can be seen to be weaker inthe interstellar opacity case (right most panel), than in any of theother cases, and the spiral arms most diminished.The diminution of the spiral arms with increasing opacity canbe further seen in Figs. 18 and 19, which are surface density anddensity weighted temperature maps respectively of discs with em-bedded protoplanets, once again modelled by point masses withsurfaces of radius 0.03 R H . From left to right the protoplanet massesare 33, 100, and 333 M ⊕ respectively, and from top to bottom themodels are locally-isothermal, and RHD with 1, 10, and 100 percent interstellar grain opacities respectively; we have omitted the10 M ⊕ cases because the disturbance caused by the protoplanet isbarely visible. Moving from the locally-isothermal density distribu-tion (top-row of figure 18) to the interstellar opacity case (bottomrow) it is evident that where a gap forms, it is less clear at higheropacities, and that the density waves propagating outwards are less Figure 16.
Migration timescales for protoplanets in radiation hydrodynam-ics calculations, modelled by point masses with surface radii of 0.03 R H ,with three di ff erent opacities: 100 per cent (plus symbols), 10 per cent(asterisks), and 1 per cent (diamonds) interstellar grain opacities. The re-duced opacity models yield shorter migration timescales than their locally-isothermal equivalent models (c.f. Fig. 14), with only the full opacity caseleading to slower migration. Increasing the grain opacity of the disc leadsto slower migration rates, with the e ff ect becoming smaller for higher massprotoplanets. The analytic Type I model of Tanaka et al. (2002) is plotted asa dashed line. peaked. Fig. 11 shows midplane densities and velocity vectors inthe vicinity of the protoplanets. The velocity vectors make plain thespiral shocks, horseshoe orbits, and the transition to circumplane-tary flow deep within the Hill sphere. 20 shows density weightedtemperature maps for the same models but concentrating on just thearea surrounding the protoplanet.The diminution of the propagating density waves with increas-ing opacity maybe in part due to the process of wave channelling. Ina vertically isothermal disc, these density perturbations propagateas plane waves (Lubow & Pringle 1993). However, once radiativetransfer is introduced, a vertical temperature profile is established,in which the midplane of the disc is hottest due to the e ff ects ofcompression and viscous heating. This leads to wave channelling,causing the energy of the wave to quickly become confined to thesurfaces of the disc, where it is dissipated (Korycansky & Pringle1995; Lubow & Ogilvie 1998; Ogilvie & Lubow 1999; Bate et al.2002). This results in the more limited propagation of the spiralarms seen in the non-isothermal models, and which is clearest forlow mass planets where the weak spirals barely last an entire plan-etary orbital period. It is also possible to see in Fig. 11 that thebreadth of the spiral arms increases with opacity, as was discussedearlier.In Fig. 19 and Fig. 20 it is possible to see the hot gasthat develops around the protoplanets, particularly those of highermass, when modelled with RHD and interstellar grain opacities.In Fig. 19 it is also possible to see the heat produced in the spiralshocks. The spatial extent of the hot gas surrounding a protoplanet c (cid:13) , 000–000 lanetary migration and growth Figure 17.
Torque distributions around a 33 M ⊕ protoplanet modelled by a point mass with a surface of radius 0.03 R H after 50 orbits. From left to right themodels are locally-isothermal, and then RHD with 1 per cent, 10 per cent, and 100 per cent interstellar grain opacity respectively, all including self-gravity.The Roche lobe is marked on in each case using a white line. All the torques above the white dashed line ( y =
0) are positive, whilst all those below arenegative. The spiral arms are particularly poorly defined in the interstellar grain opacity case (right most panel), and the torques acting from around the Rochelobe are weakest in this case. reduces with the reduction of the opacity, which allows the lessdense outer reaches to cool more easily. It is also possible to seethe higher temperatures in the spiral shocks, even in those wavesthat have propagated far from the planet, when compared with theunderlying temperature distribution shown in the top row.
It seems that when planets form, they should also migrate. Thismigration poses problems for the survival of gas giant cores. Con-versely, such migration may aid in planet growth by opening upundepleted regions of the disc for these protoplanets to accrete (Al-ibert et al. 2004). It is known from the ever increasing number ofexoplanet detections that a large population of planets do grow andsurvive around many di ff erent types of star. The work discussedin this paper, like the work of many authors, was undertaken totry and ascertain the influence of di ff erent physical processes uponmigration. The models we present in section 3.5.1 are the most self-consistent to date, including self-gravity, radiative transfer, and pro-toplanets modelled with surfaces to produce realistic envelopes andcircumplanetary discs. We first made comparisons with the results of a grid code to ensurethat our SPH model was suited to tackling the problem of migra-tion, a capacity in which it has been seldom used before. The abil-ity of the SPH method to model protoplanet migration was consid-ered by de Val-Borro et al. (2006), who made comparisons between many grid codes and two SPH codes. They found reasonable agree-ment in disc structures around embedded protoplanets, though theSPH calculations under-resolved some disc features. This di ff er-ence is mitigated to an extent in our models by the increased num-ber of particles, more than 6 times greater than the number usedin these previous comparisons. Our migration rates obtained us-ing non-self-gravitating locally-isothermal calculations are in goodagreement with the analytic model of Tanaka et al. (2002), andso also with many previous similar numerical models which showsimilar agreement to the analytic model. We introduced disc self-gravity to examine its impact upon mi-gration. Pierens & Hur´e (2005) examined the inclusion of discself-gravity analytically. They broke down the analysis into twocompeting influences, that of the circumstellar disc-planet inter-action, and the circumstellar disc-circumplanetary disc interaction(gas self-gravity). All of our SPH models include the first of thesetwo interactions, which Pierens & Hur´e (2005) found to acceler-ate Type I migration. So by introducing self-gravity, we mean theaddition of the second interaction. Pierens & Hur´e found this toasymmetrically shift the outer and inner Lindblad resonances suchthat the migration rate reduces. More recently, numerical modelswere performed by Baruteau & Masset (2008b) to further explorethe impact of including self-gravity. They ran models to explore theimpact of both interactions discussed here, and found, in agreementwith Pierens & Hur´e (2005) and the previous numerical models ofNelson & Benz (2003), that the disc-disc interaction reduced theType I migration rates. The consistent, though small, increases in c (cid:13)000
It seems that when planets form, they should also migrate. Thismigration poses problems for the survival of gas giant cores. Con-versely, such migration may aid in planet growth by opening upundepleted regions of the disc for these protoplanets to accrete (Al-ibert et al. 2004). It is known from the ever increasing number ofexoplanet detections that a large population of planets do grow andsurvive around many di ff erent types of star. The work discussedin this paper, like the work of many authors, was undertaken totry and ascertain the influence of di ff erent physical processes uponmigration. The models we present in section 3.5.1 are the most self-consistent to date, including self-gravity, radiative transfer, and pro-toplanets modelled with surfaces to produce realistic envelopes andcircumplanetary discs. We first made comparisons with the results of a grid code to ensurethat our SPH model was suited to tackling the problem of migra-tion, a capacity in which it has been seldom used before. The abil-ity of the SPH method to model protoplanet migration was consid-ered by de Val-Borro et al. (2006), who made comparisons between many grid codes and two SPH codes. They found reasonable agree-ment in disc structures around embedded protoplanets, though theSPH calculations under-resolved some disc features. This di ff er-ence is mitigated to an extent in our models by the increased num-ber of particles, more than 6 times greater than the number usedin these previous comparisons. Our migration rates obtained us-ing non-self-gravitating locally-isothermal calculations are in goodagreement with the analytic model of Tanaka et al. (2002), andso also with many previous similar numerical models which showsimilar agreement to the analytic model. We introduced disc self-gravity to examine its impact upon mi-gration. Pierens & Hur´e (2005) examined the inclusion of discself-gravity analytically. They broke down the analysis into twocompeting influences, that of the circumstellar disc-planet inter-action, and the circumstellar disc-circumplanetary disc interaction(gas self-gravity). All of our SPH models include the first of thesetwo interactions, which Pierens & Hur´e (2005) found to acceler-ate Type I migration. So by introducing self-gravity, we mean theaddition of the second interaction. Pierens & Hur´e found this toasymmetrically shift the outer and inner Lindblad resonances suchthat the migration rate reduces. More recently, numerical modelswere performed by Baruteau & Masset (2008b) to further explorethe impact of including self-gravity. They ran models to explore theimpact of both interactions discussed here, and found, in agreementwith Pierens & Hur´e (2005) and the previous numerical models ofNelson & Benz (2003), that the disc-disc interaction reduced theType I migration rates. The consistent, though small, increases in c (cid:13)000 , 000–000 B.A. Ayli ff e & M.R. Bate migration timescales found for the SPH calculations including self-gravity (Fig. 7) are in agreement with these previous results. Whilstthe e ff ect is small in our models, in discs with masses several timesthe minimum mass solar nebula self-gravity might slow migrationby a more considerable factor (Baruteau & Masset 2008b). Introducing radiative transfer has allowed us to investigate migra-tion beyond the more common locally-isothermal models. Previousstudies concerned with the impact of non-isothermality, have foundthat the inclusion of more complicated thermodynamics leads to in-creased migration timescales (Jang-Condell & Sasselov 2005), andeven the reversal of the Type I migration direction (Paardekooper& Mellema 2006, 2008b; Kley & Crida 2008; Kley et al. 2009). Weare also able to obtain outward migration for a 10 M ⊕ protoplanetwhen modelling it with an Accreting sink particle of radius 0.03 R H with radiation hydrodynamics in an interstellar grain opacity disc,as shown in Fig. 21.However, once we progress to using our point mass with sur-face model, which is the most realistic case, we find no cases ofoutward migration. The 10 M ⊕ case modelled in this way, with asurface radius of 0.03 R H and using interstellar grain opacities, isshown in Fig. 21 for comparison with the Accreting sink model.Kley et al. (2009) find outward migration to result from a densitypeak in the planet’s immediate vicinity, just ahead of it, which pullsit around in its orbit, increasing its angular momentum and so caus-ing it to migrate outwards. To test the reliability of this result theyexperiment with di ff erent resolutions and various gravitational soft-ening lengths, the shortest being 0.5 R H , to explore their impact onthe resulting migration. The resulting outward migration was foundto be invariant to these changes. The migration rates are based oncalculations of torques acting on planet at a fixed orbital radius of5.2 AU, where the torques are scaled depending upon their distancefrom the protoplanet, d , according to a fraction obtained from thetapering function f b ( d ) = (cid:34) exp (cid:32) − d / R H − bb / (cid:33) + (cid:35) − . (5)This function is such that at b , the torque cut-o ff radius in units of R H , the torques are halved (where they typically use b = . / R H .Changing the value of b to 0.6 altered their total torques by 10per cent in locally-isothermal models, and 30 per cent in the radi-ation hydrodynamics models, but the resolution limits their abilityto usefully reduce the value of b further. In Fig. 10, where our un-derlying model resolution is ≈
20 times higher than that of Kleyet al. (2009), it can be seen that the greatest torques are actingfrom within the Roche lobe. In fact the peak torques are actingfrom a radius of ≈ . R H in the radiation hydrodynamics casewhere gas is still passing through the Hill sphere, and able to carryaway angular momentum. The ability of the SPH model to resolvethese torques clearly, removes the need to use a tapering function ofthe kind above, which otherwise would have scaled these torquesdown to 73 per cent of their true values. Furthermore, using ourplanet surface treatment yields a self-consistent envelope around a20 M ⊕ protoplanet with densities that are an order of magnitudehigher than those obtained using gravitational softening in Kleyet al.’s otherwise similar disc model. Such a change in density issignificant given the source of the positive torques suggested byKley et al.. Given the improved resolution and realism of the planet treatment employed in our SPH models, our results suggest thatoutward migration is not an inevitable consequence of introducingmore complex thermodynamics.We find that introducing radiative transfer to our most real-istic models of low mass protoplanets in high opacity discs canslow migration rates by up to an order of magnitude. The reduc-tion of Type I migration rates may be a result of enhanced corota-tion torques due to entropy gradients across the co-orbital region(Paardekooper & Papaloizou 2008a; Baruteau & Masset 2008b;Kley & Crida 2008; Kley et al. 2009; Paardekooper et al. 2010).For high mass protoplanets, those capable of forming circumplan-etary discs (Ayli ff e & Bate 2009b), the dominance of gravity overthermal e ff ects leaves the migration rates very little altered uponthe introduction of radiative transfer when compared with locally-isothermal conditions. This competition between thermal and grav-itational influences was also found to play an important role in de-termining accretion rates onto growing protoplanets by Ayli ff e &Bate (2009a). As with migration, the accretion rates in RHD mod-els were found to be most dependent upon the grain opacity forlow mass protoplanets, where the gravitational forces at work aresmaller. Ayli ff e & Bate (2009a) found that a 33 M ⊕ protoplanet, inan interstellar grain opacity disc, is expected to double in mass in ≈ × years. We find the migration timescale for the same pro-toplanet to be ≈ × years. Under the same conditions a 100 M ⊕ protoplanet will double in mass in ≈ × years, and has a mi-gration timescale of ≈ × years. Whilst a 333 M ⊕ protoplanetdoubles in mass in ≈ × years and migrates into its centralstar in ≈ years. These results suggest that is should be commonfor a protoplanet embedded in an interstellar grain opacity disc toreach a Jupiter mass. However, to reach larger masses may requirehigher disc masses or some additional process to either slow a pro-toplanet’s migration, or accelerate its accretion.In the low mass regime the increase in migration timescalesthat we find using RHD may help remedy the problem of proto-giant core depletion that was identified by Alibert et al. (2004),Alibert et al. (2005), Ida & Lin (2008) and Mordasini et al. (2009).Ida & Lin showed that the e ffi ciency of Type I migration must bereduced by at least an order of magnitude to enable the formationof a population of planets that match current observed prevalences.We find that the inclusion of radiative transfer, in conjunction witha realistic treatment of accretion on to the protoplanet can providejust such a reduction in the low mass regime.As was discussed in the results section, in the Type II regimethe introduction of radiative transfer has two significant e ff ectswhich have opposing impacts on migration rates. The reductionin the peak density of the spiral arms due to heating results inlower torques from these features, slowing migration, whilst thepoorer evacuation of the disc gap leads to a migration process thatis intermediate between Types I and II. In all our high mass caseswhich are modelled with either Accreting sinks, or with surfaces,in locally-isothermal or radiation hydrodynamics models, the rateof migration is considerably faster than the discs viscous migrationtimescale. Only using the unrealistic Killing sinks do our Type IIrates approach the analytic model for Type II migration of Ward(1997), most likely because they can most rapidly clear a gap ofmaterial. Some component of the rapid high mass migration ratesfound here may result from the models starting with an unperturbeddisc, which over 50 orbits fails to establish a gap of the depth thatmight be expected if the planet had evolved to its initial mass insitu. However, radiative transfer might well be expected to hinderthe transition from Type I to the slower Type II migration due to in- c (cid:13) , 000–000 lanetary migration and growth Figure 18.
Surface density maps of protoplanetary discs after 50 orbits of the embedded protoplanets. From left to right the protoplanet masses are 33, 100,and 333 M ⊕ . From top to bottom the models are locally-isothermal, and then radiation hydrodynamical with 1, 10, and 100 per cent interstellar grain opacities.c (cid:13) , 000–000 B.A. Ayli ff e & M.R. Bate
Figure 19.
Density weighted temperature maps of protoplanetary discs after 50 orbits of the embedded protoplanets. From left to right the protoplanet massesare 33, 100, and 333 M ⊕ . From top to bottom the models are locally-isothermal, and then radiation hydrodynamical with 1, 10, and 100 per cent interstellargrain opacities. c (cid:13) , 000–000 lanetary migration and growth Figure 20.
Magnified sections of the whole disc midplane temperature maps shown in Fig. 19, focussing on the region surrounding the embedded protoplanets.Again, from left to right the protoplanet masses are 33, 100, and 333 M ⊕ . From top to bottom the models are locally-isothermal, and then radiation hydrody-namical with 1, 10, and 100 per cent interstellar grain opacities. It is particularly apparent for the 333 M ⊕ protoplanet that the energy released as a result ofaccretion leads to a larger region of hot gas, which is less able to cool as the opacity is increased. It is also possible to see the broadening of the spiral arms asthe opacity is increased.c (cid:13) , 000–000 B.A. Ayli ff e & M.R. Bate
Figure 21.
The orbital evolution of a 10 M ⊕ protoplanet modelled threeways, twice as an Accreting sink particle and once as a point mass with sur-face. The Accreting sink models are di ff erentiated by the chosen equation ofstate, one locally-isothermal (marked isothermal above), and one using ra-diation hydrodynamics (marked RHD) which migrates outwards. The pointmass with surface is modelled using RHD, and shows a slow but inwardsmigration. creased thermal pressure that acts to oppose gap formation, as seenin our models. We have conducted models of protoplanet migration with a rangeof physics and a raft of protoplanet treatments. Beginning withlocally-isothermal, non-self gravitating discs, we incrementallyadded physics to reach self-gravitating radiation hydrodynamic cal-culations. We progressed from protoplanets modelled by crude sinkparticles, to those modelled by gravitating point masses with plan-etary surfaces. The final models we produce are the most self-consistent to date, and turn up some interesting results. In partic-ular, • We establish that the ‘inertial mass problem’, su ff ered by pro-toplanets that develop circumplanetary discs in models withoutself-gravity, is an e ff ect of order 10 per cent for minimum massmodels. • We find that Type I migration rates, those that a ff ect low massplanets (10-33 M ⊕ ), can be slowed by up to an order of magni-tude by the inclusion of radiative transfer in protoplanetary discswith interstellar grain opacities. Such a reduction helps to justifythe arbitrary reductions in the Type I migration rate used in planetsynthesis models (Alibert et al. 2004, 2005; Ida & Lin 2008; Mor-dasini et al. 2009) where such a reduction is required to ensureplanets survive beyond the disc lifetime. • The Type I migration timescale that we obtain for a 10 M ⊕ protoplanet, modelled using RHD in an interstellar grain opacitydisc, is of similar order to its mass doubling time, as found inAyli ff e & Bate (2009a), whilst for a 33 M ⊕ protoplanet the migra-tion timescale is an order of magnitude longer than the mass dou-bling time. Similarly for a 100 M ⊕ protoplanet the mass doublingtime is shorter than the migration timescale, but for such high massprotoplanets the migration timescale is largely unchanged between locally-isothermal and RHD models. These accretion and migra-tion timescales suggest that a protoplanet embedded in an interstel-lar grain opacity disc might grow to a high mass within its availablemigration time. • We find that the introduction of radiative transfer to a modelof a 10 M ⊕ protoplanet, represented by an Accreting sink particleembedded in an interstellar grain opacity disc, causes it to migrateoutwards. However the inclusion of radiative transfer does not leadto outward migration when we employ a well resolved and morerealistic model of the protoplanet with a surface. ACKNOWLEDGMENTS
We would like to thank the referee for his perceptive comments,and useful suggestions. The calculations reported here were per-formed using the University of Exeter’s SGI Altix ICE 8200 su-percomputer. Many visualisations were produced using SPLASH(Price 2007), a visualisation tool for SPH that is publicly availableat http: // / people / dprice / splash. MRB is gratefulfor the support of a Philip Leverhulme Prize and a EURYI Awardwhich also supported BAA. This work, conducted as part of theaward “The formation of stars and planets: Radiation hydrodynam-ical and magnetohydrodynamical simulations” made under the Eu-ropean Heads of Research Councils and European Science Foun-dation EURYI (European Young Investigator) Awards scheme, wassupported by funds from the Participating Organisations of EURYIand the EC Sixth Framework Programme. REFERENCES
Alexander D. R., 1975, ApJS, 29, 363Alibert Y., Mordasini C., Benz W., 2004, A&A, 417, L25Alibert Y., Mordasini C., Benz W., Winisdoer ff er C., 2005, A&A,434, 343Ayli ff e B. A., Bate M. R., 2009a, MNRAS, 393, 49Ayli ff e B. A., Bate M. R., 2009b, MNRAS, 397, 657Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Baruteau C., Lin D. N. C., 2010, ApJ, 709, 759Baruteau C., Masset F., 2008a, ApJ, 672, 1054Baruteau C., Masset F., 2008b, ApJ, 678, 483Bate M., 1995, PhD thesis, PhD thesis, Univ. Cambridge, (1995)Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362Bate M. R., Lubow S. H., Ogilvie G. I., Miller K. A., 2003, MN-RAS, 341, 213Bate M. R., Ogilvie G. I., Lubow S. H., Pringle J. E., 2002, MN-RAS, 332, 575Benz W., 1990, in Buchler J. R., ed., Numerical Modelling ofNonlinear Stellar Pulsations Problems and Prospects SmoothParticle Hydrodynamics - a Review. pp 269– + Benz W., Cameron A. G. W., Press W. H., Bowers R. L., 1990,ApJ, 348, 647Boley A. C., Hartquist T. W., Durisen R. H., Michael S., 2007,ApJ, 656, L89Bryden G., Chen X., Lin D. N. C., Nelson R. P., PapaloizouJ. C. B., 1999, ApJ, 514, 344Crida A., Baruteau C., Kley W., Masset F., 2009, A&A, 502, 679Crida A., Morbidelli A., 2007, MNRAS, 377, 1324D’Angelo G., Bate M. R., Lubow S. H., 2005, MNRAS, 358, 316D’Angelo G., Henning T., Kley W., 2002, A&A, 385, 647D’Angelo G., Kley W., Henning T., 2003, ApJ, 586, 540 c (cid:13) , 000–000 lanetary migration and growth D’Angelo G., Lubow S. H., 2008, ApJ, 685, 560de Val-Borro M., et al., 2006, MNRAS, 370, 529Fouchet L., Mayer L., 2008, ArXiv e-prints, 806Goldreich P., Tremaine S., 1980, ApJ, 241, 425Ida S., Lin D. N. C., 2008, ApJ, 673, 487Ivanov P. B., Papaloizou J. C. B., Polnarev A. G., 1999, MNRAS,307, 79Jang-Condell H., Sasselov D. D., 2005, ApJ, 619, 1123Klahr H., Kley W., 2006, A&A, 445, 747Kley W., 1999, MNRAS, 303, 696Kley W., Bitsch B., Klahr H., 2009, A&A, 506, 971Kley W., Crida A., 2008, ArXiv e-prints, 806Korycansky D. G., Pollack J. B., 1993, Icarus, 102, 150Korycansky D. G., Pringle J. E., 1995, MNRAS, 272, 618Li H., Lubow S. H., Li S., Lin D. N. C., 2009, ApJ, 690, L52Lodato G., Price D., 2010, ArXiv e-printsLubow S. H., Ogilvie G. I., 1998, ApJ, 504, 983Lubow S. H., Pringle J. E., 1993, ApJ, 409, 360Masset F., Kley W., 2006, Disk-planet interaction and migration.Cambridge University Press, pp 216– + Masset F. S., 2002, A&A, 387, 605Masset F. S., D’Angelo G., Kley W., 2006, ApJ, 652, 730Mayor M., Queloz D., 1995, Nature, 378, 355Monaghan J. J., 1992, ARA&A, 30, 543Monaghan J. J., 2002, MNRAS, 335, 843Monaghan J. J., Gingold R. A., 1983, Journal of ComputationalPhysics, 52, 374Mordasini C., Alibert Y., Benz W., 2009, A&A, 501, 1139Morohoshi K., Tanaka H., 2003, MNRAS, 346, 915Morris J. P., Monaghan J. J., 1997, J. Comput. Phys., 136, 41Nelson A. F., Benz W., 2003, ApJ, 589, 578Nelson R. P., 2005, A&A, 443, 1067Nelson R. P., Papaloizou J. C. B., 2003, MNRAS, 339, 993Nelson R. P., Papaloizou J. C. B., 2004, MNRAS, 350, 849Nelson R. P., Papaloizou J. C. B., Masset F., Kley W., 2000, MN-RAS, 318, 18Ogilvie G. I., Lubow S. H., 1999, ApJ, 515, 767Paardekooper S., Baruteau C., Crida A., Kley W., 2010, MNRAS,401, 1950Paardekooper S., Papaloizou J. C. B., 2008a, A&A, 485, 877Paardekooper S., Papaloizou J. C. B., 2009, MNRAS, 394, 2283Paardekooper S.-J., Mellema G., 2006, A&A, 459, L17Paardekooper S.-J., Mellema G., 2008a, A&A, 478, 245Paardekooper S.-J., Mellema G., 2008b, A&A, 478, 245Paardekooper S.-J., Papaloizou J. C. B., 2008b, A&A, 485, 877Pepli´nski A., Artymowicz P., Mellema G., 2008a, MNRAS, 386,179Pepli´nski A., Artymowicz P., Mellema G., 2008b, MNRAS, 387,1063Pierens A., Hur´e J.-M., 2005, A&A, 433, L37Pollack J. B., McKay C. P., Christo ff erson B. M., 1985, Icarus,64, 471Price D. J., 2007, Publications of the Astronomical Society ofAustralia, 24, 159Price D. J., Bate M. R., 2007, MNRAS, 377, 77Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337Springel V., Hernquist L., 2002, MNRAS, 333, 649Syer D., Clarke C. J., 1995, MNRAS, 277, 758Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257Ward W. R., 1986, Icarus, 67, 164Ward W. R., 1997, Icarus, 126, 261Whitehouse S. C., Bate M. R., 2006, MNRAS, 367, 32 Whitehouse S. C., Bate M. R., Monaghan J. J., 2005, MNRAS,364, 1367Yu C., Li H., Li S., Lubow S. H., Lin D. N. C., 2010, ArXiv e-prints c (cid:13)000