© T. Baumann and B. Bitsch 2020

1 Introduction

The detection of the first exoplanet around a solar-type star (Mayor & Queloz 1995) has revealed a new class of planets that we do not have in our own solar system: hot Jupiters. These planets are similar in mass to our own Jupiter, but their orbits around their central star last just shy of a few hours or days. However, at positions so close to the star, the sublimation temperatures for silicates can be exceeded, making in-situ formation impossible, which raises the question of how these planets got into this position.

Gravitational interactions between planets and the protoplanetary disk lead to a change in the semi-major axis of the planets, that is, migration. Spiral density waves, launched at Lindblad resonances, lead to fast inwards migration (Ward 1997). Studies in recent years have also found so-called corotation torques (Tanaka et al. 2002; Baruteau & Masset 2008; Paardekooper & Papaloizou 2009), which can slow down migration or foster outwards migration (Paardekooper & Mellema 2006; Kley & Crida 2008; Kley et al. 2009) in certain disk conditions and regions (e.g., Bitsch et al. 2013).

Planetary migration is an important ingredient for planet formation simulations (Ida & Lin 2004; Mordasini et al. 2009; Bitsch et al. 2015a; Ndugu et al. 2018). However, planet formation simulations cannot incorporate 2D or 3D hydrodynamical simulations to accurately calculate the torques over millions of years. Thus, the planet formation simulations rely on torque formulae derived from hydrodynamical simulations, such as Masset & Casoli (2010), Paardekooper et al. (2011), or Jiménez & Masset (2017).

Two such models of torque formulae that were derived from fits of hydrodynamical simulations are studied and compared in this work. These are Paardekooper et al. (2011), which was derived from a two-dimensional disk, and Jiménez & Masset (2017), which is updated research in a three-dimensional disk. A new development, which is also investigated here, is the thermal torque (Benítez-Llambay et al. 2015; Masset 2017), which can lead to outwards migration due to heat released by the planet.

The migration rate of planets crucially depends on the radial gradients of the gas surface density and the entropy. This implies that the migration rates depend significantly on the underlying disk structure and disk model. As the planet grows, it can migrate across several ice lines, which would change the composition of the material the planet can accrete (Öberg et al. 2011a; Madhusudhan et al. 2017; Cridland et al. 2017).

Bitsch et al. (2015a) presents a study of planet formation in the pebble accretion scenario in the Bitsch et al. (2015b) disk model, with planets migrating according to the Paardekooper et al. (2011) formula. One aim of our work here is to determine if the differences between the new formulae from Jiménez & Masset (2017) are large compared to the migration rates by Paardekooper et al. (2011). Additionally, we investigate the influence of the heating and cooling torque (Masset 2017) on the final mass and position of planets. Here, the relevant quantity is the solid accretion rate onto the planet, which releases the heat needed for the heating torque to work (Benítez-Llambay et al. 2015).

A previous study by Guilera et al. (2019) also studied the heating torque in planet formation simulations and found it had a strong influence on planet migration and thus formation. Here, we expand on their approach and include a self-consistent calculation of the heating torque within the pebble accretion scenario. In particular, we investigate the influence of the pebble flux on outwards migration triggered by the heating torque.

That said, our results are only applicable to simulations founded on similar assumptions. A fundamental simplification of this paper is the modeling of single planet evolution only. N-body simulations (e.g., Cossou et al. 2013; Izidoro et al. 2019; Bitsch et al. 2019a), on the other hand, exhibit many additional influence factors which are not captured in our approach. For example, the trapping in resonances depends crucially on the relative migration speed between the planets. So if planets suddenly migrate slightly slower or even outwards, trapping in resonances might operate in a different way compared to situations where planets only migrate inwards.

The structure of this paper is as follows: the methods are listed in Sect. 2, followed by a comparison of the torque formulae from Paardekooper et al. (2011) and Jiménez & Masset (2017) in Sect. 3. The influence of the thermal torque is analyzed in Sect. 4. Finally, the analysis of the previous sections is repeated in the disk model from Ida et al. (2016) in Sect. 5. A discussion of the results and the conclusion can be found in Sects. 6 and 7.

2 Planetary growth and migration models

The model of planetary growth implemented here has already been described in Bitsch et al. (2015a), so we simply repeat the essentials of the growth model. We also use the disk model of Bitsch et al. (2015b), which originates from 3D hydrodynamical simulations, including stellar and viscous heating, as well as radiative cooling by micrometer-sized dust grains.

In this paper, we assume planetary cores to grow via pebble accretion and migrate in type I and type II fashion. Pebbles are rocky or icy objects of mm to cm size that form by coagulation or condensation of dust particles (Brauer et al. 2008; Birnstiel et al. 2012; Ros & Johansen 2013; Drążkowska & Alibert 2017; Ormel et al. 2017). By accreting these pebbles, planets can grow cores efficiently (Ormel & Klahr 2010; Johansen & Lacerda 2010; Lambrechts & Johansen 2012), that is, on timescales that are multiple orders of magnitude shorter than planetesimal accretion (Johansen & Lambrechts 2017; Johansen & Bitsch 2019). The size of the pebbles relates to the Stokes number (e.g., Brauer et al. 2008), which is set to 0.1 in these simulations.

In our simulations, planets are inserted at pebble transition mass, meaning they accrete pebbles within the Hill radius right away. Pebble transition mass is determined following Lambrechts & Johansen (2012, 2014).

Pebble accretion is modeled in two regimes. The planet accretes pebbles in the fast 2D scenario, when its Hill radius is larger than the pebble scale height, $$H_{\text{peb}}=H_{\text{gas}}\sqrt{\alpha /{\tau }_{\text{f}}},$$(1)

with the gas scale height H_gas, the α-viscosity parameter (Shakura & Sunyaev 1973) and the Stokes number τ_f (Youdin & Lithwick 2007; Morbidelli et al. 2015). The accretion rate in the 2D regime is given by: $${\dot{M}}_{\text{c},2\text{D}}=2r_{\text{H}}{v}_{\text{H}}{\Sigma }_{\text{peb}}.$$(2)

Here, r_H is the Hill radius, v_H = r_HΩ and Σ_peb is the pebble surface density at the planet’s location. Otherwise, the planet accretes in the 3D scenario. A reduction factor from Morbidelli et al. (2015) is applied, resulting in an accretion rate of: $${\dot{M}}_{\text{c},3\text{D}}={\dot{M}}_{\text{c},2\text{D}}\left(\frac{\pi {\left({\tau }_{\text{f}}/0.1\right)}^{1/3}{r}_{\text{H}}}{2\sqrt{2\pi }{H}_{\text{peb}}}\right).$$(3)

The amount of pebbles that are available to the planet is determined by the pebble flux through the disk, which is set to ${\dot{M}}_{\text{peb}}=2\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$ in most simulations. This corresponds to a total of about 380 M_E (Earth-masses) of pebbles in the disk throughout its lifetime of 3 Myr, equivalent to roughly a 1% dust-to-gas ratio if the initial disk mass was 10% of the solar mass.

During core growth, the planet accretes 90% in solids and 10% in gas, building a gaseous envelope alongside the core. When the planet reaches pebble isolation mass (Lambrechts et al. 2014; Bitsch et al. 2018; Ataiee et al. 2018), core growth is halted and the gaseous envelope is accreted onto the planet. Following Bitsch et al. (2018), the pebble isolation mass is calculated via: $$\begin{array}{lll}\hfill & \hfill & M_{\text{iso}}=25f_{\text{fit}}{M}_{\text{E}}+\frac{{\Pi }_{\text{crit}}}{\lambda }{M}_{\text{E}},\hfill \\ \hfill & \hfill & f_{\text{fit}}={\left[\frac{H/r}{0.05}\right]}^3\left[0.34{\left(\frac{\log \left({\alpha }_3\right)}{\log \left({\alpha }_{\text{disc}}\right)}\right)}^4+0.66\right]\left[1-\frac{\frac{\partial \ln P}{\partial \ln r}+2.5}{6}\right],\hfill \end{array}$$

with $\lambda \approx 0.00476/f_{\text{fit}},{\Pi }_{\text{crit}}=\frac{{\alpha }_{\text{disc}}}{2{\tau }_{\text{f}}}$ and α₃ = 0.001.

After the pebble isolation mass is reached, the heating of the envelope through in-falling pebbles stops and the planetary envelope can contract. Here, we follow the contraction rates from Piso & Youdin (2014). Once the envelope mass is similar to the core mass, runaway gas accretion can start and for this, we follow the rates of Machida et al. (2010).

Planetary growth is halted when the planet reaches the inner edge of the disk at 0.1 AU or when the gas disk dissipates after 3 Myr. We stop the growth of planets at 0.1 AU because gas accretion is thought to be inefficient in this region due to strong recycling flows that penetrate deeply into the Hill sphere and prevent gas contraction around young planets (e.g., Cimerman et al. 2017; Lambrechts & Lega 2017).

Migration is modeled in two regimes, called type I and II. Type I migration is a consequence of gravitational interaction between the planet and the protoplanetary disk that it grows in. The two models for type I migration that are used here are described in Sect. 3, with the formulae in Appendices A–C. Type I migration is valid for lowmass planets. As planets grow more massive, they start opening a gap in the disk and transition to type II. In this regime, migration is determined by stellar accretion and is usually slower than type I and directed towards the star. For a recentreview on type I and type II migration, see Kley & Nelson (2012) and Baruteau et al. (2014). For information about the calculation of type II migration timescales in this paper, see Appendix. D.

A simplification in these simulations is that the planets move in circular orbits with vanishing inclinations. This is a sensible approximation, as eccentricity and inclination are damped quite quickly by the gas disk (Bitsch & Kley 2010, 2011).

3 Comparison of the torque formulae from Paardekooper et al. (2011) and Jiménez & Masset (2017)

Paardekooper et al. (2011) and Jiménez & Masset (2017) both performed hydrodynamical simulations in non-barotropic disks with finite viscosity and constant thermal diffusion to derive their formulae. The major difference is that Paardekooper et al. (2011) made two-dimensional simulations, while Jiménez & Masset (2017) updated their findings based on three-dimensional disks.

Generally, the two formulae agree in the physical principles of the torque exchange, but show different numerical factors. For example, the saturation mechanisms respond similarly to viscosity and thermal diffusivity changes, but for individual values of viscosity, the behavior can be significantly different.

The Paardekooper et al. (2011) formula incorporates fewer components. The components that they included are: the Lindblad torque, which drives inwards migration due to spiral density waves at Lindblad resonances, and the barotropic and entropy related corotation torques. The corotation torques typically slow down or even set inwards migration in reverse. The entropy torque is based on the entropy gradient, that is, a combination of the temperature and surface density gradients in the disk. The second contribution is the barotropic torque, resulting from the surface density gradient.

All of the above components are incorporated in Jiménez & Masset (2017) as well, although the barotropic torque is called vortensity torque in their formulation. Additionally, there is the temperature torque as another contribution of the midplane temperature gradient. Finally, Jiménez & Masset (2017) introduced the viscous coupling term, a result of viscous production of vortensity at the density jumps at the seperatrices of the horseshoe region.

We ought to keep in mind that these formulae describe type I migration. Once the planet transitions to type II, its migration is determined by the disk properties and differences in these formulae no longer impact migration.

Type II migration is also an active area of research, where it is currently debated if the type II migration speed follows the viscous flow (e.g., Dürmann & Kley 2015, Robert et al. 2018) or is a scaling of the type I torques with the gap depth (e.g., Kanagawa et al. 2015, 2018; Ida et al. 2019; Johansen et al. 2019). Here, we use the viscous speeds for planets in type II migration.

Fig. 1 Individual torque components of the Paardekooper et al. (2011) and Jiménez & Masset (2017) formulae for planets at 3 AU and 1 Myr, depending on planetary mass. The Jiménez & Masset (2017) formula adds new contributions, but the total torque looks similar. This is just one arbitrarily chosen set of disk evolution time and position in the disk, which manifests many general trends. Outwards migration occurs when the total torque is positive. This plot shows outwards migration caused by the heating torque (see Sect. 4) for small planets as well as outwards migration caused by the corotation torques for slightly heavier planets. The pebble flux through the disk is ${\dot{M}}_{\text{peb}}=2\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$, which is relevant for the thermal torque as discussed in Sect. 4. The viscosity is set to α = 5.4 × 10−3 during these simulations.

3.1 Differences in the formulae

As the formulae are based on the same physics, we expect similar outcomes. This is true for planet formation, but migration maps can show significant differences. Figure 1 shows the individual torque components of the two formulae for planets of different masses at an orbital radius of 3 AU and 1 Myr of disk evolution time. These plots already include the thermal torque from Masset (2017), which is discussed in Sect. 4. The plot is an arbitrary example of disk evolution time and position in the disk, but serves to show many general trends.

The torques shown in Fig. 1 are normalized to: $${\Gamma }_0={\left(\frac{q}{h}\right)}^2{\Sigma }_{\text{P}}{r}_{\text{P}}^4{\Omega }_{\text{P}}^2,$$(6)

where the index P indicates that the quantities are evaluated at the planets’ position. In this formulation the Lindblad torque does not depend on the planetary mass, as this dependence is hidden in Γ₀. The torque is always negative in both formulae, but we find up to 35% stronger Lindblad torques with the Paardekooper et al. (2011) formula, where the largest differences occur in the inner disk.

The barotropic or vortensity torque component is very comparable between the formulae, but the entropy torque shows a large difference. It is the main source for outwards migration in the Paardekooper et al. (2011) formula, but is much smaller in magnitude and often even negative in the Jiménez & Masset (2017) formula.

The temperature torque in Jiménez & Masset (2017), which does not exist in Paardekooper et al. (2011), gives an additional positive contribution to the torque. However, it is usually not large enough to make up for the difference in entropy torques compared to Paardekooper et al. (2011). The viscous coupling term, the other new component in the Jiménez & Masset (2017) formula, is found to be negligible in our simulations.

The horseshoe width, a region where the material moves in horseshoe orbits around the planet, is the source of the corotation torque and is, consequently, an important parameter. This width, however, is calculated differently in Jiménez & Masset (2017) and Paardekooper et al. (2011). In fact, Jiménez & Masset (2017) have argued that the Paardekooper et al. (2011) formula is only valid for planets of up to a few Earth-masses. Jiménez & Masset (2017) find slightly smaller horseshoe widths at higher masses. This should lead to smaller corotation torques, but we find that this difference is not that large due to the transition into type II migration for larger planets, which changes the torque in any case.

To summarize, the Paardekooper et al. (2011) torque formula allows a larger positive contribution from the entropy related corotation torque compared to the Jiménez & Masset (2017) torque formula. This leads to larger regions of outwards migration using the Paardekooper et al. (2011) formula compared to the Jiménez & Masset (2017) model (see Fig. 4).

Fig. 2 Planet formation maps with different pebble fluxes in the disk. The axes display where and when a planetary embryo has been inserted in the disk. The color codes how massive the planets grow, while the solid lines outline its final orbitalradius. Planets that are below the dashed blue line have reached pebble isolation mass and accrete gas. Planets that are interior to the dashed white line have reached runaway gas accretion and can grow very massive. Ṁpeb marks the pebble flux through the disk, which decreases exponentially over time. Differences between the formulae in terms of impact on planet formation are small. As the pebble density in the outer disk is small, planets that are inserted there do not grow much.

3.2 Impact on planet formation

As hinted at before, differences in planet evolution between the two formulae are small. Figure 2 shows the final mass and orbital radius of planets, depending on where (r₀) and when (t₀) they have been inserted in the disk. When comparing the four plots, it becomes apparent that the two formulae produce planets of similar mass and final orbital radius, regardless of the amount of pebbles in the disk. We only observe a difference in the very inner disk, where planets migrating with the Paardekooper et al. (2011) formula are less massive. This is caused by the different migration pattern that moves the planets faster into the inner disk, where their pebble isolation mass is small, resulting in smaller planets and thus prevents them from reaching runaway gas accretion.

It is noteworthy that this result only holds for single planet formation. This is because the migration maps do show significant differences (Fig. 4), which might be of larger influence in an N-body simulation, as they introduce other sources of migration. The main ones are gravitational interactions between the planets and trapping in mean motion resonances.

The trapping in mean motion resonances depends crucially on the relative migration speed between the planets, where a slower relative migration speed allows trapping in wider resonances compared to faster migration speeds. The resonant configuration of the system then has important consequences for the stability of the formed planetary system and thus its final configuration (e.g., Izidoro et al. 2017, 2019).

The mutual interactions between the planets could increase the eccentricity of the planets, which quenches the entropy-driven corotation torque and stops outwards migration (e.g., Bitsch & Kley 2010; Fendyke & Nelson 2013). The impact of forming multiple planets simultaneously is to be tested in future research.

Fig. 3 Planet formation maps with different viscosities for planetary migration using the Jiménez & Masset (2017) formula. Like in Fig. 2, planets that are below the dashed blue line havereached pebble isolation mass and accrete gas. Planets that are interior to the dashed white line have reached runaway gas accretion and can grow very massive. The pebble flux is ${\dot{M}}_{\text{peb}}=2\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$ in all plots. A lower choice of α mostly leads to more massive planets due to the quicker transition to type II migration. In addition, the inwards migration of planets is reduced due to the slow type II migration.

3.3 Influence of viscosity

Viscosity has proven to be a large influence factor on migration, because viscosity is responsible to keep the horseshoe region unsaturated (e.g., Masset 2001; Baruteau & Masset 2008; Paardekooper & Papaloizou 2009; Bitsch et al. 2013). At low viscosity, the corotation torque saturates and outwards migration from the corotation torques ceases to exist.

Our simulations incorporate the α-viscosity model from Shakura & Sunyaev (1973). However, we keep the disk structure fixed in our model to show the differences on the torques, even though the disk structure is a combination of viscous and stellar heating as well as radiative cooling (Bitsch et al. 2015b), which would change if the viscosity changes. When the viscosity was changed in these simulations, only migration timescales and the gap parameter (this parameter models the transition to type II migration, see Eq. (D.1)) were calculated using different values of α. In addition, we assume that there is no influence on the accretion rate by the reduced viscosity, in order to study the effects of lower viscosity on the migration rates separately.

The influence of viscosity is mainly due to two effects. Type I migration is influenced, as corotation torques saturate more quickly at low viscosity. This leads to less outwards migration at lower α-values. On the other hand, the gap opening is also changed, meaning planets transition to type II at lower masses for lower viscosity.

The impact of a viscosity change is shown in Fig. 3, which shows final mass and orbital radius maps using the Jiménez & Masset (2017) formula. Clearly, a choice of lower α leads to more massive planets overall, which can be explained by the earlier transition to type II migration. The α = 10⁻³ plot also shows planets that do not grow as large when inserted into the inner disk, which is a consequence of the saturation of the corotation torques. The planets, in this case, grow and migrate to the inner edge of the protoplanetary disk where they stop their migration and growth in our model. Thus, these planets do not grow into gas giants.

The trends for changes in viscosity are applicable also to the Paardekooper et al. (2011) formula so we only show the results of the Jiménez & Masset (2017) formula here. Individual migration maps, however, show slightly different reactions to a change in viscosity. This is shown in Fig. 4, where migration maps with both formulae at 1 Myr are shown for different values of α. As mentioned before, the regions of outwards migration shrink with lower viscosity, until they disappear for α = 10⁻⁴, and the planets transition into type II migration at lower masses at low viscosity, since they can open gaps more easily.

The saturation functions of the corotation torque of the two models are based on the same physical requirements, but are slightly different numerically, similar to the torque components. For a plot of the Paardekooper et al. (2011) saturation function see Figs. 2 and 3 in Paardekooper et al. (2011).

Fig. 4 Migration maps after 1 Myr with different α-viscosity values for the calculation of migration timescales and the gap parameter. The left plots show the Paardekooper et al. (2011) formula and the right ones the Jiménez & Masset (2017) formula. Hatched regions mark the transition to type II migration, while planets in the cross-hatched areas are migrating completely in the type II regime. Regions of outwards migration that are caused by the thermal torque, which is calculated with an accretion rate of ${\dot{M}}_{\text{peb}}\approx 1.4\times {10}^{-4}\frac{M_{\text{E}}}{\text{year}}$, are outlined by dashed blue lines, while regions of outwards migration inside solid blue lines are due to the original torque formulae. Outwards migration in the original formulae is caused by the torques originating from the corotation region (Paardekooper et al. 2011) or by the temperature torque (Jiménez & Masset 2017). As type II migration depends on the disk properties rather the type I torque formula, the cross-hatched regions are identical between the corresponding left and right plots.

4 Thermal torque from Masset (2017)

A recent development in type I migration is the thermal torque, which is a torque that acts on the planet due to its own heat release. Two manifestations of this phenomenon have been discovered in the simulations of Lega et al. (2014) and Benítez-Llambay et al. (2015).

Thermal diffusion leads to two cold and dense lobes in the gas disk around the planet, which are asymmetric (Masset 2017). The resulting torque acting on the planet is called cooling torque or sometimes “cold finger effect” and was found in Lega et al. (2014). It is negative and accelerates inwards migration.

The other part of the thermal torque, which has been found in Benítez-Llambay et al. (2015), is the heating component. Here, hot and underdense lobes form in the gas profile around luminous planets due to heat release. This torque, again arising from asymmetry of the lobes, is positive and can lead to outwards migration.

The thermal torque acts on planets of up to intermediate mass only. It gets cut off by two factors: When the planet reaches pebble isolation mass, solid accretion stops and the associated heat release vanishes. In our model, we assume no planetesimal accretion which could still be ongoing after pebble isolation mass is reached. However, the amount of impacting planetesimals at this stage is rather small (e.g., Alibert et al. 2018), so that not enough heat can be generated to fuel the heating torque.

The other cut-off mass is called thermal critical mass, which describes the point at which the planet stops generating excess internal energy outside the Bondi sphere, meaning the gas density is no longer perturbed. This is estimated in Masset (2017) as: $$M_{\text{crit},\text{thermal}}={\chi }_{\text{P}}{c}_{\text{s},\text{P}}/G,$$(7)

where χ_P is the thermal diffusivity, c_s,P = h ⋅ Ω_P ⋅ r_P is the local sound speed and G is the gravitational constant. The index P indicates that the quantities are evaluated at the planets’ position.

Our simulations use the description of the thermal torque developed by Masset (2017), which incorporates both the heating and cooling torque. The luminosity for the heating torque is calculated as the solid accretion luminosity of the planet following Chrenko et al. (2017). This paper focuses on how the heating torque reacts to pebble accretion as the core growth model. For details about calculations and formulae related to the thermal torque that are used here, see Appendix C.

4.1 Impact on migration

The contribution of the thermal torque is already shown in Fig. 1, which shows the individual torque components and in Fig. 4, which shows the migration maps. Figure 1 shows the typical course of the thermal torque (thick blue line). As it is a composition of the negative cooling component and the positive heating part, it acts on planets with both signs. Low-mass planets typically have a small accretion luminosity, so the cooling effect is the larger contribution. Since the Hill radius scales with the third root of the core mass, planets can accrete enough pebbles for a positive thermal torque only when they get massive enough. Finally, the thermal torque vanishes at thermal critical mass or at pebble isolation mass. The thick purple lines show how the thermal torque changes the total torque according to the discussed behavior (the solid line includes the thermal torque, while the dashed one does not).

The migration maps in Fig. 4 are only influenced by thermal torque below the two red lines which mark the cut-off. Outwards migration due to the thermal torque is marked by dashed blue lines. However, these new regions of outwards migration are small and restricted to low mass planets between 2 and 5 AU.

The thermal torque is the same across all choices of α, because the thermal component is independent of alpha, as we do not change the pebble scale height when we change alpha. Instead, accretion rates were always calculated using α = 5.4 × 10⁻³. In a more realistic model, a lower viscosity reduces the pebble scale height, which results in larger accretion rates and thus a larger heating torque, potentially leading to larger regions of outwards migration due to the heating torque. When we tested this dependence,however, we found that the impact on planet formation is still small.

4.2 Impact on planet formation

Similarly to the comparison between Paardekooper et al. (2011) and Jiménez & Masset (2017), migration maps show a significant change by the thermal torque, while the impact on planet formation is essentially negligible.

Figure 4 shows one of the main reasons why the thermal torque did not lead to major changes in planet formation inour model: the thermal critical mass cuts the thermal torque off at low masses in the inner disk. This means that while outwards migration due to the heating torque does occur, most planets do not experience it during their formation.

However, this is only one part of the story, as the cut-off is much higher in the outer disk (≈ 10-100 M_E). However, the thermal torque is inversely proportional to the thermal diffusivity, which typically increases when moving away from the star. This means that while it is active, the thermal torque is small in the outer disk. Additionally, pebble accretion is less efficient there, as the pebble surface density decreases and the pebble scale height increases with radius. This means that the cooling torque outweighs the heating torque for planets that are far away from the star for the accretion rates used in this model, resulting in a low overall contribution of the heating torque to planet formation. We discuss implications of different accretion rates in the next section.

4.3 Thermal torque and pebble accretion

As the heating torque is directly proportional to the solid accretion rate (see Eqs. (C.5) and (C.6)), the corresponding model for core growth plays a crucial role in determining the thermal torque. The efficiency of pebble accretion depends on the amount of pebbles that are available to the planet. This is regulated via the pebble flux Ṁ_peb through the disk, which is typically set to ${\dot{M}}_{\text{peb}}=2\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$ (i.e., 380 M_E in total over 3 Myr) in our model. As discussed before, this is insufficient to yield a significant impact of the thermal torque on planet formation.

Research is typically directed towards slowing down migration to the star. Therefore, it is interesting to look for positive contributions of the thermal torque. Figure 5 shows the minimum solid accretion rate and corresponding pebble fluxes that are required for a positive thermal torque at a given planetary mass and orbital distance using our disk model. Keep in mind that this does not necessarily lead to outwards migration, as the Lindblad torque still needs to be overcome. The plot extends beyond the thermal critical mass, as this is only an estimate in Masset (2017) and regions of larger masses might become important in the future.

Clearly, Fig. 5 shows that unreasonably high pebble fluxes are required to achieve a meaningful positive contributionof the thermal torque. However, it is still interesting to look at studies with increased amounts of pebbles, which can be observed in Fig. 6. Here, planet evolution maps are shown with a pebble flux of ${\dot{M}}_{\text{peb}}=5\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$ (i.e., about 950 M_E in total, corresponding to a factor of 2.5 increase in pebbles to our nominal value), with and without the thermal torque. While the contribution is still not large, it is noticeable at this pebble flux. Only planets that are inserted in the inner disk are affected, however. This is due to the thermal torque being small in the outer disk and being cut-off for larger planets in the inner disk, meaning that planets that are inserted in the outer disk avoid regions of outwards migration caused by the heating torque during their growth.

When including the thermal torque, more planets reach runaway gas accretion. This happens, for example, when they form at an early stage (e.g., before 500 kyrs). This is caused by some outwards migration during the growth of the planets which has two effects: (i) as the planets move outwards, the pebble isolation mass increases so they can form larger cores which are more efficient at attracting gaseous envelopes; and (ii) as the planets are further away from the central star they need more time to migrate to the inner edge of the disk and thus have more time to grow.

Figure 7 shows select growth tracks with and without the thermal torque. The planets have been inserted after 1 Myr in a disk with a pebble flux of ${\dot{M}}_{\text{peb}}=5\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$ at various distances to the star. Comparing it to Fig. 6, we can see that this insertion scenario corresponds to planets that do not grow as much when including the thermal torque. This is the case despite significant outwards migration and serves to show that planetary growth is subject to many variables. Here, the planets spend less time in the regions of outwards migration created by corotation torques, since the initial outwards migration changes their trajectory in the migration maps and increases their pebble isolation mass, leading to fast inwards migration, where the planets have not enough time to contract their envelope and grow into gas giants before reaching the inner edge of the disk.

Fig. 5 Minimum core accretion rate and pebble flux that are required for a positive thermal torque in a disk that is 1 Myr old. The flux of pebbles needed to make a positive thermal torque also depends on the choice of Stokes number, which is set to 0.1 here. For a higher Stokes number, lower pebble fluxes are required for a positive thermal torque. The thermal critical mass used here is only an estimate, so values above might be of interest in the future. Below the thermal critical mass and in the outer disk, unreasonably high pebble fluxes are required for a positive thermal torque. Due to the exponential decrease in pebble flux that is employed in these simulations, about 72% of the initial pebble flux are left after 1 Myr, which corresponds to $1.4\times {10}^{-4}\frac{M_{\text{E}}}{\text{year}}$ in the standard case. Keep in mind that a positive thermal torque does not automatically correspond to outwards migration. For different parameters of the disk age and how it influences the minimal pebble flux needed for a positive thermal torque, see Figs. E.1 and E.2.

Fig. 6 Planet formation maps with and without the thermal torque from Masset (2017) using the Jiménez & Masset (2017) formula for the other torque components. Simulations were performed with an increased pebble flux of ${\dot{M}}_{\text{peb}}=5\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$, i.e. about 950 ME over 3 Myr. The plot should be read just like Fig. 2. Planets in the inner disk are significantly affected by the thermal torque due to the larger amount of pebbles. Planets that are inserted in the outer disk are not noticeably influenced by the thermal torque.

5 Evaluation in the Ida et al. (2016) disk model

As discussed before, the disk model massively impacts planetary growth and migration. This means that statements about migration inthis paper are valid only in the disk models discussed here, rather than being generally applicable. Therefore, it is useful to repeat the study for different disk environments. Here, we discuss the implications of the heating torque when using the disk model from Ida et al. (2016). Different components of the Paardekooper et al. (2011) torque formula and their influence on super-Earth formation has been discussed in Brasser et al. (2018).

In order to adapt our code to the different disk, the temperature has been recalculated according to Ida et al. (2016), but all other variables remained the same. In particular, stellar luminosity, viscosity and pebble accretion are modeled identically to the Bitsch et al. (2015b) disk simulations.

Compared to the Bitsch et al. (2015b) disk, this is a more simple model that employs power laws. The temperature is seperated into two regimes: The inner disk, which is dominated by viscous heating and the outer disk, where stellar irradiation dominates the temperature profile. As the disk evolves, the transition to stellar irradiation moves inwards. The temperature is then calculated using: $$\begin{array}{lll}\hfill & \hfill & T_{\text{vis}}\simeq 200M_{*0}^{3/10}{\left(\alpha \times {10}^3\right)}^{-1/5}{\left({\dot{M}}_{*}\times {10}^8\right)}^{2/5}{\left(\frac{r}{1\text{AU}}\right)}^{-9/10}\text{K},\hfill \\ \hfill & \hfill & T_{\text{irr}}\simeq 150L_{*0}^{2/7}{M}_{*0}^{-1/7}{\left(\frac{r}{1\text{AU}}\right)}^{-3/7}\text{K},\hfill \\ \hfill & \hfill & T=\max \left(T_{\text{vis}},T_{\text{irr}}\right).\hfill \end{array}$$

We used M_*0 = 1.989 × 10³³ g for the stellar mass and the values from Table 1 in Bitsch et al. (2015b) for the luminosity of the star. The gas surface density is then calculated through the gas accretion rate Ṁ, which we assume to follow the same evolution as in the Bitsch et al. (2015b) disk model. The impact of the different disk model in terms of planet formation is quite significant, as we show in the following.

Fig. 7 Growth tracks of planets that have been inserted at different radii after 1 Myr with the Jiménez & Masset (2017) formula (no thermal torque) in black and the Jiménez & Masset (2017) + Masset (2017) formula (includes thermal torque) in red. The lines become dashed when the planets reach pebble isolation mass and are accreting gas. An increased pebble flux of ${\dot{M}}_{\text{peb}}=5\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$ was used to achieve a larger contribution from the thermal torque. After performing significant outwards migration, planets thatgrow with the thermal torque migrate to the star quickly and grow less massive than the planets that were evolved without thermal torque, which reach runaway gas accretion. This might be due to the fact that heavier planets migrate faster, which gives them less time to contract an envelope before reaching the inner edge of the disk.

5.1 Changes in migration maps

Migration maps after 1 Myr of disk evolution in the Ida et al. (2016) disk model are shown in Fig. E.3. In the Ida et al. (2016) disk, the Jiménez & Masset (2017) formula produces larger regions of outwards migration than the Paardekooper et al. (2011) formula at the standard viscosity with α = 5.4 × 10⁻⁴. This stands in contrast to the Bitsch et al. (2015b) disk, where the opposite is true for most cases.

In fact, outwards migration is generally more scarce in the Ida et al. (2016) disk than in the Bitsch et al. (2015b) model. This is particularly noticeable for the thermal torque, which yielded small regions of outwards migration in the Bitsch et al. (2015b) disk, but only leads to essentially negligible amounts of outwards migration here. The transition to type II migration is similar in the inner disk, where the temperature is set by viscous heating, but larger planetary masses are required in the outer disk.

5.2 Changes in planet formation

Figure E.4 presents planet evolution maps in the Ida et al. (2016) disk model. Planets that evolve in the Ida et al. (2016) model generally grow less massive than planets in the Bitsch et al. (2015b) disk when considering the same pebble flux. This is caused by the reduced region of outwards migration by the entropy related corotation torque. In the Bitsch et al. (2015b) model, the growing planets can be parked in this region of outwards migration and stay there until the end of the disk lifetime or they then slowly contract their envelope and grow to gas giants from this region. Both options are not possible in the Ida et al. (2016) model and planets instead migrate to the inner edge of the disk as super-Earths with a small gaseous envelope. Apart from that, behavior is similar between the two models.

When comparing the Paardekooper et al. (2011) and the Jiménez & Masset (2017) torque formulae, really similar planets are formed in these migration models and the thermal torque does not impact planet growth significantly at these pebble fluxes. However, as the thermal torque produces even less outwards migration in this disk model, the impact on planet formation is also smaller when considering increased pebble fluxes.

6 Discussion

Here, we discuss the shortcomings and implications of our model and compare it with a previous study that incorporated the heating torque in planet formation simulations (Guilera et al. 2019).

6.1 Comparison with Guilera et al. (2019)

As we mentioned earlier in this paper, our results are only valid for the exact model that we tested. This becomes clear when we compare our results to Guilera et al. (2019), where the authors performed a similar study, but based on a different disk model. Their results are very different, with the most important factor being that they found a very large level of impact on the part of the thermal torque. In fact, they found very significant outwards migration of small planets due to the heating torque (see Fig. 11 in Guilera et al. 2019), which strongly impacts their final mass.

Their model is different in two main ways: as mentioned before, they operate in a different disk model, which goes beyond power laws for temperature (etc.) as well, but is very different from the Bitsch et al. (2015b) model. In their disk, the Jiménez & Masset (2017) formula produces significantly less outwards migration, while migration maps employing the Paardekooper et al. (2011) formula look similar between the two disk models, that is, their work and ours. The second major difference is that Guilera et al. (2019) employ a fixed core accretion rate of ${10}^{-5}\frac{M_{\text{E}}}{yr}$ for the calculation of the heating torque. We, on the other hand, use the solid accretion rate that follows from pebble accretion, which is usually significantly lower for small planets.

Guilera et al. (2019) come to the conclusion that the heating torque changes planetary growth and migration significantly. We argue that they overestimated the impact of the heating torque dramatically when they employed the fixed core accretion rate of ${10}^{-5}\frac{M_{\text{E}}}{yr}$ for planetary embryos. This can be seen in Fig. 5, where we show the minimal pebble flux and core accretion rate required to achieve a positive contribution of the thermal torque. In fact, in order to achieve the accretion rates used in Guilera et al. (2019), we would need a pebble flux that is more than a factor of 10 higher than our highest used value. We deem this high pebble flux unrealistic for most stellar systems.

6.2 Limitations of our model

The planet growth model is very important for determining the planets’ evolution. As the heating torque is proportional to the solid accretion rate, the growth mechanism gains additional importance in our study of the thermal torque. We only tested pebble accretion here, whereas a more realistic model would allow the planet to accrete planetesimals as well. That said, for planetary embryos heavier than 10⁻² M_E, pebble accretion is a lot more efficient than planetesimal accretion (Johansen & Lambrechts 2017), so the inclusion of planetesimals is unlikely to change our results significantly. Johansen & Bitsch (2019) showed that planetesimal accretion can only be efficient in the inner disk, if planetesimals are small and stirring is low. This implies that planetesimal accretion can not be the main contribution to core accretion of giant planets at a few AU, as we investigate here.

Additionally, the fact that the solid accretion rate depends on the α-viscosity via the pebble scale height (see Eqs. (1) and (3)) was neglected. At lower α, the pebble scale height decreases and planets enter the more efficient 2D-pebble accretion regime more quickly. Additionally, the reduction factor in the three dimensional accretion rates scales with α^−1∕2. That means at the lowest value that we tested (α = 10⁻⁴), solid accretion rates in the 3D-regime are larger by a factor of 10 compared to the highest value that we tested (α = 10⁻²). However, when we implemented this dependence we found that while planets grow significantly more massive and outwards migration from the heating torque was stronger, the influence of the thermal torque on planet evolution was still small.

The main reason why the thermal torque did not change migration much in our simulations was that the thermal critical mass cuts it off at low masses in the inner disk. It is important to note that this thermal critical mass was only an estimate of when the planet stops generating excess internal energy outside the Bondi sphere (Masset 2017). Should future research show that the thermal critical mass is in fact larger than we assumed here, our study of the thermal torque would have to be repeated.

6.3 Influence in N-body simulations

A more impactful simplification is the evolution of single planets only. As mentioned before, planet–planet interactions have a huge impact on planetary migration and are not captured here.

In particular, even though the heating torque has only minimal consequences on the overall formation of planets in our systems, the migration directions and speeds can change significantly (Fig. 7). This has important consequences for the trapping of planets in resonances, because the relative migration speed between the planets determines in which resonance the planets will be trapped. It is extremely difficult to asses the general implications of the heating torque for N-body simulations with our current simulations. In order to determine the influence of the torque formulae and the thermal torque on super-Earth systems, where the trapping in resonances plays an important role (e.g., Izidoro et al. 2017, 2019; Lambrechts et al. 2019), N-body simulations featuring these different prescriptions have to be undertaken.

6.4 Influence on the chemical composition of planets

When a planet crosses an ice line in the disk, the chemical composition of the accreted material changes (e.g., Öberg et al. 2011b; Madhusudhan et al. 2017). In particular, Bitsch et al. (2019b) showed that the direction of migration close to the ice lines has important consequences for the chemical composition of formed planets.

The heating torque has the potential to change the direction of migration if the accretion rate onto the planet is large. This implies that planets forming around metal-rich stars could have different compositions due to the larger influence of the heating torque compared to planets forming around metal poor stars. Thus, we think that the investigation of the heating torque in high metallicity environments could be crucial if planetary compositions can be linked to the planet formation pathways.

7 Conclusion

Here, we present a brief summary of our results. When comparing the type I torque formulae from Paardekooper et al. (2011) and Jiménez & Masset (2017), we found that they produce very similar planets in a single planet formation scenario (Fig. 2), despite significant differences in the migration maps (Fig. 4).

Adding the thermal torque following Masset (2017) introduces new regions of outwards migration to the migration maps due to the heating component (Fig. 4). The addition of the thermal torque to either of the two type I torque formulae investigated here does not yield a significant impact on single planet formation (Figs. 1 and 6), however, as outwards migration caused by the heating torque is limited to very small planets and a small interval of distances to the star around 3 AU (Fig. 4).

The specific influence on the planetary migration rates with the thermal torque crucially depends on the underlying disk model and on the accretion rate, where we followed an accretion recipe based on pebble accretion (e.g., Lambrechts & Johansen 2014; Bitsch et al. 2015a). If the accretion rates onto the planet are larger, outwards migration can be achieved, even thoughit is difficult to achieve these accretion rates with realistic dust-to-gas ratios (Fig. 5). However, even then the influence on planet formation seems minimal in our simulations (Fig. 6). Summing up, the main conclusion of this paper is that thermal torque only has a minimal impact on planetary growth.

The core growth model that we used here is pebble accretion. Other models will result in different accretion rates, which changes not only the trajectory in the migration maps and, hence, planet evolution in general, but also the heating torque.

While the difference between the type I torque formulae is small, it is still advisable to employ the Jiménez & Masset (2017) type Itorque formula as it is derived in a three-dimensional disk, which is a more realistic approach than the two-dimensional simulations from Paardekooper et al. (2011). Hence, the Jiménez & Masset (2017) formula should produce more realistic migration rates than the Paardekooper et al. (2011) formula.

The heating torque starts affecting planet formation at high pebble fluxes so it should be included in simulations of planets forming around high metallicity stars. Planets formed around low-metallicity stars, on the other hand, will barely exhibit any impact fromthe heating torque.

Acknowledgements

T.B. and B.B. thank the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support and the referee for their helpful comments that improved the manuscript.

Appendix A Type I migration model from Paardekooper et al. (2011)

The torque formulae of Paardekooper et al. (2011) are the result of an analysis of 2D hydrodynamical simulations in non-barotropic disks with finite viscosity and constant thermal diffusion. Furthermore, they use simple power laws in radius for surface density Σ, temperature T and viscosity ν. They set the background temperature profile as an equilibrium of viscous heating and radiative cooling. They found that the formulae they derived agreed with their simulations up to a maximum error of 20%.

All torque components are normalized to Γ₀ (Eq. (6)). We note that Γ₀ depends quadratically on the planet mass via q, so small planets migrate more slowly than more massive ones. As a result, the main part of the migration happens just before the planet becomes massive enough to open a gap which slows the migration down again as the planet enters type II migration. They also used an effective adiabatic index γ_eff which modifies the adiabatic index according to the planets position, the disk aspect ratio and thermal diffusivity according to $${\gamma }_{\text{eff}}=\frac{2Q\gamma }{\gamma Q+\frac{1}{2}\sqrt{2\sqrt{{\left({\gamma }^2{Q}^2+1\right)}^2-16Q^2(\gamma -1)}+2{\gamma }^2{Q}^2-2}},$$(A.1)

where, $$Q=\frac{2{\chi }_{\text{p}}}{3h^3{r}_{\text{p}}^2{\Omega }_{\text{p}}},$$(A.2)

where χ is the thermal diffusivity. This effective adiabatic index aims to account for regions in the disk that can cool efficiently, that is, the outer disk, producing isothermal results, while regions that can not cool efficiently in the inner disk produce adiabatic results. The width of the horseshoe region is determined as: $$x_{\text{s}}=\frac{1.11}{{\gamma }_{\text{eff}}^{1/4}}\sqrt{\frac{q}{h}}\cdot r_{\text{P}}.$$(A.3)

In this model, the saturation of the corotation torque is tackled via a modification of Masset (2001). When the viscosity gets very large, the viscous diffusion timescale is no longer much shorter than it takes for material that has performed a horseshoe turn to enter back into the horseshoe region. In this case, the horseshoe drag gets replaced by a linear corotation torque, which is typically smaller. Since this scenario requires high viscosity, saturation no longer occurs and the torque becomes fully linear. The saturation is modeled via a function F in Eq. (A.4), which can be seen in Figs. 2 and 3 in Paardekooper et al. (2011). Generally, as the viscosity increases, the horseshoe drag becomes smaller and the linear corotation torque gets larger. These two processes take place independently on different timescales and they are modeled by two different functions G ≥ 0 and K ≤ 1 such that: $${\Gamma }_{\text{c}}=G(p)\cdot F(p)\cdot {\Gamma }_{\text{hs}}+(1-K(p))\cdot {\Gamma }_{\text{c},\text{lin}},$$(A.4)

where p is a saturation parameter that depends on the planets’ position, its horseshoe width, and the viscosity. G and K follow the same general formula, but with different time scale parameters such that towards low viscosities the linear torque decreases slower than the horseshoe drag increases.

The saturation function F for the horseshoe drag is: $$F(p)=\frac{8I_{4/3}(p)}{3p{I}_{1/3}(p)+\frac{9}{2}{p}^2{I}_{4/3}(p)},$$(A.5)

where I_4∕3 is a modified Bessel function. The parameter that is entered into this function is: $$\begin{array}{lll}\hfill & \hfill & p=2\sqrt{k{x}_{\text{s}}^3}/3\hfill \\ \hfill & \hfill & k=\frac{r_{\text{p}}^2{\Omega }_{\text{p}}}{2\pi {\nu }_{\text{p}}}.\hfill \end{array}$$

The entropy torque saturates not only with viscosity but also with thermal diffusion, in this case the saturation parameter is: $$\begin{array}{lll}\hfill & \hfill & p_{\chi }=2\sqrt{k_{\chi }{x}_{\text{s}}^3}/3\hfill \\ \hfill & \hfill & k_{\chi }=\frac{r_{\text{p}}^2{\Omega }_{\text{p}}}{2\pi {\chi }_{\text{p}}}.\hfill \end{array}$$

The blending functions G and K follow $$\begin{array}{lll}\hfill & \hfill & G(p)=\{\begin{array}{cc}\frac{16}{25}{\left(\frac{45\pi }{8}\right)}^{3/4}{p}^{3/2}& p<\sqrt{\frac{8}{45\pi }}\\ 1-\frac{9}{25}{\left(\frac{8}{45\pi }\right)}^{4/3}{p}^{-8/3}& p\ge \sqrt{\frac{8}{45\pi }}\end{array}\hfill \\ \hfill & \hfill & K(p)=\{\begin{array}{cc}\frac{16}{25}{\left(\frac{45\pi }{28}\right)}^{3/4}{p}^{3/2}& p<\sqrt{\frac{28}{45\pi }}\\ 1-\frac{9}{25}{\left(\frac{28}{45\pi }\right)}^{4/3}{p}^{-8/3}& p\ge \sqrt{\frac{28}{45\pi }}\end{array}\hfill \end{array}$$

Lindblad torque

The formula for the Lindblad torque found by Paardekooper et al. (2011) is $${\gamma }_{\text{eff}}\frac{{\Gamma }_{\text{L}}}{{\Gamma }_0}=-2.5-1.7\beta +0.1\alpha ,$$(A.12)

where γ_eff is the effective adiabatic index from Eq. (A.1), − β is the temperature exponent and − α is the surface density exponent. Usually, both temperature and surface density decrease with the distance to the star so α and β are typically positive. For the disk model employed here, which is the one from Bitsch et al. (2015b), that is not always the case, however. The Lindblad torque is typically the main negative contribution and drives fast inward migration, as we can see in Fig. 1.

Barotropic torque

The barotropic torque is the result of the horseshoe trajectories which the material in the corotation region follows. The planet pushes material forwards that comes from the outer disk and pushes material backwards that comes from the inner disk. As a result of the surface density gradient in the disk, the planet typically pushes more mass backwards and migrates outwards as a consequence. Naturally, this torque component depends on the surface density gradient. The unsaturated barotropic torque from the horseshoe drag follows $${\gamma }_{\text{eff}}\frac{{\Gamma }_{\text{hs},\text{baro}}}{{\Gamma }_0}=1.1(\frac{3}{2}-\alpha )$$(A.13)

and the linear barotropic corotation torque is given by $${\gamma }_{\text{eff}}\frac{{\Gamma }_{\text{c},\text{lin},\text{baro}}}{{\Gamma }_0}=0.7(\frac{3}{2}-\alpha ).$$(A.14)

This part of the torque does not depend on thermal diffusion. As the surface density generally drops with increasing radius, α is positive, and the barotropic torque gets positive when the gradient is shallow.

Entropy torque

In this torque model the entropy torque is the main component responsible for outwards migration. It relies on the horseshoe trajectories of the corotating material like the barotropic torque, but is a consequence of the temperature difference between the inner and outer disk. Material that gets pushed from the outer disk to the inner disk typically gets warmer and contracts, leading to a greater gravitational pull in front of the planet than the material behind the planet, which cools down and expands. The entropy related horseshoe drag is: $${\gamma }_{\text{eff}}\frac{{\Gamma }_{\text{hs},\text{ent}}}{{\Gamma }_0}=7.9\frac{\xi }{\gamma }$$(A.15)

and the linear entropy corotation torque is $${\gamma }_{\text{eff}}\frac{{\Gamma }_{\text{c},\text{lin},\text{ent}}}{{\Gamma }_0}=(2.2-\frac{1.4}{\gamma })\xi ,$$(A.16)

with ξ = β − (γ − 1)α as the entropy gradient. Here, we note that γ is used, and not γ_eff, because the temperature gradient of the disk is unrelated to the cooling properties of the planet. The entropy torque needs both thermal andviscous diffusion to remain unsaturated. Thermal diffusion is required to restore the entropy gradient and viscous diffusion is required to diffuse the vortensity since the entropy torque originates from a streamline near the seperatrix of the horseshoe region. Consequently, Eq. (A.4) needs to be modified for the entropy torque as F(p) → F(p_ν)F(p_χ), with viscosity, ν, and thermal diffusivity, χ.Values for p_ν and p_χ depend on the planet’s position, horseshoe width, and on the viscosity or thermal diffusivity, respectively. Similarly, the transition functions G and K in Eq. (A.4) only get modified by choosing new arguments: $G(p)\to \sqrt{G(p_{\nu })G(p_{\chi })}$ and $1-K(p)\to \sqrt{(1-K(p_{\nu }))(1-K(p_{\chi }))}$. Since a higher thermal diffusivity keeps advection of entropy at the separatrix at bay, it increases the torque initially. However, when χ becomes toolarge, the disk essentially behaves locally isothermal where the corotation torque is given by the barotropic horseshoe drag and the linear part of the entropy torque. As the non-linear part of the entropy component is larger than the linear one, this leads to a smaller overall torque.

Appendix B Type I migration model from Jiménez & Masset (2017)

Jiménez & Masset (2017) claim that the formulae from Paardekooper et al. (2011) are not accurate enough for planet population synthesis. They argue that although the key mechanisms have been captured, many constants in the formulae are off. The main update compared to Paardekooper et al. (2011) is that these considerations are based on three-dimensional hydrodynamical simulations rather than just two-dimensional ones. They also update the horseshoe width to remain valid for more intermediate mass planets of several Earth-masses. The new formula for the horseshoe width x_s is: $$x_{\text{s}}=\frac{1.05{(q/h)}^{\frac{1}{2}}+3.4q^{\frac{7}{3}}/h^6}{1+2q^2/h^6}{r}_{\text{P}}.$$(B.1)

They do, however, still find some discrepancies between their formulae and simulations. They base their reasoning on multiple simplifications, including the neglect of the influence of the planet’s torque on the disk. Before the planet opens a full gap, it starts perturbing the density profile around it, which can alter the effect of the Lindblad torque and can change the thermal diffusivity, but this has not been taken into consideration here.

Rather than simplifying the formulae as far as possible, Jiménez & Masset (2017) kept one torque component for each source of torque. This results in four corotation components (Γ_C = Γ_V + Γ_S + Γ_T + Γ_VCT) and one Lindblad component Γ_L, all displayedin Fig. 1. Comparing the thick dashed purple lines in Fig. 1, which mark the total torque acting on the planet (without thermal torque), we can see that both follow similar curves, but are different in terms of the masses at which they act and in their magnitude. This can have a significant influence on migration.

Most torque components are still normalized to Γ₀ from Eq. (6). Also, as before, the corotation torque components are divided into two subcomponents, the linear one and the horseshoe drag which can saturate. Rather than developing a new saturation mechanism, they use the one from Masset & Casoli (2010). A notable difference to Paardekooper et al. (2011) is that the blending of the subcomponents happens at the same pace, whereas Paardekooper et al. (2011) assume the linear part to decrease slower than the horseshoe part increases as the viscosity decreases.

Lindblad torque

The Lindblad torque was taken from Tanaka et al. (2002) and reexamined in a locally isothermal setup. As it was only the coefficient for the temperature dependence that was reworked, the idea was to test the torque on a planet that was in the disk long enoughfor the corotation torque to saturate, leaving only the Lindblad torque and to fit the result. The new formula is: $$\frac{{\Gamma }_{\text{L}}}{{\Gamma }_0}=-(2.34-0.1\alpha +1.5\beta )\cdot f\left(\frac{\chi }{{\chi }_{\text{c}}}\right),$$(B.2)

with $$f(x)=\frac{\sqrt{\frac{x}{2}}+\frac{1}{\gamma }}{\sqrt{\frac{x}{2}}+1}.$$(B.3)

The factors 2.34 and −0.1 in Eq. (B.2) were taken from Tanaka et al. (2002), ${\chi }_{\text{c}}=r_{\text{P}}^2{h}^2{\Omega }_{\text{P}}$ is the critical thermal diffusivity from Masset & Casoli (2010). The function $f(\frac{\chi }{{\chi }_{\text{c}}})$ in Eq. (B.3) acts like the effective adiabatic index γ_eff in Eq. (A.12).

Entropy torque

Like in Paardekooper et al. (2011), the entropy torque is the result of an entropy gradient in the disk. However, this part has been updated in a meaningful way. Similar to Paardekooper et al. (2011), the entropy torque is mixed between the linear and horseshoe part by both the viscosity and the thermal diffusivity so there is an ɛ_ν and an ɛ_χ to determine the blending. The formula is: $${\Gamma }_{\text{S}}={\varepsilon }_{\nu }{\varepsilon }_{\chi }{\Gamma }_{\text{S}}^{\text{hs}}+(1-{\varepsilon }_{\nu }{\varepsilon }_{\chi }){\Gamma }_{\text{S}}^{\text{lin}}.$$(B.4)

The blending coefficients are taken from Masset & Casoli (2010) and amount to: $$\begin{array}{lll}\hfill & \hfill & {\varepsilon }_{\nu }={\left(1+{\left(6h{z}_{\nu }\right)}^2\right)}^{-1},\hfill \\ \hfill & \hfill & {\varepsilon }_{\chi }={(1+15h{z}_{\chi })}^{-1},\hfill \\ \hfill & \hfill & z_{\chi }=\frac{r_{\text{P}}\chi }{{\Omega }_{\text{P}}{x}_{\text{s}}^3},\hfill \\ \hfill & \hfill & z_{\nu }=\frac{r_{\text{P}}\nu }{{\Omega }_{\text{P}}{x}_{\text{s}}^3}.\hfill \end{array}$$

The linear and unsaturated horseshoe parts follow: $$\begin{array}{lll}\hfill & \hfill & \gamma \frac{{\Gamma }_{\text{S}}^{\text{lin}}}{{\Gamma }_0}=0.8{\xi }^{\prime },\hfill \\ \hfill & \hfill & {\Gamma }_{\text{S}}^{\text{uhs}}=3.3{\xi }^{\prime }{\Sigma }_0{\Omega }_{\text{P}}{x}_{\text{s}}^4,\hfill \\ \hfill & \hfill & {\xi }^{\prime }=\beta -0.4\alpha -\mathrm{0.64.}\hfill \end{array}$$

We note that ${\Gamma }_{\text{S}}^{\text{uhs}}$, as in the following unsaturated horseshoe torques, is not given in terms of Γ₀. The quantity ξ′ is called entropy gradient by Jiménez & Masset (2017) because it shows resemblance to the actual radial entropy gradient. The horseshoe drag is the product of the unsaturated formula and a saturation function F_S. $$\begin{array}{lll}\hfill & \hfill & {\Gamma }_{\text{S}}^{\text{hs}}=F_{\text{S}}{\Gamma }_{\text{S}}^{\text{uhs}},\hfill \\ \hfill & \hfill & F_{\text{S}}=1.2\cdot \overline{1.4z_{\chi }^{1/2}}\cdot \overline{1.8z_{\nu }^{1/2}},\hfill \end{array}$$

where $\overline{x}=\text{min}(1,x)$.

Temperature torque

This torque component is a result of the temperature gradient in the midplane of the disk and was incorporated into a type I torque recipe for the first time in Jiménez & Masset (2017). It is found in locally isothermal simulations with vanishing entropy and surface density gradients but β≠0. The temperature gradient produces vortensity in isothermal calculations which is concentrated at the downstream seperatrices. It is reasonable to assume that the thus generated torque depends on viscosity in the same way as the entropy torque. The formula is: $${\Gamma }_{\text{T}}={\varepsilon }_{\nu }{\Gamma }_{\text{T}}^{\text{hs}}+(1-{\varepsilon }_{\nu }){\Gamma }_{\text{T}}^{\text{lin}},$$(B.14)

with the mixing coefficient ɛ_ν from Eq. (B.5). Once again, the temperature component uses a different saturation function and linear component: $$\begin{array}{lll}\hfill & \hfill & \gamma \frac{{\Gamma }_{\text{T}}^{\text{lin}}}{{\Gamma }_0}=1.0\beta ,\hfill \\ \hfill & \hfill & {\Gamma }_{\text{T}}^{\text{hs}}=F_{\text{T}}(z_{\nu }){\Gamma }_{\text{T}}^{\text{uhs}},\hfill \\ \hfill & \hfill & {\Gamma }_{\text{T}}^{\text{uhs}}=0.73\beta {\Sigma }_0{\Omega }_{\text{P}}{x}_s^4,\hfill \\ \hfill & \hfill & F_{\text{T}}=1.2\cdot \overline{1.8z_{\nu }^{1/2}}.\hfill \end{array}$$

Vortensity torque

The vortensity torque was not updated by Jiménez & Masset (2017) since the horseshoe drag has the same expression in two and three dimensions in a barotropic disk. This component corresponds to the barotropic torque in Paardekooper et al. (2011). To simplify, they assume a vortensity gradient of $\frac{3}{2}-\alpha$, which is only valid for disks in which the surface density and angular frequency are power laws with respect to the distance to the star. The specific formula is: $${\Gamma }_{\text{V}}={\varepsilon }_{\text{b}}{\Gamma }_{\text{V}}^{\text{hs}}+(1-{\varepsilon }_{\text{b}}){\Gamma }_{\text{V}}^{\text{lin}}.$$(B.19)

With the blending coefficient: $${\varepsilon }_{\text{b}}={(1+30h{z}_{\nu })}^{-1}$$(B.20)

from Masset & Casoli (2010). The linear component is: $$\gamma \frac{{\Gamma }_{\text{V}}^{\text{lin}}}{{\Gamma }_0}=0.976-0.640\alpha .$$(B.21)

The horseshoe drag component is given by: $${\Gamma }_{\text{V}}^{\text{hs}}=\frac{8\pi }{3}{z}_{\nu }F(z_{\nu }){\Gamma }_{\text{V}}^{\text{uhs}},$$(B.22)

where ${\Gamma }_{\text{V}}^{\text{uhs}}$ is the unsaturated horseshoe drag and F(z_ν) in combination with the prefactors is the saturation function. The remaining quantities are: $${\Gamma }_{\text{V}}^{\text{uhs}}=\frac{3}{4}(\frac{3}{2}-\alpha ){\Sigma }_0{\Omega }_{\text{P}}^2{x}_s^4,$$(B.23) $$F(x)=\{\begin{array}{ll}1-x^{\text{n}rac12},\hfill & \text{if }x<\frac{4}{9}\hfill \\ \frac{4}{27x},\hfill & \text{otherwise}\text{.}\hfill \end{array}$$

Viscous coupling torque (VCT)

This part results from viscous production of vortensity at the density jumps at the seperatrices of the horseshoe region. As this phenomenon largely depends on the entropy gradient it is assumed to scale with the quantity ξ′ from Eq. (B.11), which takes the role of an entropy gradient. However, they normalize ξ′ by β = 1 for this term. Unlike the other corotation torques, this component does not have a linear equivalent. Masset & Casoli (2010) have proposed that Γ_VCT decays like the vortensity torque with the mixing coefficient ɛ_b. The formula is: $${\Gamma }_{\text{VCT}}=\frac{4\pi \xi }{\gamma }{\Sigma }_0{\Omega }_{\text{P}}{x}_{\text{s}}^4{\varepsilon }_{\text{b}}{z}_{\nu }\frac{z_{\nu }F(z_{\nu })-z_{\chi }F(z_{\chi })}{z_{\nu }-z_{\chi }}.$$(B.24)

In our the simulations, the viscous coupling term is found to have an almost negligible influence on the total torque.

Appendix C Thermal torque from Masset (2017)

Here, the thermal torque consists of two components: a cooling torque which is negative and a heating torque which is positive. Both act as a result of density perturbations in the gas profile around the planet due to heat exchange with the disk. In this paper, the planets luminosity is its accretion luminosity caused by the pebble flux onto the planet. The typical course of the thermal torque is visible in Fig. 1 (thick blue line).

One source of the cut-off effect on the thermal torque, the thermal critical mass, is plotted in Fig. C.1b at different times in the disk from Bitsch et al. (2015b). The evolution over time is dominated by the thermal diffusivity, which changes by several orders of magnitude (Fig. C.1a). The disk aspect ratio also changes over time, but the influence is smaller. It should be noted that the thermal critical mass (Eq. (7)) is only an estimate in Masset (2017) and needs to be investigated in more detail in the future, especially as it becomes very small at later times of disk evolution, which has important consequences for the thermal torque.

Fig. C.1 Disk thermal diffusivity and thermal critical mass over distance to the star at various times with a pebble flux of $2\times {10}^{-4}\cdot \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{yr}}$. The change in thermal diffusivity dominates the change in thermal critical mass. The other changing variable in the thermal critical mass is the aspect ratio.

An important ingredient for both components is a finite thermal diffusivity and the simulations can not be adiabatic as well. Masset (2017) derive their formulae from investigating low mass planets in a disk with thermal diffusivity with linear perturbation theory in a three-dimensional shearing sheet. Their results for the cooling torque are: $$\begin{array}{lll}\hfill & \hfill & \gamma \frac{{\Gamma }_{\text{thermal}}^{\text{cold}}}{{\Gamma }_0}=-1.61(\gamma -1)\frac{x_{\text{p}}}{{\lambda }_{\text{c}}}\hfill \\ \hfill & \hfill & {\lambda }_{\text{c}}=\sqrt{\frac{\chi }{q^{\prime }{\Omega }_{\text{P}}\gamma }},\hfill \end{array}$$

where x_p is the distance from a planet to its corotation, λ_c is the length scale over which a luminous planet releases heat, and q′ = 1.5 is the shear parameter for a Keplerian disk. The heating torque on the other hand is: $$\begin{array}{lll}\hfill & \hfill & \gamma \frac{{\Gamma }_{\text{thermal}}^{\text{heating}}}{{\Gamma }_0}=1.61(\gamma -1)\frac{x_{\text{p}}}{{\lambda }_{\text{c}}}\frac{L}{L_{\text{c}}}\hfill \\ \hfill & \hfill & L_{\text{c}}=\frac{4\pi GM\chi {\rho }_0}{\gamma },\hfill \end{array}$$

where L is the luminosity of the planet, ρ₀ is the unperturbed disk density and L_c is used to normalize the density perturbation caused by the heat release. See Fig. C.2 for a plot of the critical luminosity at various times. The luminosity here is the accretion luminosity of a planetary embryo, which, according to Chrenko et al. (2017) follows: $$L=\frac{GM}{R}\dot{M},$$(C.5)

where R is the radius of the planet and Ṁ is the solid accretion rate onto the planet. The radius of the planet was determined from its mass by assuming the planet is spherical with a constant density of $5\frac{\text{g}}{{\text{cm}}^3}$, which corresponds to an Earth-like density. Considering the density of ice giants, which is roughly $2\frac{\text{g}}{{\text{cm}}^3}$, we would get a smaller luminosity, meaning the heating torque would be reduced. The sum of heating and cooling torque can be simplified tothe final form: $$\begin{array}{lll}\hfill & \hfill & {\Gamma }_{\text{thermal}}^{\text{total}}=1.61\frac{\gamma -1}{\gamma }\eta \left(\frac{H}{{\lambda }_{\text{c}}}\right)\left(\frac{L}{L_{\text{c}}}-1\right){\Gamma }_0,\hfill \\ \hfill & \hfill & \eta =\frac{\alpha }{3}+\frac{\beta +3}{6},\hfill \end{array}$$

where H is the height of the disk. As a consequence, when the planet’s luminosity is smaller than the critical value from Eq. (C.4), the cooling torque dominates over the heating torque and the overall thermal torque is negative.

Fig. C.2 Critical luminosity (Eq. (C.4)) for a planet mass of 1 ME (right) at various times. The kinks in the critical luminosity are caused by the thermal diffusivity (see Fig. C.1a).

Appendix D Calculation of type II migration timescales

Type II migration describes the movement of massive planets in the disk. Massive here means that the planet opens a gap in the disk. A gap corresponds to an area in which the surface density is smaller than 10% of the unperturbed value. The opening of a gap by the planet or synonymously the transition into type II migration is modeled following Crida et al. (2006). In particular, a gap parameter g is calculated as: $$g=\frac{3H}{4r_{\text{H}}}+\frac{50\nu H^2}{{(q{r}_{\text{P}})}^2},$$(D.1)

where $r_{\text{H}}=r{\left[M_{\text{c}}/\left(3M_{\star }\right)\right]}^{1/3}$ is the planet’s Hill radius, that is, the radius in which the planet dominates accretion. When the gap parameter is smaller than 3, the migration timescale gets calculated as a linear interpolation between type I and type II migration, similar as in Dittkrist et al. (2014). This corresponds to the hatched regions in the migration maps. When this gap parameter is smaller than 1, the planet has opened a full gap and migrates completely in the type II regime, which is marked by cross hatched regions in the migration maps. Larger values of the gap parameter indicate that the planet does not perturb the disk significantly and migrates completely in the type I fashion. At later times, it is easier for a planet to open a gap because the disks aspect ratio decreases in time.

A recent study by Kanagawa et al. (2018) suggests that the type II migration timescale gets calculated just like the type I migration timescale, but gets damped by the lower surface density in the gap. The usual model, which is incorporated in this paper, however, suggests that the type II migration timescale is equal to the accretion time scale τ_visc of the disk (e.g., Dürmann & Kley 2015; Robert et al. 2018), with $${\tau }_{\text{visc}}=r_{\text{P}}^2/\nu .$$(D.2)

Depending on the viscosity in the disk, the planetary migration rate changes. In low-viscosity environments planets in type-II migration move slower and can thus accrete more material and grow bigger. Baruteau et al. (2014) have studied how the planet can slow down type II migration even further when growing much more massive than the gas surrounding the gap, which slows down the viscous accretion. When $M_{\text{P}}>4\pi {\Sigma }_{\text{g}}{r}_{\text{P}}^2$, the migration gets slowed down, resulting in a formula for the type II migration timescale: $${\tau }_{\text{II}}={\tau }_{\text{visc}}\cdot \max \left(1,\frac{M_{\text{P}}}{4\pi {\Sigma }_{\text{g}}{r}_{\text{P}}^2}\right).$$(D.3)

Appendix E Additional plots

Fig. E.1 Minimum core accretion rate that is required for the thermal torque to become positive in ME /year for different times. It generally drops in time and is high in the inner and outer disk. At a few AU there is a sweet spot where the thermal torque can become positive more easily, caused by the critical luminosity (see Fig. C.2). For a more comprehensive review of these plots see Fig. 5.

Fig. E.2 Minimum pebble flux that is required for the thermal torque to become positive in ME /year for different times. It generally drops in time and is high in the inner and outer disk. At a few AU there is a sweet spot, caused by the critical thermal diffusivity. As the thermal critical mass is only an estimate as of yet, values above are still plotted here. For a more comprehensive review of these plots see Fig. 5.

Fig. E.3 Migration maps in the Ida et al. (2016) disk model after 1 Myr at α = 5.4× 10−3. Like in Fig. 4, hatched regions mark the transition to type II migration, while planets in the cross-hatched areas are migrating completely in the type II regime and regions of outwards migration that are caused by the thermal torque are outlined by dashed blue lines, while regions of outwards migration inside solid blue lines are due to the original torque formulae. Compared to the corresponding plots in the Bitsch et al. (2015b) disk (Fig. 4c), less outwards migration takes place. This is particularly true for the Paardekooper et al. (2011) formula and outwards migration due to the heating torque.

Fig. E.4 Planet evolution maps in the Ida et al. (2016) disk model at α = 5.4 × 10−3. As in Fig. 2, planets that are below the dashed blue line have reached pebble isolation mass and accrete gas. Planets that are interior to the dashed white line have reached runaway gas accretion and can grow very massive. The pebble flux is ${\dot{M}}_{\text{peb}}=2\times {10}^{-4}\cdot \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$ in all plots. As in the Bitsch et al. (2015b) disk, the formulae produce really similar results, however planets grow much lighter with the same pebble fluxes in the Ida et al. (2016) model due to the region of outwards migration from corotation torques being smaller.

References

All Figures

Fig. 1 Individual torque components of the Paardekooper et al. (2011) and Jiménez & Masset (2017) formulae for planets at 3 AU and 1 Myr, depending on planetary mass. The Jiménez & Masset (2017) formula adds new contributions, but the total torque looks similar. This is just one arbitrarily chosen set of disk evolution time and position in the disk, which manifests many general trends. Outwards migration occurs when the total torque is positive. This plot shows outwards migration caused by the heating torque (see Sect. 4) for small planets as well as outwards migration caused by the corotation torques for slightly heavier planets. The pebble flux through the disk is ${\dot{M}}_{\text{peb}}=2\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$, which is relevant for the thermal torque as discussed in Sect. 4. The viscosity is set to α = 5.4 × 10−3 during these simulations.

In the text

Fig. 2 Planet formation maps with different pebble fluxes in the disk. The axes display where and when a planetary embryo has been inserted in the disk. The color codes how massive the planets grow, while the solid lines outline its final orbitalradius. Planets that are below the dashed blue line have reached pebble isolation mass and accrete gas. Planets that are interior to the dashed white line have reached runaway gas accretion and can grow very massive. Ṁpeb marks the pebble flux through the disk, which decreases exponentially over time. Differences between the formulae in terms of impact on planet formation are small. As the pebble density in the outer disk is small, planets that are inserted there do not grow much.

In the text

Fig. 3 Planet formation maps with different viscosities for planetary migration using the Jiménez & Masset (2017) formula. Like in Fig. 2, planets that are below the dashed blue line havereached pebble isolation mass and accrete gas. Planets that are interior to the dashed white line have reached runaway gas accretion and can grow very massive. The pebble flux is ${\dot{M}}_{\text{peb}}=2\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$ in all plots. A lower choice of α mostly leads to more massive planets due to the quicker transition to type II migration. In addition, the inwards migration of planets is reduced due to the slow type II migration.

In the text

Fig. 4 Migration maps after 1 Myr with different α-viscosity values for the calculation of migration timescales and the gap parameter. The left plots show the Paardekooper et al. (2011) formula and the right ones the Jiménez & Masset (2017) formula. Hatched regions mark the transition to type II migration, while planets in the cross-hatched areas are migrating completely in the type II regime. Regions of outwards migration that are caused by the thermal torque, which is calculated with an accretion rate of ${\dot{M}}_{\text{peb}}\approx 1.4\times {10}^{-4}\frac{M_{\text{E}}}{\text{year}}$, are outlined by dashed blue lines, while regions of outwards migration inside solid blue lines are due to the original torque formulae. Outwards migration in the original formulae is caused by the torques originating from the corotation region (Paardekooper et al. 2011) or by the temperature torque (Jiménez & Masset 2017). As type II migration depends on the disk properties rather the type I torque formula, the cross-hatched regions are identical between the corresponding left and right plots.

In the text

Fig. 5 Minimum core accretion rate and pebble flux that are required for a positive thermal torque in a disk that is 1 Myr old. The flux of pebbles needed to make a positive thermal torque also depends on the choice of Stokes number, which is set to 0.1 here. For a higher Stokes number, lower pebble fluxes are required for a positive thermal torque. The thermal critical mass used here is only an estimate, so values above might be of interest in the future. Below the thermal critical mass and in the outer disk, unreasonably high pebble fluxes are required for a positive thermal torque. Due to the exponential decrease in pebble flux that is employed in these simulations, about 72% of the initial pebble flux are left after 1 Myr, which corresponds to $1.4\times {10}^{-4}\frac{M_{\text{E}}}{\text{year}}$ in the standard case. Keep in mind that a positive thermal torque does not automatically correspond to outwards migration. For different parameters of the disk age and how it influences the minimal pebble flux needed for a positive thermal torque, see Figs. E.1 and E.2.

In the text

Fig. 6 Planet formation maps with and without the thermal torque from Masset (2017) using the Jiménez & Masset (2017) formula for the other torque components. Simulations were performed with an increased pebble flux of ${\dot{M}}_{\text{peb}}=5\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$, i.e. about 950 ME over 3 Myr. The plot should be read just like Fig. 2. Planets in the inner disk are significantly affected by the thermal torque due to the larger amount of pebbles. Planets that are inserted in the outer disk are not noticeably influenced by the thermal torque.

In the text

Fig. 7 Growth tracks of planets that have been inserted at different radii after 1 Myr with the Jiménez & Masset (2017) formula (no thermal torque) in black and the Jiménez & Masset (2017) + Masset (2017) formula (includes thermal torque) in red. The lines become dashed when the planets reach pebble isolation mass and are accreting gas. An increased pebble flux of ${\dot{M}}_{\text{peb}}=5\times {10}^{-4}\times \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$ was used to achieve a larger contribution from the thermal torque. After performing significant outwards migration, planets thatgrow with the thermal torque migrate to the star quickly and grow less massive than the planets that were evolved without thermal torque, which reach runaway gas accretion. This might be due to the fact that heavier planets migrate faster, which gives them less time to contract an envelope before reaching the inner edge of the disk.

In the text

	Fig. C.1 Disk thermal diffusivity and thermal critical mass over distance to the star at various times with a pebble flux of $2\times {10}^{-4}\cdot \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{yr}}$. The change in thermal diffusivity dominates the change in thermal critical mass. The other changing variable in the thermal critical mass is the aspect ratio.
In the text

	Fig. C.2 Critical luminosity (Eq. (C.4)) for a planet mass of 1 ME (right) at various times. The kinks in the critical luminosity are caused by the thermal diffusivity (see Fig. C.1a).
In the text

	Fig. E.1 Minimum core accretion rate that is required for the thermal torque to become positive in ME /year for different times. It generally drops in time and is high in the inner and outer disk. At a few AU there is a sweet spot where the thermal torque can become positive more easily, caused by the critical luminosity (see Fig. C.2). For a more comprehensive review of these plots see Fig. 5.
In the text

Fig. E.2 Minimum pebble flux that is required for the thermal torque to become positive in ME /year for different times. It generally drops in time and is high in the inner and outer disk. At a few AU there is a sweet spot, caused by the critical thermal diffusivity. As the thermal critical mass is only an estimate as of yet, values above are still plotted here. For a more comprehensive review of these plots see Fig. 5.

In the text

Fig. E.3 Migration maps in the Ida et al. (2016) disk model after 1 Myr at α = 5.4× 10−3. Like in Fig. 4, hatched regions mark the transition to type II migration, while planets in the cross-hatched areas are migrating completely in the type II regime and regions of outwards migration that are caused by the thermal torque are outlined by dashed blue lines, while regions of outwards migration inside solid blue lines are due to the original torque formulae. Compared to the corresponding plots in the Bitsch et al. (2015b) disk (Fig. 4c), less outwards migration takes place. This is particularly true for the Paardekooper et al. (2011) formula and outwards migration due to the heating torque.

In the text

Fig. E.4 Planet evolution maps in the Ida et al. (2016) disk model at α = 5.4 × 10−3. As in Fig. 2, planets that are below the dashed blue line have reached pebble isolation mass and accrete gas. Planets that are interior to the dashed white line have reached runaway gas accretion and can grow very massive. The pebble flux is ${\dot{M}}_{\text{peb}}=2\times {10}^{-4}\cdot \exp (-\frac{t}{t_{\text{f}}})\frac{M_{\text{E}}}{\text{year}}$ in all plots. As in the Bitsch et al. (2015b) disk, the formulae produce really similar results, however planets grow much lighter with the same pebble fluxes in the Ida et al. (2016) model due to the region of outwards migration from corotation torques being smaller.

In the text