Highlight
Free Access
Issue
A&A
Volume 579, July 2015
Article Number A65
Number of page(s) 8
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201525636
Published online 29 June 2015

© ESO, 2015

1. Introduction

Planets with masses larger than approximately 0.1 M significantly excite density waves in a protoplanetary disk, which in turn exerts torque on the planets, and migrate toward the central star (Type I migration; e.g., Goldreich & Tremaine 1980; Ward 1986). It is known that the inward migration timescale for an Earth-mass planet at 1 AU is ~105 yr in a locally isothermal disk (e.g., Tanaka et al. 2002), which is shorter than the disk lifetime and thus can be a problem for terrestrial planet formation in the solar system. Recent studies have shown that in an adiabatic disk, fully unsaturated positive corotation torques can compensate for negative Lindblad torques, which can change the direction and rate of migration (e.g., Paardekooper et al. 2010). If this effect is significant, a convergence zone where migration is convergent can be created in a disk and gas giant cores can form at a few tens of AU (Horn et al. 2012).

Paardekooper et al. (2011) revised the torque formula by taking into account the effects of viscous and thermal diffusion on the corotation torque, which indicates that only planets with a limited range of masses can experience the non-linear corotation torque due to saturation effects. In fact, planets with masses smaller than a few Earth masses are not affected by the non-linear corotation torque (the horseshoe torque) and they migrate inward under the influence of the Lindblad torque (e.g., Kretke & Lin 2012; Hellary & Nelson 2012; Cossou et al. 2014), which means that the solar system’s terrestrial planets still have the problem of inward migration. However, if the slope of the gas surface density becomes large enough, the positive linear corotation torque can reverse the migration.

During the evolution of protoplanetary disks, vertical winds as well as radial accretion play a significant role in the dispersal of the gas component. While photoevaporating winds by UV and X-rays from a central star have been widely discussed (e.g., Alexander et al. 2006; Ercolano et al. 2009), the importance of spontaneously-driven disk winds has also been pointed out (Suzuki & Inutsuka 2009; Suzuki et al. 2010); using magnetohydrodynamic (MHD) simulations they showed that MHD turbulence triggered by magnetorotational instability (MRI) in disks inevitably drives vertical outflows as well as the radial transport of angular momentum and consequent accretion. The detailed properties of such disk winds driven by turbulent Poynting flux have been studied (Bai & Stone 2013; Fromang et al. 2013; Lesur et al. 2013). An intriguing characteristic of the disk wind is that disk dispersal proceeds gradually in an inside-out manner because the dispersal timescale is approximately proportional to rotation time r3/2, where r is the radial distance from the star. As a result, the surface density gradient at r ≲ 1 AU becomes much shallower than the typical values of 0.5 to 1.5, and could even be positive, which will drastically affect the migration of protoplanets.

As we are interested in seeing how planet formation proceeds in such a disk, we investigate the formation of terrestrial planets using N-body simulations in a disk where the effect of a disk wind is considered. We calculate the disk evolution using a one-dimensional disk with viscous diffusion and a disk wind based on Suzuki et al. (2010). In this work, we focus on the late phase of planet formation when the effect of a disk wind is clearly visible, and calculate the accretion and orbital evolution of planets from planetary embryos with masses of 0.1 M for ~100 Myr. The N-body code includes the effect of Type I migration that depends on the slope of gas surface density and the mass of the planets using the prescription described by Paardekooper et al. (2011). We also include the influence of eccentricity on the corotation torque.

Our primary aim is to clarify the effect of a disk wind on terrestrial planet formation. In addition, we also discuss the origin of the terrestrial planets in the solar system. A number of authors (e.g., Kominami & Ida 2002; Nagasawa et al. 2005; Ogihara et al. 2007; Raymond et al. 2009; Morishima et al. 2010) have studied the origin of these planets using N-body simulations; however, no single simulation has yet convincingly reproduced the observed constraints simultaneously (e.g., eccentricities, spatial concentration around 1 AU, late moon-forming giant impact). We will propose a possible model for the origin of the inner solar system.

The plan of the paper is as follows. In Sect. 2 we describe our model and the initial conditions of the simulations; in Sect. 3 we present the results of N-body simulations; in Sect. 4 we discuss the origin of the terrestrial planets in the solar system; in Sect. 5 we give our conclusions.

2. Model description

2.1. Disk model

We calculate the time evolution of the gas component of a protoplanetary disk by taking into account the radial mass flows by the transport of angular momentum through MHD turbulence and the mass loss due to the disk wind. Although there are uncertainties in the mass flux of the disk wind, we follow the simple framework of Suzuki et al. (2010). The evolution of gas surface density Σg can be expressed in the form of a diffusion equation, (1)where Ω is the Keplerian frequency. The standard α viscosity prescription, ν = αcsh, is used, where cs indicates the sound velocity and h is the disk scale height. The second term on the right-hand side of Eq. (1) is the disk wind flux ρvz, which can be expressed as Cwρ0cs using the mid-plane density ρ0, cs, and a non-dimensional constant Cw (Suzuki et al. 2010). We note that the coefficient of the first term on the right-hand side of Eq. (1) is not exactly the same as that of Eq. (9) of Suzuki et al. (2010) owing to a different definition of ν, which gives a slightly different evolution of the disk surface density.

For the initial distribution of the gas disk, we assume (2)where fg is a dimensionless parameter and rcut = 50 AU is used as a cutoff radius. In our N-body code, we use the evolution of Σg to calculate the effect of the gas disk on planet evolution.

The temperature distribution is assumed to be that of an optically thin disk (Hayashi 1981), (3)We note that although Type I migration is also dependent on the temperature profile, we fix the temperature distribution to clarify the effect of the disk wind via surface density distribution. Then the disk scale height h is given by (4)where L and M are the luminosity and mass of the host star, respectively.

2.2. Type I migration and eccentricity damping

The orbits of planetary embryos with masses M1,M2,... and position vectors r1,r2,... relative to the host star are calculated by numerically integrating the equation of motion, (5)where Fdamp is a specific force for eccentricity and inclination damping and Fmig is a specific force for Type I migration (see Ogihara et al. 2014 for each formula). The timescale for damping of the eccentricity, te, is given by where Σg = 2400fg(r/ 1 AU)−3/2 g cm-2 is used in the second equality. Here the relative motion between gas and planets is assumed to be subsonic (evKcs). For planets with high eccentricities and inclinations, we include a correction factor according to Eqs. (11) and (12) of Cresswell & Nelson (2008).

The migration timescale, ta, is given by where β is a coefficient that determines the direction and speed of Type I migration. According to Paardekooper et al. (2010), the Type I migration torque depends on the Lindblad torque, the barotropic part of the horseshoe drag (or the non-linear corotation torque), the entropy-related part of the horseshoe drag, the barotropic part of the linear corotation torque, and the entropy-related part of the linear corotation torque. Paardekooper et al. (2011) derived the total Type I migration torque including both saturation and the cutoff at high viscosity. Thus we write the migration coefficient in the form (10)where βL,βc,baro, and βc,ent are related to the Lindblad torque, the barotropic part of the corotation torque, and the entropy-related part of the corotation torque, respectively. Each formula is given as where p and q denote the local surface density gradient (p(r) = −dlnΣg/ dlnr) and the local temperature gradient (q(r) = −dlnT/ dlnr), γ is the adiabatic index, and ξ is the local entropy gradient (ξ = q − (γ − 1)p). In this work, γ = 1.4 and q = 0.5. The function F(P) is a decreasing function with the value of [0,1] that is related to saturation, and the functions G(P) and K(P) are increasing functions with the value of [0,1] related to cutoff. Paardekooper et al. (2011) introduced viscous and thermal parameters Pν and Pχ, which are expressed by the ratio between the viscous/thermal diffusion timescales τν/τχ and the horseshoe libration timescale τlib, where the dimensionless half-width of the horseshoe region is (16)where ws is the half-width of the horseshoe region. For prescriptions of F(P), G(P), and K(P), we refer the readers to Eqs. (23), (30), and (31) in Paardekooper et al. (2011).

For the thermal diffusivity, we assume χ = ν in this work for simplicity. The thermal diffusivity determines the saturation of the entropy-related part of the corotation torque. The thermal diffusivity is considered to be large when the disk surface density decreases and the disk is optically thin, and thus radiative cooling is efficient. In this case, the entropy-related corotation torque diminishes, and one might speculate that the orbital evolution of planets would change. However, we do not expect that this can change our results in this paper. In the terrestrial planet forming region, the disk opacity is about 1−10 cm2 g-1 (e.g., Bell & Lin 1994). This means that the disk is optically thin only when the disk is significantly dispersed g ~ 0.1−1 g cm-2) and Type I migration is no longer effective. Therefore, our simplified model (χ = ν) may be justified. In addition, we observe that our qualitative results are not affected in test simulations in which the entropy-related torque is switched off from the beginning of the simulations.

Recent studies (e.g., Bitsch & Kley 2010) suggest that the corotation torque decreases as the planet eccentricity increases; therefore, we also consider this effect using the formulae (Fendyke & Nelson 2014), where ef = 0.5h/r + 0.01. We note that the Lindblad torque can also be reduced when planets have high eccentricities and inclinations (e.g., Papaloizou & Larwood 2000). In some runs of our simulations, we also add a reduction factor to Eq. (11), which is given by Eq. (13) of Cresswell & Nelson (2008). However, we found that even if this factor is ignored, the results do not change very much.

There are several works that take into account the effect of magnetic field on Type I migration (Terquem 2003; Fromang et al. 2005; Muto et al. 2008; Uribe et al. 2015; Bans et al. 2015). Among these works, Bans et al. (2015) consider the vertical transport of angular momentum by disk winds. In our paper, we assume the net vertical magnetic field is very weak and the disk wind is driven by the MRI triggered turbulent magnetic field, which is in contrast to the disk wind by relatively strong net magnetic field assumed in Bans et al. (2015). Therefore, the Type I migration formula based on pure hydrodynamics (which corresponds to the limit of weak magnetic field) by Paardekooper et al. (2011) is probably still reasonable in our setting.

2.3. Set of parameters and orbital integration

In our series of simulations, we mainly vary the following parameters: turbulent viscosity α, disk wind efficiency Cw. Although the efficiency of the disk wind should be related to the turbulent viscosity, we vary α and Cw independently to explore parameter space. The list of parameters for each model is summarized in Table 1.

Table 1

List of parameters of each model.

In model 1, the effect of the disk wind is not considered to clarify its effect on planet formation. The strength of the net vertical magnetic field is an important control parameter that determines α and Cw (Suzuki et al. 2010; Okuzumi & Hirose 2011). If the net vertical magnetic field is sufficiently weak, α and Cw stay more or less constant during the disk evolution. Model 2 corresponds to this weak regime; we adopt α = 8 × 10-3 and Cw = 2 × 10-5 according to the MHD simulations of Suzuki et al. (2010). If the net vertical magnetic field is relatively strong, α and Cw increase with time as gas materials are dispersed (Suzuki et al. 2010). In model 3 we consider this strong regime, and we use In addition, in model 4 we use α = 2 × 10-5 and Cw = 5 × 10-7.

The parameter set for model 4 may require some explanation. The slope of gas surface density basically depends on α/Cw. In other words, if the ratio of α/Cw is the same, the distribution of gas surface density takes the same form while the evolution timescale depends on α and Cw. In model 4 we consider the case of a more efficient disk wind (α/Cw = 40) than that of model 2 (α/Cw = 400). The value of α is also changed because we wish to discuss the origin of the solar system’s terrestrial planets in this paper. As we show in Sects. 3 and 4, outward migration of 0.1 Earth-mass embryos is favorable for reproducing the observed properties of the solar system and α ~ 10-5 is required to avoid the saturation of the corotation torque for these bodies. Therefore, α = 2 × 10-5 and Cw = α/ 40 = 5 × 10-7 are adopted. While the evolution of the surface density of model 4 is slower than models 2 and 3 because both α and Cw are smaller, the trend of the radial dependence of the surface density of this model represents a case more or less between model 2 and model 3. The initial amount of gas and the gas dispersal time are also changed in model 4. The scaling factor for the initial gas density fg is increased by a factor of three. In the other models, fg is set to unity. In addition, the disk would rapidly disperse (~0.1 Myr) after typical disk lifetimes of a few Myr, and additional effects are required to achieve this (e.g., photoevaporating winds; e.g., Alexander et al. 2014). In model 4, in order to mimic this disk evolution, α and Cw are increased by a factor of 100 after t = 3 Myr.

For the initial distribution of embryos, we put 60 protoplanets with mass M = 0.1 M between a = 0.1−3 AU. Their orbital separation, Δ, measured in mutual Hill radii is set to 10. We perform three runs of simulations for each model with different initial locations of protoplanets.

For numerical integration, we use a fourth-order Hermite scheme (Makino & Aarseth 1992). When the physical radii of two bodies overlap, they are assumed to merge into one body, conserving total mass and momentum assuming perfect accretion. The physical radius of an embryo R is determined by R = (3M/ 4πρ)1/3, where we assume the internal density of ρ = 3 g cm-3.

3. Results

3.1. Disk evolution

We first show the disk evolution and its contribution to Type I migration. The solid lines in Fig. 1 show the evolution of the disk surface density for models 2–4 while the dotted lines in each panel represent the surface density for model 1. Figure 2 shows the migration timescale (ta) for model 2 and model 4, where the migration rate and direction are indicated by the color (see color bar). When the migration timescale is negative, the direction of migration is inward. As we saw in Sect. 2, Type I migration depends on the eccentricity of planets. We assume e = 0.01 in these plots.

thumbnail Fig. 1

Evolution of gas surface density profile for a) model 2; b) model 3; and c) model 4. Solid lines in each panel indicate the gas surface density at t = 0,0.1 Myr,1 Myr, and 3 Myr, respectively. Dotted lines show results for model 1. In panels a) and b), the solid and dotted lines overlap for t = 0 because of the initial conditions.

Open with DEXTER

thumbnail Fig. 2

Time evolution of migration rate and direction for bodies with e = 0.01 for a) model 2 and b) model 4. The color indicates the migration timescale. The contours show migration timescales of 0.1 Myr, 1 Myr, and 3 Myr.

Open with DEXTER

Figure 1a clearly shows that the disk evolution of model 2 is altered by the disk wind in the following two aspects. First, the slope of the disk becomes shallow in the inner region. Second, the disk is depleted faster than the case without the disk wind. These two factors significantly reduce the inward migration of embryos. From the migration map of model 2 in Fig. 2a, we find that the inward migration timescale for sub-Earth mass planets is over 1 Myr at t = 0.1 Myr, which means they no longer undergo migration. We also observe no migration at t = 1 Myr and 3 Myr.

For the disk evolution of model 3 that has larger Cw and α shown in Fig. 1b, an inner cavity is created inside a certain radius at which 99% of the initial gas is dispersed according to Eqs. (19) and (20). As time passes, the inner cavity grows gradually and the disk inner edge moves outward. The inside-out evacuation of the disk can significantly affect the orbital evolution of planets, which will be shown in Sect. 3.2.

We also examine the disk evolution for model 4, in which the radial profile of the surface density in the inner region is approximately between model 2 and model 3. In Figs. 1 and 2 for model 2, we see the disk is depleted in the early stage of planet formation; therefore, a larger amount of initial gas disk is assumed (fg = 3). The evolution of the surface density in Fig. 1c shows that the slope becomes positive inside 1 AU. In order to get a positive surface density gradient inside 1 AU, we find that α/Cw ≲ 100 is required. The disk lifetime is longer than model 2, which is more consistent with the observationally inferred disk lifetime. The migration map in Fig. 2b is more complicated. Because the effect of the disk wind is stronger in the inner region, the orbital region where planets can move outward is limited to r ≲ 1 AU. The saturation of the corotation torque depends on the mass and the disk viscosity. If α = 2 × 10-5 is assumed, the saturation effect is weakest for bodies with M ≃ 0.1 M. Therefore, only sub-Earth-mass planets can undergo outward migration. It is also worth noting that planetary orbits can evolve until a few Myr because the disk is almost completely depleted after 10 Myr.

Before proceeding to the next subsection, where the results of N-body simulations are shown, here we briefly discuss how the inner disk profile is determined by viscous diffusion and the disk wind. We note that the disk wind brings about the variety in the radial dependence of the surface density in the inner region ranging from a clear hole (model 3) to mild inside-out evacuation (model 2). We suppose that the inner hole seen in model 3 is an extreme case. In a realistic situation low-density gas should remain in the inner region and a low level of accretion would continue, rather than the formation of a clear inner hole, even at later times as seen in model 4, because the origin of the energy to drive the disk wind is the gravitational energy liberated by gas accretion (Suzuki et al. 2010); the infall of a certain amount of gas is required to launch vertical outflows. Suzuki et al. (2010) demonstrated the time evolution of a global disk by a simple model taking into account such energy-limited disk winds, which shows that the radial profile of the surface density is slightly modified at ≲ 1 AU.

3.2. Orbital evolution

We show the results of N-body simulations in disks that evolved via viscous diffusion and a disk wind. Figure 3 shows the orbital evolution of planets. The size of each filled circle represents the radii of bodies. The color of lines indicates the eccentricity of the planets. Each panel shows the results for model 1, model 2, model 3, and model 4. We performed three runs for each model, and typical runs are shown in the plots. Figure 4 shows the final orbital configurations in the semimajor axis-mass plane for all runs of model 2, model 3, and model 4 (the results for model 1 are not shown), where planets that formed through the same run are connected with the same line.

thumbnail Fig. 3

Orbital evolution of planetary embryos for a) model 1; b) model 2; c) model 3; and d) model 4. The filled circles connected with lines represent each body. The color of lines indicates the eccentricity (see color bar). The size of circles is proportional to the radius of the body.

Open with DEXTER

In Fig. 3a, the results for model 1 (in which the disk wind is not considered) is shown for comparison. Planetary embryos migrate inward and 12 bodies are lost by t = 20 Myr. The computational cost is huge for this run because the bodies do not collide with each other and the orbital period of the innermost planet is short; therefore, the computation was stopped at t = 20 Myr, which took over four months of CPU time. Orbit crossings are not observed before t = 20 Myr; however, the system is expected to be destabilized at ~100 Myr according to previous studies on the stability of multi-planet systems (Chambers et al. 1996).

In the results for model 2, we find that planetary embryos do not migrate inward any more as expected from Figs. 1 and 2. In this model, as the disk dissipates, the eccentricity of the embryos is excited and they exhibit orbit crossings, leading to giant impacts between the bodies. At t = 100 Myr, the number of planets has decreased to 13 and their eccentricities are between about 0.01 and 0.1. The final orbital configurations are almost the same among the three different runs (see Fig. 4a). The averaged mass of the largest body in each run is ≃0.9 M. The typical orbital separation between the planets is ≃10−40rH.

For model 3, the orbital evolution is significantly different from that without a disk wind. As seen in Fig. 1b, the disk’s inner edge moves outward. The migrating inner edge sweeps up embryos and almost all bodies migrate outward. As a result, no planets are left inside 0.4 AU and a large planet forms at 5 AU. The eccentricities of the planets are ~0.04−0.3 (some bodies are still exhibiting orbit crossings). In the other two runs, the results are almost the same (see Fig. 4b). The mass of the largest body is ≃3.8 M, while the second largest planet is located at around 1 AU with a mass of ≃0.5 M.

We observe an interesting orbital evolution in model 3 in which most protoplanets are swept up by the disk edge before t ≃ 0.3 Myr but afterward several bodies drift through the edge. We briefly discuss how bodies follow the migration of the disk edge or not. First, just as there are for torques on planets, there are two different positive torques, namely the Type I migration torque and the edge torque. As seen in Sect. 2, the former can become positive owing to the corotation torque, which depends on the slope of the gas surface density. The latter is a torque recently found by Ogihara et al. (2010). When a body with non-zero eccentricity straddles a sharp disk inner edge, the body gains a net positive torque from the disk depending on the sharpness of the edge and on the eccentricity. Whether planets can follow the movement of the disk edge depends on the migration speeds of the planets and the edge. If the outward migration rate of the planets due to the positive torque (which is the sum of the Type I migration torque and the edge torque) is smaller than the migration rate of the edge, the planets go through the disk edge and stay in place. Through a series of test calculations, we find that the positive torque exerted on the planets is usually large for this model and 0.1 Earth-mass planets can follow the disk edge. Therefore, there is another reason why some planets drift through the disk, which is that small embryos are scattered inside the disk edge by larger bodies. Larger bodies grow at the disk edge by sweeping up smaller embryos. At t ≃ 1 Myr, the planets at the edge already have masses of ≃3.5 M and scatter embryos inside the edge.

The results for model 4 are also interesting. Before t ≃ 0.02 Myr, all bodies undergo inward migration, but afterward ~0.1 Earth-mass embryos in close-in orbits move outward while bodies in distant regions migrate inward, leading to a convergence of planets at around 1 AU. For α = 2 × 10-5, the positive corotation torque can be saturated for bodies larger than ≃0.3 Earth masses; as a result, several planets that undergo collisional growth migrate inward and are lost to the central star. As the disk disperses after t = 3 Myr, all the bodies cease their migration. At t = 100 Myr, a system with six planets with a maximum mass of 0.7 M and e ≃ 0.01 had formed. The orbital separations between the planets are relatively small (down to ≲ 10rH), and giant impacts can occur after t = 100 Myr. Thus we continued simulation until t = 200 Myr for this run. The system underwent a giant impact at t = 145 Myr, which is discussed in Sect. 4. In the other two runs, although the mass distributions are slightly different from each other, the results are qualitatively the same (Fig. 4c).

4. Origin of the solar system’s terrestrial planets

Here we discuss our results in terms of the origin of terrestrial planets in the solar system. The important characteristics are the radial mass concentration and small eccentricities of planets, which we characterize using two statistical measures, Sc and Sd. The radial mass concentration statistic is (Chambers 2001) (21)where j = 1,2,...,N. The angular momentum deficit, which represents the deviation from circular and coplanar orbits, is (Laskar 1997) (22)where ij is the orbital inclination of the body j.

The radial mass concentration of the solar system’s terrestrial planets is Sc = 89.9, which indicates high radial mass concentration because the masses of Mercury and Mars are small and the orbital distance between Venus and the Earth is also relatively small. Previous studies have shown that it is difficult to reproduce such a high mass concentration (e.g., Raymond et al. 2009; Morishima et al. 2010); a way to account for this statistic in N-body simulations is to assume a mass distribution initially confined to a narrow annulus (Morishima et al. 2008; Hansen 2009), or in other words to assume an initially large Sc. In the grand tack model (e.g., Walsh et al. 2011; Jacobson & Morbidelli 2014), Jupiter sweeps up to ~1.5 AU and embryos are eventually confined to a relatively compact region. However, they assumed an inner cutoff at ≃0.7 AU, which means that Sc is already relatively large in the initial state.

thumbnail Fig. 4

Final orbital configurations for a) model 2; b) model 3; and c) model 4. Only for model 4a is the configuration at t = 200 Myr shown. Results of three different runs are plotted in each panel.

Open with DEXTER

The angular momentum deficit for the terrestrial planets in the solar system is Sd = 0.0018, which means that the eccentricities of the terrestrial planets in the solar system are relatively small (e ≃ 0.01−0.1). We note that Jupiter and Saturn increase Sd after their formation, which is ignored in this paper. Therefore, the necessary condition for our results to explain terrestrial planets in the solar system is that Sd ≲ 0.0018.

Table 2

Statistical measures for final planetary systems.

Table 2 shows the statistical measures for all runs. In the cases of model 1 (1a, 1b, 1c) and model 2 (2a, 2b, 2c), Sc is kept small from the beginning of the simulations because the planets do not undergo convergent migration. The initial value of Sc,0 is 5.4. In the results of model 3 (3a, 3b, 3c), Sc is slightly increased to ≃9. The final system is more concentrated than the initial distribution, but the value is still not consistent with the value of the current solar system. In the results of model 4 (4a, 4b, 4c), Sc is as large as the value of the inner solar system (Sc is between 30 and 70). The mass concentration at around 1 AU is clearly seen in Fig. 3d.

Regarding the angular momentum deficit, Sd for model 1 is small because planets do not exhibit orbit crossings. In contrast, in model 2 and model 3, planets experience orbit crossings and the gas disk is depleted early, thus Sd is larger than that of model 1. We neglect the forced eccentricity by Jupiter and Saturn such that additional damping may be required for model 2 and model 3. In the results of model 4, Sd is sufficiently small (≃ 4 × 10-5).

The timing of a moon-forming impact is also a key constraint on formation models of the solar system. According to recent work, the time of the last giant impact on Earth could have been significantly late (50150 Myr after the formation of the first solids in the solar system; e.g., Touboul et al. 2007; Allégre et al. 2008; Jacobson et al. 2014). The timing of the latest impact for model 2 is t ~ 50 Myr, and that for model 3 is t ≃ 50−100 Myr, consistent with the timing of the last giant impact on Earth. Particularly interesting is that the timing of the last giant impacts is delayed after 100 Myr for model 4. In fact, we continue simulation for 4a until t = 200 Myr and find that the last giant impact occurs at t = 145 Myr.

In conclusion, we find that a disk wind can increase the radial mass concentration. In particular, Sc can be as large as that of the solar system even if we assume that the initial embryo distribution is extended to 0.1 AU. In the result of 4a, the orbital properties including the timing of the moon-forming impact are similar to those of the solar system. We have only performed simulations in a limited number of cases; however, we can identify a significantly large parameter space that may allow us to find a best-fit model for the observed properties of the solar system. To assess the possibility, we have to explore a broad range of parameter space to find a condition that can create a solar system analog.

We also discuss one other problem in the formation of the solar system. Recent numerical simulations on the evolution of the snow line in the solar system that solve the detailed radiative energy transfer have shown that the snow line moves in the disk and comes inside Earth’s orbit to ≃0.5 AU (e.g., Oka et al. 2011; Davis 2005; Graud & Lin 2007) during a phase with a disk accretion rate of ~10-10M yr-1. If planetesimals formed during this stage, the terrestrial planets in the solar system should have formed from icy planetesimals and contain significant amounts of water, which is not consistent with the current water content of the terrestrial planets. We can provide a pathway to avoid this problem. Protoplanets that formed from rocky planetesimals inside the snow line within ≲ 0.5 AU migrate outward in disks that are affected by a disk wind, resulting in water-poor planets at around 1 AU.

5. Conclusions

We performed N-body simulations that include the effects of a disk wind and found that the orbital evolution of terrestrial planets and their precursor protoplanets can significantly differ from the results of previous simulations. Because of the disk wind, the efficiency of which depends approximately on r3/2, gas materials blow out from the surface of the disk from the inside out, leading to a shallower surface density profile in the inner region (≲ 1 AU). Planets with masses less than a few Earth masses would no longer migrate inward. The orbital evolution depends on the parameters (α, Cw), and in fact, 0.1 Earth-mass embryos inside 1 AU never undergo significant inward migration for all the parameters in our simulations. In an extreme case in which the vertical magnetic field is strong, we observed the inner edge of the disk move outward and sweep up embryos.

In a certain parameter range, we found that planetary embryos move outward. When α/Cw is less than ~100, the surface density slope becomes positive in the inner region, which can induce convergent migration to around 1 AU. In order for outward migration to take place, the corotation torque should not be saturated, which depends on diffusion. For 0.1 Earth-mass embryos, the saturation is inhibited when α is ≃10-5. Thus, we performed additional simulations with α = 2 × 10-5 and Cw = 5 × 10-7 and found that the resultant systems are radially concentrated at around 1 AU. This is an important result. In this case, the observed constraints of the inner solar system (e.g., eccentricity, radial mass concentration, late moon-forming impact) may be reproduced. This expectation requires that the turbulent viscosity be quite low (α ~ 10-5) in the terrestrial planet-forming region in the protosolar nebula. Although we focus on the properties of the terrestrial planets, other observed properties in the solar system (e.g., the dynamically excited asteroid belt) should be reproduced in future work. A combination of our model with the grand tack model (Walsh et al. 2011), which explains asteroid excitation, may be a possible scenario. In addition, a disk wind can play other important roles. For example, if the snow line had been located inside 1 AU, the Earth would have accreted a significant amount of icy planetesimals. However, this problem can be overcome if the Earth formed from rocky protoplanets that migrated outward owing to the effect of the disk wind.

The disk wind can also apply to exoplanet formation models. A large number of close-in low-mass planets have been discovered, and thus we can statistically compare observational

results with theoretical models. To do so, we need to explore parameter space (α,CW), which we leave for future work. If disk winds operate efficiently in exoplanet systems, super-Earths would tend to pile up at around 1 AU rather than 0.1 AU.

Acknowledgments

We thank the anonymous referee for helpful comments. M.O. thanks Shigeru Ida and Alessandro Morbidelli for helpful discussions. We also thank Jennifer M. Stone for valuable comments. Numerical computations were in part conducted on the general-purpose PC farm at the Center for Computational Astrophysics, CfCA, of the National Astronomical Observatory of Japan. H.K. gratefully acknowledges support from Grant-in-Aid for Scientific Research (B) (26287101). S.I. is supported by Grant-in-Aid for Scientific Research (23244027, 23103005).

References

  1. Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006, MNRAS, 369, 229 [NASA ADS] [CrossRef] [Google Scholar]
  2. Alexander, R. D., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VI, eds. H. Beuther, R. Klessen, C. Dullemond, & Th. Henning (Tucson: Univ. Arizona Press) [Google Scholar]
  3. Allègre, C. J., Manhès, G., & Göpel, C. 2008, Earth Planet. Sci. Lett., 267, 386 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bai, X.-N., & Stone, J. M. 2013, ApJ, 767, 30 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bans, A., Königl, A., & Uribe, A. 2015, ApJ, 802, 55 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bitsch, B., & Kley, W. 2010, A&A, 523, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Chambers, J. E. 2001, Icarus, 152, 205 [NASA ADS] [CrossRef] [Google Scholar]
  9. Chambers, J. E., Wetherill, Q. W., & Boss, A. P. 1996, Icarus, 119, 261 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  10. Cossou, C., Raymond, S. N., Hersant, F., & Pierens, A. 2014, A&A, 569, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Cresswell, P., & Nelson, R. P. 2008, A&A, 482, 677 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Davis, S. S. 2005, ApJ, 620, 994 [NASA ADS] [CrossRef] [Google Scholar]
  13. Ercolano, B., Clarke, C. J., & Drake, J. J. 2009, ApJ, 699, 1639 [NASA ADS] [CrossRef] [Google Scholar]
  14. Fendyke, S. M., & Nelson, R. P. 2014, MNRAS, 437, 96 [NASA ADS] [CrossRef] [Google Scholar]
  15. Fromang, S., Terquem, C., & Nelson, R. P. 2005, MNRAS, 363, 943 [NASA ADS] [CrossRef] [Google Scholar]
  16. Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606 [NASA ADS] [CrossRef] [Google Scholar]
  18. Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  19. Hansen, B. M. 2009, ApJ, 703, 1131 [NASA ADS] [CrossRef] [Google Scholar]
  20. Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  21. Hellary, P., & Nelson, R. P. 2012, MNRAS, 419, 2737 [NASA ADS] [CrossRef] [Google Scholar]
  22. Horn, B., Lyra, W., Mac Low, M. -M., & Sándor, Zsolt 2012, ApJ, 750, 34 [NASA ADS] [CrossRef] [Google Scholar]
  23. Jacobson, S. A., & Morbidelli, A. 2014, Phil. Trans. R. Soc. A, 372, 20130174 [NASA ADS] [CrossRef] [Google Scholar]
  24. Jacobson, S. A., Morbidelli, A., Raymond, S. N., et al. 2014, Nature, 508, 84 [NASA ADS] [CrossRef] [Google Scholar]
  25. Kominami, J., & Ida, S. 2002, Icarus, 157, 43 [NASA ADS] [CrossRef] [Google Scholar]
  26. Kretke, K. A., & Lin, D. N. C. 2012, ApJ, 755, 74 [NASA ADS] [CrossRef] [Google Scholar]
  27. Laskar, J. 1997, A&A, 317, L75 [NASA ADS] [Google Scholar]
  28. Lesur, G., Ferreira, J., & Ogilvie, G. I. 2013, A&A, 550, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141 [NASA ADS] [Google Scholar]
  30. Morishima, R., Schmidt, M. W., Stadel, J., & Moore, B. 2008, ApJ, 685, 1247 [NASA ADS] [CrossRef] [Google Scholar]
  31. Morishima, R., Stadel, J., & Moore, B. 2010, Icarus, 207, 517 [NASA ADS] [CrossRef] [Google Scholar]
  32. Muto, T., Machida, M. N., & Inutsuka, S. 2008, ApJ, 679, 813 [NASA ADS] [CrossRef] [Google Scholar]
  33. Nagasawa, M., Lin, D. N. C., & Thommes, E. 2005, ApJ, 635, 578 [NASA ADS] [CrossRef] [Google Scholar]
  34. Ogihara, M., Ida, S., & Morbidelli, A. 2007, Icarus, 188, 522 [NASA ADS] [CrossRef] [Google Scholar]
  35. Ogihara, M., Duncan, M. J., & Ida, S. 2010, ApJ, 721, 1184 [NASA ADS] [CrossRef] [Google Scholar]
  36. Ogihara, M., Kobayashi, H., & Inutsuka, S. 2014, ApJ, 787, 172 [NASA ADS] [CrossRef] [Google Scholar]
  37. Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 1413 [NASA ADS] [CrossRef] [Google Scholar]
  38. Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65 [NASA ADS] [CrossRef] [Google Scholar]
  39. Paardekooper, S.-J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950 [NASA ADS] [CrossRef] [Google Scholar]
  40. Paardekooper, S.-J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293 [NASA ADS] [CrossRef] [Google Scholar]
  41. Raymond, S. N., O’Brien, D. P., Morbidelli, A., & Kaib, N. A. 2009, Icarus, 203, 644 [NASA ADS] [CrossRef] [Google Scholar]
  42. Suzuki, T. K., & Inutsuka, S. 2009, ApJ, 691, L49 [NASA ADS] [CrossRef] [Google Scholar]
  43. Suzuki, T. K., Muto, T., & Inutsuka, S. 2010, ApJ, 718, 1289 [NASA ADS] [CrossRef] [Google Scholar]
  44. Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [NASA ADS] [CrossRef] [Google Scholar]
  45. Terquem, C. E. J. M. L. J. 2003, MNRAS, 341, 1157 [NASA ADS] [CrossRef] [Google Scholar]
  46. Touboul, M., Kleine, T., Bourdon, B., Palme, H., & Wieler, R. 2007, Nature, 450, 1206 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  47. Uribe, A., Bans, A., & Königl, A. 2015, ApJ, 802, 54 [NASA ADS] [CrossRef] [Google Scholar]
  48. Walsh, K. J., Morbidelli, A., Raymond, S. N., O’Brien, D. P., & Mandell, A. M. 2011, Nature, 475, 206 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  49. Ward, W. R. 1986, Icarus, 67, 164 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

List of parameters of each model.

Table 2

Statistical measures for final planetary systems.

All Figures

thumbnail Fig. 1

Evolution of gas surface density profile for a) model 2; b) model 3; and c) model 4. Solid lines in each panel indicate the gas surface density at t = 0,0.1 Myr,1 Myr, and 3 Myr, respectively. Dotted lines show results for model 1. In panels a) and b), the solid and dotted lines overlap for t = 0 because of the initial conditions.

Open with DEXTER
In the text
thumbnail Fig. 2

Time evolution of migration rate and direction for bodies with e = 0.01 for a) model 2 and b) model 4. The color indicates the migration timescale. The contours show migration timescales of 0.1 Myr, 1 Myr, and 3 Myr.

Open with DEXTER
In the text
thumbnail Fig. 3

Orbital evolution of planetary embryos for a) model 1; b) model 2; c) model 3; and d) model 4. The filled circles connected with lines represent each body. The color of lines indicates the eccentricity (see color bar). The size of circles is proportional to the radius of the body.

Open with DEXTER
In the text
thumbnail Fig. 4

Final orbital configurations for a) model 2; b) model 3; and c) model 4. Only for model 4a is the configuration at t = 200 Myr shown. Results of three different runs are plotted in each panel.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.