# Migration of accreting planets in radiative discs from dynamical torques

###### Abstract

We present the results of hydrodynamical simulations of the orbital evolution of planets undergoing runaway gas accretion in radiative discs. We consider accreting disc models with constant mass flux through the disc, and where radiative cooling balances the effect of viscous heating and stellar irradiation. We assume that 20-30 giant planet cores are formed in the region where viscous heating dominates and migrate outward under the action of a strong entropy-related corotation torque. In the case where gas accretion is neglected and for an viscous stress parameter , we find evidence for strong dynamical torques in accreting discs with accretion rates . Their main effect is to increase outward migration rates by a factor of typically. In the presence of gas accretion, however, runaway outward migration is observed with the planet passing through the zero-torque radius and the transition between the viscous heating and stellar heating dominated regimes. The ability for an accreting planet to enter a fast migration regime is found to depend strongly on the planet growth rate, but can occur for values of the mass flux through the disc of . We find that an episode of runaway outward migration can cause an accreting planet formed in the 5-10 AU region to temporarily orbit at star-planet separations as large as 60-70 AU. However, increase in the amplitude of the Lindblad torque associated with planet growth plus change in the streamline topology near the planet systematically cause the direction of migration to be reversed. Subsequent evolution corresponds to the planet migrating inward rapidly until it becomes massive enough to open a gap in the disc and migrate in the Type II regime. Our results indicate that a planet can reach large orbital distances under the combined effect of dynamical torques and gas accretion, but an alternative mechanism is required to explain the presence of massive planets on wide orbits.

###### keywords:

accretion, accretion discs – planet-disc interactions– planets and satellites: formation – hydrodynamics – methods: numerical^{†}

^{†}pagerange: Migration of accreting planets in radiative discs from dynamical torques–LABEL:lastpage

^{†}

^{†}pubyear: 2014

## 1 Introduction

Radial velocity surveys indicate that giant planets exist around of Sun-like stars (Cumming et al. 2008; Mayor et al. 2011),
with
a frequency that is strongly correlated to the metallicity of the host star. There exists a great diversity of orbital
architectures among this population of
exoplanets, which span a range of eccentricities and periods from days to several thousands of days. The abundance
of giant planets is found to depend strongly on the star-planet orbital separation. For example, hot Jupiters with semi major axes AU are found to orbit around only of main sequence stars, whereas there is an excess of Jupiter-mass planets located between 1-2 AU
(Cumming et al. 2008). Giant planets orbiting beyond the snow-line at 5-10 AU from their star also appear to be common, and
microlensing surveys (Shvartzvald et al. 2016) even suggest that half of the stars may host such planets.
Giant planets with mass have also been observed at larger orbital distances 10-100 AU by direct imaging (Brandt et al. 2014). The
system HR8799 (Marois et al. 2008, 2010) provides an example of such a planetary system that has been directly imaged. It contains four giant planets with mass
in the range and located at orbital distances AU. Interestingly, it has been suggested that these planets may be engulfed
in 8:4:2:1 Laplace resonances (Gozdziewski & Migaszewski 2014; Maire et al. 2015) that guarantee the stability of the system.

It has been proposed that the planets
observed by direct imaging may have formed by gravitational instability
(Dodson-Robinson et al. 2009) or through pebble accretion (Lambrechts et al. 2014). Due to the difficulty for forming cores that are massive enough to undergo runaway gas accretion at distances AU, explaining the in-situ
formation of massive
planets on wide orbits through core accretion does indeed appear challenging. In the context of this scenario, however, it can not be excluded
that these planets formed at distances AU where the accretion timescales are short compared to the lifetime
of the disc, and then migrated outward. Such a process arises for example when two planets are in mean motion resonance and evolve within
a common gap. Provided that the inner planet is more massive than the outer one and that the disc viscosity is small enough, Crida et al. (2009a)
have shown that a pair of planets formed in the 5-20 AU region may reach distances as large as 100 AU through outward migration.

A single planet can also reach large semi major axes due to the action of dynamical corotation torques (Paardekooper 20114), which
operate in protoplanetary discs with low viscosity and a few times more
massive than the Minimum Mass Solar Nebula (hereafter MMSN; Hayashi 1981). The dynamical corotation torque results from the torque exerted by the region of
the disc
that is bound to the planet plus the torque due to the gas material that flows across the orbit as the planet migrates. In isothermal
discs, the dynamical torque scales with the vortensity gradient inside the disc and the planet drift rate (Paardekooper 2014), with the possibility of
runaway migration in the case where the planet migrates in the direction set by the corotation torque. A similar effect can arise in
non-isothermal discs, where the entropy-related corotation torque is the main agent for driving outward migration. Provided the latter is strong enough
to allow for outward migration, the region trapped to the planet tends to be gas depleted and the dynamical torque restricts to the orbit-crossing torque which
has a positive feedback on migration. This can make the planet reaching high drift rates and eventually migrating through zero-torque radii where
analytical formulae for the disc torque (Paardekooper et al. 2011) predict stalling of migration. For example, Pierens (2015) found that a Neptune-mass
planet embedded in a 5 MMSN radiative disc can migrate from 5 AU up to AU by this process, while analytical estimations would predict
stopping of migration at AU.

In this paper, we extend our previous study by including the effect of gas accretion onto the planet. The aims are to i) investigate how
the effects of dynamical torques on the planet orbital evolution are impacted by this process and ii) examine whether massive planets evolving on wide orbits can be
formed by the combined effect of dynamical torques and gas accretion. With respect to our previous work, we also consider accreting disc models
with constant mass flow through the disc (Bitsch et al. 2014), in order to estimate the stages of disc evolution during which dynamical torques can be important.
Our results suggest that gas accretion has two opposite effects on dynamical torques. It tends to make runaway outward migration easier due to the
onset of non-linear effects, but planet growth ultimately leads to a reversal of migration due to the continuous increase of the Lindblad torque and change of the
streamline topology near the planet. This makes it difficult for massive planets at large distances from their star to be formed by this process.

This paper is organized as follows. In Sect. 1, we present the numerical setup. In Sect. 2, we describe the structure of the accreting
disc models we employ. In Sect. 3, we present the results of simulations for non-accreting cores and discuss the
effects of gas accretion in Sect. 4. Finally, we draw our conclusions in Sect. 5.

## 2 Model set-up

### 2.1 Numerical method

Hydrodynamical simulations were performed using the GENESIS numerical code which solves the governing equations for the disc evolution on a polar grid. This code employs an advection scheme based on the monotonic transport algorithm (Van Leer 1977) and includes the FARGO algorithm (Masset 2000) to avoid time step limitation due to the Keplerian velocity at the inner edge of the disc. The energy equation that is implemented in the code includes the effect of viscous heating, stellar irradiation, and radiative cooling. It reads:

(1) |

where is the thermal energy density, the velocity, the adiabatic index which is set to , and is the disc scale height. In the previous equation, is the viscous heating term, and is the local radiative cooling from the disc surfaces, where is the Stephan-Boltzmann constant and the effective temperature which is given by (Menou & Goodman 2004):

(2) |

Here, is the midplane temperature and is the vertical optical depth, where is the surface density and the Rosseland mean opacity which is taken from Bell & Lin (1994) and computed assuming a constant solid/gas ratio . is the irradiation temperature which is computed from the irradiation flux (Menou & Goodman 2004):

(3) |

where is the disc albedo and the stellar luminosity. The factor is set to
be (Chiang & Goldreich 1997), which implies that self-shadowing effects are not taken into account in this study.
Provided that an embedded planet does not carve a deep gap in the disc, test simulations have shown that this leads to temperature and surface density structure of accreting disc models that are in very good agreement with
the 3D hydrodynamical models of Bitsch et al. (2014).

In Eq. 1, is the radiative flux which is treated in the flux-limited diffusion approach and which
reads (e.g. Kley & Crida 2008):

(4) |

where is the midplane density and where is a flux-limiter (e.g. Kley 1989).

In this work, we do not take into account the effects of the disc self-gravity. For the most massive discs we consider, this could artificially alter the torque felt by a migrating planet. In the case where self-gravity is not included, disc gravity indeed causes a change in the angular velocity of the planet only, resulting in a spurious shift of Lindblad resonances (Pierens & Huré 2005). To force both the disc and the planet to evolve in the same gravitational potential, we therefore subtract the axisymmetric component of the disc gravity prior to evaluating the disc torque exerted on the planet. Such a procedure results in migration rates that are consistent with fully self-gravitating calculations (Baruteau & Masset 2008), and has been already employed in previous simulations of migrating planets undergoing gas accretion (D’angelo & Lubow 2008).

When calculating the torque exerted by
the disc on the planet, we also exclude the material contained within a distance from the planet, where is the
planet Hill radius. Such a value is consistent with the one recommended by Crida et al. (2009b) to obtain
an accurate estimation of the migration rate. Moreover,
the potential of the planet is smoothed over a distance , where is the disc scale height at the location of the planet.
Using leads to linear torques that are in good agreement with those derived from a 3D linear theory (Paardekooper et al. 2010). Also,
we note that it has been shown that for planets with planet-to-star mass ratios similar to those considered here, 2D simulations that
employ such a value for the softening length provide a correct estimation of the half-width of the horseshoe region derived
from 3D calculations (Fung et al. 2015), and therefore reasonably estimate the non-linear horseshoe drag for such planets.

Gas accretion onto the protoplanet is modelled by removing a fraction of the gas located within half of the Hill sphere during each time step,
and by adding this mass to the planet. The gas removal rate is chosen such that the accretion timescale onto the planet is , where
is the orbital period of the planet and is a free parameter for which we considered values of
. The case corresponds to the maximum rate at which the planet can accrete gas material from the disc (Kley 1999).

The computational domain is covered by radial grid cells uniformly distributed between and AU typically, and azimuthal grid cells. For simulations in which gas accretion is included, the disc extends up to AU and we employ a number of radial grid cells twice as large.

### 2.2 Initial conditions

We consider accreting disc models with constant mass flow through the disc , with given by:

(5) |

where is the viscosity. Here, the effective kinematic viscosity is modeled using the standard alpha prescription , with and where
is the disc aspect ratio. In the outer parts of the disc where stellar heating dominates over viscous heating, we expect
(Chiang & Goldreich 1997) and this leads to, using Eq. 5, for a disc with constant
(Bitsch et al. 2014). These profiles for the aspect ratio and surface density of the disc correspond to the ones that are adopted as initial conditions in our simulations.

The surface density at is determined according to the chosen , for which we consider values of . For our reference disc with , at AU. For this value of the accretion rate, the stellar luminosity is set to where is the solar luminosity. However, to account for the change in luminosity due to stellar evolution, we follow Bitsch et al. (2015) and consider the stellar luminosity to depend on the value for we impose. The value for as a function of is taken from Bitsch et al. (2015).

### 2.3 Boundary conditions

To obtain a disc model with steady accretion flow through the disc, we allow inflow of disc material at the inner boundary while we impose the value for the accretion rate at the outer edge of the computational domain (Bitsch et al. 2015). This is done by employing a wave-killing zone at the outer boundary (De Val-Borro et al. 2006) where both the surface density and the radial velocity are relaxed toward their initial value (Kley & Haghighipour 2014). Since our initial condition for the surface density is such that , this implies that must be large enough so that the disc structure in the wave-killing region is indeed determined by stellar irradiation.

## 3 Structure of the accreting disc models

The upper panel of Fig. 1 shows the disc aspect ratio as a function of radius for the different values of we consider.
Not surprisingly, the aspect ratio increases with in the inner parts of the disc where viscous heating dominates over
stellar irradiation. This is simply a consequence of the viscous heating rate scaling linearly with . In the outer parts of the disc, however,
stellar irradiation is the main source of heating and this leads to a flared structure with and which is almost
independent of the value for . As expected, the transition
from the viscous heating dominated regime to the stellar heating dominated regime occurs closer to the central star for low . The same holds for the bump in H/R at AU that results from an opacity transition. Both the
location and amplitude of the bump get smaller as decreases, consistent with the results of Bitsch et al. (2014).

We plot in the lower panel of Fig. 1 the surface density profile of the disc for each value of . As the disc tries to
maintain a constant accretion flow with , changes in the disc aspect ratio are correlated with variations in the disc surface density (Bitsch et al. 2014). Because the accretion flow is constant throughout the disc, the strong positive gradient in in the very inner parts of the disc is compensated
by a strong negative gradient in . Also, it is clear by comparing the upper and lower panels of Fig. 1 that the bumps in in the inner disc correspond to jumps in
the surface density. Regarding the outer flared parts of the disc, these give rise to a much shallower surface density
since there. Again, these results are in perfect agreement with the results of Bitsch et al. (2014).

In Fig. 2 we present for each value of the corresponding migration map, that shows the disc torque as a function of planet to star mass ratio and orbital radius. These were calculated using the analytical formulae of Paardekooper et al. (2011) for Type I migration that include saturation effects for the corotation torque but assume that the planet evolves on fixed circular orbit. Possibly strong dynamical torques are therefore not accounted for when computing these migration maps. Because a positive gradient in tends to give rise to a positive entropy gradient and consequently to a negative corotation torque, Type I migration is directed inward in the outer parts of the disc where stellar irradiation is the main source of heating. In the inner parts of the disc where viscous heating dominates, however, migration can be directed outward due to a positive corotation torque resulting from a negative entropy gradient. Looking at the lower panel of Fig. 1, it is clear that assuming , the surface density index in the region where outward migration can occur is such that , resulting in a positive vorticity-related horseshoe drag. However, we find that this component of the corotation torque is typically smaller than the one related to the entropy gradient by a factor of , which is therefore the main engine for driving outward migration. Of course, outward migration can only occur provided that the corotation torque is close to its optimal value, which is reached when the viscous/thermal timescale across the horseshoe region is approximately equal to half the libration timescale (Baruteau & Masset 2013). At thermodynamical equilibrium, where (Masset et al. 2006) denotes the half-width of the horseshoe region, with the planet semi major axis and the disc aspect ratio at the planet location. Given that , where is the planet angular velocity, it is straightforward to show that planets that can experience outward migration have planet-to-star mass ratios such that:

(6) |

For and this gives , consistent with the migration maps shown in Fig. 2. Because
the corotation torque for planets with mass ratios given by Eq. 6 is close to its fully unsaturated value, these are good candidates to experience
dynamical torques provided the disc is massive enough. For the disc models we consider here, it is also interesting to note that this mass
range corresponds to the typical mass range for runaway gas accretion in the classical core accretion
scenario (Pollack et al. 1996). We will discuss the effect of gas accretion on dynamical torques in Sect. 5.

In a previous study (Pierens 2015), we showed that a necessary condition for dynamical torques to come into action is that the Toomre parameter at the planet location satisfies the following relation:

(7) |

where is the dimensionless torque that takes into account possible saturation effects of the corotation torque, and , with the disc surface density at the planet location. When the condition given by Eq. 7 is met during the course of migration, the drift timescale across the horseshoe region becomes shorter than the libration timescale, with the implication that the horseshoe region contracts into a tadpole-like region. Fig. 3 shows the Toomre parameter as a function of radius for each value of we consider. As expected, decreases with radius and can reach values below the gravitational stability limit for our most massive disc, with accretion rates . Comparison with Fig. 1, however, indicates that the outward-migrating region that we will focus on in this study is gravitationally stable for all values of .

From the above discussion, we expect the effect of dynamical torques to be stronger for large values of , i.e., for more massive discs. This means that the influence of dynamical torques should be most pronounced early in a disc’s evolution. There are reasons for this: i) decreases as increases in the outer regions of our accretion disc models, and ii) the outward-migrating region extends outward as increases (Bitsch et al. 2014), so that an outward migrating protoplanet is more prone to reaching disc regions where the constraint given by Eq. 7 is satisfied if is high enough.

## 4 Dynamical torques on giant planet cores

In this section we provide a quantitative
analysis of the dependence of dynamical torques experienced by protoplanets on , with the goal of estimating the critical accretion rate
above which dynamical torques strongly affect the orbital evolution of protoplanets. We focus on protoplanets
with masses of 20-30 , since these are more sensitive to dynamical torques due to their corotation torque being
close to their fully unsaturated value. We neglect gas accretion in this section but will consider this process in the next section.

These planets are released immediately after insertion in the disc and the initial planetary orbital radius is located in the region where outward migration is expected to occur from the migration maps in Fig. 2. It takes different
values that depend on the value for that is considered. It is set to AU for , AU for
, and AU for . For each , we plot in the upper panel of
Fig. 4 the evolution of the semimajor axis relative to of a planet as a function of time (measured in units of planet orbits at ). The corresponding drift rates versus normalized semi major axis are shown in the lower panel of Fig. 4, and
compared with the migration rate expected from static torques only. The migration rate assuming a static torque has been simply obtained by measuring, at different
orbital radii, the torque experienced by a planet evolving on a fixed circular orbit. We see that for , there is decent agreement between the actual migration rates of the planet and those derived assuming a static torque. For
, however,
these two estimations can significantly differ from each other, with the actual planet drift rates being larger than those obtained from static torques
by a factor of typically. For , both estimates show almost similar values initially, but these diverge at later times. This is a clear indication that in this case, dynamical torques can play an important role in the planet’s orbital evolution. As the
planet migrates outward, the Toomre parameter at the planet location continuously decreases until a point in time where the condition given
by Eq. 7 becomes eventually satisfied, resulting in the planet horseshoe region to no longer extend to the full in azimuth.
This is illustrated in Fig. 5, which shows contours of
the perturbed surface density at (as well as a few streamlines) for the case with . These clearly show that the horseshoe has contracted into a tadpole-like region (dark streamline). Since this region moves with the planet as it migrates, it has a negative feedback on migration, whereas disc material flowing across the planet orbit (represented as blue streamline) has a positive feedback on migration. From Fig. 5,
however, it is clear that the
tadpole-like region is underdense relative to the ambient disc, which arises because the planet migrates in the direction set by the
entropy-related corotation torque. As discussed in Pierens (2015), this results in a net positive feedback on migration and the possibility for the planet
to reach high drift rates, consistently with what is observed in the lower panel of Fig. 4.

The orbital evolution of planets embedded in our four accreting discs is presented in Fig. 6. Again, we find strong
dynamical torques provided that . For , migration is stopped close
to the zero-torque radius where the
differential Lindblad torque counterbalances the saturated corotation torque. For higher accretion rates, however, the planet is observed to migrate beyond the zero-torque radius and the transition between the viscous heating dominated and stellar
heating dominated regimes. In the run with for example, outward migration stalls at
( AU) whereas
static torques predict inward migration for ( AU). Moreover, the transition between the viscous heating dominated and stellar
heating dominated regimes is located at ( AU) for this model. For , outward migration is expected to stop at
(AU) whereas dynamical torques make the outermost location that can be reached by the planet be
( AU). Again, this is outside the transition between the viscous heating dominated and stellar
heating dominated regimes which is located at ( AU) for this run.

## 5 Effect of gas accretion

### 5.1 Conditions for triggering fast outward migration

As mentioned above, are expected to undergo runaway gas accretion in the standard core accretion scenario (Pollack et al. 1996). In this section, we examine how gas accretion affects the orbital evolution of a planet that is subject to strong dynamical corotation torques. Fig. 7 presents the results of simulations of the dynamical evolution of an accreting protoplanet with initial mass of , and that employ the routine for gas accretion described in Sect. 2.2.

Fig. 7 shows that gas accretion strengthens dynamical torques at early times. We indeed observe that an
accreting protoplanet can experience large excursions in the outer disc, depending on the accretion rate and the accretion parameter
. For example, in the run with and , an accreting core starting at AU
rapidly migrates outward to AU under the combined effects of dynamical torques and gas accretion. In presence of
gas accretion, there is also evidence for runaway outward migration in accretion discs models with whereas runs in which gas accretion was neglected showed fast outward migration for only.

It is plausible that the positive effect of gas accretion on the strength of dynamical torques is related to the onset of non-linearities as the planet
mass ratio approaches . Results from previous work
(Paardekooper 2014; Pierens 2015) indeed indicate that dynamical torques can easily make intermediate mass planets with enter the
outward runaway migration regime. This occurs because such planets tend to have a larger horseshoe region than low-mass cores (Masset et al. 2006),
resulting in a boost of the corotation torque associated with the gas material that flows across the planet orbit.

As discussed in Pierens (2015), fast outward migration is triggered whenever a planet’s drift rate is
faster than the critical drift rate , such that the planet migrates a distance greater than the horseshoe
width over one horseshoe libration time. The migration rate of the planet is given by:

(8) |

where corresponds to the mass of gas material contained inside the planet orbit. For low-mass planets, the criterion for runaway migration given by Eq. 7 can be recovered by simply substituting an estimation for the half-width of the horseshoe region in the condition . For higher mass planets with (Masset et al. 2006), however, it can be shown that the condition is equivalent to:

(9) |

Interestingly, we see that the condition for the planet to enter the fast migration regime now depends on the planet mass ratio and that it is
easier to get runaway migration as the planet grows. In principle, the term also depends on the planet
mass ratio since the disc surface density profile is expected to be altered if the planet mass is high enough. In a recent work,
however, Lega et al. (2015) found that even for planets that are subject to significant non-linear effects, there is still quite a good agreement between the torque obtained from 3D hydrodynamical simulations and the one
predicted by existing analytical formula. This indicates that in accreting disc models, the term is almost independent of for intermediate planet masses with . We also note that the condition given by Eq. 9 implicitly depends on the softening
parameter that is adopted for the planet potential. Since the entropy-related horseshoe drag scales as (e.g. Baruteau
& Masset 2013), more massive
discs (or equivalently higher accretion rates) would be required to trigger runaway migration for values of larger than the one
considered here.

For a given value of the Toomre parameter , Eq. 9 provides an estimation of the minimum value for the mass ratio above which fast
migration is expected to occur. It is given by:

(10) |

Considering the model with , the middle panel of Fig. 7 indicates that runaway migration is triggered once the planet is located at AU and its mass ratio has reached . Assuming that the corotation torque is unsaturated we have (Paardekooper et al. 2010):

(11) |

where and are the negative of the power-law indices associated with the temperature and entropy profiles respectively. For
, we obtain , , at AU and consequently . Moreover, examination of Fig. 1 shows that and at AU. Inserting these values
into Eq. 10 leads to
, which is in reasonable agreement with the results of the calculation corresponding to .

Fig. 7 shows that the ability of an accreting core to undergo runaway outward migration is a sensitive function of the accretion parameter . In the case where the growth rate is maximum, namely for , the resulting strong Lindblad torque rapidly leads to a reversal of the migration direction. In the opposite limit of slow gas accretion, however, corotation torque saturation can prevent the planet from entering the fast migration regime. A crude estimation of the minimum value for that is required to avoid saturation of the corotation torque as the planet accretes gas material from the disc can be obtained from the argument that the doubling-mass timescale needs to be shorter than the libration timescale . For our accretion procedure (see Sect . 2.2), the doubling mass timescale is given by:

(12) |

Substituting this previous expression in the criterion and using , we find that a necessary condition for the planet to avoid the effect of corotation torque saturation as it grows is given by:

(13) |

If we consider our nominal disc model with , at AU, which leads to the prediction that a protoplanet with initial mass ratio will undergo strong dynamical torques provided that . This is slightly smaller than what we report from the simulations, which resulted in fast outward migration in the runs with , as can be observed in Fig. 7.

### 5.2 Long term evolution

Although gas accretion can, under certain circumstances, boost the effect of dynamical torques, it appears from
Fig. 7 that an episode of fast outward migration
cannot be sustained over long timescales. This arises because of the following effects:

i) The negative Lindblad torque increases in strength with the planet mass. Because the mass contained
within a planet’s feeding zone scales as with (see Sect. 2.2) in the outer parts of the disc, outward migration promotes rapid growth of the planet, such that the Lindblad torque
comes to dominate dynamical torques. We note that for the highest accretion rates that are considered, the mass growth
at large planet masses should be in reality reduced due to a high accretion luminosity (Stamatellos 2015).

ii) The trapped underdense region cools due to thermal diffusion, resulting in a slight increase in surface density. The
implication is that the negative feedback effect from that region on migration is enhanced as the planet migrates outward.

iii) Part of the surface density deficit within the trapped region is progressively
lost as outward migration proceeds. This is illustrated in Fig. 8 where we display contour plots of the perturbed surface density at different times for the run with and
,. At early
times, the streamlines that are overplotted show gas flowing across the planet’s orbit, together with the underdense librating material that is
responsible for the strong dynamical torques experienced by the planet. At yr, however, we see that the gas material that is bound to the planet
now librates around the Lagrange point, and this tends to remove the surface density deficit in the underdense region. Interestingly, it has been shown that a similar effect is responsible for stalling
of runaway migration in the context of the classical Type III migration regime (Masset & Papaloizou 2003;
Peplinski et al. 2008; Lin & Papaloizou 2010).

The third panel of Fig. 8 shows the perturbed surface density at yr, just after the direction of migration has been reversed. At that time, the planet mass ratio is (see Fig. 7) and its
semi major axis AU. Because the planet evolves in a region of the disc whose thermodynamical state is close to the isothermal limit, the trapped region which is now located at the leading side of the planet has almost the same surface density as the ambient disc. Moreover, we see that the
planet has not opened a gap at that time. The gap opening criterion of Crida et al. (2006) is given by:

(14) |

which is not satisfied for due to the high value for the aspect ratio at the planet’s location. Gap clearing
is therefore not expected for , consistent with what we observe in the simulation. This is further illustrated in
Fig. 9, which shows the disc surface density and aspect ratio profiles at different times for the same run .

The planet subsequently migrates back inward because of the strong tidal torque resulting from the large planet and inner disc masses. At
yr (the fourth panel of Fig. 8) the planet has opened a gap in the disc,
which is confirmed by looking at the corresponding disc surface density and aspect ratio plots in Fig. 9. We note that
these profiles should be taken with care since self-shadowing effects, which are not included in this study (see Sect. 2.1 ), are expected to be important
once the planet opens a deep gap in the disc. From this time onward,
Fig. 7 shows that the planet slowly migrates inward, with a drift rate corresponding to the Type II regime. We note that
such a mode
of evolution is similar to what is seen in the other runs where the planet experienced a large excursion in the outer disc.

Compared to the prediction from the criterion given by Eq. 14, a common feature of these runs is that we find
that gap opening occurs at slightly larger masses. For the reference run with and , for example, the planet is actually found
to open a gap once its mass ratio has reached whereas Eq. 14 predicts that this should occur for
. As demonstrated by Malik et al. (2015), such a discrepancy is due to the fast inward migration of the planet and arises because the planet migrates across the region
where it is expected to open a gap on a timescale shorter than the gap opening timescale.

## 6 Discussion and conclusion

In this paper, we have presented the results of hydrodynamical simulations of accreting protoplanets embedded in radiative
discs. The goal of this study was to examine the influence of gas accretion onto the orbital evolution of outward
migrating planets subject to strong dynamical corotation torques. In relatively massive protoplanetary
discs, dynamical torques can make a low/intermediate mass planet enter a fast migration regime and experience large excursions in the outer disc.

We considered accreting disc models with constant mass flow through the disc and where radiative cooling is balanced by
viscous and stellar heating, and studied how dynamical torques depend on the accretion rate . For an
viscous stress parameter , we find evidence of strong dynamical torques
for . In that case, cores with fixed mass of 20-30 Earth masses migrate outward with drift rates that are significantly larger than those
predicted from analytical formulae for the disc torque, and are subsequently trapped at or just beyond the predicted zero-torque radius where the
Lindblad and corotation torques cancel each other. Protoplanets with initial mass of and undergoing runaway gas accretion, however, are
found to enter a runaway outward migration regime under some circumstances. Depending on the values for the mass flux through the
disc and the accretion timescale onto the planet, protoplanets formed in the 5-10 AU region may migrate well beyond the the zero-torque radius out to separations of .

We find that gas accretion has two opposite effects. It is easier for an accreting protoplanet to enter a fast migration regime due to the
onset of non-linearities as its mass ratio approaches , where is the disc aspect ratio at the planet location. Planet growth, however, is accompanied by an increase of the Lindblad torque which
tends to counteract the effect of dynamical torques for an outward migrating planet. This, combined with the fact that the structure of the region that
is bound to the planet is significantly altered as a result of mass accretion and diffusion processes, systematically leads to a change in the direction
of migration. Interestingly, it is worthwhile to notice that similar results have been reported from simulations of accreting planets undergoing outward-directed Type III
migration in isothermal discs (Peplinski et al. 2008).

Reversal of migration typically occurs once the planet mass has reached . Due to the high value for the disc aspect ratio at
the planet location, a Jupiter mass
planet that stops its outward migration does not open a gap in the outer parts of our accreting disc models. This resuts in a
very fast inward migration due to the large planet and inner disc masses, which further prevents the planet from opening a gap. Such a
process can however arise at later times, once the planet is massive enough for the gap opening timescale to be shorter than the migration
timescale. This corresponds to the transition to Type II migration, which typically occurs in our simulations at planet masses of and orbital distances of AU.

At first glance, its therefore seems difficult for such a mechanism to explain the presence of massive planets on wide orbits, and the population of giant planet residing beyond a few AU as well. The work of Coleman & Nelson (2014) indeed suggests that for a Jovian planet to survive at a few AU, it must have initiated runaway gas accretion and Type II migration beyond AU typically, and late in the disc lifetime (see also Bitsch et al 2015b). This is because giant planet cores that accrete gas tend to undergo rapid inward migration due to corotation torque saturation before they open a gap. It is clear that the mechanism presented here could potentially allow a growing core that undergoes runaway accretion to reach large distances in the outer disc and therefore avoid corotation torque saturation. Our results indicate that this could lead to the formation of Jupiter-mass planets at orbiting at AU. However, since this occurs for typical accretion rates higher than only, it would seem difficult for such planets to survive type II migration before the disc is dispersed through photo evaporation and viscous evolution.

It should be noted, however, that the simulations presented here contain a number of simplifications. Although
dynamical corotation torques play an important role on the planet orbital evolution for typical values of the
Toomre stability parameter of , we did not consider
the effect of disc self-gravity. Self-gravity is expected to cause a shift of Lindblad resonances (Pierens & Huré 2005; Baruteau & Masset 2008),
resulting in a slightly stronger Lindbad torque in comparison with the case where
self-gravity is not considered. For highest accretion rates we considered, can even fall below the gravitational stability limit
in the very outer parts of the disc, possibly resulting in gravitoturbulence if the gas cooling timescale is not too short (Baruteau et al. 2011). Moreover,
we have only considered the single planet case. In the multiple planet case, however, planet-planet interactions may prevent the fast inward migration
episode that is observed in the simulations. We will examine this issue in a future paper.

## Acknowledgments

Computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Universite de Bordeaux and by HPC resources of Cines under the allocation c2016046957 made by GENCI (Grand Equipement National de Calcul Intensif). We thank the Agence Nationale pour la Recherche under grant ANR-13-BS05-0003 (MOJO).

## References

- Baruteau & Masset (2008) Baruteau, C., & Masset, F. 2008, ApJ, 672, 1054
- Baruteau et al. (2011) Baruteau, C., Meru, F., & Paardekooper, S.-J. 2011, MNRAS, 416, 1971
- Baruteau & Masset (2013) Baruteau, C., & Masset, F. 2013, Lecture Notes in Physics, Berlin Springer Verlag, 861, 201
- Bell & Lin (1994) Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987
- Bitsch et al. (2014) Bitsch, B., Morbidelli, A., Lega, E., & Crida, A. 2014, A&A, 564, A135
- Bitsch et al. (2015) Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28
- Bitsch et al. (2015) Bitsch, B., Lambrechts, M., & Johansen, A. 2015b, A& A, 582, A112
- Brandt et al. (2014) Brandt, T. D., McElwain, M. W., Turner, E. L., et al. 2014, ApJ, 794, 159
- Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368
- Coleman & Nelson (2014) Coleman, G. A. L., & Nelson, R. P. 2014, MNRAS, 445, 479
- Crida et al. (2006) Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587
- Crida et al. (2009) Crida, A., Masset, F., & Morbidelli, A. 2009a, ApJL, 705, L148
- Crida et al. (2009) Crida, A., Baruteau, C., Kley, W., & Masset, F. 2009b, A&A, 502, 679
- Cumming et al. (2008) Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531
- D’Angelo & Lubow (2008) D’Angelo, G., & Lubow, S. H. 2008, ApJ, 685, 560
- de Val-Borro et al. (2006) de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529
- Dodson-Robinson et al. (2009) Dodson-Robinson, S. E., Veras, D., Ford, E. B., & Beichman, C. A. 2009, ApJ, 707, 79
- Fung et al. (2015) Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101
- Goździewski & Migaszewski (2014) Goździewski, K., & Migaszewski, C. 2014, MNRAS, 440, 3140
- Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
- Kley (1989) Kley, W. 1989, A&A, 208, 98
- Kley (1999) Kley, W. 1999, MNRAS, 303, 696
- Kley & Crida (2008) Kley, W., & Crida, A. 2008, A&A, 487, L9
- Kley & Haghighipour (2014) Kley, W., & Haghighipour, N. 2014, A&A, 564, A72
- Lambrechts et al. (2014) Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35
- Lega et al. (2015) Lega, E., Morbidelli, A., Bitsch, B., Crida, A., & Szulágyi, J. 2015, MNRAS, 452, 1717
- Lin & Papaloizou (2010) Lin, M.-K., & Papaloizou, J. C. B. 2010, MNRAS, 405, 1473
- Maire et al. (2015) Maire, A.-L., Skemer, A. J., Hinz, P. M., et al. 2015, A&A, 576, A133
- Malik et al. (2015) Malik, M., Meru, F., Mayer, L., & Meyer, M. 2015, ApJ, 802, 56
- Masset (2000) Masset, F. 2000, A&AS, 141, 165
- Masset & Papaloizou (2003) Masset, F. S., & Papaloizou, J. C. B. 2003, ApJ, 588, 494
- Masset et al. (2006) Masset, F. S., D’Angelo, G., & Kley, W. 2006, ApJ, 652, 730
- Marois et al. (2008) Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348
- Marois et al. (2010) Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010, Nature, 468, 1080
- Mayor et al. (2011) Mayor, M., Marmier, M., Lovis, C., et al. 2011, arXiv:1109.2497
- Menou & Goodman (2004) Menou, K., & Goodman, J. 2004, ApJ, 606, 520
- Paardekooper et al. (2010) Paardekooper, S.-J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950
- Paardekooper et al. (2011) Paardekooper, S.-J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293
- Paardekooper (2014) Paardekooper, S.-J. 2014, MNRAS, 444, 2031
- Pepliński et al. (2008) Pepliński, A., Artymowicz, P., & Mellema, G. 2008, MNRAS, 386, 179
- Pierens & Huré (2005) Pierens, A., & Huré, J.-M. 2005, A&A, 433, L37
- Pierens (2015) Pierens, A. 2015, MNRAS, 454, 2003
- Shvartzvald et al. (2016) Shvartzvald, Y., Maoz, D., Udalski, A., et al. 2016, MNRAS, 457, 4089
- Stamatellos (2015) Stamatellos, D. 2015, ApJL, 810, L11
- Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
- van Leer (1977) van Leer, B. 1977, Journal of Computational Physics, 23, 276