Mapping the Complex Kinematic Substructure in the TW Hya Disk

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, 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 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.

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

2.2. Imaging

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 ∼15, with a more diffuse outer disk. Two gaps are observed at ∼15 (originally reported in Teague et al. 2017) and ∼26 (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 ¹⁰ 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 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 42 (252 au) and CS is detected out to 39 (235 au), while the continuum edge is around 165 (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 ≈15 (90 au) and ≈26 (156 au).

Zoom In Zoom Out Reset image size

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,

$$\begin{eqnarray}&&I(v)={I}_{0}\exp \left(-{\left(\displaystyle \frac{v-{v}_{0}}{{\rm{\Delta }}V}\right)}^{2}\right),\end{eqnarray} \tag{ 1 }$$

with v₀ describing the line center, ΔでるたV is the Doppler width of the line (where FWHM = 2 $\sqrt{\mathrm{ln}2}\,{\rm{\Delta }}V$) and I₀ is the peak intensity, the CO emission profile was modeled as

$$\begin{eqnarray}&&I(v)={I}_{0}\cdot \displaystyle \frac{1-\exp \left(-\tau (v)\right)}{1-\exp (-{\tau }_{0})},\end{eqnarray} \tag{ 2 }$$

with the optical depth τたう(v) varying as a Gaussian with a width of ΔでるたV and a peak optical depth of τたう₀. Note that τたう₀ 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 χかい² 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 χかい² 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₀, 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.

Zoom In Zoom Out Reset image size

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₀ maps, we use the eddy Python package (Teague 2019a) and fit a Keplerian rotation profile of the form

$$\begin{eqnarray}&&{v}_{\mathrm{kep}}=\sqrt{\displaystyle \frac{{{GM}}_{\star }{r}^{2}}{{\left({r}^{2}+{z}^{2}\right)}^{\tfrac{3}{2}}}},\end{eqnarray} \tag{ 3 }$$

where M_⋆ 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₀, y₀). Thus, for this model, a total of five free parameters are needed: Θしーた_kep = {x₀, y₀, PA, M_⋆, v_LSR}. For face-on disks, there is an extreme degeneracy between M_⋆ and i, so we fix i = 58, 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 = 50; 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_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: Θしーた_prs = {Θしーた_kep, M_disk, γがんま, r_in}. r_in describes the cylindrical radius where v_ϕ begins to deviate from Keplerian rotation, while M_disk dictates the magnitude of this deviation at the outer edge of the disk (for this fitting, we set r_out = 4''). In essence, v_ϕ is a Keplerian profile around a star with a dynamical mass of M_* inside r_in, and approaches a Keplerian profile around a star with a dynamical mass of M_* + M_disk (remembering that M_disk is a negative value in this case) at r_out. How quickly the effective dynamical mass of the system changes with radius is governed by γがんま, with γがんま = 0 resulting in a linear progression and γがんま > 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_ϕ profile changes under various choices of Θしーた_prs are also shown.

For each map, we consider two models: pure Keplerian rotation, v_kep, requiring five free parameters, Θしーた_kep; and pressure support rotation, v_prs, requiring eight free parameters, Θしーた_prs. We model the posterior distributions through eddy, which wraps the emcee MCMC sampling package (Foreman-Mackey et al. 2013). Only regions between 04 (roughly twice the beam FWHM) and out to 375 for CO and 35 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₀. 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 = 1516, M_⋆ = 0.82 M_⊙, and v_LSR = 2843 m s⁻¹. The width of the posteriors (typically used as a measure of the statistical uncertainty) were exceptionally narrow, resulting in fractional uncertainties of ≲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 ∼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_disk = −0.12, γがんま = 1.44, and r_in = 171, 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₀ map calculated from CS emission, likely due to the low S/N of the emission in the regions where the deviation is strongest (r ≳ 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_kep, and a profile including a sub-Keplerian outer disk, v_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_ϕ 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_kep and v_prs as the background model. This is as expected because the model parameters for the pressure support were unable to converge.

Zoom In Zoom Out Reset image size

The most striking feature is a large, tightly wound spiral positive residual, starting at roughly (15, 00), and extending out to the disk edge after almost a full 2πぱい 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 ∼06, and a second in a southwesterly offset from the star at an offset of ∼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 ∼ 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 ∼ 237°) at a distance of 52 au (087).

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,

$$\begin{eqnarray}&&{v}_{0}={v}_{\mathrm{LSR}}+{v}_{\phi ,\mathrm{proj}}+{v}_{r,\mathrm{proj}}+{v}_{z,\mathrm{proj}},\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}&&{v}_{\phi ,\mathrm{proj}}=+{v}_{\phi }\cos (\phi )\sin (\left|i\right|),\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{v}_{r,\mathrm{proj}}=-{v}_{r}\sin (\phi )\sin (i),\end{eqnarray} \tag{ 5b }$$

$$\begin{eqnarray}&&{v}_{z,\mathrm{proj}}=-{v}_{z}\cos (i),\end{eqnarray} \tag{ 5c }$$

where ϕ is the polar angle in the disk, defined to be ϕ = 0° along the redshifted major axis, and increasing in a counterclockwise 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 with strong v_ϕ components crosses the minor axis of the disk, ϕ = ±90°, then the projected velocity v_ϕ,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,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).

Zoom In Zoom Out Reset image size

Around the northeast disk minor axis (PA ∼ 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 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 ∼136 (∼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.

Zoom In Zoom Out Reset image size

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, 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. ¹¹ 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_Nep), 0.3 M_Jup (1 M_Sat), and 0.5 M_Jup, from left to right, convolved with a Gaussian kernel with an FWHM of 018 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_Jup planet, in line with the inferred planet mass based on the size of the gap. Perturbations in the azimuthal direction, v_ϕ , are not shown, as their morphology is inconsistent with the observed structures. However, we do note that, as the v_ϕ 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.

Zoom In Zoom Out Reset image size

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_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 ≳ 0.1. These regions would be therefore traced by CO emission, while the CS emission probes deeper in the disk (z / r ≲ 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 ΔでるたV traced by CO or CS, as shown in Figure 9.

Zoom In Zoom Out Reset image size

A tentative enhancement in the background Δでるた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 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 ¹³CO and C¹⁸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_⊙ and 0.816 M_⊙, respectively. As discussed in Yu et al. (2021), if emission is expected to trace an elevated region parameterized by χかい = 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 ${\left(1+{\chi }^{2}\right)}^{3/2}$. If two different molecules are used to measure a dynamical mass, then the ratio of the inferred dynamical masses, M_⋆,a and M_⋆,b can be used to relate their two emission surfaces, χかい_a and χかい_b , through

$$\begin{eqnarray}&&{\chi }_{b}={\left[\left(1+{\chi }_{a}^{2}\right){\left(\displaystyle \frac{{M}_{\star ,a}}{{M}_{\star ,b}}\right)}^{2/3}-1\right]}^{1/2}.\end{eqnarray} \tag{ 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, χかい_CO ∼ 0.3, while CS is believed to trace closer to the disk midplane, χかい_CS ≲ 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 ≈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., χかい_CO = 0.3 and χかい_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 12 (70 au) and 265 (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.

Zoom In Zoom Out Reset image size

What is particularly interesting is that the radial profile of τたう₀ 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, a subtle bump is seen at 12, 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 12 (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 265 (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.

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 21 (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 135 (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_Jup) 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 (135). 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 Telescope 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).

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 ¹²CO emission presented in this paper. The methodology is as follows. First, we assume the underlying true spectrum is described by

$$\begin{eqnarray}&&I(v)={I}_{0}\cdot \displaystyle \frac{1-\exp \left(-\tau (v)\right)}{1-\exp (-{\tau }_{0})},\end{eqnarray} \tag{ B1 }$$

where I₀ and τたう₀ are the intensity and optical depth at the line center, and where τたう(v) has a Gaussian form,

$$\begin{eqnarray}&&\tau (v)={\tau }_{0}\exp \left(-{\left[\displaystyle \frac{v-{v}_{0}}{{\rm{\Delta }}V}\right]}^{2}\right),\end{eqnarray} \tag{ B2 }$$

where v₀ is the line center and ΔでるたV is the Doppler width (not the typical standard deviation of a Gaussian, but rather a factor of $\sqrt{2}$ larger, such that $\mathrm{FWHM}=2\sqrt{\mathrm{ln}2}\,{\rm{\Delta }}V$). Then, for a given {R, S/N, τたう₀} set, a model spectrum is generated on a v-axis sampled at integer values, with the line center, v₀, 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 σしぐま = I₀ /[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 χかい² 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 ≳ 10³, 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⁴ for any {R, S/N, τたう₀} set.

We then sample the {R, S/N, τたう₀} parameter space. Four values of τたう₀ were chosen: 0 (a pure Gaussian model to compare with Lenz & Ayres 1992), 1, 3, and 5. For τたう₀ > 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, τたう₀} set is shown in Figure 13. Both the line center, v₀, first column, and the line width, Δでるた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})/{\rm{\Delta }}{V}_{{\rm{chan}}}=0.69\,\sqrt{R}\,/\,[{\rm{S}}/{\rm{N}}]$ when averaged over all four τたう₀ 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 ≈ 6.7. Assuming that the S/N ∼ 10 for a given pixel, then δでるた(v₀) ≈ 7 m s⁻¹, or about a fifth of a channel. For the other three parameters, there is a strong τたう₀ dependency, which limits the utility in defining a similar relationship.

Zoom In Zoom Out Reset image size

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_ϕ 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₀ 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) · 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, ${\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{0}\,{\left(r/{r}_{c}\right)}^{-\gamma }$, then the disk mass interior to r is given by

$$\begin{eqnarray}&&{f}_{d}(r)=\displaystyle \frac{{r}^{2-\gamma }-{r}_{\mathrm{in}}^{2-\gamma }}{{r}_{\mathrm{out}}^{2-\gamma }-{r}_{\mathrm{in}}^{2-\gamma }},\end{eqnarray} \tag{ C1 }$$

where r_in and r_out describe the inner and outer disk edges. For γがんま = 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 γがんま on v_ϕ , respectively.

Zoom In Zoom Out Reset image size

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_disk. In this form, shown in panels (c) and (f) in Figure 14, M_disk and γがんま can be used to introduce sub-Keplerian motion in regions beyond r_in. Although a direct relation cannot be made to the true pressure support term,

$$\begin{eqnarray}&&{v}_{\mathrm{prs}}^{2}=\displaystyle \frac{r}{\rho }\displaystyle \frac{\partial P}{\partial r},\end{eqnarray} \tag{ C2 }$$

the value of r_in 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).