# Molecules with ALMA at Planet-forming Scales (MAPS): Complex Kinematics in the AS 209 Disk Induced by a Forming Planet and Disk Winds

Maria Galloway-Sprietsma<fnref n="1" />, Jaehan Bae<fnref n="1" />, Richard Teague<fnref n="2" />, Myriam Benisty<fnref n="3" />,<fnref n="4" />, Stefano Facchini<fnref n="5" />, Yuri Aikawa<fnref n="6" />, Felipe Alarcón<fnref n="7" />, Sean M. Andrews<fnref n="8" />, Edwin Bergin<fnref n="7" />, Gianni Cataldi<fnref n="9" />, L. Ilsedore Cleeves<fnref n="10" />, Ian Czekala<fnref n="11" />,<fnref n="12" />,<fnref n="13" />,<fnref n="14" />, Viviana V. Guzmán<fnref n="15" />,<fnref n="16" />, Jane Huang<fnref n="7" />,<fnref n="22" />, Charles J. Law<fnref n="8" />, Romane Le Gal<fnref n="4" />,<fnref n="17" />, Yao Liu<fnref n="18" />, Feng Long<fnref n="19" />,<fnref n="22" />, François Ménard<fnref n="20" />, Karin I. Öberg<fnref n="8" />, Catherine Walsh<fnref n="21" />, and David J. Wilner<fnref n="8" />

Received 2023 January 23; revised 2023 April 3; accepted 2023 April 4; published 2023 June 19

## Abstract

We study the kinematics of the AS 209 disk using the $J=2-1$ transitions of $^{12}$CO, $^{13}$CO, and C$^{18}$O. We derive the radial, azimuthal, and vertical velocity of the gas, taking into account the lowered emission surface near the annular gap at $\simeq 1''7$ (200 au) within which a candidate circumplanetary-disk-hosting planet has been reported previously. In $^{12}$CO and $^{13}$CO, we find a coherent upward flow arising from the gap. The upward gas flow is as fast as $150 \text{ m s}^{-1}$ in the regions traced by $^{12}$CO emission, which corresponds to about 50% of the local sound speed or 6% of the local Keplerian speed. Such an upward gas flow is difficult to reconcile with an embedded planet alone. Instead, we propose that magnetically driven winds via ambipolar diffusion are triggered by the low gas density within the planet-carved gap, dominating the kinematics of the gap region. We estimate the ambipolar Elsässer number, Am, using the $HCO^+$ column density as a proxy for ion density and find that Am is $\sim 0.1$ at the radial location of the upward flow. This value is broadly consistent with the value at which numerical simulations find that ambipolar diffusion drives strong winds. We hypothesize that the activation of magnetically driven winds in a planet-carved gap can control the growth of the embedded planet. We provide a scaling relationship that describes the wind-regulated terminal mass: adopting parameters relevant to 100 au from a solar-mass star, we find that the wind-regulated terminal mass is about one Jupiter mass, which may help explain the dearth of directly imaged super-Jovian-mass planets.

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

## 1. Introduction

Detecting exoplanets during their formation stages allows for a deeper understanding of planet formation processes. However, although there are more than 5000 confirmed exoplanets, only a few of them have been directly detected at a stage when they are still forming (Keppler et al. 2018; Haffert et al. 2019;

Currie et al. 2022). The Atacama Large Millimeter/submillimeter Array (ALMA) has revolutionized our ability to probe for young, forming planets. ALMA has revealed detailed substructures in continuum emission of protoplanetary disks, such as rings, gaps, and spirals (e.g., Andrews et al. 2018; Long et al. 2018; Cieza et al. 2021). These substructures provide compelling evidence that planets could be present in the disks, although we cannot rule out other origins (see reviews by Andrews 2020; Bae et al. 2022a).

In addition to continuum observations, by probing the kinematics of the protoplanetary disk gas via molecular line observations, ALMA provides a unique and powerful means to search for young planets. Molecular line observations are

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

<footnote n="1">¹Department of Astronomy, University of Florida, Gainesville, FL 32611, USA; mgallowayspriets@ufl.edu</footnote>

<footnote n="2">²Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="3">³Laboratoire Lagrange, Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, F-06304 Nice, France</footnote>

<footnote n="4">⁴Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), Université Grenoble Alpes, CNRS, F-38000 Grenoble, France</footnote>

<footnote n="5">⁵Università degli Studi di Milano, via Celoria 16, I-20133 Milano, Italy</footnote>

<footnote n="6">⁶Department of Astronomy, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan</footnote>

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

<footnote n="8">⁸Center for Astrophysics, Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA</footnote>

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

<footnote n="10">¹⁰Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA</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">¹⁵Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile</footnote>

<footnote n="16">¹⁶Núcleo Milenio de Formación Planetaria (NPF), Chile</footnote>

<footnote n="17">¹⁷Institut de Radioastronomie Millimétrique (IRAM), 300 rue de la piscine, F-38406 Saint-Martin d'Hères, France</footnote>

<footnote n="18">¹⁸Purple Mountain Observatory & Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, Nanjing 210023, People's Republic of China</footnote>

<footnote n="19">¹⁹Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="20">²⁰Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

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

<footnote n="22">²²NASA Hubble Fellowship Program Sagan Fellow.</footnote>
capable of discerning subtle localized kinematic perturbations, the so-called velocity kinks, caused by embedded planets (Perez et al. 2015; Pinte et al. 2018b, 2019, 2020). With observations of this nature, one can constrain the surface of the disk in different molecular tracers and therefore understand the 3D velocity structure of the disk. This method is particularly powerful because one can infer the location and mass of the planet (e.g., Izquierdo et al. 2021). Molecular line observations can also probe global-scale dynamics of the protoplanetary disk gas, such as radial changes of the gas velocity (Teague et al. 2018b, 2019a) and velocity variations along large-scale spirals (Teague et al. 2019b, 2021; Wölfer et al. 2022), which can be related to the perturbations created by yet-unseen planets. When multiple molecular lines probing different heights in a disk are used together, one can also probe coherent flows from the surface to the midplane (e.g., Yu et al. 2021; Teague et al. 2022). In addition, circumplanetary disks (CPDs) can be detected with molecular lines, providing unique and strong constraints on their physical and kinematic properties (Bae et al. 2022b).

Here we study the kinematics of the AS 209 protoplanetary disk using the $J=2-1$ transitions of $^{12}$CO, $^{13}$CO, and C$^{18}$O obtained as part of the ALMA Large Program Molecules with ALMA at Planet-forming Scales (MAPS; 2018.1.01055.L; Öberg et al. 2021). AS 209 is a 1–2 Myr old T Tauri star (Andrews et al. 2009, 2018) and is located 121 pc away in the Ophiuchus star-forming region (Gaia Collaboration et al. 2021). Previous continuum observations revealed multiple sets of concentric rings and gaps that extend out to ~140 au (Guzmán et al. 2018; Huang et al. 2018; Sierra et al. 2021), which are theorized to be caused by one or multiple giant planets (Fedele et al. 2018; Zhang et al. 2018). Molecular line observations also revealed rich annular substructures (Huang et al. 2016; Teague et al. 2018b; Law et al. 2021a). In particular, Teague et al. (2018b) kinematically identified a pressure minimum at ~1.9 (230 au) in $^{12}$CO, which was identified and spatially resolved previously by Guzmán et al. (2018). The previous work by Teague et al. (2018b) used the $^{12}$CO $J=2-1$ transition to measure the rotational velocity of the AS 209 disk and found deviations from Keplerian rotation. More recently, Bae et al. (2022b) reported a CPD candidate detected in $^{13}$CO $J=2-1$ emission, at a radial separation of 1.7 (200 au) from the star. With these gas substructures, along with a young, forming planet candidate in the disk, the AS 209 disk warrants a detailed study of its kinematics.

In this paper, we decompose the line-of-sight velocity into three orthogonal velocity components, namely radial, rotational (or azimuthal), and vertical velocities, for three CO isotopologues, $^{12}$CO, $^{13}$CO, and C$^{18}$O $J=2-1$. As we will show, this allows us to have a more complete 3D view of the kinematic structure of the disk.

This paper is organized as follows. We outline the observations in Section 2. In Section 3, we describe the analysis of the data, including the emission surfaces and the velocity profiles, and present the results. In Section 4, we discuss the results, focusing on the origin of the velocity structure in the AS 209 disk and its implications. We summarize our findings and discuss future directions in Section 5.

## 2. Observations

All data used in this work were obtained as part of the ALMA Large Program MAPS.<fnref n="23" /> For the observational setup and calibration process, we refer readers to Öberg et al. (2021). The imaging process is described in Czekala et al. (2021). As part of the MAPS data release, all images have been post-processed using the Jorsater & van Moorsel (1995; JvM) correction. For all analysis in this work, we use the robust = 0.5 weighted, JvM-corrected images.<fnref n="24" /> The synthesized beam size is 134 mas × 100 mas for $^{12}$CO $J=2-1$ with a position angle (PA) of 90.83, 140 mas × 104 mas for $^{13}$CO $J=2-1$ with a PA of 90.44, and 141 mas × 105 mas for C$^{18}$O $J=2-1$ with a PA of 91.37. The rms noise measured in a line-free channel is 0.562, 0.471, and 0.339 mJy beam$^{-1}$ for each data cube, respectively. The data were imaged with a channel spacing of 200 m s$^{-1}$, set by the MAPS Program.

In the top panels of Figure 1, we present peak brightness temperature maps for $^{12}$CO, $^{13}$CO, and C$^{18}$O $J=2-1$ lines calculated using bettermoments (Teague & Foreman-Mackey 2018a). The $^{12}$CO $J=2-1$ peak brightness temperature map clearly shows the annular gap at about 1.7 (≈200 au), which is the main feature we focus on in this paper. Additionally, the AS 209 disk suffers from foreground cloud contamination on the western side of the disk, visible in the $^{12}$CO brightness temperature map. Teague et al. (2018b) estimated that the cloud absorbs ~30% of the $^{12}$CO emission along the western side of the disk and showed that this level of perturbation does not impact the kinematic analyses (see their Appendix A.2).

## 3. Analysis and Results

### 3.1. Emission Surface and Disk Geometric Properties

To begin the characterization of disk kinematics, we first constrain the emission surface for $^{12}$CO, $^{13}$CO, and C$^{18}$O. For our base model, we adopt a power-law emission surface with an exponential taper, given by

$$z(r) = z_0 \left(\frac{r}{1''}\right)^\psi \exp\left(-\left[\frac{r}{r_t}\right]^{q_t}\right), \quad (1)$$

where $z(r)$ describes the height $z$ at a given radius $r$, $\psi$ is the power-law exponent, $r_t$ is the characteristic radius for the exponential taper, and $q_t$ is the exponent of the taper term, following Law et al. (2021b, 2022a). We note that Law et al. (2021b) already inferred the $^{12}$CO and $^{13}$CO emission surfaces of the AS 209 disk using the same data set as the one we use in this paper. However, Law et al. (2021b) limited the outer bound of the fit to 1.98 for $^{12}$CO, which does not cover the full radial extent of the $^{12}$CO disk (≈2.5), and to 1.35 for $^{13}$CO, which does not cover the gap around the CPD at ≈1.7. Because the main goal of this study is to study the kinematics within and around the gap, we opt to fit the emission surfaces adopting

<footnote n="23">²³ Data used for this project can be downloaded at the MAPS webpage: https://alma-maps.info/.</footnote>

<footnote n="24">²⁴ We repeated the analysis using data cubes with a 0.15 taper and confirmed that the inferred emission surfaces and velocity profiles presented in Section 3 do not change significantly. Likewise, we obtain consistent results with JvM-uncorrected cubes as we show in Appendix B.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top panels: peak brightness temperature maps of ¹²CO (left), ¹³CO (middle), and C¹⁸O J = 2–1 (right) lines. The ¹²CO emission from the western side of the disk experiences foreground cloud contamination as previously reported in independent data sets (Öberg et al. 2011; Huang et al. 2016; Guzmán et al. 2018; Teague et al. 2018b; Law et al. 2021a). Middle panels: emission surface heights above the midplane, taking into account the lowered emission surface across the gap (Model 2; see Section 3.1). For C¹⁸O, our model is consistent with a flat surface located at the disk midplane. Dashed ellipses and lines show constant radius and azimuth in the disk frame, with intervals of 0.′′5 and 30°, respectively. Bottom panels: centroid velocity maps $v_0$ (see Section 3.2). Synthesized beams are shown in the lower left corner of each panel.</figcaption></figure>

larger outer bounds of 2.′′5 in ¹²CO, 2.′′0 in ¹³CO, and 1.′′6 in C¹⁸O.

To fit the emission surface, we use *disksurf*<fnref n="25" /> (Teague et al. 2021), which implements the method outlined in Pinte et al. (2018a), who used the asymmetry of the line emission above the disk midplane to infer an emission height. This method allows us to locate emission arising from specific locations in the disk. We then use that information to construct the 3D structure of the emission layer. Following Law et al. (2021b), we use *disksurf*'s *get emission surface* function to extract the deprojected radius *r*, emission height *z*,

surface brightness $I_v$, and channel velocity *v* for each pixel associated with the emission. We do not exclude channels that suffer from foreground contamination, as including the contaminated channels is shown to have no significant effects on the retrieved surface (Teague et al. 2018b). For the initial geometric properties used to fit the surface, we assume the disk-center offsets $x_0$ and $y_0$ to be zero and adopt PA = 85.°8, inclination $i = 35°$, and stellar mass $M_* = 1.2 M_\odot$ from Öberg et al. (2021). We refit these parameters later on and confirm that the values we initially adopted describe the data well. For the individual pixels inferred from this procedure, we apply two constraints before we fit the emission surface. First, for all three isotopologues, we implement a minimum *z* value equal to minus half of the beam semimajor axis. This choice follows the

<footnote n="25">²⁵ https://disksurf.readthedocs.io/en/latest/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Colored circles represent individual pixels inferred by *disksurf*. Gray filled circles present radially binned pixels, binned by ~0.45 of the beam. Note that the radially binned emission surfaces are shown for visualization purposes only, and the emission surface is fit with individual pixels. We explore three surface models for ¹²CO and ¹³CO (see Section 3.1): Model 1 (smooth surface; solid curves) assumes a tapered power law described by Equation (1), Model 2 (Gaussian gap; dashed curves) adds a Gaussian gap to Model 1, and Model 3 (midplane gap; dotted curves) adopts the same gap center and width as in Model 2, but the gap extends to the midplane (i.e., $A_{gap} = 1$ in Equation (2)). For C¹⁸O, we assume a flat surface. The semimajor axis of the synthesized beam is shown in the upper right corner of each panel. Points above one beam size from the Model 1 curve have been removed for fitting the gap parameters.</figcaption></figure>

<table><caption>**Table 1**
Emission Surface Parameters Derived for ¹²CO and ¹³CO $J = 2-1$ Lines</caption><thead><tr><th></th><th>z<sub>0</sub><br>(au)</th><th>ψ</th><th>r<sub>taper</sub><br>(au)</th><th>q<sub>taper</sub></th><th>A<sub>gap</sub></th><th>r<sub>gap</sub><br>(au)</th><th>σ<sub>gap</sub><br>(au)</th></tr></thead><tbody><tr><td>¹²CO</td><td>22.99<sup>+1.21</sup><sub>-1.21</sub></td><td>1.47<sup>+0.13</sup><sub>-0.11</sub></td><td>217.8<sup>+7.26</sup><sub>-9.68</sub></td><td>3.69<sup>+0.53</sup><sub>-0.51</sub></td><td>0.6</td><td>216.6</td><td>24.2</td></tr><tr><td>¹³CO</td><td>9.68<sup>+2.42</sup><sub>-2.42</sub></td><td>4.53<sup>+0.34</sup><sub>-0.58</sub></td><td>152.5<sup>+9.68</sup><sub>-12.1</sub></td><td>4.10<sup>+0.64</sup><sub>-0.92</sub></td><td>0.6</td><td>216.6</td><td>24.2</td></tr></tbody></table>

Note. The errors represent statistical uncertainties and do not account for systematic ones.

methods from Law et al. (2021b), where large negative $z/r$ values were removed, but some negative values were allowed to remain to avoid positively biasing the resulting surface. Additionally, for ¹³CO and C¹⁸O, we remove the individual pixels that are above the ¹²CO emission surface because ¹³CO and C¹⁸O must be optically thinner than ¹²CO. Figure 2 shows the individual pixels after data cleaning. We then use the Markov Chain Monte Carlo (MCMC) method from *disksurf*, which wraps *emcee* (Foreman-Mackey et al. 2013), adopting 128 walkers, 500 burn-in steps, and 1000 steps to obtain $z_0$, $\psi$, $r_t$, and $q_t$. We confirmed the convergence of the MCMC fit by checking the posterior distribution. Throughout the paper, the emission surface obtained by this process is referred to as Model 1. Table 1 presents the fitted parameters.

Although Equation (1) describes the overall emission surface well, it cannot describe fine features, such as annular gaps. In particular, the gap at 1"/7 within which a candidate CPD-hosting planet is found (Bae et al. 2022b) cannot be described by Equation (1). To infer more accurate velocity structures within/around the gap, we add a Gaussian gap to the emission surface obtained in Model 1, adopting the following functional form:

$$z(r) = z_0 \left(\frac{r}{1''}\right)^\psi \exp\left(-\left[\frac{r}{r_t}\right]^{q_t}\right) \times \left(1 - A_{\text{gap}} \exp\left[-\frac{(r - r_{\text{gap}})^2}{2\sigma_{\text{gap}}^2}\right]\right). \quad (2)$$

Here $A_{gap}$, $r_{gap}$, and $\sigma_{gap}$ describe the depth of the gap, radial location of the center of the gap, and radial width of the gap, respectively. After obtaining the tapered power-law parameters using the aforementioned methods, we fit for only the gap parameters using *scipyeronimize*, *curve_fit*. For this process, we remove individual pixels above one beam from the Model 1 emission surface for a better convergence of the fit. The pixels removed through this procedure are less than 10% of all the pixels. We note that removing these individual pixels at high altitude estimates a deeper gap than would otherwise be found if these pixels were included. However, as we show below, the inferred velocity profiles are insensitive to the depth of the gap. As for Model 1, we sample the posterior distributions using an MCMC approach, adopting 128 walkers, 500 burn-in steps, and 1000 steps.

From now on, we refer to this surface with a Gaussian gap as Model 2, and this model is the main model we will use for our analysis. We do not fit the gap in ¹³CO separately because the ¹³CO emission is weak beyond ≈1"/5 and does not probe the full extent of the gap. Instead, we adopt the best-fit gap parameters from the ¹²CO data. As we found that the C¹⁸O emission surface is consistent with a flat surface at the disk midplane, we do not introduce a gap in the C¹⁸O surface (see Law et al. 2022b, for flat C¹⁸O emission surfaces in other disks). As such, throughout this paper we adopt a single model with zero emission height for C¹⁸O. Figure 2 shows emission surfaces from all the models. Table 1 presents the best-fit gap parameters.
<caption>**Table 2**</caption>

Geometric Properties Derived for the ¹²CO, ¹³CO, and C¹⁸O J = 2–1 Emission Assuming a Tapered Power Law with a Gaussian Gap (Model 2)

<table><thead><tr><th></th><th>x<sub>0</sub><br>(au)</th><th>y<sub>0</sub><br>(au)</th><th>PA<br>(deg)</th><th>M<sub>*</sub><br>(M<sub>☉</sub>)</th><th>v<sub>LSR</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td><sup>12</sup>CO</td><td>-3.87<sup>+0.08</sup><sub>-0.08</sub></td><td>1.67<sup>+0.12</sup><sub>-0.12</sub></td><td>84.88<sup>+0.07</sup><sub>+0.06</sub></td><td>1.24<sup>+0.003</sup><sub>-0.003</sub></td><td>4.64<sup>+0.001</sup><sub>-0.001</sub></td></tr><tr><td><sup>13</sup>CO</td><td>-1.67<sup>+0.36</sup><sub>-0.36</sub></td><td>0.73<sup>+0.29</sup><sub>-0.28</sub></td><td>86.01<sup>+0.21</sup><sub>+0.21</sub></td><td>1.25<sup>+0.007</sup><sub>-0.008</sub></td><td>4.65<sup>+0.003</sup><sub>-0.004</sub></td></tr><tr><td>C<sup>18</sup>O</td><td>-1.57<sup>+0.76</sup><sub>-0.72</sub></td><td>-0.57<sup>+0.69</sup><sub>-0.76</sub></td><td>86.15<sup>+0.50</sup><sub>-0.51</sub></td><td>1.26<sup>+0.01</sup><sub>-0.01</sub></td><td>4.66<sup>+0.008</sup><sub>-0.009</sub></td></tr></tbody></table>

**Note.** The errors represent statistical uncertainties and do not account for systematic ones.

Finally, we allow the Gaussian gap to reach the midplane by setting $A_{gap} = 1$, which we denote as Model 3. The purpose of having this hypothetical model is to allow the emission surface to reach the disk midplane and examine the effect of the gap depth in the derived velocity profile.

Once the emission surfaces are fitted, we take the best-fit values to infer the geometric properties of the disk using eddy<fnref n="26" /> (Teague 2019). We fit the disk-center offset $x_0$ and $y_0$, disk PA, dynamical stellar mass $M_*$, and LSR velocity of the target $v_{LSR}$, while the disk inclination is fixed to 35°, a value constrained by high-resolution continuum data (Huang et al. 2018). We use an MCMC method with the same setup mentioned previously. The geometric properties obtained using the Model 2 emission surface are listed in Table 2, while those derived using Model 1 and Model 3 are listed in Tables 3 and 4 in Appendix A. The geometric properties obtained via this method are broadly consistent with the dust-based values obtained in Huang et al. (2018), who find a PA of 85.76 ± 0.16. Our PA values for ¹³CO (86.01<sup>+0.21</sup><sub>-0.21</sub>) and C¹⁸O (86.148<sup>+0.50</sup><sub>-0.51</sub>) are closer to the value obtained via continuum fitting by Huang et al. (2018) likely because they trace closer to the midplane. These geometric properties are also consistent with those from Öberg et al. (2021).

## 3.2. Velocity Profiles

To infer the velocity profiles, we first make maps of the line centers $v_0$, using the quadratic method from bettermoments (Teague & Foreman-Mackey 2018b). The resulting $v_0$ maps are shown in the bottom panels of Figure 1. Then, with the derived emission surface and disk geometric properties, we decompose $v_0$ into the radial, rotational, and vertical velocities, following Teague et al. (2018a, 2018b, 2019a).

This is done by breaking apart the following equation:

$$v_0 = v_{\phi} \cos(\phi) \sin(|i|) + v_r \sin(\phi) \sin(i) - v_z \cos(i) + v_{LSR}, \quad (3)$$

assuming that $v_\phi$ and $v_r$ are azimuthally symmetric, where $v_\phi$ is the rotational velocity, $v_r$ is the radial velocity, $v_z$ is the vertical velocity, $i$ is the inclination of the disk,<fnref n="27" /> and $\phi$ is the azimuthal angle in the frame of reference of the disk. In practice, we use the getVelocity_profile module from eddy (Teague 2019) with 20 iterations to obtain the stacked spectra, each of which uses a random sample of independent pixels; a

<caption>**Table 3**</caption>

Geometric Properties Calculated Assuming a Tapered Power Law for the Emission Surfaces (Model 1)

<table><caption>**Table 4**</caption><thead><tr><th></th><th>x<sub>0</sub><br>(au)</th><th>y<sub>0</sub><br>(au)</th><th>PA<br>(deg)</th><th>M<sub>*</sub><br>(M<sub>☉</sub>)</th><th>v<sub>LSR</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td><sup>12</sup>CO</td><td>-3.86<sup>+0.08</sup><sub>-0.08</sub></td><td>1.66<sup>+0.12</sup><sub>-0.12</sub></td><td>84.88<sup>+0.07</sup><sub>+0.07</sub></td><td>1.24<sup>+0.003</sup><sub>-0.003</sub></td><td>4.64<sup>+0.001</sup><sub>-0.001</sub></td></tr><tr><td><sup>13</sup>CO</td><td>-3.82<sup>+0.35</sup><sub>-0.36</sub></td><td>0.71<sup>+0.29</sup><sub>-0.29</sub></td><td>86.03<sup>+0.21</sup><sub>+0.21</sub></td><td>1.25<sup>+0.008</sup><sub>-0.007</sub></td><td>4.65<sup>+0.004</sup><sub>-0.004</sub></td></tr></tbody></table>

Geometric Properties Calculated Assuming a Tapered Power Law with a Gaussian Gap That Reaches Down to the Midplane (Model 3)

<table><thead><tr><th></th><th>x<sub>0</sub><br>(au)</th><th>y<sub>0</sub><br>(au)</th><th>PA<br>(deg)</th><th>M<sub>*</sub><br>(M<sub>☉</sub>)</th><th>v<sub>LSR</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td><sup>12</sup>CO</td><td>-3.87<sup>+0.08</sup><sub>-0.08</sub></td><td>1.67<sup>+0.12</sup><sub>-0.12</sub></td><td>84.88<sup>+0.07</sup><sub>-0.06</sub></td><td>1.24<sup>+0.003</sup><sub>-0.003</sub></td><td>4.64<sup>+0.001</sup><sub>-0.002</sub></td></tr><tr><td><sup>13</sup>CO</td><td>-1.67<sup>+0.36</sup><sub>-0.36</sub></td><td>0.73<sup>+0.28</sup><sub>+0.28</sub></td><td>86.01<sup>+0.21</sup><sub>+0.21</sub></td><td>1.25<sup>+0.007</sup><sub>-0.008</sub></td><td>4.65<sup>+0.004</sup><sub>-0.004</sub></td></tr></tbody></table>

weighted average is then taken over these 20 samples to calculate $v_r$ and $v_\phi$. We choose this number of iterations based on Yu et al. (2021), who found that the gradient of the average standard deviation of the results flattens after about 20 iterations.

As in Teague et al. (2018b), we model the stacked spectrum with a Gaussian process, which allows for a more flexible and robust model (Foreman-Mackey et al. 2017). As can be seen in Equation (3), the vertical velocity has no dependence on $\phi$ and is thus not directly calculated by shifting and stacking spectra. Instead, to calculate $v_z$, we exploit Equation (3) and subtract projected radial and rotational velocities, along with $v_{LSR}$, from the $v_0$ map, following Yu et al. (2021).

The resulting velocity profiles for ¹²CO, ¹³CO, and C¹⁸O are shown in Figure 3. Looking at the rotational velocity first, we find evidence of super- and sub-Keplerian rotation in ¹²CO on the order of ±5% of the background Keplerian rotation. The sub-Keplerian rotation is most significant at ~150–240 au and has a double-peaked profile. At ~80–130 au and ≳230 au, the disk rotation is super-Keplerian, up to about 2% of the background Keplerian rotation. Overall, the ¹²CO rotational velocity profile is consistent with what was previously inferred by Teague et al. (2018b). The ¹³CO emission shows a rotational velocity profile that is broadly consistent with ¹²CO: the disk at ~90–170 au has super-Keplerian motion. Additionally, we find a rapid transition to sub-Keplerian rotation beyond 190 au. We conjecture that this is due to lower signal-to-noise ratio (S/N). We find a similar rapid transition to sub-Keplerian rotation in C¹⁸O beyond 170 au, likely due to low S/N.

Next, the radial velocity profile in ¹²CO shows a change in sign, from about −100 to 50 m s⁻¹, around 200 au. There is not a similar trend in ¹³CO, and the magnitude of the radial velocity is much smaller than that of ¹²CO, within ±20 m s⁻¹. The radial velocity of C¹⁸O is consistent with zero within uncertainties.

Lastly, we find a large upward vertical velocity flow in ¹²CO. This upward vertical motion is persistent from 140 to 220 au and has a maximum velocity of about 150 m s⁻¹ at a radius of ≈177 au, which corresponds to about 6% of the local Keplerian speed or 50% of the local sound speed adopting the 2D r-z gas temperature distribution inferred by Law et al. (2021b). The vertical velocity in ¹³CO emission also shows evidence of large coherent upward motions from 160 to 220 au, with a maximum

<footnote n="26">²⁶ https://eddy.readthedocs.io/en/latest/</footnote>

<footnote n="27">²⁷ In Equation (3), positive $i$ represents a disk that is rotating in a counterclockwise direction, while negative $i$ describes a clockwise rotation (Pinte et al. 2022). Because the AS 209 disk rotates clockwise, we adopt $i = -35^\circ$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Velocity profiles for ¹²CO (left column, blue), ¹³CO (middle column, red), and C¹⁸O (right column, purple). For ¹²CO and ¹³CO, we adopt the emission surfaces described by a tapered power law with a Gaussian gap (i.e., Model 2; Equation (2)), while a flat emission surface at the midplane is assumed for C¹⁸O. For each isotopologue, the first rows show the rotational velocity $v_\phi$; the second row shows the perturbations in $v_\phi$ from Keplerian rotation $v_{Kep}$, $\delta v_\phi / v_{Kep} \equiv (v_\phi - v_{Kep}) / v_{Kep}$; the third rows show the radial velocity $v_r$; and the fourth row shows the vertical velocity, $v_z$. The vertical dark-gray dotted line shows the radial location of the CPD. The shaded regions represent the uncertainties on the velocities; these are statistical uncertainties and do not account for systematic ones.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Gas flows in the r-z plane for ¹²CO (blue), ¹³CO (red), and C¹⁸O (purple). Note the strong upward gas flows arising from the gap. The emission surfaces are shown with thick curves. The arrows are scaled to the local sound speed $c_s$, where we use the 2D r-z gas temperature distribution inferred by Law et al. (2021b). The size of the semimajor axis of the synthesized beam (¹²CO) is shown in the upper left corner. An arrow showing 10% of the sound speed is also presented in the upper left corner. The vertical dark-gray dotted line shows the location of the CPD.</figcaption></figure>

velocity of 85 m s⁻¹ at 193 au. In C¹⁸O, we see a much smaller upward motion, reaching a maximum of about 10 m s⁻¹, but note that C¹⁸O emission is weak and does not probe the radial regions where strong upward motions are seen in ¹²CO or ¹³CO. These velocity profiles are broadly consistent with what are found by Izquierdo et al. 2023, where the authors carried out an independent kinematic analysis on the same data obtained by the MAPS program. In Section 4.1 we discuss the potential origin of these coherent, large-scale upward flows.

We examined how (in)sensitive the inferred velocity profiles are to the assumed emission surface by repeating the analysis and deriving velocity profiles adopting Model 1 (tapered power-law emission surface without a gap) and Model 3

(tapered power-law emission surface with a Gaussian gap that reaches the midplane). As we show in Figure 7 in Appendix A, varying the emission surfaces does not have a significant impact on the velocity profiles. For the rest of the paper, we thus opt to use Model 2 for our discussion. To help visualize the inferred velocity flows along with the emission surfaces, in Figure 4 we depict the gas flows in the r-z plane. As shown, it is apparent that the large upward motions in ¹²CO and ¹³CO coincide with the gap in the disk.

In addition to searching for kinematic structures in azimuthally averaged radial profiles of the velocity, we investigate structure within the deprojected residual velocity maps. To do so, we calculate a best-fit Keplerian model with
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Deprojected velocity residual maps for $^{12}$CO, $^{13}$CO, and C$^{18}$O, calculated by subtracting best-fit Keplerian models ($v_{\text{mod}}$) from the centroid velocity maps ($v_0$; see Figure 1, bottom panels). Note that in $^{12}$CO and $^{13}$CO the velocity structures within the gap at ≃1′′5–2′′ are present around the entire azimuth. The dashed circles show lines of constant radius at 1′′3 and 1′′7, highlighting the velocity structure associated with upward vertical motions arising from the gap. Synthesized beams are shown in the lower left corner of each panel.</figcaption></figure>

eddy, adopting the emission surfaces and the disk geometric properties constrained as in Section 3.1. This produces a model map, $v_{\text{mod}}$, which we subtract from the line centroid map, $v_0$ (shown in Figure 1). Figure 5 shows the resulting residual maps in the deprojected disk plane for all three CO isotopologues. As seen in $^{12}$CO and $^{13}$CO, the velocity structure in the residual maps is mostly azimuthally symmetric, indicating that the velocity perturbation contributions are largely from the vertical component (Teague et al. 2019b). We find no clear asymmetric features associated with the planet candidate proposed by Bae et al. (2022b). However, as we will discuss in Section 4.1, the kinematics of the gap is likely dominated by disk winds, not the planet, and we emphasize that lack of asymmetric features in the residual velocity maps does not dispute the presence of a planet.

# 4. Discussion

## 4.1. Origin of the Vertical Flows

The most prominent kinematic structure found in the AS 209 disk is the upward flow in $^{12}$CO and $^{13}$CO, arising from the gap at 1′′2–1′′8 (145–218 au). In this section, we explore several possibilities to explain the upward flows.

### 4.1.1. Giant Planet

Giant planets are expected to perturb the velocity structure of the disk. When a giant planet opens a gap, steep density gradients develop on the sides of the gap, driving sub-/super-Keplerian rotation (Kanagawa et al. 2015), as well as a downward flow back into the midplane (Kley et al. 2001; Gressel et al. 2013; Morbidelli et al. 2014; Szulágyi et al. 2014; Fung & Chiang 2016). However, this is not the trend we see in the AS 209 disk. As shown in Figure 4, the $^{12}$CO and $^{13}$CO vertical velocity patterns reveal upward motion at the center of the gap that tapers off at the gap edges—a meridional *fountain*. This upward motion is seen across the gap over a broad range of azimuth (only a small section of azimuth in the 4 o’clock direction does not exhibit the upward flows within the gap; see Figure 5), so it is unlikely that the observed upward flows are associated with a jet or outflow arising locally from the embedded planet.

Alternatively, one might ask whether we are seeing downward flows toward the midplane from the back side of the disk. This may be possible when the front side of the disk is sufficiently optically thin; however, $^{12}$CO remains optically thick within the gap, supported by the fact that the $^{12}$CO brightness temperature within the gap is >20 K (Law et al. 2021b) and that the CPD is visible only in $^{13}$CO and not in $^{12}$CO (Bae et al. 2022b). Overall, the upward flows seen in AS 209 are not straightforward to reconcile with the presence of a giant planet alone.

### 4.1.2. Ambipolar-diffusion-driven Winds

To explain both the presence of the CPD-hosting planet previously reported in Bae et al. (2022b) and the azimuthally symmetric upward gas flows found in this paper, we propose a scenario where the low density within the planet-carved gap triggers magnetically driven winds via ambipolar diffusion. Ambipolar diffusion is the dominant nonideal magnetohydrodynamic (MHD) effect when the underlying gas density and ionization levels are low (Wardle 2007). When ambipolar diffusion dominates the gas dynamics, ions that are coupled to the magnetic fields can drag neutral molecules/atoms, driving winds (Bai & Stone 2013; Gressel et al. 2015; Béthune et al. 2017; Suriano et al. 2018; Hu et al. 2022; see also the review by Lesur et al. 2022).

To examine this possibility more quantitatively, we estimate the ratio of the ion-neutral drift time to the dynamical time by calculating the ambipolar Elsässer number Am, given by

$$ Am \equiv \frac{v_A^2}{\eta_A \Omega_K}, \qquad (4) $$

where $v_A$ is the Alfvén velocity, $\eta_A$ is the ambipolar diffusivity, and $\Omega_K$ is the Keplerian frequency (Bai 2011). Using $v_A = B/\sqrt{4\pi\rho}$ and $\eta_A = B^2/(4\pi\gamma\rho\rho_i)$, where B is the magnetic field strength, $\rho$ is the density of the neutral gas, and $\rho_i$ is the density of the ionized gas, Equation (4) turns into

$$ Am = \frac{\gamma \rho_i}{\Omega_K}. \qquad (5) $$
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Ambipolar Elsässer number, Am, calculated for the AS 209 disk (green) and the HD 163296 disk (pink). The triangles show upper limits on Am when HCO⁺ *J* = 1–0 is not detected. The shaded regions represent the uncertainty on the measurements, which are derived from the uncertainty in the *N*(HCO⁺) measurement. The green dashed lines and arrows show the radial region where winds are found in AS 209, coincident with the annular gap in ¹²CO and the CPD. The gray dashed line shows the radial location of the CPD in AS 209 (Bae et al. 2022b). The pink dashed line and arrow show the radial region where winds are found in HD 163296 (Teague et al. 2019a).</figcaption></figure>

Here $\gamma \equiv \langle \sigma v \rangle_i / (m_n + m_i)$, where $\langle \sigma v \rangle_i$ is the momentum transfer rate coefficient for an ion-neutral collision, given by

$$ \langle \sigma v \rangle_i = 2.0 \times 10^{-9} \left( \frac{m_H}{\mu} \right)^{1/2} \text{cm}^3 \text{s}^{-1}, \quad (6) $$

where $m_n$ and $m_i$ are the mass of neutrals and ions, respectively, and $\mu \equiv m_n m_i / (m_n + m_i)$ is the reduced mass (Draine 2011).

Thermochemical models of protoplanetary disks suggest that HCO⁺ is the most abundant molecular ion in the warm molecular layer where CO gas is abundant (e.g., Aikawa et al. 2015). In fact, in all five disks observed in MAPS, the HCO⁺ column density is greater than that of N₂H⁺ and N₂D⁺ (two other ions that are believed to be abundant in protoplanetary disks) by more than an order of magnitude (Aikawa et al. 2021). Assuming that HCO⁺ and H₂ are the dominant ions and neutrals in protoplanetary disks, we obtain $\gamma = 2.82 \times 10^{13} \text{cm}^3 \text{g}^{-1} \text{s}^{-1}$. To compute the ion density $\rho_i$, we use the column density of HCO⁺ observationally constrained by Aikawa et al. (2021)<fnref n="28" /> and divide it by the gas pressure scale height, $\rho_i = m_i n_i = m_{HCO^+} \times N(\text{HCO}^+)/H(r)$, motivated by thermochemical models where HCO⁺ forms a layer having an approximately constant volume density (Aikawa et al. 2021). We calculate the scale height using a power law with a flaring index determined by Zhang et al. (2021) for the AS 209 disk. Figure 6 shows the derived ambipolar Elsässer number as a function of disk radius. Within the inner ~100 au Am is about 10, but it drops to ~0.1 at 200 au owing to the low HCO⁺ density. Recent nonideal MHD simulations have shown that when Am drops to ~1 ambipolar diffusion starts to quench the MRI (Bai & Stone 2011) and launches winds (Bai & Stone 2013; Gressel et al. 2015; Suriano et al. 2018). The inferred ambipolar Elsässer number of ~0.1 at the radial region of emerging vertical flows is thus broadly consistent with these

numerical simulations. In the case in which the outer disk is transparent to UV radiation and C⁺ dominates the ion density instead of HCO⁺, our Am estimates would provide a lower limit.

In the HD 163296 disk, Teague et al. (2019a, 2022) found upward<fnref n="29" /> meridional flows, most prominently at 240 au, the radial location of a kinematically inferred planet (Pinte et al. 2018b; Teague et al. 2018a), and radially outward disk winds beyond ~380 au. Are these findings consistent with the picture we propose for AS 209? In order to answer this question, we compute the ambipolar Elsässer number in the HD 163296 disk using the same methods we applied to the AS 209 disk, using the HCO⁺ column density from Aikawa et al. (2021) and the scale height from Zhang et al. (2021). The resulting radial profile of the ambipolar Elsässer number is shown in Figure 6. As in AS 209, the ambipolar Elsässer number is ~10 in the inner ~100 au of the disk. At 240 au in HD 163296, the ambipolar Elsässer number is ~0.5, and beyond ~380 au the ambipolar Elsässer number drops to ~0.1. Within the two disks, we find an overall trend that winds appear in the radial regions having a low ambipolar Elsässer number of ~0.1–0.5, suggesting that winds driven by ambipolar diffusion may be common in low-density regions of protoplanetary disks.

Besides a low density, ambipolar diffusion is more efficient in the presence of strong magnetic fields. We estimate the required magnetic field strength in the AS 209 disk by using the magnetic diffusion numbers from Wardle (2007), who defines the ambipolar regime as being dominant when $1 \ll \beta_i \ll \beta_e$, where the magnetic diffusion numbers, $\beta_i$ and $\beta_e$, are given by

$$ \beta_i \approx 4.6 \times 10^{-3} \left( \frac{B}{1 \text{G}} \right) \left( \frac{n_H}{10^{15} \text{cm}^{-3}} \right)^{-1} \quad (7) $$

and

$$ \beta_e = 3.5 \left( \frac{B}{1 \text{G}} \right) \left( \frac{n_H}{10^{15} \text{cm}^{-3}} \right)^{-1} \left( \frac{T}{100 \text{K}} \right)^{-1/2} \quad (8) $$

In the above equations, *B* is the magnetic field in gauss and $n_H$ is the number density of hydrogen nuclei in units of cm⁻³. A lower limit on the magnetic field strength that satisfies $\beta_i \gg 1$ is obtained using the hydrogen column density. Note that the second part of the inequality, $\beta_i \ll \beta_e$, is always satisfied with gas temperatures of tens to hundreds of kelvin, as can be seen from Equations (7) and (8). We compute the hydrogen column density $N(H_2)$ using the gas surface density derived in Zhang et al. (2021; see their Section 4.1 and Figure 16). At 200 au, $N(H_2) \approx 1.8 \times 10^{21} \text{cm}^{-2}$. Adopting the scale height at 200 au of $H = 14.3$ au from Zhang et al. (2021), the number density of hydrogen nuclei at the midplane is $n_H = (1/2) \times N(H_2)/\sqrt{2\pi H} \approx 1.7 \times 10^6 \text{cm}^{-3}$. Inserting this hydrogen nuclei number density into Equation (7), we find that a weak magnetic field strength of $B \gg 0.36 \mu\text{G}$ is sufficient for ambipolar diffusion to dominate at 200 au. Note also that the weak required magnetic field strength is consistent with nondetection of magnetic fields in the AS 209 disk (3σ upper limits of a few mG) via observations of Zeeman splitting of the CN *N* = 21 line (Harrison et al. 2021).

<footnote n="28">²⁸ Available to download at the MAPS webpage: https://alma-maps.info/.</footnote>

<footnote n="29">²⁹ The sign of the vertical velocity extracted in previous papers (Teague et al. 2019a, 2022) was incorrect and needs to be flipped.</footnote>
### 4.1.3. Vertical Shear Instability

Vertical shear in the rotational velocity of the disk gas can lead to an instability that can produce vertical flows when saturated (Nelson et al. 2013). Barraza-Alfaro et al. (2021) showed that vertical flows driven by the vertical shear instability (VSI) can manifest as nearly concentric rings of upward and downward flows in the Keplerian-subtracted centroid velocity maps of molecular line emission. However, we conjecture that the VSI is less likely to be the origin of the vertical flows seen in the AS 209 disk because the radial extent of the vertical flow in the AS 209 disk is much larger than what is typically seen in numerical simulations of the VSI. The radial width of the VSI-induced vertical flows in numerical simulations is about a gas scale height (Nelson et al. 2013; Barraza-Alfaro et al. 2021). On the other hand, the upward flow in the AS 209 disk spans about $0.5^{\prime\prime} \approx 60.5$ au, which corresponds to about 4 scale heights at $1.7^{\prime\prime} \approx 206$ au adopting the midplane scale height of $0.12^{\prime\prime} \approx 14.7$ au from Zhang et al. (2021).

In summary, we conclude that ambipolar-diffusion-driven winds from a planet-carved gap are the most viable origin for the observed vertical flows in the AS 209 disk.

### 4.2. Can Winds Stop the Growth of Giant Planets?

In the traditional meridional circulation picture without winds, the rate at which a planet would grow depends on the rate of the circumstellar disk gas flowing into the planet-carved gap (e.g., Morbidelli et al. 2014). In this picture, planets can continuously grow in mass until the circumstellar disk loses most of its mass. Indeed, hydrodynamic simulations showed that the mass-doubling time for a Jovian-mass planet is of order of 100–1000 orbital times, which can be much shorter than the lifetime of protoplanetary disks depending on the radial location of the planet (e.g., Kley 1999; Lubow et al. 1999).

In our modified picture considering ambipolar-diffusion-driven winds, the mass outflow rate via winds can exceed the mass inflow rate toward the gap, in which case the growth of the embedded planet can be limited or even ceased. In the AS 209 disk, we can estimate the mass-loss rate from the annular gap using the following equation:

$$ \dot{M}_{\text{wind}} = \int_{r_{\text{in}}}^{r_{\text{out}}} 2\pi r \rho v_{\text{wind}} dr, \quad (9) $$

where $r_{\text{in}} = 1.3^{\prime\prime} (\sim 160$ au) and $r_{\text{out}} = 2^{\prime\prime} (\sim 240$ au) are the inner and outer boundaries of the wind-launching region, respectively, and $\rho$ and $v_{\text{wind}}$ are the gas density and speed of the wind, respectively. For simplicity, we opt to use the midplane density, $\rho = \rho_{\text{mid}}$, and the vertical velocity of $^{12}$CO, $v_{\text{wind}} = v_z(^{12}\text{CO})$. To calculate $\rho_{\text{mid}}$, we use the gas surface density derived by Zhang et al. (2021) (see Section 4.1.2) divided by $\sqrt{2\pi H(r)}$: $\rho_{\text{mid}}(r) = \Sigma(r)/\sqrt{2\pi H(r)}$. With this, we estimate the mass-loss rate via winds to be $4.4 \times 10^{-8} M_{\odot} \text{yr}^{-1}$. In reality, the gas density of the wind can be smaller than the midplane density. In AS 209, the $^{12}$CO emission surface lies within $\approx 2$ scale heights from the midplane, so, assuming vertical hydrostatic equilibrium, the mass-loss rate can be reduced by a factor of $e^{-2} \approx 0.14$. Taking this into account, the mass-loss rate via winds is $6.2 \times 10^{-9} M_{\odot} \text{yr}^{-1}$. Calculating the total mass within the gap using the surface density from Zhang et al. (2021), we find that the gap would be depleted in a

minimum of $1.4 \times 10^4$ yr if the mass-loss rate is maintained and there is no gas radially advected into the gap.

Next, we estimate the mass inflow rate to the gap assuming an absence of winds using

$$ \dot{M}_{\text{in}} = 2\pi r v_r \Sigma_{\text{in}}, \quad (10) $$

where $v_r$ is the radial velocity of the inflowing gas and $\Sigma_{\text{in}}$ is the surface density of the gas that falls into the gap. For a steady-state viscous disk, the radial velocity $v_r$ can be described by $v_r = \alpha H \Omega_{\text{K}}$, where $\alpha$ is the coefficient characterizing the efficiency of the accretion (regardless of the origin), $H$ is the scale height, and $\Omega_{\text{K}}$ is the Keplerian orbital frequency. Adopting $H = 14.3$ au and $\Sigma_{\text{in}} = \Sigma_{\text{gas}} \approx 0.006 \text{ g cm}^{-2}$ at 200 au (Figure 16 of Zhang et al. 2021), we estimate a mass inflow rate of $3.0 \times 10^{-11} \times (\alpha/10^{-3}) M_{\odot} \text{yr}^{-1}$, which is smaller than $\dot{M}_{\text{wind}}$ unless $\alpha \gtrsim 0.2$. This means that the strong upward flows seen in AS 209 can lead to mass loss from the gap, possibly halting the growth of the embedded planet and depleting the gas inside the gap.

Up to this point, our discussion has been focused on AS 209, but we can use the theory discussed to make a general scaling relation for wind-regulated terminal mass of giant planets. The depth of a gap opened by a planet can be described by

$$ \frac{\Sigma_{\text{gap}}}{\Sigma} = \frac{1}{1 + 0.04K}, \quad (11) $$

where $\Sigma_{\text{gap}}$ is the surface density at the center of the gap, $\Sigma$ is the unperturbed surface density, and $K \equiv (M_p/M_*)^2(H/R)^{-5}\alpha^{-1}$ (Kanagawa et al. 2015). For $K \gg 1$, applicable for planets opening a deep gap, we can write Equation (11) as

$$ \frac{\Sigma_{\text{gap}}}{\Sigma} \approx \frac{1}{0.04K} = 0.25 \left( \frac{M_p/M_*}{10^{-3}} \right)^{-2} \left( \frac{H/R}{0.1} \right)^5 \left( \frac{\alpha}{10^{-3}} \right). \quad (12) $$

We can then relate the surface density at the gap center to the ambipolar Elsässer number as

$$ \text{Am} = \frac{\gamma \rho_i}{\Omega_{\text{K}}} = \frac{\gamma \chi_i \rho}{\Omega_{\text{K}}} = \frac{\gamma \chi_i \Sigma_{\text{gap}}}{\Omega_{\text{K}} \sqrt{2\pi H}}, \quad (13) $$

where $\chi_i \equiv \rho_i/\rho$ is the ionization fraction. Inserting $\Sigma_{\text{gap}}$ from Equation (12) into Equation (13) and reorganizing the equation, we obtain the terminal mass of a giant planet as follows:

$$ \begin{aligned} \frac{M_p}{M_*} = {}& 9.7 \times 10^{-4} \left(\frac{\text{Am}}{0.1}\right)^{-1/2} \left(\frac{\chi_i}{10^{-9}}\right)^{1/2} \left(\frac{H/R}{0.1}\right)^4 \\ & \times \left(\frac{\alpha}{10^{-3}}\right)^{1/2} \left(\frac{\Sigma}{1 \text{ g cm}^{-2}}\right)^{1/2} \left(\frac{R}{100 \text{ au}}\right)^{1/4} \\ & \times \left(\frac{M_*}{1 M_{\odot}}\right)^{-1/4} . \end{aligned} \quad (14) $$

Using the fiducial parameters used in Equation (14), the terminal mass of a giant planet around a solar-mass star is about a Jupiter mass at 100 au.

Despite the prevalence of substructures in protoplanetary disks, attempts to search for young, forming planets through direct imaging resulted in a low detection rate (see review by Benisty et al. 2022, and references therein). The properties of observed substructures suggest that the majority of the young planet population has (sub-)Jovian mass (Bae et al. 2018, 2022a; Zhang et al. 2018; Lodato et al. 2019). The wind-regulated terminal mass we estimated above coincides with the planet mass inferred from
substructure properties, potentially helping to explain the dearth of directly imaged super-Jovian-mass young planets. Future kinematic studies of a larger sample of protoplanetary disks will enable us to test whether wind-regulated growth of young planets is ubiquitous.

In the discussion above, we simplified the picture by assuming that there is no mass being fed to the CPD in the presence of large-scale winds. In hydrodynamic simulations, it is shown that the circumstellar disk gas can be supplied to the CPD through nonaxisymmetric flows (Lubow et al. 1999). Whether the same can happen in the presence of large-scale magnetically driven winds needs to be tested in the future, using nonideal magnetohydrodynamic simulations with an embedded planet. If mass can still be supplied to the CPD in the presence of large-scale magnetically driven winds, the wind-regulated terminal mass in Equation (14) would provide a lower limit to the final mass of the planet.

## 5. Summary

We have used $^{12}$CO, $^{13}$CO, and C$^{18}$O $J=2-1$ emission to carry out the detailed analysis of the kinematics within the AS 209 disk. We found significant perturbations in the rotational velocity in $^{12}$CO, up to 5% of the Keplerian rotation, which is consistent with previous findings by Teague et al. (2018b). In addition to the perturbations in the rotational velocity, we found a strong meridional fountain (coherent upward flows) in $^{12}$CO and $^{13}$CO at $\simeq 1''7$ (200 au). The upward flows are as fast as $150 \text{ m s}^{-1}$ in $^{12}$CO, corresponding to about 50% of the local sound speed or 6% of the local Keplerian speed. Interestingly, these upward flows are colocated with an annular gap within which a candidate CPD is recently reported (Bae et al. 2022b).

The observed upward flows are in the opposite direction to collapsing, downward flows within planet-carved gaps seen in hydrodynamic planet-disk interaction simulations and are difficult to explain with an embedded planet alone. Instead, we propose a scenario in which the low density within the planet-carved gap has triggered magnetically driven winds via ambipolar diffusion. To support this idea, we estimated the ambipolar Elsässer number using the $HCO^+$ column density. At the radial location of the upward flows, we found that the ambipolar Elsässer number is about 0.1, broadly consistent with the value at which ambipolar diffusion drives strong winds in numerical simulations. In this scenario, we hypothesize that magnetically driven winds from a planet-carved gap can limit/cease the growth of the planet embedded in the gap. This may be the explanation for the dearth of detections of gas giant planets in disks with observed dust substructure with ALMA. We also provided a scaling relationship that describes the wind-regulated terminal mass. Using parameters generally applicable to protoplanetary disks, we found that the wind-regulated terminal mass around a solar-mass star is about a Jupiter mass at 100 au, which can explain the dearth of directly imaged super-Jovian-mass young planets at large orbital distances.

These results show compelling kinematic evidence of disk winds arising from the gap opened by a forming planet. In the future, constraining the ion density beyond $HCO^+$ will help better constrain the environment under which ambipolar-diffusion-driven winds are launched. Observations constraining the morphology and strength of the magnetic fields in the AS 209 disk would help better understand the complex interplay between a forming planet and disk winds. Observations of species that can probe the warm outflowing gas from the low-density, higher regions, such as CI (Gressel et al. 2020; Alarcón et al. 2022), could help further characterize the nature

of the winds in the AS 209 disk. Kinematic studies for a larger sample of protoplanetary disks will help assess whether winds launched from planet-carved gaps are common or the AS 209 disk is a unique case. Additionally, nonideal magnetohydrodynamic planet-disk interaction simulations can prove (or dispute) that the activation of magnetically driven winds within planet-carved gaps can regulate the growth of embedded planets. Finally, numerical studies of orbital migration in a disk with active winds will allow us to infer whether the CPD-hosting planet in the AS 209 disk had formed at the current radial location or had formed at a different radial location but experienced inward/outward migration.

## Acknowledgments

We thank Xiao Hu, Geoffroy Lesur, and Gaylor Wafflard-Fernandez for helpful discussion. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01055.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. V.V.G. gratefully acknowledges support from FONDECYT Regular 1221352, ANID BASAL projects ACE210002 and FB210003, and ANID—Millennium Science Initiative Program—NCN19_171. Support for J.H. was provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51460.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. R.L.G. acknowledges support from the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP cofunded by CEA and CNES. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (PROTOPLANETS, grant agreement No. 101002188). C.W. acknowledges financial support from the University of Leeds, the Science and Technology Facilities Council, and UK Research and Innovation (grant Nos. ST/T000287/1 and MR/T040726/1). This work was supported by a grant from the Simons Foundation (686302, KIÖ) and an NSF AAG grant (No. 1907653, KIÖ). Support for F.L. was provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51512.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Y.A. acknowledges support by NAOJ ALMA Scientific Research Grant code 2019-13B.

Facility: ALMA.

Software: bettermoments (Teague & Foreman-Mackey 2018b), CASA (McMullin et al. 2007), disksurf (Teague et al. 2021), eddy (Teague 2019), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), SciPy (Virtanen et al. 2020).

## Appendix A

### Results with Additional Emission Surface Models

In Tables 3 and 4, we list $x_0$, $y_0$, PA, $M_*$, and $v_{LSR}$ fitted with Models 1 and 3, respectively. Figure 7 compares $^{12}$CO and $^{13}$CO velocity profiles for Models 1, 2, and 3. Note that the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison of velocity profiles found for ¹²CO (left) and ¹³CO (right) adopting different emission surface models. Orange: emission surfaces without a Gaussian dip (Model 1). Blue (¹²CO) and red (¹³CO): a tapered power law with a Gaussian gap (Model 2). Green: emission surface with Gaussian dip going to the midplane (Model 3). Note that the derived velocity profiles are insensitive to the emission surface models we adopt. The vertical dark-gray dotted line shows the radial location of the CPD (1.77 ± 0.2 au).</figcaption></figure>

derived velocity profiles are insensitive to the emission surface models we adopt.

## Appendix B

### Results with JvM-uncorrected Cubes

In this appendix we repeat the velocity analysis presented in Section 3 but with JvM-uncorrected data cubes. For

consistency, we use the emission surfaces and geometric properties derived from the JvM-corrected cubes. The resulting velocity profiles are shown in Figure 8. As shown in the figure, the inferred velocity profiles with the JvM-uncorrected data are broadly consistent with what we obtained with the JvM-corrected data. Most importantly, the upward motions at ~150–200 au are recovered.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Same as Figure 3, but with cubes that have not undergone JvM correction. Models 1, 2, and 3 represent the three emission surface models shown in Figure 2. The black curves are the velocity profiles for the JvM-corrected cubes for comparison. We do not present C¹⁸O because the S/N of the JvM-uncorrected data is not sufficient for this velocity analysis.</figcaption></figure>

## ORCID iDs

Maria Galloway-Sprietsma https://orcid.org/0000-0002-5503-5476

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

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

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

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

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

Felipe Alarcón https://orcid.org/0000-0002-2692-7862

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

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

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

L. Ilsedore Cleeves https://orcid.org/0000-0003-2076-8001

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

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

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

Yao Liu https://orcid.org/0000-0002-7616-666X

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

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

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

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

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

## References

Aikawa, Y., Cataldi, G., Yamato, Y., et al. 2021, ApJS, 257, 13

Aikawa, Y., Furuya, K., Nomura, H., & Qi, C. 2015, ApJ, 807, 120

Alarcón, F., Bergin, E. A., & Teague, R. 2022, ApJL, 941, L24

Andrews, S. M. 2020, ARA&A, 58, 483

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

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502

Bae, J., Isella, A., Zhu, Z., et al. 2022a, arXiv:2210.13314

Bae, J., Pinilla, P., & Birnstiel, T. 2018, ApJL, 864, L26

Bae, J., Teague, R., Andrews, S. M., et al. 2022b, ApJL, 934, L20

Bai, X.-N. 2011, ApJ, 739, 51

Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144

Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76

Barraza-Alfaro, M., Flock, M., Marino, S., & Pérez, S. 2021, A&A, 653, A113

Benisty, M., Dominik, C., Follette, K., et al. 2022, arXiv:2203.09991.

Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75

Cieza, L. A., González-Ruilova, C., Hales, A. S., et al. 2021, MNRAS, 501, 2934

Currie, T., Lawson, K., Schneider, G., et al. 2022, NatAs, 6, 751

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

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton, NJ: Princeton Univ. Press)

Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220

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

Fung, J., & Chiang, E. 2016, ApJ, 832, 105

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 650, C3

Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59

Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126

Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84

Guzmán, V. V., Huang, J., Andrews, S. M., et al. 2018, ApJL, 869, L48

Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, NatAs, 3, 749

Harrison, R. E., Looney, L. W., Stephens, I. W., et al. 2021, ApJ, 908, 141

Hu, X., Li, Z.-Y., Zhu, Z., & Yang, C.-C. 2022, MNRAS, 516, 2006

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

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

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

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

Izquierdo, A. F., Testi, L., Facchini, S., et al. 2023, arXiv:2304.03607

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

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

Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44

Kley, W. 1999, MNRAS, 303, 696

Kley, W., D'Angelo, G., & Henning, T. 2001, ApJ, 547, 457

Law, C. J., Crystian, S., Teague, R., et al. 2022a, ApJ, 932, 114

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

Law, C. J., Teague, R., Loomis, R. A., et al. 2021b, ApJS, 257, 4
Law, C. J., Teague, R., Öberg, K. I., et al. 2022b, ApJ, 948, 60

Lesur, G., Ercolano, B., Flock, M., et al. 2022, arXiv:2203.09821

Lodato, G., Dipierro, G., Ragusa, E., et al. 2019, MNRAS, 486, 453

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

Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001

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

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

Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610

Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Öberg, K. I., Qi, C., Fogel, J. K. J., et al. 2011, ApJ, 734, 98

Perez, 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., Teague, R., Flaherty, K., et al. 2022, arXiv:2203.09528

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

Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14

Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2018, MNRAS, 477, 1239

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

Teague, R. 2019, JOSS, 4, 1220

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

Teague, R., Bae, J., Andrews, S. M., et al. 2022, ApJ, 936, 163

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. 2018a, RNAAS, 2, 173

Teague, R., & Foreman-Mackey, D. 2018b, Bettermoments: A Robust Method To Measure Line Centroids, v1.0, Zenodo, doi:10.5281/zenodo.1419754

Teague, R., Law, C. J., Huang, J., & Meng, F. 2021, JOSS, 6, 3827

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

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

Wardle, M. 2007, Ap&SS, 311, 35

Wölfer, L., Facchini, S., van der Marel, N., et al. 2022, A&A, 670, A154

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

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

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