DYNAMICS AND ECCENTRICITY FORMATION OF PLANETS IN OGLE-06-109L SYSTEM

We assume that the masses and orbital elements of the two planets obey normal distributions,

$$\begin{equation} P(x)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left[-\frac{(x-x_0)^2}{2\sigma ^2}\right], \end{equation} \tag{ 1 }$$

with x₀ the nominal value of x and σ the observational error. The distribution of eccentricity is obtained by $e=\sqrt{\vphantom{A^A}\smash{\hbox{$x_1^2+x_2^2$}}}$, where x₁, x₂ are two independent Gaussian (1) with x₀ = 0 and different dispersions σ₁, σ₂, and it is the standard Rayleigh distribution (Zhou et al. 2007) when σ₁ = σ₂. We choose values of σ₁ and σ₂ so that the most probable eccentricity of m_c is 0.11, with 1σ confidence interval (0.07, 0.28). Also we assume that e_b follows a Rayleigh distribution with the most probable value $({\bar{e}}_b)$. The distributions of the inclinations are the same with those of e/2. The remaining angles are randomly generated.

We carry out 15,000 runs of simulations by integrating the full motion of the three bodies (m_*, m_b, m_c) in the three-dimensional physical space up to 0.1 Myr. By checking whether any one of the six resonant angles (3λ_c − λ_b − iϖ_b − jϖ_c − kΩ_b − lΩ_c with i, j, k, l non-negative integers and i + j + k + l = 2) librates, we find the probabilities that m_b, m_c in 3:1 MMR are 0.82%, 2.52%, 1.88% for ${\bar{e}}_b=0.04, 0.1, 0.2$, respectively. Thus, the probability for m_b and m_c in the 3:1 MMR is small within observational errors.

To investigate the secular dynamics of the two planets, we adopt the method of representative plane of initial conditions (Michtchenko & Malhotra 2004). Due to the existence of four center-of-mass integrals, the planar three-body (the star and two planets) system is a Hamiltonian one with four degrees of freedom, with the Hamiltonian function (e.g., Laskar & Robutel 1995)

$$\begin{equation} H=-\sum _{i=b}^c\frac{\mu _i^2m^{\prime 3}_i}{2L_i^2}-G\frac{m_b m_c}{|{\bf r}_b-{\bf r}_c|} + \frac{m_bm_c}{m_\ast }(\dot{x}_b\dot{x}_c+\dot{y}_b\dot{y}_c), \end{equation} \tag{ 2 }$$

where μ_i = G(m_* + m_i), m'_i = m_im_*/(m_* + m_i), $L_i =m^{\prime }_i\sqrt{\mu _ia_i}$, with a_i, e_i, r_i being the semimajor axis, eccentricity, relative position vector of m_b or m_c, respectively. An averaged system is obtained by eliminating the short periodic terms

$$\begin{equation} H_{\rm {sec}}=-\frac{1}{(2\pi)^2}\int _0^{2\pi }\int _0^{2\pi } Hd\lambda _b\lambda _c, \end{equation} \tag{ 3 }$$

where λ_b and λ_c are the longitudes of mean motion of m_b and m_c, respectively. Due to the D'Alembert's rule, only Δϖ_bc = ϖ_c − ϖ_b appears in the average Hamiltonian, thus the secular system is integrable with one degree of freedom in Δϖ_bc and its conjugate momentum $K_b= L_b(1-\sqrt{\vphantom{A^A}\smash{\hbox{$1-e_b^2$}}})$, and Δϖ_bc either circulates or liberates around 0 or π. So the information of secular motion with different initial conditions can be presented in the representative plane (e_bcos Δϖ_bc, e_c), where Δϖ_bc is fixed at either 0 or π.

Figure 1 shows the various motions in the representative plane of the OGLE-06-109L system, with the semimajor axes of the planets being the nominal values (Table 1). Two families of the equilibriums of the secular system (3), with Δϖ_bc(t) ≡ 0 or π are plotted as dashed lines. Nearby orbits are with Δϖ_bc(t) librating around 0 or π. Between the two libration regions are orbits with Δϖ_bc-circulating, bounded by thin solid lines. We also plot the orbital crossing time (T_c) of orbits originating from the representative plane by full three-body integrations (with other angles randomly chosen). T_c is defined as the minimum time that either a_b > a_c occurs or the separation of m_b and m_c is smaller than 1 mutual Hill radii. From Figure 1, we see that the apsidal alignment (Δϖ_bc ≈ 0) between m_b and m_c can stabilize the interacting planets.

Zoom In Zoom Out Reset image size

Some typical motions originating from Figure 1 are plotted in Figure 2. Orbits (a) and (b) are from the libration region so that Δϖ_bc librates around 0 or π, respectively. Orbits (c) are with Δϖ_bc circulating and e_c is in a near-separatrix motion, i.e., oscillates between 0 and a finite value (∼0.1) periodically, which is very similar to the Jupiter–Saturn system in panel (d).

Zoom In Zoom Out Reset image size

We integrate the three-body system with a second-order WHM code (Wisdom & Holman 1991) from the SWIFT package (Levison & Duncan 1994). The time step is set as the 1/12 of the period of the innermost orbit. Taking initial a_b and a_c with present nominal values, we carry out 2500 runs of integrations with different initial eccentricities. All the angles including ϖ_b and ϖ_c are randomly chosen. The orbital crossing time T_c, defined as the minimum time that a_b > a_c or distance of m_b and m_c being smaller than one mutual Hill radius, is presented in Figure 3. According to the results, the two-planet (m_b and m_c) system is stable (T_c > 10⁸ yr) only if

$$\begin{equation} e_b^2+e_c^2<0.3^2 \end{equation} \tag{ 4 }$$

holds approximately. Let us compare this to another commonly used stability. Hill stability requires that the ordering of the two planets remains unchanged for all the time, which allows the outer planet escaping to infinity. Topological studies of three-body systems give the following sufficient criterion for two planets being Hill stable (Marchal & Bozis 1982; Gladman 1993),

$$\begin{equation} -\left(\frac{2m_{\rm {tot}}}{G^2m^3_{\rm {pair}}}\right)L^2E> 1+\frac{3^{4/3}m_bm_c}{m_\ast ^{2/3}(m_b+m_c)^{4/3}}+\cdots, \end{equation} \tag{ 5 }$$

where m_tot = m_* + m_b + m_c, m_pair = m_*m_b + m_*m_c + m_bm_c, and L and E are the total angular momentum and energy of the three-body system, respectively. This criterion is also plotted in Figure 3. Note that the Hill stability criterion gives a larger region as that of T_c > 10⁸ yr. This implies that the stability defined by the orbital crossing time T_c > 10⁸ yr is stronger than the Hill stability, in the sense that our stability in terms of orbital crossing requires the mutual distance of two planets to be more than 1 mutual Hill radius all the time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We numerically integrate the orbits of a few hundred test particles in the planetary system. Such simulations enable us to identify regions where low-mass companions can have stable orbits (Rivera & Haghighipour 2007; Haghighipour 2008). The test particles are initially located in circular orbits coplanar with two planets, with their mean longitudes randomly set. The orbital evolution time is set as 100 Myr. A particle is removed from the simulation when its stellar distance exceeds 30 AU, or enters the Hill sphere of either planet. We study two cases that either m_b and m_c are in 3:1 MMR or not. In the non-resonance cases, m_b and m_c are initially located at the nominal elements (Table 1) with e_b = 0.06 (a value inferred from formation scenario in the following section). To show the dependence of stability on Δϖ_bc, we initially let ϖ_b = 0, M_b = 0, and ϖ_c = 0, π/2, π. M_c is randomly set. In the 3:1 MMR cases, m_b, m_c are set with nominal parameters except a_b = 2.11445 AU, e_b = 0.06, ϖ_b = ϖ_c = M_b = 0, M_c = π, so that m_b and m_c are initially in 3:1 MMR, and the three corresponding resonance angles, 3λ_c − λ_b − 2ϖ_b, 3λ_c − λ_b − ϖ_b − ϖ_c, 3λ_c − λ_b − 2ϖ_c, liberate around π, 0, π, respectively.

Inner region. Figure 4(a) shows the survival time of test particles in the inner system with different initial semimajor axes. When m_b and m_c are not in 3:1 MMR, test particles with initial a ⩽ 1.5 AU are stable in the sense of T_c > 100 Myr, except at 5:1, 3:1, 5:2, 2:1 MMRs with m_b. The locations of these MMRs depend on Δϖ_bc. However, if m_b and m_c are in 3:1 MMR, the stable region is reduced a bit to a ⩽ 1.4 AU.

Outer region. Figure 4(b) shows the survival time of test particles in the outer system. When m_b and m_c are not in 3:1 MMR, the test particles are stable as long as a ⩾ 9.7 AU. However, if the two planets are in 3:1 MMR, the stable region is enlarged to a ⩾ 7.5 AU except the 1:3 MMR with m_c. It is interesting to note that the 3:1 MMR between m_b and m_c increases the stable region in the outer regions.

Habitable zone. Kasting et al. (1993) estimated the width of the habitable zone (HZ), where an Earth-like planet can have liquid water on its surface, around main-sequence stars. For the star with mass 0.5 M_☉, the most conservative position of HZ is 0.25–0.36 AU, and the actual HZ could be much wider than this extension, e.g., the outer edge of HZ could be 0.36–0.47 AU. From Figure 4(a), if the two giant planets are out of 3:1 MMR, test particles with initial semimajor axes in HZ are stable, although their eccentricities will be excited by secular resonance due to m_b and m_c (Malhotra & Minton 2008). To see whether the planet can maintain its orbit in the HZ, we put a planet with Earth mass in different initial locations. Figure 5 shows the largest and smallest stellar distances (Q and q) of the Earth-mass planets. The shaded area is the conservative estimation of the HZ. The great variation of Q and q at a = 0.28–0.32 AU is due to the secular resonance of m_b and m_c (Malhotra & Minton 2008). Interesting to note that the variations of Q and q are small when Δϖ_bc = π, i.e., the eccentricity excited by the secular resonance of m_b, m_c is small (Migaszewski et al. 2009). Also, when m_b, m_c is in 3:1 MMR, HZ is not in their secular resonance region, so the variations of Q and q are also small. No matter in which type of Δϖ_bc, most of orbits are in HZ except a ∈ [0.28 AU, 0.32 AU] with Δϖ_bc = 0, π/2. Considering the actual HZ is wider than this, OGLE-06-109L is a hopeful candidate system for hosting a habitable terrestrial planet. Simulations in the following section indeed show evidence for the formation of super-Earth planets in its HZ.

Zoom In Zoom Out Reset image size

According to the conventional core-accretion scenario of planet formation, the planet is formed through planetesimals coagulation by means of runaway growth and became protoplanetary embryos by oligarchic growth in the protoplanetary disk (Safronov 1969; Kokubo & Ida 1998). To model the masses of embryos formed in the disk, we adopt the empirical minimum-mass solar nebula (MMSN; Hayashi 1981) so that the surface density of gas disk at stellar distance a is given as

$$\begin{equation} \Sigma _{g}=2.4\times 10^3 f _g f_{\rm dep}\left(\frac{a}{\rm 1\,AU}\right)^{-3/2}\;{\rm g\;cm}^{-2}, \end{equation} \tag{ 6 }$$

where f_g is the gas enhancement factor, f_dep = exp(−t/τ_dep) is the gas depletion factor due to disk accretion, photoevaporation or planet formation with a timescale of τ_dep ∼ 3 Myr (Haisch et al. 2001), and t is the evolution time. The surface density of solid disk is given as

$$\begin{equation} \Sigma _{d}=10 f_d \gamma _{\rm ice} \left(\frac{a}{\rm 1\,AU}\right)^{-3/2}\;{\rm g\;cm}^{-2}, \end{equation} \tag{ 7 }$$

where f_d is the solid enhancement factor, γ_ice is the volatile enhancement with a value of 4.2 or 1 for material exterior or interior to the snow line (0.68 AU for the OGLE-06-109L system), respectively. In such a disk, the embryos will grow under cohesive collisions in a timescale of (Kokubo & Ida 2002; Ida & Lin 2004)

$$\begin{eqnarray} \tau _{\rm acc}&\simeq & 1.6\times 10^5\gamma _{\rm ice}^{-1}f_d^{-1}f_g^{-2/5}\nonumber\\ &&\times\left(\frac{a}{\rm 1\,AU}\right)^{27/10} \left(\frac{M_c}{M_{\oplus }}\right)^{1/3}\left(\frac{M_*}{M_{\odot }}\right)^{-1/6} \;{\rm yr}, \end{eqnarray} \tag{ 8 }$$

where M_c is the core mass. The core growth will continue until it accretes all the dust material round its feeding zone (∼10 Hill radii) so that an isolation body is achieved with a mass of (Ida & Lin 2004)

$$\begin{equation} M_{\rm iso} = 0.12 \gamma _{\rm ice}^{3/2} f_d^{3/2} \left(\frac{a}{\rm 1\,AU}\right)^{3/4} \left(\frac{M_*}{M_\odot }\right)^{-1/2}\;M_\oplus. \end{equation} \tag{ 9 }$$

For the OGLE-06-109L system, the isolation mass at 4 AU with f_d = 2 is around 12 M_⊕, above the critical mass (∼10 M_⊕) for the onset of efficient gas accretion to form giant planets (Pollack et al. 1996), and the core growth timescale is ∼1.2 Myr.

Interactions between embryos and the gas disk may damp the eccentricities of the embryos (Goldreich & Tremaine 1980). The timescale of the eccentricity damping for an embryo with mass m can be described as (Cresswell & Nelson 2006)

$$\begin{eqnarray} \left(\frac{e}{\dot{e}}\right)_{\rm emb} &=& \frac{Q_e}{0.78}\left(\frac{M_*}{m}\right) \left(\frac{M_*}{a^2\Sigma _g}\right)\left(\frac{h}{r}\right)^4\nonumber\\ &&\!\! \times\Omega ^{-1}\left[1+\frac{1}{4}\left(e\frac{r}{h}\right)^3\right]{\rm yr }, \end{eqnarray} \tag{ 10 }$$

where r, e, h, and Ω are the stellar distance, eccentricity of the embryo, scale height of the disk, and the Kepler angular velocity, respectively, Q_e = 0.1 is a normalization factor to fit with hydrodynamical simulations.

As an embryo grows to a massive planet (⩾30 M_⊕), it will induce strong tidal torques on the disk to open a gap around it (Lin & Papaloizou 1993). Then the planet will be embedded in the viscous disk to undergo type II migration. The timescale of type II migration for a planet with mass m_p can be modeled as (Ida & Lin 2004)

$$\begin{eqnarray} \tau _{\rm II}=\frac{a}{|\dot{a}|}&=&0.8 \times 10^6 \;{\rm yr}\; f_g^{-1} \left(\frac{m_p}{M_{\rm Jup}}\right)\left(\frac{M_*}{M_{\odot }}\right)^{-1/2}\nonumber\\ &&\times\left(\frac{\alpha }{10^{-4}}\right)^{-1}\left(\frac{a}{\rm 1\,AU}\right)^{1/2}, \end{eqnarray} \tag{ 11 }$$

where α is a dimensionless parameter to adjust the effective viscosity, and we set α = 10⁻⁴ as a standard value in our simulation. When a giant planet is embedded in a gas disk, tidal interaction of disk may damp its eccentricity if it is not massive enough. As we mentioned in the abstract, the situation is quite elusive for different mass regime of the giant planets, so we adopt an empirical formula (Lee & Peale 2002)

$$\begin{equation} \left(\frac{\dot{e}}{e}\right)_{\rm pl}=-K \left| \frac{\dot{a}}{a} \right| \end{equation} \tag{ 12 }$$

to describe the eccentricity-damping rate of giant planets in the OGLE-06-109L system, where K is a positive constant with a value ranging 10–100 (Sándor et al. 2007). After some tests, we choose K = 10 in this paper to let e_b and e_c have reasonable convergent values.

In this section, we simulate the configuration formation for the OGLE-06-109L system with N-body models. We assume that the two giant planets have formed with the observed masses at 4 AU and 8–9 AU, respectively. The physical epoch corresponds to this assumption is ∼1 Myr after the formation of star, so that the gas disk is still present. The two giant planets have opened gaps around them and will undergo type II migration according to Equation (11) in the viscous disk. At the initial stage of our simulation, there are some leftover embryos in inner orbits that have obtained their isolation masses. To mimic the formation of Earth-like planets in inner orbits, we put 18 embryos, with masses ranging from 0.17 M_⊕ to 9 M_⊕ derived from Equation (9) and initial locations from 0.25 AU to 3 AU. The mutual distances among the embryos are set as 10 Hill radii. An additional embryo with in situ isolation mass will be put between m_b and m_c in model 2. All the planets and embryos are initially located in near-coplanar and near-circular orbits (e = 10⁻³, inclination i = e/2), their phase angles (mean motion, longitude of perihelion, longitude of ascending node) are randomly chosen. The acceleration of the planet (embryo) with mass m_i is given as

$$\begin{eqnarray} \frac{d}{dt}\textbf {V}_i &=& {-}\frac{G(M_*+m_i)}{{r_i}^2}\left(\frac{\textbf {r}_i}{r_i}\right) +\sum _{j\ne i}^N {Gm_j}\left[\frac{(\textbf {r}_j-\textbf {r}_i)}{|\textbf {r}_j-\textbf {r}_i|^3}-\frac{\textbf {r}_j}{r^3_j}\right]\nonumber\\ && +\,\textbf {F}_{\rm edamp}(+\textbf {F}_{\rm migII}), \end{eqnarray} \tag{ 13 }$$

where r_i and V_i are the position and velocity vectors of m_i in the stellar-centric coordinates,

$$\begin{equation} \textbf {F}_{\rm edamp} = -2\frac{(\textbf {v}_i \cdot \textbf {r}_i)\textbf {r}_i}{r_i^2\tau _e} \end{equation} \tag{ 14 }$$

is the damping acceleration, effective for all embryos and gas giants but with different τ_e in Equations (10) and (12), respectively, and the acceleration that causes the type II migration,

$$\begin{equation} \textbf {F}_{\rm migII} = -\frac{{\bf V}_i}{2 \tau _{\rm II}} \end{equation} \tag{ 15 }$$

is adopted for two giant planets. We numerically integrate the evolution of Equation (13) with a time-symmetric Hermit scheme (Aarseth 2003). The simulation is performed up to 10 Myr. As we assume that the gas disk depletes exponentially in a timescale τ_dep = 1 Myr, the gas disk almost disappears at the end of simulation.

During the earlier stage when gas disk is present (the coming models 1–2), embryos in outer disk will also induce a damping of giant planets' eccentricities through dynamic friction, which has similar effect by the gas disk. However, as the gas disk dominates before gas depletion, we did not consider the presence of embryos in the outer disk in models 1–2, except in model 3 where embryos in the outer disk are included, after the gas disk depletes.

Model 1: smooth and convergent migration. In this model, we vary the initial locations of two giant planets, a_b and a_c. The 18 embryos with in situ isolation masses in inner orbits are put so that the outermost one is in an orbit 3.5 Hill Radii away from m_b. During the evolution of the typical run R1 (see Table 2 for initial parameters), m_b and m_c are captured into 3:1 MMR at t ≈ 1 Myr, with three resonant angles liberate around either 0 or π with amplitudes around ∼0.4π (Figure 6). Their eccentricities are excited to about 0.1 inside the resonance. Eccentricities of the embryos in the inner orbits are also excited due to secular perturbations from two giant planets, which results in their inward migration in the gas disk. At the end of simulation, they merge into four planets at [0.22 AU, 0.77 AU], with masses of 4.86 M_⊕, 8.76 M_⊕, 5.79 M_⊕, and 16.06 M_⊕ in the order of increasing semimajor axes. Noticeably, the inner two are in the edge of the HZ ([0.25 AU, 0.36 AU]) of the system.

Zoom In Zoom Out Reset image size

The mechanism that the migration of giant planets triggers the merge of inner embryos in MMRs or secular resonances had been already discussed in literature, e.g., Zhou et al. (2005) and Fogg & Nelson (2005). However, in this case, the eccentricities of embryos are excited by the secular perturbations of outside giant planets. The configuration of 3:1 MMR between two planets are kept to the end of the simulation, when they are stalled at 2.53 AU and 5.26 AU due to the severe depletion of gas disk. We did eight runs in this model with a_b varying in [3.5 AU, 3.8 AU], and a_c in [8.2 AU, 8.5 AU]. In all the simulations, the two giant planets are trapped in 3:1 MMR, with their eccentricities osculating around ∼0.1 at the end of evolution as in Figure 6(c).

Model 2: planetary scattering during migration. Unlike the previous model that m_b and m_c undergo smooth migration, now, besides the 18 embryos in inner orbits, we put an additional embryo with local isolation mass (denoted as m_e1, slightly varies for different locations) between the orbits of m_b and m_c. Figure 7 shows the results of a typical run (R2), with m_e1 = 16.8 M_⊕. At t ≈ 5.2 Myr, a close encounter between m_c and m_e1 occurred, which scatters m_e1 out of the system (Figure 7). The encounter excites e_c up to 0.1, which in turn excites e_b from 0 to ∼0.08, resulting in a configuration that e_b passing 0 in the eccentricity plane (e_be_ccos Δϖ_bc, e_be_csin Δϖ_bc) almost periodically, the so-called near-separatrix (of libration and circulation of Δϖ_bc) motion (Barnes & Greenberg 2006, 2008). The eccentricities of embryos in inner orbits are also excited due to the sudden increase of e_b and e_c, which cause strong mergers among the embryos into a planet of 14.38 M_⊕ at 0.89 AU.

Zoom In Zoom Out Reset image size

If close encounters between m_e1 and one of the giant planets occur much earlier, tidal interaction between the planets (embryos) and the gas disk will eventually damp their eccentricities. In run R3, the scattering process occurred at t = 1.5 Myr when the eccentricity damping induced by the gas disk is still strong (Figure 8). As a result, e_b and e_c that excited during planetary scattering are damped to less than 0.01 quickly. In this case, 11 planets are left with their masses from 0.21 M_⊕ to 16.81 M_⊕ at [0.15 AU, 1.10 AU]. Among them, two embryos are in the HZ.

Zoom In Zoom Out Reset image size

We did 20 runs of simulations with different initial locations of m_e1. Among them, six runs have close encounter events during the evolution, like run R2, resulting in similar configurations of m_b and m_c as in Figure 7. Five runs have earlier close encounters so that e_b and e_c are damped in the end of simulations. The rest nine runs do not suffer close encounters, but have some milder encounter events occasionally. The resulting e_b and e_c have only slight changes up to ∼0.01.

In a compact configuration with a_b/a_c > 0.48, semimajor axes shift to a_b/a_c < 0.48 by planetary scattering is a possible routine that could lead to the trap of 3:1 MMR between m_b and m_c. To investigate this probability, we perform additional 900 runs of simulations with a simplified four-body (star–two giants–one embryo) model. In this model, m_b and m_c are located initially at 3.7 AU and 7.6 AU so that (a_b/a_c > 0.48), with one embryo m_e1 between their orbits. The initial semimajor axis of m_e1 is chosen at [5.0 AU, 5.4 AU]. The mass of m_e1 slightly varies at different locations, in the range of 13.8 M_⊕ − 14.6 M_⊕. Both F_edamp and F_migII in Equations (14) and (15) are included. Numerical simulations show that 185 runs out of 900 ones with close encounter events occurred between m_e1 and one of the planets. The epoch for close encounters occurred do not show a clear correlation with the relative distance between the embryo and one of the giant planets (Figure 9). Among the 185 runs with scattering events, 21 runs (2.3% of 900 runs) lead to the trap of m_b and m_c into 3:1 MMR, 11 runs (1.2%) result in the trap near the boundary of 3:1 MMR, the rest 153 runs (17.0%) do not lead to the trap of 3:1 MMR. We also observe 62 runs with m_b and m_c being trapped in 8:3 MMR. Figure 10 shows the final e_b, e_c at the end of 900 runs' simulations. As we can see, besides those being trapped into 8:3 MMR, e_b and e_c are excited significantly only in those 185 runs with planetary scattering, with the average values e_b ∼ 0.058 and e_c ∼ 0.085 (Figure 10).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Model 3: divergent migration in the presence of planetesimal disk. To study the effect of planetesimal disk in outside orbits after the gas disk is depleted, we perform simulations by including the embryos in outside orbits but discarding those in inner orbits, since they may be in hot orbits and have less affections on the outer system. We set two types of embryos in the outer region, those with masses of 5 M_⊕ and of 0.2 M_⊕. After some test simulations, we find that the solid disk out of 10 AU has little effect on the evolutions of the two giant planets at nominal location (see Figure 11 for a typical run). So we set the outer edge of solid disk within 10 AU in the following simulations. Two groups of simulations are made according to different initial a_b and a_c.

Zoom In Zoom Out Reset image size

Group 3a. We initially put m_b and m_c at a_b = 3 AU and a_c = 4 AU, so m_c is inside the 2:1 MMR location (at 4.76 AU) of m_b. The separation is about 3.3 times of their mutual Hill's radii (R_Hbc), above the threshold (2.4R_Hbc), so they are Hill stable (Gladman 1993) if there are no other perturbations. We put 63 embryos evenly at [5.5 AU, 9.5 AU], including 5 × 5 M_⊕ and 58 × 0.2 M_⊕ ones, corresponding to a solid disk of f_d = 2 with a total mass of 36.6 M_⊕ in [5 AU, 10 AU].

Figure 12 shows the evolution of m_b and m_c in a typical run, with the innermost embryos a_in = 5.5 AU. Under the perturbation of outer embryos, m_b (m_c) becomes unstable and undergoes inward (outward, resp.) migration quite soon (Figures 12(a) and (b)). The migration results in the quick crossing of 5:3 and 2:1 MMR at t ≈ 0.05 Myr, and 0.12 Myr (see Figure 12(c)), the two most strong resonances between 3:2 and 3:1 MMRs, and the system seems to be stable until m_b (m_c) reaches around 2.82 AU (4.90 AU, resp.), with a drift extension of 6% (22%, resp.). Their eccentricities are excited after the crossing of 5:3 and 2:1 MMR, with maximum values of ∼0.2–0.3 (Figure 12(d)). Finally, e_b and e_c oscillate in [0.03, 0.15] and [0.04, 0.19], respectively. The 5:3 and 2:1 MMR crossings of m_b and m_c lead to the strong scattering of embryos in outer orbits (Figure 13(a)), which results in the escape of all the embryos except one with mass 5 M_⊕ at the orbit of 15 AU and eccentricity of 0.45.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We did 15 runs of simulations in this group, by changing the initial locations of the 63 embryos so that a_in of the innermost one varies in [5.2 AU, 8 AU]. 5:3 and 2:1 MMR crossings occurred in three runs, with eccentricities e_b ∈ [0.03, 0.24] and e_c ∈ [0.04, 0.28] at the end of our simulations (t = 10 Myr). We did not observe resonance crossings in another three runs up to 10 Myr evolution, with final eccentricities e_b ∈ [0, 0.05] and e_c ∈ [0, 0.06]. In the rest nine runs, m_b and m_c have strong close encounter so that m_c is scattered out of the system, with the eccentricity of the only survival giant planet m_b ∼ 0.3. The values of a_in corresponding to these three types of outcomes (MMR crossing, non-MMR crossing, m_c being ejected) do not show clear correlation, which indicates the chaotic states of m_b and m_c under the perturbation of outer embryos.

Group 3b: m_b, m_c are initially put at a_b = 2.3 AU and a_c = 4.6 AU, the nominal locations. We put 56 embryos evenly at [6.4 AU, 10.4 AU], including 5 × 5 M_⊕ and 51 × 0.2 M_⊕ ones, corresponding to a solid disk of f_d = 2 with a total mass of 35.2 M_⊕ in [5.5 AU, 10.5 AU]. Figure 14 shows the evolution of m_b and m_c in a typical run, with the innermost embryos a_in = 6.4 AU. m_b, m_c cross 3:1 MMR during the divergent migration at t ≈ 0.06 Myr and t ≈ 2 Myr, causing the increase of e_b and e_c, At the end of simulation, m_b and m_c locate at 2.27 AU and 4.71 AU, with e_b ∈ [0.01–0.07] and e_c ∈ [0.03–0.1]. One 5 M_⊕ embryo (at 11.4 AU with e = 0.024) and 4 × 0.2 M_⊕ ones (at [16 AU, 22 AU] and e ∈[0.08, 0.53]) are left (Figure 13(b)).

Zoom In Zoom Out Reset image size

We did a total of 41 runs (13 runs with f_d = 2 and 28 runs with lower f_d) in this group. For the 13 runs of f_d = 2, five out of these eight runs with a_in < 7.8 AU are observed having m_b, m_c's 3:1 MMR crossing, with final eccentricities e_b ∈ [0, 0.08], e_c ∈ [0, 0.1]. The rest five runs with a_in > 7.8 AU do not have 3:1 MMR crossing, with final eccentricities e_b ∈ [0, 0.04] and e_c ∈ [0, 0.04]. For the 28 runs with lower solid disks (f_d = 1, 1.2, 1.5, 1.6, 1.8) and different a_in, no 3:1 MMR crossing is observed in runs with f_d < 1.6, maybe due to the small mass of solid disk (M < 22 M_⊕). While for f_d = 1.6, 1.8, the probability of 3:1 MMR crossing is ∼33%. The extensions of e_b and e_c are similar to that of f_d = 2.