Saving Super-Earths: Interplay between Pebble Accretion and Type I Migration

The topic of planet formation has vexed modelers for many decades because it encompasses a range of physical processes that are tightly interlinked, yet occur at completely different length- and timescales. Planet formation takes place in protoplanetary disks that surround young newborn stars. Gas giants have to form within the lifetime of the protoplanetary disk, which is to say, within a few million years (Mamajek 2009), while terrestrial planets, especially in our own solar system, take an order of magnitude longer until they are fully assembled. For example, hafnium–tungsten dating indicates that the Earth's core formed in about 10 Myr (Yin et al. 2002), and the terrestrial planets themselves in less than 100 Myr, consistent with early dynamical models (Chambers 2001). On the other hand, the timescale and formation mechanisms for super-Earth planets, i.e., planets that are a few times more massive than the Earth but whose bulk composition is still predominantly composed of rock and possibly ice, are not well understood. Progress has been hampered by inadequate numerical approximations to planet migration, gas accretion, and tidal damping forces from the gas disk. The long timescale, short orbital period and high number of bodies make N-body simulations extremely challenging and time-consuming (Coleman & Nelson 2014; Cossou et al. 2014). Efforts to combine N-body simulations and chemical composition variation have been applied with moderate success to the terrestrial planets of the solar system (Bond et al. 2010; Matsumura et al. 2016), but applications to exoplanets have so far been a more formidable task (Fischer & Ciesla 2014; Moriarty et al. 2014).

In the core-accretion scenario of giant planets, a planetary core of about 10 Earth masses forms first and subsequently accretes an envelope of gas as soon as accretion of leftover solids is very low (Pollack et al. 1996). However, forming a planetary core of a few Earth masses through planetesimal accretion at large orbital distances in a disk with solar metallicity is not possible during the disk lifetime. The planetary cores tend to scatter planetesimals in their vicinity rather than accrete them (Tanaka & Ida 1999). In addition, planetesimal-driven migration will isolate some cores and prevent further accretion (Levison et al. 2010). However, recent models that consider the accretion of millimeter- to centimeter-sized objects—dubbed pebbles—have solved the timescale problem for the core formation of giant planets (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Levison et al. 2015a). These works demonstrated that efficient core growth at large orbital distances is possible well within the lifetime of the protoplanetary disks.

The reason for the increased accretion rate of pebble accretion versus planetesimal accretion is the effect of gas drag. As pebbles lose their angular momentum due to friction with the sub-Keplerian gas, they drift sunward. The rapid timescale for pebble migration and the strong headwinds they endure results in the pebble accretion cross-section of a planetary embryo being roughly equal to the radius of its Hill sphere rather than the gravitational cross-sectional radius. However, this rapid accretion only works when the time it takes for the pebble to drift past the planet's orbit is longer than the stopping time (Ormel & Klahr 2010). Otherwise, the pebbles just drift past the planet relatively untouched. Planetesimals, on the other hand, are generally scattered rather than accreted because the drag force that binds pebbles to the planet and results in accretion is virtually absent for typical-sized planetesimals.

The growth of a planetary core via pebbles crucially depends on the pebble size in a complicated manner (Ida et al. 2016). The typical pebble size is determined by an interplay between pebble growth and drift. Pebble accretion ends as soon as the planet starts to carve a partial annulus in the gas of the protoplanetary disk (Lambrechts et al. 2014), which changes the pressure gradient outside of the planetary orbit and thereby opens a corresponding annulus in the dust distribution of the disk (Paardekooper & Mellema 2006). This ends the supply of pebbles to the planet. At this point the planet has reached its pebble isolation mass, which strongly depends on the disk aspect ratio. When the planet completes accreting pebbles, it can contract a gaseous envelope, which allows the transition into runaway gas accretion as soon as the planetary envelope becomes more massive than the planetary core. The contraction timescale of the envelope is highly dependent on the mass of the planetary core, with a heavier core accelerating contraction (Ikoma et al. 2000).

During their growth, the planets also migrate through the protoplanetary disk. Small planets that barely influence the structure of the disk undergo type I migration, while gap-opening planets undergo type II migration (Goldreich & Tremaine 1980; Lin & Papaloizou 1986; Ward 1997). Type I migration is a combination of the Lindblad torque—caused by the spiral waves generated in the disk by the planet—and of the barotropic and entropy-related corotation torques. The Lindblad torque is generally negative and thus causes starward (inward) migration, while the corotiaton torque can be positive and therefore cause outward migration if it is strong enough to overcome the negative Lindblad torque.

The corotation torque arises from material orbiting on horseshoe orbits, which are corotating with the planet. Material circulating on orbits interior to the planet orbit the star slightly faster and thus catch up to the planet. Once they meet, the material recedes from the planet again because the planet deflects it to exterior orbits. If the disk interior to the planet is hotter than the exterior, the material cools and expands adiabatically when it is transported to outer orbits, generating a density difference across the planetary orbit. Similarly, material circulating on exterior orbits travels more slowly than the planet and is deflected inward, where the gas heats up and adiabatically compresses. The density gradient causes the corotation torque on the planet, which can drive outward migration.

However, this torque is prone to saturation, where the effects described above can saturate and outward migration stops. Two effects are important here: (i) heating and cooling, and (ii) viscosity. Without heating and cooling effects, the material orbiting on the horseshoe orbits would become fully mixed and the temperature gradient would disappear, resulting in the destruction of the over- and underdensity in front of and behind the planet. The corotation torque thus saturates, and outward migration is no longer possible. Heating and cooling is thus necessary to keep the temperature gradient in the corotation region, where a minimum temperature gradient is required to give this effect enough power to overcome the Lindblad torque and drive outward migration.

Viscosity, on the other hand, is needed to replenish the material inside the corotation region because this is where the planet extracts angular momentum from this region and the angular momentum of the corotation region is finite. An influx of new material inside the corotation region is thus needed to sustain outward migration. In order to supply the planet with more angular momentum at each U-turn, the viscous timescale needed to replenish the corotation region must be shorter than the horseshoe libration timescale—which is the time taken for material to flow through the entire corotation region—because viscous replenishment must be faster than the time needed for the material to complete a full cycle. Additionally, the viscous timescale must be longer than the U-turn timescale—which is the time the material needs to make a U-turn close to the planet—because a shorter viscous timescale would replace the material before it completes its U-turn and prevent it from transferring angular momentum to the planet.

Note here that the viscous timescale depends on the width of the horseshoe region⁴ squared, so it scales linearly with the planet mass, while the libration and U-turn timescales are proportional to the inverse of the planet mass (Baruteau & Masset 2008). This indicates that the torques on more massive planets saturate, unless the viscosity is high.

The strength of the torques also essentially depends on the radial gradients of surface density, temperature, and entropy (Paardekooper et al. 2011). Generally, a strong radial gradient in temperature (and thus entropy, because the two are related) is needed to sustain outward migration because it will cause a steeper density gradient across the planet orbit and thereby increase the power of the corotation torque. The delicate interplay between the torques can lead to regions in the protoplanetary disk where the total torque is zero and planets do not migrate. Such outward migration regions typically occur when the relative disk scale height h = H/r drops with radial distance (Bitsch et al. 2013). The planets could then, in principle, stay at these locations until the disk dissipates, if they do not grow to masses higher than the maximum mass at which outward migration can be sustained. Clearly, both planet growth and planet migration are related to the disk structure, so that the latter is of fundamental importance in understanding the mechanisms behind planet formation.

Bitsch et al. (2015b) combined planet growth via pebble accretion with recent models of planet migration to study the formation of planets in evolving protoplanetary disks with a prescribed evolution and viscosity. In that work the survival of super-Earth planets at orbital distances greater than 1 au was only possible when the planets remained in the region of outward migration, but at the same time did not outgrow it (Coleman et al. 2017). Here we expand on this work by studying the interplay between planetary growth via pebble accretion and type I migration for super-Earths (i.e., low-mass planets) for a variety of disk parameters. The aim of this study is to explore a finite set of disk parameters and determine which set allows a suite of low-mass planets formed through pebble accretion to survive the migration instead of colliding with the central star.

Our work is structured as follows. In Section 2 we review the foundations of type I migration and study the interplay between diffusion and viscosity on planet migration and growth. In Section 3 we study the influence of different temperature gradients and viscosity on the migration and survival of super-Earth planets. We then engage in a short discussion in Section 4 and present our conclusions in Section 5, followed by an appendix summarizing pebble accretion.

In this work we focus on a simple commonly used model that relies on the α viscosity to be constant (Shakura & Sunyaev 1973). More complex disk models involving MRI and/or disk winds (Bai 2016; Suzuki et al. 2016) would deviate from the single-α approach used here and are reserved for future studies.

We assume a steady accretion rate of the disk gas onto the Sun. The gas accretion rate is related to the gas surface density and scale height via (Pringle 1981)

$$\begin{eqnarray}&&{\dot{M}}_{* }=3\pi {\rm{\Sigma }}\nu =3\pi \alpha {\rm{\Sigma }}{H}^{2}{{\rm{\Omega }}}_{{\rm{K}}},\end{eqnarray} \tag{ 1 }$$

where Σ is the gas surface density, ν is the radially dependent viscosity, H is the disk scale height, and ${{\rm{\Omega }}}_{{\rm{K}}}$ is the Kepler frequency. The α-viscosity is assumed to be constant so that $\nu =\alpha {c}_{s}^{2}{{\rm{\Omega }}}_{{\rm{K}}}^{-1}$ (Shakura & Sunyaev 1973). The disk scale height is related to the temperature via H = c_s/Ω_K, where c_s = (k_BT/μm_p) is the isothermal sound speed, k_B is the Boltzmann constant, and μ = 2.3 is the mean atomic mass of the gas.

2.1. Basic Equations

For this study we investigate the migration of a low-mass planet according to the formulation of Paardekooper et al. (2011). The normalized torque on the planet is generally a function of four parameters, but in the derivation below, we can narrow this down to three. The normalized torque is given by

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Gamma }}}{{{\rm{\Gamma }}}_{0}}={{\rm{\Gamma }}}_{{\rm{c}}}+{{\rm{\Gamma }}}_{{\rm{L}}}={{\rm{\Gamma }}}_{{\rm{c}},\mathrm{baro}}+{{\rm{\Gamma }}}_{{\rm{c}},\mathrm{ent}}+{{\rm{\Gamma }}}_{{\rm{L}}}\end{eqnarray} \tag{ 2 }$$

where Γ_c and Γ_L are the corotation and Lindblad torques, respectively, and ${{\rm{\Gamma }}}_{0}={({m}_{p}/{m}_{\odot })}^{2}{(H/r)}^{-2}{\rm{\Sigma }}{{\rm{\Omega }}}_{{\rm{K}}}^{2}$ is a normalization constant. The corotation torque is further subdivided into contributions from the barotropic and entropy parts. The corotation torque in a non-isothermal disk with thermal diffusion becomes

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{c}},\mathrm{baro}}=F({p}_{\nu })G({p}_{\nu }){{\rm{\Gamma }}}_{\mathrm{hs},\mathrm{baro}}+[1-K({p}_{\nu })]{{\rm{\Gamma }}}_{{\rm{c}},\mathrm{lin},\mathrm{baro}},\end{eqnarray} \tag{ 3 }$$

for the barotropic part, while the entropy contribution is

$$\begin{eqnarray}{{\rm{\Gamma }}}_{{\rm{c}},\mathrm{ent}} & = & F({p}_{\nu })F({p}_{\chi }){[G({p}_{\nu })G({p}_{\chi })]}^{1/2}{{\rm{\Gamma }}}_{\mathrm{hs},\mathrm{ent}}\\ & & +\{[1-K({p}_{\nu })][1-K({p}_{\chi })]\}{}^{1/2}{{\rm{\Gamma }}}_{{\rm{c}},\mathrm{lin},\mathrm{ent}}.\end{eqnarray} \tag{ 4 }$$

The "lin" and "hs" subscripts indicate the Lindblad and horseshoe (corotation) parts of the torque, respectively (see below). The functions F(p), G(p), and K(p) are defined as (Paardekooper et al. 2011)

$$\begin{eqnarray}F(p) & = & \displaystyle \frac{8{I}_{4/3}(p)}{3{{pI}}_{1/3}(p)+\tfrac{9}{2}{p}^{2}{I}_{4/3}(p)}\approx \displaystyle \frac{1}{1+{(p/1.3)}^{2}}\\ G(p) & = & \displaystyle \frac{16}{25}{\left(\displaystyle \frac{45\pi }{8}\right)}^{3/4}\left[1-{\rm{\Theta }}\left(p-\sqrt{\displaystyle \frac{8}{45\pi }}\right)\right]{p}^{3/2}\\ & & +{\rm{\Theta }}\left(p-\sqrt{\displaystyle \frac{8}{45\pi }}\right)\left[1-\displaystyle \frac{9}{25}{\left(\displaystyle \frac{8}{45\pi }\right)}^{4/3}{p}^{-8/3}\right],\\ K(p) & = & \displaystyle \frac{16}{25}{\left(\displaystyle \frac{45\pi }{28}\right)}^{3/4}\left[1-{\rm{\Theta }}\left(p-\sqrt{\displaystyle \frac{28}{45\pi }}\right)\right]{p}^{3/2}\\ & & +{\rm{\Theta }}\left(p-\sqrt{\displaystyle \frac{28}{45\pi }}\right)\left[1-\displaystyle \frac{9}{25}{\left(\displaystyle \frac{28}{45\pi }\right)}^{4/3}{p}^{-8/3}\right].\end{eqnarray} \tag{ 5 }$$

In the above equations, I_ν(p) are modified Bessel functions of the first kind and Θ(x) is the Heaviside step function. The remaining torque terms are functions of the temperature and surface density gradients $q=-\tfrac{d\mathrm{ln}T}{d\mathrm{ln}r}$ and $s=-\tfrac{d\mathrm{ln}{\rm{\Sigma }}}{d\mathrm{ln}r}$ so that T(r) ∝ r^−q and Σ ∝ r^−s. In steady state, $q+s=\tfrac{3}{2}$ and only one is independent. It follows that $\tfrac{d\mathrm{ln}H}{d\mathrm{ln}r}=\tfrac{3}{2}-\tfrac{1}{2}q$. The remaining contributions to the torque are then given by (Paardekooper et al. 2011)

$$\begin{eqnarray}{{\rm{\Gamma }}}_{{\rm{L}}} & = & -\,2.5-1.7q+0.1s=-2.35-1.8q\\ {{\rm{\Gamma }}}_{\mathrm{hs},\mathrm{baro}} & = & 1.1(3/2-s)=1.1q\\ {{\rm{\Gamma }}}_{{\rm{c}},\mathrm{lin},\mathrm{baro}} & = & 0.7(3/2-s)=0.7q\\ {{\rm{\Gamma }}}_{\mathrm{hs},\mathrm{ent}} & = & \displaystyle \frac{7.9\xi }{\gamma }=5.6\left(\displaystyle \frac{7}{5}q-\displaystyle \frac{3}{5}\right)\\ {{\rm{\Gamma }}}_{{\rm{c}},\mathrm{lin},\mathrm{ent}} & = & \left(2.2-\displaystyle \frac{1.4}{\gamma }\right)\xi =0.8\left(\displaystyle \frac{7}{5}q-\displaystyle \frac{3}{5}\right),\end{eqnarray} \tag{ 6 }$$

where $\xi =-\tfrac{d\mathrm{ln}S}{d\mathrm{ln}r}=q-(\gamma -1)s=\tfrac{7}{5}q-\tfrac{3}{5}$ is the negative entropy gradient and γ is the ratio of specific heats, which is set to 7/5. When q = 3/7, which occurs in the outer parts of the disk where the disk temperature is dominated by stellar irradiation, the entropy gradient is zero, implying that outward migration is only possible in the inner parts of the disk that are dominated by viscous heating (Bitsch et al. 2013).

The parameter p_ν is defined as (Paardekooper et al. 2011)

$$\begin{eqnarray}&&{p}_{\nu }=\displaystyle \frac{1}{3}{{rx}}_{{\rm{s}}}\sqrt{\displaystyle \frac{2{{\rm{\Omega }}}_{K}{x}_{{\rm{s}}}}{\pi \nu }},\end{eqnarray} \tag{ 7 }$$

with ${x}_{{\rm{s}}}={\left(\tfrac{{m}_{p}}{{m}_{* }}\right)}^{1/2}{h}^{-1/2}$ being the horseshoe half-width. Ultimately,

$$\begin{eqnarray}{p}_{\nu } & = & \displaystyle \frac{2}{3}{(2\pi \alpha )}^{-1/2}{\left(\displaystyle \frac{{m}_{p}}{{m}_{* }}\right)}^{3/4}{h}^{-7/4}\\ & = & \displaystyle \frac{2}{3}{(2\pi \alpha )}^{-1/2}{\left(\displaystyle \frac{{m}_{p}}{{m}_{* }}\right)}^{3/4}{h}_{0}^{-7/4}{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{-7/8(1-q)},\end{eqnarray} \tag{ 8 }$$

where h = H/r is the reduced disk scale height, and we used the fact that in steady state $\tfrac{d\mathrm{ln}h}{d\mathrm{ln}r}=\tfrac{1}{2}(1-q).$ The value h₀ is the disk scale height at 1 au (or any other unit of distance). When q = 1, p_ν has no dependence on distance.

The parameter p_χ in Equation (5) is given by ${p}_{\chi }=\tfrac{3}{2}{\left(\tfrac{\nu }{\chi }\right)}^{1/2}{p}_{\nu }$ and χ is the coefficient of thermal diffusion. It is related to the temperature, surface density, and opacity of the disk via

$$\begin{eqnarray}&&\chi =\displaystyle \frac{16\gamma (\gamma -1){\sigma }_{\mathrm{SB}}{T}^{4}}{3\kappa {\rho }_{{\rm{g}}}^{2}{H}^{2}{{\rm{\Omega }}}_{{\rm{K}}}^{2}}=\displaystyle \frac{32\pi \gamma (\gamma -1){\sigma }_{\mathrm{SB}}{T}^{4}}{3\kappa {{\rm{\Sigma }}}^{2}{{\rm{\Omega }}}_{{\rm{K}}}^{2}}.\end{eqnarray} \tag{ 9 }$$

Here κ is the Rosseland mean opacity, ${\rho }_{{\rm{g}}}=\tfrac{{\rm{\Sigma }}}{\sqrt{2\pi }H}$ is the gas midplane density, and ${\sigma }_{\mathrm{SB}}$ is the Stefan-Boltzmann radiation constant, and we inserted a missing factor 4 from Paardekooper et al. (2011). Since χ and ν have the same dimensions, one may think of χ as a diffusive viscosity and then subsequently define a parameter $\chi ={\alpha }_{\chi }{H}^{2}{{\rm{\Omega }}}_{{\rm{K}}}$, analogous to the α parameter $\nu =\alpha {H}^{2}{{\rm{\Omega }}}_{{\rm{K}}}$. Then $\tfrac{\nu }{\chi }=\tfrac{\alpha }{{\alpha }_{\chi }}$, which is mostly a function of the opacity. For the sake of simplicity, in the rest of the paper we usually assume ν = χ, which is valid in optically thick regions, unless stated otherwise.

2.2. Conditions for Outward Migration

At low eccentricity, outward migration occurs when the positive corotation torque overcomes the negative Lindblad torque and thus depends on the values of q, p_ν, and p_χ. However, since p_χ is a (non-fixed) multiple of p_ν, we can analyze the above functions to find when the outward torque is a maximum as a function of p_ν and q only if we fix $\tfrac{\nu }{\chi }.$ In Figure 1 we show the various combinations of F, G, and K functions with ν = χ as they appear in the corotation torque. The $1-K(p)$ functions are monotonically decreasing, while the product F(p)G(p) can often be approximated as a log-normal function. The products F(p_ν)G(p_ν) and F(p_χ)G(p_χ) have the same underlying distribution and maximum value when p_ν = p_χ, but the mode is shifted toward lower values of p as $\tfrac{{p}_{\chi }}{{p}_{\nu }}$ increases. The torque is maximal when p_ν ∼ 0.25 because there the curves F(p_ν)F(p_χ)[G(p_ν)G(p_χ)]^1/2 and $\{[1-K({p}_{\nu })][1-K({p}_{\chi })]\}{}^{1/2}$ intersect, although the precise location of this intersection depends on the relative value of ν versus χ.

Zoom In Zoom Out Reset image size

The term F(p_ν)F(p_χ)[G(p_ν)G(p_χ)]^1/2 determines the saturation and the strength of the entropy part of the horseshoe torque (see Equation (4)) and is mostly responsible for the outward migration of the planet. The behavior of the product $F({p}_{\nu })F({p}_{\chi }){[G({p}_{\nu })G({p}_{\chi })]}^{1/2}$ is shown in Figure 2 as a function of p_ν and p_χ. The maximum occurs when $\chi =\tfrac{9}{4}\nu$. The product $\{[1-K({p}_{\nu })][1-K({p}_{\chi })]\}{}^{1/2}$ is monotonically decreasing in both p_ν and p_χ so that its global maximum is near the origin and thus is not displayed here. From Figure 2 one may infer that the outward torque is maximal when both p_ν and p_χ are in the range 0.3–0.7, setting limitations on the ratio $\tfrac{\nu }{\chi }$ and thus on the disk opacity for which outward migration is feasible.

Zoom In Zoom Out Reset image size

The condition for outward migration does not only depend on p_ν and p_χ, but also on the temperature gradient, q, through Equation (6). A contour plot of the influence of the temperature gradient on the total torque is shown in Figure 3 for χ = ν. Outward migration is only possible inside the contour of zero torque, and only occurs for a narrow range in p_ν and when q ≳ 0.87.

Zoom In Zoom Out Reset image size

In the inner disk, the temperature gradient is determined by viscous heating balanced with radiative cooling, and the cooling is inversely proportional to the opacity. The opacity gradient thus determines the cooling gradient and with it, the temperature structure of the disk. Outward migration is related not only to the temperature gradient (where a steeper slope favors outward migration), but also to the surface density gradient, where a shallow slope is preferred for outward migration. In a steady-state accretion disk, where the negative density and temperature gradients add up to 3/2, q ≳ 0.87 is needed to sustain outward migration, which is slightly lower than the typical temperature gradient of a viscously heated and irradiated disk (Garaud & Lin 2007; Chambers 2009; Bitsch et al. 2013). In disk models with heating and cooling, q can be as high as 1.2 in the viscous part through self-shadowing (Bitsch et al. 2015a), so we expect in the viscous region that 0.85 ≲ q ≲ 1.2. The outer part of the disk is dominated by stellar heating, so that as a consequence of the heating/cooling balance with stellar irradiation q = 3/7 always (Chiang & Goldreich 1997), although self-consistent models give a slightly steeper gradient of 4/7 in the outer parts of the shadowed region, that then become shallower and reach the 3/7 profile at larger orbital distances. (Bitsch et al. 2013). Therefore, outward migration is unlikely to occur in the outer region of the disk.

Increasing χ with respect to ν will shift the region of positive torque slightly toward lower p_ν and q, but once $\chi \gt \tfrac{9}{4}\nu$, the region of outward migration is confined to higher values of q because the maximum of $F({p}_{\nu })F({p}_{\chi }){[G({p}_{\nu })G({p}_{\chi })]}^{1/2}$ is lower and the entropy contribution is stronger than the barotropic contribution (see Equation (6)). Thus, the region(s) of outward migration is obtained from solving Γ = 0 for p_ν given q and p_χ, and then applying Equation (8).

So far, we have only parametrized the region of outward migration for a fixed ratio of $\tfrac{\nu }{\chi }$. However, the opacity of the disk is not a constant value as a function of the distance to the star, and thus we expect the above ratio to vary with distance and time. Since the maximum outward torque occurs when $\chi =\tfrac{9}{4}\nu$, it is obvious that variations in the opacity do not necessarily help in increasing the planetary mass and distance range for which outward migration is possible, although the same cannot be said for the opacity gradient. We return to the role of opacity later.

In the next subsection we focus on the effect of pebble accretion and how to effectively prevent planets from spiralling all the way to the star.

2.3. Outward Migration at the Pebble Accretion Isolation Mass

Except for inside the snow line, where the pebbles may lose their volatiles quickly and they become coupled to the gas, pebble accretion is capable of producing planets of a few Earth masses in less than a million years (Lambrechts et al. 2014; Ida et al. 2016); generally much shorter than the lifetime of the gas disk. It is thought that pebble accretion ceases to be effective when the planet has reached what is termed the pebble isolation mass (Lambrechts et al. 2014). Owing to angular momentum transfer between the planet and the disk, the material around the planetary orbit is pushed aside, which generates a pressure bump outside of the planetary orbit, trapping all pebbles. The pebbles can thus no longer reach the planet, and pebble accretion self terminates. The pebble isolation mass is approximately given by (Lambrechts et al. 2014)

$$\begin{eqnarray}&&{m}_{\mathrm{iso}}\sim 20{\left(\displaystyle \frac{h}{0.05}\right)}^{3}\,{M}_{\oplus }.\end{eqnarray} \tag{ 10 }$$

Lambrechts et al. (2014) briefly discuss whether the pebble isolation mass has an α dependence, and while they did not investigate this in great detail, they conclude that a weak dependence is possible. In what follows below, we assume no such dependence for the sake of simplicity, but we return to a possible weak dependence in Section 4.

The dependence of the pebble isolation mass on the (reduced) scale height as h³ is steeper than that of p_ν, which goes as h^7/3—see Equation (8). Since the total torque depends on p_ν (and q and p_χ), which itself is a function of both the planet mass and the disk scale height (and the viscosity), it is possible to determine where the contour of zero torque as a function of h and m_p intersects that of m_iso for a given value of α and q. Equating the isolation mass to the maximum mass for which the torque is zero and solved for h using Equation (8) results in

$$\begin{eqnarray}&&{h}_{\mathrm{crit}}=\displaystyle \frac{125\pi \sqrt{3}}{16}\alpha {p}_{\nu }^{2},\end{eqnarray} \tag{ 11 }$$

where the value of p_ν comes from Figure 3 and we used the approximation that m_⊕/m_⊙ = 3 × 10⁻⁶. This is a powerful diagnostic from which one can quickly narrow down which regions of the disk support outward migration and which disk conditions are necessary to trap planets, preventing them from reaching the central star. For example, the typical maximum value of p_ν for which outward migration occurs is p_ν ∼ 0.6, which results in a maximum scale height of ${h}_{\mathrm{crit},\max }\sim 0.015$ when α = 10⁻³. A corresponding minimum value occurs when p_ν ∼ 0.15 and ${h}_{\mathrm{crit},\min }\sim 0.001$ when α = 10⁻³. A potential weak dependence of m_iso on α of the form ${m}_{\mathrm{iso}}\propto {\left(\tfrac{\alpha }{{10}^{-3}}\right)}^{\beta }$ with β ≪ 1 would lower ${h}_{\mathrm{crit},\max }$ when α > 10⁻³ and increase it when α < 10⁻³.

We have plotted a number of cases in Figure 4 for various values of α and q. The left column has α = 10⁻³ and the right column has α = 5 × 10⁻³. The top row has q = 9/10 and q = 6/5 in the bottom row. We set χ = ν. The white curve depicts the pebble isolation mass.

Zoom In Zoom Out Reset image size

The growth and migration of a typical low-mass planet in these plots is similar to the following. With α fixed, planets initially grow at fixed h and at some point in time begin to migrate inward. When this happens, h decreases while m_p still increases. Eventually, the planet will reach the line where the torque is zero and becomes trapped at the outer edge of the positive torque region. This entrapment could occur before or after it has reached its isolation mass. As time evolves, the disk cools (see Equation (12) below), which implies that h decreases. Meanwhile, the planet mass remains fixed or still increases. This combination of disk evolution and possible continued planet growth implies that the planet eventually evolves out of the positive torque zone and migrates toward the central star (unless h or α somehow increases).

It is clear from the plots on the left side that stalling planets in the outward region of migration when α = 10⁻³ is only possible at very low isolation masses and scale heights (${h}_{\mathrm{crit},\max }\lesssim 0.015$). On the other hand, when α = 5 × 10⁻³, trapping planets in the outward migration region is possible even when q = 9/10, provided the disk is thin (h ≲ 0.03). From the plot it is evident that the planet isolation mass at which outward migration occurs increases for an increasing temperature gradient. This trend suggests there is a critical value of α where it is possible to prevent planets from reaching the central star for a given temperature gradient. We determine this threshold value of α and whether we can trap planets in the outward migration region in the next section.

In this section we present the results of numerical simulations of single-planet formation with pebble accretion as well as which disk conditions are needed to trap planets in the outward migration region.

3.1. Constant Opacity

From Figure 4 it is easy to discern which disk parameters will be able to prevent low-mass planets from migrating to the star. To test this hypothesis, we ran numerical simulations of single-planet formation with pebble accretion in both static and dynamic disks. We use the disk model from Ida et al. (2016), which is based on the works of Garaud & Lin (2007) and Oka et al. (2011), for simplicity, but the analysis below can also be applied to more complex disk models, e.g., Kretke & Lin (2010). The disk is assumed to be in a steady state and the temperature and surface density are power laws of the distance. The best fit for the temperature profile for a solar-type star is given by (Garaud & Lin 2007; Ida et al. 2016)

$$\begin{eqnarray}{T}_{\mathrm{vis}} & \simeq & 200{\alpha }_{3}^{-1/5}{\dot{M}}_{* 8}^{2/5}{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{-9/10}\ {\rm{K}},\\ {T}_{\mathrm{irr}} & \simeq & 150{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{-3/7}\ {\rm{K}},\end{eqnarray} \tag{ 12 }$$

where ${\dot{M}}_{* 8}={\dot{M}}_{* }/{10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and α₃ = α/10⁻³. The inner region of the disk is viscously heated, while the outer region is heated by the stellar flux. For Equation (12) we assumed that the disk is vertically optically thin, but radially optically thick. When the disk is so depleted that it is also optically thin in the radial direction, the temperature profile becomes $T\simeq 280{(r/1\mathrm{au})}^{-1/2}\,{\rm{K}}$ (Hayashi 1981). The radially optically thin condition occurs only when ${\dot{M}}_{* }\leqslant {10}^{-10}$ M_⊙ yr⁻¹. However, the condition for the disk to be radially optically thin only occurs at very low stellar accretion rate and accordingly a very low disk gas surface density. We do not consider this here because photoevaporation would destroy the disk before such low accretion rates can be attained (Alexander et al. 2014).

From the above temperature relation the reduced scale height becomes

$$\begin{eqnarray}{h}_{\mathrm{vis}} & \simeq & 0.027{\alpha }_{3}^{-1/10}{\dot{M}}_{* 8}^{1/5}{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{1/20},\\ {h}_{\mathrm{irr}} & \simeq & 0.024{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{2/7}.\end{eqnarray} \tag{ 13 }$$

The surface density then follows from the steady-state accretion. In the viscous region $d\mathrm{ln}{\rm{\Sigma }}/d\mathrm{ln}r=-3/5$ and in the irradiative region $d\mathrm{ln}{\rm{\Sigma }}/d\mathrm{ln}r=-15/14$. The boundary between the viscous and irradiation regimes given by T_vis = T_irr, which occurs at

$$\begin{eqnarray}&&{r}_{\mathrm{vi}}\simeq 1.8{\alpha }_{3}^{-14/33}{\dot{M}}_{* 8}^{28/33}\,\mathrm{au}.\end{eqnarray} \tag{ 14 }$$

The snow or ice line, which is the distance at which water ice condensates, is given by (Ida et al. 2016)

$$\begin{eqnarray}&&{r}_{\mathrm{snow}}\sim \max (1.2{\alpha }_{3}^{-2/9}{\dot{M}}_{* 8}^{4/9},0.75)\quad \mathrm{au}.\end{eqnarray} \tag{ 15 }$$

We employ this disk model and numerically simulate the evolution of a planet that forms through pebble accretion and migrates through the disk. We computed the evolutionary tracks according to Bitsch et al. (2015b), which is slightly different from the method of Ida et al. (2016). A summary of this method is given in the appendix. We have modified the code to employ the disk model discussed above; the pebble accretion model stays the same. We first assume a constant accretion rate of ${\dot{M}}_{* 8}=1,$ corresponding to a disk age of 1 Myr, and set the inner edge of the disk at 0.1 au because it is unclear how reliable the employed disk model is closer to the Sun where MRI effects play an increasingly important role. For simplicity we also set ν = χ and γ_eff = 7/5 (Paardekooper et al. 2011). Since we are interested in super-Earth formation and keeping the calculations feasible, we place planetary seeds from 0.1 to 5 au in succession for a range of initial values of α to determine the disk conditions that are required to prevent the planets from migrating to the star. Placing planetary seeds at greater orbital distances would allow these seeds to grow to larger planetary cores because of the increased scale height of the disk. These planetary cores would then efficiently accrete gaseous envelopes and form gas giants (Ikoma et al. 2000), hence our decision to place planetary seeds in the inner disk only.

In Figure 5 we plot the normalized torque as a function of distance to the Sun and planetary mass for our chosen static-disk model. The red region encompassed by the black line is where outward migration occurs, although the maximum mass for outward migation is always below 3 M_⊕. The region of outward migration is smaller than that found in Bitsch et al. (2015a) because we chose a different disk model. Bitsch et al. (2015a) found a smooth transition between the viscous and irradiative regimes of the disk, but it had a complicated structure, while in our simplified model the change is abrupt: the temperature and surface gradients suddenly change at r_vi, which is generally closer than 2 au to the star.

Zoom In Zoom Out Reset image size

When we perform single-planet pebble accretion simulations with migration for a range of initial distances to the Sun and viscosity, we end up with the result displayed in Figure 6. The black line indicates the critical value of α above which the planets end up on orbits farther than 0.1 au from the Sun; all other planets are lost to the star. It is clear that for a temperature slope of 9/10 a critical α ∼ 4 × 10⁻³ is necessary to prevent the planets from migrating to the star.

Zoom In Zoom Out Reset image size

Increasing the temperature gradient will decrease the critical value of α that is required to trap planets far from the star, while it also enhances the overall phase space where this trapping becomes possible. This is shown in Figure 7, which is identical to Figure 6, but now the temperature gradient is increased to q = 6/5 in the viscously heated part. As mentioned earlier, the typical range of the temperature gradient in the viscous region is 0.85 ≲ q ≲ 1.2, so Figures 6 and 7 represent approximate end member states. In Figure 7 we assumed that the whole viscous region of the disk has this steep slope, which is not observed in the simulations of Bitsch et al. (2015a). Thus the outcome should be treated as an extreme case. With the steeper temperature gradient the area of phase space in Figure 7, where planet trapping now becomes possible, has increased substantially. This enlarged region is caused by the elevated strength of the corotation torque with temperature gradient, as shown in Figure 3 and discussed in Section 1. In addition, a change in the temperature gradient also results in an opposite variation in the surface density gradient, which becomes shallower as the temperature gradient grows steeper. The shallower surface density gradient also supports outward migration.

Zoom In Zoom Out Reset image size

Thus, it seems that it is certainly possible to prevent low-mass planets from reaching the star under the right conditions, provided that a threshold value of α is surpassed; this critical value is a strong function of the temperature gradient, however.

3.2. Constant Opacity with Time Evolution

Circumstellar disks evolve with time, with the surface density and temperature readily decreasing with time, while α and q are often assumed to remain constant. The cases displayed above are for a disk with a fixed value of ${\dot{M}}_{* }$, while in reality, this quantity decreases with time according to (Hartmann et al. 1998)

$$\begin{eqnarray}&&\mathrm{log}{\dot{M}}_{* 8}=-1.4\mathrm{log}\left(\displaystyle \frac{t}{\mathrm{Myr}}+0.1\right).\end{eqnarray} \tag{ 16 }$$

The extra factor of 0.1 was added by Bitsch et al. (2015a) to avoid the singularity when t = 0. The relation in Equation (16) is based on observations and it is not very well constrained, but it is widely used by the community. From Equation (16) and a fixed value of α and a temperature relation, one can compute the surface density through Equation (1).

Very young disks have an accretion rate upward of 10⁻⁷ M_⊙ yr⁻¹, which results in a very hot inner disk with a high aspect ratio. Any planets that form early in this disk quickly reach a high isolation mass and are (eventually) trapped in the region of outward migration. However, as the disk evolves with time, the region of outward migration shrinks, the planets are released from it, and they fall into the central star (Bitsch et al. 2015b). In this study we only focus on cases where the disk age is 1 Myr or older; the maximum disk age is 3 Myr.

Figure 8 shows the region where planets may be trapped farther than 0.1 au as a function of α and initial formation location for an evolving disk. In the top panel the temperature gradient in the inner disk is 9/10, while it is 6/5 in the bottom panel. It is clear that the disk evolution increases the critical value of α and also shrinks the total area of the plot where the planets may be trapped.

Zoom In Zoom Out Reset image size

Increasing the initial disk lifetime to 2 Myr yields similar results, which are shown in Figure 9. The final masses of trapped planets are generally lower than for the younger disk because the aspect ratio is now diminished (the disk is cooler). However, the planets can be trapped more easily in the region of outward migration as the region itself barely evolves beyond this point as the disk cools. The region where planets survive in the disk is therefore larger than it is for the younger (and warmer) disk. Additionally, the region in r₀–α space where planets survive extends all the way to the outer disk. These planets have not had enough time to migrate all the way to the central star as the disk dissipates at 3 Myr. This can best be seen in the top panel of Figure 9, where the two regions of surviving planets appear separated, while in the bottom panel, where the temperature gradient is steeper, the variation in the critical value of α as a function of r₀ is less pronounced.

Zoom In Zoom Out Reset image size

3.3. Variable Opacity

Thus far, we have avoided discussing the influence of the disk opacity because it adds another dimension to the problem, increasing its complexity. Unfortunately, there are various popular opacity laws (Bell & Lin 1994; Bell et al. 1997; Semenov et al. 2003), so the model choice will affect the outcome. Here we use the opacity profile of Bell & Lin (1994) because this was used by Bitsch et al. (2015a) and it makes for easier comparison with the full-disk model discussed in the next subsection.

The temperature structure of the protoplanetary disk in the viscously dominated part is determined by the balance between heating and cooling, where the cooling function is inversely proportional to the opacity of the disk. A steep gradient in the opacity function will thus result in a steep gradient in the temperature profile. However, the cooling also depends on the density of the protoplanetary disk, where a higher density blocks outgoing radiation and decreases the cooling rate. The heating, on the other hand, also depends on the density and the viscosity. The only way to safely determine the temperature structure of protoplanetary disks are thus 2D simulations, including heating and cooling (Bitsch et al. 2015a). However, these simulations are beyond the scope of this work because we wish to disentangle the interplay between the different factors responsible for growth and migration of super-Earth planets. For this purpose, we rely on simple power-law disks, which gives us an easier understanding of the interactions of the parameters, even if the temperature and opacity gradients may not be fully consistent everywhere.

Put more simply, the role of the opacity is to change the value of χ and thereby the cooling rate of the disk, effective ratio of the specific heats of the disk, γ_eff, and the value of p_χ versus p_ν (Paardekooper et al. 2011). Here we wish to show one case of how the opacity likely affects the outcome, keeping in mind that there is a larger variety of outcomes depending on the opacity relation that is employed and assumptions about how the opacity scales with temperature and density (Garaud & Lin 2007).

In Figure 10 we employ a static disk with ${\dot{M}}_{* 8}=1$, corresponding to an age of 1 Myr, but the opacity is allowed to vary. Specifically, we use the opacity relations of Bell & Lin (1994). It is clear that allowing the opacity to vary is not of much help in the top panel, for which the temperature gradient in the inner disk is 9/10. The critical α is higher than in the static disk with constant opacity, and the region encompassed by the black line is smaller. However, when the temperature gradient is 6/5, as in the bottom panel, the critical α is reduced below 10⁻³ and the region encompassed by the black line is much greater than in the static disk of Figure 7.

Zoom In Zoom Out Reset image size

3.4. Variable Opacity with Time Evolution

The migration map that depicts the scenario when we add time evolution to the simplified disk with a variable opacity but a static temperature gradient is depicted in Figure 11. It is clear that this is substantially different than what was displayed in Figure 5. Indeed, when q = 9/10 (top panel), there is no region of outward migration, precisely because the opacity changes χ/ν, which potentially lowers the effect of the corotation torque. In the bottom panel, where q = 6/5, the region of outward migration is substantial, as expected. The final masses and area where planets may be trapped in the outward migration region is displayed in Figure 12. Note that in the top panel, almost no planets are saved from falling into the star, while in the bottom panel most planets are fine.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, the inclusion of a varying opacity only appears to be beneficial for steep temperature profiles, but is a hindrance for shallow temperature profiles because outward migration seems to be more sensitive to the value of χ/ν.

3.5. Full Model

Here we briefly touch upon using a full-disk model such as the model in Bitsch et al. (2014).

In Figure 13 the migration map for the full-disk model of Bitsch et al. (2015a) is displayed a the point when the disk age is 1 Myr. Note that there are two regions of outward migration: one close to the star, and one at moderate distance. These are caused by transitions in the opacity. The outer region was not present in the simplified model, which assumes that the disk is fully irradiative beyond r_vi, in contrast to the disk model of Bitsch et al. (2015a), where a shadowed region exists for r > r_ice allowing outward migration. The inner migration region is caused by the silicate line, and the outer one is cause by the water ice line.

Zoom In Zoom Out Reset image size

Unfortunately, we cannot display a figure similar to Figure 6 because the 2D disk model simulations performed so far used α = 0.054 to model the viscous heating. Computing the evolution of disks with varying α is computationally too expensive. Thus the migration map should be used as an indication of the outcome. Instead, in Figure 14 we present the final positions and masses of planets as a function of their initial birth distance for this single full-disk model and the simplified power-law model we used earlier for the same α and with the viscous temperature gradient set to 9/10.

Zoom In Zoom Out Reset image size

It is clear that planets in the power-law disk only survive in the inner part, where the temperature slope is steep enough to trigger outward migration. In between the two regions of outward migration in the disk model of Bitsch et al. (2015a), the planets grow too massive (about 3.5 M_⊕) to be contained in the inner region of outward migration and thus move all the way to the inner disk edge. Note that the final migration map at t = 3 Myr (not shown) allows only even smaller planets to stay in the region of outward migration compared to Figure 13.

Thus the full-disk model certainly helps in stalling planets farther from the star. The main reason is that there are two regions of outward migration (see Figure 13). The first is inside 0.5 au and is the result of viscous heating of the disk. The second region between 2 and 4 au is the result of shadowing of the disk due to a drop in the opacity profile. In the shadowed region h shrinks with distance and thus the temperature gradient is sufficiently steep for the corotation torque to dominate the migration.

Several topics require further discussion, but go beyond the scope of the paper and could be the subject of future studies.

In our models the temperature profile does not evolve in the outer parts, beyond r_vi, and thus the aspect ratio there is constant in time. Only the surface density decreases because ${\dot{M}}_{* 8}$ decreases. However, this is unlikely to affect the overall result because in the outer part of the disk type I migration is almost always inward (see Section 2.2).

The evolution of the inner part of the disk is more complicated than we assumed here. Even though the temperature and surface density gradients remain the same and the steep temperature gradient allows for outward migration, the temperature itself decreases, and thus the disk aspect ratio also shrinks. Since the temperature in the viscous region declines, stellar heating will eventually take over and the irradiative regime moves inwards. This implies that the zero-torque region must also shift sunward because the amount of viscous heating subsides. Thus all planets then follow the zero-torque line and only the most massive ones would pop out of the outward migration region and subsequently fall into the star.

Is there a way to prevent this outcome? One way could be to assume that the disk dissipates very rapidly once the planets have halted in the outward migration region. In order for this to work, the disk dissipation timescale would have to be shorter than the migration timescale. The migration time is normally a few hundred thousand years for an Earth-mass planet, which is much shorter than the dissipation time of the disk, especially when the planet is at a few au. Only when the stellar accretion rate drops to 10⁻¹⁰ M_⊙ yr⁻¹ is the clearing of the disk by photoevaporation much faster (Bitsch et al. 2015a) and comparable to the migration timescale. In summary, this scenario is probably implausible most of the time.

In this work, we focused on planets that grow to their pebble isolation masses. However, it is also feasible to imagine that the final masses of planets are m_p < m_iso. Such planets could form late. Alternatively, they could appear after the pebble flux dries up because planets farther out have reached their isolation mass before the inner ones have, and thereby prevent them from growing any further. Such planets could also be grown by planetesimal accretion, whose "planetesimal" isolation mass (Lissauer 1987) is generally lower than the pebble isolation mass. In this case, the same issues apply as for a planet that has reached the isolation mass. Whether such a planet remains trapped in the outward migration region depends on the precise evolution and properties of the disk. It is not possible to determine this a priori.

The exact formation pathway of super-Earths is still a mystery. Do they form in situ or via migration? In one of the earliest studies of its kind, Ogihara & Ida (2009) simulated rocky planet formation around an M-dwarf star starting with a disk of planetesimals and evolving it all the way past a possibly giant impact phase. They studied two end cases: one with the full strength of type I migration included, and one where the migration was artificially slowed down. They discovered that in the slow-migration case, almost all planets ended up in a multiple resonant chain. With full migration they only formed a few planets close to the star, generally fewer than with the slower migration.

Later work by Hansen & Murray (2013) combined N-body and Monte Carlo simulations. They suggest that in situ formation of super-Earth systems is possible and provides a reasonably good match to existing Kepler data. On the other hand, Raymond & Cossou (2014) argue against in situ formation and suggest that all close-in super-Earth planets must have reached their current orbits as a result of migration. Recently, Ogihara et al. (2015) studied the formation pathway of super-Earths in disks that do not allow outward migration. In their study, they let planetary embryos migrate toward the inner edge of the disk where they collided and formed larger objects. The resulting planetary systems were basically all tightly packed and with a hierarchical mass order (the most massive planet was closest to the central star), in contradiction to observations. This suggests that planet migration in a disk with no zones of outward migration may not be able to explain the distribution of observed super-Earth systems. Thus a change in the migration mechanism may be necessary to explain the observed distribution of super-Earths. The results of Cossou et al. (2014) and Coleman & Nelson (2014) are somewhat more promising, but all of these studies rely on planetesimal accretion and on migration to pile up the planets close to the central star. The formation of multiple planets undergoing pebble accretion is expected to have a great variety of outcomes, with the mass gradient of the planets sensitively depending on the temperature gradient of the disk (Ida et al. 2016). A full N-body study is beyond the scope of this paper.

Zones of outward migration in the protoplanetary disk can be caused by strong enough radial gradients in temperature, as we studied here, or by an inversion of the radial surface density profile, where the surface density increases with orbital distance (Masset et al. 2006). These changes in the disk profile can be achieved due to photoevaporation, which carves a hole in the inner disk, stopping planet migration, even for gas giants (Alexander & Pascucci 2012). However, these changes of the surface density gradients only appear in the very late stages of the disk evolution and are typically farther than 1 au, which severely limits the ability of this mechanism to halt planet migration and prevent super-Earth planets from reaching the central star.

Another topic that requires discussion is the effect of whether the pebble isolation mass depends on α; this was alluded to in Section 2.3. Assume that ${m}_{\mathrm{iso}}\propto {\alpha }_{3}^{\beta }$ and that this relation is valid for ${\alpha }_{3}\gtrsim 0.01$. We repeat that α₃ = α/10⁻³ and we explore the scaling law around this value. Then we have

$$\begin{eqnarray}&&{h}_{\mathrm{crit}-\alpha }=\displaystyle \frac{\pi \sqrt{3}}{128}{p}_{\nu }^{2}{\alpha }_{3}^{1-3/2\beta }.\end{eqnarray} \tag{ 17 }$$

We assume β ≪ 1. This choice is justified because Lambrechts et al. (2014) observed no appreciable difference in the isolation mass when they varied α by a factor of five. In an independent study, Zhu et al. (2014) conclude that in inviscid disks the estimate of the isolation mass is no more than a factor of two lower than the nominal value; a value β ≫ 1 appears inconsistent with these previous works.

As an example, we consider β = 1/5, then ${h}_{\mathrm{crit}}\propto {\alpha }_{3}^{7/10}$. As α₃ is varied from 0.1 to 10, the value of h_crit varies by as much as a factor five from the case where the isolation mass has no α dependence. Dividing the modified value of ${h}_{\mathrm{crit}-\alpha }$ in Equation (17) by the original value h_crit in Equation (11) yields ${h}_{\mathrm{crit}-\alpha }/{h}_{\mathrm{crit}}={\alpha }_{3}^{-3/2\beta }$. Thus the modified value ${h}_{\mathrm{crit}-\alpha }$ is higher than the original value in Equation (11) for α₃ < 1 and lower when α₃ > 1. This α dependence hinders trapping planets in the outward migration region: when α₃ > 1, the pebble isolation mass becomes higher than without an α dependence, and the trapping needs to occur at lower h (see Figure 4). When α₃ < 1, the trapping should in theory be more efficient because the zero-torque line is crossed at higher h. In Figure 15 we show the outcome of simple pebble accretion simulations in an evolving disk with varying opacity, but now the pebble isolation mass scales as ${m}_{\mathrm{iso}}\propto {\alpha }_{3}^{1/5}$. In the top panel no planets survive, which is the same outcome as in Figure 12. In the bottom panel there is a clear dependence on α, and the higher pebble isolation mass for larger α makes trapping more difficult in the evolving disk than in the case with no α dependence (Figure 12). Note also the slightly different color scale and maximum isolation mass, confirming our expectations.

Zoom In Zoom Out Reset image size

If there is no α dependence on m_iso, then it becomes nearly impossible to trap low-mass planets in the simplified model that we have considered if the temperature gradient is shallow. A viscous region with a steeper temperature gradient works much better, as does the full model because of the large region of outward migration near 3 au. On the other hand, if the weak α dependence is confirmed in future studies, then disks having α < 10⁻³ are still not expected to be able to trap planets far from the star.

A recent study of the global evolution of protoplanetary disks indicates that the midplane is laminar, with α likely to be lower than 10⁻³ (Bai 2016). These newer models also do not appear to agree very well with traditional constant α disk models driven by isotropic MRI turbulence, however. That said, much information is still to be had about planet migration with the more traditional approach and models that we have adopted here. An alternative model, relying on disk winds (Suzuki et al. 2016), could be explored in the future. In addition, we are not yet aware of a full treatment of the torques on the planets from such a disk.

Another possibility for trapping planets would be the inner edge of the dead zone, where changes in viscosity cause a positive gradient in the surface density, blocking the influx of dust and also the inward motion of planets (Masset et al. 2006; Flock et al. 2017). These changes in viscosity are caused by changes in the ionization fraction of the disk, which is a function of the disk density. However, as disks evolve, so will the inner edge of the dead zones, washing away the gradient of the surface density, releasing trapped planets (Bitsch et al. 2014), which is why we here instead focus on trapping planets purely by the corotation torque.

Unless there is an as yet unknown and undiscovered disk mechanism that acts on migrating planets, the entropy-driven corotation torque is the only way of stopping, and possibly reversing, the inward migration of super-Earth mass planets. For this to happen, the disk profile needs to have a radial temperature slope steeper than 0.87 and a viscosity higher than a few promille. Lower temperature slopes and viscosities in disks will only allow inward migration, leading to the problematic super-Earth formation scenario from Ogihara et al. (2015). It is not guaranteed that steeper temperature slopes and higher viscosities in disks will solve the formation mystery of super-Earths, but it should be investigated in future studies, especially those relying on pebble accretion.

We investigated the conditions that are required for low-mass planets to execute outward migration in steady-state protoplanetary disks. This outward migration is accomplished by the horseshoe torque, which is a complex function of the mass, scale height, viscosity, and thermal diffusivity of the disk. We restricted ourselves mostly to the disk model employed by Ida et al. (2016) because of its simplicity, but we allowed the temperature gradient and opacity to vary. We ultimately conclude the following:

• 1.  
  the corotation torque is at a maximum when the thermal diffusion constant $\chi =\tfrac{9}{4}\nu ;$
• 2.  
  the negative temperature gradient needs to exceed approximately 0.87 for outward migration to occur;
• 3.  
  there exists a critical value of the disk scale height, h_crit, below which planets at the pebble isolation mass are trapped in the outward migration region, preventing them from reaching the central star;
• 4.  
  planets at the pebble isolation mass can be prevented from reaching the star as long as α exceeds a critical value that depends on the temperature profile and age of the disk. A typical rule of thumb is that α_crit ≳ 4 × 10⁻³;
• 5.  
  opacity variations in the disk do not necessarily increase the region of phase space where planets can be trapped, but it may create additional regions where outward migration becomes possible.

We predict that systems with a high number of low-mass planets, such as TRAPPIST-1 (Gillon et al. 2017), whose planets may be in a multiple resonant chain, underwent slower migration than those with just a few planets (Ogihara & Ida 2009). We attribute the variation in migration speed from one system to the next to a different temperature gradient and opacity in the viscously heated part of their circumstellar disks.

The authors thank Michiel Lambrechts for a positive and constructive review that substantially improved this paper. R.B. is grateful for financial support from JSPS KAKENHI (16K17662).

Here we briefly summarize our implentation of pebble accretion.

The timescale of growing planetary seeds from planetary embryos all the way to planetary cores can be greatly enhanced by considering the accretion of leftover pebbles in the disk (Ormel & Klahr 2010; Lambrechts & Johansen 2012). In contrast to planetesimals, small pebbles feel gas drag and migrate through the disk. If these pebbles enter within the Hill sphere of the planet, they spiral toward the planet due to their loss of angular momentum caused by the friction with the gas and are subsequently accreted.

Very small bodies (100–200 km in diameter) accrete pebbles in the inefficient Bondi accretion branch, see Lambrechts & Johansen (2012), where the accretion timescale is quite long (depending on the disk structure, this can be in the range of 1 Myr) until the pebble transition mass is reached

$$\begin{eqnarray}&&{M}_{{\rm{t}}}=\displaystyle \frac{\sqrt{3}}{3}\displaystyle \frac{{(\eta {v}_{{\rm{K}}})}^{3}}{G{{\rm{\Omega }}}_{{\rm{K}}}},\end{eqnarray} \tag{ 18 }$$

where G is the gravitational constant, v_K = Ω_Kr, and

$$\begin{eqnarray}&&\eta =-\displaystyle \frac{1}{2}\displaystyle \frac{\partial \mathrm{ln}P}{\partial \mathrm{ln}r}{h}^{2}.\end{eqnarray} \tag{ 19 }$$

Here, $\partial \mathrm{ln}P/\partial \mathrm{ln}r$ is the radial pressure gradient in the disk; generally, $\partial \mathrm{ln}P/\partial \mathrm{ln}r=3-\tfrac{1}{2}q$. These masses are typically in the range of 0.01 M_⊕ and it marks the mass at which planets can accrete within the framework of the efficient Hill accretion. The accretion rate is additionally a function of the pebble scale height given by

$$\begin{eqnarray}&&{H}_{\mathrm{peb}}={H}_{{\rm{g}}}\sqrt{\displaystyle \frac{\alpha }{{ \mathcal S }}},\end{eqnarray} \tag{ 20 }$$

where α is the viscosity parameter and ${ \mathcal S }$ the Stokes number. When the Hill radius of the planet is larger than the pebble scale height, the accretion rate of the planet is given by (Lambrechts & Johansen 2014)

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{c}},2{\rm{D}}}=2{\left(\displaystyle \frac{{ \mathcal S }}{0.1}\right)}^{2/3}{r}_{{\rm{H}}}{v}_{{\rm{H}}}{{\rm{\Sigma }}}_{\mathrm{peb}},\end{eqnarray} \tag{ 21 }$$

where ${r}_{{\rm{H}}}=r{[{M}_{{\rm{c}}}/(3{M}_{\star })]}^{1/3}$ is the Hill radius, v_H = Ω_Kr_H the Hill speed, and Σ_peb the pebble surface density. If the Stokes number of the particles exceeds 0.1, the accretion rate is limited to

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{c}},2{\rm{D}}}=2{r}_{{\rm{H}}}{v}_{{\rm{H}}}{{\rm{\Sigma }}}_{\mathrm{peb}},\end{eqnarray} \tag{ 22 }$$

because the planetary seed cannot accrete particles from outside its Hill radius (Lambrechts & Johansen 2012). However, if the planetary seeds are small (as at the beginning of growth), the planet's Hill radius is smaller than the pebble scale height of the disk, meaning that the planet is fully embedded in the flow of pebbles and the accretion rate changes to (Morbidelli et al. 2015)

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{c}},3{\rm{D}}}={\dot{M}}_{{\rm{c}},2{\rm{D}}}\left(\displaystyle \frac{\pi {({ \mathcal S }/0.1)}^{1/3}{r}_{{\rm{H}}}}{2\sqrt{2\pi }{H}_{\mathrm{peb}}}\right).\end{eqnarray} \tag{ 23 }$$

The transition from 3D to 2D pebble accretion is then reached when

$$\begin{eqnarray}&&\displaystyle \frac{\pi {({ \mathcal S }/0.1)}^{1/3}{r}_{{\rm{H}}}}{2\sqrt{2\pi }}\gt {H}_{\mathrm{peb}}.\end{eqnarray} \tag{ 24 }$$

This transition depends on the scale height of the disk and on the particle size, indicating that a higher planetary mass is needed to reach the faster 2D accretion branch in the outer parts of the disk where H_peb is larger. The Stokes number of the dominant particle size given by

$$\begin{eqnarray}&&{ \mathcal S }=\displaystyle \frac{\sqrt{3}}{8}\displaystyle \frac{{\epsilon }_{{\rm{P}}}}{\eta }\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{peb}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}.\end{eqnarray} \tag{ 25 }$$

Here we set the parameters _P to 0.5 and _D to 0.05 as in Lambrechts & Johansen (2014). The final Stokes number obtained is constrained by the equilibrium between growth and drift to fit advanced coagulation models (Birnstiel et al. 2012). The pebble surface density depends on the gas surface density Σ_g and the semimajor axis r_p through

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{peb}}=\sqrt{\displaystyle \frac{2{\dot{M}}_{\mathrm{peb}}{{\rm{\Sigma }}}_{{\rm{g}}}}{\sqrt{3}\pi {\epsilon }_{{\rm{P}}}{r}_{{\rm{P}}}{v}_{{\rm{K}}}}},\end{eqnarray} \tag{ 26 }$$

where the pebble flux is

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{peb}}=2\pi {r}_{{\rm{g}}}\displaystyle \frac{{{dr}}_{{\rm{g}}}}{{dt}}({Z}_{\mathrm{peb}}{{\rm{\Sigma }}}_{{\rm{g}}}).\end{eqnarray} \tag{ 27 }$$

Here Z_peb denotes the fraction of solids (metallicity) in the disk that can be transformed into pebbles at the pebble production line r_g at time t

$$\begin{eqnarray}&&{r}_{{\rm{g}}}={\left(\displaystyle \frac{3}{16}\right)}^{1/3}{({{GM}}_{\star })}^{1/3}{({\epsilon }_{{\rm{D}}}Z)}^{2/3}{t}^{2/3},\end{eqnarray} \tag{ 28 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{{dr}}_{{\rm{g}}}}{{dt}}=\displaystyle \frac{2}{3}{\left(\displaystyle \frac{3}{16}\right)}^{1/3}{({{GM}}_{\star })}^{1/3}{({\epsilon }_{{\rm{D}}}Z)}^{2/3}{t}^{-1/3},\end{eqnarray} \tag{ 29 }$$

where M_⋆ is the stellar mass, which we set to 1 M_⊙.

This pebble growth mechanism in combination with planet migration and disk structure evolution was used in Bitsch et al. (2015b) to constrain the growth of planets in disks, and the same code is used here as well. We just adapt it to the different disk models and viscosities, but the pebble accretion model stays the same.