Classification of the Circumstellar Disk Evolution During the Main Accretion Phase
We have carried out hydro-dynamical simulations to investigate the formation and evolution of protostar and circumstellar disks from the prestellar cloud. As the initial state, we adopt the molecular cloud core with two non-dimensional parameters representing the thermal and rotational energies. With these parameters, we make 17 models and calculate the cloud evolution years after the protostar formation. We find that early evolution of the star-disk system can be qualitatively classified into four modes: the massive disk, early fragmentation, late fragmentation, and protostar dominant modes. In the ’massive disk mode’ to which the majority of models belong, the disk mass is greater than the protostellar mass for over years and no fragmentation occurs in the circumstellar disk. The collapsing cloud shows fragmentation before the protostar formation in the ’early fragmentation mode’. The circumstellar disk shows fragmentation after the protostar formation in the ’late fragmentation mode’, in which the secondary star substantially gains its mass from the circumstellar disk after fragmentation and it has a mass comparable to that of the primary star. The protostellar mass rapidly increases and exceeds the circumstellar disk mass in the ‘protostar dominant mode’. This mode appears only when the initial molecular cloud core has a considerably small rotational energy. Comparison of our results with observations indicates that a majority of protostar has a fairly massive disk during the main accretion phase: the circumstellar disk mass is comparable to or more massive than the protostar. It is expected that such a massive disk promotes gas-giant formation by gravitational instability in a subsequent evolution stage.
keywords:star formation – circumstellar disk – methods: hydrodynamics – smoothed particle hydrodynamics – protoplanetary disk – planet formation
There are two major mechanisms for gas-giant planet formation: one is the core accretion mechanism in which a massive solid core forms first and the disk gas accretes onto the core (Safronov, 1969; Goldreich & Ward, 1973; Pollack et al., 1996), and the other is the gravitational instability (GI) mechanism in which the circumstellar disk directly fragments into gas-giant planets via GI (Cameron, 1978). Recent discovery of extra-solar planets at a great distance from the central star such as HR8799b, c, d and e (Marois et al., 2008, 2010) and GJ579b (Thalmann et al., 2009) creates a new problem for the planet formation. It is difficult to form planets in the regions far from central stars according to the core accretion mechanism, because massive solid core formation before the dissipation of gaseous disk seems to be difficult (Dodson-Robinson et al., 2009).
The gravitationally instability mechanism may be more plausible for the formation of these planets. Many studies of disk fragmentation have been done using either an analytic approach (Rafikov, 2005) or numerical simulations (e.g., Stamatellos & Whitworth, 2007; Cai et al., 2008; Boley et al., 2006; Mejia et al., 2005; Pickett et al., 2003; Meru & Bate, 2010). These efforts, however, seem to build a consensus that the planet formation by GI within AU is highly difficult when the disk-to-stellar mass ratio is whose ratio is suggested by observations (see, e.g., Kitamura et al., 2002).
On the other hand, Inutsuka et al. (2010) showed that the circumstellar disk is comparable to or more massive than the protostar () during the (early) main accretion phase (i.e., Class 0 or Class I stages) and is highly gravitationally unstable. Recently, such massive disks were observed around very young protostars (e.g., Eisner & Carpenter, 2006; Enoch et al., 2009). Fragmentation in such massive disks may potentially account for the formation of gas-giant planets with a wide separation such as HR 8799 or GJ579. Thus, we may have to focus on fragmentation and gas-giant planet formation in a massive disk during the main accretion phase. There are several unsolved problems: how long such a massive disk exists, or whether the massive disk can fragment to form sub-stellar objects.
So far, a few authors have investigated the early evolution of the circumstellar disk and its fragmentation from prestellar core stage. Using the thin disk approximation, Vorobyov & Basu (2010a) calculated the formation and evolution of the circumstellar disk in the molecular cloud core with a barotropic equation of state (or including cooling and heating effects, Vorobyov & Basu, 2010b). They found that a gravitationally unstable disk frequently appears in the star formation process and tends to show fragmentation during the main accretion phase. With isothermal equation of state, Kratter et al. (2010) investigated the formation of the circumstellar disk in their three dimensional calculation. They also showed massive disk formation during the main accretion phase and concluded that the disk fragments if the disk-to-stellar mass ratio is greater than . Using an adiabatic equation of state and radiative cooling, Walch et al. (2009) studied the circumstellar disk formation with a somewhat coarse spatial resolution. They showed that a massive disk becomes hot during the early stage of the main accretion phase and fragmentation is suppressed.
Machida et al. (2010) also investigated the disk formation with a finer spatial resolution than Walch et al. (2009) using a barotropic equation of state. In the collapsing cloud, the adiabatic heating dominates the radiative cooling at cm and the first (adiabatic) core with a size of 1 AU forms before the protostar formation (Larson, 2001; Masunaga & Inutsuka, 2010). Machida et al. (2010) showed that the first core formed before the protostar formation directly evolves into the circumstellar disk after the protostar formation. Thus, to investigate the (early) formation of the disk, we need to resolve the first core that has a size of AU. In their calculation, however, they implicitly assumed the point-symmetric structure because they fixed the protostellar location to the center of the computational domain. Michael & Durisen (2010) pointed out that the stellar motion weakens GI activity in some degree, when the disk-to-stellar ratio is . Thus, a more adequate treatment may be necessary for the central protostar.
All previous studies of circumstellar disk formation have shown that the circumstellar disk can become more massive than a protostar during the mass accretion phase. However, in such studies, the disk evolution was investigated in a limited parameter range (Walch et al., 2009; Machida et al., 2010; Machida & Matsumoto, 2011). Thus, we cannot conclude whether such a massive disk that is a favorable site for gas-giant planet formation generally appears in the star formation process. To determine this, we need to investigate the disk formation resolving 1AU structure in a wider parameter range.
The initial cloud for star formation is conventionally characterized by two parameters: the ratio of thermal and rotational energy toward the gravitational energy (for detailed description, see §2.2). With these parameters, many studies have calculated and classified the cloud evolution in the isothermal (Boss, 1993; Miyama, 1984; Tsuribe, 2002) and adiabatic (Cha & Whitworth, 2003; Matsumoto & Hanawa, 2003) gas contraction phases to investigate fragmentation before the protostar formation. However, the evolution and fragmentation of the circumstellar disk that is formed after the protostar formation has not been investigated with these parameters.
In this study, with parameters of and , we simulate the formation and evolution of circumstellar disk from molecular cloud core until years after the protostar formation using a smoothed particle hydrodynamics code with sufficient spatial resolution. Parameters and significantly influence the disk formation, because the accretion rate onto the circumstellar disk and the size and mass of the circumstellar disk are determined by these parameters (see, §4). We calculate the cloud evolution for 17 models using the barotropic equation of state.
2 Numerical Method and Initial Condition
Our simulations are carried out using an SPH code newly developed for this work. The code includes an individual time-steps technique and uses Barnes-Hut tree algorithm to calculate self-gravity with opening angle . We use adaptive softening length according to Price & Monaghan (2007). Artificial viscosity is included according to the prescription by Monaghan (1997) with plus the Balsara switch (Balsara, 1995). Our code is parallelized with MPI.
To mimic the thermal evolution of molecular cloud, we use the barotropic equation of state as
where and is adopted. As the initial state, we take a uniform density sphere with axisymmetric density perturbation. The density perturbation is given as . The initial cloud is rigidly rotating. We parameterized the ratio of thermal and rotational energy to the gravitational energy of the initial cloud. They are characterized by two non-dimensional parameters
where , and are the initial cloud radius, angular velocity and mass of cloud core, respectively. We vary parameters and in the range of and . In this study, total cloud mass is fixed to in all models to compare the evolution of cloud with the same mass. In addition, we fixed the initial cloud temperature to 10 K. Thus, parameter and initial cloud density are determined when the initial cloud radius is given. Note that such treatment changes the gravitational energy of the cloud, After the cloud radius is fixed with arbitrary , parameter is determined when the initial angular velocity is given. The model name, parameters and , initial cloud radius, initial angular velocity and initial cloud density are listed in Table 1. The initial cloud is modeled with about 500 000 SPH particles.
To calculate the disk formation several years after the protostar formation, we adopt a sink particle technique according to the prescription by Bate et al. (1995). Starting from the prestellar core stage, the cloud evolution is calculated without sink particle. Then, we assume the protostar formation and dynamically introduce a sink particle when the gas particle density exceeds the threshold density, . The threshold density roughly corresponds to the density at which the second collapse begins. As shown in Larson (2001) and Masunaga & Inutsuka (2010), the second collapse begins when the gas density reaches at which the gas temperature exceeds K and molecular hydrogen begins to dissociate. The protostar forms immediately after the second collapse. Thus, in this paper, we safely define the protostar formation epoch as that at which the gas density exceeds .
To treat the gas accretion onto the sink particle after the creation of the sink particle, we set the accretion radius . Then, within the accretion radius, we allow the gas accretion onto the sink particle when the following condition is fulfilled: (1) the gas particle density exceeds the accretion density , (2) it is gravitationally bound to the sink particle and (3) its specific angular momentum is less than that required for it to form a circular orbit at . We did not implement a boundary condition for sink particle.
|Model||R (AU)||Mode||number of star||
”M”, ”EF”, ”LF”, ”P” mean ”Massive disk mode”, ”Early fragmentation mode”, ”Late fragmentation mode” and ”Protostar dominant mode”, respectively.
”Massive disk era” is defined as the period during which the disk mass is greater than the protostar mass. Some of simulations are terminated during the massive disk era due to computational limits and we represent it with .
As listed in Table 1, we calculated the cloud evolution for 17 models in total, and classified them into four modes as follows:
Massive disk mode: the disk mass dominates the protostellar mass for over years after the protostar formation.
Early fragmentation mode: fragmentation occurs in the collapsing cloud before the protostar formation.
Late fragmentation mode: fragmentation occurs in the circumstellar disk after the protostar formation.
Protostar Dominant mode: the protostellar mass rapidly dominates the disk mass within years after the protostar formation.
The calculation results and its classification are summarized in Figure 1.
3.1 Massive Disk Mode
Massive disk mode is indicated by the circles in Figure 1. Figures 2 and 3 show the time evolution of the center of the cloud for model 8 ( and ) that is a typical model for massive disk mode. The upper left and upper middle panels in these figures show snapshots before the protostar formation. A strong bar-like structure appears in the central high-density gas region in the top left panel in Figure 2. Then, the bar structure effectively transfers the angular momentum outward and the high-density gas region shrinks. As a result, a bimodal structure composed of the central high-density core and its surrounding disk appears as seen in the top middle panel. These figures clearly show the disk formation before the protostar formation. The top right panel shows the snapshot just when the second collapse occurs (i.e., the protostar forms). The time when the protostar forms is years from beginning. This roughly corresponds to the free-fall timescale of the initial cloud. By this epoch, the spiral arm has developed around the central object. The bottom left, middle, and right panels show the snapshots and years after the protostar formation. These figures show that the disk gradually increases its size with time retaining the global spiral structure.
Figure 4 shows the mass evolution of the disk and protostar. For this model, the disk mass is greater than the protostellar mass for more than years after the protostar formation. Such a disk is expected to be gravitationally unstable. To investigate the disk stability, the contours of Toomre’s Q parameter at and years after the protostar formation (the same epoch of the bottom left and middle panels of Fig. 2) are plotted in Figure 5. The Toomre’s Q parameter is described as
where , and are the sound velocity and epicyclic frequency and surface density of the disk, respectively. For this model, the circumstellar disk did not show fragmentation even though the circumstellar disk has the region of at years after the protostar formation. This is because GI cannot grow sufficiently fast in the circumstellar disk. The characteristic timescale of Toomre’s analysis, , is comparable to the orbital period of the disk (both are several thousand years). In this case, there are non-linear stabilize mechanisms against GI. As seen in the right panel of Figure 5 or bottom middle panel of Figure 2, the spiral arm globally redistribute the mass and angular momentum in a short duration, years ().
Furthermore, the disk-star configuration also dynamically changes in a short duration. Figure 6 shows the trajectory of the protostar for years after the protostar formation. The asterisk indicates the position where the protostar forms. The central protostar drifts toward the dense part of the spiral arm due to the gravitational interaction between the protostar and spiral arms. Since the strength of Keplerian shear is proportional to , its radial dependency is very strong and the drift motion increases the shear stress of the dense region and suppresses fragmenting of the disk. As described in Gammie (2001), the thermal evolution of the disk also may play an important role to suppressing/promoting fragmentation. They showed that no fragmentation occurs when the cooling time is much longer than the orbital period. In our simulation, since we adopted the barotropic equation of state (see, §2.2), the disk evolves adiabatically (practically all the time, the disk mid-plane density is greater than ). With these effects, the disk is stabilized as seen in the right panel of figure 5. For disk fragmentation, the mass in-fall must be sufficiently fast to overcome these non-linear back reactions.
3.2 Early Fragmentation Mode
The open triangles in Figure 1 indicate ”Early fragmentation mode,” in which the collapsing cloud undergoes fragmentation before the protostar formation. Figure 8 shows the time evolution of the center of the cloud for model 11 ( and ). The figure shows that the ring pattern arises (top left panel) and it fragments to form four clumps (or protostars). These protostars have roughly the same mass when they form. After fragmentation, they interact with each other. Finally, a quadruple stellar system composed of two close binary appears as seen in the lower right panel of Figure 8. The condition of the ring like structure was substantially investigated in Cha & Whitworth (2003). Since they fixed parameter (but changed the rotational low), we cannot quantitatively compare their results with ours. However, our result is not contradict with them.
In models showing fragmentation before the protostar formation (i.e., early fragmentation mode), fragmentation pattern qualitatively changes as the rotation energy of the initial cloud decreases. Figure 9 shows the time evolution of the center of the cloud for the model 12 ( and ). A clear bar-like structure arises before the protostar formation for model 12. The bar fragments into two clumps. Then, a binary system appears around the center of the collapsing cloud (upper middle panel of Fig. 9). Furthermore, the remnant of the bar fragments and form two extra protostars about years after the first fragmentation, as seen in the upper right panel of Figure 9. Then, one protostar is ejected from the central region by gravitational interaction among protostars, and a triple stellar system remains around the center of the cloud at the end of the calculation.
Fragmentation and subsequent evolution of fragments are very complicated in models belonging to the early fragmentation mode. However, since fragmentation before the protostar formation was well investigated in previous studies (see, Tsuribe 2002, Goodwin et al. 2007 and Bodenheimer et al. 2000 for review), we do not comment it any more in this paper.
3.3 Late Fragmentation Mode
The filled triangles in Figure 1 indicates the ”Late fragmentation mode,” in which fragmentation occurs in the circumstellar disk after the protostar formation. Models 13 and 14 belong to this mode. Figure 10 shows the snapshots for model 14 ( and ). As seen in the massive disk mode, after a strong bar develops (top left panel), bimodal structure composed of the central high-density gas region and disk-like structure appears even for model 14 (top middle panel). However, fragmentation occurs in the disk before the spiral arm sufficiently develops (top right panel) in this model, while no fragmentation occurs in the models belonging to the massive disk mode (see, Figs. 2 and 3). After fragmentation, the secondary object appears in the region AU far from the primary star (bottom left panel) years after the primary star formation. Thus, Figure 10 indicates that fragmentation can occur in the circumstellar disk during the main accretion phase if the mass in-fall is sufficiently fast (with smaller , see §5).
The mass evolution of two protostars is plotted in Figure 7, in which the secondary star increases its mass in a way similar to that of the primary star. At the end of the calculation, the mass of secondary star is comparable to that of the primary star. Thus, this system evolves into a binary stellar system, not a star-planet system.
3.4 Protostar Dominant Mode
The crosses in Figure 1 indicate ”protostar dominant mode”, in which the protostar rapidly increases its mass and exceeds the circumstellar disk mass within years. For this mode, the protostellar mass exceeds the disk mass within years after the protostar formation. Figures 11 and 12 show the time evolution of the center of the cloud for model 16 ( and ) that is a typical model for the protostar dominant mode. The top left and middle panels show the snapshots before the protostar formation for model 16. Unlike the massive disk mode, the first core keeps an almost axisymmetric structure without developing the non-axisymmetric perturbation. The second collapse occurs at years from the beginning. Even after the protostar formation, the axisymmetric structure is maintained until the disk radius sufficiently grows ( AU).
To investigate the disk stability for model 11, the Toomre’s Q parameter at years after the protostar formation (the same epoch as the bottom middle of Fig. 11) is plotted in Figure 13. At this epoch, the disk mass is still greater than the protostellar mass (see, Fig. 14). Nevertheless, the Q parameter is greater than unity in the whole region of the disk due to compactness of the disk.
The mass of the protostar and circumstellar disk is plotted against time after the protostar formation in Figure 14. The figure indicates that the accretion rate onto the protostar is higher than that of the circumstellar disk, and the protostellar mass exceeds the circumstellar disk mass at years after the protostar formation.
4.1 Effect of Radiative Cooling
In this study, we used the barotropic equation of state (eq. ), in which the gas in the disk behaves adiabatically when the disk mid-plane density is greater than cm. However, in reality, the disk radiatively cools with time. According to Rafikov (2005), we can roughly estimate the cooling time of the disk as
where , , T, and are the surface density, Stephan-Boltzmann constant, mid-plane temperature, optical depth of the disk and mean molecular weight, respectively. In equation (5), we estimate the vertical optical depth of the disk as In addition, we assume the optical depth as , where is the opacity. We assume that the disk is optically thick, and derive RHS of eq (5) with .
Equation (5) indicates that the cooling timescale of the circumstellar disk is typically about years. Thus, the radiative cooling may play an important role in investigating the thermal evolution of the disk and its fragmentation over years. In this study, however, we showed that the evolution of the circumstellar disk can be qualitatively classified into four modes during very early stages of the star formation ( years). Thus, we expect that our classification does not change qualitatively even when the radiative effects are included. However, the radiative cooling of disk is important to investigate the disk evolution and fragmentation for years.
4.2 Treatment of Protostar and Sink
In this study, no model shows the formation of a star-planet system during the main accretion phase because the secondary object continues to increase its mass and finally exceeds the hydrogen-burning limit ( M). On the other hand, Vorobyov & Basu (2010a) and Machida et al. (2010) showed the formation of a star-planet system during the main accretion phase. The difference is thought to be caused by treatment of the protostar and sink.
In Vorobyov & Basu (2010a) and Machida et al. (2010), the protostar (or sink cell) is fixed at the center of the computational domain. It is expected that such treatment promotes fragmentation. In reality, the density fluctuation arising around the protostar can cancel out by movement of the protostar. Thus, fragmentation tends to occur when the protostar is fixed.
In addition, Vorobyov & Basu (2010a) and Machida et al. (2010) did not impose the sink on fragments formed in the circumstellar disk. Instead, they suppressed further collapse of fragments with adiabatic equation of state. Such treatment decreases the mass accretion onto fragments in some degree. On the other hand, our sink treatment may overestimate the mass accretion onto fragments or protostar.
4.3 Comparison with Previous Works
So far, few authors have investigated the evolution of circumstellar disk from molecular cloud core. With an isothermal equation of state, Kratter et al. (2010) showed that fragmentation occurs in the circumstellar disk with a wide parameter space during the early stage in the main accretion phase and claimed that fragmentation frequently occurs when the disk-to-stellar mass ratio is greater than unity. However, their assumption of isothermality seems not to be valid for the disk evolution of the early main accretion phase, because the gas becomes opaque and behaves adiabatically when the gas density exceeds the critical density of (e.g., Masunaga & Inutsuka 2010). In addition, the radiative cooling can affect the disk evolution years after the protostar formation as described in §4.1. Since the adiabatic equation of state stabilizes the circumstellar disk, fragmentation barely occurs in our calculation even when the disk-to-stellar mass ratio exceeds unity. Thus, Kratter et al. (2010) may overestimate the fragmentation condition, while our calculation may underestimate it because of lack of radiative cooling.
Walch et al. (2009) also studied the circumstellar disk formation with the adiabatic equation of state and radiative cooling. They showed no fragmentation in the circumstellar disk, because the disk becomes very hot during the early stage of main accretion phase. Their spatial resolution is, however, somewhat coarse; the minimum smoothing length of their study is AU. On the other hand, AU in our simulations. In addition, they restricted their initial conditions to rapidly rotating cases () to investigate disk evolution before the central density becomes high. The observation suggested that molecular cloud cores have the rotational energy of with a typical value of (Caselli et al., 2002). Thus, they studied the cloud evolution in the limited parameter range.
We have carried out hydro-dynamical simulation to investigate the evolution of the circumstellar disk with two non-dimensional parameters representing the thermal and rotational energy of the initial cloud. The thermal energy is related to the mass accretion rate onto the circumstellar disk as (see, Machida et al. 2011). Thus, smaller provides a high accretion rate onto the circumstellar disk and vice versa. On the other hand, the initial rotational energy that is represented by parameter is related to the disk radius. The centrifugal radius of the initial cloud is related to as , where is the initial cloud radius. Thus, with larger , the cloud forms a larger disk in the main accretion phase. In other words, with larger , a large fraction of the in-falling matter accretes onto the disk, rather than directly onto the primary protostar. As a result, smaller and larger increase disk surface density and makes a gravitationally unstable disk.
On the other hand, the non-axisymmetric structure arose in such unstable disk can stabilize the disk, because it can redistribute the angular momentum and promote mass accretion onto the protostar. Thus, no fragmentation occurs when a strong non-axisymmetricity grows and transfers sufficient angular momentum outward. By contrast, the disk becomes highly gravitationally unstable and shows fragmentation when non-axisymmetric structure does not grow sufficiently or when the growth timescale of the non-axisymmetricity is much longer than the disk growth timescale. Thus, fragmentation condition depends also on the growth of the non-axisymmetriciy, which is closely related to parameters and because they determine the evolution of the mass and angular momentum of the disk.
With parameters and , we found that the disk evolution is qualitatively classified into four modes: protostar dominant, massive disk, early fragmentation and late fragmentation modes. The schematic classification of the circumstellar disk is shown in Figure 1 that covers a wide parameter range of the rotational energy supported by observations (; Caselli et al., 2002). We describe each mode in the following.
Protostar Dominant Mode: For protostar dominant mode, the protostellar mass exceeds the circumstellar disk mass within years after protostar formation. This mode appears in the range of when the protostar forms, the circumstellar disk is more massive than the protostar. However, for this mode, since the mass increase rate of the protostar is larger than that of the circumstellar disk, the protostellar mass exceeds the circumstellar disk mass in years after the protostar formation, as described in §3.4.
Massive Disk Mode: For the massive disk mode, the circumstellar disk is more massive than the protostar over years after protostar formation. As seen in Figure 1, this mode appears when the initial cloud has larger rotational and thermal energies. In our calculation, a majority of models (8 out of 17 models) belongs to this mode. This indicates that the circumstellar disk that is more massive than the protostar frequently appears in the star formation process. The circumstellar disk is self-regulated for this mode. When the circumstellar disk becomes massive, it decreases its mass because the non-axisymmetric structure develops and promotes the mass accretion onto the protostar. On the other hand, when the disk is relatively less massive, the disk mass increases because a relatively stable disk acquires its mass from the in-falling envelope without effective angular momentum transfer.
Early Fragmentation Mode: For early fragmentation mode, fragmentation occurs before the protostar formation. As seen in Figure 1, this mode appears when the initial cloud has large rotational but small thermal energies. If the initial rotational energy is large enough, the gas cannot condense to the center, and a ring like structure appears as seen in figure 8. This structure is unstable against perturbation and it fragments into protostars. As the rotational energy decreases, the gas can condense to the center and bar like structure develops. However, unlike in the massive disk mode, since the bar cannot transfer angular momentum outward fast enough, fragmentation occurs and binary system appears as seen in Figure 9.
Late Fragmentation Mode: For late fragmentation mode, fragmentation occurs after the protostar formation. The rotational energy for this mode is smaller than that for early fragmentation mode. Thus, without fragmentation before the protostar formation, the central bar like structure shrinks transferring angular momentum outward and form a single protostar and a disk.
In this study, since we did not include proper radiation treatment, we cannot discuss whether a gas-giant planet can form in such a massive disk. However, such a massive disk is a plausible site for the gas-giant formation by GI. To determine whether a gas-giant planet forms in the circumstellar disk by GI, we need to perform more long-term calculation including radiative effects in the future.
We thank T. Tsuribe, S. Inutsuka, K. Osuga, K. Tomida and J. Makino for fruitful discussions. We also thank anonymous referee for helpful comments. The snapshots were produced by SPLASH (Price, 2007). This work was supported by the Grants-in-Aid from MEXT (21740136). The computations were performed on a parallel computer, XT4 system at CfCA of NAOJ. Y.T. is financially supported by Research Fellowships of JSPS for Young Scientists.
- Bate M. R., Bonnell I. A., Bromm V., 1995, MNRAS, 277, 362
- Balsara D. S., 1995, J. Comput. Phys., 121, 357
- Boss A.P., 1993, ApJ., 410, 157
- Bodenheimer, P., Burkert, A., Klein, R. I., & Boss, A. P. 2000, Protostars and Planets IV, 675
- Boley A.C. et al, 2006, ApJ, 651, 517
- Cai K. et al, 2008, ApJ, 673, 1138
- Cameron A. G. W., 1978, M&P, 18, 5
- Caselli P. et al, 2002, ApJ, 572, 238
- Cha, S. H.,,Whitworth, A. P., 2003, MNRAS, 340, 91
- Dodson-Robinson, S. E. et al, 2009, ApJ, 707, 79
- Eisner J. A., Carpenter J. M., 2006, ApJ, 641, 1162
- Enoch, M. L., Corder, S., Dunham, M. M., & Duchêne, G. 2009, arXiv:0910.2715
- Gammie C. F, 2001, ApJ, 553, 174
- Goldreich P.,Ward, W. R., .,1973, ApJ, 183, 1051
- Goodwin, S. P., Kroupa, P., Goodman, A., & Burkert, A. 2007, Protostars and Planets V, 133
- Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347
- Inutsuka, S., Machida, M. N., & Matsumoto, T. 2010, ApJ, 718, 58
- Kratter K. M. et al, 2010, ApJ, 708, 1585
- Kitamura Y. et al, 2002, ApJ, 581, 357
- Larson R. B., 1969, MNRAS, 145, 271
- Lin C. C., Lau Y. Y., 1979, Studies in Applied Mathematics, 60, 97
- Machida, M. N., & Matsumoto, T. 2011, MNRAS, 284
- Machida M. N. et al, 2010, ApJ, 724, 1006
- Machida M. N. et al, 2011, ApJ, 729, 42
- Marois et al. 2008, Science, 322, 1348
- Marois et al. 2010, Natur, 468, 1080
- Masunaga H., Inutsuka S., 2000, ApJ, 531, 350
- Matsumoto T., Hanawa T., 2003, ApJ, 595, 913
- Mejia A.C. et al, 2005, ApJ, 619, 1098
- Michael S., Durisen R. H.,2010, MNRAS, 406, 279
- Miyama S., 1984, ApJ, 279, 621
- Monaghan J, J., 1997, J. Comput. Phys., 136, 298
- Meru F., Bate R.B., 2010, MNRAS, 406, 2279
- Pickett B. K.,2003 ApJ, 590, 1060
- Pollack et al. 1996, Icar, 124, 62
- Price D. ,2007, PASA, 24, 159
- Rafikov R. R., 2005, ApJ, 621, 69
- Safronov V. S., 1969, QB, 981, 26
- Saigo K., Tomisaka K., 2006, ApJ, 645, 381
- Stamatellos D., Whitworth A.P., 2008, A&A, 480, 879
- Tsuribe T, 2002, PThPS, 147, 155
- Thalmann et al., 2009, ApJ, 707, 123
- Vorobyov E. I., Basu S, 2010, ApJL, 714, 133
- Vorobyov E. I., Basu S, 2010, ApJ, 719, 1896
- Walch S. et al., 2010, MNRAS, 402, 2253