{'fraction_non_alphanumeric': 0.03509440600445643, 'fraction_numerical': 0.02018588014542043, 'mean_word_length': 4.324877058777613, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
From planetesimals to terrestrial planets: N-body simulations including the effects of nebular gas and giant planets

We present results from a suite of N-body simulations that follow the accretion history of the terrestrial planets using a new parallel treecode that we have developed. We initially place 2000 equal size planetesimals between 0.5--4.0 AU and the collisional growth is followed until the completion of planetary accretion (>100 Myr). All the important effect of gas in laminar disks are taken into account: the aerodynamic gas drag, the disk-planet interaction including Type I migration, and the global disk potential which causes inward migration of secular resonances as the gas dissipates. We vary the initial total mass and spatial distribution of the planetesimals, the time scale of dissipation of nebular gas, and orbits of Jupiter and Saturn. We end up with one to five planets in the terrestrial region. In order to maintain sufficient mass in this region in the presence of Type I migration, the time scale of gas dissipation needs to be 1-2 Myr. The final configurations and collisional histories strongly depend on the orbital eccentricity of Jupiter. If today's eccentricity of Jupiter is used, then most of bodies in the asteroidal region are swept up within the terrestrial region owing to the inward migration of the secular resonance, and giant impacts between protoplanets occur most commonly around 10 Myr. If the orbital eccentricity of Jupiter is close to zero, as suggested in the Nice model, the effect of the secular resonance is negligible and a large amount of mass stays for a long period of time in the asteroidal region. With a circular orbit for Jupiter, giant impacts usually occur around 100 Myr, consistent with the accretion time scale indicated from isotope records. However, we inevitably have an Earth size planet at around 2 AU in this case. It is very difficult to obtain spatially concentrated terrestrial planets together with very late giant impacts.


Introduction
Accretion models for the formation of the terrestrial planets need to reproduce the following properties of our solar system. (1) The spatial mass distribution of the terrestrial planets is radially concentrated. Namely, the Earth and Venus contain most of the mass inside of Jupiter's orbit, whereas the masses of Mars and Mercury are an order of magnitude smaller than the Earth's mass (M ⊕ ) and the total mass of asteroids is only 5 × 10 −4 M ⊕ . The orbital separation between Earth and Venus is only 0.3 AU. (2) The orbits of the Earth and Venus are nearly circular and coplanar even though orbital excitements are expected in the late stage of accretion where large protoplanets experience mutual collisions. This probably indicates that there were some damping mechanisms to reduce their orbital eccentricities and inclinations near the end of their formation processes. (3) Isotope records such as W/Hf and U/Pb for the Earth and Moon indicate that the Moon-forming giant impact occurred ∼ 100 Myr after the beginning of the solar system (Touboul et al., 2007;Allégre et al., 2008). Some accretion models can reproduce two of the three properties but there is no accretion model that satisfies all three properties.
In particular, the dynamical shake-up model of Nagasawa et al. (2005) and Thommes et al. (2008) is the most successful model so far in reproducing both the small orbital eccentricities and the radial mass concentration. Their model initially places a few tens of protoplanets embedded in a gas disk and takes into account the effect of the gas giant planets with orbital eccentricities comparable to or slightly larger than the present values (the present eccentricities of both Jupiter and Saturn are ∼ 0.05). The secular resonance with Jupiter migrates from ∼4.0 to 0.6 AU with dissipation of the gas disk and sweeps protoplanets toward the current location of the Earth. Mutual impacts between protoplanets occur even when a large amount of gas remains, and the remnant gas disk reduces the orbital eccentricities of formed Earth-like planets.
If the formation of the Moon occurred very late, as indicated by recent isotope measurements (Touboul et al., 2007), significant damping due to remnant gas may be problematic for the following reasons. Significant depletion of gas already at ∼ 10 Myr is suggested from direct measurements of line emissions of circumstellar disks (Pascucci et al. 2006) and from Hα emission profiles of central stars due to accretion shock (Muzerolle et al., 2000), although we cannot rule out the possibility that our protosolar nebular was exceptionally long-lived. Another potential problem for long-lived gas disks is depletion of the solid component due to Type I migration, which is neglected in the model of Nagasawa et al. (2005) and Thommes et al. (2008). Once a protoplanet becomes as massive as Mars in a gas disk, the gravitational torque from spiral density waves launched by the protoplanet causes rapid inward migration of the protoplanet, referred to as Type I migration (Goldreich and Tremaine, 1980;Ward, 1986Ward, , 1997Tanaka et al., 2002). McNeil et al. (2005) and Daisaka et al. (2006) conducted N-body simulations taking Type I migration into account and found that the decay time scale for a gas disk needs to be ∼ 1 Myr in order to retain sufficient solid mass for the formation of Earth-size planets.
Some remnant planetesimals are likely to have existed even at 100 Myr particularly in the asteroidal region. 1 Several authors have conducted N-body simulations starting with protoplanets and remnant planetesimals (Chambers, 2001;Chambers and Cassen, 2002;O'Brien et al. 2006;Raymond et al. 2006Raymond et al. , 2009. Their results show that final eccentricities of planets tend to be lower with increasing total mass of remnant planetesimals. Even with the same total mass of planetesimals, the damping is more effective when the size of an individual planetesimal is smaller and planetesimals are placed uniformly over both the terrestrial and asteroidal regions rather than in the asteroidal region only (O'Brien et al., 2006;Raymond et al., 2006). An eccentricity of Jupiter comparable to or slightly larger than the present value is preferred in regard to the radial mass concentration (Chambers and Cassen, 2002;O'Brien et al., 2006;Raymond et al., 2009), although none of simulations have succeeded in reproducing the very high radial mass concentration of the current terrestrial planets. The spatial and size distributions of remnant planetesimals are arbitrarily assumed in these simulations. In order to clarify these distributions one needs to follow the growth of planetesimals and protoplanets in the runaway growth stage. In addition, the above simulations neglect the effects of gas. The protosolar nebular must have existed at least until the formation of Jupiter, and probably survived somewhat longer. Even a rapid dissipation of gas of the order of 0.1 Myr in the presence of Jupiter strongly affects the orbital evolution of asteroidal bodies (Nagasawa et al., 2000).
In the present study we perform N-body simulations starting with a wide planetesimal disk and continue until the end of the accretion phase. Since we directly calculate the runaway growth of the protoplanets, we can obtain reasonable spatial and size distributions of remnant planetesimals in the late stage. We take into account the gas giant planets and all the possible effects of gas in laminar disks: aerodynamic gas drag, disk-planet interaction including Type I migration, and global disk potential. In a previous study, Morishima et al. (2008) also conducted N-body simulations that followed the evolution of a protoplanetary disk starting with planetesimals only until the end of accretion. However, their simulations started from radially narrow annuli (widths of 0.3-0.5 AU) and neglected any external forces. In Sec. 2, we introduce our new N-body code, review the effects of nebular gas, and outline our choice of simulation parameters. In Sec. 3, we show time evolution of two contrasting examples, present the results of all simulations and compare them with the solar system. In principle, none of our simulations succeed in satisfying all the three constraints listed in the first paragraph of this section. In Sec. 4, we discuss possible solutions for this problem. In Sec. 5, we give our conclusions.

N-body code
Any N-body code consists of a gravity solver to calculate the mutual gravity force and an integrator to update positions and velocities for a next time step. Existing N-body codes applied for the runaway growth stage have good gravity solvers, such as tree methods (Richardson et al., 2000;Stadel, 2001) or special hardware such as HARP or GRAPE (Kokubo and Ida, 1996). On the other hand, N-body codes applied for the late stage have good integrators, such as SyMBA (Duncan et al., 1998) and Mercury (Chambers, 1999), so that one can take a large time step size, while solving the mutual gravity by direct N 2 calculation. The new N-body code we have developed has both a good gravity solver and a good integrator, allowing us to follow accretion stages of planetary systems with large N and over a large number of orbits. The code is still being upgraded and details of the complete version will be reported elsewhere. Here we briefly describe the version of the code used for simulations in this paper.
The gravity solver we use is a parallel-tree method, PKDGRAV (Stadel, 2001). The original version of PKDGRAV uses a leap-frog integration scheme, where individual particles are updated in hierarchical time step blocks, each block with a factor of two smaller time-step according to the local dynamical time of its particles. Richardson et al. (2000) added some functions for planetary accretion to the code and conducted simulations for the runaway growth stage. In the tree method, all particles are divided into hierarchal tree-structured groups. Then, the gravity forces from distant particle groups are calculated using their multipole expansions. The zeroth-order multipole expansion corresponds to the gravity from the mass center of a particle group. In order to keep good accuracy, we include multipole expansions up to the fourth order. The computer time of gravity calculation using this gravity solver is only proportional to N log N. We find a significant speed up in the gravity calculation with our tree method as compared with the direct N 2 calculation when N > 200. In addition, the tree method is suitable for parallel computing, and this results in further speed up when N > 1000. Since the number of particles decreases with growth of planets, we switch calculation methods from parallel-tree to serial-tree, and to serial-direct N 2 calculations. The error in the total energy due to the tree method appears as a random walk, and thus is proportional to the square root of the number of time steps. Therefore, with some short additional simulations, we can estimate the accumulated error after a long term evolution of a system. The accuracy of the method is controlled by the opening angle for tree cells and we set it to be 0.5 in the present paper. We confirm that the expected error due to the gravity calculation is sufficiently smaller than the error due to the integrator's, which rises most notably during close encounters.
The integrator implemented in the code is a mixed-variable symplectic (MVS) integrator, SyMBA (Duncan et al., 1998). In MVS integrators such as SyMBA and Mercury, the Hamiltonian is split into a part which describes the Keplerian motion around the central star and another part which describes the perturbation from other bodies. The former part corresponds to updates of positions and velocities along Kepler orbits. This can be done nearly analytically. The latter part corresponds to velocity changes due to the perturbation from other bodies, which we calculate using the tree method described above. These types of integrators are very accurate and allow us to take very large time steps for long term simulations. We use a time step size, ∆t, of 6 days. This is comparable to those used in previous high-resolution N-body simulations using MVS integrators (O'Brien et al., 2006;Raymond et al., 2006). Practically, the time step size needs to be chosen so that a fraction of particles in close encounters defined below is small enough.
In close encounter search, we first search for 10-20 nearest neighbors of each particle using the same algorithm in the tree SPH code, Gasoline (Wadsley et al., 2004). With the relative velocities and positions at the beginning and end of a time step, the closest distance of a pair during the step is calculated using a third order interpolation (Eq. (11) of Chambers, 1999). In order to judge whether a particle is in a close encounter or not, we first define the critical distance for each particle. The critical distance is defined as the larger of (1) n 1 R H or (2) n 2 v ep ∆t, where n 1 and n 2 are numerical factor (we use n 1 = 4.5 and n 2 = 1.0), R H is the Hill radius, and v ep = v kep √ e 2 + i 2 is the epicyclic velocity of the particle with v kep , e, and i being the Kepler velocity, eccentricity, and inclination, respectively. If the closest distance of a pair is smaller than the sum of the critical distances, we regard this pair to be in a close encounter. In the original SyMBA code (Duncan et al., 1998), only the first condition is used. However, the second condition is necessary if the relative velocity is much larger than the mutual escape velocity; such encounters usually occur near resonances of the giant planets or between small planetesimals stirred by large protoplanets. In the Mercury code (Chambers, 1999), the second condition uses the largest v kep in all particles instead of each body's v ep which we apply. Although this criterion of the Mercury code gives better energy conservation (for a same n 2 ), it also results in too many pairs being identified as involved in close encounters during the runaway growth stage, where the number density of particles is large. If the number density is not large, as is the case in the late stage, use of the Mercury's criterion is recommended (although we do not use it even in the late stage in this paper).
All the above computations are performed in parallel (if N > 1000), while particles in close encounters are collected on a single cpu during the domain decomposition; the procedure which distributes the workload, and hence the remaining particles, amongst processors. Then, multistep time stepping is done until all particles proceed to the end of the original time step, following the procedure described in Duncan et al. (1998). The closest mutual distance in a substep is always estimated by interpolation using the information of both the beginning and the end of the step. This guarantees time symmetry of the integrator which is essential to avoid secular errors. The computational time for the substep procedure is small enough as it is simply proportional to the number of particles in close encounters. We set the ratio of the cut-off radius to that of the next order to be 2.08 and the number of next order substeps in the current substep to be 3 after Duncan et al. (1998).
We assume that two bodies merge when they partially overlap. We assume the density of particles to be 2 g cm −3 . Since we do not apply the artificial enhancement of radii Ida, 2002, Morishima et al., 2008), the accretion time scale we will show is realistic. A particle is assumed to collide with the Sun when either the heliocentric distance is less than 0.1 AU or the semimajor axis is less than 0.25 AU. A particle is assumed to escape from the system when the heliocentric distance is larger than 100 AU. The energy error after completion of accretion is typically ∆E/E = 10 −5 -10 −3 . Note that energy changes due to merging or gas forces can be explicitly calculated and are subtracted in the above error. The largest portion of the error usually accumulates after the runaway growth stage where a large number of excited small particles experience mutual encounters. Therefore, the energy error primarily depends on how long small particles survive.

Effects of gas
In this section, we review the forces due to the nebular gas. The acceleration due to the gas is incorporated in the integrator by placing the half kicks before and after the Kepler drift operation. This symmetric placement once again avoids secular errors coming from the inclusion of gas forces.

Aerodynamic drag
The aerodynamic drag force per unit mass (or acceleration) is given by (Adachi et al., 1976) where m and r p are the mass and the radius of a particle, c D is the numerical coefficient (assumed to be 2 in the present paper), ρ gas is the gas density, and u rel = u − u gas is the velocity of a particle relative to the gas velocity. In the cylindrical coordinate (r, θ, z), the gas velocity in the θ component is obtained by the force balance in the r component as where c is the isothermal sound velocity, T is the gas temperature, G is the gravitational constant, and M ⊙ is the mass of the central star (assumed to be the solar mass). The left hand side of Eq.
(2) represents the centrifugal force while the first, second, and third terms of the right hand side represent the pressure gradient, the gravity of the central star, and the self-gravity of the gas disk (Eq. (16)), respectively. Assuming the gas is isothermal in the z direction and neglecting the self gravity of the gas, the gas density is approximately given by where Σ gas is the surface density of the gas and h = c/Ω kep is the scale height, with the Keplerian frequency Ω kep . If we assume Σ gas ∝ r −p , T ∝ r −q , we obtain, We adopt the temperature profile of an optically thin disk, T = 280 (r/1 AU) −1/2 K (q = 1/2), after Hayashi et al. (1981). We will explain how Σ gas and p are given in Sec. 2.3.
The growth rate of a protoplanet is proportional to e −2 , where e is the orbital eccentricity of planetesimals, and the balance between viscous stirring by protoplanets and damping due to gas drag gives e ∝ m 1/18 (Eq. (25) of Kokubo and Ida, 2000). The initial planetesimal mass m 0 in our simulations is ∼ 10 25 g. On the other hand, the size of planetesimals m GI formed by gravitational instability is estimated to be several order of magnitude smaller (Goldreich and Ward, 1973). Therefore, the growth time scale of protoplanets will be significantly overestimated in our model. In order to reproduce the growth rate of protoplanets surrounded by small planetesimals (∼ m GI ), we artificially enhance the gas drag force on planetesimals by a factor of (m/m ′ ) 2/3 in N-body simulations if the particle mass m is smaller than m 1 (> m 0 ). Here the mass m ′ is given as With this scheme, a bunch of small planetesimals of mass m ′ are represented by a single super planetesimal of mass m (e.g., McNeil et al., 2005). The mass m 1 should be as small as possible in order to avoid artificial clean up of large planetesimals by protoplanets. We use m 1 = 0.01 M ⊕ as a compromise, and m GI = 10 19 g. The migration and damping time scales with the artificial enhancement of the gas drag are shown in Fig

Tidal interaction between a protoplanet and a nebular disk
We use the formulation of Tanaka et al. (2002) for Type I migration rate. For the damping rates of orbital eccentricity e and inclination i, we use the formulation given in Tanaka and Ward (2004). We make corrections in these formulations for the case of large e and i after Papaloizou and Larwood (2000) and Cresswell et al. (2007). Tanaka et al. (2002) derived the torque due to tidal interaction between a protoplanet and a three-dimensional baroclinic disk, for the case of e = i = 0. The decay time scale of the semimajor axis a is given as where f tid1,θ is the acceleration in the azimuthal direction. Tanaka and Ward (2004) extended their work to the case of nonzero e and i (≪ h/a) and obtained the rate for damping of e and i due to the tidal interaction. The r, θ, and z components of the damping force per unit mass are given as (a typo was corrected; see Ogihara et al., 2007) where characteristic time of the orbital evolution, τ wave , is given by In Eqs. (7)-(9), the velocity components of a particle are given as u kep − r f glob,r (r, 0)) 1/2 is the Keplerian velocity of solid bodies including the gravity of the gas disk. Note that all the components of the damping force become zero when e = i = 0.
The formulation of Tanaka and Ward (2004) is limited to small e and i (≪ h/a), or that the planet motion relative to the gas remains subsonic. From a two dimensional linear theory, Papaloizou and Larwood (2000) found that the time scales τ tid1 and τ wave (Eqs. (6) and (10)) increase with increasing e by factors of g tid1 and g wave , respectively: Here the original x used in Papaloizou and Larwood (2000) is x PL = ae/h. On the other hand, we use a modified value as in order to reproduce similar results of simulations of Cresswell et al. (2007) (see the next paragraph). The evolution of the orbital eccentricity with the corrections for large e is given as This equation corresponds to Eq. (45) of Tanaka and Ward (2004), but the correction factor is multiplied by τ wave . From Eqs. (6), (11), and (14), the evolution of the semimajor axis due to the tidal interaction is given by, where the first term in the right hand side comes from exchange of orbital angular momenta while the second term represents the migration rate due to energy dissipation. Cresswell et al. (2007) conducted hydrodynamic simulations for migration of a 20 Earth mass planet with nonzero e and i embedded in a gas disk. They found good agreement with Tanaka et al. (2002) and Tanaka and Ward (2004) for small e and i (< 0.1), while the damping time scale e/ė is proportional to e 3 for large e as predicted by Papaloizou and Larwood (2000) (Eq. (12)). They also showed the same dependence of the inclination decay rate on i, so we first simply replace e by (e 2 + i 2 ) 1/2 in Eq. (13). With the original x PL , Eq. (15) gives outward migration for ae/h > 1.1. Cresswell et al. (2007) also found that the sign of the torque changes to positive (g tid1 < 0) when ae/h > 1.6, but did not find any indication of outward migration. This means that the absolute value of the second term in the right hand side of Eq. (15) is always larger than the first term for large e. It is unclear whether outward migration really occurs or not and we need to wait for future studies. In the present paper, we use the modified x (Eq. (13)) to reproduce similar dependences ofė andȧ on e and i to those obtained in Cresswell et al. (2007). The migration and damping time scales (e/ė and a/ȧ) with various e are shown in Fig. 1. Although there are some differences between the time scales shown in Fig. 1 and those from Cresswell et al. (2007), these do not result in significant differences in our N-body simulations unless the time scales become too large or negative. Practically, large protoplanets obtain large e when they are trapped in secular resonances. If outward migration really occurs, protoplanets can be removed from secular resonances which migrate inward.

Global nebular force and secular resonances
We calculate the gravity of the gas disk after Nagasawa et al. (2000). The r and z components of the nebular gravitational force are given as where K(k) and E(k) are elliptic integrals of the first and second kind, and k is given by These components of the force are initially tabulated in a (r, z) grid, and interpolated in simulations. A secular resonance occurs when precession rates of two orbiting bodies coincide (Heppenheimer, 1980;Ward, 1981;Nagasawa et al., 2000). We analytically calculate precession rates of planetesimals and the gas giant planets for e, i ≪ 1, following Appendix A of Nagasawa et al. (2000). Then, we obtain the locations of the secular resonances with Jupiter (ν 5 ) and Saturn (ν 6 ) (e.g., Fig. 2).

Input parameters
We adopt four different initial spatial distributions of planetesimals, four different dissipation time scales of the nebular gas, and four different types of orbits of Jovian planets. In total, therefore, we made 64 runs with different combinations of parameters. Each run takes one to three cpu months until completion of accretion.

Initial spatial distribution of planetesimals
We initially place 2000 equal-mass planetesimals between 0.5 and 4.0 AU. The initial surface density distribution of planetesimals is given by where Σ solid,0 is Σ solid at 1 AU. We adopt p = 1 and 2 in simulations. The initial total solid mass is given as We adopt M T = 5 and 10 M ⊕ . In the case of M T = 5 M ⊕ , Σ solid,0 = 6.1 and 10.2 gcm −2 for p = 1 and 2, respectively. For comparison, the minimum mass solar nebula (MMSN) by Hayashi et al. (1981) has Σ solid,0 = 7.1 gcm −2 , p = 1.5, and M T = 4.3 M ⊕ .

Dissipation time scale of nebular gas
In the present paper, we limit ourselves to the case where the surface density of gas dissipates exponentially in time and uniformly in space: where Σ gas,0 is Σ gas at 1AU and t = 0 and τ decay is the time scale for gas decay. We adopt four different decay times: τ decay = 1, 2, 3, and 5 Myr. We assume that Σ gas,0 = 2000 and 3366 gcm −2 for p = 1 and 2, respectively. The gas-to-solid ratio is constant between 0.5 AU and 4.0 AU in the beginning of simulations and is 329 and 165 for M T = 5 and 10 M ⊕ , respectively. For comparison, MMSN has Σ gas,0 = 1700 gcm −2 (Hayashi et al., 1981). We do not consider a gap around Jupiter's orbit in the present paper. The formation of a gap does not significantly alter the time of the sweeping secular resonance's passage in the gaseous region (Nagasawa et al., 2005). However, orbital evolutions of asteroidal bodies in a gas-free gap will be different, in particular, if a gap is very wide.

Orbits of Jovian planets
We introduce Jupiter and Saturn after τ decay from the beginning of a simulation, assuming that they form via core accretion (Bodenheimer and Polack, 1986;Ikoma et al., 2000;Lissauer et al., 2009). Rapid formation of Jovian planets via disk instability is unlikely for our solar system because efficient cooling for clump formation is only possible at r > 100 AU (Boley, 2009).
We define a J and e J to be the semimajor axis and the orbital eccentricity of Jupiter, and a S and e S to be the corresponding orbital elements for Saturn. We adopt four types of orbits for Jupiter and Saturn at their formation and use the same notation introduced in O' Brien et al. (2006) and Raymond et al. (2009), although we have slightly different orbits: • EJS ("Eccentric Jupiter and Saturn"). Jupiter and Saturn are placed on their current orbits: a J = 5.20 AU, e J = 0.048, a S = 9.55 AU, and e S = 0.056.
• CJS ("Circular Jupiter and Saturn"). These are the initial conditions used in the Nice model . Jupiter and Saturn are radially more confined with circular orbits: a J = 5.45 AU, e J = 0.0, a S = 8.18 AU, and e S = 0.0. The mutual inclination is 0.5 degrees.
• EEJS ("Extra Eccentric Jupiter and Saturn"). The same as EJS except a higher eccentricity of Jupiter: e J = 0.1.
• CJSECC. Jupiter and Saturn have the CJS semimajor axes with higher eccentricities: e J = 0.05 and e S = 0.05.
The EJS and EEJS simulations assume that Jupiter and Saturn formed at their current locations. On the other hand, the Nice model suggests that Jupiter migrated inward and Saturn migrated outward from their original CJS orbits via interactions with Kuiper Belt objects leading to a sudden change to the EJS orbits when Jupiter and Saturn crossed their 1:2 orbital resonance . Asteroidal and cometary bodies strongly perturbed by this event are considered to have caused the Late Heavy Bombardment at ∼ 3.85 billion years ago; well after formation of the terrestrial planets Strom et al., 2005). Initially higher e J and e S in the EEJS and CJSECC simulations than those in the EJS and CJS simulations, respectively, are considered because e J and e S decrease via interaction with planetesimals (Sec. 3).
We only take into account the gravity of other bodies and the global nebular force on motion of Jupiter and Saturn. We neglect any further interaction between a gas disk and the giant planets. Details of orbital evolution of a giant planet, which opens a gap in a gas disk, have a complicated dependence on disk properties, such as shape of a gap and disk mass, for which explicit and reliable formulae are not yet available. Nevertheless, some basic trends, which should be included in future studies, started to be clarified as follows. Analytic (Goldreich and Sari, 2003;Ogilvie and Lubow, 2003) and numerical ) studies indicate that eccentricity growth for a Jupiter-mass planet occurs together with the saturation of the corotation resonance in a low-viscosity disk. The migration time scale of a Jovian planet is much longer than the viscous diffusion time scale and even outward migration occurs, either when the gap is clean and the orbit is eccentric  or when the gap is not deep enough . A long migration time scale of Jupiter due to disk-sculpturing by Saturn is also pointed out by Morbidelli and Crida (2007).

Case with EJS orbits
Figures 2-4 show an example of the time evolution for the case with the present orbits of Jupiter and Saturn: snapshots on the plane of the semimajor axis versus the orbital eccentricity (Fig. 2), evolution of the total masses and the mean eccentricities of small and big particles (Fig. 3), and evolution of the masses and orbital excursions of surviving planets (Fig. 4). The gas [Figs. 2-4] decay time is t decay = 1 Myr, the initial total mass is 5 M ⊕ , the power-law index for the surface density is p = 2. In the present paper, we distinguish between a big particle, or (proto)planet, defined as a body with m > 2 × 10 26 g (∼ 0.03 M ⊕ ), and smaller bodies or planetesimals. This border mass is slightly smaller than Mercury's mass (3.3 × 10 26 g), and one order of magnitude smaller than the isolation mass at 1AU (Kokubo and Ida, 2002).
We follow growth of protoplanets without Jupiter and Saturn, until t decay from the beginning of the simulation. We first observe the mass loss (∼ 0.2 M ⊕ at 10 5 yr; see Fig. 3) by the migration of the smallest planetesimals due to gas drag at the inner most region. Protoplanets grow rapidly in the inner regions while many smaller bodies remain unaccreted in the outer regions. Protoplanets first tend to collide with smaller planetesimals because their orbital eccentricities are smaller. At 1 Myr, the smallest particles significantly deplete at r < 2 AU, while some larger planetesimals still survive there (Fig. 2). Once a planet reaches a certain size, it starts rapid inward migration (Type I migration). Since the growth time is shorter for smaller r, significant mass depletion first occurs in the inner disk. With increasing time, the mass depletion propagates outward (see Ida and Lin (2008) for more details). The mass loss due to Type I migration is easily recognized in the upper panel of Fig. 3, as we see stepwise depletion in the total mass. More than half of the total mass is lost by this mechanism until 3 Myr. The half of the total mass corresponds to the initial total mass inside of ∼1.5 AU. The three planets, which survive till the end, are only lunar-size at 1 Myr and located outside of 1.2 AU (Fig. 4). There are other larger bodies located inside of 1.2 AU at 1 Myr, but they will not survive Type I migration.
At 1 Myr, we introduce Jupiter and Saturn with their present orbital elements. Their formation causes a sudden increase of orbital eccentricities of nearby particles, particularly those in mean and secular resonances, while such a jump does not appear for inner planets (see the lower panel of Fig. 3). The locations of the ν 5 and ν 6 resonances at t = 1 Myr are 3.10 and 3.34 AU, respectively. These resonances, in particular ν 5 , effectively enhance orbital eccentricities of planetesimals and protoplanets, and these bodies migrate inward as the gas drag and the tidal interaction damp their eccentricities. At 3-5 Myr, the mean eccentricities of big and small bodies are comparable (Fig. 3) as many protoplanets are also involved in the ν 5 resonance. Near the resonance, the surface density of solid material is highly enhanced and relative velocities of solid bodies are small due to orbital alignments (Nagasawa et al., 2005). Consequently, protoplanets grow very quickly (Fig. 4). Collisions between protoplanets also commonly occur at this stage.
When the gas density becomes so low that inward migration of planetesimals and protoplanets due to the gas drag or the tidal interaction is slower than the migration rate of the secular resonance, these bodies are no longer swept by the resonance (Nagasawa et al., 2005). After the passage of the ν 5 resonance, orbital eccentricities of planets are reduced first by tidal interaction with a remnant gas disk. The damping due to gas is effective as long as the damping time is shorter than the gas decay time (τ damp < τ decay ). For an Earth-size planet at 1 AU, this condition is satisfied until ∼ 7-8 τ decay . Indeed, we see a rapid decrease of the mean eccentricity of big bodies at 5-7 Myr in Fig. 3. Later, the damping due to dynamical friction of planetesimals further reduces the eccentricities of planets. Its effect is determined by the total mass and the mass distribution of planetesimals (Morishima et al., 2008). Small bodies contain ∼ 20% of the total mass at 10 Myr and this amount is likely to be sufficient to achieve energy equipartition between big and small bodies as long as the planet-planet interaction is negligible. The typical size of remnant planetesimals is the lunar size or slightly smaller. There are only a few remnant planetesimals with the initial size because the smallest planetesimals had small orbital eccentricities due to strong gas drag and efficiently merged with larger bodies. The equilibrium eccentricity of Earth-size bodies due to dynamical friction of lunar-size planetesimals is estimated to be ∼ 0.03 (Morishima et al., 2008), and we obtain a quite consistent value after most of the remnant planetesimals merge with planets.
Through the above processes, three planets form near 1 AU at the end of this simulation. The largest planet has 1.4 M ⊕ while both of other two have 0.3 M ⊕ . The locations of these three planets resemble those of the current terrestrial planets (except we do not have a Mercury analog), although the masses are somewhat different. Orbital eccentricities are as small as those for the current planets. The most serious issue for this simulation is that the last of the potentially Moon-forming impacts occurs too early (at 13.5 Myr; see Sec. 3.2 for the definition) when compared with that suggested from isotope records (Touboul et al., 2007).  (Figs. 2-4). If the orbits of the gas giant planets are nearly circular, the effects of the secular resonances are much weaker than the case of the EJS orbits. Consequently, a large fraction of particles stay in the asteroidal region for a long period of time. Nevertheless, the orbital eccentricities of the gas giant planets can not be completely zero when there are two (or more) gas giant planets. Therefore, resonances gradually sculpture orbits of particles in the asteroidal region. Some enhancements at the 3:1 mean motion resonance (2.6 AU) and the ν 6 secular resonance (3.2 AU) are seen at 3 Myr (Fig. 5). These enhancements gradually propagate to non-resonant locations via mutual gravitational interactions (Raymond et al., 2006). The mass transfer from the asteroidal region to the terrestrial region slowly occurs as excited asteroidal bodies encounter and collide with inner planets.

Case with CJS orbits
At 3 Myr, there are about 15 Mars-size protoplanets in the terrestrial region. The fraction of small particles is very small in the terrestrial region whereas small particles are still dominant in mass in the asteroidal region. This situation is similar to some of the initial conditions adopted in previous N-body simulations (e.g., Raymond et al., 2006). In our simulation, however, the mean orbital separation between protoplanets is ∼ 20 Hill radius at 3 Myr while previous late stage simulations usually assume the orbital separation to be ∼ 10 Hill radius, based on the oligarchic growth model of Kokubo and Ida (1998). The larger separation in our simulation is due to Type I migration, which makes inner protoplanets start to migrate inward earlier than outer planets (McNeil et al., 2005). With the large orbital separation and the orbital stabilization due to the nebular gas, mutual collisions between protoplanets are rare until ∼ 10 Myr.
After 10 Myr, mutual collisions between protoplanets start to occur, and the orbital separation increases as the number of protoplanets decreases. Strictly speaking, the orbital separation normalized by the mutual Hill radius, and averaged over pairs of planets, increases nearly logarithmically with time (from ∼ 20 at 3 Myr to ∼ 60 at 800 Myr), and this is probably because the stable time scale of a planetary system exponentially increases with the normalized separation (Chambers et al., 1996;Ito and Tanikawa, 1999;Yoshinaga et al., 1999). Through giant impacts, an Earth-size planet forms slightly inside of 2 AU. This planet highly enhances orbital eccentricities of asteroidal bodies and their number gradually decreases by either hitting this planet or escaping from the system after encounters with Jupiter.
There are still a large number of remnant planetesimals at 100 Myr, but not at 300 Myr (Figs. 5 and 6). Therefore, damping due to planetesimals is expected until giant impacts complete at around 100 Myr. In this simulation, the last potentially Moon-forming impact occurs at 139.5 Myr on the outermost planet (Fig. 7). Nevertheless, the mean orbital eccentricity in the end of the simulation (∼ 800 Myr) is not small, as the two middle planets are still strongly interacting (Fig. 7). These two planets are stirred by the innermost and outermost planets which have larger masses than the middle planets. This simulation succeeded in reproducing the very late giant impact, but failed to reproduce the nearly circular orbits and the radial mass concentration of the current system.

Comparison with constraints
In this section, we examine the properties of planetary systems formed in all runs and compare them with those for our own system. We consider the following three properties of systems (see also Raymond et al., 2009). (Chambers, 2001) S c = max m j m j log 10 (a/a j ) , where m j and a j are the mass and the semimajor axis of each planet, and the subscript j is given in the order of a j . The function in the parenthesis is calculated for a throughout the terrestrial planet region, with S c being the maximum value. This quantity increases with radial mass concentration. We use two additional parameters as constraints. The first one is the minimum orbital separation between neighboring planets b min normalized by their mutual Hill radius r H, The second being the mass weighted mean semimajor axis of planets where M T,final = m j is the final total mass of planets. For the terrestrial planets, S c = 89.9, b min = 26.3 (the normalized orbital separation between the Earth and Venus), and a m =0.90 AU.
2. The deviation from circular and coplanar orbits. The orbital excitation of terrestrial planets is characterized by the angular momentum deficit, S d (Laskar, 1997): where e j and i j are the orbital eccentricity and inclination of each planet with respect to the fiducial plane of an initial planetesimal disk. We also measure the mass weighted mean eccentricity: We use the orbital elements averaged over a few Myr. For the terrestrial planets, S d = 0.0018 and e m = 0.038 (Quinn et al., 1991).
3. The timing of the Moon-forming impact. Based on simulations for giant impacts and accretion of moons (Cameron and Benz, 1991;Ida et al., 1997;Canup et al., 2001;Canup, 2004;2008), we roughly regard an impact as potentially Moon-forming, if a) the total mass of the impactor and the target is > 0.5 M ⊕ , b) the impactor's mass (< the target's mass) is > 0.05 M ⊕ , and c) the impact angular momentum is > 0.5 L EM , where L EM is the current total angular momentum of the Earth-Moon system (i.e., the Earth's spin angular momentum and the Moon's orbital angular momentum). We define the time of the last potentially Moon-forming impact in a run to be t imp . The time of the actual Moon-forming impact is estimated as t imp = 50-150 Myr (Touboul et al., 2007;Allégre et al., 2008). 2 The small content of iron in the Moon is explained when the impact velocity is as slow as the mutual escape velocity (Canup, 2004(Canup, , 2008; this constraint will be also discussed ( Figs. 13 and 16). Raymond et al. (2009) considered two additional constraints. The first one is the absence of planetary bodies (> 0.05 M ⊕ ) in the asteroidal region (> 2 AU). They showed that bodies which are too massive destroy the observed radial variation of the asteroidal taxonomic types (Gradie and Tedesco, 1982). If massive bodies exist in the asteroidal region, this system usually has too large a m or too small S c , because we do not exclude a massive body in the asteroidal region to be defined as a planet. Therefore, this constraint is automatically included in our first constraint. Their second additional constraint is the Earth's water content. They showed that the mass supply from asteroidal bodies to terrestrial planets decreases with increasing eccentricity of Jupiter. As a result, the water content of an Earth-analog is insufficient in their EJS and EEJS simulations. On the other hand, in our EJS and EEJS simulations with the nebular gas, a large fraction of mass of an Earth-analog comes from the asteroidal region, while protoplanets which formed from planetesimals originally in the terrestrial region spiral into the Sun by Type I migration (Sec. 3.1.1). Probably, the problem here is rather that we need to consider how to remove excess water from our Earth-analogs.

Comparison between EJS and CJS simulations
The properties of the planetary systems from the EJS and CJS simulations are summarized in Table 1. Also, Figures 8 and 9 show the snapshots on the a-e plane at the end of all EJS and [Table 1] CJS simulations, respectively. For both EJS and CJS simulations, the total mass and the typical [Figs. 8 and 9] size of planets decrease with increasing gas decay time τ decay because of Type I migration. When the surface density is large, protoplanets rapidly grow and then spiral into the Sun. Therefore, the final total mass M T,final depends on the initial surface density only weakly for a certain τ decay . This trend is particularly evident when τ decay ≥ 2 Myr, because the mass depletion due to Type I migration propagates to the outer asteroidal region for large τ decay . The number of planets N p is usually 3 or 4. In the EJS simulations with τ decay = 5 Myr, however, N p is smaller; in these simulations, the ν 5 resonance effectively sweeps most of the remnant mass to its current location (0.61AU), and a planet with dominant mass (0.3 − 0.5 M ⊕ ) forms there.
We plot the final planets on the plane of semimajor axis versus mass (Fig. 10). The crosses [Fig. 10] on the figure represent the current terrestrial planets. The EJS simulations clearly have better radial mass concentration than the CJS simulations. In the CJS simulations, the effect of the ν 5 resonance is negligible so that we inevitably have a massive planet around 2 AU, quite in contrast to the EJS simulations. In both EJS and CJS simulations, sizes of outer planets are smaller when the initial surface density in the asteroidal region is small, as is the case for small M T or large p. In addition, for the EJS simulations, the radial mass concentration is indirectly affected by the initial spatial distribution of planetesimals via their interactions with Jupiter. Asteroidal particles reduce e J either by angular momentum exchange at the ν 5 resonance or by energy exchange in close encounters with Jupiter; in both cases, the decrease in e J is given as ∆e J ∼ m/(m J e J ), where m J is Jupiter's mass. With decreasing e J , the effect of the ν 5 resonance decreases (Nagasawa et al., 2005). If the total mass in the asteroidal region is small due to small M T , large p, or large τ decay , the decrease in e J is small (see final e J in Table 1). Such systems can obtain high radial mass concentration since Jupiter's perturbation remains effective. This accounts for the large difference in the radial mass concentration between cases with M T = 5 and 10 M ⊕ in the EJS simulations (Fig. 10). The ν 6 resonance is also important for the radial mass concentration, and it also works effectively when e S is large (we will discuss its effect at Fig. 17 more in detail). Since Saturn's eccentricity e S is primarily determined by the eccentricity forced by Jupiter, e S is usually comparable to or slightly larger than e J . Figure 11 shows the plot of the final planets on the plane of mass versus orbital eccentricity. The small orbital eccentricities are well reproduced in the EJS simulations with M T = 5 [ Fig. 11] M ⊕ but not in those with M T = 10 M ⊕ . This difference is closely related to the radial mass concentration. In the case with M T = 5M ⊕ , the ν 5 resonance effectively sweeps most of the bodies and giant impacts between protoplanets complete early while there remains a sufficient amount of gas and planetesimals. On the other hand, the systems starting with M T = 10 M ⊕ radially expand even after the secular resonance passage and Jupiter continuously perturbs the outer planets. These planets eventually trigger orbital instability and cause giant impacts after both gas and planetesimals significantly deplete. Giant impacts tend to occur very late in the CJS simulations. However, small asteroidal bodies remain for a long time here as well and tend to damp eccentricities of the planets. This is why planets in CJS simulations have moderately small eccentricities both for M T = 5 and 10 M ⊕ . For the CJS simulations, we found a simple correlation between the mass weighted mean orbital eccentricity, e m , and the mass weighted mean semimajor axis, a m (Fig. 12). The longer the gas decay time, the larger a m , because mass [ Fig. 12] depletion due to Type I migration propagates outward with time. Since Jupiter perturbs outer planets more strongly than inner planets, we have this simple correlation. The CJS simulations always produce too large a value for a m , while the EJS simulations well reproduce the current value of a m except for runs with τ decay = 5 Myr where a m becomes too small. Figure 13 shows angular momenta and impact velocities of all potentially Moon-forming impacts as functions of time. The most important result is that giant impacts occur most com- [Fig. 13] monly at ∼ 10 Myr in the EJS simulations whereas at ∼ 100 Myr in the CJS simulations. In the EJS simulations, some early giant impacts occur between protoplanets trapped in the ν 5 resonance. It is more common, however, that giant impacts occur slightly later when orbital eccentricities of planets have sufficiently damped after passage of the resonance so that mutual gravitational attraction is enhanced. In the CJS simulations, giant impacts usually occur after gas completely depletes and planetesimals in the asteroidal region also significantly (but not completely) deplete. The time of the last Moon-forming impact in each simulation is shown in Table 1. In fact, last giant impacts in some of the EJS simulations occur very late. These systems, however, inevitably have either large eccentricities or low radial mass concentrations.
When we see a single simulation in Fig. 13, the impact velocity, v imp , is as small as the mutual escape velocity v esc at early time, but v imp /v esc increases with time so that its value for the last impact is notably larger than unity. Some simulations, however, have final impacts with v imp /v esc ≃ 1, which is favorable for the small iron content of the Moon (Canup, 2004(Canup, , 2008. When the orbital separation of two planets, normalized by the mutual Hill radius, slowly decreases due to their growth or weak perturbation from other planets, an impact between them with v imp /v esc ≃ 1 can occur after the Jacobi integral in Hill's equation slightly exceeds the potential at the Hill's sphere. Orbits of other protoplanets are affected only weakly during the time that two colliding protoplanets interact. On the other hand, if three or larger number of planets interact strongly, their relative velocity increases to the escape velocity. In this case, the typical impact velocity is v imp /v esc ≃ √ 2. Further higher v imp /v esc are possible between small protoplanets stirred by larger planets or if Jupiter's perturbation is strong. The former type of slow giant impacts usually occur in early stages while the latter type of high-velocity giant impacts occur after depletion of both gas and planetesimals (Fig. 13). Systems with S d and e m as small as those for the current system do not usually experience the latter type of impacts. Therefore, not only the small orbital eccentricities of the Earth and Venus but also the small iron content of the Moon are likely to support the presence of some damping forces at the time of the Moon-forming impact.

Cases with larger eccentricities of giant planets: EEJS and CJSECC simulations
As shown in Table 1, e J (and e S ) decreases via interactions with particles in the asteroidal region. Remnant planetesimals orbiting outside of Jupiter's orbit may also reduce e J , although we do not include these bodies in our simulations. Therefore, the initial e J and e S are likely to have been larger than those adopted in the EJS simulations, if their initial semimajor axes were (nearly) the same as the present values. On the other hand, if their initial semimajor axes were those given in the Nice model, e J and e S might have been nearly zero throughout the terrestrial planet formation, but larger initial e J and e S are not excluded. The summary of simulations with larger initial e J and e S , namely the EEJS and CJSECC simulations, is given in Table 2. The [ Table 2] snapshots of planets in the end of simulations in the a-e plane are also shown in Figs. 14 (the EEJS runs) and 15 (the CJSECC runs).
[Figs. 14 and 15] In the EEJS simulations, the ν 5 resonance is so strong that a m usually becomes too small. The number of planets is typically one or two, because protoplanets experience large orbital excursions. None of the simulations have S c as large as the present value, because the semimajor axis of the innermost planet is usually too small. The innermost planet is pushed toward the Sun by outer planets which are swept by the ν 5 resonance. The potentially Moon-forming impacts usually occur too early and impact velocities are too high (Fig. 16). Clearly, the initial e J of our [ Fig. 16] EEJS simulations is too high. If e J is somewhat smaller but still larger than the present value, the ν 5 resonance may work much better (Thommes et al., 2008).
The initial orbital eccentricities of Jupiter and Saturn in the CJSECC simulations are similar to those in the EJS simulations. The orbital separation between Jupiter and Saturn in the CJSECC simulations is narrower than that in the EJS simulations. This makes the ν 5 and ν 6 resonances in the CJSECC simulations more distant from the Sun than those in the EJS simulations. Because of this, the averaged a m and t imp in the CJSECC simulations become larger than those in the EJS simulations (Tables 1 and 2). Figure 17 shows a comparison between a CJSECC and an EJS simulation with the same parameters except for the orbits of Jupiter and Saturn. The ν 5 resonance sweeps planetesimals effectively at 3 Myr but no longer works at 5 [Fig. 17] Myr in both cases. In the EJS simulation, however, the ν 6 resonance (∼ 2 AU) still sweeps planetesimals toward the terrestrial region. On the other hand, the location of the ν 6 resonance in the CJSECC simulation is too far (∼ 3 AU) to contribute to radial mass concentration. In the CJSECC simulation, accretion of outer bodies, which survived the sweeping by the ν 5 resonance, proceeds slowly and this makes the final giant impact late. Surprisingly, 7 out of 10 CJSECC simulations in which potentially Moon-forming impacts occur satisfy the time of the last giant impact suggested from isotope records (Table 2). However, v imp /v esc 's of most of these late impacts are too high (Fig. 16) and also S c 's of these systems are too small ( Table 2).
As we have seen, simulation outcomes strongly depend on Jupiter's eccentricity, e J . Through interactions with small bodies, e J decreases. Clearly, the effect of Jupiter's perturbation is more important if e J in the late stage is larger regardless of initial values of e J . In Fig. 18, we plot the angular momentum deficit S d , the radial mass concentration parameter S c /a 2 m , and the time of the last one of potentially Moon-forming impacts t imp versus the value of e J at the end of a simulation, e J,final (averaged over a few Myr). We find that S c /a 2 m more appropriately represents [ Fig. 18] the radial mass confinement than does S c because systems consisting of small planets with large a tend to obtain large S c while b min is large (particularly in the CJS simulations). Simple linear χ 2 fits suggest that S c /a 2 m and t imp increases/decreases with e J,final , although the scatter of the simulation data are very large. The problem is that S c /a 2 m can be as large as the present value when e J,final is large whereas t imp can be as large as that suggested from isotope records only when e J,final is small enough. For S d , we fit the data with a quartic function as two competing effects are expected: large e J means stronger perturbation while quick accretion due to large e J means stronger damping forces due to remnant gas and planetesimals. We find that S d is nearly flat for small e J,final and increases with increasing e J,final when e J,final is larger than the present value. In an ensemble average sense, S d is usually much larger than the current value, although some simulations have S d consistent or even smaller than the present day value. The dependence of all three properties on e J,final is stronger for the CJS+CJSECC simulations than for the EJS+EEJS simulations. It is unclear to us, however, whether this is really physically meaningful because fits for the CJS+CJSECC simulations are not very reliable with only a few points at large e J .

Discussion
None of our simulations succeeded in satisfying all of the three constraints listed in Sec. 3.2. The most serious issue we found is a trade-off between the radial mass concentration parameter, S c /a 2 m , and the time of the potentially Moon-forming impact, t imp (the upper panel of Fig. 19). The large S c /a 2 m in the current system is due to the small masses of Mars and Mercury and the [ Fig. 19] small orbital separation between the Earth and Venus. It is hard to satisfy both of these conditions in a system with large t imp . Since the evolution time scale is longer in the outer regions, giant impacts can occur late when a large amount of mass stays in the outer region for a long time. Such a system usually has outer planets more massive than Mars. Also, if a planetary system slowly evolves with weak perturbations from Jupiter (e.g., the CJS simulations), the orbital separation between neighboring planets normalized by the mutual Hill radius inevitably becomes larger than that between the Earth and Venus (26.3) owing to the mutual orbital repulsion. If e J is as large as the present value or larger, S c /a 2 m can be as large as today's value because Jupiter tends to radially compress a system. However, accretion completes too rapidly in such a case.
There is also a trade-off between the angular momentum deficit S d and t imp (the lower panel of Fig. 19). Namely, in order to obtain small S d , a final large impact needs to occur while some damping forces work effectively. Otherwise, S d ends up with much larger than the current value. This trend is particularly clear for the EJS and CJS simulations (filled and open circles). Small S d even without any damping forces seems to be possible but such a case is statistically rare.
Our failure indicates that there may be some missing physics in our simulations or flaws with the assumptions that we have adopted. In the following, we suggest three possible scenarios/solutions for future studies.
In the case with an eccentric Jupiter, a possible solution is that τ decay is much longer and Type I migration is much slower than we assumed. This situation actually corresponds to that adopted in the original dynamical shake up model of Nagasawa et al. (2005) and Thommes et al. (2008). In their simulations (τ decay = 3 and 5 Myr), the giant impacts usually occur at t ≃ 4 − 5 τ decay when the ν 5 resonance reaches near 1 AU. Therefore, τ decay > 10 Myr is required for the late Moon-forming impact indicated from isotope records (t imp > 50 Myr; Touboul et al. 2007). If τ decay = 10 Myr, the gas density is still ∼ 1% of the minimum mass solar nebular at 50 Myr, and observations suggest that such a gas disk is exceptionally long-living (Pascucci et al. 2006). Slightly smaller t imp (≃ 30 Myr) is not excluded from the terrestrial W/Hf records only (Jacobsen et al., 2009), but this still requires a long-surviving disk (τ decay ∼ 6 Myr). Recent theories indicate significant changes in Type I migration rate due to the non-isothermal effect (Baruteau and Masset, 2008;Paardekooper and Mellema, 2008) and the non-linear effect of the co-rotation resonance (Masset et al., 2006;Paardekooper and Paparoizou, 2009). These mechanisms cause slow down of inward migration or even outward migration if the radial entropy gradient is negative for the former mechanism or if the radial surface density gradient is nearly flat or positive for the latter one. Therefore, if a local density maximum exists, planetary migration slows down or protoplanets even tend to migrate toward this point (Kretke et al., 2009). A local density maximum can exist at the boundary between the optically thin inner region and the optically thick outer region (i.e., the inner boundary of the so-called dead-zone; Kretke and Lin, 2007;Brauer et al., 2008). Protoplanets might have survived in the protosolar nebular for a long time with such mechanisms.
An alternative possibility for the case with an eccentric Jupiter is that the nebular gas inside of Jupiter might have dissipated very rapidly after the formation of Jupiter. This situation corresponds to previous N-body simulations for the late stages without gas (Chambers, 2001;O'Brien et al., 2006;Raymond et al., 2009). Without gas, t imp tends to be larger than those in our simulations, because the effect of ν 5 is less signifiant. Non-uniform gas dissipation models, such as inside-out and gap-opening models discussed in Nagasawa et al. (2000), probably have similar outcomes, because ν 5 passes a certain location after the gas dissipates there. Raymond et al. (2009) showed that some of their EEJS simulations nearly satisfy all the three constraints listed in Sec. 3.2. In fact, these simulations satisfy the small mass of Mars, but not the small orbital separation between the Earth and Venus. This is why the values of S c in all of their simulations are smaller than the current value. This problem might be resolved, however, if they adopt even higher resolution with smaller planetesimals (e.g., Morishima et al., 2008) and inelastic rebounds for high-speed or grazing collisions (Agnor and Asphau, 2004;Asphau et al., 2006). Because these effects are likely to enhance the damping effects and reduce the orbital excursions of planets, the orbital separation between neighboring planets will probably be smaller. Very quick or non-uniform dissipation of gas seems to be favorable, at least, in the asteroidal region in the case with large e J . As we have shown, the ν 5 resonance is so strong in the asteroidal region, except the CJS simulations, that the radial compositional structure in the asteroidal region is destroyed even for τ decay = 1 Myr. The gas dissipation of the order of 0.1 Myr or shorter seems to be necessary for the uniform dissipation model (Nagasawa et al., 2000). Jupiter might have opened a wide gap around its orbit, extending to the entire asteroidal region immediately after its formation (however, τ decay = 1 Myr or larger may be still possible if e J is small, when ν 5 passes the asteroid region, and increases later).
In the case of a circularly orbiting Jupiter, the initial planetesimals needed to be radially highly concentrated in order to satisfy a large S c . In other words, the depletion of mass in the asteroidal region needed to occur before the formation of planetesimals. Such a concentration of planetesimals may be possible if planetesimal formation occurs in radially limited regions after the increase of the surface density due to dust migration (Youdin and Shu, 2002;Youdin and Chiang, 2004;Kretke and Lin, 2007;Brauer et al., 2008). Morishima et al. (2008) performed N-body simulations starting from highly concentrated disks in the absence of gas, and found that some of the resulting planetary systems were similar to the current solar system. The problem of their simulations is that accretion usually completes very quickly ( < 10 Myr), as well as our EJS simulations. Nevertheless, one simulation (run 1 of Morishima et al., 2008) has a very late giant impact. This system unfortunately has slightly too large an orbital separation between the two largest planets. Hansen (2009) also made simulations similar to Morishima et al. (2008) and found the mean time of last giant impacts to be 45 Myr. However, the lower limit to the impactor mass (0.02 M ⊕ ) used for his definition of giant impacts is probably too small for formation the Moon. Further explorations for planetary accretion starting from radially concentrated disks are necessary.

Conclusions
We have performed N-body simulations which follow the growth and accretion histories of the terrestrial planets, including the effects of the giant planets and the interaction with gas in laminar disks. We vary initial total mass and spatial distribution of planetesimals, the time scale of dissipation of a gas disk, and the orbits of Jupiter and Saturn. It is a remarkable success of numerical simulations, that planetary systems with global properties similar to the terrestrial planets can be obtained starting with a proto-planetary disk. We are now at the stage where fine details of the models can be compared with solar system constraints. In order to retain a final total mass as massive as the present total mass (∼ 2 Earth mass) in the presence of standard Type I migration, the time scale of dissipation of gas needs to be 1-2 Myr.
Using the results from a large number of simulations we compared the evolutionary histories of the resulting planetary systems with the following three properties of the terrestrial planets: (1) the high radial mass concentration, (2) the small orbital eccentricities and (3) the very late Moon-forming impact (∼ 100 Myr). The orbital eccentricities of the planets can be as small as their present values if giant impacts finish whilst there is still remnant gas or planetesimals, as expected. The most serious issue we found was a trade off between the radial mass concentration and the late Moon-forming impact. If the orbital eccentricities of Jupiter are as large as the present value, most of bodies in the asteroidal region are swept up to the terrestrial region owing to inward migration of the secular resonance. Although this mechanism helps planetary systems to have high radial concentrations in the end, giant impacts between protoplanets most commonly occur too early (∼10 Myr). On the other hand, if the orbital eccentricity of Jupiter is close to zero as suggested in the Nice model, the effect of the secular resonance is negligible and a large amount of mass stays for a long period of time in the asteroidal region. Although giant impacts usually occur around 100 Myr, we inevitably have an Earth size planet at around 2 AU in this case. It is very difficult to obtain spatially concentrated terrestrial planets with very late giant impacts, as long as we include all the effects of gas and assume initial disks similar to the minimum mass solar nebular.  Table 1: Simulation parameters and outcomes for the EJS and CJS simulations. Table columns are the orbits of Jupiter and Saturn, the gas decay time scale τ decay , the power-law index for the radial gradient of the initial surface density of planetesimals p, the initial total mass of planetesimals M T , the number of planets N p , the final total mass of planets M T,final , the radial mass concentration statistics S c , the mass weighted mean semimajor axis a m , the minimum orbital separation between neighboring planets b min normalized by the mutual Hill radius, the mass weighted mean eccentricity e m , the angular momentum deficit S d , the time of the last potentially Moon-forming impacts t imp , and the final orbital eccentricity of Jupiter e J,final . Any orbital elements averaged over a few Myr are used. Each value is boldly marked when 3 ≤ N p ≤ 5, 1.0 ≤ M T,final /M ⊕ ≤ 3.0, S c > 45.0, 0.75 AU < a m < 1.05 AU, b min < 35, e m < 0.076, S d < 0.0036, or 50 Myr < t imp < 150 Myr. SS represents the solar system. -Time scales for radial migration τ mig and damping of eccentricity τ damp due to gas drag and tidal interaction (τ mig = a[(da/dt) drag + (da/dt) tid ] −1 and τ damp = e[(de/dt) drag + (de/dt) tid ] −1 ) at 1 AU with Σ gas = 2000 g cm −2 and i = 0. The migration and damping rates due to gas drag, (da/dt) drag and (de/dt) drag , are derived after Adachi et al. (1976) with some corrections given in Inaba et al. (2001). The labels to the lines represent e. For e = 0.05, we also plot the time scales with the enhancement of the gas drag force for m 0 = 0.0025 M ⊕ , m 1 = 0.01 M ⊕ , and m GI = 10 19 g (see Sec. 2.2.1). The masses of solar system bodies are also plotted in the upper panel. The size of each body is proportional to the radius, except Jupiter's size is modified to the Earth size on this plot. The color represents spin rate ω: ω > ω crit , ω crit ≥ ω > 0.5 ω crit , 0.5 ω crit ≥ ω > 0, and ω = 0 for red, blue, green, and black bodies, respectively, where ω crit is the spin rate with which the centrifugal force balances with the gravity at the surface. The longdashed and short-dashed lines are locations of the ν 5 and ν 6 resonances, respectively. The black dashed curve is given by a(1 + e) = a J ; particles in the right side of this curve experience orbital crossing with Jupiter. Note the difference in scales of e between the top panels and others. and the mass weighted mean eccentricities of small and big particles (the lower panel) for an EJS simulation. We define mass of a small (big) particle to be below (above) 2 × 10 26 g.           are results from the EJS, EEJS, CJS, and CJSECC simulations, respectively. The color represents the gas decay time: τ decay = 1, 2, and 3 Myr for black, red, and blue symbols. A fit to a quartic function is applied to log S d , while linear fits are applied to S c /a 2 m and log t imp for sets of the EJS+EEJS and CJS+CJSECC simulations, respectively. Horizontal and vertical dotted lines indicate the current values (for t imp , a possible range indicated from isotope records is between two dotted lines) Fig. 19.-Radial mass concentration parameter S c /a 2 m (upper panel) and angular momentum deficit S d (lower panel) vs. time of the Moon-forming impact t imp . The data are the same as those in Fig. 18.