Trapping low-mass planets at the inner edge of the protostellar disc
The formation of multiple close-in low-mass exoplanets is still a mystery. The challenge is to build a system wherein the outermost planet is beyond 0.2 AU from the star. Here we investigate how the prescription for type I planet migration affects the ability to trap multiple planets in a resonant chain near the inner edge of the protostellar disc. A sharp edge modelled as a hyperbolic tangent function coupled with supersonic corrections to the classical type I migration torques results in the innermost planets being pushed inside the cavity through resonant interaction with farther planets because migration is starward at slightly supersonic eccentricities. Planets below a few Earth masses are generally trapped in a resonant chain with the outermost planet near the disc edge, but long-term stability is not guaranteed. For more massive planets the migration is so fast that the eccentricity of the innermost resonant pair is excited to highly supersonic levels due to decreased damping on the innermost planet as it is pushed inside the cavity; collisions frequently occur and the system consists one or two intermediate-mass planets residing closer to the star than the disc’s inner edge. We found a neat pileup of resonant planets outside the disc edge only if the corotation torque does not rapidly diminish at high eccentricity. We call for detailed studies on planet migration near the disc’s inner edge, which is still uncertain, and for an improved understanding of eccentricity damping and disc torques in the supersonic regime.
Preventing low-mass planets from migrating to their host star is a long-standing problem. Mergers with their host star can be prevented if the protostellar disc has a sharp inner edge at a few stellar radii (Masset et al., 2006; Ogihara et al., 2010). Terquem & Papaloizou (2007) showed that migrating protoplanets usually end up in resonances. Some of these resided inside the disc’s inner cavity. Ogihara & Ida (2009) also build a resonant chain of low-mass protoplanets that stalled near the disc’s inner edge when the migration was artificially slowed down and the reduction of the corotation torque was ignored. Recently Izidoro et al. (2017) trapped a high number of low-mass planets outside the disc edge in a resonant chain, that subsequently needed to break to account for the currently-observed exoplanet period distribution. On the other hand, Matsumura et al. (2017) had trouble trapping multiple planets near the disc’s inner edge, even though their migration prescription was very similar to that of Izidoro et al. (2017): both include supersonic corrections to the migration and eccentricity damping timescales. Terquem & Papaloizou (2007) also included such corrections, but they followed the prescription of Papaloizou & Larwood (2000) while Izidoro et al. (2017) and Matsumura et al. (2017) followed Coleman & Nelson (2014). The simulations of Matsumura et al. (2017) usually resulted in one or two hot Neptune planets rather than a multiplet of smaller planets. The disparity between all of these results warrants further study.
2 Disc model and planet migration
2.1 Disc parameters
We assume steady-state accretion onto the Sun. The gas accretion rate is
where is the gas surface density, is the disc scale height and is the orbital frequency. The -viscosity is assumed to be constant (Shakura & Sunyaev, 1973). The disc scale height is , where with , is the Boltzmann constant, the proton mass and is the mean atomic mass of the gas. Now evolves as (Hartmann et al., 1998)
The extra 0.1 Myr avoids the logarithmic singularity (Bitsch et al., 2015b).
2.2 Disc inner edge implementation
Near the star the surface density of the disc is assumed to smoothly decrease to zero. Cossou et al. (2014) suggested that
where is the planet trap location from the star; here the surface density is maximal. The trap is at AU for solar-type stars; it is unclear how reliable the employed disc model is closer to the star where MHD effects become important. Trapping planets requires a sharp edge (Masset et al., 2006; Ogihara et al., 2010). Therefore we set the inner edge of the disc at , which is approximately 2 scale heights inside of . The surface density slope is computed as
where . Now as and all
planets in the type I regime cease migrating at or near .
To avoid the divergence of at , which could possibly cause numerical artefacts, we also tested a linear decrease in , for which . The resulting discontinuity of at was made smooth via a linear connection over length .
2.3 Planet migration
The gas disc exerts torques and tidal forces on the embedded planets which result in a combined effect of radial
migration and the damping of the eccentricity and inclination. For low-mass planets the migration is of type I
(Tanaka et al., 2002) while massive planets that are able to clear the gas in their vicinity experience type II migration
(Lin & Papaloizou, 1986). Here we are only interested in the former.
We follow Coleman & Nelson (2014) for computing the torque and the direction of migration. Their formulae are based on Paardekooper et al. (2011) for the torque and on Fendyke & Nelson (2014) and Cresswell & Nelson (2008) for the eccentricity damping, including corrections to the damping timescale and corotation torque in the supersonic regime when the eccentricity . We restrict ourselves to planar orbits. The normalised torque is (Paardekooper et al., 2011)
where and are the corotation and Lindblad torques respectively and is a normalisation constant. The corotation torque in a non-isothermal disc with thermal diffusion becomes
Here and depend on , and (Paardekooper et al., 2011). The functions , and determine the amount of torque saturation and are dependent on the planet mass and disc scale height (Paardekooper et al., 2011). In steady state , where . The remaining contributions are then (Paardekooper et al., 2011)
where is the negative entropy gradient. The
, and functions disallow writing the explicit dependence of on , but generally and .
The factors and are
where and . The eccentricity damping timescale is (Cresswell & Nelson, 2008)
where the wave timescale is (Tanaka & Ward, 2004)
The ’migration timescale’ (Cresswell & Nelson, 2008) is , and is
Despite its name, is not the actual migration timescale. By definition (Cresswell et al., 2007) and , so that when . Here is the timescale for semi-major axis evolution. In the supersonic regime hydrodynamical simulations show that the torque reverses direction when and the torque is maximal (and positive) when (Cresswell & Nelson, 2008). However, Cresswell & Nelson (2008) show that the planet always migrates inwards, despite the torque reversal at high eccentricity. The analytical approach of Muto et al. (2011) agrees with this result: inward migration persists at high eccentricity. Therefore, it is incorrect to use to compute the evolution of the semi-major axis; should be used instead (Matsumura et al., 2017). This is
Izidoro et al. (2017) adopted for the migration timescale, which led to the aforementioned difference. All timescales are positive for inward migration and negative for outward migration.
2.4 Planet migration near the disc edge
The reciprocal migration timescale, , and its two components, and
, are plotted in the top panels of Fig. 1 for the fiducial case (all
supersonic corrections enabled) as a function of planetary eccentricity. The middle panels show the individual
timescales and the bottom panels depict the normalised torques. The left column pertains to a 1
planet at 1 AU while the to the right the planet is at 0.099 AU. The outcome for the linear edge prescription is very
At 1 AU when we generally have and angular momentum transfer between the disc
and the planet usually leads to inward migration (); the example in Fig. 1, however,
has outward migration when because for this mass and the strength of the corotation
torque slightly exceeds that of the Lindblad torque (Paardekooper et al., 2011; Brasser et al., 2017); see bottom of Fig. 1. When we have and the inward migration is the result
of eccentricity damping at constant angular momentum. These two cases are common to all the combinations of the torque
formulae. However, the relations for , and are different when for the various
At 0.099 AU (near the disc’s inner edge), these two cases are also applicable, except that the difference between
and when is smaller than at 1 AU. When at 0.099 AU, however, the relations
for , and are different from that at 1 AU.
These fundamental differences are clearly seen between the two panels, assuming that equation (7)
and the supersonic corrections are applicable at the disc edge.
Suppose there is a single planet at . When a second planet approaches the first one, it is likely to become
trapped in resonance (Petrovich et al., 2013). The equilibrium eccentricity of both planets is the result of the balance of
migration and damping, and is (Goldreich & Schlichting, 2014). At this eccentricity , but when , and the planets migrate inwards: the innermost planet has its eccentricity
excited as the outer planet pushes the pair starward until the latter reaches ; the inner planet is
now parked inside the disc cavity. Should a third, fourth and additional planets approach the inner pair the mechanism
will likely repeat itself until the outermost planet is at the disc edge and all the other planets are inside of it, or
until the outermost planets can no longer push the inner chain deeper into the cavity due to the chain’s inertia.
The above arguments assume that the planet is surrounded by a gas disk on both sides. Near the edge, the torque is one-sided: at the inner disk edge there is only the outer Lindblad torque that pushes the planet inwards, but no inner Lindblad torque. Liu et al. (2017) implemented simplified one-sided corotation and Lindblad torques and showed that small planets could be trapped near the edge, but their implementation only works for an infinitely sharp edge.
3 Numerical methods
To study the behaviour of planets near the disc edge we performed a set of numerical N-body simulations consisting of a
solar-mass star and four equal-mass planets. These integrations used the symplectic -body code SyMBA (Duncan et al., 1998),
which was heavily modified to include the effects of eccentricity and inclination damping as well as planet migration
by the gas disc according to the formulation above (Matsumura et al., 2017).
We initially place the planets beyond their 2:1 mean-motion resonances. The planet’s masses are all either 0.1, 0.5,
2, 3, 5 or 8 Earth masses (). We compute the torques at each time step for each body; we apply eccentricity
damping only when . There is a surface density maximum at AU from the star. Closer to the
star than the disc edge at there is no migration nor damping.
Simulations are run for 2 Myr with a time step of 0.146 d. Bodies are removed when they are closer than 0.02 AU or
farther than 100 AU from the star, or when they collide. We assume perfect accretion during collisions. Initially
and . For simplicity we keep fixed despite its potential to change close
to the star where MRI effects (Bai & Stone, 2013) and disc ionisation are important (Gammie, 1996). The value of
affects the torque in a complicated manner as described in equation (8) (Paardekooper et al., 2011; Brasser et al., 2017). As we show
below, the most important factor in altering the torque is the supersonic corrections, in particular the corotation
reduction, , which is independent of and .
We test the dependence of the migration on eccentricity by selectively enabling or disabling the supersonic corrections to the corotation torque, the Lindblad torque and the eccentricity damping timescale. Apart from Masset et al. (2006) we are not aware of any systemic hydrodynamical studies of embedded planets near the disc edge. Motivated by this, we enable or disable (’on’ or ’off’) the corotation reduction term ( in equation (10), ’CR’) as well as the supersonic corrections to the eccentricity damping timescale (the term inside the parentheses in equation (11), ’ED’) and to the Lindblad torque ( in equation (10), ’EL’). This simple approach should suffice in the absence of more detailed torque prescriptions near the edge.
We show an example of the evolution of four planets of near the disc edge in Fig. 2. The
evolution is very different depending on the migration prescription. In the broadest sense the planets become trapped
in resonances outside the disc edge when the corotation reduction is not applied, and preferably when the
supersonic correction to the Lindblad torque, , is applied. If either or both of these are applied, multiple
planets will end up inside the disc cavity.
This outcome is in contrast with Izidoro et al. (2017) because they use to compute the semi-major
axis evolution of the planets instead of . Since at the disc edge for any value of , their
innermost planet can often stall any additional incoming planets even if it or the others are supersonic, though
the exact evolution is mass dependent. For example, in their Figure 4, the innermost planet is 10 and all
planets are trapped beyond 0.1 AU, while in their Figure 5 the innermost planet (initially) has a few and a
chain of planets is pushed inward as more massive planets migrate in. In resonance the approximation breaks down because the planets are supersonic; equation (14) should be used instead.
Our results are also in disagreement with Ogihara et al. (2010), who found that a sharp edge was able to prevent the planets from
falling into the cavity. Ogihara et al. (2010) concluded that the imbalance of increased drag on the planet at aphelion versus
little to no drag at perihelion caused the planet to be stationary at the edge with a non-zero eccentricity. However,
they did not include any supersonic corrections to their migration formulae.
In our approach the edge is fairly sharp, but the the tanh function quickly flattens beyond . Thus if the planet is near the drag at aphelion and perihelion is within 25% if
the eccentricity is . Only when the planet is roughly halfway between and is the
drag at aphelion (where the tanh function is 1) much stronger than at perihelion (where it is 0) and do we
possibly recover the situation from Ogihara et al. (2010), but only if the supersonic reduction of the corotation torque is
4.1 Different torque prescriptions yield different outcomes
What aspects of the torque prescription are responsible for the different behaviour at the edge? At low eccentricity
outward migration is caused by the corotation torque while inward migration is caused by the Lindblad torque; the
latter only depends on the temperature and surface density gradients (Paardekooper et al., 2011). When supersonic corrections are
considered the total torque also depends on . In general with our prescription is not a
monotonically decreasing function of , but instead has a peak near and a trough at
(Cresswell & Nelson, 2008). This behaviour is caused by how varies with eccentricity, and also because the maxima of
and do not coincide (top-right panel of Fig. 1). The migration
rate peaks for all planetary masses near , reaches a minimum when and a further maximum when . This non-monotonic behaviour is inconsistent with the analytical results of (Muto et al., 2011).
However, the behaviour is different when the planet is inside of the trap. At very low eccentricity, and
for nearly all planet masses, is negative and low, implying slow outward migration. Migration is
inward when and outward again at higher eccentricity. The top row of Fig. 3 is a contour map
showing at 1 AU (left), near the disc’s inner edge (middle), together with near the
The region of inward migration near moderate prevents the trap from stalling migrating planets in resonances, and
any planet situated at is pushed deep inside the cavity, together with any planets interior to it.
Therefore, the outcome of numerical simulations and the ability to trap planets near the disc edge depends on the exact
migration prescription employed. When using the migration is always outwards.
Figure 2 suggests that eliminating the corotation reduction, , and weakening the Lindblad torque
by applying provides the best prescription to trap multiple planets in resonance outside of ,
assuming the current torque formulae hold near the edge (cf. Liu et al. (2017)). The bottom row of Fig. 3
shows similar contour maps but now i.e. there is no corotation reduction. The behaviour is qualitatively
different everywhere: at the edge migration is always outward, but the strength is a complicated function of both the
planetary mass and the eccentricity.
Hydrodynamical simulations show that the corotation torque weakens as the eccentricity increases and mostly disappears
at (Cresswell et al., 2007). It thus appears to be unphysical to remove the corotation reduction far from the disc
edge, but it is unclear if this removal is applicable near the edge. The exponential reducion of Fendyke & Nelson (2014) does not
appear to hold for low values of ; the corotation reduction also depends on how the torque is measured.
Their Figs. 4 and 9 clearly show torque maxima near so that the exponential reduction may not be
universally applicable. Clearly more work is needed, both on the reduction itself but also how it behaves near the disc
In Fig. 2 the bottom-middle panel with all the supersonic effects enabled was able to temporarily trap the planets in a resonance even though the innermost planets were pushed into the cavity. Increasing the corotation torque kept the planets outside of the cavity (CR: OFF). This structure does not hold for higher-mass planets because they migrate faster and therefore excite themselves to higher eccentricities once the innermost planet is in the cavity and eccentricity damping is weak or non-existent. An example is shown in Fig. 4, which is the same as Fig. 2 but now the planets are 5.
The ability to trap multiple low-mass planets in a resonant chain outside the inner edge of the protostellar disc has
been investigated. These low-mass planets execute type I migration which pulls them invariably towards the star. In
the absence of a barrier these would all collide with the star. The disc’s inner edge could provide a trapping
mechanism (Masset et al., 2006). We have tested two types of sharp inner edges of the disc: a hyperbolic tangent and a linear
function, along with different migration prescriptions.
We find that a neat pileup of resonant planets outside the disc edge is established if the corotation torque does not rapidly diminish at high eccentricity. The expectation is that if the resonant chain of the planets remains outside the inner disc edge they eventually start orbit crossing and instigate a phase of giant impacts. This may account for formation of similar-sized, regularly spaced, non-resonant low-mass planets that are found to be common in relatively close-in regions by Kepler observations. However, the eccentricity damping and disc torques in the supersonic regime remain uncertain near the discâs inner edge. Due to resonant interactions, eccentricity is generally excited to values for which the migration is generally inward. Therefore we call for detailed studies on eccentricity damping and disc torques in the supersonic regime and near the disc edge. Such a study will play an important role in understanding the common architecture of compact systems.
- Bai & Stone (2013) Bai, X.-N., & Stone, J. M. 2013, ApJ, 767, 30
- Bitsch et al. (2013) Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & Dobbs-Dixon, I. 2013, A&A, 549, A124
- Bitsch et al. (2015a) Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28
- Bitsch et al. (2015b) Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112
- Brasser et al. (2017) Brasser, R., Bitsch, B., & Matsumura, S. 2017, AJ, 153, 222
- Coleman & Nelson (2014) Coleman, G. A. L., & Nelson, R. P. 2014, MNRAS, 445, 479
- Cossou et al. (2014) Cossou, C., Raymond, S. N., Hersant, F., & Pierens, A. 2014, A&A, 569, A56
- Cresswell et al. (2007) Cresswell, P., Dirksen, G., Kley, W., & Nelson, R. P. 2007, A&A, 473, 329
- Cresswell & Nelson (2008) Cresswell, P., & Nelson, R. P. 2008, A&A, 482, 677
- Duncan et al. (1998) Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067
- Fendyke & Nelson (2014) Fendyke, S. M., & Nelson, R. P. 2014, MNRAS, 437, 96
- Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355
- Garaud & Lin (2007) Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606
- Goldreich & Schlichting (2014) Goldreich, P., & Schlichting, H. E. 2014, AJ, 147, 32
- Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385
- Ida et al. (2016) Ida, S., Morbidelli, A. & Guillot, T. 2016. A& A 600, 154
- Izidoro et al. (2017) Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750
- Lin & Papaloizou (1986) Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846
- Liu et al. (2017) Liu, B., Ormel, C. W., & Lin, D. N. C. 2017, A&A, 601, A15
- Masset et al. (2006) Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478
- Matsumura et al. (2017) Matsumura, S., Brasser, R., & Ida, S. 2017, A&A, 607, A67
- Muto et al. (2011) Muto, T., Takeuchi, T., & Ida, S. 2011, ApJ, 737, 37
- Ogihara et al. (2010) Ogihara, M., Duncan, M. J., & Ida, S. 2010, ApJ, 721, 1184
- Ogihara & Ida (2009) Ogihara, M., & Ida, S. 2009, ApJ, 699, 824
- Ogihara et al. (2015) Ogihara, M., Morbidelli, A., & Guillot, T. 2015, A&A, 584, L1
- Oka et al. (2011) Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141
- Paardekooper et al. (2011) Paardekooper, S.-J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293
- Papaloizou & Larwood (2000) Papaloizou, J. C. B., & Larwood, J. D. 2000, MNRAS, 315, 823
- Petrovich et al. (2013) Petrovich, C., Malhotra, R., & Tremaine, S. 2013, ApJ, 770, 24
- Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
- Tanaka et al. (2002) Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257
- Tanaka & Ward (2004) Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388
- Terquem & Papaloizou (2007) Terquem, C., & Papaloizou, J. C. B. 2007, ApJ, 654, 1110