Eccentricity Distribution of Gas Giants

Toward a Deterministic Model of Planetary Formation VII: Eccentricity Distribution of Gas Giants


The ubiquity of planets and diversity of planetary systems reveal planet formation encompass many complex and competing processes. In this series of papers, we develop and upgrade a population synthesis model as a tool to identify the dominant physical effects and to calibrate the range of physical conditions. Recent planet searches leads to the discovery of many multiple-planet systems. Any theoretical models of their origins must take into account dynamical interaction between emerging protoplanets. Here, we introduce a prescription to approximate the close encounters between multiple planets. We apply this method to simulate the growth, migration, and dynamical interaction of planetary systems. Our models show that in relatively massive disks, several gas giants and rocky/icy planets emerge, migrate, and undergo dynamical instability. Secular perturbation between planets leads to orbital crossings, eccentricity excitation, and planetary ejection. In disks with modest masses, two or less gas giants form with multiple super-Earths. Orbital stability in these systems is generally maintained and they retain the kinematic structure after gas in their natal disks is depleted. These results reproduce the observed planetary mass-eccentricity and semimajor axis-eccentricity correlations. They also suggest that emerging gas giants can scatter residual cores to the outer disk regions. Subsequent in situ gas accretion onto these cores can lead to the formation of distant (AU) gas giants with nearly circular orbits.

planetary systems: formation – solar system: formation – stars: statics

1 Introduction

The accelerating pace of observational discovery (with radial velocity, transit, microlensing, and direct imaging surveys) provide rich data sets on the statistical properties of extra solar planets and planetary systems. In order to extract, from these data, useful information on the dominant physical processes and appropriate range of physical condition associated with planetary formation and evolution, population synthesis models have been developed. In a series of papers (Ida & Lin, 2004a, b, 2005, 2008a, 2008b, referred to as Paper I-V), we and others (Mordasini et al., 2009a, b; Alibert et al., 2011) have outlined prescriptions for various physical processes and described attempts to reproduce observational data with simulated models. Some of the relevant effects have been taken into account so far include the coagulation of planetesimals into embryos, gas accretion onto sufficiently massive cores, type I and II migration of individual protoplanets due to their interaction with their natal disks.

Recent transit search by the Kepler space telescope indicates that many, if not most, super-Earths have siblings, i.e., they are members of multiple-planet systems. In these systems, gravitational interaction between planets affects not only the formation but also their evolution and destiny. There have been several attempts to describe the dynamical evolution of these systems with N-body simulations. Due to the time consuming nature of such approaches, only limited range of initial parameters have been explored. It is not clear whether these models adequately represent the distribution of multiple planet systems.

The population synthesis method is based on a statistical mechanics approach to evaluate the most likely outcome of dynamical interactions between planets. In order to describe the statistical outcome of gravitational interactions and collisions between rocky planetary embryos, we have upgraded our population synthesis method to generate a series of many simulations. In paper VI (Ida & Lin, 2010), we described our analytic prescription for the collisions and close scatterings between modest-mass planets. We applied this method to show that super-Earths may be the in situ merger products of a population of embryos with Earth-mass or less which converged towards the proximity of their host stars through type I migration.

Radial velocity survey show that modest to long-period ( days) gas giants also have significant eccentricity. These planets most likely have formed in their natal disks with nearly circular orbits as in the case of Jupiter and Saturn (Pollack et al 1996, Ida & Lin 2004). Their eccentricity may have been excited by the gravitational perturbation of their planetary siblings (Rasio & Ford 1996, Zhou et al 2007, Jaric & Tremaine 2008, Chatterjie et al 2008). Indeed, many gas giant planets are members of known multiple-planet systems.

In order to simulate the outcome of dynamical interaction between gas giants and their siblings including other gas giants as well as less-massive protoplanetary embryos, we construct in this paper, construct a ”multiple-planet-in-a-disk model”. This approach is a natural continuation of the method presented in Paper VI. The main technical challenge is that the gravitational perturbations from gas giants significantly affect orbital configuration of the entire planetary system.

The method we have developed to approximate the evolution of systems of N gas giants is presented in §2. This prescription takes into account gravitational interactions between all the planets around their common host stars. It reproduces the general results of direct N-body simulations with a large reduction (by many orders of magnitude) in the computational cost. A detailed description of the method is given in the Appendix.

We incorporate this prescription into the latest version of the population synthesis models. In §3, we briefly recapitulate the methodology used in our population synthesis models including prescriptions for: 1) disk structure and evolution, 2) planetesimal growth and the planetary embryos’ dynamical isolation, 3) migration due to tidal interaction between protoplanets and their natal disks, and 4) resonant interaction between protoplanets.

In §4, we apply these prescription to simulate some sample multiple-planet systems. These results show that that 1) the most favorable location for gas giant formation in the core-accretion scenario is at a few AU’s, 2) the emergence of gas giants leads to the scattering of residual protoplanetary embryos, some of which may subsequently undergo gas accretion and form later generation gas giants, 3) as a consequence of dynamical instability, some gas giants in multiple-planet systems may be scattered to attain the observed eccentricity distribution.

An advantage of the population synthesis model is its computational efficiency. With it, we are able simulate in §5 rich sets of planetary systems based on the observed range of disk properties. We compare the results obtained from the population synthesis models with the observational data. Our models reproduce the observed planetary mass - eccentricity and semimajor axis - eccentricity correlations. These comparisons are important not only to calibrate some uncertainties in the model parameters but also to highlight the dominant physical processes.

Finally in §6, we summarize our results and discuss their implications.

2 Modeling of eccentricity excitation and ejection of giant planets as a result of gravitational instability

The objective of this section is to introduce a simple-to-use and computationally-efficient prescription to approximate the dynamical interaction between multiple planets around common host stars.

2.1 An overview of our technical approach

Before outlining the technical details of our prescription, we first qualitatively describe the principles of our analytic approximation for the gravitational interactions between several gas giants. This approach is guided by the results generated from various N-body simulations of the dynamical interaction between giant planets. Many such simulations have been carried out for the purpose of studying the statistical properties associated with close encounters between multiple planets especially with regard to the origin of the eccentric as well as short-period Jupiter-mass planets (e.g. Rasio & Ford, 1996; Weidenschilling & Marzari, 1996; Ford et al., 2000; Zhou et al., 2007; Ford & Rasio, 2008; Jurić & Tremaine, 2008; Marzari & Weidenschilling, 2002; Nagasawa et al., 2008; Chatterjee et al., 2008).

Most of the N-body experiments have been carried out with idealized initial conditions. They show that 1) the onset of dynamical instability is sensitively determined by the planets’ initial normalized (in terms of their Hill’s radius) separation and 2) dynamically unstable systems undergo orbit crossing and relaxation which generally lead to a Rayleigh distribution in the planets’ asymptotic eccentricity (Ida & Makino, 1992; Zhou et al., 2007). In order to check the validity of our prescription (to be presented below), we perform calculations with our analytic prescriptions for the same initial conditions as those adopted by the previous N-body simulations.

There are some analytic approximation for the stability of planetary systems. Two planets with masses and on initially circular orbits around a common host star with a mass () would immediately cross each other’s path if the difference between their semimajor axes () is smaller than a critical value,


where is the average of their semimajor axis and


is there mutual Hill’s radius. But, orbital crossings in such a system would not occur if initially is slightly larger than (Gladman, 1993). In §2.2, we present a prescription to evaluate the asymptotic semimajor axis and eccentricity distributions of two-planet systems which undergo orbital crossings.

In systems with three or more gas giants, the orbit-crossing condition is significantly different. For a finite duration of time, systems with initial larger than may be maintained in a quasi-stable state with a limited variation in the amplitude of eccentricity. However, in due course, such systems may undergo a transition with a rapid eccentricity increase which is followed by orbital crossings (Chambers et al., 1996; Lin & Ida, 1997; Marzari & Weidenschilling, 2002; Zhou et al., 2007). The time scale for the onset of this transition depends sensitively on the normalized initial orbital separation , the planet-to-star mass ratio, and weakly on the number of planets.

In the population synthesis model (see §3), we consider the possibility that the giant planets and rocky/icy planetary embryos may form, coexist, and evolve contemporaneously in a gaseous environment. During advanced stages of their formation, gas giants grow through runaway gas accretion, i.e., the time scale for them to double their mass decreases with . Their Hill’s radius and the width of their feeding zone also increases with their . The eccentricity of the nearby planetary siblings and residual embryos (with ) would be excited over their synodic periods with respect to the orbits of the gas giants.

While multiple planets’ eccentricity is excited by their secular interaction with each other, it is also damped by their tidal interaction with their natal disks (Artymowicz, 1993; Ward, 1993; Tanaka & Ward, 2004). It is appropriate to take into account these competing effects (Iwasaki et al., 2002) (see §3). Planet-disk interaction also leads to type I and II migration (for modest-mass embryos and high-mass gas giants respectively). Idealized prescriptions for isolated planets’ migration have already been incorporated into early versions of population synthesis models (Papers I & V). In a subsequent paper, we will develop and apply a new prescriptions for embryos’ type I migration (Paardekooper et al., 2011) and consider the feedback on the disk structure by multiple gas giants and their interference on each other’s direction and speed of migration.

In this paper, we focus on dynamical interaction between multiple gas giants in the limit that the difference between their semimajor axis may be reduced during the course of their migration. Two gas giants would capture each other onto their mutual low-order mean-motion resonances provided the convergent speed is slow. This effect has already been analyzed with our population synthesis models for interaction between rocky and icy embryos (Paper VI). In §2.2, we modify our prescription to take into the nonlinear interaction and determine the resonant capture condition between multiple gas giants. We also consider the critical converging speed above which the gas giants may overcome the resonances’ barrier and intrude into each other’s Hill’s sphere and undergo close encounters (see §3.6).

The intensity of planet-disk interaction reduces with the depletion of the disk gas. Observational signatures of protostellar disks diminish on a time scale of 3-5 Myr. The accretion rate from the disk to the central stars (inferred from the UV veiling) also declines on a similar time scale (Hartmann, 1998). Our population synthesis models take into account disk evolution and adjust the effects of gas accretion, eccentricity damping, and orbital migration accordingly. After the disk gas is severely depleted, both gas damping and accretion cease while secular interaction between the planets persists.

In Paper VI, we presented prescriptions for eccentricity excitations and merging events (giant impacts) among the rocky/icy planetary embryos. For these relatively low-mass objects, we construct an analytic approximation to determine the outcome of the dynamical interaction between planets based on celestial mechanics. This prescription quantitatively reproduce the results of N-body simulations. We applied this prescription to simulate the evolution of systems of 10–20 embryos with a total mass of a few . This initial condition corresponds to the the advanced stage of oligarchic growth when the embryos become dynamically isolated. We showed that close encounters generally lead to a velocity dispersion which is a significant fraction of the embryos’ surface escape speed. Within a few AU’s, embryos are unlikely to be ejected because is smaller than the local Keplerian speed . After the gas depletion, repeated close encounters lead to mergers and the emergence of one or two Earth-mass planets with few residual embryos. Thereafter, the magnitude of increases beyond the main sequence life span of solar type stars (ie yrs). We also used this prescription to show that some super-Earths may be the merger products of a population of embryos which migrated to the proximity of their host stars.

In §2.3 of this paper, we generalize our prescription to simulate the onset of dynamical instability and resultant orbital crossings among multiple gas giants in the outer regions of an evolving disk environment. The recoil motion excites the eccentricity of the retained planets to establish a Rayleigh’s distribution similar to that observed among eccentric gas giants (Rasio & Ford, 1996; Weidenschilling & Marzari, 1996; Lin & Ida, 1997; Marzari & Weidenschilling, 2002). In contrast to the embryos close to their host stars, the surface escape speed of typical gas giants exceeds the Keplerian speed in the outer regions of their natal disks where they emerge. Close encounters with gas giants at few AU from the host stars may lead to ejection rather than merger events.

The crossing timescale increases with the semimajor axis separations and decreases with the orbital eccentricities. As the population of gas giants declines through merger and ejection events or migration and stellar consumption, their normalized semimajor axis separation () is enlarged. Energy dissipation during cohesive collisions also leads to eccentricity damping. These combined effects lengthen by several orders of magnitude after each merger or ejection event.

Recent numerical simulations show that there is a chance for orbital crossings to eventually excite one of the gas giants to attain nearly parabolic orbit by Kozai effect (Nagasawa et al., 2008).1 The orbit of an inwardly scattered gas giant may be circularized if its pericenter distance is sufficiently close to the host star (Ivanov & Papaloizou, 2004). This scenario has been suggested as one of the mechanism for the origin of “hot Jupiters” (Rasio & Ford, 1996; Nagasawa et al., 2008). We have not incorporated tidal circularization and secular perturbations including the Kozai effect into the current version of the population synthesis models, although the overall effects of secular eccentricity excitation is implicitly taken into account through the assumption that the orbital crossings start on timescale .

2.2 Two gas giants case

The above description qualitatively outlines the overall methodology of our population synthesis model. We now provide some technical details on its latest upgrade, i.e., the construction of an analytic plus Monte Carlo prescription.

We begin with the simplest case of dynamical interaction between two gas giants. As we have indicated above, in this case, orbit crossing are expected to occur between two planets with nearly circular co-planar orbits which are initially separated in semimajor axis by . Extensive numerical simulations show that close encounters between these planets generally leads to the expansion of their semimajor axis separation and eccentricity excitation. We construct a method to compute the outcome of these close encounters. This prescription is calibrated with the statistical results generated from numerical simulations (Ford et al., 2000; Ford & Rasio, 2008). Detailed prescriptions are described in the Appendix A1.

The procedure of our model for two-planet encounters is structured in the following manner.


We first identify two interacting gas giants, with mass and .


We compute the maximum relative eccentricity attainable during a close scattering in terms of the two-body surface escape velocity () scaled by the local Keplerian velocity (), such that,


A probable eccentricity is generated for individual bodies with


where refers to the perturbing companios and


is a random number with a Rayleigh distribution of unit root mean square dispersion ,


If, at least one of the trial eccentricities (given by Eq. 4) exceeds unity, the body with the larger eccentricity would be identified as a planet to be ejected. The remaining planet is assumed to be retained. For it, a new eccentricity is generated (with Eqs. 6 and 7) together with an associated semimajor axis (with Eq. 45) (see §2.2.1).


If the trial eccentricities of both planets are less than unity, they would be adopted as the expected values for the two planets after their close encounters with some damping (see step 3-b in Appendix A1). The associated semimajor axis changes are evaluated in accordance with the prescription in §2.2.2.

This procedure is applied throughout the disk evolution. Prior to severe gas depletion, the damping of eccentricity due to the planets’ tidal interaction with their natal disks is taken into account. The combined effect of planet-planet scattering and planet-disk tidal interaction induce the two planets to enlarge their semimajor axis separation until they become dynamically isolated with . The efficiency of gas damping declines during gas depletion. In the absence of any significant amount of residual disk gas, two planets with would continue to scatter each other until either they merge with each other or one of them is ejected.

Planetary escape after close encounters between two planets

We consider the possibility that one member of a planet pair may attain sufficient recoil velocity to escape from the gravitational potential of their host star. We compare the results generated with our analytic plus Monte Carlo prescription with those obtained from N-body simulations which were computed without taking into account of any eccentricity damping effect (due to either planet-disk interaction or dissipative collisions.).

Since both interacting planets were initially bound to their host stars, conservation of energy implies that at least one of them must be retained. We label the retained and ejected planets as body 1 and 2 respectively. In Appendix A1, we derived the differential distribution function of the post-scattering eccentricity of the retained body 1 as




and and are masses of retained and ejected planets respectively. Since it takes less energy to eject the lighter body, we assume which is also consistent with the results of N-body simulations.

In Eqs. (6) and (7), the eccentricity distribution function depends on the ratio rather than the individual values of and . The results of numerical simulation of close encounters between two equal mass planets Ford & Rasio (2008) also show a negligible dependence of (in the range between ). Note from Eq. (3), ejection cannot be produced from close encounters between planets with ).

Figure 1: Dynamical instability and asymptotic properties of systems with two closely separated gas giant planets. All the models start with two planets with a total mass around a common host star. Their initial orbits are nearly circular and coplanar with a unstable semimajor axis separation (ie ). In the absence of gas drag, one planet in each system is eventually ejected while the other is retained with an changed orbit. The eccentricity distributions of the retained planets are plotted for three sets of models (with and 0.5 from left to the right). The solid and dotted lines represent results obtained from our N-body simulations and our analytic plus Monte Carlo prescription.

In order to verify the validity of the analytic plus Monte Carlo prescription, we compare the results generated with it and those obtained from a series of N-body simulations with a set of similar initial conditions as those presented by Ford & Rasio (2008). In all of these models, we consider the interaction between two planets with small initial eccentricities and inclinations. We set their initial semimajor axis to be and is randomly specified to be in the range of so that they generally undergo orbital crossings.

We consider three series with and 0.5. In all cases, . For each value of , we carried out 300 N-body simulations with different random number seeds for initial relative orbital phases between the two planets. In the absence of gas drag, one planet is eventually ejected. The differential eccentricity distribution function of the retained planet is represented by the solid lines in Figure 1. For the cases with and , the more massive planet 2 is generally retained. But, this mass preference is reduced for the cases.

We also adopted the analytic plus Monte Carlo prescription to generated models for the same sets of values with the distribution in Eqs. (6) and (7). The results of the Monte Carlo calculations are represented by the dotted lines in Figure 1. They are in good agreement with those of the N-body simulations.

Depending on the difference between their semimajor axes, two planet systems either rapidly undergo orbital crossings (if ) or completely avoid close encounters. Some systems may evolve into unstable states through disk migration (which leads to a reduction in ) or gas accretion (which widens ). In their natal disks, the eccentricities of gas giants are also damped. In our population synthesis simulations, we include the effect of eccentricity damping for the retained planet. The efficiency of this effect depends on the mass ratio between the planet and local disk (see Appendix A1).

Finally, the semimajor axis of the retained planet () is obtained from the conservation of energy, such that


where is the semimajor axis before scattering. For computational simplicity, we neglect 1) the decrease in the semimajor axis change associated with the eccentricity damping of the retained planet, and 2) the finite kinetic energy of the ejected planet.

Retention of both gas giants after close encounters between them

If the trial eccentricities (generated at step 3) of both planets are less than unity, they would both be adopted as the asymptotic values. The effect of the eccentricity damping due to the post-encounter disk-planet interactions is applied to both planets with the same prescription as that described in the previous section (i.e., in §2.2.1). In §5.4, we also discuss the eccentricity damping due to disk gas accretion onto a planet.

Close encounters launch both planets into eccentric orbits from nearly circular initial orbits through changes in their energy and semimajor axis. We assume that a less massive body (represented by the subscript “2”) is always scattered outward such that its semimajor axis prior to the close encounter become its new periastron distance, such that


The new semimajor axis of the inwardly scattered planet () is given by the energy conservation such that


2.3 Three gas giants case

In contrast to the stability criterion in Eq. (1) for two planet systems, three planet systems may maintain relatively low-eccentricity states until gas in their natal disk is severely depleted and then undergo close encounters after a time scale . The magnitude of is a sensitive function of their initial separation in semimajor axis .

In order to develop an easy-to-use and robust prescription for close encounters among three planets, we consider the situation that is smaller than the expected main sequence lifespan of their host stars. In these systems, eccentricities of the planets are always excited during the close encounters. Since these events take place in the absence of residual gas drag, the planets’ excited eccentricity can only be damped through 1) dissipative merger events, 2) tidal dissipation in the proximity of their host stars, or 3) dynamical friction by a residual population of planetesimals. In this paper, we neglect the effect of the tidal dissipation and residual planetesimals.

Following a similar procedure from the previous section, we construct an analytic plus Monte Carlo prescription for the outcomes of close encounters and compare it with the result generated from N-body simulations of three-giant-planet systems which undergo orbital crossings in a gas-free environment (Marzari & Weidenschilling, 2002; Nagasawa et al., 2008; Chatterjee et al., 2008). In this case, two or more planets ensure repeated encounters with a range of impact parameter. Physical collisions occur when the impact parameter is less than the sum of the planets’ radii. Such merger events reduce the number of residual planets bound to the host star and consequently stabilize the remaining systems.

Slightly wider scatterings, with impact parameter marginally larger than the sum of planets’ radii, lead to recoil velocities comparable to their surface escape speed. At a distance of several AU from their solar-type host stars, such strong orbital deflection generally leads to the escape of one planet from the gravitational potential of the host star. The remaining two planets are retained with modest to high eccentricities and widely separated semimajor axes such that the residual systems are stabilized.

We present below our analytic plus Monte Carlo prescriptions for the outcomes of close encounters in three-planet systems. Based on N-body simulations, we assume that this intense interaction leads to 1) mostly the ejection of one planet and 2) a significant fractional () incidents of physical collisions and merger events. In both outcomes, the residual systems are left with two remaining planets. The main objective of our prescription is to evaluate the asymptotic eccentricities and semimajor axes of these two retained planets.

We justify the validity of our prescription based on the direct comparison of our results with those generated from the N-body simulations. The technical procedure essentially follows that constructed for the two planet case in §2.2, with some modifications applied to the inwardly and outwardly scattered retained planets. Their detailed description and physical justification are given in Appendix A2. We outline the following basic steps.


We first set initial conditions of the three giants and assign the identity of the most massive planet to be ,


We assign a probability () that two planets physically collide and merge with each other (assumed to always occur after each collision). A merged planet is formed from a randomly selected pair of planets which are involved in the orbital crossings. Its resultant semimajor axis and eccentricity are calculated under the assumed conservation of total mass, orbital energy and angular momentum.


We also assign a probability () that one of the three planets is ejected while the other two residual planets are retained.

For the ejection events,


we compute the maximum eccentricity for all three () planets, in accordance with the procedure described in step 2 in Appendix A2. The magnitude of the most massive body tends to be less than those of the other bodies.


We select the planet with the largest value of of to be the ejector.


We determine asymptotic eccentricities of the two retained planets to be


where is the physical radius of planet and is a random number following a Rayleigh distribution of unit root mean square. The above expression is similar to that in Eq. (4) with the consideration of incomplete stirring (see Appendix A2 for justification).


Among the two retained planets, we select the inwardly scattered body with a probability distribution function which is weighted by the square of their masses.


Based on the nature of close scattering, the new semimajor axis of the outwardly scattered body is prescribed to be


where is the excited eccentricity of the outer planet, and and are the maximum and minimum semimajor axes of the three planets in their initial state prior to orbital crossings.


The semimajor axis of the inwardly scattered body is determined by


where is the mass of the inner planet and is the total energy calculated from the initial semimajor axes of the three planets.

In order to compare the results of our prescription with those obtained from previous N-body simulations, we adopt a set of models with three planets in nearly circular and coplanar orbits as initial conditions in step B1. For these test cases, we consider the situation that orbital crossings have already been initiated after a timescale of has elapsed, so that the initial orbital separations are arbitrarily set. In the actual population synthesis simulations, is calculated and compared with the system time at each timestep (see Appendix A3).

Figure 2: The eccentricity and semimajor axis distributions of surviving planets in the case of three Jupiter-mass planets with AU, AU, AU and initially. The dotted and solid lines represent the distributions of inner and outer planets in the top and middle panels. The left panels show the results of models generated with our analytic plus Monte Carlo prescription and the right ones are the results of 75 models obtained from Nagasawa et al. (2008)’s N-body simulations.

The comparison of our prescribed method and N-body simulations is shown in Figure 2 where we illustrate the eccentricity and semimajor axis distributions of the surviving planets for a test case with three Jupiter-mass planets (). We adopt a set of initial conditions with AU, AU, AU and , in accordance with some existing N-body simulations by Marzari & Weidenschilling (2002) and Nagasawa et al. (2008) (their set N without tidal damping).

The left panels show the results of models generated with our analytic plus Monte Carlo prescription and the right panels illustrate the results of 75 models of Nagasawa et al. (2008)’s N-body simulations. The initial orbital phase distributions of planets are modified in each of the N-body simulations. In the analytic plus Monte Carlo simulations, the seed for random number generation to determine eccentricity is changed in each simulation.

The N-body simulations are carried out in 3D, i.e., the orbital inclinations of each planet are also calculated. In our analytic plus Monte Carlo simulations, although the planets’ inclinations are not explicitly computed, their distributions can be evaluated by (in radian) (Nagasawa et al., 2008).

In this equal-mass example, the top panel of Figure 2 shows that the analytic plus Monte Carlo prescription produces a broad differential eccentricity distribution which is similar for the inner and outer planets. We prescribed (in step B5) a Rayleigh’s distribution through which has the form (see Eq. 11). The peak of this distribution function occurs at .

The middle panel of Figure 2 shows a bimodal semimajor axis distribution because inwardly scattered planets are generally well separated from outwardly scattered planets. The middle panel of Figure 2 shows a correlation between semimajor axis () and eccentricity () of the outwardly scattered planets. This correlation arises because we assume (in step B7) that these planets are scattered from the region near the initial locations and their periastron radii are close to their original semimajor axes despite some outward diffusion before the scattering (see Appendix A2). These results are all consistent with the results obtained by the previous N-body simulations.

Figure 3: The eccentricity and semimajor axis distributions of the retained planets in the case of . All other initial conditions are the same as those in Fig. 2.

Marzari & Weidenschilling (2002) also carried out simulations with and . They found that in 54 runs of 64 runs (corresponds to probability), the more massive planet (with ) is scattered inward and the inner retained planets have a broad eccentricity distribution peaked around . Here, we carry out a similar simulation with the analytic plus Monte Carlo prescription to generate models. These results show that the massive planet is scattered inward in 8382 models with a peak eccentricity (in accordance with eq. (11)).

For comparison purpose, we also simulated 100 N-body models. For the massive planets’ inward scattering probability and their average peak eccentricities, the results obtained with the analytic-plus Monte Carlo prescription are almost consistent with those generated with the N-body simulation, which provides support for the validity of the mass-weighted probability function we have assumed. The middle panel in Fig. 3 show two sharp peaks in the semimajor axis distribution for the inwardly scattered planets. These peaks are also found in the N-body simulations. They are generated from the requirement for energy conservation such that


If the inner planet has a mass AU whereas AU if .

2.4 General case

In the population synthesis simulations, multiple giant planets and rocky/icy planetary embryos co-exist. Here we describe a summary of the prescriptions for dynamical interactions in the general case.

During the early stages of their dynamical evolution, newly formed planets are embedded in and interact with their natal disks. In the presence of the disk gas, we only consider scattering between two planets when their orbital separation becomes less than the critical value with the prescription given in §2.2.

We calculate the planets’ eccentricity damping timescale (), due to their gravitational interaction with the residual disk, and the orbit-crossing timescale () for all pairs of existing planets in gas-free environment. As the disk gas is depleted, eventually becomes larger than and the planets’ eccentricity grows due to secular excitation (Iwasaki et al., 2002). It then becomes possible for planets with initial separation greater than to cross each other’s orbits.

Orbital crossings and close encounters among gas giants strongly affect the asymptotic global structure of planetary systems. The gas giants’ eccentricities are highly excited during their orbital crossings. The order of their semimajor axis is occasionally swapped. Their secular resonances may sweep across the entire systems multiple times. Most of the residual small planets are either ejected or collide with their host stars (Matsumura et al., 2013).

For the evaluation on the outcome of planetary scattering after the severe disk-gas depletion, we apply the following prescription. The details are presented in Appendix A3. Note that the order of the prescription presented in this section is different from those in Appendix A3, for the purpose of easier understanding.


We identify the ”giant planets” from a list of all the planets in the system with the criteria that they have (i) a mass and (ii) (Eq. 51).


We apply the procedures for orbital crossings among the small planets (as outlined in Paper VI) if there is only one (or less) giant planet.


In systems with two (and only two) widely separated (with ) giant planets, they retain their initial kinematic properties (because they do not significantly perturb each other). Their dynamical influence on the less-massive planets is computed independently, following the procedures in step C2.


In systems with two (and only two) closely separated (with ) giant planets, the outcome of their dynamical interaction is computed following steps A2-A3 in §2.2. We also remove all the small planets under the assumption that they are ejected by the gas giants’ sweeping secular resonances.


In systems with 3 or more giant planets around the same host stars, we first evaluate for all pairs. For this case,


if for all giant planet pairs is longer than the system time (), we would neglect any interaction between them and compute their influence on the low-mass planets following step C2-a.


If for one or more pairs of giant planets is less than but there are no triples (overlapping pairs) with , we would determine the outcome of interaction between the giant planets following step C2-b or C2-c.


Any groups of three giant planets around the same host stars would be labelled as interacting triples if for at least two planet pairs. In this case, we follow steps B2-B8 for the triple giants interaction case (see §2.3). All the small planets would be removed under the assumption that they are ejected by the triple giants’ sweeping secular resonances.


After the kinematic properties of interacting triple systems are modified to a set of new configuration in the previous step, we re-calculate the for each pairs of giant planets. Steps C2-C3 are repeated until either 1) there are only two giant planets remaining around the host stars or 2) when yrs.


For systems with more than 3 giant planets, we follow the same procedure as C3-a and C4 and eject all the small planets. But, we let only one giant planet be scattered inward and all other giant planets are assumed to be scattered outward.

3 Prescriptions of other physical processes in the population synthesis models

We follow Paper VI for prescriptions of a) planetesimals’ growth through cohesive collisions, b) the evolution of planetesimal surface density, c) embryos’ type I migration, and their stoppage at the disk inner edge. For the gas giants, d) the onset, rate, and termination (through gap opening or and global depletion) of efficient gas accretion, and e) their type II migration, we follow prescriptions in Paper IV, except for a slight modification on gas accretion termination due to global depletion. Here we present a brief summary.

3.1 Disk models

We adopt the minimum mass solar nebula (MMSN) model (Hayashi, 1981) as a fiducial set of initial conditions for planetesimal surface density () and introduce a multiplicative factor (). For gas surface density , although we adopt the -dependence of steady accretion disk with constant viscosity (), rather than that of the original MMSN model (), we scale by that of the MMSN at 10AU with a scaling factor (). Following Paper IV, we set


where normalization factors and correspond to 1.4 times of and at 10AU of the MMSN model, and the step function inside the ice line at (eq. [17]) and 4.2 for [the latter can be slightly smaller () (Pollack et al., 1994)]. We specify an inner disk boundary where vanishes and planetesimals’ type I migration is arrested (see §3.5).

Neglecting the detailed energy balance in the disk (Chiang & Goldreich, 1997; Garaud & Lin, 2007), we adopt the equilibrium temperature distribution of optically thin disks prescribed by Hayashi (1981) such that


where and are stellar and solar luminosity. We set the ice line to be that determined by an equilibrium temperature (eq. [16]) in optically thin disk regions,


The effects of evolution of the disk temperature and associated migration of the ice line adopting more realistic thermal disk structure (Garaud & Lin, 2007; Oka et al., 2011; Min et al., 2011) will be studied in a subsequent paper.

Dependence of disk metallicity is attributed to distribution of , where and are initial values of and , respectively. Due to viscous diffusion and photoevaporation, decreases with time. For simplicity, we adopt


where is disk lifetime (for detailed discussion, see Paper IV). The constant self-similar solution obtained by Lynden-Bell & Pringle (1974) is expressed by with an asymptotic exponential cut-off at radius of the maximum viscous couple. In the region at , decreases uniformly independent of as the exponential decay does, although the time dependence is different. In the self-similar solution, at decays as . In the exponential decay model that we adopt, decays more rapidly as becomes larger at . If the effect of photoevaporation is taken into account, decays rapidly after it is significantly depleted, so that the exponential decay partially mimics the effect of photoevaporation. However, the final planet distributions hardly depend on whether or . In this paper, we adopt . The parameter values of , , and [Fe/H] are specified for each run.

3.2 From Oligarchic Growth to Isolation

On the basis of oligarchic growth model (Kokubo & Ida, 1998, 2002), growth rate of embryos/cores at any location and time in the presence of disk gas, is described by where, after correcting some typos in Paper IV,


where is the mass of the embryo (core), and we set the mass of the typical field planetesimals to be g.

If type I migration is not effective, embryos would formed through oligarchic growth and attain a local isolation mass, where the full width of feeding zone of an embryo with a mass is given by (Kokubo & Ida, 1998, 2002)


and is the Hill radius for two bodies with comparable masses. From equations (15) and (20) we find


However, type I migration may induce embryos to migrate before they acquire all the residual planetesimals within their feeding zone and acquire their isolation mass. Type I migration may also be stalled near some trapping radius. At these locations, the congregation of planetesimals may increase the magnitude of and enlarge the magnitude of .

We compute the evolution of distribution due to accretion by all the emerging embryos in a self-consistent manner. The growth and migration of many planets are integrated simultaneously with the evolution of the -distribution. We set up linear grids for across the disk with typical width of AU. We introduce a population of seed embryos, all with an initial mass g (i.e., that of the residual planetesimals) and compute their mass accretion rate. In the inner disk region, we set their initial separation of the full feeding-zone width () of embryos with local asymptotic isolation mass .

The planetesimals’ growth time scale increases rapidly with their distance from the central stars and embryos in the outer disk region of the disk are unlikely to attain a local isolation masses within the life span of their host stars. There, we place seed embryos which are separated by feeding zone width of embryos with masses evaluated (using Eq 19) for the local after Gyr. We follow the growth of the seed embryos due to planetesimal accretion in accordance with Eqs. (19) and (20). The planetesimals’ mass accreted by the embryos is uniformly subtracted within the embryos’ current feeding zone. We follow the evolution of throughout the disk and use its values to evaluate embryos’ local accretion rates. We also use to estimate the strength of planetesimals’ dynamical friction on the embryos.

During the early phase of evolution, embryos are embedded in their natal disks. Despite their mutual gravitational perturbation, the embryos preserve their circular orbits due to the gravitational drag from disk gas(Artymowicz, 1993; Ward, 1993) and dynamical friction from the residual planetesimals(e.g., Stewart & Ida, 2000). However, after the disk gas is severely depleted, the efficiency of eccentricity damping mechanism is reduced. The embryos’ eccentricity grow until they cross each others’ orbits. We use the prescription in Paper VI and Appendix A3 to compute the occurrence and consequence of giant impacts between embryos.

3.3 Type I migration

Type I migration of an embryo is caused by the sum of tidal torque from the disk regions both interior and exterior to the embryos. The rate and direction of embryos’ migration are determined by the differential Lindblad and corotation torques. While a conventional formula of type I migration assuming locally isothermal disks (e.g., Tanaka et al., 2002) shows that the migration is always inward, recent developments of type I migration of isolated embryos in non-isothermal disks (e.g., Paardekooper et al., 2011) show the magnitude and sign of the tidal torque, especially that due to corotation resonances is a sensitive function of the surface density and temperature distribution of the disk gas.

Since the properties of migration depend on detailed thermal/dynamical structure of disks and ”saturation” degree of the corotation torque (e.g., Kretke & Lin, 2012), substantial discussions are required for the effects of type I migration in non-isothermal disks. The pace of individual embryos’ type I migration in dense multiple planet systems remains uncertain because mutual perturbation between nearest neighbors may modify their Lindblad and especially corotation torque. On the other hand, as shown below, the eccentricity distributions that we are primarily interested in here are not sensitive to a change in the formula of type I migration. Thereby, in order to highlight the validity of our treatment of dynamical interaction between multiple planets in this paper, we postpone the detailed discussions on the non-isothermal type I migration to a subsequent paper. While we also show a result with a non-isothermal formula, in most of the results we present here, we use a conventional formula of type I migration in isothermal disks derived by Tanaka et al. (2002) with a scaling factor :


The expression of Tanaka et al. (2002) corresponds to , and for slower migration, . We study a dependence of the eccentricity distribution by changing a value of , rather than using variations of non-isothermal formula.

We assume type I migration ceases inside the inner boundary of the disk, because is locally zero there. For computational convenience, we set the disk inner boundary to be the edge of the magnetospheric cavity at AU.

3.4 Formation of gas giant planets

Prescriptions for formation of gas giant planets are the same as those used in Paper IV, except for slight modification for reduction and termination of gas infall. Embryos are surrounded by gaseous envelopes when their surface escape velocity becomes larger than the sound speed of the surrounding disk gas. When their mass grow (through planetesimal bombardment) above a critical mass


both the radiative and convective transport of heat become sufficiently efficient to allow their envelope to contract dynamically(Ikoma et al., 2000).

In the above equation, we neglected the dependence on the opacity in the envelope (see Paper I and Hori & Ikoma (2010)). In regions where the cores have already acquired isolation mass, their planetesimal-accretion rate would be much diminished (Ikoma et al., 2000; Zhou et al., 2007) and can be comparable to an Earth mass . But, gas accretion also releases energy and its rate is still regulated by the efficiency of radiative transfer in the envelope such that


where is the planet mass including gas envelope.

In Paper I, we approximated the Kelvin-Helmholtz contraction timescale of the envelope with


where is the contraction timescale for . Since there are uncertainties associated with dust sedimentation and opacity in the envelope(Pollack et al., 1996; Helled et al., 2008; Hori & Ikoma, 2011), we adopt a range of values years and –4. Here we fix and change (In a nominal case, we use years).

Equation (24) shows that rapidly increases as grows. But, it is limited by the global gas accretion rate throughout the disk and by the process of gap formation near the protoplanets’ orbits. The disk accretion rate is


During the advanced stage of disk evolution, we assume both and evolves where is the gas depletion time scale. The rate of accretion onto the cores cannot exceed .

An (at least partial) gap is formed when the planets’ tidal torque exceeds the disk’s intrinsic viscous stress (Lin & Papaloizou, 1985). This viscous condition for gap formation is satisfied for planets with


Then, type I migration is switched to type II migration. Fluid dynamical simulations (e.g., D’Angelo et al., 2003; D’Angelo & Lubow, 2008) show that some faction of gas still flows into the gap. According to this, we allow a protoplanet to continue accreting the residual gas which flows past it, as shown below. These fluid dynamical simulations, however, also show that the accretion rate rapidly decreases with after exceeds a Jupiter mass. According to this result and the analysis by Dobbs-Dixon et al. (2007), we completely terminate gas accretion, when the planet’s Hill radius () becomes larger than 2 times of disk scale height (), which corresponds to (extended) thermal condition (Lin & Papaloizou, 1985), that is,


In general, our prescriptions for gas accretion rates onto the cores are


where is that in the absence of any feedback on the disk structure, i.e., that without the effect of gap opening,


and is a reduction factor due to gap opening,


The formula for is constructed to avoid any abrupt truncation. Some gas may leak through the gap. But, for planets with (D’Angelo et al., 2003), the diffuse gas flows along horseshoe streamlines without being accreted by the planet (Dobbs-Dixon et al., 2007). Further discussion on this issue will be discussed in a future paper.

In Papers I-III, we incorporated the effect of the global depletion by setting a limiting value for the planet’s mass to be . This prescription was modified in Paper IV, in which we constrained . These previous prescriptions do not take into account of the viscous diffusion of gas from other disk regions. In order to consider the possibility that gas accretion may also occur in the inner disk regions where the local gas content is limited, we evaluate in this paper both and using the instantaneous values of . The quenching of gas accretion due to the disk’s global depletion is taken into account without any additional specification. This modification enables us to compute the gas accretion rate onto cores with relatively small semimajor axis.

3.5 Type II migration

During the gap formation, embedded gas giants adjust their positions in the gap to establish a quasi equilibrium between the torque applied on them from the regions of the disk both interior and exterior to their orbits. Subsequently, as the disk gas undergoes viscous diffusion, this interaction leads to type II migration.

We assume that planets undergo type II migration after they have accreted a sufficient mass to satisfy the viscous () condition for gap formation, as commented in section 3.4.

While increases, the disk mass declines due to stellar and planetary accretion and photoevaporation. While the disk mass exceeds (during the disk-dominated region), planets’ type II migration is locked with the viscous diffusion of the disk gas. During the advanced stages of the disk evolution when its mass become smaller than (during the planet-dominated region), the embedded planets carry a major share of the total angular momentum content. A significant fraction of the total (viscous plus advective) angular momentum flux transported by the disk gas is absorbed by the planet in its orbital evolution (for detailed discussion, see Paper IV).

For the disk-dominated regime, the migration timescale is given by


For the temperature distribution given by Eq. (16), and reduces to


where is global disk viscous-evolution time scale and is the disk radius where the viscous couple attains a maximum value. (The magnitude of is approximately the characteristic disk size). In general, is a fraction of the disk lifetime except for the planets which formed near .

Observationally inferred disk lifetime is 1-10Myrs. If this lifetime is comparable to the disk viscous-evolution time scale, AU. In this limit, Eq. (33) would imply that many gas giants form just beyond the snow line migrate to the proximity of their host stars and become hot Jupiters unless they emerged during the advanced (planet-dominated) stage of disk evolution or the type II tidal torque is significantly reduced by the disk flow through the gap or by the mutual interaction between multiple planets and their natal disks.

In Paper IV, we showed that the type II migration time scale in the planet-dominated regime can be evaluated in terms of


where is an efficiency factor associated with the degree of asymmetry in the torques between the inner and outer disk regions. If the inner disk is severely depleted, . But, in evolving disks with comparable surface density on either sides of multiple planets, the degree of torque asymmetry has large uncertainties in both the disk and planet dominated regimes. Although the reduction factors for the disk and planet-dominated regimes may be independent from each other, we adopt the same reduction factor for both regime and treat the factor as a model parameter. In the results shown in this paper, we set . The kinematic distribution of the planets are weakly affected by the choice of , as long as . The eccentricity distribution of the emerging gas giant planets is not strongly affected by the magnitude of .

3.6 Resonant capture

The pace of both type I and II depends on both the disk structure and . As different-mass planets migrate, the separation between their semimajor axis may either increase or decrease through divergent and convergent migration respectively. Multiple planets may also converge near some special disk locations such as the edge of the magnetospheric cavity near the inner disk boundary. As the migration of first-born planets is stalled at these trapping radii, later-generation planets may continue to arrive.

The planets gradually converging orbits capture each other on their mutual mean-motion resonances provided their migration time scale across the resonances’ width is shorter than the resonances’ libration time scales. After they enter into their mean motion resonances, these planets have a tendency to migrate together while maintaining the ratio of their semimajor axes (McNeil et al., 2005; Ogihara & Ida, 2009). In their simulations, Ogihara & Ida (2009) and Ogihara et al. (2010) found that outside the disk inner boundary at AU, it is possible for a few dozen rocky planetary embryos to be trapped in adjacent embryos’ mean-motion resonances. During and after the depletion of the disk gas, their orbits become unstable, they undergo orbit crossing and cohesive collisions to merge into a few super-Earths.

In Paper VI, we implemented a prescription for resonant capture between rocky/icy planets into our sequential planet formation model. With this prescription, we constructed population synthesis models for multiple super-Earth systems. These simulations showed a rich population of embryos accumulate near the inner disk boundary. Although differential migration leads to embryos’ orbital convergence, stochastic secular perturbation between multiple embryos induces orbital diffusion. Instead of locking into the most widely separated 2:1 mean motion resonances, the nearest neighbors are packed into low-order main-motion resonance with semimajor axis separations comparable to several . The Hill radius is defined to be where and are the masses of the closest pairs, is the host star’s mass, and is mean semimajor axis of the pair.

For computational simplicity, we adopted in Paper VI, a prescription that, as a consequence of converging migration, the closest embedded embryo pairs enter into and are locked in mean motion resonances with a separation of between their semimajor axes. After capture each other, the resonant embryos migrate in steps while maintaining their semimajor axis ratios and sharing angular momentum loss associated with the gravitational interactions between the disk and individual embryos.

In the present paper, we consider the possibility that the converging differential migration may proceed so rapidly that the neighboring embryos may not have time to capture each other into mean motion resonances. We evaluate the condition for a dynamical equilibrium in which convergence time scale (due to differential migration) matches with the mean motion libration time scale (see Paper VI). If this equilibrium is established with an orbital separation which is less than , we assume the converging embryos will collide and coalesce (see §2.2).

We further extend the resonant capture condition for planet pairs including both super-Earths and gas giants. Based on the rate of orbital decay due to gas drag, Shiraishi & Ida (2008) derived a criterion for non-resonant planetesimals to enter the feeding zone of a growing gas giant. Here we adopt their criterion for the embryos, the migration of which is determined by their tidal interaction with the disk gas. We replace the migration rate due to gas drag with differential speed between type II migration of a giant and type I migration of rocky/icy embryos.

In the case that a planet (an embryo or another giant) enters a feeding zone of a gas giant, we apply the prescription in §2.1 to determine the outcome of their encounters. Although the treatment in §2.1 is constructed for scattering between two gas giants, the same prescriptions are also applicable for the case where an embryo is scattered by a giant. In this limit, the perturbations from the giant is so strong and dominant that these events tend to result in the ejection of the embryos or in widely separated orbits rather than cohesive collisions.

In Paper VI, we assumed that the disk inner edge is a rigid wall for type I migration of embryos. This assumption is based on a strong ”eccentricity trapping” effect(Ogihara et al., 2010). The torque exerted onto the innermost embryo in an eccentric orbit is so strong that it can halt the type I migration and offset the angular momentum loss of several outer embryos. The magnitude of torque associated with embryos-planet interaction depends on the planet’s mass. The condition for a single planet (with a mass ) to frustrate the inner migration of other embryos (with masses ) at the inner boundary (Ogihara et al., 2010) is


Leakage of innermost embryo across the disk inner boundary would be imposed if the above condition is not satisfied. We also allow the leakage of embryos across the disk inner boundary if they lie in the path of an inwardly migrating gas giant planet.

4 Simulated individual systems

Before we enter into an comprehensive discussion on the statistical properties and population synthesis of multiple-planet systems, we first highlight the dynamical interaction among several bodies around common host stars with a few sample simulations of individual systems. Such studies have not been possible with the ”one-planet-per-disk-model” prescription in which planet-planet interactions are neglected (e.g., Papers I-V and the models of Mordasini et al. (2009a, b)).

Figure 4: Growth and migration of planets in a system with and . In the color version, the green, blue and red lines represent planets with their main component is rock, ice and gas, respectively.

For illustration purpose, we choose for Model 1, , ([Fe/H]), years, years around a solar-mass star (). Figure 4a and b show time and mass evolution of planets for Model 1. The green, blue and red lines represent rocky, icy, and gas giant planets with their main component being rock, ice and gas, respectively. Although rocky and icy planets are formed respectively interior and exterior to the snow line, their migration and dynamical interaction may lead to their spatial diffusion. The main composition of the planets may change through gas or planetesimal accretion or embryo-embryo collisions (giant impacts).

The results in Panel b show that at 0.1–1AU small embryos grow in situ until their masses reach and then undergo type I migration to accumulate near the inner boundary of their natal disks. We can determine the critical embryo mass for transition from local mass growth to efficient type I migration by matching Eqs. (19) and (22) such that


A few small () embryos are also induced to undergo orbital decay. The orbital evolution of these low-mass embryos is induced by the perturbation from the migrating planets.

Many resonant embryos accumulate in the vicinity of the disk inner boundary which is set to be 0.04 AU (Panel a, Figure  4). These planets are preserved in the disk region provided the condition (35) is satisfied. But when the total mass of the trapped planets in the proximity of the inner disk edge exceeds the critical mass given by eq. (35), innermost planets are driven to cross the disk inner edge until the retention condition in Eq . (35) is restored. In this model, more than half of the embryos arrived in the inner disk region eventually cross its inner boundary. We assume that after planets migrate inside the central magnetospheric cavity in the disk, their tidal and magnetic interaction with their host stars provide an effective torque to induce them to undergo further orbital decay. Note that in the weak-field limit, the size of the magnetospheric cavity is much smaller and the disk inner boundary is adjacent to the stellar surface. In this limit the innermost planets are unlikely to survive (Paper VI).

In the intermediate disk region, a migrating core attains just outside the ice line at AU. This core starts its runaway gas accretion without any significant type I migration. When it evolves into a gas giant with a surrounding gap, it undergoes type II migration which is much slower than type I migration. The emerging gas giant scatter and eject nearby embryos (they are represented by the lines which are discontinued at yrs.)

Figure 5: Planetary systems as functions of and . The radii of circles are proportional to one cube of planetary masses; legends are shown in upper right of the figure. In the color version, the green, blue and red lines represent planets with their main component is rock, ice and gas, respectively.

In order to explore the dependence of various model parameters, we simulated a series of models with a range of different and values. We use the same values (as Model 1) for [Fe/H], yrs, and . The asymptotic distributions of mass, semimajor axis and eccentricity of planets for these models are summarized in Figure 5.

For high values of , type I migration is so efficient that gas giants are rarely formed. In runs with , rocky/icy planets are not retained at AU except the vicinity of the disk edge at which type I migration is halted.

Figure 6: Planetary systems with as functions of . In order to test the dispersion of the simulated results, we choose three random number seeds for the scattering prescription. Models in the middle column are identical to those in the middle column in Figure 5.

For models with and 0.01, the planet distributions are sensitively determined by the total planetesimal mass in the disk, i.e., the value of . In our procedures for close scattering (see §2.2) depend on the choice of random number seeds. In order to test the robustness and dispersion of our models, two additional series of models are carried out for with different random number seeds for the scattering prescription (Figure 6).

All these models show similar planetary systems’ dependence on values of . With relatively high values of , the emerging planetary systems consist mostly of gas giants. In these massive disks, the amount of rocky/icy materials is sufficiently large to enable the formation of multiple cores of gas giants. Dynamical instability leads to giant planets to undergo close encounters which excite their eccentricity. Some giant planets may collide and coalesce into massive and eccentric merger products (Lin & Ida, 1997). Other planets attain eccentric orbits including those which are ejected from the system or scattered to endure close encounters with their host stars. During these violent dynamical relaxation, the orbits of rocky/icy planets are strongly perturbed by the giant planets and they are generally not retained.

In contrast, disks with moderate masses (or values of ) generally nurture fewer and more widely separated cores. Gas accretion leads to the emergence of either single or widely separated and sparsely populated giant planets. The mass of these planets is determined by the gap formation condition and is generally smaller than that of merger products. The mutual interaction between single gas giant planets and residual embryos or between sparsely populated gas giant planets is generally too weak for them to become dynamical unstable. They then retain relatively low eccentricities.

Comparison of Figures 5 and 6 indicates relatively massive and eccentric gas giant planets are formed out of massive disks whereas disks with modest masses generally produce gas giant planets with lower masses and eccentricities. Through the population synthesis models (next section), we show that this correlation can be used to infer initial formation condition of mature planets long after their natal disks have been depleted.

5 Population Synthesis of Planetary Systems

This set of prescription enable us to incorporate gravitational interactions among planets into our population synthesis models. The main objective of simulating ”multi-planets-in-a-disk” models is to generate mass and eccentricity distributions of the emerging multiple-planet systems. These quantities can be directly compared with observational data.

5.1 Initial conditions

We adopt a range of disk model parameters (see §3.1) and assign each their values with appropriate statistical weights to match the observed distribution of disk properties. We assume that and have log normal distributions in ranges of 0.1–10 and years, respectively. We use a normal distribution of [Fe/H] in a range of -0.2 to 0.2 which corresponds to that of typical target stars in various on-going radial velocity surveys. Since the target stars of radial velocity surveys are mostly G dwarfs, we use a log normal distribution, in a range of 0.8-, for the host stars’ mass ().

The initial distributions of seed planetesimals are described in §3.1 and their growth rate and migration time scales are specified in §3.2–3.5. Prescriptions for dynamical interactions among multiple planets are presented in §2.1–2.3 and 3.6.

5.2 Distributions of mass, semimajor axis and eccentricities

The asymptotic , , and distributions of a population of multiple-planet systems are shown in Figures 7. Similar to previous figures, the green, blue and red symbols represent rocky, icy and gas giant planets, respectively.

Models with similar parameters have been simulated in Paper IV with the ”one-planet-in-a-disk” prescription. In these previous models, the formation of gas giant planets require a substantial reduction in the pace of their cores’ type I migration (with ). This result reflects the requirement for not only the formation but also retention of sufficiently massive cores (§3.6) to efficiently accrete gas. The present prescription bypasses this competition between type I migration and mass growth because it takes into consideration the possibility that multiple embryos may undergo convergent migration, dynamical instability, orbit crossing and close encounters. Coalescence of two or more embryos (with relatively low masses and long type I migration time scales) promote sporadic, substantial mass increases. The merger products can impulsively acquire a critical mass for the onset of gas accretion before they undergo extensive migration.

A direct comparison between Figure 7c and Figure 5 in Paper IV show that, with , the new prescription indeed leads to much more prolific production of gas giants. Some of these gas giant planets form during the advanced stages of disk evolution after 90% of the disk gas been depleted. In these low density environment, the eccentricity of embryos may be excited by their perturbation on each other but not effectively damped. Eventually they cross each other’s orbit, merge into sufficiently massive cores to accrete residual disk gas and evolve into gas giant planets.

In contrast to the “one-planet-in-a-disk” models, we can determine, with the new prescription, planets’ velocity dispersion from their mutual interaction in multiple-planet systems. From this kinematic information, we obtain the and distributions in Figure 7a and c. The ”V” shape centered at a few AU in the distribution is a signature of strong scatterings by the most massive planets in the system which preferentially form near the snow line (at a few AU’s). Less massive planets are either scattered inward, outward, or collide with the gas giant planets. The left/right branch of the “V” corresponds to the loci of planets which have been scattered inward/outward. These inwardly/outwardly scattered planets preserve their apo/peri astrons to be orbital radii of the gas giant planets. Finally, the distribution shows that eccentricity distribution is uniform or rather increases with . If we consider only massive planets that are currently observable by radial velocity survey, eccentricity distribution clearly increases with as shown below.

Figure 7: The asymptotic kinematic distributions of a simulated population of planets in multiple-planet systems. For these fiducial models, we adopt yr and and . The eccentricity versus semimajor axis distribution is plotted in panel (a). The distribution is plotted in panel (b). The distribution is plotted in panel (c). The green, blue and red symbols represent rocky, icy and gas giant planets, respectively.

We compare our simulated population-synthesis models with observational data. In Figures 8, we plot observed (a) , (b) , and (c) distributions in the left, middle, and right columns respectively. The observational data are displayed in the top row. These data include only those planets which were first discovered by radial velocity surveys. We excluded planetary candidates which were first discovered by transit surveys, regardless whether they were confirmed by follow-up radial velocity observations. The main rationale for this selection is that the detection probability by transit surveys is highly biased to short-period planets.

The simulated population of multiple-planet systems is plotted in the bottom row. Here refers to with a distribution of inclination angle between the line of sight and the planets’ orbital angular momentum axis that assumes random orientations of orbit normals of the planets. Some of these planets have low-mass and large semimajor axis. They are below the present-day detection capability. In order to take into account the limited radial velocity precision and data base line, we selected a subset of data in which the predicted reflective motion of the host stars is larger than 1 m/s and the planets’ orbital period is shorter than 10 years. These “detectable planets” are plotted in the middle row.

Comparison between top and middle panels in Figures 8a and b shows that the theoretical models match well the observed kinematic distributions. Both observations and population synthesis models show that increases with and . In the next subsection, we will show that these correlations are robustly reproduced, almost independent of model parameters.

Most gas giants have semimajor axis AU. Their eccentricities are not affected by the dissipation of tidal perturbation induced by the host stars on the planets. Detailed comparison of the top and mid row of the middle column indicate that the fraction of high-eccentricity () planets in the theoretical model is low () compared with that () in observed data. This discrepancy may be observational uncertainties which generally tend introduce over-estimates for planets’ eccentricities. It may also be due the secular perturbations between gas giants which has not yet been incorporated in our prescription. In multiple-planet systems with significant angular momentum deficit, secular interaction can induce large-amplitude angular momentum exchange, eccentricity modulation, and long-term orbital instabilities. A treatment of secular perturbations will be constructed and presented elsewhere.

In this paper, we focus primarily on the planets’ eccentricity distributions. Comparison between top and middle panels in Figure 8c show that our prescription qualitatively reproduces the observed correlations in the distribution. Whereas the observational data clearly show an over density of gas giants at AU, their semimajor axis distribution in the simulated population synthesis models is relatively uniform on a log scale. The sensitive dependence of the distribution on the parameters of population synthesis models will be discussed in detail in the next paper.

Figure 8: Comparison of theoretical results in a nominal case with observational data: (a) , (b) , and (c) distributions. The top row are the data obtained from radial velocity surveys. The planets that were discovered transit survey and confirmed with follow-up radial velocity observation are excluded. The bottom row are the distributions of planets in systems simulated with the updated population-synthesis prescription. Panels in the middle row represent only the simulated observable planets, i.e., those with radial velocities greater than 1 m/s and periods less than 10 years.

5.3 Correlations between eccentricity and mass or semimajor axis

The origin of the correlation of eccentricity and semimajor axis can be traced back through the procedures of our theoretical model. In the observed distribution, is apparently lower at AU than at AU (the top panel in Figure 8b). This correlation has been attributed to the dissipation of stellar tidal perturbation inside the planetary envelope which leads to the orbital circularization of close-in planets. This tidal effects has not yet been implemented in the population synthesis models. Yet, this correlation is also very well established in the simulated results. In Eq. (3), we showed that the maximum eccentricity excited by close scattering between gas giants is , because the two-body surface escape velocity () is independent of whereas the Kepler velocity () is . Consequently, in the deep potential near the host stars, it is difficult for close scattering to excite high eccentricities. Based on the above prescription, the maximum we obtain at AU is 3 times smaller than that at AU in the simulated models (the middle panel in Figure 8b).

Multiple massive giants are preferentially formed in relatively massive disks (see discussions in §4). These dynamically packed systems are more prone to dynamical instabilities, orbit crossing, excitation of high eccentricities and cohesive collisions. A correlation between high eccentricity and planetary mass is naturally expected. In order to verify this scenario, we plot in Figure 9, and of gas giants as a function of in the fiducial models (in which yr and ) of Figure 7. This figure clearly shows that, in the range of , both the mass of gas giants and the fraction of high planets increase with . Based on this inference, we can infer the distribution of disk masses from the observed distribution. Raymond et al. (2010) performed N-body simulations of multiiple giant planets with different masses and found a similar correlation between high eccentricity and planetary mass, although their initial conditions are artificial. Note that Figure 9 only includes gas giants at AU. As shown in Figures 7a and c, the kinematic distribution of gas giants with AU have two (high and low ) components. We discuss this issue in the next subsection.

Figure 9: The distributions of and of gas giants with AU as a function of in the fiducial models of Figure 7. The thick and thin symbols represent giants with and , respectively (In the color version, red and orange symbols).

The above discussion indicates that the and correlations are the natural outcome of close scatterings and the formation of gas giants. It also shows that these correlations are insensitive to other physical parameters such as the rate of type I migration rate () and the magnitude of the gas accretion timescale (). These independences are reflected in Figures 10 and 11 for simulated populations with different values of and . Although the number of gas giants per host star is lower among populations with relatively high and long , the trend that is lower at AU than at AU persists.

This trend is less conspicuous in the limit of relatively small ( in Figures 10) and gas accretion time scale ( years in Figure 11). In these limits, the prolific formation of gas giants and the high frequency of their close scattering leads to a prominent ”V” feature in the distribution. This feature also affects the distribution of observable planets (displayed in the upper rows of these figures). The lack of any obvious evidences of this “V” feature in the observational data (Figure 8a), places a limit on the efficiency of gas giant formation and constraints on the magnitude of and . In the fiducial population synthesis model (with and yr), the fraction of solar-type stars that harbor gas giants with AU and (cool jupiters) is 21%. This fraction for both the models (47%) and the years models (42%) are both much higher than that inferred by radial velocity surveys.

Figure 10: Dependence on of the simulated distribution. The lower panels are distributions of systems of planets generated with our Monte Carlo simulation. The upper panels contain only the observable (with radial velocities greater than 1 m/s and periods less than 10 years) simulated planets. These observable criteria correspond to the technical limitation of current radial velocity surveys. The left, middle and right panels display results for models with type I migration coefficient and 0.3, respectively. The symbols are the same as those in Figures 8.
Figure 11: The distribution of simulated models with various values of . The left, middle and right panels include results generated with the parameter for envelope contraction timescale (for ) and years, respectively. The meaning of upper and lower panels are the same as in Figures 10.

Figures 12 and 13 show the dependence of the distribution on and , respectively. Although more massive planets, in general, tend to have larger , this trend is less conspicuous for models with and years. In these limiting cases, the formation of gas giants is affected by the rapid migration of the cores and slow gas accretion. These restrictions are more severe for formation of the relatively massive gas giants. Despite these extreme cases, the correlations in and distributions are well established in the population synthesis models for a wide range of parameters. The general reproduction of the observed correlations suggests that these properties are associated with close scattering among gas giants.

Figure 12: Dependence on of the predicted distribution. The meaning of panels are the same as in Figures 10.
Figure 13: Dependence on of the predicted distribution. The meaning of panels are the same as in Figures 11.
Figure 14: Same as Figure 7 except for use of non-isothermal type I migration formula instead of isothermal one with .

In Figures 14a, b and c, we show the , and distributions with the non-isothermal type I migration formula. Because in inner (optically thick) disk regions, migration is outward for some range of (Kretke & Lin, 2012), many gas giants survive even without any reduction in amplitude of the migration velocity. While the distribution is affected by the different migration formula, the and distributions are similar to those obtained by the isothermal migration formula and they are within a range of variations of the results with different values of and in figures  10, 11, 12 and 13. Therefore, here we do not go into the details of the non-isothermal migration formula and focus on the eccentricity distributions obtained with the isothermal migration formula and a scaling factor , as discussed in §3.3.

5.4 Formation of distant large gas giants in circular orbits

The population synthesis models also generate a population of massive () gas giants with large semimajor axis (AU) (Figure 7c). In principle, these planets can be observationally detected during their infancy with direct imaging searches. In the fiducial models, the fraction of host stars that have these planets is and most of them have low-eccentricity () orbits (Figure 7a). If these planet were grown to gas giants near the snow line and scattered to these large distances, they would preserve their periastron distance. For example, planets formed interior to 30 AU may be scattered to attain a semimajor axis AU with an . Although the effect of eccentricity damping is included in our prescription for disk-planet interaction, its efficiency is generally too low to account for the low eccentricity of distant, massive planets generated in these population synthesis models (step 3-a of Appendix A1).

A close inspection of our results yields an alternative path for the formation of distant gas giant planets. In systems which contain multiple gas giants, the rapid gas accretion of first generation gas giants destabilize the orbits of nearby residual embryos. Some embryos are scattered to large distance where the surface densities of both residual planetesimals and gas are relatively low. Because the scattered embryo masses are usually well below the local disk mass, even the relatively low surface density gas damps the eccentricity of the scattered embryos to the level below 0.1 (also see the discussion below).

A reduction in the planetesimal accretion rate at these distance also lowers the critical core mass for the onset of gas accretion (see Eq. 23). In extended protostellar disks, there is adequate supply of gas supply, even during the advanced stages of disk evolution to enable the circularized planets to acquire masses larger than that of Jupiter.

The ratio of planets’ asymptotic mass to that of their host stars is a monotonically increasing function of their natal disks’ aspect ratio which is also an increasing function of (see Eqs. 27 and 28). The simulated populations confirm that the masses of the long-period gas giants are robustly correlated with in the fiducial (Figure 7c) and other models. If this trend is confirmed by observation, it provides evidence to support the scenario that long-period, low-eccentricity, massive gas giants are the byproducts of gas accretion onto outwardly scattered cores.

During the course of mass accretion, the planets also acquire angular momentum of the disk gas. Although the scattered embryos have high eccentricity initially, they gain angular momentum through planet-disk tidal interaction. Our prescription for the eccentricity damping may be less efficient than that, due to supersonic relative velocity between the core and disk gas (Ostriker, 1999; Papaloizou & Larwood, 2000; Muto et al., 2010), although it nonetheless may lead to the circularization of planets’ orbits. In addition to the tidal torque, E. Thommes (2010, private communication) used a hybrid N-body and 2D hydrodynamical (FARGO) scheme to demonstrate that during their eccentric excursion, scattered embryos take longer time to pass through their apoastrons and they tend to accrete gas from that region which has relatively large specific angular momentum. We anticipate that as the embryos evolve into gas giants, their specific angular momentum becomes comparable to that of the disk gas near their apoastron and their orbits become circularized with a radius close to their apoastron radius. The detailed calculation will be shown in a separate paper (Kikuchi, Higuchi & Ida, in preparation).

Through this process, multiple cores can be scattered outward to initiate gas accretion. We find that in the fiducial models, the fraction of systems with two, three, and four giants with low-eccentricity orbits (such as the solar system) is 0.4%, 0.07%, and 0.04%, respectively. The orbital configuration of systems of multiple gas giants with relatively large masses, low eccentricities, and long periods are consistent with the observed data obtained from the direct imaging initiatives.

We have not introduced an outer boundary for planets’ natal disks. If most protostellar disks are truncated either by photoevaporation or magnetic braking, the outwardly scattered cores would not have access to gas beyond the truncation radius. An abrupt decline in the observed same major axis distribution of gas giant planets could be used to place constraints on the structure and evolution of their natal disks.

In the models with , the fraction of systems with massive, low-eccentricity, and long-period gas giants is similar to that of the fiducial models (with ). But this fraction is reduced to in the models with . This dichotomy suggests that type I migration of cores is a major cause for them to be scattered by a giant planet. In order to match with the frequency of planets discovered through direct imaging observation, type I migration cannot be neglected.

6 Summary and Discussions

Radial velocity and transit surveys indicate that multiple-planet systems are common around solar-type stars. We adopt a new paradigm that all planets are formed in families (as in the case of the solar system) and their asymptotic properties are determined by their dynamical interaction with each other as well as with their natal disks.

In order to incorporate this scenario into our population synthesis models, we upgrade an existing numerical population synthesis scheme which was developed in Paper I-IV for idealized “one-planet-in-a-disk” models. The original prescriptions were constructed and applied to simulate the semimajor axis and mass distributions of extrasolar planets, calculating the growth and migration of each planet, independently of any other planets. Similar approaches have also been adopted by other authors (Mordasini et al., 2009a, b; Alibert et al., 2011).

In Paper VI, we first introduced a modification to partly include the gravitational interactions and collisions between rocky planetary embryos. In this paper, we generalize an efficient and robust prescription to take into account close scattering between all planets. We calibrate our prescriptions for scattering among gas giants with the results of comprehensive N-body simulations. For example, in the aftermass of close scatterings between gas giants, ejection commonly occurs. Our prescriptions reproduce well the distributions of eccentricity and semimajor distributions of retained giants obtained from the N-body simulations.

With this upgraded tool, we can adopt a Monte-Carlo approach to simulate a large sample of population synthesis models. In this paper, we mostly investigate the effects of scattering by giant planets on asymptotic architecture of planetary systems and their eccentricity distributions. We focus on the correlations between the eccentricity of gas giants with their mass and semimajor axis.

Observational data obtained through radial velocity surveys clearly shows that the mean eccentricity increases with planetary mass (at least for ) and with semimajor axis (at least for AU). Our new population synthesis models show that in relatively massive disks, several massive giants may form and interact with each other intensely. Dynamical instabilities lead to eccentricity excitation, orbit crossing, close encounters, collisions, and ejection. But in protostellar disks with moderate masses (comparable to the surface density distribution of the minimum mass nebula model), one or two relatively small-mass giants are formed. They commonly avoid dynamical instability and preserve their low initial eccentricities. The observationally inferred correlation between gas giants’ eccentricity and mass is robustly reproduced by our models.

Both gas giants and rocky/icy planets form prolifically in massive disks. But the low-mass planets are mostly ejected or accreted by the gas giants. Moderate planetary systems in which rocky/icy planets preserve their modest eccentricity and coexist with one or two sibling gas giants are preferentially formed in disks with modest masses (analogous to the formation of the solar system in a minimum mass solar nebula). Multiple, short-period super-Earths, which are commonly found around solar-type stars, are also retained only in the systems formed from such modest-mass disks but not in the massive disks.

The observed correlation between eccentricity and semimajor axis is also reproduced by the population synthesis models. After the onset of dynamical instability and orbit crossing, planets’ velocity dispersion is determined by close scatterings. The maximum recoil speed is comparable to the surface escape velocity of the dominant perturber, which is independent of semimajor axis. Since mean eccentricity is given approximately by velocity dispersion scaled by the local Keplerian velocity, the mean eccentricity is proportional to square root of semimajor axis. Although we have not included the effect of tidal circularization, which preferentially damps the eccentricity of short-period planets, the correlation of eccentricity and semimajor axis can also be produced by close scatterings alone.

Recently, direct imaging surveys discovered several systems which contain massive () gas giants in distant (AU) orbits with small eccentricities (). Some of these planets are members of multiple planets systems. We found analogous planets among a non-negligible fraction (a few %) of systems in our population synthesis simulations. In our models, the cores of these gas giants were originally formed well interior to their present-day semimajor axis (with AU). They were scattered outwardly by first generation gas giants to the outer regions of their natal disks. The surface density of field planetesimals and their collision frequency with the cores rapidly decreases with the disk radius. A reduction in the planetesimal bombardment rate also decreases the critical core mass needed for the onset of efficient gas accretion. In the outer regions of relative massive disks, there is an adequate supply of gas to enable the cores to accrete gas even during the advanced stage of disk evolution. The high eccentricity of the scattered cores is damped by disk-planet interactions and accretion of gas in the outer disk where its specific angular momentum is relatively high.

The emergence of single massive giants can scatter multiple nearby cores, we also found systems with multiple massive gas giants in nearly circular distant orbits. Therefore, the core accretion model can indeed produce the systems that are similar to those discovered by direct imaging surveys. This origin predicts a clear correlation that the mass increases with semimajor axis for these giants by the following reason. The aspect ratio of the disks increases with their radius. At these large distances, gas giants continue to accrete gas until they acquire sufficient mass to open gaps in the disks near their orbits. Because the critical mass for gap opening is higher for larger orbital radius, the correlation is established. It is of interest to check if observed data also indeed exhibits such a trend.

ACKNOWLEDGMENTS. This work is partially supported by NASA (NNX08AM84G), NSF (AST-0908807), JSPS (23103005) and a University of California Lab Fee grant.

CORRESPONDENCE should be addressed to S. I. (

Appendix. Prescription for eccentricity excitation and ejection of giant planets as a result of orbital instability

In our simulations, giant planets and rocky/icy planetary embryos co-exist. We separate the treatment of planets’ dynamical interactions into those between (i) giant planets, (ii) embryos, and (iii) a giant planet and an embryo.

We consider two-giants scatterings in a disk environment with residual gas. We adopt an assumption that three-giants scatterings mostly occur after the disk has been depleted in a gas-free environment (see §2). The prescriptions for scattering between two gas giants and three gas giants are described in A1 and A2 respectively. In A3, we present prescriptions for dynamical interactions in general cases in which giant planets and rocky/icy planetary embryos co-exist.

A1. Two giants

The prescriptions for scattering of two giants with mass and are as follows:


Specify initial conditions: We assume that orbital instability between two gas giants occurs when their orbital separation becomes smaller than