# Mapping the Vertical Gas Structure of the Planet-hosting PDS 70 Disk

Charles J. Law, Myriam Benisty, Stefano Facchini, Richard Teague, Jaehan

Bae, Andrea Isella, Inga Kamp, Karin I. Öberg, Bayron Portilla-Revelo, Luna

Rampinelli

► To cite this version:

Charles J. Law, Myriam Benisty, Stefano Facchini, Richard Teague, Jaehan Bae, et al. Mapping the Vertical Gas Structure of the Planet-hosting PDS 70 Disk. The Astrophysical Journal, 2024, 964, 10.3847/1538-4357/ad24d2. insu-04836885

HAL Id: insu-04836885

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

Submitted on 13 Dec 2024

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

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

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

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

# Mapping the Vertical Gas Structure of the Planet-hosting PDS 70 Disk

Charles J. Law<sup>1,2,10</sup>①, Myriam Benisty<sup>3,4</sup>①, Stefano Facchini<sup>5</sup>①, Richard Teague<sup>6</sup>①, Jaehan Bae<sup>7</sup>①, Andrea Isella<sup>8</sup>①, Inga Kamp<sup>9</sup>①, Karin I. Öberg<sup>2</sup>①, Bayron Portilla-Revelo<sup>9</sup>①, and Luna Rampinelli<sup>5</sup>

Received 2023 July 31; revised 2023 November 30; accepted 2024 January 3; published 2024 March 29

## Abstract

PDS 70 hosts two massive, still-accreting planets and the inclined orientation of its protoplanetary disk presents a unique opportunity to directly probe the vertical gas structure of a planet-hosting disk. Here, we use high-spatial-resolution ($\approx$0.1; 10 au) observations in a set of CO isotopologue lines and $HCO^+$ $J=4-3$ to map the full 2D ($r, z$) disk structure from the disk atmosphere, as traced by $^{12}$CO, to closer to the midplane, as probed by less abundant isotopologues and $HCO^+$. In the PDS 70 disk, $^{12}$CO traces a height of $z/r \approx 0.3$, $^{13}$CO is found at $z/r \approx 0.1$, and $^{C^{18}O}$originates at, or near, the midplane. The $HCO^+$ surface arises from $z/r \approx 0.2$ and is one of the few non-CO emission surfaces constrained with high-fidelity in disks to date. In the $^{12}$CO $J=3-2$ line, we resolve a vertical dip and steep rise in height at the cavity wall, making PDS 70 the first transition disk where this effect is directly seen in line-emitting heights. In the outer disk, the CO emission heights of PDS 70 appear typical for its stellar mass and disk size and are not substantially altered by the two inner embedded planets. By combining CO isotopologue and $HCO^+$ lines, we derive the 2D gas temperature structure and estimate a midplane CO snowline of $\approx 56-85$ au. This implies that both PDS 70b and 70c are located interior to the CO snowline and are likely accreting gas with a high C/O ratio of $\approx 1.0$, which provides context for future planetary atmospheric measurements from, e.g., JWST, and for properly modeling their formation histories.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Planet formation (1241); CO line emission (262); High angular resolution (2167)*

*Supporting material: data behind figures, figure set*

## 1. Introduction

Planets assemble and acquire their compositions from gas and dust in their natal protoplanetary disks, while the planet formation process is also expected to simultaneously alter the disk physical and chemical structure (e.g., Cleeves et al. 2015; Facchini et al. 2018; Favre et al. 2019). In particular, vertical gas flows have been identified in multiple disks with suspected planets, which suggests that perturbations in the vertical gas structure of disks may be common planetary signposts (e.g., Teague et al. 2019; Yu et al. 2021; Galloway-Sprietsma et al. 2023; Izquierdo et al. 2023). However, the direct detection of planets embedded in disks remains difficult and only a handful of such systems have thus far been robustly identified (Keppler et al. 2018; Haffert et al. 2019; Currie et al. 2022; Hammond et al. 2023; Wagner et al. 2023), making it difficult to conclusively assess the potential influence of giant planets on disk structure. Of these systems, PDS 70 is the only source whose protoplanetary disk is at a favorable inclination (51°; Keppler et al. 2019) to allow for a direct view of its

vertical gas distribution and thus provides a unique opportunity to understand if and how forming planets alter disk vertical structure.

PDS 70 is an ~5 Myr-old, K7 star (Müller et al. 2018) at a distance of 112 pc (Gaia Collaboration et al. 2021) in the Upper Centaurus Lupus association (Pecaut & Mamajek 2016) that hosts a gas-rich protoplanetary disk (Long et al. 2018; Facchini et al. 2021) with two, still-accreting giant planets PDS 70b and 70c that have been directly imaged at multiple wavelengths (Keppler et al. 2018; Müller et al. 2018; Wagner et al. 2018; Christiaens et al. 2019; Haffert et al. 2019; Mesa et al. 2019; Follette et al. 2023). PDS 70b and 70c are located at $\approx 22$ and $\approx 34$ au, respectively, in a near 2:1 mean motion resonance (Bae et al. 2019; Wang et al. 2021). Although the precise planet properties remain considerably uncertain, due in part to the presence of circumplanetary dust (e.g., Isella et al. 2019; Stolker et al. 2020; Benisty et al. 2021), they must be sufficiently massive (i.e., a few $M_{Jup}$) to have carved a central gas and dust cavity seen in molecular gas, millimeter continuum, and NIR/scattered light (Hashimoto et al. 2012, 2015; Keppler et al. 2019; Facchini et al. 2021; Portilla-Revelo et al. 2022, 2023).

The relative proximity and inclined orientation of the PDS 70 disk makes it possible to spatially resolve emission arising from elevated regions above and below the disk midplane using high-angular-resolution observations from the

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

<footnote n="1"><sup>1</sup>Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA</footnote>

<footnote n="2"><sup>2</sup>Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

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

<footnote n="4"><sup>4</sup>Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="5"><sup>5</sup>Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy</footnote>

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

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

<footnote n="8"><sup>8</sup>Department of Physics and Astronomy, Rice University, 6100 Main Street, MS-108, Houston, TX 77005, USA</footnote>

<footnote n="9"><sup>9</sup>Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands</footnote>

<footnote n="10"><sup>10</sup>NASA Hubble Fellowship Program Sagan Fellow.</footnote>
<table><caption>**Table 1**
Details of ALMA LB Observations of PDS 70</caption><thead><tr><th rowspan="2">UT Date</th><th rowspan="2">Ants.</th><th rowspan="2">Int. (min)</th><th rowspan="2">Baselines (m)</th><th rowspan="2">Res. (arcsec)</th><th rowspan="2">M.R.S. (arcsec)</th><th rowspan="2">PWV (mm)</th><th colspan="2">Calibrators</th></tr><tr><th>Flux/Bandpass</th><th>Phase</th></tr></thead><tbody><tr><td>2021-07-14</td><td>44</td><td>75.9</td><td>15.0-3396.4</td><td>0.1</td><td>1.8</td><td>1.7</td><td>J1427-4206</td><td>J1352-4412</td></tr><tr><td>2021-07-15</td><td>29</td><td>75.6</td><td>15.0-3396.4</td><td>0.1</td><td>1.8</td><td>1.4</td><td>J1427-4206</td><td>J1352-4412</td></tr><tr><td>2021-07-15</td><td>43</td><td>74.6</td><td>15.0-3396.4</td><td>0.1</td><td>1.8</td><td>1.1</td><td>J1427-4206</td><td>J1352-4412</td></tr><tr><td>2021-07-16</td><td>46</td><td>76.4</td><td>15.0-3638.2</td><td>0.1</td><td>2.0</td><td>1.2</td><td>J1427-4206</td><td>J1352-4412</td></tr></tbody></table>

Atacama Large Millimeter/submillimeter Array (ALMA) (e.g., Rosenfeld et al. 2013). This then allows for the direct extraction of emission surfaces of vertically extended molecular lines using similar techniques that have now been applied to a substantial number of mid-inclination disks (e.g., Pinte et al. 2018; Teague et al. 2019; Izquierdo et al. 2021, 2023; Law et al. 2021a, 2022, 2023; Rich et al. 2021; Paneque-Carreño et al. 2022, 2023; Stapper et al. 2023). Using these methods for the PDS 70 disk, we can map out the line-emitting heights and conduct a detailed search for vertical perturbations, such as those driven by the embedded planets PDS 70b and 70c.

By combining multiple optically thick molecular lines, which trace different heights in the disk, we can also empirically derive the full two-dimensional (2D) temperature distribution (e.g., Dartois et al. 2003; Rosenfeld et al. 2013; Pinte et al. 2018; Law et al. 2021a, 2023; Leemker et al. 2022). Having a well-constrained thermal structure is crucial as it allows for the determination of important volatile snowline locations, which have a direct impact on the expected atmospheric composition of embedded planets (e.g., Öberg et al. 2011; Mordasini et al. 2016; Mollière et al. 2022). Moreover, disk gas temperatures are vital inputs to numerical simulations and thermochemical models (e.g., Bae & Zhu 2018; Calahan et al. 2021; Zhang et al. 2021) and are required to infer the radial and vertical locations of molecular species, e.g., large organic molecules, that are otherwise too faint to directly map out (e.g., Ilee et al. 2021). This is especially critical to determine what types of chemical reservoirs are accessible to be accreted by PDS 70b and 70c and to more generally connect planetary atmospheric composition with that of the disk environment.

Here, we extract line-emission surfaces in the PDS 70 disk using high-spatial-resolution observations of CO isotopologues and HCO$^{+}$. In Section 2, we describe the observations, imaging, and emission-surface-extraction process. We present the derived emission surfaces along with radial and vertical temperature profiles in Section 3. In Section 4, we explore the origins of the observed disk vertical structure and estimate the CO snowline location based on the derived thermal structure. We summarize our conclusions in Section 5.

## 2. Observations and Analysis

### 2.1. Observational Details

We present long-baseline (LB) observations associated with the ALMA program (2019.1.01619.S; PI: S. Facchini) and Table 1 lists a summary of observation dates and details. Here, we restrict our analysis to a subset of the lines presented in the chemical inventory of Facchini et al. (2021) based on the short-baseline (SB) data from the same program, namely, the $J=2-1$ line of the $^{12}\text{CO}$, $^{13}\text{CO}$, and $\text{C}^{18}\text{O}$ isotopologues, which are

sufficiently bright and vertically extended to allow for extraction of their vertical structure. Full details about these SB data, including the spectral setup, are presented in Facchini et al. (2021) and the LB data for the remaining molecules will be published in forthcoming papers. We also make use of existing archival observations of the $^{12}\text{CO}$ $J=3-2$ and HCO$^{+}$ $J=4-3$ lines in the PDS 70 disk compiled from ALMA programs 2015.1.00888.S (PI: E. Akiyama) and 2017.A.00006. S (PI: M. Keppler), which were presented in Long et al. (2018), Keppler et al. (2019), and Facchini et al. (2021). We also include unpublished LB data of the same lines from 2018. A.00030.S (PI: M. Benisty). See Benisty et al. (2021) for additional details about each of these programs. We include HCO$^{+}$ in our analysis, due to its vertically elevated emission and the availability of existing high-angular-resolution line data. Since HCO$^{+}$ typically shows bright, optically thick emission in disks (e.g., Booth et al. 2019; Huang et al. 2020), it is thus also useful for mapping the thermal structure of PDS 70 in combination with CO isotopologues.

### 2.2. Self-calibration and Imaging

The three $J=2-1$ CO isotopologue lines are included in one spectral setup of the Facchini et al. (2021) program. For these Band 6 data, we followed the same self-calibration procedure, in line with Andrews et al. (2018) and Benisty et al. (2021). The self-calibration was performed with CASA v5.8 (McMullin et al. 2007; CASA Team et al. 2022). We first realigned the data of the individual Execution Blocks (EBs) to a common disk center with the fixvis and fixplanets tasks. The center was computed by fitting an ellipse to the outer ring (also clearly discernible in the SB data, see Facchini et al. 2021). This approach had previously performed well for the self-calibration of the Band 7 data of the PDS 70 disk by Benisty et al. (2021). After flagging the lines listed in Facchini et al. (2021) within $\pm 15$ km s$^{-1}$ from the line center, we spectrally averaged the data into 250 MHz channels to create pseudo-continuum ms tables. Individual EBs were rescaled in amplitude to provide the same disk flux density, after verifying that phase de-coherence was not artificially reducing it. The SB data were self-calibrated first, and then concatenated to the LB data before running a new self-calibration. The tclean model was used as the source model for each step, where tclean was performed over an elliptical mask with a $1''5$ radius and inclination and position angle as in Facchini et al. (2021). Four rounds of phase-only and one round of amplitude self-calibration were performed on the SB data, and three rounds of phase-only and one of amplitude were used on the LB data, respectively. Before amplitude calibration, we applied tclean to the data with a deep $1\sigma$ threshold. Using Briggs weighting with robust = 0.5, the synthesized beam was $0''15 \times 0''12$ (PA = $-82^\circ9$). The resulting 233 GHz continuum image had a flux density of
<table><caption>**Table 2**
Image Cube Properties</caption><thead><tr><th>Transition</th><th>Beam<br>(arcsec × arcsec, deg)</th><th>JvM ϵ<sup>a</sup></th><th>robust</th><th>Chan. δν<br>(km s<sup>-1</sup>)</th><th>rms<br>(mJy beam<sup>-1</sup>)</th><th>R<sub>size</sub><sup>b</sup><br>(au)</th><th>ALMA<br>Project Code(s)</th></tr></thead><tbody><tr><td><sup>12</sup>CO J = 2–1</td><td>0.13 × 0.10, -83.7</td><td>0.61</td><td>0.25</td><td>0.30</td><td>0.93</td><td>209 ± 4</td><td>P2019</td></tr><tr><td><sup>13</sup>CO J = 2–1</td><td>0.15 × 0.12, -80.0</td><td>0.46</td><td>0.5</td><td>0.20</td><td>0.62</td><td>167 ± 5</td><td>P2019</td></tr><tr><td>C<sup>18</sup>O J = 2–1</td><td>0.16 × 0.13, -82.4</td><td>0.46</td><td>0.5</td><td>0.40</td><td>0.42</td><td>142 ± 7</td><td>P2019</td></tr><tr><td><sup>12</sup>CO J = 3–2</td><td>0.07 × 0.06, 87.1</td><td>0.38</td><td>1.0</td><td>0.43</td><td>0.26</td><td>210 ± 2</td><td>P2015, P2017, P2018</td></tr><tr><td>HCO<sup>+</sup> J = 4–3</td><td>0.05 × 0.04, 69.3</td><td>0.30</td><td>1.0</td><td>0.43</td><td>0.31</td><td>165 ± 3</td><td>P2015, P2018</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The ratio of the CLEAN beam and dirty beam effective area used to scale image residuals to account for the effects of non-Gaussian beams. See Section 2.2 and Jorsater & van Moorsel (1995) and Czekala et al. (2021) for further details.

<sup>b</sup> Disk size ($R_{\text{size}}$) was computed as the radius that encloses 90% of the total disk flux.

58.3 mJy within the CLEAN mask (in agreement with Facchini et al. 2021), with an rms of 12.5 µJy beam⁻¹. We finally applied the gain solutions to the full spectral data.

We obtained archival Band 7 observations of the <sup>12</sup>CO J = 3–2 and HCO<sup>+</sup> J = 4–3 lines from the data sets presented in Benisty et al. (2021), for which we applied the same gain solutions to the line spectral windows. For both the Band 6 and Band 7 lines, we continuum subtracted the data with the contsub task using a first-order polynomial.

We imaged the five lines discussed in this paper with CASA v6.2. The channel maps were restored using tclean with a threshold of 3.5σ. Table 2 lists the imaging parameters of all cubes. We used Keplerian masks for the channels generated from the keplerian_mask code (Teague 2020), where we assumed the disk to be extended 3" in radius, and we conservatively assumed that all emission is as elevated as <sup>12</sup>CO. Given the focus of the paper, extracting accurate brightness temperatures is imperative. We thus applied the so-called JvM correction to the final images (Jorsater & van Moorsel 1995; Czekala et al. 2021) to more accurately recover the low surface brightness intensities in the disk outer regions and of weaker lines (e.g., C<sup>18</sup>O J = 2–1). The effective rms, which can be underestimated when applying the JvM correction (Casassus & Cárcamo 2022), is not relevant for the analysis of this paper. The full image cube channel maps are available in Appendix A. We also imaged the non-continuum-subtracted line data by adopting the same imaging parameters as the line-only emission image cubes. The non-continuum-subtracted image cubes are required for the calculation of gas temperatures (Section 3.4). Overall, the ability to extract line-emitting heights depends on having sufficiently high angular resolution, spectral resolution, and line sensitivities. Considering this, we selected the imaging parameters that best suited surface extraction for this work. Table 2 lists the properties of the image cubes selected for analysis of the vertical line-emission structure.

Figure 1 shows an overview of the CO isotopologue and HCO<sup>+</sup> line-emission velocity-integrated intensity, or “zeroth moment,” maps, and the 344 GHz continuum emission. The continuum image is taken from Isella et al. (2019). We generated zeroth moment maps of line emission from the image cubes using bettermoments (Teague & Foreman-Mackey 2018) with the same Keplerian masks employed during CLEANing and with no flux threshold for pixel inclusion to ensure accurate flux recovery.

## 2.3. Emission-surface Extraction

We derived vertical emission heights from the line-emission image cubes using the disksurf (Teague et al. 2021) Python code, which is based on the methodology originally presented in Pinte et al. (2018). We closely followed the methods outlined in Law et al. (2021a), which we briefly summarize below.

We first adopted a PA = 160°.4 and inc = 51°.7 for the PDS 70 disk (Keppler et al. 2019; Facchini et al. 2021) and then extracted a deprojected radius *r*, emission height *z*, surface brightness $I_v$, and channel velocity *v* for each pixel associated with the emitting surface. We then employed an iterative fitting approach, which further refines the extracted surfaces using one to five subsequent iterations depending on the line. To ensure the robustness of the extracted points, we also filtered those pixels with unphysical *z/*r values based on the expected disk structure and removed those points with low surface brightnesses (2–3 × rms). At each step, we visually confirmed the fidelity of the derived emission surfaces. To further reduce scatter in the extracted surfaces, we also performed each of the same two types of binning (radially binned and moving averages) as in Law et al. (2021a) with half-beam spacings.

We then fitted parametric models to all line-emission surfaces using an exponentially tapered power law, which captures the inner flared surfaces and the plateau and turnover regions at larger radii (e.g., Teague et al. 2019). We adopted the following functional form:

$$z(r) = z_0 \times \left( \frac{r - r_{\text{cavity}}}{l''} \right)^{\phi} \times \exp \left( - \left[ \frac{r}{r_{\text{taper}}} \right]^{\psi} \right) \quad (1)$$

This is similar to the form employed by Law et al. (2021a) but here we include an additional term $r_{\text{cavity}}$ to describe the inner gas cavity of the PDS 70 disk.

We used the Monte Carlo Markov Chain (MCMC) sampler implemented in emcee (Foreman-Mackey et al. 2013) to estimate the posterior probability distributions for: $z_0$, $\phi$, $r_{\text{taper}}$, $r_{\text{cavity}}$, and $\psi$. Each ensemble used 64 walkers with 1000 burn-in steps and an additional 500 steps to sample the posterior distribution function. The posteriors were approximately Gaussian with no significant degeneracies between parameters. Table 3 shows the median values of the posterior distribution, with uncertainties given as the 16th and 84th percentiles, as well as the radial range, $r_{\text{fit, max}}$, considered for each surface.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** CO isotopologue and HCO⁺ zeroth moment maps and continuum images (from top left to bottom right) of the PDS 70 disk. Panels for each disk have the same field of view. The locations of the planets PDS 70b and 70c are marked in each panel (Wang et al. 2021). Color stretches were individually optimized and applied to each panel to increase the visibility of outer disk structure. The synthesized beam and a scale bar indicating 50 au is shown in the lower left and right corner, respectively, of each panel. The 344 GHz continuum is from Isella et al. (2019). Table 2 lists details about each image.</figcaption></figure>

<table><caption>**Table 3**
Parameters for CO Isotopologue and HCO⁺ Emission-surface Fits</caption><thead><tr><th rowspan="2">Line</th><th rowspan="2">r<sub>fit,max</sub> (arcsec)</th><th colspan="5">Exponentially Tapered Power Law</th><th rowspan="2">Char. z/r<sup>a</sup></th></tr><tr><th>z<sub>0</sub> (arcsec)</th><th>φ</th><th>r<sub>cavity</sub> (arcsec)</th><th>r<sub>taper</sub> (arcsec)</th><th>ψ</th></tr></thead><tbody><tr><td><sup>12</sup>CO J = 2–1</td><td>1.80</td><td>0.41<sup>+0.02</sup><sub>-0.02</sub></td><td>0.51<sup>+0.08</sup><sub>-0.06</sub></td><td>0.37<sup>+0.02</sup><sub>-0.05</sub></td><td>1.27<sup>+0.03</sup><sub>-0.02</sub></td><td>5.74<sup>+0.97</sup><sub>-0.84</sub></td><td>0.32<sup>+0.01</sup><sub>-0.06</sub></td></tr><tr><td><sup>13</sup>CO J = 2–1</td><td>1.45</td><td>0.41<sup>+1.30</sup><sub>-0.24</sub></td><td>1.29<sup>+1.16</sup><sub>-0.87</sub></td><td>0.28<sup>+0.21</sup><sub>-0.21</sub></td><td>0.74<sup>+0.52</sup><sub>-0.37</sub></td><td>1.61<sup>+1.28</sup><sub>-0.75</sub></td><td>0.09<sup>+0.04</sup><sub>-0.02</sub></td></tr><tr><td><sup>12</sup>CO J = 3–2</td><td>1.90</td><td>0.43<sup>+0.03</sup><sub>-0.05</sub></td><td>0.63<sup>+0.05</sup><sub>-0.13</sub></td><td>0.35<sup>+0.00</sup><sub>-0.02</sub></td><td>1.34<sup>+0.02</sup><sub>-0.03</sub></td><td>5.40<sup>+4.09</sup><sub>-1.80</sub></td><td>0.32<sup>+0.02</sup><sub>-0.04</sub></td></tr><tr><td>HCO<sup>+</sup> J = 4–3</td><td>1.55</td><td>0.70<sup>+0.20</sup><sub>-0.23</sub></td><td>1.48<sup>+0.89</sup><sub>-0.57</sub></td><td>0.29<sup>+0.11</sup><sub>-0.18</sub></td><td>0.65<sup>+0.14</sup><sub>-0.10</sub></td><td>3.00<sup>+0.91</sup><sub>-0.70</sub></td><td>0.18<sup>+0.02</sup><sub>-0.09</sub></td></tr></tbody></table>

**Note.**

## 3. Results

### 3.1. Emission Surfaces in the PDS 70 Disk

Figure 2 shows the derived emission surfaces in the PDS 70 disk. All of the CO isotopologue lines show characteristic emission-surface profiles, i.e., a sharply rising inner component followed by subsequent plateau and turnover phases toward larger radii (e.g., Teague et al. 2019; Law et al. 2021a). Due to the high sensitivity and angular resolution of the observations combined with the iterative surface-fitting procedure of *disksurf*, we are also able to confidently track emission heights even at large disk radii.

Both the $J=3-2$ and $J=2-1$ lines of $^{12}\text{CO}$ show altitudes of $z/r \approx 0.3$ and reach a maximum height of approximately 40 au at a radius of 150 au. The $^{12}\text{CO}$ $J=3-2$ heights are generally consistent with those derived in Keppler et al. (2019), who fit a single power-law surface profile to the Keplerian rotation map but did not extract individual heights. The $^{12}\text{CO}$

$J=3-2$ emission surface shows evidence for vertical variation due to the cavity wall at $\approx 45$ au. From a radius of $\approx 15-40$ au, the surface is quite flat ($z/r \approx 0.1$) and located near the midplane, but then abruptly increases height at a radius of $\approx 45$ au up to $z/r \approx 0.3$. We are unable to extract emission heights for the $^{12}\text{CO}$ $J=2-1$ line in the inner 50 au (square box in Figure 2) despite $^{12}\text{CO}$ emission being clearly detected in the inner disk (Figure 1). This is likely due to the larger beam size of the $J=2-1$ observations (three times larger in beam area) rather than an intrinsic difference in emitting heights, especially since the $J=2-1$ and $J=3-2$ $^{12}\text{CO}$ lines have typically been observed to originate from the same heights in other disks (Law et al. 2023). The $^{13}\text{CO}$ $J=2-1$ line originates from a deeper layer in the disk ($z/r \approx 0.1$). We do not resolve the vertical structure of the $^{18}\text{O}$ emission, which suggests that $^{18}\text{O}$ arises at, or near, the disk midplane. In addition to CO isotopologues, we also measured the HCO$^+$ $J=4-3$ emitting surface, which lies at $z/r \approx 0.2$. Beyond $\approx 100$ au, the HCO$^+$

<footnote>a Characteristic $z/r$ is computed as the mean and 16th to 84th percentile range within 80% of the $r_{cutoff}$ as in Law et al. (2022, 2023).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** CO isotopologue and HCO⁺ emission surfaces for the PDS 70 disk. Large gray points show radially binned surfaces and small, light gray points represent individual measurements. The orange lines show the exponentially tapered power-law fits from Table 3. The solid lines show the radial range used in the fitting, while the dashed lines are extrapolations. Lines of constant $z/r$ from 0.1 to 0.5 are shown in gray. Vertical substructures and the locations of PDS 70b and 70c are marked in each panel (Wang et al. 2021). The lack of $^{12}$CO J = 2–1 heights in the inner 50 au (square box) is due to insufficient angular resolution rather than the absence of emission at these radii. The FWHM of the major axis of the synthesized beam is shown in the bottom right corner of each panel. The emission surfaces shown in this figure are available as Data behind the Figure.</figcaption></figure>

(The data used to create this figure are available.)

surface shows an extended plateau of flat, near-midplane emission ($z/r \lesssim 0.1$), which is not observed in any other line. At $< 100$ au, the HCO⁺ arises from higher altitudes than $^{13}$CO, but at $> 100$ au, HCO⁺ is coming from lower altitudes than $^{13}$CO.

In addition to the cavity wall seen in $^{12}$CO J = 3–2, each of the $^{12}$CO and $^{13}$CO surfaces shows a shallow vertical dip between $\approx 90$ and 120 au. We confirmed that each of substructures are also detected in emission surfaces extracted using the non-continuum-subtracted image cubes and with no filtering of low-surface-brightness points, and thus, are not artifacts of the continuum subtraction or surface-extraction process. To quantify the significance of these features, including the cavity wall, we follow the fitting procedure outlined in Law et al. (2021a). In brief, we remove a local baseline around each vertical substructure and then fit for its radial position, width, and depth. Table 4 shows the properties of each substructure, which are labeled with a Z to indicate that they are vertical variations followed by their radial location rounded to the nearest whole number in astronomical units. The cavity wall at Z42 is the deepest vertical feature, followed by the Z95 dip in $^{13}$CO, while the dips at Z119 and Z113 in the $^{12}$CO surfaces are the most shallow.

In Figure 3, we overlay the inferred emission surfaces on peak line intensity maps to better illustrate their 3D geometries. We generated all peak-intensity maps with the “quadratic” method of better moments using the full Planck function.

<table><caption>**Table 4**
Properties of Vertical Substructures</caption><thead><tr><th>Line</th><th>Feature</th><th>r<sub>0</sub> (au)</th><th>Width (au)</th><th>Δz (au)</th><th>Depth</th></tr></thead><tbody><tr><td><sup>12</sup>CO J = 2–1</td><td>Z119</td><td>119 ± 1</td><td>14 ± 7</td><td>2 ± 1</td><td>0.07 ± 0.06</td></tr><tr><td rowspan="2"><sup>12</sup>CO J = 3–2</td><td>Z42</td><td>42 ± 0.1</td><td>21 ± 4</td><td>6 ± 0.2</td><td>0.48 ± 0.07</td></tr><tr><td>Z113</td><td>113 ± 1</td><td>25 ± 6</td><td>3 ± 1</td><td>0.10 ± 0.08</td></tr><tr><td><sup>13</sup>CO J = 2–1</td><td>Z95</td><td>95 ± 2</td><td>20 ± 0.3</td><td>4 ± 1</td><td>0.28 ± 0.14</td></tr></tbody></table>

## 3.2. Radial Disk Chemical Substructure

While several molecular lines considered here have been previously observed toward the PDS 70 disk (e.g., Long et al. 2018; Facchini et al. 2021), these new data represent improvements in both angular resolution and sensitivity of more than a factor of two. Thus, in addition to mapping the vertical gas structure, we can also constrain the presence of small-scale radial substructures at the highest spatial resolution to date. To do so, we computed azimuthally averaged radial profiles using the radial_profile function in the GoFish Python package (Teague 2019) to deproject the zeroth moment maps. During deprojection, we incorporated the derived emission surfaces listed in Table 3. This is particularly important for highly elevated surfaces, e.g., $^{12}$CO, to derive accurate radial locations of line-emission substructures (e.g., Law et al. 2021b; Rosotti et al. 2021). Figure 4 shows the resultant radial profiles.

The location of line-emission substructures are labeled according to their radial location rounded to the nearest au
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Peak-intensity maps with overlaid contours showing the fitted emission surfaces, as listed in Table 3. No surface could be derived for the C<sup>18</sup>O J = 2−1 line, which is consistent with emission originating from near the midplane. The synthesized beam and a scale bar indicating 50 au is shown in the lower left and right corner, respectively, of each panel.</figcaption></figure>

following standard nomenclature (e.g., Huang et al. 2018; Law et al. 2021b). The <sup>12</sup>CO lines show centrally peaked profiles with a dip at 21 au, followed by an emission peak between 40 and 50 au and an extended plateau out to ≈ 250 au. The higher resolution <sup>12</sup>CO J = 3−2 profile resolves an inner ring at ≈5 au, indicating the possible presence of small-scale inner-disk substructure. The <sup>13</sup>CO and C<sup>18</sup>O J = 2−1 lines instead show a central gap and broad ring at 67 au. The HCO<sup>+</sup> profile is more similar to that of the <sup>12</sup>CO lines with a dip at 25 au and ring at 50 au. However, HCO<sup>+</sup> is much more sharply centrally peaked with no indications of a central dip in the inner few astronomical units, despite the HCO<sup>+</sup> having the highest spatial resolution. Overall, the radial morphologies are consistent with existing lower resolution observations (Long et al. 2018; Facchini et al. 2021), namely that <sup>12</sup>CO and HCO<sup>+</sup> peak interior to the millimeter dust ring and <sup>13</sup>CO and C<sup>18</sup>O are colocated with the millimeter continuum peak. As discussed in Facchini et al. (2021), these differences are likely due to varying line optical depths and the presence of high temperatures and illumination at the edge of the cavity wall.

We also computed the disk size of each emission line by determining the radius in which 90% of the total flux is contained (e.g., Tripathi et al. 2017; Ansdell et al. 2018). Table 2 lists the derived values. The <sup>12</sup>CO lines have have sizes of ≈ 200 au, while <sup>13</sup>CO and HCO<sup>+</sup> have sizes of ≈ 165 au. The C<sup>18</sup>O is the most compact line at ≈ 140 au. The millimeter continuum edge of the PDS 70 disk is approximately 100 au. This implies gas-to-dust sizes of ≈ 1.4−2, which are typical, albeit on the lower end, for large, resolved disks (e.g., Law et al. 2021b; Sanchis et al. 2021; Long et al. 2022).

### 3.3. Comparison with NIR Scattering Surface

The PDS 70 disk has been extensively observed at NIR wavelengths, which revealed an elliptical ring, inner-disk component, and asymmetrical feature to the northwest of the star due to either a double ring or planet-disk interactions (e.g., Hashimoto et al. 2012; Keppler et al. 2018; Müller et al. 2018; Mesa et al. 2019; Juillard et al. 2022). Here, we only focus on the scattering surface in the outer disk, which is tracing the small dust-grain population, and compare this to the derived line-
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Deprojected radial integrated intensity profiles lines and the 344 GHz continuum in the PDS 70 disk. Shaded regions show the $1\sigma$ uncertainty. Solid lines mark emission rings and dotted lines indicate gaps. The FWHM of the synthesized beam is shown by a horizontal bar in the upper right corner of each panel.</figcaption></figure>

emission surfaces. Only a handful of disks have independent measurements of both NIR and gas heights, and as a result, the general relationship between the vertical distribution of gas and small dust in disks is not yet well established.

Figure 5 shows the power-law NIR surface derived in Keppler et al. (2018) versus the extracted molecular gas heights. We also directly measured the height of the scattered light ring at ≈54 au in the PDS 70 disk by fitting an ellipse to the peak flux using the 2016 SPHERE IRDIS *J*-band, $Q_\phi$ polarimetric image from Keppler et al. (2018). To do so, we followed the same approach outlined in Rich et al. (2021) and described in detail in Appendix C of Law et al. (2023). We derived a height of $4.5 \pm 0.1$ au at a radius of $53.7 \pm 0.1$ au, which agrees well with the overall power-law trend from Keppler et al. (2018) and is shown as a diamond marker in Figure 5. The NIR surface is located at the same vertical height as the $^{13}$CO surface over the entire radial range of the NIR disk. Existing measurements in other disks show that the NIR surface typically lies between the $^{12}$CO and $^{13}$CO layers, but with considerable diversity in the relative heights of the small dust grains (e.g., Rich et al. 2021; Law et al. 2022, 2023; Paneque-Carreño et al. 2023). In particular, the well-studied IM Lup and HD 163296 disks show comparable heights in $^{12}$CO gas and small dust grains within 100 au, but significantly more elevated $^{12}$CO emitting heights at larger radii (Rich et al. 2021; Paneque-Carreño et al. 2023). Given that NIR emission in PDS 70 only extends to ≈115 au (Keppler et al. 2018), it is thus possible that the PDS 70 disk shows slightly lower than expected small dust heights compared to that of $^{12}$CO. However, we caution that only a small number of disks have had their CO isotopologue and small dust heights jointly mapped and additional disk observations are required to better place PDS 70 in context.

### 3.4. Disk Thermal Structure

We extracted radial (Section 3.4.1) and vertical (Section 3.4.2) gas temperatures for the PDS 70 disk using $^{12}$CO, $^{13}$CO, and HCO$^+$. To do so, we followed the procedures of Law et al. (2021a). We assume that all lines are optically thick and in LTE,

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Emission surfaces of $^{12}$CO and $^{13}$CO *J* = 2–1 (top) and $^{12}$CO *J* = 3–2 and HCO$^+$ *J* = 4–3 (bottom) vs. the NIR scattering surface from Keppler et al. (2018) in the PDS 70 disk. The red marker shows the individual height measurement of the PDS 70 NIR ring derived in this work. The lines show the moving-average surfaces and gray shaded regions show the $1\sigma$ uncertainty. The FWHM of the major axis of the synthesized beam is shown in the bottom right corner of each panel.</figcaption></figure>

which is expected to be the case at the typical densities and temperatures of protoplanetary disks (e.g., Weaver et al. 2018). We also assume that the line emission fills the beam due to the high angular resolution of the observations. We then repeated the
<figure><img src="image_9.png" /><figcaption>**Figure 6.** CO, ¹³CO, and HCO⁺ radial brightness temperature profiles of the PDS 70 disk. These profiles represent the mean temperatures computed by radially binning the individual measurements, similar to the procedure used to compute the radially binned surfaces (see Section 3.4). Vertical lines show the 1σ uncertainty, given as the standard deviation of the individual measurements in each bin. The solid gray lines show the power-law fits from Table 5, while the dashed lines are extrapolations. The FWHM of the major axis of the synthesized beam is shown in the bottom right corner of each panel.</figcaption></figure>

surface extractions as in Section 2.3 on the non-continuum-subtracted image cubes, which ensures that we do not underestimate the line intensity along lines of sight containing strong continuum emission (e.g., Boehler et al. 2017). We then converted the peak surface brightness $I_ν$ of each extracted pixel to a brightness temperature using the full Planck function and assume the resulting brightness temperature is equal to the local gas temperature.

All subsequent radial and 2D gas temperature distributions are taken directly from individual surface measurements, rather than from mapping peak brightness temperatures back onto derived emission surfaces (i.e., Figure 3) or radially deprojecting peak-intensity maps (see Appendix B). Hence, we only consider the brightness temperature of those pixels with derived emission heights.

### 3.4.1. Radial Temperature Profiles

Figure 6 shows the radial temperature distributions along the emission surfaces. The ¹²CO lines trace the warmest temperatures, followed by HCO⁺ and then ¹³CO. For all lines, the temperatures generally increase toward the central star, except in the inner ≈ 60–70 au where the temperatures of the CO lines plateau or decrease. This is unsurprising since the CO lines all either show a central cavity or deep gap in the inner disk, which means that at these radii, the lines are likely no longer optically thick and thus do not trace the true gas temperature. We also confirmed that the radial temperatures extracted directly from the emission surfaces in Figure 6 show excellent agreement with the peak-intensity profiles generated from Figure 3, which are shown in Appendix B.

The $J=2-1$ and $J=3-2$ lines of ¹²CO are generally comparable in temperature, but the $J=3-2$ line has a peak temperature of ≳ 40 K that is a few K warmer than that of the $J=2-1$ line, which is likely due to beam dilution at smaller radii, since the $J=2-1$ beam is more than twice as large than that of the $J=3-2$ line. While the ¹³CO and HCO⁺ temperatures profiles appear smooth, both ¹²CO profiles show

<table><caption>**Table 5**

Radial Temperature Profile Fits</caption><thead><tr><th>Line</th><th>r<sub>fit,in</sub><br>(au)</th><th>r<sub>fit,out</sub><br>(au)</th><th>T<sub>100</sub> (K)</th><th>q</th><th>Feat.<sup>a</sup></th></tr></thead><tbody><tr><td>¹²CO J = 2–1</td><td>72</td><td>204</td><td>32 ± 0.3</td><td>0.57 ± 0.03</td><td>D115</td></tr><tr><td>¹³CO J = 2–1</td><td>75</td><td>164</td><td>18 ± 0.5</td><td>0.44 ± 0.10</td><td></td></tr><tr><td>¹²CO J = 3–2</td><td>70</td><td>216</td><td>31 ± 0.2</td><td>0.79 ± 0.02</td><td>D127</td></tr><tr><td>HCO<sup>+</sup> J = 4–3</td><td>46</td><td>161</td><td>22 ± 0.2</td><td>0.43 ± 0.02</td><td></td></tr></tbody></table>

**Note.**

the presence of a dip at ≈ 120 au. This dip occurs exterior to the continuum edge and corresponds to a small vertical dip in molecular emission height at the same radius (Section 3.1). This temperature decrease could either be due to locally reduced gas surface density, in which the emission surface traces a deeper and thus cooler layer, or due to a change in heating or radiation properties related to the continuum edge. To better quantify these temperature dips, we use the same method as for the vertical substructures in Section 3.1. The dip in temperature for ¹²CO $J=2-1$ occurs at 115 au, while in ¹²CO $J=3-2$, it is at 127 au. Both features have comparable widths (≈ 20–30 au) and depths (≤10%). A full listing of temperature substructure properties is provided in Appendix B.

To better characterize the radial temperature profiles, we fitted each line with a power-law profile as:

$$T = T_{100} \times \left( \frac{r}{100 \text{ au}} \right)^{-q}, \qquad (2)$$

where $q$ is the slope and $T_{100}$ is the brightness temperature at 100 au. Table 5 shows the fitted parameters. Overall, the profiles have slopes of $q \approx 0.4–0.8$, which is generally consistent with other disks (e.g., Law et al. 2022), with ¹²CO $J=3-2$ having the steepest radial temperature profile.

<footnote>a Local temperature dips (D) labeled according to their approximate radial location in au. See Appendix B for further details.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** 2D temperature profiles of the ¹²CO and ¹³CO *J* = 2–1 (left) and ¹²CO *J* = 3–2 and HCO⁺ *J* = 4–3 (right) lines in the PDS 70 disk. Points are those from the binned surfaces and error bars are the 1σ uncertainties in *z*. Temperature measurements with radii near the emission gap (Figure 6) are marked with hollow markers and not used in the 2D temperature fits. The 2D temperature profiles shown in this figure are available as Data behind the Figure. (The data used to create this figure are available.)</figcaption></figure>

### 3.4.2. 2D Temperature Profiles

Figure 7 shows the thermal structure of the CO isotopologues and HCO⁺ emitting layers as a function of (*r*, *z*). By combining all lines with derived surfaces, which trace different heights in the disk, we can construct a 2D model of the temperature distribution. Following Law et al. (2021a), we fit a two-layer temperature model in which disk midplane and atmosphere temperature are power laws smoothly connected with a hyperbolic tangent function (Dartois et al. 2003; Dullemond et al. 2020):

$$T_{\text{mid}}(r) = T_{\text{mid,0}}(r/100 \text{ au})^{q_{\text{mid}}} \quad (3)$$

$$T_{\text{atm}}(r) = T_{\text{atm,0}}(r/100 \text{ au})^{q_{\text{atm}}} \quad (4)$$

$$T^4(r, z) = T_{\text{mid}}^4(r) + \frac{1}{2} \left[ 1 + \tanh \left( \frac{z - \alpha z_q(r)}{z_q(r)} \right) \right] T_{\text{atm}}^4(r), \quad (5)$$

where $z_q(r) = z_0(r/100 \text{ au})^\beta$. The $\alpha$ parameter defines the height at which the transition in the tanh vertical temperature profile occurs, while $\beta$ describes how the transition height varies with radius.

We used the MCMC sampler in emcee (Foreman-Mackey et al. 2013) to estimate the following seven parameters: $T_{\text{atm,0}}$, $q_{\text{atm}}$, $T_{\text{mid,0}}$, $q_{\text{mid}}$, $\alpha$, $z_0$, and $\beta$. We used 256 walkers, which take 500 steps to burn in and then an additional 5000 steps to sample the posterior distribution function. Figure 8 shows the fitted 2D temperature profiles and Table 6 lists the fitted parameter values and uncertainties, which are given as the 50th, 16th, and 84th percentiles from the marginalized posterior distributions, respectively.

## 4. Discussion

### 4.1. Influence of Embedded Planets PDS 70b and 70c on Disk Vertical Structure

The presence of embedded massive protoplanets may be expected to either locally or perhaps globally alter the vertical gas structure of protoplanetary disks, which in turn, would affect the properties of the material available to be accreted by those planets. The PDS 70 system provides an ideal environment to search for such effects. Here, we first discuss

the potential presence of any local perturbations due to PDS 70b and 70c and then compare the overall vertical structure of the PDS 70 disk with other disks in the literature.

#### 4.1.1. Local Perturbations and Cavity Wall

In the ¹²CO *J* = 3–2 emission surface, we identify relatively flat (*z*/*r* ≈ 0.1) heights that overlap with the planet locations, i.e., within the central ≈ 45 au. This vertically flat region extends to the edge of the cavity wall, where the heights sharply rise again to *z*/*r* ≈ 0.3 (Figure 2). This is the same region that shows a deep drop in brightness temperature, reaching as low as ≈ 15 K (Figure 6).

The flat emitting heights observed within ≈ 45 au are consistent with expectations of a transition disk with a deep gas- and dust-depleted central cavity, which in the case of the PDS 70 disk has been carved out by PDS 70b and 70c (Bae et al. 2019; Portilla-Revelo et al. 2023). Line-emission surfaces are tracing regions where optical depths reach unity. Thus, in this region near the planets, the overall line optical depth decreases due to the reduced gas surface density and the surface is pushed deeper into the disk, i.e., close to the midplane, which naturally explains both the dip in vertical heights and the lower brightness temperatures. This is the first transition disk where this effect is directly observed in line-emitting heights (see Section 4.2 for more details). For all other lines besides ¹²CO *J* = 3–2, no estimates of vertical heights within the central cavity were possible due to insufficient angular resolution or intrinsically flat emitting heights likely due to the low gas surface densities.

#### 4.1.2. Outer Disk Vertical Structure of PDS 70 in Context

Although we do not see any clear local perturbations around PDS 70b or 70c beyond that of the cavity wall, the overall structure of the gas disk may still be influenced by the presence of both planets in the form of, e.g., dynamical or planet-disk interactions. In particular, here we focus on the outer disk vertical gas structure. One way to assess this is to compare the PDS 70 disk with others in the literature with similar constraints on their line-emitting heights.

To do this, we computed the characteristic *z*/*r* heights of each line in the PDS 70 disk following the procedure outlined in Law et al. (2022), i.e., mean *z*/*r* within the steeply rising portion of the surfaces, for consistent comparison with
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison of the measured temperatures (points) with the fitted 2D temperature structures (background) for the PDS 70 disk, as listed in Table 6. The <sup>12</sup>CO and <sup>13</sup>CO *J* = 2–1 (left) and <sup>12</sup>CO *J* = 3–2 and HCO<sup>+</sup> *J* = 4–3 (right) are shown in separate panels for visual clarity but the 2D temperature was fit to both sets of lines simultaneously. The same color scale is used for the data and fitted model in each panel. Points excluded from the fits are shown as hollow markers and contours show constant temperatures. The uncertainty of the temperature measurements, which is not shown here, can be found in Figure 7.</figcaption></figure>

<table><caption>**Table 6**
Summary of 2D Temperature Structure Fits</caption><thead><tr><th><i>T</i><sub>atm,0</sub> (K)</th><th><i>T</i><sub>mid,0</sub> (K)</th><th><i>q</i><sub>atm</sub></th><th><i>q</i><sub>mid</sub></th><th><i>z</i><sub>0</sub> (au)</th><th><i>α</i></th><th><i>β</i></th></tr></thead><tbody><tr><td>33<sup>+0.2</sup><sub>-0.2</sub></td><td>16<sup>+0.3</sup><sub>-0.3</sub></td><td>-0.95<sup>+0.02</sup><sub>-0.02</sub></td><td>-0.01<sup>+0.04</sup><sub>-0.04</sub></td><td>13<sup>+0.2</sup><sub>-0.2</sub></td><td>2.13<sup>+0.03</sup><sub>-0.03</sub></td><td>0.14<sup>+0.03</sup><sub>-0.03</sub></td></tr></tbody></table>

literature values. Table 3 lists the computed values. We first consider the CO isotopologues and then HCO<sup>+</sup>.

### 4.1.2.1 CO Isotopologues

Figure 9 shows the characteristic $z/r$ emission heights of <sup>12</sup>CO and <sup>13</sup>CO versus the stellar mass and disk size compiled from consistently derived heights (Law et al. 2023). We also include several expected scaling relations, i.e., assuming line-emission heights scale with gas-pressure scale heights, and a previously identified positive <sup>12</sup>CO-R<sub>12CO</sub> trend (for more details, see Law et al. 2022, 2023). With respect to the literature sample, the PDS 70 disk appears quite typical in terms of its <sup>12</sup>CO and <sup>13</sup>CO-emitting heights compared to disks with similar host stellar masses and disk sizes. This is true for all disks, or if we only consider the subsample of transition disks. Thus, the two planets in the PDS 70 disk do not seem to be substantially altering the heights of the line-emission surfaces via, e.g., planet-disk interactions. However, we note the caveat that many of the literature disks used as the comparison sample also have indirect evidence for the presence of embedded planets. Thus, the comparison of emitting heights in the PDS 70 disk versus those disks without planets is not necessarily as straightforward as Figure 9 suggests, but such a comparison nonetheless establishes that the PDS 70 disk has a vertical gas structure that is typical of the large, resolved disks thus far studied in detail with ALMA.

The presence of embedded planets could, for instance, have resulted in flatter emitting surfaces due to the removal or redistribution of gas via vertical flows (e.g., Teague et al. 2019; Yu et al. 2021; Galloway-Sprietsma et al. 2023), or alternatively, driven dynamical interactions that would inflate the disk gas vertical distribution (e.g., Montesinos et al. 2021; Kuo et al. 2022). However, the close-in radial locations of the PDS 70 planets may provide an explanation, in which the vertical disk gas distribution is decoupled from planet formation occurring in the inner few 10 s of astronomical units.

Several previous works (Izquierdo et al. 2023; Paneque-Carreño et al. 2023) have identified vertical variations in line emitting heights that are coherent across several tracers at the radial locations of suspected planetary companions or kinematic deviations. In the PDS 70 disk, we find vertical variations at radii between ≈ 95 and 120 au in both <sup>12</sup>CO and <sup>13</sup>CO. However, no evidence for kinematic deviations or additional embedded planets exists in the outer disk of PDS 70. Thus, it is more likely that these vertical substructures are instead related to the edge of the millimeter continuum disk due to, e.g., changing radiation fields at the dust edge.

### 4.1.2.2 HCO<sup>+</sup>

Only a few measurements of HCO<sup>+</sup> emitting heights exist in protoplanetary disks. Paneque-Carreño et al. (2023) derived heights of the *J* = 1–0 line in the MAPS disks (Öberg et al. 2021) and for the *J* = 4–3 line in the WaOph 6 disk, while Huang et al. (2020) constrained the *J* = 3–2 emitting surface in the GM Aur disk. The *J* = 1–0 lines showed typical heights of $z/r \lesssim 0.1$, but cannot be directly compared to PDS 70, because the difference in excitation properties between the *J* = 4–3 ($E_u \approx 43$ K) and *J* = 1–0 ($E_u \approx 4$ K) lines mean that we do not necessarily expect both lines to trace the same disk vertical layers. The WaOph 6 measurement, although using the same transition as in PDS 70, suffered from high uncertainties with a possible $z/r$ ranging from 0.1 to 0.5 due to the coarse beam size ($\approx 0.3''$) and lower SNR of the observations. This range of heights is consistent with what we measure in PDS 70 but does not permit a detailed comparison of the emitting surfaces. The most similar previous measurement is in the transition disk GM Aur, where Huang et al. (2020) used high-angular-resolution observations to fit a single power-law profile to the HCO<sup>+</sup> *J* = 3–2 emission surface and found an approximate height of $z/r \gtrsim 0.1$. Although not identical, the *J* = 3–2 ($E_u \approx 26$ K) line has more comparable excitation properties to that of the *J* = 4–3 line.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Comparison of the characteristic ¹²CO (top row) and ¹³CO (bottom row) emission heights of the PDS 70 disk (in pink) vs other disks in the literature (Law et al. 2023, and see references therein). Transition disks are shown as symbols with a hollow center. All stellar masses are dynamical, *z/r* values are computed homogeneously via directly extracted emission surfaces, and disk sizes represent the radius enclosing 90% of the total flux. The dynamical stellar mass of PDS 70 is taken from Keppler et al. (2019). Several scaling relations, i.e., line-emission heights scale with gas-pressure scale height and a previously identified positive ¹²CO height trend with increasing disk size, are labeled (see Law et al. 2022, 2023, for more details).</figcaption></figure>

Moreover, both PDS 70 and GM Aur are T Tauri stars of similar spectral types and host transition disks with central cavities of comparable size with the main difference being that GM Aur is at least twice as large in overall disk size in both its millimeter dust and ¹²CO line emission. For consistency, we refit the HCO⁺ *J* = 3–2 data in the GM Aur disk from Huang et al. (2020) using an exponentially tapered power-law with the same procedure as in Section 2.3.

Figure 10 shows the resultant emission surface, with best-fit parameters of $z_0 = 0.22_{-0.02}^{+0.02}$, $\phi = 0.99_{-0.17}^{+0.18}$, $r_{\text{taper}} = 1.93_{-0.10}^{+0.09}$, $\psi = 4.01_{-0.73}^{+0.61}$, and $r_{\text{cavity}} = 0.08_{-0.06}^{+0.08}$. Our emission surface is higher than the one derived in Huang et al. (2020) at all radii within 250 au, which is due to our surface-extraction method being able to more accurately track individual pixels in the channels.

HCO⁺ is emitting from similar disk vertical regions ($z/r \approx 0.2$) in both the GM Aur and PDS 70 disks, and thus, similar to the CO isotopologues, PDS 70 does not appear to show any significant differences, which could be ascribed to the influence of PDS 70b or 70c. As noted in Section 3.1 and shown in

Figure 10, the relative heights of the ¹³CO *J* = 2–1 and HCO⁺ *J* = 4–3 surfaces cross over beyond ~100 au. We find the same trend in the GM Aur disk, with the HCO⁺ surface located at higher altitudes than ¹³CO in the inner ~200 au, while at large radii, the ¹³CO is coming from a higher disk layer than HCO⁺. The cross-over point of these surfaces occurs directly exterior to the bulk of the millimeter continuum emission in both disks, which suggests that this effect may be due to a change in illumination at the continuum edge. When combined with the lower densities expected at larger disk radii, this, in turn, would likely alter how far ionizing photons can penetrate in the disk and could thus naturally explain the lower heights of the HCO⁺ surfaces in the outer disk. Although intriguing, with only two sources, it is difficult to make any additional conclusions about the vertical distribution of HCO⁺ in protoplanetary disks more generally. We stress that high-angular-resolution HCO⁺ observations in additional disks, especially those with different properties than PDS 70 and GM Aur, are urgently needed to set stringent constraints on disk ionization structure.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** (Top) Comparison of HCO$^+$ surfaces in the PDS 70 (left) and GM Aur (right) disks. The GM Aur data were obtained from Huang et al. (2020) and we extracted the emission surface as in Section 2.3. The single power-law fit from Huang et al. (2020) is also shown for reference. Otherwise, as in Figure 2. (Bottom) Zoom-in on the HCO$^+$ (orange) and $^{13}$CO (gray) emission surfaces. The millimeter continuum disk size is marked by a vertical line. The $^{13}$CO surface in the GM Aur disk is from Law et al. (2021a).</figcaption></figure>

## 4.2. Comparison to the LkCa 15 Transition Disk

Here, we compare the PDS 70 disk to that of LkCa 15, which is one of the few transition disks whose vertical emission structure has been mapped in detail with multiple molecular lines (Leemker et al. 2022; Law et al. 2023). LkCa 15 has a similar stellar spectral type (K5) and age (~5 Myr) (Donati et al. 2019) and while the LkCa 15 disk is much larger than that of PDS 70 in both millimeter dust continuum and CO isotopologue emission by a factor of ≈2–3, both disks have nearly identical central cavity sizes (e.g., Piétu et al. 2006; Facchini et al. 2020). Moreover, both disks have comparable gas temperatures within 200 au and similar $^{12}$CO and $^{13}$CO emission layer $z/r$ heights (Leemker et al. 2022; Law et al. 2023). Thus, they are ideal systems to compare how line-emitting heights change across their central cavities.

Figure 11 shows line-emission surfaces in each disk zoomed into the location of the central cavity, as indicated in the corresponding azimuthally averaged millimeter continuum radial profile. The PDS 70 disk shows an increase in emitting heights in $^{12}$CO $J=3-2$ associated with the edge of the cavity wall, while no such effect is seen in LkCa 15. While in Figure 11, we show the $J=2-1$ lines of $^{12}$CO and $^{13}$CO for the LkCa 15 disk, emitting height measurements also exist for $^{12}$CO $J=3-2$ and $^{13}$CO $J=6-5$ (Leemker et al. 2022), but neither of these surfaces show a similar height increase as in PDS 70.

However, the absence of this cavity-wall effect in the LkCa 15 disk is likely an observational limitation, rather than indicating any fundamental differences in inner-disk structure. For instance,

in Figure 11, the $^{12}$CO $J=2-1$ beam is ≈0.′′3, which at the distance of LkCa 15 (157 pc; Gaia Collaboration et al. 2021), corresponds to a physical scale of ≈50 au, which is approximately the size of the central cavity. In comparison, our $^{12}$CO $J=3-2$ observations of PDS 70 trace physical scales of ≈11–15 au, a factor of several times better. We note that while the $^{13}$CO $J=6-5$ surface in LkCa 15 from Leemker et al. (2022) was derived from observations of comparable angular resolution, emitting heights interior to ≈40 au could not be extracted. Moreover, $^{13}$CO originates at lower elevations than that of $^{12}$CO, which means that this height increase at the cavity wall would be intrinsically more difficult to observe using $^{13}$CO.

Emission from $^{12}$CO $J=2-1$, 3–2 and $^{13}$CO $J=6-5$ have been detected in the central cavity of LkCa 15 (Jin et al. 2019; Leemker et al. 2022; Law et al. 2023), which suggests that emission height constraints here would be possible with sufficiently sensitive and high-angular-resolution data. Leemker et al. (2022) found that the gas column density in the dust cavity (≈45 au) of LkCa 15 drops by a factor of >2 compared to the outer disk, with an additional order-of-magnitude decrease inside 10 au. In comparison, the PDS 70 cavity shows a higher gas depletion (≥10×) over its full cavity (Portilla-Revelo et al. 2023). Thus, while we still expect flatter emission surfaces in the LkCa 15 cavity, this cavity-wall feature may be less pronounced in the LkCa 15 disk due to its more modest gas density decrease. Nonetheless, observations from the ongoing exoALMA Large Program will probe similar physical resolutions as in the PDS 70 disk in several transition disks,
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Zoomed-in emission surfaces of transition disks around PDS 70 and LkCa 15 (top row) compared to azimuthally averaged millimeter continuum radial profiles (bottom row). Lines of constant z/r from 0.1 to 0.5 are shown in gray. Emission surfaces and continuum profiles in LkCa 15 are from Leemker et al. (2022), Law et al. (2023), and Facchini et al. (2020), respectively. Only the PDS 70 disk shows an increase in line-emitting heights at the location of the cavity wall.</figcaption></figure>

including LkCa 15, and should be able to confirm if such
changes in emitting heights at cavity walls are commonplace or
if the PDS 70 disk structure is unique.

4.3. CO Snowline Location

The elemental carbon-to-oxygen (C/O) ratio is of vital
importance to understand the accretion history of protoplanets
and connect disk chemistry with planetary organic compositions
and atmospheres (e.g., Madhusudhan 2019; Öberg & Ber-
gin 2021). While chemical evolution may alter the C/O ratio
(e.g., Eistrup et al. 2018; Cridland et al. 2019; Krijt et al. 2020),
the condensation of major volatile species across the disk has the
most dramatic influence on the elemental ratio of gas-phase
material (e.g., Öberg et al. 2011). These condensation fronts,
often referred to as snowlines, occur at specific temperatures and
thus, with detailed knowledge of the disk temperature structure,
it is possible to estimate the location of these snowlines directly.
Here, we focus on the CO snowline, which is readily
observationally accessible at spatial scales probed by ALMA.

Figure 12 shows the midplane gas temperature as inferred
from the 2D empirical fit (Figure 8). It is important to note that
the midplane temperature from this fit is an extrapolation from
higher disk heights since we lack lines that directly trace the

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Midplane gas temperature as a function of radius in the PDS 70 disk as inferred from the 2D empirical temperature fit. The blue horizontal region indicates a range of CO freeze-out temperatures (18–22 K), which corresponds to a midplane snowline radius of ≈56–85 au. The peak brightness temperature profile of C<sup>18</sup>O J = 2–1 extracted in a narrow ±5° wedge along the disk major axis is shown in orange.</figcaption></figure>
disk midplane interior to a radius of 100 au. As a result, we also show the C¹⁸O J = 2–1 peak brightness temperature profile in Figure 12, which we expect to be emitting at, or near, the disk midplane (Section 3.1). Beyond the central cavity, the C¹⁸O brightness temperature is generally consistent with our estimated midplane temperature and is ≈5–10 K colder at larger radii (≈100 au). As C¹⁸O is not expected to be fully optically thick, especially at large radii, it only provides a lower limit on the gas temperature, which is consistent with our warmer midplane temperature estimates.

If we adopt a reasonable range of possible CO freeze-out temperatures from 18–22 K (e.g., Collings et al. 2004; Öberg et al. 2005; Martín-Doménech et al. 2014; Facchini et al. 2017), we estimate an approximate snowline radius of ≈56–85 au. Despite this range, we can confidently conclude that both PDS 70b and 70c are located interior of the CO snowline. Thus, the nearby gas from which both planets are accreting likely has a super-stellar<fnref n="11" /> C/O ratio. Depending on the exact location of the CO₂ snowline, this gas likely has a C/O of ≲1.0 (Öberg et al. 2011), with the exact value depending on the details of radial transport and chemical processing (e.g., Krijt et al. 2018, 2020). From our midplane temperature estimate, we would estimate a CO₂ snowline location of ≈21–26 au (assuming freeze-out temperatures of 42–52 K), which is consistent with C/O ≈1, with the caveat that these are extrapolations to radii within the inner cavity where we lack direct temperature measurements.

Overall, this is consistent with Facchini et al. (2021), who found that the PDS 70 molecular layer hosts a high C/O ratio inferred from bright C₂H line fluxes and their lower limit on the CS/SO column density ratio (e.g., Bergin et al. 2016; Semenov et al. 2018; Miotello et al. 2019; Le Gal et al. 2021). The fact that both planets in the PDS 70 disk are likely accreting gas with a high C/O ratio provides important context for interpreting future atmospheric planetary measurements from, e.g., JWST, and for properly modeling their formation histories within the PDS 70 disk (e.g., Cridland et al. 2023).

Further observations of temperature tracers, such as additional *J* lines of CO isotopologues (e.g., Leemker et al. 2022), especially those closer to the midplane and in the inner disk would allow for a better constrained temperature model and lead to a more accurate prediction of the snowline location. Complementary to this, observations of N₂H⁺, which has been shown to be an accurate observational tracer of the CO snowline (Qi et al. 2019), would provide an independent estimate of the CO snowline location. If combined with our temperature-based estimate, this would provide a unique opportunity to measure the CO binding energy in an astrophysical context and provide important and independent context to laboratory-based estimates.

## 5. Conclusions

We present observations of a suite of CO isotopologue and HCO⁺ lines toward the PDS 70 disk at high angular resolution (≈0.′1). We extracted line-emission surfaces and conclude the following:

1. ¹²CO emission surfaces originate in elevated disk regions ($z/r \approx 0.3$), while the less abundant ¹³CO-emitting heights trace deeper disk layers ($z/r \approx 0.1$). We also

derived emitting heights for the HCO⁺ *J* = 4–3 line, which arises from $z/r \approx 0.2$ in the inner 100 au and shows an extended, flat ($z/r \lesssim 0.1$) component at larger radii.

2. In the ¹²CO *J* = 3–2 line, we clearly resolve the vertical dip and steep rise in emitting heights at ≈15–45 au due to the cavity wall in the PDS 70 disk. This is the first transition disk where this effect has been directly seen in emitting heights.

3. The emitting heights of the CO isotopologues in the outer disk of PDS 70 appear typical for its stellar mass and disk size and do not appear to be substantially altered by the presence of its two embedded planets. Overall, this suggests that the outer disk is largely decoupled from the planet formation occurring in the inner few 10 s of astronomical units.

4. By combining CO isotopologue and HCO⁺ lines, we derive an empirical 2D temperature structure for the PDS 70 disk. Using this, we estimate an approximate CO midplane snowline radius of ≈56–85 au, which implies that both PDS 70b and 70c are located well interior to the CO snowline and suggests they are both accreting from gas with a C/O ratio of ≈1.0.

## Acknowledgments

The authors thank the anonymous referee for valuable comments that improved both the content and presentation of this work. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.00888.S, 2017.A.00006.S, 2017.1.01151.S, 2018.A.00030.S, 2018.1.01230.S, and 2019.1.01619.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Support for C.J.L. was provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51535.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. S. F. is funded by the European Union under the European Union's Horizon Europe Research & Innovation Programme 101076613 (UNVEIL). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. 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).

Facility: ALMA

Software: Astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018; Astropy Collaboration et al. 2022), bettermoments (Teague & Foreman-Mackey 2018), CASA (McMullin et al. 2007; CASA Team et al. 2022), cmasher (van der Velden 2020), disksurf (Teague et al. 2021),

<footnote n="11">¹¹ Based on its photospheric abundances of carbon and oxygen, PDS 70 has a C/O ratio of 0.44 ± 0.19 (Suárez-Andrés et al. 2018; Cridland et al. 2023).</footnote>
GoFish (Teague 2019), keplerian_mask (Teague 2020), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), SciPy (Virtanen et al. 2020).

## Appendix A
Line-emission Channel Maps

A complete gallery of channel maps for all CO isotopologue lines and HCO$^+$ $J = 4-3$ is shown in Figure Set 13, which is available in the electronic edition of the journal.

## Appendix B
Radial Temperatures versus Peak Brightness Temperatures

We extracted radial profiles from the peak-intensity maps (Figure 3) by adopting the inferred emission surfaces and following the same procedure as in Section 3.2. In Figure 14, we compare them to the radial temperatures inferred directly from the emission surfaces. Overall, we find excellent agreement between the two curves. The most notable differences occur in the inner 100 au of the $^{13}$CO $J = 2-1$

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Channel maps of the $^{12}$CO $J = 2-1$ emission of the PDS 70 disk. The synthesized beam is shown in the lower left corner of each panel and the LSRK velocity in km s$^{-1}$ is printed in the upper right.</figcaption></figure>

(The complete figure set (five images) is available.)
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Peak brightness (blue curve) vs. radial temperature profiles (points) computed along the emission surfaces. The azimuthally averaged peak-intensity profiles are generated by deprojecting along the derived emission surfaces.</figcaption></figure>

<table><caption>**Table 7**
Properties of Radial Temperature Dips</caption><thead><tr><td>Line</td><td>Feature</td><td>r<sub>0</sub> (au)</td><td>Width (au)</td><td>ΔT (K)</td><td>Depth</td></tr></thead><tbody><tr><td><sup>12</sup>CO J = 2–1</td><td>D115</td><td>115 ± 0.3</td><td>22 ± 11</td><td>2.8 ± 0.9</td><td>0.09 ± 0.003</td></tr><tr><td><sup>12</sup>CO J = 3–2</td><td>D50<sup>a</sup></td><td>50 ± 2</td><td>13 ± 0.0</td><td>2.3 ± 0.3</td><td>0.07 ± 0.002</td></tr><tr><td><sup>12</sup>CO J = 3–2</td><td>D127</td><td>127 ± 1</td><td>33 ± 12</td><td>2.1 ± 0.7</td><td>0.08 ± 0.01</td></tr></tbody></table>

**Note.**

<sup>a</sup> Computed from the peak-intensity radial profile.

profile, where the peak intensities have temperatures that are at most ≲5 K larger than the surface-extracted ones. We also note that the dips identified at ≈ 120 au in both <sup>12</sup>CO J = 2–1 and J = 3–2 from Figure 6 are also present in the peak intensities.

We also find a small dip in the peak-intensity profiles of <sup>12</sup>CO J = 3–2 at ~50 au that is not apparent in the radial brightness temperature profiles. This dip was identified previously in integrated intensity by Keppler et al. (2019) and was attributed to continuum absorption of the back side of the disk, since it appears coincident with the millimeter continuum gap. However, given its appearance in the peak-intensity profile, which should be less affected by continuum absorption, this suggests that it is a real feature and may hint at

a true gas and dust gap. While we do not see this feature in any of the other CO lines, this is perhaps not surprising given the small width of this feature and the larger beam sizes of the other CO lines. Additional observations at high angular resolution of other CO isotopologues are necessary to confirm the nature of this feature.

As noted in Section 3.4.1, we also determined the properties of the temperature dips seen in Figures 6 and 14 using the same method as for the vertical substructures, i.e., via removing a local baseline (Law et al. 2021a). Table 7 shows the properties of each of these features, including the 50 au dip in <sup>12</sup>CO J = 3–2 identified in the peak-intensity radial profiles.
ORCID iDs

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

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

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

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

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

Andrea Isella https://orcid.org/0000-0001-8061-2207

Inga Kamp https://orcid.org/0000-0001-7455-5349

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

Bayron Portilla-Revelo https://orcid.org/0000-0002-6278-9006

References

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

Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

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

Bae, J., Zhu, Z., Baruteau, C., et al. 2019, ApJL, 884, L41

Benisty, M., Bae, J., Facchini, S., et al. 2021, ApJL, 916, L2

Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101

Boehler, Y., Weaver, E., Isella, A., et al. 2017, ApJ, 840, 60

Booth, A. S., Walsh, C., & Ilee, J. D. 2019, A&A, 629, A75

Calahan, J. K., Bergin, E. A., Zhang, K., et al. 2021, ApJS, 257, 17

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

CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501

Christiaens, V., Casassus, S., Absil, O., et al. 2019, MNRAS, 486, 5819

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

Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133

Cridland, A. J., Eistrup, C., & van Dishoeck, E. F. 2019, A&A, 627, A127

Cridland, A. J., Facchini, S., van Dishoeck, E. F., & Benisty, M. 2023, A&A, 674, A211

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

Dartois, E., Dutrey, A., & Guilloteau, S. 2003, A&A, 399, 773

Donati, J. F., Bouvier, J., Alencar, S. H., et al. 2019, MNRAS, 483, L1

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

Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 613, A14

Facchini, S., Benisty, M., Bae, J., et al. 2020, A&A, 639, A121

Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 605, A16

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

Facchini, S., Teague, R., Bae, J., et al. 2021, AJ, 162, 99

Favre, C., Fedele, D., Maud, L., et al. 2019, ApJ, 871, 107

Follette, K. B., Close, L. M., Males, J. R., et al. 2023, AJ, 165, 225

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

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1

Galloway-Sprietsma, M., Bae, J., Teague, R., et al. 2023, ApJ, 950, 147

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

Hammond, I., Christiaens, V., Price, D. J., et al. 2023, MNRAS, 522, L51

Hashimoto, J., Dong, R., Kudo, T., et al. 2012, ApJL, 758, L19

Hashimoto, J., Tsukagoshi, T., Brown, J. M., et al. 2015, ApJ, 799, 43

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

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2020, ApJ, 891, 48

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

Ilee, J. D., Walsh, C., Booth, A. S., et al. 2021, ApJS, 257, 9

Isella, A., Benisty, M., Teague, R., et al. 2019, ApJL, 879, L25

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, A&A, 674, A113

Jin, S., Isella, A., Huang, P., et al. 2019, ApJ, 881, 108

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

Juillard, S., Christiaens, V., & Absil, O. 2022, A&A, 668, A125

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

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

Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134

Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78

Kuo, I. H. G., Yen, H.-W., Gu, P.-G., & Chang, T.-E. 2022, ApJ, 938, 50

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

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

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

Law, C. J., Teague, R., Öberg, K. I., et al. 2023, ApJ, 948, 60

Le Gal, R., Öberg, K. I., Teague, R., et al. 2021, ApJS, 257, 12

Leemker, M., Booth, A. S., van Dishoeck, E. F., et al. 2022, A&A, 663, A23

Long, F., Andrews, S. M., Rosotti, G., et al. 2022, ApJ, 931, 6

Long, Z. C., Akiyama, E., Sitko, M., et al. 2018, ApJ, 858, 112

Madhusudhan, N. 2019, ARA&A, 57, 617

Martín-Doménech, R., Muñoz Caro, G. M., Bueno, J., & Goesmann, F. 2014, A&A, 564, A8

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

Mesa, D., Keppler, M., Cantalloube, F., et al. 2019, A&A, 632, A25

Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A&A, 631, A69

Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74

Montesinos, M., Cuello, N., Olofsson, J., et al. 2021, ApJ, 910, 31

Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, 832, 41

Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2

Öberg, K. I., & Bergin, E. A. 2021, PhR, 893, 1

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

Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJL, 743, L16

Öberg, K. I., van Broekhuizen, F., Fraser, H. J., et al. 2005, ApJL, 621, L33

Paneque-Carreño, T., Miotello, A., van Dishoeck, E. F., et al. 2022, A&A, 666, A168

Paneque-Carreño, T., Miotello, A., van Dishoeck, E. F., et al. 2023, A&A, 669, A126

Pecaut, M. J., & Mamajek, E. E. 2016, MNRAS, 461, 794

Piétu, V., Dutrey, A., Guilloteau, S., Chapillon, E., & Pety, J. 2006, A&A, 460, L43

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

Portilla-Revelo, B., Kamp, I., Facchini, S., et al. 2023, A&A, 677, A76

Portilla-Revelo, B., Kamp, I., Rab, C., et al. 2022, A&A, 658, A89

Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123

Qi, C., Öberg, K. I., Espaillat, C. C., et al. 2019, ApJ, 882, 160

Rich, E. A., Teague, R., Monnier, J. D., et al. 2021, ApJ, 913, 138

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

Rosotti, G. P., Ilee, J. D., Facchini, S., et al. 2021, MNRAS, 501, 3427

Sanchis, E., Testi, L., Natta, A., et al. 2021, A&A, 649, A19

Semenov, D., Favre, C., Fedele, D., et al. 2018, A&A, 617, A28

Stapper, L. M., Hogerheijde, M. R., van Dishoeck, E. F., & Paneque-Carreño, T. 2023, A&A, 669, A158

Stolker, T., Marleau, G. D., Cugno, G., et al. 2020, A&A, 644, A13

Suárez-Andrés, L., Israelian, G., González Hernández, J. I., et al. 2018, A&A, 614, A84

Teague, R. 2019, JOSS, 4, 1632

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

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

Teague, R., & Foreman-Mackey, D. 2018, 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

Tripathi, A., Andrews, S. M., Birnstiel, T., & Wilner, D. J. 2017, ApJ, 845, 44 van der Velden, E. 2020, JOSS, 5, 2004

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

Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJL, 863, L8

Wagner, K., Stone, J., Skemer, A., et al. 2023, NatAs, 7, 1208

Wang, J. J., Vigan, A., Lacour, S., et al. 2021, AJ, 161, 148

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

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