Low mass planet migration in magnetically torqued dead zones – II. Flow-locked and runaway migration, and a torque prescription
We examine the migration of low mass planets in laminar protoplanetary discs, threaded by large scale magnetic fields in the dead zone that drive radial gas flows. As shown in Paper I, a dynamical corotation torque arises due to the flow-induced asymmetric distortion of the corotation region and the evolving vortensity contrast between the librating horseshoe material and background disc flow. Using simulations of laminar torqued discs containing migrating planets, we demonstrate the existence of the four distinct migration regimes predicted in Paper I. In two regimes, the migration is approximately locked to the inward or outward radial gas flow, and in the other regimes the planet undergoes outward runaway migration that eventually settles to fast steady migration. In addition, we demonstrate torque and migration reversals induced by midplane magnetic stresses, with a bifurcation dependent on the disc surface density. We develop a model for fast migration, and show why the outward runaway saturates to a steady speed, and examine phenomenologically its termination due to changing local disc conditions. We also develop an analytical model for the corotation torque at late times that includes viscosity, for application to discs that sustain modest turbulence. Finally, we use the simulation results to develop torque prescriptions for inclusion in population synthesis models of planet formation.
keywords:planet–disc interactions – protoplanetary discs – planets and satellites: dynamical evolution and stability
As a planet grows in a protoplanetary disc, angular momentum exchange between the planet and disc due to tidal torques begin to become important in driving migration when it exceeds the mass of Mars and approaches the Earth’s mass. The theory of low mass planet migration has explicitly or implicitly been formulated in the context of viscous discs, by assuming that the underlying disc structure near the planet is smooth and unperturbed, and it has been supposed that the effective viscosity would be supplied by turbulence, which would also provide the accretion-driving stresses in the disc to account for mass flow onto the central star.
Our current understanding of the non-ideal magnetohydrodynamics of protoplanetary discs, however, suggests that they are likely not sufficiently turbulent to explain observed accretion rates, and are essentially laminar over large regions that incorporate areas of the disc traditionally associated with the sites of planet formation (i.e. between –10 AU). The picture that emerges instead is one where disc accretion is driven by magnetised winds that are launched from thin ionised surface layers, and possibly also by large scale horizontal magnetic fields near the midplane, with the disc remaining laminar throughout the region described above (Bai & Stone, 2013; Bai, 2013, 2014a, 2014b; Lesur et al., 2014; Gressel et al., 2015; Bai, 2015; Simon et al., 2015a; Bai et al., 2016; Bai, 2016; Béthune et al., 2017). In those discs where the vertical magnetic field is parallel to the angular momentum vector of the disc, the Hall effect, acting in the Ohmic dead zone (Wardle & Ng, 1999; Pandey & Wardle, 2008) can generate large scale radial and azimuthal fields, and a radial Maxwell stress, through the Hall-shear instability acting in concert with the Keplerian shear (Kunz, 2008; Kunz & Lesur, 2013), even in a laminar flow (Bai & Stone, 2013; Lesur et al., 2014; Béthune et al., 2017). In this paper, which follows on from Paper I, we aim to derive the properties of planetary migration torques in such discs, where we continue to focus on the Hall-modified region that drives an accretion flow throughout the vertical column of the disc and not just in its surface layers.
For planets which do not open a gap in the disc (i.e. the Type-I migration regime)111Although all planets eventually alter the surface density profile of the disc to some extent, the approximation that the surface density profile is unperturbed is a useful analytical approximation in understanding this limit., the migration torque consists of two basic components: the Lindblad torque, that arises due to the spiral wake driven by the planet, and the corotation torque, that originates from material close to the planet’s orbit which undergoes horseshoe turns (or U-turns) in front of and behind the planet. The Lindblad torque (Goldreich & Tremaine, 1979) is insensitive to the presence or absence of viscous accretion stresses in the disc, and drives inward migration for most disc models. In the absence of other effects, the migration can be rapid, with migration times at 1 AU being years for an Earth mass planet and years for a 10 Earth mass body.
The corotation torque (e.g. Goldreich & Tremaine, 1979; Ward, 1991), however, has much richer dependencies than the Lindblad torque, that in turn lead to a plethora of migration phenomena, some of which are the focus of this paper. The corotation torque can be viewed as resulting from the exchange of angular momentum with gas parcels undergoing horseshoe turns when encountering the planet. Clearly, there must be an asymmetry between the fluid elements making the turns in front of and behind the planet in order for a net torque to arise, and this is provided by a radial vortensity gradient in the disc, and perhaps also an entropy gradient if present (Masset, 2001; Paardekooper & Mellema, 2008; Baruteau & Masset, 2008a). Furthermore, phase-mixing of material trapped on librating horseshoe orbits can remove the gradients that give rise to a net corotation torque, causing it to saturate and switch off. Viscous diffusion, however, can reestablish the vortensity gradient by promoting mixing between the corotation region and the surrounding disc, and hence can enable the corotation torque to be sustained. Under favourable disc conditions, with negative entropy and vortensity gradients and appropriate levels of viscous and radiative diffusion, a strong and positive corotation torque can be maintained that balances or exceeds the Lindblad torque, and hence prevents rapid migration of a planet into the central star (Paardekooper et al., 2010, 2011; Masset & Casoli, 2010). We note that this corotation torque can be estimated at any point in time according to the instantaneous conditions in the disc at the location of the planet, and does not depend on the motion of the planet or on the history of the torques that have been applied to the planet or the disc. As such, we can describe it as a static corotation torque.
A corotation torque can arise in an inviscid disc, but in this case the motion of the planet relative to the disc gas must be taken into account, resulting in a dynamical corotation torque whose magnitude depends on the planet’s motion (Paardekooper, 2014). Here, the radial motion of the planet with respect to the disc causes distortion of the horseshoe region222In addition to the possibly small asymmetry arising from a vortensity gradient (Casoli & Masset, 2009)., which takes on a teardrop shape, with one wider and one narrower horseshoe turn, and the orbits narrowing in between. On the horseshoe turn where the outermost librating streamline is narrower, some background disc material will flow across the corotation region due to the radial motion of the planet, making only a single horseshoe turn before continuing on as part of the background disc. The mismatch between this material and the librating material, and the asymmetry of the librating streamlines, can lead to a corotation torque even in the absence of a turbulent viscosity to mix the vortensity of the librating region with the surrounding disc. For an inward migrating planet driven by Lindblad torques, this corotation torque can grow with time as the planet migrates and can eventually cause the inward migration to stall (Paardekooper, 2014).
In McNally et al. (2017) (hereafter Paper I) we showed that in the non turbulent, very magnetically dead inner regions of a protoplanetary disc where radial flow of gas may be driven by stresses due to a Hall-generated laminar magnetic field, the resulting corotation torque can be described as a generalisation of the dynamical corotation torque. This generalisation involves accounting for the relative radial motion of the gas and planet. In Paper I, we treated only torques acting on a planet on a fixed circular orbit with the gas flowing past the planet due to large scale magnetic accretion stresses. However, we were able to predict that there should be four distinct regimes of planet migration, that depend on the relative motion of the gas and the planet. In this paper, we introduce live (fully dynamical) planets into our previous models and follow the longer term evolution of the net (Lindblad plus corotation) torque, allowing us to probe the four regimes of the dynamical corotation torque. Solutions where the planet is driven inwards, outwards, and where the net torque reverses are found. As with Paper I, we consider disc models with a simple barotropic isothermal equation of state, and hence we do not consider the influence of entropy gradients on the dynamical corotation torque, as have been considered previously by Pierens (2015) and Pierens & Raymond (2016) in the context of viscous disc models. Hence, in this work we retain our focus on purely vortensity-related effects.
In summary, we expect that in all reasonable disc models, an embedded low mass planet will experience Lindblad torques that attempt to drive inwards migration. In a disc where the accretion stresses are primarily provided by turbulence, we expect that the corotation torques will be provided by the static torques described above and in Paardekooper et al. (2011). In a laminar magnetised disc, for which Hall EMFs are ineffective and accretion is driven only in narrow surface layers by the launching of a magnetised wind, such that there are no radial gas flows near the midplane, then we would expect the corotation torque to behave as the dynamical torque discussed in Paardekooper (2014). In a disc where the Hall effect plays an important role, such that accretion flows occur near the midplane and in the surface layers, we would expect the corotation torque to behave according to dynamical torques described in this paper and in Paper I, with the long term evolution having a strong dependence on the relative motion between the planet and disc gas at the point when migration is initiated.
The paper is organised as follows. In Section 2 we discuss the basic analytical properties of the problem, define basic quantities and recapitulate relevant results from Paper I. The first numerical experiments involve enforcing specified initial migration to provoke the four basic regimes of migration described in Paper I in Section 3. We then explore reversing inward migration by magnetic disc torques in Section 4, including a comparison between a magnetically driven inflow and a viscously driven one in Section 5. In Section 6 we analyse the fast migration regime both analytically and numerically, and present an interpolation table for the late-time torque in inviscid discs that we use in Appendix B to construct a migration torque prescription for inclusion in N-body simulations of planet formation and population synthesis models. We propose and test an analytical model for steady-state torques which can exist when a residual viscosity is present in the disc in Section 7 and in Appendix C. We offer discussion of the results and their possible implications for planet formation scenarios in Section 8, and draw our conclusions in Section 9.
2 Fundamentals of Inviscid Type-I Migration
This paper is concerned with the dynamical motion of a planet in an inviscid disc that is subject to a torque that drives a radial laminar flow through the full column density of the disc. As discussed in Paper I, such a torque can arise from horizontal magnetic fields that are generated by the Hall Shear Instability when the Hall effect is included in disc evolution calculations. We work in the limit where Type-I migration theory is valid, such that any change to the local surface density induced by the planet is small, and is not a primary factor in understanding the origin of the migration torque. Furthermore, we assume that the presence of the planet does not alter the background torque that acts on the disc. When considering a torque of magnetic origin, this is equivalent to assuming that the Ohmic resistivity is large enough to rapidly diffuse any magnetic field perturbations that are induced by the planet (as discussed below in Section 2.1).
At the end of Paper I we gave a formula for the corotation torque acting on a low mass planet in such a disc, assuming slow migration (the precise meaning of which will be discussed later) as
where is the characteristic inverse vortensity of the librating streamlines trapped on horseshoe orbits, is the inverse vortensity of the background disc at the planet location, is the disc surface density at the planet position, is the planet’s radial position (such that is the migration speed), is the half-width of the corotation region, is the Keplerian orbital frequency at the planet position, and is the radial velocity of the disc gas flow. The vortensity, which is sometimes called the specific vorticity, can be written as . As the flow is two dimensional, the curl of the two dimensional velocity only has a -component so without ambiguity the vortensity can be treated as a scalar. Even if the migration is fast, in that the planet moves radially quickly with respect to the gas, then we expect the torque to have the same sign as suggested by this form. Areas of parameter space close to where the torque changes sign are in the slow migration regime.
In Paper I we identified four regimes of behaviour in this torque expression separated by sign changes of two quantities. These regimes correspond to qualitative differences in the expected migration behaviour of a planet. We restrict our attention in this paper to Keplerian discs with radially decreasing vortensity profiles. Since vorticity , where is the Keplerian orbital angular velocity, our power-law disc models have where . As equation (1) applies formally when the relative radial velocity of the gas and planet are small and the vortensity contrast is small, these parameters determining the borders between regimes are the disc gas radial flow velocity and the ‘initial’ planet migration velocity close to transitions between regimes. Given this, the four regimes and the expected long-term migration behaviours are:
- (i) and :
The disc accretion flow and planet migration are inwards, and the planet initially migrates faster than the disc flow The planet migrates inwards, close to but slightly faster than the gas inflow speed.
- (ii) and :
The disc accretion flow and planet migration are inwards, and the planet initially migrates slower than the disc flow The planet’s inward migration reverses and runs away outwards.
- (iii) and :
The disc flow is outward and the planet initially migrates inwards The planet’s inward migration reverses but the outward migration speed cannot exceed that of the gas.
- (iv) and :
The disc flow and planet migration are outwards, and the planet is initially migrating faster than the disc flow The planet’s outward migration can run away.
Regimes (i) and (iii) result in slow migration, by which we mean that the planet asymptotically migrates with the radial motion of the gas. Hence, the planet’s motion is locked to the radial gas flow. Regimes (ii) and (iv), however, result in the planet moving fast with respect to the gas. The meaning of fast and slow radial motion is defined in terms of the relevant flow timescale across the planet’s corotation region.
To describe the speed of the radial flow of gas driven by the background torque induced by a laminar magnetic field, we use the ratio of the flow crossing time across the corotation region to the libration time
where is the flushing timescale, and is the libration timescale. In turn, the flushing timescale is
where is the width of the corotation region, and is the radial velocity of the disc gas flow, and the libration timescale is given as
Hence, a positive value of corresponds to a negative background torque driving gas inflow. Finally, in these expressions the half-width of the corotation region can be approximated as
when the disc is two-dimensional and the planet potential is smoothed as a Plummer sphere with smoothing length , with the planet-star mass ratio , being the radial position of the planet, the disc scale height and , following Masset et al. (2006).
The analytical model of Paardekooper (2014) was constructed in the slow migration regime, where the time taken for the planet to migrate across its own corotation regime is long compared to the libration period, that is
In Paper I the analogous condition for slow radial gas flow was (we are slightly more careful here to explicitly write the absolute value bars when discussing this condition). Thus, the unified form equation (1) is only strictly valid when
which gives the earlier promised definition of slow migration. However, freely migrating planets will often migrate at rates which violate this condition, particularly in regime (ii) or regime (iv). We address this fast migration in Section 6. It is also convenient to define an appropriately generalised version of the parameter taking into account the radial motion of the planet as
Again corresponds to slow migration.
A second property that immediately becomes apparent from an examination of Equation (1) is that there exists a bifurcation in the behaviour of planets attempting to change migration from regimes (i) to (ii). This could occur, for example, if a planet is in regime (i) with relatively slow gas inflow, but moves into a region of the disc where the inflow is faster. Or if the gas inflow changes with time and speeds up such that it becomes faster than the planet’s migration speed. The bifurcation occurs because changes to the factor in Equation (1) that occur during the regime (i) phase of migration are reversed as the disc torque and inflow speed of the gas are increased and the planet tries to enters regime (ii). This unwinding of the relative inverse vortensity perturbation causes the corotation torque, , to pass through zero, and at that moment the migration speed of the planet is only determined by the Lindblad (wake) torque. Hence, for a particular value of , a critical surface density separates scenarios where planet migration can transition from regime (i) to regime (ii), since the migration speed at the point of change over depends linearly on the surface density. To be able to make the transition from regime (i) to regime (ii) the Lindblad torque must drive the planet radially inwards slower than the disc gas inflow velocity , or
The migration rate can be expressed in terms of the migration torque as
where is the planet mass, is the gravitational constant, is the mass of the central star, and is a scale used to nondimensionalized the torque. A linear calculation of the Lindblad torque in a more general adiabatic disc yields
for the rest of the dependence where is the adiabatic index of the gas ( for isothermal gas), is the negative power law slope of the radial temperature gradient, and is the negative power law slope of the radial surface density profile of the disc, is the aspect ratio of the disc, and is the smoothing length of the planet potential (for a Plummer sphere softening) (Paardekooper & Papaloizou, 2008). The only dependence of the migration rate on surface density at the planet position is contained in the parameter nondimensionalizing the Lindblad torque
where is the angular frequency of the planet’s orbit. Thus in a globally isothermal disc with a laminar magnetic inflow torque acting on the gas there exists a above which a planet cannot cross from regime (i) to regime (ii). This phenomena will be demonstrated in Section 4, where the critical surface density sets a divide between planets which can reverse their migration and migrate outwards, and those which are driven inwards in the same disc conditions.
2.1 Characteristics of the problem established in Paper I
Here we discuss other aspects of the problem established previously. Importantly, in [sections 2.1.1 and 3.2]2017MNRAS.472.1565M, by means of ionisation chemistry calculations, we found that in the Ohmic resistivity dominated dead zone of a protoplanetary disc, the Ohmic diffusion timescale for the magnetic field is much faster then the U-turn flow timescale for low mass planets. This implies that the deviations from Keplerian flow due to the presence of the planet are not expected to result in significant changes to the magnetic field configuration, and that hence simulations including the live evolution of the magnetic field will be equivalent to one including only a fixed field or simply the equivalent Lorentz force from the static magnetic field given in the initial condition.
As a consequence of this fact, though the models in this work are constructed with a magnetic field strength and a resistivity value, these are degenerate, and the only physical parameter is the radial flow velocity driven by the resulting magnetic field configuration. Significant midplane flow velocities can be included within the total accretion rate typically observed onto protoplanetary disc hosting stars (Paper I). Thus our models are stated in terms of this radial flow velocity, and not the underlying magnetic field strengths and resistivities. This allows our two dimensional models to be compared to three dimensional configurations with the same radial flow velocity.
As in Paper I the models in this work are two dimensional. The corotation region for a low mass planet has been previously found to to have a columnar flow structure in three dimensions (Masset & Benítez-Llambay, 2016). Thus, the results for two dimensional models of the corotation torque should be similar to three dimensional models. In this work we again consider only an isothermal disc to isolate the effects of the vortensity related corotation torques. We further discuss the possible effects of realistic thermodynamics in Section 8.
3 Demonstrating Four Migration Regimes
To demonstrate that the four migration regimes described above and predicted in Paper I exist, in this section we present customised experiments with forced initial planet migration and disc parameters such that reasonably clean examples of all four regimes are produced. In a disc with a fixed radial inflow, we force the initial inward migration rate, either a little slower than the disc inflow, or much faster. This excites regime (i) or regime (ii) behaviour once the planet is released and allowed to migrate according to the disc torque. A second set of experiments with disc outflow excites regime (iii) by forcing initial inward migration, which results in the planet turning around and migrating outwards, with the planet asymptotically migrating at close to the disc gas outflow speed. Regime (iv) is excited by initially specifying rapid outward migration, and outward fast migration is maintained and increases in speed slightly when the planet is released.
We adopt the two dimensional disc model defined in Paper I. This is a globally isothermal disk in cylindrical coordinates with surface density where is a reference radius and is a constant determining the power-law slope of the disc surface density profile. Our disks are always globally isothermal, with the vertical scale height aspect ratio at radius , which also specifies the constant isothermal sound speed . Throughout this paper, when times are given in units of orbits, they are orbits at . All models are presented with a central stellar mass of and gravitational constant . The planet mass is denoted but usually specified in terms of the planet-star mass ratio . Our simulations are performed with FARGO3D (Benítez Llambay & Masset, 2015; Benítez-Llambay & Masset, 2016) in two dimensions on a cylindrical grid with linear spacing in radius, using a body force to drive radial gas flow instead of a magnetic field. The domain extends in azimuth and in radius from to , and the fiducial resolution used throughout is . As in Paper I, modified de Val-Borro et al. (2006) style damping zones at the boundary are used with the inner one having a radial width of and the outer a radial width of . The body force added to the momentum equation is the Lorentz force which would be produced by a spiral magnetic field given by
where is a constant with units of the magnitude of magnetic field, is the Keplerian angular velocity at , is the magnetic permeability of free space, and is the Ohmic resistivity of the disc gas. The velocity considered in calculating this force is the Keplerian velocity, as this is the thin disc approximation used in constructing the spiral magnetic field in Paper I, and the presence of the planet will not significantly alter the magnetic field configuration as demonstrated in Paper I. This force produces a radial flow in the two-dimensional disc with velocity
where , is the midplane density, and we take to be a constant in these models. As in Paper I, using just this force and eliminating the induction equation from the model is possible as the disc has sufficient Ohmic diffusivity that the shape of the magnetic field lines is not significantly altered by the flow perturbations due to the planet. We initially specify the disc gas velocity with , so the model is run for orbits to relax the radial flow to its equilibrium state before adding the planet potential. Self-gravity of the disc gas is neglected, and the axisymmetric component of gravity is thus neglected when calculating the planetary migration torque to ensure consistency between the angular velocity of the planet and the disc gas (Baruteau & Masset, 2008b). The entire disc mass, including that in the Hill sphere (which is not bound to a low mass planet) is included in the torque calculation.
The first experiment which demonstrates regime (i) is given in Figure 1. In this model the disk has surface density and . The disk surface density at is such that the coorbital region mass ratio
is and the planet mass ratio . The body force drives an inflow with speed . For this experiment the parameters must be chosen so that the Lindblad torque driven migration rate is faster than the gas inflow rate, otherwise the planet will eventually enter regime (ii) migration. The planet potential is turned on after orbits, and the migration from is forced until orbits, which is , and the planet is released at times the inwards radial offset that a gas parcel also starting at would have. In accordance with expectations from equation (1) the planet, having started migrating inwards faster than the gas, slows down under the influence of the dynamical corotation torque when it is released, but is trapped migrating inwards faster than the gas. How closely to the ideal, inviscid prediction of equation (1) the planet trajectory falls depends on the numerical diffusion induced by the finite grid resolution of the simulation. In addition to the fiducial resolution of , simulations at one half and one quarter of this resolution are shown, and the effect of the numerical diffusion limiting the vortensity contrast in the corotation region, and hence the size of the corotation torque can be seen. As the resolution is lowered, the planet migration rate diverges from the gas inflow rate in accordance with the expected effects of numerical diffusion mixing the libration region vortensity with the background disc.
The second experiment demonstrates regime (ii), and the results are given in Figure 2. Here the disk has a flat profile and , and the inflow is set to , such that the radial gas flow is on the order of . The planet potential is turned on after orbits, and the migration from is forced until orbits, which is . The planet is released at of the inwards radial offset that a gas parcel also starting at would have. In accordance with expectations from equation (1) the planet, having started migrating inwards slower than the gas, experiences an increasing dynamical corotation torque after it is released, and this eventually turns the planet around and causes it to migrate outwards as a runaway. At late times the planet migrates outwards at a roughly constant rate, dependent on the surface density of the disc. Measured outward migration rates for the three surface density values are given in Table 1, taken as the average over the radial interval where the motion is roughly constant in all cases. Sustained fast outward migration rates are observed to scale with surface density, but the proportionality is sub-linear. This suggests complicated dependencies apply, and the physics which determine regime (ii) fast migration rates is discussed in detail in Section 6.
The third experiment demonstrates regime (iii) and the results are given in Figure 3. Regime (iii) is the outward mode of migration with the planet asymptotically locked to the disc flow, a counterpart to regime (i). This disc has with , but the radial flow is now outwards with . The planet potential is turned on after orbits, and the migration from is forced until orbits, which is . For regime (iii), the planet should be moving outwards more slowly than the gas, so it is initially forced inwards to of the outward radial offset that a gas parcel also starting at would have. When the planet is released, the dynamical corotation torque causes its migration to rapidly reverse, but in accordance with expectations from equation (1) it is limited to moving outwards somewhat more slowly then the disc gas.
The fourth experiment demonstrates regime (iv) and the results are given in Figure 4. The disc used has and . A slower radial outward flow is driven with . The planet potential is turned on after orbits, and the migration from is forced until orbits, which is . To produce regime (iv) behaviour the planet is released at times the outward radial offset that a gas parcel also starting at would have. It was found that for regime (iv) migration to continue once the planet was released, that the torque on the planet at the release time needed to be approximately consistent with outward migration at the forced outward migration velocity. Even then, the outward regime (iv) migration displays a slowly growing oscillation, possibly indicative of an overstability.
Having demonstrated that the four regimes predicted by equation (1) exist and can be achieved in simulations, we proceed to demonstrate interesting consequences in some more physically plausible situations.
4 Reversing Migration
One interesting scenario suggested by equation (1) is where a planet initially migrating inwards in an inviscid disc, that is also not initially subject to an accretion-driving magnetic torque, can experience a migration torque reversal if the protoplanetary disc’s magnetic field undergoes secular evolution, and establishes a Hall-enabled laminar magnetic field torque with an associated midplane gas inflow. In other words, the planet transitions from regime (i) to regime (ii) because a magnetic torque is switched on after being absent during earlier evolution. To model this scenario, we let a planet migrate inwards in an untorqued, inviscid disc, and then ramp up the body force representing a magnetic field. Here it is ramped up over orbits, driving gas inflow. During the initial migration phase the corotation region acquires a vortensity deficit. Here, the planet migrates inwards in the disc from a lower to a higher vortensity region, while keeping the vortensity of librating coorbital material constant. To make the transition and induce the planet to run away outwards in regime (ii) migration, the vortensity of the corotation region will need to reverse from having a deficit to having an excess. During this transition, the dynamical corotation torque will pass though zero as the vortensity deficit passes through zero. As discussed in Section 2, this point is a bifurcation in the behaviour of the system. At this time, the migration torque will be entirely specified by the Lindblad (wave) torque. If the disc surface density is low enough that the Lindblad torque does not drive the planet faster than the gas flow at this most vulnerable point in the transition, then the planet can enter regime (ii) and experience a runaway corotation torque, eventually resulting in sustained fast outward migration. If, however, the surface density is high enough such that the inward migration driven by the Lindblad torque alone is faster than the disc inflow, then the planet will be trapped in regime (i) and will be driven inwards, faster than the gas. Thus, there is a bifurcation at a critical surface density, below which regime (ii) migration should result, and above which regime (i) migration should occur at a speed close to the gas inflow speed. Here we demonstrate this transition via customised numerical experiments.
We tune the disc parameters with the use of the torque formula from Paardekooper et al. (2011) to predict the Lindblad torque driven migration rate for a given planet mass and disc surface density, and choose two values on either side of the critical surface density. The disk has surface density with a radial flow , and a planet with mass is used. These experiments are shown in the left panel of Figure 5. When , below the critical surface density, the planet initially in regime (i) inward migration transitions into regime (ii) outward migration, and eventually ends up migrating outwards at a sustained fast speed. However, with double the surface density, that is , the planet is unable to transition into regime (ii) migration and is driven inwards with the disc gas in regime (i) migration.
In the first case, when the magnetic inflow torque is turned on, the vortensity of the corotation region is slowly reversed before runaway occurs. The characteristic timescale for this to occur can be expressed by use of a simple model for the driving of the inverse vortensity of the corotation region as:
To lowest order in this gives the timescale to reverse a inverse vortensity contrast of
We further discuss the sustained fast migration which occurs after this runaway phase, and conditions that can terminate fast migration, in Section 6.
In the right panel of Figure 5, we also demonstrate that the introduction of an outflow torque with to the disc has completely different results from the use of an inflow torque. If an outflow torque is introduced, both of these surface density choices result in the planet smoothly adopting regime (iii) outward migration, as the vortensity contrast of the corotation region does not need to reverse for the planet to move from regime (i) to (iii). In this regime (iii) simulation, the planet in the higher surface density disc migrates closer to the disc outflow speed, as is expected from the effects of numerical diffusion in the computation.
5 Viscous versus Laminar Accretion Flows
Given the contrasting behaviours that can be induced in the migration of a low mass planet by the introduction of a laminar magnetic disc torque, it is worthwhile to contrast these to the effects of a viscous torque driving exactly the same accretion rate. In this section, we repeat the experiment in the previous section that resulted in the planet reversing its migration, but instead of adding a laminar magnetic field torque to the disc we ramp up the disc viscosity to drive an equivalent inward accretion flow. The difference in the resulting migration of the planet serves to illustrate how the migration effects due to a laminar magnetic torque are not equivalent to those generated by a viscous accretion stress.
The experiment of Section 4 is repeated here with a ramped up viscosity with . This produces a steady-state accretion flow that can be characterised as having an inflow at . The resulting planet migration trajectory is shown along with the equivalent result with a magnetic torque in Figure 6. In the viscous accretion disc case, the planet continues to migrate inwards after the viscous inflow is established. This is the opposite of what happens in the magnetically driven inflow case. Why this happens is illustrated in the vortensity maps shown in Figure 7. The vortensity in the region near the planet is shown at time orbits, where in both cases the planet is at the same radial position, but is moving inwards in the viscous case and outwards in the magnetically torqued disc case. In the viscous case, the corotation torque is fully unsaturated by the viscous diffusion of vortensity, and this is able to maintain the background vortensity profile across the librating streamlines in the corotation region (Paardekooper et al., 2011; Paardekooper, 2014). Hence, in this case the corotation torque can be thought of as a static corotation torque which does not depend on this history of the planet’s migration, or on the history of the torques that have been applied to the corotation region, or on the planet’s migration rate. It operates in addition to the Lindblad torque, and steady inwards migration results. In the inviscid case, with a magnetic torque acting on the disc, the librating streamlines have a strongly enhanced vortensity and -asymmetrical shape that leads to a strong dynamical corotation torque that does depend on the history of the disc-planet system, at least until the moment when fast migration at an asymptotically constant rate sets in at late times (Paper I). This illustrates why the nature of corotation torques in viscous and laminar torqued discs are fundamentally different.
6 Fast migration
In Sections 3 and 4 experiments were shown where a planet in regime (ii) reversed its inward migration, underwent a period of runaway outward migration, during which the migration speed accelerated, and then settled into a period of sustained migration at constant speed. The runaway phase is simply understood as outlined when this regime was predicted in Paper I: the magnetic torque causes to decrease with respect to , and as the relative disc-planet velocity increases the torque and the rate of change of increases. However, this does not predict the existence of a sustained and stable outward migration velocity.
Sustained fast migration is possible when the torque stops increasing as the planet moves further and faster. From the phenomenology of slow migration, it is reasonable to expect that a process limiting the vortensity contrast between the libration region and the disc is important in setting this rate, and as such numerical diffusion has an important role in setting the fast migration rates observed in Sections 3 and 4. Thus, we instead explore this flow regime with simulations of static planets and an imposed body force driving the disk past the planet, which greatly reduces the numerical diffusion across the corotation region since this no longer migrates across the computational grid.
For predicting the torque in the fast migration regime, we find that for parameters where the vortensity gradient in the disc is small, an extension of the analytical model used for slow migration in Paper I performs well at early times. This establishes that the key parameter in setting the torque is still the vortensity contrast between the libration region and surrounding disc. However, for the late-time torque, which gives steady fast outward migration, this simple model breaks down in several ways, and it does not itself predict a limit on the (inverse) vortensity contrast in an inviscid disc. High resolution static planet simulations do, however, reveal a saturation mechanism, as described below. Having reached an understanding of why runaway migration leads eventually to steady fast migration, we decided that the most productive way to treat the resulting, somewhat complicated situation, is to seek a fit for the steady torque in the fast regime from a suite of numerical simulations. We have produced such a fit, and this is also presented and discussed below.
Before presenting the models that indicate what the saturated torque and migration rates should be in the steady fast regime, we first present representative simulations to elucidate the features of the flow in this case.
6.1 Streamlines, angular momentum exchange and ram-pressure stripping
For these simulation we use a static planet, as in Paper I, but also employ a radial mesh with refinement focussed on the corotation region. The radial grid spacing is described in Appendix A. The use of a planet on a fixed circular orbit with a fast radial disc flow is justified by the fact that the important parameter in the problem is the relative migration rate between the planet and disc, not the intrinsic migration rate of the planet. Furthermore, the phenomenon that leads to the transition from runaway migration to steady fast migration is observed in simulations with a static planet and rapid gas flow.
A first example, with , has a very shallow vortensity gradient, set by using . We remind the reader that corresponds to an inward radial gas flow velocity for which the time to cross the horseshoe libration region is just 30% of the horseshoe libration period, and hence corresponds to fast migration because . As we are interested in understanding what sets the corotation torque in the fast migration regime, it is worth recalling how this is determined in the slow migration limit under quasi-steady conditions before discussing how the pictures changes under fast migration.
Angular momentum is exchanged between the planet and fluid elements that undergo U-turns. In a simple symmetric switch model, applied to a disc that is flowing inwards past the planet, a fluid element that orbits outside of the planet’s orbit, at a distance (defined by equation 5 in the slow migration limit), is forced by the planet’s gravity to undergo a U-turn in front of the planet, and it jumps symmetrically from an exterior orbit to one interior to the planet. Clearly the maximum specific angular momentum exchange associated with these switches is defined by the value of . These fluid elements are of two types. Those orbiting closer to the planet’s orbit are on librating horseshoe streamlines, and their vortensity is defined by the vortensity of the corotating material, which is modified continuously by the magnetic torque that drives the radial accretion flow. Those fluid elements that orbit slightly further from the planet pass directly from the outer disc to the inner disc, and never librate. Their vortensity is defined by that in the background disc near the planet. The positive torque due to this material can be defined by summing over the inverse vortensity of all streamlines that undergo U-turns in front of the planet, as shown by equation (18) below. There are also fluid elements that undergo U-turns behind the planet as gas is propelled from interior to exterior orbits, and all of these elements are trapped on librating horseshoe streamlines, such that their vortensity is defined by that in the corotation region. The negative torque due to this material can also be determined by summing the inverse vortensity of all streamlines that undergo U-turns behind the planet, and the radial width of the region that is included in this sum corresponds to , whose value is defined as described above in the slow migration regime.
The situation is more complicated in the fast migration regime. First, the fluid elements that undergo U-turns in front of the planet all flow through directly from the outer to the inner disc. Second, the librating region and associated streamlines become confined to a “bubble" or island that sits behind the planet, as shown in Figure 8, which displays results from the simulation described above. The fluid elements that undergo U-turns behind the planet are all contained in this bubble. Third, one needs to be careful when defining the value of associated with the U-turns because the fast relative flow between planet and gas causes a relatively wide region of the disc exterior to the planet to flow from the outer to the inner disc, past the orbital radius of the planet, in one synodic period, even if the planet’s gravity is switched off.
When determining the value of for the streamlines that undergo U-turns in front of the planet and exchange angular momentum with it, one needs to identify the streamline in the flow when the planet’s gravity is switched off that corresponds to the streamline that passes closest to the planet and in front of it when its gravity is switched on. The value of is then the radial distance between this unperturbed streamline and the planet at the planet’s azimuthal location, as illustrated by Figure 8. Similarly, the value of for the streamlines that undergo U-turns behind the planet should be determined by examining the streamline in the flow with the planet’s gravity switched off that just passes through the tail of the librating bubble when the planet’s gravity is switched on (because the tail of the librating bubble defines the outermost librating streamline). The value of corresponds to the radial distance from the planet of this unperturbed streamline when it crosses the planet’s azimuthal location, as shown in Figure 8. Analysis of the streamlines plotted in Figure 8 gives an estimate of for the fluid elements that U-turn in front of the planet and for those that U-turn behind it, which should be compared with the value obtained in the slow migration regime. Hence, we see that for a disc with a small vortensity contrast between the librating island and background disc, even in the fast migration regime the value of closely matches that in the slow migration regime.
However, when the vortensity contrast between the libration island and the disc becomes large, the streamline geometry is modified compared to the simple picture described above. Figure 9, left panel, shows the vortensity map and critical streamline analysis for a run with and a large background vortensity gradient set by after 800 orbits. At this time the torque has saturated to a steady value. First, in contrast to Figure 8, it is clear that the libration island has narrowed significantly, so that the libration island U-turn behind the planet cannot be approximated by a symmetric switch with the same radial width as the flow through U-turn in front of the planet. This can be understood in terms of the local modification in the shear due to the change in vortensity of the libration region. As the vortensity increases, the shear in the libration region increases. The material in the libration region is thus making a horseshoe turn in a disk with a steeper rotation law that is hence narrower. Indeed, if the vortensity perturbation in the libration region is negative, the libration island becomes wider, due to the local shear in the libration region being smaller than the Keplerian value, as shown in the right panel of Figure 9 with a negative value of . This is in essence an extreme version of the alteration of the width of the libration region phenomenon described by Casoli & Masset (2009), due in their case to the gradient of vortensity across the libration region in a viscous barotropic disc.
The above discussion does not explain why runaway outward migration is observed to level out to a constant migration rate. To understand this we note that at late times in simulations like that shown in Figure 9, left panel, with a positive , we observe that the vortensity of the libration island has stopped growing, and a stream of high-vortensity material is continually stripped off from the front of the libration island and swept downstream into the disc. This interaction is shown in detail in Figure 10. Hence, it is the ram-pressure stripping of vortensity from the head of the librating island that explains why the runaway changes to a constant migration speed. This ram-pressure stripping also appears to be the main source of the vortensity stripes seen downstream of the planet in Figure 9. That there exists a resolvable mechanism which limits the growth of the vortensity contrast of the libration island in inviscid simulations allows us to construct an approximation for the torque as a function of relative gas-planet radial motion, that in turn gives a prediction for the sustained steady-state fast outward migration speed. Before doing so, we will develop an analytical model applicable in the low vortensity contrast regime, which further elucidates the nature of the fast migration regime torque beyond the discussion given above.
6.2 Analytical model
where and are values of inverse vortensity on two radial cuts located just in front of and behind the two horseshoe turns, and is the radial distance between the planet and fluid streamline that makes a switch. Formally the derivation of this expression invokes conservation of vortensity. To use it here we only require conservation over timescales of the U-turn, which is shorter than the libration time . Thus, we apply the argument from the Appendix A of Paper I, but to faster migration, up to . On the libration region turn, we know from the numerical simulations that is constant with respect to due to phase-mixing. The flow-through turn in front of the planet processes unmodified disc material, so to first order in the inverse vortensity on the integral contour drawn in front of the leading horseshoe turn is
It is clear only the leading order term is needed as is -independent, so we can write for the fast migration torque:
which evaluates to
The use of equation (18) shows how the transition from the slow migration regime of Paper I occurs - here the two horseshoe turns become completely disconnected so no sections of the two turns cancel out. This occurs because the turns in front of the planet consist of fluid elements that pass directly from exterior to interior orbits when the relative radial motion between the planet and disc material is fast. The turns behind the planet consist only of librating fluid elements.
In this case we analytically predict fast outward migration due to the fast torque as
With a static planet, the time evolution of inverse vortensity is given by
This form agrees well with simulations when the flow stays close to the simple model for the libration region. For example, with a planet mass , density power law and , very good agreement is found with the analytical model, as shown in Figure 11. This model, and the excellent agreement with numerical experiment, make it clear that the dynamical corotation torque continues to apply even as the planet-disc relative radial velocity becomes fast (). As the vortensity contrast in the libration island grows, however, we find that the flow deviates from the analytical model so much as to render the model inapplicable, as shown in Figure 9. Hence, we resort to a numerical parameter study of migration torques in this area of parameter space that is presented in the next section.
6.3 Numerical parameter study and model
We find in our simulations that as the vortensity contrast between the libration island and the disc grows, the width of the the libration island narrows, going beyond the regime of small vortensity contrasts where the simple model used in the previous section applies. Moreover, at late times, the vortensity contrast saturates due to ram-pressure stripping of high vortensity material off the head of the libration island, as shown in Figure 10. Hence, so that progress might be made in understanding the consequences of regime (ii) migration, we present here a simple numerical analysis of the late-time torques in a inviscid disc for a range of disc and planet parameters.
To produce a model for the late-time dynamical corotation torque on a planet in regime (ii) migration, we have run a parameter scan of static planet simulations with varying planet mass, disc flow rate, and surface density gradient. Planet masses , , , disc flow rates , , and surface density radial power law indices , , . Expending on the order of a million CPU hours on evolving these simulations, we found late-time convergence to a constant torque, although in two cases (, , and ) we were not able to evolve the simulation for a long enough physical time. However, for the cases which did evolve well, we found no significant dependence of the dimensionless total torque, , on the planet mass ratio . Therefore, we averaged the values obtained for different , and to derive the final measured corotation torque we subtracted from the total torque measured in our simulations the Lindblad torque as predicted by the formulas in Paardekooper et al. (2010). We include a reasonable limit of zero corotation torque when (zero vortensity gradient) in the fit, although it is not part of the simulated parameter scan. The numerical values are reported in Figure 12, and in Appendix B we provide a prescription for the torque as function of the surface density power-law index, , and the disc accretion flow parameter, , based on a simple bi-linear interpolation of the values given in Figure 12. We also give a more extended discussion there about how one may practically model a magnetically torqued disc with an embedded planet in the context of N-body simulations of planet formation or population synthesis models.
6.4 Termination of fast migration
As the previous section has established the mechanisms which allow sustained fast regime (ii) outward migration to occur, we explore here one mechanism by which such an episode can end. Experiments like in Section 4 have explored the migration behaviour of a planet when the disc inflow torque is increased, so here we explore the effect of decreasing this torque with time on a planet in regime (ii) migration. This can cause the planet to terminate its fast migration, and return to regime (i).
Thus, our experiment is an extension in time of the run in Section 4 with , , and an inflow torque , but now is ramped to zero by decreasing in equation (14) over or orbits starting at orbits. The planet trajectories for this relevant part of the simulation are shown in Figure 13. The vortensity profile as shown in Figure 14 is smooth, due to numerical diffusion as the planet moves radially across the grid. However, it is still clearly much narrower than the canonical value of which would be expected in the left panel, as judged by the width of the inner white contour. As shown in Figure 13, although the ramping timescale varies by a factor of five, the timescale for the planet to transition to inward regime (i) migration does not vary by such a large factor. Details of the flow shown in Figure 14 offer an explanation for this. As the planet-gas relative velocity slows, an amount of low-vortensity-contrast material that was originally outside of the libration island begins making the horseshoe turn behind the planet along with the preexisting libration region material. Although the planet is constantly moving outwards in the disc, in the sequence of three panels the outer white contour of vortensity moves steadily rightwards to the outside of the planet’s position, indicating the flow of this low vortensity material around the horseshoe turn. Hence, the material making the rear horseshoe turn becomes more symmetrical with that making the front horseshoe turn, decreasing the corotation torque expected from equation (18). Given this phenomenology, one expects that the timescale for the drop from regime (ii) to regime (i) migration should occur on timescales related to the horseshoe turn. For the narrowed libration region observed at with half width of the associated is orbits, whereas for the full horseshoe turn width as applies with small vortensity contrasts of the associated timescale is orbits. The termination of fast regime (ii) migration evidentially involves motions of gas across this range of timescales.
When planets in regime (ii) reach a radial position such that the magnetic field driven disc inflow significantly decreases, the most simple scenario is that they will drop to regime (i) migration and asymptotically stall with respect to the gas. This forms a basic conjecture for the fate of planets in regime (ii) as they reach the outer edge of the dead zone, or enter a region of the disc where the magnetic torque acting on material near the midplane undergoes a significant decrease.
7 Viscous steady states
Our discussion of corotation torques in a magnetically torqued disc to this point has focussed on inviscid discs, where there is no mixing of vortensity between streamlines, only evolution of the vortensity on those streamlines due to forces exerted by ordered large scale magnetic fields. In principle, however, some level of turbulent diffusion may also be present in a protoplanetary disc, either because the magneto-rotational instability operates intermittently at a low level, as described by Simon et al. (2015b), or because a hydrodynamic instability, such as the vertical shear instability (e.g. Nelson et al., 2013), operates in the presence of the magnetic torques.
This issue is relevant because our torque model expressed by equation (1) does not predict a steady-state for a live planet beyond asymptotically locking the planet to the disc radial flow (no relative radial motion between the disc gas and planet). However, when a finite viscosity is introduced to mix the vortensity of the libration region with the surrounding disc flow, the vortensity contrast and corotation torque can settle into a steady-state even when there is relative motion between the planet and gas. Thus, we have extended the approach taken in Paardekooper (2014) for accounting for a finite turbulent viscosity in the disc to the situation at hand. This gives rise to a steady-state corotation torque.
Our analysis of this problem, which includes an analytic study and a suite of numerical simulations that are compared with the analysis, is presented in Appendix C. Generally speaking, decent agreement between the two is obtained in the parameter domain for which the analytical theory applies, hence providing a means of predicting the corotation torque in a magnetically torqued disc in which diffusion of vortensity also occurs.
Our analytical calculations and simulations show a diversity of behaviour for low mass planets in laminar magnetically torqued discs, including migration that is essentially locked to the disc flow, and runaway outward migration that saturates and transitions to steady fast migration. A significant limitation in our attempts to model dynamical corotation torques in the inviscid limit has been the effect of purely numerical diffusion in mixing the vortensity of the libration island and background disc flow, which inherently leads to resolution-dependent migration rates as demonstrated in Figure 1. In regimes (i) and (iii) we can understand from these simulations the underlying mechanisms, and extrapolate to a concept that the planet migration should be asymptotically locked to the radial motion of the disc. In the simulations of fast regime (ii) migration in Sections 3, 4 and 6.4, however, the analysis predicts continued runaway migration whereas the simulations demonstrate saturation of this runaway. Although our adoption of a refined grid and a non-migrating planet allowed us to explore fast migration under idealised conditions with reduced influence of numerical diffusion, and resulted in a natural saturation mechanism for runaway migration being identified, it is probable that the computed migration rates are slower than would be found in a fully inviscid disc. In these latter calculations the low numerical diffusion has come at the price of not having a fully dynamic planet, and hence the radial remapping technique of Benítez-Llambay et al. (2016) would be of particular utility in this respect. Furthermore, it is noteworthy that we have adopted a very simple equation of state in this study that removes baroclinic effects, and it is possible that under more general conditions instabilities might come into play as the vortensity contrast between the corotation region and surrounding disc grows. Such an instability was identified in a disc with a radial temperature gradient by Paardekooper (2014), and this effect might play a role in regulating the vortensity contrast and corotation torque under more realistic disc conditions.
Both inward and outward runaway migration can occur for higher mass (e.g. Saturn-mass) planets in viscous discs (Masset & Papaloizou, 2003), so some discussion of the differences between that phenomenon and the runaway migration discussed in this paper are warranted. The drift rate during regime (ii) outward migration is fast, in the sense defined by Masset & Papaloizou (2003), and is driven by corotation torques. Unlike in the classical analyses of so-called Type-III migration (Papaloizou et al., 2007; Pepliński et al., 2008a, b, c), however, vortensity-related dynamical corotation torques can have either a negative or positive feedback on planet migration, depending on the surface density power law (Paardekooper, 2014), whereas Type-III migration torques have only a positive feedback on planet migration. For inviscid discs with surface density power laws , runaway migration of low mass planets only exists in an outward-migrating form. Additionally, Type-III migration requires, for a given planet mass, a disc surface density above a critical value, whereas the torque reversal and runaway of our regime (ii) requires a disc surface density below a critical value (Section 4). Indeed, the mechanisms by which these two migration phenomena can be triggered and terminated also appear to be different. How fast regime (ii) phenomenology transitions to Type-III as the planet mass under consideration is increased remains as an interesting and as yet unexplored problem.
Our analysis in this paper applies to the classical Type-I migration regime in which the planet mass is small enough to not perturb the disc local surface density significantly. For an inviscid disc without an accretion flow near the midplane, however, nonlinear steeping of the spiral density waves launched at Lindblad resonances can lead to gap formation (Goodman & Rafikov, 2001) for relatively low mass planets, and to the phenomenon of the inertial limit discussed originally by Hourigan & Ward (1984), and considered more recently by Rafikov (2002), for which migration stalls completely when a planet forms a radially asymmetric gap (Li et al., 2009). The question of how things change in a laminar disc with a significant midplane accretion flow has not yet been addressed, but will almost certainly lead to a modification of previous estimates of the planet masses for which gap formation and migration stalling occur. Planets migrating in regimes (i) and (iii), for which the planet approximately locks to the radial drift of the disc, are of particular interest because, in a frame that migrates with the planet, there is essentially no relative radial motion between the disc and planet, and hence no flow to oppose the gap forming torques. A relative flow will develop, however, if the onset of gap formation leads to a slowing of planet migration, so understanding the interplay between migration, gap formation and laminar accretion flows for low mass planets is an interesting issue for future study, and one that can address important issues in planet formation such as the pebble isolation mass (Lambrechts et al., 2014).
Having adopted a very simple equation of state, we have not considered the role of entropy-related dynamical corotation torque effects as have been discussed for viscous discs by Pierens (2015) and Pierens & Raymond (2016). We expect that such entropy-related torques do not have the Galilean-invariance properties of the vortensity-related dynamical corotation torque because the rate that they build up is not proportional to the relative planet-disc gas radial velocity, but only the radial motion of the planet with relation to the disc’s entropy structure. However, the asymmetry of the libration island itself does have the same Galilean-invariant property. The basic expectation would be that an entropy contrast (or entropy deficit in the terminology of Pierens, 2015) retained by the trapped librating material will lead to a extra component of the dynamical corotation torque. For a non-viscous disc that is passively heated by the star, we would generally expect the entropy gradient to be positive, with entropy increasing as one moves away from the star. In the quasi-adiabatic limit, with inefficient heat exchange between the corotation region and the surrounding disc, a planet migrating inwards in regime (i) would migrate inwards at a faster rate due the influence of the entropy. A planet that is in regime (ii), but which starts to migrate inwards before turning around and migrating outwards would display more complicated behaviour. The entropy contrast would first slow the inwards migration, and would make the initial phase of the outward runway migration faster. Once the planet passes its point of origin, however, then the entropy contrast would change sign and would act as a drag on the outward migration. This scenario would be modified if instead the radial entropy gradient in the disc was negative, possibly due to magnetic dissipation accompanying the laminar magnetic stress or related to the dissipation in a disc wind. We note that, as opposed to being relaxed by viscous mixing and/or ram-pressure stripping, the entropy contrast of the libration island can be relaxed by radiative transfer effects. Unlike turbulent viscosity, this may not have a large variation between wind-driven and turbulent disc models. Clearly, the joint analysis of entropy-related and vortensity-related corotation torques in appropriate passively-heated wind-driven disc models is an additional area of interest for further work. The new rich behaviour introduced by the consideration of magnetically driven radial flows in laminar dead zones can be expected to produce a multitude of new scenarios for planet formation. Some possibilities which warrant exploration include:
A low mass planet at a small radius enters regime (ii) fast outward migration, exits that regime where the disc torques change significantly, and undergoes a regime (i) migration stalling at large radius where it is able to increase its mass by accreting locally.
Multiple planets in a disc evolve via regime (i) migration, such that their drift rates are controlled by the radial gas flow, altering the expectations for convergent Type-I migration and the production of resonant pairs and chains of planets.
Changing disc ionisation conditions, due to the evolution of the disc surface density and dust, lead in turn to changing non-ideal MHD conditions, magnetic field structures, and hence radial flows and coupling regimes, leaving a signature on planet formation.
Secular changes in the magnetic field configuration of the inner disc can be expected, so the motions of low mass planets and cores changes at different points in the disc’s lifetime.
Planets that form in discs where the vertical magnetic field is parallel to the disc rotation axis have their migration histories strongly influenced by Hall-effect-induced radial gas flows near the midplane. Planets that form in discs where the field and rotation vector are anti-parallel do not experience radial gas flows near the midplane, as Hall EMFs are ineffective, and hence can stall their migration completely due to dynamical corotation torques.
As full global non-ideal MHD simulations with sufficient dynamical range in space and time to both resolve the evolution of the disc, magnetic field, thermodynamics and chemistry, and the evolution of planets and their orbits and interactions are still computationally speaking, some ways off, we expect that the exploration of many of these scenarios will need to be conducted using the types of customised hydrodynamic simulations that we have presented in this paper, and using reduced N-body models, using parameterised versions of the migration torques that are presented in appendix B.
In Paper I we predicted that there are four regimes of migration for a low mass planet in a laminar magnetically torqued disc, and in this paper we have demonstrated that all four regimes can occur. We have demonstrated that inwardly migrating planets can undergo migration torque reversals and migrate outwards due to the introduction of laminar magnetic disc torques. When magnetic torques drive rapid gas inflow, we have shown that for a given planet mass a critical disc surface density exists, below which torque reversal and outward runaway migration is possible. We have demonstrated that this runaway migration eventually saturates, and transitions to steady fast migration, and we have identified the process that causes this transition. Our simulations show that if the magnetic torques acting on the disc decrease over time, the planet drops out of the fast migration regime and migrates inwards at approximately the flow rate of the gas. This is a general result for slower gas flows, as we have shown that in this case planet migration can approximately lock to the drift of the gas, such that the planet migration rate is determined by the radial velocity of the gas.
Our study has also examined why the corotation torque in a disc where viscosity drives radial gas flows differs from that in a laminar magnetically torqued disc, and we have demonstrated the critical role of viscous diffusion in smoothing vortensity contrasts such that the different behaviour is established. We have also examined the role of viscous diffusion in establishing steady state corotation torques in magnetically torqued discs that also support low level turbulent diffusion, and have developed expressions for these steady torques as a function of system parameters.
Finally, in Appendix B, we have provided a prescription for including the dynamical corotation torques discussed in this paper in N-body simulations of planet formation and population synthesis models.
Future work suggested by these results includes the application of novel techniques for reducing numerical diffusion, extension of the parameter regime to include less weakly coupled magnetic fields and transitions between disc parameter regimes, non-isothermal discs with stellar irradiation and internal radiative transfer, three dimensional models, and N-body models of planetary system formation applying the torque formulae given in this paper.
This research was supported by STFC Consolidated grants awarded to the QMUL Astronomy Unit 2015-2018 ST/M001202/1 and 2017-2020 ST/P000592/1, and in part by the National Science Foundation under Grant No. NSF PHY17-48958. This research utilised Queen Mary’s Apocrita HPC facility, supported by QMUL Research-IT (King et al., 2017); and the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/K00087X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. SJP is supported by a Royal Society University Research Fellowship.
- Bai (2013) Bai X.-N., 2013, ApJ, 772, 96
- Bai (2014a) Bai X.-N., 2014a, ApJ, 791, 73
- Bai (2014b) Bai X.-N., 2014b, ApJ, 791, 137
- Bai (2015) Bai X.-N., 2015, ApJ, 798, 84
- Bai (2016) Bai X.-N., 2016, ApJ, 821, 80
- Bai & Stone (2013) Bai X.-N., Stone J. M., 2013, ApJ, 769, 76
- Bai et al. (2016) Bai X.-N., Ye J., Goodman J., Yuan F., 2016, ApJ, 818, 152
- Baruteau & Masset (2008a) Baruteau C., Masset F., 2008a, ApJ, 672, 1054
- Baruteau & Masset (2008b) Baruteau C., Masset F., 2008b, ApJ, 678, 483
- Benítez Llambay & Masset (2015) Benítez Llambay P., Masset F., 2015, FARGO3D: Hydrodynamics/magnetohydrodynamics code, Astrophysics Source Code Library (ascl:1509.006)
- Benítez-Llambay & Masset (2016) Benítez-Llambay P., Masset F. S., 2016, ApJS, 223, 11
- Benítez-Llambay et al. (2016) Benítez-Llambay P., Ramos X. S., Beaugé C., Masset F. S., 2016, ApJ, 826, 13
- Béthune et al. (2017) Béthune W., Lesur G., Ferreira J., 2017, A&A, 600, A75
- Casoli & Masset (2009) Casoli J., Masset F. S., 2009, ApJ, 703, 845
- Goldreich & Tremaine (1979) Goldreich P., Tremaine S., 1979, ApJ, 233, 857
- Goodman & Rafikov (2001) Goodman J., Rafikov R. R., 2001, ApJ, 552, 793
- Gressel et al. (2015) Gressel O., Turner N. J., Nelson R. P., McNally C. P., 2015, ApJ, 801, 84
- Hourigan & Ward (1984) Hourigan K., Ward W. R., 1984, Icarus, 60, 29
- King et al. (2017) King T., Butcher S., Zalewski L., 2017, Technical report, Apocrita - High Performance Computing Cluster for Queen Mary University of London. Queen Mary University of London, doi:10.5281/zenodo.438045
- Kunz (2008) Kunz M. W., 2008, MNRAS, 385, 1494
- Kunz & Lesur (2013) Kunz M. W., Lesur G., 2013, MNRAS, 434, 2295
- Lambrechts et al. (2014) Lambrechts M., Johansen A., Morbidelli A., 2014, A&A, 572, A35
- Lesur et al. (2014) Lesur G., Kunz M. W., Fromang S., 2014, A&A, 566, A56
- Li et al. (2009) Li H., Lubow S. H., Li S., Lin D. N. C., 2009, ApJ, 690, L52
- Masset (2001) Masset F. S., 2001, ApJ, 558, 453
- Masset & Benítez-Llambay (2016) Masset F. S., Benítez-Llambay P., 2016, ApJ, 817, 19
- Masset & Casoli (2010) Masset F. S., Casoli J., 2010, ApJ, 723, 1393
- 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
- McNally et al. (2017) McNally C. P., Nelson R. P., Paardekooper S.-J., Gressel O., Lyra W., 2017, MNRAS, 472, 1565
- Nelson et al. (2013) Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610
- Ogihara et al. (2017) Ogihara M., Kokubo E., Suzuki T. K., Morbidelli A., Crida A., 2017, preprint, (arXiv:1710.01240)
- Paardekooper (2014) Paardekooper S.-J., 2014, MNRAS, 444, 2031
- Paardekooper & Mellema (2008) Paardekooper S.-J., Mellema G., 2008, A&A, 478, 245
- Paardekooper & Papaloizou (2008) Paardekooper S.-J., Papaloizou J. C. B., 2008, A&A, 485, 877
- 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
- Pandey & Wardle (2008) Pandey B. P., Wardle M., 2008, MNRAS, 385, 2269
- Papaloizou et al. (2007) Papaloizou J. C. B., Nelson R. P., Kley W., Masset F. S., Artymowicz P., 2007, Protostars and Planets V, pp 655–668
- Pepliński et al. (2008a) Pepliński A., Artymowicz P., Mellema G., 2008a, MNRAS, 386, 164
- Pepliński et al. (2008b) Pepliński A., Artymowicz P., Mellema G., 2008b, MNRAS, 386, 179
- Pepliński et al. (2008c) Pepliński A., Artymowicz P., Mellema G., 2008c, MNRAS, 387, 1063
- Pierens (2015) Pierens A., 2015, MNRAS, 454, 2003
- Pierens & Raymond (2016) Pierens A., Raymond S. N., 2016, MNRAS, 462, 4130
- Rafikov (2002) Rafikov R. R., 2002, ApJ, 572, 566
- Simon et al. (2015a) Simon J. B., Lesur G., Kunz M. W., Armitage P. J., 2015a, MNRAS, 454, 1117
- Simon et al. (2015b) Simon J. B., Lesur G., Kunz M. W., Armitage P. J., 2015b, MNRAS, 454, 1117
- Ward (1991) Ward W. R., 1991, in Lunar and Planetary Science Conference.
- Wardle & Ng (1999) Wardle M., Ng C., 1999, MNRAS, 303, 239
- de Val-Borro et al. (2006) de Val-Borro M., et al., 2006, MNRAS, 370, 529
Appendix A Radially refined grid
We choose a core region of with to each side of the planet, and grid it at a constant resolution , we then gradually increase the grid spacing outside that region, following a forth-order polynomial for the grid positions, thus matching the grid spacing at the join to the inner high resolution region and having zero first, second, and third derivatives of the grid spacing at the join.
Then, for a list of evenly spaced points in the interval the mapping which produces the grid is
which specifies a linearly spaced core grid of half-width transitioning into a variably spaced grid where the generating polynomial coefficients are
which is a polynomial with slope where it joins the core grid, and slope at the edge (). We use and . The final grid point positions are the mapping of this half-interval to half of the domain with
where and , so that the above specifies the grid points at . The outer half of the grid where is specified symmetrically about , as shown in Figure 15.
Appendix B Migration torque prescription for laminar magnetically torqued discs
Here we outline a simple approach to incorporating the corotation torque in a laminar, magnetically torqued disc in N-body simulations of planet formation or population synthesis calculations. We assume that the Lindblad torque will be included using a torque formula such as may be obtained from Paardekooper et al. (2010). The two most likely regimes for a planet to be in are the regimes (i) and (ii) discussed in Section 2, which both apply to a disc that sustains an inwards accretion flow onto the star, and here we focus on modelling the corotation torque in these two regimes.
We assume that the disc model has a power-law surface density profile of the form
and that the time dependent mass accretion rate is constant at each radius and is given by
where is the radius dependent radial velocity associated with the magnetically driven accretion flow. The accretion rate given in equation (33) is the component driven by the Hall-effect-induced horizontal magnetic fields, and does not include the contribution due to the launching of a magnetised wind near the disc surface. This surface component does not contribute to the evolution of the corotation torque in our model, but it would be relatively trivial to include its effects on the evolution of the disc’s mass budget if so desired. For such other unspecified accretion sources we include a second mass accretion term in the following.
If the disc has fixed values of the inner and outer radii, and , such that we ignore any evolution of these due to the torques that may be applied to them, then at each moment in time we can write
We assume that a model is computed by stepping forward in time using discrete time steps, , such that . In the absence of any matter sources, the disc accretes onto the star and the total disc mass evolves according to
In principle, one has a choice about how to treat the global, time dependent mass accretion rate through the disc. One choice could be to specify how and vary with time and to obtain by rearranging equation (33). Another choice could be to assume that is time independent, , and allow the mass accretion rate to be determined by equation (33). Here, for simplicity, we adopt the second of these choices such that decreases with time as and decrease, since we assume that no mass is added to the disc during its life time.
For a Keplerian disc we have the following relation between the radial velocity associated with the accretion flow and the azimuthal acceleration, , applied to the disc due to the Lorentz force
We note that the Lorentz acceleration acting on a two dimensional disc can be written
where is the current density, which may be expressed as , and is the magnetic field. The angled brackets indicate that the Lorentz force has been vertically averaged. When performing simple disc modelling, however, knowledge of the magnetic field strength and current density are not required. All that is needed to define how the disc and corotation torque evolves is knowledge of the torque that is applied to the disc, and not its mathematical form or physical origin333If one wants to include the dissipation of magnetic fields via Joule heating in the energy balance of the disc, then clearly a more sophisticated model that includes knowledge of the magnetic fields and currents would be required.. By defining an initial mass accretion rate and an initial surface density profile, the required value of can be determined from equation (33), and this then defines the value of that must operate in the disc. From equation (36) we have
This is used below to calculate the evolution of the corotation torque.
b.1 Corotation torque in regime (i)
We recall that regime (i) applies when both planet and disc gas are drifting inwards, and the planet is drifting faster than the gas. Quoting from equation (1) in Section 2, the corotation torque in regime (i) should be computed according to
where and are the inverse vortensities in the corotation region and in the background disc at the planet’s location, respectively, and subscript ‘p’ indicates that a quantity should be evaluated at the planet’s location. As demonstrated in Paper I (see equation (29) in Section 2.2.2 of that paper), the accretion-driving magnetic torque causes to evolve according to
For a model that is evolving via discrete time steps, the solution to equation (40) that can be used to evolve can be written
Integration of equation (41) can be initiated at time by setting , such that the inverse vortensity in the corotation region initially has the same value as the surrounding disc. The motion of the planet through the disc causes to evolve since for a Keplerian disc we have
In summary, when in regime (i) one can use equation (39) to calculate the corotation torque as a function of time, using equation (41) to update , and equation (42) to update as the planet migrates and the disc is torqued. At , the corotation torque because , and then is determined by the Lindblad torque. In regime (i), when , the migration will be slowed as the corotation torque evolves until .
b.2 Corotation torque in regime (ii)
We recall that regime (ii) corresponds to both planet and disc drifting inwards initially, but with the disc drifting in faster than the planet. Here, we expect the planet to slow down, stop, reverse its migration, undergo a period of accelerating outward migration, and to then settle into fast outward migration at a steady speed. As discussed in Section 6, we do not have an analytical prediction for what the torque should be during the fast steady phase of outward migration, since saturation of the runaway is caused by complex, nonlinear gas flows. Hence, to obtain an expression for this torque we adopt a simple formula that interpolates between the steady torque values obtained in the suite of runs described in Section 6.3 and presented in figure 12. The torque, as a function of the surface density power law in the disc and the flow-through parameter, , is given by
where, in accordance with conventional mathematical notation the indices of the matrix are in the order row, column (which differs from the array index notation of some programming languages). We recall that
We would expect that in a simulation of planet migration in a torqued inviscid disc that uses torque formulae such as those presented here to compute a planet’s orbital evolution, at the beginning of the calculation and the planet will migrate inwards at a rate determined by the Lindblad torque. If the planet is in regime (ii), then the corotation torque determined by equation (39) will increase, and will cause the planet to reverse its migration and enter a runaway phase. The runaway will saturate, however, when the corotation torques predicted by equations (39) and (B.2) are equal, and at this point one should calculate the steady torque given by equation (B.2).
b.3 Evolution with changing planet mass
So far we have only discussed how to calculate torques for planets of fixed mass, but to be useful in planet formation simulations we need to adjust the corotation torque when the planet gains mass because this widens the horseshoe region, which allows gas from the background disc enter the corotation region and start librating. The vortensity of this newly librating gas will clearly differ from that which was originally in the horseshoe region, and this will affect the corotation torque. To account for this, it is reasonable to assume that there will be mixing of material in the horseshoe region, and the value of should then be modified each time the planet mass is incremented such that it becomes an area weighted average of the original value of and .
Appendix C Model for Regimes (i) and (iii) with viscosity
The evolution of the corotation region in a laminar magnetically torqued disc leads to a sharp top-hat distribution of vortensity in the libration island, as shown in Paper I. Here we extend the approach taken in Paardekooper (2014) that accounts for a finite turbulent viscosity in the disc, and examine the steady state corotation torques that can that arise in the presence of vortensity diffusion. We first develop an analytical approximation, and then verify its application with a set of simulations. A similar final torque formula has been conjectured by Ogihara et al. (2017), although credited to Paardekooper (2014). In this section we make it clear how the model which generates this torque formula differs from the one given in Paardekooper (2014).
c.1 Analytical theory with viscosity
Before presenting the derivation, we note that this section develops a model for regime (i) and (iii) with a finite viscosity only, as the canonical model for the inviscid case in these regimes at late times is that the radial migration of the planet is simply the same as the radial motion of gas in the disc.
Our model for the case of finite viscosity follows from the ones of Paardekooper et al. (2011) and Paardekooper (2014). The evolution of the inverse vortensity, , in the corotation region is modelled with a one-dimensional viscous disc equation. Here, the radial scaling of the viscosity in assumed to be such that it does not drive accretion () but as the model only applies locally at the corotation region on the length scale (the half-width of the corotation region) the radial gradient does not significantly affect the results, so that it can be used for general viscosity power laws.
In terms of the nondimensional radial coordinate , the evolution of the inverse vortensity in the disc is
where is the rectangular function ( and otherwise), is the accretion torque exerted on the gas by Maxwell stresses, and is the kinematic viscosity of the disc at the planet radius. The body force, or Lorentz force of the laminar magnetic field, is included as the final term in equation (45). Equation (45) can thus also be read as an extension to the model for the inviscid corotation region vortensity given in equation (29) of Paper I. Note that this model neglects all radial advection of vortensity, and thus consistently neglects the magnetic torque body force outside of the corotation region. As in Paardekooper (2014) we can seek a steady-state solution for in the corotation region.
The final driving term of equation (45) can be rewritten in terms of the steady-state inflow velocity driven by the magnetic field
where again . Using this form in equation (45) we then change to the variables and , and seek steady-state solutions for , which yields:
where is . It is now apparent where Galilean invariance of the torque will arise from in the expression in square brackets. Now, expanding to leading (zeroth) order in , as the corotation region width is small compared to the semimajor axis of the planet’s orbit, gives:
We can then write the solution as where is the unperturbed value of at the planet’s position and expand in (as the perturbation in vortensity due to the planet is small) to leading order which gives
This is essentially the same equation as Paardekooper (2014) is solving to get his equation (27) (except that our boundary condition is ), so we only give here the new solution in the corotation region, which is
Taking leads to
Substituting this model of the inverse vortensity of the corotation region into equation (1) gives a steady-state horseshoe torque for regime (i) and regime (iii) migration as
This form again encodes the Galilean invariance property of the dynamic corotation torque, in that the final velocity term is the relative disc gas-planet radial velocity. We now proceed to test this model with numerical simulations.
c.2 Numerical tests of viscous regime (i) and (iii)
To demonstrate the level of agreement between the steady-state analytical model equation (52) and numerical simulations we present here a series of tests run in the same the fiducial domain as in Section 3 and resolution . Again as in Section 5 we use the surface density scaling and viscosity radial scaling . First we ran a series of tests with a magnetic inflow torque producing , surface density , and viscosities , , and . The resulting migration rates are shown in Figure 16 along with the steady-state migration rate predicted by the total toque modelled as:
where is the Lindblad torque given by Paardekooper et al. (2011) their equation (3), is the barotropic corotation torque given by Paardekooper et al. (2011) their equation (32), and is given by equation (52). At the largest viscosity, the dynamical corotation torque contribution is very small, and the corotation torque is largely unsaturated, so that the total torque is essentially the linearly predicted Lindblad and unsaturated corotation torques. As the viscosity is decreased the relative contribution of the dynamical corotation torque increases. At the lower viscosities one can also see the impact of libration oscillations on the planet motion; an initial transient which dies away as the libration island becomes well mixed.
A second scan of viscosities was run with a midplane outflow torque, producing with disc surface density . The resulting migration rates and predictions of equation (53) for viscosities , , and are shown in Figure 17. Again at the lowest viscosity the libration oscillations are very apparent, but equation (53) does a quite reasonable job at predicting the migration rate and the transition from an unsaturated linear corotation torque to a dynamical corotation torque.
The theory leading to equation (52) has formally limited applicability when the corotation region vortensity evolution is not dominantly a diffusion problem with time-constant source terms. To quantify this one can evaluate the ratio between the crossing time of the radial flow across the corotation region (in the planet’s frame) and the diffusion timescale across the corotation region. For convenience we define this ratio as
and tabulate the characteristic values for each set of parameters in Figures 16 and 17 in Table 2. For the lower values of viscosity, is less than unity, indicating the the torque model would not be expected to preform well. However, the torque model generally follows the correct trend, suggesting that it does a reasonable job of including the relevant effects. We also list in Table 2 the values of the parameter controlling viscous unsaturation of the corotation torque as defined in Paardekooper et al. (2011), where for the corotation torque is unsaturated by viscosity. When this is the case, the dynamical corotation torque effect is also negligible, as the vortensity in the corotation region simply acquires a gradient matching the background disk flow. In this case the migration rate from equation (53) is essentially as given by the Paardekooper et al. (2011) formulae, with increasingly little contribution from dynamical torques.