© ESO 2018

1 Introduction

The properties of the terrestrial planets of our solar system are well understood compared to other systems. A number of N-body simulations of planetary accretion were performed to reproduce these properties, such as small eccentricities and late moon formation (e.g., Nagasawa et al. 2005; Ogihara et al. 2007; Raymond et al. 2009; Morishima et al. 2010).

One of the major issues facing the theory of planet formation is reproducing the current orbital configuration of the solar system’s terrestrial planets. Their orbital locations are confined to a region between 0.38 and 1.5 au from the Sun. Several models have been proposed to reproduce the outer boundary of the localized distribution which includes the grand tack model (e.g., Walsh et al. 2011). In the meantime, it is quite hard to explain why there is no large planets in close-in orbits. If planetesimals form when the condensation line for silicates is located inside r = 0.7 au, protoplanets/planets should form in close-in orbits. However, no planets exist inside the orbit of Mercury in our solar system. Regarding this problem, Ogihara et al. (2015) demonstrated that protoplanets larger than Mars that form in close-in orbits can undergo outward type I migration in disks with positive surface density profiles, which would help to reproduce the deficit of close-in planets.

Based on three-dimensional (3D) magnetohydrodynamic (MHD) simulations, Suzuki et al. (2016) derived the long-term global disk evolution including effects of magnetically driven disk winds. They considered the mass loss due to disk winds (Suzuki & Inutsuka 2009) and the wind-driven accretion (Bai & Stone 2013; Gressel et al. 2015; magnetic braking) in addition to the standard viscous evolution. They found that the disk profile in the close-in region (r ≲ 1 au) can be significantly altered from the minimum-mass solar nebular (MMSN) model (Weidenschilling 1977; Hayashi 1981). In particular, the slope of the gas surface density is almost flat inside r ≃ 1 au for MRI-active disk, and the slope can even be positive for MRI-inactive disk. In a separate paper (Ogihara et al. 2018, hereafter Paper II), we investigated the orbital evolution of planetary embryos larger than Mars in such disks under various conditions, and found that type I migration of protoplanets can be significantly suppressed in many cases compared to the type I migration rate estimated for the MMSN model. On the other hand, it is expected that smaller particles (dust/planetesimals) can undergo significant radial drift in MRI-inactive disks of Suzuki et al. (2016).

In this work, we examine another path to form the localized configuration of the solar system’s terrestrial planets using the recent disk evolution model of Suzuki et al. (2016). We examine whether or not the following formation scenario can take place by numerical simulations: planetesimals undergo convergent radial drift to form localized planetesimal disks, then planetary embryos grow to terrestrial planets without undergoing significant type I migration.

The paper is organized as follows. In Sect. 2, we present the numerical model of disk evolution. In Sect. 3, we show the outcomes of simulations that consider the radial drift of planetesimals. In Sect. 4, we present results of N-body simulations from planetary embryos. In Sect. 5, we make concluding remarks.

2 Disk model

For evolution of disk surface density, we numerically solve the diffusion equation that includes effects of disk winds based on MHD simulations (the same as used by Suzuki et al. 2016; see also Simon et al. 2015; Hasegawa et al. 2017): $$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{2}{r\Omega }\left\{\frac{\partial }{\partial r}(r^2{\Sigma }_{\text{g}}\overline{{\alpha }_{r,\varphi }}{c}_{\text{s}}^2)+r^2\overline{{\alpha }_{\varphi ,z}}\frac{{\Sigma }_{\text{g}}H{\Omega }^2}{2\sqrt{\pi }}\right\}\right]-C_{\text{w}}\frac{{\Sigma }_{\text{g}}\Omega }{\sqrt{2\pi }},$$(1)

where Σ_g, r, Ω, c_s, and H are the gas surface density, the radial distance, the Keplerian frequency, the sound speed, and the disk scale height, respectively. An effective turbulent viscosity $\overline{{\alpha }_{r,\varphi }}$ is based on a Shakura & Sunyaev (1973) α viscosity parameterization. Parameters $\overline{{\alpha }_{\varphi ,z}}$ and C_w describe measures of the angular momentum loss due to the wind torque and mass loss due to disk winds, respectively. According to results of MHD simulations, we chose the parameters as follows. We adopt two values of the effective turbulent viscosity; namely, $\overline{{\alpha }_{r,\varphi }}=8\times {10}^{-3}$ (MRI-active case) and $\overline{{\alpha }_{r,\varphi }}=8\times {10}^{-5}$ (MRI-inactive case). Other parameters are fixed as $\overline{{\alpha }_{\varphi ,z}}={10}^{-5}{({\Sigma }_{\text{g}}/{\Sigma }_{\text{g},\text{ini}})}^{-0.66}$ (Bai 2013), C_w = 2 × 10⁻⁵ (MRI-active case; Suzuki & Inutsuka 2009) and C_w = 10⁻⁵ (MRI-inactive case; Suzuki et al. 2010). As described in Sect. 2.3 of Suzuki et al. (2016), the mass-loss rate due to disk winds C_w is also constrained by energetics.

The initial gas profile is given by ${\Sigma }_{\text{g},\text{ini}}=1.7\times {10}^4{(r/1 \text{au})}^{-3/2}\exp (-r/30 \text{au}) \text{g} {\text{cm}}^{-2}$. In order to consider the early stage of disk evolution, which is connected to star formation, the initial gas disk is about ten times more massive than the MMSN model. We note that planetesimals and planetary embryos take some time to form, therefore we start simulations some time after the disk formation. In the same way as in our previous works, we start simulations at t = 10⁵ yr in the disk evolution of Suzuki et al. (2016), and refer to this time as “initial”.

Figure 1 shows the time evolution of the gas surface density for MRI-active disks and MRI-inactive disks. The evolution is the same as in Suzuki et al. (2016). We see that disk profiles are altered in the close-in region (r ≲ 1 au). An important point is that the disk obtains a positive surface density slope in a broad close-in region in MRI-inactive disks. As discussed in Suzuki et al. (2016), the angular momentum loss due to the wind torque is inversely correlated with the gas density. In addition, the mass loss rate is proportional to the Keplerian time. Therefore, the gas density is carved out from inside based on the model of Suzuki et al. (2016). It is likely that in the late stage of disk evolution, the disk is MRI active due to smaller gas surface density. For the evolution of the disk temperature, the stellar irradiation and viscous heating are considered. We refer to Sect. 2.4 of Suzuki et al. (2016) for a detailed description.

Fig. 1 Time evolution of gas surface density based on Suzuki et al. (2016). Blue and red lines represent MRI-active disk and MRI-inactive disk, respectively.

3 Radial drift of planetesimals

We first examine the radial drift of planetesimals before the formation of planetary embryos. Particles with the size of planetesimals or smaller undergo radial migration due to gas drag. As stated in Sect. 2, it is considered that the MRI activity increases with decreasing gas surface density. Therefore, here we use the MRI-inactive disk for the orbital evolution phase of planetesimals. As discussed in Suzuki et al. (2016), disk winds affect the radial drift of planetesimals. As the radial pressure gradient force and hence the rotation velocity of gas are altered, planetesimals can feel a tail wind in the region of a positive pressure gradient, leading to outward drift.

Here we perform simulations of the orbital evolution of planetesimals in MRI-inactive disks with mass M = 1 × 10¹⁶ g (R ≃1 km assuming the internal density of 3 g cm⁻³) that are placed between r = 0.3 and 3 au. Mutual interactions between planetesimals (i.e., gravity, collision) are ignored. For aerodynamical gas drag force, we use the formula of Adachi et al. (1976). The drift timescale (see Appendix A) is inversely proportional to a scaling factor of pressure gradient force $\eta =-1/2{(c_{\text{s}}/v_{\text{K}})}^2(\text{dln\hspace{0.17em}}P/\text{dln\hspace{0.17em}}r)$. See Fig. A.1 for the evolution of η. Figure 2 shows the result of our simulation, in which the evolution of solid surface density is indicated. We find that planetesimals undergo convergent radial drift. As the disk evolves, the location of the pressure maximum moves outward. After t = 0.3 Myr, planetesimals accumulate in a narrow ring-like region because all planetesimals inside ~0.7 au undergo outward drift. Therefore, it is possible that a localized configuration is attained by radial drift of planetesimals before formation ofplanetary embryos. We note that typical sizes of planetesimals that form via streaming instability would be larger than 1 km (e.g., Simon et al. 2016). Larger planetesimals suffer less efficient radial drift. Thus, if planetesimals larger than 10 km are distributed between 0.3 and 3 au, the strong convergent drift seen in Fig. 2 does not occur.

Although it is reasonable to assume that the disk is MRI inactive in the early/middle phase of planet formation, we also comment onradial drift of planetesimals in MRI active disks. When the surface density distribution exhibits flat profiles (e.g., MRI-active disk), planetesimals undergo slow inward drift instead of outward drift. In this case, narrow ring-like regions would not form; however, as the drift velocity is slowed down significantly, it would be useful to overcome the drift barrier of dust/planetesimals.

We note that we adopt several simplifications that affect the realism of the simulation. First, a single size population of planetesimals is considered, and their growth is also ignored. In reality, planetesimals have size distribution, which changes with time. Second, mutual interactions between planetesimals and the back reaction of planetesimals on the gas disk are neglected. When solid particles accumulate in a narrow ring, where the dust-to-gas ratio is close to unity, the back reaction on the gas disk should be important. The evolution in Fig. 2 is therefore an idealized result. We will investigate these effects in future works.

Fig. 2 Time evolution of solid surface density. Planetesimals with R ≃ 1 km are placed between r = 0.3 and 3 au. The disk is considered to be MRI-inactive.

4 Orbital evolution of planetary embryos via type I migration

Now we examine the late stage of planet formation from planetary embryos. As we showed in the previous subsection, planetesimals can accumulate in a narrow-ring region. We therefore consider the subsequent evolution for the initial embryo distribution. Hansen (2009) adopted a localized configuration for the initial embryo distribution and performed N-body simulations. We use the same initial conditions as Hansen (2009). As an initial condition of N-body simulations, 400 planetary embryos with mass M = 0.005 M_⊕ are randomly placed between 0.7 and 1 au according to a uniform surface density. Their initial eccentricities and inclinations are set to be small (≃10⁻²). Although Hansen (2009) considered the presence of Jupiter, we do not include any giant planets in our simulations.

For the disk evolution, as we stated in the previous section, it is reasonable to consider that the disk is MRI-inactive in the early phase, which evolves to the MRI-active disk. However, the MRI activity is uncertain when embryos form; therefore we perform N-body simulations of embryos both in MRI-active and -inactive disks. The disk evolution model is already shown in Sect. 2.

For the formulae of type I migration, we use those described in Paardekooper et al. (2011). The unsaturated torque is given by Γ = (−0.6 + 2.3p − 2.8q)Γ₀, where p = dln Σ_g∕dln r is the local surface density slope and q = dln T∕dln r is the local temperature slope (Paardekooper et al. 2010). In our disk model, the temperature slope is roughly given by q = −1∕2 (see Fig. 5 of Suzuki et al. 2016), therefore the rate and direction of type I migration is mainly determined by the surface density slope q. The actual torque is calculated by considering the saturation effect. Here the saturation is determined by the viscous diffusion coefficient ν, which is proportional to $\overline{{\alpha }_{r,\varphi }}$ (see Sect. 4.1 of Ogihara et al. 2017).

Figure 3a shows snapshots of the system for MRI inactive disks (see also Paper II for the migration map). No significant inward migration is seen. Although some protoplanets move outward to r ≃ 2 au, the innermost larger protoplanets (M ≳ 0.1 M_⊕) are located at r = 0.56 au at t = 100 Myr.

Figure 3b shows results of the N-body simulation for MRI-active disks. We find no clear tendency for planets to migrate inward or outward. In the final state (t = 100 Myr), four larger protoplanets remain between r = 0.46 and 1.1 au. The region where planets exist spreads compared to the initial ring-like region. This is not because planets exhibit migration but because they undergo scattering between protoplanets.

For comparison, we perform N-body simulations for the MMSN model. We find that after planetary embryos grow to Mars-mass protoplanets(M ≃ 0.1 M_⊕), they rapidly migrate inward and are lost to the central star with a migration timescale of ≃0.1 Myr. At t = 10 Myr, the total mass of remaining objects is only 0.445 M_⊕. Therefore we can confirm that disk winds play a role in suppressing type I migration, which leads to localized orbital configurations. We note that from geochemical constraints, Mars grew to half of its present size in about 1.8 Myr (Dauphas & Pourmand 2011). If Mars grew on this timescale and the MMSN decays exponentially with a timescale of 1 Myr, Mars would migrate marginally and not fall onto the star.

We have also performed additional simulations that start from protoplanets with an isolation mass instead of embryos with M = 0.005 M_⊕. For the calculation of the isolation mass, see Eq. (13) of Kokubo & Ida (2002) for example. We chose a relatively small initial orbital separation of 5 r_H, where r_H represents the mutual Hill radius of neighboring protoplanets. The initial and final orbital configuration is overplotted in Fig. 3 by circles. We find that our conclusion does not change. No planets remain in the close-in region (a < 0.4 au) and planetary orbits are confined to a localized region. Detailed results of N-body simulations of isolation-mass protoplanets are presented in Paper II.

Fig. 3 Snapshots of the system for MRI-inactive disk (a) and MRI-active disk (b). The filled circles represent the size of particles. When particles grow to more massive than 0.1 M⊕, they are indicated by filled blue circles. The red lines indicate the gas surface density (right axis). The red circles represent the initial and final orbital configurations of simulations from isolation-mass protoplanets.

5 Conclusions

We investigate the formation of the solar system’s terrestrial planets in disks evolving with disk winds. We demonstrate that planetesimals formed in close-in orbits can undergo outward radial drift in MRI-inactive disks and the inner edge of the planetesimal ring is created at r ~ 0.7 au. Subsequently, we examine the orbital evolution and growth of embryos from a narrow ring by N-body simulation, and find that migration of protoplanets can be significantly suppressed. As a result, we confirm that no large planets remain inside ≃0.5 au. Hansen (2009) showed that the current orbital configuration of the solar system’s terrestrial planets can be reproduced if a narrow solid ring is formed in a gas-free environment. Here, in this work, even if we start simulations from an earlier stage when gas is present, we show that the localized configuration can be maintained. In a previous work, Ogihara et al. (2015) showed that protoplanets can undergo outward migration to ≃1 au in a region with a positive surface density slope. In the present work, we showed that even if planetary embryos do not undergo convergent migration, the localized configuration can be created due to convergent radial drift of planetesimals. The reason why all planetesimals that formed in the close-in region undergo outward radial drift, which creates the inner edge of planetesimal ring, is that the gas disk obtains a positive pressure gradient in a broad close-in region, which expands outwards in time. This results in outward drift of all planetesimals in the close-in region, leaving no planetesimals distributed in the broad close-in region. This is qualitatively different from other narrower pressure traps (e.g., the snow line). Such traps can stop inward drift, but cannot drag planetesimals outwards over the full close-in region. If there was no clear disk inner edge, the narrow-trap case would give a ring-like planetesimal distribution because planetesimals inside the pressure trap would fall onto the star. However, if there existed an inner edge, inward drifting planetesimals would be accumulated there; the planetesimal distribution could be very different from our model. We note also that the outer boundary is determined by the initial distribution of planetesimals/embryos. Other models are required (e.g., the grand tack model) to explain the origin of the outer boundary in the solar system’s terrestrial planets.

In Sect. 3, planetesimals are initially distributed between 0.3 and 3 au. Although we do not specify their origin, they would be formed via radial drift of dust and pebbles, which undergo rapid drift. We note here that no pressure bump exists in the early phase of disk evolution of Suzuki et al. (2016), and it takes some time (~0.1 Myr) to form positive slope in the MRI-inactive disk. Therefore, in the earlier stage, pebbles can drift inwards with no limit and planetesimals can form everywhere in the close-in region, which justifies our choice of initial planetesimal distribution. On the other hand, if a pressure bump exists before the planetesimal formation, dust can accumulate at the pressure maximum and planetesimals only form in a local region at the pressure bump. In this case,the planetesimal distribution would differ from that used in Sect. 3. However, our conclusion would still be valid because this situation is consistent with our assumption in Sect. 4 that embryos form in a narrow ring-like region. In this case, the radial drift of planetesimals is not needed and big (>10 km) planetesimals could form. In either case, as we stated in Sect. 3, a simulation of growth and radial drift of dust particles in a consistent manner is required. We note also that an alternative mechanism to create the narrow planetesimal ring by drifting pebbles and their backreaction to the disk has been proposed (Dra̧żkowska et al. 2016).

In addition, disk winds would contribute to a solution of the “snow line problem”. According to radiative transfer calculations, the ice condensation line (snow line) at the midplane moves inside 1 au when the accretion rate of the disk drops below ≃ 10⁻⁹ M_⊙ yr⁻¹ (e.g., Oka et al. 2011). If planetesimals form at this time and the Earth forms in situ, the Earth should accrete a large amount of icy material, and it is hard to explain the origin of the amount of water (0.023 wt%) on the Earth. Here we apply the disk-wind model to this problem. When the snow line moves inside 1 au, we can assume that the pressure bump is located outside 1 au. If planetary embryos or planetesimals that grow to form the Earth were to form at 1 au by the time of the passage of the snow line through 1 au, the pressure maximum outside 1 au would possibly prevent the inward drift of icy pebbles in the vicinity of the Earth. This could be another mechanism of a fossilization of the snow line proposed by Morbidelli et al. (2016). This point should be explored in future works.

Acknowledgements

We thank Cornelis Dullemond for constructive comments. Numerical computations were conducted on the general-purpose PC farm at the Center for Computational Astrophysics, CfCA, of the National Astronomical Observatoryof Japan. This work was supported by JSPS KAKENHI Grant Numbers 16H07415 and 17H01105.

Appendix A: Radial drift due to aerodynamical gas drag

We use a force formula of aerodynamical gas drag on planetesimals developed by Adachi et al. (1976); $$F_{\text{aero}}=-\frac{1}{2}{C}_{\text{D}}\pi R^2{\rho }_{\text{g}}\Delta u\Delta u,$$(A.1)

where C_D, ρ_g, and Δu are the gas drag coefficient, the gas density, and the relative velocity of the planetesimal to the disk gas. The relative velocity is expressed by $$\Delta u\simeq {\left(\frac{5}{8}{e}^2+\frac{1}{2}{i}^2+{\eta }^2\right)}^{1/2}{v}_{\text{K}}.$$(A.2)

Here we use the random velocity relative to the local Keplerian motion and the gas velocity (v_gas = [1 − η]v_K). The stopping time is given by t_stop = MΔu∕|F_aero|. Then the drift timescale is $$\begin{array}{lll}{t}_{\text{drift}}\hfill & \simeq \hfill & \frac{t_{\text{stop}}}{2\left|\eta \right|}\hfill \\ \hfill & \simeq \hfill & 2\times {10}^5{\left(\frac{\eta }{{10}^{-3}}\right)}^{-1}{\left(\frac{C_{\text{D}}}{2}\right)}^{-1}{\left(\frac{M}{{10}^{16} \text{g}}\right)}^{1/3}{\left(\frac{\rho }{3 \text{g\hspace{0.17em}}{\text{cm}}^{-3}}\right)}^{2/3}\hfill \\ \hfill & \hfill & \times {\left(\frac{{\Sigma }_{\text{g}}}{{10}^3 \text{g\hspace{0.17em}}{\text{cm}}^{-2}}\right)}^{-1}\left(\frac{H/r}{0.05}\right){\left(\frac{r}{1 \text{au}}\right)}^{3/2}\hfill \\ \hfill & \hfill & \times {\left(\frac{\sqrt{\frac{5}{8}{e}^2+\frac{1}{2}{i}^2+{\eta }^2}}{0.01}\right)}^{-1}\text{yr},\hfill \end{array}$$

where ρ is the internal density of planetesimals. In Sect. 3, e = 0.01 and i² = e²∕4 are used. To obtain a more accurate orbital evolution of planetesimals, we need to calculate the evolution of e and i, which should be considered in a future work.

Figure A.1 shows the time evolution of η for an MRI-inactive disk. A similar plot is shown in Fig. 10 of Suzuki et al. (2016). The direction of radial drift is determined by the sign of η. Although Eq. (A.3) can be a little more complicated for |η| < e² (see Eq. (4.21) of Adachi et al. 1976), outward drift is clearly seen for η ≲ −1 × 10⁻⁴. Thick lines represent the region of negative η. There exists a region with positive pressure gradient in the inner region (< 1 au), which move outwards in time. Thepositive pressure gradient is seen in region between 0.05 and 0.4 au at the initial time, while the region is between 0.2 and 0.6 au at t = 0.1 Myr. Therefore, all planetesimals inside ~1 au can undergo outward radial drift before t = 1 Myr in results of Sect. 3.

Fig. A.1 Time evolution of the absolute value of the scaling factor of pressure gradient force η for an MRI-inactive disk based on Suzuki et al. (2016). Thick lines correspond to the region where η is negative.

References

All Figures

	Fig. 1 Time evolution of gas surface density based on Suzuki et al. (2016). Blue and red lines represent MRI-active disk and MRI-inactive disk, respectively.
In the text

	Fig. 2 Time evolution of solid surface density. Planetesimals with R ≃ 1 km are placed between r = 0.3 and 3 au. The disk is considered to be MRI-inactive.
In the text

Fig. 3 Snapshots of the system for MRI-inactive disk (a) and MRI-active disk (b). The filled circles represent the size of particles. When particles grow to more massive than 0.1 M⊕, they are indicated by filled blue circles. The red lines indicate the gas surface density (right axis). The red circles represent the initial and final orbital configurations of simulations from isolation-mass protoplanets.

In the text

	Fig. A.1 Time evolution of the absolute value of the scaling factor of pressure gradient force η for an MRI-inactive disk based on Suzuki et al. (2016). Thick lines correspond to the region where η is negative.
In the text