{'fraction_non_alphanumeric': 0.03612570635471399, 'fraction_numerical': 0.018558540695945275, 'mean_word_length': 4.300683131897005, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Dynamics and Evolution of Gravitational Instability-Dominated Disks

We derive the evolution equations describing a thin axisymmetric disk of gas and stars with an arbitrary rotation curve that is kept in a state of marginal gravitational instability and energy equilibrium due to the balance between energy released by accretion and energy lost due to decay of turbulence. Rather than adopt a parameterized alpha prescription, we instead use the condition of marginal gravitational instability to self-consistently determine the position- and time-dependent transport rates. We show that there is a steady-state configuration for disks dominated by gravitational instability, and that this steady state persists even when star formation is taken into account if the accretion rate is sufficiently large. For disks in this state we analytically determine the velocity dispersion, surface density, and rates of mass and angular momentum transport as a function of the gas mass fraction, the rotation curve, and the rate of external accretion onto the disk edge. We show that disks that are initially out of steady state will evolve into it on the viscous timescale of the disk, which is comparable to the orbital period if the accretion rate is high. Finally, we discuss the implications of these results for the structure of disks in a broad range of environments, including high redshift galaxies, the outer gaseous disks of local galaxies, and accretion disks around protostars.


INTRODUCTION
In the past few years, observational and theoretical advances in many areas have led to intense study of the role of accretion and gravitational instability in determining the structure and rates of transport through disks. In the high redshift universe, clumpy star-forming galaxies at redshifts z ∼ 2 − 3 (e.g. Elmegreen et al. 2004Elmegreen et al. , 2005Genzel et al. 2008;Förster Schreiber et al. 2009;Cresci et al. 2009) appear to be undergoing rapid accretion, and also have velocity dispersions that are much larger than those present in local galaxies. Numerical simulations (e.g. Dekel et al. 2009a;Ceverino et al. 2009;Agertz et al. 2009;Bournaud & Elmegreen 2009;McNally et al. 2009) suggest that the large velocity dispersion and the massive clump morphology are both produced by a combination of gravitational instability and rapid external accretion. Around active galactic nuclei, radiative cooling pushes thin accretion disks into a state of gravitational instability (Shlosman et al. 1990;Goodman 2003), and in this state their accretion rates and structures may be determined by gravitationallydriven turbulence (Gammie 2001). Closer to home, accretion disks around protostars of mass 1 M ⊙ are expected to experience strong gravitational instability for a significant part of their lives (Kratter et al. 2008) due to a combination of rapid accretion and strong radiative cooling. Numerical simulations indicate that the noncircular motions produced by this instability provide the dominant mechanism for mass and angular momentum krumholz@ucolick.org burkert@usm.uni-muenchen.de transport in the disk (e.g. Krumholz et al. 2007bKrumholz et al. , 2009a. Dozens of simulations of gravitational instability in disks have been published, both for disks undergoing external accretion (Vorobyov & Basu 2007Kratter et al. 2010;Machida et al. 2010, e.g.) and for those in isolation (e.g. Lodato & Rice 2004Kim & Ostriker 2007;Cai et al. 2008;Cossins et al. 2009;Agertz et al. 2009; see Durisen et al. 2007 for a review of earlier work). Based on these, several authors have presented one-dimensional time-dependent disk evolution models in which the effects of gravitational instability are approximated by an α prescription, with α obtained by fits to simulation results or by general energy arguments (Goodman 2003;Hueso & Guillot 2005;Kratter et al. 2008;Rice & Armitage 2009;Rice et al. 2010). Similar order-of-magnitude energy arguments have been extended to the case of galactic disks by Dekel et al. (2009b), Klessen & Hennebelle (2009), and Elmegreen & Burkert (2010). In the realm of purely analytic work, Bertin & Lodato (1999) present steady-state solutions for self-gravitating disks with decaying turbulence, while Rafikov (2009) and Clarke (2009) derived steady-state accretion rates for disks in balance between radiative cooling and accretion-driven heating in protostellar disks.
While the analytic and one-dimensional models have provided a good understanding of the basic mechanism of gravitationally-driven turbulence and transport in disks, they also suffer from significant weaknesses. No analytic models published to date consider the case of gravitational instability-dominated disks that are time-dependent rather than in steady state. Most previous work has been limited to a particular rotation curve (e.g. Keplerian or flat), to a disk of pure gas without stars, and to a disk that is vertically supported by thermal pressure rather than supersonic turbulence. As a result of these limitations it is not clear under what circumstances disks that are not in equilibrium can be expected to evolve into it, and it is not even clear what the equilibrium state is for a star-forming, supersonically turbulent disk such as a galactic disk.
Our goal is to improve this situation by developing a first-principles theory for the evolution of a thin, supersonically turbulent disk of star-forming gas in a state of marginal gravitational stability, driven by a specified rate of external accretion. The only assumptions we make are (1) that the disk maintains a state of marginal gravitational instability (Q = 1) at all times, and (2) the rate of energy loss due to radiative cooling can be parameterized as a certain fraction of the energy per crossing time of a disk scale height. We do not assume that the disk is in steady state or that it is characterized by any particular rotation curve. From these assumptions, in Section 2 we derive equations describing the instantaneous, position-dependent rates of mass, energy, and angular momentum transport, and the time evolution of the disk surface density and velocity dispersion. In Section 3 we show that these equations admit an exact steady-state solution, and we derive the steady-state profiles of surface density, velocity dispersion, and transport of mass, energy, and angular momentum. In Section 4 we show that disks that are not in the steady state will evolve toward it, and that for high accretion rates this evolution occurs on an orbital timescale. Finally, in Section 5 we discuss the implications of our findings, and we summarize in Section 6.

EVOLUTION EQUATIONS
2.1. Mass, Angular Momentum, and Energy Transport 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: Here ρ, v, e, and p, are the gas density, velocity, specific internal energy, and pressure,ρ * 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. 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 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), (2), and (3) imply where is the viscous torque, j = rv φ is the specific angular momentum, Ω = v φ /r is the angular velocity, β = ∂ ln v φ /∂ ln r, and G = Γ dz and L = Λ dz are the vertically-integrated rates of non-fluid energy gain and loss. Equations (4), (5), and (6) are the standard equations of mass, angular momentum, and energy conservation 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 2 φ /r and ∂j/∂t ≈ 0. Using these conditions in equations (4), (5), and (6), we arrive at the evolution equations for Σ and σ: 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 § 5.5.

Star Formation and Radiative Gain and Loss
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Σ * , the rates of radiative gain and loss G and L, and the torque 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 2009;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 ratė 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 midplane density in the disk. Even with this correction, however, we can be confident that ǫ ff 0.1. Moreover, if a significant fraction of the 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Σ * = 0.0067ΣΩ; in comparison, Kennicutt (1998) reports an observed star formation rateΣ * = 0.011ΣΩ, identical to within the errors. 1 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)Σσ 2 . The condition that the disk lose this amount of energy per crossing time of the disk scale height then reduces to 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 G = 0. We return to the question of the real value of G in Section 5.3.

Gravitational Stability and the Torque Equation
We now turn to the central hypothesis of our model, which is that the 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. (2009). For a disk of gas plus stars, the parameter that describes its stability is (Rafikov 2001) and I 0 is the Bessel function of order zero. Here Σ * and σ * are the surface density and velocity dispersion of stars, κ = [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 § 5.2 for a further discussion of this point. For these disks we therefore adopt σ * = σ. In this case the expression is accurate to better than 7%. With this approximation, the condition for a disk to remain marginally stable at Q = 1 becomes Plugging in our expressions for the temporal derivatives of σ, Σ, and Σ * , we obtain an equation that describes the torque required to maintain Q = 1: 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Ṁ =Ṁ ext (so ∂T /∂r = −v φ (β +1)Ṁ 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Ṁ =Ṁ 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 > 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 0 , 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). 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 T = 0 is a solution forṀ 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.

Non-Dimensional Equations
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 T = τṀ ext v φ (R)R, then we can rewrite equation (17) as subject to the boundary condition τ ′ = −β − 1 at x = 1.
Here primes indicate differentiation with respect to x, and we have defined The form of equation (22) is instructive. The coefficients h 0 and h 1 appearing on the left-hand side depend on the current state of the disk without reference to external accretion or turbulent dissipation. Those 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 dis- we have defined Σ * = S * Ṁext /(v φ (R)R) in analogy with S, and the Q = 1 condition in dimensionless form is In these units the orbital period is 1, and the fraction of the disk mass that accretion adds per orbit is S + S * −1 , where the angle brackets indicate an average over the disk area.

The Steady-State Solution
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 ∂Ṁ /∂r = 0, or in dimensionless form d dx 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 those 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 and the steady-state condition τ = −x then implies 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), and the corresponding dimensionless viscosity parameter Shakura & Sunyaev (1973) and viscous evolution timescales are 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 τ = − √ 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 ησ 2 Ω, 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 balances this radiative loss. For a flat rotation curve, the decrease in potential involved in moving from radius r 0 to radius r is v 2 φ ln(r/r 0 ), so inward drift at a velocity v r causes the potential energy per unit mass to decrease at a rate v 2 φ (v r /r). Equating the rates of dissipation and energy increase gives ησ 2 Ω = v 2 φ (v r /r), and it follows immediately that v r = ησ 2 /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 , gravitationallyinduced turbulence produces an effective viscosity α = [γ(γ − 1)(9/4)Ωτ c ] −1 . 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 the 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).

Conditions for Steady-State
In deriving the steady solution we have ignored the change in total disk mass due to accretion. More subtly, our steady solution has constantṀ all the way in to r = 0, so mass effectively vanishes through the origin. These approximations are only reasonable for time scales over which the total disk mass changes little. We can define the accretion timescale for our steady solution as (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 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: (41) For our fiducial ǫ ff = 0.01, this gives 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 Comparing the viscous time to the timescales relevant for accretion and star formation, we have t visc t SF = 162 π 2 1/4 f 7/6 g ǫ ff η 1/3 χ 2/3 = 0.08f 7/6 g χ −2/3 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 −3 , or longer if the disk is non-star-forming (ǫ ff = 0). Thus we expect that our steady state model is reasonable under these conditions.

TIME-DEPENDENT SOLUTIONS
We now consider a disk that is initially out of steady state, with a specified initial value of χ and f g . For the reasons discussed in the previous section, we take ǫ ff = 0 and χ and f g as constant. We start each calculation from an initial velocity dispersion specified on a grid of N x cells, logarithmically spaced. The grid runs from x = x 0 to 1, where x 0 is close to zero. For the boundary conditions, we cannot simultaneously require that τ obey the regularity condition derived in Appendix B at the inner boundary and that τ ′ = −(β + 1)u at the outer boundary; this would amount to applying three boundary conditions to a second-order ODE, and for general choices of s no solution for τ exists that satisfies all three conditions. Instead, we continue to fix τ ′ (and thus the accretion rate) at the outer boundary, and at the inner boundary we require that the torque be τ = −x 0 . This choice amounts to requiring that the work done by gas in the region x < x 0 on the computational domain at x ≥ x 0 approaches 0 as x 0 → 0, while the mass flux through the inner boundary at x = x 0 is allowed to vary freely. This is a good approximation to a disk being fed from the outside a fixed rate and that has an inner bulge region or an inner boundary layer that is stress free, but which can accept mass at varying rates. Note that we do not set τ = 0 at x = x 0 exactly, because this is inconsistent with the steady state solution.
With this setup, we evolve the system according to equation (28) using the algorithm given in Appendix C 2 . Note that we find that the unmodified evolution equation (28) for s is numerically unstable to the growth of small oscillations on the grid scale. We damp these by adding a small amount of viscosity to the disk evolution, implemented in a manner that maintains exact energy conservation.
In Figure 1 we show the evolution of disks with β = 0 (flat rotation curve), f g = 0.5, η = 3/2, x 0 = 0.1 and N x = 500. The left panels show runs with χ = 0.1, and the right panels show runs with χ = 0.01. The top row shows disks with an initial velocity dispersion s 0 equal to either twice the equilibrium value s eq given by equation (34), and the bottom row shows disks with an initial velocity dispersion that is half this value. As the figure shows, all disks evolve toward the equilibrium solution found in Section 3 very rapidly, regardless of whether they start with initial velocity dispersions smaller or larger than s eq . The runs with χ = 0.1 halve their distance from the equilibrium solution in less than a quarter of an outer orbital period, and converge to within 10% of the equilibrium solution within one full outer orbit, after which point they are essentially static. 3 in fact, the time required to reach equilibrium seems to be a factor of ∼ 2 − 3 less than our naive estimate that it should the viscous time t visc . This may be because t visc is an estimate of the time for the material to reach zero radius, while in our case the material need not travel as far Here seq = [χ/(ηfg )] 1/3 / √ 2 (equation 34) is the analytically-computed equilibrium velocity dispersion, and these runs use η = 3/2, fg = 1/2. In each plot, the black curve labeled τ = 0 is the initial velocity dispersion, the red curve is the velocity dispersion at τ = 0.25 (i.e. after 0.25 outer orbits), and each subsequent curve after that represents a factor of 2 increase in time: τ = 0.5 (green), τ = 1 (blue), τ = 2 (purple), and, for the χ = 0.01 runs, τ = 4 (aqua). The thick dashed black line is s = seq.
to set up an equilibrium velocity dispersion and surface density profile.
Thus the time required to reach equilibrium is well below the accretion time t acc = 1.6t orb (for f g = 1/2) we computed in Equation (40), and is far less than the star formation time t SF = 20t orb given by Equation (42), in accord with our analytic expectations. We therefore conclude that disks with χ = 0.1 converge to the equilibrium configuration on a timescale short compared to both the accretion and star formation timescales. For χ = 0.01 the convergence to equilibrium is slightly slower, but the runs are within 10% of equilibrium by 2 outer orbits, and are within ∼ 1% of equilibrium by 4 outer orbits. Since the accretion timescale is 7.6 orbital times for χ = 0.01, and the star formation time is 20t orb , these runs too reach equilibrium fast compared to t acc or t SF .
Thus we have demonstrated that the time-independent solution we obtained in Section 3 is not only an exact equilibrium, it is an attractor toward which initially outof-equilbrium disks will converge. As long as χ 10 −3 (or for arbitrarily small χ in non-star-forming disks), this convergence occurs on a timescale short compared to either the either the star formation timescale or the accretion timescale over which the disk mass changes appreciably. This means that our earlier decision to neglect both star formation and changes in the rotation curve due to accretion is reasonable, and that a disk that is out of steady state will converge to its time-independent equilibrium state faster than either star formation or accretion can alter that equilibrium. It also implies a vast simplification in comparing to observations: since convergence to equilibrium is fast, we can generally assume that the velocity dispersions and surface densities of ob- Fig. 2.-Disk velocity dispersion σ versus redshift z for halos of mass M h,12 = 0.1, 0.3, and 1.0 (blue, red, and black lines) and gas fraction fg = 1/2 or 1 (dashed and solid lines). The plot uses our fiducial η = 3/2 and assumes f b,0.18 = 1/2, i.e. half the infalling baryons are gas and half are stars. The fg = 1 case is appropriate for systems where there no stars or where σ * ≫ σ, such as presentday galactic disks, while the case fg = 1/2 is appropriate for highredshift disks where σ * ≈ σ. served disks reflect the instantaneous equilibrium state dictated by their current external accretion rates and gas fractions. Of course we still have not included star formation feedback, a topic we approach in Section 5.3.

Cosmological Evolution of the Velocity Dispersion of Galactic Disks
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. (2009) extend this to give an estimate of the gas accretion rate onto the disk at the center of the halo, which we adopt: 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 12 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 where f (z) is a function that is linear in time and varies from unity at z = 2.2 to 0.5 at z = 1. 4 InsertingṀ g from equation (46) forṀ ext in 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.
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 −1 . While this is very slightly higher than the value of σ ≃ 8 km s −1 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 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 −1 found in galactic disks found at redshifts ∼ 1.5 − 3 (Cresci et al. 2009). The observed galaxies likely correspond to ∼ 10 12 M ⊙ halos. For redshifts in this range and f g ∼ 1/2, typical of galaxies at that redshift (Daddi et al. 2009;Fig. 3.-Ratio of disk maximum circular velocity Vmax to velocity dispersion σ as a function of redshift z for halos of mass M h,12 = 0.1, 0.3, and 1.0 (blue, red, and black lines, top to bottom) and gas fraction fg = 0.1, 0.5, or 1 (dot-dashed, dashed, and solid lines). All parameters are the same as for Figure 2. Tacconi et al. 2010) we predict typical velocity dispersions of 30 − 50 km s −1 , 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.
Interestingly, we find that there is little dependence of V max /σ on halo mass. Instead, the primary dependence is in 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 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.

High-Redshift Galaxies
It is particularly interesting to apply our models to the z ∼ 2 − 3 galaxies observed by Elmegreen et al. (2004Elmegreen et al. ( , 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 12 M ⊙ , thought to be typical of the observed systems, the models shown in Figures 2 and 3 give χ = 8 × 10 −3 − 1.1 × 10 −2 . 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Ṁ * ≈Ṁ ext . Since the star formation rates are observed (and have typical values ∼ 100 M ⊙ yr −1 ), we can plug them into our model in place ofṀ ext in order to predict disk properties.
Doing so, we find that the redshift 2 − 3 disks should have velocity dispersions (from equation 34) whereṀ * ,100 =Ṁ * /100 M ⊙ yr −1 , independent of their maximum rotation velocities V max . Thus galaxies of 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) where R 10 = R/10 kpc and V 200 = V max /200 km s −1 , and the gas mass is (from equation 35) The ratio of baryonic to dynamical mass within the disk region R is 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 over which the heat is where we have taken the radial scale length to be half of R. In contrast, the time required to double the stellar mass is (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.

Effects of Stellar Feedback
In our idealized models, we have neglected the influence of stellar feedback by setting 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 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. 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 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 L remains larger than 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 G appears on the right hand side of the torque equation with the opposite sign as L, but their functional dependence on other disk parameters (surface density, velocity dispersion, etc.) are 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.

Protostellar Disks
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. (2008Kratter et al. ( , 2010, who introduce the dimensionless numbers: where c s,d is the sound speed in the disk, M * d is the total mass of the disk and star and Ω k,in is the Keplerian orbital period of the infalling material. Physically, ξ represents the ratio of the external accretion rate to the maximum rate (∼ c 3 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 outer disk orbit. Indeed, since v φ (R) = RΩ k,in = GM * 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 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 M ∼ 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Ṁ ext ∼ 10 −5 M ⊙ yr −1 for typical out 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 transsonic 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.

On the Validity of a Local Viscous Approximation
for Gravitational Instability-Induced Transport 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 with 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 2004Boley 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 ( § 5.1) and for non-fragmenting protostellar disks ( § 5.4). It is marginally violated for the observed disks at z ∼ 2 ( § 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).

SUMMARY
In this paper we derive the basic evolution equations for a disk of gas and stars kept in a state of marginal gravitational instability by a combination of external accretion, inward migration of gas, and decay of turbulent motions due to radiative shocks. In such a disk, we use the equations of conservation of mass, angular momentum, and energy to derive an equation (22) that characterizes the instantaneous rates of mass and angular momentum transport required to maintain the state of marginal stability, and we show that this equation has an analytic steady-state solution in which the disk velocity dispersion (Equation 34), surface density (Equation 35), and rates of transport (Equations 36 and 37) through the disk are determined by the rate of external infall onto the disk and the gas mass fraction within it. We show that disks converge to this steady state on timescales of order the orbital time, much less than the time over which either the rotation curve or the gas mass fraction changes significantly.
Based on our analytic solution for the properties of a gravitational instability-dominated disk and their dependence on the gas mass fraction and the infall rate, we are able to gain new insight into several processes. We show that the velocity dispersions of both the outer H i disks of present day galaxies and the main disks of redshift ∼ 2 galaxies can be understood naturally if they are in a state of gravitational instability-regulated equilibrium. Moreover, we can understand the general progression of galactic disks from low values of rotation speed to velocity dispersion ratio, V max /σ, at high redshift to much higher values today. This progression is driven primarily by a falloff in galaxy accretion rates and secondarily by the development of disks with stellar velocity dispersion much lager than the gas velocity dispersion, reducing the importance of stars in setting the gravitational instability condition. We also predict that the observed range of V max /σ values seen at z ∼ 2 is primarily a sequence in gas mass fraction. Finally, we use the same model to study the velocity dispersions of protostellar disks. We show that our results are in good agreement with numerical simulations of gravitational instability in disks, and we predict that velocity dispersions should be transsonic in class 0 and I protostars, dropping to subsonic for class II and III sources.
Although our attention in this paper is focused on cases that can be solved analytically or nearly so, we close by pointing out that our model, as a result of its grounding in the basic equations of fluid dynamics, is also amenable to a more general numerical treatment. One can easily relax our assumptions of constant gas fraction, negligible influence from stellar feedback, and a fixed relationship between gas and star velocity dispersion. The resulting equations are identical to the ones we have already solved, except that they would need to be solved numerically. There is no fundamental barrier to doing so however, and the result would be a new method for simulating the evolution of marginally unstable star-forming disks that is intermediate between purely analytic models such as those we have pursued here and full numerical simulations that can be extremely costly. We plan to pursue this avenue in future work.