ROTATIONAL INSTABILITY IN THE OUTER REGION OF PROTOPLANETARY DISKS

We consider an axisymmetric rotating disk with hydrostatic equilibrium in the vertical direction. The equation for radial force balance in cylindrical coordinates is

$$\begin{eqnarray} \frac{j^2}{R^3} &=& \frac{GM}{R^2}+\frac{1}{\rho } \frac{\partial P}{\partial R}, \nonumber \\\ms &=& \frac{GM}{R^2}+\frac{H}{\Sigma } \frac{\partial \left(c^2_s \Sigma /H\right)}{\partial R}, \nonumber \\\ms &=& \frac{GM}{R^2}+\frac{1}{\Sigma } \frac{\partial \left(c^2_s \Sigma\right)}{\partial R}- \frac{c^2_s}{H} \frac{\partial H}{\partial R}, \end{eqnarray} \tag{ 2 }$$

where G is the gravitational constant, M the mass of the central star, P the pressure, ρ the density at midplane, Σ the surface density, H ≡ c_s/Ω_K the scale height, c_s the sound speed of the gas, and Ω_K ≡ (GM/R³)^1/2 the Keplerian angular velocity. We adopt the isothermal equation of state, $P=c_s^2\rho$. We assume that the disk is geometrically thin H/R ≪ 1, and that the vertical structure is isothermal, i.e., $\Sigma =\sqrt{2 \pi } \rho H$.

For generalization, we define normalized non-dimensional parameters for the disk radius, r, and the surface density, σ, as follows,

$$\begin{equation} R = R_0 r \ \ {\rm and} \ \ \Sigma = \Sigma _0 \sigma, \end{equation} \tag{ 3 }$$

where R₀ and Σ₀ are the normalization constants of the disk radius and the surface density, respectively. Also, we assume that the radial temperature profile has a power-law distribution with an index, β,

$$\begin{equation} T=T_0 r^{-\beta }. \end{equation} \tag{ 4 }$$

In this case, the radial dependences of the sound speed and the scale height become

$$\begin{equation} c_s = c_0 r^{-\beta /2} \ \ {\rm and} \ \ H = H_0 r^{(3-\beta)/2}, \end{equation} \tag{ 5 }$$

where T₀, c₀, and $H_0=c_0/({GM}/R_0^3)^{1/2}$ are the temperature, the sound speed, and the scale height, respectively, at R = R₀.

Substituting Equations (3) and (5) into Equation (2), we obtain

$$\begin{equation} \frac{j^2}{{GMR}_0}=r+\left(\frac{H_0}{R_0}\right)^2 r^{(2-\beta)} \left\lbrace \frac{\partial (\mathrm{ln }\; \sigma)}{\partial (\mathrm{ln }\; r)} -\frac{3+\beta }{2} \right\rbrace. \end{equation} \tag{ 6 }$$

Substituting Equation (6) into Equation (1), we obtain

$$\begin{equation} (2-\beta) \frac{\partial (\mathrm{ln }\; \sigma)}{\partial (\mathrm{ln }\; r)} +\frac{\partial ^2 (\mathrm{ln }\; \sigma)}{\partial (\mathrm{ln }\; r)^2} +\left(\frac{R_0}{H_0} \right)^{2}r^{(\beta -1)}-\frac{1}{2}(3+\beta)(2-\beta) >0, \end{equation} \tag{ 7 }$$

where we use r > 0. The disk is rotationally stable when the surface density profile satisfies this condition (see also Yang & Menou 2010).

When the disk is marginally stable for rotational instability, the specific angular momentum is constant, ∂j/∂R = 0. From Equation (7), this condition is equivalent to

$$\begin{equation} (2-\beta) \frac{\partial (\mathrm{ln }\; \sigma)}{\partial (\mathrm{ln }\; r)} +\frac{\partial ^2 (\mathrm{ln }\; \sigma)}{\partial (\mathrm{ln }\; r)^2} +\left(\frac{R_0}{H_0} \right)^{2}r^{(\beta -1)}-\frac{1}{2}(3+\beta)(2-\beta) = 0. \end{equation} \tag{ 8 }$$

The solution of Equation (8) is σ_ms,

$$\begin{eqnarray} \sigma _{{\rm ms}} &=& \exp \left[\frac{C_1}{(\beta -2)}r^{(\beta -2)} -\frac{1}{\beta -1}\left(\frac{R_0}{H_0} \right)^{2}r^{(\beta -1)}\right. \nonumber\\ &&\left. +\frac{3+\beta }{2}\mathrm{ln }\; r +C_2 \right], \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \frac{\partial (\mathrm{ln }\; \sigma _{{\rm ms}})}{\partial (\mathrm{ln }\; r)}&=& C_1r^{(\beta -2)}-\left(\frac{R_0}{H_0} \right)^{2}r^{(\beta -1)}+\frac{3+\beta }{2}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \frac{\partial ^2 (\mathrm{ln }\; \sigma _{{\rm ms}})}{\partial (\mathrm{ln }\; r)^2} &=& (\beta -2) C_1 r^{(\beta -2)}-(\beta -1) \left(\frac{R_0}{H_0} \right)^{2} r^{(\beta -1)},\nonumber\\ \end{eqnarray} \tag{ 11 }$$

where C₁ and C₂ are constants of integration. Equation (9) represents the marginally stable surface density profile for rotational instability.

Here, we examine the critical radius where rotational instability occurs in the similarity solution of the LBP74 model. The surface density evolution of viscous disks rotating with the Keplerian velocity is given by (e.g., Pringle 1981)

$$\begin{equation} \frac{\partial \Sigma }{\partial t} = \frac{3}{R}\frac{\partial }{\partial R} \left[R^{1/2} \frac{\partial }{\partial R} (R^{1/2} \nu \Sigma) \right], \end{equation} \tag{ 12 }$$

where t is the time and ν is the coefficient of kinematic viscosity. In the α-viscous disk model (Shakura & Sunyaev 1973), the viscosity ν is given by

$$\begin{equation} \nu = \alpha c_s H \propto r^{(3/2-\beta)}, \end{equation} \tag{ 13 }$$

where α is the standard viscous parameter and Equation (5) is used. The similarity solution of Equation (12) is given by Lynden-Bell & Pringle (1974; see also Hartmann et al. 1998)

$$\begin{eqnarray} \Sigma &=& \frac{C}{3\pi \nu _1} {\bf T^{-5/(1+2\beta)}}r^{(\beta -\frac{3}{2})}\exp {[ -r^{(\beta +1/2)} ]}, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} R_0 &\equiv & T^{2/(1+2\beta)} R_1, \end{eqnarray} \tag{ 15 }$$

where R₁ is an initial radial scale factor, C is a scaling constant, and ν₁ is the coefficient of viscosity at R₁. We note that the disk radius is normalized by R₀, which is equal to R₁ at t = 0 and increases with time. While Hartmann et al. (1998) used the normalization constant R₁, we use R₀ instead since it is more convenient for comparison with the marginally stable state discussed in Section 2.2. The non-dimensional time T is

$$\begin{equation} T \equiv \frac{t}{T_\nu } +1, \end{equation} \tag{ 16 }$$

where the viscous scaling time, T_ν, is

$$\begin{equation} T_\nu \equiv \frac{4}{3(1+2\beta)^{\bf 2}}\frac{R_1^2}{\nu _1}. \end{equation} \tag{ 17 }$$

In Equation (15), we set R₀ increasing with time. We set the surface density normalization constant as

$$\begin{equation} \Sigma _0= \frac{C}{3\pi \nu _1} {\bf T^{-5/(1+2\beta)}}. \end{equation} \tag{ 18 }$$

Then the non-dimensional surface density (σ_s) becomes

$$\begin{eqnarray} \sigma _s &=&r^{(\beta -\frac{3}{2})}\exp {[ -r^{(\beta +1/2)} ]}, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} \frac{\partial (\mathrm{ln }\; \sigma _s)}{\partial (\mathrm{ln }\; r)}&=& -\left(\beta +\frac{1}{2}\right)r^{(\beta +1/2)}+ \left(\beta -\frac{3}{2}\right). \end{eqnarray} \tag{ 20 }$$

The radius at which the disk becomes rotationally unstable, r_m, is given by substituting Equation (19) into Equation (8)

$$\begin{equation} \left(\frac{R_0}{H_0} \right)^{2}=\frac{5}{2}\left(\beta +\frac{1}{2}\right)r_m^{3/2}+\frac{1}{2}(6-\beta)(2-\beta)r_m^{(1-\beta)}. \end{equation} \tag{ 21 }$$

In the outermost part of the disk, r > r_m, the similarity profile cannot be maintained because of rotational instability. Figure 1 shows the relation between r_m and H₀/R₀ given by Equation (21) with β = 0 (red solid line), β = 1/2 (green dashed line), and β = 3/4 (blue dotted line). As discussed in Section 4.1 below, H₀/R₀ of typical protoplanetary disks ranges from 0.1 to 0.3. The critical radius r_m is about 10, 3, and 1 for H₀/R₀ = 0.1, 0.2, and 0.3, respectively (see Figure 1).

Zoom In Zoom Out Reset image size

In this paper, we simply assume that the similarity profile connects smoothly to the marginally stable profile, i.e., σ_ms = σ_s and ∂σ_ms/∂r = ∂σ_s/∂r at r = r_m (e.g., Tanigawa & Ikoma 2007). Then, from Equations (9), (10), and (19)–(21), the constants of integration, C₁ and C₂, are obtained as

$$\begin{eqnarray} C_1 &=& \frac{3}{2}\left(\beta +\frac{1}{2}\right)r_m^{5/2} +\frac{1}{2}(1-\beta)(6-\beta)r_m^{(2-\beta)}, \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} C_2 &=&-\frac{C_1}{\beta -2} r_m^{(\beta -2)}+\frac{\beta -6}{2} \mathrm{ln }\; r_m\nonumber\\ &&+\frac{3(2\beta +3)}{4(\beta -1)}r_m^{(\beta +\frac{1}{2})} +\frac{(2-\beta)(6- \beta)}{2(\beta -1)}. \end{eqnarray} \tag{ 23 }$$

This assumption, however, must be examined and we will discuss the more realistic evolution of the surface density profile in Section 4.3.

3.1.1. Disks Irradiated by the Central Stars (β = 1/2)

The similarity profile of the surface density of the disks is given by Equation (19)

$$\begin{equation} \sigma _s =r^{-1}\exp {\left(-r \right)}. \end{equation} \tag{ 24 }$$

The critical radius where the profile becomes rotationally unstable, r_m, is given by Equation (21)

$$\begin{equation} \left(\frac{R_0}{H_0} \right)^{2}=\frac{5}{2}r_m^{3/2}+\frac{33}{8}r_m^{1/2}. \end{equation} \tag{ 25 }$$

The marginally stable surface density profile for β = 1/2 in Equation (9) is

$$\begin{equation} \sigma _{{\rm ms}} = \exp {\left[-\frac{2}{3}C_1r^{-3/2} +2\left(\frac{R_0}{H_0} \right)^{2}r^{-1/2}+\frac{7}{4}\mathrm{ln }\; r +C_2 \right]}. \end{equation} \tag{ 26 }$$

The smooth connection of the similarity profile to the marginally stable profile at r = r_m, gives the constants of integration (Equations (22) and (23)),

$$\begin{eqnarray} C_1 &=& \frac{3}{2}r_m^{5/2} +\frac{11}{8}r_m^{3/2}, \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} C_2 &=&\frac{2}{3}C_1 r_m^{-3/2}-\frac{11}{4} \mathrm{ln }\; r_m -6r_m -\frac{33}{4}. \end{eqnarray} \tag{ 28 }$$

Figure 2 shows the surface density profile σ as a function of radial distance from the central star r in disks with β = 1/2. The black solid line shows the similarity profile. Other lines show the marginally stable profiles, which smoothly connect to the similarity profile at r_m, for H₀/R₀ = 0.1 (red solid line), H₀/R₀ = 0.2 (green long dashed line), and H₀/R₀ = 0.3 (blue dotted line). Each square point represents r_m. It is apparent that for r > r_m the marginally stable profile deviates from the similarity profile. It is also possible to construct marginally stable profiles for r < r_m, as shown by the dashed lines, but the similarity profile (black solid line) should be realized in r < r_m in the actual disks. For H₀/R₀ = 0.1, the marginally stable profile appears only at the extremely low density region (see the red line). This is because the critical radius r_m, the location of which is shown by the red square, is much further than the disk radius, r = 1, of the similarity solution, i.e., r_m ≫ 1 (see also Figure 1). In such disks, deviation from the similarity profile would hardly be detectable. However, for thicker disks (for larger H₀/R₀), the critical radius r_m decreases, and the marginally stable profile appears at inner part of the disk where the gas density is higher (see green long dashed and blue dotted lines). It is apparent that for H₀/R₀ = 0.3 the marginally stable profile (the blue dotted line) significantly deviates from the similarity profile, keeping relatively high density. In such geometrically thick disks, it is expected that the deviation from the similarity profile will be detected by future observations of high sensitivity.

Zoom In Zoom Out Reset image size

3.1.2. Isothermal Disks (β = 0)

The similarity profile of the surface density of the disks is given by Equation (19):

$$\begin{equation} \sigma _s =r^{-3/2}\exp {(-\sqrt{r})}. \end{equation} \tag{ 29 }$$

The critical radius is given by Equation (21):

$$\begin{equation} \left(\frac{R_0}{H_0} \right)^{2}=\frac{5}{4}r_m^{3/2}+6r_m. \end{equation} \tag{ 30 }$$

The marginally stable surface density profile of the disks is

$$\begin{equation} \sigma _{{\rm ms}} = \exp {\left[-\frac{C_1}{2}r^{-2} +\left(\frac{R_0}{H_0} \right)^{2}r^{-1}+\frac{3}{2}\mathrm{ln }\; r +C_2 \right]}. \end{equation} \tag{ 31 }$$

The constants of integration are

$$\begin{eqnarray} C_1 &=& \frac{3}{4}r_m^{5/2} +3r_m^{2}, \end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray} C_2 &=&\frac{C_1}{2} r_m^{-2}-3 \mathrm{ln }\; r_m -\frac{9}{4}r_m -6. \end{eqnarray} \tag{ 33 }$$

Figure 3 shows the surface density profile σ as a function of radial distance from the central star r in isothermal disks (β = 0). As with Figure 2, the black solid line shows the similarity profile and the other lines show the marginally stable profiles, which smoothly connect to the similarity profile at r_m, for H₀/R₀ = 0.1 (red solid line), H₀/R₀ = 0.2 (green long dashed line), and H₀/R₀ = 0.3 (blue dotted line). Each square point represents r_m. Compared to the β = 1/2 case (Figure 2), in isothermal disks the marginally stable profile outside r_m more significantly deviates from the similarity solution. Even if H₀/R₀ = 0.1 (red solid line), the marginally stable region appears in high surface density region. This is because the similarity profile for β = 0 (black solid line in Figure 3) drops more slowly in the outer region than for β = 1/2 (black solid line in Figure 2). In disks with shallower temperature profiles, the deviation from the similarity profile is more easily detectable. We note that geometrically thin approximation is broken in the region where the surface density increases (see Section 4.2).

Zoom In Zoom Out Reset image size

In the marginally stable region (r ⩾ r_m), the specific angular momentum is uniform (j(r) = j(r_m)) because ∂j/∂R = 0. The rotation velocity (v_ϕ) in r ⩾ r_m is given by

$$\begin{equation} v_{\phi } \equiv \frac{j}{R} = \sqrt{\frac{{GM} r_m}{R_0}} r^{-1}, \end{equation} \tag{ 34 }$$

which decreases more steeply than the Keplerian profile $v_K \equiv \sqrt{{GM}/R_0} r^{-1/2}$, and thus is slower than the Keplerian velocity. In reality, the specific angular momentum j would increases gradually with r (see Section 4.3) but (∂j/∂R) is suppressed. If (∂j/∂R) is significantly reduced, v_ϕ should be much smaller than v_K, and the deviation of rotation velocity from the Keplerian velocity profile can be observable in this region (r > r_m).

The deviation of rotation velocity also causes strong drag forces between dust grains and the gas. Dust grains migrate toward the central star and the drag forces make a difference between the surface density profile of the gas and dust grains. Observations indicate the difference (e.g., Panic et al. 2009) and there are some theoretical works (e.g., Takeuchi et al. 2005; Birnstiel & Andrews 2014) that explain the observations. These theoretical works, however, have not taken into account rotational instability and assume the similarity profiles of the disk gas, that is, they have taken into account the effect of non-Keplerian rotation of the gas for the dust evolution but not for the gas evolution (see Section 4.3 for details). The gas disk profile deviates from the similarity profile in non-Keplerian disks, which affects the dust disk profile. Because future observations of the outer structure of disks, by, e.g., ALMA, are expected to unveil both the gas and dust disk structures, theoretical investigations on the gas and dust evolution in non-Keplerian disks is an important topic.

In this section, we discuss the possibility of the observational detection of non-similarity profiles in protoplanetary disks, which depends on the value of H₀/R₀. In typical T Tauri disks, the scale height is roughly estimated as

$$\begin{equation} \frac{H}{R}\approx 0.1\left (\frac{T}{28\,{\rm K}}\right)^{1/2}\left (\frac{M_*}{M_{\odot }}\right)^{-1/2}\left (\frac{R}{100\,{\rm AU}}\right)^{1/4}, \end{equation} \tag{ 35 }$$

where we assume T = 280 K at R = 1 AU and adopt β = 1/2, considering that the temperature profile is controlled by the central star. In the outer regions where the temperature is as low as that of the surrounding molecular clouds (∼10–30 K, depending on low- or high-mass star forming regions), the temperature profile approaches isothermal (β = 0; see below). Observationally, R₀ ranges from ∼15 AU to 200 AU (Andrews et al. 2009, 2010), and then H₀/R₀ ∼ 0.1–0.18 from Equation (35) if we adopt M_* = 0.5 M_☉. In this case, r_m = R_m/R₀ is roughly 3–10 (see Figure 1). For disks with a relatively less massive central star and relatively high temperature, H₀/R₀ would be close to 0.2. For such disks, the deviation of surface density from the similarity profile is so large near R_m that the deviation will be detected by future observations. If H₀/R₀ is as low as 0.1, it will be difficult to observe the deviation (see Figure 2). We note that Andrews & Williams (2007) reported some deviation from the similarity profiles in their observations of surface density profiles of protoplanetary disks. In typical Herbig Ae disks, the temperature is higher (T ≈ 100 K at R = 10 AU; e.g., Dullemond & Dominik 2004) but the stellar mass is larger (M_* ≈ 2 M_☉). Therefore, H₀/R₀ becomes smaller and the deviation seems more difficult to observe.

For disks irradiated by nearby massive stars in young clusters, the external irradiation makes β ≃ 0 in the outer regions as in the case of protoplanetary disks in the Trapezium Cluster in the Orion Nebula (e.g., Robberto et al. 2002; Walsh et al. 2013). In this case, the scale height is roughly estimated as

$$\begin{equation} \frac{H}{R}\approx 0.2 \left (\frac{T}{60\,{\rm K}}\right)^{1/2}\left (\frac{M_*}{0.5\,M_{\odot }}\right)^{-1/2}\left (\frac{R}{100\,{\rm AU}}\right)^{1/2}, \end{equation} \tag{ 36 }$$

where we assume T = 60 K in the outer region and adopt β = 0. H₀/R₀ approaches about 0.3 if the mass of the central star is as low as ∼0.2 M_☉ (e.g., Hillenbrand et al. 1998). In this case, r_m = R_m/R₀ is roughly 1 (see Figure 1), and the deviation from the similarity profile is possibly observable (see Figure 3). In disks near a massive star, however, the gas in the outer region escapes due to photoevaporation, which should be taken into account (e.g., Johnstone et al. 1998; Richling & Yorke 2000). On the other hand, our result also suggests that the radial pressure gradient force may affect the evolution process due to photoevaporation by surrounding high-mass stars (see Section 4.3).

To summarize, it is expected that deep observations in the near future can detect deviation from the similarity profile of surface density in the outer regions of the disks whose central stars are less massive and temperature is relatively high. In the disks heated by irradiation from a nearby massive star in young star clusters, the deviation will be more easily observable, but we should take into account the effect of photoevaporation on the surface density profile in the outer disks. Observations by facilities with high sensitivity, such as ALMA, will reveal the surface density profiles of the outer regions of protoplanetary disks, which tell us the physical properties of disks with sharp surface density profiles.

In this paper, we assume that disks are geometrically thin. The assumption, however, breaks in the outer region of the disks. From Equations (3) and (5), the scale height is given by

$$\begin{equation} \frac{H}{R}=\frac{H_0}{R_0} r^{(1-\beta)/2}. \end{equation} \tag{ 37 }$$

Figure 4 shows the scale height, H_m, at R_m ≡ R₀r_m as a function of H₀/R₀, obtained by Equations (3), (5), and (21) for β = 0 (red solid line), β = 1/2 (green dashed line), and β = 3/4 (blue dotted line). Because H_m/R_m ∼ 0.4 at most, the figure indicates that geometrically thin approximation is appropriate around R_m. Beyond r_m, H/R exceeds unity around the minimum point of the marginally stable surface density profile in Figures 2 and 3 (for β = 1/2, r ∼ 10⁴, 625, and 100, and for β = 0, r ∼ 100, 25, and 10 for H₀/R₀ = 0.1, 0.2, and 0.3, respectively) and the geometrically thin approximation is invalid beyond the radius. If the gas of the surrounding envelope continues to infall to the disk, the density profile in the outer edge may be more shallower than the self-similar profile. In order to consider this situation, let the density profile be a simple power-law in the whole disk. In this case, R_m would be larger compared to the disk of similarity profiles, resulting in a breakdown of the geometrically thin approximation around R_m (see Equations (8) and (37)). The steep drop of the density at the exponential tail is essential to cause rotational instability.

Zoom In Zoom Out Reset image size

In the previous sections, we simply assume that the surface density profiles connect smoothly from the similarity profile to the marginally stable profile at r_m. In reality, the surface density profile evolves so that it avoids being rotationally unstable. In this subsection, we discuss the time evolution of the surface density profile. Here, we use only three approximations that disks are axisymmetric and geometrically thin and that the self-gravitation of disks is negligible. Now, the equation of continuity is

$$\begin{equation} R\frac{\partial \Sigma }{\partial t}+ \frac{\partial }{\partial R}\left(R\Sigma v_R \right)=0, \end{equation} \tag{ 38 }$$

the equation of motion in the radial direction is

$$\begin{equation} \frac{\partial v_R}{\partial t} +v_R \frac{\partial v_R}{\partial R} -\frac{j^2}{R^3}= -\frac{1}{\rho }\frac{\partial P}{\partial R}-\frac{GM}{R^2}, \end{equation} \tag{ 39 }$$

and the equation of angular momentum transfer is

$$\begin{equation} 2\pi R \Sigma \frac{\partial j}{\partial t} -\dot{M} \frac{\partial j}{\partial R}= \frac{\partial W}{\partial R}, \end{equation} \tag{ 40 }$$

where v_R is the velocity of the gas in the radial direction, $\dot{M}$ is the inward mass flux, and W is the viscous torque caused by turbulent motion. $\dot{M}$ is defined as

$$\begin{equation} \dot{M} \equiv -2\pi R \Sigma v_R. \end{equation} \tag{ 41 }$$

W is defined as

$$\begin{equation} W \equiv 2\pi R^2 T_{R \varphi } =2\pi R^3 \nu \frac{\partial \Omega }{\partial R}, \end{equation} \tag{ 42 }$$

where φ is the azimuthal component in the cylindrical coordinate system and T_Rφ is the Rφ component of the viscous stress tensor. We assume that the other components of the viscous stress tensor are negligible. Using Equations (41) and (42), Equations (38)–(40) are

$$\begin{eqnarray} \frac{\partial \Sigma }{\partial t} &=&\frac{1}{2\pi R} \frac{\partial \dot{M}}{\partial R}, \end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray} \frac{\partial \dot{M}}{\partial t} &=& {-}2\pi R \Sigma \left(\frac{j^2}{R^3} -\frac{GM}{R^2} -\frac{1}{\rho } \frac{\partial P}{\partial R} \right)\nonumber\\ \ms\ms &&{-}\frac{\dot{M}^2}{2\pi R^2 \Sigma }\left(1+ \frac{\partial \mathrm{ln} \; \Sigma }{\partial \mathrm{ln} \; R} -2\frac{\partial \mathrm{ln} \; \dot{M}}{\partial \mathrm{ln} \; R} \right), \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray} \frac{\partial j}{\partial t} &=& \frac{1}{2\pi R\Sigma } \left(\dot{M} \frac{\partial j}{\partial R} +\frac{\partial W}{\partial R} \right). \end{eqnarray} \tag{ 45 }$$

Equations (43)–(45) are differential equations of Σ(R, t), $\dot{M}(R, t)$, and j(R, t).

In the previous sections, we ignore the terms that contain second-order of $\dot{M}$ and the terms of $(\partial \dot{M}/ \partial t)$ and (∂j/∂t). In general, we cannot ignore these terms when ∂j/∂R approaches 0. Here, in order to qualitatively comprehend how the surface density profile evolves when ∂j/∂R approaches 0, we retain only the term which contains second-order of $\dot{M}$, ignoring the terms $(\partial \dot{M}/ \partial t)$ and (∂j/∂t). Then Equations (43)–(45) are rewritten as

$$\begin{eqnarray} \frac{\partial \Sigma }{\partial t} &=&\frac{1}{R}\frac{\partial }{\partial R}\left[ \frac{1}{\left(\frac{\partial j}{\partial R} \right)} \frac{\partial }{\partial R}\left\lbrace \nu j \Sigma \left(2- \frac{\partial \mathrm{ln} \; j}{\partial \mathrm{ln} \; R} \right) \right\rbrace \right],\quad \end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray} j^2&=&{GMR} +\frac{R^3}{\rho }\frac{\partial P}{\partial R}-\frac{1}{\Sigma ^2 \left(\frac{\partial j}{\partial R} \right)^2}\nonumber\\\ms\ms &&\times\left[\frac{\partial }{\partial R} \left\lbrace \nu j \Sigma \left(2- \frac{\partial \mathrm{ln} \; j}{\partial \mathrm{ln} \; R} \right) \right\rbrace \right]^2\nonumber\\\ms\ms &&\times\left(1+ \frac{\partial \mathrm{ln} \; \Sigma }{\partial \mathrm{ln} \; R} -2\frac{\partial \mathrm{ln} \; \dot{M}}{\partial \mathrm{ln} \; R} \right), \end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} \dot{M} &=& \frac{2\pi }{\left(\frac{\partial j}{\partial R} \right)}\frac{\partial }{\partial R} \left\lbrace \nu j \Sigma \left(2- \frac{\partial \mathrm{ln} \; j}{\partial \mathrm{ln} \; R} \right) \right\rbrace. \end{eqnarray} \tag{ 48 }$$

Equation (46) is the evolution of the surface density of a geometrically thin viscous disk (similar to Equation (12), but without an assumption of the Keplerian rotation) and Equation (47) is the equation for radial force balance (similar to Equation (2)). Equations (46) and (48) indicate that when (∂j/∂R) approaches 0 around R ∼ R_m, the surface density diffuses fast and $\dot{M}$ becomes large, and then the third term on right side of Equation (47) becomes non-negligible. However, since the diffusion velocity should be subsonic, the accretion rate is limited to $|\dot{M}|=2\pi R\Sigma c_s$. In this case, (∂j/∂R) is roughly estimated from Equation (48) as

$$\begin{eqnarray} \frac{\partial j}{\partial R}&=&\frac{2\pi }{\dot{M}}\frac{\partial }{\partial R} \left\lbrace \nu j \Sigma \left(2- \frac{\partial \mathrm{ln} \; j}{\partial \mathrm{ln} \; R} \right) \right\rbrace \nonumber \\\ms &\sim & \frac{2\pi }{2\pi R\Sigma c_s}\frac{1}{R}(2\alpha c_s H j\Sigma) \nonumber \\\ms &\sim & \alpha c_s\Omega /\Omega _K, \end{eqnarray} \tag{ 49 }$$

where we adopt ν = αc_sH, H = c_s/Ω_K and j = ΩR². Thus, around R ∼ R_m, (∂j/∂R) is much smaller than that of the disk rotating with the Keplerian velocity, ∂j/∂R ∼ v_K, since c_s < v_K and Ω < Ω_K. Therefore, the surface density profile around R ∼ R_m will be close to the marginally stable profile, σ_ms, obtained in Section 3.1. However, it is not exactly the same as σ_ms since ∂j/∂R ≠ 0. In order to obtain the actual surface density profile, we need to solve Equations (43)–(45) in future work.