# Viscosity prescription for gravitationally unstable accretion disks

###### Abstract

Gravitationally unstable accretion disks emerge in a variety of astrophysical contexts — giant planet formation, FU Orioni outbursts, feeding of AGNs, and the origin of Pop III stars. When a gravitationally unstable disk is unable to cool rapidly it settles into a quasi-stationary, fluctuating gravitoturbulent state, in which its Toomre remains close to a constant value . Here we develop an analytical formalism describing the evolution of such a disk, which is based on the assumptions of and local thermal equilibrium. Our approach works in the presence of additional sources of angular momentum transport (e.g. MRI), as well as external irradiation. Thermal balance dictates a unique value of the gravitoturbulent stress driving disk evolution, which is a function of the local surface density and angular frequency. We compare this approach with other commonly used gravitoturbulent viscosity prescriptions, which specify the explicit dependence of stress on Toomre in an ad hoc fashion, and identify the ones that provide consistent results. We nevertheless argue that our approach is more flexible, robust, and straightforward, and should be given preference in applications. We illustrate this with a couple of analytical calculations — locations of the snow line and of the outer edge of the dead zone in a gravitoturbulent protoplanetary disk — which clearly show the simplicity and versatility of the approach.

###### Subject headings:

accretion, accretion disks — instabilities — (stars:) planetary systems: protoplanetary disks — (galaxies:) quasars: general## 1. Introduction.

Astrophysical accretion disks can be prone to gravitational instability (hereafter GI) when their temperature is low and surface density is high (Safronov 1960; Toomre 1964). In particular, GI is conceivable in the outer parts of protoplanetary disks (Cameron 1978; Boss 1998) at the early stages, when the disks are still massive because of ongoing infall. Accretion outbursts of FU Orioni stars can be driven by the gravito-magnetic cycle in the dense, gravitationally unstable parts of the protoplanetary disk (Audard et al. 2014). Outer parts of AGN disks accreting at high are also expected to become gravitationally unstable far from the black hole (Paczynski 1978a, 1978b; Kozlowski et al. 1979; Kolykhalov & Sunyaev 1980; Goodman 2003). Turk et al. (2009), Clark et al. (2011), Greif et al. (2012) find gravitationally unstable disks around young Population III stars in their simulations of star formation at redshifts of .

Possibility of the GI in a gaseous disk is characterized by the so-called Toomre parameter defined as (Toomre 1964)

(1) |

where and are the surface density and sound speed of the disk, and is the angular frequency (assuming Keplerian rotation around a central object with mass ). In the linear regime GI sets in as from above, where is the threshold value suggested by numerical experiments (Boss 2002; Cossins et al. 2009, 2010).

Operation of the GI is accompanied by enhanced angular momentum transport in the disk driven by the non-axisymmetric gravitational torques. As a result, GI results in mass redistribution in the disk, driving its evolution. Energy dissipation caused by the application of gravitational stresses heats up the disk and tends to oppose the GI.

Because of this feedback, the nonlinear outcome of the GI sensitively depends on another dimensionless parameter — the product of the local cooling time in the disk and . Gammie (2001) demonstrated that when cooling is fast and , the disk disintegrates into a number of bound self-gravitating structures. Such fragmentation has been invoked by Boss (1998) to explain the origin of giant planets (see Rafikov (2005, 2007) for constraints on this scenario). Goodman & Tan (2004) and Levin (2007) suggested that fragmentation of gravitationally unstable quasar disks should result in formation of massive stars migrating through the disk. Clark et al. (2011), Greif et al. (2012), Latif & Schleicher (2014) propose that fragmentation of protostellar disks around Pop III stars can produce low mass extremely metal poor stars that can survive until present days.

In the opposite limit of long cooling time it was shown by Gammie (2001) that the disk settles into a quasi-stationary state of gravitoturbulence. In this regime surface density fluctuates in time, but the disk state averaged over the period longer than the dynamical timescale remains the same. The time-averaged value of the Toomre Q parameter in a gravitoturbulent disk is around , i.e. the disk maintains itself in a marginally stable state. This general picture has been confirmed with global simulations by different groups (Rice et al. 2003, 2005; Durisen et al. 2007; Cossins et al. 2009,2010).

The critical value of the cooling time corresponding to the transition between the gravitoturbulent and fragmenting regimes, depends on a variety of factors. One of the key determinants is the equation of state (EOS) of the gas, with softer EOS promoting fragmentation and resulting in higher (Rice et al. 2005; Jiang & Goodman 2011). Opacity transitions, e.g. due to dust grain evaporation, also affect critical (Johnson & Gammie 2003), as well as other forms of the temperature dependence of opacity (Cossins et al. 2010). External irradiation (Rice et al. 2011) and the details of the disk structure (Meru & Bate 2011a) may also affects the value of .

Some concerns have been raised regarding the convergence of gravitoturbulent disk simulations and the existence of the well-defined , as its value has been claimed (Meru & Bate 2011b; Paardekooper 2012) to vary with the grid resolution. However, later this non-convergence has been traced mainly to the numerical effects (Paardekooper 2011; Meru & Bate 2012; Rice et al. 2014).

It has also been debated whether fragmentation ensues because rapid cooling facilitates collapse of unstable fragments, or because the disk can withstand only certain amount of stress and fragments when the dimensionless stress parameter (Shakura & Sunyaev 1973) exceeds a threshold value . Since in thermal equilibrium in the absence of external energy inputs (Gammie 2001)

(2) |

where is the adiabatic index of gas, the values of and are directly connected. Comparing simulations with different (different EOS) Rice et al. (2005) noted that is essentially independent of and is about 0.06. This led them to suggest that the primary reason for fragmentation is the maximum stress that can be sustained by the disk. However, simulations with external irradiation (Rice et al. 2011) suggest that does depend on the level of irradiation and is thus non-universal.

The main goal of this work is to explore disk characteristics in the gravitoturbulent state, which can persist for a long time over a significant range of radii. For that reason it is important to understand how the disk evolves in this regime under the action of the non-axisymmetric gravitational torques. Currently direct 2D and 3D simulations are too numerically expensive to permit such exploration, and one often has to resort to azimuthally-averaged, one-dimensional (in radius) disk models. To evolve them properly one must provide a description of the angular momentum transport by the gravitoturbulence. Formulating such a description is the focus of the present work.

Balbus & Papaloizou (1999) argued that due to the long-range nature of the gravitational interaction the angular momentum transport due to the non-axisymmetric gravitational perturbations is inherently non-local. However, Gammie (2001) showed that in cold, geometrically thin disks angular momentum transport by the gravitational torques can still be described as a local process. Lodato & Rice (2004, 2005) and Cossins et al. (2009) confirmed this conclusion numerically, certainly for the low-mass disks.

In this work we adopt the latter point of view and will explicitly assume the gravitoturbulent transport to be local and characterized by the effective viscosity parameter . Different explicit and implicit prescriptions for have been proposed, which can be classified into two general categories.

The first class of prescriptions relies on the fact that in local dynamical and thermal equilibrium the angular momentum transport in the disk is intimately related to its thermal state. This allows one to directly express via the disk temperature and density. Using constant disk without external energy inputs as an example, Rafikov (2009) has demonstrated that this property, when coupled with the defining characteristic of the gravitoturbulence — the condition , allows one to directly relate to providing a complete description of the disk evolution. This approach has also been implicitly featured in several numerical studies (Terquem 2008; Clarke 2009; Zhu et al. 2009a).

A different way of describing the gravitoturbulent disk evolution does not explicitly assume . Instead, it specifies the explicit dependence of on (Kratter et al. 2008; Zhu et al. 2009b, 2010a, 2010b; Martin & Lubow 2011, 2013, 2014). This approach dates back to the work of Lin & Pringle (1987), who suggested that based on a phenomenological description of transport in gravitationally unstable disks. However, in most cases such prescriptions do not naturally follow from physical arguments. Instead, they are designed to replicate certain qualitative features of behavior.

The goal of this work is to formulate, analyze and compare different approaches to characterizing . First, in §2 we discuss equations used to describe evolution of gravitationally unstable disks. Then, in §3 we highlight the differences between the two aforementioned approaches to closing the system of evolution equations. After considering the steady, constant models in §4, we contrast the performance of different gravitoturbulent viscosity prescriptions in §5. We discuss our results in §6, where we show, in particular, how the gravitoturbulent prescription based on condition can be naturally used to obtain simple analytical expressions for the locations of the ice lines and dead zone edges in the gravitoturbulent protoplanetary disks.

## 2. Basic equations.

Provided that gravitoturbulent stress can be described as local -viscosity, it is convenient (Rafikov 2013) to characterize the disk structure via the angular momentum flux defined (in a Keplerian disk) as

(3) |

where is the specific angular momentum and is the kinematic viscosity expressed through the dimensionless parameter as (Shakura & Sunyaev 1973)

(4) |

As always, we will assume the gas sound speed to be determined by the midplane temperature of the disk. Angular momentum transport is effected both by gravitoturbulence and by any other background sources of effective viscosity, such as MRI.

Mass accretion rate through the disk (defined to be positive for inflow) is related to via a simple relation (Rafikov 2013)

(5) |

Coupled with the continuity equation, this results in the following evolution equation:

(6) |

Expressing via equation (3) one reproduces the classical equation for the viscous disk evolution (Lightman & Eardley 1974; Lin & Papaloizou 1996)

(7) |

Solving either of these equations requires the knowledge of the thermodynamical properties of the disk since both and are functions of gas temperature , see equations (3) & (4).

We now address the thermal state of the disk. We assume the disk to be in thermal equilibrium, so that energy sources and sinks are in local balance. The former include viscous dissipation and external irradiation intercepted by the disk. The latter is effected by radiative cooling. For local angular momentum transport the rate of viscous energy dissipation per unit radial distance in the disk is . Then the rate of energy production by viscous stresses per unit area of the disk is

(8) |

Gravitationally unstable disks are often studied in the limit of a self-luminous disk (Gammie 2001; Rafikov 2009), in which viscous heating of any nature dominates over the energy input due to external irradiation. This assumption should be reasonable for e.g. Class 0 T Tauri stars, which are characterized by intense mass accretion and relatively low luminosity of a very young central star, which is also heavily obscured.

However, in general, the disk, especially its outer parts, may also be heated appreciably by the external radiation field with temperature (Rice et al. 2011; Zhu et al. 2012); can be a function of if irradiation is due to the central object. To account for this possibility we describe the local thermal balance, which ultimately determines the midplane disk temperature , via the following approximate relation:

(9) |

Details of the vertical transport of radiation in the disk are specified here via an explicit function of the optical depth

(10) |

To zeroth order the disk optical depth is determined predominantly by the midplane value of the temperature ( can also be a function of gas density).

The explicit form of depends on whether the disk is optically thick or thin and on the mechanism of the vertical energy transport (Rafikov 2007). Here we assume radiative energy transfer (Rafikov 2007 has derived for convective transport of energy) with dust-dominated opacity. In this case can be reasonably well approximated (up to constant factors of order unity) by

(11) |

This expression does not discriminate between the Rosseland and Planck mean opacities and is approximately valid (up to constant factors of order unity) both in the optically thick and thin regimes (Goodman 2003). It is easy to see that the condition (9) supplemented with equation (11) properly describes the disk thermal balance in different asymptotic regimes of . Indeed, in the optically thick case () one finds , while in the optically thin regime () one has . Both expressions are valid up to constant factors in the appropriate limits.

For simplicity in this work we will consider the following behavior of appropriate for dust opacity in certain temperature regimes (Bell & Lin 1994):

(12) |

with , being constants. For example, at low temperatures, below K the opacity is dominated by ice grains and is characterized by (Bell & Lin 1994)

(13) |

Our results can be trivially extended to other, more complicated forms of .

Equations (9)-(13) provide the sought relation between the disk temperature and surface density, which is used to determine the viscosity behavior in equation (7).

Our final note on disk thermodynamics concerns the case when cooling is not due to dust emission. In particular, Latif & Schleicher (2014a,b) considered the structure of gravitationally unstable disks around Pop III stars. Such disks are believed to be metal-free (i.e. no dust opacity), with cooling provided by molecular hydrogen. In this case, the right-hand side of the energy balance equation (9) should be modified to , where is the disk scale height and is the volumetric cooling rate. This expression is valid in the optically thin regime, with certain modifications required when line cooling switches to the optically thick regime (Latif & Schleicher 2014a; Ripamonti & Abel 2004).

## 3. Closure in the gravitoturbulent state.

Definition (3) and thermal balance condition (9) supplemented with equations (10) and (11) allow one to uniquely relate to if the behavior of is known. This situation is typical for the disk regions dominated by the background viscosity (e.g. due to MRI), where one usually knows (or, more often, postulates) some behavior of . Knowing the dependence one can express and as a function of only, leaving as the only dependent variable in the equation (7). This provides a closure of the system of equations for the density, angular momentum, and thermal balance, and allows one to self-consistently evolve in time and space.

In gravitoturbulent state the behavior of is not known a priori and this approach does not apply. To close the system of evolution equations one needs to make additional assumptions. Two conceptually different approaches to this problem are described next.

### 3.1. Closure via the condition.

Rafikov (2009) pointed out that a necessary closure is naturally provided by the the basic property of a gravitoturbulent disk following directly from simulations — that the disk hovers on the margin of gravitational stability with

(14) |

where .

Indeed, with definition (1) this condition predicts an explicit relation between and in the form

(15) |

Plugging this expression into the thermal balance condition (9) gives the following expression for the gravitoturbulent stress

(16) |

The dependence on surface density arises here because , which enters this expression both explicitly and also through

(17) |

is a function of .

Combining equations (3)-(16) one also comes up with the explicit expression for the effective -parameter due to gravitoturbulence:

(18) |

The dependence on comes only through .
This result was first obtained by Rafikov (2009) in the
non-irradiated case () when studying constant
gravitoturbulent disks. However, it clearly
also holds more generally, for evolving and irradiated
gravitoturbulent disks, because thermal balance gets
established faster than the viscous evolution proceeds.
Note that Rice et al. (2011) found a different (linear)
scaling with^{1}^{1}1Note that the expression in parentheses in
equation (18) can be written as
, with
(Rice et al. 2011). in their expression for the
gravitoturbulent because they used a different
cooling model, namely assuming a constant cooling time.
Expression (18) for completes
the closure and makes evolution equation (7)
self-consistent.

On the other hand, gravitoturbulent disk evolution can be described entirely without introducing -parameter. Indeed, substituting given by the expression (16) into equation (6) one finds

(19) |

where, again, , and behavior is specified. This (in general nonlinear) equation is a closed-form, fully self-contained evolution equation for as a function of and . It represents one of the main results of this work.

In the parts of the disk where external irradiation plays the dominant role this equation adopts a simple time-independent solution , so that (Rafikov 2009).

#### 3.1.1 Additional sources of viscosity.

Angular momentum transport in the disk may be effected not only by gravitoturbulence but also by additional stresses. This would be the case, for example, if the disk is both sufficiently ionized for the MRI to operate and massive enough for being gravitationally unstable.

To account for the possibility of additional, non-gravitoturbulent viscosity parametrized by -parameter , one can simply write

(20) |

Here is still given by the expression (16) with the switch function introduced as described below. The non-gravitational viscous angular momentum flux is

(21) |

and the behavior of is specified. Thermal balance relation (9) now gives

(22) |

from which we determine as a function of and . Then equation (20) yields as a function of and , allowing the evolution of to be followed using equation (6).

A subtle point in this prescription is that equation (20) always includes the gravitoturbulent transport in the form (16), with independent of the actual disk temperature and the value of . As a result, our prescription formally has non-zero gravitoturbulent stress even in the gravitationally stable part of the disk, which is not expected.

For this reason we introduce a switch function in equation (20), which takes care of this issue. There are several different ways in which this can be done. One method would be to adopt which quickly goes to zero as soon as . However, in the absence of a microscopic theory of gravitoturbulence such a factor would necessarily be introduced in an ad hoc fashion.

Instead, we have chosen to use

(23) |

where is the Heaviside step function ( for ; for ). With this approach we simply keep in the form (16) in equation (20), as long as it is positive. This is not a problem, since, as demonstrated in Rafikov (2009), as soon as exceeds the threshold value , the transport is guaranteed to be dominated by the viscous stress, i.e. . In other words, in gravitationally stable disks, even if . And vice versa, and when the disk is gravitationally unstable and . This situation is further discussed in §4 and is illustrated in Figure 5, where we display the actual run of the viscous and gravitoturbulent transport contributions for a particular disk model.

Under certain circumstances one may find that in the gravitationally stable part of the disk, so that according to equation (16). In this case our choice (23) of simply guarantees that the gravitoturbulent part of the angular momentum flux vanishes and exactly.

Evolution equation (6) based on the prescription (20) with given by equation (23) allows us to explore the transition between the disk regions dominated by gravitoturbulence and background viscosity. In the appropriate limits its solution reduces to either the gravitoturbulent disk solution (for large ) following from equation (19), or the conventional viscous disk solution (for small , where ). Alternatively, one can evolve the disk structure using equation (7) with determined by

(24) |

It is easy to show that this approach results in almost the same disk properties, except for the region where and the disk is close to marginal gravitational stability.

### 3.2. Closure via the explicit dependence.

An alternative closure scheme that has been broadly used in the literature postulates some dependence of the gravitoturbulent viscosity on . In this case closure is analogous to a regular viscous disk anzatz: equations (1), (3) with , and (9)-(11) are combined to yield a unique relation. It is then plugged into the relation (4) for viscosity, resulting in the expression for and allowing equation (7) to be evolved in time.

Note that the condition (14) is not used explicitly in this approach and thus in general there is no guarantee that this prescription would result in a truly gravitoturbulent disk structure, with the disk hovering at the edge of instability with .

In the absence of a microscopic model of gravitoturbulence that could motivate a possible dependence of on , several ad hoc prescriptions for have been suggested. Their main feature is the rapid increase of as approaches some critical value from above.

In particular, Zhu et al. (2010a,b)
used^{2}^{2}2Zhu et al. (2009b) used
a very similar prescription .

(25) |

Martin & Lubow (2011,2013,2014) adopted

(26) |

with and . The dependence is motivated by the work of Lin & Pringle (1987). Note that this prescription explicitly relates the value of gravitoturbulent viscosity to the background viscosity . This recipe was also used in Owen & Jacquet (2014) to explore the chemical evolution of a protoplanetary disk undergoing accretion outbursts.

Kratter et al. (2008) have used

(27) |

One can see that reduces to if one uses and instead of and . Also, saturates at a constant value for , unlike , which can be arbitrarily large for small .

## 4. Constant disks.

One is often interested in knowing the structure of a steady-state accretion disk, which necessarily has independent of . This limiting case provides a nice point of comparison between the different types of viscosity prescriptions and will be used in §5.

If is independent of (and, correspondingly, ) and there are no external torques acting on the disk, then equation (5) naturally yields (Rafikov 2013). We now illustrate the difference in approaches between the various closure philosophies when calculating the constant disk structure.

### 4.1. Constant disk: closure

When we explicitly demand a constant disk to be marginally unstable with the condition (14), the equation (20) implies a relation between , , in the form

(28) |

The second algebraic relation between these variables is provided by equations (9)-(12), thus fully specifying the behavior of and . Rafikov (2009) and Clarke (2009) used this approach to understand the properties of gravitoturbulent disks with constant .

### 4.2. Constant disk: explicit closure

When one uses an explicit prescription, the approach to determining disk structure is different from that in §4.1. Assuming a non-zero background viscosity, equation (3) now yields for the constant disk

(29) |

Here both and explicitly depend on temperature. Thus, is a function of in this approach, unlike the case considered in §3.1.

## 5. Comparison of different prescriptions

We now compare the performance of different prescriptions for described in §3.1-3.2. To that effect we construct steady state (constant ) models of gravitoturbulent protoplanetary disks using the results of §4. We then compare the behaviors of various disk characteristics obtained with different recipes.

Our calculations adopt the opacity behavior (12)-(13) typical for low temperatures K (Bell & Lin 1994). For simplicity we will assume this behavior to extend also to higher temperatures; this should not be a problem as out present goal is to explore the differences between the various prescriptions. We take the mass of the central star to be and set the threshold for gravitational stability at . In this and subsequent calculations we assume that fragmentation ensues when exceeds (rather than unity, as e.g. is assumed in equations (A1)-(A3)). Numerical studies suggest that a lower value of is more realistic (Rice et al. 2005; Paardekooper et al. 2011).

Rafikov (2009) showed that gravitoturbulent state results in a set of fiducial values of various disk variables: surface density , midplane temperature and sound speed and , mass accretion rate , and angular frequency . They are chosen such that at the radius corresponding to the angular frequency a steady gravitoturbulent disk with the mass accretion rate equal to simultaneously (1) has midplane temperature and sound speed equal to and , (2) fragments, and (3) has optical depth . In Appendix A we provide explicit expressions and numerical estimates for these variables for the case of interest to us (opacity ). These fiducial quantities provide characteristic values of the important parameters of gravitoturbulent disk and will be used in our subsequent comparisons.

For a given mass of a central object characteristic angular frequency determines a fiducial distance according to the formula

(30) |

where we took the numerical value of from equation (A2) in Appendix A. This is the radius beyond which an optically thick disk with opacity inevitably fragments (Matzner & Levin 2005; Rafikov 2009; Clarke 2009).

In all our models we assume that the disk is immersed in a uniform radiation field with the radially-independent temperature K. We have checked that all our conclusions remain the same for other values of , in particular for the self-luminous disks with .

In Figure 2 we show the radial profiles of the midplane temperature , surface density , optical depth , -parameter and Toomre for a gravitoturbulent disk that has a high mass accretion rate, M yr (leaving aside the question of how realistic such is). We also use a relatively high value of the background viscosity, , which is often adopted in the literature (Zhu etal 2009a; Martin & Lubow 2011, 2013). Different curves correspond to prescriptions given by equations (18), (25)-(27), as indicated in the Figure.

Shaded region outside of AU corresponds to a fragmenting part of the disk where , and the behavior of disk parameters in this region is irrelevant (many curves are very discrepant there). As expected (Rafikov 2009), at the fragmentation edge our disk model has , , and is optically thick ().

One can see that all prescriptions reliably reproduce the radius at which the transition from the gravitationally stable to the gravitoturbulent state occurs (around the vertical line). Interior to this radius background viscosity dominates, , and different curves overlap. But outside this radius, in the region where , some differences between prescriptions start emerging. In particular, at the fragmentation edge (boundary of the shaded region) computed through closure is different from of Martin & Lubow (2011) by about a factor of . Similar differences at the same AU are seen in the values of and . Other explicit prescriptions show better agreement (within tens of per cent) for , , and with the results of closure. On the other hand, the behavior of the midplane temperature (Figure 2e) shows good agreement between different prescriptions. As expected, never drops below about 10 K as it is limited by far from the star; however, for this disk model this is true only beyond the fragmentation edge so that in practice irradiation with the adopted is irrelevant here.

To investigate these differences further we recomputed the disk structure for the same high mass accretion rate M yr but decreasing the background viscosity by two orders of magnitude, to . This lower value of may be more realistic for the cold protoplanetary disks, where the MRI is weakened by the non-ideal effects such as ambipolar diffusion (Bai & Stone 2011), or may not be operating at all, e.g. in a dead zone (Gammie 1996; Fleming & Stone 2003).

Figure 3 shows the results of this calculation. Its inspection reveals good quantitative agreement between the models computed via closure (black solid), which rely on given by the equation (18), and those using explicit viscosity in the form proposed by Kratter et al. (2008), , and Zhu et al. (2010a,b), . They typically agree with several tens of per cent for all variables in the whole gravitoturbulent region of the disk, between the fragmentation edge and the line.

It is also obvious that the use of suggested by Martin & Lubow (2011) and this low leads to rather poor quantitative agreement with all other prescriptions: , and show a discrepancy of more than an order of magnitude in the outer disk, at the edge of the fragmentation zone. The value of at this radius plunges to for , in contrast to maintained by all other prescriptions in agreement with simulations. Lower value of used in this calculation pushes gravitoturbulent zone down to AU, accentuating the deviations which were already visible at higher , with less extended gravitoturbulent region (see Figure 2).

This general situation holds for other disk models. For example, in Figure 4 we compare different recipes in a low- disk with M yr, while keeping . This disk is optically thin in its outer parts, with transition happening interior to . In this case, as shown in Rafikov (2009), fragmentation can be pushed out to very large radii and the disk can extend beyond AU in the optically thin, gravitoturbulent regime. Moreover, external irradiation is very important in this case: K sets the temperature floor outside 100 AU, in the part of the disk which is stable against fragmentation, see Figure 4e.

Again, with such low disk models using suggested by Martin & Lubow (2011) are considerably different from all other models, which tend to agree with each other. For the discrepancy can be as high as an order of magnitude. From this we can conclude that prescription underperforms at low values . We discuss the reasons for this next.

## 6. Discussion.

Our interpretation of the behavior found in models using is that it is due to the rather gradual dependence of the former on . Indeed, Figure 1 shows that when is between and the critical value for fragmentation ( in our case) grows with decreasing much slower than either or . This is especially obvious in the case, when has to go down to for to reach . This explains a significant drop in towards the fragmentation edge in Figures 3b and 4b. Thus, we conclude that does not properly capture the main feature of the gravitoturbulent disk — almost constant and close to unity value of the Toomre . The explicit dependence of on may be another feature that causes it to underpeform compared to and .

At the same time, we must emphasize that we are not aware of the existing gravitoturbulent disk calculations using low — most of them assume . Thus, Figures 3 and 4 are not providing comparison with any published results, but merely urging caution in choosing recipe when working with low . And at higher calculations using can still be acceptable (see Figure 2), depending on the required accuracy.

Explicit prescriptions suggested in Kratter et al. (2008) and Zhu et al. (2010a,b) ( and ) do rather well at reproducing the calculation using . Based on this we expect that effectively any explicit scheme that exhibits very sharp rise in the vicinity of as decreases, should be suitable for practical calculations of the gravitoturbulent disk properties. In fact, the sharper is the better and something as simple as with will maintain in the gravitoturbulent state quite well (conditional upon ).

Nevertheless, we strongly believe that there are significant benefits to using the approach outlined in §3.1. First, the only assumption that it employs is that the disk is able to maintain itself in a state of marginal gravitational stability . This stipulation is in agreement with the results of direct numerical simulations (Cossins et al. 2009, 2010). Explicit prescriptions also make this assumption, even though implicitly. But in addition they have to postulate some actual functional form of the dependence, which does not have physical justification and is always introduced in an ad hoc fashion. This makes the explicit approach rather arbitrary, which may result in certain problems, as we have demonstrated in §5.

Second, deep in the gravitoturbulent regime, when one can neglect the background viscosity compared to the gravitoturbulent contribution , the explicit approach results in a self-contained equation (19) for the surface density evolution. This is allowed by the fact that the expression (18) already implicitly accounts for the thermal balance of the disk, even in presence of external irradiation. Disk temperature is uniquely defined in this case by the relation (15) after the disk surface density has been solved for. This is not the case for the explicit recipes, which always need to explicitly relate to via the equation of thermal balance (even deep in the gravitoturbulent state, when ), because both enter the definition (1) of the Toomre . Thus, our prescription allows a more efficient numerical exploration of the gravitoturbulent disk structure.

Third, when the background viscosity cannot be neglected compared to , the explicit prescription (18) does not make the task of computing the disk properties more complicated than the explicit prescriptions, see §3.1.1. Dealing with this case using our prescription (24), which results in a non-zero even in the gravitationally stable part of the disk, is not a problem.

Indeed, in Figure 5 we demonstrate
that computed in our disk models via equation
(18) rapidly becomes subdominant as
the disk transitions into the regime dominated by the background
viscosity . Thus, even though in general
is not precisely zero^{3}^{3}3Graviturbulent contribution does
vanish inside of 25 AU in Figure 5 for
because of our use of in the form
(23) and the fact that drops below
inside this radius. in this part of the disk,
it still does not affect the disk properties. The only
place where our approach may be
quantitatively less accurate is around the transition
. However, this is the location
where all prescriptions have
some degree of arbitrariness.

Fourth, prescription allows a much more transparent analytical derivation of certain gravitoturbulent disk characteristics. We demonstrate this next when we compare the derivations of the snow line location (§6.1) and of the dead zone edge (§6.2) in gravitoturbulent disks using the two classes of prescriptions.

### 6.1. Snow line location

Snow (or ice) lines for different volatile materials (HO, CO, etc.) arise in protoplanetary disks at the locations where these materials undergo a transition from solid to gaseous phase. It is conventional to identify snow lines with the disk radii where the midplane temperature is equal to some characteristic value for sublimation of a particular material at the typical midplane pressure. Following common wisdom we will adopt this as a working definition of a snow line.

Martin & Livio (2013) derived an analytical expression
for the snow line radius in a gravitoturbulent,
optically thick, self-luminous protoplanetary disk with
constant . Here we compare their calculation with
the derivation of using our favored
recipe (18) based on
the condition . To facilitate comparison we
employ a common set of assumptions: (1) neglect background
viscosity and external irradiation
(i.e. ), (2) describe optically thick radiation
transfer by^{4}^{4}4This is what equations (4) & (5) of
Martin & Livio (2013) effectively reduce to.

(31) |

in place of our equation (11), with (instead of our definition (17)), and (3) take opacity at the icy grain sublimation radius to be given by equation (12) with and (in appropriate CGS units).

#### 6.1.1 via the prescription

In our approach we first express via and by setting in equation (15). Then we plug this into the expression (18) for . Finally, inserting all this into equation (28), written in the form with all variables expressed through and , and resolving it for we get

(32) | |||||

Note that this derivation does not use thermodynamic balance condition (9) as it has already been used in deriving the expression (18) for .

#### 6.1.2 via the explicit

Now we outline the steps used in Martin & Lubow (2013)
to derive . They adopted an explicit
with^{5}^{5}5Martin & Livio (2013) used rather
than to denote in their
equation (9) for ,
but this notation would cause confusion with the background viscosity
in our paper. . This expression is similar to
that in Zhu et al. (2009b; see our equation (25))
in spirit, but not in detail, as the maximum value of
is modulated by a small factor
in the Martin & Livio (2013) case. Note that this anzatz
would never allow a gravitationally unstable disk to fragment
unless the critical value of for fragmentation
were very low, below .

With this prescription the constant assumption, embodied in the equation (29) allowed Martin & Livio (2013) to relate at the snow line to the disk midplane temperature (assumed equal to ) and . Because of the adopted form of this relation depended on and involved a non-elementary function (Lambert function), complicating the analysis. The latter problem was overcome by using an asymptotic relation for the Lambert function, valid for relatively low values of M yr.

Additional (algebraic) relation between , , and at the snow line comes from the thermodynamical relation (9) with . Eliminating with the aid of these two relations Martin & Lubow (2013) obtained a scaling for in terms of powers of , , and the central star mass , given by the equation (19) of their paper.

#### 6.1.3 Comparison of approaches

Comparing the result for in Martin & Lubow (2013) with equation (32) shows that the two are essentially identical, with all the scalings being the same and the numerical estimates different at level for . However, it is clear that our derivation of the expression (32) is much more straightforward and flexible.

First, it involves only elementary functions. Second, it does not involve ad hoc factors, such as used by Martin & Lubow (2013) in their definition of , which are confusing and do not appear in the final result anyway. Third, our result (32) is clearly valid regardless of the value of , while the asymptotic representation of the Lambert function used by Martin & Lubow (2013) works only for relatively low M yr.

Analogous problems would arise if one were to use other explicit prescriptions for analytical calculations. This provides strong motivation for adopting the approach advocated in this work rather than the explicit recipes, when studying gravitoturbulent disks.

### 6.2. Outer edge of the dead zone

In a very similar vein we can estimate the location of the outer edge of a dead zone (Gammie 1996) in a gravitoturbulent protoplanetary disk. This edge should be located at the radius , where the disk surface density is twice the surface density of the active layer g cm, down to which external ionizing radiation can penetrate into the disk and keep it fully ionized.

For this estimate the equation (15) of the approach allows one to express the midplane temperature through and , while is already a function of the two, see equation (18). Then, repeating the steps in §6.1.1 and eliminating one obtains

(33) | |||||

In making the estimate we again adopted the opacity in the form (13), which is perfectly reasonable since and ice grains dominate opacity.

It is easy to see that deriving via some explicit prescription would again require going through the rather contrived procedure described in §6.1.2, involving non-elementary functions, their asymptotic expansions, and so on. Thus, similar to §6.1.3, we can again argue that the approach provides the shortest and clearest path to finding .

## 7. Conclusions.

One of the main outcomes of this work is the development of a self-consistent analytical formalism for the evolution of a gravitoturbulent accretion disk, including the possibility of external irradiation. It explicitly enforces Toomre in the disk to stay close to the value needed for marginal gravitational stability. Implicit numerical implementations of this approach do exist (e.g. Terquem 2008; Clarke 2009; Zhu et al. 2009a) but a complete analytical formalism was lacking until now. Our work fills this gap with the prescriptions described in §3.1, in particular the explicit expression (18) for the gravitoturbulent viscosity , which depends only on and and fully accounts for the thermodynamic equilibrium of the disk. Our formalism includes the possibility of the non-zero background viscosity (§3.1.1) provided by sources other than the gravitoturbulence (e.g. MRI), and external irradiation of the disk. This work thus generalizes the models of steady, constant gravitoturbulent disks explored in Rafikov (2009).

In fully gravitoturbulent state, when we derive an explicit equation (19) for the surface density evolution, which is fully self-contained — it involves only and — and is non-linear in general. This equation provides a powerful tool for exploring gravitoturbulent disk evolution and is thus an important result of our work.

We then contrasted our approach with the calculations specifying an explicit dependence of the gravitoturbulent viscosity on the Toomre , which is often done in the literature. We clearly demonstrated that our approach provides a much faster and more natural route to deriving analytical expressions for various important characteristics of the gravitoturbulent disks, such as the locations of the snow lines (§6.1) and the edge of the dead zone (6.2). We also compared disk models computed using different recipes (§5) and found that some explicit prescriptions are able to reproduce the results obtained using approach reasonably well. This is the case only for those prescriptions that cause to increase very steeply as decreases in the vicinity of (e.g. Kratter et al. 2008; Zhu et al. 2010a,b). Certain disagreement may arise if this requirement is violated (Martin & Lubow 2010), especially when the background (non-gravitoturbulent) viscosity in the disk is low (see our discussion in §6). Thus, although it is often stated in the literature that the precise form of the dependence used for building gravitoturbulent disk models is unimportant, we find this to be only partly true.

Finally, we emphasize that our formalism has the minimal number of imposed assumptions — namely, that is approximately maintained in the gravitoturbulent state, in agreement with simulations. It thus provides a natural and robust pathway to building efficient time-dependent models of gravitoturbulent disks around objects such as young stars, quasars, and Pop III stars.

## References

- (1)
- (2) Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, arXiv:1401.3368
- (3) Bai, X.-N. & Stone, J. M. 2011, 736, id. 144
- (4) Balbus, S. A. & Papaloizou, J. C. B. 1999, ApJ, 521, 650
- (5) Bell, K. R. & Lin, D. N. C. 1994, ApJ, 427, 987
- (6) Boss, A. P. 1998, ApJ, 503, 923
- (7) Boss, A. P. 2002, ApJ, 576, 462
- (8) Cameron, A. G. W. 1978, Moon Planets, 18, 5
- (9) Clark, P. C., Glover, S. C. O., Smith, R. J., Greif, T. H., Klessen, R. S., & Bromm, V. 2011, Science, 331, 1040
- (10) Clarke, C. J. 2009, MNRAS, 396, 1066
- (11) Cossins, P., Lodato, G., & Clarke, C. J. 2009, MNRAS, 393, 1157
- (12) Cossins, P., Lodato, G., & Clarke, C. J. 2010, MNRAS, 401, 2587
- (13) Durisen, R. H., Boss, A. P., Mayer, L., Nelson, A. F., Quinn, T., & Rice, W. K. M. 2007, in Protostars and Planets V, B. Reipurth, D. Jewitt, and K. Keil (eds.), University of Arizona Press, Tucson; p. 607
- (14) Fleming, T. & Stone, J. M. 2003, ApJ, 585, 908
- (15) Gammie, C. F. 1996, ApJ, 457, 355
- (16) Gammie, C. F. 2001, ApJ, 553, 174
- (17) Greif, T. H., Bromm, V., Clark, P. C., Glover, S. C. O., Smith, R. J., Klessen, R. S., Yoshida, N., & Springel, V. 2012, MNRAS, 424, 399
- (18) Goodman, J. 2003, MNRAS, 339, 937
- (19) Goodman, J. & Tan, J. C. 2004, ApJ, 608, 108
- (20) Jiang, Y.-F. & Goodman, J. 2011, ApJ, 730,id. 45
- (21) Johnson, B. M. & Gammie, C. F. 2003, ApJ, 597, 131
- (22) Kolykhalov, P. I. & Syunyaev, R. A. 1980, Soviet Ast. Lett., 6, 357
- (23) Kozlowski, M., Wiita, P. J., & Paczynski, B. 1979, Acta Astronomica, 29, 157
- (24) Kratter, K. M., Matzner, C. D., & Krumholz, M. R. 2008, ApJ, 681, 375
- (25) Latif, M. A. & Schleicher, D. R. G. 2014a, arXiv:1411.0096
- (26) Latif, M. A. & Schleicher, D. R. G. 2014b, arXiv:1411.5902
- (27) Levin, Y. 2007, MNRAS, 374, 515
- (28) Lightman, A. P., & Eardley, D. M. 1974, ApJL, 187, L1
- (29) Lin, D. N. C., & Papaloizou, J. 1996, ARA&A, 34, 703
- (30) Lin, D. N. C. & Pringle, J. E. 1987, MNRAS, 225, 607
- (31) Lodato, G. & Rice, W. K. M. 2004, MNRAS, 351, 630
- (32) Lodato, G. & Rice, W. K. M. 2005, MNRAS, 358, 1489
- (33) Martin, R. G. & Livio, M. 2013, MNRAS, 434, 633
- (34) Martin, R. G. & Lubow, S. H. 2011, ApJL, 740, id. L6
- (35) Martin, R. G. & Lubow, S. H. 2013, MNRAS, 432, 1616
- (36) Martin, R. G. & Lubow, S. H. 2014, MNRAS, 437, 682
- (37) Matzner, C. D. & Levin, Y. 2005, ApJ, 628, 817
- (38) Meru, F. & Bate, M. R. 2011a, MNRAS, 410, 559
- (39) Meru, F. & Bate, M. R. 2011b, MNRAS, 411, L1
- (40) Owen, J. E. & Jacquet, E. 2014, arXiv:1410.6844
- (41) Paardekooper, S.-J. 2012, MNRAS, 421, 3286
- (42) Paardekooper, S.-J., Baruteau, C., & Meru, F. 2011, MNRAS, 416, L65
- (43) Paczynski, B. 1987a, Acta Astronomica, 28, 91
- (44) Paczynski, B. 1987b, Acta Astronomica, 28, 241
- (45) Rafikov, R. R. 2005, ApJ, 621, L69
- (46) Rafikov, R. R. 2007, ApJ, 662, 642
- (47) Rafikov, R. R. 2009, ApJ, 704, 281
- (48) Rafikov, R. R. 2013, ApJ, 774, id. 144
- (49) Rice, W. K. M., Armitage, P. J., Bate, M. R., & Bonnell, I. A. 2003, MNRAS, 339, 1025
- (50) Rice, W. K. M., Armitage, P. J., Mamatsashvili, G. R., Lodato, G., & Clarke, C. J. 2011, MNRAS, 418, 1356
- (51) Rice, W. K. M., Paardekooper, S.-J., Forgan, D. H., & Armitage, P. J. 2014, MNRAS, 438, 1593
- (52) Rice, W. K. M., Lodato, G., & Armitage, P. J. 2005, MNRAS, 364, L56
- (53) Ripamonti, E. & Abel, T. 2004, MNRAS, 348, 1019
- (54) Safronov, V. S. 1960, Annales dâAstrophys., 23, 979
- (55) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
- (56) Terquem, C. E. J. M. L. J. 2008, ApJ, 689, 532
- (57) Toomre, A. 1964, ApJ, 139, 1217
- (58) Turk, M. J., Abel, T., & O’Shea, B. 2009, Science, 325, 601
- (59) Zhu, Z., Hartmann, L., & Gammie, C. 2009a, ApJ, 694, 1045
- (60) Zhu, Z., Hartmann, L., & Gammie, C. 2010a, ApJ, 713, 1143
- (61) Zhu, Z., Hartmann, L., Gammie, C., & McKinney, J. C. 2009b, ApJ, 701, 620
- (62) Zhu, Z., Hartmann, L., Gammie, C., McKinney, J. C., Book, L. G., Simon, J. B., & Engelhard, E. 2010b, ApJ, 713, 1134
- (63) Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. 2012, ApJ, 746, id. 110

## Appendix A Characteristic disk parameters

For our adopted opacity behavior (12)-(13), assuming that fragmentation happens at , the fiducial values of the surface density , midplane temperature and sound speed and , mass accretion rate , and angular frequency are

(A1) | |||||

(A2) | |||||

(A3) |

The numerical estimates are for . More general expressions for these variables for different opacity behavior () and , as well as the coefficient in the expression (18) for , can be found in Rafikov (2009).