# Mapping the Complex Kinematic Substructure in the TW Hya Disk

Richard Teague, Jaehan Bae, Sean M. Andrews, Myriam Benisty, Edwin A. Bergin, Stefano Facchini, Jane Huang, Cristiano Longarini, David Wilner

► To cite this version:

Richard Teague, Jaehan Bae, Sean M. Andrews, Myriam Benisty, Edwin A. Bergin, et al.. Mapping the Complex Kinematic Substructure in the TW Hya Disk. The Astrophysical Journal, 2022, 936, 10.3847/1538-4357/ac88ca. insu-03860299

HAL Id: insu-03860299

https://insu.hal.science/insu-03860299

Submitted on 19 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Mapping the Complex Kinematic Substructure in the TW Hya Disk

Richard Teague¹,²,*, Jaehan Bae³,†, Sean M. Andrews²,‡, Myriam Benisty⁴,†, Edwin A. Bergin⁵,†, Stefano Facchini⁶,†, Jane Huang⁵,<fnref n="7" />,†, Cristiano Longarini⁶,†, and David Wilner²,†

¹Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA; rteague@mit.edu

²Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

³Department of Astronomy, University of Florida, Gainesville, FL 32611, USA

⁴Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁵Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, USA

⁶Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16, I-20133 Milano, Italy

Received 2022 June 8; revised 2022 July 29; accepted 2022 August 9; published 2022 September 12

## Abstract

We present ALMA observations of CO $J=2-1$ and CS $J=5-4$ emission from the disk around TW Hydrae. Both molecules trace a predominantly Keplerian velocity structure, although a slowing of the rotation velocity is detected at the outer edge of the disk beyond ≈140 au in CO emission. This was attributed to the enhanced pressure support from the gas density taper near the outer edge of the disk. Subtraction of an azimuthally symmetric background velocity structure reveals localized deviations in the gas kinematics traced by each of the molecules. Both CO and CS exhibit a “Doppler flip” feature, centered nearly along the minor axis of the disk (PA ~ 60°) at a radius of 1.′′35, coinciding with the large gap observed in scattered light and mm continuum. In addition, the CO emission, both through changes in intensity and its kinematics, traces a tightly wound spiral, previously seen with higher-frequency CO $J=3-2$ observations. Through comparison with linear models of the spiral wakes generated by embedded planets, we interpret these features in the context of interactions with a Saturn-mass planet within the gap at a position angle of PA = 60°, consistent with the theoretical predictions of Mentiplay et al. The lack of a corresponding spiral in the CS emission is attributed to the strong vertical dependence on the buoyancy spirals, which are believed to only grow in the atmospheric of the disk, rather than those traced by CS emission.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Exoplanet formation (492); Radio interferometry (1346)*

## 1. Introduction

There is mounting evidence that protoplanetary disks exhibit a similar complexity of substructures in their gas density and velocity structures as they do in their dust distributions (Disk Dynamics Collaboration et al. 2020; Pinte et al. 2022). Identifying and characterizing these structures promises to yield unique insights into the planet formation process, as young, recently formed planets can be detected through the subtle dynamical perturbations they impart on the gas through which they move (Perez et al. 2015; Pinte et al. 2018, 2019; Casassus & Pérez 2019). Such detections provide opportunities to study planet formation in action, and facilitate the detection of planets that evade detection at optical and/or infrared wavelengths.

Beyond the detection of planets, kinematic-focused observations help localize variations in the rotation speed of the gas, facilitating the inference of the radial gas pressure profile (Teague et al. 2018a, 2018b; Yu et al. 2021) and enabling direct tests of grain trapping in pressure maxima (Rosotti et al. 2020; Yen & Gu 2020; Boehler et al. 2021). With the Atacama Large Millimeter/submillimeter Array (ALMA) delivering unprecedented spatial resolutions and sensitivities, these analyses have

been extended to also probe vertical motions (Teague et al. 2019a; Yu et al. 2021), revealing large-scale flows that are likely examples of meridional flows (Morbidelli et al. 2014; Szulágyi et al. 2014). Such structures provide a key dynamical link between the chemically rich atmospheric layers of the disk with the relatively inert midplanes, potentially transporting volatile rich gas to accreting planets (Cridland et al. 2020). Mapping these structures is thus of immense importance in understanding the delivery and formation of exoplanetary atmospheres.

It is clear, therefore, that dedicated studies of the gas velocity structure of disks provide new constraints for the key dynamical processes and the presence (or absence) of still-forming planets. Previous studies have focused predominantly on CO emission out of a necessity to achieve a good signal-to-noise ratio across a narrow velocity (or frequency) bandwidth (although there are a few recent exceptions that have been able to use CO isotopologue emission, e.g., Teague et al. 2021; Yu et al. 2021). We present new observations designed specifically for the kinematical analyses of both CO and CS emission from the disk around TW Hya. The sensitivities of the observations were chosen to demonstrate the utility of investing in longer integrations that allow the study of the dynamics of the gas traced by molecules less abundant than CO, which exists closer to the disk midplane.

As one of the closest protoplanetary disks to Earth ($d = 60.1$ pc; Gaia Collaboration et al. 2018), TW Hya has been extensively studied across a range of frequencies. At submm wavelengths, the continuum emission has been mapped at an exquisite ~20 mas resolution (1.2 au; Andrews et al. 2016; Macías et al. 2021), revealing multiple concentric rings. A more recent study,

<figure><img src="image_3.png" /></figure>

<footnote n="7">⁷NASA Hubble Fellowship Program Sagan Fellow.</footnote>
achieving an unprecedented sensitivity, reports the detection of mm continuum emission to larger radii of ~100 au (Ilee et al. 2022). An intriguing, slightly elongated feature was reported by Tsukagoshi et al. (2019), which Zhu et al. (2022) argue could be the spatially resolved envelope of a 10–20 $M_\odot$ planet. In terms of the gas, CO emission was originally reported by Zuckerman et al. (1995), with Qi et al. (2004) resolving the characteristic emission morphology indicative of Keplerian rotation. At shorter wavelengths, the sub-µm dust in the disk atmosphere has been mapped through scattered light observations from both space- (Roberge et al. 2005; Debes et al. 2013, 2016) and ground-based facilities (Akiyama et al. 2015; Rapson et al. 2015; van Boekel et al. 2017). Gaps coincident with those detected in the mm continuum were identified, in addition to a large, broad gap at ~90 au, likely associated with a depletion in CS emission reported by Teague et al. (2017), and aligning with the recently reported mm continuum gap (Ilee et al. 2022). A kinematic analysis of the high spatial resolution, but low spectral resolution, images of CO emission presented in Huang et al. (2018) found a tightly wound spiral feature in the gas velocity structure on top of the background rotation (Teague et al. 2019b). Such spiral features were proposed by Bae et al. (2021) to be buoyancy resonances driven by an embedded planet within the outer gap.

As such, TW Hya represents an ideal candidate for a thorough kinematic analysis to fully characterize the possibilities of ongoing planet–disk interactions. In Section 2, the observations and the subsequent calibration and imaging are described, including the generation of moment maps. A study of the velocity structure is performed in Section 3, looking at bulk background motion, localized deviations, and deprojecting disk-frame velocities. The interpretations of these velocity structures are discussed in Section 4, with the main findings summarized in Section 5.

## 2. Observations

Observations were acquired as part of project 2018.A.00021. S (PI: Teague). Teague et al. (2022) described the observational setup and calibration of these data in detail. We also combined this with data from project 2013.1.00387.S (PI: Guilloteau), originally presented in Teague et al. (2016), which has an almost identical correlator setup. In this paper, we focus on the CO and CS emission.

### 2.1. Post Processing

We start with the calibrated measurement sets from Teague et al. (2016; “2013 data”) and Teague et al. (2022; “2018” data). The continuum fluxes were compared using the tools<fnref n="8" /> provided by the DSHARP Large Program (Andrews et al. 2018). The two continuum fluxes were found to be within ≲2% of one another, so no flux rescaling was applied. As the 2013 data were already aligned, the phase center was updated to match that of the 2018 data using the `fixplanets` command in CASA to account for any proper motion of the source.

### 2.2. Imaging

As in Teague et al. (2022), Briggs weighting with a robust value of 0.5 and a channel width of 40 m s⁻¹ was adopted for the imaging of both CO and CS emission. This resulted in a

synthesized beam of 0.′′19 × 0.′′17 (89°) and 0.′′18 × 0.′′17 (87°) for CO and CS, respectively. A Keplerian mask was generated for both molecules for the CLEANing,<fnref n="9" /> making sure that all disk emission was contained within the mask. Following Teague et al. (2018b), a stellar mass of $M_\star = 0.6 M_\odot$ and a viewing geometry described by $i = 5°.8$ and $PA = 151°$ were adopted for this mask. After the imaging, a correction was applied to the image to account for the non-Gaussian synthesized beams, due to the combination of several different array configurations (see the discussion in Czekala et al. 2021; Casassus & Cár-camo 2022), as proposed by Jorsater & van Moorsel (1995). The resulting rms in a line-free channel was measured to be 0.6 mJy beam⁻¹ and 0.7 mJy beam⁻¹ for CO and CS, respectively. Integrating over the Keplerian mask for CO and CS, we recover integrated intensities of $18.31 \pm 0.03$ Jy km s⁻¹ and $1.69 \pm 0.03$ Jy km s⁻¹. These uncertainties do not include the systematic uncertainty of ~10% associated with the flux calibration of the data.

### 2.3. Channel Maps

Channel maps from the images are presented in Figures 1 and 2 for CO and CS, respectively. Only the central channels where there is extended emission are shown, although emission close to the disk center is detected out to velocities of ±~3 km s⁻¹ relative to the systemic velocity of 2.84 km s⁻¹ for CO (comparable to that reported by Rosenfeld et al. 2012) and ±~1 km s⁻¹ for CS. Both lines exhibit a high level of substructure that is seen across multiple channels. The CO emission morphology is dominated by large arcs such that, in the central channels (those close to the systemic velocity of 2.84 km s⁻¹), emission can be seen extending around the full azimuth of the disk. In these central channels, two rings are seen bridging the two lobes of emission across the major axis of the disk, suggesting regions of higher density. These arcs can also be seen in channels offset from the line center by ~100 m s⁻¹ as small protrusions from the bulk emission.

The CS emission shows a similar level of substructure, but with a different morphology, likely because the emission is optically thin and is thus more sensitive to changes in the background gas density or abundance structure than CO emission. The inner regions are dominated by a bright ring, out to a radius of ~1.′′5, with a more diffuse outer disk. Two gaps are observed at ~1.′′5 (originally reported in Teague et al. 2017) and ~2.′′6 (156 au; previously reported by Teague et al. 2018c), and can be traced in many channels around the systemic velocity.

### 2.4. Moment Maps

To aid in the analysis of the data, we collapse all data along the spectral axis into 2D maps using the methods described below. All collapsing was done with the `bettermoments`<fnref n="10" /> Python package.

#### 2.4.1. Integrated Intensity

For the total integrated intensity, or zeroth-moment maps, we integrated the emission along the spectral axis using the Keplerian CLEAN mask to circumvent any sigma clipping. The CO zeroth-moment map has already been published in Teague

<footnote n="8">⁸ https://bulk.cv.nrao.edu/almadata/lp/DSHARP/scripts/reduction utilis.py</footnote>

<footnote n="9">⁹ github.com/richteague/keplerian_mask</footnote>

<footnote n="10">¹⁰ github.com/richteague/bettermoments</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1. CO (2–1) channel maps, imaged at 40 m s⁻¹ spacing. The beam size is shown in the bottom left of the panels. The channel velocity is shown in the top left of each panel. Each panel spans ±4''5 (±270 au) across. Note that the color scaling has been adjusted to bring out the detail in the line wings, and the peak intensity is 176.4 mJy beam⁻¹.</figcaption></figure>

et al. (2022), where a shadow was detected in the outer disk. These moment maps are shown in Figure 3 and compared to the 233 GHz continuum originally presented in Macías et al. (2021) to highlight the difference in the radial extent of the molecular gas and the mm-sized grains. We use the Python package GoFish (Teague 2019) to produce azimuthally averaged radial profiles of these zero moment maps. Emission from CO is detected out to 4''2 (252 au) and CS is detected out to 3''9 (235 au), while the continuum edge is around 1''65 (100 au; Ilee et al. 2022). In these radial profiles, the substructures identified in the CO and CS channel maps are more apparent, most notably the inner cavity and two gaps in the CS emission at ≈1''5 (90 au) and ≈2''6 (156 au).

### 2.4.2. Analytical Spectral Profiles

To calculate maps of line profile properties, such as the line center, width, and peak intensity, fitting of analytical profiles with bettermoments was found to yield the best results,

owing to the high spectral resolution of the data. This is in contrast to the results presented in Teague et al. (2019b), where a quadratic curve was fit to the velocity pixel of peak intensity and its two adjacent velocity pixels (see, e.g., Teague & Foreman-Mackey 2018), as the data were taken at a much lower spectral resolution. At higher spectral resolutions the quadratic method is limited by the low level of curvature at the line center for a well-resolved Gaussian profile. This was particularly an issue for the CO emission, which was found to have a highly saturated core, due to the large optical depths (as discussed in Teague et al. 2016).

While the CS was well-reproduced by a Gaussian profile,

$$I(v) = I_0 \exp \left( - \left( \frac{v - v_0}{\Delta V} \right)^2 \right), \quad (1)$$

with $v_0$ describing the line center, $\Delta V$ is the Doppler width of the line (where FWHM = $2\sqrt{\ln 2}\Delta V$) and $I_0$ is the peak
<figure><img src="image_5.png" /><figcaption>Figure 2. CS (5–4) channel maps, imaged at 40 m s⁻¹ spacing. The beam size is shown in the bottom left of the panels, and the channel velocity is shown in the top left of each panel. Each panel spans ±4''5 (±270 au) across. Note that the color scaling has been adjusted to bring out the detail in the line wings. The peak intensity is 26.4 mJy beam⁻¹.</figcaption></figure>

intensity, the CO emission profile was modeled as

$$I(v) = I_0 \cdot \frac{1 - \exp(-\tau(v))}{1 - \exp(-\tau_0)}, \quad (2)$$

with the optical depth $\tau(v)$ varying as a Gaussian with a width of $\Delta V$ and a peak optical depth of $\tau_0$. Note that $\tau_0$ should not be taken as an absolute measure of the column density, but rather a measure of how saturated is the core of the line, which is related to the optical depth of the line, with larger values relating to more saturated cores. It was tested whether a Gauss–Hermite quadrature expansion, as is often used for studies of galactic dynamics (e.g., van der Marel & Franx 1993), provided a better fit than an “optically thick Gaussian” to the CO emission. However, the Gauss–Hermite expansion failed to describe the data as well as the “optically thick Gaussian,” in addition to requiring an additional free parameter.

A straightforward minimization of $\chi^2$ was achieved with the `curve_fit` routine in `scipy` (Virtanen et al. 2020) to yield an

initial fit to the data. Although fast, this approach often found local minima in the $\chi^2$ surface resulting in moment maps with a large variance over small spatial ranges. Instead, an MCMC sampling approach, implemented with the `emcee` package (Foreman-Mackey et al. 2013), was found to be more forgiving in such cases. Taking the median of the posterior distributions as the “best-fit” parameter produced far more spatially coherent maps. For this fitting, the Keplerian `CLEAN` masks were adopted and only pixels with at least as many unmasked voxels as there were free parameters in the model (four for the CO emission, three for the CS) were included in the fit. It was tested whether smoothing with a simple top-hat kernel with a width of four channels (120 m s⁻¹) would improve the results, but little difference was found between maps made with or without the smoothing. For this manuscript, we focus primarily on $v_0$, shown in Figure 4, while the moment maps for each additional parameter can be found in Appendix A. An exploration of how accurately the parameters are recovered in the presence of noise can be found in Appendix B.
<figure><img src="image_6.png" /><figcaption>Figure 3. Comparison of the archival 233 GHz continuum from Macías et al. (2021) and the velocity-integrated intensity (zeroth-moment) maps of CO and CS, shown in panels (a), (b), and (c), respectively. Note that the CO data has a “square-root” stretch applied in order to bring out the morphology. In each panel, the synthesized beam is shown in the lower left and a 100 au scale bar in the lower right. The right-most panel, (d), shows the azimuthally averaged radial profiles of the CO and CS emission. The shaded regions show uncertainties in steps of 1σ.</figcaption></figure>

### 3. Mapping the Velocity Structure

In this section, we explore how we can decompose the observed line-of-sight velocity into its constituent parts in order to map the velocity structure of the disk and isolate localized velocity perturbations.

#### 3.1. Background Motions

It is possible to model the background rotation and subtract this component to reveal any non-Keplerian motions within the disk (e.g., Casassus & Pérez 2019; Teague et al. 2019b; Casassus et al. 2021; Wölfer et al. 2021). There are multiple approaches that can be employed for such a goal, each with their own benefits and limitations. One main consideration is whether to try to model the disk as a whole, thereby imposing some functional form to $v_φ$, such as a purely Keplerian rotation profile (e.g., Teague et al. 2019b, 2021; Rosotti et al. 2020; Wölfer et al. 2021; Izquierdo et al. 2021, 2022), or to instead split the disk into concentric annuli and fit each annulus independently, allowing for a far more flexible profile (e.g., Teague et al. 2018a, 2018b, 2019a; Casassus & Pérez 2019; Casassus et al. 2021).

One consideration is that, with improvements in the precision at which rotation curves can be measured, systematic deviations associated with the radial pressure profile (Dullemond et al. 2020; Teague et al. 2021) or disk self-gravity (Veronesi et al. 2021) are becoming apparent, resulting in large residuals at the outer edge of the disk. This issue can be circumvented with the independent annuli approach. However, as a model comprised of independent annuli is highly flexible, it can run the risk of overfitting the data and confusing the disentangling of various velocity components, in particular if the background disk is not azimuthally symmetric (see the discussion regarding the choice of a background velocity profile in Teague et al. 2019a). As there is substantial structure in the TW Hya disk that is not azimuthally symmetric, we focus here on just the analytical models for $v_φ$.

##### 3.1.1. Analytical Models

To perform the fitting of the background rotation profile to the $v_0$ maps, we use the eddy Python package (Teague 2019a)

<figure><img src="image_7.png" /><figcaption>Figure 4. Line-of-sight velocity maps, $v_0$, for CO (left) and CS (right). The synthesized beam is shown in the bottom left of each panel. The $v_0$ maps are masked where the peak intensities fall below a 5σ level for CO and a 3σ level for CS.</figcaption></figure>

and fit a Keplerian rotation profile of the form

$$v_{\text{kep}} = \sqrt{\frac{GM_{\star}r^2}{(r^2 + z^2)^{\frac{3}{2}}},} \quad (3)$$

where $M_\star$ is the dynamical stellar mass and $(r, z)$ are the disk-frame cylindrical radius and height, related to the on-sky pixel through the disk inclination, $i$, position angle, PA, and a source center offset, $(x_0, y_0)$. Thus, for this model, a total of five free parameters are needed: $\Theta_{\text{kep}} = \{x_0, y_0, \text{PA}, M_\star, v_{\text{LSR}}\}$. For face-on disks, there is an extreme degeneracy between $M_\star$ and $i$, so we fix $i = 5^{\circ}8$, following Teague et al. (2019b). We note that this is different from other recent works studying TW Hya, such as Huang et al. (2018) and Macías et al. (2021), which adopt a lower $i = 5^{\circ}0$; however, we adopt the higher inclination for a more direct comparison of the kinematical features described in Teague et al. (2019b). We have verified that, at such low inclinations, the structures and features described in this manuscript are all robust to the choice of inclination.

As previously mentioned, there is mounting evidence that many disks exhibit non-Keplerian rotation due to the contributions of disk self-gravity (Veronesi et al. 2021) or the radial pressure gradient (Dullemond et al. 2020; Teague et al. 2021). To mimic the effect of these on the rotation curve, we have
included a pseudo-disk-mass term, $M_{\text{disk}}$, where positive values hasten the rotation speed and negative values slow the rotation. For the case of a sub-Keplerian outer disk, three free parameters are included to the modeling: $\Theta_{\text{prs}} = \{\Theta_{\text{kep}}, M_{\text{disk}}, \gamma, r_{\text{in}}\}$. $r_{\text{in}}$ describes the cylindrical radius where $v_{\phi}$ begins to deviate from Keplerian rotation, while $M_{\text{disk}}$ dictates the magnitude of this deviation at the outer edge of the disk (for this fitting, we set $r_{\text{out}} = 4''$). In essence, $v_{\phi}$ is a Keplerian profile around a star with a dynamical mass of $M_{*}$ inside $r_{\text{in}}$, and approaches a Keplerian profile around a star with a dynamical mass of $M_{*} + M_{\text{disk}}$ (remembering that $M_{\text{disk}}$ is a negative value in this case) at $r_{\text{out}}$. How quickly the effective dynamical mass of the system changes with radius is governed by $\gamma$, with $\gamma = 0$ resulting in a linear progression and $\gamma > 0$ describing a quick initial change that then slowly approaches the final value. The full parameterization used in the eddy fitting is described in more detail in Appendix C. Plots of how the $v_{\phi}$ profile changes under various choices of $\Theta_{\text{prs}}$ are also shown.

For each map, we consider two models: pure Keplerian rotation, $v_{\text{kep}}$, requiring five free parameters, $\Theta_{\text{kep}}$; and pressure support rotation, $v_{\text{prs}}$, requiring eight free parameters, $\Theta_{\text{prs}}$. We model the posterior distributions through eddy, which wraps the emcee MCMC sampling package (Foreman-Mackey et al. 2013). Only regions between 0"/4 (roughly twice the beam FWHM) and out to 3"/5 for CO and 3"/5 for CS, were fit, where the inner boundary was to chosen to avoid effects from beam smearing, while the outer boundary removed regions with low-S/N measurements of $v_0$. The cubes were spatially downsampled such that only spatially independent pixels were considered. One hundred and twenty-eight walkers were initialized in a cluster around previous literature values; they were then given 5000 steps to burn in and an additional 5000 steps to sample the posterior distribution. We find consistent posterior distributions for all parameters between both lines and for those from previous works (Huang et al. 2018; Teague et al. 2019a). Neglecting the nuisance source-center offset parameters, the median values for the Keplerian model were found to be PA = 151°/6, $M_* = 0.82 M_\odot$, and $v_{\text{LSR}} = 2843 \text{ m s}^{-1}$. The width of the posteriors (typically used as a measure of the statistical uncertainty) were exceptionally narrow, resulting in fractional uncertainties of $\lesssim 0.1\%$ for all parameters. Such low statistical uncertainties are likely a symptom of a relatively inflexible model rather than a realistic constraint on those parameters, and as such they are not quoted, in order to avoid misrepresentation.

For the CO, the outer $\sim 1''$ was found to be overpredicted with a Keplerian model suggesting a sub-Keplerian outer disk. With the pressure support rotation model, Gaussian posteriors were found for the three additional parameters, with $M_{\text{disk}} = -0.12$, $\gamma = 1.44$, and $r_{\text{in}} = 1''/71$, again with the widths of the posteriors being below 1% of these quoted values. The five other parameters were consistent with the values found for the purely Keplerian fit, including the dynamical mass of the star. No similar constraints could be made for the $v_0$ map calculated from CS emission, likely due to the low S/N of the emission in the regions where the deviation is strongest ($r \gtrsim 3''$).

## 3.2. Non-Keplerian Motions

Figure 5 shows the residuals after subtracting the background velocity structures inferred in the previous subsection.

Panels (a) and (b) compare the CO velocity residuals after subtracting a purely Keplerian profile, $v_{\text{kep}}$, and a profile including a sub-Keplerian outer disk, $v_{\text{prs}}$, respectively. The sub-Keplerian outer edge is clear in panel (a), with negative and positive residuals along the redshifted and blueshifted major axis of the disk, which are mostly removed when allowing for the more flexible $v_{\phi}$ profile. This asymmetry is removed with the pressure-supported model in (b), with little structure seen in the outer regions of the disk. No difference can be seen between the CS residuals when changing between $v_{\text{kep}}$ and $v_{\text{prs}}$ as the background model. This is as expected because the model parameters for the pressure support were unable to converge.

The most striking feature is a large, tightly wound spiral positive residual, starting at roughly (1"/5, 0"/0), and extending out to the disk edge after almost a full $2\pi$ winding. The start of this spiral was previously detected by Teague et al. (2019b) using CO (3-2) emission, but due to the lower spectral resolution of that data, it was only considered an arc, rather than a full spiral. At smaller radii, two arcs of positive residuals are seen: one directly south of the disk center at an offset of $\sim 0''/6$, and a second in a southwesterly offset from the star at an offset of $\sim 1''$. Both these features were seen in the CO (3-2) data, suggesting that these are real features and that these observations are hitting the kinematical equivalent of the “confusion limit.”

The velocity structure traced by CS emission shows far less structure, due to the lower S/N of the emission. However, a large positive feature is detected coincident with the peak residual near the start of the spiral arm detected in CO. A weaker negative residual is seen ahead (the disk is rotating in a clockwise direction), up to a PA $\sim 0°$, before turning into a positive residual again. Such a feature is highly reminiscent of the “Doppler flip” found in HD 100546 (Casassus & Pérez 2019), proposed to be due to the spiral shocks from an embedded planet (e.g., Perez et al. 2015). This feature is marked in Figure 5(d) with a dashed contour. Guided by the presence of this feature in the CS map, a similar feature is potentially traced by the CO emission at the same location.

No obvious kinematic deviations are observed around the continuum feature reported by Tsukagoshi et al. (2019) close to the southeastern minor axis of the disk (PA $\sim 237°$) at a distance of 52 au (0"/87).

## 3.3. Deprojecting Velocity Components

The measured line-of-sight velocity shown in Figure 4 is the sum of the projection of all local velocity components,

$$v_0 = v_{\text{LSR}} + v_{\phi,\text{proj}} + v_{r,\text{proj}} + v_{z,\text{proj}}, \quad (4)$$

where the projection of each component along the line of sight is given by

$$v_{\phi,\text{proj}} = +v_{\phi} \cos(\phi) \sin(|i|), \quad (5a)$$

$$v_{r,\text{proj}} = -v_{r} \sin(\phi) \sin(i), \quad (5b)$$

$$v_{z,\text{proj}} = -v_{z} \cos(i), \quad (5c)$$

where $\phi$ is the polar angle in the disk, defined to be $\phi = 0°$ along the redshifted major axis, and increasing in a counter-clockwise direction (i.e., increasing with PA). As the projection of each velocity component has a different azimuthal dependence, it is possible to leverage this difference to disentangle the velocity residuals. For example, if a feature
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The $v_0$ residuals for CO, top, and CS, bottom, after subtracting a purely Keplerian background model, $v_{\text{kep}}$, left column, or a pressure-supported Keplerian model, $v_{\text{prs}}$, right column. In each panel, the synthesized beam is shown in the lower left, while a 100 au scale bar is shown in the lower right. The dotted contours show radial steps of 0.5 and azimuthal steps of 30°. The colorbar shows the residuals in units of $\sigma$, with absolute residuals less than $1\sigma$ left white. The dashed contour in (d) highlights the Doppler flip.</figcaption></figure>

with strong $v_\phi$ components crosses the minor axis of the disk, $\phi = \pm 90^\circ$, then the projected velocity $v_{\phi,\text{proj}}$, which is observed, will flip sign. Conversely, a $v_r$ component that crosses the minor axis of the disk will not flip sign; the observed $v_{r,\text{proj}}$ will remain constant. Therefore, under the assumption that the velocity residuals are driven by a single velocity component, we can deproject the residuals in Figure 5 to recover the true, intrinsic velocity structure. If these are found to exhibit multiple changes in sign coincident with disk axes, an unlikely physical situation, then a strong argument can be made against that velocity component being responsible for the observed residuals.

To apply this technique, we deprojected the residuals after subtracting the best pressure-supported Keplerian rotation model, assuming that the perturbations were fully azimuthal, radial, or vertical in nature through Equations (5a), (5b), and

(5c), respectively. Figure 6 shows the results for CO in the top row, and CS in the bottom row. For CO, it is clear that a majority of the structure cannot be azimuthal or radial in nature, due to the numerous changes in sign observed across both axes. Indeed, the azimuthal coherence of the spiral-shaped residuals point to a vertical component, with gas moving from the atmosphere toward the disk midplane, as argued for in Teague et al. (2019b).

Around the northeast disk minor axis (PA $\sim$ 60°), the interpretation is less straightforward. Both CO and CS show residuals that, under the assumption of rotational velocity components, yield the most azimuthally coherent features, as shown in Figures 6(b) and (d). A region of hastened rotation is expected at the outer edge of a gap, due to the large positive radial pressure gradient, as would be expected based on the substantial surface density perturbation at that location (Teague et al. 2017; van Boekel et al. 2017). However, such a feature
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The deprojected velocity residuals for CO (top row) and CS (bottom row), assuming all velocities are azimuthal (left column), radial (central column), or vertical (right column). For the azimuthal and radial velocities, wedges along the minor and major axes have been masked, as here the observations are insensitive to those components. In each panel, the synthesized beam and a 50 au scale bar are shown in the bottom left and right corners. Note that the field of view for the bottom row is smaller than the top row.</figcaption></figure>

would be expected to be azimuthally symmetric and not limited to just the northeast half of the disk. In this scenario, one could argue that the large-scale spiral traced by CO may hide the super-Keplerian rotating gas along the southeast side of the disk; however, the lack of a similar spiral feature traced by the CS emission would therefore require a different explanation for the deeper regions of the disk.

An alternative interpretation is to consider the case of when the underlying velocity structure is not constant in direction as a function of azimuth. For example, if there is in fact a “Doppler flip,” a feature predicted to be indicative of an embedded planet (Casassus & Pérez 2019), then a change in velocity direction would be expected at the location of the planet. Calcino et al. (2022), following the analytical prescription of spiral wakes from Rafikov (2002), demonstrate that the

radial velocity components will be the dominant velocity component, except in the very immediate vicinity of the planet, where a radially inward (negative $v_r$) component should be leading the planet and a radially outward (positive $v_r$) component should be trailing the planet. This morphology is indeed what is seen for both CO and CS (Figure 6(b) and (e)), with the “Doppler flip” center aligning almost exactly with the disk minor axis. In this scenario, an embedded planet at a radial location of ~1.36 (~~82 au) would naturally explain the presence of the surface density depletion at ~85 au, as well as providing a mechanisms for launching the large spiral structure that appears to emanate from a similar location. Figure 7 demonstrates this proposed scenario, which is discussed more in Section 4.1.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Polar projections of the velocity residual maps from Figure 5 with CO after subtracting a purely Keplerian background profile in panel (a), after subtracting a pressure-supported background in panel (b), and CS after subtracting a Keplerian background profile in panel (c). In all panels, the residuals are normalized by the local uncertainty on the line center, $\delta v_0$. The two horizontal lines in each show the gap edges of the continuum gap at 82 au ($1''/35$) reported by Ilee et al. (2022), coincident with the gap in scattered light (Debes et al. 2013; Rapson et al. 2015; van Boekel et al. 2017) and molecular line emission (Teague et al. 2017). The dotted boxes highlight the “Doppler flip” feature. The spirals traced by the CO are annotated. The circle marks the center of the “Doppler flip.” The dotted square marks the location of the continuum excess reported in Tsukagoshi et al. (2019). The disk rotates from right to left.</figcaption></figure>

# 4. Discussion

We have shown that the velocity structure of the protoplanetary disk around TW Hya, traced by CO and CS emission, exhibits substantial perturbations in the velocity structure. What is particularly interesting about this set of observations is that these molecules are believed to trace two distinct vertical regions in the disk: CO traces the disk atmosphere while CS traced closer to the disk midplane (e.g., Dutrey et al. 2017). This combination of lines therefore offers a

glimpse of the vertical dependence of kinematic features within the disk. In this section, we discuss the implication of these findings.

## 4.1. An Embedded Planet

The search for planets in TW Hya has had a long history. Multiple works have modeled various properties of the outer gap probed by different observations, be it total scattered light (Debes et al. 2013), polarized scattered light (van Boekel et al. 2017; Mentiplay et al. 2019), dust continuum (Ilee et al. 2022), or molecular line emission (Teague et al. 2017), to infer the mass of the putative planet responsible for opening the gap. These works have reported planet masses spanning between 10–100 $M_\oplus$, all well below the current limits of direct detection methods, which are sensitive down to a mass of ~2 $M_{Jup}$ (Asensio-Torres et al. 2021) at 90 au. The detection of localized velocity perturbations coincident with this gap therefore offers a new quantitative probe of the putative planet’s mass.

### 4.1.1. “Doppler Flip”

The most intriguing aspect of the velocity residuals is the feature that bears a striking similarity to the “Doppler flip” morphology discussed in Casassus & Pérez (2019), where it was argued this it would be due to an embedded planet. This interpretation is particularly enticing for the case of TW Hya, due to coincidence of the “Doppler flip” and the outer gap location at ~82 au, as shown in Figure 7. Such a coincidence between a localized velocity perturbation and a continuum gap has been previously argued to be a strong indication of an embedded planet (e.g., Pinte et al. 2019). Under the assumption that the motions associated with the “Doppler flip” are dominated by in-plane velocity components (i.e., either azimuthal or radial in direction) then Figures 6(b) and (e) show that CO and CS emission are tracing perturbations on the order of 50 m s⁻¹ and 100 m s⁻¹, respectively. A stronger perturbation closer to the midplane is exactly what is predicted by numerical simulations, as the spirals arising due to Lindblad resonances grow weaker with height in the disk (e.g., Pinte et al. 2019). Indeed, the variation in amplitudes of the perturbations found in the simulations of Pinte et al. (2019) are consistent with the difference in vertical heights expected to be traced by CO and CS emission.

While a full hydrodynamical modeling of these observations is left for future work, the observed velocity perturbations can be compared with predictions from linear models (e.g., Bollati et al. 2021) to estimate the putative planet’s mass.<fnref n="11" /> As these linear models only consider a flat, 2D disk, we limited this analysis to the CS emission, for which the assumption of a 2D system is a more appropriate approximation than for the highly elevated region traced by CO emission. We assumed that the temperature and surface density of the disk over the region of interest are well-described by single power laws with exponents $p = -0.75$ and $q = -0.47$, respectively (e.g., Calahan et al. 2021; Huang et al. 2018), which resulted in a local pressure scale height of $h = 0.073r$ at the radius of the planet when adopting the midplane temperature profile from Huang et al. (2018).

Figure 8 shows the $v_r$ components driven by embedded planets with respective masses of 0.05 $M_{Jup}$ (1 $M_{\text{Nep}}$), 0.3 $M_{Jup}$

<footnote n="11">¹¹ https://github.com/DanieleFasano/Analytical_Kinks</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Radial velocity perturbations driven by planets of different masses: 1 $M_{\text{Nep}}$, 1 $M_{\text{Sat}}$, and 0.5 $M_{\text{Jup}}$, shown in panels (a), (b), and (c), respectively. Each panel is centered on the putative planet's location, with the dotted line showing the orbital radius of the planet, 82 au. The solid contours mark 100 m s⁻¹, the velocity traced by CS emission. The cross marks the center of the disk. Each panel spans a field of view of 2′′5 (150 au). The analytical models were calculated following Bollati et al. (2021).</figcaption></figure>

(1 $M_{\text{Sat}}$), and 0.5 $M_{\text{Jup}}$, from left to right, convolved with a Gaussian kernel with an FWHM of 0.18 to mimic the observations. Although the morphologies of the perturbations do not perfectly match the deprojected observations in Figure 6, the magnitudes of the perturbations, ~100 m s⁻¹, are best matched by the 0.3 $M_{\text{Jup}}$ planet, in line with the inferred planet mass based on the size of the gap. Perturbations in the azimuthal direction, $v_\phi$, are not shown, as their morphology is inconsistent with the observed structures. However, we do note that, as the $v_\phi$ motions for each planet mass are comparable in magnitude to the $v_r$ perturbations, regardless of the velocity component these observations are tracing, the magnitudes of the perturbations are consistent with a Saturn-mass planet.

### 4.1.2. Spiral

A Saturn-mass planet is also consistent with the sub-Jovian mass planet put forward in Bae et al. (2021) to account for the large-scale velocity spiral. As argued in Section 3.3, the azimuthal coherence of the spiral structure suggests that this must be dominated by vertical motions moving from the disk atmosphere toward the disk midplane. Bae et al. (2021) found that, at large azimuthal distances from the planet, Lindblad spirals that drive the “Doppler flip” described above have negligible vertical motions, while buoyancy spirals were dominated by such vertical motions.

Comparing the velocity residuals shown in Figure 5 to the numerical simulations of Bae et al. (2021), it is clear that buoyancy resonances can qualitatively explain the morphology of the observed residuals, namely the tightly wound spirals with substantial vertical motions. The magnitude of the residuals from the least massive model considered—a 0.5 $M_{\text{Jup}}$ planet—is larger than what is observed in TW Hya (~10 m s⁻¹, if the perturbations are believed to be in the vertical direction). As the perturbation scales linearly with planet mass, the observed kinematic perturbations are therefore consistent with the inference of a Saturn-mass planet based on the “Doppler flip” described above. Furthermore, the strong vertical dependence on the strength of the buoyancy spirals naturally explains the lack of a spiral in the CS data. As buoyancy resonances, the driving mechanism of the buoyancy

spirals, require regions with slow thermal relaxation (i.e., a reduced rate of gas-particle collisions), they are expected to only grow in the atmospheric regions of disks, $z/r \gtrsim 0.1$. These regions would be therefore traced by CO emission, while the CS emission probes deeper in the disk ($z/r \lesssim 0.1$; Dutrey et al. 2017) where the thermal relaxation is sufficiently fast that buoyancy spirals cannot grow.

A spiral in brightness temperature was reported in Teague et al. (2019b) through CO (3–2) emission presented in Huang et al. (2018). A direct comparison between the data presented here and these archival data is hindered by the detection of a large-scale shadow believed to be tracing the shadow observed in scattered light (Debes et al. 2017; Teague et al. 2022). Although Teague et al. (2022) did report the presence of a spiral arm in the outer disk consistent with that from Teague et al. (2019b), a more comprehensive analysis would be limited by the different spectral response functions applied to each data set, due to the differing spectral resolutions used when taking the data. Different spectral response functions will result in the intrinsic line profiles being modified in different ways, leading to systematic differences in line widths or peak intensities that could be mistaken for true physical variations.

### 4.1.3. Line Widths

Although this scenario paints a compelling picture of an embedded planet in TW Hya, having the “Doppler flip” feature aligned with the disk minor axis limits the unambiguous decomposition of the velocity residuals. A search for enhanced line widths due to spatially unresolved motions around the putative embedded planet (either from unresolved circumplanetary motion or planet-driven turbulent motions, e.g., Perez et al. 2015; Dong et al. 2019) would offer a method to potentially disentangle these scenarios. Unfortunately, no evidence for localized enhancements was found in the maps of $\Delta V$ traced by CO or CS, as shown in Figure 9.

A tentative enhancement in the background $\Delta V$ is found around the gap location, highlighted by the vertical dashed line in Figure 9(c). This is counter to the reduced line widths found in both ¹²CO and ¹³CO in the gaps of HD 163296 (Izquierdo et al. 2022), which were explained as due to the drop in local
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Residuals in the local line width, $\Delta V$, when an azimuthally averaged model, $\Delta V_{\text{avg}}$, is subtracted for CO (panel (a)) and CS (panel (b)), normalized by the uncertainty of the linewidth measurement, $\delta\Delta V$. The radial profile used as the background model, $\Delta V_{\text{avg}}$, is shown in panel (c). In both (a) and (b), the center of the proposed “Doppler flip” is shown by the white point, while the ellipse in the bottom left shows the synthesized beam. Dotted contours show steps of 0.5 in radius and 30° in azimuth. In (c), the shaded regions show the 1σ standard deviation in each azimuthal bin. The vertical dotted line shows the location of the gap center as reported by Ilee et al. (2022).</figcaption></figure>

density limiting the pressure broadening of the line. One potential solution to this discrepancy is that, for the face-on orientation of TW Hya, observations are far more sensitive to the vertical velocity dispersions within the gap driven by the embedded planet (e.g., Dong et al. 2019). For a moderately inclined disk like HD 163296, the line of sight will only intersect a narrow region of the gap, tracing a smaller column of material, and thus a smaller velocity dispersion, than a line-of-sight directly into the gap.

A precise measurement of the width of a line requires substantially better data than what is needed for a similar precision on the line center (Teague 2019b). As such, the most direct way to search for these unresolved motions would be to move to higher-frequency observations, for example, those tracing the $J=6-5$ transitions of CO isotopologues, which can achieve far better velocity resolutions, improving on those reported here by almost a factor of 4. However, the larger upper state energies of these transitions may result in emission not being observable to the outer edge of the disk (as was found for $^{13}$CO and $^{18}$O emission in Schwarz et al. 2016).

## 4.2. Emission Heights

A particularly interesting point is that a similar dynamical stellar mass was found for both CO and CS: 0.814 $M_{\odot}$ and 0.816 $M_{\odot}$, respectively. As discussed in Yu et al. (2021), if emission is expected to trace an elevated region parameterized by $\chi = z/r$, and a 2D Keplerian model is used to infer a stellar mass, then the true stellar mass would be underestimated by a factor of $(1+\chi^2)^{3/2}$. If two different molecules are used to measure a dynamical mass, then the ratio of the inferred dynamical masses, $M_{\star,a}$ and $M_{\star,b}$ can be used to relate their two emission surfaces, $\chi_a$ and $\chi_b$, through

$$ \chi_b = \left[ (1 + \chi_a^2) \left( \frac{M_{\star,a}}{M_{\star,b}} \right)^{2/3} - 1 \right]^{1/2} . \quad (6) $$

For the two stellar masses measured for CO and CS, this relationship would suggest CS traces a similar height to that of CO. This is seemingly at odds with what is observed in the Flying Saucer (Dutrey et al. 2017) and found more generally from the MAPS program (Law et al. 2021), where CO is observed to trace atmospheric disk regions, $\chi_{CO} \sim 0.3$, while CS is believed to trace closer to the disk midplane, $\chi_{CS} \lesssim 0.1$. To add to this, Teague et al. (2018c) used multiband observations of CS in TW Hya to constrain the temperature of the gas, finding values close to 40 K in the inner 1" of the disk, and dropping to $\approx 20$ K beyond 2". This is considerably lower than the temperature traced by CO emission by between 10 or 20 K (e.g., Huang et al. 2018; Calahan et al. 2021), suggesting that CO traces a more elevated, and thus hotter, region. Assuming then that CO and CS do trace different vertical heights, e.g., $\chi_{CO} = 0.3$ and $\chi_{CS} = 0.1$, then the two dynamical masses should vary by 10%. It is therefore likely that the large spiral structure observed in CO biases the inferred dynamical mass by roughly a similar level.

## 4.3. CO Optical Depth

The arcs and spurs that are annotated in Figure 1 are found to be coincident in radius with the two rings of elevated optical depth at 1.2" (70 au) and 2.65" (160 au), as demonstrated in Figure 10. This suggests that, at these radii, there is a considerable increase in the local CO column density, enhancing the emission in the line wings, which results in the observed protrusions that are symmetric about the line center.

What is particularly interesting is that the radial profile of $\tau_0$ exhibits the two aforementioned peaks separated by a broad plateau. This appears to bear little resemblance to the radial profile of the CO integrated intensity, shown as the dark blue line in Figure 10, which appears relatively smooth. However,
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparing the azimuthally averaged radial profile of the CO optical depth, τ₀, in black, with the azimuthally averaged radial profiles of the integrated intensities of the CO and CS emission, dark and light blue, respectively. The uncertainty shown on τ₀ is dominated by the statistical uncertainty associated with the fitting of the line profile rather than azimuthal scatter.</figcaption></figure>

a subtle bump is seen at 1.2, as previously reported by Huang et al. (2018), consistent with a local enhancement in the CO column density. Comparison with the radial profile of ¹³CO (Figure 7 of Zhang et al. 2019), which is less optically thick than ¹²CO, does show peaks in integrated intensity at broadly the same radial locations at peaks in τ₀, further evidencing that these peaks are indeed tracing annuli of enhanced CO column density. These observations therefore demonstrate that, with sufficient spectral resolution, the subtle change in emission line morphology—namely the saturation of the core—enables a characterization of column density variations even in highly optically thick lines.

The comparison with the CS integrated intensity, shown in light blue in Figure 10, however, is more complex. The inner peak at 1.2 (70 au) does appear to align with the outer edge of the CS ring discussed in Section 2.4, signaling a localized enhancement in the total gas surface density rather than just CO column density. This is also consistent with the rings observed in CN emission (Teague & Loomis 2020) and DCN emission (Öberg et al. 2021), which align with the edge of the pebble disk. The outer peak in τ₀ at 2.65 (160 au), however, instead aligns with a drop in both CS and CN integrated intensity. Teague et al. (2018b) used a multi-transition excitation analysis to demonstrate that this emission deficit was due to a drop in CS column density rather than a change in temperature. Similarly, Teague & Loomis (2020) showed that no change in the excitation temperature of CN was found, requiring instead a drop in CN column density to account for the dip. The difference in the behavior of the CO emission to CS and CN emission therefore points toward chemical processing that could lead to a preferential increase in CO column densities relative to other species. A more thorough exploration of chemical scenarios that could lead to this difference is left for future work.

## 5. Summary

We have presented new, high velocity resolution observations of CO and CS emission from the disk around TW Hya. We have used these to map and characterize the complex kinematic structure of the disk. Line-of-sight velocity maps were calculated by fitting analytical forms to the spectrum in

each pixel, and the bulk background motion was found to be well-described by a Keplerian profile. Regions beyond 2.1 (126 au) traced by CO emission were found to be rotating at slightly sub-Keplerian speeds, likely due to the influence of a strong radial pressure gradient at the edge of the disk (Dullemond et al. 2020).

After subtracting the best-fit rotation model from the data, complex kinematic substructure was found in the gas velocities traced by both molecules, most notably a large, open spiral traced by CO, and a “Doppler flip” feature traced by both CO and CS. The location of the “Doppler flip,” at a radial offset from the central star of 1.35 (82 au) and position angle PA = 60°, lies along the northeast minor axis of the disk and coincides with a gap in the disk. This coincidence in location is strong evidence of an embedded planet opening the gap (Pinte et al. 2019).

A Saturn-mass planet (0.3 M<sub>Jup</sub>) was found to reproduce the magnitude of the velocity perturbations associated with the “Doppler flip,” while also falling within the range of masses inferred based on the width and depth of the outer gap (e.g., Debes et al. 2013; Teague et al. 2017; van Boekel et al. 2017; Mentiplay et al. 2019; Ilee et al. 2022). Furthermore, the modeling of Bae et al. (2021) showed that a sub-Jovian mass planet is capable of launching buoyancy spirals in the disk atmosphere that exhibit morphologies and magnitudes similar to those traced by the CO emission. The interpretation of this spiral as due to buoyancy resonances, which are found only in the atmospheres of disks, naturally explains why it is traced by CO emission but not CS emission.

Although no elevated velocity dispersions were found around the “Doppler flip” location, a tentative enhancement in the background line width was found at the radius of the gap. Such motions could be interpreted as highly localized vertical motions excited by the embedded planet within the gap (e.g., Dong et al. 2019).

In sum, these results show compelling evidence of an embedded, Saturn-mass planet in the disk around TW Hya, carving out the gap at 82 au (1.35). They demonstrate the power of observations designed to trace the kinematical properties of disks, while also emphasizing the utility of coupling molecules known to trace distinct vertical regions in order to extract critical information about the 3D structure of the features found.

We thank the anonymous referee for a thorough reading of the manuscript and for providing helpful comments. This paper makes use of the following ALMA data: ADS/JAO.ALMA# 2013.1.00387.S and ADS/JAO.ALMA#2018.A.00021.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. These observations are associated with program #16228. This project has received funding
from the European Research Council (ERC) under the
European Union's Horizon 2020 research and innovation
program (grant agreement No. 101002188). Support for J.H.
was provided by NASA through the NASA Hubble Fellowship
grant \#HST-HF2-51460.001-A awarded by the Space Tele-
scope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., for NASA, under
contract NAS5-26555.

Software: emcee (Foreman-Mackey et al. 2013), GoFish (Teague 2019), bettermoments (Teague & Foreman-Mackey 2018), CASA (v5.8.0; McMullin et al. 2007), keplerian_mask.py (Teague 2020).

## Appendix A
### Moment Maps

In this appendix, we present the moment maps made from the “optically thick Gaussian” and Gaussian analytical fits discussed in Section 2.4 for CO and CS, shown in Figures 11 and 12, respectively. The maps were generated with the bettermoments package. We note that, although the optical depth is accounted for with the τ₀ variable, the spectral response function of the ALMA correlator lessens the strength of the optically thick core, due to the convolution of the signal with the spectral response function, resulting in an underestimation of the true optical depth.

<figure><img src="image_14.png" /><figcaption>Figure 11. Moment maps of the CO emission. Only regions where $I_0 > 5\sigma$, where $\sigma = 0.7$ mJy beam$^{-1}$, are shown.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 12. Moment maps of the CS emission. The $v_0$ and $\Delta V$ panels are shown on the same scaling as CO in Figure 11. Only regions where $I_0 > 3\sigma$, where $\sigma = 0.3$ mJy beam$^{-1}$, are shown.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 13. The accuracy to which the parameters of a Gaussian profile can be fit to noisy data for a given sampling rate, R, and signal-to-noise ratio, S/N. Each column represents a different parameter, each row represents a different τ₀: 5, 3, 1, and 0, from top to bottom.</figcaption></figure>

# Appendix B

## Theoretical Limits in Extracting Velocity Perturbations

It is useful to consider what the theoretical limits are in terms of detecting velocity deviations, in order to know what deviations we should be able to detect in a given data set. In practice, this is answering the question: “How well can we measure the line centroid?” In this appendix, we build upon the work of Lenz & Ayres (1992) who were able to define a relationship between the sampling rate, R, defined as the number of samples over the FWHM of the line, and the signal-to-noise ratio of the spectrum peak, S/N, to the accuracy to which parameters describing a Gaussian profile could be constrained. From this relationship, it was possible to estimate what level of velocity deviation one could expect to detect for a data set.

Here, we extend this methodology to an opacity-broadened Gaussian profile, more suited to the optically thick $^{12}$CO emission presented in this paper. The methodology is as follows. First, we assume the underlying true spectrum is described by

$$I(v) = I_0 \cdot \frac{1 - \exp(-\tau(v))}{1 - \exp(-\tau_0)}, \qquad (B1)$$

where $I_0$ and $\tau_0$ are the intensity and optical depth at the line center, and where $\tau(v)$ has a Gaussian form,

$$\tau(v) = \tau_0 \exp\left(-\left[\frac{v - v_0}{\Delta V}\right]^2\right), \qquad (B2)$$

where $v_0$ is the line center and $\Delta V$ is the Doppler width (not the typical standard deviation of a Gaussian, but rather a factor of $\sqrt{2}$ larger, such that FWHM = $2\sqrt{\ln 2} \Delta V$). Then, for a given $\{R, S/N, \tau_0\}$ set, a model spectrum is generated on a v-axis sampled at integer values, with the line center, $v_0$, randomly drawn from a normal distribution. This defines the “true” spectrum. Gaussian noise is then added to the model spectrum, with a standard deviation of $\sigma = I_0/[S/N]$, creating the “observed” spectrum. The assumption of white noise is not perfect, but it is a reasonable approximation for ALMA data that has been imaged with a channel spacing larger than the spectral resolution.

Using the Python package lmfit, a model profile is fit to the “observed” spectrum by minimizing $\chi^2$ using the Levenberg-Marquardt algorithm. This returns optimized model parameters and their associated uncertainties (while the full covariance matrix is returned, we only take the diagonal components and ignore any parameter covariances for this analysis). The difference between the optimized parameter and the true parameter is recorded as a measure of the accuracy of the fit. This process is then repeated N times, with each iteration using a new, random spectrum. It was found that, for $N \gtrsim 10^3$, the distribution of (true-optimized) for a given parameter followed a Gaussian distribution, with a mean equal to the average uncertainty returned from the fitting. To ensure this distribution was well sampled, we adopted $N = 10^4$ for any $\{R, S/N, \tau_0\}$ set.
We then sample the {$R$, S/N, $\tau_0$} parameter space. Four values of $\tau_0$ were chosen: 0 (a pure Gaussian model to compare with Lenz & Ayres 1992), 1, 3, and 5. For $\tau_0 > 5$, the line is almost entirely saturated and very little difference in the resulting spectra can be seen. The sampling rate, $R$, was chosen to span from 2 to 20, covering the range of possible rates based on typical line widths (between 100 m s⁻¹ and 250 m s⁻¹, the widths of thermally broadened lines at temperatures between 20 K and 100 K) and spectral resolutions afforded by ALMA (roughly between 30 m s⁻¹ and 250 m s⁻¹). While many observations result in $R < 2$, there is little utility in fitting an analytical model to data so poorly sampled, and thus these values are ignored. The S/N was chosen to range between 3 and 40, again typical of ALMA observations of molecular line emission from protoplanetary disks.

The standard deviation of the (true-optimized) distribution for each parameter and {$R$, S/N, $\tau_0$} set is shown in Figure 13. Both the line center, $v_0$, first column, and the line width, $\Delta V$, second column, are normalized relative to the sampling rate. As can clearly be seen, and as expected, the accuracy increases with both an increased $R$ and S/N. Adopting a similar functional form to that used in Lenz & Ayres (1992), we can model the relationship between the accuracy to which the line center can be determined and the data characteristics. Again using lmfit, we find $\delta(v_0)/\Delta V_{chan} = 0.69 \sqrt{R} / [S/N]$ when averaged over all four $\tau_0$ values. Applying this relationship to the data presented here, the channel spacing is 40 m s⁻¹. The typical Doppler widths of the lines are 160 m s⁻¹, resulting in $R \approx 6.7$. Assuming that the S/N ~ 10 for a given pixel, then $\delta(v_0) \approx 7$ m s⁻¹, or about a fifth of a channel. For the other three parameters, there is a strong $\tau_0$ dependency, which limits the utility in defining a similar relationship.

## Appendix C: Beyond Keplerian Rotation

With improvements in both the quality of observations and the analysis techniques, it is no longer reasonable to assume a purely Keplerian background rotation profile, as the contribution from either the disk self-gravity (e.g., Veronesi et al. 2021) or the background pressure gradient (e.g., Dullemond et al. 2020) can now be discerned. As discussed in Section 3.1, methods that split the observations into independent annuli can circumvent explicitly dealing with these contributions, as the azimuthally average $v_\phi$ profile can be easily reconstructed by interpolating between each annulus. However, as also discussed in Section 3.1, there is still utility in developing a parameterized profile for situations where strong azimuthal variations in $v_0$ prevent the annulus-by-annulus approach. In this appendix, we discuss a parameterization that was included in the eddy code (Teague 2019a), originally developed to search for the impact of disk self-gravity on the rotation curve, but later found to provide a reasonable description of the sub-Keplerian rotation due to the radial pressure gradient.

To mimic the effect of disk self-gravity on the rotational velocity profile, we have included a disk mass term, $M_{disk}$, such that the total dynamical mass of the system is given by $M_{tot} = M_* + f_d(r) \cdot M_{disk}$, where $f_d(r)$ describes the fraction of the disk mass within a cylindrical radius $r$. Assuming that the disk surface density is described by a power-law profile, $\Sigma(r) = \Sigma_0 (r/r_c)^{-\gamma}$, then the disk mass interior to $r$ is given by

$$f_d(r) = \frac{r^{2-\gamma} - r_{\text{in}}^{2-\gamma}}{r_{\text{out}}^{2-\gamma} - r_{\text{in}}^{2-\gamma}}, \qquad (C1)$$

where $r_{in}$ and $r_{out}$ describe the inner and outer disk edges. For $\gamma = 1$, this is simply a linear interpolation between 0 and 1 from $r_{in}$ to $r_{out}$. Panels (a) and (b) in Figure 14 demonstrate the effect of changing $M_{disk}$ and $\gamma$ on $v_\phi$, respectively.

<figure><img src="image_17.png" /><figcaption>Figure 14. Rotation profiles using the new parameterization to model non-Keplerian components. The top row shows the $v_\phi$ in km s⁻¹, while the bottom row shows them relative to a purely Keplerian rotation profile. Panels (a) and (d) show the impact of changing $M_{disk}$, while panels (b) and (e) show the impact of changing $\gamma$. Panels (c) and (f) demonstrate how a negative $M_{disk}$ in this parameterization can be used to model sub-Keplerian rotation in the outer disk.</figcaption></figure>
While this parameterization is not as accurate as the full form presented in Veronesi et al. (2021) (notably, it does not take into account mass beyond the cylindrical radius *r*), it allows for a simple way to mimic the sub-Keplerian rotation due to pressure support in the outer disk by considering negative values of *M*<sub>disk</sub>. In this form, shown in panels (c) and (f) in Figure 14, *M*<sub>disk</sub> and *γ* can be used to introduce sub-Keplerian motion in regions beyond *r*<sub>in</sub>. Although a direct relation cannot be made to the true pressure support term,

$$v_{\text{prs}}^2 = \frac{r}{\rho} \frac{\partial P}{\partial r}, \qquad (\text{C2})$$

the value of *r*<sub>in</sub> can be broadly interpreted as where the exponential tail of the gas surface density kicks in, while *γ* can be used to quantify how rapidly the rotation curve deviates from Keplerian. For a full exploration of the impact of pressure support, full forward modeling is required (e.g., Dullemond et al. 2020; Yu et al. 2021).

## ORCID iDs

Richard Teague https://orcid.org/0000-0003-1534-5186

Jaehan Bae https://orcid.org/0000-0001-7258-770X

Sean M. Andrews https://orcid.org/0000-0003-2253-2270

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Edwin A. Bergin https://orcid.org/0000-0003-4179-6394

Stefano Facchini https://orcid.org/0000-0003-4689-2684

Jane Huang https://orcid.org/0000-0001-6947-6072

Cristiano Longarini https://orcid.org/0000-0003-4663-0318

David Wilner https://orcid.org/0000-0003-1526-7587

## References

Akiyama, E., Muto, T., Kusakabe, N., et al. 2015, ApJL, 802, L17

Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJL, 820, L40

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

Asensio-Torres, R., Henning, T., Cantalloube, F., et al. 2021, A&A, 652, A101

Bae, J., Teague, R., & Zhu, Z. 2021, ApJ, 912, 56

Boehler, Y., Ménard, F., Robert, C. M. T., et al. 2021, A&A, 650, A59

Bollati, F., Lodato, G., Price, D. J., & Pinte, C. 2021, MNRAS, 504, 5444

Calahan, J. K., Bergin, E., Zhang, K., et al. 2021, ApJ, 908, 8

Calcino, J., Hilder, T., Price, D. J., et al. 2022, ApJL, 929, L25

Casassus, S., & Cárcamo, M. 2022, MNRAS, 513, 5790

Casassus, S., & Pérez, S. 2019, ApJL, 883, L41

Casassus, S., Christiaens, V., Cárcamo, M., et al. 2021, MNRAS, 507, 3789

Cridland, A. J., Bosman, A. D., & van Dishoeck, E. F. 2020, A&A, 635, A68

Czekala, I., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 2

Debes, J. H., Jang-Condell, H., & Schneider, G. 2016, ApJL, 819, L1

Debes, J. H., Jang-Condell, H., Weinberger, A. J., Roberge, A., & Schneider, G. 2013, ApJ, 771, 45

Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205

Disk Dynamics Collaboration, Armitage, P. J., Bae, J., et al. 2020, arXiv:2009.04345

Dong, R., Liu, S.-Y., & Fung, J. 2019, ApJ, 870, 72

Dullemond, C. P., Isella, A., Andrews, S. M., Skobleva, I., & Dzyurkevich, N. 2020, A&A, 633, A137

Dutrey, A., Guilloteau, S., Piétu, V., et al. 2017, A&A, 607, A130

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Huang, J., Andrews, S. M., Cleeves, L. I., et al. 2018, ApJ, 852, 122

Ilee, J. D., Walsh, C., Jennings, J., et al. 2022, MNRAS, 515, L23

Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022, ApJ, 928, 2

Izquierdo, A. F., Testi, L., Facchini, S., Rosotti, G. P., & van Dishoeck, E. F. 2021, A&A, 650, A179

Jorsater, S., & van Moorsel, G. A. 1995, AJ, 110, 2037

Law, C. J., Teague, R., Loomis, R. A., et al. 2021, ApJS, 257, 4

Lenz, D. D., & Ayres, T. R. 1992, PASP, 104, 1104

Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A, 648, A33

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser., Astronomical Data Analysis Software and Systems XVI, 376, ed. R. A. Shaw, F. Hill, & D. J Bell (San Francisco, CA: ASP), 127

Mentiplay, D., Price, D. J., & Pinte, C. 2019, MNRAS, 484, L130

Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icar, 232, 266

Öberg, K. I., Cleeves, L. I., Bergner, J. B., et al. 2021, AJ, 161, 38

Perez, S., Dunhill, A., Casassus, S., et al. 2015, ApJL, 811, L5

Pinte, C., Teague, R., Flaherty, K., et al. 2022, arXiv:2203.09528

Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJL, 860, L13

Pinte, C., van der Plas, G., Ménard, F., et al. 2019, NatAs, 3, 1109

Qi, C., Ho, P. T. P., Wilner, D. J., et al. 2004, ApJL, 616, L11

Rafikov, R. R. 2002, ApJ, 569, 997

Rapson, V. A., Kastner, J. H., Millar-Blanchaer, M. A., & Dong, R. 2015, ApJL, 815, L26

Roberge, A., Weinberger, A. J., & Malumuth, E. M. 2005, ApJ, 622, 1171

Rosenfeld, K. A., Qi, C., Andrews, S. M., et al. 2012, ApJ, 757, 129

Rosotti, G. P., Benisty, M., Juhász, A., et al. 2020, MNRAS, 491, 1335

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91

Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65

Teague, R. 2019, JOSS, 4, 1632

Teague, R. 2019a, JOSS, 4, 1220

Teague, R. 2019b, RNAAS, 3, 74

Teague, R. 2020, richteague/keplerian_mask: Initial Release, v1.0 Zenodo, doi:10.5281/zenodo.4321137

Teague, R., Bae, J., Benisty, M., et al. 2022, ApJ, 930, 144

Teague, R., Bae, J., & Bergin, E. A. 2019a, Natur, 574, 378

Teague, R., Bae, J., Huang, J., & Bergin, E. A. 2019b, ApJL, 884, L56

Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018a, ApJL, 860, L12

Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018b, ApJ, 868, 113

Teague, R., & Foreman-Mackey, D. 2018, RNAAS, 2, 173

Teague, R., & Loomis, R. 2020, ApJ, 899, 157

Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49

Teague, R., Semenov, D., Gorti, U., et al. 2017, ApJ, 835, 228

Teague, R., Henning, T., Guilloteau, S., et al. 2018c, ApJ, 864, 133

Teague, R., Bae, J., Aikawa, Y., et al. 2021, ApJS, 257, 18

Tsukagoshi, T., Muto, T., Nomura, H., et al. 2019, ApJL, 878, L8

van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132

van der Marel, R. P., & Franx, M. 1993, ApJ, 407, 525

Veronesi, B., Paneque-Carreño, T., Lodato, G., et al. 2021, ApJL, 914, L27

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Wölfer, L., Facchini, S., Kurtovic, N. T., et al. 2021, A&A, 648, A19

Yen, H.-W., & Gu, P.-G. 2020, ApJ, 905, 89

Yu, H., Teague, R., Bae, J., & Öberg, K. 2021, ApJL, 920, L33

Zhang, K., Bergin, E. A., Schwarz, K., Krijt, S., & Ciesla, F. 2019, ApJ, 883, 98

Zhu, Z., Bailey, A., Macías, E., Muto, T., & Andrews, S. M. 2022, arXiv:2204.04404

Zuckerman, B., Forveille, T., & Kastner, J. H. 1995, Natur, 373, 494