Two saturated states of the vertical shear instability in protoplanetary disks with vertically varying cooling times

Abstract

1 Introduction

The initial stage of planet formation is the formation of kilometer-sized planetesimals from micron-sized dust grains in protoplanetary disks (see Johansen et al. 2014 and Drazkowska et al. 2022 for recent reviews). Micrometer-sized grains initially grow into larger aggregates through coagulation. Sufficiently large aggregates settle onto the disk midplane and form a dense dust layer (Weidenschilling 1980; Nakagawa et al. 1981; Dullemond & Dominik 2005; Tanaka et al. 2005). Large aggregates can also concentrate radially toward local maxima of the disk gas pressure (e.g., Whipple 1972; Kretke & Lin 2007; Pinilla et al. 2012). These dust concentration processes may lead to planetesimal formation through streaming instability (Youdin & Goodman 2005; Johansen & Youdin 2007; Johansen et al. 2009; Carrera et al. 2015; Yang et al. 2017) and gravitational instabilities (Goldreich & Ward 1973; Youdin 2011; Takahashi & Inutsuka 2014; Tominaga et al. 2018, 2019, 2020; Pierens 2021). If the aggregates are sticky enough, they may form planetesimals directly by successive coagulation (e.g., Okuzumi et al. 2012; Windmark et al. 2012; Kataoka et al. 2013).

In these planetesimal formation processes, gas disk turbulence can play many important roles, including both positive and negative ones. Large long-lived vortices created by turbulence may lead to dust concentration and subsequent planetesimal formation through gravitational collapse (e.g., Barge & Sommeria 1995; Raettig et al. 2021; Lehmann & Lin 2022). On the other hand, strong turbulence induces a large impact velocity between dust particles (e.g., Ormel & Cuzzi 2007) that may inhibit planetesimal formation through coagulation (e.g., Brauer et al. 2008; Okuzumi & Hirose 2012). Turbulence can also diffuse the midplane dust layer and prevent the onset of gravitational and streaming instabilities (e.g., Dubrulle et al. 1995; Umurhan et al. 2020; Chen & Lin 2020). Therefore, the strength of the turbulence would determine how planetesimals form.

Recent radio interferometric observations with the Atacama Large Millimeter–submillimeter Array (ALMA) and other interferometers have provided important constraints on turbulence intensity in the outer region of protoplanetary disks. Molecular line emission observations point to weak or no turbulence in the upper layers (Hughes et al. 2011; Flaherty et al. 2015, 2017, 2018, 2020) and regions within one scale height above and below the midplane (Guilloteau et al. 2012; Teague et al. 2016, 2018) in the outer disk region (for a recent review, see Pinte et al. 2022). Strong nonthermal gas motion is detected in a disk around DM Tau (Flaherty et al. 2020). Observations of dust rings and gaps in massive and large disks (e.g., ALMA Partnership 2015; Andrews et al. 2018; Long et al. 2018; van der Marel et al. 2019) can also be used to indirectly constrain the turbulence intensity because the morphology of the rings and gaps can be affected by turbulent dust diffusion. The well-defined morphology of the dust gaps in the disks around HL Tau (Pinte et al. 2016) and Oph 163131 (Villenave et al. 2022) suggests a significantly low level of turbulence. On the other hand, one of the two major dust rings in the HD 163296 disk shows a high degree of vertical dust diffusion (Doi & Kataoka 2021), potentially implying strong turbulent diffusion in that location. Taken together, the observations so far suggest that the outer disk regions are mostly laminar but can become strongly turbulent under some circumstances.

To better understand when and where strong turbulence emerges, it is essential to study what mechanisms drive disk turbulence. Regarding the outer disk region, it was previously believed that the magnetorotational instability (MRI; Balbus & Hawley 1991) is the most viable mechanism for driving turbulence (e.g., Sano et al. 2000). However, recent theoretical studies have shown that the viability of the MRI in the outer disk is limited by ambipolar diffusion (Simon et al. 2013a, 2013b; Bai 2015; Béthune et al. 2017; Riols & Lesur 2018; Cui & Bai 2021). Instead, recent studies have pointed out the importance of purely hydrodynamic instabilities (for reviews, see Lyra & Umurhan 2019; Lesur et al. 2022), in particular the vertical shear instability (Urpin & Brandenburg 1998; Arlt & Urpin 2004; Nelson et al. 2013; Lin & Youdin 2015). This is a hydrodyamical instability similar to the Goldreich–Schubert–Fricke instability in differentially rotating stars (Goldreich & Schubert 1967; Fricke 1968). This instability requires a vertical gradient in the gas orbital velocity and rapid cooling of disk gas (Urpin 2003; Nelson et al. 2013; Lin & Youdin 2015; Manger et al. 2021). Turbulence driven by the VSI has a predominant vertical motion (e.g., Nelson et al. 2013; Stoll & Kley 2014) that can prevent dust settling (Stoll & Kley 2016; Flock et al. 2017, 2020). Because the VSI requires rapid gas cooling, it is most active in outer disk regions with low optical depths (Malygin et al. 2017; Pfeil & Klahr 2019; Fukuhara et al. 2021). In the outer regions, the VSI can take over the MRI (Cui & Bai 2020, 2022).

Despite its importance, it is yet to be fully understand in what conditions the VSI produces turbulence around the midplane, where dust evolution and planetesimal formation mainly take place. According to local linear analysis, the VSI tends to be suppressed in a midplane region with a high optical depth and hence a low cooling rate (Malygin et al. 2017; Pfeil & Klahr 2019). However, Pfeil and Klahr (2021) recently showed that such an optically thick midplane region can become turbulent if the VSI operates above the region. They performed hydrodynamical simulations using a vertically varying cooling rate profile, and found that a vertical gas motion generated by the VSI crosses the midplane in the final saturated state. Yet, it remains unclear how strongly their results depend on the assumed vertical profile of the gas cooling rate. In general, the disk cooling profile is determined by the size and spatial distribution of the dust grains that dominate the disk opacity (Malygin et al. 2017). Therefore, the cooling profile can change significantly as the grains grow and/or settle toward the midplane (Barranco et al. 2018; Fukuhara et al. 2021). Very recently, Dullemond et al. (2022) also showed that depletion of sub-micrometer-sized dust grains due to coagulation increases the gas cooling time and makes the disk VSI-stable. These aspects are particularly important in the context of planetesimal formation, where dust evolution is inevitable.

In this paper, we investigate systematically how the saturated state of VSI-driven turbulence depends on the vertical profile of the disk cooling rate. We perform global two-dimensional (2D) hydrodynamical simulations of an axisymmetric protoplanetary disk with a parametrized vertical profile of the disk gas cooling timescale. We show that, depending on the thickness of the VSI-stable midplane region, a vertical gas motion generated by the VSI above the midplane either does or does not penetrate into the midplane. We also provide empirical formulas that predict the saturated level of VSI-driven turbulence from a given vertical profile of the gas cooling timescale. Such empirical relations will be useful to study how the saturated state of VSI-driven turbulence varies with the long-term evolution of the background disk.

The paper is organized as follows. In section 2, we describe our numerical hydrodynamic simulation setup, our disk cooling (thermal relaxation) model, and a simulation analysis method. We then present the main results in section 3, and discuss the implications of our study in section 4. Section 5 presents a summary.

2 Numerical method

To solve the hydrodynamical equations [equations (1)–(3)], we use the Athena++ code (Stone et al. 2020) with the combination of the Harten–Lax–Van Leer (HLLC) approximate Riemann solver (Mignone & Bodo 2005), the reconstruction scheme of a second-order piecewise linear method (van Leer 1974), and second-order Runge–Kutta time integrator. The Courant–Friedrichs–Lewy (CFL) number is set to 0.3.

2.1 Simulation setup

At all computational boundaries, we fix both density and pressure to the initial values. For the velocity components normal to the boundaries, we apply the outflow boundary conditions preventing inflow at the inner and outer radial boundaries and reflecting boundary conditions at the upper and lower meridional boundaries.

Grid cells have logarithmically and linearly uniform spacings in the radial and meridional directions, respectively. The radial and meridional domains cover 0.5 ≤ r ≤ 2.5 and π/2 − 5H₀/R₀ ≤ θ ≤ π/2 + 5H₀/R₀. Because a resolution of 100 cells or more per scale height is necessary to resolve VSI-driven turbulence (Flores-Rivera et al. 2020), we adopt a resolution of 128 cells per scale height in both radial and vertical directions. Therefore, our simulations use 4160 × 1280 grid cells.¹

We adopt the code units M = G = R₀ = 1. In this unit system, the orbital period at R = R₀ is P_in = 2π.

2.2 Cooling model

In this study we treat a and b as free parameters. In reality, these values can vary with dust growth and settling. The value of a can increase or decrease depending on dust growth and settling. This is because the optically thick area around the midplane widens for a certain amount of dust concentration due to the presence of regions that change from a high to low optical depth, and narrows beyond that point. On the other hand, dust growth and settling cause dust depletion and increase cooling time at high altitudes (Barranco et al. 2018; Fukuhara et al. 2021), which corresponds to a decrease in b.

We perform 46 simulations with different values of a, β₁, and b, with β₀ = 2 × 10⁻³, as summarized in table 2 of appendix  1. We also perform one simulation with the locally isothermal equation of state. We take β to be constant in time. Figure 1 illustrates the vertical profile of β for four runs. The β profiles for all our runs are shown in figure 11 in appendix  1.

Fig. 1.

Vertical β profile for runs of (a, β₁) = (2.10, 0) and (0.75, 0) (upper panels) and (a, β₁, b) = (2.1, 0.5, 0.6) and (2.1, 0.5, 0.9) (lower panels) at R = 1.0. The gray regions are marked as “unstable layer” and other white regions are marked as “stable layer,” which are determined by applying the linear criterion [equation (13)]. The dotted lines represent β = β_gc. We refer to the thicknesses of the unstable layer (all panels) and stable midplane layer (only lower panels) as ΔL_u and ΔL_s, respectively.

Open in new tabDownload slide

2.3 Linearly unstable and stable layers

All our simulations have linearly unstable layers at some heights (see figure 1 for illustrative examples). Runs with β₁ = 0 have a cooling time monotonically decreasing toward the midplane, yielding a single unstable layer at |z| < z_u, where z_u is the height of the unstable layer’s upper boundary. Runs with β₁ > 0 have two linearly unstable layers sandwiching a midplane region where equation (13) breaks down. We call this midplane region the linearly stable layer and denote the height of its boundary by z_s. When the linearly stable layer is absent, we set z_s = 0.

We define the thicknesses of the linearly stable and unstable layers as ΔL_s = 2z_s and ΔL_u = 2z_u − ΔL_s. When z_s > 0, ΔL_u accounts for the thicknesses of the two separated unstable layers lying at z < 0 and z > 0 (see figure 1). In general, ΔL_u decreases with decreasing a, and ΔL_s increases with increasing b. The values of ΔL_u and ΔL_s for all runs are summarized in table 2 of appendix  1. As we show in section 3, ΔL_u and ΔL_s are key quantities that dictate the saturated state of VSI-driven turbulence.

2.4 Turbulence diagnostics

We quantify the strength of VSI-driven turbulence using the time average of the squared vertical velocity |$\langle v_z^2\rangle$|⁠, where v_z = v_rcos θ − v_θsin θ is the vertical velocity. The bracket 〈·〉 denotes the time average. In our simulations, the time averaging is performed after the system relaxes into a quasi-steady state. The mean squared vertical velocity is related to the vertical diffusion coefficient for gas and small dust particles (Fromang & Papaloizou 2006). We discuss this in more detail in subsection 4.1.

Turbulence also transports the disk’s radial angular momentum. The efficiency of angular momentum transport is measured by the Reynolds stress 〈ρ_gδv_rδv_φ〉, where δv_r = v_r − 〈v_r〉 and δv_φ = v_φ − 〈v_φ〉 are the dispersions of the radial and azimuthal velocities, respectively.

The time that it takes for the system to reach a quasi-steady state differs from one run to another (see figure 3), and may depend on the unstable layer thickness. Therefore, we stop a run at 400 orbits if the quasi-steady state has already been reached by 250 orbits; otherwise, we continue the run until 800 orbits. The stopping times for all runs are summarized in table 2 of appendix  1. For simulations with shorter and longer runtimes, time averaging is performed over periods of 250–400 and 650–800 orbits, respectively.

Fig. 2.

Vertical velocity v_z/c_s for runs with (a, β₁, b) = (2.1, 0.5, 0.6) (left panel) and (a, β₁, b) = (2.1, 0.5, 0.9) (right panel), which relaxes into the T state and pT state, respectively, as a function of R and z at the end of the simulations. The dashed lines indicate the unstable layer’s upper boundaries at z = ±z_u (uppermost and lowermost lines), and the stable midplane layer’s upper boundaries at z = ±z_s (two middle lines).

Open in new tabDownload slide

Fig. 3.

Vertical velocity v_z/c_s for runs with (a, β₁, b) = (2.1, 0.5, 0.6) (T state; left panel) and (a, β₁, b) = (2.1, 0.5, 0.9) (pT state; right panel) as a function of time and z/H_g at R = 1.0. The dashed lines indicate z = ±z_u (top and bottom lines) and z = ±z_s (two middle lines).

Open in new tabDownload slide

3 Results

In this section we present our simulation results to study how the vertical profile of the cooling time affects VSI-driven turbulence in protoplanetary disks. We begin by defining two saturated states of VSI-driven turbulence and then analyze how the saturated state depends on the thicknesses of the unstable and stable layers in subsection 3.1. We construct empirical formulas of |$\langle v_z^2\rangle |_{\rm mid}$|⁠, |$\langle \overline{v_z^2}\rangle$|⁠, and |$\overline{\alpha _{r\phi }}$| as a function of the thicknesses of the unstable and midplane stable layers, i.e., ΔL_u and ΔL_s, in subsection 3.2.

3.1 Two saturated states of VSI-driven turbulence

We find that the width of the linearly VSI-stable layer at the midplane determines the vertical structure of VSI-driven turbulence in a steady state. Specifically, we identify two possible saturated states of turbulence. In the first class of saturated states, which we call the T (turbulent) states, the vertical gas motion generated in the linearly VSI-unstable layers penetrates into the linearly VSI-stable midplane. This state was already seen in the simulations by Pfeil and Klahr (2021). In the second class, which we call the pT (partially turbulent) states, the vertical gas motion is well confined in the VSI-unstable layers, leaving the VSI-stable midplane layer only weakly turbulent.

Figure 2 illustrates the two saturated states. Here, we present two-dimensional maps of the vertical velocity v_z, normalized by the sound speed c_s, at the end of the simulations for (a, β₁, b) = (2.1, 0.5, 0.6) and (2.1, 0.5, 0.9). The run with (a, β₁, b) = (2.1, 0.5, 0.6) relaxes into the T state with a vertically uniform gas motion through |z| < z_u = 3.9H_g. This state is not expected from the linear stability analysis, which predicts that the midplane region of |z| < z_s = 0.9H_g is linearly VSI stable. The run with (a, β₁, b) = (2.1, 0.5, 0.9) relaxes into the pT state with little vertical gas motion inside the linearly VSI-stable layers (z_s = 1.3H_g for this case).

Our simulations do not clearly reproduce a finer spatial profile of vertical velocity caused by the secondary parametric instability, which can be related to the nonlinear saturation process of the VSI. This is because our simulation resolution is insufficient to resolve the parametric instability, which may require ∼300 cells per gas scale height in the radial direction (Cui & Latter 2022), while it is sufficient to resolve the dominant VSI modes (Flores-Rivera et al. 2020).

To see how the final saturated states are reached, we plot in figure 3 the vertical profiles of v_z at R = 1.0 for the two cases displayed in figure 2. In both cases, the vertical gas motion starts to develop near the unstable layer’s upper boundary at |z| = z_u, where the vertical shear is strong. For (a, β₁, b) = (2.1, 0.5, 0.6) (T state), the vertical flows that develop in the VSI-unstable layers overshoot the stable midplane layer and eventually form a unified flow (see also Pfeil & Klahr 2021). For (a, β₁, b) = (2.1, 0.5, 0.9) (pT state), one can see that the gas vertical motion is well confined in the unstable layers. The vertical profile of |$\langle v_z^2\rangle$| is steady over t = 400–800 orbits, suggesting that our simulation captures the final saturated state. Figure 12 in appendix  2 presents the time evolution of turbulence diagnostics for the two runs presented here.

We define the T and pT states more quantitatively by using the time-averaged |$v_z^2$| at the midplane, |$\langle v_z^2 \rangle |_{\rm mid}$|⁠, and |$v_z^2$| averaged both in time and in the vertical direction, |$\langle \overline{v_z^2}\rangle$|⁠. The ratio between the two quantities reflects how strongly the vertical gas flow penetrates into the midplane region. Simulations exhibiting strongly overshooting vertical gas flows (T states) yield |$\langle v_z^2 \rangle |_{\rm mid} \approx \langle \overline{v_z^2}\rangle$|⁠, whereas those with a less turbulent midplane region than in the unstable layers (pT states) yield |$\langle v_z^2 \rangle |_{\rm mid} \ll \langle \overline{v_z^2}\rangle$|⁠. Figure 4 shows |$\langle v_z^2 \rangle |_{\rm mid}$| versus |$\langle \overline{v_z^2} \rangle$| for all the runs presented in this study. The majority of our simulations result in either |$\langle v_z^2 \rangle |_{\rm mid} \approx \langle \overline{v_z^2}\rangle$| or |$\langle v_z^2 \rangle |_{\rm mid} \approx 0.01$|–|$0.1\langle \overline{v_z^2}\rangle$|⁠. In the following, we refer to the T and pT states as saturated states with |$\langle v_z^2 \rangle |_{\rm mid} > 0.1\langle \overline{v_z^2}\rangle$| and |$\langle v_z^2 \rangle |_{\rm mid} < 0.1\langle \overline{v_z^2}\rangle$|⁠, respectively. The states of all our simulation runs are summarized in table 2 of appendix  1.

Fig. 4.

Time and vertical mean squared vertical velocity, |$\langle \overline{v_z^2}\rangle$|⁠, vs. time mean squared vertical velocity at the midplane, |$\langle v_z^2 \rangle |_{\rm mid}$|⁠, at R = 1.0 for all the runs presented in this study. The symbols correspond to the T states (circles) and pT states (crosses). The dashed lines show |$\langle v_z^2 \rangle |_{\rm mid} = C_1 \times \langle \overline{v_z^2}\rangle$| with C₁ = 1, 0.1, and 0.01.

Open in new tabDownload slide

From the examples shown in figures 2 and 3, one can expect that the thickness of the linearly VSI-stable layer, ΔL_s, determines the final saturated state. To test this hypothesis, we map in figure 5 the values of the turbulence diagnostics |$\langle v_z^2 \rangle |_{\rm mid}$| and |$\langle \overline{v_z^2}\rangle$| from all our simulations against ΔL_u and ΔL_s. We find that pT states (⁠|$\langle v_z^2 \rangle |_{\rm mid} < 0.1\langle \overline{v_z^2}\rangle$|⁠) are realized when the linearly stable midplane layer is as wide as ΔL_s ≳ 2H_g.

Fig. 5.

Time mean squared vertical velocity at the midplane |$\langle v_z^2 \rangle |_{\rm mid}$| (upper panel) and its vertical average |$\langle \overline{v_z^2}\rangle$| (lower panel) at R = 1.0 from all the simulation runs, mapped in the ΔL_u–ΔL_s plane. The dashed line indicates the simulations relaxing to pT states (⁠|$\langle v_z^2 \rangle |_{\rm mid} < 0.1\langle \overline{v_z^2}\rangle$|⁠).

Open in new tabDownload slide

Another important finding from figure 5 is that the thickness of the linearly unstable layer, ΔL_u, determines the vertically averaged saturation level |$\langle \overline{v_z^2}\rangle$|⁠. Turbulence is largely suppressed at all heights in the cases of ΔL_u ≲ 2H_g. One can see that |$\langle \overline{v_z^2}\rangle$| as well as |$\langle v_z^2 \rangle |_{\rm mid}$| decrease sharply from |${\sim }2\times 10^{-3}c_0^2$| to |$\ll 10^{-3}c_0^2$| as ΔL_u falls below 2H_g.

For both the T and pT states, the vertical averaged Reynolds stress |$\overline{\alpha _{r\phi }}$| is tightly correlated with |$\langle \overline{v_z^2}\rangle$|⁠. Figure 6 plots |$\langle \overline{v_z^2}\rangle$| versus |$\overline{\alpha _{r\phi }}$| at R = 1.0 for all the runs presented in this study. This figure shows that most simulations produce |$\overline{\alpha _{r\phi }} = 0.01$|–|$0.1\langle \overline{v_z^2}\rangle /c_0^2$|⁠.

Fig. 6.

Time and vertical mean squared vertical velocity, |$\langle \overline{v_z^2}\rangle$|⁠, vs. vertical mean Reynolds stress, |$\overline{\alpha _{r\phi }}$|⁠, at R = 1.0 for all runs. The symbols correspond to the T states (circles) and pT states (crosses). The dashed lines show |$\overline{\alpha _{r\phi }} = C_2\times \langle \overline{v_z^2}\rangle /c_0^{\,2}$| with C₂ = 1, 0.1, and 0.01.

Open in new tabDownload slide

Our simulation results may change if we perform a three-dimensional (3D) simulation. However, the 3D effecs may be minor because Pfeil and Klahr (2021) already showed that the average vertical velocity in a 3D simulation is comparable to that in 2D simulations.

3.2 Empirical formulas for turbulent quantities

The results presented in the previous subsection show that important turbulent quantities |$\langle v_z^2 \rangle |_{\rm mid}$|⁠, |$\langle \overline{v_z^2}\rangle$|⁠, and |$\overline{\alpha _{r\phi }}$| are all tightly correlated with ΔL_u and ΔL_s. This suggests that one can predict these quantities for general cases as a function of ΔL_u and ΔL_s without having to perform further simulations. Here, we construct such formulas based on our simulation results.

The upper panel of figure 7 compares the values of |$\langle v_z^2 \rangle |_{\rm mid}$| from our simulations with those from the empirical formula. This formula reproduces the simulation results of |$\langle v_z^2 \rangle |_{\rm mid}/c_0^2 \approx 10^{-3}$| (T state) and ∼10⁻⁵ (pT state). The sigmoid functions in f_T and f_pT are also useful to represent the sharp decrease in ΔL_u ≲ 2H_g and ΔL_s ≈ 2H_g. Furthermore, this formula also represents the simulation results of Pfeil and Klahr (2021). To quantify the errors between the simulation results and formula, we show in figure 8 the values of |$\langle v_z^2 \rangle |_{\rm mid}$| for some simulations and their corresponding values for the empirical formula as a function of ΔL_s. This figure indicates that the formula is accurate to less than an order of magnitude in |$\langle v_z^2 \rangle |_{\rm mid}$| at ΔL_s ≲ 2H_g. Furthermore, the formula replicates the sharp drop at ΔL_s > 2H_g and |$\langle v_z^2 \rangle |_{\rm mid}/c_0^2 \approx 10^{-5}$| for the pT state’s simulations of |$\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s} = 10.0$|⁠.

Fig. 7.

Upper panel: Comparison of |$\langle v_z^2 \rangle |_{\rm mid}$| from simulations (points) and the empirical formula (equations (16)–(19); background) on the ΔL_u–ΔL_s plane. The triangles show the simulation results of Pfeil and Klahr (2021). The dashed lines are contours of |$\langle v_z^2 \rangle |_{\rm mid} = 10^{-3}$|⁠, 10⁻⁴, 10⁻⁵, and 10⁻⁶ from the formula. The dotted lines show |$3.5\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s}-4.8 = 0$|⁠, |$0.07\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s}+1.8 = 0$|⁠, and |$\ln (\max \lbrace \Delta \tilde{L}_{\rm u}-2.5,\,\,0\rbrace )-\Delta \tilde{L}_{\rm s}+1.8 = 0$|⁠. Lower panel: Same as the upper panel, but comparing |$\langle \overline{v_z^2}\rangle$| from the simulations and from equations (20)–(22).

Open in new tabDownload slide

Fig. 8.

Comparison of |$\langle v_z^2\rangle _{\rm mid}$| from simulations and the empirical formula (equations (16)–(19); dashed lines) as a function of ΔL_s. The circles, crosses, and triangles are from simulations with |$\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s} = 10.0$|⁠, 5.6, and 3.7, respectively.

Open in new tabDownload slide

To verify the accuracy of the formula, we calculate the root-mean-squared error of the formula with respect to |$\log _{10}\langle v_z^2 \rangle |_{\rm mid}$|⁠. For simulations with |$\langle v_z^2 \rangle |_{\rm mid} > 10^{-5}$| the error is 0.48 dex, which means that the formula has an accuracy of less than an order of magnitude. The error increases to 0.95 dex if we include all simulations.

Fig. 9.

Ratio of |$\overline{\alpha _{r\phi }}$| and |$\langle \overline{v_z^2}\rangle /c_0^2$|⁠, C₂, as a function of ΔL_u. The points and crosses plot the simulation results for all the runs presented in this study and for Pfeil and Klahr (2021), respectively. The dashed line shows the fitting function in equation (24).

Open in new tabDownload slide

Fig. 10.

Comparison of |$\overline{\alpha _{r\phi }}$| from simulation results (circles, crosses, and triangles) and empirical formula (equation (23); dashed lines) as a function of ΔL_s. The circles, crosses, and triangles are from simulations with |$\Delta \tilde{L}_{\rm u}-\Delta \tilde{L}_{\rm s} = 10.0$|⁠, 5.6, and 3.7, respectively.

Open in new tabDownload slide

4 Discussion

4.1 Estimating the dust vertical diffusion coefficient

The predominantly vertical gas motion in VSI-driven turbulence causes strong vertical dust diffusion (Flock et al. 2017, 2020). Our simulations show that the gas velocity dispersion at the midplane varies with the thicknesses of the VSI-stable and -unstable layers. Qualifying how the dust diffusion coefficient varies will be useful for studying planetesimal formation (e.g., Johansen et al. 2009) and testing theory with millimeter observations of dust rings and gaps (e.g., Pinte et al. 2016). A direct measurement of the dust vertical diffusion coefficient requires a calculation of dust grain motion in hydrodynamical simulations, which is beyond the scope of this study. Here, we indrectly estimate the dust vertical diffusion coefficient using the gas vertical velocity dispersion measured in our simulations.

We are particularly interested in the value of α_z at the midplane, |$\alpha _{z,\rm mid}$|⁠, because the thickness of the dust sedimentary layer at the midplane scales as |$\alpha _{z,\rm mid}^{-1/2}$| (Dubrulle et al. 1995; Youdin & Lithwick 2007). Our simulations show |$\langle v_z^2 \rangle |_{\rm mid}/c_{\rm s}^2 \approx 2\times 10^{-3}$| (subsection 3.1) for fully developed VSI-driven turbulence in T states with ΔL_u ≳ 2H_g. For this case, equation (27) predicts a dimensionless vertical diffusion coefficient of |$\alpha _{z, \rm mid} \approx 4 \times 10^{-4\cdots -2}$|⁠, with the uncertainty originating from that of τ_corr. This predicted value is higher than the dimensionless Reynolds stress in VSI-driven turbulence (α_rφ ≈ 2 × 10⁻⁴), reflecting the predominantly vertical motion of VSI-driven turbulence. For the cases of ΔL_u < 2H_g and ΔL_s > 2H_g, equation (27) predicts much smaller diffusion coefficients of |$\alpha _{z, \rm mid} \ll 10^{-4}$| and |$\alpha _{z, \rm mid} \approx 4 \times 10^{-6\cdots -4}$|⁠, respectively. The implications of the suppressed VSI-driven turbulence for dust evolution and disk observations are discussed in the following subsection.

4.2 Implications for dust evolution and observations of protoplanetary disks

The suppression of VSI-driven turbulence at the midplane in the cases of ΔL_u < 2H_g or ΔL_s > 2H_g has important implications for dust settling and planetesimal formation in the outer regions of protoplanetary disks. Weak turbulence yields low relative velocities of dust particles, which is preferred for dust growth through coagulation without collisional fragmentation and erosion (e.g., Brauer et al. 2008; Okuzumi & Hirose 2012). Furthermore, the weak turbulent diffusion (α_z ≲ 10⁻⁴) in suppressed turbulence would promote planetesimal formation through the streaming and gravitational instabilities (e.g., Sekiya 1998; Youdin & Shu 2002; Johansen et al. 2009; Gole et al. 2020; Umurhan et al. 2020; Chen & Lin 2020). Therefore, disk regions with ΔL_u < 2H_g or ΔL_s > 2H_g would be preferential sites for planetesimal formation.

The suppression of VSI-driven turbulence may also explain the low level of vertical dust diffusion inferred from ALMA observations of some protoplanetary disks. The dust rings around HL Tau and Oph 163131 exhibit well-separated morphology in millimeter images, indicating that the large dust particles in the rings have settled onto the midplane. Assuming that millimeter-sized particles dominate the millimeter emission, these observations point to small vertical diffusion coefficients of α_z around a few 10⁻⁴ for HL Tau (Pinte et al. 2016) and of α ≲ 10⁻⁵ for Oph 163131 (Villenave et al. 2022). These estimated values of α_z are consistent with suppressed VSI-driven turbulence in the cases of ΔL_u < 2H_g or ΔL_s > 2H_g (see subsection 4.1). Therefore, we hypothesize that VSI-driven turbulence is indeed suppressed in the outer regions of these disks. Testing this hypothesis requires detailed modeling of these disks’ cooling structure.

4.3 Need for self-consistent modeling of dust and VSI evolution

Because dust particles control disk cooling, whether the condition ΔL_u < 2H_g or ΔL_s > 2H_g for the suppression of VSI-driven turbulence is realized would depend on the size and amount of the particles. Dust growth and settling lead to a VSI-unstable region that has smaller ΔL_u, i.e., that is more confined around the midplane (Fukuhara et al. 2021). Depletion of small grains that dominate the gas disk cooling can also increase the cooling time (Dullemond et al. 2022) and thereby make the VSI-unstable region smaller. These effects may result in suppression of VSI-driven turbulence. On the other hand, an increase in the dust surface density leads to a wider optically thick region around the midplane. This may make the cooling in the region around the midplane less efficient and consequently suppress VSI-driven turbulence at the midplane. Assessing whether ΔL_u < 2H_g or ΔL_s > 2H_g can be realized under realistic conditions requires self-consistent modeling of dust growth, dust surface density evolution, and disk cooling that takes into account gas–dust thermal coupling.

Moreover, because a change in the saturated state of VSI-driven turbulence would also affect dust evolution, the evolution of dust and VSI should be coupled. For instance, suppression of VSI-driven turbulence by dust growth, if it really occurs, would reduce the collision velocity between the dust particles and thus promote their further growth. Similar positive feedback can also be expected for suppression of VSI-driven turbulence by dust settling. These positive feedback effects can be important for understanding planet formation and turbulence in outer disk regions.

To quantify these effects, the empirical formulas presented in subsection 3.2 will be useful. These formulas represent the correlation between the thicknesses of the VSI-unstable and -stable layers and the VSI-driven turbulence intensity. The size of the dust particles as well as their spatial distribution controls gas cooling and determines the thicknesses of the VSI-unstable and -stable layers. Our formulas can be used to predict how the saturated level of VSI-driven turbulence would evolve with the long-term evolution of dust. This will be considered in future work.

4.4 Limitations of our simulations

Our simulations are subject to two important limitations that should be addressed in future work. First, the cooling rates adopted in our simulations are vertically varying but constant in time. In reality, because VSI-driven turbulence alters the dust distribution (e.g., Flock et al. 2020), the cooling time distribution determined by dust can change with turbulence. Therefore, turbulence and cooling time would evolve simultaneously until the turbulence and dust profile are in equilibrium states. If the dust evolution timescale is longer than the timescale for VSI-driven turbulence saturation, we can investigate the co-evolution of dust and turbulence using the empirical formulas presented in subsection 3.2. This is because it is reasonable to assume that the size and spatial distribution of the dust grains will not change while the VSI develops turbulence. On the other hand, if this assumption breaks down, the stability of this system should be studied in hydrodynamical simulations that include the thermal coupling between gas and dust. We plan to address these open issues in future work.

Second, our simulations neglect the effects of dust and magnetic fields on gas disk dynamics. Dust would increase the effective buoyancy frequency of the gas (Lin & Youdin 2017) and weaken the VSI (Lin 2019). Magnetic fields threading the disk may also suppress the VSI either directly through magnetic tension or indirectly through MRI turbulence (Nelson et al. 2013; Latter & Papaloizou 2018; Cui & Bai 2020). The roles of magnetic fields in VSI suppression can be positive or negative depending on non-ideal magnetohydrodynamical effects (ambipolar diffusion, Ohmic resistivity, and the Hall effect; Cui & Bai 2020, 2022; Cui & Lin 2021; Latter & Kunz 2022). We plan to quantify these effects using simulations including dust feedback and magnetic fields in the future.

5 Summary

We have investigated how the saturated state of VSI-driven turbulence depends on the vertical profile of the disk cooling rate. We have performed global two-dimensional hydrodynamical simulations of an axisymmetric protoplanetary disk with vertically varying cooling times. Our key findings are summarized as follows.

• The thickness of the linearly VSI-stable layer at the midplane determines the vertical structure of VSI-driven turbulence in a steady state (figures 2 and 3). We have identified two final saturated states of turbulence. In the first state, the vertical gas motion generated in the linearly VSI-unstable layers penetrates the VSI-stable midplane layer (T states). In the second state, the vertical gas motion is well confined in the unstable layers (pT states), leaving the stable midplane layer only weakly turbulent. Using the time averaged squared vertical velocity at the midplane, |$\langle v_z^2 \rangle |_{\rm mid}$|⁠, and |$v_z^2$| averaged both in time and in the vertical direction, |$\langle \overline{v_z^2}\rangle$|⁠, we refer to the T and pT states with |$\langle v_z^2 \rangle |_{\rm mid} > 0.1\langle \overline{v_z^2}\rangle$| and |$\langle v_z^2 \rangle |_{\rm mid} < 0.1\langle \overline{v_z^2}\rangle$|⁠, respectively (figure 4).

• The pT states are realized when the thickness of the VSI-stable midplane layer ΔL_s is larger than two gas scale heights (figure 5). When the thickness of the VSI-unstable layer ΔL_u is thinner than 2H_g, VSI-driven turbulence is also largely suppressed at all heights. The turbulence diagnostic value |$\langle \overline{v_z^2}\rangle$| as well as |$\langle v_z^2 \rangle |_{\rm mid}$| decreases sharply from |${\sim }2\times 10^{-3}c_0^2$| to |$\ll 10^{-3}c_0^2$|⁠, where c₀ is the sound speed, as ΔL_u falls below 2H_g (figure 5).

• For both the T and pT states, the vertical averaged Reynolds stress |$\overline{\alpha _{r\phi }}$| and |$\langle \overline{v_z^2}\rangle$| are connected by |$\overline{\alpha _{r\phi }} = 0.01$|–|$0.1\langle \overline{v_z^2}\rangle /c_0^2$| (figure 6).

• We propose empirical formulas for the turbulence diagnostics |$\langle v_z^2 \rangle |_{\rm mid}$|⁠, |$\langle \overline{v_z^2}\rangle$|⁠, and |$\overline{\alpha _{r\phi }}$| as a function of ΔL_u and ΔL_s [equations (16), (20), and (23), respectively]. These formulas reproduce the tight correlation of turbulence diagnostics with ΔL_u and ΔL_s (figure 7) and have an accuracy of less than an order of magnitude in strong turbulence (figures 8 and 10). These formulas will be useful for predicting how the states of VSI-driven turbulence vary with the long-term dust evolution.

Our results suggest that the suppression of VSI-driven turbulence at the midplane in the cases of ΔL_u < 2H_g or ΔL_s > 2H_g can lead to strongly vertical settling of dust particles. This effect may promote planetesimal formation through dust coagulation and the gravitational and streaming instabilities in outer disk regions. This effect may also explain the low level of vertical dust diffusion observed in the dust rings of some protoplanetary disks.

For simplicity, the present study has modeled the vertical cooling rate profile with a parametrized analytic function. In reality, gas cooling is regulated by dust particles, and therefore the disk cooling structure should depend on the particle size and spatial distribution. This implies that the evolution of dust and VSI can be coupled because VSI-driven turbulence can affect dust evolution. For instance, dust growth would make the VSI-unstable region confined around the midplane and consequently suppress VSI-driven turbulence. This turbulence suppression by dust growth may reduce the collision velocity between dust grains and thus promote their further growth. To quantify this positive feedback effect, the empirical formulas presented in this study may be useful. Verification of planetesimal formation by this positive feedback requires future study of how the VSI and dust co-evolve.

Acknowledgements

We thank Mario Flock for a discussion about the correlation time of VSI-driven turbulence. We also thank the anonymous referee for helpful comments. This work was supported by JSPS KAKENHI Grant Numbers JP18H05438, JP20H00182, JP20H01948, JP20J01376, and JP22J22593. Numerical computations were carried out on a Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.

Appendix 1 Parameter choices and final saturated states for all runs

We perform 46 hydrodynamical simulations with different values of a, β₁, and b that determine the vertical profile of the cooling time β [equation (12)] and thicknesses of the unstable and midplane stable layers [equation (13)]. The values of a, β₁, b, ΔL_u, and ΔL_s for all the runs presented in this study are summarized in table 2. Figure 11 illustrates the vertical β profiles for all runs.

We stop a run at a different time for each run, because the time until the system reaches a quasi-steady state differs from one run to another. The runtimes and final saturated states of all our simulations are also summarized in table 2.

Fig. 11.

Vertical profiles of the dimensionless cooling time β as a function og z/H_g at R = 1.0 for all the runs presented in this study.

Open in new tabDownload slide

Appendix 2 Examples of turbulence time evolution and vertical structure for runs with T and pT states

In figure 3, we display the two-dimensional maps and vertical profile’s time evolution of the vertical velocity for runs with (a, β₁, b) = (2.1, 0.5, 0.6) and (2.1, 0.5, 0.9) that reach the T state and pT state, respectively. These two runs also differ in their time evolution of |$v_z^2|_{\rm mid}$| and |$\overline{v_z^2}$|⁠. Figure 12 shows the time evolution of |$v_z^2|_{\rm mid}$| and |$\overline{v_z^2}$| for the two runs presented in subsection 3.1. For (a, β₁, b) = (2.1, 0.5, 0.6) (T state), the vertical flows develop at ≲ 50P_in, corresponding to the growth rate Γ_VSI being ∼10⁻²Ω_K. On the other hand, for (a, β₁, b) = (2.1, 0.5, 0.9) (pT state), the vertical gas motion grows at ≲ 300P_in with Γ_VSI ≈ 10⁻³Ω_K. The values of |$v_z^2|_{\rm mid}$| and |$\overline{v_z^2}$| are steady over t = 250–400 orbits and t = 650–800 orbits for runs with the T state and pT state, respectively, which are consistent with the vertical velocity’s time evolution shown in figure 3.

Fig. 12.

Time evolution of the squared vertical velocity at the midplane |$v_z^2|_{\rm mid}$| and |$v_z^2$| averaged in the vertical direction |$\overline{v_z^2}$| for runs with (a, β₁, b) = (2.1, 0.5, 0.6) (T state; upper panel) and (a, β₁, b) = (2.1, 0.5, 0.9) (pT state; lower panel) at R = 1.0.

Open in new tabDownload slide

These two runs reach different quasi-steady states, resulting in different turbulence profiles (see figures 2 and 3). To evaluate this quantitatively, we plot in figure 13 vertical profiles of |$\langle v_z^2 \rangle$| and α_rφ for these two cases. In both cases, |$\langle v_z^2 \rangle$| and α_rφ have a vertical profile with peaks near the unstable layer’s upper boundary at z = ±z_u. The run with (a, β₁, b) = (2.1, 0.5, 0.6) (T state) has almost the same turbulent structure as the run of the isothermal state, which is the fundamental ideal state for the VSI. In contrast, for (a, β₁, b) = (2.1, 0.5, 0.9) (pT state), both |$\langle v_z^2 \rangle$| and α_rφ decrease sharply near the stable midplane layer’s upper boundary at z = ±z_s. This vertical distribution of |$\langle v_z^2 \rangle$| yields the difference of more than an order of magnitude between |$\langle v_z^2 \rangle |_{\rm mid}$| and |$\langle \overline{v_z^2} \rangle$| (see figure 4). The increase in |$\langle v_z^2 \rangle$| around the midplane may be related to the linear VSI local criterion in appendix  3. Additionally, figure 13 indicates that α_rφ is an order of magnitude smaller than |$\langle v_z^2 \rangle$| for both the T and pT states. This is consistent with the correlation between |$\overline{\alpha _{r\phi }}$| and |$\langle \overline{v_z^2}\rangle$| shown in figure 6.

Fig. 13.

Vertical profile of time mean squared vertical velocity |$\langle v_z^2 \rangle$| (upper panel) and Reynolds stress α_rφ (lower panel) for runs with the isothermal state, (a, β₁, b) = (2.1, 0.5, 0.6) (T state), and (a, β₁, b) = (2.1, 0.5, 0.9) (pT state) at R = 1.0. The dashed dotted lines represent the height of the unstable layer’s upper boundaries z = ±z_u, which is the same for these two runs. The dotted lines represent the height of the midplane stable layer’s upper boundaries z = ±z_s for each run.

Open in new tabDownload slide

Appendix 3 Local criterion for linear VSI

So far, we have used the global criterion of the linear VSI in equation (13) to determine the thicknesses of the unstable and stable layers. However, the slight rise of |$\langle v_z^2 \rangle$| near the midplane of the run with (a, β₁, b) = (2.1, 0.5, 0.9) (pT state) shown in the upper panel of figure 13 may be due to a slightly thin unstable midplane layer determined by a local criterion of the VSI.

The presence of the thin unstable midplane layer can be confirmed by looking at the vertical dependence of the local criterion. Using the ideal gas law P = k_Bρ_gT/m_g, where k_B is the Boltzmann constant and m_g is the mean molecular mass of the gas, and assuming vertically isothermal and vertical hydrostatic equilibrium [equation (6)], we get |$N_z^2 \propto (\partial _z \ln \rho _{\rm g})^2 \propto z^2$|⁠. With this assumption and ∂_z(RΩ) ∝ z [see equation (10)], the local dimensionless critical cooling time can be analytically performed, resulting in β_lc∝|z|⁻¹. Therefore, in regions where z is sufficiently small (but not zero), β_lc increases, and thus the local criterion in equation (A1) is satisfied for any finite β. As a result, the thin VSI-active layer determined by the local criterion always exists near the midplane,⁴ although the VSI does not operate in the midplane (z = 0) where the vertical shear vanishes.

This thin unstable layer near the midplane may be weakly turbulent by the VSI. This may produce the slight vertical shading visible near the midplane in the right panel of figure 3 and the slight increase of |$\langle v_z^2 \rangle$| near the midplane in the upper panel of figure 13. Therefore, to discuss exactly where the VSI operates turbulence, it may be necessary to apply the local criterion as well as the global criterion.

Footnotes

References