Structure Formation in a Young Protoplanetary Disk by a Magnetic Disk Wind
Structure formation in young protoplanetary disks is investigated using a one-dimensional model including the formation and the evolution of disks. Recent observations with ALMA found that a ring-hole structure may be formed in young protoplanetary disks, even when the disk is embedded in the envelope. We present a one-dimensional model for the formation of a protoplanetary disk from a molecular cloud core and its subsequent long-term evolution within a single framework. Such long-term evolution has not been explored by numerical simulations due to the limitation of computational power. In our model, we calculate the time evolution of the surface density of the gas and the dust with the wind mass loss and the radial drift of the dust in the disk. We find that the MHD disk wind is a viable mechanism for the formation of ring-hole structure in young disks. We perform a parameter study of our model and derive condition of the formation of ring-hole structures within years after the start of the collapse of the molecular cloud core. The final outcome of the disk shows five types of morphology and this can be understood by comparing the timescale of the viscous diffusion, the mass loss by MHD disk wind and the radial drift of the dust. We discuss the implication of the model for the WL 17 system, which is suspected to be an embedded, yet transitional, disk.
Subject headings:protoplanetary disks, stars: formation
Protoplanetary disks are thought to be the birth place of planets. The density and temperature structures strongly affect the formation processes of planets. Recent observations reveal the detailed structures of protoplanetary disk. They have found that spiral structure (e.g. Muto et al. 2012; Grady et al. 2013; Benisty et al. 2015; Pérez et al. 2016), non-axisymmetric structure (e.g. van der Marel et al. 2013b; Casassus et al. 2013; Fukagawa et al. 2013; Muto et al. 2015), and ring-like structure (e.g. Geers et al. 2007; Isella et al. 2010, 2012, 2013, 2016; Andrews et al. 2011; Hashimoto et al. 2011, 2012; Mathews et al. 2012; Mayama et al. 2012; ALMA Partnership et al. 2015; Andrews et al. 2016; Tsukagoshi et al. 2016; Fedele et al. 2017a, b; Loomis et al. 2017) are formed in protoplanetary disks.
The observed disk structures will give us some clues to reveal the planet formation scenario. Especially, the disk structure formed in the early evolutionary phase is important to understand the disk evolution process and the initial condition of the planet formation. Recently, Sheehan & Eisner (2017) have found that the young protoplanetary disk around WL 17 in Ophiuchus star forming region exhibits a ring structure. The dust continuum emission from the disk shows a ring-like structure with a central hole. This is similar to the structure of transitional disks. However, one remarkable difference between standard transitional disks and the WL 17 system is the disk age. Observations of Sheehan & Eisner (2017) suggest that WL 17 is still covered by the envelope. It is regarded as class I YSO van Kempen et al. (2009); Enoch et al. (2009), whose age is suggested to be Myr Evans et al. (2009), despite large uncertainty. WL 17 appears like transitional disks, which are one or two orders of magnitude older. The ring-hole structure formation mechanism has not been well studied for such young disks. One possibility is the gap formation by (an) unseen planet(s), but it may be very difficult to form planets at such an early stage. The observations of WL 17 suggest that there may be another mechanism that results in the ring structure formation at the early evolutionary phase.
In this work, we investigate the early phase of disk formation and evolution to explore the possibility of forming small scale structures in young disks. The important processes for the early evolution of the disk are the gravitational collapse of the cloud core, the angular momentum transfer due to the gravitational instability within the disk, the growth and the radial drift of dust particles, and the effect of the magnetic fields.
The disk formation through the gravitational collapse of the cloud core and the angular momentum transfer due to the gravitational instability are investigated well by using three-dimensional numerical simulations (e.g. Bate 1998; Machida et al. 2007; Tomida et al. 2010; Inutsuka et al. 2010; Tsukamoto & Machida 2011). The subsequent long-term evolution of the star-disk systems have been investigated separately from the star formation phase by assuming some initial disk model. This is partly due to the limitation of computational power of solving all the way from the initial star and disk formation from the molecular cloud core to the final dispersal. However, the evolution of young disks that we focus on in this work may strongly depend on the final outcome of the star-disk formation, which is the initial condition of the disk evolution.
The pioneering theoretical work on the formation and evolution of protoplanetary disks is done by Cassen & Moosman (1981) (see also Cassen & Summers 1983). They have developed the model for the evolution of a one-dimensional viscous accretion disk including the effect of the infall from the cloud core onto the disk. Following their work, Stahler et al. (1994) have done the detailed investigation on the trajectory of accreting gas onto the disk in the early stage of the star-disk formation. Hueso & Guillot (2005) have performed the model calculations with large parameter space of the viscosity, temperature and rotation rate of the cloud cores and compared the results with the observed disk structures. The model calculations show that the disk becomes gravitationally unstable in their formation phase Nakamoto & Nakagawa (1994). Zhu et al. (2010) have investigated the disk evolution due to the angular momentum transfer caused by the gravitational instability and MRI using two-layered disk model. These previous studies used the mass accretion rate from the cloud core onto the disk obtained from the self-similar solution of singular isothermal sphere Shu (1977), which gives time-independent mass accretion rate. In Takahashi et al. (2013), we have constructed the analytical model providing the time-dependent mass accretion rate from the cloud core with arbitrary radial density profile.
The importance of dust growth and radial drift at the initial disk formation stage has been pointed out recently. Two-dimensional numerical simulations show that dust particles grow and drift inward resulting in the small dust-to-gas mass ratio during the disk formation Vorobyov et al. (2018). The dust growth and the decrease of dust-to-gas ratio in the disk formation stage are also indicated by the steady accretion disk model Tsukamoto et al. (2017). Such behavior of dust particles in the disk formation phase will affect the following disk evolution. To investigate the early evolution of disks, we need to deal with the gravitational collapse of the cloud core, the disk evolution caused by the angular momentum transfer, and the dust evolution in the disk comprehensively.
In addition, we investigate the effects of the disk wind that can be driven by the magneto-rotational instability Suzuki & Inutsuka (2009); Suzuki et al. (2010) on the formation of small scale structures in young disks. Suzuki et al. (2010) investigated the mass loss rate due to the disk wind. They parameterized the mass loss by , where . According to the three-dimensional local MHD simulations, . The timescale of wind mass loss scales with the local Kepler time of the disk, which is faster at inner radii. Therefore, if is constant in the disk, the wind mass loss is efficient in the inner region and the inner hole is naturally formed by the wind. In such disks, it is expected that the dust concentrates around the inner edge of the disk and the ring-hole structure in dust distribution will be formed. The long-term evolution of the disks including the disk wind is investigated by using one-dimensional models treating isolated disks that have already been formed (e.g. Suzuki et al. 2016; Pinilla et al. 2016). These work, however, did not calculate the disk formation phase.
In this paper, we calculate the formation and the evolution of protoplanetary disks within a single framework. We extend the model provided in Takahashi et al. (2013), which takes into account the time-dependent mass accretion, to include the effect of the dust and the wind mass loss on the disk evolution. We show that various gas and dust distributions can be formed by the MRI disk wind in young protopranetary disks.
In this work, we calculate the formation and the evolution of protoplanetary disks to investigate the ring structure formation in young disks. Since we focus on young disks whose ages are similar to the disk formation timescale ( yr), we cannot assume the already formed protoplanetary disks, for example minimum mass solar nebulae, as initial stage. Instead, we adopt molecular cloud cores as initial conditions and calculate the gravitational collapse of the cloud core and the disk evolution simultaneously.
Figure 1 shows the schematic picture of the formation and evolution of protoplanetary disks. We calculate the gravitational collapse of the cloud core by a one-dimensional semi-analytic model (Panel (a) of Figure 1). The gas around the center of the core has small angular momentum so that it makes the protostar through the gravitational collapse. The gas in the outer region of the cloud core has large angular momentum. Thus, it cannot directly falls onto the protostar but makes the disk around the protostar. Here, we assume that the gas falls at the centrifugal radius of the disk. Using the semi-analytic model, we obtain the mass infall rate per unit area on the protoplanetary disk Takahashi et al. (2013). We also calculate the evolution of the disk with the mass infall from the cloud core (Panel (b) of Figure 1). After the end of the collapse of the cloud core, we take into account the wind mass loss from the disk (Panel (c) of Figure 1). In this section, we show the equations to calculate the gravitational collapse of the cloud core (Panel (a) of Figure 1) and the evolution of the protoplanetary disks (Panel (b) and (c) of Figure 1).
2.1.1 Gravitational collapse of cloud cores and gas infall onto disks
First, we describe the one-dimensional semi-analytic model to calculate gravitational collapse of cloud cores (Panel (a) of Figure 1). Here, we use the models given in Takahashi et al. (2013) (see also Takahashi et al. (2016a) for a detailed derivation). The validity of the model is tested by the comparison with a three-dimensional numerical simulation in Takahashi et al. (2013, 2016a) and we found that the model reproduces the results of the numerical simulation well. In this model, we first evaluate the total mass accretion rate onto the central star and the disk to obtain the increase rate of the disk surface density per unit time due to the mass infall from the cloud core. Since the typical scale of the cloud cores ( pc) is much larger than that of the disk ( au), we neglect the centrifugal force and assume the spherical symmetry in the modeling of the collapsing cloud core. We adopt the mass accretion rate onto the center of the cloud core as the total accretion rate onto the star and the disk. We assume that molecular cloud cores and infalling envelopes are isothermal (10 K). The initial density profile is given so that the initial gravitational force is larger than the initial pressure gradient force by a factor of . Because of the pressure gradient force, the time when the collapsing gas reaches the disk is larger than the free-fall time (Equation (6) in Takahashi et al. (2016a)):
Since the free-fall time depends on the initial radius of the cloud core, the infall time also depends on the initial radius of the gas in the core. The mass accretion rate onto the disk is given by
where is the initial radius of the infalling gas and is the initial density of the cloud core at radius .
We assume that cloud cores initially rotate rigidly and the specific angular momentum of the infalling gas is conserved. To derive the increase of the surface density of the disk per unit time , we assume the envelope gas falls at the centrifugal radius (see also Equation (10) in Takahashi et al. 2013):
where is the initial angular velocity of the cloud core, is the disk radius, and is the specific angular momentum distribution of the disk, which is related with the disk angular velocity by . This equation is the same as Equation (6) in Hueso & Guillot (2005) if . In this work, we assume that the total molecular cloud core mass . We stop the infall when the all of gas in the cloud core falls onto the central star or the disk.
2.1.2 Disk evolution
Due to the gas infall from the core, protoplanetary disks are formed. In order to investigate the small structure formation in young disks, we calculate the time evolution of the surface densities of the gas and the dust (Panel (b) and (c) of Figure 1). We assume that the disk evolution timescale is much larger than the Kepler timescale. The radial motion of dust due to the interaction with the gas has been well studied (cf. Whipple 1972; Weidenschilling 1977; Birnstiel et al. 2016, and the references therein). The one-dimensional model of gas and dust evolution is obtained for both inviscid and viscous disks (cf. Nakagawa et al. 1986; Kretke et al. 2009; Kanagawa et al. 2017). In this work, we extend the model to the disk with the mass accretion from the cloud core and the wind mass loss. The equations for the surface densities and the radial velocities of the gas and the dust are given as follows (a derivation of the equations are given in Appendix),
where and are the surface density and the radial velocity of the gas, and are those of the dust, is the mass loss rate per unit area due to the disk wind, is the dust-to-gas mass ratio in the infalling envelope, is the angular velocity of gas and dust, is the enclosed mass within the radius . The specific torque is given by
where is the coefficient of the kinematic viscosity. Here, we assume that the viscosity is caused by the gravitational instability and magnetorotational instability (see Equation (15)). The factor represents the effect of the self-gravity of the disk of gas,
The Stokes number St is
where is the radius of the dust, is the internal density of the dust and is the gas density at the midplane of the disk. We assume . is the modified Stokes number,
We evaluate as follows,
where is the pressure at the disk midplane and is the isothermal sound speed with mean molecular weight . The mass input/loss rate due to the infall and disk wind within the radius is
We assume the centrifugal balance in the disk and that the gravitational force is approximately given by and the pressure gradient force, the frictional force and the viscosity are small compered with the centrifugal force and the gravitational force. Thus, the angular velocity distribution in the disk is given by
As shown in Equation (8), we evaluate the torque by using the effective viscosity represented by , where is a dimensionless measure of turbulent intensity Shakura & Sunyaev (1973). We assume the locally isothermal in the envelope and the disk, which are given in Section 2.3.
In this work, we take into account the angular momentum transport caused by gravitational instability and magnetorotational instability (MRI, Balbus & Hawley 1991),
The is given by a function of the Toomre’s parameter
This effective viscosity becomes efficient when the disk becomes gravitationally unstable Zhu et al. (2010). This model can mimic the disk evolution due to the angular momentum transfer by spiral arms formed in gravitationally unstable disks Takahashi et al. (2013). We assume that is constant in the disk.
We use the simple model to evaluate the mass loss rate due to the disk wind . We assume that the disk wind becomes efficient after the infall is finished. We discuss the validity of this treatment in Section 4.3. We introduce the efficiency parameter for the wind mass loss as follows
where is the mass that has already fallen onto the central star and the disk. Since the gas density decreases with increasing height and that at the launching point of the wind is a factor of smaller than that on the midplane, the coupling between the gas and the dust in the wind is much weaker than that on the midplane. The dust grains are blown out by the wind only when they are small enough to be coupled with the gas in the wind. Miyake et al. (2016) have shown that the dust grains with can be blown out by the disk wind. Since such grains are much smaller than that we adopt in this work (see Table 2), the dust mass loss rate is assumed to be zero in our model. 111The small dust blown out by the wind will be the origin of the continuum emission observed in the outflow(cf. Gueth et al. 2003). When we calculate the time evolution of the dust size distribution and wind mass loss of the dust, we can compare the model with the observations and test the validity of the model of the dust wind mass loss.
The list of symbols used for equations in this paper is shown in Table 1.
|time that the envelope reaches center, Equation (1)|
|free fall time of the envelope|
|mass infall rate into the disk, Equation (2)|
|initial radius of the infalling gas|
|initial density of the infalling gas|
|specific angular momentum|
|surface density of gas|
|radial velocity of gas|
|surface density of dust|
|radial velocity of dust|
|mass infall rate per unit area|
|wind mass loss rate per unit area|
|coefficient of the kinematic viscosity|
|enclosed mass within the radius|
|specific torque, Equation (8)|
|internal density of dust|
|gas density at the midplane of the disk|
|St’||modified Stokes number, Equation (11)|
|pressure at the disk midplane|
|isothermal sound speed|
|mass infall/loss rate within|
|Toomre’s parameter, Equation (17)|
|total infalled mass|
|temperature of the disk, Equation (19)|
|mass enhancement factor of the core|
|initial angular velocity of the core|
|total mass of the cloud core|
|dust-to-gas mass ratio in the envelope|
|St||Stokes number, Equation (10)|
|mean molecular weight|
|dimensionless measure of MRI turbulent intensity|
|efficiency parameter for wind mass loss|
|initial central density of the core|
|temperature of the core|
2.2. Boundary Conditions
In the calculations, we set the outer boundary at au and the inner boundary at 1 au. At the outer boundary, we assume . For the inner boundary, we assume at the center to evaluate . We also assume is constant at the inner boundary.
2.3. Model setup
Since the calculations start from gravitational collapse of cloud cores, we adopt the density and rotation profiles of cloud cores as initial conditions. The density distributions of cloud cores are given by the Bonner-Evert sphere with the central density is and the temperature is K. In order to make cores gravitationally unstable, we increase the density of the BE sphere by a factor of 1.4. The rotation velocity of the cores is given by the rigid rotation with the angular velocity 0.3 (the dependence of the disk evolution on is discussed in Section 4.1).
We assume that the temperature is given by the equilibrium temperature that is obtained from the balance between the irradiation heating from the central star and the radiation cooling at the disk surface Chiang & Goldreich (1997),
The parameters used in this work are shown in Table 2.
Due to dust growth in the disk, the typical dust size depends on the disk radius. In the inner region, the small dust grows quickly. The resulting dust size distribution is roughly given by the constant Stokes number of Okuzumi et al. (2012). We therefore assume that the Stokes number is the same at all radii of the disk and set the fiducial value of the Stokes number to be 0.1, but we explore the parameter space of St. The disk evolution in the case of the constant dust radius is discussed in Section 4.2. The smallest value of St= corresponds to the tight coupling between gas and dust resulting in almost constant dust-to-gas mass ratio in the disk. We set the largest Stokes number to be 10 times larger than the fiducial value: St=1. The largest value corresponds to the size for which the interaction between the gas and the dust is the most efficient. Note that the fiducial value of the Stokes number corresponds to the mm size grains in the ring structure, which is the size probed by mm- and submm- wavelengths observations (see Section 4.4). Since the strength of the turbulence and the wind mass loss is quite uncertain, we also explore the parameter space of and obtained in the numerical simulations on the MRI and disk wind (e.g. Suzuki et al. 2010; Okuzumi & Hirose 2011). The fiducial parameters that we involve in this work is summarized in Table 2.
3.1. Disk Evolution and Ring Formation
First of all, we present the result of a model with one set of parameters, where we observe the formation of a ring structure. The parameters are and St=0.1.
We first look at the mass infall phase where a protoplanetary disk is formed around a forming central star. Figure 2 shows the mass infall rate from the core onto the disk and the central star.
The infall rate does not depend on the parameters and St but depends only on the density and temperature profiles of the cloud core. The infalling gas reaches the center of the cloud core at Myr, which is the time of the protostar formation. The gas infall continues until Myr. The time evolution of the mass infall rate per unit area is shown in Figure 3.
In this case, the angular momentum of the cloud core is small and all the gas infalls within au. The size of the disk depends on the angular velocity of the core, which will be discussed in Section 4.1. The top panel of Figure 4 shows the time evolution of the surface density of the gas and dust before the gas infall from the cloud is finished at Myr. The disk is kept gravitationally unstable during this phase since gas is continuously supplied from the envelope. Thus, the angular momentum is redistributed mainly due to the effect of the gravitational instability (Equation (16)). The gas accretes inward and/or expands outward resulting in the decrease of the surface density while the infall from the cloud increases it. As a result, the disk with infall from cloud core sustains (cf. Lodato & Rice 2004; Vorobyov & Basu 2007; Tsukamoto et al. 2015b; Takahashi et al. 2016b). Using the central star mass and temperature distribution given in Equation (19), the radial distribution of the surface density with is
This surface density distribution is shown as a red dotted line in the top panel of Figure 4. Note that the gas disk radius is larger than the maximum centrifugal radius of the infalling gas. The radius of the disk expands due to the outward mass flux caused by the angular momentum transfer in the disk. On the other hand, dust disk radius is almost same as the maximum centrifugal radius of the infalling gas. This is because the dust decouples from the gas and does not expand outward due to the large Stokes number St=0.1. Moreover, the radial drift of the dust causes the smaller dust-to-gas mass ratio than that of the infalling envelope .
We now turn our attention after the mass infall from the cloud core is finished at Myr. We note that our model calculations continue from the infall phase, although we present the results before and after the infall stops separately. The bottom panel of Figure 4 shows the evolution of the surface density of the gas and the dust. The time Myr is just after the infall stops and the wind mass loss stars. In Figure 5, we show the time evolution of the mass loss rate.
The wind mass loss make the ring-hole structure in the gas disk. The dust accumulates at the radius of the pressure maximum around the inner edge of the gas ring structure. The dust ring structure at the radius of au with the width of au is formed at Myr, which corresponds to the protostar age (calculated from 0.1 Myr) of Myr.
Figure 6 shows the ring radius and dust-to-gas mass ratio at the center of the ring with , and St=0.1.
The ring radius increases with time because the radius of the pressure maximum increases due to the wind mass loss of gas. The dust-to-gas mass ratio is about 0.3 for Myr. For Myr, the dust-to-gas mass ratio increases because significant amount of gas is depleted due to the disk wind while the dust particles remain in the disk. The expansion of the ring slows down as dust particles condensate. To move the dust outward in the disk, the gas need to give the angular momentum to the dust. However, when the dust-to-gas mass ratio increases, the gas has no enough inertia to move the dust. Thus, the ring radius does not increase after the time when the dust-to-gas mass ratio become larger than unity. This behavior appears only when we use the equations taking into account the back reaction from dust to gas.
3.2. Parameter Study
The gas and dust structures formed in disks depend on the parameters , , and St. In this section, we show other disk structures formed in the disk. Figure 7 shows four typical outcomes of the disk structures obtained form our model calculations: the dust gap disk, the filled dust disk, the dust poor disk, and the normal filled disk are formed.
In the case of the dust gap disk (, , St=0.03), the dust ring is formed at au. However, the interior of the ring ( au) is filled by dust. Thus, this structure is observed as a gap structure in the dust emission. In the case of the filled dust disk (, , St=), the surface density of the dust has no remarkable structure. On the other hand, the gas is quickly removed by the wind mass loss resulting in no or quite small amount of gas in the disk. In contrast to this case, there is no structure in the gas surface density and only few amount of dust remains in the case of the dust poor disk (, , St=1). This small dust surface density is due to the efficient radial drift of the dust. Finally, in the case of normal filled disk (, , St=), there is no small scale structure in both gas and dust surface density. In this case, the dust-to-gas mass ratio is about 0.01, which is the same as that of the infalling envelope. The diversity of the disk structure can be attributed to the difference of timescales of viscosity, wind mass loss, and dust radial drift, as described in the next section.
3.3. Comparison of Timescales
The formation of various disk structures can be understood by comparing three important timescales: the timescales of the viscous diffusion, the wind mass loss and the radial drift of the dust. We estimate the viscous timescale by using the standard viscous disk with Kepler rotation. The evolution of the surface density only due to the viscosity is given by
where and . We now consider the location close to the local maximum of the pressure, where dust particles are expected to accumulate. At such radius, const, where is the disk scale height. Thus, we obtain and . This gives and the viscous timescale
The wind timescale is
The radial drift timescale is
where we assume St and dust-to-gas mass ratio are much smaller than unity and the gas surface density is given in such a way that Toomre’s -parameter is unity (Equation (20)).
Figure 8 shows the timescales of the viscous diffusion, the wind mass loss and the radial drift of the dust for five parameter sets explained above.
Dust ring disk: The viscous timescale is larger than the wind timescales . Thus, the time evolution of the gas surface density is mainly determined by the disk wind after the wind mass loss starts. The wind timescale around 10-20 au is about yr. Thus, the wind mass loss makes the pressure maximum at 10 au at Myr, which is Myr after the time when the wind mass loss starts. Since the wind timescale is larger for larger radius, the radius of the pressure maximum moves outward with time. This makes the ring move outward as shown in Figure 4 and 6. The drift timescale is comparable to or smaller than the wind timescale. Thus, the dust can concentrate on the pressure maximum formed due to the disk wind and the ring structure is formed.
Dust gap disk: In this case, the time evolution of the gas disk is dominated by the wind mass loss because . In the inner part of the disk ( au), the gas in the disk is removed before dust particles concentrate at the location of pressure maximum since . As a result, the dust inner disk is formed in 20 au. On the other hand, is satisfied at outer radii ( 20 au). This is similar to the case of the ring structure formation described above. Therefore, the ring like structure is also formed in au in this case. The gap structure consists of the inner dust disk and the outer ring structure.
Filled dust disk: In the case that is satisfied in the whole disk, the dust disk without gas is formed since the gas is removed by the wind very rapidly. This structure corresponds to the inner region of the gap structure.
Dust poor disk: When , the hole structure of the gas disk is not formed by the wind mass loss. Since the radial drift is faster than the viscous evolution (), the dust surface density is much smaller than that of the gas, resulting in the dust poor disk.
Normal filled disk: In this case, neither the gas nor the dust hole structures are formed since . However, in contrast to the case of dust poor disk, is satisfied in the disk. Thus, the radial drift of the dust is not significant and the dust-to-gas mass ratio is almost same as that of the infall envelope.
The relation between the timescales and the disk structure is summarized in Figure 9.
3.4. Parameter dependence of ring structure
In this section, we present the results of parameter study, where we have varied , , and St. We investigate how these parameters are related to the disk morphology. We focus on the ring radius, the ring width, the dust surface density at the ring, existence of the inner dust disk interior to the ring (cf. the dust disk in au in left top panel of Figure 7 (dust gap disk)), and dust-to-gas mass ratio in the ring structure. The ring radius is defined by the radius of the local maximum of the dust surface density. The ring width is given by the full width of the half maximum of the dust surface density distribution. We regard the local maximum of the dust surface density as a ring when the inner radius of the ring is larger than the radius of the inner boundary. We only investigate the ring whose radius is larger than 3 au, because the surface density distribution at the innermost region of the computational domain may be affected by the inner boundary condition.
Figure 10 shows the ring radius at Myr for , and .
The green dashed lines indicate the parameters where and satisfy at the inner boundary. The region under this line corresponds to . In this region, the viscous diffusion prevents the inner hole formation in the gas disk by the wind mass loss. As a result, the ring structure is not formed. This disk structure corresponds to the dust poor disk or the normal filled disk. The red dashed lines indicate the parameters where and St satisfy . The parameters right of the red lines and above the green lines correspond to the case where the timescale of the wind mass loss is smaller than the radial drift of the dust (). If , ring structures do not form since dust particles cannot accumulate to the gas pressure maximum before gas is lost by the wind. Therefore, this parameter space corresponds to the filled dust disk. For , the radial drift timescale is of order of yr in the disk. Thus, the ring structure is not formed in most of the parameters of and . It is noted that the case when , and St=0.03 is categorized as no ring structure, despite we see ring structures in other parameters around them. This is partly due to our setting of threshold values to determine “ring structure”. We require that the dust surface density contrast in the ring is larger than factor of two. Around the parameter region, dust surface density contrast is about factor of two, which is close to our threshold value. We do see weak ring structures are formed when , and St=0.03. The same is true with parameters around , and St=0.01.
In Figure 10, the “ring radius” is measured in two kinds of dust distribution. One is the “dust ring disk” and the other is the “dust gap disk”. For the dust gap disk, the bump of the surface density at the outer edge of the gap (e.g, the dust structure in au in the top left panel of Figure 7) is classified as the ring structure. To distinguish the ring structure and the dust gap disk, we show the ratio of the dust mass inside the ring and the total ring dust mass in Figure 11.
For the dust with small Stokes number, the radial drift timescale is large. When the drift timescales is larger than the wind timescale, the expansion of the inner hole of the gas disk due to wind mass loss is faster than the redial drift of the dust. As a result, the inner hole is not formed in the dust distribution for small St. The critical value of as to whether the dust inner disk remains or not can be estimated from the comparison of and as follows. The ring radius is estimated by the wind timescale given by Equation (23). Since the wind mass loss start at Myr and we focus on the ring structure at Myr here, the relation of the ring radius and is obtained from Equation (23) by substituting 0.3 Myr for ,
Then, substituting this radius in Equation (24), we obtain the drift timescale at the ring radius. To make the inner hole in the dust disk, this timescale is required to be smaller than Myr. As a result, we obtain the critical
The red dashed lines in Figure 11 show this critical . This figure shows that the critical estimated here roughly traces the mass ratio . The mass ratio corresponds to that is an order of magnitude smaller than the critical . The disk structure with large amount of dust inside the ring () corresponds to the dust gap disk.
The dust-to-gas mass ratio is the other important value because it is strongly related to the planetesimal formation in the ring. Figure 12 shows the dust-to-gas mass ratio at the ring radius.
This figure shows that the ratio is large for large because the gas surface density decreases rapidly when the wind mass loss is efficient. The ring radius is large when is large, and therefore, dust-to-gas mass ratio is large when the ring radius is large. The dust-to-gas mass ratio is large for small Stokes number. This is because the radial drift of the dust with small Stokes number is not efficient so that the dust does not drift towards the central star before the wind mass loss starts. As a result, the total mass of the dust in the disk is large for small Stokes number, resulting in the large dust-to-gas ratio at the ring.
Finally, we investigate how the ring radius, width, and the maximum surface density depends on input parameters. Figure 13 shows these values for various and with St=0.1 at Myr. The top panel of Figure 13 shows the parameter dependence of the ring radius. For , the ring radius depends only on because the disk evolution is dominated by the disk wind. For , the viscous diffusion is efficient and is satisfied only in the inner region. As a result, the ring radius is small for large . The middle panel of Figure 13 shows the width of the ring. The ring width is of the order of 10 au when and . For , the ring width is of order of 1 au because the large wind timescale allows gas to remain in the disk for long time and the dust particles concentrate strongly at the location of pressure maximum. The bottom panel of Figure 13 shows the maximum surface density of the dust in the ring. For , the wind timescale is larger than the drift timescale so that only weak gas pressure maximum is formed and the dust surface density decreases due to the radial drift to the central star. As a result, the dust surface density with is much smaller than those of the other cases.
Figure 14 shows the ring radius, the width and the maximum surface density of dust for various and St with at Myr. The top panel of Figure 14 is the ring radius. For , the dust strongly couples with the gas and cannot concentrate into the pressure maximum. Thus, no ring structure is formed in the disk. For , the radius of the ring does not strongly depend on St because the radius of the pressure maximum depends mainly on the wind efficiency and viscous diffusion timescale. The middle panel of Figure 14 shows the ring width. For , the ring width decreases with St increases. This tendency is understood as a result of the rapid concentration of the dust in the pressure maximum. The bottom panel of Figure 14 shows the maximum surface density in the ring. For large Stokes number, the radial drift timescale is small and the dust surface density decreases before the wind mass loss stats. As a result, the total dust-to-gas mass ratio and the ring surface density is small for the dust with large Stokes number.
4.1. Dependence on the angular velocity of the cloud cores
In this section, we discuss the effect of the angular velocity of the cloud core, which we have fixed at . The angular momentum affects mainly the dust mass in the disk at the time when the wind mass loss starts. Figure 15 shows the evolution of the surface density with and , which are same as those of Figure 4. The top and bottom panels are the result with and , respectively. For , the maximum centrifugal radius is about 10 au and the infalling envelope accrete onto the inner radius of the disk. Thus, the radius of the gas and dust is smaller than that with . As a result, the total dust mass in the disk is also small. Thus, the resultant ring surface density is smaller than that with . On the other hand, the ring radius is almost the same as that with , because it is mainly determined by the wind timescale. For the case with , the large gas disk is formed. The total dust mass is also larger than that with . Because of the large dust mass, dust-to-gas mass ratio becomes about unity at Myr. In this case, the back reaction from dust to gas prevents the radial drift of the dust. As a result, the width of the ring is larger than that of . In this case, the gas surface density distribution is changed by the back reaction from dust to gas. Thus, the location of pressure maximum moves inward and the ring radius is slightly smaller than that of .
4.2. Constant dust radius
In this paper, we have fixed the Stokes number of dust particles within the disk for simplicity, as motivated by previous study of grain growth. The actual dust properties within the disk is largely unknown, and therefore we have also tried another extreme, where the size of dust particles does not change within the disk. Figure 16 shows the time evolution of the surface density of the gas and the dust with and mm.
In this case, the dust particles are tightly coupled with the gas. In the inner region, the Stokes number is . Thus, the dust-to-gas mass ratio is about in the disk at Myr. The wind mass loss makes the inner hole in the gas disk in Myr. However, the dust inner hole is not formed because the dust is coupled too tightly with the gas to drift in the disk and to concentrate at the pressure maximum of the gas disk. This is even in the case with the large dust radius. Figure 17 shows the evolution of the surface densities of the gas and the dust with and cm, which corresponds to St at 1 au at Myr. In this case, the coupling between the gas and the dust is weaker than the previous case since the dust size is larger. However, the inner hole of the dust is not formed even in this case although the dust surface density in au is smaller than the previous case.
We have checked for other dust sizes and found that the dust inner hole is not formed when the protostar age is Myr if we assume the constant dust radius cm.
4.3. Treatment of the disk wind in Class 0, I YSOs
In our model, the wind mass loss starts after the mass accretion from the cloud core onto the disk terminates for simplicity. In reality, however, the existence of the outflow in Class 0, I phase is reported by both observations (e.g. Arce & Sargent 2006; van der Marel et al. 2013a) and numerical simulations (e.g. Tomisaka 1998; Inutsuka 2012, and references therein). Since the outflow transports the mass and angular momentum of the disk, the outflow in the early star formation phase will decrease the total mass and angular momentum of the Class II star-disk system.
The effect of the outflow in Class 0, I phase on the disk surface density distribution is not yet fully understood. The numerical simulations for the star and disk formation have shown that the disk with outflow is gravitationally unstable (cf. Inutsuka et al. 2010). The surface density profile of such a disk is expected to satisfy . Thus, the disk structure before the infall stops obtained from our model (e.g. upper panel of Fig 4) might mimic that with the outflow. Even when the outflow changes the surface density profile of the gas and the dust at the end of the infall, we expect that the classification of the disk morphology at the early phase of disk formation may be applicable at least qualitatively because it depend only on the timescales of the viscous diffusion, the wind mass loss and the radial drift of the dust.
To confirm the validity of our model and obtain the realistic processes of the disk structure formation, we need to develop the model treating the outflow in Class 0, I phase consistently. This involves detailed comparison between the 1D model and the state-of-the-art 3D numerical simulations, which we shall explore in future work.
4.4. Comparison with the observation of WL 17
In this subsection, we discuss the parameter set of , , and suitable to explain the observed structure in the disk around WL 17. Sheehan & Eisner (2017) presented observations as well as a simple model for this object. In the following discussion, we quote their values as demonstration. The observations show that the ring radius is 10-20 au. According to Figure 10 (see also Figure 13), it corresponds to the result with . To make the ring structure with , it is also required that . Since dust emission from WL 17 shows a ring-like structure, the dust and gas should be reasonably decoupled and the Stokes number larger than about 0.03 is preferred. However, the Stokes number should not be too large considering the disk mass (see Figure 14). Since the disk mass of WL 17 is about 0.05 and the disk area is , the mean surface density of dust is . Comparing this surface density with results of model calculations, Stokes number maybe lower than 0.3. The outflow of the WL 17 is also observed by the CO line emission in van der Marel et al. (2013a). The estimated mass and timescale of the observed outflow are and yr. Thus, the mass loss rate is . As shown in Figure 5, the mass loss rate is of the order of at Myr when , , St=0.1. This mass loss rate is also consistent with the value obtained form observations. Based on the above discussion, the disk structure around WL 17 may be interpreted as the disk with , , and St. In the ring at Myr, Stokes number St=0.1 corresponds to the dust size of a few mm, which is consistent with the fact that WL 17 is observed in the mm-wavelength. These parameters correspond to the viscous timescale yr, the wind timescale yr, and the radial drift timescale is yr. According to Figure 12, the dust-to-gas mass ratio is larger than 0.1 if the structure is formed by the disk wind.
4.5. Evolution of the dust in the disk
In this work, we treat the Stokes number (or dust size) as a parameter. In principle, however, it should be solved simultaneously with the disk evolution. As discussed in Section 4.2, it is difficult to form the ring structure with an inner hole in the dust disk by using the models with a constant dust radius. Thus, dust growth might be important to explain the ring structure like WL 17. As zeroth order approximation, we take into account the growth of the dust by using the models with constant Stokes number according to the one-dimensional calculation of the dust growth Okuzumi et al. (2012). Since the growth timescale of the dust is about yr at 10 au (cf. Takeuchi et al. 1996), we consider that the assumption that the dust grows to St=0.1 within the disk evolution timescale 0.1 Myr is not too unrealistic. However, there might be some problem in this treatment. The most serious problem will be that the constant Stokes number corresponds to the decrease of the dust radius when the surface density of gas decreases due to the wind mass loss. The constant Stokes number assumption may be valid if the destruction of the dust due to the collision or supplying the small dust from the outer radius by the drift of the dust makes the dust radius decrease. If the dust radius does not change when the gas surface density decreases, Stokes number of the dust increases. This makes the dust concentration on the pressure maximum faster than the results presented in this paper. Thus, the condition of the ring structure formation may not be too affected, or it may be possible to form ring structures with smaller initial Stokes number.
4.6. Effects of magnetic field on disk formation and evolution
In our model, we assume the conservation of the angular momentum in the infalling envelope for simplicity. When the cloud core is magnetized, however, the angular momentum is transferred outward in the envelope. This is called “magnetic braking”. The initial configuration of a cloud core and non-ideal MHD processes both affect the angular momentum transfer in the infalling envelope (cf. Li et al. 2014; Machida et al. 2014; Tsukamoto et al. 2015a), while we still do not have quantitative models for these effects. Since the angular momentum transfer in the envelope affects the disk radius, we need further investigation of modeling of the gravitational collapse of the cloud core taking into account the magnetic braking.
The strength of the turbulence and the efficiency of the wind mass loss depend on the magnetic field in the disk. Thus, in principle, these parameters cannot be chosen independently but related to each other. The assumption that and are constant corresponds to the smooth distribution of the magnetic field in the disk. However, it depends on the radius in reality (cf. Suzuki & Inutsuka 2014; Bai & Stone 2014; Flock et al. 2015). Although the and depend on the radius of the disk, the condition for the ring structure formation based on the comparison of the timescales obtained in this work will be available. According to Suzuki et al. (2010), is about an order of magnitude larger than (Figure 8 in Suzuki et al. (2010), see also Okuzumi & Hirose (2011)). We have observed the inner hole formation when and . These values are roughly consistent with or the wind mass loss is slightly (a factor of a few) more effective than those obtained in the local simulations.
The non-ideal MHD effects are also important for the disk structure formed by the disk wind. In Suzuki et al. (2010), only the Ohmic resistivity is taken into account. As a result, the layered structure consisting of the MRI-inactive midplane and the MRI-active surface is realized. On the other hand, the MRI at the surface layer is suppressed by the ambipolar diffusion Bai & Stone (2013); Bai (2013); Gressel et al. (2015) and the accretion flow is laminar in the inner region of the disk. In such a case, the disk accretes due to the wind torque.
Moreover, the Hall effects also changes the angular momentum transfer and wind mass loss rate Bai (2014, 2015); Lesur et al. (2014); Simon et al. (2015). In the inner disks ( au), the disk evolution depends on the orientation of the magnetic field to the rotation axis of the disk, . If , the azimuthal magnetic field is amplified leading to the efficient angular momentum transfer and wind mass loss. On the other hand, if , the horizontal magnetic field is reduced and the Maxwell stress in the disk and wind mass loss rate decreases. Thus, the ring structure formation in a young disk will be difficult for the disk with because the wind mass loss timescale may be longer than the age of the disk. In the outer disk ( au), the surface layer of the disk is ionized with FUV irradiation and MRI turbulence is sustained there. As a result, efficient mass accretion is expected at outer radii. If the gas from the outer disk piles up at the boundary of the inner region, the additional ring structure might be formed.
Although the ring structure formation is affected by the complicated non-ideal MHD effects, the effects have not been fully understood quantitatively. Further investigations of global, 3D simulations on magnetorotational instability are necessary to construct more realistic models of early disk evolution.
In this work, we investigate the ring structure formation in young disks by the wind mass loss caused by the MRI turbulence. To calculate the ring formation in young disks, we use a one-dimensional disk model that treats the formation and evolution of disks in a single framework. In this model, the strength of the turbulence, the mass loss rate by the disk wind, and dust size (Stokes number) are treated as parameters, and the dependence of the disk evolution on these parameters are investigated. Main results obtained in this work are summarized as follows:
We find five types of disk structures as a result of disk formation and evolution model including the effects of wind mass loss: the ring structure, the dust gap disk, the filled disk, the dust poor disk, and the normal filled disk are obtained from our model calculations.
The disk evolution is characterized by the timescales of the viscous diffusion , the wind mass loss , and the radial drift of the dust . The formation of various disk structures can be understood by comparing these timescales.
The relation between the timescales and the disk structure is summarized in Figure 9. When , the inner hole structure is formed in the gas disk. In this case, the dust can concentrate at the pressure maximum of the gas disk and the ring structure is formed when . In the case or in the entire disk, the gas inner hole expands faster than the radial drift of the dust. As a result, the dust disk remains in the inner hole of the gas disk and the filled dust disk is formed. If the disk satisfies or only in the inner region and the outer region satisfies , the dust gap disk is formed. When , the inner hole structure of the gas disk is not formed. In this case, the dust ring structure is not formed either. If , the dust in the disk drifts inward faster than the gas disk evolution. As the result, the dust poor disk is formed. On the other hand, the dust-to-gas mass ratio is almost same as the initial value (0.01) if or . This case is categorized as the normal filled disk.
When the ring structure is formed, the ring radius and the dust-to-gas mass ratio at the radius of the dust surface density maximum increase with time. When the dust-to-gas mass ratio becomes larger than about unity, the back reaction from the dust to the gas becomes efficient and the ring radius does not increase after that.
To explain the ring structure observed around WL 17, , and St are suitable based on our model calculations. The resultant structures obtained by using these parameters are consistent with the ring radius, the dust surface density of the ring, formation of the inner hole in the dust disk, and mass loss rate estimated from observations. Our model calculations suggest that the dust-to-gas mass ratio at the ring is larger than 0.1 and the outer radius of the gas disk is larger than that of the dust ring.
In reality, the dust radius, strength of the turbulence, and the mass loss rate due to the wind are given by the result of the dust growth and the evolution of the magnetic field in the disk. The calculation of these values with the disk evolution is required to obtain more realistic structure of the disk, which will be the subject of our future work.
We thank Takeru K. Suzuki for fruitful discussions and his valuable comments. This work was supported by NAOJ ALMA Scientific Research Grant Numbers 2016-02A. TM is supported by JSPS KAKENHI Grant Nos. 17H01103, 15H02074, and 26800106.
Appendix A Derivation of the velocity of the gas and dust in the disk
The equations for the gas and the dust in the disks are given as follows:
where , , are the surface density, the radial velocity, and the azimuthal velocity of the gas, , , are those of the dust, is the mass loss rate due to the disk wind, is the vertically integrated pressure, is the enclosed mass of the gas within the radius , is the stopping time of the dust, is the coefficient of the kinematic viscosity, and is the dust-to-gas mass ratio in the infalling envelope.
We assume that the effect of the pressure gradient force, the frictional force between the gas and the dust, and the viscosity is small compared with the centrifugal force and the gravitational force. If we neglect these forces, we obtain and . We calculate the deviation of the velocities from these values in the first order. Here we treat that is the same order as . We define and as the deviations of the azimuthal velocities. From Equation (A2) and (A5), we obtain
Equation (A7) gives
where . To obtain the radial drift of the dust in the disk midplane, we evaluate the first term by using defined in the midplane, where and is the pressure in the midplane. Equation (A8) gives
Thus, we obtain