GRAVITO-TURBULENT DISKS IN THREE DIMENSIONS: TURBULENT VELOCITIES VERSUS DEPTH

The simulations with t_cool = 10Ω⁻¹ (dubbed tc=10 at standard resolution and tc=10.hi at high resolution in Table 1) start from the hydrostatic equilibrium solution described in Section 2.3 and evolve for 30t_cool = 300Ω⁻¹. The evolution of the disk as a whole can be followed by constructing an effective Toomre's Q:

$$\begin{equation} Q\,{\equiv}\,\frac{\langle c_{\rm s} \rangle _\rho \Omega }{\pi G\langle \Sigma \rangle }, \end{equation} \tag{ 16 }$$

where 〈Σ〉 is the vertical surface density averaged horizontally and $\langle c_{\rm s} \rangle _\rho\,{\equiv}\,\langle \sqrt{\gamma P / \rho } \rangle _\rho$ is the density-weighted, box-averaged sound speed.

Figure 2 displays the time evolution of Q. The disk settles into steady gravito-turbulence after an initial transient phase lasting ∼50Ω⁻¹. The transient phase is violent: the cooling disk collapses under its own weight during the first 10Ω⁻¹; rebounds vertically as gas becomes strongly heated by shocks from t = 10Ω⁻¹ to 25Ω⁻¹; and relaxes from t = 25Ω⁻¹ to 50Ω⁻¹ into a quasi-steady state that lasts the remainder of the simulation. During the rebound phase, about 10% of the total mass in the box is lost through the vertical boundaries. This initial mass loss should be of no consequence as we are interested in the final dynamical equilibrium attained by the disk, i.e., the steady self-regulated state in which heating driven by gravitational instability closely balances the imposed cooling.

Zoom In Zoom Out Reset image size

The final equilibrium density profile is broader than our initial profile, as illustrated in Figure 1. In the final state, mass continues to be lost through the top and bottom of the box, but at a rate so slow (≲10% from t = 50–300Ω⁻¹) that we discern no obvious long-term trend in any other statistical property that we measure. This is demonstrated in Figure 2, which attests that in self-regulated gravito-turbulence (for β = 10), Q hovers around 1.33; density fluctuations

$$\begin{equation} \delta \rho = \rho (x,y,z) - \langle \rho (z) \rangle \end{equation} \tag{ 17 }$$

are on the order of unity when normalized to the average local density; velocity fluctuations

$$\begin{equation} \delta v = \sqrt{v_x^2 + \delta v_y^2 + v_z^2} \end{equation} \tag{ 18 }$$

are mildly sonic (here δv_y = v_y + 3Ωx/2 is the non-Keplerian azimuthal velocity); and the Shakura–Sunyaev stress-to-pressure parameter

$$\begin{equation} \alpha\,{\equiv}\,\frac{{\langle w_{xy} \rangle }_\rho }{\langle P \rangle _\rho }\,{\equiv}\,\frac{\langle g_{x} g_{y}/4\pi G+ \rho v_{x} \delta v_{y}\rangle _\rho }{\langle P \rangle _\rho } \end{equation} \tag{ 19 }$$

fluctuates about 0.055, with the gravitational stress (g_xg_y/4πG) exceeding the Reynolds stress (ρv_xδv_y) by roughly a factor of three (here $g_x \mathbf {\hat{x}} + g_y \mathbf {\hat{y}} + g_z \mathbf {\hat{z}}$ is the local gravitational acceleration). Our definition for α weights by density; it tries to "follow the mass" in our stratified simulations and should be more accurate, e.g., when applied to 2D studies that evolve Σ instead of ρ. In the literature—which reports on simulations that are commonly unstratified—the stress-to-pressure ratio is more often calculated as a volume average without weighting by density (e.g., Equation (19) in Gammie 2001):

$$\begin{equation} \alpha ^{\prime }\,{\equiv}\,|\frac{d\ln {\Omega }}{d\ln {R}}|^{-1}\!\! \frac{\int \! w_{xy} dxdydz}{\int \!\rho c_{\rm s}^2 dxdydz} = \frac{2}{3\gamma }\, \frac{\langle w_{xy}\rangle }{\langle P \rangle }, \end{equation} \tag{ 20 }$$

where R is disk radius. This conventional α' is 0.5–0.8 times our α in constant cooling time simulations. We calculate both measures of stress in this paper; see, e.g., Table 1. Our results appear robust insofar as doubling the resolution in every direction changes α by just a few percent (compare red and black curves in Figure 2, and see also Table 1).

All simulations were run for long enough that we can recover steady curves at late times like those shown in Figure 2, enabling us to evaluate meaningful time averages. In only one constant cooling run, tc=3.hi (β = 3), did the disk fragment instead of settling into steady gravito-turbulence (see also our analogous optically thin run b=100.hi for which β < 3, which also fragmented). Thus we establish (for γ = 5/3) that the collapse criterion in 3D local disks is t_cool ≲ 3Ω⁻¹. This is consistent with the fragmentation boundary reported previously in 2D local (e.g., Gammie 2001) and 3D global studies (e.g., Clarke et al. 2007; but see Meru & Bate 2011, Meru & Bate 2012, and Rice et al. 2014 for cautionary remarks).

Figures 3 and 4 show how the density-weighted 〈α〉_t and the volume-averaged 〈α'〉_t—and their gravitational and hydrodynamic components—vary with t_cool for our high-resolution runs (time averages are taken over intervals lasting at least 100Ω⁻¹; see Table 1). Gravitational stresses are observed to exceed hydrodynamic stresses by factors of ∼2–5. We find that α and especially α' scale closely with 1/(Ωt_cool), as expected from energy conservation (Gammie 2001). Figure 4 for α' reveals that our simulation results match the analytic expectation

$$\begin{equation} \alpha ^{\prime } = \frac{4}{9\gamma (\gamma -1)} \frac{1}{\Omega t_{\rm {cool}}} \end{equation} \tag{ 21 }$$

nearly perfectly for our chosen γ = 5/3, indicating excellent energy conservation in our code.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A similar scaling of stress with t_cool applies at all heights, as demonstrated in Figure 5. Note how the local stress—which is dominated by gravitational forces—diminishes by factors of ∼6–16 from |z| = H to |z| = 4H, or roughly as |z|⁻ⁿ with n = 1.2–2, by contrast to the density ρ which drops exponentially with height (Figure 1). The comparatively slow drop-off of stress with height arises because gravity is a long-range force: the stress at altitude arises from gravitational forces exerted by density fluctuations near the midplane. We illustrate this in Figure 6 by plotting the contribution to the total gravitational stress exerted at z = 4H from all heights. About 80% of the stress at z = 4H comes from the gravity of the disk between −H < z < H.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Total α values generally change by 10% or less when the resolution is doubled in all directions; when the imposed density floor is lowered from 10⁻⁴ρ₀ to 10⁻⁵ρ₀; when the box size is doubled horizontally; or when the box height varies from 12H to 10H or 14H (Table 1). Increasing the spatial resolution tends to increase the Reynolds stress and lower the gravitational stress. Velocity fluctuations are better resolved with finer grids, but the gains are modest; doubling the resolution in all directions changes non-circular velocities by ∼30% at most and typically by only several percent; at the same time, the increased Reynolds stress should be offset by a decrease in gravitational stress, since the total heating rate must balance the imposed cooling rate (see Figure 4). Note that simulations with larger β (longer cooling times) have smaller amplitude fluctuations and are more computationally demanding—this may be the reason why, e.g., our β = 80 runs exhibit more sensitivity to resolution than our β < 80 runs (see Table 1). Unless indicated otherwise, the data plotted in all our figures and discussed in the text are drawn from our high-resolution simulations with a standard density floor of 10⁻⁴ρ₀ (i.e., the .hi models in Table 1).

The variation with height of the rms value of the turbulent velocity δv (i.e., the non-Keplerian bulk motion) is shown in Figure 7 for β = 10. Turbulent velocities vary only weakly with altitude, increasing by less than a factor of two from midplane to surface as the overlying column density drops by more than three orders of magnitude. The drop of turbulent velocities from |z| ≈ 5H to 6H in our standard box is a numerical artifact, as we have discovered by varying the box height to 10H (run tc=10.hi.lz10) and 14H (tc=10.hi.lz14)—see Figure 8. The roll-over might be due to the enforced outflow boundary conditions that artificially eliminate any inward motion. However, our results away from the vertical boundaries appear to have converged with box size.

That turbulent motions are just as vigorous at altitude as they are at depth might at first glance seem surprising, since densities at altitude fall below the Toomre density Ω²/2πG ∼ 0.5 (Shi & Chiang 2013). In other words, gas at altitude is not self-gravitating and does not generate turbulence in and of itself. However, in hindsight, the relatively flat velocity profile is understandable, because as we noted in Section 3.1.2, gravity is long-range: density fluctuations at the midplane—where gas is of Toomre density and is self-gravitating—exert gravitational forces at altitude and strongly perturb the rarefied gas there.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In fact, non-Keplerian motions are, if anything, greater at altitude than at depth. Vertical velocities v_z vary more strongly with height (factor of five increase from midplane to surface) than do in-plane velocities v_p, which is sensible because the vertical gravitational acceleration from disk self-gravity tends to increase away from the midplane: at the midplane, vertical forces from the upper and lower disk tend to cancel. In addition, fluid motions in all directions should amplify at altitude from the steepening of waves launched upward from the midplane and which propagate down density gradients. Simon et al. (2011) suggested from their numerical experiments that increasing velocity dispersions with height are generic to turbulent disks, and our results are consistent with this proposal. At |z| > H, motions are strongest in the x–z plane, driven by self-gravity in those directions; azimuthal perturbations are suppressed by Keplerian differential rotation which shears overdensities into nearly axisymmetric structures. We see evidence for a kind of x–z circulation, whereby gas compresses radially from self-gravity and spurts out vertically, as illustrated in Figure 9.

Zoom In Zoom Out Reset image size

We also calculate velocity dispersions more closely related to observables. We forego a full spectral line analysis and instead compute a line-of-sight (los) turbulent velocity in the disk plane, averaged over azimuth ϕ:

$$\begin{equation} |v_{\rm p}| \equiv \frac{1}{2\pi }\int ^{2\pi }_{0} d\phi |v_{\rm x} \cos \phi - \delta v_{\rm y} \sin \phi| \end{equation} \tag{ 22 }$$

(Simon et al. 2011; Forgan et al. 2012). Here we use δv_y for the turbulent velocity in the azimuthal direction, not to be confused with the bulk velocity v_y which includes the background shear velocity. The quantity |v_p| is a proxy for the actual los turbulent linewidth for disks seen edge-on. For disks seen face-on, the corresponding proxy velocity is |v_z|. Figure 7 plots the vertical profiles for 〈|v_p(z)|〉 and 〈|v_z(z)|〉 (each horizontally averaged), while the top left panel of Figure 10 displays the full probability density functions for |v_p|/c_s and |v_z|/c_s for several cuts in height, all for β = 10. There is little variation in these proxy los velocities with height, particularly for |v_p|. The most probable value of |v_p|/c_s increases by a factor of only 1.2 from z > 0 to z > 4H (for β = 10), while |v_z|/c_s increases by a factor of 2.5. We obtain similarly flat velocity profiles for all β > 3 runs; the only difference is a systematic shift toward smaller velocities as t_cool increases, e.g., for β = 10(40), the most probable value of |v_p|/c_s is 0.35(0.22), sampled over all z > 0; see Figure 10).⁵

Zoom In Zoom Out Reset image size

Note that the variations in Mach numbers between simulations with different cooling times are almost all due to variations in turbulent velocities, not to variations in the sound speed. Table 1 indicates that the sound speeds of various simulations are all about the same—as expected, since the input surface densities and self-regulated Toomre Q-values are all about the same, irrespective of cooling time. Prolonging the cooling time reduces turbulent velocities but does not alter background temperatures (at fixed surface density).

Figure 10 also makes a head-to-head comparison between the velocity distributions of our gravito-turbulent disks and those of a disk made turbulent by the MRI. Our MRI data are extracted from the standard run (STD32) of Shi et al. (2010), who studied the MRI using a stratified, non-self-gravitating, ideal MHD shearing box with zero net magnetic flux (for details, see their Sections 2 and 3.1). Turbulent velocities vary more strongly with height in MRI-turbulent disks than in gravito-turbulent disks; as z increases from the midplane and the overlying column density drops by three orders of magnitude, both |v_p|/c_s and |v_z|/c_s increase in the MRI-turbulent disk by factors of ∼15, from ∼0.06 to ∼1. This variation characterizes ideal MHD disks. Similar results, including supersonic velocities at z > 3H, were obtained by Simon et al. (2011). In non-ideal disks (those affected by Ohmic dissipation, ambipolar diffusion, and/or the Hall effect), the variation of magnetically driven velocities with altitude tends to be even steeper (Simon et al. 2011, 2013a; see also Figure 7 of Lesur et al. 2014).

Forgan et al. (2012) suggested that the ratio |v_p|/|v_z| could be used to differentiate between gravito-turbulent and MRI-turbulent disks. They found that |v_p| ∼ 6|v_z| near the midplanes of their gravito-turbulent disks, whereas MHD turbulence is more isotropic (even at the midplane of a dead zone, |v_p| differs from |v_z| by a factor ≲ 3 according to the top panel of Figure 4 of Simon et al. 2011). Our results are nominally of higher resolution, and show that |v_p| is indeed higher than |v_z|, but only by factors of 2–3 at the midplane (Figure 10). Unfortunately, even this difference appears to vanish at altitude. Thus the usefulness of |v_p|/|v_z| in distinguishing between gravito-turbulence and MRI turbulence appears limited.