Thermal torque effects on the migration of growing low-mass planets

TL;DR

This paper investigates how thermal effects and planetary luminosity influence type I migration of low-mass planets in protoplanetary discs, revealing that luminosity can cause significant outward migration, impacting planet formation models.

Contribution

It compares updated non-isothermal migration rates with the effects of planetary luminosity, highlighting their combined impact on planetary migration directions.

Findings

Luminosity can reverse inward migration to outward migration.

Non-isothermal effects significantly alter migration rates.

Outward migration influences planet formation outcomes.

Abstract

As planets grow the exchange of angular momentum with the gaseous component of the protoplanetary disc produces a net torque resulting in a variation of the semi-major axis of the planet. For low-mass planets not able to open a gap in the gaseous disc this regime is known as type I migration. Pioneer works studied this mechanism in isothermal discs finding fast inward type I migration rates that were unable to reproduce the observed properties of extrasolar planets. In the last years, several improvements have been made in order to extend the study of type I migration rates to non-isothermal discs. Moreover, it was recently shown that if the planet's luminosity due to solid accretion is taken into account, inward migration could be slowed down and even reversed. In this work, we study the planet formation process incorporating, and comparing, updated type I migration rates for non-isothermal discs and the role of planet's luminosity over such rates. We find that the latter can have important effects on planetary evolution, producing a significant outward migration for the growing planets.

Equations

$$\frac{\partial P}{\partial z}$$

$$\frac{\partial F}{\partial z}$$

$$\frac{\partial T}{\partial z}$$

$$P_s$$

$$F_s$$

$$0$$

$$T_{\text{irr}}=T_{\star }{\left[\frac{2}{3\pi }{\left(\frac{R_{\star }}{R}\right)}^3+\frac{1}{2}{\left(\frac{R_{\star }}{R}\right)}^2\left(\frac{H}{R}\right)\left(\frac{d\log H}{d\log R}-1\right)\right]}^{0.5},$$

$$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}=\frac{3}{R}\frac{\partial }{\partial R}\left[R^{1/2}\frac{\partial }{\partial R}\left(\bar{\nu }{\Sigma }_{\text{g}}{R}^{1/2}\right)\right],$$

$${\Sigma }_{\text{g}}(R)={\Sigma }_{\text{g}}^0{\left(\frac{R}{R_{\text{c}}}\right)}^{-b}{e}^{-(R/R_{\text{c}}{)}^{2-b}},$$

$${\Sigma }_{\text{p}}(R)=z_0{\eta }_{\text{ice}}{\Sigma }_{\text{g}}(R).$$

$$\frac{\partial }{\partial t}\left[{\Sigma }_{\text{p}}(R,r_{p_j})\right]$$

$$\dfrac{{\rm d}\text{M}_{\text{C}}}{{\rm d}t}\bigg{|}_{\text{psimal}}=\Sigma_{\text{p}}\,(R_{\text{P}})\,\text{R}_{\text{H}}^{2}\,P_{\text{coll}}\,\Omega_{\text{P}},$$

$$\displaystyle\frac{{\rm d}\text{M}_{\text{C}}}{{\rm d}t}\bigg{|}_{\text{pebbles}}=\begin{cases}2\beta\,\Sigma_{\text{p}}(R_{\text{P}})\,\text{R}_{\text{H}}^{2}\,\Omega_{\text{P}},~{}\text{if }~{}0.1\leq\text{$S_{t}$}<1,\\ \\ 2\beta\,\left(\frac{\text{$S_{t}$}}{0.1}\right)^{2/3}\,\Sigma_{\text{p}}(R_{\text{P}})\,\text{R}_{\text{H}}^{2}\,\Omega_{\text{P}},~{}\text{if }~{}\text{$S_{t}$}<0.1,\end{cases}$$

$${\text{M}}_{\text{crit}}=10{\left(\frac{{\dot{\text{M}}}_{\text{C}}}{10^{-6}{\text{M}}_{\oplus }{\text{yr}}^{-1}}\right)}^{0.25}{\text{M}}_{\oplus },$$

$$\frac{\mathrm{d}{\text{M}}_{\text{g}}}{\mathrm{d}t}=\text{min}[{\dot{\text{M}}}_{\text{KH}},{\dot{\text{M}}}_{\text{disc}},{\dot{\text{M}}}_{\text{GP}}],$$

$${\dot{\text{M}}}_{\text{disc}}=3\pi \nu {\Sigma }_{\text{g,P}},$$

$${\dot{\text{M}}}_{\text{GP}}=\dot{A}{\Sigma }_{\text{acc}},$$

$$\dot{A}=0.29{\left(\frac{H_{\text{g,P}}}{R_{\text{P}}}\right)}^{-2}{\left(\frac{{\text{M}}_{\text{P}}}{{\text{M}}_{\star }}\right)}^{4/3}{R}_{\text{P}}^2{\Omega }_{\text{P}}$$

$$\Sigma(x)=\begin{cases}\Sigma_{\text{g,P}}\text{exp}\left[-\left(\dfrac{x}{l}\right)^{-3}\right]&\text{ if $x>x_{\text{m}}$},\\ \\ \!\begin{aligned} &\Sigma_{\text{g,P}}\text{exp}\Bigg{[}-\dfrac{1}{2}\left(\dfrac{x}{H_{\text{g,P}}}-\dfrac{5}{4}\dfrac{x_{\text{m}}}{H_{\text{g,P}}}\right)^{2}\\ \\ &+\dfrac{1}{32}\left(\dfrac{x_{\text{m}}}{H_{\text{g,P}}}\right)^{2}-\left(\dfrac{x_{\text{m}}}{l}\right)^{-3}\Bigg{]}\\ \end{aligned}&\text{if $x\leq x_{\text{m}}$},\end{cases}$$

$$l=0.146{\left(\frac{\nu }{10^{-5}{R}_{\text{P}}^2{\Omega }_{\text{P}}}\right)}^{-1/3}{\left(\frac{{\text{M}}_{\text{P}}}{10^{-3}{\text{M}}_{\star }}\right)}^{2/3}{R}_{\text{P}},$$

$$x_{\text{m}}=0.207{\left(\frac{H_{\text{g,P}}}{0.1R_{\text{P}}}\frac{{\text{M}}_{\text{P}}}{10^{-3}{\text{M}}_{\star }}\right)}^{2/5}{\left(\frac{\nu }{10^{-5}{R}_{\text{P}}^2{\Omega }_{\text{P}}}\right)}^{-1/5}{R}_{\text{P}},$$

$$\dfrac{{\rm d}R_{\text{P}}}{{\rm d}t}\bigg{|}_{\text{migI}}=-2R_{\text{P}}\dfrac{\Gamma}{\mathcal{L}_{\text{P}}},$$

$$\Gamma ={\Gamma }_{\text{L}}+{\Gamma }_{\text{C}}.$$

$${\Gamma }_{\text{L}}=(-2.5-1.7{\beta }^{\prime }+0.1{\alpha }^{\prime })\frac{{\Gamma }_0}{{\gamma }_{\text{eff}}}$$

$${\Gamma }_0={\left(\frac{q}{h}\right)}^2{\Sigma }_{\text{g,P}}{R}_{\text{P}}^4{\Omega }_{\text{P}}^2,$$

$${\gamma }_{\text{eff}}=$$

$$\frac{2\gamma Q}{\gamma Q+\frac{1}{2}\sqrt{2\sqrt{({\gamma }^2{Q}^2+1{)}^2-16Q^2(\gamma -1)}+2{\gamma }^2{Q}^2-2}},$$

$${\chi }_{\text{P}}=\frac{16(\gamma -1)\sigma T_{\text{P}}^4}{3{\kappa }_{\text{P}}{\rho }_{g,\text{P}}{P}_{\text{P}}},$$

$${\Gamma }_{\text{C}}={\Gamma }_{\text{C,bar}}+{\Gamma }_{\text{C,ent}}.$$

$$\Gamma ={\Gamma }_{\text{Type I}}+{\Gamma }_{\text{thermal}},$$

$${\Gamma }_{\text{L}}=(-2.34-1.5{\beta }^{\prime }+0.1{\alpha }^{\prime }){\Gamma }_{0}f({\chi }_{\text{P}}/{\chi }_{\text{C}}),$$

$$f(x)=\frac{\sqrt{x/2}+1/\gamma }{\sqrt{x/2}+1},$$

$${\Gamma }_{\text{C}}={\Gamma }_{\text{C,vor}}+{\Gamma }_{\text{C,ent}}+{\Gamma }_{\text{C,temp}}+{\Gamma }_{\text{C,vv}}.$$

Full text

††pubyear: 2019††pagerange: Thermal torque effects on the migration of growing low-mass planets–A

1 Introduction

In the classical core accretion mechanism (Safronov, 1969; Mizuno, 1980; Stevenson, 1982; Bodenheimer & Pollack, 1986; Greenzweig & Lissauer, 1992), the formation of both terrestrial and giant planets occurs through subsequent stages. Immediately after the collapse of the primordial gas cloud a circumstellar disc is formed. The dust of the protoplanetary disc quickly settles into the disc mid-plane, where it coagulates to form larger solid structures known as planetesimals. Then planetesimals coaggregate to form planetary embryos, which can then grow to form planets.

At the moment, there are two main scenarios for planetesimal formation: collisional growth models, wherein a gradual growth by continuous collisions is invoked; or gravoturbulent models, wherein the gravitational collapse of dust accumulation is invoked (see Blum, 2018, and references therein for an update review). The first scenario predicts that most of the dust particles coagulate to form planetesimals, which continue growing by mutual accretion until planetary embryos of about the mass of the Moon form. At this stage, these are the only bodies that continue growing by the accretion of planetesimals in a process often called oligarchic growth regime (Ida & Makino, 1993; Kokubo & Ida, 1998; Ormel et al., 2010). In the second scenario, dust particles are concentrated in high pressure regions in sufficiently high amount to be able to trigger the streaming instability. In this case, the planetesimals are directly formed by the gravitational collapse of the dust aggregates (Youdin & Goodman, 2005; Johansen et al., 2007).

Once planetesimals are formed they can efficiently accrete small particles and quickly form massive solid embryos (Ormel & Klahr, 2010; Johansen & Lacerda, 2010; Lambrechts & Johansen, 2012; Lambrechts et al., 2014; Chambers, 2016). There are two necessary conditions for this process: massive planetesimals combined to a high concentration of small particles.

Dust particles and small planetesimals experience a fast radial-drift. This is because the gaseous component of the disc rotates at a sub-Keplerian velocity, while dust and planetesimals orbit at Keplerian velocities. This difference causes a drag force on the solid bodies, which results in a loss of angular momentum. As bodies grow and become planetary embryos, the drag force becomes negligible. However, planetary embryos are able to locally modify the gas density structure. This effect leads to an exchange of angular momentum, which results in a net torque that makes the planets migrate within the disc. Different regimes of migration have been studied during the last decades, starting with the pioneer works of Lin & Papaloizou (1986) and Ward (1997). Planet–disc interactions are key in the orbital evolution of forming planets, leading to potentially large variations in their semi-major axis. For planets not massive enough to open a gap in the gaseous disc, the migration regime is usually known as type I migration. Early models of type I migration were developed for isothermal discs (Tanaka et al., 2002). Such models found very high inward type I migration rates, so planets migrate too fast towards the central star. Therefore, in order to reproduce observations, many authors had to reduce by up to a factor of 0.001 these previous migration rates (Ida & Lin, 2004; Alibert et al., 2005; Miguel et al., 2011).

Numerical studies including further physical processes — magnetic fields (Guilet et al., 2013) and non-isothermal protoplanetary discs for instance — have shown that type I migration can be slowed down, and even reversed. In hydrodynamical simulations of non-isothermal discs — in 2D and 3D — it was found that the outward type I migration strongly depends on the planetary mass and semi-major axis, but also on the disc structure and thermodynamics (Paardekooper & Mellema, 2006; Baruteau & Masset, 2008; Paardekooper & Papaloizou, 2008; Kley et al., 2009). Later, Masset & Casoli (2010) and Paardekooper et al. (2010, 2011) derived semi-analytical prescriptions for type I migration rates in non-isothermal discs, which can be incorporated in planet formation models. In fact, the migration rates derived by Paardekooper et al. (2010, 2011) are the most often employed in planet formation models that consider non-isothermal discs (Fortier et al., 2013; Alibert et al., 2013; Coleman & Nelson, 2014; Dittkrist et al., 2014; Bitsch et al., 2015a, b; Coleman & Nelson, 2016; Guilera et al., 2017; Izidoro et al., 2017; Kobayashi & Tanaka, 2018; Ndugu et al., 2018). It is worth noting that the type I migration rate prescriptions given by Masset & Casoli (2010) and Paardekooper et al. (2010, 2011) were derived from the results of 2D hydrodynamical simulations. Although they account for most of the mechanisms that contribute to the total torque over the planet, these prescriptions do not match the more physical 3D radiative hydrodynamic simulations with enough accuracy (Kley et al., 2009; Bitsch & Kley, 2011).

More recently, however, Jiménez & Masset (2017) updated the torque formula of Masset & Casoli (2010) by means of 3D numerical simulations. This resulted in an improvement in the derived type I migration prescriptions for both the Lindblad and the corotation torques with respect to the aforementioned works. These new type I migration recipes were compared with the results from the 3D radiative hydrodynamical simulations from Bitsch & Kley (2011) and Lega et al. (2015), showing good agreement.

Using 3D radiative hydrodynamical simulations, Lega et al. (2014) found that in the proximity to a low-mass planet the gas is cooler and denser than it would be if it behaved adiabatically. Such effect is asymmetric, generating an additional negative cold torque, which generates an inward migration. Additionally, Benítez-Llambay et al. (2015) found that the heat released by a planet due to accretion of solids is diffused in the nearby disc, generating two asymmetric hot and low-density lobes. This effect produces a heating torque that can become positive, depending on the amount of heat released by the planet. Then, Chrenko et al. (2017) also performed radiative hydrodynamical simulations considering growing planets immersed in the disc. They also found planet outward migration if the heat released into the disc by the planets due to the accretion of solids is taken into account. However, neither of these works derived analytical prescriptions that reproduced such phenomena.

It is important to note that the effects found by Lega et al. (2014) and Benítez-Llambay et al. (2015) were not considered by Jiménez & Masset (2017) for the derivation of their migration recipes. However, these two phenomena were studied analytically in a novel work by Masset (2017), who proposed a semi-analytical prescription for the thermal torque, defined by the author as the sum of the two aforementioned effects. This torque can be incorporated in planet formation models.

In this work, we update our planet formation model PLANETALP (Ronco et al., 2017; Guilera et al., 2017) by implementing the type I migration rates derived by Jiménez & Masset (2017), along with the thermal torque derived by Masset (2017). We aim to study the impact of these new prescriptions on planet formation. In doing so, we pay particular attention to the thermal torque. This work is organised in the following way: in Sect. 2, we describe our planet formation model and the improvements to the type I migration rates. Readers not interested in the technical details of the simulations can directly jump to Sect. 3, where we present the migration maps obtained with these new prescriptions and the resulting impact on planet formation. Finally, we discuss our findings and draw our conclusions in Sect. 4.

2 Description of PLANETALP

In a series of previous works (Guilera et al., 2010, 2014, 2017; Guilera & Sándor, 2017; Ronco et al., 2017), we developed a planet formation code called PLANETALP. It models the formation of a planetary system immersed in a protoplanetary disc that evolves in time. The protoplanetary disc is modelled by two components: a gaseous and a solid one. Planets immersed in the disc start to grow by the accretion of solid material. Once they are massive enough, they are also able to accrete the surrounding gas. The growing planets are allowed to migrate under type I and type II migration along the disc. Even though PLANETALP allow us to calculate the simultaneous formation of several planets, in this work we will only calculate the formation of only one planet per disc.

Below, we describe the most relevant features of PLANETALP, especially the new type I migrations rates incorporated in this work.

2.1 The gaseous disc

In the latest version of PLANETALP (Guilera et al., 2017), the disc gaseous component is modelled by a standard 1D+1D viscous accretion disc. We consider an axisymmetric, irradiated disc in hydrostatic equilibrium. To calculate the vertical structure of the disc, we solve for each radial bin the structure equations (Papaloizou & Terquem, 1999; Alibert et al., 2005; Migaszewski, 2015)

[TABLE]

Where $P$, $\rho$, $F$, $T$ and $z$ represent the pressure, density, radiative heat flux, temperature, and vertical coordinate of the disc, respectively. We note $\Omega$ the Keplerian frequency at a given radial distance and $\nu=\alpha c_{s}^{2}/\Omega$ the viscosity (Shakura & Sunyaev, 1973), with $\alpha$ a free parameter (we varied $\alpha$ between $10^{-4}$ and $10^{-2}$) and $c_{s}^{2}=P/\rho$ the square of the local isothermal sound speed. In our model, convection sets in when the standard Schwarzschild criterion is fulfilled. In convective regions, the vertical temperature gradient $\nabla=d\log T/d\log P$ is calculated considering the mixing length theory, following the prescription of Kippenhahn et al. (2012); while in purely radiative regions, the temperature gradient $\nabla_{\rm rad}$ is given by $\nabla_{\text{rad}}=3\kappa\rho F/16\sigma\Omega^{2}zT^{4}$, where $\sigma$ is the Stefan-Boltzmann constant and $\kappa$ is the local Rosseland mean opacity (Bell & Lin, 1994). We adopt the equation of state of an ideal diatomic gas $P=\rho kT/\mu m_{\text{H}}$, where $k$ is the Boltzmann constant, $\mu$ the mean molecular weight, $m_{\text{H}}$ the mass of the hydrogen atom, and $\nabla_{\text{ad}}$ the adiabatic gradient. Since the main component of the gas disk is molecular hydrogen, it is reasonable to assume as in previous works that $\mu=2$ and $\nabla_{\text{ad}}=2/7$.

To solve the system of Eqs. 1, we consider that the boundary conditions at the surface of the disc $H$ as in Papaloizou & Terquem (1999), Alibert et al. (2005) and Migaszewski (2015). Namely, $P_{s}=P(z=H)$, $F_{s}=F(z=H)$, and $T_{s}=T(z=H)$. These quantities are given by:

[TABLE]

where $\tau_{\text{ab}}=10^{-2}$ is the optical depth, $T_{\text{b}}=10$ K is the background temperature, and $\dot{M}_{\text{st}}$ is the equilibrium accretion rate. The temperature associated with the stellar irradiation is given by:

[TABLE]

where $R$ is the radial coordinate, and $d\log H/d\log R=9/7$ (Chiang & Goldreich, 1997). We note that here we are not considering the self-shadowing of the disc (Baillié et al., 2015, 2016; Bitsch et al., 2015a; Migaszewski, 2015). The central object is assumed to be a $1\,M_{\odot}$ protostar with $R_{\star}=2~{}R_{\odot}$ and $T_{\star}=4000$ K (Baraffe et al., 2015). We note that the stellar luminosity does not evolve in time in our model. In addition to the surface boundary conditions, in the mid-plane of the disc the heat flux should vanish: $F(z=0)=F_{0}=0$. The system of Eqs. 2 is solved by a multidimensional Newton-Raphson algorithm. The vertical structure of the disc is solved following Alibert et al. (2005) and Migaszewski (2015), using a shooting method with a Runge-Kutta-Fehlberg integrator.

Then, the time evolution of the gas surface density $\Sigma_{\text{g}}$ is calculated solving the classical diffusion equation of Pringle (1981):

[TABLE]

where $\bar{\nu}$ is the mean viscosity. Therefore, to solve Eq. 4, we need to know the radial profile of $\bar{\nu}$ at each time-step as a function of the corresponding $\Sigma_{\text{g}}$, i.e. $\bar{\nu}=\bar{\nu}(\Sigma_{\text{g}},\,R)$. To do this, we solve the vertical structure of the disc only at the beginning of the simulation, for a fixed value of the $\alpha$-parameter and varying $\dot{M}_{\text{st}}$ from $10^{-16}~{}\text{M}_{\odot}/\text{yr}$ to $10^{-4}~{}\text{M}_{\odot}/\text{yr}$ using a dense grid between these values. For each $\dot{M}_{\text{st}}$ –solving the vertical structure of the disc– we calculate the radial profiles of $\Sigma_{\text{g}}$ and $\bar{\nu}$ assuming that at each radial bin the disc is in a steady state (i.e. $\dot{M}_{\text{st}}=3\pi\nu\Sigma_{\text{g}}$). Then, from the initial gas surface density radial profile we interpolate at each radial bin the corresponding value of $\dot{M}_{\text{st}}$. With this information we can then interpolate the value of the corresponding mean viscosity at each radial bin. Given the initial gas surface density and mean viscosity radial profiles we solve for the first time-step the diffusion equation, obtaining a new radial profile for the gas surface density, and repeating the aforementioned mechanism. Finally, we note that Eq. 4 is solved using an implicit Crank-Nicholson method considering zero torques as boundary conditions, which is equivalent to consider zero density as boundary condition. We also remark that the radial profiles time evolution of the other relevant thermodynamic quantities at the disc mid-plane (temperature, pressure, opacities, etc.), needed to calculate the planet migration rates, are calculated in the same way as the mean viscosity radial profiles (see the extended version on the arXiv of Migaszewski, 2015, for details).

Throughout this work, the initial surface density of the disc is assumed to be equal to (Andrews et al., 2010):

[TABLE]

where $R_{\text{c}}=39$ au is the characteristic disc radius, $b=1$, and $\Sigma_{\text{g}}^{0}$ is a normalisation constant. The latter is a function of the disc mass $\text{M}_{\text{d}}$. In all our simulations, the characteristic disc radius and the exponent $b$ are kept fixed. For our fiducial model, we assume a disc with a mass of $0.05~{}\text{M}_{\odot}$, and we adopt a value of $\alpha=10^{-3}$ for the Shakura-Sunyaev $\alpha$-viscosity parameter (Shakura & Sunyaev, 1973). As we are interested in the early stages of planet formation, we only analyse the first million years of the disc evolution.

2.2 The solid component of the disc

In PLANETALP, the solid component of the disc is characterised by a particle population (planetesimals or smaller particles as pebbles), represented by a solid surface density. The initial solid surface density is given by

[TABLE]

Here $z_{0}$ is the initial dust to gas ratio (in this work we assume $z_{0}=0.0153$, the initial solar metallicity Lodders et al., 2009), and $\eta_{\text{ice}}$ represents the sublimation of volatiles inside the ice-line. $\eta_{\text{ice}}$ is equal to $1$ and $1/3$ outside and inside the ice-line, respectively. In our model, the ice-line is located at the distance to the central star where the initial temperature at the disc mid-plane drops to 170 K (for the fiducial disc it is located at $\sim 3$ au). We note that the ice-line does not evolve in time, and it is only used to define the initial solid surface density (which decreases as we move far away from the central star).

This population evolves according to the drift due to the nebular gas considering the Epstein, Stokes and quadratic regimes (Rafikov, 2004; Chambers, 2008), and according to the accretion by the embryos (see next section). In general, small particles as pebbles are always in the Epstein regime along the disc; small planetesimals are either in the Stokes or quadratic regimes; and big planetesimals ($\gtrsim 10$ km) are in the quadratic regime along the disc. For each radial bin, the drift velocities of the solid particles are calculated corresponding to the drag regime. In the case of planetesimals, we also take into account ejection/dispersion processes (see Ronco et al., 2017, for the approach employed).

As in previous works, we take the particle eccentricities and inclinations out of equilibrium considering the planet gravitational excitations (Ohtsuki et al., 2002), and the damping due to the nebular gas drag (Rafikov, 2004; Chambers, 2008). The numerical treatment of these phenomena is described in detail in Guilera et al. (2014) and Guilera & Sándor (2017). Although PLANETALP allows us to model the collisional evolution of the particle population (Guilera et al., 2014; San Sebastián et al., 2018), due to its high computational cost, this effect is not taken into account in the present work.

Finally, the time evolution of the solid surface density is calculated by using the following continuity equation:

[TABLE]

where $\dot{\Sigma}^{\text{tot}}_{\text{p}}(R,r_{p_{j}})$ represents the sink terms due to the accretion (and ejection/dispersion in the case of planetesimals) by the planets (see Ronco et al., 2017, for numerical details) and $v_{\text{mig}}(R,r_{p_{j}})$ is the particle drift velocity. The $r_{p_{j}}$ dependence emphasises that Eq. 7 is solved independently for each particle size for a given planetesimal distribution. In this case, the total planetesimal surface density is given by: $\Sigma_{\text{p}}(R)=\sum_{j}\Sigma_{\text{p}}(R,r_{p_{j}})$. For simplicity, we adopt a single size distribution in this work and we solve Eq. 7 using an implicit Donor cell algorithm considering zero density as boundary conditions.

2.3 Growth of the planets

In PLANETALP, planets grow by the continuous and simultaneous accretion of solid material and the surrounding gas (see Ronco et al., 2017, for details).

On the one hand, in the case of planetesimals (solid particles with Stokes number111The Stokes number of a particle is defined as $S_{t}=\Omega~{}t_{\text{stop}}$ (Stokes, 1851; Whipple, 1972; Weidenschilling, 1977), being $t_{\text{stop}}$ the stopping time that depends of the drag regime (see Guilera et al., 2014, for the adopted stopping times). greater than unity), the solid cores of the planets grow according to the accretion rate for oligarchic growth derived by Inaba et al. (2001),

[TABLE]

where $\text{M}_{\text{C}}$ is the mass of the core, $\Sigma_{\text{p}}(R_{\text{P}})$ the surface density of solids at the planet location $R_{\text{P}}$, $\text{R}_{\text{H}}$ the planet Hill radius, $P_{\text{coll}}$ the collision probability, and $\Omega_{\text{P}}$ the Keplerian frequency. $P_{\text{coll}}$ is a function of the planet capture radius, considering the drag force experienced by the planetesimals when entering the planetary envelope following Inaba & Ikoma (2003), the planet Hill radius, and the relative velocity between the planetesimals and the planet (assuming a coplanar and circular orbit for the latter). There are three regimes for $P_{\text{coll}}$ depending on whether the relative velocities are low, intermediate or high. This in turn depends on the planetesimal eccentricities and inclinations (see Guilera et al., 2010, 2014, for details). In general, small planetesimals ($r_{p}\lesssim 1$ km) are in the low relative velocity regime, while big planetesimals ($r_{p}\gtrsim 10$ km) are in high relative velocity regime. Thus, small planetesimals have greater $P_{\text{coll}}$ and so, greater accretion rates (typically of the order of $10^{-5}~{}\text{M}_{\oplus}/\text{yr}$–$10^{-4}~{}\text{M}_{\oplus}/\text{yr}$) with respect to big planetesimals (which typically have accretion rates of the order of $10^{-6}~{}\text{M}_{\oplus}/\text{yr}$)222See Guilera et al. (2010, 2011, 2014) and Fortier et al. (2007, 2013) for classical planetesimal accretion rates in the framework of the oligarchic growth depending on planetesimal sizes and disc masses.. The Eq. 8 is radially integrated along the planetesimal feeding zone, which extends $4.5~{}\text{R}_{\text{H}}$ at each side of the planet location, using a normalisation function (Guilera et al., 2014). We do not consider some potentially relevant phenomena that could change the planetesimals accretion rates, such as planetesimal resonant trapping or gap opening in the planetesimal disc (Tanaka & Ida, 1999; Shiraishi & Ida, 2008). These phenomena are difficult to incorporate without considering N-body interactions between the planets and the planetesimals.

On the other hand, if the solid particles have Stokes numbers $S_{t}$ lower or equal than unity we consider them as pebbles, and the accretion rates are given by Lambrechts et al. (2014):

[TABLE]

where $\beta=\text{min}(1,\text{R}_{\text{H}}/H_{\text{p}})$ takes into account a reduction in the pebble accretion rates if the scale-height of the pebbles ($H_{\text{p}}$) becomes greater than the planet Hill radius. This scale-height is defined as $H_{\text{p}}=H_{\text{g}}\sqrt{\alpha/(\alpha+S_{t})}$ (Youdin & Lithwick, 2007), where $H_{\text{g}}$ is the scale-height of the gaseous disc. Eq. 9 is radially integrated along the pebble feeding zone, which extends $1~{}\text{R}_{\text{H}}$ at each side of the planet location (Ormel & Klahr, 2010). While planetesimals can be accreted by a fraction of the Hill sphere of the planet, pebbles can be accreted by the full Hill sphere producing hence higher accretion rates. We do not take into account the pebble isolation mass (Lambrechts et al., 2014; Bitsch et al., 2018) or the possible pebble ablation in the planet envelope (Brouwers et al., 2018). The first phenomenon stops the accretion of pebbles at high planetary masses, while the second could strongly inhibit the formation of massive cores by pebble accretion.

Regarding the gas accretion, we follow the approach used in Ronco et al. (2017) and Guilera et al. (2017). We consider that if the planet core mass reaches a critical mass equal to:

[TABLE]

where $\dot{\text{M}}_{\text{C}}$ is the solid accretion rate, then the gas accretion rate onto the planet is given by:

[TABLE]

where $\dot{\text{M}}_{\text{KH}}=\text{M}_{\text{P}}/\tau_{\text{g}}$, being $\text{M}_{\text{P}}$ the total mass of the planet and $\tau_{\text{g}}=8.35\times 10^{10}\,(\text{M}_{\text{P}}/\text{M}_{\oplus})^{-3.65}\,\text{yr}$ the characteristic Kelvin-Helmholtz growth timescale of the envelope (Ronco et al., 2017). Additionally, the maximum rate at which the gas can be delivered by the disc onto the planet is equal to:

[TABLE]

where $\Sigma_{\text{g,P}}$ is the local gas surface density (hereafter by local we mean at the planet location). Finally, if the planet is able to open a gap in the gaseous disc (Tanigawa & Ikoma, 2007), the gas accretion onto the planet $\dot{\text{M}}_{\text{GP}}$ reads as follows:

[TABLE]

where

[TABLE]

with $\Sigma_{\text{acc}}=\Sigma(x=2\text{R}_{\text{H}})$ and,

[TABLE]

wherein $l$ and $x_{\text{m}}$ are defined as

[TABLE]

and

[TABLE]

being $\text{M}_{\star}$ the stellar mass, and $H_{\text{g,P}}$ the local gaseous disc scale height.

2.4 Planetary orbital evolution

As we mentioned before, as a planet grows, the exchange of angular momentum with the gaseous disc produces a torque that results in the planet’s migration along the disc. If the planet is not able to open a gap in the disc, then we consider that it migrates under type I migration. Its stellocentric distance varies as follows:

[TABLE]

where $\Gamma$ is the total torque on the planet, $\mathcal{L}_{\text{P}}$ the planet angular momentum, and $R_{\text{P}}$ the planet semi-major axis (assuming that it always remains in a nearly circular orbit).

When the planet becomes massive enough to open a gap in the gaseous disc, it enters in the type II migration regime. We follow the criterion derived by Crida et al. (2006) for the transition between type I and type II migration. The treatment of the latter regime in our model is described in detail in Ronco et al. (2017).

2.4.1 Previous type I migration rates

In Guilera et al. (2017), we implemented the type I migration rates derived by Paardekooper et al. (2010, 2011) to study the formation of giant planets in wide orbits. The total torque is given by:

[TABLE]

where $\Gamma_{\text{L}}$ and $\Gamma_{\text{C}}$ are the Lindblad and corotation torques, respectively. The former is equal to:

[TABLE]

where $\alpha^{\prime}=-{\rm d}\ln\Sigma_{\text{g,P}}/{\rm d}\ln R$ and $\beta^{\prime}=-{\rm d}\ln T_{\text{P}}/{\rm d}\ln R$, with $T_{\text{P}}$ the local disc temperature at the mid-plane. The normalised torque $\Gamma_{0}$ is given by:

[TABLE]

where $q=\text{M}_{\text{P}}/\text{M}_{\star}$ and $h=H_{g,\text{P}}/R_{\text{P}}$. Finally, $\gamma_{\text{eff}}$ reads:

[TABLE]

where $\gamma=7/5$ is the specific heat ratio for the diatomic gas, and $Q=2\chi_{\text{P}}\Omega_{\text{P}}/3hc_{s,\text{P}}^{2}$. $c_{s,\text{P}}$ and $\chi_{\text{P}}$ are the local sound speed and the local thermal diffusion coefficient at the disc mid-plane, respectively. The latter is computed as follows333As it was noted by Bitsch & Kley (2011), there is a typo in the Eq. 34 in Paardekooper et al. (2011), the factor in the numerator has to be 16 instead of 4.:

[TABLE]

where $\kappa_{\text{P}}$, $\rho_{g,{\text{P}}}$, and $P_{\text{P}}$ are the local opacity, local volumetric gas density, and local pressure of the disc also at the mid-plane (respectively).

The corotation torque has two main contributions: a barotropic component and an entropy one. Therefore,

[TABLE]

The detailed expressions of these two terms are described in Appendix A. It is important to note that the type I migration rates from Paardekooper et al. (2011) are the most used rates in planet formation models that consider the evolution of non-isothermal discs.

2.4.2 Updated type I migration rates

In this work we incorporate for the first time the new type I migration rates recently proposed by Jiménez & Masset (2017) and Masset (2017). We consider that the total torque on the planet is given by:

[TABLE]

where $\Gamma_{\text{Type~{}I}}$ is the sum of the Lindblad and corotation torques (possibly saturated) on the planet, while $\Gamma_{\text{thermal}}$ is the additional torque on the planet that arises from thermal effects (the perturbation of the flow in the planet’s vicinity due to a finite thermal diffusivity and that arising from the heat released in its surroundings by the accreting planet). Jiménez & Masset (2017) recently studied the first term of the right side of Eq. 25 through 3D hydrodynamical simulations, improving the previous study of Masset & Casoli (2010). They provide the following expression: $\Gamma_{\text{Type~{}I}}=\Gamma_{\text{L}}+\Gamma_{\text{C}}$, where the Lindblad torque is now given by:

[TABLE]

where $\alpha^{\prime}$, $\beta^{\prime}$, $\Gamma_{0}$ and $\chi_{\text{P}}$ are the same quantities aforeintroduced. We see that while the dependence of the Lindblad torque on the radial surface density gradient is the same as in Paardekooper et al. (2011), the dependence on the temperature radial gradient is slightly smaller. The function $f(\chi_{\text{P}}/\chi_{\text{C}})$ plays the role of the inverse of an effective adiabatic index $\gamma_{\text{eff}}$, as in Paardekooper et al. (2011) and it reads as follows:

[TABLE]

wherein $x=\chi_{\text{P}}/\chi_{\text{C}}$, and where $\chi_{\text{C}}=R_{\text{P}}^{2}h^{2}\Omega_{\text{P}}$ is the critical diffusivity. The corotation torque is now given as the sum of four contributions: three of them associated with disc radial gradients (vortensity, entropy, and temperature) and a fourth contribution which arises from a viscous production of vortensity. The sum of the three contributions arising from radial gradients may seem surprising at first glance, as the disc radial profile is determined by two quantities (e.g. the density and temperature) and two indexes (e.g. $\alpha^{\prime}$ and $\beta^{\prime}$). While the first three components of the corotation torque could be compacted into two components, each has different saturation properties, and the result would obfuscate the physical origin of each of the components. This has been discussed by Jiménez & Masset (2017) in section 4.3. Thus,

[TABLE]

These components are also described in detail in Appendix A. We remark that these new type I migration prescriptions are derived for low-mass planets ($q<0.2h^{3}$) and intermediate-mass planets ($0.2h^{3}\lesssim q\lesssim 2h^{3}$). For more massive planets able to deplete significantly the coorbital region this migration recipe could not be accurate enough. But such massive planets should be near the transition to the type II regime.

The second term of the right side of Eq. 25 was recently studied by Masset (2017) from the analytical perspective. This work considers luminous and non-luminous planets and derives a semi-analytical prescription for $\Gamma_{\text{thermal}}$, valid for both regimes. If the planet does not release energy into the disc, two cold and dense lobes on either side of the orbit appear due to the thermal diffusion of the disc. This effect exerts a torque on the planet, comparable in magnitude to the Lindblad torque. Masset (2017) defines this as a cold torque and argues that it corresponds to the cold fingers reported by Lega et al. (2014). Alternatively, if the planet releases energy into the disc, hot and low-density lobes are generated in the surroundings. In this case, a heating torque — of opposite sign with respect to the cold torque — is exerted on the planet (Benítez-Llambay et al., 2015). Here, we implement the thermal torque in the same way as proposed by Masset (2017), i.e. as the sum of the cold and the heating torque. The former is given by:

[TABLE]

where $x_{p}=\eta h^{2}R_{\text{P}}$, with $\eta=\alpha^{\prime}/3+\beta^{\prime}/6+1/2$, and $\lambda_{c}=\sqrt{\chi_{\text{P}}/q\Omega_{\text{P}}\gamma}$. Alternatively:

[TABLE]

where the critical luminosity $\text{L}_{\text{C}}$ is given by:

[TABLE]

and the luminosity released by the planet, $\text{L}_{\text{P}}$, due to the accretion of planetesimals reads:

[TABLE]

being $\text{R}_{\text{C}}$ the radius of the planet’s core. We remark that the total luminosity of the planet is given by the accretion of solids (which is a reasonable approximation until the planet accrete a significant envelope) and that the radius of the planet’s core increases as planet grows by the accretion of solid material considering a constant density for the core of $3.2~{}\text{g}/\text{cm}^{3}$. Thus, the thermal torque can be written as follows444The normalised torque ($\Gamma_{0}$) defined in Masset (2017) is not equal that the one defined in this work, so our Eq. 33 does not have a factor $h$ multiplying $\Gamma_{0}$ as in Eq. 146 from Masset (2017).:

[TABLE]

Finally, we note that the thermal torque is significant while the heat released by the planet generates an excess of internal energy outside of the Bondi sphere. Masset (2017) argues that this condition is satisfied when the time-scale for heat diffusion across the Bondi radius is shorter than the acoustic time. This condition translates into $\text{M}_{\text{P}}<\text{M}_{\text{crit}}^{\text{thermal}}$, where the critical thermal mass is given by

[TABLE]

If the planet mass is greater than this critical mass, it is not guaranteed that the internal energy injected in the gas near the planet emerges as an excess of internal energy outside of the Bondi sphere. In such case, a cut-off of the thermal torque is expected. From a numerical point of view, when $\text{M}_{\text{P}}>\text{M}_{\text{crit}}^{\text{thermal}}$ we consider that $\Gamma_{\text{thermal}}=0$.

3 Results

The aim of this work is to analyse how the updated type I migration rates described in Sec. 2.4.2, and especially the thermal torque, impacts the process of planet formation. To do so, we first construct migration maps and then model the formation of a planet at different initial locations and for different protoplanetary disc models.

3.1 Migration maps for the updated type I migration rates

3.1.1 Comparison between Paardekooper et al. (2011) and Jiménez &

Masset (2017) type I migration rates

Here we compare the type I migration rates derived by Paardekooper et al. (2011) and Jiménez & Masset (2017) constructing migration maps for both cases. For the calculation of the torques of the disc on the planet we require the following quantities: the gas surface density, the mean viscosity, the scale height of the disc, the temperature, the volumetric gas density, the sound speed, the pressure, and the opacity at the mid-plane of the disc. In Fig. 1 we plot the radial profiles of the aforementioned quantities from our fiducial simulation555Recall from Sect. 2.1 that for our fiducial model considers $R_{\text{c}}=39$ au, $\gamma=1$, $\text{M}_{\text{d}}=0.05~{}\text{M}_{\odot}$ and $\alpha=10^{-3}$., at four different stages of the disc early evolution (50 Kyr, 250 Kyr, 500 Kyr and 1 Myr).

In Fig. 2 we show the migration maps calculated using the type I migration rates by both Paardekooper et al. (2011) (left) and Jiménez & Masset (2017) (right panels). Each row corresponds to a different time, from 50 Kyr (top) to 1 Myr (bottom). These maps were constructed using $1000\times 1000$ radial and mass bins, both logarithmic equally spaced. The disc extension and the mass range are comprised between 0.1 au and 1000 au and between $0.1~{}\text{M}_{\oplus}$ and $1000~{}\text{M}_{\oplus}$, respectively. The green and blue regions of the maps indicate inward and outward migration of the planet, respectively. The white region in the migration maps indicates that the migration rate is very low, reaching the bottom of the colour scale. The shadowed region for each panel at high masses represents the planet mass at which a gap is opened in the disc, following Crida et al. (2006). This implies that planets above this line migrate in our simulations according to type II migration rates. The dashed line across each panel represents the planet mass for which the condition $q=2h^{3}$ is satisfied. We observe that this curve is close to the transition between type I and type II migration regimes. We see in Fig. 2 that the migration maps generated with the two prescriptions considered are quite different, especially the location of the outward migration (blue) regions. Using the migration recipes from Paardekooper et al. (2011), two outward migration zones appear. One at $\sim 3$ au for planet masses between $1\lesssim\text{M}_{\text{P}}/\text{M}_{\oplus}\lesssim 10$, and an inner one at $R<1$ au for low-mass planets ($1\lesssim\text{M}_{\text{P}}/\text{M}_{\oplus}\lesssim 5$). Moreover, these regions of outward migration remain roughly constant over the first $1$ Myr of the disc evolution. On the other hand, the model by Jiménez & Masset (2017) does not present blue regions of outward migration during the first $1$ Myr. The comparison between the two sets of migration maps shows that the corotation torques predicted by the migration recipes from Paardekooper et al. (2011) are higher than the corresponding predicted by Jiménez & Masset (2017).

In Fig. 3, we plot the normalised difference between the migration maps for the four different times, defined as:

[TABLE]

where $\dot{a}_{\text{JM2017}}$ and $\dot{a}_{\text{P2011}}$ represent the migration rates from Jiménez & Masset (2017) and Paardekooper et al. (2011), respectively. In these panels, the yellow-black (yellow-red) colour palette represents a negative (positive) difference. We clearly see that there are significant differences between both migration prescriptions for the four different times. These are particularly important in the inner disc regions at $R\lesssim 20$ au, where planet formation is expected to take place. These differences also correspond to a wide range of planetary masses, from $0.1~{}\text{M}_{\oplus}$ up to tens of Earth masses. Moreover, the normalised difference becomes extremely high in the regions where both migration recipes predict opposite migration directions. As a consequence, each prescription is expected to produce a different early planetary evolution. As we will show in next sections, these differences in the migration rates have a dramatic impact on planet formation.

It is well known that the Shakura-Sunyaev $\alpha$-viscosity parameter and the disc mass (which defines the initial gas surface density) play an important role in the time and radial disc evolution. This in turn affects the torques exerted on the planet. In Fig. 4, we plot the dependence of the normalised difference between the migration maps at 250 Kyr with the $\alpha$-viscosity parameter. We compare three different migration maps, one using the fiducial disc (panel b), and the other two using the same parameters as the fiducial disc except for the $\alpha$-viscosity parameter. In particular, we consider lower ($10^{-4}$, panel a) and higher ($10^{-2}$, panel c) values for $\alpha$. In Fig. 5 we plot, also at 250 Kyr, the dependence of the normalised difference with the disc mass. The panel b) represents The fiducial case is shown in panel b), while the disc model with lower ($\text{M}_{\text{d}}=0.01~{}\text{M}_{\odot}$) and higher ($\text{M}_{\text{d}}=0.15~{}\text{M}_{\odot}$) masses are shown in panels a) and c), respectively. Figures 4 and 5 show that –independently of the values of $\alpha$ and $\text{M}_{\text{d}}$– the migration maps are substantially different for the two recipes considered. This result is robust since it applies for a wide range of disc parameters. Additionally, in Fig. 4 and Fig. 5 we see that the black region moves to higher masses as the $\alpha$-viscosity parameter and the disc mass increase (i.e. as the viscosity increases), being more evident the dependence with the $\alpha$-viscosity parameter. This black region presents the major discrepancy between both models and it is related with the unsaturated horseshoe component of the corotation torque. For low masses, the corotation torque is linear and its value is similar from both recipes, while at high masses the corotation torque is saturated and it tends to vanish, also for both recipes (see Fig 13). The transition mass between the linear regime and the horseshoe drag regime of the corotation torque depends on the viscosity (the transition mass increases as the viscosity increases). In the same way, the mass for the saturation of the horseshoe drag also increases as the viscosity increases. This explains why the black region moves to higher masses as the viscosity increases.

3.1.2 Comparison between Jiménez &

Masset (2017) type I migration rates with and without Masset (2017) thermal torque

Here we compare the migration maps obtained for 3D adiabatic discs (Jiménez & Masset, 2017), with and without the inclusion of the thermal torque derived by (Masset, 2017). It is important to note that we do not add the thermal torque to the type I migration prescriptions derived by Paardekooper et al. (2011). As Masset (2017) pointed out, the thermal torque is derived taking into account how the disc thermal diffusion affects the torques over the planet. Paardekooper et al. (2011), derived their migration prescriptions based on 2D radiative hydrodynamics simulations wherein thermal diffusion was already included. Thus, the migration recipes derived by Paardekooper et al. (2011) could have partially incorporated the effect of the cold torque, although we can anticipate that it would have a very small value in their two-dimensional experiments with a finite smoothing length. This does not occur with the heating torque because Paardekooper et al. (2011) did not considered the heat released by the planet due to solid accretion. As mentioned in Section 1, we stress that the migration recipes derived by Jiménez & Masset (2017) do not take into account the effects of the cold and heating torques.

In Fig. 6, we plot the migration maps constructed employing the type I migration rates combined from Jiménez & Masset (2017) with the ones corresponding to the thermal torque from Masset (2017). As in Fig. 2, each panel corresponds to a different time, from 50 Kyr to 1 Myr. The green (blue) regions of the maps indicate an inward (outward) planet migration. Again, the shadowed region in each panel corresponds to the transition to the type II migration regime, and the dashed line corresponds to the condition $q=2h^{3}$. The dotted line in each panel corresponds to the critical thermal mass. These migration maps are constructed considering a constant solid accretion rate of $10^{-5}~{}\text{M}_{\oplus}/{\rm yr}$ for the thermal torque computation. We observe that the inclusion of the thermal torque introduces two regions of outward migration at low masses up to $\sim 2.5~{}\text{M}_{\oplus}$. The small blue island between 1 au and 3 au disappears as the disc evolves. However, the larger one between 4 au and 15 au, remains practically unchanged during the first Myr of the disc evolution. We note that in this region of outward migration, the transition to inward migration is smooth at the right side of the island. This means that in this specific region the Lindblad and cold torques dominate over the heating torque. On the other hand, at the left side of the blue island, the transition is more abrupt. This is due to the fact that in this case, as opposed to what happens on the right side, the planet reaches the critical mass for the thermal torque. The presence of this region of outward migration plays an important role in planetary evolution (see formation tracks in Sect. 3.2).

In Fig. 7, we plot the normalised difference between the type I migration rate prescriptions combined from Jiménez & Masset (2017) and Masset (2017) minus the ones derived only from Jiménez & Masset (2017). In the region where the planet mass is greater than the critical thermal mass (in bright yellow at the left of each panel), the type I migration prescriptions for both models are the same. The regions in red represent a positive normalised difference. These are regions of outward migration where the heating torque dominates the total torque. At some point, the combination between the Lindblad and the cold torque compensates the heating torque and an interface of small or null torque is generated (yellow transition between the red and black regions). This situation occurs due to the fact that, for a fixed solid accretion rate and a fixed mass of the planet, the heating torque decreases with distance to the central star (see Eqs. 30 and 31). The black region represents a negative normalised difference. This means that the total torque derived from the sum of Jiménez & Masset (2017) and Masset (2017) is more negative than the one derived only from Jiménez & Masset (2017). Here the cold thermal torque dominates over the heating torque. In other words, the planet’s luminosity is smaller than the critical value given by Eq. 31 and hence a more negative total torque is obtained. As we move far away from the central star, the contribution of the cold torque decreases (see Eq. 29) and the planet enters the Type I migration regime again (strong yellow).

In Fig. 8, we plot the migration maps at 250 Kyr for three different constant solid accretion rates, $10^{-6}~{}\text{M}_{\oplus}/{\rm yr}$ (left), $10^{-5}~{}\text{M}_{\oplus}/{\rm yr}$ (centre), and $10^{-4}~{}\text{M}_{\oplus}/{\rm yr}$ (right). As the heating torque is proportional to the luminosity of the planet generated by the accretion of the pebbles or planetesimals, the blue region of outward migration increases as the solid accretion grows. For a constant solid accretion rate of $10^{-4}~{}\text{M}_{\oplus}/{\rm yr}$, the blue outward migration region extends up to  30 au and up to a planet mass greater than $\sim 10~{}\text{M}_{\oplus}$. As before, at the right side of this region the transition to the inward migration regime is soft, but at the left side is abrupt because the planet reaches the critical mass for the thermal torque. As a planet grows, the solid accretion is expected to increase as well as long as there is enough available material in the planet feeding zone.

3.2 Effects on planet formation

In the previous sections, we constructed migration maps and showed that the migration rates predicted by different authors have significant differences, which can be of more than an order of magnitude for low-mass and intermediate-mass planets. In this section, we study the differences that arise from the adoption of either migration prescription in our planet formation model.

3.2.1 Planet formation affected by Paardekooper et al. (2011) and Jiménez &

Masset (2017) type I migration rates

In Fig. 9, we plot the planet formation tracks considering different initial locations and different initial pebble and planetesimal sizes. The black lines represent the formation tracks using the migration prescriptions from Paardekooper et al. (2011), while the red lines represent the tracks using the migration prescriptions from Jiménez & Masset (2017). It is important to note that we performed one simulation per formation track, i.e. we do not calculate the simultaneous formation of several planets. In our simulations, the proto-planet begins as a rocky core of $0.01~{}\text{M}_{\oplus}$ immersed in our fiducial disc. The simulations end when the planet either reaches the disc inner radius or after 3 Myr of evolution.

When we considered pebbles of 1 cm radius (panel a of Fig. 9), the planets always reach the inner part of the disc experiencing both type I migration recipes. For planets with initial locations lower than 9 au, the planet formation tracks are very different despite the fact that all of them end with a similar mass (of $\sim 25~{}\text{M}_{\oplus}$). We note that if the planet is initially located at 11 au, the planet formation tracks and masses are dramatically different. For planets at further distances, the formation tracks and final masses remain similar.

For the case of initial planetesimals of 100 m of radius (panel b of Fig. 9), the formation tracks and final masses are very different using both type I migration rates (in this case for all the initial planets locations). For initial planet locations beyond 7 au, the type I migration rates derived using the prescription of Jiménez & Masset (2017) seem to be higher than the corresponding rates derived using the prescription of Paardekooper et al. (2011). In all the simulations, the planets reach the inner radius of the disc, except for the planet initially located at 5 au using the migration recipes from (Paardekooper et al., 2011).

For the case of initial planetesimals of 1 km (panels c of Fig. 9), planet formation tracks present significant differences between both cases when planets grow up to a a few Earth masses.

Finally, for 100 km planetesimals (panels d of Fig. 9), the formation tracks are similar using the two migration prescriptions considered. This is due to the fact that planets only grow up to about the mass of Mars at most.

In Fig. 10, we plot the planet formation tracks using both migration prescriptions for initial planetesimals of 100 m of radius and the planet initial locations at 3 au and 5 au. These two cases are the ones showing the largest differences in the planet formation tracks presented in the panel b) of Fig. 9. We emphasise the role of the Lindblad and corotation torques for both migration planetary tracks. The light-blue and orange curves correspond to parts of the evolution in which either the Lindblad or the corotation torques dominate, respectively. It is clear that both prescriptions predict very different torques over the evolution of the planet. Adopting the migration recipes from Paardekooper et al. (2011), the corotation torque becomes dominant at planet masses between $\sim 2~{}\text{M}_{\oplus}$ – $\sim 5.5~{}\text{M}_{\oplus}$, and between $\sim 2.5~{}\text{M}_{\oplus}$ – $\sim 3.5~{}\text{M}_{\oplus}$ for planets initially located at 3 au and 5 au, respectively. This is a natural consequence of the blue islands in the left column of Fig. 2, located at $R\sim 3$ au. On the other hand, adopting the migration recipes from Jiménez & Masset (2017), the corotation torque never becomes dominant, and planets always migrate inward. That is the reason of the radically different formation tracks.

3.2.2 Planet formation including Masset (2017) thermal torque

In this section, we analyse the role of the thermal torque over the formation of a planet. In Fig. 11, we plot the formation tracks of the planets considering different initial locations and different initial planetesimal sizes. The blue lines represent the formation tracks using the migration prescriptions from Jiménez & Masset (2017) and Masset (2017), while the red lines only the ones from Jiménez & Masset (2017) (i.e. without the thermal torque). Again, we perform one simulation per formation track, where the planets are initially rocky cores of $0.01~{}\text{M}_{\oplus}$. The disc is the fiducial one and the simulations end when the planets reach the inner radius of the disc or at 3 Myr of evolution.

For the case where the solid disc is formed by pebbles of 1 cm (panel a of Fig. 11), and when the thermal torque is considered, the heating torque produces a significant outward migration. For the case of the planets initially located at $a_{\text{P}}\lesssim 3$ au, this migration is not significant, due to the fact that the planet reaches the critical thermal mass at a very low-mass and the thermal torque is then suppressed. Afterwards, the formation track is similar to the pure type I case (as the case of the planet initially located at 1 au). However, for planets initially located beyond 3 au the heating torque becomes dominant (see Fig. 12 that illustrates a similar situation) and the planets migrate outward until $\sim 24$ au. At this distance from the central star, the planets achieve the critical thermal mass (of about $10~{}\text{M}_{\oplus}$), and quickly migrate inward until reaching the disc inner radius. In these cases, the final planet masses are very different with respect to the case where the thermal torque is not considered. We note that not considering the pebble isolation mass does not affect the planet formation tracks before the planets reach the critical thermal mass since pebble isolation masses are normally higher than the critical thermal masses. However, this assumption could lead to an overestimation of pebble accretion rates at large masses.

When planetesimals of 100 m are considered (panel b of Fig. 11), only for planets initially located at 5 au and 7 au the heating torque becomes dominant666This happens because pebble accretion rates are higher than the planetesimal accretion rates, even for such small planetesimals. and planets migrate outward until they reach the critical thermal mass of about $10~{}\text{M}_{\oplus}$ at $\sim 23$ au. Then, the planets migrate inward and they reach the disc inner radius with final masses significantly greater than the ones corresponding to the case where the thermal torque is not included. For the planets initially located at 3 au and 9 au, the heating torque only acts for very low-masses, so planets do not migrate outward efficiently. In the former case, the planet quickly achieves the critical thermal mass and the thermal torque is suppressed. In the latter, the planetesimal accretion rate is not high enough to produce a heating torque greater than the combination of the cold and the Lindbland torques. This results into inward migration. This is because in the fiducial case, as we mentioned before, the snowline is located near 3 au. Thus, as we move far away from the central star the planetesimal surface density decreases and then the planetesimal accretion rate also decreases. For planets located beyond 9 au we obtain a faster inward migration compared to the case where the thermal torque is not considered. This is due to the fact that, as we mentioned before, planetesimal accretion rates are lower as we move far way from the central star, and thus the heating torque is negligible. So, the cold torque dominates the thermal torque. When added to the Lindbland one, the net torque is then more negative. Similar situations occur for the cases where planetesimals of 1 km and 100 km are considered (panels c and d of Fig. 11, respectively). However, in these cases the planets that efficiently migrate outward only migrated until $\sim 19$ au for 1 km planetesimals and $\sim 13$ au for 100 km planetesimals. In these cases, since planetesimal accretion rates decrease as planetesimal sizes increase, the cold torque plus the Lindblad torque reverse the migration direction of the planet at lower distances with respect to the case of smaller planetesimals or pebbles.

Finally, in Fig. 12 we plot the formation tracks for both migration models, for 100 m planetesimals and for a planet initially at 5 au. In this case, we compare the Lindblad (light-blue) and the thermal torques (yellow) by considering their absolute magnitudes. In this case, the corotation torque never becomes dominant. We observe that initially the thermal torque, and specifically the heating torque, dominates the migration track of the planet. This is due to the fact that planetesimal accretion rates are high for small planetesimals of 100 m at 5 au, while the cold and Lindblad torques are proportional to the mass of the planet, which is low initially. Thus, the planet migrates outward until it reaches $\sim 10~{}\text{M}_{\oplus}$ at $\sim 23$ au. When this happens, the Lindblad torque becomes dominant and the migration direction of the planet is reversed until it reaches the disc inner radius with a final mass of $\sim 50~{}\text{M}_{\oplus}$. For the case where the thermal torque is not considered, as we mentioned before, the Lindblad torques always dominate and the planet always migrates inward until reaching the inner disc radius with a final mass of $\sim 10~{}\text{M}_{\oplus}$.

4 Summary and discussion

In this work, we incorporated the most updated type I migration rate prescriptions to analyse their role in the process of planet formation. These prescriptions have two components: one related to the classical type I torque, in which the role of thermal diffusion within the disc is limited to its non-linear effect on the saturation of the corotation torque, and another related to the thermal torque. The latter takes into account the perturbation of the disc’s flow arising from thermal diffusion and the heat released by the planet due to the accretion of solid material, in the linear limit. The first ones are derived from 3D hydrodynamical simulations (Jiménez & Masset, 2017), improving the calculation of the corotation torque respect to previous 2D calculations (Masset & Casoli, 2010; Paardekooper et al., 2011). The thermal torques are derived from the analytical study by Masset (2017), motivated by previous hydrodynamical simulations (Lega et al., 2014; Benítez-Llambay et al., 2015).

We first compared the migration maps constructed using type I migration prescriptions from 2D and 3D hydrodynamical simulations for a fiducial case. For the 2D hydrodynamical simulations we adopted the commonly used type I prescriptions of Paardekooper et al. (2011), while for the 3D hydrodynamical simulations we adopted the new prescriptions derived by Jiménez & Masset (2017). We showed that the corresponding maps present significant differences for the regions of outward and inward migration, which are related to the corotation torques. We also analysed the dependence of the migration maps with the $\alpha$-viscosity parameter and the disc mass finding large differences between both type I migration prescriptions in a wide range of such parameters. This is particularly striking in the disc inner regions, at $R\lesssim 20$ au, where planet formation is expected to take place.

Then, we added the calculation of the thermal torque derived by Masset (2017) to the updated prescriptions by Jiménez & Masset (2017), for the first time to the best of our knowledge. We find that the inclusion of the thermal torque introduces a new region of outward migration, whose size directly depends on the planet solid accretion rate because the heating torque is proportional to this quantity. We also showed that at some point the sum of the Lindbland torque and the cold torque counteracts the total positive torque generated by the heating torque. Then, due to the cold torque, the total torque becomes more negative compared to the pure type I migration case. It is important to stress that, as in Masset (2017), the critical thermal mass given by Eq. 34 is only an order of magnitude estimate. More detailed studies on the mass limit for which the planet stops to efficiently release heat on the disc are needed. In addition, the heating torque is calculated assuming that the heat released in the disc is only due to the solid accretion. This is certainly true for a low-mass or intermediate-mass planet. However, if the planet has a significant envelope, the heat released due to its contraction could also play an important role. We note that in this version of PLANETALP, the equations of transport and structure of the envelope are not calculated, and we implemented a parametric gas accretion model as in Ronco et al. (2017). Thus, we are not able to calculate in an accurate way such additional production of heat. This will be studied in detail in a future work.

After the comparison between the migration maps, we studied the planet formation affected by the different type I migration rates. We first analysed the planet formation employing Paardekooper et al. (2011) and Jiménez & Masset (2017) prescriptions. To do so, we calculated the planet formation tracks adopting the fiducial disc for planets located at different initial locations, and for different pebble and planetesimal sizes. We showed that the planet formation tracks are systematically different, especially in the models where the planets grow up to several Earth masses or more. This is precisely when the corotation torque becomes significant.

Then, we calculated the formation tracks considering the inclusion of the thermal torque. For the case of pebbles, we found that all the planets initially located beyond 3 au suffered a significant outward migration due to the thermal torque, and specifically to the heating torque, because of the high accretion rates onto them. In contrast to Bitsch et al. (2015b), we found that the planet formation tracks are remarkably different compared to the case where the thermal torque is not considered. As a result, the planets end with significantly different final masses. We associate this discrepancy with the fact that, in this early work published before the release of an analytic expression for the thermal torque, they used an empirical fit of the heating torque measured by Benítez-Llambay et al. (2015), which was obtained for planets with masses comparable to or larger than the critical mass. These could had been therefore subjected to thermal torques smaller than the ones predicted by Masset (2017). Other discrepancies can arise from the solid accretion rates involved in both models: the higher the accrection rate, the stronger the heating torque. More recently, Chrenko et al. (2017) computed the simultaneous formation of four super-Earths (initially located between 5 au and 10 au) by the accretion of pebbles through 2D radiative hydrodynamical simulations. These calculations also included the heat released by the planets due to the pebble accretion onto the nearby gas disc. They analysed the orbital evolution and growth in three different scenarios: without pebble accretion, with pebble accretion but without heating, and with accretion plus the heating torque. For the latter case, their results are different with respect to the first two cases regarding for the outward migration due to heating torque.

For planetesimals of 100 m of radius, the significant planet outward migration only occurs for planets initially located at 5 and 7 au. This is because the snowline of the fiducial disc is slightly beyond 3 au, so at larger distances from the central star the planetesimal accretion rates become smaller. For planets with initial locations beyond 7 au, we note that these migrated faster compared to the case where the thermal torque was not considered. This is because the cold (negative) torque is added to the Lindblad torque resulting in a more negative total torque. Similar situations occur when planetesimals of 1 km and 100 km are considered, with the particularity that in these cases the planets did not reach high masses as in the previous cases.

There are a number of shortcomings in our analysis. While we take into account the effects of the planet’s luminosity on the time evolution of its semi-major axis, we neglect its impact on the eccentricity and inclination, which is potentially quite strong (Eklund & Masset, 2017; Chrenko et al., 2017). Sizeable eccentricities and inclinations may significantly alter the migration paths of planetary embryos, and may also affect the accretion rate on the embryos (Chrenko et al., 2017) and the rate of collisions between embryos. As we consider only one planetary embryo at a time, such effect cannot be captured by this work. Our torque expressions assume that the torque is a function only of the underlying disc properties, but not of the planet’s migration rate itself. The existence of a migration feedback on the corotation torque, initially considered for sub-giant planets (Masset & Papaloizou, 2003), has been extended to lower mass planets (Paardekooper, 2014; Pierens, 2015; McNally et al., 2018). To the best of our knowledge, prescriptions for the dynamical corotation torque have not been used yet in planetary migration paths studies such as the one presented here. In spite of these shortcomings, our work shows that the new type I migration prescriptions by Jiménez & Masset (2017) and Masset (2017) are very different from the ones by Masset & Casoli (2010) and Paardekooper et al. (2011) and that the planet formation tracks are highly sensitive to these differences.

This calls for a reanalysis of the predictions and conclusions of previous work based on previous type I migration rates. We also emphasise that the thermal torque, and especially the heating torque, can lead to very different planet formation tracks when the solid accretion rates are high. In this case, the planet can experience a significant outward migration. As a first step, we studied the formation of only one planet in a fiducial disc. The consequences of considering a population of embryos and the impact of these new type I migration prescriptions on a population synthesis will be explored in a future work.

Acknowledgements

We thanks the anonymous referee who help us to improve the work with his/her comments. OMG and M3B are supported by the PICT 2016-0053 from ANPCyT, Argentina. OMG also acknowledges the hosting as invited researcher from IA-PUC. NC acknowledges financial support provided by FONDECYT grant 3170680. MPR and JC acknowledge financial support from the Iniciativa Científica Milenio (ICM) via the Núcleo Milenio de Formación Planetaria Grant. NC, MPR and JC acknowledge support from CONICYT project Basal AFB-170002. MM acknowledge financial support from the Chinese Academy of Sciences (CAS) through a CAS-CONICYT Postdoctoral Fellowship administered by the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile.

Appendix A Corotation torque prescriptions from 2D and 3D hydrodynamical simulations

As we mentioned in Sec. 2.4.1, the corotation torque derived by Paardekooper et al. (2011) from 2D hydrodynamical simulations is composed by two components, a barotropic and an entropy one, and thus $\Gamma_{\text{C}}=\Gamma_{\text{C,bar}}+\Gamma_{\text{C,ent}}$. The barotropic contribution is given by

[TABLE]

where $\Gamma_{\text{hs,bar}}$ is the fully unsaturated horseshoe drag part given by

[TABLE]

and $\Gamma_{\text{C,lin,bar}}$ is the linear barotropic corotation part given by

[TABLE]

and where the function $F$ governs the saturation, while the functions $G$ and $K$ govern the cut-off at high viscosity (see Paardekooper et al., 2011, for the explicit form of such functions). $p_{\nu}$ is the viscosity saturation parameter given by

[TABLE]

where $\nu_{\text{P}}$ is the viscosity at the planet location and $x_{s}$ is the non-dimensional half-width of the horseshoe region, given by

[TABLE]

and where Paardekooper et al. (2011) considered the softening parameter $b/h=0.4$ in order to mimic 3D effects. Finally, the entropy contribution to the corotation torque is given by

[TABLE]

where $\Gamma_{\text{hs,ent}}$ is the fully unsaturated entropy part, given by

[TABLE]

with $\xi=\beta^{\prime}-(\gamma-1)\alpha^{\prime}$, and $\Gamma_{\text{C,lin,ent}}$ is the linear entropy part given by

[TABLE]

The parameter $p_{\chi}$ is the saturation parameter associated with thermal diffusion, and it is given by

[TABLE]

On the other hand, Jiménez & Masset (2017) improved the expression of the corotation torque using 3D hydrodynamical simulations. In their expression, the latter has four components $\Gamma_{\text{C}}=\Gamma_{\text{C,vor}}+\Gamma_{\text{C,ent}}+\Gamma_{\text{C,temp}}+\Gamma_{\text{C,vv}}$, associated with the radial gradient of vortensity, entropy, and temperature, and with a viscously created vortensity. The different contributions for the corotation torque are given by

[TABLE]

where $\epsilon_{b}=1/(1+30hz_{\nu})$, with $z_{\nu}=R_{\text{P}}\nu_{\text{P}}/\Omega_{\text{P}}x_{s}^{3}$. A new result from Jiménez & Masset (2017) is the improvement in the calculation of the half width of the horseshoe region $x_{s}$ given now by

[TABLE]

where $h^{\prime}=h\sqrt{\gamma}$. $\Gamma_{\text{vor,lin}}$ is the linear part of the vortencity component and is given by

[TABLE]

while $\Gamma_{\text{vor,hs}}$ is the horseshoe drag part given by

[TABLE]

where $\mathcal{F}_{v}$ is the saturation function for the vortensity component of the horseshoe drag (see Jiménez & Masset, 2017, for the explicit form of $\mathcal{F}_{v}$), and

[TABLE]

is the vortensity unsaturated horseshoe drag component. The entropy component of the corotation torque is given by

[TABLE]

where $\epsilon_{\nu}=1/(1+6hz_{\nu})^{2}$ and $\epsilon_{\chi}=1/(1+15hz_{\chi})^{2}$, with $z_{\chi}=R_{\text{P}}\chi_{\text{P}}/\Omega_{\text{P}}x_{s}^{3}$. The linear part is given by

[TABLE]

where $\xi=\beta^{\prime}-0.4\alpha^{\prime}-0.64$. The horseshoe drag term of the entropy component is given by

[TABLE]

where $\mathcal{F}_{S}$ is saturation function of the entropy corotation torque (see Jiménez & Masset, 2017, for its explicit form) and the unsaturated horseshoe drag contribution is given by

[TABLE]

The temperature component is given by

[TABLE]

being the linear contribution

[TABLE]

and the temperature horseshoe drag

[TABLE]

where $\mathcal{F}_{T}$ is the saturation function (Jiménez & Masset, 2017), and the unsaturated horseshoe drag component is given by

[TABLE]

Finally, the last contribution to the corotation torque is given by

[TABLE]

where the explicit form of $g$ can be found in Jiménez & Masset (2017, Eq. 64).

In Fig. 13, we plot the normalised Lindblad, corotation and total torques as a function of the planet mass for a fixed planet location of 3.5 au, using both torque recipes for the case of the fiducial disc at 50 Kyr. We see that while the Lindblad torques remain very similar between both models, the corotation torque given by Paardekooper et al. (2011) overestimates the corresponding one from Jiménez & Masset (2017). Finally, in Fig. 14 we plot the components of the corotation torque given by Jiménez & Masset (2017). The functional forms of such components are very similar to the one presented by Jiménez & Masset (2017) (right panel of their figure 8). The differences in the values that reach the different components with respect to Jiménez & Masset (2017) arise from the difference in the disc model adopted.