© ESO, 2017

1. Introduction

According to numerous recent studies, global disk structures and evolution are more complicated than previously considered (e.g., Fromang et al. 2011; Flock et al. 2011, 2012; Dzyukevich et al. 2013; Parkin & Biknell 2013; Suzuki & Inutsuka 2014; Béthune et al. 2017). Recently, magnetically-driven disk winds have been shown to have a crucial role in the disk evolution.

Suzuki & Inutsuka (2009) performed local three-dimensional magnetohydrodynamic (MHD) simulations and showed the existence of turbulence-driven disk winds, in which vertical gas outflows are triggered by magnetorotational instability (MRI) turbulence. Suzuki et al. (2010) further performed MHD simulations and found that turbulence-driven disk winds are observed even in MRI-inefficient disks, because disk winds can be driven near the surface region where magnetic fields are more turbulent. They also calculated the global disk evolution with mass loss (vertical gas outflows) due to disk winds, and found that mass loss due to disk winds alter the surface density distribution of gas disks, which can affect planet formation process in several ways. The dynamics and evolution of dust grains can be affected; small grains that are well coupled to the gas component are blown upward and away by the gas drag force, while larger grains that are loosely coupled to the gas remain in the disk (Miyake et al. 2016).

One of the most remarkable effects is on the migration of low-mass planets (type I migration). Ogihara et al. (2015a,b) performed N-body simulations of protoplanets and/or planetesimals in disks evolving via disk winds, which is based on Suzuki et al. (2010), and found that type I migration is suppressed over the whole close-in region when the mass loss due to disk winds is relatively strong. They also demonstrated that type I migration can be even reversed in some cases. In addition, results of N-body simulations indicate that disk winds would help in reproducing observed distributions of terrestrial planets including extrasolar super-Earths.

In recent disk evolution models (Bai et al. 2016; Suzuki et al. 2016), the effect of disk wind torque (magnetic braking) is also investigated. The disk winds carry away the angular momentum from the disk, which causes mass accretion. This is investigated in magneto-centrifugal disk wind model (e.g., Blandford & Payne 1982; Bai & Stone 2013; Simon et al. 2013; Gressel et al. 2015). Suzuki et al. (2010) considered the mass loss due to disk winds alone, and ignored the effect of disk wind torque. Suzuki et al. (2016) took this into account and calculated the global disk evolution. They found that the wind torque can alter the structure of gas surface density and flow, and in some cases, the radial mass accretion is dominated not by turbulence-driven accretion but by wind-driven accretion. Taking into account the wind-driven accretion, we can explain the observed mass accretion rate onto the central star. The observationally inferred mass accretion rate is larger than 10^-9M_⊙ yr^-1 at t = 10⁶ yr, which was not explained by MRI-inactive disks. However, it is shown that the accretion rate via magnetic braking may reproduce the observed rate (Suzuki et al. 2016) even in MRI-inactive disks.

As an advance on our previous papers (Ogihara et al. 2015a,b), here we incorporate two effects of disk winds. First, the evolution of gas surface density was calculated according to Suzuki et al. (2016), which considered the effect of wind torque. Second, effects of global gas flows on type I migration were considered. Global gas flows can be induced, for example, by wind-driven accretion, which may alter the picture of type I migration. According to recent studies of type I migration, the direction and rate of migration depend on saturation level of the corotation torque (e.g., Paardekooper et al. 2011). When there exists an accretion flow, differences in radial velocity between planets and disk gas would change the contribution of corotation torque (e.g., Masset & Papaloizou 2003; Paardekooper 2014). If there are high-velocity gas flows, the corotation torque may change. However, previous models did not examine these effects.

In this paper, we have investigated the effects of global gas flows on type I migration. In particular, we examined the case in which gas flow structures are developed by wind-driven accretion. We note that one can use the same recipe for more general gas flows. For example, the accretion flows would be vertically non-uniform in a protoplanetary disk with an MRI-dead zone, namely a layered accretion takes place (e.g., Gammie 1996; Turner & Sano 2008; Flock et al. 2015). In doing so, we kept the problem as simple as possible and develop a general analysis. Computationally expensive and challenging hydrodynamic simulations are not performed in this work. The plan of the paper is as follows. In Sect. 2 we describe the evolution of protoplanetary disks under the influence of wind-driven accretion; in Sect. 3 we examine type I migration using dynamical torque model; in Sect. 4.1 we discuss effects of global gas flows on type I migration using a different approach (namely, a viscous diffusion approximation); in Sect. 4.2 we discuss the outcomes of our simulations and give our conclusions.

Fig. 1 Time evolution of gas surface density of MRI-inactive disks (, panel a) and MRI-active disks (, panel b). The black lines represent disk evolution without wind-driven accretion (wind torque), while red lines indicate that with wind-driven accretion.

2. Effects of wind-driven accretion on mass accretion at midplane

In our previous papers (Suzuki et al. 2010; Ogihara et al. 2015a,b), the evolution of protoplanetary disks is investigated with effects of mass loss due to disk winds without considering the wind torque (wind-driven accretion). In Fig. 1 we show the gas surface density evolution of MRI-inactive disks () and MRI-active disks () without considering wind torque, where indicates an effective turbulent viscosity, which is commonly referred to as the α parameter based on Shakura & Sunyaev (1973). The disk evolution is the same as those shown in Figs. 1 and 5 in Suzuki et al. (2016). Their initial disk profile is set as a power-law distribution with a power-law index of –3/2 (Eq. (29) in Suzuki et al. 2016), in which an exponential cutoff at r = 30 au is considered. We note that we are interested in orbital migration of planets in the type I regime and consider the stage when planets have grown to the isolation mass (e.g., Kokubo & Ida 1998). Taking into account the growth time of isolation-mass planets, we refer to t = 10⁵ yr in Suzuki et al. (2016) as initial (t = 0). We see that the surface density exhibits shallow negative slopes or even flat in inner region. We also found that type I migration can be suppressed in such region (e.g., Ogihara et al. 2015a,b).

Suzuki et al. (2016) developed a one-dimensional model of protoplanetary disks that evolve through magnetically driven disk winds. They also considered effects of the wind torque, in which wind carries away some angular momentum from the disk (also known as magnetic braking). This causes the radial mass accretion. Evolution of gas surface density including the effects of the wind torque is also shown in Fig. 1. We also refer the reader to Sect. 2 of Suzuki et al. (2016) for details. The surface density slope is significantly different from the initial index of –3/2 in the close-in region (r ≲ 1 au) both for active and inactive disks. The slope can even be positive. Thus we can say that the wind torque increases the surface density slope. In MRI-inactive disks with smaller , the positive surface density slope is more likely to be realized. In the region where the mass transport due to disk winds dominates over that due to viscous transport, the surface density slope tends to be positive values.

The surface density profile is altered only in the inner region, because the mass transport and mass loss due to disk winds are stronger in regions closer to the star. This is due to the fact that disk evolution time is proportional to the dynamical time at each location r, that is r^{− 3/2}, and shorter for the inner region.

Figure 2 shows the temperature evolution, which is the same as that in Figs. 1 and 5 of Suzuki et al. (2016). A detailed description of the temperature model is given in Sect. 2.4 of Suzuki et al. (2016), but here both viscous heating and stellar irradiation are considered. When the wind torque is not taken into account, the gas surface density is high in the inner region during the early phase of disk evolution. Thus, the disk temperature is determined by the viscous heating at r < 10 au and t < 0.1 Myr, and the temperature slope is different from that of the radiative equilibrium for the minimum mass solar nebula (Hayashi 1981). In the late phase of disk evolution without the wind torque and in almost all phases for a model with the wind torque, the temperature slope is considered to be −1/2.

Figure 3 indicates the evolution of , which gives a measure of the angular momentum loss due to the wind torque. This is calculated according to Bai (2013) and Eq. (30) in Suzuki et al. (2016). That is, (1)where Σ_g,int is the initial gas surface density given by Eq. (29) in Suzuki et al. (2016). We see in Fig. 3 that decreases with increasing r, which means that the global gas flow can be rapid in the inner disk¹. We used Eq. (1) up to by extending the original expression introduced by Bai (2013) for .

Fig. 2 Time evolution of temperature of MRI-inactive disks (, panel a) and MRI-active disks (, panel b). The black lines represent disk evolution without wind-driven accretion (wind torque), while red lines indicate that with wind-driven accretion.

Fig. 3 Time evolution of of MRI-inactive disks (, red) and MRI-active disks (, blue).

We cautiously note that the slopes of the surface densities that were obtained based on our simple model may be affected by various effects that are not considered here. First, we adopted the zero-torque-gradient condition, , at the inner boundary, which forces the slope to be set to Σ_g ∝ r^{− 3/2} in the inner region. However, if a magnetospheric cavity is formed around the central star, the condition, Σ_g = 0, would be more realistic for the inner boundary, as is adopted by, for example, Matter et al. (2016). In such circumstances, the positive slope of Σ_g would be enhanced in the inner region (Eq. (25) and Fig. 4 of Lynden-Bell & Pringle 1974).

Second, we assumed a spatially constant , and that depends only on Σ_g (Eq. (1)). If the ionization effect by the radiation heating from the central star is included, both and are supposed to increase with r for a “usual” negative slope of Σ_g because the radiation can penetrate toward the midplane for smaller Σ_g, which is considered in the original expression of derived by Bai (2013). If the radial dependence of , and is taken into account, the obtained positive slope of Σ_g would be suppressed at early time when the negative slope of Σ_g is still preserved. However, once the gas in the inner region is cleared out, we expect that the positive dependence of and on r will be suppressed or even reversed.

There is an important caveat in the one-dimensional disk model; the vertical dependence of mass accretion is not known. Although the wind torque should operate at the disk surface, we do not know whether it affects the mass accretion at the midplane or not. We considered two cases. The first case is that the wind torque affects only the disk near the surface. At the midplane, the disk accretes solely through the turbulent viscosity (v_r = v_r,turb). Hereafter we refer this situation as “case A”. The second case is that the wind torque changes the accretion at the midplane and the mass accretion occurs over the entire region between midplane and surface (v_r = v_r,turb + v_r,wind). Hereafter we refer this situation as “case B”. Case A provides a rough lower limit on the mass accretion at the midplane, while case B corresponds to an upper limit. Figure 4 shows a schematic picture of the vertical dependence of mass accretion for each case. The reality would lie between these two cases. To determine which cases represent the reality, detailed three-dimensional MHD simulations are required. Recently, three-dimensional global MHD simulations have been performed by various groups. Simulations with a net vertical magnetic flux show the accretion takes place near the surfaces, which can be modeled by our case A, even without dead zones (Suzuki & Inutsuka 2014; Zhu & Stone 2017). On the other hand, simulations without a net vertical flux show the opposite radial flow pattern (Flock et al. 2011; Fromang et al. 2011). Vertical streams such as meridional circulation further perturb the radial flows and make the flow pattern more complicated (Béthune et al. 2017). At the moment, the situation is far from conclusive. Instead of taking into account the effect of complicated flow structure, in this paper we consider each end-member case and examine planet migration in that framework.

Fig. 4 Schematic illustration of the vertical distribution of radial mass accretion for case A and case B. We note that the vertically integrated accretion rate is the same for each case. Thus, the surface radial velocity in case A is generally much higher than throughout the disk’s height in case B.

The effect of global gas flows on type I migration have not been examined. Although we especially focused on the wind-driven accretion in this study, our model can be applied to more general flows. For example, case A resembles a case with layered accretion. Again, we needed hydrodynamical simulations to assess this effect. As a first step, in this study, we adopted a simple recipe. We considered that streamlines are altered by the mass accretion onto the central star, which causes an imbalance of corotation torques, and investigate type I migration. We note that we considered type I migration of planets with small inclinations (≪1) and focus on motion near the midplane.

3. Dynamical torque model

3.1. Model

We treat the gas drift as the feedback effect of the relative motion of the gas and planet. Masset & Papaloizou (2003) considered the relative motion between the gas and planet and found that a positive feedback exists as the corotation torque is proportional to the relative velocity, possibly leading to a runaway migration (type III migration). They suggested that the runaway migration happens for intermediate-mass planets. Paardekooper (2014) proposed that a similar feedback effect can occur for low-mass planets, which they call dynamical corotation torque, and developed a formula of the torque. The type I migration torque is expressed by superposition of the static torque and dynamical corotation torque: (2)where the static torque is the commonly used type I migration torque. For detailed description of the torque formulae, see Eqs. (50)–(53) in Paardekooper et al. (2011). The dynamical torque is given by (3)where x_s and Ω are the half-width of the horseshoe region and the Keplerian frequency, respectively, and v_rel = v_r,planet−v_r,gas is the relative velocity between the gas and planet. The quantity w_c indicates the inverse of the specific vorticity at the initial location of the planet, while w(r) is the inverse of the specific vorticity at r. The specific vorticity is defined as the vertical component of the flow vorticity divided by the gas surface density: (∇ × v)_z/ Σ_g. As stated in Paardekooper (2014), the factor 1−w_c/w(r) plays the same role as the coorbital mass deficit in classical type III migration. This factor can be calculated by (Paardekooper 2014) (4)where p = (dlnΣ_g/ dlnr) and ν represents the viscosity.

Using this approach, special care must be taken. First, Paardekooper (2014) developed the dynamical torque formula assuming that the drift timescale across the horseshoe region is longer than the libration time (x_s/v_rel >τ_lib). The libration time is longer for smaller planets, so that this condition is not satisfied for such planets as will be seen in Sect. 3.2. When this condition is violated, it is suggested that planets undergo runaway migration (Masset & Papaloizou 2003; Paardekooper 2014; Pierens 2015). However, it is uncertain whether such low-mass planets (≲1 M_⊕) can undergo runaway migration. Instead, a cut-off of the corotation torque may be expected when the radial drift time across the horseshoe region is shorter than the libration time (Masset 2001, 2002). Thus we set Γ_dynamic → 0 for x_s/v_rel <τ_lib. Second, the dynamical torque formula (Eq. (3)) is derived for locally isothermal disks. So we do not know the actual dynamical torque in radiative disks. In this study, the barotropic part of the corotation torque has a large value; thus the torque formula can be a good approximation. However, this would have to be taken into account in future work.

The relative velocity between the gas and planet v_rel, which was used in the calculation of the dynamical corotation torque, for cases A and B are indicated below.

3.1.1. Case A: no wind-driven accretion at midplane

In case A, the velocity of gas flow is determined by the turbulence-driven accretion: (5)where H is the disk scale hight.

3.1.2. Case B: wind-driven accretion at midplane

In case B, the radial gas velocity at the midplane is expressed by the sum of turbulence-driven accretion and wind-driven accretion (v_r,gas = v_r,turb + v_r,wind). The velocity due to the turbulence-driven accretion is given by Eq. (5). The radial velocity due to the wind torque is (Eq. (34) in Suzuki et al. 2016): (6)where c_s indicates the sound speed.

3.2. Migration map

Figure 5 shows the migration map when the midplane gas flow is solely determined by (case A; Sect. 3.1.1), while Fig. 6 indicates the map including the wind-driven accretion at the midplane (case B; Sect. 3.1.2). In Figs. 5a and 6a, is assumed (MRI inactive), while in Figs. 5b and 6b, is used (MRI active). Here we assumed v_r,planet = 0 when drawing the map. The color indicates the migration timescale with inward migration, outward migration, and no migration regions all shown. As described in Eq. (2), the type I migration torque consists of the static torque and dynamical corotation torque. The static type I migration torques is expressed by the superposition of Lindblad and corotation torques. The Lindblad torque is basically negative, but the corotation torque can be positive depending on the vortensity gradient and the entropy gradient of disk gas in the horseshoe region. As shown in Eqs. (4)–(7) in Paardekooper et al. (2011), the barotropic part of corotation torque depends on the surface density slope, while the entropy-related corotation torque depends on the temperature slope and the surface density slope. In general, when the surface density slope is positive and the temperature slope is negative, the gas in the horseshoe region provides a positive torque, which can overwhelm the negative Lindblad torque. Panels for t = 100 yr indicates the initial state. Table 1 shows the summary of the model.

Fig. 5 Migration maps of MRI-inactive disks () and MRI-active disks () based on case A prescription. The color indicates the migration timescale. In the blue region, planet migration is inward, while in the red region planets undergo outward migration. Snapshots of systems of N-body simulations are also shown by open circles (Σ0 = 30) and filled circles (Σ0 = 75), respectively.

For comparison, we plot the migration map for cases with mass loss due to disk winds but without wind-driven accretion due to the wind torque in Fig. 7. The migration map shows that inward migration is expected in almost the whole region. The inward migration timescale is longer than in the case with power-law density slope of –3/2 (see also Ogihara et al. 2015a,b). In the inner region (r ≲ 1 au), the surface density profile is shallower than −3/2 and the barotropic part of the corotation torque partially compensate for the Lindblad torque. In addition, Fig. 2 shows that the temperature profile has a steep slope near r ~ 1 au for t < 0.1 Myr in the case without wind torque, especially for inactive disks. The entropy-related corotation torque has a positive value near r ~ 1 au, which also partially compensates for the Lindblad torque. The total torque is still negative, leading to inward migration. Comparing Figs. 5 and 7, effects of the wind torque on surface density slope (and hence on type I migration) are seen. When the wind torque is considered, the surface density is carved out from inside (see Fig. 1), so that the corotation torque compensates for the Lindblad torque and the disk dissipation timescale is shortened in the inner region, leading to a further suppression of inward migration. That is, the inward migration region in Fig. 7 changes to no migration or outward migration region in Fig. 5. The temperature profile is altered for t < 0.1 Myr in the case without wind torque, which gives positive entropy-related corotation torque. In the case with wind torque, however, the temperature profile can be considered as a power-law distribution with index −1/2. Although the temperature profile is different, the main differences in migration map (see outward migration region in Fig. 5) can be attributed to the difference in the surface density slope.

Comparing Figs.5a and b, we see that the region of outward migration is wider in the first case (MRI inactive disk). This is because, as seen in Fig. 1, the region where the surface density slope is positive is wider in MRI-inactive disks than in MRI-active disks. As stated above, the positive surface density slope gives positive corotation torque. Thus the outward migration region is wider in Figs. 5a and 6a than in Figs. 5b and 6b. As seen in Fig. 2, the temperature distribution can be considered as a power-law distribution, so that the temperature slope does not play a major role in changing the corotation torque. We note that in MRI active disks, the surface density slope is almost flat between r = 0.1–1 au (see Fig. 1b), which results in substantial suppression of both inward and outward migration. We note also that in disks with no effects of disk winds, the outward migration region moves toward the lower-left region (e.g., Bitsch et al. 2015). However, we see in migration maps that the outward migration region moves toward the right of the figure. This is because the region with a positive surface density slope expands to the outer region as time increases. We also note that in case A, the direction and rate of type I migration are almost solely determined by the static torque, because the radial gas velocity at the midplane is not large enough to change the migration map.

Figure 6 shows the migration map, in which wind-driven accretion at the midplane is considered (case B). That is, v_rel = v_r,planet−(v_r,turb + v_r,wind). We see that the outward migration region extends toward the upper-left region (smaller r and larger M). As seen in Fig. 3, the radial velocity due to wind-driven accretion is large in inner region, leading to a high v_rel. The magnitude of dynamical corotation torque is proportional to , thus the torque is large in inner region. In the lower-left corner of the panels of Fig. 6a, the condition x_s/v_rel >τ_lib is not satisfied and Γ_dynamic is set to zero. Another important point concerns the sign of the dynamical torque. As discussed in Paardekooper (2014), the sign of the dynamical torque depends on the background vortensity gradient. If the surface density slope is positive, the dynamical torque has a positive value (and fully unsaturated corotation torque tends to be positive). We see no significant change in Fig. 6b from Fig. 5b before t = 1 Myr. This is because the disk evolution is dominated by the turbulence-driven accretion in MRI-active disks ().

Fig. 6 Same as Fig. 5 but based on case B prescription.

Fig. 7 Migration map of MRI-inactive disks () and MRI-active disks (), in which effects of wind torque are not considered for calculation of gas surface density eovlution. Corresponding disk evolution is shown in Fig. 1 by black lines.

Fig. 8 Orbital evolution of planets for Σ0 = 75 in MRI-inactive disks. Snapshots of the system for panels a and b are shown by filled circles in Fig. 5a and Fig. 6a, respectively. The filled circles connected with solid lines represent the sizes of planets. The color of the lines indicates the eccentricity (color bar).

3.3. Orbital evolution

In order to demonstrate the effects on planet formation, we perform N-body simulations of orbital evolution of isolation-mass planets. Planets undergo the type I migration and the tidal damping of eccentricities and inclinations (see Eq. (6) of Ogihara et al. 2015a, for e-damping timescale). We note that, although v_r,planet = 0 is used when drawing migration maps in Figs. 5 and 6, the actual radial velocity of the planet is used in calculating the dynamical torque. For numerical integration, we use a fourth-order Hermite scheme with a hierarchical individual time step. When the physical radii of two spherical bodies overlap, they are assumed to merge into one body, assuming perfect accretion. The initial planetary mass (M_iso = 2πaΔaΣ_d) is calculated according to power-law solid distributions: (7)where Σ₀ = 30 or 75 is used. Then (8)The initial separation between oligarchs is set by Δa = 10r_H, where r_H is the mutual Hill radius. Snapshots of each system are overplotted in Figs. 5 and 6. Open circles represent the case in which Σ₀ = 30, while solid circles represent the case for Σ₀ = 75.

As seen in Figs. 5b and 6b, results are similar between cases A and B in MRI-active disks because migration maps are similar. The surface density slope is almost flat in the inner region (Fig. 1b) and migration is significantly suppressed in MRI-active case.

From Figs. 5a and 6a, we see that orbital evolution and final orbital configurations differ in each model. In case B, no planets migrate inside r = 0.08 au (Fig. 6a). In case A, planets of a few Earth masses pass over the outward migration zone (red region in Fig. 5a). This does not occur in case B, because the outward migration zone extends to high-mass regions (red region in Fig. 6a), which prevents the inward migration of super-Earth mass planets. Thus wind-driven accretion has an effect of preventing inward migration. In addition, in Fig. 6a, super-Earth mass planets undergo outward migration to outer region (r > 5 au) (filled circles). For more detailed time evolution of semimajor axis for Σ₀ = 75 in MRI-inactive disks, see Fig. 8a (case A) and b (case B). As stated above, no inward migration is observed in Fig. 8b.

According to Paardekooper (2014), planets may undergo runaway migration when the k-coefficient defined by Eq. (31) of Paardekooper (2014) is larger than 0.5. This condition for runaway migration is satisfied only for r < 0.1 au and M ≳ 5 M_⊕ at t = 0.1 Myr in Fig. 6a. In our N-body simulations, no planets satisfy this condition.

4. Discussion and conclusions

4.1. A different approach: diffusion model

Before we make concluding remarks, we briefly see effects of global gas flows on migration based on a different approach. According to Paardekooper et al. (2011), the corotation torque for type I migration is expressed by using the following parameters: where , ν and χ are dimensionless half-width of the horseshoe region, viscous and thermal coefficients, respectively. Values of effective viscosity ν is presented below. For thermal diffusion coefficient χ, we use the same form as Eq. (34) of Paardekooper et al. (2011), in which a coefficient is corrected as described in Bitsch & Kley (2011). The corotation torque has a maximum for P_ν ≃ 0.3 (see Fig. 3 in Paardekooper et al. 2011). A simple physical interpretation for the dependence on diffusivities is as follows. The corotation torque comes about because of the gas in the planet’s horseshoe region making inward and outward U-turns relative to the planet at different vortensities and specific entropies. Both these quantities evolve in time because of their advection-diffusion-creation within the horseshoe region, and the magnitude of the corotation torque thus depends on how the timescales for viscous and thermal diffusions across the horseshoe region compare with the horseshoe U-turn and libration timescales (Baruteau et al. 2014).

We calculated an effective “viscosity” from radial accretion and examine the consequent de-saturation of corotation torque. To incorporate effects of gas flows due to wind-driven accretion in the formula of type I migration, we derived an “effective” viscosity from the radial gas flow.

4.1.1. Case A: no wind-driven accretion at midplane

In case A, in which the viscosity at the midplane is solely expressed by , the diffusion at the midplane is solely due to MRI turbulence (ν = ν_turb): (11)Then we calculated the corotation torque by using (12)In this case, we confirmed that the migration map is almost identical to that in Fig. 5 both for MRI-inactive and active disks.

4.1.2. Case B: wind-driven accretion at midplane

In case B, using an approximate relation between ν and v_r of (13)the effective “viscous” coefficient due to wind-driven accretion can be expressed as (14)where Eq. (6) is used for v_r,wind. Then we calculated the corotation torque using (15)Figure 9 shows the migration map for case B. Regarding the effect of wind-driven accretion on desaturation of corotation torque, the outward migration region is more extended towards the upper-left region in Fig. 9a. This is because the effective viscosity is larger in closer region due to high . The corotation torque has a maximum for P_ν ≃ 0.3, and this corresponds to planets with M ≃ 6 M_⊕ (t = 100 yr) and M ≃ 30 M_⊕ (t = 1 Myr) at r = 0.1 au in Fig. 9a. However, for MRI-active disks, we see no significant differences between Figs. 5b and 9b before t = 1 Myr. This is because in MRI-active disks, is large and ν is determined primary by except for the very late stage (t > 1 Myr).

We see that the migration map of Fig. 6a and 9a show similar trends. The similarity might be explained by the fact that in the region of the disk with a positive surface density gradient the de-saturated corotation torque is positive (Masset et al. 2006); likewise, the dynamical torque is positive, because the gas flows from a high vortensity region to a low vortensity region (Paardekooper 2014).

One concern may be whether the gas flow promoted by the wind torque can be interpreted as a flow of the angular momentum into the corotation region, as presented in Sect. 4.1. We think hydrodynamical simulations are required to determine if this approach is valid or not. However, the rapid gas flow at the surface due to wind-driven accretion induces a large velocity difference between the surface and the midplane, which would cause some diffusion processes (e.g., turbulence). In this case, our viscosity approach may represent reality. We note that it is not clear how likely it is that this turbulence can be interpreted as a source of diffusion, because this depends on the statistical properties of turbulence (e.g., Baruteau & Lin 2010). This should be investigated in future studies.

Fig. 9 Same as Fig. 6 but based on a diffusion approach.

4.2. Conclusions

We have examined effects of wind-driven accretion on type I migration. As seen in Suzuki et al. (2016), the disk wind torque has an effect of changing the surface density slope. The surface density slope can be carved out from inside. In addition, we considered effects of gas flows on type I migration and found that gas flows at the midplane has an effect of preventing the infall of super-Earths onto the central star. In MRI-inactive disks, the surface density slope can be positive due to the wind torque and an outward migration region appears in the migration map. The outward migration region extends for case B, in which the gas flow at the midplane is induced due to wind-driven accretion. This prevents the inward migration of massive super-Earths in the close-in region. We simulated orbital evolution of planets by N-body simulations and observed that no planets fall onto the star in this case. We also see super-Earth mass planets migrate outward. In the MRI-active disks, we usually obtained a disk with flat slope in inner region due to the effect of wind torque, which results in substantial suppression of both inward and outward migration. In MRI-active case, the effect of gas flows due to wind-driven accretion on desaturation of the corotation torque is not important, because the turbulence-driven accretion primarily controls the surface density slope and corotation torque.

We note that although similar trends are observed in migration maps using different approaches, they are not perfectly matched with each other. One may also wonder that if planets open a density gap, the formulae of type I migration would not be valid and planets would migrate in the type II regime². So, we checked the gap-opening criterion of Crida et al. (2006) and found that even massive super-Earths do not open a gap. This is due to high mass accretion rate due to wind-driven accretion.

We note that the torque formula of Paardekooper et al. (2011) would be accurate only when (M/M_∗)/(H/r)³ < 1. For example, planets with M > 5 M_⊕ at r = 0.1 au and with M > 40 M_⊕ at r = 1 au do not satisfy this condition. Although no planets in our N-body simulations break this condition, the torque formula should be treated with some caution for more massive planets.

In this paper, we focused on wind-driven accretion; however one can apply the model to more general gas flows. For example, case A of MRI-inactive disks would correspond to a dead zone region with an active surface layer.

As another application of this study, it is interesting to investigate how planets with high inclination undergo migration. The relative velocity between gas and planet changes with z, which may change migration. Special care would have to be taken in this case because the standard formula for type I migration is derived for i <H/r.

As Suzuki et al. (2016) cautioned, the disk evolution with effects of magnetically-driven disk winds depends on the evolution of the net vertical magnetic field. It is very important to reveal the evolution of the vertical magnetic field by MHD simulations. From another perspective, effects of the magnetic filed on the corotation torque have also been investigated (e.g., Guilet et al. 2013; Uribe et al. 2015), which also depends on the evolution of the magnetic field.

It is very important to address two issues by performing hydrodynamical simulations in future work. First, we need to determine the vertical dependence of mass accretion when the wind-driven accretion is considered. Second, we need to examine how planets undergo type I migration when the radial mass accretion is considered. We would encourage (M)HD simulations that explore these issues.

Acknowledgments

We thank the anonymous referee for valuable comments. We also thank Bertram Bitsch, Hidekazu Tanaka and Takanori Sasaki for useful comments and suggestions. Numerical computations were conducted on the general-purpose PC farm at the Center for Computational Astrophysics, CfCA, of the National Astronomical Observatory of Japan. This work was supported by JSPS KAKENHI Grant Number 16H07415.

References

All Tables

All Figures

	Fig. 1 Time evolution of gas surface density of MRI-inactive disks (, panel a) and MRI-active disks (, panel b). The black lines represent disk evolution without wind-driven accretion (wind torque), while red lines indicate that with wind-driven accretion.
In the text

	Fig. 2 Time evolution of temperature of MRI-inactive disks (, panel a) and MRI-active disks (, panel b). The black lines represent disk evolution without wind-driven accretion (wind torque), while red lines indicate that with wind-driven accretion.
In the text

	Fig. 3 Time evolution of of MRI-inactive disks (, red) and MRI-active disks (, blue).
In the text

	Fig. 4 Schematic illustration of the vertical distribution of radial mass accretion for case A and case B. We note that the vertically integrated accretion rate is the same for each case. Thus, the surface radial velocity in case A is generally much higher than throughout the disk’s height in case B.
In the text

	Fig. 5 Migration maps of MRI-inactive disks () and MRI-active disks () based on case A prescription. The color indicates the migration timescale. In the blue region, planet migration is inward, while in the red region planets undergo outward migration. Snapshots of systems of N-body simulations are also shown by open circles (Σ0 = 30) and filled circles (Σ0 = 75), respectively.
In the text

	Fig. 6 Same as Fig. 5 but based on case B prescription.
In the text

	Fig. 7 Migration map of MRI-inactive disks () and MRI-active disks (), in which effects of wind torque are not considered for calculation of gas surface density eovlution. Corresponding disk evolution is shown in Fig. 1 by black lines.
In the text

	Fig. 8 Orbital evolution of planets for Σ0 = 75 in MRI-inactive disks. Snapshots of the system for panels a and b are shown by filled circles in Fig. 5a and Fig. 6a, respectively. The filled circles connected with solid lines represent the sizes of planets. The color of the lines indicates the eccentricity (color bar).
In the text

	Fig. 9 Same as Fig. 6 but based on a diffusion approach.
In the text