Geoscience Frontiers 8 (2017) 215e222

H O S T E D BY

Contents lists available at ScienceDirect
China University of Geosciences (Beijing)

Geoscience Frontiers
journal homepage: www.elsevier.com/locate/gsf

Research paper

Population synthesis of planet formation using a torque formula withdynamic effects
Takanori Sasaki a, *, Toshikazu Ebisuzaki b
ab

Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502, JapanRIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japana r t i c l e i n f o

a b s t r a c t

Article history:
Received 24 February 2016
Received in revised form
2 April 2016
Accepted 11 April 2016
Available online 2 May 2016

Population synthesis studies into planet formation have suggested that distributions consistent withobservations can only be reproduced if the actual Type I migration timescale is at least an order ofmagnitude longer than that deduced from linear theories. Although past studies considered the effect of
the Type I migration of protoplanetary embryos, in most cases they used a conventional formula basedon static torques in isothermal disks, and employed a reduction factor to account for uncertainty in themechanism details. However, in addition to static torques, a migrating planet experiences dynamic
torques that are proportional to the migration rate. These dynamic torques can impact on planet
migration and predicted planetary populations. In this study, we derived a new torque formula for Type I
migration by taking into account dynamic corrections. This formula was used to perform population
synthesis simulations with and without the effect of dynamic torques. In many cases, inward migration
was slowed signiﬁcantly by the dynamic effects. For the static torque case, gas giant formation was
effectively suppressed by Type I migration; however, when dynamic effects were considered, a substantial fraction of cores survived and grew into gas giants.
Ó 2016, China University of Geosciences (Beijing) and Peking University. Production and hosting by
Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/
licenses/by-nc-nd/4.0/).

Keywords:
Planetary formation
Population synthesis
Type I migration

1. Introduction
The population of extrasolar planets (Baruteau and Papaloizou,2013), and perhaps even the solar system (Walsh et al., 2011),
provide strong evidence that migration has played a role in shaping
planetary systems. Low-mass planets (i.e., those with masses up tothat of Neptune) migrate through the excitation of linear density
waves in the disk, and through a contribution from the corotation
region (i.e., Type I migration). Early analytical work (Tanaka et al.,2002) focused on isothermal disks, in which the temperature wasprescribed and ﬁxed. These studies found that migration was always directed inward for reasonable disk parameters, and that
migration time scales were much shorter than the disk life time;
therefore, according to migration theory, all the planets should end
up very close to the central star.
While Type I migration has always been linked to linear interactions with the disk, Paardekooper and Papaloizou (2009)* Corresponding author.
E-mail address: takanori@kusastro.kyoto-u.ac.jp (T. Sasaki).
Peer-review under responsibility of China University of Geosciences (Beijing).

showed that corotation torque (or horseshoe drag) in isothermal
disks show nonlinear behavior and can be much larger than previous linear estimates, which works against fast inward migration.
However, the corotation tends to be prone to saturation and fail toprevent rapid inward migration for most of the cases, since in theabsence of a diffusive process, the corotation region is a closed
system; therefore, it can only provide a ﬁnite amount of angularmomentum to a planet.
Several well-established theoretical models of planet formation
based on the core accretion scenario adopted a population synthesis approach (e.g., Ida and Lin, 2004, 2008; Mordasini et al.,2009a,b; Ida et al., 2013). Ida and Lin (2004) focused on the inﬂuence of Type I migration on planetary formation processes andfound that when the effects of Type I migration are taken into account, planetary cores have a tendency to migrate into their hoststars before they acquire adequate mass to initiate efﬁcient gas
accretion. In order to preserve a sufﬁcient fraction of gas giants
around solar-type stars, they introduced a Type I migration
reduction factor, where factor magnitudes of smaller than unity
work to lengthen the Type I migration timescale relative to thosededuced from linear theories. With a range of small factors (w0.01),
it was possible to produce a planetary Mpa distribution that was

http://dx.doi.org/10.1016/j.gsf.2016.04.002
1674-9871/Ó 2016, China University of Geosciences (Beijing) and Peking University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NCND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

216

T. Sasaki, T. Ebisuzaki / Geoscience Frontiers 8 (2017) 215e222

qualitatively consistent with observations from a radial velocity
survey. While several suppression mechanisms for Type I migration
under various circumstances have been suggested (e.g.,
Paardekooper et al., 2011), the origin of the extremely small
reduction factor values remains unknown.
Recently, it was proposed that dynamic corotation torque canalso play a role for low-mass planets, especially where static
corotation torques saturate. Paardeokooper (2014) presented ananalysis of the torques on migrating, low-mass planets in locally
isothermal disks. They found that planets experience dynamic
torques whenever there is a radial gradient in vortensity in addition
to static torques, which do not depend on the migration rate. These
dynamic torques are proportional to the migration rate and canhave either a positive or a negative feedback on migration,
depending on whether the planet is migrating with or against thestatic corotation torque. Moreover, they showed that in disks a few
times more massive than the minimum mass solar nebula (MMSN),
the effects of dynamic torques are signiﬁcant to reduce inward
migration.
In this study, we deduced a torque formula for Type I migration
by taking into account dynamic corrections. Using this formula, we
performed population synthesis simulations with and without the
dynamic corrections in order to evaluate the migration velocity
quantitatively. We found that the effective torques with dynamic
correction were much smaller than the simple static torques whenapplied to disks of the MMSN model. We used dynamic torques
based on the theory of Paardeokooper (2014) and estimated actual
Type I migration, and then simulated various sets of planetarysystems based on the observed range of disk properties. Finally, we
compared the simulated results with observational data.

vp þ


2qd 
Gstatic þ Gdynamic ¼ 0
3pqrp US

According to Paardekooper et al. (2011; hereafter Pa11) andColeman and Nelson (2015; hereafter CN14), new static torqueformula for Type I migration can be derived (see Appendix). Whendeveloping the dynamic correction formula, we considered thework of Paardeokooper (2014), who showed that Gdynamic, the term
proportional to drp/dt, must be included in the torque formula, or inother words:

(1)

x2

2.1. Inviscid case 6rsp nvp > 1
We derived vp by substituting Eqs. 2 and 3 into Eq. 4 as:Ginviscid ¼  

(2)where S is the surface density of the disk, rp is the semimajor axis of
the protoplanet, xs is the thickness of the horseshoe region,
U ¼ ðGM =rp3 Þ1=2 (G is the gravitational constant, and M* is thestellar mass), and vp ¼ drp/dt is the radial velocity of the protoplanet. Here, 1wc/w(rp) was calculated by (Pa14 Eq. 28, modiﬁedby TE):



  
x2
1  wc w rp ¼ ð3=2 þ pÞmin 1; s vp
6rp n

(3)

where p ¼ d lnS/d ln r and n is the viscosity of the disk. Assuming acircular orbit of the protoplanet, vp can be calculated by theequation:

slib



dvp
2qd
G
¼ vp þ
þ Gdynamicdt
pqrp3 US static1G
1  ð3=2 þ pÞmc static(6)

where, mc is given by (Pa14 Eq. 20):mc ¼ 4qd xs =q(7)where xs ¼ xs =rpx2

2.2. Viscid case 6rsp nvp < 1
We derived a quadratic equation of vp by substituting Eqs. 2 and3 into Eq. 4 as:

Av2p þ Bvp þ C ¼ 0

(8)

where

A ¼





2qd x3s rp 3
sn3
þ p ¼ mc
þp
2
2
3qn
6rp

C ¼

(9)

(4)

where slib ¼ 4prp/(3Uxs) is the liberation timescale of gas in thedisk, q ¼ Mp/M*, and qd ¼ prp2 S=M . For the case of slow migration

(10)

rp
2qd q
g
r Ug
¼
smig static
ph2 p static

where

h

is

the

scale

(11)
height

of

the

disk,

gstatic ¼ Gstatic =G0 ðG0 ¼ ðq=hÞ2 Sr4 U2 Þ, and sn and smig are thetimescales of diffusion and migration, respectively, as given by(Pa14 Eqs. 10 and 23):sn ¼

rp2 x2s(12)nwhich is given by (Pa14 Eq. 18):

Gdynamic ¼ 2p 1  wc w rp Srp2 xs Uvp(5)

B ¼ 1

2. Dynamic correction of Type I migration formula

G ¼ Gstatic þ Gdynamic(i.e., slib4prp/(3Uxs)rp/vp), we were able to assume a steadystate for Eq. 4 to determine vp, or in other words:

smig ¼ph2
2qd qU(13)

The quadratic formula gives the total torque after dynamic
correction:

Gviscid ¼ QðkÞGstatic(14)

where the function Q(k) is deﬁned by (Pa14 Eq. 30):QðkÞ ¼

1

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1  2kk(15)

where k is the coefﬁcients given by (Pa14 Eqs. 31e32):

k ¼





gstatic q2d x3s rp U
8 3
3
mc sn gstatic
þp
þ
p
¼
n
3p 2
2
6smig
h2

(16)

The function Q(k) takes a critical value of 2 at k ¼ 1/2, but fork > 1/2 it does not take any value, since the inside of the square rootof Q(k) becomes negative. This suggests that runaway migration
takes place for the case k > 1/2. Paardeokooper (2014) suggested

T. Sasaki, T. Ebisuzaki / Geoscience Frontiers 8 (2017) 215e222

that the time scale of migration for runaway case would be mcsmig,
and in such a case, torque for the runaway migration would be:

Grw ¼

q
G
4qd 0

(17)

The results are consistent with the numerical results ofPaardekooper et al. (2011). In their simulation, vp rapidly converged
to the values obtained here, after a short transient phase (Figs. 7, 8,
and 9 in Paardeokooper (2014)).
In summary, our new torque formula of Type I migration, taking
into account dynamic effects, is given as:1
GI ¼ Gstatic min; QðkÞ
1  mc ð3=2 þ pÞ
¼

q
G
4qd 0

k>0:5

k < 0:5

(18)

(19)

Fig. 1 shows the dynamic correction factor at each semimajor
axis for embryos with masses of 0.01 M4, 0.1 Mp4ﬃﬃﬃﬃﬃﬃ, 1.0 M4, and 10
M4. The massespﬃﬃﬃﬃﬃﬃ
of the disks were (a) 1= 10  MMSN, (b)1.0  MMSN, (c) 10  MMSN, and (d) 10  MMSN. Except for aclose-in small protoplanet, most protoplanets had correction factors of signiﬁcantly less than 0.1; therefore, Type I migration was
generally signiﬁcantly slowed by dynamic effects.
3. Planet formation and migration modelIn our model, we adopted the models of Ida and Lin (2004, 2008)
and Ida et al. (2013) for (1) planetesimals’ growth through cohesive

217

collisions, (2) the evolution of planetesimal surface density, (3)embryos’ Type I migration and their stoppage at the disk inner edge
(except for a modiﬁcation of the Type I migration formula toinclude dynamic correction; see Section 2), and for the gas giants,(4) the onset, rate, and termination (through gap opening and/or
global depletion) of efﬁcient gas accretion, and (5) their Type II
migration.
3.1. Disk models
We adopted the MMSN model (Hayaschi, 1981) as a ﬁducial setof initial conditions for planetesimal surface density (Sd) andintroduced a multiplicative factor (fd). For the gas surface density(Sg), we adopted the r-dependence of steady accretion disk withconstant viscosity ðSg fr 1 Þ scaled by that of the MMSN at 10 AU
with a scaling factor (fg). Following Ida and Lin (2008), we set:(

Sd ¼ Sd;10 hice fd ðr=10AUÞ1:5
Sg ¼ Sg;10 fg ðr=10AUÞ1:0(20)

where normalization factors Sd,10 ¼ 0.32 g/cm2 andSg,10 ¼ 75 g/cm2, and the step function was hice ¼ 1 inside theice line at aice and 4.2 for r > aice.
Neglecting the detailed energy balance in the disk (Chiang andGoldreich, 1997), we adopted the equilibrium temperature distribution of optically thin disks given by Hayaschi (1981), such that: r 1=2  L 1=4
T ¼ 280
K
1AU
L1

(21)

Figurepﬃﬃﬃﬃﬃﬃ 1. Dynamic correction factor atpeach
ﬃﬃﬃﬃﬃﬃ semimajor axis for 0.01 M4 (M4 ¼ ME), 0.1 M4, 1.0 M4, and 10 M4 from top to bottom in each panel. The mass of the disks are (a)1= 10  MMSN, (b) 1.0  MMSN, (c) 10  MMSN MMSN, and (d) 10  MMSN.

218

T. Sasaki, T. Ebisuzaki / Geoscience Frontiers 8 (2017) 215e222

where L* and L1 are stellar and solar luminosity. We set the ice line
to be that determined by an equilibrium temperature in opticallythin disk regions:



1=2



Mp Sg r 2 r UK 2dr
xC1  1:08ðp þ 0:80q  2:52Þ
r UK
dt
M M
cs

(26)

where fg,0 is the initial value of fg and sdep is the gas depletion
timescale.

where p ¼ d log Sg/d log r, q ¼ d log T/d log r, cs is the sound speed,
and UK is the Keplerian angular velocity. The expression of Tanaka
et al. (2002) corresponds to C1 ¼ 1, and for slower migration, C1 < 1.
While we derived a new torque formula for Type I migrationthat included dynamic corrections (see Section 2), for comparison
we also used the conventional formula with the scaling factor
C1 ¼ 1.0.
We assumed that Type I migration ceases inside the innerboundary of the disk, because at this point fg is locally zero. For
computational convenience, we set the disk inner boundary to be
the edge of the magnetospheric cavity at 0.04 AU.

3.2. From oligarchic growth to isolation

3.4. Formation of gas giant planets

On the basis of the oligarchic growth model (Kokubo and Ida,
1998, 2002), the growth rate of embryos/cores at any location, a,and time t, in the presence of disk gas, was described by:Models for the formation of gas giant planets were the same as
those used in Ida et al. (2013). Embryos were surrounded bygaseous envelopes when their surface escape velocities became
larger than the sound speed of the surrounding disk gas. Whentheir mass grew (through planetesimal bombardment) above a
critical mass:

aice ¼ 2:7

L
L1

AU

(22)

Owing to viscous diffusion and photoevaporation, fg decreases
with time. For simplicity, we adopted:

fg ¼ fg;0 exp 

t

!
(23)

sdep

dMc
Mc
¼
sc;acc
dt

(24)

where
1 2=5
sc;acc ¼ 3:5  105 h1
ice fd fg

 a 5=2  M 1=3  M 1=6
c

yr
1AU
M4
M1
(25)

where Mc is the mass of the embryo/core. Furthermore, we set themass of typical ﬁeld planetesimals to be m ¼ 1020 g.
We computed the evolution of Sd distribution due to accretionby all emerging embryos in a self-consistent manner. The growthand migration of many planets were integrated simultaneouslywith the evolution of the Sd-distribution.
During the early phase of evolution, embryos are embeddedin their natal disks. Despite their mutual gravitational perturbation, embryos preserve their circular orbits owing to gravitational drag from disk gas (Ward, 1993) and dynamic friction fromresidual planetesimals (Stewart and Ida, 2000). After the disk gas
is severely depleted, the efﬁciency of the eccentricity damping
mechanism is reduced, and the embryos’ eccentricity grow until
they cross each other’s orbits (i.e., giant impact). However, in
this study, growth via the giant impact process was notconsidered. Moreover, we also ignored dynamic interaction between planets, with the growth of individual planets integrated
independently.

3.3. Type I migration
Type I migration of an embryo is caused by the sum of tidal
torque from disk regions that are both interior and exterior to theembryos. The rate and direction of embryos’ migration are determined by the differential Lindblad and corotation torques. While aconventional formula of Type I migration, which assumes locally
isothermal disks (Tanaka et al., 2002), shows that the migration isalways inward, recent developments have shown Type I migrationof isolated embryos in non-isothermal disks; therefore, themagnitude and sign of tidal torque can be changed. Ida and Lin(2008) used the conventional formula of Type I migration in
isothermal disks derived by Tanaka et al. (2002) with a scaling
factor C1 of:

M_ cMc;hydro x10
106 M4 yr1

!0:25
M4

(27)

both the radiative and convective transport of heat became sufﬁciently efﬁcient to allow their envelope to contract dynamically(Ikoma et al., 2000).
In the above equation, we neglected the dependence on opacity
in the envelope (Hori and Ikoma, 2010). In regions where cores
have already acquired isolation mass, their planetesimal-accretion
rate ðM_ c Þ would be much diminished (Ikoma et al., 2000) andMc,hydro would be comparable to an Earth-mass, M4. However, gas
accretion also releases energy and its rate is still regulated by the
efﬁciency of radiative transfer in the envelope, such that:dMp Mp
x
sKH
dt

(28)

where Mp is the planet mass including gas envelope. According to
Ida and Lin (2008), we approximated the Kelvin-Helmholtz
contraction timescale, sKH, of the envelope using:sKH x109



Mp
M4

3
(29)

Eq. 28 shows that dMp/dt rapidly increases as Mp grows; however, this is limited by the global gas accretion rate throughout thedisk and by the process of gap formation near the protoplanets’
orbits. The disk accretion rate can be expressed as:

M_ disk x3  109 fg


M1 yr1
103

a

(30)

where a is a parameter of alpha prescription for turbulent viscosity(Shakura and Sunyaev, 1973). During the advanced stage of disk
evolution, we assumed that both M_ disk and Sg evolved in proportion to exp(t/sdep). The rate of accretion onto the cores cannot
exceed M_ disk .
A gap, or at least a partial gap, is formed when a planet’s tidal
torque exceeds the disk’s intrinsic viscous stress (Lin and
Papaloizou, 1985). This viscous condition for gap formation issatisﬁed for planets with:

T. Sasaki, T. Ebisuzaki / Geoscience Frontiers 8 (2017) 215e222


Mp >Mg;vis x30



a

103

 
a 1=2 L 1=4
M4
1AU
L1

(31)

In this case, Type I migration transitions to Type II migration.
Fluid dynamic simulations (D’Angelo and Kley, 2003; D’Angelo and
Lubow, 2008) show that some fraction of gas still ﬂows into the gap.
Following the results of Dobbs-Dixon et al. (2007), we completely
terminated gas accretion when a planet’s Hill radius became larger
than two times the disk scale height, which corresponded to thethermal condition of Lin and Papaloizou (1985), that is:

Mp >Mg;th x0:95103 a 3=4  L 3=8  M 1=2


M4
1AU
L1
M1

(32)

In general, our models for gas accretion rates onto the cores
were:

dMp
¼ fgap M_ p;nogap
dt

(33)

when in the absence of any feedback on the disk structure.
Therefore, without the effect of gap opening:



Mp _
; Mdisk
M_ p;nogap ¼ min

(34)

sKH

and fgap is a reduction factor due to gap opening:

81
for Mp < Mg;vis
>
>
>
>
>
h
i
< logMp  logMg;vis
for Mg;vis < Mp < Mg;thfgap ¼ logMg;th  logMg;vis
>
>
>
>
i
h
>
:
0
for Mp >Mg;th(35)

3.5. Type II migration
During gap formation, embedded gas giants adjust their positions in the gap to establish a quasi equilibrium between the torqueapplied on them from disk regions both interior and exterior totheir orbits. Subsequently, as the disk gas undergoes viscous
diffusion, this interaction leads to Type II migration.
We assumed that planets undergo Type II migration after they
have accreted a sufﬁcient mass to satisfy the viscous condition(Mg,vis < Mp) for gap formation.
While Mp increases, the disk mass declines owing to stellar andplanetary accretion and photoevaporation. While disk mass exceeds Mp (during the disk-dominated regime), planets’ Type II
migration is locked, with the viscous diffusion of the disk gas.
During the advanced stages of disk evolution, when the mass becomes smaller than Mp (during the planet-dominated regime),
embedded planets carry a major share of the total angular momentum content.
For the disk-dominated regime, the migration timescale is givenby:smig2;disk x0:7  105



1 

a

103



a  M 1=2
yr
1AU M1(36)

For the planet-dominated regime, the migration timescale isgiven by:smig2;pl x5  105 fg1



C2 a
104

1  Mpa 1=2 M 1=2
yr
1AU
MJ
M1
(37)

219

where C2 is an efﬁciency factor associated with the degree ofasymmetry in the torques between the inner and outer disk regions. If the inner disk is severely depleted, C2 ¼ 1. We treated thefactor C2 as a model parameter, and we set C2 ¼ 0.1.
4. Population synthesis of planetary systems
Using our new torque formula for Type I migration, we modeled
the formation of planetary systems using Monte Carlo simulations.
The predicted mass and period distributions were compared withthose from a conventional Type I migration model.
4.1. Numerical settings
We ﬁrst generated a set of 1000 disks with various values of fg,0(the initial value of fg) and sdep. We adopted a range of disk model
parameters that represented the observed distribution of disk
properties and assigned them to each model with an appropriate
statistical weight. For the gaseous component, we assumed that fg,0
had a lognormal distribution centered on fg,0 ¼ 1 with a dispersionof 1 and an upper cutoff at fg,0 ¼ 30, independent of the stellar
metallicity. For heavy elements, we choose fd;0 ¼ 10½Fe=Hd fg;0 ,where [Fe/H]d is the metallicity of the disk. We assumed that thesedisks had the same metallicities as their host stars. We also
assumed that sdep had log-uniform distributions in the range
106e107 yr.
For each disk, 15 values of a for the protoplanetary seeds wereselected from a long-uniform distribution in the range 0.05e30 AU,
assuming that the mean orbital separation between planets was 0.2
on a logarithmic scale. Constant spacing in the logarithm corresponded to the spacings between the cores, which were proportional to a. This represented the simplest choice and a natural
outcome of dynamic isolation at the end of oligarchic growth.
In all simulations, the values a ¼ 103 and M* ¼ 1M1 wereassumed. Since ongoing radial velocity surveys are focused onrelatively metal-rich stars, we presented our results with [Fe/H] ¼ 0.1.
We artiﬁcially terminated Type I and Type II migration near the
disk inner edge at 0.04 AU. We did not specify a survival criterion
for the close-in planets because we lacked adequate knowledge
about planets’ migration and their interaction with host stars near
the inner edge of their nascent disks. Hence, we recorded all of theplanets that migrated to the vicinity of their host stars. In reality, alarge fraction of the giant planets that have migrated to small disk
radii were either consumed (Sandquist et al., 1998) or tidally disrupted (Trilling et al., 1998) by their host stars. Cores that migrate to
the inner edge of the disk may also coagulated and form superEarths (Ogihara and Ida, 2009); however, this was not consideredin our simulations.
4.2. Simulated individual systems
We compared the time evolution of planetary masses andsemimajor axes for the new torque formula model (Figs. 2a and 3a)and the conventional torque formula model with C1 ¼ 1.0 (Figs. 2band 3b). We choose a disk a few times more massive (fg,0 ¼ 6.0 and
8.0) than the minimum solar nebula, and with sdep ¼ 3  106 yr. The
results showed that inward migration of planet embryos was
slowed signiﬁcantly by dynamic effects. When dynamic effects
were considered, some cores survived and grew into gas giants;
however, when considering only static torque, all cores migrated to
the vicinity of their central star before growing enough to accrete
the nebula gas.
Fig. 4 shows the evolution of the dynamic correction factor (seeEq. 18),

220

T. Sasaki, T. Ebisuzaki / Geoscience Frontiers 8 (2017) 215e222

Figure 2. Growth and migration of planets for scaling factor fg,0 ¼ 3.0. Units of mass(Mp) and semimajor axis (a) are Earth masses (M4 ¼ ME) and AU. (a) Mass evolutions
obtained from simulations with the new torque formula for Type I migration. (b)Evolutions using the conventional formula (C1 ¼ 1.0).

Figure 4. Evolution of the dynamic correction factor for a scaling factor of (a) fg,0 ¼ 3.0,
and (b) fg,0 ¼ 8.0. Units of mass (Mp) are Earth masses (M4 ¼ ME).

GI =Gstatic ¼ min
1; QðkÞ
1  mc ð3=2 þ pÞ

(38)

with the mass of each planet embryo. The correction factors
remained small (0.1) throughout the simulation; therefore, thedynamic correction of Type I migration effectively prevented embryos from migrating to the central star.
4.3. Distributions of mass and semimajor axes

Figure 3. Growth and migration of planets for scaling factor fg,0 ¼ 8.0. Units of mass(Mp) and semimajor axis (a) are Earth masses (M4 ¼ ME) and AU. (a) Mass evolutions
obtained from simulations with the new torque formula for Type I migration. (b)Evolutions using the conventional formula (C1 ¼ 1.0).

We compared the predicted Mpea distributions using the new
torque formula (Fig. 5a) and the conventional formula with C1 ¼ 1.0
at t ¼ 2  107 yr. In order to directly compare the theoretical predictions with the observed data, we plotted values of Mp that were1.27 times the values of Mpsini, as determined from radial velocity
measurements (Fig. 5c). This correction factor corresponded tomean values of 1/(sini) ¼ 4/p for a sample of planetary systems withrandomly oriented orbital plants. To compare the theoretical results with M* ¼ 1M1, we plotted only the data of planets aroundstars with M* ¼ 0.81.2M1 that have been observed by radial velocity surveys.1
For the conventional models, the formation probability of gas
giants dramatically changed with C1 (Ida and Lin, 2008). Within thelimits of Type I migration with an efﬁciency comparable to that
deduced from the traditional linear torque analysis (i.e., with
C1 ¼ 1; Fig. 5b), all cores were cleared prior to gas depletion, suchthat gas giant formation was effectively suppressed. However,when considering the dynamic effects, a substantial fraction of thecores survived and grew into large gas giants (Fig. 5a). We carry out
a Kolmogorov-Smirnov (K-S) test for statistical similarity between
the predicted Mpea distributions and the observed data for the1

See http://exoplanet.eu/.

T. Sasaki, T. Ebisuzaki / Geoscience Frontiers 8 (2017) 215e222

221

contrast, simulations using a conventional formula (C1 ¼ 1.0;
Fig. 5b) cannot explain the observed data. These results show thatplanet populations consistent with observations can be reproduced
naturally (i.e., without considering the reduction factor) if we take
into account dynamic corrections for Type I migration torque.
5. Conclusions
We derived a new torque formula for Type I migration by taking
into account dynamic corrections. Using this formula, we performed
population synthesis simulations with and without the effects ofdynamic torques. In most cases, inward migration was signiﬁcantly
slowed by the dynamic effects. Considering just static torques, gas
giant formation was effectively suppressed by Type I migration ofcores; however, when dynamic effects were considered, a substantial fraction of cores survived and grew into gas giants.
Acknowledgments
We thank an anonymous reviewer for the helpful comments. This
research was supported by Grant-in-Aid for Scientiﬁc Research on
Innovative Areas from the Ministry of Education, Culture, Sports,
Science and Technology (MEXT; Grant No. 26106006). T.S. was
supported by a Grand-in-Aid for Young Scientists (KAKENHI B) fromthe Japan Society for the Promotion of Science (JSPS; Grant No.
24740120).
Appendix. Static Torque formula
Static torque (Pa11 and CN14) is given by

Gstatic ¼ FL GLR þ GVHS Fðpn ÞGðpn Þ þ GEHS Fðpn ÞF pc

rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
 ﬃ
 Gðpn ÞG pc þ GLVCT ð1  Kðpn ÞÞrﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

 ﬃ
Fe Fi
þ GLECT ð1  Kðpn ÞÞ 1  K pc
Figure 5. Planetary mass and semimajor axis distribution. Units of mass (Mp) andsemimajor axis (a) are Earth masses (M4 ¼ ME) and AU. (a) Distribution obtained fromMonte Carlo simulations with the new torque formula for Type I migration. (b) Distributions using the conventional formula (C1 ¼ 1.0). (c) Observational data of extrasolar planets around stars with M* ¼ 0.8  1.2M1 detected by radial velocity surveys.
The determined value of Mpsin i is multiplied by 1/(sin i) ¼ 4/p x 1.27, where a randomorientation of the planetary orbital planes is assumed.

(A1)
where GLR, GVHS, GEHS, GLVCT, and GLECT are the Lindblad torque,vortensity and entropy related horseshoe drag torques, and linear
vortensity and entropy related corotation torques, respectively, as
given by Eqs. 3e7 in Paardekooper et al. (2011):.

GLR ¼ ð2:5  1:7b þ 0:1aÞG0 geffparameter domain of 0.1 AU < a < 5 AU and Mp > 100M. While theconventional model produces no giant planets (Fig. 5b), the predicted Mp  a distribution using the new torque formula (Fig. 5a) is
statistically similar to the observed data (Fig. 5c) within a signiﬁcance level of p-value > 0.05 for both the semimajor axis and mass
cumulative distribution functions.
Without considering the dynamic correction for Type I migration, when a planet’s mass exceeded that of Earth, the corotation
torque became smaller owing to saturation. At this point, static
torque affected the planet more efﬁciently so that its inward
migration was rapid. However, when the dynamic correction was
included, the migration timescale was short, and the saturation of
the corotation torque was less effective. Under these conditions,inward migration slowed, which allowed for the formation of gas
giants before migration to the central star.
In summary, population synthesis simulations using our new
torque formula with dynamic correction (Fig. 5a) can explain thegas giants (>100M4) observed in exoplanetary systems (Fig. 5c). In.

GVHS ¼ ½1:1ð3=2  aÞG0 geff
 .



.

GEHS ¼ 7:9 x geff G0 geff.

GLVCT ¼ ½0:7ð3=2  aÞG0 geff
GLECT ¼hi .
.
2:2  1:4 geff G0 geff

(A2)
(A3)
(A4)
(A5)
(A6)

where a ¼ d ln S/d ln r, b ¼ d ln Tm/d ln r, and x ¼ b  (geff  1)a.
Here, G0 ¼ ðq=hÞ2 Sr 4 U2 .
The functions F(pn), F(pc), G(pn), F(pc), K(pn), and K(pc) are related
to the ratio between the viscous/thermal diffusion time scale andhorseshoe liberation/horseshoe U-turn time scales, given by Eqs.
21, 23 and 30 in Paardekooper et al. (2011):

222

T. Sasaki, T. Ebisuzaki / Geoscience Frontiers 8 (2017) 215e222

1FðpÞ ¼GðpÞ ¼KðpÞ ¼FL ¼


16 45p 3=4 3=2
p
25
8



9
8 4=3 8=31p
25 45pp<rﬃﬃﬃﬃﬃﬃﬃﬃﬃ
845p

rﬃﬃﬃﬃﬃﬃﬃﬃﬃ
8
p>
45p



16 45p 3=4 3=2
p
25 28



9
28 4=3 8=31p
25 45pp<rﬃﬃﬃﬃﬃﬃﬃﬃﬃ
28
45p

rﬃﬃﬃﬃﬃﬃﬃﬃﬃ
28
p>
45p

sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
r 2 Ux3s
2pn

sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
r 2 Ux3s
pc ¼
2pc

(A9)

C ¼

1:1

g1=4
eff



(A10)

(A11)

(A12)

(A14)


0:4 1=4
b=h

(A15)

and (Pa11 Eq. 34)

c ¼

4gðg  1ÞsT 4
3kr2 H 2 U2(A16)

The effective adiabatic index geff is given by (Pa11 Eqs. 45, 46):

geff ¼

2Q g
rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ﬃgQ þ 12 2 g2 Q 2 þ 1 2  16Q 2 ðg  1Þ þ 2g2 Q 2  2(A17)taking into account of the photon diffusion in a disk. Here,Q ¼2c3h3 r 2 U

(A18)

Fe and Fi are the reduction factor due to eccentricity and inclination of the planets, which are give by (CN15 Eqs. 16 and 20):e
Fe ¼ exp 
ef

!#
 4
 
 e  i 2 1ii
þ 0:0850:07
 0:08
h
e
h h
(A22)

Pe ¼

1þ



1=2

e
2:25h

1

(A13)

qﬃﬃﬃﬃﬃﬃﬃﬃ
q=h

Pe

jPe j

where Pe is deﬁned aswhere (Pa11 Eqs. 48, 49) xs ¼ xs =r is given byxs ¼ CPe þ(A8)The pn and pc are given by (Pa11 Eq. 19 and 40):

2
pn ¼
3

"

(A7)

1 þ ðp=1:3Þ2

!
(A19)

where e is the plane’s eccentricity and ef is deﬁned as:ef ¼ h=2 þ 0:01

(A20)

Fi ¼ 1  tanhði=hÞ

(A21)

where i is the inclination of the planet. The factor FL is the reductionin Lindblad torques when planets are on eccentric or inclined orbits, and is given by Cresswell and Nelson (2008):



6

e
þ 2:84h
4

(A23)

e
2:02h

References
Baruteau, C., Papaloizou, J.C.B., 2013. Disk-planets interactions and the diversity of
period ratios in Kepler’s multi-planetary systems. Astrophysical Journal 778 (7), 15.
Chiang, E.I., Goldreich, P., 1997. Spectral energy distributions of T Tauri stars with
passive circumstellar disks. Astrophysical Journal 490, 368e376.
Coleman, G.A.L., Nelson, R.P., 2015. On the formation of planetary systems via
oligarchic growth in thermally evolving viscous discs. Monthly Notices of theRoyal Astronomical Society 445, 479e499.
Cresswell, P., Nelson, R.P., 2008. Three-dimensional simulations of multiple protoplanets
embedded in a protostellar disc. Astronomy & Astrophysics 482, 677e690.
D’Angelo, G., Kley, W., 2003. Orbital migration and mass accretion of protoplanets inthree-dimensional global computations with nested grids. Astrophysical Journal 586, 540e561.
D’Angelo, G., Lubow, S.H., 2008. Evolution of migrating planets undergoing gas
accretion. Astrophysical Journal 685, 560e583.
Dobbs-Dixon, I., Li, S.L., Lin, D.N.C., 2007. Tidal barrier and the asymptotic mass ofproto gas-giant planets. Astrophysical Journal 660, 791e806.
Hayaschi, C., 1981. Structure of the solar nebula. Growth and decay of magnetic
ﬁelds, and effect of magnetic and turbulent viscosities on the nebula. Progressof Theoretical Physics Supplement 70, 35e53.
Hori, Y., Ikoma, M., 2010. Critical core masses for gas giant formation with grain-freeenvelopes. Astrophysical Journal 714, 1343e1346.
Ida, S., Lin, D.N.C., 2004. Toward a deterministic model of planetary formation. I. A
desert in the mass and semimajor axis distributions of extrasolar planets.
Astrophysical Journal 604, 388e413.
Ida, S., Lin, D.N.C., 2008. Toward a deterministic model of planetary formation. IV.
Effects of type I migration. Astrophysical Journal 673, 487e501.
Ida, S., Lin, D.N.C., Nagasawa, M., 2013. Toward a deterministic model of planetary formation. VII. Eccentricity distribution of gas giants. Astrophysical Journal 775 (42), 25.
Ikoma, M., Nakazawa, K., Emori, H., 2000. Formation of giant planets: dependenceson core accretion rate and grain opacity. Astrophysical Journal 537, 1013e1025.
Kokubo, E., Ida, S., 1998. Oligarchic growth of protoplanets. Icarus 131, 171e178.
Kokubo, E., Ida, S., 2002. Formation of protoplanet systems and diversity of planetary systems. Astrophysical Journal 581, 666e680.
Lin, D.N.C., Papaloizou, J.C.B., 1985. On the dynamical origin of the solar system. In:
Black, D.C., Matthew, M.S. (Eds.), Protostars and Planets II. Univ. Arizona Press,
Tucson, pp. 981e1072.
Mordasini, C., Alibert, Y., Benz, W., 2009a. Extrasolar planet population synthesis. I.
Method, formation tracks, and mass-distance distribution. Astronomy & Astrophysics 501, 1139e1160.
Mordasini, C., Alibert, Y., Benz, W., Naef, D., 2009b. Extrasolar planet population
synthesis. II. Statistical comparison with observations. Astronomy & Astrophysics 501, 1161e1184.
Ogihara, M., Ida, S., 2009. N-body simulations of planetary accretion around M
dwarf stars. Astrophysical Journal 699, 824e838.
Paardekooper, S.-J., Papaloizou, J.C.B., 2009. On corotation torques, horseshoe dragand the possibility of sustained stalled or outward protoplanetary migration.
Monthly Notices of the Royal Astronomical Society 394, 2283e2296.
Paardekooper, S.-J., Baruteau, C., Kley, W., 2011. A torque formula for non-isothermal
type I planetary migration e II. Unsaturated horseshoe drag. Monthly Notices of
the Royal Astronomical Society 410, 293e303.
Paardeokooper, S.-J., 2014. Dynamical corotation torque on low-mass planets.
Monthly Notices of the Royal Astronomical Society 444, 2031e2042.
Sandquist, E., Taam, R.E., Lin, D.N.C., Burkert, A., 1998. Planet consumption andstellar metallicity enhancements. Astrophysical Journal 506, L65eL68.
Shakura, N.I., Sunyaev, R.A., 1973. Black holes in binary systems. Observational
appearance. Astronomy & Astrophysics 24, 337e355.
Stewart, G.R., Ida, S., 2000. Velocity evolution of planetesimals: uniﬁeld analyticalformulas and comparisons with N-body simulations. Icarus 143, 28e44.
Tanaka, H., Takeuchi, T., Ward, W.R., 2002. Three-dimensional interaction between aplanet and an isothermal gaseous disk. I. Corotation and Lindblad torques andplanet migration. Astrophysical Journal 565, 1257e1274.
Trilling, D.E., Benz, W., Guillot, T., Lunine, J.I., Hubbard, W.B., Burrows, A., 1998.
Orbital evolution and migration of giant planets: modeling extrasolar planets.
Astrophysical Journal 500, 428e439.
Walsh, K.J., Morbidelli, A., Raymond, S.N., O’Brien, D.P., Mandel, A.M., 2011. A low
mass for Mars from Jupiter’s early gas-driven migration. Nature 475, 206e209.
Ward, W.R., 1993. Density waves in the solar nebula: planetesimal velocities. Icarus
106, 274e287.