© ESO, 2014

1. Introduction

The huge diversity found in the properties of extrasolar planets is challenging to reproduce with global theoretical planet formation models. The goal of such a model is to explain all the different planet types, which range from low-mass rocky planets such as Kepler-10 b (Batalha et al. 2011) and multiplanet systems like our solar system to high-mass planets orbiting far from their star, such as NR 8977 (Marois et al. 2008).

The only way to study this problem is to use the results of global formation and evolution models and to compare them statistically with the steadily increasing number of known planets and their physical properties. This is done in planet populations synthesis calculations, in which the evolution of one or several planets and the harboring protoplanetary disk is calculated at the same time in Monte Carlo simulations to create whole populations of planets. Several groups presented studies based on this method, for example, Ida & Lin (2004, 2010), Thommes et al. (2008), Miguel & Brunini (2008), Hellary & Nelson (2012), or by our group, Mordasini et al. (2009a,b, 2012a,b), and Alibert et al. (2011).

One general result found in all these models is that giant planets close to the star (“hot” Jupiters) do not form in situ: the extrapolation of disk properties found at larger distances to small distances indicates that there is probably not enough solid material close-in to form a sufficiently large core that whould be able to accrete gas. The amount of material a core can accrete locally is given by the isolation mass M_iso. According to the empirical minimum-mass solar nebula model (MMSN), the isolation mass is only a fraction of the Earth mass (M_⊕) inside of 1 AU (Ida & Lin 2004). Therefore an increase of solid matter by roughly two orders of magnitude compared with the MMSN would be needed (Ida & Lin 2004). This means that to explain the close-in “hot” Jupiters, they would have to form initially at larger separations from the star and move inward by some mechanism (e.g., planet-planet scattering Rasio & Ford 1996; Kozai mechanism Nagasawa et al. 2008, or tidal interactions with the gas disk Goldreich & Tremaine 1980; Tanaka et al. 2002). In the planet formation model used in this work, we consider only one core per disk, which means that only migration caused by tidal interactions can be studied. This migration is generally described by two different regimes that depend on the mass of the protoplanet. The first is type I migration for low-mass planets, which are too small to form a gap in the disk, and the second is type II migration for planets that open a gap (D’Angelo et al. 2002).

In our previous work (e.g., Mordasini et al. 2009a; Alibert et al. 2011), we used the results obtained for isothermal disks reported in Tanaka et al. (2002) to calculate type I migration rates. The migration rate in this model only depends on the disk surface density profile and not on the temperature profile of the disk because it assumes a globally isothermal disk. Migration rates obtained with this model always lead to rapid inward migration in disks with profiles similar to the MMSN. We showed (Mordasini et al. 2009a) that to obtain a synthetic population compatible with observations, one needs to artificially reduce the isothermal type I migration rate by a large factor. Ida & Lin (2004) used a similar type I migration prescription and found necessary reduction factors of ≲0.1 for the migration rate.

Recent studies of type I migration in 2D or 3D hydrodynamical simulations also found outward migration for some masses or semimajor axes, depending on the disk temperature structure (Masset et al. 2006; Paardekooper & Mellema 2008; Kley et al. 2009). More analytical work derived a formulation that could be used in planet population synthesis calculations (Casoli & Masset 2009; Masset & Casoli 2009, 2010). Finally, Paardekooper et al. (2010) derived a formalism for type I migration in adiabatic or locally isothermal disks and improved this even more in Paardekooper et al. (2011).

In the present work we describe a semi-analytical type I migration model that can be applied to a wider range of planet and disk properties than that of Paardekooper et al. (2011). For this we used the adiabatic and locally isothermal migration equations from Paardekooper et al. (2010) and ratios of relevant time scales to determine the transition between different regimes. We then studied the global consequences of the physics included in the migration model in new sets of population synthesis calculations.

As an overview of this work, the new semi-analytical model we created is introduced in Sect. 2, where we also compare torques obtained with this model with torques obtained with the model of Paardekooper et al. (2011) and data obtained in 3D radiative hydrodynamic simulations of Kley et al. (2009) and Bitsch & Kley (2011). We discuss in Sect. 3 a reference synthetic planet population calculated with the nominal model, and in Sect. 4 we study the effect of different migration models on planetary populations. In Sect. 5 we draw our conclusions and summarize our results.

2. Migration model

The migration module in our planet formation model (Alibert et al. 2004) distinguishes three main regimes, type I, disk-dominated type II, and planet-dominated type II migration (Armitage & Rice 2005). Low-mass planets up to typically a few 10 Earth masses migrate in the type I regime, followed by migration in disk-dominated type II regime for more massive planets, which finally pass into planet-dominated type II migration when they reach masses of typically 100−200 M_⊕. We first describe our improvements to the description of the type II regime and introduce the new type I migration model afterwards.

2.1. Type II model and outward migration

In the old model, the type II migration rate was calculated using the equilibrium flux of gas in the disk, which was always assumed to be directed inward (Mordasini et al. 2009a). Now the direction and rate of type II migration is numerically calculated by considering the radial velocity v_gas of the (nonequilibrium) flux of gas at the position of the planet. It therefore allows outward migration if the planet is in a part of the disk where the gas is flowing outward (Veras & Armitage 2004). The planetary migration rate is then $\mathit{ȧ}\mathrm{p}\mathit{,}\mathrm{T}\mathrm{2}\mathrm{=}{\mathit{v}}_{\mathrm{gas}}\mathrm{\times }\mathrm{Min}\mathrm{(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\Sigma }{\mathit{a}}_{\mathrm{p}}^{\mathrm{2}}\mathit{/}{\mathit{M}}_{\mathrm{p}}\mathrm{)}\mathit{,}$(1)where a_p is the semimajor axis of the planet, M_p is its mass, and Σ is the gas surface density.

This mechanism has been invoked to explain the formation of exoplanets with semimajor axes larger than 20 AU, which cannot have formed in situ via the core accretion model (Veras & Armitage 2004). However, we find in the syntheses using the nonequilibrium model that outward migration in type II is seldom important, and no large-scale net outward migration over more than ~1 AU typically occurs due to it. The reason is twofold:

The radius of velocity reversal (or of maximum viscous couple) R_MVC (Lynden-Bell & Pringle 1974) is relatively close-in only early in the disk evolution. But, at these early times, the planets have usually not yet grown to a mass regime in which they migrate via type II. The evolution of the disk leads to the subsequent recession of R_MVC to larger radii. This occurs faster than the growth of the planets. Therefore type II outward migration is a very rare event during the spreading phase of the disk: at the moment planets have grown massive enough to migrate in type II, they are most of the time already located inside the R_MVC.

Another chance of outward migration exists towards the end of the disk lifetime, when parts of the disk flow outward because mass is removed at the outer border due to external photoevaporation. The gas surface densities, however, are typically already quite low at the position of the planet at this moment, so that the reduction factor of the planet’s migration rate relative to the viscous velocity ∝$\mathrm{\Sigma }{\mathit{a}}_{\mathrm{planet}}^{\mathrm{2}}\mathit{/}{\mathit{M}}_{\mathrm{planet}}$ (see Alibert et al. 2005) is low as well, leading again to only modest amounts of outward migration. As in the original model, we assumed a linear reduction factor if ${\mathit{M}}_{\mathrm{planet}}\mathit{>}\mathrm{2}\mathrm{\Sigma }{\mathit{a}}_{\mathrm{planet}}^{\mathrm{2}}$ (“fully suppressed” planet-dominated type II), because this agrees a better with hydrodynamical simulations than a square-root dependence on the planetary mass (Alexander & Armitage 2009). Outward migration during effective photoevaporation is additionally limited because the remaining disk-lifetime is short.

2.2. Type I migration regimes

Here we discuss the migration of low-mass planets below a few tens of Earth masses. As mentioned above, one of the problems with the original description of orbital migration of low-mass planets is the short time scale found in isothermal type I migration (Tanaka et al. 2002), which resulted in too many close-in planets in planet population synthesis calculations (see also the comparison in Sect. 4.1.1). These rates had to be artificially reduced by correction factors to produce enough “cold” giant planets at larger semimajor axes (Mordasini et al. 2009b; Ida & Lin 2008) to fit the observations.

On the other hand, Paardekooper & Mellema (2006) and Kley et al. (2009) showed that migration of small-mass bodies is not always inward because in a nonisothermal disk migration can be directed outward for some masses. Outward migration in the Type I regime can also occur in isothermal disks if full MHD turbulence is considered (see also Uribe et al. 2011; Guilet et al. 2013), but these new effects are not considered here because they need to be studied in more detail more studies before they can be parametrized.

Recently, Paardekooper et al. (2010) derived semi-analytical formulas for migration in the limiting case of adiabatic disks. Here we combine different formulas that are valid in different thermodynamical regimes into a model that can be applied to the wide range of planet masses and thermodynamic properties of the disk that are needed for our population-synthesis models.

2.2.1. Type I fit formula

Nonisothermal migration rates are more complex than isothermal rates. Generally, the gravitational interactions of the planet and the gas disk can lead to three characteristic flow regions that produce different types of torques:

• Lindblad torques: gas sufficiently far from the planet is onlyslightly perturbed and orbits the star on nearly circular orbitsinside or outside the planet’s orbit. The gravitational interactiongenerated by the planet in the regions outside and inside thecorotation region produces the Lindblad torques.

• Corotation torque: if the orbit of a gas parcel is closer to the planet, its flow is more and more deflected, until it passes the planet’s orbit in front or behind the planet. Thus it forms so-called horseshoe orbits (see Fig. 1). The gas in the horseshoe orbits produces by deflection the so-called horseshoe drag or corotation torque. The strength of this torque depends on the thermodynamic regime of the interaction between the planet and the disk and the mass of the planet.

For typical properties (radially decreasing density and temperature) of the disk the Lindblad torques lead to inward migration, whereas the isolated torques from the corotation region can result in either inward or outward migration. For certain thermal and surface density profiles in the disk, these torques can be stronger than the Lindblad torques. Thus the combination of Lindblad and corotation torques can lead to either inward or outward migration, depending on their relative strength, which is determined by the disk properties.

The total torque Γ can be expressed in the following way, as shown in Tanaka et al. (2002), Paardekooper et al. (2010) or Masset & Casoli (2010): $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ \mathrm{\Gamma }\end{array}& \mathrm{=}& \frac{\mathrm{\Gamma }}{{\mathit{M}}_{\mathrm{p}}}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}{\mathit{a}}_{\mathrm{p}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{p}}^{\mathrm{2}}\mathit{,}\mathrm{where}\\ \begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}& \mathrm{=}& \end{array}$In these equation, $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathrm{\Gamma }\end{array}$ is the specific torque (torque per unit mass). With a small p we denote all quantities of or at the position of the planet. M_p is the total mass of the planet and a_p the semimajor axis, while M_∗ is the mass of the star, h = H/a is the aspect ratio of the disk with a vertical scale height H. Ω_p is the Keplerian frequency of the planet and Σ_p the gas surface density. The dimensionless factor C in the second part of Eq. (3) gives the direction and strength of the migration and is discussed below for different thermodynamical regimes (Sect. 2.2). The dimensionless number $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{C}\end{array}$ is proportional to the migration rate ȧ_p and inversely proportional to the number of orbits needed for a planet to migrate over the distance of its semimajor axis N_orb. It is also proportional to the ratio of the two time scales for migration τ_mig and orbital motion τ_orb

Additional parameters in the following equations are β_Σ as the local power-law exponent of the gas surface density (Σ ∝ r^β_Σ), β_T the local power-law exponent of the temperature (T ∝ r^β_T), and γ the adiabatic index (ratio of the heat capacities) of the gas.

Depending on the properties of the disk, C is a combination of various torque contributions of variable importance. Before we describe how we combined the contributions, we introduce expressions for their individual strength.

Paardekooper et al. (2010) derived that the Lindblad torque in an adiabatic disk is proportional to (their Eq. (47), part 1) ${\mathit{C}}_{\mathrm{Lind}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\gamma }}\left(\mathrm{-}\mathrm{2.5}\mathrm{+}\mathrm{1.7}{\mathit{\beta }}_{\mathit{T}}\mathrm{-}\mathrm{0.1}{\mathit{\beta }}_{\mathrm{\Sigma }}\right)\mathit{.}$(4)They also found that the horseshoe drag in the adiabatic case is proportional to (their Eq. (47), part 2) ${\mathit{C}}_{\mathrm{HS}\mathit{,}\mathrm{adia}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\gamma }}\left(\mathrm{1.65}\mathrm{+}{\mathit{\beta }}_{\mathrm{\Sigma }}\mathrm{(}\mathrm{9}\mathrm{-}\mathrm{7.9}\mathit{/}\mathit{\gamma }\mathrm{)}\mathrm{-}\mathrm{7.9}{\mathit{\beta }}_{\mathit{T}}\mathit{/}\mathit{\gamma }\right)\mathit{.}$(5)The coefficient C in the adiabatic regime due to the combination of the Lindblad and corotation torques is ${\mathit{C}}_{\mathrm{adia}}\mathrm{=}{\mathit{C}}_{\mathrm{Lind}}\mathrm{+}{\mathit{C}}_{\mathrm{HS}\mathit{,}\mathrm{adia}}\mathit{.}$(6)Paardekooper et al. (2010) also found that the total torque in a locally isothermal regime, where the temperature T is constant in time but not with semimajor axis, is proportional to (their Eq. (49)) ${\mathit{C}}_{\mathrm{loc}}\mathrm{=}\mathrm{-}\mathrm{0.85}\mathrm{+}\mathrm{0.9}{\mathit{\beta }}_{\mathit{T}}\mathrm{+}{\mathit{\beta }}_{\mathrm{\Sigma }}\mathit{.}$(7)Subtracting from this the Lindblad torque in the adiabatic regime, but setting γ = 1 (compare Paardekooper et al. 2010 Sect. 5.4¹), one finds the horseshoe drag part in the locally isothermal regime as ${\mathit{C}}_{\mathrm{HS}\mathit{,}\mathrm{loc}}\mathrm{=}\mathrm{1.65}\mathrm{-}\mathrm{0.8}{\mathit{\beta }}_{\mathit{T}}\mathrm{+}\mathrm{1.1}{\mathit{\beta }}_{\mathrm{\Sigma }}\mathit{.}$(8)Compared with Eq. (4), Masset & Casoli (2010) derived a partially different Lindblad torque. We study the effect of this weaker Lindblad torque in Sect. 4.2.1. ${\mathit{C}}_{\mathrm{Lind}\mathit{,}\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\gamma }}\left(\mathrm{-}\mathrm{2.5}\mathrm{+}\mathrm{0.5}{\mathit{\beta }}_{\mathit{T}}\mathrm{-}\mathrm{0.1}{\mathit{\beta }}_{\mathrm{\Sigma }}\right)\mathit{.}$(9)

2.3. Time scales

The proper mix of the above described torque contributions can be determined by investigating the relevant time scales. For instance a disk behaving adiabatically produces a different torque than a locally isothermal one. Here the relevant time scales are the cooling time in comparison to the dynamic time.

To decide in which subtype the planet belongs to, we compared four characteristic time scales in total. The different time scales are schematically shown in Fig. 1. In all the following estimates of the time scale, the important characteristic length scale is the width of the horseshoe region x_s given as (Masset et al. 2006; Baruteau & Masset 2008; Paardekooper et al. 2010) ${\mathit{x}}_{\mathrm{s}}\mathrm{=}\mathrm{1.16}{\mathit{a}}_{\mathrm{p}}\sqrt{\frac{\mathit{q}}{\mathit{h}\sqrt{\mathit{\gamma }}}}\mathrm{·}$(10)In this equation q is the ratio of the planet mass to the stellar mass.

The first two time scales we compared are the cooling time and the u-turn time to distinguish between the locally isothermal and the adiabatic regime. The u-turn time is the time a gas particle needs to undergo one turn in front or behind the planet. Its value is approximately given as (Baruteau & Masset 2008) ${\mathit{\tau }}_{\mathit{u}\mathrm{-}\mathrm{turn}}\mathrm{=}\frac{\mathrm{64}{\mathit{x}}_{\mathrm{s}}{\mathit{h}}^{\mathrm{2}}}{\mathrm{9}\mathit{q}{\mathit{a}}_{\mathrm{p}}{\mathrm{\Omega }}_{\mathrm{p}}}\mathrm{·}$(11)The cooling time of a gas blob undergoing a turn is calculated by solving the 1D equation (Kley et al. 2009) $\frac{\mathrm{d}\mathit{T}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }{\mathit{C}}_{\mathit{V}}}\frac{\mathit{\partial }}{\mathit{\partial a}}\mathrm{(}{\mathit{F}}^{\mathrm{)}}\mathit{,}$(12)where F is the heat flux, ρ the gas density, and C_V the heat capacity at constant volume. We assumed a cooling over the length l_cool, which is the minimum of H, and x_s, which corresponds to either horizontal cooling over the width of the horseshoe region or vertical cooling through the disk. The flux F in the optically thin case, when $\mathit{\rho \kappa }{\mathit{l}}_{\mathrm{cool}}\mathit{<}\sqrt{\mathrm{1}\mathit{/}\mathrm{8}}$ is F = τσT⁴, and in the optically thick case in the diffusion description $\mathit{F}\mathrm{=}\frac{\mathrm{4}\mathit{ac}{\mathit{T}}^{\mathrm{3}}}{\mathrm{3}\mathit{\rho \kappa }}\frac{\mathit{\partial T}}{\mathit{\partial a}}$ (Kley et al. 2009). This means that we have two different types of the cooling time scale with a smooth transition when the optical depth is 1/$\sqrt{\mathrm{8}}$${\mathit{\tau }}_{\mathrm{cool}}\mathrm{=}\frac{{\mathit{l}}_{\mathrm{cool}}\mathit{\rho }{\mathit{C}}_{\mathit{V}}}{\mathrm{8}\mathit{\sigma }{\mathit{T}}^{\mathrm{3}}}\left(\mathrm{8}\mathit{\rho \kappa }{\mathit{l}}_{\mathrm{cool}}\mathrm{+}\frac{\mathrm{1}}{\mathit{\rho \kappa }{\mathit{l}}_{\mathrm{cool}}}\right)\mathrm{·}$(13)In a similar way we compared the viscous time scale and the libration time scale to find out whether the horseshoe drag is saturated or not (Masset & Casoli 2010). The viscous time scale is ${\mathit{\tau }}_{\mathrm{visc}}\mathrm{=}\frac{{\mathit{x}}_{\mathrm{s}}^{\mathrm{2}}}{\mathit{\nu }}\mathit{,}$(14)while the libration time is given as (Baruteau & Masset 2008) ${\mathit{\tau }}_{\mathrm{lib}}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }{\mathit{a}}_{\mathrm{p}}}{\mathrm{3}{\mathrm{\Omega }}_{\mathrm{p}}{\mathit{x}}_{\mathrm{s}}}\mathrm{·}$(15)We denote the ratios of the relevant time scales as ${\mathit{s}}_{\mathrm{1}}\mathrm{=}{\mathit{f}}_{\mathrm{cool}}\frac{{\mathit{\tau }}_{\mathrm{cool}}}{{\mathit{\tau }}_{\mathit{u}\mathrm{-}\mathrm{turn}}}$(16)and ${\mathit{s}}_{\mathrm{2}}\mathrm{=}{\mathit{f}}_{\mathrm{visc}}\frac{{\mathit{\tau }}_{\mathrm{visc}}}{{\mathit{\tau }}_{\mathrm{lib}}}\mathrm{·}$(17)An equivalent approach can also be found in Casoli & Masset (2009) and Masset & Casoli (2009).

These time scales are typically order-of-magnitude estimates. Therefore, we introduced two factors, f_cool and f_visc, of order unity to adjust the point of transition between different regimes to obtain a better agreement with radiative hydrodynamic simulations. In this work we set f_cool = 1.0, since this proved toagree well with numerical results of migration rates, as we show in Sect. 2.5. There, we also study the influence of various values of f_visc. In general, we find that in nonirradiated α disks, planets transit from the locally isothermal into the adiabatic regime before the corotation torque saturates (Sect. 3.5).

We also note that for a given h, the ratio of the u-turn time and libration time is constant: ${\mathit{\tau }}_{\mathit{u}\mathrm{-}\mathrm{turn}}\mathrm{=}\frac{{\mathrm{1.16}}^{\mathrm{2}}\mathrm{8}}{\mathrm{3}\mathit{\pi }}\frac{\mathit{h}}{\sqrt{\mathit{\gamma }}}{\mathit{\tau }}_{\mathrm{lib}}\mathrm{\approx }\mathrm{1.14}\frac{\mathit{h}}{\sqrt{\mathit{\gamma }}}{\mathit{\tau }}_{\mathrm{lib}}\mathit{.}$(18)

2.4. Type I migration formula

To combine the migration rates of these different regimes we defined an arbitrary transition function z to obtain a smooth shift from the locally isothermal to the adiabatic regime as a function of the variable s₁: $\mathit{z}\mathrm{\left(}{\mathit{s}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathit{s}}_{\mathrm{1}}^{\mathit{b}}}\mathit{,}$(19)which has the properties that for s₁ → 0,z(s₁) → 1 and for s₁ → ∞,z(s₁) → 0. Furthermore, depending on the value of b, the transition from 1 to 0 occurs more or less quickly around s₁ = 1. Therefore, uncertainties in overlap of different regimes can be approximated with a lower b value. Additionally, the continuity of the transition function allows one to use longer numerical time steps in a simulation when the planet is close to the transition from one regime into another.

For individual transitions (e.g., locally isothermal to adiabatic, or type I into type II), other studies (Paardekooper et al. 2011; Masset & Casoli 2010) derived physically motivated transition functions. But the comparison of physically motivated transition functions Paardekooper et al. (2011) with our simple transition function defined above shows little difference in population synthesis models (see Sect. 4.1.2). We set b = 4.0, but the actual value of b is again not very important for the global outcome seen in a population if 1.5 ≲ b ≲ 100 (see Appendix B).

We multiplied the horseshoe part with min(1/s₂,1) to account for the saturation of the torque that originates in the horseshoe region. As shown in previous studies (Masset 2002; Masset & Casoli 2010), even when the cooling time is shorter than the u-turn time, and thus it is also much shorter than the libration time, viscosity is the driving factor for the saturation of the horseshoe region. This yields the following type I migration coefficient C_T1 (in the second term of Eq. (3)): ${\mathit{C}}_{\mathit{T}\mathrm{1}}\mathrm{=}{\mathit{C}}_{\mathrm{Lind}}\mathrm{+}\begin{array}{c}\min \end{array}\mathrm{(}\mathrm{1}\mathit{/}{\mathit{s}}_{\mathrm{2}}\mathit{,}{\mathrm{1}}^{\mathrm{)}}\left(\mathit{z}\mathrm{\left(}{\mathit{s}}_{\mathrm{1}}\mathrm{\right)}{\mathit{C}}_{\mathrm{HS}\mathit{,}\mathrm{loc}}\mathrm{+}\left(\mathrm{1}\mathrm{-}\mathit{z}\mathrm{\left(}{\mathit{s}}_{\mathrm{1}}\mathrm{\right)}\right){\mathit{C}}_{\mathrm{HS}\mathit{,}\mathrm{adia}}\right)\mathit{.}$(20)This assumes that even when cooling is much faster than the libration time scale the horseshoe drag will saturate without sufficient viscous coupling of the horseshoe region to the rest of the disk. The horseshoe part itself shifts between the adiabatic and the locally isothermal value depending on the ratio s₁ of the cooling time scale to the u-turn time scale (if h is constant also with respect to the libration time scale (Eq. (18)).

The reduction of the surface density at the planet’s location due to the beginning gap formation for more massive planets will also modify the migration behavior. Crida & Morbidelli (2007) derived a formula to estimate the depth of the gap relative to the unperturbed disk. In their definition, a gap is formed when the surface density is reduced to 10% of the unperturbed value. They defined a factor (their Eq. (12)) ${\mathit{P}}_{\mathrm{\Sigma }}\mathrm{=}\frac{\mathrm{3}\mathit{h}}{\mathrm{4}\mathrm{3}\sqrt{\mathit{q}\mathit{/}\mathrm{3}}}\mathrm{+}\frac{\mathrm{50}\mathit{\nu }}{\mathit{q}{\mathit{a}}_{\mathrm{p}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{p}}}$(21)as a combination of two criteria (the thermal and the viscous condition) and used this factor to approximate the depth of the gap as ${\mathit{f}}_{\mathrm{\Sigma }}\mathrm{\left(}{\mathit{P}}_{\mathrm{\Sigma }}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}{\mathit{P}}_{\mathrm{\Sigma }}\mathit{<}\mathrm{2.4646}\mathit{,}\\ & \mathrm{otherwise}\mathrm{.}\end{array}$(22)We multiplied Σ_p in Eq. (3) with this factor to reduce the surface density in our torque calculations. This means that the type I migration rate is reduced by a factor of up to 10 at the point of transition to type II migration. This reduced type I migration rate is still about one order of magnitude higher than the type II migration rate for typical disk conditions (Sect. 2.5), which necessitates using a smooth transition function necessary.

Eventually, when the planet mass reaches the gap opening mass M_gap, that is, the mass at which f_Σ = 0.1, the torque transitions to type II: $\mathit{C}\mathrm{=}\mathit{z}\left(\frac{{\mathit{M}}_{\mathrm{p}}}{{\mathit{M}}_{\mathrm{gap}}}\right){\mathit{C}}_{\mathit{T}\mathrm{1}}\mathrm{+}\left(\mathrm{1}\mathrm{-}\mathit{z}\left(\frac{{\mathit{M}}_{\mathrm{p}}}{{\mathit{M}}_{\mathrm{gap}}}\right)\right){\mathit{C}}_{\mathit{T}\mathrm{2}}\mathrm{·}$(23)We used the same functional dependency as in (Eq. (19)) for a smooth transition from the type I to the type II torque. Here, we set the smoothing parameter b = 10.0 (fast transition), because the reduction of the surface density already partly smoothes the transition (see Sect. B). The transition factor s_{type II} is in this case the ratio of the planet mass to the mass at which the planet would open a gap in the disk: ${\mathit{s}}_{\mathrm{type}\hspace{0.17em}\mathrm{II}}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{p}}}{{\mathit{M}}_{\mathrm{gap}}}\mathrm{·}$(24)The factor C_T2 is the corresponding type II migration coefficient in the same dimensions as C_T1, given as ${\mathit{C}}_{\mathit{T}\mathrm{2}}\mathrm{=}\frac{\mathit{ȧ}\mathrm{p}\mathit{,}\mathrm{T}\mathrm{2}{\mathit{M}}_{\mathrm{\ast }}{\mathit{h}}^{\mathrm{2}}}{\mathrm{2}{\mathit{a}}_{\mathrm{p}}^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{p}}\mathit{q}{\mathrm{\Omega }}_{\mathrm{p}}}\mathrm{·}$(25)

2.5. Comparison with other models

In this section we compare our model with that of Paardekooper et al. (2011) for a specific choice of parameters. The global consequences in a full planetary population are shown in Sect. 4.1.2. Here, we use a simple, nonevolving 1D-disk and compare the torques predicted by the two models in such a disk with results from the 3D radiation-hydrodynamic simulations of Kley et al. (2009) and Bitsch & Kley (2011).

Paardekooper et al. (2011) developed a more complex migration model trying to use first principles to derive in particular the transition functions compared with our model described in Sect. 2. They derived their torques by calculating the effect of viscosity and thermal diffusion onto the horseshoe region itself. From this they found transition functions between different parts of the corotation torque, using thermal diffusion time scales and viscous time scales as transition criteria between barotropic and entropy-related parts of the horseshoe drag and the linear corotation torque and also the saturation of both. On the other hand, we altered our torque according to the ratio of these time scales using an arbitrary transition. They called their ratio of these time scales p_ν and p_χ with ${\mathit{p}}_{\mathit{\nu }}^{\mathrm{2}}$, the ratio scaling the saturation due to viscosity, being equal to our s₂ ratio if f_visc = 0.75. For p_χ one can show that $\frac{\mathit{\pi }}{\mathrm{3}}{\mathit{p}}_{\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\frac{{\mathrm{1.16}}^{\mathrm{2}}\mathrm{8}}{\mathrm{9}}\frac{\mathit{h}}{\sqrt{\mathit{\gamma }}}{\mathit{s}}_{\mathrm{1}}\mathrm{·}$(26)The numerical factors in front of both sides of the equation are almost equal 1. The transition functions in Paardekooper et al. (2011) and our work are different. Therefore we considered the turnover points to compare the transitional behavior. In our model it occurs when s₁ = 1. For typical values of h and γ this leads to a value of p_χ ≈ 0.2. While not as easy to characterize due to the multiplication of several factors (see their Eqs. (23), (30), (31), and (53)), this is in the range where the transition in the entropy torques occurs also in the model of Paardekooper et al. (2011). The main difference in the description is the separate treatment of the “entropy” and “barotropic” parts in Paardekooper et al. (2011) with three slightly different transition functions between different parts of the torque. We only used one transition function between the locally isothermal and adiabatic regime to attribute the dependence of the torque on cooling processes. Finally, they used only one free parameter, the smoothing factor for the planetary potential, which we set to 0.4 in our comparison here. They also stated that their model is mostly applicable to planets with a mass of a few Earth masses, while our model is intended to be applicable to a wider mass range.

To quantify the impact of the differences, we next compared the torques predicted by our simple model with those in 3D radiation hydrodynamic simulations, and also with the results found with the Paardekooper model. As will become clear, one finds that despite the simplifications in our model, we can fit the numerical results relatively well.

The 1D disk for the comparison was set up to be as similar as possible to the disk from the 3D simulations of Kley et al. (2009) and Bitsch & Kley (2011). The surface density Σ is a power-law in semimajor axis with a fixed slope of −0.5, and the temperature T is a power-law in semimajor axis of −1.5 inward of 2.5 a_Jup, where a_Jup is the semimajor axis of Jupiter. Farther out, when the temperature almost reaches 10 K, its power-law exponent β_T increases and approaches 0. We used the same EOS (Saumon et al. 1995) and opacity tables (Bell & Lin 1994) as in our full model. The viscosity was set to the same value as in Kley et al. (2009). All other disk parameters, such as disk height H or density ρ, were calculated from Σ and T. In this disk we then calculated the specific torques on planets of either different mass or different semimajor axes. With β_Σ = − 0.5 and β_T = − 1.5 there will be no type II torque since there is no radial gas flow in this disk. We stress that our disk model for the planet population synthesis calculations is in contrast a time-evolving 1+1D α model and is summarized in Sect. 3.1.

In Fig. 2 we plot torques from our model as a function of planetary mass and compare them with simulations from Kley et al. (2009), while in Fig. 3, we compare our model with Bitsch & Kley (2011), where torques at different semimajor axes for a 20 M_⊕ planet were calculated. In both figures we also show the torque as predicted from the model of Paardekooper et al. (2011; for their suggested optimal free parameter b = 0.4).

We plot the migration rates for our model with and without the reduction of the surface density due to gap formation (solid and dashed lines respectively). This reduction was not only applied to our own model but also to the Paardekooper et al. (2011) prediction. For the dependency of planet masses (Fig. 2), the curves with the gap reduction remain closer to the data obtained by Kley et al. (2009). At around 200 M_⊕ when the transition to type II occurs, this factor causes a 10 times lower type I torque than without the reduction. Therefore with the zero type II torque in this disk-setup, the overall torque is also 10 times smaller. The red lines correspond to our model described above (with f_visc = 1: hereafter the STD “standard”-model), while the green lines are our model with f_visc = 0.55 (hereafter the BMF “best mass fit”-model). The latter value for f_visc was chosen to increase the saturation mass in a way so that the mass of zero torque in our own semi-analytical model agrees with the data of Kley et al. (2009). The associated reduction of transition parameter s₂ could be interpreted as a more efficient viscous injection of angular momentum into the horseshoe region from the rest of the disk than in the STD-model. Only a factor ~2 is needed to bring our model in agreement with the hydrodynamic simulations, meaning that the simple estimate of the transition using the characteristic time scales leads to relatively accurate results. Setting f_visc = 0.55 also increases the maximum specific torque, which follows the data from Kley et al. (2009) more closely than our model with f_visc = 1, or the model of Paardekooper et al. (2011), at masses larger than 15 M_⊕.

In Fig. 3 we compare torques on a planet of a fixed mass (20 M_⊕) as a function of semimajor axis. Here we see as well that the BMF-model provides the best agreement with the data of Bitsch & Kley (2011). The semimajor axis at which the torque vanishes is with 2.7 a_Jup relatively close to the 2.5 a_Jup found in the hydrodynamical simulations. The model of Paardekooper et al. (2011) places the position of zero torque at 1.8 a_Jup, while in the STD-model it occurs at 1.4 a_Jup. Thus, for planets closer than 2.5 a_Jup the BMF-model seems to be the best representation of both sets of numerical simulations.

At larger distance from the star the BMF-model and Paardekooper et al. (2011) again yield similar results and compare relatively well to the data. All curves also indicate that including the reduction of the surface density in Eq. (3) due to gap formation leads to a better agreement with the torques found in the hydrodynamic simulations, especially in the outer parts of the disk. In summary, we see that by setting f_visc = 0.55, the model agrees fairly well with the data provided from 3D simulations.

Some diskrepancy exists for the innermost point considered in the hydrodynamic simulation. It is most likely caused by effects of the change of the opacity due to ice evaporation, which is handled in a different way in our model vs. the 3D full models, so that the gradient of temperature and surface density differ.

Note that this comparison in principle only applies for a fixed value of the Prandtl number. The latter depends on γ and the optical depth τ (Paardekooper et al. 2011), which varies, for instance with the distance from the star or temperature. In the limit of high optical depth, the Prandtl number approaches $\mathit{Pr}\mathrm{=}\frac{\mathrm{9}}{\mathrm{4}}\mathit{\gamma }\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}$.

Nevertheless, we conclude that our BMF-model can quite adequately reproduce the current understanding of planet migration based on (semi-) analytical models for the modeling purpose of planetary populations.

2.6. Model versions

In the last section, we have calibrated our semi-analytical model with the results of a specific set of radiation-hydrodynamic simulations. With these results, we define three versions of the model that are used below to study the (global) effects of these different variants of our migration model:

• STD model. In the standard version we simply setf_visc = 1.0 so that the plain time scales as defined in Sect. 2.3 are employed to calculate the point of saturation.

• BMF model. In the best mass fit version we multiply s₂ by 0.55 (f_visc = 0.55). This model compares best with the 3D migration simulations as shown above.

• RED model. In the reduced version we reduce s₂ by multiplying it with 0.125 (f_visc = 0.125). This results in a four times higher saturation mass than for the STD case. We use this even larger reduction to study the effect of different reduction factors.

3. Formation models and migration tracks

In addition to the comparison shown above for the migration rates alone, we are interested in the global effects of the new migration model onto a planet population compared with our older isothermal migration model. The planet formation model used for this is based on the paradigm of core accretion (Perri & Cameron 1974; Mizuno et al. 1978; Pollack et al. 1996; Alibert et al. 2005). The model is described in detail in Alibert et al. (2005) with later modifications shown in Mordasini et al. (2009a), Mordasini et al. (2010, 2012a,b), and Fouchet et al. (2012). It consists of different computational modules, namely the planet accretion module, the disk module, and the migration module. The migration module is described above and was used in recently published work of Fortier et al. (2013), Mordasini et al. (2012a, b). In the following sections we give a short overview of the disk and accretion modules. Then we present calculations of a few specific planets, show the associated formation tracks and the important features found in them. Finally, we present results from planet population synthesis calculations. Since we concentrate in this paper on the direct effects of the new migration model, we present here only simulations with just one embryo per protoplanetary disk. The interplay of migration and multiple concurrently forming planets (Alibert et al. 2013) will be addressed in future work.

3.1. Disk model

As a model for the protoplanetary disk, we used a 1+1D model as in Papaloizou & Terquem (1999) or Bell & Lin (1994), and present only a short summery here. The protoplanetary disk is described as a viscously evolving α disk, where α is assumed to be constant throughout the disk. The viscosity is given as ν = αHc_s with H the disk scale height and c_s the sound speed (Shakura & Sunyaev 1973). Irradiation effects from the host star can be included in the calculations of the vertical structure (Fouchet et al. 2012). For the evolution of the gas surface density Σ over time t and distance a from the star we solved the standard viscous evolution equation from Lynden-Bell & Pringle (1974). For most of our simulations in this paper we neglected stellar irradiation and used only viscous heating if not otherwise mentioned: $\frac{\mathrm{d\Sigma }}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{3}}{\mathit{a}}\frac{\mathit{\partial }}{\mathit{\partial a}}\left[{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{\mathit{\partial }}{\mathit{\partial a}}\mathrm{\left(}\mathit{\nu }\mathrm{\Sigma }{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\right)}\right]\mathrm{+}\mathrm{\Sigma ̇}\mathrm{w}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{.}$(27)As a sink term we included the photoevaporation rate as given in Mordasini et al. (2009a). Together with α and the initial disk mass it determines the disk lifetime. At the start, the solid surface density of the planetesimals is equal to the gas density times the dust-to-gas ratio f_{D / G}. It is further reduced inside of the iceline at a temperature T = 170 K by a factor of 4. Other than by accretion onto and ejection by the planet, we did not evolve the solid surface density.

3.2. Model of accretion and internal structure

We simulated the growth of one planet per disk. For this, we inserted a 0.6 M_⊕ seed embryo at a random position in the disk. The core has a constant density of 3.2 g / cm² and contains all the heavy matter the planet accretes, that is, we assumed that all planetesimals sink to the core (Pollack et al. 1996). The seed will initially start to accrete mostly planetesimals, which leads to a growth of the core. The amount of energy released from the infalling planetesimals is high at this point and the core mass is low, therefore initially only small amounts of gas are bound. As the core grows, it binds an increasingly massive envelope. The accretion rate of gas is found by solving the standard equations of the structure of planet interiors, but with the simplification that the luminosity is constant throughout the envelope.

In the original model of Mordasini et al. (2009a), the luminosity of the planet is due to the accretion of planetesimals alone. This is usually the dominant source of luminosity at smaller masses (Pollack et al. 1996). Here, we adopted a simplified version of the model described in Mordasini et al. (2012b) to take also into account the luminosity generated by the accretion of gas. In the original model it was sufficient to consider only the luminosity due to the planetesimals, because the migration was always directed inward. This means that the cores always migrate into regions with a high solid surface density. With the new migration model, this is no longer the case: because the positive torques act at certain masses (see Sect. 2.5), it is possible that a core migrates through parts of the disk with a very low solid surface density content. There the luminosity of gas accretion becomes important.

3.3. Convergence zones

These positive torques lead to regions in a disk where a planet within the adiabatic migration regime migrates outward instead of inward (Lyra et al. 2010; Mordasini et al. 2011). Figure 4 shows the values of C_adia (which represents the strength of the migration in the adiabatic regime, see Eq. (6)) at time equal 0.02, 0.1 and 0.5 Myr in a nonirradiated evolving α-disk with Σ₀ = 420 g / cm² and Ṁ_w = 6.7 × 10^-9M_⊙/ yr. These parameters results in a lifetime of the disk of 2.8 Myr.

In Fig. 4 one can see two regions of outward migration, i.e., with positive values of C_adia, from 0.46 AU to 2.6 AU and farther out from 3.6 AU to 11.5 AU at 0.02 Myr (red curve). At later times these regions are closer to the star as the disk evolves. To further illustrate the time evolution of the regions in Fig. 5 we show the direction of migration as a function of time for this disk. Blue indicates regions of inward migration while green shows outward migration. One can see that these regions slowly move inward as the disk evolves. This occurs on the viscous time scale of the disk (Lyra et al. 2010). Owing to the existence of the outward and inward regimes, there are special locations in the disk. For example, at an age of 0.5 Myr, the migration changes from outward to inward at 10 AU (and at about 1.2 AU) when moving to a larger distance (i.e., these are points of zero torque where the derivative of the torque with distance is additionally negative). Such a point is called a convergence zone (CZ). It is called this because planets in its vicinity converge on this zone from both inside and outside (see Lyra et al. 2010). The domain in orbital distances from which planets migrate to the convergence zone is the associated convergence region. At 0.5 Myr, the inner convergence region extends from 0.2 to 1.8 AU for example. Similar results for two convergence zones for various disk models were also recently presented by Kretke & Lin (2012) and Yamada & Inaba (2012). After reaching the convergence zone, a planet remains slightly outward, but close to it, so that the net torque pushes the planet inward at the same migration rate as the zone itself.

Once captured in a convergence zone, the planet migrates on a time scale that is at least an order of magnitude longer than typical type I migration time scales. For example, while captured in the inner convergence zone, the planet discussed below has an equivalent migration coefficient C (Eq. (3)) of 0.034 while C is during normal type I migration on the order of 1 as shown for instance in Fig. 4 or Paardekooper et al. (2010).

However, for some conditions, a low-mass planet cannot migrate at a sufficiently high rate to remain close to the CZ. Instead, the planet leaves the CZ and falls behind it (the planet still migrates inward, but is overtaken by the CZ). This occurs especially for the inner convergence zone. Typically, after leaving the inner zone, the planet transitions into the outer one, where it is again captured. The reason is that the type I migration rate is proportional to the planet mass. For a sufficiently low mass of the planet, the type I migration time scale thus becomes longer than the viscous time scale of the disk, which sets the speed at which the CZ moves. This characteristically occurs at the end of the disk lifetime, when the gas surface density is low, so that the type I migration time scale becomes even longer.

The position of the inner CZ is close to the distance of the local minimum in the opacity at a temperature of ≈200 K (Lyra et al. 2010), which is the temperature where ice grains are completely molten in the opacity law of Bell & Lin (1994). The change in the slope of the opacity at this point leads to a change in the temperature power-law exponent β_T, which leads to the change of the sign of C_adia and therefore to a convergence zone. The reason for the outer CZ is the change of the temperature slope due to the convergence onto the background temperature in the outer part of the disk.

The convergence zones only exist for a certain range of planet masses (Kretke & Lin 2012). A low-mass planet will migrate in the locally isothermal regime since the thermal processes in the disk are fast enough to regulate the temperature during a u-turn of the gas. On the other hand, when a planets becomes massive enough, the angular momentum flux through the horseshoe region will be too low to support the unsaturated horseshoe drag, so this part of the torque saturates and in general rapid inward migration sets in.

3.4. Migration and formation tracks

The new migration model leads to changes in the evolution of a protoplanet that are different from those described in Mordasini et al. (2009a). We discuss the new behavior in this section where we show the migration and formation tracks of a typical protoplanet of our synthesis simulation. Not every evolving protoplanet in our calculations shows every behavior described below, but most evolutionary tracks of planets feature at least some part of the behavior. The consequences of the new migration model for a whole synthetic population are studied after this.

Further to what was discussed before, Fig. 5 also shows two tracks of an evolving protoplanet starting at 0.57 AU at 2000 yr after the start of disk evolution. Figure 6 shows the corresponding formation tracks of this protoplanet in the semimajor axis mass plane. In both diagrams, the simulation represented by the solid line uses the BMF migration model, where f_visc = 0.55, while the dashed line is used for the RED model with f_visc = 0.125 (see Sect. 2.6), simulating a larger saturation mass. The beginning of the evolution of the protoplanet is the same in both models. It starts to migrate inward in the locally isothermal migration regime (blue part of the track in Fig. 6) and accretes solids from its surrounding, depleting this part of the disk of planetesimals. When it reaches a mass of approximately 2.5 M_⊕, its horseshoe region is broad enough that cooling cannot keep the gas that is on horseshoe orbits in a locally isothermal state. The planet enters the unsaturated adiabatic regime, shown in red in Fig. 6, and starts to migrate outward (at the position marked by a circle in the figures). This occurs at 0.1 Myr as seen in the migration track in Fig. 5. While migrating outward, the planet does not grow at first because there is almost no solid material left because of its previous passage through this part of the disk. The protoplanet is also still too small to accrete significant amounts of gas. After crossing its initial position, growth by accretion of planetesimals sets back in and the planet migrates outward until it enters the convergence zone at 0.6 Myr (diamond symbol). Here the direction of migration again changes to inward as the planet is now bound to the CZ. While the disk evolves, this zone and the captured planet move slowly inward. The planet again moves through a region that is depleted of almost all planetesimals. But, with a 10 M_⊕ core, the planet is now massive enough to bind nebular gas in its envelope. This is especially true when there is not much solid accretion, which means that the core luminosity is low, and therefore the gas accretion rate is high. The planet grows by accreting gas, until at about 1.25 Myr in the BMF model (solid line), it reaches the mass where saturation sets in (square symbol). The positive horseshoe drag is reduced with increasing mass and is not sufficient to counterbalance the negative Lindblad torques. The planets therefore leaves the CZ and rapidly migrates inward (saturated adiabatic migration regime, magenta in Fig. 6).

The planet continues to accrete gas until the inner edge of its feeding zone reaches the distance of the previous closest approach to the star (upward-facing triangular symbol). Again in reach of planetesimals, solid accretion restarts and increases the core luminosity. The envelope expands, heated by this process, and because it is still connected to the disk at this time (attached phase, see Mordasini et al. 2012b), the heating leads to the loss of envelope mass and a corresponding reduction of the total planetary mass (see Sect. 3.2). This reduces the level of saturation, which in return increases the horseshoe drag, which can again balance the Lindblad torque resulting in an almost complete stop of migration: if the planet were to migrate outward again, solid accretion would stop and rapid gas accretion would start to increase the mass. This would reduce the horseshoe drag and push the planet back inward. In the opposite case of further inward migration, solid accretion would become stronger. This would increase the core luminosity and remove more envelope mass, and the total mass decreases. This would increase the outward-directed horseshoe drag because of the reduced saturation of the corotation torques. The planet would thus tend to migrate outward. The combination of these two points means that the planet has reached a quasi-stable point due to the interplay of accretion and migration.

With the evolution of the disk, the saturation mass at a given point in the disk decreases over time. Thus the mass where the partly saturated horseshoe drag balances the Lindblad torques is also reduced with time. Because of the interaction of accretion and migration, the planet’s mass remains exactly at this zero-torque mass.

The planet is forced to loose more mass and therefore remains at this semimajor axis for a given moment in time, while, as a result, the planet moves in slowly over longer time scales. It just “nibbles” on the edge of the planetesimal disk that was depleted up to this distance, while reducing its mass. The gas envelope mass is decreased until a level is reached that the planet can support with the full solid accretion rate of normal inward migration.

In case of the BMF model, this is a slow process. The protoplanet remains for 0.5 Myr at 0.5 AU while it looses 6 M_⊕ of the envelope it has accreted while moving inward in the CZ and later in the saturated adiabatic migration regime. At this time it has only 0.1 M_⊕ of envelope left. This is the same amount as at the time it became bound to the CZ. The RED model and its increase of the saturation mass leads to a slightly different behavior of the protoplanet. It remains for a longer time in the CZ and can thus accrete gas for a longer time than the BMF model protoplanet. This results in a larger envelope mass when it reaches the distance of the previous closest approach (left-facing triangle). But since in both cases the protoplanet migrated through the same part of the disk, both have the same amount of solids accreted and thus the same core-mass. Therefore, the increase of the solid accretion rate also leads to a mass loss of the envelope in the RED model case. But here, the planet is still in the unsaturated migration regime while loosing most of its envelope mass. Therefore, the mass loss does not lead to a change in the migration rate. The protoplanet is still bound to the CZ and is pushed by it into the remaining planetesimal disk. Therefore, the mass loss occurs much faster here, the planets looses 11 M_⊕ in only 0.1 Myr.

This illustrates the strong inter-dependence of migration and accretion. The solid accretion rate is set by the amount of solids reachable by the planet, and therefore by the migration that brings the protoplanet into new regions of the disk. But this behavior is only true if the availability of planetesimals themselves, and not the collision rate, is the limiting factor for the solid accretion rate Ṁ_core. In the planet-envelope structure model used here, the solid accretion rate, the associated core luminosity, and the mass of the core itself define the envelope structure and thus the envelope mass of a planet.

The small loops in the track that occur while the planet in the BMF model looses mass (between the up- and down-facing triangles) are caused by a finite time-step. Some larger loops are visible in the formation tracks in Fig. 7. We show here the data obtained during the population synthesis calculations discussed in Sect. 3.5 and shown in that figure. We also separately simulated the same case with a much shorter time-step and obtained the same results without these small loops.

We did not consider an increase in the random velocities of the planetesimals due to the presence of a planet. They remain small, as reported in Pollack et al. (1996). A more realistic oligarchic growth model as described in Fortier et al. (2007) and Fortier et al. (2013) leads to higher random velocities and thus a smaller focusing factor in the calculation of the solid accretion rate Ṁ_core. This would lead to less envelope loss if Ṁ_core remains small enough.

Another point to consider, in addition to the planetesimal accretion rate, is the treatment of the protoplanetary disk in our model. We did not let the disk evolve except for depletion due to accretion onto the protoplanet. Therefore there is a sharp edge into which the planet can migrate. A more realistic treatment (e.g., planetesimal drifting or diffusion, scattering) would lead to a gradual increase of Ṁ_core and thus to a slower loss of the envelope. But we expect the final outcome to be similar, only the track in the a − M plane would be smoother (mass loss setting in earlier and more gradually).

In both cases shown here, the protoplanet starts again to accrete solids and nebula gas and migrates inward after it looses almost all its envelope mass (down- respectively right-facing triangles in Figs. 5 and 6). In case of the BMF model it migrates in the faster saturated adiabatic migration regime until in reaches 0.1 AU and a final mass of 17 M_⊕ (0.9 M_⊕ in its envelope) and the simulation stops. In the RED model the simulation ends with the disappearance of the disk at 2.8 Myr and the planet migrated to 0.43 AU with a mass of 10.5 M_⊕ (only 0.1 M_⊕ in the envelope).

3.5. Reference population synthesis calculation

After studying the single case, we now look at a population synthesis calculation with 10 000 different initial conditions that we used as our reference when we investigated the effects of different migration models on a synthetic planet population. The important parameters of the synthesis can be found in Table 1. We used the BMF-migration model and an α parameter of 7 × 10^-3 for the nonirradiated disk.

Figure 7 shows the tracks of 250 cases, the different migration regimes color coded. The meaning of the colors is described in the caption of Fig. 7. One can see that most planets start in the locally isothermal migration regime, changing into the unsaturated adiabatic migration regime before the horseshoe drag saturates and they migrate inward. Some are large enough to transition into type II migration while others end up at 0.1 AU. Two groups corresponding to the inner (inside of 1AU) and outer (outside of 1 AU) convergence zone are visible in Fig. 7.

When a planet migrates outward through the iceline, the migration rate does not change, but its solid accretion rate does change by a factor of 4, because the increase in the solid surface density. This results in a bend and a much steeper slope in some of the tracks in the outer zone. Horizontal parts in the formation tracks that indicate migration without growth are also lacking in the tracks of the outer planets. In the outer parts, especially outside of the iceline, the amount of solid matter is too large to be completely accreted onto the planet at the given accretion rates Ṁ_core. Thus, planets can also collect material on their second pass through a part of a disk and grow.

The inner and outer CZ are the reason for three groups of planets in the final position a − M diagram. They can be seen in Fig. 8. A fourth possible cluster is not visible because the BMF-migration model leads to a situation where all those planets migrate inward of 0.1 AU (see Sect. 4.2.2):

• The first cluster, which is the one with the highest number ofplanets, lies approximately between 0.15 and0.8 AU and 1 and 20 M_⊕ (red box in Fig. 8). It consists of planets captured by the inner convergence zone. Most of the planets are directly attached to the CZ when the disk disappears and the simulation stops. Their mass is too small for a transition into saturated migration and departure of the CZ. A strip of planets extends inward from this cluster. These planets saturated shortly before the disk ended. However, most planets that transition into the saturated migration regime while being in the inner CZ eventually end up at 0.1 AU.

• The second cluster lies farther away from the star (1 AU to 4 AU) but is also in the mass range between 1 and 20 M_⊕ (blue box in Fig. 8). It consists of planets captured by the outer second CZ. It is less populated because fewer embryos start at the larger distances because the distribution of the starting positions, which are uniform in log(a). Additionally, the amount of solid material to grow is larger, thus many planets at these distances can reach masses above the saturation mass.

• The third cluster consists of planets that saturated while being in the outer convergence region (green box in Fig. 8). But they are massive enough and can accrete enough solids on their way in that most of them can transition into type II migration. Here they migrate on the time scale of the disk evolution or slower and can accrete gas until the disk ends. They reach masses between 100 and a few 10 000 M_⊕. The planets can reach such high masses because we neglected the effect of gap formation on the gas accretion rate. If the reduction of the gas accretion rate due to gap formation were included, the planet masses would be restricted to lower masses, depending on the disk viscosity and mass (cf. Bodenheimer et al. 2013).

Planets in the inner convergence region will migrate to the inner CZ and thus migrate through a large part of the inner part of the disk. But they remain completely inward of the iceline since the inner convergence region ends there due to opacity transitions. Therefore they are able to accrete most matter in the first few 0.1 AU and thus finally have at least a few Earth masses. In our model we obtain a small planet of several M_⊕ or less only when its starting time is in the last few 0.1 Myr of the disk lifetime. With more than one core per disk one planet alone cannot accrete all the matter in the inner part. The solids will be distributed in many small planets. We therefore overestimate the number of planets between 2 and 30 M_⊕ and underestimate the number of planets smaller than 2 M_⊕.

The synthesis with the new migration model also shows the desert of planets between 30 and 200 M_⊕. This is a feature of the runaway gas accretion that occurrs in the core-accretion model (Pollack et al. 1996). On the other hand, the region of close-in, low-mass planets, which was empty in Mordasini et al. (2009b), is now well populated with the new nonisothermal migration model.

The spread of the first cluster originates in the spread of the initial conditions of the photoevaporation rate . In the implementation of our disk module the photoevaporation rate determines the mass of the disk at the end of a simulation and therefore the position of the CZ at the end of a simulation. A single value of in all simulations of a synthesis would result in only one position of the CZ and therefore a high concentration of small mass planets on one radius (see also Sect. 3.3).

We used some basic statistics, namely the number of “hot” and “cold” planets and the number of massive and small planets (cf. below), to compare synthesis calculations with different migration models. While none of these numbers are compared with those of the observed population of extra-solar planets, the difference between the calculations allows us to see the importance of different parameters or parts of the migration model. Comparison with the observed population will be made in future work when other such as effects like the decrease of the disk mass due to accretion onto the planet or multiple concurrently forming planets (Mordasini et al. 2012b,a; Alibert et al. 2013) are also included. Out of 10 000 initial conditions we obtained 6850 planets more massive than our starting mass of 0.6 M_⊕ in the synthesis calculation described above. The remaining 3150 initial conditions have starting times (time when we insert the planet embryo into the disk) longer than the lifetime of the corresponding disk. Out of these 6850 planets 54,4% migrated to 0.1 AU, the inner border of the computational disk, and are called “hot” planets. The other 45,6% ended further out and will be called “cold” planet.

Finally, there are 912 (13.3%) massive planets in total with M> 100 M_⊕. While the majority of all the planets in the synthesis are “hot”, the massive planets split into 276 “hot” and 636 “cold” planets. Thus there are more “cold” giant planets than “hot” ones.

4. Impact of different migration prescriptions

In this section we compare results from population synthesis calculations where we changed one aspect of the migration model relative to the calculation above (Sect. 3.5), but otherwise used the same initial conditions and settings. We therefore refer to the synthesis above as the reference synthesis. We first compare it with migration models described in earlier studies before we change some physics of the model itself.

4.1. Earlier prescriptions

We use two earlier prescriptions, the isothermal migration model of Tanaka et al. (2002) used in our earlier work and the migration prescription from Paardekooper et al. (2011).

4.1.1. Isothermal migration model of Tanaka et al. (2002)

The first comparison was made with calculations preformed with the nonreduced isothermal migration model (f₁ = 1) based on the results of Tanaka et al. (2002). While we used the same type I migration model as in Mordasini et al. (2009a), the results here still differ from those published in Mordasini et al. (2009a), since we here use the same transition criteria for the transition from type I into type II migration as in the reference synthesis. This is different from the one in Mordasini et al. (2009a). There only the thermal condition (H>R_H) for the transition into type II was used. This leads to much smaller transition masses than here.

Figure 9 shows the final position of the planets in the distance-mass diagram in red, right-facing triangles. Blue, left-facing triangles show the reference synthesis. While the range in mass and distance covered by the planets is the same, there is no clustering caused by the migration into a CZ at small masses (0.6 to 30 M_⊕). The total number of “cold” planets (a> 0.1 AU) is only ≈45% of the number found in the reference calculation. From the 6850 synthetic planets more massive than 0.6 M_⊕ we find 79.2% “hot” planets with 2.7% “hot” massive planets and only 20.8% planets outside of 0.1 AU. Of these, 346 are more massive than 100 M_⊕ (5%). The new nominal nonisothermal migration model in comparison gives twice as many planets that remain outside of 0.1 AU and also almost twice as many massive planets. Some preliminary calculations indicate that this ratio can increase even more with lower values of α. The nominal migration model still leads to a loss of more than half of the planets with more than 0.6 M_⊕ into the inner part of the disk inside of 0.1 AU and therefore potentially into the star. On the other hand, the new model doubles the number of planets outside of 0.1 AU compared with the isothermal migration model without any arbitrary reduction factor.

4.1.2. Paardekooper et al. (2011) migration model

Paardekooper et al. (2011) developed a migration model that is similar to the one described in Sect. 2, but more sophisticated, because it uses thermal-diffusion time scales and viscous time scales as transition criteria between barotropic and entropy-related parts of the horseshoe drag and the saturation of both. The difference between their and our model are discussed in Sect. 2.5, where we compared torque curves from our model with those of this model. We also made two synthesis calculations using this migration model. The final a − M distribution of the two simulations is shown in Fig. 10. Green dots represent the final positions of the synthetic planets obtained with the Paardekooper et al. (2011) model and a fixed γ = 1.4. The second calculation is shown in red, here γ was determined with the EOS we usually use (Saumon et al. 1995). The blue symbols show the reference synthesis.

For planet masses larger than three Earth masses there is no real difference between the three plotted syntheses. The simulation with the free γ produces nearly the same number of “cold” planets, the simulation with the fixed γ slightly (5 − 10%) less relative to the nominal BMF model. The situation is slightly different for the number of massive planets with M_p> 100 M_⊕ outside of 1 AU: the free γ synthesis has only 85% of the number of “cold” massive planets of the BMF-model synthesis and the fixed-gamma simulation about 60% of the nominal model.

The reason for the small differences in the overall amount of “cold” planets but the larger differences for the “cold” massive planets is, as visible in Figs. 2 and 3, that overall both models produce the same general migration behavior: first planets migrate inward, then outward to a convergence zone, and after saturation inwards again. But the point of crossover from outward to inward migration is closer in and at lower masses for the Paardekooper et al. (2011) model relative to the nominal BMF model. The planets saturate at lower masses and therefore start to rapidly migrate inward earlier in their evolution and fewer planets can reach the type II migration regime and become massive.

Overall, there are no large difference between the different migration models. On the other hand, our simpler BMF model seems to agree better with the torques obtained with the radiative hydrodynamical simulations of Kley et al. (2009) and Bitsch & Kley (2011).

4.2. Different input physics

We now study the effect of different Lindblad torques and of different saturation masses onto the synthetic planet population.

4.2.1. Lindblad torques

As stated in Sect. 2.2, there are two different formulas for the Lindblad torques. We changed the Lindblad torque to the one used in Masset & Casoli (2010) for one synthesis. We show the final positions of the planets in the distance-mass diagram in Fig. 11. The color-coding is the same as before. The resulting Lindblad torque is weaker with the equation of Masset & Casoli (2010), thus the migration in the saturated regime is slower. This results in moving the high-mass, third clusters farther out. In the reference synthesis the outer cluster is located between 0.1 AU and 0.6 AU. Here, it is located between 0.2 AU and 1 AU (red solid box). Inside of 0.2 AU one can see some planets of a fourth cluster. The planets in this inner group originate in part from fast growing planets, that is, those large dust-to-gas-ratios, of the inner convergence zone. The rest are planets from the outer zone, which saturated but grow rapidly enough during their fast inward-migration phase to transition into type II migration. In contrast to the reference synthesis, where all planets of this group migrated to the inner border of the computational domain at 0.1 AU, the planets are now at larger distances and can be seen in the calculation (red dashed box). The third cluster (outside of 1AU) completely consists of planets in the outer convergence zone.

The smaller torque also results in a smaller and weaker region of inward migration between the CZ. It is small and weak enough that in some cases, planets can drift from the inner into the outer convergence region, because the migration rate is lower than the movement rate of the CZ due to the evolution of the disk (see Sect. 3.3).

For low-mass planets below the saturation mass the differences are smaller between the two syntheses. There are still the two clusters, associated with the two convergence zones.

We can thus conclude that the weaker Lindblad torque results in an shift of the final position of a planet in the distance-mass diagram to the right that is, to larger distances. Because of this we also see more “cold” (increase from 45.6% to 51.2%) and more massive planets relative to the reference synthesis (increase from 9.3% to 12.1%). The calculated migration rate in saturation is only about a factor of 2−3 smaller with the formula of Masset & Casoli (2010) than that of Paardekooper et al. (2010). This seems small compared with the previous reduction of type I migration by a factor of 10−1000. But still, this small change in the description of a part of the torque by a factor of three has observable influence on the distribution of massive planets by up to a factor of 2 in semimajor axis for some cases.

4.2.2. Saturation mass

We furthermore conducted population synthesis calculations with two different values of f_visc. One with a larger f_visc = 1 (STD model) and one with a smaller f_visc = 0.125 (RED model) than the reference synthesis (BMF model, f_visc = 0.55) (see also 2.6). This only affects planets that are massive enough to undergo the transition into a saturated migration regime. An eight times smaller f_visc will result in a four times larger saturation mass (Eq. (A.1)). Figure 12 shows the final positions of the planets in the distance-mass diagram for the two nonnominal calculations and for the reference synthesis. The blue (red) points are the RED (STD) case. For most points of the two clusters of low-mass planets there are only small differences in between these three simulations. They result from the onset of saturation at different masses. For example, some of the planets in these clusters will start to saturate with the STD model while they are still in the unsaturated adiabatic migration regime at the same mass with the RED model and thus migrate faster and also accrete differently because of the different migration rates.

The high-mass clusters are shifted farther out in the synthesis with f_visc = 0.125 than in the reference calculation. Moreover, the fourth group mentioned in Sect. 4.2.1 is visible for planets obtained from the RED case. In this situation the time planets spend in the saturated migration regime is shorter because of the higher saturation mass, therefore the distance they migrate inward is smaller and the planets end up farther out. This has an effect on the content of heavy elements in the planet. Less migration means that the planet can reach less amount of planetesimals in our disk. The difference is the amount of solids in the annulus of the planetesimal disk that is not visited by the planet with the larger saturation mass. But the difference in the final mass for most massive planets (>1000 M_⊕) is small, lower than 1%. The reason is that they grow the most while they are in type II migration and the runaway accretion phase. In this phase the accretion is dominated by gas accretion. Here the growth is set by the remaining disk lifetime, and this is almost the same for the different cases.

Note that we made one simplifying assumption: we included the eccentric instability (Kley & Dirksen 2006). Therefore gap formation does not lead to a reduction of the planetary gas accretion rate. If this effect were not included, the maximal planet masses would be on the order of 10 Jovian masses (Lubow et al. 1999; Armitage 2007).

Planets beyond 2 AU and between 60 and 400 M_⊕ are closer to the star with a higher saturation mass level than the planets mentioned above. The disk temperature is lower in the outer parts and the slope of the temperature profile begins to become flatter (β_T → 0) as the temperature approaches the assumed background temperature of 10 K. This change in the slope leads to a strong change in the horseshoe drag at a radius around 10 AU for early disk times and farther in at later times. It dictates the position of the outer convergence zone (Kretke & Lin 2012). The change of the slope causes the horseshoe drag itself to become negative not far outside of the outer CZ and pushes the planet inward, as does the Lindblad torque. This means that an increase in the saturation mass will result in faster overall inward migration for the planets that start outside the outer CZ and are saturated in this part of the disk.

When comparing the numbers of “cold” planets in the STD and the RED case we see that the increase of the saturation mass by a factor of 4 results in 64% more “cold” planets and also 156% more “cold” massive planets. The number of “hot” massive planets is reduced from 7.5% to 1.4% as the slower overall migration shifts the clusters outward. Even if we cannot directly compare our results with observations, as a reference, Mayor et al. (2011) stated an observational value of ≈1%.

5. Irradiated vs. nonirradiated disks

We recently included irradiation of the host star into the disk model assuming an equilibrium flaring angle (Fouchet et al. 2012). In all calculations presented above, viscous heating only determined the thermal structure of the disk. This is sufficient in the inner parts of the disk at the beginning of the simulations, but leads to too low temperatures in the outer parts of the disk and in the later phases of the disk evolution, when the rate of gas flow through the disk becomes low and the irradiation is dominant. This effect leads to a different temperature gradient, which is important for the strength of the torques (Lyra et al. 2010). The increased heating in the outer parts makes the temperature profile less steep throughout the disk. The temperature is still around 20 K at 50 AU and decreases with distance, while in the nonirradiated case the temperature dropped to the background temperature of 10 K at 20 − 30 AU and became nearly constant. The temperature structure in the whole disk is set only by the irradiation when the disk is nearly gone and viscous heating is unimportant for the temperature structure of the entire disk. The different profile also affects the shape of the convergence regions. The parts of inward migration in the disk become smaller with time and vanish after a few million years, but still before the disk dissolves (Kretke & Lin 2012). This affects first the outer and then the inner CZ. Thus, the CZ is no longer a stopping point for type I migration throughout the complete lifetime of a disk. Therefore the confinement of planets in a small radius (see Appendix A) even when only one value of the photoevaporation rate is used is not the case for irradiated disks.

Figure 13 shows the formation tracks of 250 initial conditions and Fig. 14 shows the final positions of the planets in the distance-mass diagram. In both figures color shows the different migration regimes a planet is in in the same ways as in Fig. 7. The different symbols in Fig. 14 indicate if a planet migrated in adiabatic regime during its evolution (empty triangles or circle) or not (filled triangles or circle). The new temperature structure and increase in temperature results in a switch into the adiabatic migration regime at a higher mass than in the nonirradiated disk calculations. Moreovre, as seen in these figures, almost all small planets end in a locally isothermal migration regime even when they went into the adiabatic regime for some part of their evolution (as shown by the red part of the tracks in Fig. 13) and the unfilled cyan triangles in Fig. 14, they end in a locally isothermal regime again when the disk is almost gone. The higher temperature than in the one in nonirradiated disks, especially at the end when almost all gas is gone and only minor viscous heating occurs, leads to a shorter cooling time and therefore to the transition from adiabatic migration into locally isothermal migration.

The irradiation also leads to higher disk accretion rates in the outer regions of the disk, since ν = αc_s/ Ω ∝ T, and therefore to the faster depletion of the outer disk regions. The surface density is therefore lower in the irradiation case than in the nonirradiated at same time of the simulation. This leads to lower migration rates even in locally isothermal migration. This compensates to some degree the larger inward migration zones in our disks and the transition into locally isothermal migration at the end of the disk lifetime.

We used the same initial conditions as in all other syntheses. The different disk model leads to different disk lifetimes, therefore the results are only partially comparable. We considered 7700 planets more massive than 0.6 M_⊕. From these planets 68.8% migrated to the inner border of our computational domain. From the planets with a final distance larger than 0.1 AU from the star ≈20% are more massive than 100 M_⊕, which is similar to the ratio in our reference synthesis. The larger amount of “hot” planets is the result from the disappearing CZs at toward end of the disk evolution.

6. Summary and conclusions

We have compiled a prescription for type I migration based on the latest hydrodynamic simulations of planet-disk interactions. We tested the influence of the prescription on the outcome of the planetary population synthesis calculations. Our migration model is based on the combination of results from Paardekooper et al. (2010) and different time scales to distinguish between different migration regimes. These time scales reflect the thermodynamical behavior of the interaction between the planet and the surrounding disk like the horseshoe drag. We first compared the migration torques of this model with a model of Paardekooper et al. (2011) and with radiative-hydrodynamical simulations from Kley et al. (2009) and Bitsch & Kley (2011). The comparison of the theoretical torque curves with data from Kley et al. (2009) and Bitsch & Kley (2011) suggests that with an adjustment of the viscous time scale in the calculation of the saturation mass (the mass when the corotation torque starts to vanish) by a factor of f_visc ≈ 0.55 our model reproduces the torque better.

We also showed the global effects of different parameters of various migration models in a number of population synthesis calculations. Here, the comparison of our nominal BMF model (f_visc ≈ 0.55) and the migration model of Paardekooper et al. (2011) indicates similar results (Sect. 4.2.2) even with the difference shown in the torque curves.

Like Lyra et al. (2010), Kretke & Lin (2012), and Hellary & Nelson (2012), we also find with our prescription that nonisothermal type I migration leads to convergence zones (CZ), that is, points in a disk to which planets migrate to from the inside and outside. As in previous studies (Hellary & Nelson 2012), we find that the convergence zones move inward as the disk evolves and take the captured planets with it. This occurs on the slower time scale of disk evolution and therefore the captured planets are trapped in it and also only migrate on this time scale. The planets leave this zone when their horseshoe drag saturates.

The difference between the migration rate of planets captured in a CZ and the migration rate in the saturated regime is significant. This means that an increase of the mass where saturation occurs by a modest factor of 2, for example, the time spent in rapid-saturated type I migration is significantly shortened and therefore also the extent of migration. The level of saturation at a given mass, and the mass at which saturation begins, are among the critical aspects for the evolution of a giant planet since small changes by a factor of 2 in f_visc can change the final distribution in the mass semimajor axis diagram by a measurable degree. But a similar degree of change in the final semimajor axis-mass distribution was seen when we changed the description of the Lindblad torque to the formula found by Masset & Casoli (2010).

Finally, we determined the formation tracks of a planet, illustrating that under certain conditions, a planet can loose almost all of its gas mass during its evolution. This mass loss can lead to a stop of migration for a few 10⁵ years when the mass loss is strong enough to desaturate the horseshoe drag. The reason for the envelope mass loss is a jump in solid accretion rate, which is caused by the migration of the planet from a solid-depleted into a solid-rich region of the disk. This behavior strongly depends on the treatment of the planetesimal disk. Here we did not evolve the disk or changed the random velocities because of a giant planet. Both change the accretion rate and therefore whether or how this mass loss occurs.

Neither did we consider random variations in the torque due to turbulent density variations in the disk. Recent studies showed that in some cases random migration due to turbulence can dominate the migration behavior for low-mass planets (Papaloizou et al. 2004; Uribe et al. 2011). Depending on the strength of the random torque, it could push planets from the inner into the outer convergence zone even for the stronger Lindblad torque of Paardekooper et al. (2011).

We made no detailed comparison with observed extrasolar planet populations because the migration model is only a small part of the improvements to the overall model, and defer such comparisons to future work. But we quantified the impact of different model settings by studying the fractional yield of different planet types in a synthesis.

With the new nonisothermal migration model described here or that of Paardekooper et al. (2011) there are still about 50% of all low-mass planets lost in the innermost part of the disk. However, it is about a factor two better than the isothermal migration model without any artificial factors. One way to reduce the number of “lost” planets is an increase of the saturation mass. A lower critical mass for runaway gas accretion can also help due to the faster transition into type II migration (Hori & Ikoma 2011). Furthermore, there are hints that there are additional effects leading to outward migration of planets in 3D simulations of magneto-hydrodynamic disks (Uribe et al. 2011).

We also note that up to now there is a shortcoming of all analytic torque predictions: they all neglect the fact that the horseshoe region is over-wide compared with the prediction used here when q/h³ ≈ 1 (Masset et al. 2006). The fast growth of the width of the horseshoe region in that mass range (a few ten Earth masses to a hundred Earth masses) yields a boost of the corotation torque, which is a strong effect. In the same mass range gap opening and the transition into type II migration occur. And as shown in Sect. 2.5, saturation plays a main role in the change of the direction of the torque in that mass range as well. Moreover, our fit of the torques data from Kley et al. (2009) and Bitsch & Kley (2011) depends on the progression of the torque curve in that mass range. This means that even when our fit produced the torque data well, a closer study of the torque in the mass range might uncover new effects, that may change the outcome of population synthesis calculations.

Nevertheless, the CZ or a similar effect might explain the concentration of planets in clusters one finds in distance-mass diagrams of observed extrasolar planets for high- and low-mass planets (Mayor et al. 2011). In particular, the absence of close-in, low-mass planets in Mordasini et al. (2009b) is not seen any longer with the updated migration model.

Additional comparisons of our results with new data from radial velocity measurements and Kepler data will be important, especially when we combine the new migration model with new improvements of our model, i.e., the long term evolution of giant planets (Mordasini et al. 2012b) and the concurring evolution of multiple planets per disk (Ida & Lin 2010; Alibert et al. 2013). Multiple cores in one disk might be collected into one CZ and form one larger core. Therefore the convergence region may function as a large feeding zone of solid matter onto a core captured in a convergence zone and thus enhance the solid accretion rate (Hellary & Nelson 2012; Sándor et al. 2011; Horn et al. 2012). Moreover, in particular the low-mass planets of the inner cluster only migrated through the inner convergence region. The outer boundary of this region is due to the transition in the gas opacity at the iceline. This means that all these planets moved only through the ice-free part of the protoplanetary disk. We will also study the global effects of different viscous α values and deadzones on the behavior of the convergences zones with population synthesis calculations (cf. Hasegawa & Pudritz 2011). This brings us closer to a description of orbital migration without ad hoc efficiency factors.

Online material

Appendix A: Saturation and disk evolution

From the positions of planets at the time of saturation in the formation tracks in Fig. 7 in Sect. 3.5 one sees that most positions lie on two lines, one for the planets of the inner and one for the outer CZ. We study here the reasons for this feature.

Figure A.1 shows the positions of the planets at the transition from unsaturated to saturated adiabatic migration, the most common transition into a saturated migration regime. The positions of all planets undergoing this transition are shown in black, while the colored points show that those with a distance between their semimajor axis and the position of the CZ at that time are fewer than 10% of the semimajor axis. There are two large groups, again one of the inner and one of the outer convergence region. Almost all planets of the inner group are in the CZ when they saturate, while the outer ones are much more spread out. But the planets that are in the CZ form a line here as well. There is a third CZ inside of 0.3 AU in the most massive disks in the synthesis, which leads to a third minor group. However, this CZ evolves quickly and disappears after the first few 0.01 Myr and all associated planets end at 0.1 AU. The dozen points inside of 0.2 AU in Fig. A.1 correspond to planets in this small CZ.

This behavior results from the interaction of the following points:

One can calculate the saturation mass as a function of the orbital distance a and time t by setting s₂ = 1 in Eq. (17). There, f_visc is 1, 0.55 or 0.125 in each of the model versions. ${\mathit{M}}_{\mathrm{sat}}\mathrm{\left(}\mathit{a,t}\mathrm{\right)}\mathrm{=}\frac{\mathit{h}\mathrm{\left(}\mathit{a,t}\mathrm{\right)}}{\mathrm{1.16}}{\left(\frac{\mathrm{8}\mathit{\pi }{\mathit{M}}_{\mathrm{\ast }}\mathit{\nu }\mathrm{\left(}\mathit{a,t}\mathrm{\right)}}{\mathrm{3}{\mathit{f}}_{\mathrm{visc}}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}\frac{\mathrm{1}}{\mathrm{(}\mathit{Ga}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}}\mathrm{·}$(A.1)

For a fixed orbital distance and stellar mass it only depends on the disk aspect ratio h and the viscosity ν. Both quantities are decreasing with time, as the disk mass decreases, and therefore, the saturation mass also becomes smaller as the disk evolves.

While the photoevaporation rate is important for the lifetime of a disk, the constant value of α in all simulations of one synthesis results in a similar disk structure in the part of the disk where viscosity is dominant. There, the disks go through the same series of disk states (radial profile of temperature, surface density, etc.) and only the speed with which the disks go through the states is different and depends on the photoevaporation rate. Given one semimajor axis a, there is only one disk state where the inner (our outer) CZ lie at this position. Therefore also h and ν are fixed for this semimajor axis of the CZ. Therefore, the semimajor axis of the CZ corresponds to only one saturation mass. And as the CZ moves inward while the disk evolves, the saturation mass decreases. Both processes approximately follow power-laws and thus we see a line-like structure for the planets that saturate while they are in one CZ. Planets that saturate early in the disk evolution do so at a higher mass and farther out than planets that saturate in later times of disk evolution.

Finally, the spread in the outer group results from planets that saturated before they reached the CZ. These are planets in disks with high solid surface densities where the planet cores can grow fast. Compared with the inner group, the outer group also contains more planets that saturate outside of the CZ. The larger amount of solids outside of the iceline leads to higher accretion rates. The scatter is reduced when the saturation mass is increased by reducing f_visc. The higher saturation mass gives the planets more time to migrate to the CZ and to saturate there.

Appendix B: Impact of numerical parameters

We made several population syntheses calculation to test the effects of different numerical parameters and comment on the effects here.

We made calculations with the STD model and exponents b = 2.0 and b = 10.0 (nominal value b = 4.0) in the transition function between locally isothermal and adiabatic migration regime (Eq. (20)). This only has a small effect for low-mass planets, while for massive planets the final semimajor axis and mass is almost the same for different values of b. At smaller masses no clear pattern can be seen. An increase of b can lead to either more or less massive planets and to either larger or smaller distances from the star of a few per cent.

We also made calculations with b = 4.0 and a hard jump for the transition function between type I and type II migration (Eq. (23), nominal value b = 10). This only affects massive planets and gives only a change in distance from the star of a few percent, with a larger b leading to planets farther away from the star but almost no change in the final mass.

Overall, the range of the different parameters studied here only leads to minor changes in the overall distribution of planets in semimajor axis and mass.

---

Acknowledgments

We thank Betram Bitsch for providing important data for the radiation hydrodynamic calculations, and Willy Benz, Kees Dullemond, Willy Kley and Sijme-Jan Paardekooper for enlightening discussions. This work was supported in part by the Swiss National Science Foundation and the European Research Council under grant 239605. We thank an anonymous referee whose comments lead to a substantial clarification of the manuscript. Christoph Mordasini thanks the Max-Planck-Gesellschaft for the Reimar-Lüst Fellowship. Kai-Martin Dittkrist thanks the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg.

References

All Tables

All Figures