# Molecules with ALMA at Planet-forming Scales (MAPS). II. CLEAN Strategies for Synthesizing Images of Molecular Line Emission in Protoplanetary Disks

Ian Czekala, Ryan A. Loomis, Richard Teague, Alice S. Booth, Jane Huang, Gianni Cataldi, John D. Ilee, Charles J. Law, Catherine Walsh, Arthur D. Bosman, et al.

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

## ► To cite this version:

Ian Czekala, Ryan A. Loomis, Richard Teague, Alice S. Booth, Jane Huang, et al., Molecules with ALMA at Planet-forming Scales (MAPS). II. CLEAN Strategies for Synthesizing Images of Molecular Line Emission in Protoplanetary Disks. The Astrophysical Journal Supplement Series, 2021, 257, 10.3847/1538-4365/ac1430. insu-03672379

HAL Id: insu-03672379

https://insu.hal.science/insu-03672379v1

Submitted on 19 Sep 2025

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_2.png" /></figure>

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

# Molecules with ALMA at Planet-forming Scales (MAPS). XVIII. Kinematic Substructures in the Disks of HD 163296 and MWC 480

Richard Teague<fnref n="1" />, Jaehan Bae<fnref n="2" />,<fnref n="3" />,<fnref n="23" />, Yuri Aikawa<fnref n="4" />, Sean M. Andrews<fnref n="1" />, Edwin A. Bergin<fnref n="5" />, Jennifer B. Bergner<fnref n="6" />,<fnref n="23" />, Yann Boehler<fnref n="7" />, Alice S. Booth<fnref n="8" />,<fnref n="9" />, Arthur D. Bosman<fnref n="5" />, Gianni Cataldi<fnref n="4" />,<fnref n="10" />, Ian Czekala<fnref n="11" />,<fnref n="12" />,<fnref n="13" />,<fnref n="14" />,<fnref n="15" />,<fnref n="23" />, Viviana V. Guzmán<fnref n="16" />, Jane Huang<fnref n="1" />,<fnref n="5" />,<fnref n="23" />, John D. Ilee<fnref n="9" />, Charles J. Law<fnref n="1" />, Romane Le Gal<fnref n="1" />,<fnref n="7" />,<fnref n="17" />,<fnref n="18" />, Feng Long<fnref n="1" />, Ryan A. Loomis<fnref n="19" />, François Ménard<fnref n="7" />, Karin I. Öberg<fnref n="1" />, Laura M. Pérez<fnref n="20" />, Kamber R. Schwarz<fnref n="21" />,<fnref n="23" />, Anibal Sierra<fnref n="20" />, Catherine Walsh<fnref n="9" />, David J. Wilner<fnref n="1" />, Yoshihide Yamato<fnref n="4" />, and Ke Zhang<fnref n="5" />,<fnref n="22" />,<fnref n="24" />

## Abstract

We explore the dynamical structure of the protoplanetary disks surrounding HD 163296 and MWC 480 as part of the Molecules with ALMA at Planet-forming Scales (MAPS) large program. Using the *J* = 2–1 transitions of ¹²CO, ¹³CO, and C¹⁸O imaged at spatial resolutions of ~0.′′15 and with a channel spacing of 200 m s⁻<fnref n="1" />, we find perturbations from Keplerian rotation in the projected velocity fields of both disks (≤5% of the local Keplerian velocity), suggestive of large-scale (tens of astronomical units in size), coherent flows. By accounting for the azimuthal dependence on the projection of the velocity field, the velocity fields were decomposed into azimuthally averaged orthogonal components, *v*<sub>ϕ</sub>, *v*<sub>r</sub>, and *v*<sub>z</sub>. Using the optically thick ¹²CO emission as a probe of the gas temperature, local variations of ≈3 K (≈5% relative changes) were observed and found to be associated with the kinematic substructures. The MWC 480 disk hosts a suite of tightly wound spiral arms. The spirals arms, in conjunction with the highly localized perturbations in the gas velocity structure (kinematic planetary signatures), indicate a giant planet, ~1 *M*<sub>Jup</sub>, at a radius of ≈245 au. In the disk of HD 163296, the kinematic substructures were consistent with previous studies of Pinte et al. and Teague et al. advocating for multiple ~1 *M*<sub>Jup</sub> planets embedded in the disk. These results demonstrate that molecular line observations that characterize the dynamical structure of disks can be used to search for the signatures of embedded planets. This paper is part of the MAPS special issue of the Astrophysical Journal Supplement.

*Unified Astronomy Thesaurus concepts: Interferometry (808); Millimeter astronomy (1061); Exoplanet formation (492); Protoplanetary disks (1300)*

## 1. Introduction

It is clear that the Atacama Large Millimeter/submillimeter Array (ALMA) has driven a revolution in our understanding of the physical and chemical structure of protoplanetary disks, the birthplaces of planets. Not only has it revealed a stunning variety and ubiquity of substructures in the millimeter continuum, such as gaps, rings, and vortices (e.g., ALMA

Partnership et al. 2015; Andrews et al. 2018; Long et al. 2018), but also a similar level of complexity in the morphology of molecular line emission (e.g., Christiaens et al. 2014; Teague et al. 2017, 2019b; Huang et al. 2020). While the former is likely due to the grains being shepherded by changes in the global gas physical structure (Dullemond et al. 2018), identifying the cause of changes in molecular emission requires distinguishing between physical and chemical effects (Facchini et al. 2018; van der Marel et al. 2018).

The unparalleled sensitivity of ALMA has enabled an alternative approach to detailing the physical structure of the

<footnote n="1">¹Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA; richard.d.teague@cfa.harvard.edu</footnote>

<footnote n="2">²Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA</footnote>

<footnote n="3">³Department of Astronomy, University of Florida, Gainesville, FL 32611, USA</footnote>

<footnote n="4">⁴Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan</footnote>

<footnote n="5">⁵Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, USA</footnote>

<footnote n="6">⁶University of Chicago, Department of the Geophysical Sciences, Chicago, IL 60637, USA</footnote>

<footnote n="7">⁷Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="8">⁸Leiden Observatory, Leiden University, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="9">⁹School of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, UK</footnote>

<footnote n="10">¹⁰National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="11">¹¹Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="12">¹²Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="13">¹³Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="14">¹⁴Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="15">¹⁵Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA</footnote>

<footnote n="16">¹⁶Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile</footnote>

<footnote n="17">¹⁷IRAP, Université de Toulouse, CNRS, CNES, UT3, 31400 Toulouse, France</footnote>

<footnote n="18">¹⁸IRAM, 300 rue de la piscine, F-38406 Saint-Martin d'Hères, France</footnote>

<footnote n="19">¹⁹National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA</footnote>

<footnote n="20">²⁰Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="21">²¹Lunar and Planetary Laboratory, University of Arizona, 1629 E. University Blvd., Tucson, AZ 85721, USA</footnote>

<footnote n="22">²²Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St., Madison, WI 53706, USA</footnote>

<footnote>Received 2021 February 15; revised 2021 May 20; accepted 2021 June 14; published 2021 November 3</footnote>

<footnote n="23">²³NASA Hubble Fellowship Program Sagan Fellow.</footnote>

<footnote n="24">²⁴NASA Hubble Fellow.</footnote>
gas: the gas dynamics. A significant advantage of a dynamical study of the planet formation environment is that it circumvents the need to account for the chemical effects when relating measured fluxes to column densities of H₂ gas. This is most clearly seen with the pressure modulation of the rotation velocity,

$$ \frac{v_{\phi}^2}{r} = \frac{GM_*r}{(r^2 + z^2)^{3/2}} + \frac{1}{\rho_{\text{gas}}} \frac{\partial P}{\partial r}, \quad (1) $$

where r and z are the cylindrical radius and vertical height above the midplane, respectively; $M_*$ is the stellar mass; $\rho_{\text{gas}}$ is the gas density, dominated by H₂; and P is the gas pressure. Changes in the local gas pressure, or gas surface density, will therefore manifest as changes in $v_{\phi}$, irrespective of the abundance of the observed molecule. The advantage of this approach is clear: variations in the total gas column can be inferred without the need for expensive astrochemical models. For example, Teague et al. (2018a, 2018b) detected subtle variations in $v_{\phi}$ for the disks around HD 163296 and AS 209, which were used to model the perturbed gas surface density profile.

It is well understood that giant planets will open up gaps in the gas distribution of their parental protoplanetary disks when they reach a sufficient mass (e.g., Papaloizou & Lin 1984). While this is an attractive, and frequently invoked, scenario to explain perturbations in the gas surface density, other dynamical processes, such as the magnetorotational instability (MRI; Flock et al. 2015) or the vertical shear instability (VSI; Flock et al. 2017), have been demonstrated to result in comparable perturbations. A way to distinguish between these scenarios is to leverage the fact that planets should only drive localized perturbations, as the shocks from their spirals are strongest in their immediate vicinity. Such localized perturbations (~10% of the background Keplerian rotation) in the velocity structure of HD 163296 and HD 97048 were reported by Pinte et al. (2018b, 2019), which the authors showed were consistent with giant, ~2 $M_{\text{Jup}}$ planets.

Embedded planets are also expected to drive large velocity perturbations in all three components of the gas velocity, $v_{\phi}$, $v_r$, and $v_z$ (e.g., Perez et al. 2015; Pérez et al. 2018; Pinte et al. 2019; Disk Dynamics Collaboration et al. 2020). These can be probed by decomposing the line-of-sight velocity $v_0$, which is the sum of the systemic velocity, $v_{\text{LSR}}$, and the projected disk velocity field, $v_{\text{disk}}$, into the three cardinal velocity components,

$$ v_{\text{disk}} = v_{\phi} \cos(\phi)\sin(i) + v_r \sin(\phi)\sin(i) + v_z \cos(i), \quad (2) $$

where $\phi$ is the azimuthal angle of the disk (measured such that $\phi=0$ is the redshifted major axis) and $i$ is the inclination of the disk. Teague et al. (2019a, 2019b) demonstrated that these components can be disentangled if we know their geometric projection and assume azimuthal asymmetry. A similar approach is frequently used in studies of galactic dynamics (e.g., Krajnović et al. 2006).

There has been much recent success in using the high spatial and spectral observations from ALMA to probe the gas dynamics and infer the presence of unseen planets by the detection of kinematic planetary signatures (KPSs; Disk Dynamics Collaboration et al. 2020). This term encompasses all localized features in the gas dynamical structure that are driven by an unseen planetary perturber, such as “kinks” or “wiggles” in the emission morphology or “Doppler-flips” and spiral patterns in rotation map residuals. Pinte et al. (2018b)

found localized velocity deviations in the outer disk of HD 163296, consistent with a ~2 $M_{\text{Jup}}$ planet at a radius of ~260 au. In the same source, Teague et al. (2019a) used $^{12}$CO emission to measure the azimuthally averaged 3D velocity structure, revealing a large gas pressure minimum at the orbital radius of the planet and evidence for meridional flows: the motion of gas from the atmosphere of the disk into (likely planet-opened) gaps (Morbidelli et al. 2014; Szulágyi et al. 2014). Other KPSs have been found in AB Aur (Tang et al. 2017), HD 97048 (Pinte et al. 2019), HD 100546 (Casassus and Pérez 2019; Pérez et al. 2020), and TW Hya (Teague et al. 2019b), along with several tentative features in sources targeted by the DSHARP survey (Pinte et al. 2020). These results demonstrate not only that we are on the verge of painting a six-dimensional view of the planet formation process but also that we have the potential to transform ALMA into a (sub) millimeter planet-hunting instrument. Observations of molecular line emission will enable a unique quantification of the planet-disk interactions in a less ambiguous way than with (sub)millimeter continuum alone.

The Molecules with ALMA at Planet-forming Scales (MAPS) large program (Öberg et al. 2021) produced some of the deepest and highest angular resolution observations of molecular line emission from protoplanetary disks to date with a view to characterize the chemical complexities associated with the planet formation process. These observations are ideal for also characterizing the dynamical structure of the planet-hosting disks. In this paper, we use these data to analyze the velocity structures of the disks around HD 163296 and MWC 480, the two Herbig Ae stars in the program. These two sources were selected because the other three MAPS sources contained complex $^{12}$CO emission morphologies: GM Aur was found to have highly complex extended emission, presented in Huang et al. (2021); AS 209 has substantial cloud contamination over a large region of the disk (e.g., Huang et al. 2016); and IM Lup has extensive diffuse $^{12}$CO emission, potentially associated with a photoevaporative wind (e.g., Cleeves et al. 2016). In Section 2 we describe the observations used and measure the line-of-sight velocities for the bright $^{12}$CO emission in Section 3. In Section 4 these are disentangled into azimuthally symmetric orthogonal components. Section 5 describes the local changes in the gas temperature that are associated with the perturbations in the gas dynamics. Finally, these perturbations are discussed in the context of embedded planets in Section 6, with the results summarized in Section 7.

## 2. Observations

As part of the MAPS program, both HD 163296 and MWC 480 were observed with the correlator tuned to cover the $J=2-1$ rotational transitions of $^{12}$CO, $^{13}$CO, and C$^{18}$O so as to trace a range of vertical heights in the disk (Law et al. 2021b). We refer the reader to Öberg et al. (2021) for a full description of the observational setup and calibration process and Czekala et al. (2021) for a thorough description of the imaging process.

For this paper we use the robust = 0.5 weighted, $JvM$-corrected$^{25}$ images, as these produced the smallest synthesized beam. This resulted in beam sizes of approximately

<footnote>$^{25}$ A correction related to the mismatch in CLEAN beam and dirty beam sizes when combining a CLEAN image and the residuals (Jorsater & van Moorsel 1995; Czekala et al. 2021).</footnote>
$0''14 \times 0''11$ at a position angle of $104^\circ$ for the CO isotopologue emission in HD 163296 and $0''17 \times 0''12$ at a position angle of $6^\circ$ for isotopologue emission in MWC 480. All images were produced with a channel spacing of 200 m s$^{-1}$. The rms was measured as the standard deviation of pixel values in a circular area of 1" in radius centered on the phase center of the first and last channel ($\pm 12$ km s$^{-1}$ from the systemic velocity) in the image cube. For HD 163296, the rms was measured to be 0.58, 0.54, and 0.37 mJy beam$^{-1}$ for $^{12}$CO, $^{13}$CO, and $^{C^{18}O}$, respectively. In MWC 480, these were 0.73, 0.68, and 0.36 mJy beam$^{-1}$. A set of channel maps for the $^{12}$CO, $^{13}$CO, and $^{C^{18}O}$emission can be found in Appendix D.

### 3. Rotation Maps

In this section, we collapse the image cubes along the spectral axis to study their projected velocity fields and search for any deviations relative to a background Keplerian model.

#### 3.1. Method

We use the “quadratic” method implemented in better-moments to produce maps of the line center, $v_0$ maps, including a statistical uncertainty for the measured $v_0$ as described in Teague & Foreman-Mackey (2018). The maps are then masked to only include regions where the peak intensities are greater than 5 times the rms value measured in a line-free channel in order to remove noisy values at the outer edge of the disk. We explored whether an additional smoothing along the spectral axis prior to the collapse of the data cube would improve the fidelity of the rotation map (as advocated by Teague & Foreman-Mackey 2018 for low signal-to-noise ratio spectra); however, we found that this did not significantly improve the resulting $v_0$ maps. Appendix A discusses this in more detail.

To each of these maps we fit a Keplerian rotation model, $v_{\text{kep}}$, given by Equation (1) when assuming $\partial P/\partial r = 0$. To find the best-fit Keplerian velocity field, $v_{\text{kep}}$, we use the Python package `eddy` (Teague 2019a) to explore the posterior distributions of the model parameters. For these models there is a strong degeneracy between the disk inclination and the stellar mass that can usually be broken with high-resolution observations such as those from the MAPS project. However, if there are perturbations in the velocity structure, these will dominate the residuals and thus preclude a tight constraint on the inclination. As such, incorporating prior information on the disk inclination, from the continuum emission, for example, is preferable. We adopt the inclination and binned $z(r)$ profiles from Huang et al. (2018a) and Law et al. (2021b), respectively, for HD 163296. For MWC 480, we found that the velocity structure in the outer disk is more consistent with an inclination of 33°, rather than the 37° based on continuum fitting (Liu et al. 2019) and that was adopted for other MAPS papers. As systematic residuals from a misspecified inclination dominate our residual map (see, e.g., Yen & Gu 2020), we choose to use the 33° inclination for the entire MWC 480 disk and recalculate the $z(r)$ profiles following the method described in Law et al. (2021b). Appendix B provides more details on this. Other geometrical properties, including source center ($x_0, y_0$), position angle PA, dynamical mass $M_*$, and systemic velocity $v_{\text{LSR}}$, are left free. Appendix C discusses the model-fitting procedure in more detail.

#### 3.2. HD 163296

The $v_0$ maps of HD 163296 for all three CO isotopologues are shown in the left column of Figure 1, with each isotopologue on a separate row. The color scale has been saturated in the central regions to bring out the detail along the minor axes of the disk. Solid contour lines show isovelocity contours in steps of 500 m s$^{-1}$ in this region. The middle column shows the residual after subtracting a best-fit Keplerian mode, $v_{\text{kep}}$. To suppress the noise in the residual map, these are smoothed by convolution with a circular Gaussian kernel with an FWHM of 2 pixels (see Appendix A for a comparison of how the smoothing affects the residuals). Overlaid are contours of constant radius and azimuth in the disk frame of reference, in steps of $0''5$ and $30^\circ$, respectively. In the inner $\approx 0''3$ ($\approx 30$ au), beam convolution effects preclude a good fit of the data and are therefore masked in the residual maps. In the right column, the residual map has been deprojected into a face-on view, taking account of the disk inclination, position angle, and emission surface. The north-facing minor axis is oriented along the positive y-axis. In all these panels, the disk is rotating in a clockwise fashion. Note that for HD 163296, continuum emission is only detected out to a radius of $\sim 1''7$ (Isella et al. 2016), such that the gas is over 3 times larger in radius.

We compare the rotation maps for $^{12}$CO emission made with MAPS data to those made with data from DSHARP (Andrews et al. 2018; Isella et al. 2018) in Figure 1 for HD 163296. The DSHARP data were taken at a slightly higher angular resolution ($\approx 100$ mas vs. $\approx 125$ mas for MAPS), but at a lower effective velocity resolution (640 m s$^{-1}$ vs. the binned down 200 m s$^{-1}$ of the MAPS data). There are two striking features in this comparison. First, there is the extraordinary similarity in the observed residuals, despite being entirely independent data sets. This strongly suggests that with ALMA we are hitting the “confusion limit” for rotation maps: that is, our residuals are now dominated by real variations in the gas kinematics rather than thermal noise. Second, there is considerable “feathering” (radial features that result in a beating pattern as a function of azimuth) observed in the DSHARP data that is due to the channelization of the data. This is exacerbated by the low effective velocity resolution of the data (640 m s$^{-1}$) and is mostly removed with the finer effective velocity resolution achieved by MAPS (Öberg et al. 2021).

We can identify several dominant features in the residual map of $^{12}$CO emission, as annotated in Figure 2. First, there are large positive and negative residuals along the blueshifted and redshifted major axes, respectively, labeled “pressure.” These can be attributed to the sub-Keplerian rotation of the outer disk, driven by the large pressure gradient as the gas density drops at the edge of the disk (Rosenfeld et al. 2013; Dullemond et al. 2020). There is also mild contamination from the rear side of the disk, along the top of the plot, spanning between 230° and 310°.

We also see three large perturbations, marked by the arc-shaped regions in Figure 2 centered at $(r, \phi) = (2''85, 305^\circ)$, $(1''65, 315^\circ)$, and $(0''9, 220^\circ)$, labeled A, B, and C, respectively. Feature A is the manifestation of the “kink” KPS in channel maps reported by Pinte et al. (2018b), which the authors attribute to a $\approx 2 M_{\text{Jup}}$ planet. The innermost perturbation, C, is coincident with the “kink” at 86 au reported by Pinte et al. (2020), while the extended arm closely traces the azimuthal asymmetry detected in the millimeter continuum (Isella et al. 2018) in azimuthal extent but appears to be radially
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Gallery of $v_0$ maps for different $J = 2-1$ CO isotopologue emission from the disk around HD 163296 using both MAPS and DSHARP data (Andrews et al. 2018; Isella et al. 2018). The left column shows the raw $v_0$ map, the middle column shows the residuals after subtracting a best-fit Keplerian rotation model, and the right column shows the residuals deprojected into a face-on orientation accounting for the elevated emission surface. Solid colored contours in the left panel show isovelocity contours of $v_{LSR} \pm \{2, 2.5, 3, \dots\}$ km s$^{-1}$. Dashed black contours in the middle and right columns show contours of constant radius and azimuth in the disk frame, with intervals of 0.5 and 30°, respectively. The beam size is shown in the lower left corner of each panel. For the deprojected maps in the right column, the beam size is only representative, as the deprojection results in a nonuniform spatial resolution across the image. The disk rotates clockwise in all panels.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

**Figure 2.** A zoom-in of Figure 1(c). A representative beam size is shown in the lower right; however, due to the elevated emission surface, the true angular resolution varies across the map. The figure extends out to 4′′4, and the color scaling ranges between residuals of ±180 m s⁻¹. Residuals due to the pressure gradient slowing the rotation of the disk are observed along the major axis, marked as “pressure.” Three regions of large velocity deviations are marked by black boxes, centered at (r, φ) = (2′′85, 305°) (A), (1′′65, 315°) (B), and (0′′9, 220°) (C), all of which have extended structure shown by the dashed lines. The outer feature, A, is associated with the planet identified by Pinte et al. (2018b), while the inner two, B and C, are likely associated with the potential “kink” reported by Pinte et al. (2020), or due to the asymmetry in the continuum emission (Isella et al. 2018). Two arc-like structures are also seen at radii of 1′′15 (116 au), and 1′′8 (182 au), shown by the solid lines, labeled D and E, respectively, in addition to a larger arc feature at 4′′1 (414 au), shown by the thick solid line with a faded end, annotated as F. It is unclear where this outermost arc ends, as the perturbation extends into a region dominated by pressure-modulated deviations.

offset to larger orbital radii by ~0′′1. This is discussed more in Section 6.1. The middle perturbation, B, has no previous detection but lies at a radius that is broadly consistent with the 145 au (1′′45) gap in the millimeter continuum (Isella et al. 2016, 2018). Both these features appear to have a large azimuthal extent, traced by the dashed lines. In particular, this demonstrates that the previously detected “kink” KPS feature is far larger and more complex than previously thought from previous observations, the cause of which is discussed more in Section 6.1.

In addition to these larger deviations, we find evidence of three arc structures: two well-defined arcs in the inner disk at 115 and 181 au (1′′15 and 1′′8), marked with solid lines labeled D and E, respectively, and one large spiral along the southwest edge of the disk at 414 au (4′′1), shown by the solid line with a faded end, labeled F. As the latter arc blends into the perturbation from the pressure gradient, we are unable to constrain its azimuthal extent.

The ¹³CO $v_0$ map, residual from a Keplerian rotation pattern, and deprojected residual map are shown in Figures 1(g), (h), and (i), respectively. Large systematic residuals are seen along the northwest of the disk, similar in morphology to those seen in ¹²CO emission, but with much larger velocities. The morphology of these residuals is similar to that proposed by Yen & Gu (2020) for a misspecified emission surface. As with the ¹²CO emission, this arises because the emission starts to become thin in the outer edge of the disk and thus traces both the front and rear side of the disk. The average height traced by the emission is then underestimated by the emission surface,

giving rise to these residuals. The pressure gradient in the outer disk also introduces small perturbations along the major axis of the disk, as for ¹²CO.

Beyond these features, there is one notable residual: a localized perturbation at a radial offset of ≈157 au (1′′55), coincident with the large velocity residual annotated as B in Figure 2. Tracing the feature in two isotopologues points toward a real velocity deviation rather than noise. The coincidence of this feature with the gap in the dust continuum, D145, is suggestive of a planetary origin. However, the similar azimuthal angle to the previously discussed planet at ≈260 au raises the question of whether these features are related to the same origin, a scenario discussed more in Section 6.1.

Maps showing $v_0$, the residuals from Keplerian rotation, and spatially deprojected residuals from Keplerian rotation for C¹⁸O emission are plotted in Figures 1(j), (k), and (l), respectively. The residual maps are dominated by large systematic residuals, similar in morphology to those seen in the northern half of the disk traced by ¹²CO and ¹³CO emission, but also reflected in the southern half of the disk. Yen & Gu (2020) argue that residuals with this morphology indicate an incorrect inclination. However, for the C¹⁸O emission this is rather due to emission arising from very close to the midplane such that there is little difference in the brightness of emission arising from the underside of the molecular layer (i.e., from the backside of the disk), or the top side of the molecular layer (i.e., from the front side of the disk). For highly abundant species, such as ¹²CO, the difference in gas temperature at the top of the molecular layer (set by the dissociation of CO by UV radiation; $T \gtrsim 50$ K) and the bottom of the molecular layer (set by the freeze-out of CO; $T \sim 21$ K) is sufficiently large that separating emission from these two regions is simple. For less abundant species, such as C¹⁸O, the top of the molecular layer is much deeper in the disk and thus spans a much narrower temperature range, the result being that disentangling between near-side and far-side emission is much tougher, leading to the line emission tracing an average height much closer to the midplane. This is very similar to the optically thin emission argument made for ¹²CO and ¹³CO emission above; however, we stress that while the outcome is the same in terms of residuals, the reason for the bias in traced height is different.

### 3.3. MWC 480

The $v_0$ maps, residuals from Keplerian rotation, and spatially deprojected residual maps for CO isotopologue emission from the disk of MWC 480 are shown in Figure 3. In all panels, the disk rotation is anticlockwise. In this source, all velocity deviations are smaller in magnitude than those observed in HD 163296, both in terms of absolute value and relative to the local Keplerian rotation. As with HD 163296, continuum emission from the disk around MWC 480 has been detected out to ≈1′′7, meaning that the gas disk, as traced by ¹²CO, is up to three times larger in radius. Note that the millimeter continuum edge reported in MAPS III (Law et al. 2021a) provides a more conservative estimate of the outer disk as the region that encompasses 90% of the continuum flux.

As MWC 480 was not a target of DSHARP, we only consider the MAPS data. An annotated version of Figure 3(c) is shown in Figure 4. Similar to HD 163296, there are systematic residuals that can be attributed to large negative pressure gradients; however, there is no obvious contamination from the backside of the disk, which would show up along the bottom of
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Same as Figure 1, but for MWC 480. As this disk was not targeted by the DSHARP survey, only the MAPS data are shown. The disk rotates in an anticlockwise direction in all frames. The beam size is shown in the lower left corner of each panel. For the deprojected maps in the right column, the beam size is only representative, as the deprojection results in a nonuniform spatial resolution across the image.</figcaption></figure>

the panel. Instead, a systematic residual with a positive value in the eastern quarter of the disk and a negative value in the northern quarter (i.e., either side of the minor axis in the NE half of the disk in Figure 3(b), or either side of 90° in the top of the disk in Figures 3(c) and 4) is seen. Yen & Gu (2020) demonstrate that such a morphology can be attributed to a misspecified emission surface. As discussed in Law et al. (2021b), the ¹²CO emission surface in MWC 480 is only well constrained inward of 2'', so a misspecified surface is the most likely the cause for this feature and annotated as such.

Aside from these systematic structures, the kinematic substructure in MWC 480 is dominated by annular structures, with the most striking ring of redshifted emission at 1''5, labeled as box A, and then two broad arcs slightly inside and outside this feature at 1'' and 2''2, labeled as B and C, respectively. There are noisy features observed at the same

radius but extending further in azimuth, suggesting that these arcs may extend over a broader azimuth than annotated; however, deeper observations would be required to fully constrain their extent. The velocity perturbations appear to be strongest around $\phi \sim 5^\circ$, as marked by the two boxed regions in Figure 4.

The kinematics traced by ¹³CO emission are shown in Figures 3(d)–(f). There is a substantial difference in the magnitude of the perturbations seen in the velocity traced by ¹³CO and those traced by ¹²CO emission. There appears to be a subtle positive-to-negative residual across the minor axis of the disk, suggesting a slightly sub-Keplerian rotation of the gas (as discussed with regard to the pressure modulation for ¹²CO).

Subtle deviations that appear to trace the same regions as those found in the ¹²CO emission are detected, most notably, three arc-like structures in the northwest of the disk at the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Annotated version of Figure 3(c), showing the features observed in the residual rotation map from ¹²CO in the disk of MWC 480. The disk rotates in an anticlockwise fashion. Residuals due to the pressure gradient slowing the rotation of the disk are observed along the major axis, marked as “pressure,” with similar systematic features likely due to a misspecified emission surface marked as “emission surface.” The most prominent is an arc extending almost the full azimuth of the disk at 1/5 with a peak marked by a solid annulus, labeled A, with an inner arc at 0/95, labeled B. An additional tentative arc is annotated as arc C at 2/1. A representative beam size is shown in the lower right corner; however, due to the elevated emission surface, the true beam size varies across the map.</figcaption></figure>

locations of the peak velocity perturbations in ¹²CO. The relationship between these structures and those seen in ¹²CO is discussed in Section 6.2.

The bottom row of Figure 3 shows the C¹⁸O $v_0$ map, residuals from Keplerian rotation, and a deprojected map of the residuals, respectively. Little, if any, structure is observed in the residuals. A very low-level “three-spur” pattern is observed, suggesting a misspecified inclination (Yen & Gu 2020). As discussed above, and in more detail in Appendix B, the inclination of the MWC 480 disk is hard to constrain, and the lack of residuals seen for ¹²CO and ¹³CO, which predominantly trace the outer regions of the disk than C¹⁸O, may suggest a radially varying inclination.

# 4. Velocity Profiles

While a $v_0$ map provides an excellent summary of the dominant velocity fields in a disk, the spectrally resolved nature of the observations provides a wealth of additional dynamical information that is lost when collapsing the data cube. In this section we use the full image cube to derive radial profiles for azimuthally averaged rotational, radial, and vertical motions.

## 4.1. Method

To extract the rotational and radial velocities, $v_φ$ and $v_r$, respectively, we use the Python package GoFish (Teague 2019b). The method is fully described in Teague et al. (2018b) and extended to three dimensions in Teague et al. (2019a). Nonetheless, a brief overview of the method is provided below.

First, the disk is split into concentric annuli with a width of a quarter of the beam major FWHM, taking into account the emission surface for the ¹²CO as defined in Law et al. (2021b). For each annulus, a random sample of spatially independent pixels (i.e., those at least a beam major FWHM apart) are

selected, which, under the assumption of an azimuthally symmetric disk, should all be identical in shape and only differ by their line center, which varies owing to the projected velocity component of the disk. The velocity axis for each spectrum is then corrected for the disk velocity structure using Equation (2), assuming $v_z = 0$. The reason for not including a $v_z$ component in this correction is that the projection of this term has no azimuthal dependence and thus only results in an overall shift of the spectral axis, making it degenerate with $v_{LSR}$. The correct $v_φ$ and $v_r$ values are defined as those that best align the individual spectra, as quantified by fitting a Gaussian processes mean-model to the averaged spectra. In essence, this minimizes the standard deviation between the sample of shifted spectra without having to make assumptions about the underlying line profile, i.e., assuming an analytical profile. There will be some intrinsic variation in the line profile as a function of azimuth owing to projection effects; however, such variations are expected to be negligible compared to thermal noise of the data (see the discussion in Teague & Loomis 2020, for example).

This process is repeated 20 times for each annulus, each time taking a new, random selection of pixels in order to minimize features in the radial profiles that can be attributed to noisy pixels (e.g., Keppler et al. 2019). The posterior distributions from each sample are combined, and the results are shown in the top three rows of Figure 5. As $v_φ$ is dominated by the Keplerian rotation, a background model is subtracted to highlight the small-scale structure, as in the second row of Figure 5. As discussed in Teague et al. (2019a), while the choice of background model does not strongly affect these residuals, different background models will lead to small changes in the overall gradient and thus the measured “peak-to-trough” values (to measure spatial properties of pressure perturbations, for example; Yun et al. 2019). Such an approach will only highlight the high spatial frequency (radially localized) perturbations, rather than larger-scale perturbations, such as the sub-Keplerian rotation identified in the outer disk of HD 163296. This is demonstrated by the inset axes in Figure 5 (a), which compares a purely Keplerian rotation with a black dotted line to the observed $v_φ$ profile in red and the 10th-order polynomial fit used for Figure 5(c), shown by a gray dashed line. As such, these panels should only be used for a qualitative analysis, and any quantitative analysis must be done with the $v_φ$ profiles directly.

To measure the vertical velocity component, an analytical Gaussian fit to the aligned and averaged spectra must be performed in order to recover the line center. The difference between the measured line center and the systemic velocity is subsequently attributed to any vertical velocity components. When aligning and averaging the spectra for each annulus, emission from the backside of the disk was found to contaminate the line wings, resulting in a broad, pedestal-like component. To remove this from the fitting, only velocities within 500 m s⁻¹ of the $v_{LSR}$ were considered in the fit, which was verified to not include any wing emission. To improve the precision of the fit, the averaged spectrum is binned onto a velocity axis with a channel spacing of 50 m s⁻¹, a factor of 4 finer than the native data. This is possible because the velocity shifts due to the rotation of the disk result in subchannel shifts, meaning that the average spectrum is sampled at a much higher rate than the channel spacing. As discussed in Teague & Loomis (2020), this approach does not remove any systematic
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Azimuthally averaged radial profiles of the gas velocities in HD 163296 (left) and MWC 480 (right), traced by ¹²CO. Panels (a) and (b) show $v_φ$, while panels (c) and (d) show the residual of $v_{\phi\text{phi}}$ after subtracting a baseline $v_{\text{mod}}$. The radial and vertical velocities are shown in panels (e) and (g) and panels (f) and (h) for HD 163296 and MWC 480, respectively. All velocities have been corrected for the inclination of the disk. Following Teague et al. (2019a), the baseline model, $v_{\text{mod}}$, used to highlight localized deviation in $v_φ$ is a 10th-order polynomial fit to $v_φ$. In panel (a), an inset highlights the difference between a purely Keplerian rotation profile (black dotted line) and the observed rotation profile and 10th-order polynomial fit (gray dashed line) used for panel (c). The shaded regions show the 1σ uncertainties. Locations of the gaps and rings in the continuum are annotated in the lower three rows as vertical dashed lines, with the gray and orange shading showing their width (Huang et al. 2018a). Note that the sign convention for radial velocity has been inverted relative to Teague et al. (2019b) such that positive $v_r$ represents gas moving away from the star.</figcaption></figure>

broadening effects, for example, due to the spectral response of the correlator, but does allow for the overall structure of the line to be more precisely characterized.

It is important to note that while both $v_φ$ and $v_r$ are measured absolutely, $v_z$ is only a relative measurement and depends on the accuracy of the $v_{\text{LSR}}$ measurement. For both sources, we estimate the uncertainty on $v_{\text{LSR}}$ by fitting $v_0$ maps made from ¹³CO and C¹⁸O emission described in Appendix C. In both cases, the best-fit $v_{\text{LSR}}$ values vary by ~10 m s⁻¹, and this scatter is used as the uncertainty on

$v_{\text{LSR}}$. This uncertainty dominates the uncertainty plotted in Figures 5(g) and (h).

## 4.2. HD 163296

The azimuthally averaged velocity profiles traced by ¹²CO emission are shown in the left panel of Figure 5. We find that the HD 163296 velocity profiles are consistent with those previously published using the DSHARP data when adopting the same emission surfaces for ¹²CO (but noting the swap in
sign of $v_r$ to be consistent with the cylindrical coordinate definition; Teague et al. 2018a, 2019b).

There is a strong correspondence between variations in the gas velocity structure and the continuum gaps and rings. In particular, $v_φ$ will be modulated by changes in the gas pressure gradient due to perturbations in the disk physical structure. Such velocity deviations have been previously used to measure the depth of the gaps (and hence the mass of planets required to open these gaps; Teague et al. 2018a) and the efficiency of grain trapping in the continuum rings (Rosotti et al. 2020). Beyond the edge of the millimeter continuum there are far fewer variations in $v_φ$.

Both radial and vertical components of the velocity, $v_r$ and $v_z$, respectively, appear to have more structure in the outer disk. A radial outflow, previously reported by Teague et al. (2019b), is recovered and suggests that $^{12}$CO may be tracing the base of a wind in the outer disk. Intriguingly, HD 163296 hosts jets and a large-scale disk wind (Klaassen et al. 2013), which may be associated with these smaller-scale dynamical features. Booth et al. (2021) present a detailed analysis of the large-scale disk wind as seen with MAPS data. In that paper, Booth et al. suggest a wind launch point at ~4 au; in this case it is not clear as to the relationship between the large-scale wind and flow arising from the outer disk. The most striking kinematic feature in the outer disk is the large negative $v_z$ component at $r \approx 2''3$ (232 au) associated with meridional flows around an embedded planet (Pinte et al. 2018b; Teague et al. 2019b), believed to be driving the “kink” KPS discussed in Section 3.2.

### 4.3. MWC 480

Velocity profiles for the MWC 480 disk are shown in the right column of Figure 5 and show a different behavior from those for the HD 163296 disk. Despite the rich substructures observed in the dust continuum for this source (Long et al. 2018; Liu et al. 2019; Sierra et al. 2021), there appears to be a much weaker relation between where velocity perturbations are found in the disk and the continuum substructure; however, the large uncertainties preclude any definitive statement. In particular, the inner $1''2$ (194 au) appears to be dominated by azimuthally asymmetric features, very clearly seen in Figure 4, which result in a poor constraint of $v_φ$ and the resulting large uncertainties on both $v_φ$ and $v_r$. In the outer disk, however, there are significant deviations around ~1"5, the same radius as the large arc shown in Figure 4, in addition to potentially related structure extending out to the edge of the disk, with the radial width of the features increasing with radius.

The disk of MWC 480 appears to have negligible $v_r$ components that are azimuthally symmetric, suggesting the absence of a disk wind like HD 163296. However, this does not rule out more localized radial flows, particularly in the inner disk, where the uncertainties are large. The largest velocity deviation is seen in the $v_z$ component, tracing a large vertical flow at $r \approx 1''5$, seen clearly as the ring in Figure 4. Smaller vertical flows are also seen and correspond to the weak arcs observed in the $v_0$ map residuals.

## 5. Gas Temperature

For optically thick emission, such as $^{12}$CO, we can use the brightness temperature, $T_B$, of the emission as a proxy for the gas temperature. This has been previously suggested as a method to detect embedded planets via their heating of their

immediate surroundings (Cleeves et al. 2015), with recent studies finding spiral structure likely associated with possible planets (e.g., Teague et al. 2019a; Wölfer et al. 2021). In this section, we look for variations in the temperature structures probed by CO isotopologue emission that correlate with the velocity structures we have detected in the disks of HD 163296 and MWC 480.

### 5.1. Method

To measure an azimuthally averaged $T_B$ profile for each disk, we use the $T_B$ maps produced by bettermoments when calculating the $v_0$ maps described in Section 3. Using GoFish, these maps were binned into annuli with a width of 1/4 of the beam major axis, taking into account the elevated emission surface of $^{12}$CO by using the analytical fits described in Law et al. (2021b) and then averaged. To estimate the uncertainty, the standard deviation of each annular bin was divided by the square root of the number of independent beams that fit in that bin. These were converted to units of kelvin using the full Planck law. The resulting profiles are shown in Figure 6, with the gaps in the continuum emission included. The inner $0''4$ of all profiles are affected by beam dilution and continuum subtraction, resulting in a drop in $T_B$. Beam dilution will occur at radial separations larger than the beam, as the projected Keplerian rotation leads to a narrow azimuthal extent of the emission that remains unresolved (Horne & Marsh 1986).

The radial profiles show a level of substructure that is not as clearly seen in the radial profiles of the integrated intensities for strong lines (aside from HCO$^+$; see Law et al. 2021a). This is because while integrated intensity and peak brightness are related for a single emission source, the backside of the disk will contribute to the integrated intensity in a spatially dependent manner, while the peak brightness temperature should only be affected by the front-side emission (at least for optically thick lines).

These azimuthally averaged profiles are used to create a 2D background model, which can be subtracted from the $T_B$ maps to reveal subtle perturbations in the background temperature structure (e.g., Teague et al. 2019b). Figures 7 and 8 show, for HD 163296 and MWC 480, respectively, the $T_B$ map in panel (a), the residual after subtracting the azimuthally averaged model $T_{B,avg}$ in panel (b), and the residuals deprojected to account for the inclination and projection of the emission surface in panel (c). For both sources, an enhanced $T_B$ is found along the major axis of the disk extending ~1" from the disk center. This arises owing to the Keplerian rotation pattern. Gas tracing smaller radii, which is therefore hotter and brighter, is detected primarily along the major axis at small radial offsets owing to the projection effects from the inclination of the disk (see, e.g., Figure 1 in Horne & Marsh 1986).

### 5.2. HD 163296

The radial $T_B$ profiles for HD 163296 are shown in the left column of Figure 6. Small deviations relative to the mean background are observed in all three isotopologues, likely associated with the continuum substructure in the inner $1''5$. However, the outer disk appears fairly featureless. The residuals of the $T_B$ profiles for the three CO isotopologues after subtracting a background model (a 10th-order polynomial fit, as done in Section 3.1) are shown in Figure 6(c). However, the radial variations appear to be associated with the ring and
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Azimuthally averaged $T_B$ profiles for CO $J = 2-1$ isotopologue emission in HD 163296 (panels (a) and (c)) and MWC 480 (panels (b) and (d)). In all panels, the shaded regions show the $1\sigma$ uncertainty, while the vertical dashed lines show the locations of gaps in the millimeter continuum. Panels (c) and (d) show a zoom-in of panels (a) and (b), respectively, after subtracting a smooth background model, $T_{B,\text{mod}}$, shown as a dotted black line in panels (a) and (b). The background model is a 10th-order polynomial fit to the profile between 0.1"3 and 4" for $^{12}\text{CO}$ and $^{13}\text{CO}$ and between 0.1"3 and 2" for $^{18}\text{O}$. The shaded regions in panels (a) and (b) show the extent of panels (c) and (d), respectively.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Gas temperature perturbations probed by $^{12}\text{CO}$ in HD 163296. Panel (a) shows the $T_B$ map, tracing the gas temperature of the top surface. Panel (b) shows the residual after subtracting the azimuthally averaged $T_B$ profile. Panel (c) shows the same residuals, but deprojected into a face-on orientation. A representative beam size is shown in the lower right corner; however, due to the elevated emission surface, the true beam size varies across the map. In panels (b) and (c), dashed contours show lines of constant radius and azimuth in steps of 0.1"5 and 30°. In all frames, the disk rotates in clockwise direction. The region where Pinte et al. (2018a) detected a KPS, labeled box A in Figure 2, is shown in both panels (b) and (c).</figcaption></figure>

gap structure and point toward variations in the thermal structure of the disk driven by changes in the dust opacity in gapped regions (e.g., Teague et al. 2017; Facchini et al. 2018; van der Marel et al. 2018; Alarcón et al. 2020; Rab et al. 2020). It is important to note that optically thick continuum emission has been shown to impact the emergent radiation, either through the absorption of emission arising from behind it or through continuum subtraction (e.g., Boehler et al. 2018; Weaver et al. 2018), potentially introducing spurious features in molecular line emission. This is unlikely the case here, as the same radial structure is seen when using data with and without

continuum subtraction. Furthermore, with the molecular emission arising from an elevated surface, the continuum emission does not trace the same annuli when projected on the sky, limiting the impact of the absorption.

There are subtle substructures seen in the $T_B$ residual map, shown in Figure 7(b). A tentative arc is seen along the southwestern quadrant at a radius of ~3.5". With the current sensitivity, it is hard to constrain its azimuthal extent; however, it appears to have an annular rather than a spiral morphology. There is only tentative structure observed around the location of the planet, at (0.1"0, 1.1"5) in the middle panel and (-1.1"5,
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Same as Figure 7, but for $^{12}$CO in the disk of MWC 480. In all frames, the disk rotation is in an anticlockwise direction.</figcaption></figure>

2′′5) in the projected panel. The “kink” KPS feature seen in Figure 1(b) is bounded to the east (thus trailing the hypothesized planet) by slightly hotter gas, while at the location of the planet we find a slightly lower gas temperature. Deeper observations, facilitating a comparable mapping of the $^{13}$CO and $^{18}$O emission, are required to fully map the temperature structure around the planet and thus quantify the significance of any localized heating from accretion onto the planet.

### 5.3. MWC 480

In contrast, MWC 480 displays a more complex range of temperature profiles with ringed features seen across the inner 3" of the disk, as shown in Figure 6(b). The substructures seen in the $T_B$ profiles, Figure 6(d), appear not to strongly correlate with the location of the continuum substructures, with the exception of the large dip in $^{12}$CO seen at 245 au. Of particular note is that the features traced by different isotopologues are not merely offset in radius, which could be explained by an incorrect emission surface used in the projection, but they appear to have different radial scale wavelengths such that their peak-to-peak distances increase with height traced in the disk. This phenomenon is discussed more in Section 6.2.

MWC 480 displays tentative substructure when looking at the full 2D map, characterized by multiple, tightly wound spirals. The top panel of Figure 9 shows an annotated version of Figure 8(c), with the spiral arms shown by colored dotted lines. The spirals are only identified by eye, with the quality of the data precluding a more quantitative fit. Deeper observations are required to more accurately map these subtle $T_B$ variations. Figure 8(b) shows a polar deprojection of this panel. These spiral features align with the two spirals observed in the gas velocity (see Figure 4), where the temperature perturbations are found to be both leading and trailing the velocity perturbations. We find no analogous perturbation in the velocity structure to the outer spirals. The pitch angles of the spirals, $\beta$, are defined as $\tan(\beta) = -dr/rd\phi$ and plotted in the bottom panel of Figure 9, with the gaps in the millimeter continuum shown for reference. The significance of the radial dependence on $\beta$ is discussed in Section 6.2.2.

## 6. Discussion

We have shown that the disks around both HD 163926 and MWC 480 exhibit a substantial amount of substructure in the gas traced by $^{12}$CO emission, in both the gas velocity structure and temperature structure. These structures are predominantly at radii beyond the continuum edge of the disks, demonstrating either that the processes that are sculpting the large grains traced by (sub)millimeter continuum emission also influence the outer disk or that there are additional processes present in the outer disk. In this section, we interpret the observed substructure as the signatures of embedded planets.

### 6.1. HD 163296

With these data, we are able to confirm many of the previously detected kinematic signatures of embedded planets in HD 163296. In particular, we recover the “kink” reported by Pinte et al. (2018b) at the outer edge of the disk, finding that it is much more azimuthally extended than previously inferred. In addition, we find azimuthally averaged velocity profiles that agree with those previously reported (Teague et al. 2018a, 2019a), which showed pressure-modulated rotational velocities due to perturbations in the gas surface density (Rosotti et al. 2020) and localized radially converging flows, interpreted as the tops of meridional flows (Morbidelli et al. 2014; Szulágyi et al. 2014).

#### 6.1.1. A Planet at 260 au

With observations at higher spectral resolution and achieving a higher sensitivity, the KPS previously detected (“kink” in the emission morphology) by Pinte et al. (2018b) is clearly seen in the rotation map residuals (Figure 2, feature A) and associated with a much more azimuthally extended structure. This is also reflected in the range of velocities (i.e., number of channels) the feature appears in channel maps, shown in Appendix D, where the perturbation is seen to change morphology as a function of velocity.

As these residuals extended over large azimuthal ranges, it is likely that at large azimuthal extents from the planet there is little azimuthal dependence on their projection and thus that the velocities in these regions are dominated by vertical flows (Equation (2), and see the discussion in Teague et al. 2019b).
<figure><img src="image_12.png" /></figure>

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

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Panel (a) shows an annotated version of Figure 8(c), highlighting the spirals, show with colored dotted lines. Panel (b) is a polar deprojection of panel (a) to better demonstrate the spiral structure. The pitch angles for the spirals are shown in panel (c), with the locations of the continuum gaps annotated.</figcaption></figure>

Although Pinte et al. (2018b) only originally considered perturbations in $v_\phi$ when modeling HD 163296, Pinte et al. (2019) reported more numerical simulations that showed that embedded planets drive perturbations in all three directions, with vertical motions starting to become present at higher altitudes, while radial and rotational perturbations are stronger

in midplane regions. Bae et al. (2021) also showed that when using more realistic thermal structures for the background disk, planets will excite buoyancy resonances, driving tightly wound spirals with large vertical motions dominating the dynamics. A potential scenario is that perturbations traced by $^{12}$CO should be dominated by in-plane motions, $v_\phi$ and $v_r$, close to the planet, and predominantly in the vertical direction at larger azimuthal offsets from the planet.

However, Bae et al. (2021) also predict that for buoyancy resonances the velocity deviations relative to the background Keplerian rotation will be largest at the location of the planet. This is consistent with observed velocity deviations and the predicted location of the planet proposed by Pinte et al. (2018b). For HD 163296, unfortunately this location is also where the projected components of $v_\phi$ and $v_r$ can add coherently ($\sin(\phi) \approx \cos(\phi)$, so $\phi \approx 45^\circ$ or $225^\circ$; see Equation (2)), meaning that we have the greatest sensitivity to these motions at this particular location. Detecting how the velocity perturbations vary as a function of height and comparing with numerical-simulation-based predictions may help distinguish between spiral launching scenarios. However, the only true way to distinguish scenarios will require the orbital rotation of the planet such that the projection of the velocity components changes significantly. At 230 au the orbital velocity around a $2 M_\odot$ star is $\sim 5.8$ mas yr$^{-1}$, meaning that with current instrumentation it will be around a decade until a discernible difference can be detected (the limiting factor here being how well a peak or point of reference can be determined in each image).

A planet can readily account for the observed substructure in the disk of HD 163296, without the need for the fine-tuning of any of the disk properties (e.g., Pinte et al. 2018b; Teague et al. 2018a). A full comparison with numerical simulations of planet–disk interactions will be the focus of a future article. The inference of an embedded planet associated with meridional flows represents a unique opportunity to probe the delivery of volatile materials to the still-forming atmosphere. Future observations that achieve a much greater sensitivity are required if we are to compare the chemical complexity observed in the inner regions of the disk to that in the immediate vicinity of wide-separation, embedded planets.

### 6.1.2. Additional Localized Velocity Deviations

A Second Outer Planet?—In Section 3.2, we showed the presence of a localized velocity deviation broadly coincident in radius with the D145 gap in continuum emission (feature B in Figure 2). Applying the same method of $v_0$ map making and subtracting a Keplerian velocity model from the $^{13}$CO (2–1) data, we find a similar velocity deviation at the same location, as shown in Figure 1. The lower signal-to-noise ratios of the $^{13}$CO images result in a residual map that is far nosier, but the coincidence in location with the $^{12}$CO feature strongly suggests that this is a real perturbation and has a vertical extent that spans between the layers probed by $^{12}$CO and $^{13}$CO emission: $z/r \approx 0.3$ and $0.15$, respectively (Law et al. 2021b). It is suspicious, however, that this feature shares the same azimuthal angle as that of the planet at 260 au, and it could suggest that a large planet can generate kinematic features across a large radial extent. Simulations of planet–disk interactions employing the most realistic thermal and density structures will be essential in distinguishing between a second planet, or perturbations from a single planet influencing a much more
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Channel maps of $^{12}$CO J = 2 − 1 emission from the disk around HD 163296, centered on the proposed kinematic planetary signature reported by Pinte et al. (2020). Panels (a) and (b) used DSHARP data (Andrews et al. 2018; Isella et al. 2018), while panels (c)–(f) use MAPS data. No clear “kink” is visible at the spatial resolution of the MAPS data. The dotted circle shows the region of interest.</figcaption></figure>

extended radial range than previously thought (e.g., Dong et al.
2017; Bae & Zhu 2018a, 2018b).

*A Secondary KPS in the Inner Disk.*—Pinte et al. (2020) identified a second “kink” KPS in the DSHARP set of images. The authors argued that if this feature was driven by an embedded planet, it would require the planet to sit within the D86 gap. Figure 10 shows zoomed-in channel maps of the $^{12}$CO emission at the same location as the proposed velocity deviation and in comparison with the DSHARP data (Andrews et al. 2018; Isella et al. 2018). As the spatial resolution of the MAPS data is slightly coarser than that achieved by the DSHARP observations, we are unable to verify the feature in the channel maps. However, residuals in the $^{12}$CO $v_0$ maps made with both the MAPS and DSHARP data are found at the proposed location of the planet at ≈86 au, as seen in Figure 11(b), in addition to far more tentative features seen in the $^{13}$CO and $C^{18}$O $v_0$ maps, shown in Figures 11(c) and (d). A blue/red feature is seen around the proposed planet location, reminiscent of the “Doppler flip” found in HD 100546 (Casassus and Pérez 2019), but the orientation of the feature appears inconsistent with the clockwise rotation of the disk.

This “kink” KPS also lies very close to the azimuthal
asymmetry detected in the continuum emission at r ≈ 0.54
(55 au, Figure 11(a) and feature C in Figure 4; Isella et al.
2018). It is unknown what would cause such a feature;
however, it would require a local pressure maximum to radially
and azimuthally capture the large grains. Recent works looking

at the kinematics around the vortices, particularly in
HD 142527, have suggested that such a pressure maximum
has observable effects on the local dynamical structure of the
disk (Huang et al. 2018b). For the dust trap in HD 142527, this
manifested as super-Keplerian rotation along the outer edge of
the feature (a positive v₀ residual; Yen & Gu 2020; Boehler
et al. 2021), consistent with the velocity residuals seen here.

The features in the $v_0$ map residuals support the hypothesis of localized velocity deviations. However, given the spatial resolution, it is hard to distinguish between the source of the velocity deviations and those potentially arising from radiative transfer effects due to the large intensity gradients associated with gaps in the continuum (e.g., Keppler et al. 2019; Rab et al. 2020; Boehler et al. 2021). Follow-up observations achieving spatial resolutions of <0.1 are absolutely required if we are to probe the velocity structure of gas on spatial scales equivalent to the substructure observed in the dust continuum.

6.2. MWC 480

MWC 480 exhibits an extensive variety of substructures as
traced by the $^{12}$CO emission, in both the velocity and the
temperature. Within the inner ~1.2, there are substantial
velocity perturbations with narrow azimuthal extents that
preclude a reasonable inference of the $v_\phi$ and $v_r$ radial profiles,
as evidenced by the large uncertainties for these profiles, shown
in Figure 5. Outside this radius, there are significant deviations
in $v_\phi$ and $v_z$, while perturbations in the radial velocity appear to
be minimal. In turn, these velocity perturbations are associated
with large-scale variations in the gas temperature, a combina-
tion of annular structures that have variations on the order of
≈3 K (≈5% relative to the background), with spiral perturba-
tions detected on top of these, also with variations around
≈3 K. As the disk rotates anticlockwise on the sky, these
spirals are consistent with trailing spirals.

6.2.1. *Ordered Substructure*

What is particularly striking about the features in MWC 480
is how ordered the deviations are: features appear incredibly
regular as a function of radius, particularly when compared to
the structure observed in the disk of HD 163296 (see Figures 5
and 6). To measure the size of these features, we use a
Savitzky–Golay filter (with a window size of 0.45 and a third-
order polynomial fit) to extract the second derivative of the $v_\phi$
and $T_B$ profiles. This approach was found to yield a good
compromise between smoothing the data to suppress noise and
obtaining an accurate measure of radial location of features in
the profile. We identify peaks and troughs in the profiles as
local minima and maxima in the second-derivative profiles,
respectively. The radial width of these features, either the peak-
to-peak or trough-to-trough distance, and their radial location
are plotted in Figure 12. A linear fit to each of these is shown
by the dashed lines, finding a radial dependence of ~0.39r for
the gas velocity and ~0.15r for the gas temperature.

It is well known that hydrodynamical processes that create
structure in disks will result in features that scale with the local
pressure scale height, i.e., increasing with radius. This is
routinely seen where the width of gaps and rings formed by a
planet with a given mass scales with the local scale height (e.g.,
Kanagawa et al. 2015a; Zhang et al. 2018; Yun et al. 2019),
which in turn increases with radius. The strong correlation
shown in Figure 12 between the size of features detected in
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Kinematic structures around the continuum vortex and proposed KPS in the inner disk of HD 163296. Panel (a) shows the spatially deprojected 1.25 mm continuum data from DSHARP (Andrews et al. 2018; Isella et al. 2018) with the vortex and KPS region, box C in Figure 2, outlined. Panels (b), (c), and (d) show similar zoom-ins of the $v_0$ residuals for $^{12}$CO, $^{13}$CO, and C$^{18}$O, respectively. In all panels, the proposed planet location from Pinte et al. (2020) is plotted. All $v_0$ residual maps share the same color scaling. A representative beam size is shown in the lower left corner of each panel, taking into account the deprojection.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 12.** Radial dependence on the size of gaseous substructures in the disk of MWC 480. The figure compares the radial size of the perturbations (maximum-to-maximum or minimum-to-minimum of the second derivative of the radial profile) as a function of radius. The dashed lines show a linear fit to each set of points.</figcaption></figure>

both the gas temperature and velocity structure and their orbital radius suggests that there appears to be a dominant hydro-dynamical process sculpting the structure of the disk.

This is counter to the scenario observed in HD 163296, where the gaseous substructures appear to be somewhat independent of one another. A possible cause could be the difference in the number of planets in each disk. HD 163296 is believed to host (at least) three Jupiter-mass planets (inferred through gas and dust surface density depletions, Isella et al. 2016; and through gas dynamics, Pinte et al. 2018b; Teague et al. 2018a), such that the influence of each planet will be felt most strongly in radial regions close to the perturber. If MWC 480 were to only host a single planet, a scenario explored in the next subsection, then the observed deviations would be characterized only by the single perturber.

### 6.2.2. A Possible Planet at 245 au

An embedded planet would provide an enticing mechanism to explain the observed substructures in the gas and dust. To summarize, the main observed features in the gas are as follows:

1. Two tightly wound spiral structures seen in the $v_0$ residuals at 162 au (1") and 245 au (1".5). See Figure 4.

<figure><img src="image_18.png" /><figcaption>**Figure 13.** Azimuthally averaged 242 GHz continuum $T_B$ profile for MWC 480 (Sierra et al. 2021). The shaded region shows the $1\sigma$ uncertainty, estimated by the standard deviation in each radial bin. Three gap centers, D73, D143, and D245, are marked by orange squares, and two rings, B95 and B163, are marked by orange circles. The locations of the two annular structures in $T_B$ are marked by $\delta T_B$ (see Figure 9), while the two dominant spirals in $v_\phi$ (see Figure 4) are marked as $\delta v_\phi$.</figcaption></figure>

2. Rings of large ($\approx 50 \text{ m s}^{-1}$) vertical velocities, with the most prominent at 245 au (1".5; feature A in Figure 4). See Figures 4 and 5.

3. Tightly wound trailing spirals in the gas temperatures, with all pitch angles increasing with smaller radii. See Figure 9.

4. A series of temperature minima and maxima extending across the radius of the disk, with the radial width of the features increasing with radius. See Figure 12.

In addition to the complex structures in the gas, there are gaps detected in the millimeter continuum distribution (Law et al. 2021a; Sierra et al. 2021), as shown in Figure 13. Extending the list of features, we include the following:

1. Two gaps at 73 au (0".46) and 141 au (0".87), with the latter being a broader structure (Long et al. 2018; Law et al. 2021a). See Figure 13.

In the following, we identify the most likely location for a planet able to account for these observed features. Without the use of numerical simulations and radiative transfer modeling,
we cannot say whether the planet can drive all the observed features, or only a subset of them; however, that is the focus of a future paper. The aim of this section is to qualitatively understand which features are consistent with a planet, and where that planet would need to be located.

*Gaps in the Continuum*.—Any planet of sufficient mass to drive the substructure observed in the gas would also be massive enough to open a gap in the continuum (e.g., Papaloizou & Lin 1984; Paardekooper & Mellema 2004). This suggests that the planet could sit within one of the two continuum gaps, 73 au or 143 au, or outside the millimeter continuum at ≳290 au. There is a tentative dip in the continuum at 243 au (1"5), as shown in Figure 13; however, deeper and higher angular resolution observations are needed to verify the presence of this feature. Note that this feature is not reported in MAPS III (Law et al. 2021a) owing to the feature lying outside the radius that encloses 90% of the disk flux, which was taken as the disk edge. To distinguish between these potential locations, we must look to the gaseous substructures.

*Spirals*.—There are two mechanisms that can launch spiral arms from an embedded planet: the traditional Lindblad spirals (Goldreich & Tremaine 1979; Ogilvie & Lubow 2002), and the more recently discussed spirals arising from buoyancy resonances (Zhu et al. 2012; Bae et al. 2021). The former drives spirals that have a maximum pitch angle at the location of the planet, that is, at larger radial separations from the planet, the spirals become more tightly wound. In addition, these are usually large pitch angles, reaching ≳20° at the location of the planet (Zhu et al. 2015; Bae & Zhu 2018a, 2018b). In contrast, spirals driven by buoyancy resonances are much more tightly wound, with pitch angles monotonically decreasing with orbital radius and no significant changes around the location of the planet (Bae et al. 2021). A second distinction between these two mechanisms is the location at which the perturbations peak. For the case of Lindblad spirals, the most significant deviations from a smooth background are found at a radial distance of a few scale heights from the planet, while for spirals driven by buoyancy resonances the strength of the deviations is peaking at the location of the planet (Zhu et al. 2015; Bae & Zhu 2018a, 2018b; Bae et al. 2021).

If we assume the spirals to be predominantly Lindblad spirals, then the embedded planet needs to be located inside the observed spirals in order to account for the radial trend in $\beta$. This places the planet within the millimeter continuum, and thus the two gaps at 73 and 141 au prove to be likely locations. Alternatively, if the spirals are dominated by buoyancy resonances, then the launching radius cannot be constrained by the radial dependence of $\beta$, but rather on the strength of the perturbations. As the spirals are only tentatively detected, it is hard to isolate the most significant perturbation; however, it seems that the inner edge of the spirals, ~220 au (1"36), displays the largest azimuthal deviations after the subtraction of an azimuthally averaged background model.

We additionally observe two tightly wound spirals in the gas velocity structure at radii of 162 and 243 au (1"0 and 1"5), at a position angle that connects them to the large variations seen in the gas temperature. These two features have extremely small pitch angles, $\beta \lesssim 5°$, and constitute the largest perturbations in the gas velocity. It is therefore more likely that these are associated with buoyancy resonances than Lindblad spirals. In this scenario, the planet responsible for driving these features would be in the immediate vicinity, $\lesssim 0''2$, of the perturbations.

*Changes in Rotation Velocity*.—In addition to opening a gap in the distribution of large grains, a planet will perturb the local gas pressure profile, primarily by decreasing the total gas surface density, but also through changes in the local thermal structure (e.g., Teague et al. 2017; Rab et al. 2020). Pressure perturbations are detectable in the $v_\phi$ profile, with pressure minima characterized by a slowing of the rotation inside of the pressure minimum and a hastening of the rotation outside (Kanagawa et al. 2015b; Teague et al. 2018a). Using Figure 5(d) as a guide, we see that pressure minima are qualitatively consistent with the deviations in the $v_\phi$ profile at 143, 187, 245, and 339 au (0"88, 1"15, 1"52, and 2"10, respectively). We note that changes in emission height can also drive changes in $v_\phi$ (see Equation (1)); however, Law et al. (2021b) found no significant sharp changes in the emission height of $^{12}$CO (2–1) in MWC 480, aside from a depression at 66 au (0"41). Law et al. (2021b) report a tentative wave-like structure in the emission surface between 100 and 180 au; however, the locations of the proposed peaks and troughs in the emission surface do not align with the features reported here.

*Substructure in Molecular Emission*.—At these locations of pressure minima, it is unclear whether either the integrated emission (as traced by zeroth-moment maps; Law et al. 2021a) or the local gas temperature (probed by the brightness temperature for optically thick lines) will increase or decrease. Extensive thermochemical modeling has shown that these changes are highly dependent on not only the magnitude of the perturbation in the gas surface density but also the level of depletion of small grains in the gap, which are the dominant coolant for the gas, and any shadowing effects from puffed-up rims of the gap (e.g., Teague et al. 2017; Facchini et al. 2018; van der Marel et al. 2018; Alarcón et al. 2020; Rab et al. 2020). In general, however, there is a growing consensus that for shallow perturbations in the gas surface density, for example, those opened by a Saturn-mass planet, the local gas temperature will increase owing to the drop in opacity toward UV photons, locally enhancing the CO abundance and thus boosting the CO emission (van der Marel et al. 2018; Alarcón et al. 2020). Conversely, for more significant gas gaps, for example, driven by a Jupiter-mass planet, the drop in total gas column will dominate all other effects and the CO emission will drop, due to the combination of lower CO abundances and the CO emission tracing a cooler region, closer to the disk midplane (Rab et al. 2020). Readers are referred to MAPS V (Zhang et al. 2021) for a thorough analysis of the radial profiles of CO column densities in these sources.

As such, we remain agnostic about the expected change in temperature and instead correlate the locations of pressure minima with the local maxima or minima in the $T_B$ profile (Figure 6(d)). From these locations, it is only the 245 au (1"52) location that closely matches a pressure minimum and coincides with an extreme of the $T_B$ profile—here a minimum. Furthermore, this location is the only location that broadly matches a feature in the integrated intensity, D246, a dip at 246 au (1"55). As both of these features are a depression, this would suggest that a higher-mass planet, ~1 $M_{Jup}$, would be needed to drive such a perturbation.

Although there is no definitive location for a planet, it appears that the observed features are all broadly consistent with being launched by a planet at ≈245 au (1"52). Such a planet would open a gap in both the dust continuum, tentatively seen in Figure 13, and the gas, consistent with the pressure
minimum interpretation of the $v_φ$ variations, and drop in both the $^{12}$CO brightness temperature profile and integrated intensity at ≈245 au (the integrated intensity map is presented in Law et al. 2021a, and shown to be azimuthally symmetric, with radial variations on the order a few percent relative to the background value being found, similar to the peak brightness profile). Based on the location and strength of the spiral perturbations, it is likely that a planet would be at a position angle between 270° and 0°, i.e., in the northwest quarter of the sky.

While there is no clear “kink” KPS detected in the CO isotopologue emission at this location (see, e.g., Figure 19), the channel maps are dominated by “spurs,” features that Bae et al. (2021) argued were spirals from buoyancy resonances, which may hide a KPS. It is unclear why these spur features dominate in MWC 480 and could be due either to the lower inclination of the disk relative to that of HD 163296, thus increasing the sensitivity of our observations to the vertical motions, or to a steeper vertical temperature gradient, more favorable for the generation of these resonances, as suggested by the thermal structured extracted from the data and described in MAPS IV (Law et al. 2021b). High angular resolution observations able to spatially resolve these features may facilitate the identification of a KPS within these spurs around the location of the proposed planet.

## 7. Summary

In this work, we have analyzed high spatial and spectral resolution observations of CO isotopologue emission from the disks of HD 163296 and MWC 480. Both disks exhibit a stunning variety of gaseous substructures, traced by local variations in the gas temperature and velocity structure, which we summarize below.

A previously detected “kink” KPS in the $^{12}$CO channel maps of HD 163296 was shown to be much more azimuthally extended than originally thought, with velocity deviations observed to span at least half the disk in azimuth. These features tentatively extend the full azimuth of the disk, connecting to a spiral in the southern side of the disk. Unfortunately, deviations due to the large pressure gradient at the outer edge of the disk prevent a more definitive statement.

At smaller radii in the disk of HD 163296, substantial velocity deviations on the order of ~200 m s⁻¹ were observed in both $^{12}$CO and $^{13}$CO emission. With the magnitude of the perturbations peaking at a radius of 145 au, consistent with a gap in the continuum emission, a planetary origin is an enticing prospect. The coincidence in the azimuthal angle and that of the aforementioned planet at 260 au, however, may indicate that these deviations are related and will require follow-up observations to disentangle.

In addition to HD 163296, the gas disk of MWC 480 displays an exceptional level of substructure. $^{12}$CO, $^{13}$CO, and C$^{18}$O emission were found to show concentric radial variations in their $T_B$ profiles of ≈3 K, with the radial size of the features growing with radius. On top of this complex background, subtle spiral perturbations were detected in the $^{12}$CO emission at the level of ≈3 K. Linear fits to the spirals suggested small pitch angles (between 2° and 6°), broadly consistent with predictions for spirals driven by buoyancy resonances. Two large arc-shaped perturbations in the gas velocity structure were found to connect to the spiral arms traced by $T_B$, in addition to a significant vertical velocity perturbation centered at 245 au. A

planet at ≈245 au could be responsible for driving the observed perturbations, in addition to a tentative gap at 245 au in the dust continuum.

These observations represent some of the deepest and most well-resolved images of CO emission from a protoplanetary disk to date and highlight the utility of deep, spectral line observations of protoplanetary disks. We would see significant benefits in an improved characterization of these features from future follow-up observations better suited to spatially and spectrally resolve them. These sources are also ideal targets for the upcoming James Webb Space Telescope, as the predicted planets lie at large separations from the central star, where contamination from the stellar point-spread function is minimized.

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01055.L and 2016.1.00484.L. 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.

Y.B. acknowledges funding from ANR (Agence Nationale de la Recherche) of France under contract No. ANR-16-CE31-0013 (Planet-Forming-Disks). A.S.B. acknowledges the studentship funded by the Science and Technology Facilities Council of the United Kingdom (STFC). V.V.G. acknowledges support from FONDECYT Iniciación 11180904 and ANID project Basal AFB-170002. K.Z. acknowledges the support of the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin–Madison with funding from the Wisconsin Alumni Research Foundation. K.Z., K.R.S., J.H., J.B., J.B.B., and I.C. acknowledge the support of NASA through Hubble Fellowship grants HST-HF2-51401.001, HST-HF2-51419.001, HST-HF2-51460.001-A, HST-HF2-51427.001-A, HST-HF2-51429.001-A, and HST-HF2-51405.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. C.J.L. acknowledges funding from the National Science Foundation Graduate Research Fellowship under grant DGE1745303. A.D.B. and E.A.B. acknowledge support from NSF AAG grant No. 1907653. K.I.Ö acknowledges support from the Simons Foundation (SCOL No. 321183) and an NSF AAG Grant (No. 1907653). R.L.G. acknowledges support from a CNES fellowship grant. S.M.A. and J.H. acknowledge funding support from the National Aeronautics and Space Administration under grant No. 17-XRP17 2-0012 issued through the Exoplanets Research Program. J.D.I. acknowledges support from the Science and Technology Facilities Council of the United Kingdom (STFC) under ST/T000287/1. C.W. acknowledges financial support from the University of Leeds, STFC, and UKRI (grant Nos. ST/R000549/1, ST/T000287/1, MR/T040726/1). R.T. and F.L. acknowledge support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow. Y.A. and G.C. acknowledge support by NAOJ ALMA Scientific Research grant code 2019-13B and Grant-in-Aid for Scientific Research Nos. 18H05222 and 20H05847. A.S. acknowledges support from ANID/CONICYT Programa de Astronomía Fondo ALMA-CONICYT 2018
31180052. F.M. acknowledges support from ANR of France under contract ANR-16-CE31-0013 (Planet-Forming-Disks) and ANR-15-IDEX-02 (through CDP “Origins of Life”). Y.Y. is supported by IGPEES, WINGS Program, the University of Tokyo. L.M.P. acknowledges support from ANID project Basal AFB-170002 and from ANID FONDECYT Iniciación project No. 11181068.

*Facility:* ALMA.

*Software:* eddy (Teague 2019a) emcee (Foreman-Mackey et al. 2019) GoFish (Teague 2019b) matplotlib (Hunter 2007) scipy (Virtanen et al. 2020).

## Appendix A
### Creation of Rotation Maps

Prior to collapsing the data into a rotation map, it is often useful to apply a smoothing kernel along the spectral axis to remove noisy peaks. This can remove high-frequency noise and allow for a more accurate determination of the desired statistic; however, it also can result in subtle deviations being washed out. The Python package `bettermoments` (Teague & Foreman-Mackey 2018) offers two ways to smooth the data prior to collapsing the cube: convolution with a top-hat kernel of a user-defined size, or the use of the Savitzky-Golay filter with a user-defined window size and polynomial order. The Savitzky-Golay filter is frequently used in spectroscopy, as it

retains much of the *shape* of the profile, which can sometimes be lost with other forms of convolution.

In addition to spectrally smoothing the data when making the $v_0$ map, a convolution in the spatial plane after subtracting the model velocity map will remove some noise, aiding in the identification of coherent structures. One major source of structure in these residuals is the feathering due to the channelization of the data.

In the main text, we chose to use no smoothing prior to the collapse of the data, as it was found to wash out some of the subtler deviations. Figure 14 demonstrates the impact of the smoothing on the $v_0$ map of $^{12}$CO in HD 163296. The top left panel shows no spectral smoothing, the top right and bottom left panels show convolution with a top-hat kernel with a window size of two and three channels, respectively, while the bottom right panel shows the Savitzky-Golay filter with a window size of 5 and a second-order polynomial. All residual maps have been additionally smoothed with a 2D Gaussian convolution in the spatial plane to further suppress noise. This figure shows that the smoothing does not introduce any new features; however, it can highlight some of the more prominent features (e.g., the large red arc at a ~2'' offset directly north of the disk center), while simultaneously washing out some of the subtler detail (e.g., the azimuthal extension both leading and following the aforementioned arc).
<figure><img src="image_19.png" /><figcaption>**Figure 14.** Comparison of the HD 163296 <sup>12</sup>CO $v_0$ map residuals after subtracting a Keplerian rotation model when using different levels of spectral and spatial smoothing. Each row represents a different spectral smoothing: no smoothing, top row; top-hat convolutions, second and third rows; Savitzky-Golay filter, bottom row. Each column has a different level of spatial smoothing applied: no smoothing, left column; convolution with a Gaussian kernel with a standard deviation of 1 and 2 pixels, middle and right columns, respectively. The kernel size is plotted in the upper left corner of the appropriate panels.</figcaption></figure>

## Appendix B
### The Inclination of MWC 480

Within the MAPS collaboration, the inclination of MWC 480 was taken to be $i = 37^\circ$, based on the continuum fitting presented in Liu et al. (2019). However, when fitting the $v_0$ maps for MWC 480, as described in Appendix C, it was noticed that this fixed inclination resulted in a rotationally symmetric residual in the outer disk for all three CO

isotopologues considered, as shown in the right panel of Figure 15. Such residuals are indicative of a misspecified inclination.

To test this hypothesis, the same $v_0$ fitting procedure was performed for all three CO isotopologues, but letting the inclination vary. Fitting the entire disk yielded a larger inclination of $\approx 35^\circ$; however, this still resulted in systematic residuals seen in the outer disk. As the inclination was likely
<figure><img src="image_20.png" /><figcaption>**Figure 15.** Left: the $v_0$ map of $^{13}$CO (2–1) in MWC 480. The synthesized beam is shown by the ellipse in the lower left corner of the panel. Middle: residuals from a Keplerian model where the inclination was allowed to vary. The fit only considered emission outside a (projected) radius of 1.8, shown by the dashed black line. Right: residuals from a Keplerian model where the inclination was fixed to 37°.</figcaption></figure>

<table><caption>**Table 1**
Best-fit $v_{\text{kep}}$ Models from MAPS Data</caption><thead><tr><th rowspan="2">Model Parameter</th><th colspan="3">HD 163296</th><th colspan="3">MWC 480</th></tr><tr><th><sup>12</sup>CO (2–1)</th><th><sup>13</sup>CO (2–1)</th><th>C<sup>18</sup>O (2–1)</th><th><sup>12</sup>CO (2–1)</th><th><sup>13</sup>CO (2–1)</th><th>C<sup>18</sup>O (2–1)</th></tr></thead><tbody><tr><td>δx<sub>0</sub> (mas)</td><td>-20 ± 1</td><td>-16 ± 1</td><td>-11 ± 1</td><td>-19 ± 1</td><td>-13 ± 1</td><td>-9 ± 1</td></tr><tr><td>δy<sub>0</sub> (mas)</td><td>9 ± 1</td><td>10 ± 1</td><td>16 ± 1</td><td>5 ± 1</td><td>6 ± 1</td><td>10 ± 1</td></tr><tr><td>i (deg)</td><td>[46.7]</td><td>[46.7]</td><td>[46.7]</td><td>[-32.4]</td><td>[-32.4]</td><td>[-32.4]</td></tr><tr><td>PA (deg)</td><td>312.7 ± 0.1</td><td>312.7 ± 0.1</td><td>312.6 ± 0.1</td><td>327.9 ± 0.1</td><td>327.5 ± 0.1</td><td>327.5 ± 0.1</td></tr><tr><td>M<sub>*</sub> (M<sub>☉</sub>)</td><td>2.01 ± 0.01</td><td>1.93 ± 0.1</td><td>1.92 ± 0.1</td><td>2.61 ± 0.01</td><td>2.57 ± 0.01</td><td>2.60 ± 0.01</td></tr><tr><td>v<sub>LSR</sub> (km s<sup>-1</sup>)</td><td>5.76 ± 0.01</td><td>5.77 ± 0.01</td><td>5.75 ± 0.01</td><td>5.10 ± 0.01</td><td>5.10 ± 0.01</td><td>5.09 ± 0.01</td></tr><tr><td>z<sub>0</sub> (arcsec)</td><td>[0.39]</td><td>[0.12]</td><td>[0.17]</td><td>[0.30]</td><td>[0.11]</td><td>[0.07]</td></tr><tr><td>φ</td><td>[1.85]</td><td>[1.50]</td><td>[2.96]</td><td>[1.32]</td><td>[1.17]</td><td>[1.18]</td></tr><tr><td>r<sub>taper</sub> (arcsec)</td><td>[2.36]</td><td>[3.16]</td><td>[1.04]</td><td>[3.27]</td><td>[1.43]</td><td>[1.09]</td></tr><tr><td>q<sub>taper</sub></td><td>[1.18]</td><td>[5.00]</td><td>[4.99]</td><td>[2.73]</td><td>[4.11]</td><td>[3.11]</td></tr><tr><td>d (pc)</td><td>[101.0]</td><td>[101.0]</td><td>[101.0]</td><td>[161.8]</td><td>[161.8]</td><td>[161.8]</td></tr><tr><td>r<sub>fit,in</sub> (arcsec)</td><td>[0.28]</td><td>[0.29]</td><td>[0.29]</td><td>[0.33]</td><td>[0.34]</td><td>[0.34]</td></tr><tr><td>r<sub>fit,out</sub> (arcsec)</td><td>[3.75]</td><td>[2.25]</td><td>[2.25]</td><td>[4.25]</td><td>[3.00]</td><td>[3.00]</td></tr></tbody></table>

**Note.** Uncertainties represent the 16th–84th percentiles of the posterior distribution. Values in braces were held fixed during the fitting.

biased by the strong non-Keplerian features in the inner disk (see Section 3.3), this process was repeated, but fitting only the outer regions of the disk where there were no clear substructures, as shown in Figure 3. This resulted in fitting outside 2.5 for $^{12}$CO, outside 1.8 for $^{13}$CO, and outside 1.5 for C$^{18}$O. All three fits yielded an inclination of $i \approx 33^\circ$, while also removing the systematic residuals in the outer disk, as shown by the middle panel of Figure 15.

As the emission surfaces described in Law et al. (2021b) were based on an inclination of $37^\circ$, we repeat the fitting using the same method with an inclination of $i = 22^\circ 4$, with the results of the tapered power-law fit shown in Table 1. We find little difference between the emission surfaces reported in Law et al. (2021b); however, we find a slightly elevated $^{12}$CO emission surface owing to the change in projection from the inclination. For this paper, we adopt these emission surfaces. For a final check, the $v_0$ map fits are performed again, fitting only the outer edge of the disk, but using the updated emission surfaces. No significant change is found for the inclination from $\approx 33^\circ$, as the outer edge of the disk is close to being geometrically thin.

The continuum features discussed in Law et al. (2021a) were also checked to make sure the updated inclination did not significantly change their locations. Features were found to vary only by a few au, as the change in projected radius is small compared to the observed size of the feature.

## Appendix C: Rotation Map Fitting

In this appendix, we describe the rotation map fitting procedure with eddy and provide the best-fit results. eddy (Teague 2019a) is a wrapper for emcee (Foreman-Mackey et al. 2019), a Python ensemble sampler for affine-invariant Markov Chain Monte Carlo. In the modeling, five free parameters were considered, the source offset from the phase center, the disk position angle, the host star mass, and the systemic velocity, $\{\delta x_0, \delta y_0, \text{PA}, M_*, v_{\text{LSR}}\}$, while the disk inclination, $i$, and emission surface profiles, $z(r)$, were held constant. For each run, the inner $\approx 0.3$ was masked, as beam dilution limited the precision at which $v_0$ could be calculated. The outer fitting radius was dictated by the contamination from the rear side of the disk, with the adopted values provided in
Table 1. The uncertainty maps produced by bettermoments were adopted as the uncertainties.

To explore the posterior distributions of the free parameters, 32 walkers were used. The walkers took 2000 steps to burn in (it was found that walkers only needed ≈500 steps to converge), and then an additional 2000 steps to sample the posterior distribution. The posterior distributions for all parameters were found to be Gaussian, with little to no covariance between other parameters. The 16th to 84th percentile ranges were therefore used to estimate the statistical uncertainties on each model parameter, which are provided in Table 1. The $v_0$ maps for the CO isotopologues for HD 163296 and MWC 480 are presented in Figures 1 and 3, respectively. These include the residuals after subtracting the best Keplerian rotation model. Excellent agreement is found between model parameters for different CO isotopologues. Larger differences found in the dynamical masses can be attributed to the large kinematical deviations found within these disks, which will bias the fitting procedure.

We note that this procedure was performed for all sources within the MAPS sample—IM Lup, GM Aur, AS 209, HD 163296, and MWC 480—in order to derive dynamical masses for the generation of Keplerian masks (used for imaging

the data and extracted spectra). As all other MAPS papers adopted an inclination of 37° for MWC 480, the published literature value, the fit for MWC 480 was repeated with this larger inclination. This yielded a lower dynamical mass of $M_* = 2.1 M_\odot$, consistent with the dynamical mass reported in Simon et al. (2019).

## Appendix D: Channel Maps

Searching for coherent structures in the residuals when subtracting a rotation model from the line center maps is both quick and relatively straightforward to compare to maps of the dust continuum. As demonstrated in Disk Dynamics Collaboration et al. (2020), the same information is found in the channel maps,<fnref n="26" /> although presenting a harder challenge for interpretation. The reward to identifying the structures in a channel-map basis is that there is considerably more information compared to a rotation map, which has been first collapsed along the spectral axis. While this paper focused on the collapsed rotation maps, we provide the channel maps in this appendix. Figures 16–18 show the (2–1) transitions of ¹²CO, ¹³CO, and C¹⁸O in HD 163296, respectively. Figures 19–21 show the same molecular lines but for MWC 480.

<footnote n="26">²⁶ Indeed, the proposed planet in the disk of HD 163296 was originally identified in the channel maps rather than in the rotation map residuals (Pinte et al. 2018b).</footnote>
<figure><img src="image_21.png" /><figcaption>**Figure 16.** Channel maps of $^{12}$CO J = 2–1 emission in HD 163296. The color scaling has been chosen to highlight the morphology of the emission from the line wings. The velocity is given for each channel, as well as the synthesized beam shown in the lower left corner of each panel. Regions where the KPS and spiral are most visible are annotated.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 17. Same as Figure 16, but for $^{13}$CO J = 2 − 1 emission.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 18. Same as Figure 16, but for C$^{18}$O J = 2 − 1 emission.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 19. Channel maps of $^{12}$CO J = 2–1 emission in MWC 480. The color scaling has been chosen to highlight the morphology of the emission from the line wings. The velocity is given for each channel, as well as the synthesized beam shown in the lower left corner of each panel.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 20. Same as Figure 19, but for $^{13}$CO J = 2 − 1 emission.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 21. Same as Figure 19, but for C$^{18}$O J = 2 − 1 emission.</figcaption></figure>

ORCID iDs

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

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

Yuri Aikawa https://orcid.org/0000-0003-3283-6884

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

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

Jennifer B. Bergner https://orcid.org/0000-0002-8716-0482

Yann Boehler https://orcid.org/0000-0002-8692-8744

Alice S. Booth https://orcid.org/0000-0003-2014-2121

Arthur D. Bosman https://orcid.org/0000-0003-4001-3589

Gianni Cataldi https://orcid.org/0000-0002-2700-9676

Ian Czekala https://orcid.org/0000-0002-1483-8811

Viviana V. Guzmán https://orcid.org/0000-0003-4784-3040

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

John D. Ilee https://orcid.org/0000-0003-1008-1142

Charles J. Law https://orcid.org/0000-0003-1413-1776

Romane Le Gal https://orcid.org/0000-0003-1837-3772

Feng Long https://orcid.org/0000-0002-7607-719X

Ryan A. Loomis https://orcid.org/0000-0002-8932-1219

François Ménard https://orcid.org/0000-0002-1637-7393

Karin I. Öberg https://orcid.org/0000-0001-8798-1347

Laura M. Pérez https://orcid.org/0000-0002-1199-9564

Kamber R. Schwarz https://orcid.org/0000-0002-6429-9457

Anibal Sierra https://orcid.org/0000-0002-5991-8073

Catherine Walsh https://orcid.org/0000-0001-6078-786X

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

Yoshihide Yamato https://orcid.org/0000-0003-4099-6941

Ke Zhang https://orcid.org/0000-0002-0661-7517

References

Alarcón, F., Teague, R., Zhang, K., Bergin, E., & Barraza-Alfaro, M. 2020, ApJ, 905, 68

ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJL, 808, L3

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

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

Bae, J., & Zhu, Z. 2018a, ApJ, 859, 118

Bae, J., & Zhu, Z. 2018b, ApJ, 859, 119

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

Boehler, Y., Ricci, L., Weaver, E., et al. 2018, ApJ, 853, 162

Booth, A., Tabone, B., Ilee, J. D., et al. 2021, ApJS, 257, 16

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

Christiaens, V., Casassus, S., Perez, S., van der Plas, G., & Ménard, F. 2014, ApJL, 785, L12

Cleeves, L. I., Bergin, E. A., & Harries, T. J. 2015, ApJ, 807, 2

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110

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

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

Dong, R., Li, S., Chiang, E., & Li, H. 2017, ApJ, 843, 127
Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJL, 869, L46

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

Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104

Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131

Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68

Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, JOSS, 4, 1864

Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857

Horne, K., & Marsh, T. R. 1986, MNRAS, 218, 761

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018a, ApJL, 869, L42

Huang, J., Andrews, S. M., Öberg, K. I., et al. 2020, ApJ, 898, 140

Huang, J., Bergin, E. A., Öberg, K. I., et al. 2021, ApJS, 257, 19

Huang, J., Öberg, K. I., & Andrews, S. M. 2016, ApJL, 823, L18

Huang, P., Isella, A., Li, H., Li, S., & Ji, J. 2018b, ApJ, 867, 3

Hunter, J. D. 2007, CSE, 9, 90

Isella, A., Guidi, G., Testi, L., et al. 2016, PhRvL, 117, 251101

Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJL, 869, L49

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

Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2015a, ApJL, 806, L15

Kanagawa, K. D., Tanaka, H., Muto, T., Tanigawa, T., & Takeuchi, T. 2015b, MNRAS, 448, 994

Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118

Klaassen, P. D., Juhasz, A., Mathews, G. S., et al. 2013, A&A, 555, A73

Krajnović, D., Cappellari, M., de Zeeuw, P. T., & Copin, Y. 2006, MNRAS, 366, 787

Law, C., Loomis, R. A., Teague, R., et al. 2021a, ApJS, 257, 3

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

Liu, Y., Dipierro, G., Ragusa, E., et al. 2019, A&A, 622, A75

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17

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

Öberg, K., Bergin, E. A., Guzman, V. V., et al. 2021, ApJS, 257, 1

Ogilvie, G. I., & Lubow, S. H. 2002, MNRAS, 330, 950

Paardekooper, S. J., & Mellema, G. 2004, A&A, 425, L9

Papaloizou, J., & Lin, D. N. C. 1984, ApJ, 285, 818

Pérez, S., Casassus, S., & Benítez-Llambay, P. 2018, MNRAS, 480, L12

Pérez, S., Casassus, S., Hales, A., et al. 2020, ApJL, 889, L24

Pérez, S., Dunhill, A., Casassus, S., et al. 2015, ApJL, 811, L5

Pinte, C., Ménard, F., Duchêne, G., et al. 2018a, A&A, 609, A47

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

Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJL, 890, L9

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

Rab, C., Kamp, I., Dominik, C., et al. 2020, A&A, 642, A165

Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16

Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173

Sierra, A., Pérez, L. M., Guzman, V. V., et al. 2021, ApJS, 257, 14

Simon, M., Guilloteau, S., Beck, T. L., et al. 2019, ApJ, 884, 42

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

Tang, Y.-W., Guilloteau, S., Dutrey, A., et al. 2017, ApJ, 840, 32

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

Teague, R. 2019b, JOSS, 4, 1632

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

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., Bae, J., Huang, J., & Bergin, E. A. 2019b, ApJL, 884, L56

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

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

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

van der Marel, N., Williams, J. P., & Bruderer, S. 2018, ApJL, 867, L14

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

Weaver, E., Isella, A., & Boehler, Y. 2018, ApJ, 853, 113

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

Yun, H. G., Kim, W.-T., Bae, J., & Han, C. 2019, ApJ, 884, 142

Zhang, K., Booth, A., Law, C. J., et al. 2021, ApJS, 257, 5

Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJL, 869, L47

Zhu, Z., Dong, R., Stone, J. M., & Rafikov, R. R. 2015, ApJ, 813, 88

Zhu, Z., Stone, J. M., & Rafikov, R. R. 2012, ApJL, 758, L42