Conditions for Gravitational Instability in Protoplanetary Disks

Conditions for Gravitational Instability in Protoplanetary Disks

Shigeo S. Kimura and Toru Tsuribe Department of Earth and Space Science, Graduate School of Science, Osaka University,
1-1 Machikaneyama-cho, Toyonaka-shi, Osaka 560-0043
Abstract

Gravitational instability is one of considerable mechanisms to explain the formation of giant planets. We study the gravitational stability for the protoplanetary disks around a protostar. The temperature and Toomre’s Q-value are calculated by assuming local equilibrium between viscous heating and radiative cooling (local thermal equilibrium). We assume constant viscosity and use a cooling function with realistic opacity. Then, we derive the critical surface density that is necessary for a disk to become gravitationally unstable as a function of . This critical surface density is strongly affected by the temperature dependence of the opacity. At the radius AU, where ices form, the value of changes discontinuously by one order of magnitude. This is determined only by local thermal process and criterion of gravitational instability. By comparing a given surface density profile to , one can discuss the gravitational instability of protoplanetary disks. As an example, we discuss the gravitational instability of two semi-analytic models for protoplanetary disks. One is the steady state accretion disk, which is realized after the viscous evolution. The other is the disk that has the same angular momentum distribution with its parent cloud core, which corresponds to the disk that has just formed. As a result, it is found that the disks tend to become gravitationally unstable for because ices enable the disks to become low temperature. In the region closer to the protostar than , it is difficult for a typical protoplanetary disk to fragment because of the high temperature and the large Coriolis force. From this result, we conclude that the fragmentation near the central star is possible but difficult.

Kimura S. & Tsuribe T.The Condition of GI in Protoplanetary Disk

\KeyWords

accretion, accretion disks — instabilities — (stars:) planetary systems: protoplanetary disk — (stars:) planetary systems: formation

1 Introduction

Two major models are currently considered for the formation of gas giant planets. One is the core accretion model (CA model; [Goldreich & Ward (1973)]; [Hayashi (1981)]; [Pollack et al. (1996)]), and the other is the gravitational instability model (GI model; [Cameron (1978)]; [Durisen et al. (2007)]). In CA model, massive solid cores are made from dust first, and then gas giant planets are formed by the gas accretion onto these cores. In GI model, a disk around a protostar fragments into pieces by gravitational instability and these pieces become gas giant planets. Recently, extrasolar planets are discovered by direct imaging ([Marois et al. (2008)]; [Kalas et al. (2008)]). These planets are heavier than Jupiter and farther from the central star than Neptune. Unless some additional ideas are considered, it is difficult for CA model to explain these planets because it is thought that the gas disappears from the disk before the formation of the massive solid core (Dodson-Rovinson et al., 2009). On the other hand, GI model has the possibility to make gas giant planets far from the central star if the mass of a disk is greater than that of the protostar. Thus, GI model attracts attention to explain their formation.

In order to understand planet formation by GI model, it is desirable to know the physical condition for disk fragmentation and the position where the fragmentation occurs. At present, two criteria are frequently used to discuss whether a protoplanetary disk is likely to fragment. The first is Toomre’s stability criterion (Toomre, 1964)

 Q≡csκepπGΣ>1, (1)

with gravitational constant , epicyclic frequency , sound speed , surface density , and stability parameter . Toomre (1964) showed that the infinitesimally thin disk is stable if the stability criterion (1) is satisfied. It is necessary to violate this criterion in order for fragmentation to occur. The other criterion is Gammie’s cooling criterion (Gammie, 2001)

 β≡tcoolΩ=u(dudt)−1Ω≲βc, (2)

where is the angular velocity, is the cooling time, is the specific internal energy, is the total cooling rate, is cooling parameter, and is the critical cooling parameter. Gammie (2001) suggested that rapid cooling is necessary for fragmentation, in addition to violating the Toomre criterion (1). Using shearing sheet simulations with constant cooling rate and without heating except for artificial viscosity, he showed that a self-gravitating disk can fragment only if . If not, disk cannot fragment because the disk can be stabilized by internal heating due to the turbulence arising from gravitational instability. This state is called ”gravitoturbulence”. It is said that the heating owing to the gravitoturbulent dissipation automatically controls the Q value close to unity. The gravitoturbulence has the potential to stabilize the system where Toomre criterion is violated, while it does not arise in the disk that satisfies Toomre criterion.

The obtained by Gammie (2001) is for the specific case (local, 2D, etc.). Many three-dimensional simulations were carried out in order to obtain the general value of (e.g., Rice et al. (2005); Clarke et al. (2007); Cossins et al. (2010)). However, each simulation suggested the different , and general has not been clarified yet. Furthermore, it is not obvious whether or not exists as some value. Meru & Bate (2011a) calculated the disks with different initial parameters such as the disk size, the surface density profile, the mass of the disk, and the mass of the central star. From their calculations, it is found that is different among the disks with a different set of initial parameters. This result shows that it is not determined only by whether the disk can fragment or not. Moreover, Meru & Bate (2011b) calculated the with various numerical resolutions. They showed that the disk was able to fragment with larger if the resolution is better, and they had no indication of the convergence of . From these results, it is not clear that Gammie criterion is always applicable to discuss the fragmentation of a disk.

Some analytical works for disk fragmentation have been done. Rafikov (2005) argued that a protoplanetary disk is unlikely to fragment near the protostar (AU) by using two criteria described above. However, as noted above, Gammie criterion is uncertain, and thus this argument includes the same uncertainty. On the other hand, Toomre criterion is probably assured in the sense that the disk violating Toomre criterion is gravitationally unstable. Other previous analytical works for disk fragmentation (e.g., Clarke (2009); Kratter & Murray-Clay (2011)) assumed gravitoturbulent disks with . These studies do not have the ability to discuss the position where Toomre criterion is satisfied. Since Toomre criterion provides robust necessary condition for fragmentation, the studies based solely on Toomre criterion (1) is favorable to know the possibility of fragmentation.

In this study, we investigate the gravitational stability in protoplanetary disks based solely on Toomre criterion without the uncertainty associated with Gammie criterion. We do not aim to consider the time evolution of particular disks but derive an instability condition that is available for various disks. Since recent observation shows the diversity of protoplanetary disks and planetary systems, we think that it is important to derive such a widely useful formulation. In order to calculate Toomre’s Q value with given and , it is necessary to know the temperature. Thus, it is significant to consider realistic thermal processes. We construct an analytical model for temperature calculation with thermal process in section 2. In section 3, we derive the critical surface density that is necessary for a disk to become gravitationally unstable. This is available to discuss the possibility and position of gravitational instability, regardless of surface density profile or formation process of the disks. In section 4, we introduce two semi-analytic models for protoplanetary disks, and discuss the possibility and location of the gravitational instability by comparing the surface densities of these models to . In section 5, our results are discussed and compared to previous studies. Finally, we summarize this study in section 6.

2 Temperature

In this section, an analytic model for the temperature calculation is introduced. The local thermal equilibrium between viscous heating and radiative cooling are assumed. We use Keplerian rotation law, ideal gas, vertical hydrostatic equilibrium, and a geometrically-thin disk. We assume that radiation escapes only in the vertical direction. Effective viscosity with constant introduced by Shakura & Sunyaev (1973) is adopted. Since these assumptions are similar to the standard accretion disk (Pringle, 1981), we have

 Ω = √GMs/r3, (3) c2s = γkBT/¯¯¯¯¯m, (4) H = cs/Ω, (5) ρ = Σ/2H=ΣΩ/2cs, (6) Γ = 98αcsHΣΩ2, (7) Λ = ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩32σT43κΣ(τ≥1)8σT4κΣ3(τ<1) , (8) and Γ = Λ, (9)

where is specific heat, is mean molecular weight, is density, is scale height, is the mass of protostar, is Stefan-Boltzmann constant, is the Boltzmann constant, is viscous heating rate, is opacity, is cooling rate per unit area, and is optical depth of vertical direction. We approximate cooling rate by connecting optically thin and thick limit continuously. Here, we use the Rosseland mean opacity that is approximated by using power-law as

 κ=κ0ρaTb, (10)

where the values of are summarized in table 2 (Bell & Lin (1994); Cossins et al. (2010)). This table shows the characteristic temperatures at which the main component of opacity changes. For example, is the temperature where the main component of opacity changes between ices and sublimation of ices, and is between sublimation of ices and dust grains. Other characteristic temperatures depend on the density.

{longtable}

llllll Bell&Lin opacity & & & & Temperature & Temperature
Opacity regime & & & & range From(K) & range To(K) \endfirsthead\endlastfootIces & & 0 & 2 & 0 & 166.8
Sublimation of ices & & 0 & -7 & 166.8 & 202.6
Metal dust & & 0 & 1/2 & 202.6 & 2286.7
Sublimation of metal dust & & 1 & -24 & 2286.7 & 2029.7
Molecules & & 2/3 & 3 & 2029.7 & 10000.0
Hydrogen scattering & & 1/3 & 10 & 10000.0 & 31195.2

In order to calculate the temperature as a function of and , seven variables can be eliminated by using the eight equations from (3) to (10). Thus, we can obtain the equilibrium temperature analytically for a given and . Additionally, we introduce the minimum temperature in order to mimic the effect of ambient heating sources in a molecular cloud. If the above equilibrium temperature becomes lower than , we simply use instead of the equilibrium temperature. As a result, the temperature can be represented as

 (11)

where

 A = αγkB√GMs¯¯¯¯¯mσ (12) and B = 12√¯¯¯¯¯mGMsγkB .

From equation (11), we can calculate temperature as a function of and for a given , , and . In this study, we set , which is the typical value of caused by MRI (magnet rotational instability) turbulence, and assume , which is the typical value of a molecular cloud.

Figure 1 shows the temperature in the plane for the case with , and as the fiducial case in this study. The solid lines show the contours of temperature, and the dotted line shows the line. The characteristic temperature defined above is labeled as 166K, and is labeled as 202K. The line for 166K can be regarded as snow line because ices become the main component of the opacity in the area below this line. Since other characteristic temperatures depend on the density, they do not coincide with contour lines of iso-temperature. At these characteristic temperatures, the property of opacity drastically changes as shown in equation (10) and table 2. This property of opacity has much influence on the gravitational stability of protoplanetary disks.

From figure 1, it is seen that the temperature is high for large and small . Above dotted line, the temperature is high for large surface density because a large amount of matter prevents radiation from escaping the system. In a small radius, the dynamical time scale is short, and thus the heating rate is large. This is why the temperature is high for large and small . The property of cooling rate drastically changes whether or not the condition is satisfied as described in equation (8). Below the dotted line, which means that the disk is optically thin, it is seen that the temperature is independent of surface density. This is because all radiation emitted from disk can escape in the optically thin case. At the right-bottom area of the line with 10K, the temperature becomes the minimum temperature due to the small viscous heating rate and large cooling rate. We henceforth name this area “isothermal region”. Note that the dotted line exists only at low temperature area, such that . This is due to the high cooling efficiency for the disk with .

3 Critical Surface Density

Now we can calculate Toomre’s by using equations (1) and (11). The analytical expression for Q value is summarized in appendix A. In figure 2, the solid line represents the critical surface density that is necessary to satisfy for a given . As references, the dotted and dashed lines represent the surface densities that satisfy and , respectively. A disk is expected to become gravitationally unstable if the disk has a larger surface density than . This critical surface density is determined only by local thermal process and criterion for gravitational instability. Thus, one can discuss the possibility and position of gravitational instability regardless of disk formation process by comparing a given surface density profile to . We show examples of this discussion in the next section. However, one should keep in mind some cautions, for example, 1) gravitational instability is not the same as fragmentation and 2) Q value is not universal for the case with non-uniform medium or with non-axisymmetric situations.

In figure 2, changes discontinuously by one order of magnitude at . We henceforth represent this critical radius as . By comparing figure 1 and figure 2, it is confirmed that is large enough to satisfy in because of the large Coriolis force and the high temperature. Analytic formula for is derived from equations (1) and part of (11) as

 Σc=(27256κ0ABaC6−2b+ar−3(7+2a−2b)/2)1/(4−2b)   (τ>1) , (13)

where , , and are defined in equations (12), (2), and

 C=√γkBMsπ2G¯¯¯¯¯m , (14)

respectively. With the value of , , and in table 2, equation (13) provides us the property of (see appendix A for details). In the fiducial case, in , sublimation of metal dust mainly contributes to the opacity, and . In , the opacity is dominated by metal dust, and . In , the main component of opacity is sublimation of ices, and . These properties can be seen in figure 2 from the fact that the line breaks at the place where the main component of opacity changes. In , is given by the isothermal condition and as

 Σc=cs,minΩπG=21 g/cm2 (r100AU)−3/2(Tmin10K)1/2(MsM⊙)1/2 (15)

where is the sound speed for .

In figure 2, the critical surface density changes discontinuously at . This is because the value is independent of the surface density there. The reason is explained as follows. If the surface density increases, the effect of self-gravity increases to destabilize a disk. On the other hand, in the optically thick case, the surface density also has an ability to block radiative cooling. Then, the effect of thermal pressure increases to stabilize a disk as the surface density increases. If these two effects are balanced, the value of Toomre Q becomes independent of the surface density. This balance is realized when the opacity depends on temperature as in optically thick regime (see appendix A). The conditions to realize this balance are satisfied in the shaded area in figure 1 where ices become the main component of opacity in optically thick regime. Thus, the discontinuously changes at the critical radius where ices form. The analytical formula for critical radius is derived by equations (1) and part of (11) with ices opacity as

 (16)

In equation (16), is an increasing function of and . This is because the viscous heating rate is large with large and large . However, dependence on and is weak. Thus, varies only severalfold even if or is different from typical value by an order of magnitude.

In figure 2, it is seen that the become multivalued function at AU and . In this area, molecules are the main component of opacity, and there is the inner most radius of the unstable region around this area. We explain this inner most radius in appendix B.

Note that the critical surface density is derived with the constant viscous parameter . It is considered that this assumption corresponds to the case that the viscosity originates from not the self-gravity but the MRI turbulence. This case assumes that the disk is MRI active everywhere. Thus, this result is specific to the case with the disk that has no region where the MRI is inactive (generally called ”dead zones”). In other words, it is assumed that there is the significant heating of non-gravitational origin at all radii in the disk. On the other hand, the previous studies like Clarke (2009) use the viscous parameter determined by the self-regulated gravitoturbulence with the condition . In this case, the disk is assumed to have the large dead zones without effective viscosity owing to the MRI turbulence. The realistic protoplanetary disks must be bracketed by these two extreme cases. Considering such disks is beyond the scope of this paper and is left as the future work.

4 Instability Condition for Protoplanetary Disks

In this section, the condition for gravitational instability is discussed by comparing the surface density of particular disks to . We introduce two simple semi-analytic models for protoplanetary disks. One is the steady state accretion disk (the steady model). This model is realized if angular momentum is sufficiently transported. The other is the disk that has the same angular momentum distribution as the cloud core before it collapses. In this model, we assume the conservation of angular momentum distribution, and hence this model is hereafter called the AMC (Angular Momentum Conservation) model. The AMC model is realized if angular momentum is not sufficiently transported. We study the steady model in section 4.1. The AMC model is investigated in section 4.2. The interpretation of these two models is discussed in section 4.3.

4.1 The Steady State Accretion Disk Model

First, we consider the model of the steady state accretion disk. The surface density profile is calculated in the framework of the standard disk model with viscosity (Shakura & Sunyaev (1973); Pringle (1981)). The disk structure are determined by three parameters; viscous parameter , protostar mass , and mass accretion rate . By the equation of continuity and angular momentum conservation, we have the relation

 αcsHΣ=˙M3π. (17)

This equation shows that surface density is determined by mass accretion rate . We calculate the temperature in the same manner as in section 2 using the cooling function approximated as equation (8) and opacity represented by equation (10) with table 2, and introducing . Using equations (3), (4), (5), (11), and (17), we obtain the analytic profile for the surface density in this model as

 ΣSD=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩DX(27κ0ABa256r−3(a+1)/2)−15+1.5a−br−3X/2(τ≥1)DX(27A256κ0Bar−3(1−a)/2)−13+b−1.5ar−3X/2(τ<1)˙M√GMs3παc2s,minr−3/2(T=Tmin), (18)

where we use a symbols , and are defined in equations (12), (2),

 D = ¯¯¯¯¯m˙M3παγkB√GMs ,

and

 X = ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩3+a/2−b5+1.5a−b(τ≥1)3+b−a/23+b−1.5a(τ<1),

respectively. This steady state accretion disk does not have the inner and the outer radius. In other words, the disk is infinitely extended.

From equation (18) and table 2, we can calculate the . Figure 3 shows the surface densities for the case with (long-dashed line), (short-dashed line), and (dotted line) for the fiducial case defined in section 2. In figure 3, it is seen that the lines of break at the radii where the main component of opacity changes just like the line in figure 2. It is also seen that is large with large . The critical surface density derived in section 3 is superimposed in figure 3 by a solid line. The disk is unstable if the condition is satisfied. From figure 3, it is found that the disks with and are unstable in , where is the critical radius defined in equation (16). On the other hand, the case with is expected to be stable in all radius. It is also seen that all the lines including the and has the same dependence on radius in the outer region (AU). From equations (15) and (18) for , we can confirm that and has the same dependence on radius owing to isothermal state. By using this property, we can recast the instability condition as

 ˙M≥˙Mcrit≡3αc3s,minG=7.7×10−8M⊙/yr(α0.01)(Tmin10K)3/2. (20)

Note that this instability condition (20) is independent of protostar mass . It is also found that is increasing function of and . As shown in figure 3, the minimum radius of the unstable region defined as equals to the critical radius . The critical radius depends on the protostar mass as , and is independent of the mass accretion rate. Thus, in the steady model, mass accretion rate determines whether or not the disk becomes unstable, and protostar mass affects where the disk becomes unstable. This property described above is summarized in figure 4. Figure 4 shows the contours of the minimum radius of the unstable region in the plane. It is seen that the stable disk is present in the left area where mass accretion rate is low and that the minimum radius is large with large protostar mass. It is also seen that the disk becomes unstable at smaller radius than for the large mass accretion rate . This property appears where the surface density is large enough to satisfy in . From figure 3, it is found that is required to be larger than in order to become gravitationally unstable in .

4.2 The AMC Model

Next, we consider the disk that has the same angular momentum distribution with the cloud core before it collapses. This model corresponds to the disk that has just formed.

Here, we introduce a simplified model for the formation of a disk and a protostar. Figure 5 schematically shows the formation process of the disk and the protostar in this model. It is assumed that a rotating cloud core collapses to make a protostar and a disk. We divide the cloud core into two parts. The shaded part near the rotation axis is assumed to become the protostar. The disk is assumed to be made from the other part. The boundary radius between these two parts is determined so that the mass inside the shaded part equals to a given protostar mass as a parameter. The cloud core is assumed to rotate rigidly. The angular velocity of the cloud core is determined by introducing the rotation parameter , where is the rotation energy and is the gravitational energy. Note that the disk formed in this model has the inner radius and the outer radius , different from the steady model described in section 4.1. It is assumed that the protostar and the disk forms instantaneously and that the disk has the Keplerian rotation velocity with ignoring the self-gravity. We also assume the inviscid and axisymmetric formation of the disk. In other words, it is assumed that the relation between the mass and the angular momentum is kept without redistribution. As the former cloud core, two types of the density distributions are considered with the assumption of spherical symmetry. One is the Bonnor-Ebert sphere (Ebert (1955); Bonnor (1956)), and the other is a sphere of uniform density. In this study, we consider both cases, and the result for the case with Bonnor-Ebert sphere is mainly described. We discuss the difference of these two cases later.

The Bonnor-Ebert sphere is an isothermal sphere in which the thermal pressure balances the self-gravity and the external pressure. This sphere is identified by three quantities, the cloud radius , the sound speed , and the central density . It is known that a Bonnor-Ebert sphere is unstable when the condition

 RE≥Rcrit≡6.46cs,0√4πGρc (21)

is satisfied. The Bonnor-Ebert sphere whose radius equals to is called the critical Bonnor-Ebert sphere, and we use this critical Bonnor-Ebert sphere in this model. We set , which corresponds to the typical sound speed in a molecular cloud. The central density is fixed as . With these parameters, the critical Bonnor-Ebert sphere has the mass and the radius AU. Here, we approximate the density profile of the Bonnor-Ebert sphere (Tomide 2011, private communication) as

 ρ(R)=ρc(1+R2/R2c)3/2 , (22)

where is the radius in the spherical coordinate and . The equation (22) approximates the density profile of a Bonnor-Ebert sphere to the accuracy of the error within a few percent. This approximation enables us to calculate the surface density of the Bonnor-Ebert sphere analytically as

 ΣBE(rcy)=2∫(R2E−r2cy)1/20ρ(R)dz=2ρc1+r2cy/R2c ⎷R2E−r2cy1+R2E/R2c, (23)

where is the radius in the cylindrical coordinate.

In this model, the surface density of the disk is determined from the angular momentum distribution of the former cloud core as the follows. Once the rotation parameter is given, the angular momentum distribution of the cloud core is determined. Next, give the ratio of disk to protostar mass as a parameter, and the specific angular momentum distribution of the disk is determined. Suppose that a fluid element in the cloud core at cylindrical radius falls onto the disk at cylindrical radius . The disk is supported by centrifugal force and rotates Keplerian velocity at . By the conservation of angular momentum, the relation between and is found as

 rcy=(GMsrd/Ω20)1/4 . (24)

By using the relation (24) with the relation of mass conservation

 rdΣAMC(rd)drd=rcyΣBE(rcy)drcy , (25)

we derive the surface density of the disk as

 ΣAMC(rd)=12ρcRE1+r2cy/R2c(1−r2cy/R2E1+R2E/R2c)1/2r2cyr2d . (26)

Note that the depends on as in the equation (24). From equation (26), we can find that depends on radius as in the region where . It is also found that depends on radius as for . Thus, the enclosed mass in the disk does not diverge. In this study, we concentrate on the Keplerian rotating disk without self-gravity. In order to use this assumption with validity, the mass of the disk should be less than that of the protostar. If we request that the disk rotates Keplerian velocity, the boundary radius shown in figure 5 is always larger than . Thus, the depends on radius as in this paper. By the way, realistically, it is considered that a protostar does not have such a large mass just after the formation of the protoplanetary disk. In this model, we consider only the outer region of the disk and regard the inner region of the disk as a part of the protostar, at least, in the sense of gravitational field. The boundary radius between the inner region and the outer region of the disk is determined by the equation (24) and .

Now, we can calculate the surface density of the disk for given and . First, we look over the property of for the case with fixed . Figure 6 shows the surface densities for the cases with 0.005 (long-dashed), 0.001 (short-dashed), and 0.0002 (dotted) with . In the case with large , it is seen that is small and that the disk is present far from the protostar. The solid line in figure 6 represents the critical surface density derived in section 3. In the case with , the disk is stable because the stability condition is satisfied owing to large in AU, where is the critical radius defined in equation (16) with . In the case with , the stability condition is violated in AU, and thus the disk is unstable there. This is because changes discontinuously at . In the case with , the angular momentum is so large that the inner radius AU is larger than the critical radius AU. The disk can be unstable around the inner radius in AU. Note that, in the cases with and 0.005, the outer regions of the disks are stable because rapidly decreases as radius increases.

Next, we calculate for various values of as well as . Figure 7 summarizes the results in - plane. The symbols with “+” correspond to the cases in figure 6. The results are classified into five groups in the plane. The area labeled as “stable” in figure 7 indicates that the disks with these parameters are stable in all radius. The case with in figure 6 belongs to this area. These disks have small outer radius owing to the small angular momentum. Thus, the condition is satisfied due to the large in . In figure 7, the areas labeled as “a few AU”, “r”, and “30AU” indicate the unstable disks. The area labeled as “a few AU” exists in the small and large area. The disks with these parameters have small radius and large surface density enough to satisfy in . Thus, they have a possibility of fragmentation at a few AU from the protostar. The area labeled as “r” exists around middle area in figure 7. This area includes the case with in figure 6. Although the disks with this area are stable inside the critical radius due to small , they are unstable in . This is because drastically changes at . It is the formation of ices that enables the disks to become unstable. The area labeled as “30AU” exists at large area. This area includes the case with in figure 6. In this area, the inner radius of the disk is larger than the critical radius , and the disk is unstable around the inner radius. The area labeled as “inner” in figure 7 is difficult to conclude. In the present AMC model, the disks with this area are regarded as stable because is very small due to the too large inner radius AU. However, this large is artificial because the AMC model is too simplified in order to calculate analytically. Realistically, it is expected that the disk is present in . If the disk has the surface density represented as in equation (26), the disks with the area labeled as “inner” in figure 7 have a possibility to be unstable.

From above results, we found that the disk is expected to be unstable with large in the AMC model. Here, we estimate the critical rotation parameter that is necessary to become unstable for the disk with small mass (). In this model, the disk becomes unstable if the disk extends to the radius . The outer radius of the disk can be estimated by substituting in equation (24). On the other hand, the is represented in equation (16). By using the relation , we derive the critical rotation parameter as

 β0,c=4×10−4(rc24AU)(RE104AU)4(Mc1M⊙). (27)

Note that the depends on the viscous parameter and on the protostar mass and that the and depend on the central density and sound speed . In figure 7, it is seen that the disk becomes unstable for the even with small mass ratio . The minimum ratio of disk to protostar mass satisfying the unstable condition is with , which belongs to “r” area.

4.3 Interpretation of the Steady Model and the AMC Model

In previous two subsections, we studied the gravitational stability of protoplanetary disks by using two simple semi-analytic models. Here, we interpret these models by considering the formation scenario of a protoplanetary disk.

It is believed that a cloud core with angular momentum collapses to form a protostar and a protoplanetary disk. First, the central part of the cloud core begins to contract and makes a protostar and a small (typically a few AU) disk (e.g. Bate (1998)). Next, the outer envelope of the core mainly falls onto the disk, and the mass accretion from the envelope makes the disk grow up in mass and size as typically larger than 100 AU (e.g. Adams et al. (1988)). The time scale for this disk growth is estimated as according to the self-similar solution for the collapse of the rotating isothermal cloud (Saigo & Hanawa, 1998). This disk that has just formed corresponds to the AMC model introduced in section 4.2. Later, this disk experiences the viscous evolution with the transport of angular momentum and eventually forgets the initial distribution of the angular momentum. The time scale of viscous evolution is estimated as from the azimuthal component of the momentum equation. It is expected that a steady state disk, which is discussed in section 4.1, corresponds to the state after the viscous evolution. For the case with and , the time scales for disk growth and for viscous evolution are estimated as

 tvis∼1.6×103yr(α0.01)−1(r1AU)3/2 (28)

and

 tgrow∼25yr(cs,0190m/s)−1(r1AU)1. (29)

From equations (28) and (29), it is seen that the disk growth time is much shorter than the viscous evolution time for AU. Thus, it is expected that the viscous redistribution of angular momentum can be regarded as negligible during much longer time than the dynamical time scale after the formation of the disk. From equation (28), it is found that the viscous evolution time is short for a small radius. Hence, it is likely that the steady state accretion disk forms at inner region first, and it will spread outward in viscous time scale. In summary, we interpret the AMC model as the outer part of the young disk and the steady model as the old disk.

5 Discussion

5.1 Thermal Stability of Equilibrium State

We calculated the equilibrium temperature based on the approximated opacity given in equation (10) and table 2 (Bell & Lin, 1994) in section 2. This opacity depends on the temperature, and the property of the opacity varies when the main component of opacity changes. Thus, the cooling rate also depends on the temperature similarly. If the temperature dependence of the cooling rate is weaker than that of the heating rate, this state is thermally unstable (Field, 1965). In order to check the stability of the thermal equilibrium state we used, we compare the exponent of the temperature in the heating rate to the exponent in the cooling rate. In our model, as shown in section 2, the heating rate has the relation , regardless of the optical depth.

In optically thick case, the temperature dependence of the cooling rate is . Thus, the condition for thermal stability is represented as . This condition is satisfied except for the case that hydrogen scattering is the main component of the opacity. This case is realized where the temperature are larger than the critical temperature K for . Within the parameter range we are interested in, the equilibrium temperature is safely smaller than , and thus, in optically thick case, the equilibrium state is stable.

In optically thin case, the temperature dependence of the cooling rate is . The condition for thermal stability is obtained as . From the values in table 2, this stability condition is violated when the main component of the opacity is sublimation of ices () or sublimation of metal dust (). According to results in section 2, in the optically thin case, the equilibrium temperature is sufficiently smaller than above which sublimation of ices becomes the main component of the opacity. Thus, in optically thin case, the equilibrium state satisfies the condition of thermal stability. Therefore, the equilibrium state we used in this study is thermally stable.

5.2 Comparison to Radiation Hydrodynamical Simulations

Here, we compare the results of our model to the results of some radiation hydrodynamical (RHD) simulations. First example is the simulation by Boley et al. (2006). They calculated the disk evolution about 4000 yr, and the surface density of the disk at the final state was approximated as , where . This disk did not fragment in their simulation. Our model also predicts that this disk cannot fragment because this surface density satisfies the stable condition . Thus, our model is consistent with the result of their simulation.

Next, we make the comparison to the result by Meru & Bate (2010). In their study, the evolutions of disks were calculated for the case with different values of the opacities. Their result showed that the opacity that is smaller than the interstellar values promotes the fragmentation. In our model, the values of and represented in equations (13) and (16) are small with the small opacity. This means that the disk easily becomes unstable with small opacity. In this sense, the prediction by our model is consistent with them.

Finally, we compare the result of the simulation by Stamatellos et al. (2011). Their calculation was performed with nine initial conditions, and they obtained the result that only three of nine disks fragmented. From this result, they suggested that the fragmentation is unlikely to occur when the disk satisfies the condition that AU and that 0.36, where is the outer radius of the disk and is the ratio of disk to protostar mass. Our result is consistent with their result in the sense that fragmentation is difficult near the central star (AU). However, with their initial surface density profile, our model predicts that all of nine disks used in Stamatellos et al. (2011) become gravitationally unstable for the case with . These differences between their results and our predictions may be caused by the heating due to the gravitoturbulence. In this sense, we should notice that the instability condition investigated in this paper is not exactly the same as the fragmentation condition.

In section 4.1, we derived the critical mass accretion rate that is necessary in order that the steady disk becomes unstable. Hartmann et al. (1998) derived the mass accretion rate of T-Tauri stars in Taurus and Chamaeleon I, based on the optical observation. The typical value of the accretion rate is as small as . Since this value is less than the critical mass accretion rate given in equation (20), the disks around a T-Tauri stars are expected to be stable. However, in Hartmann et al. (1998), 7 of 40 T-Tauri stars in Taurus and 1 of 16 in Chamaeleon I are estimated to have a large mass accretion rate enough to become gravitationally unstable. Thus, these statistics indicate that about 18 % of the T-Tauri stars have the possibility to have an unstable disk in Taurus, whereas the fraction is less than 7 % in Chamaeleon I. By the way, around a protostar, the Keplerian disk has not been clearly observed yet. If such a disk is present, its mass accretion rate is estimated as

 ˙M∼MJtff∼c3sG∼2×10−6M⊙/yr(T10K)3/2, (30)

where is Jeans mass, and is free-fall time. Because this mass accretion rate is much larger than , the disk around a protostar is expected to be gravitationally unstable. In summary, it is expected that most of the disk around a T-Tauri star is stable and that the disk around a protostar is likely to fragment. Anyway, the relation between mass accretion rate and gravitational instability may be a useful tool to probe the fragmentation in the disk systems observationally.

Clarke (2009) (hereafter C09) also discussed the fragmentation of a steady state accretion disk. Assumptions set in C09 are similar to our model, but there are critical differences. C09 assumes the condition of gravitoturbulence , uses spatially-variable parameter, and adopts Gammie criterion as fragmentation criterion. On the other hand, we assume not to be in gravitoturbulent state, use constant viscous parameter , and adopt Toomre criterion. These differences produce qualitatively different results. For example, in C09, fragmentation is predicted even with very small mass accretion rate (), whereas, in our model, the disk is expected to be stable with such a small mass accretion rate. However, note that both studies imply that the formation of ices is an important process. In the regime of ices opacity , the temperature dependence of cooling rate becomes , which is weaker than that of other opacities (sublimation of ices, metal dust, and sublimation of metal dust). This property enables the disk to become low temperature when the heating rate becomes small. The importance of this fact does not change regardless of fragmentation criteria.

5.4 The AMC Model

In the AMC model introduced in section 4.2, the disk is unstable when the condition is satisfied. By the way, according to the results of radiation hydrodynamical calculations, the disk made from the core with the large rotation parameter is expected to fragment before the formation of the protostar (Bate (2011)). From the observation of the line emission of NH and NH, the rotation parameter of the cloud core is estimated to range in with the typical value (Caselli et al., 2002). Only the 3 of 20 cores have the rotation parameter . Thus, only 15 % of the cores are expected to make a stable disk. The cores with middle rotation parameter have the possibility to fragment after the formation of a protostar. The 8 of 20 cores is present within this parameter range. The number of the core with is the 9 of 20 cores. The cases with such a large are outside the scope of our AMC model.

Matsumoto & Hanawa (2003), Machida, Inutsuka, Matsumoto (2010), and Tsukamoto & Machida (2011) calculated the formation process of the disk from the cloud core using barotropic equation of state for various values of the rotation parameter . These studies showed that the disks become gravitationally unstable for the cases with large . This feature is consistent with the result of our AMC model. Matsumoto & Hanawa (2003) calculated the cloud collapse with the rotation parameter . They showed that the disk made from the cloud with become gravitationally unstable but not fragmenting until the final state of their calculation. Our AMC model also predicts that the disk with is unstable. They also showed that fragmentation occurred before the formation of the protostar for the case with . These cases are outside the scope of our AMC model. The result in Machida, Inutsuka, Matsumoto (2010) is the quantitatively different from our model. For example, the value in their calculation is smaller than that in the AMC model by one order of magnitude. However, it is difficult to compare the result of the AMC model to that of their calculation quantitatively because, in their calculation, fragmentation tends to occur during the phase that the mass of protostar is smaller than that of the disk. Tsukamoto & Machida (2011) is more consistent with our AMC model. Despite that they used the barotropic equation of state, which is different from us, the value in their calculation is consistent with our AMC model within the error of factor 3.

5.5 The effects of ignored processes

In this study, we ignored some physical processes that are expected to affect the gravitational instability in the disk. Here, we comment on the effects of processes ignored in this paper.

First, the self-gravity will affect the disk properties through the effect on the transport of angular momentum and energy dissipation. In the form of gravitational torque due to the non-axisymmetric mode and the gravitoturbulence, the self-gravity affects the viscous parameter . It is discussed that the becomes large if the value approaches unity (Kratter et al., 2008). However, this as a function of disk properties has not been clarified yet. In this study, in order to avoid this uncertainty, we used constant parameter and derived the condition for gravitational instability as a function of . The work considering from the self-gravity remains as a future work. The self-gravity will also affect the scale height . If considering the self-gravity in the vertical direction, the scale height with self-gravity is represented as a function of value. This for the case with is approximately half of the without self-gravity represented in equation (5). From this effect, the values of and are modified, but the amount of this modification is at most a few tens of percent. Thus, it is considered that our result does not qualitatively change even if considering the compression by self-gravity in the vertical direction.

Next, the effect of irradiation from the protostar will affect the disk properties. At the region far from the protostar, the heating by irradiation from the protostar becomes greater than the viscous heating because the viscous heating strongly decreases as the radius increases. Thus, at the far from the protostar, there is a possibility that the disk is stabilized by the irradiation heating. The heating rate by irradiation is affected by the disk flaring, and this flaring is determined by the temperature distribution of the disk. It is difficult to make the accurate treatment of irradiation heating in our model because we did not consider the radial distribution of the disk when deriving the instability condition in section 3. However, the effect of irradiation heating is weak in disks with large surface density. Thus, it is considered that our results for optically thick case will not critically change. This irradiation heating should be taken into account when investigating the properties of the particular disk.

Finally, we comment on the effect of magnetic field. The magnetic field has much influence on the angular momentum transport during all phases from the protostar formation to the disk evolution. In the phase of the protostar formation, it is considered that the magnetic field triggers the outflow and that it transports the mass and the angular momentum (Shu et al. (1994); Tomisaka (2002); Machida et al. (2008)). Thus, it is expected that the assumption used in the AMC model is violated by the effect of this outflow. As a future work, it remains to construct the model including the effect of the outflow. In the phase of the disk evolution, the turbulence arising from the magnet rotational instability (MRI) is considered to become the main source of viscosity. In the region where ionization degree is low (called as dead zone), the MRI turbulence is inactive (Sano et al., 2000). Thus, in the dead zone, the viscous parameter is expected to be much smaller than that in the MRI active region. Although we used a constant model in this study, as the future work, it will be important to consider the realistic viscous parameter based on the physical processes with the relations among other physical quantities.

6 Conclusions

In this study, the gravitational stability of protoplanetary disks is studied. The temperature in a protoplanetary disk was analytically calculated with considering thermal effects such as radiative cooling, viscous heating, and ambient heating source in a molecular cloud. It is found that the temperature becomes as low as 20K in optically thin regime. We derived a critical surface density , which is necessary for the disk to become unstable, as a function of radius. By comparing a given surface density to , the possibility and the position of the gravitational instability can be predicted. The formation of ices is important for the gravitational instability, and we found that the changes discontinuously at the critical radius 20AU, where ices form.

Two semi-analytic models of protoplanetary disks are used with above critical surface density in order to discuss the possibility and the position of gravitational instability. In the model of a steady state accretion disk, which corresponds to the evolved disk, it is found that the disk becomes unstable when the mass accretion rate is greater than the critical mass accretion rate in equation (20). From this result, it is expected that most of the disks around T-Tauri stars are gravitationally stable and that the disks around protostar are unstable in . In the model of the disk with the same angular momentum distribution as the former cloud core, which corresponds to the outer region of the young disk, the disk is expected to become unstable in for the case with large rotation parameter (). From the results of these two models, we conclude that the fragmentation near the central star () is expected to be rare to occur as long as the parameter range indicated by observations is considered.

We thank Fumio Takahara, Yutaka Fujita, and Hideyuki Tagoshi for fruitful discussion and continuous encouragement. We are also grateful to Taishi Nakamoto and Kei Tanaka who provided helpful comments and suggestions.

Appendix A The expression of Temperature, Q value, and Critical Surface Density

Here we summarize the value and dependence of the temperature, Toomre Q value, and the critical surface density .

a.1 Equilibrium Temperature

The equilibrium temperature given in equation (11) is explicitly written below. The temerature range for each formula is given in table 2

a.1.1 Optically Thick Case

(a)ices

 T=12K(r10AU)−3/2(Σ102g/cm2)2(Ms1M⊙)1/2(α0.01)1 (31)

(b)sublimation of ices

 T=1.8×102K(r1AU)−3/20(Σ102g/cm2)1/5(Ms1M⊙)1/20(α0.01)1/10 (32)

(c)metal dust

 T=8.0×102K(r1AU)−3/5(Σ103g/cm2)4/5(Ms1M⊙)1/5(α0.01)2/5 (33)

(d)sublimation of metal dust

 T=1.2×103K(r1AU)−6/55(Σ104g/cm2)6/55(Ms1M⊙)2/55(α0.01)2/55 (34)

a.1.2 Optically Thin Case

As described in section 2, the temperature become very low in optically thin case. Thus, we consider only the case with ices opacity.

(a)ices

 T=21K(r1AU)−3/10(Ms1M⊙)1/10(α0.01)1/5 (35)

a.2 The Q Value

The Q value defined in equation (1) is represented below.

a.2.1 Isothermal Region

 Q=Ωcs,minπGΣ=21(Σ1g/cm2)−1(r100AU)−3/2 (36)

a.2.2 Optically Thick Case

In optically thick case, Q value is represented as

 Q=C(27256κ0ABa)16+a−2br−327+2a−2b6+a−2bΣ2b−46+a−2b, (37)

where A, B, and C is defined in equations (12), (2), and (14), respectively. Using the value of table 2, we can see the dependence of Q value. From equation (37), it is found that the Q value is independent of surface density for the case with , which is realized when ices dominate the opacity.

a.2.3 Optically Thin Case

In optically thin case, Q value is represented as

 Q=C(27A64κ0Ba)16−a+2br322a−2b−76−a+2bΣ−6−2b6−a+2b, (38)

where A, B, and C is defined in equations (12), (2), and (14), respectively.

a.3 The Critical Surface Density Σc and the Critical Radius rc

(a)isothermal region

 Σc=cs,minΩπG=21 g/cm2 (r100AU)−3/2(Tmin10K)1/2(MsM⊙)1/2 (39)

This is available for , where is the critical radius represented as equation (16).

(b)ices

In this regime, unique is not present but is present.

 rc=(γkBM3/4sπ¯¯¯¯¯mG1/4√27κ0α256σ)4/9=24AU(α0.01)2/9(MsM⊙)1/3 (40)

(c)sublimation of ices

 Σc=3.4×103g/cm2(r10AU)−7/4(Ms1M⊙)7/12(α0.01)1/18 (41)

This is available in , where is the radius where the main component of opacity changes between sublimation of ices and ices. This is represented as

 r1=16AU(Ms1M⊙)1/3(α0.01)2/9 (42)

(d)metal dust

 Σc=6.2×103g/cm2(r10AU)−3(Ms1M⊙)1(α0.01)1/3 (43)

This is available in , where is the radius where the main component of opacity changes between metal dust and sublimation of ices. This is represented as

 r2=10AU(Ms1M⊙)48/141(α0.01)98/423 (44)

(e)sublimation of metal dust

 Σc=6.3×103g/cm2(r10AU)−171/104(Ms1M⊙)7/13(α0.01)1/52 (45)

This is available in , where is the inner most radius for the gravitational instability represented as

 redge=1.2AU(MsM⊙)1/3(α0.01)54/353 (46)

We explain this in Apendix B

Appendix B The inner most raius for the gravitational instability

In the regime where molecules mainly contribute to the opacity, which is realized at K, the temperature dependence of cooling rate becomes weak as . Then, the temperature becomes too high with large surface density, and the stabilizing effect of thermal pressure becomes larger than the destabilizing effect of surface density increasing. It means that increasing surface density stabilizes the disk. Thus, the line of become multi-valued function of radius. In figure 2, we see this feature for and AU. The upper line of slopes upwards when going from left to right. This is due to the lower heating rate at larger . Because the line slopes downwards where the main component of opacity is sublimation of dust grains, the unstable region has the inner edge at the point where the main component of opacity changes between molcules and sublimation of dust grains. In figure 2, we see that this inner edge locates at and . We name this radius . For , t