Observability of the Vertical Shear Instability in protoplanetary disk CO kinematics

Context. Dynamical and turbulent motions of gas in a protoplanetary disk are crucial for their evolution and affect planet formation. Recent observations suggest weak turbulence in the disk's outer regions. However, the physical mechanism of turbulence in these outer regions remains uncertain. The vertical shear instability (VSI) is a promising mechanism to produce turbulence in disks. Aims. Our aim is to study the observability of the gas velocity structure produced by the VSI via CO kinematics with ALMA. Methods. We perform global 3D hydrodynamical simulations of a VSI-unstable disk. We post-process the simulation results with radiative transfer calculations, and produce synthetic predictions of CO rotational emission lines. Following, we compute the line of sight velocity map, and its deviations from a sub-Keplerian equilibrium solution. We explore the detectability of the VSI by identifying kinematic signatures using realistic simulated observations. Results. Our 3D simulations of the VSI show the steady state dynamics of the gas in great detail. From the velocity structure we infer a turbulent stress value of $\alpha_{r\phi}=1.4 \times 10^{-4}$. On large scales, we observe velocity deviations of 50 m s$^{-1}$ as axisymmetric rings. We find optimal conditions at $i \lesssim 20^{\circ}$ to trace for the kinematic structures of the VSI. We found that current diagnostics to constrain gas turbulence from non-thermal broadening of the line emission are not applicable to anisotropic VSI turbulence. Conclusions. The detection of kinematic signatures produced by the VSI is possible with ALMA. Observations including an extended antenna configuration combined with the highest spectral resolution available are needed for robust detection. The characterization of the large-scale velocity perturbations is required to constrain the turbulence level produced by the VSI from gas observations.


Introduction
Understanding the origin of turbulence in protoplanetary disks is a key step towards the comprehension of 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 Member of the International Max-Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD), Germany 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) were able to derive turbulence upper limits for the disks TW Hya, HD 163296 and DM Tau of around 40, 300 and 130 m s −1 (Hughes et al. 2011;Guilloteau et al. 2012). Recently, the Atacama Large Millimeter/submillimeter Array (ALMA) has allowed 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 non-thermal 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 Teague et al. 2016;Flaherty et al. 2017Flaherty et al. , 2018Teague et al. 2018b;Flaherty et al. 2020). 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,    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). Dependent on the magnetic field orientation and strength, the dominant transport of angular moment in the dead-zone might be produced by a magnetic driven wind (Turner et al. 2014;Bai 2017). However, these winds are typically launched at heights of around 5 to 7 scale heights above the midplane (Gressel et al. 2020) and so at higher altitudes than the molecular line emitting region.
The Vertical-Shear-Instability (henceforth VSI, Nelson et al. 2013) is a hydro-instability that operates in the outer disk regions with fast cooling (Pfeil & Klahr 2019), that produces anisotropic turbulence (Stoll, Moritz H. R. 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. 2017Flock et al. , 2020. Moreover, it can also be effective in a magnetized disk where non-ideal 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, Moritz H. R. 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 −4 and 10 −3 (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 Teague et al. 2016;Flaherty et al. 2017Flaherty et al. , 2018Flaherty et al. , 2020 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;Teague et al. 2018a;Pinte et al. 2018a;Pérez et al. 2020;Casassus & Pérez 2019;Pinte et al. 2019;Teague et al. 2019;Pinte et al. 2020;Disk Dynamics Collaboration et al. 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 Section 2. The radiative transfer calculations and synthetic imaging are detailed in Sec-tion 3. In Section 4 we discuss our main results. Finally, we summarize our work and draw our conclusions in Section 5.

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 Section 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: where Ω K is the Keplerian angular velocity Ω K = GM/R 3 , H = c s /Ω K is the local disk scale height, c s the local sound speed, R 0 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. ρ 0 sets the value of the gas midplane density at r = R 0 . The gas in our simulation is described by a locally isothermal equation of state, i.e. P = ρ c 2 s . The isothermal sound speed is proportional to the temperature, c 2 s ∝ T . Therefore, q also represents the radial power-law index of c 2 s (R): The disk pressure scale height follows: where H 0 is the disc scale height at the reference radius R 0 , set such that H 0 /R 0 = 0.1. In our simulations the midplane density slope is p = −1.0 , while the temperature slope 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 equation 2, we observe that for the same cylindrical radius the upper layers of the disk have slower angular velocity than the midplane, i.e., the disk has a vertical shear, which sets the instability in our simulations.

Numerical method
We run global hydrodynamical simulations of an inviscid disk unstable to the Vertical Shear Instability. The simulations are performed with the public version of the Godunov-grid-based code PLUTO 1 (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: 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 0 /R 0 = 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 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

Simulation results
with the volume dV, and v φ and v r representing the turbulent components of the velocities are determined by subtracting the azimuthally averaged value at each r and θ. Fig. 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 −4 . 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 Figure 3. The fields have been computed re-scaling the code unit of length to 100 au, and considering a central Solarmass 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, (ρ − ρ 0 )/ρ 0 . At the midplane and upper layers we see a turbulent structure with spirals and ringlike 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 Keplerian rotation 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.

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 2 (version 0.41, Dullemond et al. 2012). We build our radmc-3d input files partly based on the fargo2radmc3d 3 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, 12 CO, 13 CO and C 18 O. We chose the J=2-1 transition for the three isotopologues, centered at ∼ 230.538 GHz for 12 CO, ∼ 220.399 GHz for 13 CO, and ∼ 219.560 GHz for C 18 O. All three transitions are observable with ALMA Band 6.

Setup
As mentioned before, we can re-scale the simulation since it follows an isothermal equation of state. We assume a central star of 1M , and re-scale our density and velocity fields of the hydrodynamical simulation to R 0 = 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 equations 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 −2 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 −3 . 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 , 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: 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 −4 , computed for one hemisphere of the disk following Eq. 4 in Stoll, Moritz H. R. 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. 2017Flock et al. , 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 9 photon packages to compute the dust temperature via thermal Monte Carlo including absorption only, while 10 8 photon packages for the image raytracing. 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 12 CO relative to H 2 of 1 × 10 −4 . To obtain the abundance of 13 (Wilson & Rood 1994). The line emission is computed assuming local-thermodynamic- 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. equilibrium (LTE). The molecular data is from the LAMDA 4 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 12 CO(2 − 1), 13 CO(2 − 1) and C 18 O(2 − 1) 4 https://home.strw.leidenuniv.nl/~moldata/ synthetic datacubes are approximately 42, 23 and 13 Jy km s −1 , respectively.
We compute the synthetic cubes to have a total bandwidth of 6 km s −1 , with 600 channels of 0.01 km s −1 . Then, we average these to obtain data cubes with a resolution of 0.05 km s −1 , i.e., 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 Section 3.4. The emission of four different velocity channels convolved by a circular Gaussian beam of 50 mas (5 au) is shown in Figure A.1, for each disk inclination. It is readily seen that the VSI produces velocity perturbations in Article number, page 5 of 17 A&A proofs: manuscript no. aanda individual channels, seen as a corrugated pattern at the edges of the emission, more prominent with decreasing disk inclination.

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 12 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 −1 . Then, we corrupt the visibilities with the task ms.corrupt() to obtain a RMS of ∼ 1.5 mJy beam −1 in each channel. This is roughly the expected RMS for a precipitable water vapour (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 mas (Briggs weighting 1.0). 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.

Velocity centroid maps and best fit sub-Keplerian disk model
We compute velocity centroid maps (v 0 ) 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 5 (Teague & Foreman-Mackey 2018; Teague 2019). 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 6 (Teague 2019a) to obtain the sub-Keplerian disk model that best fit our 12 CO(2 − 1) synthetic velocity centroid maps.
The model assumes that the disk emitting surface is parametrized by: where z 0 and R taper are in arcseconds. The disk rotation curve and radial velocity follow a power-law profile: 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: 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 0 , ψ, 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 0 ). The velocity deviations results are presented in Figure 4, and discussed in the following section.

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 Figure 4. The maps are computed for line emission data cubes with a velocity resolution of 0.05 km s −1 . In the first column, the line centroid at each pixel is presented as the disk velocity in the line of sight (v 0 ), 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 −1 . The sub-structures become more complex as we increase the disk inclination. In the second column, we present the velocity perturbations obtained when we subtract to v 0 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 here the 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 0 , 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 Section 3.3, the model corresponds to a geometrically thick disk, that follows powerlaw 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 20h longbaseline 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 sur- Second column: Residual map of subtracting to v 0 the velocity centroid map obtained from a disk following an equilibrium solution (v eq ). Third column: Residual of subtracting to v 0 the best fit disk model obtained using eddy (v mod ). The model corresponds to a geometrically thick disk, that follows power-law profiles for the velocities in the radial and azimuthal directions. 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: The 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 resulting beam size of our synthetic observation is 80 × 60 arcseconds with a PA of -86.8 • . The beam size is shown with a black circle at the bottom left of each panel. The black-dotted ellipses in panels of columns 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 R.A. and Dec. angular offset from the central star's position in arcseconds. face, with a secondary contribution from the errors on the disk center, PA, and velocity profiles (see also Figure 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 0 computed from the simulated ALMA observations described in Section 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 Figure 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 Figure 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.
Article number, page 7 of 17 A&A proofs: manuscript no. aanda 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.

Detectability of VSI turbulence
A previous study has shown that a disk with isotropic turbulence can produce observable signatures of non-thermal 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 non-thermal 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 observations for the VSI unstable disk, and a disk that follows a sub-Keplerian equilibrium solution with v turb = 0, i.e. a laminar disk with 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 −1 . 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 −1 ). 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 12 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 con-trast 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 non-thermal broadening in ALMA observations, with α viscosity values 10 −3 . We compute in our simulation a disk-averaged radial α value of α rφ = 1.4 × 10 −4 , and a value about 3.6 times larger in the vertical direction (computed for one hemisphere of the disk following Eq. 4 in Stoll, Moritz H. R. 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 non-thermal 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.

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 12 CO, 13 CO and C 18 O. The disk layers traced approximately by the different isotopologues are shown in Figure 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 Figure A.4, for different disk inclinations. We observe that the quasi axisymmetric ring structure present in the residuals of the 12 CO predictions are also recovered for the 13 CO and C 18 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 Figure 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 to also dominate the line of sight velocity at deeper layers, resulting in the recovered quasi axisymmetric ringed morphology (as discussed in Section 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 13 CO or C 18 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 12 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.

How to differenciate VSI signatures from 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 et al. 2020), vortices (Huang et al. 2018;Robert et al. 2020) or gravitational instabilities (Hall et al. 2020). A few points that can be used to separate the VSI signatures from other mechanisms are the following: -VSI signatures are quasi axisymmetric. Velocity perturbations of similar magnitude are visible in all velocity channels (see Figure A.1), and show ringed sub-structures when extracted from the Moment one map (see Figure 3.3). 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). -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 . 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 twodimensional velocity centroid maps are required to disentangle between both scenarios. -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). -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. -VSI-perturbations are expected to produce dust traffic jams in the outer disk (Flock et al. 2017Schä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 . Additionally, it can trigger Rossby Wave Instability vortices 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 is possible to assess if dust vertical mixing from the VSI is present (Flock et al. 2017Villenave et al. 2020).
In our simulation, short-lived small-scale Rossby Wave Instability vortices appear, consistent with previous results on VSI-Article number, page 9 of 17 A&A proofs: manuscript no. aanda 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.

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 −5 to 10 −3 . 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 2020) 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.

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 Figure 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 −1 ) 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 20h 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 km s −1 and∼ 0.03 km s −1 , for ALMA Bands 6 and 7, 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 −1 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 12 CO and 13 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 7 , and ConeRot 8 (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 9 (Teague 2019b) can be applied. However, as discussed in Section 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.

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 12 CO, 13 CO and C 18 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 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 non-thermal 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 −1 ) 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.  Results of the expected observable velocity perturbations from a VSI unstable disk 12 CO synthetic lines observations (v 0 −v eq ), 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 (δv 0 ). The FWHM of the circular Gaussian beams are shown at the bottom left corner of each panel. The x-and y-axes indicate the R.A. and Dec. angular offset from the position of the central star, in arcseconds.   Radial profile of the disk layers where the disk becomes optically thick (τ = 1 surface), for different CO isotopologues of the VSIunstable disk radiative transfer model. We illustrate the τ = 1 surface for the J:2-1 and J:3-2 transitions of three CO isotopologues, 12 CO, 13 CO and C 18 O. We assume a face-on disk orientation for the calculation. The grey 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.45R.