Gas and Dust Dynamics in Starlight-heated Protoplanetary Disks

The initial condition is derived by constructing 2D axisymmetric profiles of density, temperature, and rotation velocity, which are in radiation hydrostatic equilibrium. We use star and disk parameters from a model of a typical T Tauri system that has been applied to explain various spatially resolved multiband observations of protostellar disks (Wolf et al. 2003, 2008; Schegerer et al. 2008, 2009; Sauter et al. 2009; Madlener et al. 2012; Liu et al. 2012; Gräfe et al. 2013). In Table 1 we summarize the parameters that define the setup. As presented in Flock et al. (2016), we construct the initial condition by iterating between the radiation transfer and solving for hydrostatic equilibrium. The opacity is determined by the dust-to-gas mass ratio of sub-10 μm grains, which we take to be 10⁻³, the same as in model M3 of Flock et al. (2017b). This choice yields thermal equilibration timescales short enough for the VSI to operate as detailed below in Section 4.2. To obtain the model M3, Flock et al. (2017b) extended the 2D axisymmetric initial conditions to 3D with an azimuthal extent of π/8 radians. Small velocity perturbations were added and the simulation was run for 400 inner orbits.

Table 1.  Setup Parameters for the 3D Radiation HD Disk Model, Including the Gas Surface Density, the Stellar Parameters, the Opacities for Stellar Irradiation and Thermal Emission, the Dust to Gas Mass Ratio of Small Grains ($\leqslant 10\,\mu {\rm{m}}$), the Domain, and the Resolution

Surface Density	${\rm{\Sigma }}=6.0{\left(\tfrac{r}{100\mathrm{au}}\right)}^{-1}{\rm{g}}\,{\mathrm{cm}}^{-2}$
Stellar parameters	T* = 4000 K, ${R}_{* }=2.0\,{R}_{\odot },{M}_{* }=0.5\,{M}_{\odot }$
Opacities	κ* = 1300 cm2 g−1
	κd = 400 cm2 g−1
Dust-to-gas ratio	D2G = 10−3
(< 10 μm)
Grid	r:θ:ϕ
Domain	20 − 100 au:  ± 0.35 rad: 2π rad
Resolution	1024 × 512 × 2044

Download table as:  ASCIITypeset image

To save recomputing the linear growth of the VSI, we here extend and restart the reference model M3 after 400 inner orbits. We first stretch the azimuthal extent by a factor of 2 and then repeat it eight times to fill the new domain's 2π-radian azimuthal extent. The stretching keeps the number of azimuthal cells low enough that the grid fits in the available memory. We next add small cellwise random velocities with amplitude 10⁻⁴c_s to the existing velocity field to break the eightfold symmetry.

The logarithmically spaced radial and uniformly spaced meridional and azimuthal grids yield cells with aspect ratio approximately 1:1:2 and a resolution of around 70 cells per scale height. As model M3 produced features that were mostly axisymmetric, reducing the azimuthal resolution by a factor of 2 is unlikely to influence the results. The simulation with the new domain is then evolved for a further 400 inner orbits. After 250 inner orbits, when the VSI turbulence has reached a new steady state, we distribute uniformly over the midplane between 25 and 95 au a half-million 0.1 mm and a half-million 1 mm particles. Each particle is given the local Keplerian velocity. The Stokes number at the midplane for the 1 mm grains lies between 3 × 10⁻² at the inner boundary and 10⁻¹ at the outer boundary. The radial midplane profile of the Stokes number is shown in Appendix A.

After the full 2π radian simulation restarts, the VSI quickly reaches a new steady state. We gauge the turbulence level by the stress-to-pressure ratio α. We calculate the pressure-weighted⁶ stress-to-pressure ratio, using

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{\int {T}_{r\phi }{dV}}{\int {PdV}}=\displaystyle \frac{\int \rho v{{\prime} }_{\phi }v{{\prime} }_{r}{dV}}{\int {PdV}},\end{eqnarray} \tag{ 10 }$$

with differential volume dV. The turbulent components of the velocities indicated by the superscript primes are found by subtracting the azimuthally averaged value at each r and θ. The time evolution of α is presented in the top panel of Figure 2. Over the first 30 inner orbits, the value of α quickly increases and saturates to a new steady level near α = 1.5 × 10⁻⁴. Because the initial state already reflects VSI turbulence, saturation takes only about one-quarter as long as in model M3 (see Figure 3 in Flock et al. 2017b). We follow the evolution for 400 inner orbits. The fluctuations over time in the stress-to-pressure ratio remain small compared to the mean value.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The radial profile of the vertically and azimuthally averaged α value is shown in the bottom panel of Figure 2, where an additional time average between 100 and 400 orbits has been employed. The profile of α shows an increase, from 8 × 10⁻⁵ at 25 au to 2.5 × 10⁻⁴ at 35 au. Beyond 35 au the α value remains above 10⁻⁴. This increase of α from 25 to 35 au is connected to the thermal relaxation timescale relative to the threshold value for the VSI to operate. The small dip at 80 au might be related to the run's duration. The steady state might not be fully established near the outer radius, since the calculation ends before the outer boundary completes 36 orbits.

The turbulent flow speed ${v}_{\mathrm{rms}}={\left\{{({v}_{r}^{{\prime} })}^{2}+{({v}_{\theta }^{{\prime} })}^{2}+{({v}_{\phi }^{{\prime} })}^{2}\right\}}^{1/2}$ near the midplane well outside the 35 au surface density dip has the azimuthal power spectrum shown in Figure 3. This spectrum is averaged over the volume between 50 and 60 au in radius and within 0.14 radians of the midplane in θ. The spectrum is further time averaged from 100 to 400 inner orbits, and normalized by its maximum value. The high spatial resolution means we cover three decades in azimuthal wavenumber. The turbulent spectral energy peaks at the scale of the circumference, with a shallow decline toward shorter scales. For azimuthal mode numbers m of 10–100 the slope turns steeper, passing through m^−5/3. Between m = 100 and 1000 the slope turns steeper again, and the turbulent motions at the smallest resolved scales carry four orders of magnitude less spectral power than those at the largest azimuthal scales. Manger & Klahr (2018) found a similar spectral shape in their VSI turbulence models, though the steepening seen here at m = 100 occurs in their models at m = 40, likely due to the lower numerical resolution, steeper radial pressure profile, or differing thermal relaxation timescale. Manger & Klahr (2018) used a uniform value t_relax = 2 × 10⁻³/Ω, comparable to our smallest value here on the midplane.

Figure 4 shows a snapshot of the total θ-velocity at 290 inner orbits. Cuts along the midplane and a representative meridional plane show the 3D structure. The characteristic pattern of the VSI is visible, with large-scale upward and downward motions crossing the midplane and reaching the disk atmosphere. Speeds are as great as 100 m s⁻¹ in the atmosphere and the midplane motions are almost as fast, with amplitudes up to 50 m s⁻¹. The midplane cut makes clear that the motions are approximately axisymmetric, with small but noticeable deviations especially at low wavenumbers.

Zoom In Zoom Out Reset image size

Here we investigate the formation and evolution of vortices in the simulation. We identify vortices using the vertical component of the vortensity $\tfrac{{({\rm{\nabla }}\times v)}_{z}}{{\rm{\Sigma }}}$. Figure 5 shows a map of the vortensity in the r–ϕ midplane after 300 inner orbits. A single large vortex is visible as a dark region with a minimum in the vortensity around 30 au. The vortex is located radially inwards from the gas gap, which is visible as a bright ring. The vortex is fully developed after 150 inner orbits. The position and size of the vortex are listed in Table 2. We define the vortex as the region where $\tfrac{{({\rm{\nabla }}\times v)}_{z}}{{\rm{\Sigma }}}\lt {10}^{-11}\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}\,{{\rm{s}}}^{-1}$. The radial extent fluctuates between 3.6 and 4.6 au, the azimuthal extent between 35 and 45 au, and the aspect ratio between χ = 8.7 and χ = 12.6. This vortex is generated by Rossby wave instability and occurs in an annulus where the instability criterion of Li et al. (2000) is satisfied as we show in Appendix B. The detailed time evolution of the vortensity at 30 au is plotted in Figure 6. During the first 50 inner orbits, the eightfold symmetry from replicating run M3 is broken and a new steady state develops. Small vortices form and merge up to 100 inner orbits. Thereafter a single large vortex is present. To lessen this vortex's drift sideways on the plot, we shift each row in azimuth according to the orbital speed at 30 au. The vortex migrates only very slowly in radius (Table 2). A number of factors influence the speed at which a vortex migrates, including its size, strength, and aspect ratio (Paardekooper et al. 2010). For aspect ratios $\chi \geqslant 7$, Paardekooper et al. (2010) showed that vortex migration can be very slow, being two orders of magnitude slower than for a vortex with aspect ratio χ = 2, for example, and hence it is likely the large aspect ratio of the vortex formed in the simulation (χ ∼ 10) that explains its slow migration. Also important to note is that the vortex does not decay and is present for the whole remaining duration of the simulation. In contrast to our previous models of MHD turbulence arising from the magnetorotational instability (MRI), where a vortex was repeatedly formed and destroyed at the MRI dead-zone outer edge (Flock et al. 2015), here the vortex is not destroyed. Ignoring MHD effects, Lesur & Papaloizou (2009) showed that elliptical vortices are unstable with the exception of vortices with an aspect ratio between 4 and 6, which are stable. However, the growth rate of the elliptical instability that could destroy the vortex remains small. For the aspect ratio of the large vortex we predict a growth time of about 100 orbits at the vortex position (Lesur & Papaloizou 2009). Also, as this vortex is sustained by the RWI it might not be subject to the elliptical instability (Lyra 2013). Further longer-term studies are needed to test the stability of this vortex, but close inspection of Figure 5 indicates that the flow inside the vortex is not smooth, suggesting that the elliptical instability is operating to generate turbulence inside the vortex, but the tendency of this interior flow to destroy the vortex is counteracted by the continuous driving provided by the RWI due to the persistent pressure bump.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Many small, short-lived vortices present in our simulation are visible in Figure 5. These appear similar to those reported by Richard et al. (2016), and likely have the same origin, namely small-scale vorticity/vortensity perturbations generated by the VSI that undergo, local Rossby wave instability. In Appendix C, we isolate and observe the formation and evolution of such a small vortensity minimum. The life time is found to be approximately 1.4 local orbits, and this short life time is likely due to the elliptical instability acting strongly in these vortices that have small aspect ratios. As discussed in the next section, small local dust concentrations arise in these vortices, but are limited by the short life times.

In this section we describe the results from the evolution of the 100 μm and 1 mm grains. A snapshot of the 1 mm grains' 50 inner orbits after the grains were added is shown in Figures 7 and 8, close to edge-on and face-on, respectively. Most of the 1 mm grains remain within ±4 au of the midplane. In the outer regions, an oscillatory structure becomes visible due to the combination of radial drift and the vertical motions induced by the VSI. We note that in the outer disk the grains are not yet fully vertically mixed. For the 1 mm grains we measure a mean inwards drift speed of 7.5 m s⁻¹ outside of 70 au, which roughly corresponds to 0.2 au per inner orbit. This value is only slightly less than the analytical estimate of 9 m s⁻¹ using the initial profiles. Figure 8 shows the distribution in a face-on projection at the same time. The distribution of 1 mm grains clearly shows the effect of the surface density depression with a lack of particles in this region. The vortex at the inner edge of the steep drop is clearly visible as a concentration of the 1 mm grains to the southeast. The other disk regions exhibit small-scale concentrations of particles. These concentrations are connected to the short-lived vortices described in Section 3.2 and Richard et al. (2016). A fully edge-on view and comparison between both grain sizes can be found in Appendix D.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A critical ingredient for comparing 3D simulations of gas and dust with observations is the scale height of the dust grains. To measure the vertical distribution of the grains we consider all particles in the radial interval between 30 and 40 au. The particles are sampled by binning onto a high-resolution grid with 1000 cells in the θ direction. The distribution is calculated at each inner orbit and then averaged from 30 inner orbits after introducing the dust particles until the end of the run. As we have shown in Flock et al. (2017b), grains at this position have already been mixed well enough to describe a new steady profile. The results are shown in Figure 9. The plot shows that the 1 mm grains are mixed up very efficiently. We fit the dust scale height h_d = H_d/R with a value of ${h}_{d}^{1\,\mathrm{mm}}=0.034$. This value is much higher than has been inferred from observations of the HL Tau system (Pinte et al. 2016), with ${h}_{d}^{1\,\mathrm{mm}}=0.007$, and higher than what we found in our previous nonideal magnetohydrodynamical models (Flock et al. 2017b).

Zoom In Zoom Out Reset image size

We also compare our results with the well-known equilibrium solution for dust settling in a turbulent medium (Dubrulle et al. 1995)

$$\begin{eqnarray}&&{h}_{d}=\displaystyle \frac{{h}_{{\rm{g}}}}{\sqrt{1+\tfrac{{St}}{\widetilde{D}}}},\end{eqnarray} \tag{ 11 }$$

with the gas scale height in units of the radius h_g = H/R, the Stokes number St, and the dimensionless diffusion coefficient $\widetilde{D}$ defined as $D=\widetilde{D}{c}_{s}H$. The dimensionless diffusion coefficient is often expressed in units of the α parameter $\widetilde{D}=\alpha /\mathrm{Sc}$ (Dullemond & Dominik 2004; Schräpler & Henning 2004) with the Schmidt number Sc close to unity (Johansen & Klahr 2005; Turner et al. 2006). By using the values from our simulations for the 1 mm grains between 30 and 40 au, $\mathrm{St}=5\times {10}^{-2},\alpha =2\times {10}^{-4}$ we obtain h_d = 0.007, a similar value to what we found in our nonideal MHD simulations (Flock et al. 2017a) and to what has been inferred for the millimeter dust in HL Tau (Pinte et al. 2016). However, to obtain the value h_d = 0.034, which we directly fitted to the dust grain distribution, we need to assume ${\alpha }_{z}=5.4\times {10}^{-3}$, which tells us that the VSI is 27 times more efficient in mixing vertically than in transporting angular momentum radially. We therefore suggest that if the VSI is operating in the HL Tau disk, the vertical mixing might be quenched either by nonideal MHD effects by a magnetized wind layer (Cui & Bai 2020) or by dust drag due to a high dust to gas mass ratio (Pinte et al. 2016; Lin 2019). Recently Schäfer et al. (2020) obtained results showing the VSI operating together with a dust layer. The efficiency of the streaming instability in these environments remains still under investigation (Yang et al. 2018; Umurhan et al. 2020; Chen & Lin 2020; Schäfer et al. 2020). We discuss more on the dust scale height measurements in Section 4.7.

Following Fromang & Papaloizou (2006), another way to express the vertical diffusion is given by $D={\tau }_{\mathrm{eddy}}\left\langle {v}_{z}^{2}\right\rangle$. To compare with previous models by Stoll & Kley (2016) we determine the eddy timescales τ_eddy. For the 1 mm grains we determine in the annulus between 30 and 40 au ${(\langle {v}_{z}^{2}\rangle )}^{1/2}\,=0.0166{c}_{s}$, which leads to ${\tau }_{\mathrm{eddy}}=0.209{{\rm{\Omega }}}^{-1}$. A similar value τ_eddy ≈ 0.2 was found by Stoll & Kley (2016) for the vertical diffusion of particles (see Section 4.3 therein).

Here we illustrate how the grains evolve by showing four representative particles, chosen at random from among the 0.1 mm particles, and four from among the 1 mm particles. Figure 10 shows the spatial and thermal evolution of the four 1 mm grains. All four show a clear radial drift over the simulation runtime. The time is now given in years as their radial location is constantly changing. The 1 mm grains drift inwards approximately 20 au in 10,000 yr. Figure 10, top left panel, also shows that locally the radial drift can be reduced (green dashed–dotted line) or can halt inside the vortex (blue dashed line). The vertical distribution of the four grains is shown in Figure 10, top right panel. The VSI's characteristic upward and downward motions are clearly visible. At these distances, millimeter grains need around 5000–10,000 yr to travel from one hemisphere to the other, around 2 au above and below the midplane. None of the four millimeter grains make excursions up into the regions of the disk that are strongly heated by the star, as is clear from their temperature evolution. Figure 10, bottom panel, shows the temperature evolution of the grains. We assume the dust and gas temperatures to be equal. All four grains warm up as they drift inward. The grain that gets trapped in the vortex (blue dashed line) shows oscillations in temperature around 40 K, connected to the radial oscillations due to the vortex rotation. The evolution of four individual 100 μm grains is shown in Figure 11. The top left panel shows the radial locations over time. As expected the radial drift of these grains is strongly reduced compared to the 1 mm grains. Three of the grains show no substantial radial drift, while one drifts about 5 au in 20,000 years (blue dashed line). Figure 11, top right panel, shows the vertical locations of the grains over time. Due to the their smaller size they are able to reach higher altitudes in the disk, spanning a region ±5 au about the midplane, showing a similar pattern of upward and downward motions to the 1 mm grains. The 100 μm grains also remain below the disk's irradiated surface layers, as shown by the temperature evolution in Figure 11, bottom panel. The warm starlight-irradiated layers begin about 2.2 scale heights or 12 au above and below the midplane at 50 au.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To obtain the characteristic turbulent velocities of the grains we first perform a time average over the full simulation time of each velocity component. This temporal mean is subtracted from each particle's velocity at each time step. This way the radial drift is removed from the radial component. The azimuthal component changes erratically with time, such that a mean azimuthal velocity is not easily chosen. We therefore focus here on the radial and vertical turbulent velocity components and note that these provide a lower limit on the turbulent flows' speed. Figure 12 shows the time evolution of the turbulent speed for an individual 1 mm grain. The rms and standard deviation of all remaining 1 mm grains are overplotted for reference. We determine a value of v_rms^{r, θ} = 11.7 ± 6.7 m s⁻¹. The 100 μm grains show a slightly higher turbulence of 14.3 ± 8.0 m s⁻¹. The representative individual particle plotted in Figure 12 shows that local grains are accelerated quite substantially to several tens of meters per second. For the chosen millimeter grain, for a short time period the speed exceeds 30 m s⁻¹.

Zoom In Zoom Out Reset image size

With our high resolution simulations covering all 360° of azimuth we reach a turbulence level of $\alpha \cong {10}^{-4}$. We note that this value is around 4 times higher than what we have found for our previous model with a limited azimuthal extent. This result is not surprising as larger modes can now fit into the domain and it is known that the VSI drives perturbations with large azimuthal scales (see also the spectrum analysis in Figure 3). Our value of α is similar to previous radiation HD simulations covering 90° in azimuth by Stoll & Kley (2014), which also found α ≅ 10⁻⁴. Isothermal simulations by Nelson et al. (2013) and Manger & Klahr (2018) reported higher values with $\alpha \cong {10}^{-3}$. Such higher values are expected for isothermal simulations as those present the upper limit of fastest cooling times. In addition we note that the VSI strength depends on the disk's aspect ratio (Nelson et al. 2013) with stronger turbulence for larger scale heights. Our aspect ratio H/R ranges from 0.09 at the inner radial boundary to 0.14 at the outer boundary. We also highlight the different nature of VSI turbulence compared with the turbulence driven by MRI. The azimuthal Fourier power spectrum we obtain for the VSI is similar to previous results. The slope is shallow between m = 1 and m = 10, steepening to m^−5/3 between m = 10 and m = 100, and steepening further between m = 100 and m = 1000. Manger & Klahr (2018) reported a similar profile for wavenumbers up to 400.

Our simulation results show increasing VSI activity with radius out to 35 au, after which the strength of the induced turbulence levels off. The increase is connected to the fact that the relaxation timescale is a decreasing fraction of the orbital period as we move away from the star, falling further below the threshold value below which the VSI operates. Moving radially outwards, the disk becomes less optically thick, allowing faster thermal relaxation and strengthening the VSI. Nelson et al. (2013) found the threshold for the relaxation timescale to be ${t}_{\mathrm{relax}}=0.1{{\rm{\Omega }}}^{-1}$. For t_relax = 0.1 they found that the VSI does not operate in the disk while decreasing the relaxation time by only one order of magnitude was enough to show a fully evolved turbulence (compare Figure 12 in Nelson et al. 2013). In Figure 13, we compare the critical timescale of the VSI with the thermal relaxation timescale as described in Lin & Youdin (2015) and adopted by Malygin et al. (2017) and Lyra & Umurhan (2019). For our model we expect the VSI to be quenched inwards of 20 au. We find that between 20 and 30 au the turbulence strengthens, saturating at around 35 au (see Figure 2). This turbulence profile triggers a modification of the surface density profile, and the structure that develops leads to the subsequent generation of a vortex by the RWI. It is interesting to note that such a process of surface density perturbation is very similar to the vortex generation mechanism in a magnetized disk, which has a dead-zone in the inner disk but that supports magnetorotational turbulence in the outer disk. Here the change in turbulence activity also leads to the formation of a vortex at the outer edge of the dead-zone (Flock et al. 2015). We interpret our new result of a large-scale and persistent vortex forming in a nonmagnetized disk in terms of a VSI dead-zone transitioning to an active zone once the thermal equilibrium timescale reaches the critical timescale of the VSI. We note that the exact location and size of the vortex might be influenced by the initial surface density profile from which we restarted. We expect that the VSI transition region should also occur in simulations with a pure β cooling, which depends on radius and height. 2D maps of the thermal relaxation timescale were presented by Pfeil & Klahr (2019), and these could easily be implemented within the β cooling framework. Future studies should test this transition regime of the VSI and its possible impact on ring formation in protoplanetary disks.

Zoom In Zoom Out Reset image size

The position of the inner edge of VSI activity depends on the optical depth and so on the total amount of dust. For a lower dust abundance, the VSI activity is shifted inwards (Flock et al. 2017b) while for a higher dust amount it is shifted outward. For our model we estimate the inner edge of the VSI activity to be at about 45, 18, and 7 au for dust-to-gas mass ratios of small grains of 10⁻², 10⁻³, and 10⁻⁴.

The VSI activity is sensitive to the total optical depth and thus the local dust abundance. For example, a small local random fluctuation in the accretion stress may lead to the formation of neighboring disk annuli with higher and lower surface densities. The resulting changes in the local optical depths might reinforce the variation in accretion stress, allowing the changes to amplify. Such a process may lead to the formation of ringlike perturbations in the surface density profile, as demonstrated by Dullemond & Penzlin (2018). Furthermore, substantial concentration of grains may have a positive dynamical feedback effect, as the additional inertia the grains provide might further reduce the turbulent activity in the higher surface density regions relative to the surrounding lower density regions (Lin 2019).

Our model develops a large-scale, long-lived vortex at around 30 au. This vortex appears because of the surface density perturbation arising from the spatial variation in turbulent activity. The vortex aspect ratio of χ ∼ 8 is very similar to the value reported in Manger & Klahr (2018), who also saw large-scale, long-lived vortices in their simulations of the VSI operating in disks with very short cooling timescales. Unlike their vortex, the one in our model does not migrate because it is generated and maintained at a single persistent feature in the surface density profile. The vortices in Manger & Klahr (2018) appear to be generated by the RWI operating at the locations of intermittent surface density features that arise as natural fluctuations in the flow. We do not see the spontaneous emergence of large-scale vortices by this mechanism in our simulation, likely due to the fact that the cooling timescales adopted by Manger & Klahr (2018) lead to a more vigorously operating VSI with a larger α value of ∼10⁻³, compared to the smaller value of 10⁻⁴ we obtain in our simulations. Fluctuations around this larger value of α likely lead to larger fluctuations in the surface density profile, allowing the emergence of multiple large-scale vortices. Manger & Klahr (2018) also use steeper power-law radial declines in density and temperature, strengthening VSI activity and vortex formation.

In addition to the large persistent vortex, we observe the formation of many smaller vortices that have short life times, on the order of the orbital period, because of their small aspect ratios. Both Richard et al. (2016) and Manger & Klahr (2018) have reported the existence of similar small and short-lived vortices in simulations of the VSI. These could arise from secondary instability because of the RWI or Kelvin–Helmholtz instability, which acts on small-scale vortensity perturbations generated by the VSI (Latter & Papaloizou 2018). The formation of small, tightly wrapped vortices with small aspect ratios allows the elliptical instability to operate effectively, such that the vortices quickly dissolve into the background turbulent flow. Vorticity capturing schemes, such as that presented by Seligman & Laughlin (2017), might be well suited for studying such small-scale vortices. Finally, we note that the reduced resolution of 35 cells per H in azimuth should nonetheless be high enough to study the formation and evolution of vortices. For comparison Manger & Klahr (2018) used a resolution of 12 cells per H.

In a recent work Krijt et al. (2018) found that dust growth and settling can lead to a depletion of CO from the warm molecular layer located at higher altitudes in the disk. Connected to this idea we investigated whether the strong upward motions of the VSI can lift 100 μm sized grains into these warm molecular layers. However, none of the half-million 100 μm grains in our models reach these irradiated layers. Nonetheless, the idea that grains in the cold midplane can be mixed up into the warm molecular layer, and release molecules to the gas phase, could still work for smaller grains or grains that are more fluffy. The efficiency of CO depletion, as presented in Krijt et al. (2018), should be investigated for the VSI in the future.

The growth of planetesimals and planetary embryos by the accretion of pebbles is an important part of planet formation (Johansen et al. 2015). For efficient pebble accretion the pebbles' scale height should be smaller than the Hill sphere of the planet, H_d < R_H. The aspect ratio of our millimeter particles h_d = 0.034 means that a planet at r_p = 50 au has Hill radius ${R}_{{\rm{H}}}={r}_{p}{\left(\displaystyle \frac{{q}_{p}}{3}\right)}^{1/3}$ exceeding the particles' scale height if the planet's mass is greater than q_p = 1.2 × 10⁻⁴ times the stellar mass. For our 0.5 ${M}_{\odot }$ star this translates to efficient accretion of millimeter-sized pebbles by planets with more than $20{M}_{\oplus }$. Picogna et al. (2018) investigated the pebble accretion efficiency in VSI turbulent disks and found a peak for Stokes number unity grains because of their more effective settling. Dust coagulation, fragmentation, and radial drift models predict typically a maximum Stokes number of between 10⁻² and 0.1 (Testi et al. 2014; Ueda et al. 2019), similar to our 100 μm and millimeter-sized grains. For these values, pebble accretion efficiency is reduced according to Picogna et al. (2018). Further study is needed of the pebble accretion efficiency on planets embedded in VSI turbulent disks. Our results indicate that the turbulent velocities induced by the VSI are likely to lead to the fragmentation of colliding millimeter-sized grains, and the scale height of millimeter-sized grains is too large to allow for efficient pebble accretion by terrestrial-mass planetary embryos. It is difficult to avoid concluding that the VSI is not conducive to either the early or late phases of planet formation. Both problems, however, could be overcome if strong local concentrations of dust grains can arise such that the VSI is suppressed by the additional inertia provided by the particles (Lin 2019).

Good measurements of the dust scale height are so far available only for a few of the brightest protoplanetary disks, assuming the dust opacity and the relevant grain sizes are correctly modeled. One of these is HL Tau, where the dust scale height is much smaller than the expected gas scale height, with a best fit h_d = 0.007 at 100 au for millimeter-sized grains (Pinte et al. 2016). Millimeter-sized grains' abundance is confirmed by the system's spectrum at millimeter to centimeter wavelengths (Carrasco-González et al. 2019). However, to compare our results to HL Tau, we must allow for the differing Stokes number. If the grain density and size are the same, the Stokes number is proportional to ${\text{St}}\,\sim \,1/({\rho }_{{\rm{g}}}h)$.

The gas density in the disk around HL Tau has not been observationally determined, so we estimate upper and lower limits. The upper limit comes from the Toomre Q parameter, which governs when the disk is gravitationally unstable. As the disk shows no m = 2 pattern or spiral, it is likely stable against self-gravity. Following Hasegawa et al. (2017), the gas surface density at 35 au then cannot exceed 100 g cm⁻², eight times our value. The lower limit we derive from HCO⁺ observations by Yen et al. (2016). An HCO⁺-to-H₂ ratio of 10⁻⁸ (Cleeves et al. 2015) would mean the HCO⁺ column density of 10¹⁶ cm⁻² at this position translates to a gas surface density of 3 g cm⁻². The gas scale height at 35 au from HL Tau is only slightly smaller than in our model, with a ratio ${(80K/30K)}^{1/2}\times {(0.5{M}_{\odot }/1.7{M}_{\odot })}^{1/2}\approx 0.9$. Taking all this together, for the upper gas density limit we estimate that the Stokes number for millimeter grains at this position in HL Tau is 7.2 times less than in our model, which by Equation (11) implies an even greater dust scale height of 0.072. For the lower gas density limit, which is around four times less than in our model, the Stokes number would be bigger by a factor of 4.5, and so the dust scale height would be 0.017, still significantly greater than the value of 0.007 measured for HL Tau.

In both gas density limits, the strong VSI mixing would lead to a greater dust scale height than observed in HL Tau. In the lower limit, the dust-to-gas mass ratio would be near 0.1, perhaps large enough that the dust's back-reaction on the gas could hamper VSI. However, it is unclear how a globally enhanced dust-to-gas ratio would be reached in such a young system.

If the VSI is absent from these regions of protostellar disks, what could be the reason? One possibility is that magnetic activity in the disk's atmosphere could lead to turbulence that interferes with VSI (Lin 2019; Cui & Bai 2020). Nelson et al. (2013) found that VSI is prevented when the magnetic accretion stress-to-pressure ratio is several times 10⁻⁴. Whether magnetic activity is possible at these locations is worth further investigation.