EVOLUTIONARY TRACKS OF TRAPPED, ACCRETING PROTOPLANETS: THE ORIGIN OF THE OBSERVED MASS–PERIOD RELATION

Adopting the standard disk model, we can estimate the positions of disk inhomogeneities. These positions are crucial because (proto)planets that undergo rapid type I migration will get trapped there. We simply summarize the positions here and refer the reader to Paper I for the complete derivations (also see Table 1).

There are generally three kinds of disk inhomogeneities: dead zones, ice lines, and heat transitions (Paper I). Dead zones are present in the inner region of disks, where high-energy photons such as X-rays from the central stars and cosmic rays cannot penetrate (Gammie 1996; Matsumura & Pudritz 2006; Ilgner & Nelson 2006). The defining feature of the dead zones is the low amplitude of turbulence there that results from the poor coupling of the magnetic field with weakly ionized disks (so that magnetorotational instabilities are suppressed there; e.g., Balbus 2003). Ice lines at which the disk temperature is low enough to trigger condensation of molecules such as water are the most famous and indispensable of disk inhomogeneities (Jang-Condell & Sasselov 2004; Min et al. 2011). They play an important role in population synthesis models (Ida & Lin 2004; Mordasini et al. 2009). Heat transitions are also well recognized in the literature and arise at that radius at which the main heat source changes from viscous heating to stellar irradiation (D'Alessio et al. 1998; Menou & Goodman 2004; Paper I). We have recently demonstrated that the heat transition becomes a planet trap based on analytical arguments (Paper I; also see Kretke & Lin 2012). The validity of the heat transition traps has been recently confirmed by hydrodynamic simulations (Yamada & Inaba 2012).

In principle, the structure of dead zones can be specified by solving the ionization equations (Sano et al. 2000; Matsumura & Pudritz 2006; Ilgner & Nelson 2006). Nonetheless, the resultant structures depend sensitively on disk parameters that are difficult to determine through observations. Therefore, we adopt a parameterized treatment of dead zones (Kretke & Lin 2007; Ida & Lin 2008a; Matsumura et al. 2009; Paper I) in which the effective α in the layered region can be given as

$$\begin{equation} \alpha = \frac{\Sigma _A \alpha _A+(\Sigma _g-\Sigma _A) \alpha _D}{\Sigma _g}, \end{equation} \tag{ 2 }$$

where α_A and α_D are the strength of turbulence in the active and dead layers, respectively, and the surface density of the active layer Σ_A is modeled as

$$\begin{equation} \Sigma _A=\Sigma _{A0} f_{{\rm ice}}\left(\frac{r}{r_0} \right)^{s_A}, \end{equation} \tag{ 3 }$$

where r₀ = 1 AU is the characteristic radius and f_ice can contain the effects of ice lines. More specifically, f_ice represents possible reductions of Σ_A at the ice line that originate from the complex interplay between ice-coated, sticky dust grains and the absorption of free electrons by them (Sano et al. 2000; Ida & Lin 2008b; see Paper I for a more complete discussion). Thus, the structures of dead zones are controlled only by Σ_A0 and s_A in this formalism. This approach is useful because it enables one to investigate how important the structure of dead zones is for understanding the population of planets by simply varying parameters Σ_A0 and s_A (see Section 8). Assuming stationary disk models (see Equation (1)), the surface density of gas is given as

$$\begin{equation} \Sigma _g= \frac{\dot{M}}{3 \pi c_s H \alpha _D} - \Sigma _A \frac{\alpha _A - \alpha _D}{\alpha _D}. \end{equation} \tag{ 4 }$$

With Equation (4) in hand, the position of a dead zone trap is given as

$$\begin{equation} \frac{r_{{\rm dz}}}{r_0} = \left(\frac{\dot{M}}{3 \pi (\alpha _A+\alpha _D) \Sigma _{A0} H_0^2 \Omega _0} \right)^{{1}/({s_A+t+3/2})}, \end{equation} \tag{ 5 }$$

where H₀ and Ω₀ are the pressure scale height and Keplerian frequency at r = r₀, respectively. This can be derived from the assumption that the outer edge of dead zones is specified around Σ_A ∼ Σ_g/2.

The position of the ice line of molecular species k is

$$\begin{equation} \frac{r_{{\rm il}}}{r_0} = \left[ \frac{1}{T_{m,k}^{12}(r_{{\rm il}})} \frac{27 \bar{\kappa }_0 \mu _g \Omega _0^3}{64 \sigma _{SB} \alpha _D \gamma k_B} \left(\frac{\dot{M}}{3 \pi } \right)^2 \right]^{2/9} \propto \dot{M}^{4/9}, \end{equation} \tag{ 6 }$$

where T_{m, k} is the disk midplane temperature below which molecules k can condense, $\bar{\kappa }_0= 2\times 10^{16}$ is the opacity at the ice line of the molecule, μ_g is the mean molecular weight of the gas, k_B is the Boltzmann constant, and γ = 1.4 is the adiabatic index. This is given by recent results that show that viscous heating (rather than stellar irradiation) is generally dominant for determining the position of ice lines (Min et al. 2011; Paper I). This expression is applicable for any molecules. Nonetheless, we focus on water ice lines here since they are likely to be the most important molecule for understanding the observed mass–period relation (Paper I). For ice lines of water, $T_{m,\mbox{H}_2\mbox{O}}(r_{{\rm il}})=170$ K (Jang-Condell & Sasselov 2004).

For the case where ice lines are located within dead zones, the position of the trap needs to satisfy the following condition:

$$\begin{equation} \frac{r_{{\rm il}}}{r_{{\rm dz}}} > \left(h(r_{{\rm dz}}) \frac{\alpha _A + \alpha _D}{\alpha _A - \alpha _D} \right)^{{1}/({s_A+t/2+1})}. \end{equation} \tag{ 7 }$$

Finally, the position of heat transitions is written as

$$\begin{eqnarray} \frac{r_{{\rm ht}}}{r_0} & = & \left[ \frac{1}{T_{m0}} \left(\frac{r_0}{R_*} \right)^{3/7} \left[ \frac{27 \bar{\kappa }_0 \mu _g \Omega _0^3}{64 \sigma _{SB} \alpha _A \gamma k_B} \left(\frac{\dot{M}}{3 \pi } \right)^2 \right]^{1/3} \right]^{14/15} \nonumber\\ & \propto& \dot{M}^{28/45}, \end{eqnarray} \tag{ 8 }$$

where $\bar{\kappa }_0= 2 \times 10^{-4}$ is the opacity at the heat transitions, R_* is stellar radius,

$$\begin{equation} T_{m0} \simeq \left(\frac{1}{H} \right)^{2/7} \left(\frac{T_*}{T_c} \right)^{1/7} T_*, \end{equation} \tag{ 9 }$$

$$\begin{equation} T_c \equiv \frac{GM_* \mu _g}{k_B R_*}, \end{equation} \tag{ 10 }$$

T_* and M_* are stellar effective temperature and mass, respectively, and G is the gravitational constant. We have adopted analytical models of Chiang & Goldreich (1997) for the temperature of the disk midplane heated by stellar irradiation. Planet traps arising from the heat transitions are active only if r_ht > r_dz.

By comparing the positions of each disk inhomogeneity (see Equations (5), (6), and (8)), one immediately observes that disk evolution, which lowers the accretion rate $\dot{M}$, moves them inward, but at different rates (see Figure 2). This is important for understanding the observed mass–period relation.

Zoom In Zoom Out Reset image size

We estimate the characteristic surface densities at the disk inhomogeneities following the above formulation. Although the detailed structures of disk inhomogeneities remain to be simulated, a number of analytical and numerical studies based on the standard viscous disk theory clarified their characteristic structure (Menou & Goodman 2004; Matsumura et al. 2007; Ida & Lin 2008b). These surface densities are utilized for deriving the characteristic masses of planets and following evolutionary tracks of planets that grow in planet traps.

As briefly mentioned above, the outer edge of dead zones is determined by Σ_A ∼ Σ_g/2. Substituting this condition into Equation (4), we find that the characteristic surface density at r_dz is given as

$$\begin{equation} \Sigma _{g,{\rm dz}} \approx \frac{2 \dot{M}}{3 \pi (\alpha _A+\alpha _D) r^2 h^2 \Omega }, \end{equation} \tag{ 11 }$$

where h = H/r is the aspect ratio.

At the ice lines, the surface density is approximately written as

$$\begin{equation} \Sigma _{g,{\rm il}} \approx \frac{\dot{M}}{3 \pi \alpha _D r^2 h^2 \Omega }. \end{equation} \tag{ 12 }$$

We took the mean value of α as ∼α_D. As mentioned before, this assumption is based on the recent extensive studies of ice lines (Sano et al. 2000; Ida & Lin 2008b; Paper I). These studies indicate that ice lines can be regarded as a localized dead zone.

On the other hand, the magnitude of turbulence at the heat transition is expected to be high enough to assume that disks are fully turbulent. As a result, the characteristic surface density at the heat transitions is given as

$$\begin{equation} \Sigma _{g,{\rm ht}} \approx \frac{\dot{M}}{3 \pi \alpha _A r^2 h^2 \Omega }. \end{equation} \tag{ 13 }$$

Viscous turbulence is the dominant driver of disk evolution (e.g., Armitage 2011). We adopt similarity solutions (Lynden-Bell & Pringle 1974) for constraining the relation between the accretion rate $\dot{M}$ and time τ. Similarity solutions are derived from the conservation of the angular momentum of disks (also see Hartmann et al. 1998). Considering a disk that has mass M_d and a characteristic disk radius R_c, its angular momentum J_d can be written as

$$\begin{equation} J_{d} \approx M_d R_c^{1/2}. \end{equation} \tag{ 14 }$$

Following time evolution, where disk material is accreted onto the central star, Equation (14) ensures that R_c is an increasing function of time (since J_d is roughly constant and M_d steadily decreases with time). In a simplified analysis, the expansion rate of R_c can be written as

$$\begin{equation} \frac{d R_c}{dt} \approx \frac{R_c}{\tau _{{\rm vis}}}, \end{equation} \tag{ 15 }$$

where τ_vis = r²/(3ν) is the viscous timescale. Assuming a power-law structure for the disk temperature (T∝r^t), Equation (15) gives

$$\begin{equation} R_c \propto \tau ^{{1}/{(1/2-t)}}, \end{equation} \tag{ 16 }$$

and the total disk mass M_d decreases as (see Equation (14))

$$\begin{equation} M_d \propto \tau ^{-1/(1-2t)}. \end{equation} \tag{ 17 }$$

As a result, the accretion rate is related to time through the following relation:

$$\begin{equation} \dot{M} \propto \tau ^{(-{t+1})/({t-1/2})}. \end{equation} \tag{ 18 }$$

Combining the observations that show that the median accretion rate for classical T Tauri stars (CTTSs) of age ∼1 Myr is ∼10⁻⁸ M_☉ yr⁻¹ (Hartmann et al. 1998) and that $\dot{M} \propto M_*^2$ (Calvet et al. 2004; Muzerolle et al. 2005), we have the following scaling law for accretion rates:

$$\begin{eqnarray} &&\dot{M} \simeq 10^{-8} \,M_{\odot } \mbox{ yr}^{-1} f_{{\rm acc}} \left(\frac{\tau }{ 10^{6} \mbox{ yr}} \right)^{(-{t+1})/({t-1/2})} \left(\frac{M_*}{0.5\,M_{\odot }}\right)^2, \nonumber\\ \end{eqnarray} \tag{ 19 }$$

where we have assumed that the typical mass of CTTSs is ∼0.5 M_☉ and introduced a dimensionless factor f_acc. This factor can be utilized for varying $\dot{M}$ and investigating the subsequent consequences on disk evolution and planet formation.

It is still unclear how gas disks disperse in the final stages of their evolution (e.g., Armitage 2011). One of the leading mechanisms is photoevaporation, which arises from heating up gas by high-energy photons from the surrounding stars and subsequent evaporation of gas due to the thermal pressure (Hollenbach et al. 1994; Johnstone et al. 1998). In principle, photoevaporation rates are determined by the complex interplay between physical and chemical processes that take place in protoplanetary disks being irradiated by their central and nearby massive stars from far-UV (FUV) to extreme-UV (EUV) and up to X-rays (Gorti & Hollenbach 2009, and references therein). Recent extensive studies have investigated how effective photoevaporation is in the dispersal of gas disks. In our models, we adopt a simple scaling law to represent the effects of photoevaporation.

Following the treatment of Adams et al. (2004; see their Appendix for the complete derivation), photoevaporation rates can be scaled as

$$\begin{equation} \dot{M}_{{\rm pe}}=f_{{\rm pe}} N_C \mu _g c_s r_g \left(\frac{r_g}{r} \right) \exp \left(- \frac{r_g}{2r}\right), \end{equation} \tag{ 20 }$$

where f_pe is a dimensionless factor of order unity, N_C is the critical column density of gas that is heated by stellar radiation, and the gravitational radius r_g is given as

$$\begin{equation} r_g= \frac{GM_* \mu _g }{k_BT} \approx 100 \mbox{ AU} \left(\frac{T}{1000\mbox{ K}} \right)^{-1} \left(\frac{M_*}{1 \,M_{\odot }} \right). \end{equation} \tag{ 21 }$$

Utilizing some of the most advanced results of photoevaporation, we further simplified Equation (20). Recently, Gorti & Hollenbach (2009) investigated photoevaporation of gas disks by taking into account radiation of a central star that covers FUV, EUV, and X-rays and found that FUV heating plays the dominant role for inducing photoevaporation at r ≳ 3 AU. This can be understood by the fact that photoevaporation rates are determined by the product of the gas temperature and density. EUV heating that leads to ionizing atomic hydrogen results in higher gas temperatures (∼10⁴ K) than FUV heating (∼10²–10³ K). Nonetheless, the ionization front above which EUV heating dominates can only penetrate the disk atmosphere where gas density is much lower than that where FUV heating becomes dominant. As a result, FUV-induced photoevaporation rates exceed EUV-induced ones.

When photoevaporation is established mainly by FUV heating, the heated outgoing flow acts as an additional source of opacity for the FUV photons (Johnstone et al. 1998; Adams et al. 2004). Consequently, the critical column density heated up by FUV satisfies the self-regulation relation

$$\begin{equation} \tau _{{\rm FUV}} = \sigma _{{\rm FUV}} N_C \sim 1, \end{equation} \tag{ 22 }$$

where σ_FUV ≈ 8 × 10⁻²² cm² is the reasonable cross section of dust grains for the FUV photons. This enables us to specify N_C in Equation (20). Also, the peak of FUV-induced photoevaporation rates is attained around 0.1–0.2r_g, rather than r_g, which is valid for EUV-induced photoevaporation (Adams et al. 2004; Gorti & Hollenbach 2009). This can again be explained by the combination of the gas temperature and density, and it is also confirmed by Equation (20).

Collecting the above arguments, we obtain a simplified but physically motivated scaling law for photoevaporation rates:

$$\begin{eqnarray} &&\dot{M}_{{\rm pe}} \simeq 2.3 \times 10^{-9} \,M_{\odot } \mbox{ yr}^{-1} f_{{\rm pe}} \left(\frac{c_s}{3 \mbox{ km s}^{-1}} \right)^{-1} \left(\frac{M_*}{1 \,M_{\odot }} \right), \nonumber\\ \end{eqnarray} \tag{ 23 }$$

where we have used Equations (21) and (22) and set that r = 0.1r_g in Equation (20).

We note that Equation (23) can be applied for photoevaporation rates induced by both EUV and X-rays despite the fact that it is derived from the physical consideration based on FUV radiation. This can be done by adjusting the dimensionless factor f_pe that is determined by the comparison with more detailed simulations. In fact, the conclusion of Gorti & Hollenbach (2009) that FUV is the dominant source of photoevaporation is still a matter of debate in the literature. This is partly because they relied exclusively on hydrostatic solutions for quantifying winds driven by FUV radiation (although hydrodynamical models are needed for precisely estimating the winds), and partly because they adopted energy spectra that eventually reduce the effects of X-rays. As shown by Owen et al. (2010, 2011, 2012), photoevaporation rates induced by X-rays can attain ∼10⁻⁸ M_☉ yr⁻¹ for the most luminous X-ray sources. This high value is derived from employing observed Chandra spectra of TTSs and is comparable to the photoevaporation rate induced by FUV radiation. As a result, we intentionally avoid specifying the dominant source of photoevaporation of gas. Instead, we consider the general effects of photoevaporation on planet formation by treating f_pe as a free parameter.

We are now in the position to discuss the complete treatment of disk evolution. We assume that the accretion rate through the disk is constant in space and regulated in time by Equation (19). As time goes on, disk material accretes onto the central star and the accretion rate decreases. This change in $\dot{M}$ drives the movement of the disk inhomogeneities (see Equations (5), (6), and (8)). When $\dot{M}$ becomes equal to the photoevaporation rate $\dot{M}_{{\rm pe}}$, represented by Equation (23), we assume that the gas disks disperse completely. Although this treatment is somewhat idealized, it can account for the more detailed simulations. As an example, Gorti et al. (2009) investigated the evolution of viscous protoplanetary disks that are photoevaporated by the FUV, EUV, and X-ray radiation from their central star. They showed that viscous turbulence controls the early stage of disk evolution and the total disk mass gradually decreases initially, which can be formulated by power laws. Once the condition that $\dot{M}\sim \dot{M}_{{\rm pe}}$ is satisfied, the disk mass drops exponentially due to the combination of viscous evolution and photoevaporation of gas. They found that disks of initial mass 0.1 M_☉ around a ∼1 M_☉ star have a lifetime of ∼4 × 10⁶ yr.

We can derive similar results from our treatment by equating $\dot{M}$ with $\dot{M}_{{\rm pe}}$: the disk lifetime is approximately estimated as ∼6 × 10⁶ yr. Thus, our treatment is sufficient for the purpose of representing the evolution of protoplanetary disks that are regulated by both viscosity and photoevaporation.

In summary, we reduce the surface density of gas and the accretion rates, following Equation (19). This results in the movement of the disk inhomogeneities. Also, we locate the final position of each disk inhomogeneity that is determined by the condition that $\dot{M} = \dot{M}_{{\rm pe}}$.

We summarize important parameters that establish the configuration and physical state of protoplanetary disks (see Table 2). They are divided into three sets: stellar parameters (M_*, T_*, and R_*), disk mass (Σ_A0, s_A, f_acc, and t), and disk evolution (α_A, α_D, and f_pe). The disk evolution parameters can be translated into disk lifetimes. These three sets are fundamental to regulate planet formation in the disks, as confirmed in the population synthesis models (Ida & Lin 2004; Mordasini et al. 2009). We focus on disks around CTTSs and denote the set of values given in Table 2 as our fiducial model. Adopting these parameters, disk evolution proceeds from τ_int = 10⁵ yr to the time at $\dot{M}=\dot{M}_{{\rm pe}}$ that defines the disk lifetime (τ_disk). A parameter study in which some of these quantities are changed is presented in Section 8.

Table 2. Important Disk Quantities

Symbols	Meaning	CTTSs
M*	Stellar mass	0.5 M☉
R*	Stellar radius	2.5 R☉
T*	Stellar effective temperature	4000 K
ΣA0	Surface density of active regions at r = r0	20 g cm−2
sA	Power-law index of $\Sigma _A (\propto r^{s_A})$	3
t	Power-law index of the disk temperature (T∝rt)	−1/2
facc	Dimensionless factor for $\dot{M}$ (see Equation (19))	1
αA	Strength of turbulence in the active zone	10−3
αD	Strength of turbulence in the dead zone	10−4
fpe	Dimensionless factor for $\dot{M}_{{\rm pe}}$ (see Equation (20))	1/3

Download table as:  ASCIITypeset image

We draw upon our comprehensive analytical study (Paper I), which showed that the gas surface density and temperature modifications induced by the disk inhomogeneities are significant enough to reverse the direction of rapid type I migration. Therefore, these positions are indeed trapping points of rapid type I migrators.

Figure 2 shows the movement of all three disk inhomogeneities. As demonstrated by Paper I, they all move inward, but at different rates. The inward movements arise from the viscous evolution and the resultant reduction of the surface density of the disks (see Equation (19)). Also, different moving rates for the traps result in complex behaviors of the multiple inhomogeneities, such as convergence (e.g., merging of the heat transition with the dead zone), disappearance, and re-emergence of inhomogeneities (e.g., the behavior of the ice line). When $\dot{M}$ equals $\dot{M}_{{\rm pe}}$, photoevaporation quickly disperses gas in the disk, and hence the movement of the planet traps is terminated. This also determines the lifetime of the disk, which is ∼8.8 × 10⁶ yr in this configuration. The behavior of the planet traps is crucial for understanding the observed mass–period relation later.

For planets of mass smaller than the gap-opening mass M_gap (see below), type I migration is the main agent that governs their orbital distribution. In our models, rapid type I migration is halted at the planet traps. Hence, the location of type I migrators is predicted by the positions of disk inhomogeneities (Paper I; see Figure 2).

It is important to define the minimum mass of planets that will be captured by the planet traps. As demonstrated numerically by Lyra et al. (2010), planets captured in planet traps "dropout" if the following condition is satisfied:

$$\begin{equation} \frac{\tau _{{\rm mig,I}}}{\tau _{\nu }}>1, \end{equation} \tag{ 24 }$$

where τ_{mig, I} is the timescale of type I planetary migration and τ_ν is the timescale that determines the moving rates of planet traps. This relation expresses the fact that trapped planets drop out if the speed of type I migration becomes less than that of the moving traps. In general, τ_{mig, I} is scaled as

$$\begin{equation} \tau _{{\rm mig,I}}= \frac{M_p r_p^2 \Omega _p }{2 \Gamma }, \end{equation} \tag{ 25 }$$

where

$$\begin{equation} \Gamma = K_{{\rm mig}} \left(\frac{M_p}{M_*} \right)^2 \frac{\Sigma _{g,p} r_p^4 \Omega _p^2}{h_p^2} \end{equation} \tag{ 26 }$$

with K_mig = 1–10, depending on the optical thickness of the disk (Paardekooper et al. 2010). For stationary accretion disk models (see Equation (1)), the minimum mass of type I migrators that can be captured at planet traps is given as

$$\begin{equation} M_{{\rm mig,I}} = \frac{h_p^2 M_*^2}{2 K_{{\rm mig}} \Sigma _{g,p} r_p^2 \Omega _p \tau _{\nu }}. \end{equation} \tag{ 27 }$$

We set τ_ν = 10⁶ yr, because the moving rates of planet traps are eventually regulated by disk lifetimes (see Figure 2) and the observations revealed that the disk lifetime of any CTTS disk is an order of Myr.

The gap-opening mass M_gap distinguishes type I migration from type II and is well discussed in the literature (e.g., Ward 1997; Matsumura & Pudritz 2006). It arises when a planet becomes sufficiently massive that the torque it exerts on the disk opens a gap. There are two main arguments for estimating M_gap. The first one is the Hill radius analysis: the Hill radius should be larger than the pressure scale height for maintaining gap formation; otherwise, gaps are closed by the gas pressure. The second argument arises from viscous disks. Disk viscosity that controls disk evolution plays the counteractive role for gap formation. Therefore, the tidal torque of a planet on the disk opens a gap if it exceeds the viscous torque. Summarizing these arguments, M_gap is given as (e.g., Matsumura & Pudritz 2006)

$$\begin{equation} \frac{M_{{\rm gap}}}{M_*} = \mbox{min} \left[ 3 h_p^3, \sqrt{40 \alpha h_p^5} \right]. \end{equation} \tag{ 28 }$$

Planets of mass larger than M_gap open up a gap in their disks and undergo so-called type II migration. In the type II regime, we define two characteristic masses. One of them is the critical mass (M_crit) above which the inertia of type II migrators is significant enough to prevent type II migration from proceeding as disks evolve (otherwise, the timescale of type II migration is given as τ_{mig, II} ∼ τ_vis). This effect is also known as a damming effect (Syer & Clarke 1995; Ivanov et al. 1999). The critical mass M_crit is defined by the local disk mass:

$$\begin{equation} M_{{\rm crit}} = \pi \Sigma _{g,p} r_p^2, \end{equation} \tag{ 29 }$$

where Σ_{g, p} is the surface density of gas disks at the position of a planet (r = r_p).

The other characteristic mass is the maximum mass of planets. In general, gas accretion onto cores of gas giants is not fully terminated even if they form a gap in their disks (Lissauer et al. 2009). This suggests that the possible maximum mass of planets that start forming at time τ can be estimated as

$$\begin{eqnarray} M_{{\rm max}}(\tau) & \simeq & \int ^{\tau _{{\rm disk}}}_{\tau } d\tau \dot{M} \nonumber\\ & = & 5 \times 10^{-3} \,M_{\odot } f_{{\rm acc}}\left(t - \frac{1}{2} \right) \left(\frac{M_*}{0.5 \,M_{\odot }} \right)^{({2t-1})/({t-1})} \nonumber\\ && \times \left[ \left(\frac{\dot{M}(\tau _{{\rm disk}})}{10^{-8} f_{{\rm acc}}\,M_{\odot } \mbox{ yr}^{-1}} \right)^{-1/(2(t-1))} \right. \nonumber\\ && \left. - \left(\frac{\dot{M}(\tau)}{10^{-8} f_{{\rm acc}}\,M_{\odot } \mbox{ yr}^{-1}} \right)^{-1/(2(t-1))} \right], \end{eqnarray} \tag{ 30 }$$

where Equation (19) is used.

In conclusion, the trapping regime is defined by M_gap and M_{mig, I}, in which type I migrators follow the movement of the planet traps, while the type II regime is defined by M_max and M_gap, wherein the radial distribution of planets is established by the type II migration (see the bottom right panel of Figure 10 in Appendix A).

The formation of gas giants is divided mainly into three stages (Wetherill & Stewart 1989; Kokubo & Ida 1998; Pollack et al. 1996): formation of rocky cores through runaway and oligarchic growth (Stage I), the subsequent slow gas accretion of the cores and formation of envelopes surrounding them (Stage II), and collapse of the envelopes and runaway gas accretion onto their cores (Stage III). In order to model these three physical processes, we adopt the formulation of Ida & Lin (2004), who first attempted to understand the statistics of the observed exoplanets by carrying out population synthesis analyses. In this formulation, these processes are treated by simple, analytical prescriptions that are derived from the detailed numerical simulations. In Appendix B, we briefly describe our treatments which slightly modify the original formulation, and we refer the readers to Ida & Lin (2004) for a complete discussion. We also present a parameter study in Appendix C for confirming the validity of our tiny modifications.

Orbital evolution of planets is governed by planetary migration that arises from tidal interactions of the planets with the surrounding gaseous disks (Ward 1997; Tanaka et al. 2002). As discussed in Section 5, the characteristic masses will classify planetary migration to four modes, depending on planetary mass: slower type I, trapped type I, the standard type II, and slower type II migration (see Table 3). We discuss our treatments of them below.

Slower type I migration. This mode is applicable if planetary mass is smaller than M_p < M_{mig, I} (see Equation (27)). When planets satisfy this condition, the migration rate of these planets is much smaller than the moving rate of gas that is regulated by disk viscosity. This is the reason why we call this mode of migration the slower type I migration. Therefore, we assume that these planets remain in the same position with time.

Trapped type I migration. When planets are in the trapping regimes, that is, M_{mig, I} ⩽ M_p ⩽ M_gap (see Equations (27) and (28)), the radial positions of these planets follow the movement of planet traps.

Standard type II migration. When the mass of planets is in the range M_gap ⩽ M_p ⩽ M_crit (Equations (28) and (29)), they undergo type II migration that proceeds as the gas disks evolve; the type II migration timescale τ_{mig, II} equals τ_vis. Therefore, the planets move inward with the velocity, written as

$$\begin{equation} v_{{\rm mig,II}} \simeq - \frac{\nu }{r}. \end{equation} \tag{ 31 }$$

Slower type II migration. For planets with M_p ≳ M_crit, on the contrary, the type II migration rate slows down due to the inertia of the planets (Syer & Clarke 1995; Ivanov et al. 1999). That is why we refer to this mode as slower type II migration. As a result, the velocity of the planets becomes

$$\begin{equation} v_{{\rm mig, slowII}} \simeq - \frac{\nu }{r(1+f_{{\rm mig,slowerII}} M_p/M_{{\rm crit}})}, \end{equation} \tag{ 32 }$$

where we have followed Hellary & Nelson (2012) for taking into account the effects of the inertia of planets. In addition, we have introduced a new free parameter f_{mig, slowerII}. As shown below, both the trapped type I migration and slower type II migration are important agents that regulate the radial distribution of planets in our model. Furthermore, it is currently uncertain how effective the inertia of planets is in slowing down standard type II migration. Thus, it is useful to clarify the role of slower type II migration by performing a parameter study wherein the value of f_{mig, slowerII} varies (see Section 8.3). We set f_{mig, slowerII} = 1 for our fiducial model.

We adopt the disk models discussed in Section 3. More specifically, we use the characteristic surface densities at three disk inhomogeneities. In addition, the surface density of dust Σ_d is required to examine planetary growth there.

We simply assume that

$$\begin{equation} \Sigma _d=f_{{\rm dtg}} \Sigma _g, \end{equation} \tag{ 33 }$$

where f_dtg is the dust-to-gas ratio. Table 4 summarizes the values of f_dtg at each disk inhomogeneity. We took f_dtg = 0.01 at the dead zone, because the dust mass is canonically about a hundredth of the gas mass in protoplanetary disks (e.g., Dullemond et al. 2007). At the ice line, condensation of water increases the dust density there and beyond. Therefore, we used f_dtg = 0.05 at the heat transition. The reason that f_dtg = 0.01 at the ice line is that we have already taken into account the effect of the ice line on Σ_g by reducing the mean value of α (see Equation (12)). As a result, f_dtg = 0.01 is reasonable for specifying Σ_d there.

We choose a value for the initial mass of cores of ≃ 0.01 M_⊕, which is sufficiently smaller than the mass that is finally obtained by the oligarchic growth (Kokubo & Ida 1998, 2002). We confirmed that this choice does not affect our results.

The cores start growing at a position r at a time τ. In principle, core formation takes place anywhere in disks. Nonetheless, we assume that the cores will quickly end up on one of the traps in the initial setup. It is noted that the assumption does not always assure the cores to be initially captured at their traps. This is because trapping happens only if the mass of the cores is larger than M_{mig, I} (see Equation (27)). Although one may consider the assumption of the initial τ and r to be somewhat artificial, this is not the case. As shown by Ida & Lin (2008a), planetary cores that undergo rapid type I migration do not contribute to the population of gas giants (since they plunge into their central star within the disk lifetime). This implies that only the cores that experience slower type I migration will play an important role for reproducing the observed gas giants. Thus, it is reasonable to focus on planet formation proceeding only in planet traps that can substantially slow down the type I migration.

Based on the assumption, it is only necessary to choose a distribution of the initial time τ (or position r) for the growth of cores to begin. The positions of disk inhomogeneities are related to the time τ through the accretion rate (see Equation (19)). Table 5 summarizes our seven choices of τ, which are selected to cover the entire disk lifetime within which planetary growth and migration take place. For reference purposes, the initial positions that are determined by Equations (5), (6), and (8) are also shown in the same table. It is noted that the multiple choices of the initial time result in forming multiple planets in each planet trap. We emphasize that the productivity of each planet trap—how many planets eventually form in each planet trap during the disk lifetime—and the relation between the number of finally formed planets and planet traps should be investigated separately.

We neglect the planet–planet interactions of the cores that grow in different planet traps—we leave this for our future work.

We may now follow the evolutionary tracks in the mass–semimajor axis diagram for planets that grow in all three planet traps. We adopt the above analytical prescriptions for planetary growth and migration. We summarize our technical procedures here. The standard treatment of mass accretion and planetary growth is given in Appendix B.

When the mass of protoplanets is less than M_{mig, I}, their mass increases with time following the standard oligarchic growth (see Equation (B2)), while their semimajor axes remain roughly the same. Time evolution also reduces the surface density of disks (gas and dust), and hence the growth rate also changes with time (see Equation (B2)). Once they acquire masses that are larger than M_{mig, I}, they start to migrate inward. When they are at their planet traps, they move inward at the same rate as their traps. If they are left behind, they quickly catch up with their traps due to the standard rapid type I migration and then follow the movement of the traps. If their planet traps disappear due to convergence with other traps, it is assumed that the planets follow new planet traps that survive the convergence. If the planets become more massive than the critical mass of cores above which their envelopes cannot maintain hydrostatic equilibrium (see Equation (B7)), then accretion of gas onto the cores begins. The gas accretion rates are regulated purely by the mass of cores (Equation (B8)). Through our experiments, we find that, for most cases, core formation is completed when they are captured in their traps.

When a planet's mass reaches the gap-opening mass, it undergoes standard type II migration. This results in "dropping out" of the trapped planet from its planet trap and happens because the planet is now so massive that it opens up a gap in the disk (which leads to different orbital evolutions between the planet and the planet trap). When the planets are within the dead zone, type II migration becomes slower through a low value of α(= α_D), while, for the planets outside the dead zone, the value of α_A is used for the migration. If the planet attains the mass of f_maxM_max, where f_max is a controllable parameter (see Appendix B), its accretion is terminated. M_max is a decreasing function of time (see Equation (30)), so that planets that need a long time to grow up to gas giants tend to be less massive, while planets that can quickly become gas giants tend to be more massive. Even when planet formation is largely complete, their disks may still have a sufficient amount of gas to drive type II migration. In this case, the type II migration is slowed down by the inertia of the planets. The accretion rate $\dot{M}$ declines with time, and at a certain time $\dot{M}$ becomes equal to photoevaporation rates $\dot{M}_{{\rm pe}}$. When this is satisfied, the positions of the planets freeze-in the mass–semimajor axis diagram.

The evolutionary track of a growing planet consists of four distinct phases in the mass–semimajor axis diagram (see Figure 3). The first phase is formation of cores of gas giants through runaway and oligarchic growth (e.g., Wetherill & Stewart 1989; Kokubo & Ida 1998). The mass of the core in this phase is high enough to keep up with the movement of the trap, while the torque of the core acting on the disk is too weak to open up a gap there. Thus, the core remains within a trapping regime and follows its movement. The timescale of this phase is of order ∼10⁵–10⁶ yr, which is much shorter than the disk lifetime (τ_disk ∼ 8.8 × 10⁶ yr in this setup) as shown in previous studies (Kokubo & Ida 2002). Hence, the protoplanet moves upward in mass while moving little in orbital radius or period.

Zoom In Zoom Out Reset image size

As the feeding zone empties, core formation is terminated and the second phase begins, wherein gas accretion onto the envelope occurs. It was well known that the timescale of this phase was problematically long for earlier models (≳ 10⁷ yr; Pollack et al. 1996). However, recent studies improved the previous models and revealed that the timescale is highly sensitive to the optical depth of the envelope. For these realistic conditions, it is significantly shorter than the disk lifetime (Lissauer et al. 2009). In our calculation, this timescale is about 2 × 10⁶ yr, and hence the core still has sufficient time to finally grow up to a gas giant within τ_disk. The core during most of this phase is trapped. As a result, its radial evolution is mainly determined by the slow movement of the dead zone trap, and the protoplanet moves to shorter radii and periods while at nearly a constant mass. Toward the end of this phase, the core becomes massive enough that the tidal torque it exerts upon the disk becomes comparable to the viscous torque that evolves gas disks, leading to gap formation in the disks and type II migration of the planetary core.

When the mass of the gaseous envelope cannot be supported by the gas pressure, runaway gas accretion onto the core takes place (Phase III). The timescale of this phase is very short (≲ 10⁵ yr), and consequently its evolutionary path is almost vertical in the mass–semimajor axis diagram. These three successive phases are the main path to forming gas giants in the core accretion scenario (Pollack et al. 1996; Lissauer et al. 2009). The massive planet opens up a gap in the disk and undergoes type II migration. This switch from type I to type II migration results in "dropping out" of the planet from the moving trap and decouples it from the movement of the planet trap.

The onset of Phase IV completes the formation of a gas giant. During this phase (≳ 10⁶ yr), type II migration moves the gas giant inward further. However, this process is minimized by the inertia of the massive planet (Syer & Clarke 1995; Ivanov et al. 1999).

Planets arrive at their final position in the mass–semimajor axis diagram when the disk is finally dissipated. Photoevaporation of the disk by high-energy radiation from the central star is likely to be the dominant mechanism of gas dispersal in the disks (e.g., Gorti & Hollenbach 2009) and will terminate type II migration. As a result, the final orbital period and mass of the planet are achieved. Thus, Figure 3 summarizes how concurrent evolution of planetary growth and migration proceeds in the mass–semimajor axis diagram: a core is formed in a dead zone trap that is initially located at ∼7 AU. Following the movement of the dead zone trap, the core is transported to ∼3 AU. Simultaneously, it undergoes the two main phases of gas giant formation. The completion of the final runaway gas accretion onto the core and subsequent type II migration involve further evolution of the planet in the diagram. When photoevaporation becomes important, the gas disk is removed and the position of the planet in the diagram "freezes-out."

Figure 4 shows the computed evolutionary tracks of planets that grow at all three disk inhomogeneities. Different lines at each planet trap correspond to different evolutionary tracks in which planetary growth starts at different times (see Table 5). Despite the difference in the starting time (and position), most planets formed at the dead zone and heat transition traps end up at r ∼ 1 AU (∼500 days) and r ∼ 0.1 AU (∼10 days), respectively. At the heat transition trap, the surface density of dust is low. Therefore, cores that grow there spend a long time in the trapping phases (Phases I and II). This maximizes the distance over which cores are transported and results in the distribution of cores that hover preferentially around ≳1 AU. Since the low-mass cores get distributed over smaller orbital radii and less time remains for the cores to grow up to gas giants, they finally remain less massive (≲ 100 M_⊕) and are located around smaller orbital radii (∼0.1 AU). The same argument is applied to planets formed in the dead zone trap. However, the surface density of dust at the dead zone is considerably higher than that at the heat transition. Consequently, the final mass of cores trapped at dead zones becomes larger, core formation completes earlier, and the distribution of cores is shifted to ∼3 AU. These combined differences result in the populations of more massive planets orbiting at ∼1 AU.

Zoom In Zoom Out Reset image size

The evolutionary tracks associated with protoplanets carried by the ice line trap show some differences. The resultant planetary population spreads out over a wider range in the mass–semimajor axis diagram (see Figure 4). Nonetheless, this can be also understood by the surface density of dust and the resultant core formation there. At the ice line, the surface densities are substantially higher than that at the dead zone and heat transition, and hence the formation of cores is most efficient. This typically results in most massive cores. At the early stage of disk evolution, therefore, the most massive cores are preferentially formed there. They can readily dropout from the moving trap and pileup around larger orbital radii (r ∼ 5 AU). These massive cores at larger orbital radii lead to the formation of more massive gas giants that finally orbit at ≳1 AU. In the later stage of disk evolution, the high dust densities at the ice line can still form cores, while at that time the other traps are not due to lower dust density. This is the physical reason of the wide spread of planetary population due to the ice line traps.

We now compare our results with observations. As already presented in Figure 4, our model shows that the superposition of all tracks for planets that grow in three planet traps constitutes a theoretical mass–period relation, wherein the final distribution of the mass of the planets is an increasing function of their periods. This is consistent with the observed mass–period relation, as the observational data scatter around the locus of end points of our tracks (see Figure 5).

Zoom In Zoom Out Reset image size

This is one of the most important findings in this paper. As discussed in Section 7.2, this arises from the fact that there are considerable differences in the properties of the planet traps that regulate planet formation and migration. As a result, different planet traps have different preferred loci at which evolutionary tracks end up in the mass–semimajor axis diagram. Thus, planet traps act as a filter for distributing cores—massive cores readily drop out from moving traps and tend to orbit farther away from the central star, while low-mass cores are trapped for a long time and tend to orbit close to the host star—and play a central role in generating the theoretical mass–period relation.

In addition, the prediction that distinct sub-populations can arise depending on the trapping mechanism has several observational consequences. For example, our model provides a physical explanation for the observed pileup of gas giants at ∼1 AU. This again relies on the argument that planet formation efficiency highly depends on the surface density of dust at planet traps. At the dead zone and ice lines, the dust density is expected to be high due to the low disk turbulence, and hence planet formation rates are high there. On the other hand, the formation rate would be low at the heat transition trap due to low dust density. This results in a general trend that more planets are readily formed at the dead zone and ice line traps that end up at r ∼ 1 AU (see Figure 5).

Furthermore, our model predicts the population of low-mass planets (≲50 M_⊕) with r ≲ 0.5 AU. This arises from planet formation that takes place in the moving ice line trap (see Figure 5). Even in the later stage of disk evolution, the highest dust density there enables the formation of low-mass planets that end up in the desert. On the contrary, the most advanced population synthesis models predict a planet desert there (Ida & Lin 2004, 2008b; also see footnote ³ in Section 1). The presence of the many observed exoplanets in the region agrees well with our findings.

Finally, our models predict the existence of planet deserts that are quite different in the mass–period space than those claimed by Ida & Lin (2004, 2008b). Figure 6 shows our deserts, denoted by hatched regions. They are produced due to trapping and subsequent transport of cores. This leads to the evacuation of the cores from these regions in which they have initially grown up. As a result, these regions are regarded as void of planets. More specifically, we define our deserts by estimating the mass ranges of planets that can be captured at the planet traps and following their movement: M_p < M_gap and τ_{mig, I} < τ_vis (see Section 5). This kind of planet desert is active only for gas disks. There are a number of possibilities to fill out our deserts: that successive formation of rocky planets after gas disks disperse may ultimately fill out the regime; that, even in the epoch of gas disks, planetary cores formed far beyond our deserts may eventually distribute there due to planetary migration; and that planet–planet scatterings induced by convergence of multiple planet traps may deliver the scattered cores into our deserts. Nonetheless, our predictions are useful in the sense that such regions are the primary target of the current and ongoing observational surveys (Mayor et al. 2011; Howard et al. 2012).

Zoom In Zoom Out Reset image size

We first focus on parameters for dead zones. We have adopted the parameterized treatment for the structures of dead zones, wherein they are represented by Σ_A0 and s_A (see Equation (3)). Even in the most recent studies, it is still somewhat uncertain what the precise structure of the dead zones is (Matsumura & Pudritz 2006; Martin et al. 2012). Therefore, we utilize our parameter study in order to discuss how sensitive our findings are to the structures of the dead zones.

Table 6 summarizes parameters we varied. For Runs A1 and A2, the value of Σ_A0 is changed, while s_A varies for Runs A3 and A4. Any other parameters remain the same as the fiducial ones for all four runs. Figure 7 shows the results of the evolution of the positions of three disk inhomogeneities (the left column), the evolutionary tracks of planets (the central column), and the trapping regimes (the right column). The top panels are for the case of Run A1, the second for Run A2, the third for Run A3, and the bottom for Run A4. One immediately observes that the results for all four cases, especially the behaviors of the evolutionary tracks, are surprisingly similar to those of the fiducial case and give all three key results. Thus, we can conclude that our findings discussed in Section 7.3 are robust even if the structures of dead zones somewhat change due to their surrounding environments and disk configurations.

Zoom In Zoom Out Reset image size

The above parameter study leads to the conclusion that disk inhomogeneities and the resultant multiple planet traps play a crucial role in reproducing the key properties of observed exoplanets. Nonetheless, there is one population of planets that is not covered by the fiducial model: gas giants orbiting at r ≳ 5 AU. We now examine whether or not this population is also predicted by our model. In order to proceed, we change the stellar mass from 0.5 to 0.9 M_☉ and keep other parameters the same as the fiducial ones. The main motivation for changing the stellar mass is that the observational data are obtained from low- to high-mass stars that cover F, G, and K stars.

Figure 8 shows the results of the movement of the disk inhomogeneities, the evolutionary tracks of planets that are formed in the traps, and the locations of the planet deserts. This figure confirms that the planets not covered in our previous setup (10² M_⊕ ≲ M_p ≲ 5 × 10³ M_⊕ and 5 AU ≲ r ≲ 10 AU) can indeed be reproduced. This is because high-mass stars result in high accretion rates (see Equation (19)), which corresponds to the situation that disks have high mass. As a result, planet formation efficiencies at all three planet traps become high, and planets readily attain high mass, which ends up with planets distributing farther away from the central star. Thus, this finding indicates that the full range of the statistical properties of exoplanets can be understood by our model.

Zoom In Zoom Out Reset image size

Based on the above discussion, a number of the fundamental statistical properties of observed exoplanets are most likely to be explained by multiple planet traps that capture and transport planetary cores, depending on planetary mass. As shown in Figure 4, however, slower type II migration also drives radial drifts in the mass–semimajor axis diagram, following evolutionary tracks of growing planets.⁴ In addition, the slower type II migration is a function of planetary mass (see Equation (32)). Thus, they imply that a combination of planet traps and slower type II migration (not only planet traps) may play a central role in reproducing the observations. In order to examine the effects of the slower type II migration on our results, we carry out a parameter study in which the value of f_{mig, slowerII} changes (see Table 7).

Figure 9 shows the results of the evolutionary tracks of planets growing in all three planet traps. The top left panel shows the results of Run B1, the top right for Run B2, the bottom left for Run B3, and the bottom right for Run B4. For the top panels, there is no mass dependency in Equation (32), while it depends on planetary mass for the bottom panels. Comparison of Run B1 (top left) with any other runs demonstrates that a process that slows down standard type II migration is clearly needed for reproducing the observed population of exoplanets even if planetary cores are saved due to planet traps. This is because standard type II migration that proceeds as a local viscous timescale leads planets to spiraling into the host stars within the disk lifetime at r ≲ 10 AU.

Zoom In Zoom Out Reset image size

When type II migration slows down sufficiently, our findings discussed above, especially the mass–period relation, are valid for various cases (see top right and two bottom panels). It is important that the overall feature of the resultant populations is very insensitive to the origins of slowing-down mechanisms (mass dependence versus independence; compare the top right panel with Figure 5), although the radial distribution varies somewhat for each case. We confirmed, through experiments, that the observed planetary population can be reproduced when f_{mig, slowerII} ≳ 1 or $f_{{\rm mig,slowerII}} = \bar{f} M_{{\rm crit}}/M_{p}$ with $\bar{f} \gtrsim 1$.

In addition to the slower type II migration, photoevaporation of gas is also important for establishing the final radial distribution of planets, since it terminates the migration. The combined effects of type II migration and photoevaporation on planetary populations were discussed in the literature (e.g., Armitage et al. 2002; Matsuyama et al. 2003; Armitage 2007; Alexander & Armitage 2009; Alexander & Pascucci 2012). Our models are more fundamental than theirs in the sense that we have incorporated the effects of type I migration that can be captured at planet traps. We will investigate in detail the combined effects in a forthcoming paper.

In summary, we can conclude that planet traps, not slower type II migration, play the primary role in reproducing the observations, provided that the standard type II migration slows down sufficiently when planetary mass exceeds the local disk mass.

Formation of cores in planetesimal disks is well understood in the literature (e.g., Kokubo & Ida 2002), and the growth timescale in the disks is given as

$$\begin{eqnarray} \tau _{c,{\rm acc}} & \simeq & 1.2 \times 10^{5} \mbox{ yr} \left(\frac{\Sigma _d}{10 \mbox{ g cm}^{-2}} \right)^{-1} \left(\frac{r}{r_0} \right)^{1/2} \nonumber\\ &&\times \left(\frac{M_c}{M_{\oplus }} \right)^{1/3} \left(\frac{M_*}{M_{\odot }} \right)^{-1/6} \nonumber\\ && \times \left[ \left(\frac{b}{10} \right)^{-1/5} \left(\frac{\Sigma _g}{2.4 \times 10^3 \mbox{ g cm}^{-3}} \right)^{-1/5}\right. \nonumber\\ &&\left.\times \left(\frac{r}{r_0} \right)^{1/20} \left(\frac{m}{10^{18} \mbox{ g}} \right)^{1/15} \right]^2, \end{eqnarray} \tag{ B1 }$$

where Σ_d is the surface density of dust, M_c is the mass of a core, b = 10 is a parameter for determining the feeding zone of the core (see below), and m = 10¹⁸ g is the mass of planetesimals that are accreted onto the core. Adopting this timescale, the growth of cores is regulated by

$$\begin{equation} \frac{d M_p}{dt} =\frac{M_p}{\tau _{c,{\rm acc}}}. \end{equation} \tag{ B2 }$$

As cores grow, their feeding zones ▵r_c empty. These zones are scaled by br_H, where r_H = (M_c/(3M_*))^1/3 is the Hill radius of a core and b ∼ 10. The decrease of planetesimals in the zones results in the reduction of the accretion rate of cores $\dot{M}_c$. In the limit of the modest to high velocity dispersion σ of planetesimals that are accreted onto the cores, $\dot{M}_c$ can be written as (Safronov 1972; Ida & Lin 2004)

$$\begin{equation} \dot{M}_c \sim 2 \pi \left(\frac{R_c}{r} \right) \left(\frac{M_c}{M_*} \right) \left(\frac{r\Omega }{\sigma }\right)^{2} \Sigma _d r^2 \Omega, \end{equation} \tag{ B3 }$$

where R_c is the radius of cores. Planetesimals within the feeding zones can reach cores when σ/Ω ∼ ▵r_c(= br_H). Assuming R_c to be similar to that of the Earth,

$$\begin{equation} R_c = 6.4 \times 10^{8} \mbox{ cm} \left(\frac{M_c}{M_{\oplus }}\right)^{1/3} \left(\frac{\rho _c}{5.5 \mbox{ g cm}^{-3}} \right)^{-1/3}, \end{equation} \tag{ B4 }$$

$\dot{M}_c$ is given as

$$\begin{eqnarray} \dot{M}_c & \sim & 3.0\times 10^{-8} \,M_{\oplus } \mbox{ yr}^{-1} \left(\frac{b}{10} \right)^{-2} \left(\frac{\rho _c}{5.5 \mbox{ g cm}^{-3}} \right)^{-1/3} \nonumber\\ &&\times \left(\frac{M_c}{M_{\oplus }} \right)^{2/3} \left(\frac{M_*}{M_{\odot }} \right)^{-1/3} \left(\frac{\Sigma _d}{10 \mbox{ g cm}^{-2}} \right) \nonumber\\ &&\times \left(\frac{r}{r_0} \right) \left(\frac{\mbox{ yr}}{1 /\Omega } \right). \end{eqnarray} \tag{ B5 }$$

When all the planetesimals in their feeding zones are consumed, the cores attain the maximum mass that is known as the isolation mass, which is defined by (Kokubo & Ida 2002; Ida & Lin 2004)

$$\begin{eqnarray} M_{c,{\rm iso}} &=& 2 \pi r \triangle r_c \Sigma _d \simeq 0.16 \,M_{\oplus } \left(\frac{b}{10} \right)^{3/2} \left(\frac{\Sigma _d}{10 \mbox{ g cm}^{-2}} \right)^{3/2} \nonumber\\ &&\times \left(\frac{r}{r_0} \right)^{3} \left(\frac{M_*}{M_{\odot }} \right)^{-1/2}. \end{eqnarray} \tag{ B6 }$$

The accretion of gas onto cores and subsequent envelope formation are initiated when the mass of the cores becomes larger than

$$\begin{equation} M_{c,{\rm crit}} \simeq 10 f_{c,{\rm crit}} \,M_{\oplus } \left(\frac{\dot{M}_c}{10^{-6}\,M_{\oplus } \mbox{ yr}^{-1}} \right)^{1/4}. \end{equation} \tag{ B7 }$$

This critical mass was originally derived from a series of numerical simulations that investigated the effect of the cores' accretion rates and opacity in the envelope on formation of gas giants (Ikoma et al. 2000). Here, we have adopted a simplified one, following Ida & Lin (2004). Recent studies, however, revealed that M_{c, crit} might be smaller than that predicted by Equation (B7) with f_{c, crit} = 1 (e.g., Hori & Ikoma 2011). In order to take this into account, we have introduced a dimensionless factor f_{c, crit}.

The gas accretion rate of cores is prescribed by (Ida & Lin 2004)

$$\begin{equation} \frac{d M_p}{dt} \simeq \frac{M_p}{\tau _{{\rm KH}}}, \end{equation} \tag{ B8 }$$

where the Kelvin–Helmholtz timescale is given as

$$\begin{equation} \tau _{{\rm KH}} \simeq 10^{c} \mbox{ yr} \left(\frac{M_p}{M_{\oplus }} \right)^{-d}, \end{equation} \tag{ B9 }$$

where c = 8 and d = 2.5. This is a simplified timescale that was originally estimated from numerical simulations (Ikoma et al. 2000). As shown by the more detailed numerical simulations (Pollack et al. 1996; Lissauer et al. 2009), this stage is slow (≳10⁶ yr).

When planets become massive enough, runaway gas accretion onto the cores starts. In the detailed numerical simulations, this stage commences when the envelope of cores becomes more massive than the cores (Pollack et al. 1996; Ikoma et al. 2000; Lissauer et al. 2009). Nonetheless, we monitor this stage by the condition that

$$\begin{equation} \frac{\tau _{{\rm KH}}}{10^5 \mbox{ yr}} < 1. \end{equation} \tag{ B10 }$$

This is because our approach is rather simple. This condition never affects our results. The growth rate of this stage is also prescribed by Equation (B8).

It is totally unclear how gas accretion onto the cores terminates and what physical process(es) determines the final mass of gas giants. Therefore, we assume that Stage III continues until planets gain the mass f_maxM_max, where f_max is an adjustable parameter.

As discussed above, planetary growth and consequent evolutionary tracks of planets are regulated by four parameters in our model (see Table 8). The parameter f_{c, crit} governs the onset of gas accretion of cores, a set of parameters c and d determine the efficiency of gas accretion onto cores, and f_max controls the final mass of planets. We denote the values given in Table 8 our fiducial model. Note that Ida & Lin (2004) adopted the values of f_{c, crit} = 1, c = 9, and d = 3 rather than our set of f_{c, crit} = 0.3, c = 8, and d = 2.5. We performed a parameter study and confirmed that the different choice of these values does not change our results very much. Therefore, our choice and results are robust in a sense that we can compare our results with the results of Ida & Lin (2004, 2008b). We present only a parameter study of f_max in Appendix C, because the choice of this value is probably the most uncertain.