© ESO, 2015

1. Introduction

Planets with masses lower than approximately 50 M_⊕ (depending on the disk scale height and viscosity) migrate toward the central star under a type I migration regime. Different from linear theories, planet formation simulations generally require the rate of type I migration to be reduced by at least a factor of ten and throughout the disk in order to reproduce the distribution of orbital distances of known exoplanets (Ida & Lin 2008; Ogihara & Ida 2009).

It has been shown that type I migration can be locally outward depending on the disk properties (e.g., Kretke & Lin 2012; Bitsch et al. 2015), which would change the picture of planet formation (e.g., Hellary & Nelson 2012; Cossou et al. 2014). The region of local outward migration arises from inhomogeneities in disks (e.g., opacity transitions). According to recent numerical simulations, it has been suggested that local inhomogeneities may help in reducing the type I migration speed (Dittkrist et al. 2014; Mordasini et al. 2015). However, weakening of type I migration over a wide region of the disk for a long time would be required to reproduce the observed distributions of exoplanets (Ida & Lin 2008; Ogihara & Ida 2009); local traps due to disk inhomogeneities would be insufficient.

Recent studies have shown that turbulence-driven disk winds, in which gas is blown away from the surface of the disk, can alter the density profile of the gas disk (Suzuki & Inutsuka 2009; Suzuki et al. 2010), which can slow down or even reverse type I migration. Ogihara et al. (2015b, hereafter OKIS15) performed N-body simulations of terrestrial planet formation in disks including disk winds and found that type I migration can be weakened or even reversed. They also demonstrated that characteristic features of the solar system’s terrestrial planets (e.g., a mass concentration around 1 au) can be reproduced by simulations with disk winds. We anticipate that disk winds play an important role in reproducing observed orbital distributions of exoplanets by slowing type I migration over a wide range of disks for a long time.

In this work we revisit the in situ formation of close-in super-Earths in disks affected by winds. In our previous work (Ogihara et al. 2015a, hereafter OMG15), we reassessed the in situ formation of close-in super-Earths using N-body simulations and observed that super-Earths undergo rapid inward migration, resulting in compact configurations near the disk inner edge, which do not match the observed distributions of super-Earths. On the other hand, we performed additional simulations in which migration is about 100 times slower. The results matched the observations much better. However, the reduction of the type I migration rate in OMG15 was just artificial and there is no physical explanation. Here we investigate whether it can be justified for disks with winds, thus providing an explanation for the observed distribution of super-Earths.

In this paper, we examine the condition for the onset of slow migration. Because of the lack of studies of disk winds, the correlation between the strength of disk winds and the resulting disk surface density slope (or type I migration rate) has not been determined. We first investigate this by numerical experiments in Sect. 2. Then we perform N-body simulations of in situ formation of close-in super-Earths in a disk that evolves through disk winds in Sect. 3. In Sect. 4 we give a summary.

2. Condition for weakening of type I migration

We first investigated the gas surface density slope and the type I migration rate for a wide range of parameters. We numerically solved the following diffusion equation using the same recipe as described in Suzuki et al. (2010) and OKIS15, $\frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{g}}}{\mathit{\partial t}}\mathrm{=}\frac{\mathrm{3}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[{\mathit{r}}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}\mathit{\nu }{\mathrm{\Sigma }}_{\mathrm{g}}{\mathit{r}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\right)}\right]\mathrm{-}{\mathit{C}}_{\mathrm{w}}\frac{{\mathrm{\Sigma }}_{\mathrm{g}}}{\sqrt{\mathrm{2}\mathit{\pi }}}\mathrm{\Omega }\mathit{,}$(1)where Ω is the Keplerian frequency and ν( = αc_sH) is the viscosity. We used the α-prescription for the viscosity, where c_s and H indicate the sound velocity and the disk scale height, respectively. The disk wind flux (ρv_z) can be expressed as ${\mathit{C}}_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{c}}_{\mathrm{s}}\mathrm{=}{\mathit{C}}_{\mathrm{w}}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\Omega }\mathit{/}\sqrt{\mathrm{2}\mathit{\pi }}$ using the mid-plane density ρ₀ and a non-dimensional constant C_w (Suzuki et al. 2010). The initial condition for the gas disk is Σ_g = 2400(r/ 1 au)^{− 3/2}exp(−r/ 50 au) g cm^-2. The temperature profile is assumed to be that of Hayashi (1981) as T = 280(r/ 1 au)^{− 1/2} K.

Fig. 1 Evolution of gas surface density profile for t = 0.01 Myr, 0.1 Myr, and 1 Myr. Dotted lines show the case of weak disk winds (Kw = 200). Solid lines indicate the case of strong disk winds (Kw = 40).

Figure 1 shows examples of the gas surface density evolution. To highlight the difference between the simulation including disk winds and those in OMG15 in the next section, a disk inner edge at r = 0.1 au is superposed to the gas surface density. It is readily seen that the disk profile is altered, especially in the close-in region. The dotted lines indicate the disk evolution for α = 10^-3 and C_w = 5 × 10^-6. The surface density slope is gentle inside r = 1 au and almost flat at r = 0.1 au. Solid lines represent the evolution for α = 10^-4 and C_w = 2.5 × 10^-6. The slope of the surface density of the gas is positive inside r = 1 au.

We here introduce a parameter K_w( ≡ α/C_w) for later discussion. The mass transport rate due to viscous transport and the mass-loss rate due to disk winds in an annulus with Δr are

$\begin{array}{ccc}& & \mathrm{\Delta }\mathit{Ṁ}\mathrm{vis}\mathrm{=}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[\mathrm{-}\mathrm{3}\mathit{\pi }\left({\mathrm{\Sigma }}_{\mathrm{g}}\mathit{\nu }\mathrm{+}\mathrm{2}\mathit{r}\frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{g}}\mathit{\nu }}{\mathit{\partial r}}\right)\right]\mathrm{\Delta }\mathit{r,}\\ & & \mathrm{\Delta }\mathit{Ṁ}\mathrm{wind}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{\pi r}{\mathit{C}}_{\mathrm{w}}\frac{{\mathrm{\Sigma }}_{\mathrm{g}}}{\sqrt{\mathrm{2}\mathit{\pi }}}\mathrm{\Omega \Delta }\mathit{r,}\end{array}$

respectively. When ΔṀ_vis> ΔṀ_wind the disk evolution is dominated by the viscous transport. Here, $\begin{array}{ccc}\frac{\mathrm{\Delta }\mathit{Ṁ}\mathrm{vis}}{\mathrm{\Delta }\mathit{Ṁ}\mathrm{wind}}\mathrm{\simeq }\frac{\mathrm{9}\sqrt{\mathit{\pi }}}{\mathrm{2}}{\left(\frac{\mathit{H}}{\mathit{r}}\right)}^{\mathrm{2}}\frac{\mathit{\alpha }}{{\mathit{C}}_{\mathrm{w}}}\mathrm{\simeq }\mathrm{0.02}{\left(\frac{\mathit{r}}{\mathrm{1}\mathrm{au}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{K}}_{\mathrm{w}}\mathit{,}& & \end{array}$(4)meaning that disk winds become significant inside 1 au when K_w ≲ 100. Disk evolution of dashed lines and solid lines in Fig. 1 correspond to K_w = 200 and 40, respectively.

Next, by performing a number of simulations for a wide range of parameters, we determined the efficiency of type I migration. The surface density slope was determined by viscous diffusion and mass loss due to disk winds; a key parameter is K_w. The total torque for type I migration is given by $\begin{array}{ccc}\mathrm{\Gamma }\mathrm{=}\frac{\mathit{\beta }}{\mathrm{2}}\left(\frac{\mathit{M}}{{\mathit{M}}_{\mathrm{\ast }}}\right)\left(\frac{{\mathrm{\Sigma }}_{\mathrm{g}}{\mathit{r}}^{\mathrm{2}}}{{\mathit{M}}_{\mathrm{\ast }}}\right){\left(\frac{{\mathit{c}}_{\mathrm{s}}}{{\mathit{v}}_{\mathrm{K}}}\right)}^{-2}\mathit{M}{\mathit{v}}_{\mathrm{K}}^{\mathrm{2}}\mathit{,}& & \end{array}$(5)where β,M,M_∗, and v_K are a coefficient that determines the direction and rate of type I migration, the mass of the planet, the mass of the host star, and the Keplerian velocity, respectively. For details of expression of β, we refer to Eqs. (11)–(13) in OKIS15, which is based on Paardekooper et al. (2011). When the surface density slope and the temperature gradient are given, β is determined by the saturation of corotation torque. The level of saturation is expressed by the parameter $\begin{array}{ccc}{\mathit{P}}_{\mathit{\nu }}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\sqrt{\frac{\mathrm{\Omega }{\mathit{r}}^{\mathrm{2}}{\mathit{x}}_{\mathrm{s}}^{\mathrm{3}}}{\mathrm{2}\mathit{\pi \nu }}}\mathit{,}{\mathit{x}}_{\mathrm{s}}\mathrm{=}\frac{\mathrm{1.1}}{{\mathit{\gamma }}^{\mathrm{1}\mathit{/}\mathrm{4}}}\sqrt{\frac{\mathit{M}}{{\mathit{M}}_{\mathrm{\ast }}}\frac{\mathit{r}}{\mathit{H}}}\mathit{,}& & \end{array}$(6)where x_s is the dimensionless half-width of the horseshoe region and γ is the adiabatic index. When the temperature distribution is fixed, P_ν is determined by viscosity and planetary mass. The saturation of the entropy-related corotation torque is determined by the thermal diffusivity ξ, and we assumed ξ = ν for simplicity.

Fig. 2 Surface density slope in a steady state at r = 0.1 au and 1 au for various values of Kw (panel a)). Migration efficiency (Γ/(−ΓTTW)) for planets with e = 0.01 at r = 0.1 au (panel b)) and r = 1 au (panel c)). The contours show a migration efficiency of −1, −0.3, −0.1, 0.1, 0.3, and 1. When Γ > 0, planets move outward.

Figure 2a represents the surface density slope at r = 0.1 au (the surface density slope just beyond the disk edge at 0.1 au) and 1 au. Figures 2b and c show migration maps of r = 0.1 au and 1 au, respectively. The color scale indicates the migration efficiency as compared to the migration rate in a locally isothermal disk with a power-law index of −3/2 derived by a three-dimensional linear analysis by Tanaka et al. (2002, hereafter TTW02), which is defined by Γ/(−Γ_TTW)¹. Γ_TTW has a negative value, so Γ > 0 indicates outward migration. This value corresponds to the negative of the efficiency parameter, −C₁, used in Ida & Lin (2008). We included the dependence of the corotation torque on the eccentricity by assuming e = 0.01 (Fendyke & Nelson 2014). We adopted a higher value for C_w( = 10^-4) in Fig. 2 to reduce the computation time. As already stated above, the surface density slope is determined only by K_w, which is confirmed by numerical experiments that adopt different values for C_w. The surface density slope relaxes to a steady state after t = t_ν = r²/ν, where t_ν is the viscous timescale, so the maps are plotted at t = 10⁵ yr.

From Fig. 2a, we find that the surface density slope decreases as K_w increases. Thus, a smaller K_w yields slower or even outward migration. At r = 0.1 au for K_w ≲ 150, we can find a range of P_ν in which the migration speed is reduced by a factor of more than ten relative to that predicted by TTW02. As r increases, the deviation of the slope from that of initial power-law disks decreases. At r = 1 au, the migration can be reduced by a factor of ten from TTW02 for K_w ≲ 70. Thus, disk winds are able to modify type I migration in a wide region inside 1 au.

The migration rate depends not only on K_w but also on P_ν. The value of P_ν increases with decreasing viscosity; the corotation torque saturates at low viscosity, and the Lindblad torque dominates the total torque. On the other hand, it is also known that there exists a cut-off for the horseshoe drag at high viscosity (small P_ν) and the corotation torque approaches its linear value (e.g., Masset 2002; Paardekooper & Papaloizou 2009a).

Using Fig. 2, we can estimate the migration rate for different sets of parameters (K_w,P_ν). According to Suzuki et al. (2010), a possible value of K_w would be ~100; however, further investigation would be required by global magnetohydrodynamics simulations that cover enough grid points in the vertical direction to constrain the range.

Our results depend strongly on the fact that K_w is constant with time and radius. As discussed by Suzuki et al. (2010), the gravitational energy is released by gas accretion, and a part of the energy is used to drive winds. Thus C_w is proportional to the kinetic energy of winds, while α is proportional to the accretion rate. Thus K_w = const is a natural assumption². We note that although we assumed that gas blown out of the disk surface escapes from the disk, some materials may return to the disk. Stellar winds can push away the lifted-up gas (Suzuki et al. 2010); however, if a large amount of gas returns to the disk, the surface density slope would be smaller than in Fig. 2a. Moreover, we assumed a weak vertical magnetic field (β_z ≳ 10⁴; β_z is the vertical component of plasma β at the midplane).

As discussed above, our model uses the same thermal profile for the disk as in OMG15. Clearly, this is a simplification that we adopted to compare the results to those of OMG15 more directly. However, we checked that our main results would hold with a more realistic temperature profile (Bitsch et al. 2015; B. Bitsch, priv. comm.). Specifically, we found that while the temperature and scale height are lower than for our model for α = 10^-4, the changes in the migration maps remain limited compared to those shown in Fig. 2. In particular, we confirm that type I migration is still suppressed for some value of P_ν for K_w< 100. Furthermore, if a density gap is opened by planets, the formulae we used would not be valid. While in our disk model only massive super-Earths (≳ 20 M_⊕) would open a gap, in the disk with a more realistic temperature profile (i.e., colder), even lower-mass planets (≳2 M_⊕) open a gap. However, in this case the planets would migrate in the type II regime, which, for the value of α = 10^-4 that we assumed, would result in a slower migration rate than in our disk model. Thus the results we present in the next section concerning the effects of a reduced migration speed are conservative in the sense that the migration speed of the most massive planets would be even lower in a disk with a more realistic temperature model.

3. In situ formation of close-in super-Earths

We now perform edN-body simulations of formation of close-in super-Earths in disks that evolve through disk winds. The simulation model is the same as that of OMG15, except for the gas disk model. According to a power-law distribution, 250 embryos with a mass of 0.1 M_⊕ and 1250 planetesimals with a mass of 0.02 M_⊕ were distributed between 0.1 and 1 au. The total mass in the system was set to 50 M_⊕. This high solid-to-gas mass ratio can be explained by invoking the radial drift of dust and pebbles into the inner part of the disk, but we place ourselves here at a stage where most of the dust has already been converted into large bodies (embryos and planetesimals). Planetesimals suffer aerodynamical gas drag assuming a physical size of 50 km in radius (Adachi et al. 1976), while planets with more than roughly 0.1 M_⊕ undergo the tidal damping of eccentricities, inclinations, and semimajor axes (see Ogihara et al. 2014 and OKIS15 for each formula).

To show the effects of disk winds in a disk model, we assumed that disk winds are relatively strong (K_w = 40), where α = 10^-4 and C_w = 2.5 × 10^-6 were used. The evolution of the gas surface density is shown by solid lines in Fig. 1. In this disk, we see from Fig. 2 that type I migration can be slower by a factor of more than ten from TTW02 or even reversed for 0.1 ≲ P_ν ≲ 1. This range roughly corresponds to M ~ 0.1−1 M_⊕ between r = 0.1 and 1 au. We assumed a rapid disk dispersal (~ 0.1 Myr) after a typical disk lifetime of 3 Myr to be consistent with observations.

Fig. 3 Time evolution of planets for a typical run. The filled circles connected with solid lines represent the sizes of planets. The smallest circle represents an embryo of 0.2 Earth-mass, while the largest circle represents a 20 Earth-mass planet. The color of the lines indicates the eccentricity (color bar).

We performed ten runs with different initial positions of solid bodies; a typical run is shown in Fig. 3. As seen in OMG15, the growth of embryos is quite rapid. However, we find that the migration speed is significantly slower than in the fiducial model (model 1) in OMG15 (see Fig. 2 in OMG15). The actual migration timescale of massive planets (~ 5 M_⊕) is ≳ 0.1 Myr between t = 0.01 and 0.1 Myr. This migration rate is about a few to ten times slower than that predicted by TTW02. At t ≃ 0.2 Myr, planets are captured in mutual mean motion resonances, making a long resonance chain (nine bodies). The chain undergoes orbital instability at about 4 Myr after gas dispersal, leading to mutual collisions and a relatively separated system of three planets out of resonances. Compared with the results of slow-migration case (model 3) in OMG15, the migration speed is faster, and hence there are fewer planets in the resonant chain in our results. However, the final phase of the planet formation process, in which planets undergo close encounters and collisions after gas dispersal, is quite similar.

Fig. 4 Comparison of cumulative period ratio distributions and eccentricity distributions. Thin solid lines show observed distributions of confirmed close-in super-Earths as of June 2015 (341 systems with 854 planets). Thick solid lines indicate results of simulations including the effects of disk winds, while thick dashed lines show the previous results for model 1 in OMG15. As discussed in OMG15, observed eccentricities can be overestimated (e.g., Shen & Turner 2008; Zakamska et al. 2011). The thin dotted line shows the eccentricity distribution, in which each eccentricity is assumed to be e−σ. Here, σ is the estimated error.

Figure 4 compares the results of simulations with the observed distributions of period ratios of adjacent pairs and eccentricities. For reference, the results of a fiducial model (model 1) in OMG15 are also plotted. We find that the period ratio distribution of simulations including disk winds matches the observations much better than that of model 1 in OMG15. The eccentricity distribution also matches the observations well. We acknowledge, however, that the observed mass distribution is very poorly matched by our simulations; the averaged slope of solid surface density in our simulation, which is deduced from the distribution of the final planets, is ~−3, which is significantly steeper than the averaged slope of close-in super-Earths (≃−1.5). This is presumably because we assumed an initial distribution of solids that is confined to a region between 0.1 and 1 au, so the slope would be improved if planets with M ≳ M_⊕ (in this case Γ < 0 at 1 au) migrated from outside 1 au.

Here we briefly discuss a successful scenario of in situ formation of close-in super-Earths. According to the results of N-body simulations presented in this and previous papers, in a successful model planets are captured in a resonance chain in a disk and then undergo orbit crossings and collisions during disk dissipation. On the other hand, unsuccessful models include the rapid migration case (e.g., model 1 in OMG15) and the no migration case (e.g., model 4 in OMG15). In the former case, planets formed in a compact resonance chain with a small number of planets (N ≃ 5). The small number of planets in the chain prevents close encounters after gas dispersal, leading to mismatches in the period ratio (too compact) and eccentricity (too low). In the latter case, planets undergo a too violent instability because they are not in a resonant chain, which also results in mismatches to observations (too separated and high eccentricities).

To realize the successful scenario described above, the number of planets in a resonance chain should be large enough to trigger orbit crossings during disk dissipation. According to a study on the orbital stability of a resonance chain (Matsumoto et al. 2012), there should be more than five to ten planets in a chain. It is difficult to precisely assess the conditions for the formation of a resonant chain with N> 5−10; however, results of N-body simulations imply that type I migration should be reduced by a factor of about ten from that predicted by the linear theory. According to Fig. 2, this condition corresponds to K_w ≲ 100 and 0.1 ≲ P_ν ≲ 1.

4. Summary

We computed the surface density profile of disks affected by winds of various strengths and the resulting type I migration rates. We confirmed that the type I migration can be slowed down in the whole close-in region (r< 1 au) if the wind is sufficiently strong. Using the migration map in Fig. 2, the migration rate for different sets of parameters can be estimated without the need for additional calculations. We also performed N-body simulations of the formation of close-in super-Earths that included the effects of disk winds. Given that the effects of disk winds are relatively strong, we demonstrated that type I migration is significantly slowed down. This is the first simulation in which the observed statistical orbital distributions of close-in super-Earths are reproduced by results of N-body simulations without applying an artificial reduction of type I migration.

Acknowledgments

We thank the anonymous referee for helpful comments and T. Suzuki for valuable discussions. This work was supported by ANR, project number ANR-13–13-BS05-0003-01 projet MOJO (Modeling the Origin of JOvian planets).

References

All Figures

	Fig. 1 Evolution of gas surface density profile for t = 0.01 Myr, 0.1 Myr, and 1 Myr. Dotted lines show the case of weak disk winds (Kw = 200). Solid lines indicate the case of strong disk winds (Kw = 40).
In the text

	Fig. 2 Surface density slope in a steady state at r = 0.1 au and 1 au for various values of Kw (panel a)). Migration efficiency (Γ/(−ΓTTW)) for planets with e = 0.01 at r = 0.1 au (panel b)) and r = 1 au (panel c)). The contours show a migration efficiency of −1, −0.3, −0.1, 0.1, 0.3, and 1. When Γ > 0, planets move outward.
In the text

	Fig. 3 Time evolution of planets for a typical run. The filled circles connected with solid lines represent the sizes of planets. The smallest circle represents an embryo of 0.2 Earth-mass, while the largest circle represents a 20 Earth-mass planet. The color of the lines indicates the eccentricity (color bar).
In the text

Fig. 4 Comparison of cumulative period ratio distributions and eccentricity distributions. Thin solid lines show observed distributions of confirmed close-in super-Earths as of June 2015 (341 systems with 854 planets). Thick solid lines indicate results of simulations including the effects of disk winds, while thick dashed lines show the previous results for model 1 in OMG15. As discussed in OMG15, observed eccentricities can be overestimated (e.g., Shen & Turner 2008; Zakamska et al. 2011). The thin dotted line shows the eccentricity distribution, in which each eccentricity is assumed to be e−σ. Here, σ is the estimated error.

In the text