© M. Barraza-Alfaro et al. 2021

1 Introduction

Understanding the origin of turbulence in protoplanetary disks is a key step towards comprehending how circumstellar disks evolve and how planets form. Globally, it can regulate the radial transport of angular momentum outwards through the generation of an effective viscosity, allowing the gas to accrete onto the central star (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974; see also PPVI review chapter in Turner et al. 2014). Locally, turbulence can be crucial for different physical processes relevant for planet formation, for example, gas accretion onto protoplanets (Picogna et al. 2018), the formation and lifetime of vortices (Fu et al. 2014), dust settling (Fromang & Papaloizou 2006), and dust growth and concentration in pressure bumps (Ormel & Cuzzi 2007; Birnstiel et al. 2010; Pinilla et al. 2012a,b). Moreover, radial variations in the strength of the turbulence can trigger the formation of pressure bumps (Lyra et al. 2009; Regály et al. 2012; Flock et al. 2015).

First observations of disk turbulence from spatially unresolved line broadening with the Submillimeter Array (SMA) enabled Hughes et al. (2011) and Guilloteau et al. (2012) to derive turbulence upper limits for the disks TW Hya, HD 163296 and DM Tau of around 40, 300, and 130 m s⁻¹. Recently, the Atacama Large Millimeter/submillimeter Array (ALMA) has allowed for a direct study of the gas turbulence in the disk’s outer regions (≳30 au) for the first time via spatially resolved observations of molecular gas lines. From the search of nonthermal line broadening, weak levels of turbulence (≲5−10% the local sound speed) have been constrained in HD 163296, TW Hya, MWC 480, and V4046 Sgr protoplanetary disks (Flaherty et al. 2015, 2017, 2018, 2020; Teague et al. 2016, 2018b). Most of the protoplanetary disk regions are expected to be weakly ionized (Dzyurkevich et al. 2013; Desch & Turner 2015). Therefore, the magneto-rotational instability (MRI, Balbus & Hawley 1991; Hawley & Balbus 1991) is unlikely to be active in most disk regions. In this MRI “dead-zone”, pure hydrodynamical instabilities are expected to dominate as the turbulence source, where the prevailing mechanisms strongly depend on the local cooling timescale (see e.g., Lyra & Umurhan 2019, and references therein). The dominant transport of the angular momentum in the dead-zone is dependent on the magnetic field orientation and strength; it might also be produced by a magnetic driven wind (Turner et al. 2014; Bai 2017). However, these winds are typically launched at heights of around 5–7 scale heights above the midplane (Gressel et al. 2020) and thus they are at higher altitudes than the molecular line emitting region.

The vertical shear instability (VSI, Nelson et al. 2013) is a hydro-instability that operates in the outer disk regions with fast cooling (Pfeil & Klahr 2019), which produces anisotropic turbulence (Stoll et al. 2017). It arises naturally from the vertical gradient of the disk’s angular velocity, where the surface layers rotate more slowly than the midplane. The VSI works more efficiently in isothermal disks (Nelson et al. 2013), but can also be active in a radiative disk (Stoll & Kley 2014; Flock et al. 2017, 2020). Moreover, it can also be effective in a magnetized disk where nonideal MHD effects suppress the MRI (Cui & Bai 2020). The VSI produces large-scale gas motions in its saturated state, which correspond to the long-wavelength modes or “body modes”. The large-scale gas motions have a characteristic meridional circulation pattern, in which the vertical motions can produce stresses significantly larger than in the radial motions (Stoll et al. 2017). A small-scale component can also be seen as turbulent eddies (Flores-Rivera et al. 2020).

The effective Shakura & Sunyaev α viscosity generated in the radial direction by the VSI has been found to be small, with values between 10⁻⁴ and 10⁻³ (Nelson et al. 2013; Flock et al. 2017; Stoll & Kley 2014; Manger & Klahr 2018). This low level of turbulence is consistent with the upper limits constrained from ALMA molecular gas line observations (Flaherty et al. 2015, 2017, 2018, 2020; Teague et al. 2016)and with the recent estimates from CO evolution models (Trapman et al. 2020). Despite these observational constraints, no direct comparison between synthetic predictions and observations has been made.

Determining the turbulence level alone by unresolved observations is not enough to confirm the responsible mechanism. In order to advance in determining the dominant physical mechanism of the disk dynamics, observations that can resolve the spatial structure of the gas kinematics are necessary. Luckily, ALMA observations have allowed to directly probe the kinematics in disks and significant deviations from the expected sub-Keplerian motion in equilibrium (Perez et al. 2015; Pérez et al. 2018, 2020; Teague et al. 2018a, 2019; Pinte et al. 2018a, 2019, 2020; Casassus & Pérez 2019; Disk Dynamics Collaboration 2020), although most of the times these are interpreted as the dynamical imprint of Jupiter mass planets due to the lack of quantitative predictions from alternative scenarios.

In this work we study the observability of the VSI-triggered gas motions by producing synthetic ALMA observations of the CO gas line emission. We address which spatial and spectral resolutions are needed, and discuss how the identification of signatures of the VSI would add constraints to the disk properties. Especially the detection of large scale motions in gas kinematics may support or discard the VSI-origin of the gas turbulence.

The disk physical model, and description of the hydrodynamical simulations are presented in Sect. 2. The radiative transfer calculations and synthetic imaging are detailed in Sect. 3. In Sect. 4 we discuss our main results. Finally, we summarize our work and draw our conclusions in Sect. 5.

2 Hydrodynamical simulations

2.1 Physical model

We adopted the physical model of the disk from Nelson et al. (2013). The physical model is described using cylindrical coordinates (R, Z, ϕ). The hydrodynamical simulations described in Sect. 2 are performed in spherical coordinates (r, θ, ϕ).

As initial conditions, we construct a set of 2D axisymmetric profiles for the density and rotation velocity that fulfill hydrostatic equilibrium.

The initial conditions for density and angular velocity are given by the equilibrium solutions: $$\rho (R,Z)={\rho }_0{\left(\frac{R}{R_0}\right)}^p\text{exp}\left(\frac{GM}{c_{\text{s}}^2}\left[\frac{1}{\sqrt{R^2+Z^2}}-\frac{1}{R}\right]\right),$$(1) $$\Omega (R,Z)={\Omega }_{\text{K}}{\left[(p+q){\left(\frac{H}{R}\right)}^2+(1+q)-\frac{qR}{\sqrt{R^2+Z^2}}\right]}^{1/2},$$(2)

where Ω_K is the Keplerian angular velocity ${\Omega }_{\text{K}}=\sqrt{GM/R^3}$, H = c_s∕Ω_K is the local disk scale height, c_s the local sound speed, R₀ is the reference radius, set to be equal to the code unit of length, and p and q are the power-law profiles for the density and temperature. ρ₀ sets the value of the gas midplane density at r = R₀. The gas in our simulation is described by a locally isothermal equation of state, that is $P=\rho c_{\text{s}}^2$. The isothermal sound speed is proportional to the temperature, $c_{\text{s}}^2\propto T$. Therefore, q also represents the radial power-law index of $c_{\text{s}}^2(R)$: $$c_{\text{s}}^2(R)=c_0^2{\left(\frac{R}{R_0}\right)}^q.$$(3)

The disk pressure scale height follows: $$H(R)=H_0{\left(\frac{R}{R_0}\right)}^{(q+3)/2},$$(4)

where H₀ is the disk scale height at the reference radius R₀, set such that H₀∕R₀ = 0.1. In our simulations the midplane density slope is p = −1.0, while the temperatureslope is q = −0.5, resulting in a flared disk with a flaring index of 0.25. Looking at the last term of the prescription given in Eq. (2), we observe that for the same cylindrical radius the upper layers of thedisk have slower angular velocity than the midplane; that is, the disk has a vertical shear, which sets the instability in our simulations.

2.2 Numerical method

We run global hydrodynamical simulations of an inviscid disk unstable to the VSI. The simulations are performed with the public version of the Godunov-grid-based code PLUTO¹ (version 4.3 Mignone et al. 2007).

We use the HD module of the PLUTO code to solve the Navier-Stokes equations of classical fluid dynamics: $$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho v)=0$$(5) $$\frac{\partial (\rho v)}{\partial t}+\nabla \cdot (\rho v{v}^{\text{T}})=-\nabla P-\rho \nabla \Phi ,$$(6)

were ρ is the gas mass density, v is the gas velocity vector, P is the pressure, and Φ is the gravitational potential. The fluid is affected by the gravitational potentials of a star only, Φ_⋆ = −GM_⋆∕r. The hydrodynamical equations were solved using a second-order accurate scheme with linear spatial reconstruction (LINEAR), second-order Runge-Kutta time-stepping, and the Harten-Lax-van Leer-Contact (HLLC) Riemann solver. The Courant number is set to 0.25.

The computational domain extends from 0.4 to 2.5 code units of length in the radial direction. In the vertical direction the model covers ~ 9 disk scale heights in total, ~4.5 at each side of the equator. In the azimuthal direction the grid covers the full 2π rad. The spherical grid is logarithmically spaced in the radial direction r, while evenly-spaced in colatitude and azimuth. The resolution of the grid is (r, θ, ϕ) = (512, 128, 1024), which gives a resolution of 14 cells per scale height and cells with an aspect ratio δr: rδθ:rδϕ of approximately 1: 1.9: 1.7. Since our simulation follows a locally isothermal equation of state, the simulation’s physical radius can be re-scaled. However, the reference aspect ratio of H₀ ∕R₀ = 0.1 set in our simulations fits best for a T Tauri disk model re-scaled to around 100 au. We display a summary of parameters used in the simulations in the upper section of Table 1.

We adopted similar boundary conditions as Flock et al. (2017). In colatitude, it corresponds to a modified outflow boundary condition which enforces zero inflow of gas into the computational domain. In the density, it extrapolates the logarithmic field along the meridional direction into the ghost cells. Additionally, we include a softening of the azimuthal velocity radial gradient, only at the meridional interface between the ghost cells and the computational domain to reduce effects from the boundary corners. In the radial direction, enforced zero inflow of gas is imposed. A buffer zone is also applied close to the inner and outer radial boundaries, whose sizes are equal to 20% of the inner and outer edge radius, respectively. In the buffer zones, the density and radial velocity are damped to the initial profiles.

Fig. 1 Time evolution of the stress-to-pressure ratio αrϕ of the VSI-unstable 3D hydrodynamical simulation. We show the time evolution in units of the orbital timescale at R = 100 au. The dashed line indicates the time-averaged value of αrϕ between 200 and 300 orbits.

2.3 Simulation results

Figure 1 shows the temporal evolution of the α_rϕ value consisting of the radial Reynolds stress-to-pressure ratio. We determine α_rϕ using the pressure-weighted stress-to-pressure ratio $${\alpha }_{r\varphi }=\frac{{\displaystyle \int T_{r\varphi }}\text{d}V}{{\displaystyle \int P}\text{d}V}=\frac{{\displaystyle \int \rho }{v^{\prime }}_{\varphi }{v^{\prime }}_r\text{d}V}{{\displaystyle \int \rho }{c}_{\text{s}}^2\text{d}V},$$(7)

with the volume dV, and ${v^{\prime }}_{\varphi }$ and ${v^{\prime }}_r$ representingthe turbulent components of the velocities are determined by subtracting the azimuthally averaged value at each r and θ. Figure 1 shows the time evolution, with a growth from very small values (we set the initial radial and vertical velocities to zero) until reaching a saturated value after around 170 orbits of evolution with α_rϕ ~ 1.4 × 10⁻⁴. For the further post-processing we use a snapshot after 300 orbits. A 3D plot showing the 3D gas dynamics of the disk is shown in Fig. 2.

As a next step, we present detailed horizontal and vertical 2D plots out of the 3D dataset using the same snapshots after 300 orbits (measured at the code unit of length). These 2D plots are shown in Fig. 3. The fields have been computed re-scaling the code unit of length to 100 au, and considering a central Solar-mass star, therefore, the snapshot corresponds to a disk evolution time of ~0.3 Myr. Each row represents snapshots of the midplane, top row, at three scale heights middle row, and a vertical cut in the R-Z plane. The four columns present density and the three velocity components. In the left column of Fig. 3, we plot the gas density perturbations, (ρ − ρ₀)∕ρ₀. At the midplane and upper layers we see a turbulent structure with spirals and ring-like structures. At the midplane the density deviations are on a level of several percent of the background value. At the upper layers of the disk the density deviations reaches levels of up to 30%. We note that even though the fluctuations are strong in the upper layers, there is no particular morphology that could be easily identified as VSI-generated.

Figure 3, right columns, show the 3D velocity structure typical of a disk in which the VSI is active in its saturated state. Strong perturbations are present in the three velocity components. In the radial direction, the velocity shows turbulent behavior, thus, not well-defined structures are seen. At the midplane the radial velocities reach values of several 10’s of meters per second, however, we note that in the upper layers these radial motions can reach over 100’s of meters per second. Such large values become of interest especially when looking later on the line emission.

Similarly, very turbulent structures are present in the azimuthal velocity deviations from a disk in sub-Keplerian rotation (expected from an unperturbed disk) at the midplane. Close to the surface, however, axisymmetric rings are visible. Interestingly, in some of these rings the rotation is sub-Keplerian in one hemisphere while super-Keplerian in the other. We also note that similar to the radial perturbations, the azimuthal ones are also weakest in the midplane with values of around 10’s meters per second while they become larger at the upper layers.

We summarize that the density perturbations, radial velocity perturbations and the azimuthal velocity deviation from a disk in Keplerianrotation are weak close to the disk midplane, however, above three pressure scale heights they can reach up to 100’s of meters per second.

Even more interesting are the meridional velocity motions, seen in Fig. 3 third column. There are prominent large scale motions through the whole vertical extent of the disk reaching velocities of around 50 meters per second. These motions are mostly axisymmetric, and therefore they show a clearer signature compared to the radial and azimuthal components.

Fig. 2 3D rendering of a VSI-unstable protoplanetary disk. The color map shows the perturbed total velocity of the gas in units of the local sound speed. Only half of our simulation domain in azimuth is shown for visualization purposes.

3 Synthetic observations

To produce synthetic images of molecular line emission of a VSI-active disk, we post-process the outputs of our simulations with the Monte-Carlo radiative transfer code RADMC-3D² (version 0.41, Dullemond et al. 2012). We build our RADMC-3D input files partly based on the FARGO2RADMC3D³ script (Baruteau et al. 2019).

We explore as observable the spatially resolved velocity centroid map (also labeled line of sight velocity map or Moment one map) computed from synthetic CO channel maps. To study the optically thin and optically thick cases, we compute the synthetic maps for three different CO isotopologues, ¹²CO, ¹³CO and C¹⁸O. We chose the J = 2–1 transition for the three isotopologues, centered at ~ 230.538 GHz for ¹²CO, ~ 220.399 GHz for ¹³CO, and ~ 219.560 GHz for C¹⁸O. All three transitions are observable with ALMA Band 6.

3.1 Setup

As mentioned before, we can re-scale the simulation since it follows an isothermal equation of state. We assume a central star of 1 M_⊙, and re-scale our density and velocity fields of the hydrodynamical simulation to R₀ = 100 au, and use it as input into RADMC-3D. To improve the performance, we volume-average the simulation data onto a coarser grid, halving the grid resolution for the radiative transfer post-processing. We also extended our disk, including an inner disk that goes from 10 au to the simulation grid’s inner edge of 40 au, that follows the equilibrium solution used as the initial condition in the simulation (see Eqs. (1) and (2)). Therefore, our disk model extends from 10 to 250 au. Additionally, we set the gas density of the model so the total gas mass of the disk in molecular hydrogen is 0.05 M_⊙, which result in a gas surface density of ~3 g cm⁻² at 100 au.

We compute the dust temperature via thermal Monte Carlo and assume that gas and dust have the same temperature. We adopt a constant gas-to-dust mass ratio of 100. The dust is composed of 60% astrosilicates (Draine & Lee 1984) and 40% amporphous carbon (Li & Greenberg 1997). The optical constant of the mixture was calculated using the Bruggeman mixing formula, with a resulting intrinsic density of 2.7 g cm⁻³. Mie theory was used to compute the dust opacities (Bohren & Huffman 1983). In our model the dust size distribution is defined adopting 10 logarithmically spaced dust size bins, with sizes ranging between a minimum of 0.01 μm, and a maximum of 1 cm. The dust size distribution follows a power-law profile n(s) ∝ s^−p with p = 3.5, in which the mass in the i-th size bin is $M_{d,i}=\chi M_{\text{gas}}\times s_i^{4-p}/{\sum }_i{s}_i^{4-p}$, with s_i the size of the i-th dust bin, χ = 0.01 the dust-to-gas mass ratio and M_gas = 0.05 M_⊙ the gas mass in the disk model (see also Baruteau et al. 2019). We adopted a simplified dust vertical settling prescription in which the dust smaller than 10 μm has the same scale height as the gas, whereas the dust with sizes larger than 10 μm is settled towards the midplane, with the dust scale height of the i-th dust size bin following: $H_{d,i}=\Lambda H_{\text{g}}{(s_i/s_{\text{min}})}^{-0.1}$, where H_g is the scale height of the gas and s_min is the minimum dust size in the model. The settling parameter Λ is set to 0.2, selected to roughly match our mm-sized dust scale height with constraints from observations (e.g., Pinte et al. 2016; Villenave et al. 2020). The adopted dust scale heights are comparable to the obtained using standard diffusion models (e.g., Dubrulle et al. 1995) for the disk-averaged vertical alpha in our simulation α_zϕ ~ 5.1 × 10⁻⁴, computed for one hemisphere of the disk following Eq. (4) in Stoll et al. (2017). Global simulations of the VSI that including dust and gas dynamics constrain larger dust scale heights due to the vertical mixing produced by the VSI vertical motions (Stoll & Kley 2016; Flock et al. 2017, 2020). Nevertheless, in our models a different vertical distribution of dust would only change the vertical temperature gradient of the disk in our radiative transfer predictions, which would have a small effect on the observable kinematic signatures.

The central object surface temperature is T_* = 7000 K with a radius of 1 R_⊙, and we assume that it emits as a perfect Black-body. The obtained midplane temperature approximately matches the used in the hydrodynamical simulations, which follows a power law ∝ R^−0.5, with R the cylindrical radius. Given the assumed central star, the midplane temperature goes from ~ 40 K at 40 au to ~16 K at 250 au. The dust is not included in the image ray-tracing to avoid optical depth effects produced by the dust continuum in the line emission. We use 10⁹ photon packages to compute the dust temperature via thermal Monte Carlo including absorption only, while 10⁸ photon packages for the image ray-tracing. A summary of relevant parameters used in the radiative transfer calculations are shown in Table 1.

For the molecular abundances, we assume a constant fraction of ¹²CO relative to H₂ of 1 × 10⁻⁴. To obtain the abundance of ¹³CO and C¹⁸O isotopologues, we scale the ¹²CO abundance with the ISM isotope ratios [¹² C]/[¹³C] ~ 77 and [¹⁶O]/[¹⁸O] ~ 560, respectively (Wilson & Rood 1994). The line emission is computed assuming local-thermodynamic-equilibrium (LTE). The molecular data is from the LAMDA⁴ database (Schöier et al. 2005). Variations of CO abundance from freeze-out or photo-dissociation are not included in our model. Although we omit the effect of CO freeze-out, the dust temperature obtained via thermal Monte Carlo shows a vertical gradient, leading to stronger emission in the upper layers. Using the molecular abundance fraction above, we generate sets of 3D grids of CO number density, and the three components of the velocity. Finally, the gas temperature is derived using the dust temperature from RADMC3D (see above).

As the gas motions produced by the VSI are roughly axisymmetric, we only study the disk emission by varying the inclination of the disk (5°, 45° and 85°), while keeping the position angle (PA) fixed to 90° (east of north). We assume a distance to the source of 100 pc. The disk near side is the south, and its rotation is counter-clockwise with respect to the observer. The total integrated flux of the ¹²CO(2−1), ¹³CO(2−1) and C¹⁸O(2−1) synthetic datacubes are approximately 42, 23 and 13 Jy km s⁻¹, respectively.

We compute the synthetic cubes to have a total bandwidth of 6 km s⁻¹, with 600 channels of 0.01 km s⁻¹. Then, we average these to obtain data cubes with a resolution of 0.05 km s⁻¹, that is the current maximum resolution of ALMA, that we use in the analysis. For a first approximation of a synthetic observation, we convolved the raw resulting synthetic channel maps by a circular-Gaussian beam, and study the addition of noise following the white noise model implemented in FARGO2RADMC3D (see also Baruteau et al. 2019). A mask of 0.42 arcseconds (42 au) was applied to all raw images to remove the emission from the disk regions close to the simulation’s inner edge. Additional data cubes varying the velocity resolution and spatial resolution are computed to explore the dependency of our results with resolution which we discuss in Sect. 3.4. The emission of four different velocity channels convolved by a circular Gaussian beam of 50 mas (5 au) is shown in Fig. A.1, for each disk inclination. It is readily seen that the VSI produces velocity perturbations in individual channels, seen as a corrugated pattern at the edges of the emission, more prominent with decreasing disk inclination.

Fig. 3 Results of the 3D hydrodynamical simulation. The panels show snapshots of the density and the velocities, taken after 300 orbital periods. From left to right: perturbed gas density, radial velocity, meridional velocity, and azimuthal velocity deviations from a sub-Keplerian velocity equilibrium solution. The velocity values are in physical units, where the simulation code unit of length is re-scaled to 100 au, and the central star is a Solar-mass star. From top to bottom: r-ϕ slices of the fields at the disk midplane, r-ϕ slices of the fields at approximately three pressure scale heights (3H) above the midplane, and r-θ slices of the fields. In the left-most column plots, the circles in dashed lines indicate the location of the inner and outer buffer zones. The dotted lines in the r-ϕ slices indicate the azimuth where the r-θ slice is taken, while in the r-θ plots the dotted line indicates the location where z ~ 3H.

3.2 ALMA synthetic observations

For a more accurate prediction of how our model would look in an interferometric ALMA observation, we simulate observations with CASA (McMullin et al. 2007) version 5.6 using the task simobserve. We use our ¹²CO(2−1) synthetic channel maps as input to simobserve. We simulate a 20 h integration in configuration C43-8 (8.5 km longest baseline) combined with a 4.4 h integration in C43-5 (1.4 km longest baseline). As described in the ALMA proposers guide, this configuration is ideal for a good uv-coverage in both short and long baselines.

We compute model visibilities using simobserve for each of the disk inclinations, with channel velocity resolution of 0.05 km s⁻¹. Then, we corrupt the visibilities with the task MS.CORRUPT() to obtain a RMS of ~ 1.5 mJy beam⁻¹ in each channel. This is roughly the expected RMS for a precipitable water vapor (PWV) of 0.9 mm and 20 h of integration time according to the ALMA sensitivity calculator. Finally, we compute CLEANed spectral cubes with CASA tclean. The restored beam has a FWHM of ~ 80 × 60 milliarcseconds (Briggs weighting 1.0), with a PA of −86.8°. The average signal-to-noise per beam and channel is ~10 at ~100 au from the central star, and larger than 5 in all the disk.

3.3 Velocity centroid maps and best fit sub-Keplerian disk model

We compute velocity centroid maps (v₀) from the ALMA simulated data cubes using a Gaussian function to fit the line emission in each pixel. For this we use the python package BETTERMOMENTS⁵ (Teague & Foreman-Mackey 2018; Teague 2019c). Alternatively, fitting a quadratic function to the line centroid could produce more precise results, however, a very high average signal-to-noise observation is needed for the method to work properly.

Following, we extract the velocity perturbations of the projected velocity image. For this purpose we use the Extracting Disk DYnamics python suite EDDY⁶ (Teague 2019a) to obtain the sub-Keplerian disk model that best fit our ¹²CO(2−1) synthetic velocity centroid maps.

The model assumes that the disk emitting surface is parametrized by: $$z(r)=z_0\times {\left(\frac{R}{1^{\prime \prime }}\right)}^{\psi }\times \exp \left(-{\left[\frac{R}{R_{\text{taper}}}\right]}^{q_{\text{taper}}}\right),$$(8)

where z₀ and R_taper are in arcseconds. The disk rotation curve and radial velocity follow a power-law profile: $$v_{\varphi }=v_{\varphi ,100}\times {\left(\frac{R}{100\text{au}}\right)}^{v_{\varphi ,q}}$$(9) $$v_R=v_{R,100}\times {\left(\frac{R}{100\text{au}}\right)}^{v_{R,q}},$$(10)

with R the cylindrical radius, and v_R,100 and v_ϕ,100 the disk radial and azimuthal velocity at 100 au from the star. The disk velocity model projected into the line of sight considers the contributions of both radial and azimuthal velocity components: $${\text{v}}_{\text{mod}}=v_{\varphi }\cdot \cos \varphi \cdot \sin i+v_R\cdot \sin \varphi \cdot \sin i+v_{\text{LSR}},$$(11)

where ϕ is the polar angle of the pixel (measured east of north relative to the red-shifted major axis) and v_LSR is the systemic velocity.

We set as free disk parameters the disk position angle, systemic velocity, disk center, the radial and azimuthal velocity at 100 au from the star, the slope of the radial and azimuthal velocity power-law profiles, and emission surface parameters (z₀, ψ, R_taper and q_taper). We fix the distance to the source to 100 pc and the disk inclination equal to the value used to compute the input data cube.

A series of MCMC chains are run to find the best fit sub-Keplerian model for a geometrically thick disk. We use 128 walkers that take 2000 burn-in steps and additional 500 steps to sample the posterior distributions for the model parameters. A delimited radial region of the disk is considered in the model fitting, set to [0.55,2.0], [0.58,1.85] and [0.6,1.7] arcseconds for inclinations of 5, 25 and 35 degrees, respectively. Finally, we extract the velocity perturbations subtracting the velocity centroid map of the best fit disk model (v_mod) to the original (v₀). The velocitydeviations results are presented in Fig. 4, and discussed in the following section.

3.4 Deviations from a sub-Keplerian disk model

The line of sight velocity maps, and extracted velocity deviations from a sub-Keplerian disk for three different disk inclinations (5°, 20° and 35°) are shown in Fig. 4. The maps are computed for line emission data cubes with a velocity resolution of 0.05 km s⁻¹. In the first column, the line centroid at each pixel is presented as the disk velocity in the line of sight (v₀), for images convolved by a circular Gaussian beam of 50 mas. Substantial deviations from the typical pattern from an inclined sub-Keplerian disk are already noticeable in the low inclination case. The following columns show the extracted velocity perturbations, which are seen as quasi axisymmetric rings that reach magnitudes of ~ 50 m s⁻¹. The substructures become more complex as we increase the disk inclination. In the second column, we present the velocity perturbations obtained when we subtract to v₀ shown in the first column, a second velocity centroid map computed from a synthetic observation of a smooth disk following an equilibrium solution (v_eq). We use herethe initial conditions from the simulation as input for the radiative transfer calculations, to compute the synthetic observations of the disk in equilibrium. The resulting residual map can be used as the expected velocity perturbations in an ideal observation with 5 au resolution. In the third column, the residual of subtracting to v₀, in this case computed for a synthetic cube convolved by a 10 au circular Gaussian beam, the best fit disk model obtained with EDDY (v_mod). As detailed in Sect. 3.3, the model corresponds to a geometrically thick disk, that follows power-law profiles for the velocities in the radial and azimuthal directions. The resulting residual maps are the observable velocity perturbations in a 10 au resolution observation, with the addition of white noise with an RMS level expected for a 20 h long-baseline ALMA observation. Recovering properly the emission surface and velocity profiles is harder for the highest inclination case presented (i =35°). We tested a higher inclination case with i = 45°, and we found we were unable to recover the velocity perturbations reliably. The extracted velocity perturbations match fairly well with the expected pattern for inclinations of 5 and 20 degrees. Additional modulations are present depending on the polar angle with respect to the disk’s major axis for the inclination of 35 degrees case. These modulations are due to systematic differences between the fitted and the true emission surface, with a secondary contribution from the errors on the disk center, PA, and velocity profiles (see also Fig. 12 in Yen & Gu 2020). In the fourth column, we applied the same procedure to extract the velocity perturbations shown in the third column, but in this case for a v₀ computed from the simulated ALMA observations described in Sect. 3.2. We observe that the perturbations are consistent with the expected VSI ringed structure, therefore, VSI-signatures are observable within ALMA capabilities.

How the different velocity components contribute to the projected line of sight velocity (LOSV) is crucial to interpret the velocity deviations. As shown by Teague et al. (2019) (see their Fig. 5), a ring of super(sub) Keplerian azimuthal motions show a sign-flip (from red- to blue-shifted, or vice versa) with respect to the PA of the disk. While for a ring of inward(outward) flow the sign flip is with respect to the line perpendicular to the disk PA. A ring of upward(downward) vertical motion, however, has the same sign for all PA. To better understand the contributions of each velocity component to the projected line of sight velocity, we computed the expected velocity perturbations considering each component separately. The synthetic data cubes were calculated using an input velocity field in which only one of the three velocity components is from a VSI unstable disk simulation and the remaining two are set to follow an equilibrium disk solution (v_eq). In Fig. A.3 we observe the contributions to the LOSV from the radial velocity v_r, meridional velocity v_θ and azimuthal velocity v_ϕ for different disk inclinations. It is clear that the meridional velocity dominates the LOSV. With increasing inclination, the radial and azimuthal velocities contribute more to the LOSV, where the radial velocity introduces asymmetric features in our highest inclination case. The contribution from the azimuthal velocity perturbations is minor for all the inclinations explored.

To summarize, the meridional flows from the VSI body modes are the most distinctive feature observable in the velocity centroid maps of a VSI-unstable disk. For low disk inclinations we have rings of red- or blue-shifted emission. Increasing the inclination we have a significant contribution of the radial velocity, in which the asymmetric nature of the radial perturbations adds complications to interpret the observable VSI signatures. Moreover, modulations of the extracted perturbations depending on the polar angle relative to the disk PA can appear due to the limitations of the method, making the characterization of the VSI-signatures difficult.

Fig. 4 Results of the line of sight velocity map and extracted velocity perturbations from a VSI unstable disk 12CO(2–1) synthetic lines observations. The velocity centroid of the line was computed at each pixel from mock data cubes with a velocity resolution of 0.05 km s−1. The input fields are shown in Fig. 3, which corresponds to a disk after ~ 0.3 Myr of evolution. From top to bottom: results for disk inclinations of 5°, 20° and 35°. First column: velocity centroid maps (v0). The images are convolved by a circular Gaussian beam of 50 mas and have no noise. Second column: residual map of subtracting to v0 the velocity centroid map obtained from a disk following an equilibrium solution (veq). Third column: residual of subtracting to v0 the best fit disk model obtained using EDDY (vmod). The images in this case are convolved by a 0.1 arcseconds circular Gaussian beam and have an RMS noise of ~ 1.5 mJy beam−1. Fourth column: same as the third column, but for a 20 h Cycle 7 ALMA simulated observation using configurations C43-8 and C43-5, with an RMS noise of ~ 1.5 mJy beam−1. The beam size is shown with a black circle at the bottom left of each panel. The black-dotted ellipses in panels of Cols. 3 and 4 are the inner and outer edge of the region considered to obtain the best fit model. The x- and y- axes indicate the RA and Dec angular offset from the central star’s position in arcseconds.

3.5 Detectability of VSI turbulence

A previous study has shown that a disk with isotropic turbulence can produce observable signatures of nonthermal gas motions in ALMA observations (Simon et al. 2015). First, it can produce a change in the shape of the spatially integrated line of emission, increasing the ratio of the peak line flux to the flux at line center. Such diagnosis has been applied to constrain disk turbulence in the outer regions of protoplanetary disks (see e.g., Flaherty et al. 2020 and references therein). We note that Simon et al. (2015) used a parametric fit to nonthermal gas motions to determine a turbulent broadening parameter as a function of radius and height. This parameter was then used as input in the line radiation transfer code LIME for the sake of computational efficiency. For our work we include the full 3D velocity field to calculate the line emission.

We explored if we can see the effect of line broadening by turbulence in our simulation data with anisotropic turbulence from the VSI. We compare directly the shape of the line of the synthetic observationsfor the VSI unstable disk, and a disk that follows a sub-Keplerian equilibrium solution with v_turb = 0, that is a laminar diskwith the radial and meridional velocities being zero. We assume the same dust temperature structure for both disks. The comparison is shown in Fig. 5 for data cubes with spatial resolution of 50 mas (5 au) and velocity resolution of 50 m s⁻¹. We observe that the effect of VSI turbulence in the line shape is negligible compared to the thermal broadening for all inclinations and CO isotopologues explored. We recovered similar results for a disk around a central star with lower temperature, analogous to a young Sun, with an effective temperature of 4300 K and radius of 2.6 R_⊙. This highlights how important is to use the full 3D velocity data to produce synthetic observation and study the line broadening in protoplanetary disks.

The second main result by Simon et al. (2015) showed that isotropic turbulence can affect the distribution of the emission in a given velocity channel. Flaherty et al. (2020) has shown that anisotropic turbulent motions can mimic such effect in spatially or spectrally unresolved observations. We observe that this effect is minor in our predictions, as the physical size and magnitude of the velocity perturbations is rather small (about 10 au and 50 m s⁻¹). The finger-like features are already indistinguishable in independent channels for a spatial resolution of about 30 au. We summarize again that a full 3D velocity field is needed when studying line-broadening. For ¹²CO, the thickness of the emission layer is smaller than a scale height. Furthermore, line broadening is only caused by velocity fluctuations with physical scales smaller than the depth of the emission layer. This fact causes the turbulent broadening to be negligible because the amplitudes of these velocity fluctuations are about one order of magnitude smaller than the large scale fluctuations.We stress that it might be more promising to determine the kinematics by spatially and spectrally resolved observations in contrast to determine the line broadening in unresolved disk observations. Our results of the α viscosity values from the hydrodynamical simulations are consistent with the upper limits from nonthermal broadening in ALMA observations, with α viscosity values ≲ 10⁻³. We compute in our simulation a disk-averaged radial α value of α_rϕ = 1.4 × 10⁻⁴, and a value about 3.6 times larger in the vertical direction (computed for one hemisphere of the disk following Eq. (4) in Stoll et al. 2017), i.e. α_zϕ ~ 3.6α_rϕ. However, we note again that these are not directly comparable as stated above, and going for higher resolution, or more sensitive observations would not result on the constraint of a turbulent α from VSI. We emphasize that spatially resolving the velocity perturbations from the VSI is required to confirm it as a source of turbulence. Moreover, it is feasible to estimate an α viscosity value of the disk by directly comparing the observed velocity structure with synthetic predictions from 3D hydrodynamical simulations. We conclude that the currently applied diagnostics to detect nonthermal turbulent motions in ALMA observations using a parametrization of the turbulence level are not reliable if the turbulence is anisotropic. Therefore, we can not discard that the VSI is active in these disks, and that turbulence is still playing an important role. Finally, despite our simulations do not resolve the smallest scales of the instability (Flores-Rivera et al. 2020), we do not expect a change in the turbulent line broadening when resolving these. The smaller turbulence scales are expected to have lower RMS velocities due to the turbulent cascade. Therefore, they do not significantly contribute to the broadening of the integrated line emission.

Fig. 5 Line profiles from synthetic observations for different CO isotopologues and disk inclinations. In each panel we show a comparison of the integrated line emission for a VSI unstable disk (solid lines) and a laminar disk (dashed lines).

4 Discussion

4.1 Effect of optical depth

Observing different CO isotopologues, allows to study different layers of the disk. To test the effect of the tracer optical depth, we run the radiative transfer calculations for ¹²CO, ¹³CO and C¹⁸O. The disk layers traced approximately by the different isotopologues are shown in Fig. A.5, where we display the layer where the disk becomes optically thick (τ = 1 surface). For simplicity, these calculations are for a disk face-on. The emissions traced in our predictions roughly match with constraints from observations (Zhang et al. 2017; Pinte et al. 2018b). The expected velocity perturbations from a VSI-unstable disk for the three different isotopologues is shown in Fig. A.4, for different disk inclinations.

We observe that the quasi axisymmetric ring structure present in the residuals of the ¹²CO predictions are also recovered for the ¹³CO and C¹⁸O predictions. Further, the remaining velocity residuals are feebler for optically thinner tracers, where the difference between isotopologues becomes more noticeable with increasing disk inclination. These findings are expected since the emission layer for these lines trace lower disk heights where the VSI perturbations are weaker. Moreover, the radial and azimuthal velocity components decrease faster than the vertical moving towards the midplane, as seen in the r-θ slices of the disk in Fig. 3, where is shown that the vertical motions from the VSI body modes penetrate through the whole vertical extent of the disk. Therefore, the meridional component of the velocity perturbations is predicted toalso dominate the line of sight velocity at deeper layers, resulting in the recovered quasi axisymmetric ringed morphology (as discussed in Sect. 3.4). Furthermore, the difference among tracers becomes larger for more inclined disks as the contribution to the line of sight velocity from the radial and azimuthal gas velocity increases with inclination. Additionally, the velocity structure is smeared for optically thinner tracers due to the emission arising from a larger range of heights (Zhang et al. 2017). Therefore, detecting VSI signatures in ¹³CO or C¹⁸O would require higher sensitivity, especially for inclined disks. Adding the difficulties to get high spatial resolution observations with enough signal-to-noise, we conclude that ¹²CO is the best tracer to detect VSI signatures. Nevertheless, observing the weak dependence with disk height of the morphology and magnitude of the VSI velocity perturbations for low inclined disks is essential to confirm its origin.

4.2 Differences between VSI signatures and other mechanisms

Recognizing VSI-driven gas motions is very important to interpret kinematics observations. Large-scale perturbations in the gas velocities can also be produced by other mechanisms. For example, massive planets (Disk Dynamics Collaboration 2020), vortices (Huang et al. 2018; Robert et al. 2020) or gravitational instabilities (Hall et al. 2020). Therefore, it is crucial to identify the features that can be used to separate the VSI signatures from other mechanisms.

Firstly, VSI signatures are quasi axisymmetric. Velocity perturbations of similar magnitude are visible in all velocity channels (see Fig. A.1), and show ringed substructures when extracted from the Moment one map (see Fig. 4). On the contrary, perturbations from the aforementioned alternative mechanisms are asymmetric, as trace spiral-arms or localized deviations. Planetary-wakes and vortices are expected to be strongest in a few velocity channels and show a localized signature in the velocity centroid maps (Perez et al. 2015; Pérez et al. 2018; Huang et al. 2018; Robert et al. 2020). Large-scale spiral arms, however, can also appear across all velocity channels, yet are distinguishable from rings when extracted from the velocity centroid maps (Pérez et al. 2018; Hall et al. 2020).

Secondly, the velocity component that contributes more to the line of sight velocity of the perturbations in a VSI-unstable disk is the vertical component. For planetary wakes and spiral-arm structures the largest contributions are from the radial and azimuthal velocities (Pérez et al. 2018; Hall et al. 2020). Meridional flows can also be generated in the gaps carved by embedded massive planets (Kley et al. 2001; Morbidelli et al. 2014; Fung & Chiang 2016; Teague et al. 2019). For Jupiter mass planets, the azimuthally averaged velocities at the gap and gap’s edges can have similar magnitude as the predicted VSI signatures (Teague et al. 2019). However, such flows are typically present together with spiral arms and perturbations localized around the planet’s location (Pérez et al. 2018; Juhász & Rosotti 2018). Therefore, resolved two-dimensional velocity centroid maps are required to disentangle between both scenarios.

Thirdly, the vertical motions of the VSI-modes are present in the whole disk vertical extent. Therefore, no significant change in the morphology of the observable VSI signatures are predicted when tracing different CO isotopologues. For planetary wakes and spiral arms, it is expected that the observable velocity perturbations change depending on disk height (Perez et al. 2015).

Fourthly, an embedded planet in a slow-cooling disk can generate additional tightly wound spiral arms due to buoyancy resonances. These spirals can have a strong vertical velocity component and produce observational signatures that can be confused with VSI signatures (Bae et al. 2021). Nevertheless, high spatial resolution observations can potentially disentangle between tightly wound spirals and VSI quasi axisymmetric rings in the line of sight velocity residuals.

Finally, VSI-perturbations are expected to produce dust traffic jams in the outer disk (Flock et al. 2017, 2020; Schäfer et al. 2020). It can also generate a dust trap at the radial transition where the cooling-timescale is short enough to sustain the VSI (Flock et al. 2020). Additionally, it can trigger Rossby Wave Instability vortices (Flock et al. 2020; Manger & Klahr 2018), which can efficiently trap dust at its center (Barge & Sommeria 1995). Similar structures in the dust distribution can be produced by embedded planets. Thus, it is difficult to disentangle between VSI and other mechanisms when studying the radial and azimuthal distribution of the dust continuum. Nevertheless, by studying the vertical distribution of the dust from edge on it ispossible to assess if dust vertical mixing from the VSI is present (Flock et al. 2017, 2020; Villenave et al. 2020).

In our simulation, short-lived small-scale Rossby Wave Instability vortices appear, consistent with previous results on VSI-unstable disks with very short thermal relaxation time scales (Richard et al. 2016; Flock et al. 2020). These vortices could produce kinematic signatures in CO line emission, as previously shown for large-scale planet-induced vortices (Huang et al. 2018; Robert et al. 2020). However, we focus on the global quasi-axisymmetric ringed structure produced by the VSI. Further analysis is needed to constrain the observability of vortices induced by the VSI in CO kinematics.

4.3 Constraints on the disk properties from a VSI-detection

Strong constraints can be done from the detection of the VSI in a protoplanetary disk. First, a VSI detection would confirm that the disk has a vertical shear and a radial gradient of the temperature. Second, a short cooling timescale is necessary for the VSI to operate in the traced layers of the disk. Via hydrodynamical simulations it has been shown that the VSI is active in disks with cooling timescales <10% of the local orbital timescale (Nelson et al. 2013; Lin & Youdin 2015; Flock et al. 2020). Hence, a VSI detection would constrain the disk cooling timescales. Finally, it also shows that the turbulence in the disk is fully or partly due to the VSI, which has been constrained from hydrodynamical simulations to α viscosity values in the range of 10⁻⁵–10⁻³. Damping of the VSI is predicted for the upper layers of the disk where the densities are low enough to allow collisional dust-gas decoupling (Pfeil & Klahr 2021) and magnetic effects dominate (Cui & Bai 2020). Therefore, studying the velocity perturbations at different layers of the disk can help us to understand changes in physical conditions at different disk heights.

4.4 Observations spatial and spectral resolution

We explored the effect of spatial and spectral resolution on the observable velocity perturbations for a disk with an inclination of 20 degrees. We display in Fig. A.2 the expected velocity perturbations from a VSI unstable disk varying the velocity resolution of the data cube, and size of the convolved circular Gaussian beam, in the first and third row, respectively. Additionally, the statistical uncertainties from the calculation of the line centroid maps derived using EDDY (Teague & Foreman-Mackey 2018) are shown in the third and fourth rows. Each data cube has added white noise with an RMS level expected from a 20 h integration ALMA observation. We observe that for beam sizes ≤ 10 au the VSI structure is well resolved. For such a high spatial resolution, all velocity resolutions explored reach an uncertainty level smaller than ~20% the channel width and the structure is recovered. However, only for the highest velocity resolution cases (channel width ≤ 50 m s⁻¹) these uncertainties are below ~20% the magnitude of the velocity perturbations. In summary, we found that high resolution observations using ALMA extended antenna configurations can spatially resolve the VSI signatures. However, the highest spectral resolution observations are needed for a robust detection.

For our predictions, with an assumed distance to the system of 100 pc, a 0.1 arcsecond resolution is enough. In ALMA band 6, the antenna configuration 8 gives the required resolution, where a 20 h integration time can get sufficient signal-to-noise ratio. Disks further away would require higher angular resolution to resolve the VSI, which would require a longer integration time. Our simulated observation setup is at the limit of the capabilities available for ALMA Cycle 8. Therefore, detection of VSI-signatures in the next ALMA cycle might be feasible for the brightest protoplanetary disks only.

Considering a standard Hanning smoothing, the highest resolutions are ~ 0.05 and ~ 0.03 km s⁻¹, for ALMA Bands 6 and7, respectively. According to our predictions, with these resolutions it is possible to resolve well the VSI signatures. Nonetheless, no spectral averaging and avoiding Hanning smoothing when preparing ALMA observations allows to achieve finer spectral resolution (down to ~0.025 km s⁻¹ for Band 7),which could ease the identification of VSI signatures.

We highlight that a good calibration is key to recover velocity perturbations from the data. For such a high spectral resolution observation to be successful, at least one of the basebands should be set up in TDM/wideband to get the best calibration. Otherwise, there may not be enough total bandwidth for best phase calibration, and this is truly essential to get the highest dynamic range images using self-calibration.

Our predictions are also valid for the CO J:3–2 transitions, observable within ALMA Band 7 in which higher velocity resolution observations are possible. Moreover, observing ¹²CO and ¹³CO within the same observation can be done for the J:3–2 transition. Obtaining the VSI signatures at different layers of the disk is important to better understand the vertical flows in the disk, and can potentially confirm the origin of the velocity perturbations.

The velocity resolution required to identify the VSI signatures could also vary depending on the method used to compute the line centroids and extract the velocity perturbations. Alternative methods to obtain the moment maps and the non-perturbed velocity centroid map are GMOMENTS⁷, and CONEROT⁸ (Casassus & Pérez 2019, Casassus et al., in prep.). The advantage of CONEROT compared to EDDY is that the perturbations can be extracted directly from observations without strong assumptions about the underlying disk model, and employing a reduced number of free parameters. Nevertheless, we do not expect significant differences in our results, but in its application to ALMA data.

Last, further study of the disk velocity structure can be performed analysing the velocity radial profiles. For this purpose CONEROT and GOFISH⁹ (Teague 2019b) can be applied. However, as discussed in Sect. 4.2, a thorough study of the three velocity components at different layers of the disk would be required to possibly disentangle between perturbations produced by planets or VSI.

5 Summary

We conducted a study of the observability of gas motions produced by the vertical shear instability (VSI) in ALMA gas kinematics observations. We explore the effect of disk inclination, spatial and spectral resolution. Furthermore, we compare the results for different CO isotopologues. We couple the results of 3D global hydrodynamical simulations of a disk unstable to the VSI, with radiative transfer calculations to obtain synthetic line emission predictions for ¹²CO, ¹³CO and C¹⁸O molecules. Next, we produce simulated observations that are directly comparable to those obtained by ALMA. Finally, we extract the velocity perturbations from our simulated observations, and recover the VSI perturbations. Our main findings are:

• We find that the VSI perturbations are seen as a corrugated pattern in the emission of individual velocity channels, and quasi axisymmetric concentric rings in the line of sight velocity residuals after subtracting a sub-Keplerian disk model fitted to the data.

• The characteristic morphology of the extracted perturbations results from the meridional velocity component being dominant in the velocity projected into the line of sight.

• With increasing inclination, the line of sight velocity has stronger addition of the radial and azimuthal velocity of the disk, adding asymmetries on top of the rings, making the VSI-perturbations challenging to characterize. Moreover, the method to extract the perturbations does not correctly recover the expected morphology for large inclinations (i ≳ 35°).

• The morphology of the extracted VSI velocity perturbations is predicted to be similar for different layers of the disk. The change in the magnitude of the expected velocity perturbations for different tracers remains relatively low for all inclinations explored. Further, the difference in magnitude increases with disk inclination.

• Our results show that the nonthermal broadening produced by the VSI in integrated line emission is negligible, consistent with current limits from ALMA radio observations. We emphasize that resolving the structure is fundamental to determine the turbulence sources in the disk.

• To recover the VSI-perturbations, spatial resolution below the disk pressure scale height might be enough to resolve the structures, which in our case corresponds to beam sizes ≲10 au. However, the spectral resolution needed to capture the perturbations reliably (≲0.05 km s⁻¹) is the highest available with the ALMA interferometer.

We highlight that our predictions are optimistic and the predicted observational signatures should be treated as upper limits. The conducted simulations are inviscid and follow a locally isothermal equation of state, giving the ideal conditions for the development of vigorous VSI motions in the disk. The inclusion of a finite cooling time or background viscosity can weaken the VSI perturbations. Moreover, the presence of velocity perturbations from a different origin (e.g., massive planets) can obstruct the identification of the VSI signatures. Further work is needed to study these scenarios.

Observing the gas kinematics at different layers of the disk is important to disentangle between the VSI and other mechanisms that could produce similar kinematic signatures. Finally, we conclude that the best cases to detect VSI signatures are gap-less disks close to face-on, in which the VSI meridional flows can be directly traced.

Acknowledgements

We thank the anonymous referee for constructive and detailed review of the manuscript. We thank F. Alarcón, N. Kurtovic and R. Teague, for helpful advice on the use of RADMC3D, CASA and bettermoments/eddy, respectively. Figures where produced by python matplotlib library (Hunter 2007), and ParaView (Hansen 2005; Ayachit 2015). M.B and M.F. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 757957). S.M. is supported by a Junior Research Fellowship from Jesus College, Cambridge. S.P. acknowledges support from ANID-FONDECYT grant 1191934. The numerical simulations were run on the HPC system COBRA at MPCDF (Max Planck Computing and Data Facility).

Appendix A Additional figures

Fig. A.1 Predictions of the 12CO disk emission of a VSI-unstable disk for different velocity channels. From left to right, the columns correspond to predictions for disk inclinations i = 5, 20 and 35 degrees. Thechannel widths are 0.05 km s−1. The images are convolved with a circular Gaussian beam of 50 mas, shown in the bottom-left corner of each panel. The disk emission shows a corrugated pattern. Only red-shifted channels are shown, as the structures are quasi axisymmetric, therefore, the symmetric blue-shifted channels show similar substructures.

Fig. A.2 Results of the expected observable velocity perturbations from a VSI unstable disk 12CO synthetic lines observations (v0 − veq), for various velocity and spatial resolution data cubes. Each channel of the cube has a RMS noise expected from a 20 h integration ALMA observation with an extended configuration (1.9, 1.5, 1.0, 0.8, 0.6 mJy from the higher to the lower velocity resolution cases). In addition, we show the uncertainties in the calculation of the velocity centroid map using  BETTERMOMENTS (δv0). The FWHM of the circular Gaussian beams are shown at the bottom left corner of each panel. The x- and y- axes indicate the RA and Dec angular offset from the position of the central star, in arcseconds.

Fig. A.3 Results of the expected observable velocity perturbations from 12CO synthetic lines observations, in which one of the velocity components is from a VSI unstable disk. The input VSI-active velocity field into the radiative transfer is modified to follow an equilibrium disk solution in two of the components. From left toright, the columns: predictions considering the contributions of the radial velocity, the meridional velocity and azimuthal velocity, respectively. From top to bottom: different rows show the results for disk inclinations of 5°, 20° and 35°. The mock data cubes have a velocity resolution of 0.05 km s−1, and the images are convolved by a circular Gaussian beam of 50 mas, shown at the bottom left corner of each panel. The x- and y-axes indicate the RA and Dec angular offset from the position of the central star, in arcseconds.

Fig. A.4 Results of the expected observable velocity perturbations from a VSI unstable disk synthetic lines observations for different CO isotopologues. From left to right, the columns show predictions for 12CO, 13CO and C18O. From top to bottom, the rows show the results for disk inclinations of 5°, 20° and 35°. The mock data cubes have a velocity resolution of 0.05 km s−1, and the images are convolved by a circular Gaussian beam of 50 mas, shown at the bottom left corner of each panel. The x- and y-axes indicate the RA and Dec angular offset from the position of the central star, in arcseconds.

Fig. A.5 Radial profile of the disk layers where the disk becomes optically thick (τ = 1 surface), for different CO isotopologues of the VSI-unstable disk radiative transfer model. We illustrate the τ = 1 surface for the J:2–1 and J:3–2 transitions of three CO isotopologues, 12CO, 13CO and C18O. We assume a face-on disk orientation for the calculation. The gray dashed lines show where Z = AR for A =0.1;0.2;0.3;0.4 and 0.5. The solid black line shows the grid boundary in colatitude of the simulated disk, located approximately at Z = 0.45 R.

References

All Tables

All Figures

	Fig. 1 Time evolution of the stress-to-pressure ratio αrϕ of the VSI-unstable 3D hydrodynamical simulation. We show the time evolution in units of the orbital timescale at R = 100 au. The dashed line indicates the time-averaged value of αrϕ between 200 and 300 orbits.
In the text

	Fig. 2 3D rendering of a VSI-unstable protoplanetary disk. The color map shows the perturbed total velocity of the gas in units of the local sound speed. Only half of our simulation domain in azimuth is shown for visualization purposes.
In the text

Fig. 3 Results of the 3D hydrodynamical simulation. The panels show snapshots of the density and the velocities, taken after 300 orbital periods. From left to right: perturbed gas density, radial velocity, meridional velocity, and azimuthal velocity deviations from a sub-Keplerian velocity equilibrium solution. The velocity values are in physical units, where the simulation code unit of length is re-scaled to 100 au, and the central star is a Solar-mass star. From top to bottom: r-ϕ slices of the fields at the disk midplane, r-ϕ slices of the fields at approximately three pressure scale heights (3H) above the midplane, and r-θ slices of the fields. In the left-most column plots, the circles in dashed lines indicate the location of the inner and outer buffer zones. The dotted lines in the r-ϕ slices indicate the azimuth where the r-θ slice is taken, while in the r-θ plots the dotted line indicates the location where z ~ 3H.

In the text

Fig. 4 Results of the line of sight velocity map and extracted velocity perturbations from a VSI unstable disk 12CO(2–1) synthetic lines observations. The velocity centroid of the line was computed at each pixel from mock data cubes with a velocity resolution of 0.05 km s−1. The input fields are shown in Fig. 3, which corresponds to a disk after ~ 0.3 Myr of evolution. From top to bottom: results for disk inclinations of 5°, 20° and 35°. First column: velocity centroid maps (v0). The images are convolved by a circular Gaussian beam of 50 mas and have no noise. Second column: residual map of subtracting to v0 the velocity centroid map obtained from a disk following an equilibrium solution (veq). Third column: residual of subtracting to v0 the best fit disk model obtained using EDDY (vmod). The images in this case are convolved by a 0.1 arcseconds circular Gaussian beam and have an RMS noise of ~ 1.5 mJy beam−1. Fourth column: same as the third column, but for a 20 h Cycle 7 ALMA simulated observation using configurations C43-8 and C43-5, with an RMS noise of ~ 1.5 mJy beam−1. The beam size is shown with a black circle at the bottom left of each panel. The black-dotted ellipses in panels of Cols. 3 and 4 are the inner and outer edge of the region considered to obtain the best fit model. The x- and y- axes indicate the RA and Dec angular offset from the central star’s position in arcseconds.

In the text

	Fig. 5 Line profiles from synthetic observations for different CO isotopologues and disk inclinations. In each panel we show a comparison of the integrated line emission for a VSI unstable disk (solid lines) and a laminar disk (dashed lines).
In the text

Fig. A.1 Predictions of the 12CO disk emission of a VSI-unstable disk for different velocity channels. From left to right, the columns correspond to predictions for disk inclinations i = 5, 20 and 35 degrees. Thechannel widths are 0.05 km s−1. The images are convolved with a circular Gaussian beam of 50 mas, shown in the bottom-left corner of each panel. The disk emission shows a corrugated pattern. Only red-shifted channels are shown, as the structures are quasi axisymmetric, therefore, the symmetric blue-shifted channels show similar substructures.

In the text

Fig. A.2 Results of the expected observable velocity perturbations from a VSI unstable disk 12CO synthetic lines observations (v0 − veq), for various velocity and spatial resolution data cubes. Each channel of the cube has a RMS noise expected from a 20 h integration ALMA observation with an extended configuration (1.9, 1.5, 1.0, 0.8, 0.6 mJy from the higher to the lower velocity resolution cases). In addition, we show the uncertainties in the calculation of the velocity centroid map using  BETTERMOMENTS (δv0). The FWHM of the circular Gaussian beams are shown at the bottom left corner of each panel. The x- and y- axes indicate the RA and Dec angular offset from the position of the central star, in arcseconds.

In the text

Fig. A.3 Results of the expected observable velocity perturbations from 12CO synthetic lines observations, in which one of the velocity components is from a VSI unstable disk. The input VSI-active velocity field into the radiative transfer is modified to follow an equilibrium disk solution in two of the components. From left toright, the columns: predictions considering the contributions of the radial velocity, the meridional velocity and azimuthal velocity, respectively. From top to bottom: different rows show the results for disk inclinations of 5°, 20° and 35°. The mock data cubes have a velocity resolution of 0.05 km s−1, and the images are convolved by a circular Gaussian beam of 50 mas, shown at the bottom left corner of each panel. The x- and y-axes indicate the RA and Dec angular offset from the position of the central star, in arcseconds.

In the text

Fig. A.4 Results of the expected observable velocity perturbations from a VSI unstable disk synthetic lines observations for different CO isotopologues. From left to right, the columns show predictions for 12CO, 13CO and C18O. From top to bottom, the rows show the results for disk inclinations of 5°, 20° and 35°. The mock data cubes have a velocity resolution of 0.05 km s−1, and the images are convolved by a circular Gaussian beam of 50 mas, shown at the bottom left corner of each panel. The x- and y-axes indicate the RA and Dec angular offset from the position of the central star, in arcseconds.

In the text

Fig. A.5 Radial profile of the disk layers where the disk becomes optically thick (τ = 1 surface), for different CO isotopologues of the VSI-unstable disk radiative transfer model. We illustrate the τ = 1 surface for the J:2–1 and J:3–2 transitions of three CO isotopologues, 12CO, 13CO and C18O. We assume a face-on disk orientation for the calculation. The gray dashed lines show where Z = AR for A =0.1;0.2;0.3;0.4 and 0.5. The solid black line shows the grid boundary in colatitude of the simulated disk, located approximately at Z = 0.45 R.

In the text