ON THE DYNAMICS AND EVOLUTION OF GRAVITATIONAL INSTABILITY-DOMINATED DISKS

We begin from the equations describing the evolution of a viscous fluid in a gravitational field that is in the process of turning its mass into collisionless stars. These are the equation of continuity, the Navier–Stokes equation, and the first law of thermodynamics:

$$\begin{eqnarray} \frac{1}{\rho } \frac{D\rho }{Dt} & = & -\nabla \cdot \mathbf {v}-\frac{\dot{\rho }_*}{\rho },\qquad\quad\ \ \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \rho \frac{D\mathbf {v}}{Dt} & = & - \nabla p - \rho \nabla \psi + \nabla \cdot \mathbf {T}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \rho \frac{De}{Dt} & = & -p\nabla \cdot \mathbf {v}+ \Phi + \Gamma -\Lambda. \end{eqnarray} \tag{ 3 }$$

Here ρ, v, e, and p are the gas density, velocity, specific internal energy, and pressure, respectively, $\dot{\rho }_*$ is the rate per unit volume at which gas mass turns into stellar mass, ψ is the gravitational potential, T is the viscous stress tensor, Φ = T^ij(∂v_i/∂x_j) is the dissipation function, and Γ and Λ are the volumetric rates of energy gain and loss due to non-fluid (e.g., radiative or chemical) processes, respectively. Note that no terms associated with star formation appear in the first law of thermodynamics or the Navier–Stokes equations because star formation does not alter the bulk velocity or specific internal energy of the gas.

We consider a thin, axisymmetric disk centered on the origin lying in the plane z = 0. At every radius r, the disk is characterized by a surface density Σ and a total thermal plus non-thermal velocity dispersion σ. The material orbits the origin with the angular velocity v_ϕ and has a radial velocity v_r ≪ v_ϕ. In Appendix A, we show that for such a star-forming disk Equations (1)–(3) imply

$$\begin{eqnarray} \frac{\partial }{\partial t} \Sigma + \frac{1}{r}\frac{\partial }{\partial r}(r\Sigma v_r) & = & - \dot{\Sigma }_*,\quad\ \ \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \Sigma \left(\frac{\partial j}{\partial t} + v_r\frac{\partial }{\partial r}j\right) & = & \frac{1}{2\pi r} \frac{\partial }{\partial r} \mathcal {T}, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &&\frac{1}{2}\Sigma \left[\frac{\partial }{\partial t}\left(v_{\phi }^2+3\sigma ^2\right) + v_r\frac{\partial }{\partial r}\left(v_{\phi }^2+3\sigma ^2+2\psi \right)\right] \nonumber\\ &&\quad + \frac{1}{r} \frac{\partial }{\partial r} (r \Sigma v_r \sigma ^2) = \frac{1}{2 \pi r}\frac{\partial }{\partial r}(\Omega \mathcal {T}) + \mathcal {G}- \mathcal {L}, \end{eqnarray} \tag{ 6 }$$

where

$$\begin{equation} \mathcal {T}= 2\pi \int r^2 T_{r\phi } \, dz \end{equation} \tag{ 7 }$$

is the viscous torque, j = rv_ϕ is the specific angular momentum, Ω = v_ϕ/r is the angular velocity, β = ∂ln v_ϕ/∂ln r, and $\mathcal {G}= \int \Gamma \, dz$ and $\mathcal {L}=\int \Lambda \, dz$ are the vertically integrated rates of non-fluid energy gain and loss, respectively. Equations (4)–(6) are the standard equations of mass, angular momentum, and energy conservation, respectively, for a thin disk (e.g., Balbus & Hawley 1998) generalized to the case of a supersonically turbulent gas that is forming stars.

If we assume that the disk is always close to radial force balance and that the potential varies slowly in time then we have ∂ψ/∂r ≈ v²_ϕ/r and ∂j/∂t ≈ 0. Using these conditions in Equations (4)–(6), we arrive at the evolution equations for Σ and σ:

$$\begin{eqnarray} \frac{\partial \Sigma }{\partial t} & = & \frac{1}{2\pi (\beta +1) r v_\phi } \left[\frac{\beta (\beta +1) + r\frac{\partial \beta }{\partial r}}{(\beta +1) r} \left(\frac{\partial \mathcal {T}}{\partial r}\right) - \frac{\partial ^2 \mathcal {T}}{\partial r^2}\right]- \dot{\Sigma }_*, \nonumber \\ \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \frac{\partial \sigma }{\partial t} & = & \frac{\mathcal {G}- \mathcal {L}}{3\sigma \Sigma } + \frac{1}{6\pi r \Sigma } \left[ (\beta -1) \frac{v_\phi }{r^2 \sigma } \mathcal {T}\right. \nonumber \\ & & {} +\frac{\beta ^2\sigma + \sigma \left(r \frac{d\beta }{dr}+\beta \right) - 5 (\beta +1) r \frac{\partial \sigma }{\partial r} }{(\beta +1)^2 r v_\phi } \left(\frac{\partial \mathcal {T}}{\partial r}\right) \nonumber \\ & & \left. -\frac{\sigma }{(\beta +1)v_\phi } \left(\frac{\partial ^2\mathcal {T}}{\partial r^2}\right) \right]. \end{eqnarray} \tag{ 9 }$$

Finally, we note that the underlying assumption of our model is that we can represent the transport processes in a disk dominated by gravitational instability using a local viscosity. This assumption is only valid in certain circumstances, and this sets limits on the applicability of our model that we discuss in Section 5.5.

Equations (8) and (9) fully specify the time evolution of the system. The only physical approximations we have made thus far are that the disk is thin and axisymmetric, and that turbulent eddies on size scales of the disk scale height provide an effective pressure proportional to the square of the turbulent velocity dispersion. However, we have not yet determined the functions describing the star formation rate $\dot{\Sigma }_*$, the rates of radiative gain and loss $\mathcal {G}$ and $\mathcal {L}$, and the torque $\mathcal {T}$. Since the physics involved in these terms is complex, we proceed in a simple parameterized manner.

For the star formation rate, we note that both observation (e.g., Krumholz & Tan 2007; Bigiel et al. 2008; Evans et al. 2009) and theory (e.g., Krumholz & McKee 2005; Krumholz et al. 2009b; Padoan & Nordlund 2010; Hennebelle & Chabrier 2009) indicate that molecular gas forms stars at a rate of ∼1% of its mass per free-fall time. We compute the free-fall time using the mid-plane density; since the scale height is H = σ/Ω, this is ρ = ΣΩ/σ. Thus, we adopt a star formation rate

$$\begin{equation} \dot{\Sigma }_* = \epsilon _{\rm ff}\Sigma \sqrt{G\rho } = \epsilon _{\rm ff} \left(\frac{G \Sigma ^3 \Omega }{\sigma }\right)^{1/2}, \end{equation} \tag{ 10 }$$

where _ff ≈ 0.01. The true value of _ff is likely to be slightly higher than this, because the clouds where stars form are generally denser than the mean mid-plane density in the disk. Even with this correction, however, we can be confident that _ff ≲ 0.1. Moreover, if a significant fraction of the interstellar medium (ISM) is atomic rather than molecular, the star formation rate is greatly reduced (e.g., Wyder et al. 2009; Blanc et al. 2009), in which case _ff can be much smaller. More complex and accurate star formation laws are possible (e.g., Krumholz et al. 2009b), but we will see below that such increased accuracy is not necessary at this stage. However, we do note that for a gaseous disk with Q = 1 (see Section 2.3, Equation (13)) and a flat rotation curve (β = 0), Equation (10) gives a star formation rate $\dot{\Sigma }_* = 0.0067 \Sigma \Omega$; in comparison, Kennicutt (1998) reports an observed star formation rate $\dot{\Sigma }_* = 0.011 \Sigma \Omega$, identical to within the errors.⁵

For the energy loss rate, we note that numerous simulations of turbulence show that it decays via radiative shocks on roughly a crossing timescale (e.g., Stone et al. 1998; Mac Low et al. 1998). For the purpose of estimating the loss rate, we take the characteristic length scale of the turbulence to be comparable to the gas scale height H, so the crossing time is 1/Ω. This is consistent with observations that show that turbulent power is dominated by the largest scales. We also limit our attention to galaxies where the total velocity dispersion is significantly in excess of the thermal value, since these are the only galaxies where one needs to explain the observed velocity dispersion by appealing to physics other than radiative balance. Thus, we take the kinetic energy per unit area to be (3/2)Σσ². The condition that the disk lose this amount of energy per crossing time of the disk scale height then reduces to

$$\begin{equation} \mathcal {L}= \eta \Sigma \sigma ^2\Omega, \end{equation} \tag{ 11 }$$

where η is a dimensionless number of order unity. As a fiducial parameter we adopt η = 3/2, which corresponds to radiating away the full kinetic energy every scale height-crossing time. In disks where the velocity dispersion is primarily thermal rather than non-thermal, the loss rate will assume a different functional form (e.g., Rafikov 2009), but the remainder of our analysis will be unchanged.

The gain rate is much more complex, since it involves turbulent motions generated by star formation. Since we are interested in systems where gravitationally driven turbulence dominates, however, we make the extreme assumption that $\mathcal {G}= 0$. We return to the question of the real value of $\mathcal {G}$ in Section 5.3.

We now turn to the central hypothesis of our model, which is that a self-gravitating disk will adjust its torque, and therefore its radial mass flux, so as to remain in a state of marginal stability. This hypothesis has also been investigated in the models of Burkert et al. (2010). For a disk of gas plus stars, the parameter that describes its stability is (Rafikov 2001)

$$\begin{equation} Q(q)^{-1} = 2Q_*^{-1} \frac{1}{q} \left[1 - e^{-q^2}I_0(q^2)\right] + 2Q_g^{-1} R \frac{q}{1 + q^2 R^2}, \end{equation} \tag{ 12 }$$

where

$$\begin{equation} Q_* = \frac{\kappa \sigma _*}{\pi G \Sigma _*} \qquad Q_g = \frac{\kappa \sigma }{\pi G \Sigma } \qquad R = \frac{\sigma }{\sigma _*}, \end{equation} \tag{ 13 }$$

and I₀ is the Bessel function of order zero. Here Σ_* and σ_* are the surface density and velocity dispersion of stars, respectively, κ = [2(β + 1)]^1/2Ω is the epicyclic frequency, and q = kσ_*/κ is the dimensionless wavenumber of the mode in question. Modes for which Q(q) < 1 are unstable. Note that Romeo et al. (2010) have proposed a generalization of this condition for the case of gas with a scale-dependent velocity dispersion, as expected for turbulence, but for simplicity we use the Rafikov (2001) criterion instead.

It is not generally possible to find the minimum value of Q(q) analytically. However, we can obtain a significant simplification if we focus on the most interesting cases for gravitationally driven turbulence. In disks at high redshift there has not been time for the stars and gas to evolve so that their velocity dispersions are very different—see Section 5.2 for a further discussion of this point. For these disks, we therefore adopt σ_* = σ. In this case, the expression

$$\begin{equation} Q = \min (Q(q)) \approx \left(\frac{1}{Q_g} + \frac{1}{Q_*}\right)^{-1} = \frac{\kappa \sigma }{\pi G (\Sigma +\Sigma _*)} \end{equation} \tag{ 14 }$$

is accurate to better than 7%. With this approximation, the condition for a disk to remain marginally stable at Q = 1 becomes

$$\begin{eqnarray} 0 & = & \frac{\partial Q}{\partial \sigma } \frac{\partial \sigma }{\partial t} + \frac{\partial Q}{\partial \Sigma } \frac{\partial \Sigma }{\partial t} + \frac{\partial Q}{\partial \Sigma _*} \frac{\partial \Sigma _*}{\partial t}\ \ \ \ \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} & = & \frac{1}{\sigma }\frac{\partial \sigma }{\partial t}-\frac{1}{\Sigma +\Sigma _*}\left(\frac{\partial \Sigma }{\partial t} + \frac{\partial \Sigma _*}{\partial t}\right). \end{eqnarray} \tag{ 16 }$$

Plugging in our expressions for the temporal derivatives of σ, Σ, and Σ_*, we obtain an equation that describes the torque required to maintain Q = 1:

$$\begin{equation} f_2 \frac{\partial ^2 \mathcal {T}}{\partial r^2} + f_1 \frac{\partial \mathcal {T}}{\partial r} + f_0 \mathcal {T}= F \end{equation} \tag{ 17 }$$

with

$$\begin{eqnarray} f_0 & = & \left[\frac{\beta -1}{6\sqrt{2(\beta +1)}}\right] \frac{G}{f_g r^2 \sigma ^3}, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} f_1 & = & -\frac{(3 f_g-1) r \sigma \frac{d\beta }{dr} - (\beta +1)\left[(3 f_g-1) \beta \sigma +5 r \frac{\partial \sigma }{\partial r}\right]}{6\sqrt{2(\beta +1)^5}},\nonumber\\ & &\times \frac{G}{f_g r v_\phi ^2 \sigma } \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} f_2 & = & \left[\frac{3f_g-1}{6\sqrt{2(\beta +1)^3}}\right] \frac{G}{f_g v_\phi ^2\sigma }, \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} F & = & \frac{\eta v_\phi }{3r}, \end{eqnarray} \tag{ 21 }$$

where f_g = Σ/(Σ + Σ_*) is the gas fraction in the disk, and we have used Q = 1 to replace any dependence on Σ with a dependence on σ and f_g. We refer to this as the torque equation. Note that in deriving Equation (17) we have implicitly assumed that stars do not migrate radially from their formation locations. If we were to make the opposite assumption that stars and gas move together, then the combined gas plus star disk would act essentially like a purely gaseous one, except that the stars would be dissipation free. The corresponding torque equation is simply Equation (17) with η replaced by ηf_g, and f_g → 1 in all other terms.

The other interesting location to consider for gravitationally driven turbulence is in outer parts of present-day disks. In these regions, star formation occurs at a negligible rate, and σ_* ≫ σ. In this case Q ≈ Q_g, and a little calculation shows that the torque equation is simply (17) with f_g = 1 everywhere. The stars simply become irrelevant for gravitational stability. Thus, Equation (17), with the appropriate choice of f_g and η, is capable of representing both a present-day outer galaxy disk and a high-redshift disk.

To finish specifying the torque, we must provide boundary conditions for Equation (17). One boundary condition comes from fixing the external accretion rate onto the galaxy to a value $\dot{M}=\dot{M}_{\rm ext}$ (so $\partial \mathcal {T}/\partial r = -v_\phi (\beta +1) \dot{M}_{\rm ext}$) at the disk outer edge, r = R. We imagine the inflow rate at this radius to be set by cosmological infall, which occurs at a rate unaffected by what happens in the disk. The other boundary condition depends on how we handle the inner boundary. If we imagine that our model applies all the way to r = 0, we must find a solution that remains regular in the vicinity of the singular point there. We show below that there does exist a steady-state solution that is regular at r = 0 and has the specified accretion rate $\dot{M}=\dot{M}_{\rm ext}$ at r = R. Outside of that steady state it is not possible to simultaneously have regularity at r = 0 and an externally imposed accretion rate at r = R, so we must instead truncate the model at some radius r₀>0, where we imagine that a stellar bulge or some other non-disky structure forms. In that case, we require that the torque have a small value at r = r₀, so that the inner bulge region does not do work on the disk.

The evolution of the system is now fully specified. At any given time, Equation (17) specifies the viscous torque. That torque in turn sets the time evolution of Σ, σ, and Σ_* via Equations (8), (9), and (10), respectively. Before moving on, however, we pause to point out some of the important physical properties embodied in Equation (17). First, the star formation rate does not appear explicitly in the torque equation. This is because for σ_* ≈ σ changing gas into stars does not significantly affect the stability of the disk, and with our approximate form for Q it does not affect the stability at all. Star formation enters the problem solely through its effects on f_g. Second, if η = 0 then clearly $\mathcal {T}= 0$ is a solution for $\dot{M}_{\rm ext} = 0$. Physically, this represents the fact that if there is neither dissipation of turbulence nor accretion, then the disk can remain marginally stable without any mass transport.

It is helpful at this point to non-dimensionalize our equations and derive some characteristic numbers. If we make a change of variables r = xR, σ = sv_ϕ(R), v_ϕ = uv_ϕ(R), and $\mathcal {T}= \tau \dot{M}_{\rm ext} v_\phi (R) R$, then we can rewrite Equation (17) as

$$\begin{equation} \tau ^{\prime \prime } + h_1 \tau ^{\prime } + h_0 \tau = H, \end{equation} \tag{ 22 }$$

with

$$\begin{eqnarray} h_0 & = & \left(\frac{\beta ^2 -1}{3f_g-1}\right) \frac{u^2}{x^2 s^2} \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} h_1 & = & -\frac{5(\beta +1) x s^{\prime }+(3f_g-1)s(\beta +\beta ^2+x\beta ^{\prime })}{(3f_g-1)(\beta +1) s x} \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} H & = & \frac{\eta }{\chi } \left(\frac{2 f_g \sqrt{2(\beta +1)^3}}{3f_g-1}\right) \frac{s u^3}{x}, \end{eqnarray} \tag{ 25 }$$

subject to the boundary condition τ' = −β − 1 at x = 1. Here primes indicate differentiation with respect to x, and we have defined

$$\begin{equation} \chi = \frac{G \dot{M}_{\rm ext}}{v_\phi (R)^3}. \end{equation} \tag{ 26 }$$

The form of Equation (22) is instructive. The coefficients h₀ and h₁ appearing on the left-hand side depend on the current state of the disk without reference to external accretion or turbulent dissipation. These affect only the inhomogeneous term H on the right-hand side, which is proportional to η/χ. We may view H as the driving term for the system, with more rapid dissipation of turbulence (i.e., larger η) and lower accretion rates (lower χ) both tending to produce larger torques.

We can similarly non-dimensionalize the evolution equation for the gas surface density and velocity dispersion by defining $\Sigma = S \dot{M}_{\rm ext}/(v_\phi (R) R)$ and t = T[2πR/v_ϕ(R)], so that the evolution equations become

$$\begin{eqnarray} \frac{\partial S}{\partial T} & = & \frac{(\beta ^2+\beta +x \beta ^{\prime })\tau ^{\prime }-x(\beta +1)\tau ^{\prime \prime }}{(\beta +1)^2 u x^2} -\frac{dS_*}{dT},\quad \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} \frac{\partial s}{\partial T} & = & \frac{1}{3 (\beta +1)^2 s S u x^3} \lbrace u^2 (\beta +1)^2(\beta -1) \tau \nonumber\\ & & + s x [s(\beta +\beta ^2+x\beta ^{\prime })-5 (\beta +1) x s^{\prime }] \tau ^{\prime } \nonumber\\ & & - (\beta +1) s^2 x^2 \tau ^{\prime \prime } - 2 \pi (\beta +1)^2 \eta s^2 S u^2 x^2 \rbrace, \end{eqnarray} \tag{ 28 }$$

where

$$\begin{equation} \frac{dS_{*}}{dT} = 2\pi \epsilon _{\rm ff} \sqrt{\frac{u S^{3} \chi }{s x}}. \end{equation} \tag{ 29 }$$

We have defined $\Sigma _{*} = S_{*} \dot{M}_{\rm ext}/(v_\phi (R) R)$ in analogy with S, and the Q = 1 condition in the dimensionless form is

$$\begin{equation} \frac{\sqrt{2(\beta +1)} u s}{\pi \chi x (S+S_{*})} = 1. \end{equation} \tag{ 30 }$$

In these units the orbital period is 1, and the fraction of the disk mass that accretion adds per orbit is 〈S + S_*〉⁻¹, where the angle brackets indicate an average over the disk area.

Having derived the basic equations that govern the system, we now search for steady-state solutions, which we can obtain analytically. A true steady state is obviously not possible in a real disk that undergoes mass accretion and star formation, but we can find solutions which are steady for time periods that are short compared to the star formation and accretion timescales. We will therefore set _ff = 0, and in Section 3.2 we will check that this is a reasonable approximation. In this case, a steady-state solution is one for which $\partial \dot{M}/\partial r = 0$, or in the dimensionless form

$$\begin{equation} \frac{d}{dx} \left[\frac{\tau ^{\prime }}{(\beta +1)u}\right] = 0. \end{equation} \tag{ 31 }$$

Combined with the boundary condition that τ' = −(β + 1) at x = 1, this implies that the steady-state solution is τ' = −(β + 1)u. One can immediately verify using Equation (27) that such a torque gives ∂S/∂T = 0.

To make further progress toward an analytic solution, we concentrate our attention on cases where β has a constant value. This is a reasonable limitation, since most galaxies have flat rotation curves (β = 0), while Keplerian disks have β = −1/2. For constant β, we have u = x^β, and we can analytically integrate the steady solution for τ' to obtain τ = −x^β+1 + c, where c is a constant to be determined by the regularity condition at the inner boundary. In Appendix B, we provide this analysis for two of the most physically important cases, β = 0 (flat rotation) and β = −1/2 (Keplerian rotation), which shows that c = 0 in both the cases. Of course, β cannot have a constant value of 0 or −1/2 all the way to r = 0, since the rotation velocity would diverge, but our approximation is appropriate in cases where β has a constant value over a large dynamic range in radius, and changes only at small radii where a bulge forms (in the case of a flat rotation-curve galaxy) or a boundary layer joins a disk to a star (for a Keplerian disk).

Given the value that τ must have in order to produce a steady state, we can plug it into the torque equation to determine the corresponding disk properties that are required. For a flat rotation curve, β = 0, we have

$$\begin{eqnarray} &&\tau ^{\prime \prime } \,{-}\, \left(\!\frac{5}{3 f_g - 1}\!\right) \frac{s^{\prime }}{s} \tau ^{\prime } - \frac{1}{(3f_g-1) s^2 x^2} \tau \,{=}\, \left(\!\frac{2^{3/2}f_g }{3 f_g-1}\!\right)\! \left(\!\frac{\eta }{\chi }\!\right) \frac{s}{x}.\nonumber\\ \end{eqnarray} \tag{ 32 }$$

and the steady-state condition τ = −x then implies

$$\begin{equation} s^{\prime } = \frac{2\sqrt{2} f_g \eta s^3 - \chi }{5 s x \chi }. \end{equation} \tag{ 33 }$$

Clearly,

$$\begin{equation} s = \frac{1}{\sqrt{2}}\left(\frac{\chi }{\eta f_g}\right)^{1/3} \quad \mbox{ or }\quad \sigma = \frac{1}{\sqrt{2}}\left(\frac{G\dot{M}_{\rm ext}}{\eta f_g}\right)^{1/3} \end{equation} \tag{ 34 }$$

is an exact solution, and numerical integration shows that all solutions converge to this value very quickly at x < 1 regardless of the value of s at x = 1. The corresponding surface density and inward velocity of the material are, from Equations (A2) and (14),

$$\begin{eqnarray} \Sigma & = & \frac{v_\phi }{\pi G r} \left(\frac{f_g^2 G \dot{M}_{\rm ext}}{\eta }\right)^{1/3}, \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} v_r & = & -\eta \frac{\sigma ^2}{v_\phi },\qquad\qquad\ \ \end{eqnarray} \tag{ 36 }$$

and the corresponding dimensionless viscosity parameter (Shakura & Sunyaev 1973) and viscous evolution timescales are

$$\begin{eqnarray} \alpha & = & \frac{G \dot{M}}{3 \sigma ^3} = \frac{2\sqrt{2}}{3}\eta f_g,\qquad \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} t_{\rm visc} & = & \frac{R}{v_r(R)} = \left(\frac{f_g^2}{\eta \chi ^2}\right)^{1/3} \frac{t_{\rm orb}}{\pi }, \end{eqnarray} \tag{ 38 }$$

where t_orb = 2πR/v_ϕ(R) is the outer disk orbital period. (We omit a factor of Q in the equation for α because it is set to unity in our model.) For the corresponding Keplerian case (β = −1/2), it is easy to verify by a similar procedure, and using the analysis of the singular point provided in Appendix B, that the steady solution is $\tau =-\sqrt{x}$, s = [3χ/(4ηf_g)]^1/3.

Our result is easy to understand intuitively. For σ ≫ c_s, the rate per unit mass at which the turbulence decays is ησ²Ω, so the turbulent decay time is of order the orbital timescale times η. To keep the velocity dispersion constant, mass must move inward at a rate such that the decrease in gravitational potential energy balances this radiative loss. For a flat rotation curve, the decrease in potential involved in moving from radius r₀ to radius r is v²_ϕln(r/r₀), so the inward drift at a velocity v_r causes the potential energy per unit mass to decrease at a rate v²_ϕ(v_r/r). Equating the rates of dissipation and energy increase gives ησ²Ω = v²_ϕ(v_r/r), and it follows immediately that v_r = ησ²/v_ϕ.

It is also worth noting that this result is very similar to that of Gammie (2001), who finds that, in steady state in a disk that cools on a timescale τ_c, gravitationally induced turbulence produces an effective viscosity α = [γ(γ − 1)(9/4)Ωτ_c]⁻¹. Our effective "cooling time" for the supersonic turbulence is τ_c = 1/(ηΩ), and the remainder of our result differs from his only in that we have modeled disks with a stellar component, and that our pressure and internal energy are appropriate for supersonically turbulent motion on scales comparable to the disk scale height, rather than for an adiabatic gas described by a polytropic equation of state. The former introduces a dependence on f_g, and the latter produces a slight change in leading coefficient. That Gammie's result can be obtained on energetic grounds, rather than by computing stresses as in his derivation, has also been pointed out by Rice et al. (2005) and Rafikov (2009).

In deriving the steady solution, we have ignored the change in the total disk mass due to accretion. More subtly, our steady solution has constant $\dot{M}$ all the way in to r = 0, so the mass effectively vanishes through the origin. These approximations are only reasonable for timescales over which the total disk mass changes a little. We can define the accretion timescale for our steady solution as

$$\begin{equation} t_{\rm acc} = \frac{\int _0^R 2 \pi r (\Sigma +\Sigma _*)\, dr}{\dot{M}_{\rm ext}} = \frac{t_{\rm orb}}{\pi (\eta f_g \chi ^2)^{1/3}} = \frac{t_{\rm orb}}{2\pi \eta f_g s^2}. \end{equation} \tag{ 39 }$$

Note that t_acc is the time required for external accretion to double the disk mass, and is distinct from the viscous accretion time defined by Equation (38). To avoid confusion, we will always refer to that quantity as the viscous timescale and use the accretion time solely to refer to the mass-doubling time produced by external infall. For our fiducial η = 3/2, we have

$$\begin{equation} t_{\rm acc} \simeq 1.3 f_g^{-1/3} \chi _{0.1}^{-2/3} t_{\rm orb}, \end{equation} \tag{ 40 }$$

where χ_0.1 = χ/0.1. Similarly, we have neglected star formation, which is only reasonable for time short compared to the star formation timescale. We can define the star formation timescale as the mean ratio of star formation rate to gas surface density in the disk:

$$\begin{equation} t_{\rm SF} = \frac{ t_{\rm orb}}{\pi } \int _0^1 (2\pi x) \frac{S}{dS_*/dT}\, dx = \frac{t_{\rm orb}}{(162\pi ^2)^{1/4} f_g^{1/2}\epsilon _{\rm ff}}. \end{equation} \tag{ 41 }$$

For our fiducial _ff = 0.01, this gives

$$\begin{equation} t_{\rm SF} \simeq 16 f_g^{-1/2} t_{\rm orb}. \end{equation} \tag{ 42 }$$

In comparison, the characteristic evolution timescale for the disk should be the viscous time defined by Equation (38), since this is the time required to drain the material in the disk and replace it with newly accreted material. For our fiducial η = 3/2 and a flat rotation curve, this is

$$\begin{equation} t_{\rm visc} = 1.3 f_g^{2/3} \chi _{0.1}^{-2/3} t_{\rm orb}. \end{equation} \tag{ 43 }$$

Comparing the viscous time to the timescales relevant for accretion and star formation, we have

$$\begin{eqnarray} \frac{t_{\rm visc}}{t_{\rm acc}} & = & f_g,\qquad\qquad\qquad\qquad\quad\qquad\ \ \ \ \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray} \frac{t_{\rm visc}}{t_{\rm SF}} & = & \left(\frac{162}{\pi ^2}\right)^{1/4} \frac{f_g^{7/6} \epsilon _{\rm ff}}{\eta ^{1/3}\chi ^{2/3}} = 0.08 f_g^{7/6} \chi _{0.1}^{-2/3} \end{eqnarray} \tag{ 45 }$$

for our fiducial values of η and _ff. Our neglect of accretion-induced changes in the disk mass and star formation-induced changes in the gas fraction is reasonable when these two ratios are ≲1. Clearly, the requirement that t_visc/t_acc ≲ 1 is always satisfied, although perhaps only marginally if f_g is large. The requirement that t_visc/t_SF ≲ 1 is satisfied as long as χ ≳ 10⁻³, or longer if the disk is non-star-forming (_ff = 0). Thus, we expect that our steady-state model is reasonable under these conditions.

We have now shown that disks dominated by gravitationally driven turbulence rapidly converge to an equilibrium state in which their velocity dispersions are determined by their gas fractions and accretion rates. Since this convergence happens on an orbital timescale, most galactic disks should be found near their equilibrium state. We can use this result, coupled to a simple model for how galaxy halos accrete mass, to study the evolution of disk velocity dispersions over cosmic time. Using Press–Schechter fits to dark matter simulations, Neistein & Dekel (2008) estimate the mean dark matter accretion rate onto halos of a given mass at a given redshift. Bouche et al. (2010) extend this to give an estimate of the gas accretion rate onto the disk at the center of the halo, which we adopt

$$\begin{equation} \dot{M}_g = 7.0\epsilon _{\rm in} f_{b,0.18} M_{h,12}^{1.1} (1+z)^{2.2}\, M_{\odot }\mbox{ yr}^{-1}, \end{equation} \tag{ 46 }$$

where z is the redshift, f_b,0.18 is the gas mass fraction of the infall divided by 0.18, the universal baryon fraction, M_h,12 is the halo mass in units of 10¹² M_☉, and _in is the fraction of gas entering the halo that reaches the galactic disk rather than being shock-heated and joining the halo. This is approximately given by

$$\begin{equation} \epsilon _{\rm in} = \left\lbrace \begin{array}{@{}ll@{}} 0.7 f(z), \quad & M_{h,12} < 1.5 \\ 0, & M_{h,12} > 1.5, \end{array} \right. \end{equation} \tag{ 47 }$$

where f(z) is a function that is linear in time and varies from unity at z = 2.2 to 0.5 and at z = 1.⁸ Inserting $\dot{M}_g$ from Equation (46) for $\dot{M}_{\rm ext}$ into Equation (34), we are able to evaluate the expected velocity dispersion of gravitational instability-dominated galactic disks as a function of halo mass and redshift. We do so in Figure 2.

Zoom In Zoom Out Reset image size

Examining the plot, we see that for a Milky Way-like halo (M_h,12 = 1; Xue et al. 2008), where σ_* ≫ σ (so that the f_g = 1 case applies), we predict a typical velocity dispersion of 10.7 km s⁻¹. While this is very slightly higher than the value of σ ≃ 8 km s⁻¹ observed in typical Milky Way-like disks today (Blitz & Rosolowsky 2004, Dib et al. 2006, and references therein), the agreement is quite good given our purely analytic model. Our results are quite similar to the numerical ones obtained by Kim & Ostriker (2007) and Agertz et al. (2009). Moreover, as Agertz et al. point out, gravitationally driven turbulence has the advantage that it can operate even in the outer H i disk where there is very little star formation, so mechanisms such as supernovae that are invoked to explain turbulence in the inner disk (e.g., de Avillez & Breitschwerdt 2007; Joung et al. 2009) are unavailable. We also predict lower velocity dispersion in smaller halos, and this appears to be consistent with the somewhat lower H i velocity dispersions seen in dwarf galaxies (Walter et al. 2008; Chung et al. 2009). Finally, however, we do note that there are alternative models to explain outer disk turbulence, including magnetorotational instability (Sellwood & Balbus 1999; Piontek & Ostriker 2007) and accretion of clumpy gas (Santillán et al. 2007).

Within the same framework we are able to explain the large velocity dispersions of 20–80 km s⁻¹ found in galactic disks found at redshifts ∼1.5–3 (Cresci et al. 2009). The observed galaxies likely correspond to ∼10¹² M_☉ halos. For redshifts in this range and f_g ∼ 1/2, typical of galaxies at that redshift (Daddi et al. 2010; Tacconi et al. 2010), we predict typical velocity dispersions of 30–50 km s⁻¹, with fluctuations at the factor of ∼1.5 level, corresponding to the expected factor of ∼3 level variations in the accretion rates of halos at the same mass and redshift. This is in good agreement with the observations.

It is also instructive to compare the velocity dispersions we predict to the expected rotation velocities of galactic disks. We compute the approximate virial velocity of a halo as a function of mass and redshift following the approximation given in Appendix A2 of Dekel & Birnboim (2006), and we take the maximum circular velocity to be 1.2 times this based on fitting the zero point of the Tully–Fisher relation (Dutton et al. 2007). With this approximation, we plot V_max/σ in Figure 3. We see that accretion-driven turbulence naturally produces the transition from disks with V_max/σ ∼ 5 found at redshifts ≳2 to disks with V_max/σ ∼ 20–25 found today.

Zoom In Zoom Out Reset image size

Interestingly, we find that there is little dependence of V_max/σ on halo mass. Instead, the primary dependence is on f_g, the gas mass fraction; analytically, V_max/σ ∝ f^−1/3_g. Thus, the most gas-dominated systems (or old galaxies that have σ_* ≫ σ) have the largest V_max/σ, while gas-poor systems have smaller V_max/σ. This suggests that the range of V_max/σ seen for galaxies at z ∼ 2 by the SINS survey (Cresci et al. 2009) represents a sequence in gas fraction. The dispersion-dominated galaxies should on average be comparatively gas-poor, while rotation-dominated ones should be gas-rich. Of course, fluctuations in the accretion rate can also cause changes in σ, so detecting this effect will require samples large enough for this noise source to be averaged out. Nonetheless, it seems likely that data to test this prediction will become available in the next few years.

It is particularly interesting to apply our models to the z ∼ 2–3 galaxies observed by Elmegreen et al. (2004, 2005), Genzel et al. (2008), Förster Schreiber et al. (2009), Cresci et al. (2009), and others, since these are thought to be examples of strongly gravitational instability-dominated disks. We first note that, in the redshift range z = 2–3 for halos of mass M_h = 10¹² M_☉, thought to be typical of the observed systems, the models shown in Figures 2 and 3 give χ = 8 × 10⁻³–1.1 × 10⁻². For these values of χ and gas fractions f_g = 1/2, using Equations (40) and (42), the ratio of star formation time to accretion time t_SF/t_acc = 1.8–2.6, so the star formation rate is roughly 1/2–1/3 of the total accretion rate. Given the uncertainties in this model and the dispersion in expected accretion rates, for simplicity we can simply adopt $\dot{M}_* \approx \dot{M}_{\rm ext}$. Since the star formation rates are observed (and have typical values ∼100 M_☉ yr⁻¹), we can plug them into our model in place of $\dot{M}_{\rm ext}$ in order to predict disk properties.

Doing so, we find that the redshift 2–3 disks should have velocity dispersions (from Equation (34))

$$\begin{equation} \sigma \approx 47\mbox{ km s}^{-1}\, f_g^{-1/3} \dot{M}_{*,100}^{1/3}, \end{equation} \tag{ 48 }$$

where $\dot{M}_{*,100} = \dot{M}_* / 100$ M_☉ yr⁻¹, independent of their maximum rotation velocities V_max. Thus, galaxies of the similar star formation rate and gas fraction should have the same σ independent of V_max. The viscous accretion timescale required for the gas at the edge of one of these disks to reach the center is (from Equation (38))

$$\begin{equation} t_{\rm visc} \approx 600\mbox{ Myr } f_g^{2/3} R_{10} V_{200}^{-1} \dot{M}_{*,100}^{-2/3}, \end{equation} \tag{ 49 }$$

where R₁₀ = R/10 kpc and V₂₀₀ = V_max/200 km s⁻¹, and the gas mass is (from Equation (35))

$$\begin{eqnarray} M_g & = & R v_{\phi } \left(\frac{f_g^2 \dot{M}_*}{\eta }\right)^{1/3}\nonumber\\ & = & 3\times 10^{10}\,M_{\odot }\; f_g^{2/3} R_{10} V_{200} \dot{M}_{*,100}^{1/3}. \end{eqnarray} \tag{ 50 }$$

The ratio of baryonic to dynamical mass within the disk region R is

$$\begin{equation} \frac{M_{\rm bar}}{M_{\rm dyn}} = 0.3\, f_g^{-1/3} V_{200}^{-1} \dot{M}_{*,100}^{1/2}. \end{equation} \tag{ 51 }$$

To the extent that these quantities have been observed, they are in good agreement with the results of our model.

It is also useful to verify that our approximation that σ_* ≈ σ is valid for these galaxies. Once stars form, transient spiral structures will dynamically heat them until the stellar disk becomes stable against further spiral patterns (Sellwood & Carlberg 1984; Carlberg & Sellwood 1985). The characteristic timescale for this heating is [(Q_lim − Q_*)/τ]t_orb, where Q_lim ≈ 2 is the limiting value at which the disk becomes stable against spiral perturbations, Q_* is the current Toomre Q parameter for the stars, and numerical experiments show that τ ∼ 4–5 for t_orb evaluated at one disk scale length. If we assume that the stars are born at Q_* = 1, then the characteristic timescale for heating is

$$\begin{equation} t_{\rm heat} \approx 750\mbox{ Myr } R_{10} V_{200}^{-1}, \end{equation} \tag{ 52 }$$

where we have taken the radial scale length to be half of R. In contrast, the time required to double the stellar mass is

$$\begin{equation} t_{\rm *} = \frac{1-f_g}{f_g}\left(\frac{M_g}{\dot{M}_*}\right) = 300\mbox{ Myr } \frac{1-f_g}{f_g^{1/3}} R_{10} V_{200} \dot{M}_*^{-2/3}. \end{equation} \tag{ 53 }$$

Thus, we see that the time required for stars to increase their velocity dispersion via spiral structure is generally comparable to or longer than the time required for a new generation of stars to form with the same velocity dispersion as the gas. Our approximation that the stellar population has the same velocity dispersion as the gas in these galaxies is therefore reasonable.

In our idealized models, we have neglected the influence of stellar feedback by setting $\mathcal {G}= 0$. This is obviously reasonable if we are concerned with the outer parts of a galactic disk where there is no star formation, or a protostellar disk where at most a few stars will form. It is not reasonable for the centers of present-day galactic disks, where supernovae are clearly important. It is questionable whether stellar feedback is important in ultra-luminous infrared galaxies (ULIRGs) or the high surface density galaxies found in the early universe. Supernovae are not effective in such environments (Thompson et al. 2005; Joung et al. 2009), but stellar radiation pressure may be. Whether it can actually drive the observed velocity dispersions in high-redshift galaxies is a matter of debate (Murray et al. 2010; Krumholz & Matzner 2009).

In those situations where feedback is significant, we can qualitatively see how it would change our results by noting that adding a non-zero $\mathcal {G}$ to our equations would have roughly the same effect as lowering η. Physically, if star formation injects turbulence into the ISM at an appreciable rate, this is equivalent to reducing the rate at which turbulence decays—we effectively increase the "cooling time" of the disk. Consulting Equations (34)–(37), we see that the effect of this is to increase the velocity dispersion and surface density in the equilibrium state, while reducing the radial velocity and the rate of angular momentum transport.

We caution that this analysis is only valid as long as the feedback is not too strong. In particular, we require that $\mathcal {L}$ remains larger than $\mathcal {G}$ for a Q = 1 disk, and that the turbulent stresses created by the feedback mechanism are significantly weaker than the stresses induced by gravitational instability-driven turbulence. If the first requirement is not met, then feedback will drive the velocity dispersion up to the point where Q>1, and the gravitational instability will shut off. If the latter condition fails, then gravitational instability will continue, but our calculation of the transport rate will not be correct because we have not included stresses induced by feedback. Even if both requirements are met, our analysis of feedback effects should be regarded as qualitative rather than quantitative. Energy injection $\mathcal {G}$ appears on the right-hand side of the torque equation with the opposite sign as $\mathcal {L}$, but their functional dependence on other disk parameters (surface density, velocity dispersion, etc.) is almost certainly different. The exact effects of feedback will depend on how energy injection varies with these quantities, which will in turn depend on the type of feedback and the physics of the ISM.

Although we have focused our discussion thus far on galactic disks, our model applies for arbitrary rotation curves, gas fractions, and infall rates, so we can apply it equally well to protostellar disks. To understand the expected levels of turbulence in protostellar disks, we use the parameterization of infall due to Kratter et al. (2008, 2010), who introduce the dimensionless numbers

$$\begin{equation} \xi = \frac{G \dot{M}_{\rm ext}}{c_{s,d}^3} \qquad \Gamma = \frac{\dot{M}_{\rm ext}}{M_{*d} \Omega _{k,\rm in}}, \end{equation} \tag{ 54 }$$

where c_s,d is the sound speed in the disk, M_*d is the total mass of the disk and star, and $\Omega _{k,\rm in}$ is the Keplerian orbital period of the infalling material. Physically, ξ represents the ratio of the external accretion rate to the maximum rate (∼c³_s,d/G) at which a stable disk can process material, while Γ represents (up to a factor of 2π) the fraction by which the disk plus star mass changes per the outer disk orbit. Indeed, since $v_{\phi }(R) = R\Omega _{k,\rm in} = \sqrt{G M_{*d}/R}$, with a little algebra it is easy to show that, in the case of a Keplerian disk consisting entirely of gas, our χ simply reduces to Kratter et al.'s Γ parameter.

With this understanding, we can explain the observation by Kratter et al. (2010) that, in their simulations, the typical velocity dispersion of disks that do not fragment is comparable to the disk thermal sound speed (see their Figure 8). For a purely gaseous Keplerian disk, our model gives s = [3χ/(4η)]^1/3, and Kratter et al. show that the disk sound speed is related to the Keplerian velocity at the disk edge by c_s,d/v_ϕ(R) = (Γ/ξ)^1/3 (their Equation (18)). Combining these two results, the expected Mach number of the accretion-driven turbulence is

$$\begin{equation} \mathcal {M} = \frac{\sigma }{c_{s,d}} = s \frac{c_{s,d}}{v_\phi (R)} = \left(\frac{3\xi }{4\eta }\right)^{1/3}. \end{equation} \tag{ 55 }$$

Since fragmentation is avoided only for disks with ξ of no more than a few, we can take ξ ∼ 1, and it immediately follows that the expected Mach number $\mathcal {M} \sim 1$.

We can apply a similar analysis to real protostellar disks: the Mach number of the turbulence in these disks should follow Equation (55). This means that disks accreting with ξ ∼ 1, corresponding to $\dot{M}_{\rm ext} \sim 10^{-5}$ M_☉ yr⁻¹ for typical outer disk temperatures T ∼ 50 K, should have disks whose turbulent velocity dispersions are roughly transonic. This state should prevail during the majority of the main accretion phase. Once the main accretion phase ends and the accretion rate drops, the turbulent velocity dispersion should drop to subsonic values. This represents another prediction from our analysis: class 0 and class I protostars should have disks with transonic turbulent velocity dispersions, while class II and class III sources should have subsonic turbulent velocity dispersions. As ALMA comes online in the next few years and provides resolved molecular line maps of protostellar disks at a variety of stages in their evolution (e.g., Krumholz et al. 2007a), we will be able to test this prediction.

The central approximation we make in our model is that transport of mass, angular momentum, and energy produced by gravitationally driven turbulence can be represented by a local viscous stress tensor. The validity of this approximation has been the subject of great debate in the past decade. Balbus & Papaloizou (1999) show that self-gravitating disks cannot in general be modeled with a viscous formalism, but that such an approximation may be reasonable for disks near Q = 1, the condition that we adopt throughout this work and that appears to apply to the galactic and protostellar disks we are interested in studying. Based on a combination of analytic arguments and local simulations, Gammie (2001) argues that a local prescription is applicable to Q = 1 disks that are sufficiently thin, s ≲ 0.12, and more recent global simulations (Lodato & Rice 2004, 2005; Boley et al. 2006; Cossins et al. 2009) generally support this result. Gammie's condition for a local transport approximation to apply is well satisfied for galactic disks at redshifts z ≲ 2 (Section 5.1) and for non-fragmenting protostellar disks (Section 5.4). It is marginally violated for the observed disks at z ∼ 2 (Section 5.2), suggesting that our model should be considered with some caution for them. At a minimum the thickness of these disks likely produces different fragmentation behavior than a standard thin disk analysis would suggest (Begelman & Shlosman 2009).

For β = 0, the torque Equation (22) reduces to

$$\begin{equation} \tau ^{\prime \prime } - \left(\!\frac{5}{3 f_g - 1}\!\right) \frac{s^{\prime }}{s} \tau ^{\prime } - \frac{1}{(3f_g-1) s^2 x^2} \tau = \left(\frac{2^{3/2} f_g}{3 f_g-1}\right) \!\left(\!\frac{\eta }{\chi }\!\right) \frac{1}{x}. \end{equation} \tag{ B1 }$$

We expand τ in a power series about x = 0, τ = ∑^∞_n=0a_nxⁿ, to obtain

$$\begin{eqnarray} &&-\frac{a_0}{(3f_g-1)s^2} x^{-2} - \left[\frac{2\sqrt{2} s^3 f_g \eta +a_1\chi }{(3f_g-1)s^2\chi }\right]x^{-1} \nonumber\\ &&- \left\lbrace \frac{a_2[1-2(3f_g-1)s^2]+5 a_1 f_g s s^{\prime }}{(3f_g-1)s^2}\right\rbrace + O(x) = 0.\qquad \end{eqnarray} \tag{ B2 }$$

Thus, the leading coefficients near x = 0 are

$$\begin{eqnarray} a_0 & = & 0\qquad\qquad\quad\ \ \ \ \end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray} a_1 & = & -2\sqrt{2} s^3 f_g \frac{\eta }{\chi }\quad\ \ \ \end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray} \qquad a_2 & = & \left[\frac{10\sqrt{2} f_g s^4 s^{\prime }}{1-2(3f_g-1) s^2}\right] \frac{\eta }{\chi }. \end{eqnarray} \tag{ B5 }$$

For β = −1/2, the torque Equation (22) reduces to

$$\begin{eqnarray} &&\tau ^{\prime \prime } + \frac{(3f_g-1) s - 10 x s^{\prime }}{(3f_g-1) s x} \tau ^{\prime } - \frac{3}{4(3f_g-1) s^2 x^3} \tau \nonumber\\ &&\quad = \left(\frac{2^{3/2}f_g s}{3 f_g-1}\right) \left(\frac{\eta }{\chi }\right) \frac{1}{x^{5/2}}. \end{eqnarray} \tag{ B6 }$$

We expand τ in a power series about x = 0, τ = x^1/2∑^∞_n=0a_nxⁿ, to obtain

$$\begin{eqnarray} &&-\frac{4 s^3 \eta f_g + 3 a_0 \chi }{4 (3f_g-1) s^2 \chi }x^{-5/2}-\frac{3 a_1}{4(3f_g-1) s^2}x^{-3/2} \nonumber\\ &&- \frac{3 a_2 + 2 s[5 a_0 s^{\prime } + 3 (3f_g-1) a_1]}{4 (3f_g - 1) s^2} x^{-1/2} + O(x^{1/2}) = 0.\nonumber\\ \end{eqnarray} \tag{ B7 }$$

Thus, the leading coefficients near x = 0 are

$$\begin{eqnarray} \ a_ 0 & = & -\frac{4}{3} s^3 \frac{\eta f_g}{\chi } \end{eqnarray} \tag{ B8 }$$

$$\begin{eqnarray} a_1 & = & 0\qquad\quad\ \end{eqnarray} \tag{ B9 }$$

$$\begin{eqnarray} \qquad a_2 & = & \left(\frac{40 f_g s^4 s^{\prime }}{9}\right) \frac{\eta }{\chi }. \end{eqnarray} \tag{ B10 }$$