# Molecules with ALMA at Planet-forming Scales (MAPS). IV. Emission Surfaces and Vertical Distribution of Molecules

Charles J. Law, Richard Teague, Ryan A. Loomis, Jaehan Bae, Karin I. Öberg, Ian Czekala, Sean M. Andrews, Yuri Aikawa, Felipe Alarcón, Edwin A. Bergin, et al.

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

► To cite this version:

Charles J. Law, Richard Teague, Ryan A. Loomis, Jaehan Bae, Karin I. Öberg, et al. Molecules with ALMA at Planet-forming Scales (MAPS). IV. Emission Surfaces and Vertical Distribution of Molecules. The Astrophysical Journal Supplement Series, 2021, 257, 10.3847/1538-4365/ac1439. insu-03672377

HAL Id: insu-03672377

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

Submitted on 19 Sep 2025

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

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

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

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

Molecules with ALMA at Planet-forming Scales (MAPS). IV. Emission Surfaces and Vertical Distribution of Molecules

Charles J. Law¹, Richard Teague¹, Ryan A. Loomis², Jaehan Bae³,4,25, Karin I. Öberg¹, Ian Czekala⁵,⁶,⁷,⁸,⁹,<fnref n="25" />,  
Sean M. Andrews¹, Yuri Aikawa¹⁰, Felipe Alarcón¹¹, Edwin A. Bergin¹¹, Jennifer B. Bergner¹²,<fnref n="25" />,  
Alice S. Booth¹³,¹⁴, Arthur D. Bosman¹¹, Jenny K. Calahan¹¹, Gianni Cataldi¹⁰,¹⁵, L. Ilsedore Cleeves¹⁶,  
Kenji Furuya¹⁵, Viviana V. Guzmán¹⁷, Jane Huang¹,¹¹,<fnref n="25" />, John D. Ilee¹⁴, Romane Le Gal¹,¹⁸,¹⁹,²⁰, Yao Liu²¹,  
Feng Long¹, François Ménard¹⁹, Hideko Nomura¹⁵, Laura M. Pérez²², Chunhua Qi¹, Kamber R. Schwarz²³,<fnref n="25" />,  
Daniela Soto¹⁷, Takashi Tsukagoshi¹⁵, Yoshihide Yamato¹⁰, Merel L. R. van 't Hoff¹¹, Catherine Walsh¹⁴,  
David J. Wilner¹, and Ke Zhang¹¹,²⁴,<fnref n="26" />

¹ Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA; charles.law@cfa.harvard.edu

² National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA

³ Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA

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

⁵ Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

⁶ Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

⁷ Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

⁸ Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA

⁹ Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA

¹⁰ Department of Astronomy, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan

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

¹² University of Chicago, Department of the Geophysical Sciences, Chicago, IL 60637, USA

¹³ Leiden Observatory, Leiden University, 2300 RA Leiden, The Netherlands

¹⁴ School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK

¹⁵ National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

¹⁶ Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA

¹⁷ Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile

¹⁸ IRAP, Université de Toulouse, CNRS, CNES, UT3, F-31000 Toulouse, France

¹⁹ Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

²⁰ IRAM, 300 rue de la piscine, F-38406 Saint-Martin d'Hères, France

²¹ Purple Mountain Observatory & Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, Nanjing 210023, People's Republic of China

²² Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile

²³ Lunar and Planetary Laboratory, University of Arizona, 1629 E. University Blvd., Tucson, AZ 85721, USA

²⁴ Department of Astronomy, University of Wisconsin–Madison, 475 N. Charter St., Madison, WI 53706, USA

Received 2021 February 22; revised 2021 July 2; accepted 2021 July 4; published 2021 November 3

Abstract

The Molecules with ALMA at Planet-forming Scales (MAPS) Large Program provides a unique opportunity to study the vertical distribution of gas, chemistry, and temperature in the protoplanetary disks around IM Lup, GM Aur, AS 209, HD 163296, and MWC 480. By using the asymmetry of molecular line emission relative to the disk major axis, we infer the emission height (z) above the midplane as a function of radius (r). Using this method, we measure emitting surfaces for a suite of CO isotopologues, HCN, and C2H. We find that ¹²CO emission traces the most elevated regions with z/r > 0.3, while emission from the less abundant ¹³CO and C¹⁸O probes deeper into the disk at altitudes of z/r ~ 0.2. C2H and HCN have lower opacities and signal-to-noise ratios, making surface fitting more difficult, and could only be reliably constrained in AS 209, HD 163296, and MWC 480, with z/r ~ 0.1, i.e., relatively close to the planet-forming midplanes. We determine peak brightness temperatures of the optically thick CO isotopologues and use these to trace 2D disk temperature structures. Several CO temperature profiles and emission surfaces show dips in temperature or vertical height, some of which are associated with gaps and rings in line and/or continuum emission. These substructures may be due to local changes in CO column density, gas surface density, or gas temperatures, and detailed thermochemical models are necessary to better constrain their origins and relate the chemical compositions of elevated disk layers with those of planet-forming material in disk midplanes. This paper is part of the MAPS special issue of the Astrophysical Journal Supplement.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Interstellar molecules (849); Astrochemistry (75); High angular resolution (2167); Circumstellar disks (235); Planet formation (1241)*

*Supporting material:* figure set

**1. Introduction**

Protoplanetary disks are highly stratified in their physical
and chemical properties (Williams & Cieza 2011) with flared

emitting surfaces set by the balance of hydrostatic equilibrium,
as first recognized in their spectral energy distributions
(Kenyon & Hartmann 1987). In particular, vertical gradients
in gas temperature, density, radiation, and ionization result in a
rich chemical structure over the height of the disk (e.g., van
Zadelhoff et al. 2001; Woitke et al. 2009; Walsh et al. 2010;

<footnote n="25">²⁵ NASA Hubble Fellowship Program Sagan Fellow.</footnote>

<footnote n="26">²⁶ NASA Hubble Fellow.</footnote>
Fogel et al. 2011). The efficiency of vertical mixing (Ilgner et al. 2004; Semenov & Wiebe 2011) and the presence of meridional flows driven by embedded planets (Morbidelli et al. 2014; Dong et al. 2019; Teague et al. 2019) also influence the vertical distribution of molecular material in disks. Vertical chemical structures have been seen in observations of highly inclined disks, which allow the emission distribution to be mapped directly (Dutrey et al. 2017; Teague et al. 2020; Podio et al. 2020; Ruíz-Rodríguez et al. 2021). In more moderately inclined disks, the excitation temperatures of different species and molecular isotopologues have instead been used to infer the properties of the vertical gas distribution (Dartois et al. 2003; Piétu et al. 2007; Öberg et al. 2021; Cleeves et al. 2021). A detailed understanding of this vertical structure is required to interpret these observations and to assess how well connected the molecular gas abundances derived from line observations are to those of the planet-forming material in disk midplanes.

The high spatial and spectral resolutions offered by ALMA allow for the direct measurement of the height at which molecular emission arises for mid-inclination disks. With sufficient angular resolution and surface brightness sensitivity, it is possible to spatially resolve emission arising from elevated regions above and below the midplane (e.g., de Gregorio-Monsalvo et al. 2013; Rosenfeld et al. 2013; Isella et al. 2018; Huang et al. 2020). In cases such as these, the emission surface of molecular lines can be directly extracted with a technique similar to that used in near-IR (NIR) observations to infer scattering surfaces (e.g., Monnier et al. 2017; Avenhaus et al. 2018). Such a method was first presented by Pinte et al. (2018), who used it to map the CO, ¹³CO, and C¹⁸O 2–1 emission surfaces in IM Lup. As the less abundant isotopologues probe deeper in the disk (i.e., closer to the midplane), this also allows for an empirical derivation of the two-dimensional gas temperature structure, an essential input for models and simulations of planet formation. A similar approach has now been employed to map CO isotopologue surfaces in a handful of disks (e.g., Teague et al. 2019; Paneque-Carreño et al. 2021; Rich et al. 2021).

As part of the Molecules with ALMA at Planet-forming Scales (MAPS) Large Program, five protoplanetary disks were observed in several molecular lines expected to emit strongly at different vertical locations. In this paper, we provide a framework for generalizing the method presented in Pinte et al. (2018) to a larger sample of disks. We use this framework to characterize line emission heights, gas temperatures, and disk vertical substructures. The layout of the paper is as follows: We present a brief overview of the observations in Section 2. In Section 3, we describe how emission surfaces were derived and fit with analytical functions. We calculate the radial and vertical temperature profiles and compare the observed vertical structures with previous millimeter and NIR observations in Section 4. In Section 5, we discuss possible origins of disk vertical substructures. We summarize our findings in Section 6. All publicly available data products are listed in Section 7.

## 2. Observations

The observations used in this study were obtained as part of the MAPS ALMA Large Program (2018.1.01055.L), which targeted the protoplanetary disks around IM Lup, GM Aur, AS 209, HD 163296, and MWC 480. An overview of the survey, including observational setup, calibration, and rationale, is provided in Öberg et al. (2021), while the imaging process is described in detail in Czekala et al. (2021). The analysis in this

work is based primarily on the ALMA Band 6 images generated with a robust parameter of 0.5, rather than the fiducial images presented in the overview paper. We opted to use these images to leverage their higher spatial resolutions (10%–40% smaller total beam area) relative to the fiducial images with a 0.′′15 circularized beam. The major and minor axes of the synthesized beam were in the range of 0.′′13–0.′′17 and 0.′′10–0.′′13, respectively. At the distances of the MAPS disks, these correspond to physical scales of ~14–27 au and 10–20 au, respectively. For extracting brightness temperatures, we also made use of the corresponding non-continuum-subtracted image cubes, which were imaged in the same way as the line-only data.

This work is based on the following five species: CO, ¹³CO, C¹⁸O, HCN, and C₂H. The primary focus for the CO isotopologues is on the Band 6 transitions, i.e., $J = 2–1$, as they possess the highest spatial resolutions. For HCN and C₂H, we only considered the brightest hyperfine components of the Band 6 transitions, i.e., C₂H $N = 3–2$, $J = \frac{7}{2} \frac{5}{2}$, $F = 4–3$ and HCN $J = 3–2$, $F = 3–2$. For simplicity, we refer to these lines as HCN 3–2 and C₂H 3–2, respectively. This set of molecules and lines was selected for this analysis because they are consistently the brightest in the MAPS sample and possessed radially extended emission that allowed for the determination of robust emission surfaces (Law et al. 2021).

Since the ¹³CO and C¹⁸O isotopologues were observed in both Bands 3 and 6, we also aimed to assess the influence of excitation on the derived surfaces, namely, if different transitions from the same molecule emit at different disk heights. To do so, we used the tapered (0.′′30) images (see Section 6.2, Czekala et al. 2021), which allowed us to match the spatial resolutions of the 2–1 (Band 6) and 1–0 (Band 3) lines for both species. This comparison, however, did not yield any conclusive results (see Appendix C).

## 3. Emission Surfaces

In the following subsections, we present an outline of how we derived emitting layers starting from the line image cubes to the final data products. We then describe how we fit an analytical function to each of the surfaces.

### 3.1. Deriving Emitting Layers

The emission surfaces derived in this work represent the mean height of the emission surface for each molecular tracer, or put simply, where the bulk of the emission arises from in each line and disk. We extracted these emission heights from the image cubes by using the asymmetry of the emission relative to the major axis of the disk. By assuming that disks are azimuthally symmetric and that the gas is on circular orbits, this allows us to infer the height of the emission above the disk midplane. A prerequisite for this approach is the ability to spatially resolve the front and back of the disk in multiple channels, as illustrated for CO 2–1, ¹³CO 2–1, and C¹⁸O 2–1 in Figure 1. For each line and disk, we determined whether both disk sides were sufficiently spatially resolved via visual inspection and then confirmed that the predicted isovelocity contours matched the spatial distribution of line emission in the channel maps (see Appendix B). We refer readers to Pinte et al. (2018) for additional details about this method.

We derived emission surfaces using the *disksurf*<fnref n="27" /> Python package, which implements this method while providing

<footnote n="27">²⁷ https://github.com/richteague/disksurf</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Representative channels for the MAPS sample, ordered from left to right by increasing stellar mass (see Table 1 in Öberg et al. 2021), for CO 2–1, $^{13}$CO 2–1, and C$^{18}$O 2–1. Both sides of the disk are visible and can be used to fit for the emitting layer. The C$^{18}$O 2–1 images in IM Lup and GM Aur have been tapered to 0.′′30, as described in Section 4.1. Kinematic local standard of rest (LSRK) velocities are shown in the upper right corner, and the plus sign indicates the disk center. Axes are angular offsets from the disk center with 1′′ tick marks. An arcsinh color stretch has been applied to highlight faint outer structures.</figcaption></figure>

additional functionality to filter the data to extract more precise emission surfaces. This series of filtering steps is described in detail in Appendix A. We used the `get_emission_` surface function to extract the deprojected radius *r*, emission height *z*, surface brightness *I*<sub>ν</sub>, and channel velocity *v* for each pixel associated with the emitting surface. These surfaces represent individual measurements, i.e., pixels, from the line image cubes.

We then use two different methods to further reduce scatter in the individual emission surface measurements and help better identify substructure. First, we radially bin the individual measurements using bin sizes of 1/4 of the FWHM of the beam major axis, i.e., the same as the radial intensity profiles in Law et al. (2021). The uncertainty is given by the standard deviation in each bin. We note that Pinte et al. (2018) included the uncertainty in the disk inclination in these uncertainties. We opted not to do this, as the disk inclination is a systematic uncertainty and results in a scaling of the vertical height axis, not a relative uncertainty between radial bins. Besides binning, we also calculated the moving average and standard deviation of the individual surface measurements. As the spacing between radial points is not uniform, we used a window size with a minimum size of one-quarter of the FWHM of the beam major axis. This window was required to contain a fixed number of points, which means that the physical size that it represents changes with radius owing to the nonuniform radial sampling from the deprojection process, i.e., in the less dense, outer disk, the window expands in order to still encompass this fixed number of points. A summary of these different data products is shown in Figure 2. While it was found that binned and moving average surfaces showed the same trends, the

binned surface benefited from a uniform radial sampling, while the moving average retained a finer radial sampling, essential for identifying subtle perturbations associated with features in the dust continuum.

All three types of line emission surfaces—individual measurements, radially binned, and moving averages—are provided as value-added data products (VADPs) and are made available to the community through the MAPS project homepage (http://www.alma-maps.info). See Section 7 for further details. Throughout this work, we sometimes bin these data products further for visual clarity, but all quantitative analysis is done using the original binning of each type of emission surface.

## 3.2. Analytical Fits

To facilitate implementing these emission surfaces in models and for comparison with other observations, we fit an exponentially tapered power law to all CO emission surfaces. This parametric fit was chosen because it describes the flared surface in the inner disk (≤200 au) and captures the expected drop in the outer disk due to decreasing gas surface density, as seen in Figure 2. We adopt the following functional form:

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

where $z_0$, $\phi$, and $\psi$ should always be nonnegative. A value of $\phi > 1$ indicates that $z$ increases with radius, while $0 < \phi < 1$ tends toward a flat $z(r)$ profile. When $r_{\text{taper}} \gg 1''$, $z_0$ represents the $z/r$ value at 1$''$. Note that some previous works (e.g.,
<table><caption>**Table 1**
Parameters for Emission Surface Fits</caption><thead><tr><th rowspan="2">Source</th><th rowspan="2">Line</th><th rowspan="2">Velocity Range<br>(km s<sup>-1</sup>)</th><th colspan="6">Exponentially Tapered Power Law<sup>a</sup></th></tr><tr><th>r<sub>fit,in</sub> (arcsec)</th><th>r<sub>fit,out</sub> (arcsec)</th><th>z<sub>0</sub> (arcsec)</th><th>φ</th><th>r<sub>taper</sub> (arcsec)</th><th>ψ</th></tr></thead><tbody><tr><td rowspan="2">IM Lup</td><td>CO 2–1</td><td>[2.8, 6.4]</td><td>0.21</td><td>3.26</td><td>4.37</td><td>3.144</td><td>0.254</td><td>0.655</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>[2.6, 6.4]</td><td>0.61</td><td>2.02</td><td>0.159</td><td>2.599</td><td>1.928</td><td>4.993</td></tr><tr><td rowspan="3">GM Aur</td><td>CO 2–1</td><td>[3.1, 7.5]</td><td>0.08</td><td>3.36</td><td>0.385</td><td>1.066</td><td>3.767</td><td>4.988</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>[3.7, 7.1]</td><td>0.2</td><td>1.89</td><td>0.113</td><td>4.539</td><td>1.496</td><td>4.989</td></tr><tr><td>C<sup>18</sup>O 2–1</td><td>[-2.9, 13.9]</td><td>0.21</td><td>0.66</td><td>0.95</td><td>3.556</td><td>0.402</td><td>3.766</td></tr><tr><td rowspan="2">AS 209</td><td>CO 2–1</td><td>[2.9, 6.5]</td><td>0.07</td><td>1.98</td><td>0.219</td><td>1.292</td><td>1.786</td><td>4.854</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>[2.9, 6.5]</td><td>0.73</td><td>1.35</td><td>0.175</td><td>2.98</td><td>1.124</td><td>2.445</td></tr><tr><td rowspan="3">HD 163296</td><td>CO 2–1</td><td>[4.3, 13.5]</td><td>0.19</td><td>4.73</td><td>0.388</td><td>1.851</td><td>2.362</td><td>1.182</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>[3.5, 13.5]</td><td>0.31</td><td>3.42</td><td>0.121</td><td>1.503</td><td>3.158</td><td>4.996</td></tr><tr><td>C<sup>18</sup>O 2–1</td><td>[3.5, 8.1]</td><td>0.39</td><td>1.43</td><td>0.174</td><td>2.956</td><td>1.043</td><td>4.994</td></tr><tr><td rowspan="3">MWC 480</td><td>CO 2–1</td><td>[2.8, 7.4]</td><td>0.13</td><td>3.69</td><td>0.261</td><td>1.35</td><td>3.098</td><td>3.074</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>[2.8, 7.2]</td><td>0.05</td><td>2.26</td><td>1.248</td><td>2.165</td><td>0.215</td><td>0.683</td></tr><tr><td>C<sup>18</sup>O 2–1</td><td>[-8.2, 18.2]</td><td>0.13</td><td>1.46</td><td>0.065</td><td>1.37</td><td>0.961</td><td>4.834</td></tr></tbody></table>

**Note.**

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison of emission heights derived for CO 2–1 in HD 163296 for individual measurements, moving average, and radially binned surfaces (from top to bottom). For increased visual clarity, we occasionally bin surfaces by an extra factor of a few, as illustrated in the bottom panel. The FWHM of the major axis of the synthesized beam is shown in the lower right corner of each panel.</figcaption></figure>

Teague et al. 2019) instead used a double power-law profile to capture the drop in emission height at large radius. It was found that this tapered form, on average, provided a better fit to the data with less manual tuning required.

To ensure the robustness of these fits, we also restricted the radial range used for fitting to locations with high densities of (r, z) measurements. The radial ranges used in each fit are given by $r_{fit,in}$ and $r_{fit,out}$ in Table 1.

We used the affine-invariant Markov Chain Monte Carlo (MCMC) sampler (Goodman & Weare 2010) implemented in emcee (Foreman-Mackey et al. 2013) to estimate the posterior distributions of these fits. We used 64 walkers that take 1000 steps to burn in and an additional 500 steps to sample the posterior distribution function. We chose an MCMC fitting approach rather than a simple chi-squared minimization to better handle the degeneracies between fitted parameters. Individual pixels are not all necessarily independent, as they may originate within a single beam, which can lead to an underestimation of the true uncertainties on the extracted heights, i.e., on how well we can extract the mean z from a sample of columns. However, this will not necessarily affect the mean height itself. This is analogous to drawing random samples from a normal distribution, where, given a sufficiently large number of samples, the standard deviation of those samples provides a good estimate of the uncertainty on the mean of the distribution. Instead, if these samples are correlated and, e.g., for every second draw the sample is biased toward the one immediately preceding it, we will oversample the central region compared to the wings. In this case, the standard deviation of the total ensemble will underestimate the true standard deviation of the distribution but will not alter the underlying mean value.

The presence of this potential spatial correlation between pixels does not affect the analytical fits, which are instead dominated by overall radial trends rather than the vertical scatter in height measurements. Fits were also performed using the individual measurements, binned, and moving average surfaces, and we confirmed that all produced consistent results.

<footnote><sup>a</sup> The statistical uncertainties on the fitted parameters were typically ≲1%, but this does not account for the systematic uncertainties associated with extracting the emission surfaces, which, although more difficult to quantify, are likely substantially larger than those related to the analytical fits.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Peak intensity maps of CO 2–1 for all MAPS sources, with overlaid contours showing the fitting emission surfaces, as listed in Table 1. The synthesized beam and a scale bar indicating 50 au are shown in the lower left and lower right corners, respectively, of each panel.</figcaption></figure>

We found that fits to the moving average surfaces were slightly more reliable at the outer disk radii. This is likely because the moving average is a compromise in terms of the number of radial points and relative signal-to-noise ratio (S/N), as the raw data have a finer grid of radial points but much larger scatter, while the binned data have coarser radial information but higher S/N.

The fitted $\phi$ (~1.5–4) for the CO emission surfaces are often several times larger than the flaring indices of the gas pressure scale heights (1.08–1.35), which are derived by fitting the observed spectral energy distributions of each star+disk system (for more details, see Zhang et al. 2021). This difference was previously noted in IM Lup by Pinte et al. (2018), who ascribed it to the sharp drop-off in UV radiation from the central star. The stellar irradiation determines the shape of the emitting surface, which follows a layer of approximately constant optical depth from the perspective of the star, rather than tracing the disk scale height. In particular, in the inner <300 au (<120 au for AS 209) where the surfaces are still steeply rising, the CO surface originates from a height that is 2.5–3.5 times the scale height, while the ¹³CO and C¹⁸O surfaces are at approximately 1–1.5 scale heights. Similar values have been reported previously for IM Lup (Pinte et al. 2018), as well as in DM Tau (Dartois et al. 2003) and the Flying Saucer (Dutrey et al. 2017).

To further illustrate the geometry of the fitted surfaces, Figure 3 shows an overlay of the inferred emission surfaces on the peak intensity maps of CO 2–1 for all disks. Isovelocity

contours generated using the surface fits from Table 1 for the CO isotopologues are also provided in Appendix B.

# 4. Results

## 4.1. Overview of Emission Surfaces

Figure 4 shows the surfaces derived for CO isotopologues in all disks. For each disk, the CO 2–1 surface lies higher than that of the ¹³CO 2–1, which, in turn, is higher than C¹⁸O 2–1. Such a progression is consistent with a line optical depth of ~1 being reached at deeper layers for rarer isotopologues. There are considerable differences in the absolute surface heights of CO and ¹³CO between disks. For instance, the CO surface reaches a peak height of $z \approx 200$ au in IM Lup and GM Aur, while they are below ~100 au in AS 209, HD 163296, and MWC 480. A similar trend is seen in ¹³CO, where IM Lup has a maximum height of $z \approx 100$ au, while AS 209 and MWC 480 peak at $z < 40$ au. This range in absolute emission heights translates into a range of peak $z/r$. CO emission is present at $z/r \gtrsim 0.5$ in IM Lup and GM Aur, while $z/r \sim 0.2$ in AS 209. ¹³CO shows less overall variation between disks than CO and is generally present at $z/r < 0.2$. C¹⁸O has $z/r \lesssim 0.1$ toward those disks where we had enough signal to estimate emission heights. Finally, the relationship among the CO, ¹³CO, and C¹⁸O emission heights within disks varies across the sample. In MWC 480, the CO emission surface is relatively elevated with $z/r \sim 0.3$, while ¹³CO and C¹⁸O are both very flat, i.e., $z/r < 0.1$. By contrast, HD 163296 shows a gradual
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Emission surfaces for CO 2–1, $^{13}$CO 2–1, and C$^{18}$O 2–1 in all MAPS sources. 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 1. The solid lines show the radial range used in the fitting, while the dashed lines are extrapolations. The C$^{18}$O 2–1 surfaces in IM Lup and AS 209 are tentative and were derived from the tapered (0 ′′ 30) images. Lines of constant z/r are shown in gray. The FWHM of the major axis of the synthesized beam is shown in the lower right corner of each panel.</figcaption></figure>

progression of z/r ~ 0.3 to 0.2 to 0.1 for CO, $^{13}$CO, and C$^{18}$O,
respectively.

No surfaces could be derived from the full-resolution images
of C$^{18}$O 2–1 in IM Lup and AS 209 owing to insufficient S/N.
However, we were able to extract surfaces from the corres-
onding tapered (0 ′′ 30) image cubes (see Section 6.2, Czekala
et al. 2021) but consider these to be tentative and did not
attempt to fit analytical functions to these surfaces. Both are
shown in Figure 4 but are otherwise omitted from subsequent
analysis.

In all CO lines, except for CO in IM Lup, we see an initial
increase of z/r with radius, i.e., flaring, a flattening, and then
eventual turnover due to decreasing gas surface densities at
large radii. IM Lup, which is known to possess extended
diffuse CO emission (Cleeves et al. 2016), does not show clear
evidence of this turnover and only shows moderate indications
of flattening. All surfaces show some degree of vertical scatter,
which is a combination of thermal noise in the images and
potential azimuthal variations in the underlying emission
surfaces. The relative contribution is, however, specific to
each disk and emission line. As this scatter increases
substantially with radius, the line S/N is likely the most
important factor in setting the vertical scatter, at least at large
radii. Due to different projections at varying azimuths in
specific channels, the height of a particular pixel can often be
easier or more difficult to determine, which provides an
additional source of uncertainty in vertical pixel positions. For
instance, channels with less favorable viewing geometries
where the two disk sides cannot be easily distinguished make it
harder to measure emission surface heights. This often occurs
at velocities either very close to or substantially offset from the
source systemic velocity. For example, in Figure 16, channels
at larger velocity offsets (≤2.6 km s⁻¹ or ≥6.2 km s⁻¹) and
those near the systemic velocity (≈4.5 km s⁻¹) show poorly
separated upper and lower disk surfaces.

Nonetheless, some surfaces appear more tightly constrained
than others, i.e., CO in HD 163296 shows considerably less
vertical dispersion compared to that of the MWC 480 disk.
This scatter in the MWC 480 disk is not just due to noise but is

the result of localized azimuthal deviations. Perturbations in
azimuthal velocity, on the order of a few percent, are located at
~240, 340, 370, and 450 au in MWC 480 (Teague et al. 2021),
which approximately align with regions of prominent vertical
scatter in its CO emission surface. Similarly detailed and disk-
specific analyses are required to discern the origins of vertical
scatter in the other MAPS sources. Several disks also show
evidence of substructure in their surfaces, e.g., C¹⁸O in
HD 163296, which we discuss in detail in Section 4.4.

Figure 5 shows the surfaces for C<sub>2</sub>H and HCN. Of the disks around T Tauri stars, only surfaces for AS 209 could be extracted and appear to be at z/r ~ 0.1. The C<sub>2</sub>H and HCN surfaces in IM Lup and GM Aur could not be reliably constrained owing to their low line optical depths and S/Ns compared to CO and <sup>13</sup>CO. The two disks around Herbig Ae stars, HD 163296 and MWC 480, also show emission at a z/r of 0.1 or less. In MWC 480, both HCN and C<sub>2</sub>H are present at z/r < 0.1, similar to the <sup>13</sup>CO and C<sup>18</sup>O surfaces. HD 163296 is the only source where the C<sub>2</sub>H and HCN lines show any structure; there is a clear gap in the surfaces corresponding to the gap between the two innermost rings in the radial emission profiles (Law et al. 2021). The first ring at 45 au is less vertically extended with z/r < 0.1, while the emission in the second ring at 110 au is more elevated at z/r ≈ 0.1. We do not attempt parametric fits for any of the HCN and C<sub>2</sub>H lines.

4.2. *Comparison with NIR Rings*

All of the MAPS sources have been observed in scattered
light (Schneider et al. 2003; Kusakabe et al. 2012; Monnier
et al. 2017; Avenhaus et al. 2018; Muro-Arena et al. 2018),
which provides valuable information about the micron-sized
dust grains in these disks. The IM Lup, AS 209,<fnref n="28" /> and
HD 163296 disks have well-defined rings in the NIR, but only

<footnote n="28">²⁸ If deprojected with a nonzero flaring angle, Avenhaus et al. (2018) found that AS 209 possesses either one (112 au) or three (78, 140, and 243 au) NIR rings depending on whether the northern side is the near or far side, respectively. Subsequent observations (Guzmán et al. 2018; Teague et al. 2018b) showed that the latter interpretation is correct.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Emission surfaces for HCN 3–2 and C$_2$H 3–2 in AS 209, HD 163296, and MWC 480. Large gray points show radially binned surfaces, and small, light-gray points represent individual measurements. Lines of constant z/r are shown in gray. The FWHM of the major axis of the synthesized beam is shown in the lower right corner of each panel. The HCN and C$_2$H surfaces for the remaining MAPS sources IM Lup and GM Aur lacked sufficient S/N for robust surface determinations.</figcaption></figure>

IM Lup and HD 163296 have direct estimates of their NIR-emitting surfaces, as measured from individual rings. The inner NIR ring in HD 163296 has a height measured from Monnier et al. (2017), while the outer ring at 330 au was recently found to have a dust scale height of 64 au in NIR/HST observations (Rich et al. 2020). All four rings in IM Lup have measured NIR heights (Avenhaus et al. 2018). Figure 6 shows these NIR heights compared to the CO and $^{13}$CO 2–1 emission surfaces. We also plot the NIR-emitting height relation identified in a sample of disks around T Tauri stars as part of the DARTTS-S program (Avenhaus et al. 2018) as a dashed red line in Figure 6.

The NIR surfaces lie between the CO- and $^{13}$CO-emitting layers in HD 163296, while in IM Lup the NIR surface appears at approximately the same height as that of the $^{13}$CO, which is roughly consistent with the findings from Pinte et al. (2018) and Rich et al. (2021). Although lacking well-defined rings, MWC 480 has been reported to have a very flat NIR surface, i.e., $z/r \sim 0.03$ (Kusakabe et al. 2012), which suggests that the micron-sized dust lies at or below the $^{13}$CO- and C$^{18}$O-emitting layers.

### 4.3. Gas Temperatures

We can use line-emitting surfaces together with line brightness temperatures to map disk temperature structures. When we extracted individual pixels from the image cubes, we also obtained a corresponding set of peak surface brightnesses. In Section 4.3.1 we describe how we converted these peak surface brightnesses into gas temperatures as a function of (r, z). Then, in Section 4.3.2 we present the radial temperature profiles, and in Section 4.3.3 we analyze the full 2D empirical temperature structure of each disk. Both the radial temperature profiles and full (r, z) temperature structures for each MAPS are provided as publicly available VADPs (see Section 7).

#### 4.3.1. Calculating Gas Temperatures

The peaks of the CO and $^{13}$CO 2–1 lines are expected to be optically thick with CO rotational levels in local thermodynamic equilibrium (e.g., Weaver et al. 2018) at the typical densities and temperatures of protoplanetary disks. Provided that the emission fills the beam, the peak surface brightness $I_ν$ provides a measure of the temperature of the emitting gas. In

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Emission surfaces of CO and $^{13}$CO 2–1 in IM Lup (top) and HD 163296 (bottom). The black lines are the moving average surfaces, and gray shaded regions show the $1\sigma$ uncertainty. The red diamonds show individual height measurements of NIR rings (Monnier et al. 2017; Avenhaus et al. 2018; Rich et al. 2020), while the red dashed line shows the inferred NIR surface using the relation found in a sample of disks in Avenhaus et al. (2018). The values from Monnier et al. (2017) have been scaled to reflect the updated Gaia distance for HD 163296. The error bars are smaller than the marker for the innermost rings in both IM Lup and HD 163296, and the 330 au ring in HD 163296 does not have reported uncertainties.</figcaption></figure>

order to not underestimate the line intensity along lines of sight containing optically thick continuum emission (e.g., Boehler et al. 2017; Weaver et al. 2018), we repeated the surface fitting, as in Section 3, using the non-continuum-subtracted image cubes (Czekala et al. 2021).

Each individual pixel (r, z) that was extracted has a peak surface brightness $I_ν$, which was then used to calculate the associated gas temperature using the full Planck function:

$$T_b = \frac{h\nu}{k_B} \left[ \ln \left( \frac{2h\nu^3}{c^2 I_\nu} + 1 \right) \right]^{-1} \quad (2)$$

In addition to CO 2–1 and $^{13}$CO 2–1, we also calculated the brightness temperatures of C$^{18}$O 2–1 in all disks and those of
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Radial brightness temperature profiles for CO 2–1, $^{13}$CO 2–1, and, when available, C$^{18}$O 2–1. The top and bottom rows are each shown on a consistent temperature scale and are grouped by whether the host star is a T Tauri or Herbig Ae star, respectively. 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 4.3.1). Vertical lines show the $1\sigma$ uncertainty, given as the standard deviation of the individual measurements in each bin. For increased visual clarity, all disks have been binned by an extra factor of two, except for IM Lup and HD 163296, which have been binned by an additional factor of three and four, respectively. The solid red line shows the fitted power-law profile from Table 2. The inner gray shaded region is the FWHM of beam major axis. The orange and gray dashes correspond to the millimeter continuum rings and gaps, respectively. The black circles mark the outer edge of the millimeter continuum. Temperature dips in GM Aur, AS 209, and HD 163296 are labeled with arrows, as is a temperature bump in IM Lup.</figcaption></figure>

HCN 3–2 and C$_2$H 3–2 in HD 163296, but as we expect these lines to be partly optically thin, their brightness temperatures will be lower limits on the gas temperatures.

The western half of the AS 209 disk suffers from foreground cloud contamination (Öberg et al. 2011) in CO 2–1. Therefore, we calculated CO 2–1 temperatures using only the eastern half of the disk, which corresponds to the velocity range of 4.90–6.90 km s$^{-1}$ (see Appendix B), to avoid underestimating the peak brightness temperatures. For all other lines, we used the same velocity/channel ranges as in Table 1 for temperature calculations.

All subsequent radial and 2D gas temperature distributions represent those derived directly from individual surface measurements, rather than radially deprojecting peak intensity maps (see Teague et al. 2021) or mapping peak brightness temperatures back onto derived emission surfaces (i.e., Figure 3). We only consider the brightness temperatures of those pixels where we were able to determine an emission height.

### 4.3.2. Radial Temperature Distributions

Figure 7 shows the radial temperature profiles for the CO isotopologues. We first reiterate that these temperatures are measurements of surface brightnesses, rather than integrated intensities that are used to identify line emission substructures (Law et al. 2021) or derive column densities (Zhang et al. 2021). As expected, in each disk, CO is the warmest, followed by $^{13}$CO and then C$^{18}$O. CO displays the largest range of measured temperatures, while $^{13}$CO and C$^{18}$O span a more

limited range. The radial temperature gradients are consistent within each disk with similar slopes across CO isotopologues with the exception of AS 209, where the $^{13}$CO is nearly flat over the entire radial range in which it was measured. This flatness in temperature structure is due to the $^{13}$CO emission rings in AS 209, and in this case, we are only able to derive the brightness temperature of the outer ring at ~120 au (Favre et al. 2019; Law et al. 2021).

The disks around T Tauri stars have brightness temperatures spanning ~10–40 K. The disks around Herbig Ae stars HD 163296 and MWC 480 are generally warmer at a given radius and have an overall larger total temperature range from ~10 to 70 K. The CO 2–1 temperatures are about 10 K higher in MWC 480 than in HD 163296, with the greatest differences occurring within <200 au. In particular, the HD 163296 and MWC 480 profiles are consistent with those presented in Teague et al. (2021), which were instead generated by deprojecting the peak intensity maps rather than direct extraction from emitting surfaces. Likewise, the CO 2–1 temperature profile of HD 163296 is approximately consistent with, although slightly cooler than, the one derived from a similar direct extraction method in Isella et al. (2018). We note that brightness temperatures less than 20 K are below the CO freezeout temperature, which suggests that the associated line emission is at least partially optically thin and thus only provides a lower limit on the true gas temperatures. This conclusion is supported by our data, where CO lines with <20 K are most common for the rarer isotopologues and at large disk radii.
The drop in brightness temperature seen within 20–40 au in all disks and lines, which is marked as a shaded region in Figure 7, is due to beam dilution as the emitting area becomes comparable to or smaller than the angular resolution of the observations. In the case of IM Lup and AS 209, the central temperature dip extends further than the beam size. This may suggest enough CO depletion for the lines to become optically thin at these innermost radii, unresolved CO emission substructure, or that a substantial fraction of the CO emission is absorbed by dust. Indeed, Cleeves et al. (2016) and Sierra et al. (2021) find that the dust is optically thick in the inner regions of the IM Lup disk, and Bosman et al. (2021) also see a large CO emission gap that is best explained by dust absorption. Dust absorption may also contribute to the low CO temperature in the inner AS 209 disk, but not out to 100 au. Optically thin emission is also an unlikely explanation: Zhang et al. (2021) find that while AS 209 has a lower CO surface density than all other MAPS disks, it is still far from the optically thin limit in CO 2–1. This leaves CO substructure as an explanation. AS 209 does present several gaps in CO emission interior to 100 au (Guzmán et al. 2018; Law et al. 2021; Zhang et al. 2021; Bosman et al. 2021), which are barely resolved and may therefore result in a low brightness temperature in the inner disk.

Figure 7 also shows the locations of millimeter continuum gaps and rings, as reported in Law et al. (2021). The radial temperature profiles are quite smooth, and hence the opportunity for coincidences between temperature substructures and other substructures is small. In two cases, the temperature substructure that is seen does line up with known disk substructures, however. HD 163296 has a 5 K drop in temperature at ~80–90 au in all three CO lines, which aligns with a gap at 85 au in the millimeter continuum. A similar drop in temperature at ~90 au is also present in CO in GM Aur, which roughly aligns with a gap–ring pair at 68 and 86 au in the millimeter continuum. In AS 209, a slightly deeper (~8 K) drop occurs at 200 au in CO 2–1 and is coincident with a CO 2–1 line emission gap at 197 au (Law et al. 2021). Low-amplitude (~2–3 K) wave-like fluctuations are seen in CO 2–1 temperature in MWC 480 (for further discussion of the features, see Teague et al. 2021). We find no association between temperature trends and the outer continuum edge in GM Aur and AS 209, but we do notice a modest flattening of the CO temperature gradient at the edge of the millimeter continuum in HD 163296 and MWC 480. Although about 50 au beyond the continuum edge, the CO 2–1 temperature in IM Lup shows a modest increase at 450 au. This may be associated with a temperature inversion in the midplane (Cleeves 2016) and is broadly consistent with the radial location of 400 au predicted in the models of Facchini et al. (2017).

We fitted the temperature profiles with power laws as

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

We first visually chose the radial range in which the temperature profiles behave like a power law and then fitted each profile using the Levenberg–Marquardt minimization implementation in `scipy. The fitting ranges and derived parameters are listed in Table 2.

As shown in Figure 7, these fits work well beyond ~100–150 au in all disks but overpredict the measured

brightness temperatures interior to this. The temperature profile of CO in MWC 480 changes slope at approximately 350 au, which complicates the choice of radial fitting range. It is possible to achieve a modestly more accurate fit if instead two power laws are used, one for the inner disk between 40 and 200 au and another for the outer disk from >135 au. However, for simplicity, we fit a single power to the maximal possible range.

In addition to the CO lines, we also derived the brightness temperature profiles of HCN and C$_2$H in HD 163296, as shown in Figure 8. The shapes of the temperature profiles are consistent but offset, as HCN is warmer by 4–6 K at all radii. The gap at ~80 au has a C$_2$H temperature of <15 K and HCN temperature of <20 K. Both lines seem to be cooler near the gap by a few kelvin relative to the outer ring and by almost 10 K versus the inner ring. However, the beam filling factor will be reduced at locations closer to the gap, likely becoming significant within 1/2–1 beam away. Thus, the line emission near this gap may become increasingly optically thin. In this case, the lower brightness temperatures would reflect reduced gas density, rather than cooler HCN and C$_2$H gas temperatures. Overall, the HCN brightness temperatures are consistent with the excitation temperatures derived in the multiline analysis of Guzmán et al. (2021). However, the C$_2$H temperatures are a factor of two lower than those reported in Guzmán et al. (2021), which suggests that the C$_2$H 3–2 line is optically thin or not in local thermal equilibrium (LTE). Thus, a non-LTE analysis of C$_2$H in HD 163296 is warranted.

### 4.3.3. 2D Temperature Structure

The advantage of having multiple CO isotopologues that trace different disk heights is access to the vertical temperature distribution. Dartois et al. (2003) were the first to demonstrate this in moderate-resolution (~1") observations of DM Tau. More recently, Pinte et al. (2018) presented a framework for directly mapping the temperature structure of each emitting layer in (~0"4) observations of IM Lup. Here, we expanded this analysis to the high spatial resolution observations of the MAPS disks. Figure 9 shows the full 2D temperature distributions.

Since we have temperature information as a function of (r, z), we can construct a full 2D model of the temperature distribution of each disk. To do so, we adopt a two-layer model similar to the one proposed by Dartois et al. (2003) but then modified by Dullemond et al. (2020) with a different connecting term. Both formulae were initially tried, but substantially better fits were obtained with that of Dullemond et al. (2020). The midplane temperature $T_{\text{mid}}$ and atmosphere temperature $T_{\text{atm}}$ are assumed to have a power-law profile with slopes $q_{\text{mid}}$ and $q_{\text{atm}}$, respectively:

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

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

Between the midplane and atmosphere, the temperature is smoothly connected using a tangent hyperbolic function

$$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 (6)$$
<table><caption>**Table 2**

Radial Temperature Profile Fits</caption><thead><tr><th>Source</th><th>Line</th><th>r<sub>fit,in</sub> (au)</th><th>r<sub>fit,out</sub> (au)</th><th>T<sub>100</sub> (K)</th><th>q</th></tr></thead><tbody><tr><td rowspan="2">IM Lup</td><td>CO 2–1</td><td>170</td><td>559</td><td>55 ± 0.9</td><td>0.58 ± 0.01</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>145</td><td>339</td><td>30 ± 0.6</td><td>0.32 ± 0.03</td></tr><tr><td rowspan="3">GM Aur</td><td>CO 2–1</td><td>135</td><td>613</td><td>52 ± 0.9</td><td>0.61 ± 0.02</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>50</td><td>314</td><td>22 ± 0.2</td><td>0.26 ± 0.01</td></tr><tr><td>C<sup>18</sup>O 2–1</td><td>30</td><td>126</td><td>14 ± 0.4</td><td>0.38 ± 0.05</td></tr><tr><td rowspan="2">AS 209</td><td>CO 2–1</td><td>95</td><td>244</td><td>42 ± 1.0</td><td>0.78 ± 0.05</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>125</td><td>163</td><td>28 ± 1.3</td><td>0.80 ± 0.13</td></tr><tr><td rowspan="3">HD 163296</td><td>CO 2–1</td><td>150</td><td>527</td><td>78 ± 1.0</td><td>0.82 ± 0.01</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>50</td><td>356</td><td>31 ± 0.2</td><td>0.37 ± 0.01</td></tr><tr><td>C<sup>18</sup>O 2–1</td><td>40</td><td>148</td><td>21 ± 0.3</td><td>0.37 ± 0.03</td></tr><tr><td rowspan="3">MWC 480</td><td>CO 2–1</td><td>100</td><td>632</td><td>70 ± 1.1</td><td>0.69 ± 0.02</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>100</td><td>388</td><td>42 ± 0.9</td><td>0.60 ± 0.03</td></tr><tr><td>C<sup>18</sup>O 2–1</td><td>80</td><td>251</td><td>26 ± 0.2</td><td>0.75 ± 0.02</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Radial brightness temperature distributions of HCN 3–2 and C$_2$H 3–2 in HD 163296. 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 4.3.1). Vertical lines show the $1\sigma$ uncertainty, given as the standard deviation of the individual measurements in each bin. The gap location seen in the line emission (Law et al. 2021) is labeled and shaded in gray. The FWHM of the beam major axis is shown in the lower right corner and also indicated by the inner gray shaded region.</figcaption></figure>

where $z_q(r) = z_0(r/100 \text{ au})^\beta$. We note that the $\alpha$ parameter defines where in height the transition in the tanh vertical temperature profile occurs and $\beta$ describes how the transition height varies over radius. In total, we fitted the following seven parameters: $T_{\text{atm},0}$, $q_{\text{atm}}$, $T_{\text{mid},0}$, $q_{\text{mid}}$, $\alpha$, $z_0$, and $\beta$.

We performed the fitting using MCMC with emcee (Foreman-Mackey et al. 2013) with 256 walkers that take 500 steps to burn in and an additional 5000 steps to sample the posterior distribution function. All available CO lines were fitted using the individual measurements, and only those points with $T_B > 20$ K were considered. Temperatures below 20 K, close to the CO freezeout temperature, are likely optically thin and thus not useful for constraining the gas temperature structure. Parameter values and

associated uncertainties are taken to be the 50th, 16th, and 84th percentiles from the marginalized posterior distributions, respectively, and are listed in Table 3.

The 2D fitted models are shown in comparison with the data in Figure 10. For all disks, the median residuals between the fitted model and measured temperatures are typically no more than 10%. The most informative fits are those with a well-sampled $(r, z)$ space, which means that we have a set of CO isotopologue lines with a diverse set of $z/r$ values, e.g., HD 163296. In contrast, IM Lup is poorly constrained over the height of the disk, since surfaces were only able to be determined for CO and $^{13}$CO and they are not widely spaced in $z/r$. The abrupt change in $z/r$ from CO to $^{13}$CO and $C^{18}$O in MWC 480 is also reflected in its inferred 2D temperature structure by its small fitted $\alpha$ and $\beta$ values (Table 3). This means that the transition in vertical temperature, as described in Equation (6), occurs close to the midplane and the transition height does not increase over radius, unlike other disks. In general, as the emitting surfaces do not provide direct constraints in the disk midplanes, we caution the use of the empirically derived $T_{\text{mid}}$, which are considerably warmer than predictions from thermochemical models (Zhang et al. 2021).

## 4.4. Substructures in Emission Surfaces

Localized vertical substructures are observed in many of the emission surfaces derived from the MAPS data. The properties of these substructures, namely, their radial locations, widths, and depths, provide important constraints that are necessary for detailed thermochemical modeling (e.g., Rab et al. 2020; Calahan et al. 2021). In the following subsections, we identify and catalog all substructures present in the derived emitting layers and compare them with the gas temperature profiles and with substructures observed in the millimeter continuum and CO line emission.

### 4.4.1. Fitting Vertical Substructures

Each substructure is labeled with its radial location rounded to the nearest whole number in astronomical units and is preceded with “Z” to indicate that these features are vertical variations. This nomenclature is also chosen to avoid ambiguity with that used to denote radial substructures in the continuum (Huang et al. 2018) and molecular line emission (Law et al. 2021) profiles, which label rings by “B” (“bright”) and gaps by “D” (“dark”).
<figure><img src="image_12.png" /><figcaption>**Figure 9.** 2D temperature distributions of CO 2–1, $^{13}$CO 2–1, and, when available, C$^{18}$O 2–1 in all MAPS sources. Points are those from the binned surfaces, and error bars are the $1\sigma$ uncertainties in z. For some of the innermost points, the uncertainty is smaller than the marker. Data for all disks and lines have been binned by a factor of three for visual clarity. The uncertainty of the temperature measurements, which is not shown here, can be found in Figure 7.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison of the measured temperatures (points) with the fitted 2D temperature structures (background), as listed in Table 3. The same color scale is used for the data and fitted model and is consistent across all panels. Contours show constant temperatures in increments of 10 K. Data for all disks and lines have been binned by a factor of three for visual clarity. The uncertainty of the temperature measurements, which is not shown here, can be found in Figure 7.</figcaption></figure>

<table><caption>**Table 3**
Summary of 2D Temperature Structure Fits</caption><thead><tr><th>Source</th><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>IM Lup</td><td>36<sup>+0.1</sup><sub>-0.1</sub></td><td>25<sup>+0.1</sup><sub>-0.1</sub></td><td>0.03<sup>+0.01</sup><sub>-0.01</sub></td><td>-0.02<sup>+0.01</sup><sub>-0.01</sub></td><td>3<sup>+0.1</sup><sub>-0.1</sub></td><td>4.91<sup>+0.17</sup><sub>-0.16</sub></td><td>2.07<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>GM Aur</td><td>48<sup>+0.3</sup><sub>-0.3</sub></td><td>20<sup>+0.2</sup><sub>-0.2</sub></td><td>-0.55<sup>+0.01</sup><sub>-0.01</sub></td><td>-0.01<sup>+0.01</sup><sub>-0.01</sub></td><td>13<sup>+0.2</sup><sub>-0.2</sub></td><td>2.57<sup>+0.03</sup><sub>-0.03</sub></td><td>0.54<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>AS 209</td><td>37<sup>+0.2</sup><sub>-0.2</sub></td><td>25<sup>+0.2</sup><sub>-0.2</sub></td><td>-0.59<sup>+0.02</sup><sub>-0.02</sub></td><td>-0.18<sup>+0.01</sup><sub>-0.01</sub></td><td>5<sup>+0.2</sup><sub>-0.2</sub></td><td>3.31<sup>+0.12</sup><sub>-0.11</sub></td><td>0.02<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>HD 163296</td><td>63<sup>+0.2</sup><sub>-0.2</sub></td><td>24<sup>+0.1</sup><sub>-0.1</sub></td><td>-0.61<sup>+0.003</sup><sub>-0.003</sub></td><td>-0.18<sup>+0.004</sup><sub>-0.004</sub></td><td>9<sup>+0.1</sup><sub>-0.1</sub></td><td>3.01<sup>+0.02</sup><sub>-0.02</sub></td><td>0.42<sup>+0.004</sup><sub>-0.004</sub></td></tr><tr><td>MWC 480</td><td>69<sup>+0.2</sup><sub>-0.2</sub></td><td>27<sup>+0.2</sup><sub>-0.2</sub></td><td>-0.7<sup>+0.004</sup><sub>-0.004</sub></td><td>-0.23<sup>+0.01</sup><sub>-0.01</sub></td><td>7<sup>+0.1</sup><sub>-0.1</sub></td><td>2.78<sup>+0.02</sup><sub>-0.02</sub></td><td>-0.05<sup>+0.01</sup><sub>-0.01</sub></td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Example of the vertical substructure fitting process in the CO 2–1 emission surface of HD 163296. The top panel shows the removal of a local baseline (solid red line) in the form of a quadratic polynomial fit. The bottom panel shows a Gaussian fit (solid blue line) to the Z46 substructure in the baseline-subtracted emission surface. The radial range for the local baseline was visually identified and is shown by the orange shaded regions.</figcaption></figure>

Feature identification was done visually and focused on the inner, rising portion of the surfaces within ~200 au, which was the most well constrained and possessed the highest S/Ns. We used the moving average surfaces to search for substructures in the form of vertical dips, i.e., we assumed that substructures represent localized decreases in $z$ in an otherwise smoothly varying emitting layer. To fit each substructure, we first visually estimated a local baseline. This baseline was then fitted with a quadratic polynomial and subtracted from the original emitting surface. The derived properties of each feature will depend on the assumed form of the local baseline, but a low-order polynomial baseline is sensible for the inner <200 au of each disk. We then fitted a Gaussian profile to characterize each feature in the baseline-subtracted surface. An example of this fitting process for CO 2–1 in HD 163296 is shown in Figure 11.

The fitted centers and FWHMs of each Gaussian are taken to be the radial location and widths of each feature, respectively. Substructure depths are defined as $\Delta z/z_{\text{baseline}}$, where $z_{\text{baseline}}$ is the vertical height of the fitted baseline and $\Delta z$ is the fitted vertical height of the emitting surface at the radial position of the substructure. Depths are subsequently referred to according to their fractional decrease in vertical height with deeper features having lower height ratios, e.g., the Z46 in HD 163296 has a $\Delta z/z_{\text{baseline}} = 0.52$, which indicates a depth of 48%. The center, width, and relative depth of each feature are listed in Table 4, and their radial locations are labeled in Figure 12.

The relative depth of each feature is sometimes more uncertain than Table 4 suggests, as depth strongly depends on the assumed baseline, but overall we find that this method works well to identify substructure radial locations and provides a preliminary characterization. These definitions also do not explicitly account for beam effects. In a few cases, the widths and depths of individual features are smaller than the minor axis of the beam FWHM. However, this is not generally

a concern, as surfaces are derived from the positional offsets of peak intensities, which are sensitive to scales smaller than the beam size.

Typical feature depths range from ~30% to 70%, and widths range from 10 to 50 au. HD 163296 has the largest total number of identified substructures, which are relatively narrow (~10–15 au), while those in IM Lup and MWC 480 have broader widths (~30–50 au). A consistent broad, bowl-shaped depression is seen in all CO lines in MWC 480 around 66 au with a width of $\geq 30–40$ au and depth of 40%–60%. The two features Z170 and Z375 associated with IM Lup are notable, as they occur at the largest radii of all identified substructures. Although having modest relative depths of 68% and 49%, they possess an absolute $\Delta z$ of 46 and 18 au, which are the largest in au by a factor of a few to an order of magnitude compared to all other substructures. Features are not always present across all CO isotopologues. For instance, Z170 in IM Lup is only seen in $^{13}$CO but not in CO, while in HD 163296 Z81 and Z83 are present in $^{13}$CO and $^{18}$O, respectively, but no corresponding feature is identified in CO.

In addition to the isolated, Gaussian-like dips we report above, we detect a few more complex trends. For instance, prominent changes in the slope of the emission surface are present at ~150 au in $^{13}$CO in GM Aur and at ~115 au in CO in IM Lup. Both AS 209 and MWC 480 also show large-scale, wave-like patterns in their CO surfaces with peak-to-trough separations of roughly 20 and 40 au, respectively and amplitudes of no more than a few au. Thus, Z56 in AS 209 and Z66 in MWC 480 may in fact be local minima associated with this larger wave rather than separate isolated substructures.

### 4.4.2. Comparison with Gas Temperature

Since we have estimates for the gas temperatures, we also searched for coincidences between vertical and temperature substructures, but we only found a single one, i.e., Z81 in $^{13}$CO in HD 163296. In IM Lup, the outer edge of the CO temperature plateau occurs around 180 au, which is roughly coincident with the Z170 feature seen in the $^{13}$CO surface. The CO 2–1 radial temperature profiles and emitting surfaces in MWC 480 also both show wave-like patterns, which are very roughly coincident in radial location. For further discussion of these features, see Teague et al. (2021). Otherwise, no other temperature trends are identified in any of the CO lines in the MAPS disks. Thus, the empirical gas temperatures do not generally seem to be sensitive to the presence of surface substructures. This may be explained, in part, by the deeper layers in the disk being more isothermal than the upper layers. Thus, a small change in the emission height would look like a bigger dip in temperature for CO than $^{13}$CO, for instance. Moreover, since the scale height is proportional to $\sqrt{T}$, this means that significant temperature differences are required for noticeable changes in the scale height.

### 4.4.3. Comparison with Millimeter Continuum and Line Emission Substructures

The majority of the MAPS disks show at least some spatial links between vertical substructures and either continuum or radial substructure in CO line emission, as shown in Figure 12. Below, we consider each of these possible spatial associations.

All MAPS disks show some degree of spatial association between surface features and millimeter continuum substructures.
<table><caption>**Table 4**
Properties of Vertical Substructures</caption><thead><tr><th>Source</th><th>Line</th><th>Feature</th><th>r<sub>0</sub> (arcsec)</th><th>r<sub>0</sub> (au)</th><th>Width (arcsec)</th><th>Width (au)</th><th>Δz (arcsec)</th><th>Δz (au)</th><th>Depth<sup>a</sup></th></tr></thead><tbody><tr><td rowspan="2">IM Lup</td><td>CO 2–1</td><td>Z375</td><td>2.38 ± 0.03</td><td>375 ± 5</td><td>0.18 ± 0.09</td><td>29 ± 14</td><td>0.29 ± 0.1</td><td>46 ± 16</td><td>0.32 ± 0.28</td></tr><tr><td><sup>13</sup>CO 2–1</td><td>Z170</td><td>1.08 ± 0.01</td><td>170 ± 2</td><td>0.18 ± 0.002</td><td>29 ± 0.4</td><td>0.11 ± 0.004</td><td>18 ± 1</td><td>0.51 ± 0.49</td></tr><tr><td rowspan="2">AS 209</td><td>CO 2–1</td><td>Z56</td><td>0.47 ± 0.02</td><td>56 ± 2</td><td>0.1 ± 0.01</td><td>13 ± 1</td><td>0.03 ± 0.02</td><td>4 ± 2</td><td>0.33 ± 0.18</td></tr><tr><td>CO 2–1</td><td>Z46</td><td>0.45 ± 0.01</td><td>46 ± 1</td><td>0.11 ± 0.01</td><td>12 ± 1</td><td>0.06 ± 0.01</td><td>6 ± 1</td><td>0.52 ± 0.13</td></tr><tr><td rowspan="5">HD 163296</td><td><sup>13</sup>CO 2–1</td><td>Z49</td><td>0.48 ± 0.02</td><td>49 ± 2</td><td>0.1 ± 0.003</td><td>11 ± 0.3</td><td>0.02 ± 0.004</td><td>2 ± 0.4</td><td>0.36 ± 0.03</td></tr><tr><td></td><td>Z81</td><td>0.8 ± 0.001</td><td>81 ± 0.1</td><td>0.11 ± 0.03</td><td>12 ± 3</td><td>0.05 ± 0.03</td><td>5 ± 4</td><td>0.43 ± 0.24</td></tr><tr><td></td><td>Z145</td><td>1.43 ± 0.01</td><td>145 ± 1</td><td>0.1 ± 0.01</td><td>10 ± 1</td><td>0.02 ± 0.03</td><td>2 ± 3</td><td>0.33 ± 0.28</td></tr><tr><td>C<sup>18</sup>O 2–1</td><td>Z83</td><td>0.82 ± 0.004</td><td>83 ± 0.4</td><td>0.16 ± 0.02</td><td>16 ± 2</td><td>0.06 ± 0.001</td><td>6 ± 0.1</td><td>0.66 ± 0.02</td></tr><tr><td>CO 2–1</td><td>Z66</td><td>0.41 ± 0.01</td><td>66 ± 1</td><td>0.19 ± 0.04</td><td>31 ± 7</td><td>0.04 ± 0.003</td><td>7 ± 1</td><td>0.42 ± 0.06</td></tr><tr><td rowspan="2">MWC 480</td><td><sup>13</sup>CO 2–1</td><td>Z66</td><td>0.41 ± 0.06</td><td>66 ± 9</td><td>0.28 ± 0.04</td><td>46 ± 7</td><td>0.03 ± 0.02</td><td>6 ± 4</td><td>0.63 ± 0.09</td></tr><tr><td>C<sup>18</sup>O 2–1</td><td>Z71</td><td>0.44 ± 0.03</td><td>71 ± 5</td><td>0.32 ± 0.06</td><td>51 ± 9</td><td>0.03 ± 0.003</td><td>4 ± 0.5</td><td>0.63 ± 0.55</td></tr></tbody></table>

**Note.**

Each dust gap in HD 163296 aligns with a surface feature in at least one of the CO isotopologue surfaces, and in the case of ¹³CO, there is a one-to-one match between millimeter gaps and surface substructures. In MWC 480, the inner dust gap D76 roughly aligns with the Z66 surface feature, and the outer dust gap at D149 approximately matches the location of the 140 au trough of the wave-like fluctuations. However, these associations in MWC 480 are considerably more tentative, and considering that the surface features are more than twice the width of continuum gaps, these may be chance alignments. Several of the wave-like surface features in AS 209 align with the radial locations of substructures in the millimeter continuum. In particular, Z56 and the troughs at ~100 and 140 au are radially coincident with dust gaps. The changes in emitting surface slope in CO 2–1 in IM Lup and ¹³CO 2–1 in GM Aur are also both colocated with the D116 and D142 dust gaps, respectively.

In a few cases, we identify features in the CO line emission profiles that are radially coincident with vertical substructures. In HD 163296, the CO line emission peaks at B49 and B81 directly align with the vertical substructures at Z49 and Z81, respectively. The CO peak at B59 is also colocated with the Z56 feature in AS 209. Both changes in slope in IM Lup and GM Aur are spatially associated with CO features at B106 and B137, respectively. In IM Lup, we also find that the D360 gap in line emission may be spatially related to the Z375 vertical dip. Given that both line emission peaks and gaps show some spatial association with surface substructures, this may point toward multiple mechanisms producing these surface features.

# 5. Discussion

## 5.1. Comparison to Previous Results

Pinte et al. (2018) were the first to demonstrate an approach to directly extract CO emission surfaces in moderate-resolution observations of IM Lup. Similar methods were then used by Teague et al. (2018a, 2018b, 2019) and Rich et al. (2021) to constrain the emitting layers in AS 209 and HD 163296. Below, we compare these previous results to the high spatial resolution MAPS observations and comment on the relative consistency and any salient differences.

### 5.1.1. IM Lup

Pinte et al. (2018) found CO and ¹³CO 2–1 emission heights of z/r ≈ 0.325 and z/r ≈ 0.125 in the inner flared disk region,

i.e., r < 300 au, using ~0.′4 resolution observations. In this same region, we find a considerably more elevated CO surface of z/r ≳ 0.5 and ¹³CO of z/r ≈ 0.2. The differences in the derived heights are likely the result of our improved angular resolution, allowing for better separation of the front and back disk sides, e.g., see Figure 17, which demonstrates the limitations of poorer spatial resolutions. Pinte et al. (2018) also found that the CO 2–1 surface flattens out beyond 300 au, while our higher-sensitivity observations reveal that this is due to a local minimum (Z375) and that globally the CO surface continues to rise out to 550 au owing to diffuse large radii CO emission (e.g., Cleeves et al. 2016). Pinte et al. (2018) do not detect vertical substructures in any of their CO surfaces, but this is likely a consequence of their modest spatial resolution.

We also find consistent brightness temperature distributions with those of Pinte et al. (2018), but with overall systematically warmer temperatures by ~5 K. This systematic offset is likely due to beam dilution in Pinte et al. (2018), whose larger beam (~4×) would have smeared out some of this emission. Overall, these comparisons illustrate the importance of high spatial resolutions in accurately constraining emission surfaces. At more moderate spatial resolution, estimates of the heights of emitting layers and gas temperatures are both underestimated. A detailed exploration of the effect of spatial resolution on extracted surfaces within the MAPS disk is found in Appendix D. In short, we find that surfaces derived from images with beam sizes between 0.′′12 and 0.′′2 are consistent, suggesting that the results presented here are not underestimated owing to insufficient angular resolutions.

Rich et al. (2021) fit the CO 2–1 emission surfaces in the IM Lup disk using the DSHARP data cubes (Andrews et al. 2018) with a similar extraction method. These images have a spatial resolution of ~0.′′12, which is comparable to the MAPS resolution. Their CO 2–1 surface is nearly identical to the one we derived with the MAPS data and shows the same slope change at ~110 au and a localized dip around 375 au.

### 5.1.2. AS 209

Teague et al. (2018b) derived the CO 2–1 emission surface in the AS 209 disk using the same approach as Pinte et al. (2018) using high spatial resolution (~0.′′2) CO line data. The authors found a surface with z/r ~ 0.2, which is the same z/r we derived. However, they find a continually rising surface out to 300 au, but our surface begins to plateau and turn over at

<footnote>a Depth of vertical substructure, defined as the ratio of Δz to z<sub>baseline</sub> at r<sub>0</sub> (see Section 4.4.1).</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison of CO 2–1, $^{13}$CO 2–1, and C$^{18}$O 2–1 (when available) surfaces with annular substructures in the millimeter continuum (left) and CO line emission (right). Substructures are labeled following the nomenclature of Huang et al. (2018) and taken from Law et al. (2021). The CO line emission substructures are from CO 2–1 for all disks, except for MWC 480, which shows those from $^{13}$CO 2–1. Substructures are labeled as follows: dust ring (solid orange line), dust gap (dashed gray line), chemical ring (solid blue line), and chemical gap (dashed purple line). Dotted black lines mark vertical substructures, as listed in Table 4. Notable changes in emitting surface slope and the suggestive wave-like features in AS 209 and MWC 480 are marked with arrows. The CO velocity kinks from Pinte et al. (2018, 2020) are labeled in IM Lup and HD 163296.</figcaption></figure>
≥150 au. This difference may be due to their factor of two coarser spatial resolution. Teague et al. (2018b) also model the CO 2–1 emission surface that best reproduces the observed deviations in rotation velocities. Unlike the directly mapped surface, which is mostly smooth, the modeled surface has wave-like vertical substructures that appear very similar to those seen in the MAPS CO 2–1 emitting surface.

### 5.1.3. HD 163296

Teague et al. (2019) mapped the CO 2–1 emission surface in the HD 163296 disk at an angular resolution of ~0"/1. The authors constrained an emission surface out to a radius of 4" that is nearly identical to the one that we derive. In the same disk but using lower-resolution observations (~0"/25), Teague et al. (2018a) modeled the emitting layer of $C^{18}O$ 2–1 as a Gaussian process and found a typical $z/r \approx 0.1$, consistent with our $C^{18}O$ surface. Teague et al. (2018a) also found slight dips in their emission surface at the millimeter gap locations, i.e., ~50, 80, and 130 au (once rescaled to the Gaia distance). The first two depressions correspond to the Z49 and Z81 dips in our surface, while the 130 au dip lies beyond the turnover region in our data. The discrepancy at large radii between these surfaces may be due in part to their larger beam size (about twice that of the MAPS beam), as well as a different approach to surface extraction. Nonetheless, they are still broadly consistent with one another.

Rich et al. (2021) also used the DSHARP images (~0"/1) to extract the CO 2–1 surface of HD 163296. The authors find a surface that is very consistent, i.e., $z/r \approx 0.35$, with that derived from the MAPS data. However, the authors do not resolve the inner dip at 46 au that is clearly seen in the MAPS surfaces.

## 5.2. Origins of Emission Layer Heights

While the majority of emitting surfaces show a consistent general behavior, i.e., are well described by an exponentially tapered power-law profile, there is considerable variation in their emitting heights. Here, we briefly explore possible mechanisms that may be important in setting the heights of disk emitting layers. Specifically, to explore the origins of the observed diversity in $z/r$ structure across disks, i.e., from $z/r \gtrsim 0.5$ to $z/r \lesssim 0.1$, we searched for trends between physical parameters among the MAPS sources. While differences between isotopologues within disks are expected, differences between disks in the same isotopologues reveal variations in radiation fields and thermal, density, or CO abundance structures. To ensure a consistent comparison, we focus on the typical $z/r$ in the inner 150 au of each disk. We only consider the CO surfaces, as they are well constrained in all disks and show a sufficiently wide range of $z/r$ values.

As irradiation from the central star plays a large role in setting the shape of the emitting layers (Dullemond et al. 2001; Dullemond & Dominik 2004a, 2004b), we first consider whether differences in $z/r$ could be explained by differences in incident stellar radiation. In Figure 13, we identify a tentative negative trend between bolometric stellar luminosity and the $z/r$ of $^{12}$CO 2–1 emission surfaces, as well as a modest positive association between X-ray luminosity and emission height. However, in both cases AS 209 is an obvious outlier with an emission surface that is substantially flatter than the other two T

Tauri sources IM Lup and GM Aur, which possess similar stellar and X-ray luminosities.

Another parameter that may set disk emitting-layer heights is the temperature of the vertically isothermal layer (e.g., Qi et al. 2019). To check this, we compared the $^{13}$CO 2–1 gas temperatures at 100 au with emission heights in Figure 13. With the exception of AS 209, we find a negative trend, where the warmer temperatures of the two disks around the Herbig Ae stars have flatter surfaces, while the cooler temperatures of those around the T Tauri stars IM Lup and GM Aur have higher $z/r$ surfaces. We find a similar association if we instead consider the midplane temperature estimates derived from the thermochemical models of Zhang et al. (2021).

We next consider the physical properties of the gas itself, i.e., total disk size and column densities. We find a tight positive trend between CO 2–1 disk size and the $z/r$ of CO surfaces, as shown in Figure 13. However, the surfaces of small disks turn over at smaller radii, which may affect the $z/r$, but this trend remains unchanged if we instead compare using the turnover radius for each disk. Assuming that $H_2$ number density scales with CO column densities, we expect more dust grains in the disk upper layers due to increased dynamical gas-grain coupling. This may, in turn, manifest as higher $z/r$ surfaces. However, if we compare peak CO column densities, i.e., <50 au, versus $z/r$, we find an inverse relation, with AS 209 being a notable outlier to this trend. Taken together, this suggests that larger disks with lower column densities may preferentially exhibit elevated emitting surfaces.

However, these source characteristics are not all independent, since the mass of the central star either sets or influences many of them. Therefore, we also compared the stellar mass and emission surface $z/r$ in Figure 13 and found a negative trend, very similar to that of the stellar luminosity. The stellar mass sets the stellar luminosity, including the X-ray luminosity (with lower-mass stars being more active), which in turn controls the disk temperature structure. Physically, warmer disks may be expected to result in increasingly flared surfaces, but this is the opposite of what we find in Figure 13. As vertical surfaces are set by the balance of pressure and gravity, disks around more massive stars should, in contrast, exhibit flatter surfaces. Since stellar mass positively correlates with disk mass$^{29}$ (e.g., Andrews et al. 2013; Ansdell et al. 2016; Pascucci et al. 2016), this scenario is consistent with the observed trends and suggests that stellar mass is the dominant factor in setting emission surface heights. Thus, the majority of observed trends may simply be tracing the impact of varying stellar masses and the effects on the surrounding disks.

Overall, however, we caution that this small and highly biased sample of disks limits generalized conclusions. A survey aimed at targeting CO lines in a large set of moderately inclined disks with sufficient resolution and sensitivity is needed to provide further constraints on the origins and distribution of the heights of disk emitting layers.

## 5.3. Origins of Vertical Substructures

The emitting surfaces show several dips in vertical heights, which may have their origins in a variety of mechanisms. They may be due to CO depletion, i.e., decreased CO column

<footnote>$^{29}$ Literature $M_*$–$M_{\text{disk}}$ correlations are typically derived in the optically thin limit, but as disk continuum emission may be partially optically thick, the estimated disk masses should be considered lower bounds (e.g., Zhu et al. 2019; Andrews 2020).</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Average $z/r$ of $^{12}$CO 2–1 emission heights vs. disk physical properties for the MAPS sources. Stellar masses and bolometric luminosities are taken from Öberg et al. (2021) and references therein, while references for X-ray luminosities are as follows: IM Lup (Cleeves et al. 2017), AS 209 (Walter & Kuhi 1981), GM Aur (Espaillat et al. 2019), HD 163296 (Günther & Schmitt 2009), and MWC 480 (Grady et al. 2010). Horizontal bars indicate observed intrinsic X-ray variability. Brightness temperatures of $^{13}$CO 2–1 at 100 au are interpolated directly from the $T_B$ radial profiles in Figure 7, so they sometimes differ from the fitted $T_{100}$ in Table 2. If we instead adopted the fitted $T_{100}$ values, the conclusions are unchanged. Peak CO column densities are from Zhang et al. (2021) (and see Figure 14), while the CO 2–1 disk sizes are from Law et al. (2021).</figcaption></figure>

density, decreases in total H₂ surface density but with constant CO abundance, or true geometrical features, e.g., warps. Here, we focus on the first two explanations and note that changes in CO abundance suggest a chemical origin, while decreases in H₂ hint at dynamical, planet-based origins.

The chemical explanation requires CO gas to be sufficiently depleted at the locations of vertical substructures. Alignments of vertical substructures with CO column density gaps, while suggestive, are not conclusive proof of a chemical origin, as gas surface density perturbations and chemical processing are often degenerate in models (e.g., Alarcón et al. 2021). However, if these features have their origins in chemistry, the depletion of gas-phase CO should lead to higher C/O ratios, causing an increase in the column density of molecules such as C₂H (Bergin et al. 2016; Alarcón et al. 2021). To test this possibility, we compare the column density profiles of CO (Zhang et al. 2021) and C₂H (Guzmán et al. 2021) with the identified vertical substructures in Figure 14. The vertical substructures Z56 in AS 209, Z46 in HD 163296, and Z66 in MWC 480 are all associated with CO column density depletions and C₂H enhancements. In contrast, Z81 and Z145 in HD 163296 and Z170 and Z375 in IM Lup are not. This

suggests that chemical conversion of CO into other species may provide at best a partial explanation for the observed vertical substructures.

The dynamical explanation, i.e., if vertical substructures are caused by forming planets, instead requires drops in total H₂ gas surface density. Perhaps, in this case, we expect vertical substructures to be associated with gaps in the millimeter-sized grains. As large grains should be concentrated in gas pressure maxima, this means that dust gaps will correspond to pressure minima and be associated with drops in the total H₂ gas surface density. As discussed in Section 4.4.3, all of the MAPS disks show some degree of spatial association between millimeter continuum gaps and vertical substructures. Below, we consider the plausibility of this interpretation for each MAPS disk.

In AS 209, the Z56 substructure is radially coincident with a deep gap in H₂ (Teague et al. 2018b). Moreover, Fedele et al. (2018) showed that either a single Saturn-mass planet at 95 au or a second ($<0.1 M_{Jup}$) planet at 57 au reproduced the observed continuum profile. Thus, a planetary origin for Z56 (and the larger wave-like structure) in AS 209 is possible, but see Alarcón et al. (2021), who constrain the mass of a putative planet at 100 au to be $<0.2 M_{Jup}$.

In HD 163296, the derived CO surfaces, and in particular the measured depths of vertical substructures, most closely match the models of Rab et al. (2020) that include deep gas gaps, i.e., similar depletion to that for the dust, at the locations of the observed millimeter continuum gaps. Similarly, the gas gap models (vs. that of CO depletion) from Calahan et al. (2021) are better able to reproduce the Z81 dip in the C¹⁸O emission surface. As HD 163296 is also believed to host three Jupiter-mass planets (Pinte et al. 2018; Teague et al. 2018a), this offers a plausible explanation for the vertical substructures in this disk.

In IM Lup, GM Aur, and MWC 480, the plausibility of vertical substructures having their origins in planets is less clear. In IM Lup, Pinte et al. (2020) reported a tentative localized deviation from Keplerian rotation at 117 au, which is thought to be due to a planetary perturber. Intriguingly, this is at the same radius where we observe a change in the slope of the CO 2–1 emitting surface. GM Aur has been suggested to have a 0.1–0.4 $M_{Jup}$ planet at 67 au based on the width of the nearby dust gap (Huang et al. 2020), but no vertical substructures are observed near this radius. The wave-like feature in MWC 480 shows broad associations with dust gaps, and Teague et al. (2021) propose a planet at 245 au, which is driving the wave-like perturbations.

Regardless of the specific mechanisms responsible for these vertical substructures, the MAPS data suggest that emitting surfaces are far from smooth. The locations, depths, and widths of such features provide important inputs to disk thermochemical models and serve as powerful probes of the planet formation. As such, they may also offer another promising mechanism to infer the existence of embedded, newly forming planets in disks.

## 6. Conclusions

We present a detailed analysis of the vertical distribution of molecules and their emitting surfaces in high-angular-resolution observations in five protoplanetary disks from MAPS. We conclude the following:
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Column density profiles for CO (Zhang et al. 2021) and C<sub>2</sub>H (Guzmán et al. 2021) for all MAPS sources vs. locations of vertical substructures in CO, <sup>13</sup>CO, and C<sup>18</sup>O 2–1 emission lines. Orange dashed lines indicate vertical dips listed in Table 4. The CO snowlines from Zhang et al. (2021) are shaded in blue. Red solid lines indicate the radial locations of changes in emission surface slope. The low-amplitude, wave-like features in AS 209 and MWC 480 are shown in black, with dashed and solid lines marking troughs and peaks, respectively.</figcaption></figure>

1. CO emission traces the most elevated regions $z/r > 0.3$, while the less abundant <sup>13</sup>CO and C<sup>18</sup>O probe deeper into the disk $z/r \approx 0.1–0.2$. These heights correspond to approximately 3 and 1 scale height, respectively.

2. In the disks around the T Tauri star AS 209 and Herbig Ae stars HD 163296 and MWC 480, C<sub>2</sub>H and HCN emission heights are also measurable, and they emit from $z/r \lesssim 0.1$, a region relatively close to the planet-forming disk midplane.

3. The NIR surfaces, which trace micron-sized dust, of HD 163296 and IM Lup are lower than the CO 2–1 emission surface and lie at or slightly above that of <sup>13</sup>CO 2–1.

4. We derive radial temperature distributions for all CO isotopologues and use them to estimate full 2D, (r, z) empirical temperature models for each disk.

5. Emission surfaces present substructures in the form of vertical dips, often seen in more than one CO isotopologue, and are detected in a majority of MAPS disks.

6. The wide range of vertical emission heights across the sample indicates a diversity in thermal, density, or CO abundance structures. Tentative trends suggest that star

+disk systems with lower stellar masses and luminosities, as well as larger CO disk sizes, exhibit the most elevated CO line-emitting surfaces. However, a larger sample of disks with well-constrained disk emitting layers is required to better understand what sets emitting layer heights in disks.

7. At least some, and possibly the majority, of vertical disk substructures have their origins in local H<sub>2</sub> surface density drops owing to embedded planets. Others may have their origin in chemical effects, namely, local reductions in CO abundance and thus CO optical depth.

Overall, we have shown an effective method for extracting the emitting layers for a sample of disks and emission lines. As disks are highly structured both radially and vertically, emission surfaces in a set of lines with varying optical depths, e.g., CO isotopologues, provide direct observational constraints on the overall 2D disk structure. Moreover, these surfaces serve as critical inputs to thermochemical models of disks, which are necessary not only to understand the true origins of vertical gas structures but also to connect observed molecular emission to midplane abundances, and therefore the chemical environment within which planets form.
## 7. Value-added Data Products

The MAPS VADPs described in this work can be accessed through the ALMA Archive via https://almascience.nrao.edu/alma-data/lp/maps. An interactive browser for this repository is also available on the MAPS project homepage at http://www.alma-maps.info.

For each combination of data processing (individual measurements, radially binned, and moving average), the following data products are available:

1. Emission surfaces.

2. Gas temperature structures, radial and full 2D (r, z) profiles.

3. Python script to generate the data products.

Each of these VADPs is provided for CO 2–1, $^{13}$CO 2–1, and, when available, C$^{18}$O 2–1 in all MAPS disks, and for C$_2$H 3–2 and HCN 3–2 in HD 163296. The naming scheme for these VADPs is as follows: [disk]_[line]_[frequency]_[resolution]_[datatype], where datatype is “individual measurements,” “radially binned,” or “moving average.” Additional data products associated with the MAPS Large Program, including line image cubes (see Section 9, Czekala et al. 2021) and radial profiles and moment maps (see Section 7, Law et al. 2021), are also available.

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

C.J.L. acknowledges funding from the National Science Foundation Graduate Research Fellowship under grant DGE1745303. R.T. and F.L. acknowledge support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow. J. B. acknowledges support by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51427.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. K.I.Ö. acknowledges support from the Simons Foundation (SCOL No. 321183) and an NSF AAG grant (No. 1907653). I.C. was supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51405.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. S.M. A. and J.H. acknowledge funding support from the National Aeronautics and Space Administration under grant No. 17-XRP17 2-0012 issued through the Exoplanets Research Program. V.V.G. acknowledges support from FONDECYT Iniciación 11180904 and ANID project Basal AFB-170002. J.H. acknowledges support for this work 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. Y.A. acknowledges support by NAOJ ALMA Scientific Research grant code 2019-13B, and

Grant-in-Aid for Scientific Research Nos. 18H05222 and 20H05847. F.A. acknowledges support from NSF AAG grant No. 1907653. E.A.B. and A.D.B. acknowledge support from an NSF AAG grant (No. 1907653). J.B.B. acknowledges support from NASA through the NASA Hubble Fellowship grant No. HST-HF2-51429.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. A.S.B. acknowledges the studentship funded by the Science and Technology Facilities Council of the United Kingdom (STFC). J.K.C. acknowledges support from the National Science Foundation Graduate Research Fellowship under grant No. DGE 1256260 and the National Aeronautics and Space Administration FINESST grant, under grant No. 80NSSC19K1534. G.C. is supported by NAOJ ALMA Scientific Research grant code 2019-13B. L.I.C. gratefully acknowledges support from the David and Lucile Packard Foundation and Johnson & Johnson's WiSTEM2D Program. J.D.I. acknowledges support from the STFC under ST/T000287/1. R.L.G. acknowledges support from a CNES fellowship grant. Y.L. acknowledges the financial support by the Natural Science Foundation of China (grant No. 11973090). F.M. acknowledges support from ANR of France under contract ANR-16-CE31-0013 (Planet-Forming-Disks) and ANR-15-IDEX-02 (through CDP “Origins of Life”). H.N. acknowledges support by NAOJ ALMA Scientific Research grant code 2018-10B and Grant-in-Aid for Scientific Research No. 18H05441. L. M.P. acknowledges support from ANID project Basal AFB-170002 and from ANID FONDECYT Iniciación project No. 11181068. K.R.S. acknowledges the support of NASA through Hubble Fellowship Program grant HST-HF2-51419.001, 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. T.T. is supported by JSPS KAKENHI grant Nos. JP17K14244 and JP20K04017. Y.Y. is supported by IGPEES, WINGS Program, the University of Tokyo. M.L.R.H. acknowledges support from the Michigan Society of Fellows. C.W. acknowledges financial support from the University of Leeds, STFC, and UKRI (grant Nos. ST/R000549/1, ST/T000287/1, MR/T040726/1). K.Z. acknowledges the support of the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin–Madison, with funding from the Wisconsin Alumni Research Foundation, and the support of NASA through Hubble Fellowship grant HST-HF2-51401.001, 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.

**Facility:** ALMA.

**Software:** Astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), bettermoments (Teague & Foreman-Mackey 2018), CASA (McMullin et al. 2007), disksurf (https://github.com/richteague/disksurf), emcee (Foreman-Mackey et al. 2013), GoFish (Teague 2019), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), SciPy (Virtanen et al. 2020).

## Appendix A
### Surface Filtering and Extraction

In this appendix we describe in detail the functionality of disksurf and how we used it to extract emission surfaces from line image cubes:

Before extracting surfaces, we applied the following two data filtering steps. First, a radially varying clip was used to
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Clipping based on disk priors and flux (left) and filtered emission surfaces (right) for CO 2–1 in IM Lup. The gray box in the left panel shows the field of view in the right panel. Light-gray points are filtered according to the steps described in detail in Appendix A.</figcaption></figure>

remove pixels that are not related to the peak of the line. This was done by calculating an azimuthally averaged profile of the peak surface brightness and then clipping values that are more than $1\sigma$ away. This clip threshold was increased to $2\sigma$ for lower-S/N lines, e.g., C$_2$H 3–2 and HCN 3–2. Then, we performed a 1D smoothing$^{30}$ to better define the peaks using a Gaussian kernel with an FWHM equal to half of the beam major-axis FWHM. These two steps were found to significantly improve the ability of the code to identify the emission peaks and minimize contamination from background thermal noise.

We used the `get_emissionSurface` function to extract the deprojected radius *r*, emission height *z*, surface brightness $I_\nu$, and channel velocity *v* for each pixel associated with the emitting surface. This requires knowledge of the disk inclination and position angle in order to correctly account for the deprojection. We adopted these values from Table 1 in Öberg et al. (2021), with disk inclinations ranging from 37°0 (AS 209) to 53°2 (GM Aur).

We then applied additional clipping based on our priors of disk physical structure. In particular, we removed extremely high $z/r$ values$^{31}$ and large negative *z* values, as the emission

must arise from at least the midplane. We allowed points with a small negative *z* value, i.e., $z/r > -0.1$, to remain to avoid positively biasing our averages to nonzero *z* values. We also filtered those points with low surface brightness (less than 5 times the image cube rms) to ensure that noise did not significantly bias the derived surfaces. However, due to the lower S/N of the C$_2$H and HCN lines, we did not perform this clipping in any disk except for HD 163296. Figure 15 shows an example of this process for CO 2–1 in IM Lup.

## Appendix B

### Full List of Isovelocity Contours

Isovelocity contours for the CO isotopologues are shown for the IM Lup disk in Figure 16. A full set of isovelocity contours for all MAPS disks are shown in the figure set, which is available online. All isovelocity contours are calculated using updated dynamical masses taken from Öberg et al. (2021), which are based on CO rotation profiles (Teague et al. 2021).

<footnote>$^{30}$ We note that this smoothing is performed prior to and as a part of the extraction process and is thus distinct from the radial binning of the extracted data used subsequently to increase the S/N.</footnote>

<footnote>$^{31}$ This threshold was initially set to clip vertical heights exceeding $z/r = 0.5$ for all disks but had to be increased to $z/r = 0.7$ for CO 2–1 in IM Lup and GM Aur, due to their very elevated surfaces.</footnote>
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Isovelocity contours of the best-fitting model, as indicated in Table 1, for the CO 2–1 (top) and $^{13}$CO 2–1 (bottom) emitting surfaces in IM Lup plotted for selected channels of the observed emission. Solid curves indicate the upper surface of the disk, and dashed curves mark the lower surface. LSRK velocities are noted in the lower right corner. The synthesized beam is shown in the lower left corner of each panel. (The complete figure set (5 images) showing isovelocity contours for each MAPS disk is available in the online journal.)</figcaption></figure>

(The complete figure set (5 images) is available.)

## Appendix C
### Excitation and Band 3 CO Surfaces

Since we also had access to $^{13}$CO 1–0, we were able to compare against the $^{13}$CO 2–1 line to see whether we could identify any excitation-related effects in the emission surfaces, i.e., differing heights (e.g., van Zadelhoff et al. 2001; Dartois et al. 2003). Due to the coarser spatial resolution and lower S/N of the 1–0 line, we did not attempt to extract the emission surfaces directly. Instead, we compared the $^{13}$CO 2–1 isovelocity contours derived from the parametric fit in Table 1 with the spatial distribution of the 1–0 line. To ensure a consistent comparison, we also included the tapered (0/30) resolution $^{13}$CO 2–1 images. We checked C$^{18}$O 1–0, which was also covered by the MAPS observations, but it did not possess sufficient S/N for this comparison, so we instead focused on $^{13}$CO 1–0. At this lower

resolution, only GM Aur, HD 163296, and IM Lup had suffi- ciently elevated $^{13}$CO 2–1 surfaces to allow for a meaningful comparison.

Figure 17 shows isovelocity contours overlaid on a representative channel of $^{13}$CO 1–0 emission that should best show the emitting layers, if resolved. In IM Lup and GM Aur, line emission tracing the backside of the disks is visible in the tapered 2–1 image, but we cannot determine whether the contours are consistent with low-S/N emission surfaces in $^{13}$CO 1–0, or if the 1–0 line is truly flatter than the 2–1 line. In the case of HD 163296, the spatial resolution is insufficient to reveal any vertical disk structure in either the $^{13}$CO 2–1 tapered or 1–0 images, and we therefore cannot compare 2–1 and 1–0 emission layer heights. Thus, in order to infer the emitting heights of $^{13}$CO 1–0, we likely require both higher spatial resolution and higher S/N.
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Representative channels for IM Lup, GM Aur, and HD 163296 for the ¹³CO 2–1 full-resolution (left column); ¹³CO 2–1, tapered to 0.′′30 (middle column); and ¹³CO 1–0 (right column) images. The ¹³CO 2–1 isovelocity contours derived using the parametric fit in Table 1 are shown in pink. Solid curves indicate the upper surface of the disk, and dashed curves mark the lower surface. LSRK velocities are noted in the upper right corner. The synthesized beam and a scale bar indicating 50 au are shown in the lower left and right corners, respectively, of each panel.</figcaption></figure>

# Appendix D

## Effects of Spatial Resolution on Derived Emission Surfaces

To extract emission surfaces using *disksurf*, the image cubes must have some minimum angular resolution, i.e., the front and back sides of the disk must be sufficiently spatially resolved to be separable. This means that emission surfaces will be sensitive to the spatial resolution of the images used to derive them. To investigate the effects of spatial resolution on our surfaces, we repeated the surface extraction for CO, ¹³CO, and C¹⁸O 2–1 using all angular resolutions, i.e., 0.′′3, 0.′′2, 0.′′15, imaged as part of MAPS (Öberg et al. 2021). We then compared them to the surfaces derived from the images generated with a robust parameter of 0.5 used throughout this work. Figures 18, 19, and 20 show the resulting surfaces.

The CO isotopologue surfaces are generally consistent across differing spatial resolutions. The lower resolutions (0.′′3, 0.′′2) occasionally underestimate the average *z/r* surface height, e.g., ¹³CO and C¹⁸O in HD 163296, and this effect is

more conspicuous in intrinsically flatter surfaces, e.g., ¹³CO and C¹⁸O in MWC 480 or those from rarer isotopologues. In highly elevated surfaces like ¹²CO 2–1, the emission structure along a given column of pixels will be two well-separated Gaussians (and another two Gaussians for the backside of the disk). At lower spatial resolutions, these Gaussians are broadened but do not overlap. Conversely, for ¹³CO or other emission lines with lower intrinsic *z/r*, these two components may overlap and thus lead to a single-peaked Gaussian with a *z/r* that approaches 0. Even for those highly flared surfaces, images with higher spatial resolutions are preferable, as they allow for the detection and characterization of vertical substructures, such as the dips, wave-like features, and slope changes seen in many of the MAPS disks (see Section 4.4). For instance, the ¹²CO surface in HD 163296 appears nearly identical between the 0.′′3 and robust = 0.5 images, with the important exception of the Z46 dip, which can only be identified in the 0.′′15 and robust = 0.5 images.
<figure><img src="image_21.png" /><figcaption>**Figure 18.** $^{12}$CO 2–1 emission surfaces derived using images with different angular resolutions (see Table 5 in Öberg et al. 2021). Large gray points show a consistent radial binning of 1/2× the beam major axis of each image, while small light-gray points represent individual measurements. Lines of constant z/r from 0.1 to 0.5 in increments of 0.1 are shown in gray. The FWHM of the major axis of the synthesized beam is shown in the upper left corner of each panel.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 19.** $^{13}$CO 2–1 emission surfaces derived using different spatial resolutions. Otherwise, same as in Figure 18.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** C$^{18}$O 2–1 emission surfaces derived using different spatial resolutions for HD 163296 and MWC 480. The C$^{18}$O 2–1 surface in GM Aur is not shown, as it was only able to be extracted from the images generated with a robust parameter of 0.5, as only this image possessed a sufficiently high spatial resolution necessary to separate the front and back disk surfaces. Otherwise, same as in Figure 18.</figcaption></figure>

## ORCID iDs

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

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

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

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

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

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

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

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

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

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

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

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

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

Jenny K. Calahan https://orcid.org/0000-0002-0150-0125

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

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

Kenji Furuya https://orcid.org/0000-0002-2026-8157

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

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

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

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
Hideko Nomura https://orcid.org/0000-0002-7058-7682

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

Chunhua Qi https://orcid.org/0000-0001-8642-1786

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

Daniela Soto https://orcid.org/0000-0001-8455-7954

Takashi Tsukagoshi https://orcid.org/0000-0002-6034-2892

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

Merel L. R. van 't Hoff https://orcid.org/0000-0002-2555-9869

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

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

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

## References

Alarcón, F., Bosman, A. D., Bergin, E. A., et al. 2021, ApJS, 257, 8

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., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129

Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46

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

Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44

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

Bosman, A., Bergin, E. A., Loomis, R. A., et al. 2021, ApJS, 257, 15

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

Cleeves, L. I. 2016, ApJL, 816, L21

Cleeves, L. I., Bergin, E. A., Öberg, K. I., et al. 2017, ApJL, 843, L3

Cleeves, L. I., Loomis, R. A., Teague, R., et al. 2021, ApJ, 911, 29

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

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

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

de Gregorio-Monsalvo, I., Ménard, F., Dent, W., et al. 2013, A&A, 557, A133

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

Dullemond, C. P., & Dominik, C. 2004a, A&A, 417, 159

Dullemond, C. P., & Dominik, C. 2004b, A&A, 421, 1075

Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957

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

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

Espaillat, C. C., Robinson, C., Grant, S., & Reynolds, M. 2019, ApJ, 876, 121

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

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

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

Fogel, J. K. J., Bethell, T. J., Bergin, E. A., Calvet, N., & Semenov, D. 2011, ApJ, 726, 29

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

Goodman, J., & Weare, J. 2010, CAMCS, 5, 65

Grady, C. A., Hamaguchi, K., Schneider, G., et al. 2010, ApJ, 719, 1565

Günther, H. M., & Schmitt, J. H. M. M. 2009, A&A, 494, 1041

Guzmán, V. V., Bergner, J. B., Law, C. J., et al. 2021, ApJS, 257, 6

Guzmán, V. V., Öberg, K. I., Carpenter, J., et al. 2018, ApJ, 864, 170

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

Ilgner, M., Henning, T., Markwick, A. J., & Millar, T. J. 2004, A&A, 415, 643

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

Kenyon, S. J., & Hartmann, L. 1987, ApJ, 323, 714

Kusakabe, N., Grady, C. A., Sitko, M. L., et al. 2012, ApJ, 753, 153

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

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

Monnier, J. D., Harries, T. J., Aarnio, A., et al. 2017, ApJ, 838, 20

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

Muro-Arena, G. A., Dominik, C., Waters, L. B. F. M., et al. 2018, A&A, 614, A24

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

Ö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

Paneque-Carreño, T., Pérez, L. M., Benisty, M., et al. 2021, ApJ, 914, 88

Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125

Piétu, V., Dutrey, A., & Guilloteau, S. 2007, A&A, 467, 163

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

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

Podio, L., Garufi, A., Codella, C., et al. 2020, A&A, 642, L7

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

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

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

Rich, E. A., Wisniewski, J. P., Sitko, M. L., et al. 2020, ApJ, 902, 4

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

Ruíz-Rodríguez, D., Kastner, J., Hily-Blant, P., & Forveille, T. 2021, A&A, 646, A59

Schneider, G., Wood, K., Silverstone, M. D., et al. 2003, AJ, 125, 1467

Semenov, D., & Wiebe, D. 2011, ApJS, 196, 25

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

Teague, R. 2019, JOSS, 4, 1632

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

Teague, R., Bae, J., & Bergin, E. A. 2019, 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., & Foreman-Mackey, D. 2018, Bettermoments: A Robust Method To Measure Line Centroids, v1.0, Zenodo, doi: 10.5281/zenodo.1419754

Teague, R., Jankovic, M. R., Haworth, T. J., Qi, C., & Ilee, J. D. 2020, MNRAS, 495, 451

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

van Zadelhoff, G. J., van Dishoeck, E. F., Thi, W. F., & Blake, G. A. 2001, A&A, 377, 566

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

Walsh, C., Millar, T. J., & Nomura, H. 2010, ApJ, 722, 1607

Walter, F. M., & Kuhi, L. V. 1981, ApJ, 250, 254

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

Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67

Woitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383

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

Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJL, 877, L18