# Molecules with ALMA at Planet-forming Scales (MAPS). V. CO Gas Distributions

Ke Zhang, Alice S. Booth, Charles J. Law, Arthur D. Bosman, Kamber R. Schwarz, Edwin A. Bergin, Karin I. Öberg, Sean M. Andrews, Viviana V. Guzmán, Catherine Walsh, et al.

► To cite this version:

Ke Zhang, Alice S. Booth, Charles J. Law, Arthur D. Bosman, Kamber R. Schwarz, et al. Molecules with ALMA at Planet-forming Scales (MAPS). V. CO Gas Distributions. The Astrophysical Journal Supplement Series, 2021, 257, 10.3847/1538-4365/ac1580. insu-03672376

HAL Id: insu-03672376

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

Submitted on 29 Jul 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Molecules with ALMA at Planet-forming Scales (MAPS) V: CO gas distributions

KE ZHANG<sup>1,2,*</sup> ALICE S. BOOTH<sup>3,4</sup> CHARLES J. LAW<sup>5</sup> ARTHUR D. BOSMAN<sup>1,2</sup> KAMBER R. SCHWARZ<sup>6,†</sup>

EDWIN A. BERGIN<sup>1,2</sup> KARIN I. ÖBERG<sup>1,5</sup> SEAN M. ANDREWS<sup>1,5</sup> VIVIANA V. GUZMÁN<sup>1,7</sup> CATHERINE WALSH<sup>1,4</sup>

CHUNHUA QI<sup>1,5</sup> MEREL L. R. VAN 'T HOFF<sup>1,2</sup> FENG LONG<sup>1,5</sup> DAVID J. WILNER<sup>1,5</sup> JANE HUANG<sup>1,2,5,†</sup>

IAN CZEKALA<sup>1,8,9,10,11,12,†</sup> JOHN D. ILEE<sup>1,13</sup> GIANNI CATALDI<sup>1,14,15</sup> JENNIFER B. BERGNER<sup>1,16,†</sup> YURI AIKAWA<sup>1,17</sup>

RICHARD TEAGUE<sup>1,5</sup> JAEHAN BAE<sup>1,18,19,†</sup> RYAN A. LOOMIS<sup>1,20</sup> JENNY K. CALAHAN<sup>1,2</sup> FELIPE ALARCÓN<sup>1,2</sup>

FRANÇOIS MÉNARD<sup>1,21</sup> ROMANE LE GAL<sup>1,5,22,23,24</sup> ANIBAL SIERRA<sup>1,25</sup> YOSHIHIDE YAMATO<sup>1,15</sup> HIDEKO NOMURA<sup>1,14</sup>

TAKASHI TSUKAGOSHI<sup>1,14</sup> LAURA M. PÉREZ<sup>1,25</sup> LEON TRAPMAN<sup>1,1</sup> YAO LIU<sup>1,26</sup> AND KENJI FURUYA<sup>1,14</sup>

<sup>1</sup>Department of Astronomy, University of Wisconsin-Madison, 475 N Charter St, Madison, WI 53706

<sup>2</sup>Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, USA

<sup>3</sup>Leiden Observatory, Leiden University, 2300 RA Leiden, the Netherlands

<sup>4</sup>School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK

<sup>5</sup>Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

<sup>6</sup>Lunar and Planetary Laboratory, University of Arizona, 1629 E. University Blvd, Tucson, AZ 85721, USA

<sup>7</sup>Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile

<sup>8</sup>Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

<sup>9</sup>Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

<sup>10</sup>Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

<sup>11</sup>Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA

<sup>12</sup>Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA

<sup>13</sup>School of Physics and Astronomy, University of Leeds, Leeds, UK, LS2 9JT

<sup>14</sup>National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

<sup>15</sup>Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

<sup>16</sup>University of Chicago Department of the Geophysical Sciences, Chicago, IL 60637, USA

<sup>17</sup>Department of Astronomy, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan

<sup>18</sup>Department of Astronomy, University of Florida, Gainesville, FL 32611, USA

<sup>19</sup>Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA

<sup>20</sup>National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA

<sup>21</sup>Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

<sup>22</sup>IRAP, Université de Toulouse, CNRS, CNES, UT3, 31400 Toulouse, France

<sup>23</sup>IPAG, Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

<sup>24</sup>IRAM, 300 rue de la piscine, F-38406 Saint-Martin d'Hères, France

<sup>25</sup>Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile

<sup>26</sup>Purple Mountain Observatory & Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, Nanjing 210023, China

ABSTRACT

Here we present high resolution (15-24 au) observations of CO isotopologue lines from the Molecules with ALMA on Planet-forming Scales (MAPS) ALMA Large Program. Our analysis employs <sup>13</sup>CO and C<sup>18</sup>O (J=2-1), (1-0), and C<sup>17</sup>O (1-0) line observations of five protoplanetary disks. We retrieve CO gas density distributions, using three independent methods: (1) a thermo-chemical modeling framework based on the CO data, the broadband spectral energy distribution, and the mm-continuum emission; (2) an empirical temperature distribution based on optically thick CO lines; and (3) a direct fit to the C<sup>17</sup>O hyperfine lines. Results from these methods generally show excellent agreement. The CO gas column density profiles of the five disks show significant variations in the absolute value and the radial
shape. Assuming a gas-to-dust mass ratio of 100, all five disks have a global CO-to-H₂ abundance of 10-100 times lower than the ISM ratio. The CO gas distributions between 150-400 au match well with models of viscous disks, supporting the long-standing theory. CO gas gaps appear to be correlated with continuum gap locations, but some deep continuum gaps do not have corresponding CO gaps. The relative depths of CO and dust gaps are generally consistent with predictions of planet-disk interactions, but some CO gaps are 5-10 times shallower than predictions based on dust gaps. This paper is part of the MAPS special issue of the Astrophysical Journal Supplement.

**Keywords:** Astrochemistry — Protoplanetary disks— Exoplanet formation — Interferometry — Millimeter astronomy

# 1. INTRODUCTION

In protoplanetary disks, H₂ is the primary mass carrier, but it does not have strong transitions in the temperature range that most of the disk mass resides. Instead, CO is the most widely used gas tracer, as it is abundant and has rotational transitions sensitive to cold, 10-50 K, gas (e.g., Bergin & Williams 2017). CO is also chemically stable (van Dishoeck & Black 1988; Visser et al. 2009), and stays in the gas phase at temperature warmer than ~20 K (Molyarova et al. 2017). Besides being a tracer of gas mass, CO is one of the primary carbon and oxygen carriers in protoplanetary disks, a stepping-stone to a series of chemical reactions that eventually lead to the formation of complex organics (e.g., Walsh et al. 2014).

Observations of (sub)mm CO ($J \le 3$) lines have been widely used to study fundamental physical properties of protoplanetary disks, including total disk gas mass (Williams & Best 2014; Ansdell et al. 2016; Miotello et al. 2017; Long et al. 2017), gas depletion in cavities and gaps (van der Marel et al. 2015; Facchini et al. 2018; Favre et al. 2019), temperature structures (Qi et al. 2006; Qi et al. 2011; Piétu et al. 2007; Schwarz et al. 2016; Kama et al. 2016; Pinte et al. 2018a), kinematics in disks (Teague et al. 2018a; Pinte et al. 2018b; Teague et al. 2019), and turbulence levels (Teague et al. 2016; Flaherty et al. 2015, 2017, 2020). Therefore, accurately measuring the CO gas distribution in disks is essential for our understanding of the physical and chemical evolution of protoplanetary disks as well as the planet formation processes occurring within them.

Over the past few years, one of the most extensive debates is whether CO is a robust tracer of H₂ gas mass in disks. The gas mass derived from CO observations hinges on the crucial assumption that the CO-to-H₂ ratio is the canonical ISM ratio of $10^{-4}$ in the warm molecular layer throughout the entire disk life-time (Bergin &

Williams 2017). However, observations are in tension with this assumption. Surveys of protoplanetary disks in nearby star-forming regions show that gas masses derived from CO lead to a gas-to-dust ratio $\le 10$ in the majority of disks (Ansdell et al. 2016; Miotello et al. 2017; Long et al. 2017). In contrast, the high gas accretion rates onto the central star suggest many disks are still gas-rich (Manara et al. 2016). The most direct evidence of CO depletion comes from independent gas mass measurements with the HD $J=1-0$ line in three disks (Bergin et al. 2013; McClure et al. 2016; Trapman et al. 2017). Their CO-to-H₂ ratios are one to two orders of magnitude lower than the canonical ISM ratio (Favre et al. 2013; Schwarz et al. 2016; Zhang et al. 2017, 2019; Calahan et al. 2021a).

To explain low CO abundances seen in protoplanetary disks, two broad types of processes have been proposed. The first type is chemical processes — CO is chemically processed into other molecules (e.g., Aikawa et al. 1998; Bergin et al. 2014; Schwarz et al. 2018; Bosman et al. 2018). The second type is physical processes — CO gas is gradually removed from the disk atmosphere, as icy dust grains grow and settle towards the mid-plane (Xu et al. 2017; Krijt et al. 2016, 2018, 2020). Recent comparisons between the CO gas mass in Class I and II disks suggest that the CO-to-H₂ abundance ratio decreases rapidly on a timescale of $\sim 1$ Myr (Zhang et al. 2020b; Bergner et al. 2020). In short, current observations and theoretical studies suggest that the CO gas abundance evolves along with the disk evolution and planet formation, rather than a simple heritage of the molecular cloud.

Another important question is whether the CO gas distribution has substructures at dust gap locations. Substructures (e.g., rings, gaps, and spirals) have been commonly seen in (sub)mm continuum emission of protoplanetary disks (e.g., Zhang et al. 2016; Isella et al. 2016; Andrews et al. 2018; Long et al. 2018; Cieza et al. 2021). Various mechanisms have been proposed as possible causes, including planet-disk interaction (Dong et al. 2015; Rosotti et al. 2016; Zhang et al. 2018), dust prop-

<footnote>* NASA Hubble Fellow</footnote>

<footnote>† NASA Hubble Fellowship Program Sagan Fellow</footnote>
erty changes at snowlines (Zhang et al. 2015; Okuzumi et al. 2016; Pinilla et al. 2017), and magnetic instability (Flock et al. 2015; Suriano et al. 2018; Riols & Lesur 2019). As gas and mm-sized dust grains have different dynamic behaviors, gas substructures can provide an important complementary insight into possible causes. In particular, planet-disk interaction models predict that giant planets can open deep gas gaps that manifest as gaps in CO line emissions (e.g., Facchini et al. 2018; van der Marel et al. 2018; Alarcón et al. 2020).

Most of previous studies of the CO gas mass distri-
butions in Class II disks have relatively low spatial res-
olutions, between 25-100 au (e.g., Schwarz et al. 2016;
Nomura et al. 2016; Isella et al. 2016; Fedele et al. 2017;
Zhang et al. 2019; Favre et al. 2019; Rosotti et al. 2021).
Although local coincidences in CO and dust gap loca-
tions have been reported in a few disks, previous obser-
vations often did not have sufficient spatial resolution to
fully characterize CO substructures. Also, studies often
focused on local perturbations and thus lacked a holistic
view of the CO gas distribution across the entire disk.

Here we study the CO gas mass distribution in five
protoplanetary disks observed as part of the Molecules
with ALMA on Planet-forming Scales (MAPS) ALMA
Large Program. We use five CO isotopologue lines to
constrain CO gas distribution at a scale of ~15-24 au.
We aim to robustly measure CO gas distributions in
these five disks, focusing on three important questions:
(1) How does the CO column density vary with ra-
dius? (2) How much does the CO column density change
across continuum substructures? (3) What are the lower
limits of gas masses in these disks?

This paper is structured as follows. In Section 2, we
briefly present the observations of CO lines in the MAPS
program and the data reduction process. We describe
the three methods used to constrain the CO column den-
sity distributions and present the results in Section 3. In
Section 4, we discuss how these results answer the three
questions above. Then we summarize our findings in
Section 5.

## 2. OBSERVATIONS

This study uses CO isotopologue line data covered in
the MAPS ALMA Large Program (2018.1.01055.L) for
five disk sources: IM Lup, GM Aur, AS 209, HD 163296,
and MWC 480. The ALMA Band 6 setting targeted
the CO, ¹³CO, C¹⁸O *J*=2-1 transitions and the Band
3 setting targeted the ¹³CO, C¹⁸O, C¹⁷O *J*=1-0 transi-
tions (Oberg et al. 2021). The analysis presented here
is based on the MAPS fiducial images, which have 0.′′15
and 0.′′30 circularized beams for lines in Bands 6 and 3,
respectively. The full details of the data calibration and

imaging are described in Czekala et al. (2021). The ob-
servational logs, line frequencies, and image statistics are
listed in Oberg et al. (2021). We use the integrated in-
tensity maps and radial emission profiles from Law et al.
(2021a). Figures 1 and 2 show integrated line intensity
maps and radial profiles used in this study.

We also present the C¹⁷O (1-0) data which are not
covered by Law et al. (2021a). The C¹⁷O (1-0) line has
three hyperfine components which are blended in the
integrated intensity map. As shown in Figure 1 and
2, the C¹⁷O (1-0) lines were robustly detected in the
HD 163296 and MWC 480 disks, and tentatively in the
IM Lup, GM Aur, and AS 209 disks. For the tentative
detections, we used the matched filter method developed
by Loomis et al. (2018) to estimate the signal-to-noise
ratio of their detections. The C¹⁷O data were searched
with a Keplerian mask extending to 100 au for all five
disks. This mask size is chosen because most of the
C¹⁷O emission is expected to be compact (≤ 100 au)
originating from within the CO snowline (Zhang et al.
2017; Booth et al. 2019). The matched filter responses
for the GM Aur, IM Lup, and AS 209 disks are presented
in the Appendix in Figure 22. The C¹⁷O (1-0) line is
robustly detected (20 σ) in the GM Aur disk, but only
tentatively (3-4 σ) in the AS 209 and IM Lup disks.

3. METHODS AND RESULTS

In this work, our goal is to measure how CO gas col-
umn density varies with radial distance from the central
star in the five MAPS disks (N<sub>CO</sub>(R)), using spatially
resolved CO isotopologue line images.

To test the robustness of the results, we employ three
independent methods to constrain N<sub>CO</sub>(R) profiles:

Method 1: we build thermo-chemical models for indi-
vidual disks to constrain the temperature structures of
gas and dust in these disks. These models are based on
fitting broad band spectral energy distribution (SED),
(sub)mm continuum images, and vertical locations of
CO emission surfaces. We then use the thermal struc-
ture to retrieve the N<sub>CO</sub> at each radius by matching the
observed radial intensity profiles of lines. The workflow
of the modeling processes is shown in Figure 3, and the
details are described in Section 3.1.

Method 2: we replace the gas temperature in thermo-
chemical models with empirically constructed tempera-
ture structures based on optically thick CO lines, and
then redo the N<sub>CO</sub> retrieval. This process is described
in Section 3.2.

Method 3: we use line ratios and absolute intensities
of C¹⁷O (1-0) hyperfine components to constrain opti-
cal depth and excitation temperature. This analysis is
described in Section 3.3.
<figure><img src="image_2.png" /><figcaption>Figure 1. Integrated intensity maps of CO, $^{13}$CO, C$^{18}$O, and C$^{17}$O lines for the MAPS disks. Axes are angular offsets from the disk center with each white tick mark representing a spacing of 2". The synthesized beam is plotted in the lower left corner of each panel. The white scale bar at the lower right corner of each panel represents 50 au scale.</figcaption></figure>

In the following subsections, we describe the detailed setup of each method. We note that all of CO data and thermochemical models used in this paper can be downloaded from the official website of the MAPS program: www.alma-maps.info.

## 3.1. Themo-chemical models

### 3.1.1. Density structure

Our disk model is axisymmetric and consists of three mass components – gas, a small dust population, and a large dust population. We assume that gas and the small

dust population are spatially coupled and thus their densities differ only by a constant factor across the disk. The mass surface distribution of gas and small dust is set as a self-similarly viscous disk (e.g., Lynden-Bell & Pringle 1974; Andrews et al. 2011).

$$ \Sigma(R) = \Sigma_c \left(\frac{R}{R_c}\right)^{-\gamma} \exp\left[-\left(\frac{R}{R_c}\right)^{2-\gamma}\right] \quad (1) $$

where $R_c$ is the characteristic scaling radius and $\gamma$ is the gas surface density exponent.
<figure><img src="image_3.png" /></figure>

**Figure 2.** Deprojected radial intensity profiles of CO lines, for the MAPS sample, ordered by increasing stellar mass from left to right. The radial profiles are azimuthally-averaged, see details in Law et al. (2021a). The only exception is the AS 209 CO (2-1) line, which is derived from a 55° wedge to avoid cloud contamination. Blue shaded regions show the standard error on the mean in the annulus or arc over which the emission was averaged.

The individual beam size is plotted at the upper right corner of each panel. The beam sizes are 0.′15×0.′15 for all (2-1) lines, and 0.′3×0.′3 for (1-0) lines, respectively. A plot of the intensity profiles on a logarithmic scale is shown in the Appendix Figure 23.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** An outline of the modeling processes. Codes used are in light blue color, input data are highlighted in solid blue boxes, and the output of each step is highlighted in red.</figcaption></figure>

<table><caption>**Table 1.** Stellar and Disk Properties</caption><thead><tr><th>Source</th><th>SpT</th><th>d<br>(pc)</th><th>incl<br>(deg)</th><th>PA<br>(deg)</th><th>T<sub>eff</sub><br>(K)</th><th>L<sub>*</sub><br>(L<sub>☉</sub>)</th><th>M<sub>*</sub><br>(M<sub>☉</sub>)</th><th>log<sub>10</sub>(M&#775;)<br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>v<sub>sys</sub><br>(km s<sup>-1</sup>)</th><th>References</th></tr></thead><tbody><tr><td>IM Lup</td><td>K5</td><td>158</td><td>47.5</td><td>144.5</td><td>4266</td><td>2.57</td><td>1.1</td><td>-7.9</td><td>4.5</td><td>1,2,3,4,5</td></tr><tr><td>GM Aur</td><td>K6</td><td>159</td><td>53.2</td><td>57.2</td><td>4350</td><td>1.2</td><td>1.1</td><td>-8.1</td><td>5.6</td><td>1,6,7,8,9</td></tr><tr><td>AS 209</td><td>K5</td><td>121</td><td>35.0</td><td>85.8</td><td>4266</td><td>1.41</td><td>1.2</td><td>-7.3</td><td>4.6</td><td>1,2,5,10,11</td></tr><tr><td>HD 163296</td><td>A1</td><td>101</td><td>46.7</td><td>133.3</td><td>9332</td><td>17.0</td><td>2.0</td><td>-7.4</td><td>5.8</td><td>1,2,5,12,13</td></tr><tr><td>MWC 480</td><td>A5</td><td>162</td><td>37.0</td><td>148.0</td><td>8250</td><td>21.9</td><td>2.1</td><td>-6.9</td><td>5.1</td><td>1,14,15,16,17,18</td></tr></tbody></table>

NOTE—Reproduced from Oberg et al. (2021), where further details about the source characteristics are available.

References are (1) Gaia Collaboration et al. (2018); (2) Huang et al. (2018); (3) Alcalá et al. (2017); (4) Pinte et al. (2018a); (5) Andrews et al. (2018); (6) Huang et al. (2020); (7) Macías et al. (2018); (8) Ingleby et al. (2015); (9) Espaillat et al. (2010); (10) Salyk et al. (2013); (11) Huang et al. (2017); (12) Fairlamb et al. (2015); (13) Teague et al. (2019); (14) Liu et al. (2019); (15) Montesinos et al. (2009); (16) Simon et al. (2019); (17) Piétu et al. (2007); (18) Mendigutía et al. (2013)
For the large dust population (pebbles), its surface
density is set by matching mm continuum images. The
details of the process are described in Section 3.1.6.

The vertical density structure is assumed to be a
Gaussian function characterized by a scale height H(R)
that is a power-law function of radius.

$$
\rho_i(R, Z) = f_i \frac{\Sigma(R)}{\sqrt{2\pi} H_i(R)} \exp\left[-\frac{1}{2} \left(\frac{Z}{H_i(R)}\right)^2\right], \quad (2)
$$

$$
H_i(R) = \chi_i H_{100}(R/100 \text{ au})^\psi \quad (3)
$$

where $f_i$ is the mass fraction of each mass component,
$H_{100}$ is the scale height at 100 au, $\psi$ is a parameter that
characterizes the radial dependence of the scale height.
The large grain population is expected to be more set-
tled compared to the gas and small grains (Nakagawa
et al. 1986; Dullemond & Dominik 2004). This settling
effect is characterized with the parameter $\chi_i$. Here we
fix $\chi=1$ for the gas and the small grain population, and
$\chi=0.2$ for the large grain population (Andrews et al.
2011).

3.1.2. Dust Opacity

We employ the dust size distribution of $n(a) \propto a^{-3.5}$ from the Mathis, Rumpl, & Nordsieck model (Mathis et al. 1977), for both small and large dust populations. The minimum grain size $a_{min}$ is fixed at 0.005 µm. The maximum size $a_{max}$ of the small and large dust populations is set to be 1 µm and 1 mm, respectively. The $a_{max}$ of the large dust population is consistent with multi-wavelength analysis of these five disks. Based on the 1.3 mm and 3 mm continuum images, Sierra et al. (2021) showed that all five MAPS disks have an $a_{max}$ between 1 mm and 3 mm beyond 25 au.

For the composition of the large dust population, we
adopt the recommendation of the DSHARP collabora-
tion (Birnstiel et al. 2018), which is a mixture of water
ice, astronomical silicates, troilite, and refractory organ-
ics. For the small dust population, we adopt a 50%-  
50% mixture of silicates and refractory organics because  
current observations and dust evolution models suggest  
that water ice is destroyed by photodesorption at higher  
disk atmosphere region (Hogerheijde et al. 2011; Krijt  
et al. 2016; Du et al. 2017). The dsharp_opac package  
from Birnstiel et al. (2018) is used to compute the wave-  
length dependence of dust opacity. We show the dust  
opacity profiles in the Appendix Figure 24.

3.1.3. Stellar Spectrum

Young stars accrete disk materials actively and pro-
duce strong X-ray and UV radiation (Hartmann et al.
2016). This high energy radiation is particularly impor-
tant for chemistry in disks (Bergin et al. 2007). Here we

make a composite spectrum for each source by combin-
ing a model stellar photosphere spectrum with UV and
X-ray spectra from observations.

The stellar photospheric spectra are computed by in-
terpolating the standard Nextgen spectral model library
in the (log₁₀ T<sub>eff</sub>, log₁₀ g) space, assuming a solar metal-
licity (Hauschildt et al. 1999). The UV and X-ray spec-
tra of individual sources are taken from the DIANA
database (Dionatos et al. 2019). As UV spectra are
highly variable with time, we use averaged spectra. The
observed UV spectra may suffer from extinction of cir-
cumstellar materials, and therefore do not represent the
UV spectra received by the disk. We apply the extinc-
tion at UV wavelength based on A<sub>V</sub> of individual sources
and the wavelength dependent dust opacity of our small
dust grain population. Some of the sources do not have
short wavelength coverage, and we replace these missing
parts by scaling the UV spectrum of TW Hya (Herczeg
et al. 2002, 2004). The UV spectra are then scaled to
fit the stellar spectra. The final composite spectra are
shown in the Appendix in Figure 24.

3.1.4. Thermo-chemical modeling

We use the thermo-chemical code RAC2D to compute
the CO abundance structure and gas temperature. The
detailed description of RAC2D can be found in Du &
Bergin (2014) and we only briefly describe the code here.

Taking a static density structure of gas and dust as
inputs, RAC2D first uses a Monte Carlo approach to
simulate absorption and scattering events of photons
in the disk. The first step produces a dust thermal
structure and radiation field (from X-ray to centimeter
wavelengths). The cosmic-ray ionization inside the disk
is simulated with an attenuation length of 96 g cm⁻²
(Umebayashi & Nakano 1981; Bergin & Williams 2017).
With the computed radiation field, the code then simul-
taneously solves the chemical evolution and gas ther-
mal structure in the disk. The chemical network has
484 species and 5830 reactions, including the full gas-
phase network from the UMIST 2006 database<fnref n="1" />, disso-
ciation of H₂O and OH by Lyα photons, adsorption and
desorption of species on the dust grain surface either
thermally or induced by cosmic rays and UV photons,
and two-body reactions on the dust grain surface. The
self-shielding effects of H₂O and CO are included in the
models.

The initial chemical abundances are listed in Table 1 of
Du & Bergin (2014). We adopt a relatively low cosmic-
ray rate of 1.36×10⁻¹⁸ s⁻<fnref n="1" /> at the disk surface, because

<footnote n="1">1 http://udfa.ajmarkwick.net</footnote>
<table><caption>**Table 2.** Best-fit disk Model Parameters</caption><thead><tr><th>Source</th><th>M<sub>g</sub><br>(M<sub>☉</sub>)</th><th>Σ<sup>g</sup><sub>c</sub><br>(g cm<sup>-2</sup>)</th><th>M<sub>mm</sub><br>(M<sub>☉</sub>)</th><th>M<sub>μm</sub><br>(M<sub>☉</sub>)</th><th>Σ<sup>μm</sup><sub>c</sub><br>(g cm<sup>-2</sup>)</th><th>R<sup>g</sup><sub>c</sub><br>(au)</th><th>γ<sub>g</sub></th><th>H<sub>100</sub><br>(au)</th><th>ψ</th><th>r<sub>in</sub><br>(au)</th><th>r<sub>out</sub><br>(au)</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th><th>(10)</th><th>(11)</th><th>(12)</th></tr></thead><tbody><tr><td>IM Lup</td><td>0.2</td><td>28.4</td><td>1.97e-03</td><td>2.02e-05</td><td>2.9e-03</td><td>100</td><td>1.0</td><td>10.0</td><td>1.17</td><td>0.20</td><td>1200</td></tr><tr><td>GM Aur</td><td>0.2</td><td>9.4</td><td>5.94e-04</td><td>1.03e-04</td><td>6.0e-03</td><td>176</td><td>1.0</td><td>7.5</td><td>1.35</td><td>1.00</td><td>650</td></tr><tr><td>AS 209</td><td>0.0045</td><td>1.0</td><td>4.50e-04</td><td>5.23e-05</td><td>1.2e-02</td><td>80</td><td>1.0</td><td>6.0</td><td>1.25</td><td>0.50</td><td>300</td></tr><tr><td>HD 163296</td><td>0.14</td><td>8.8</td><td>2.31e-03</td><td>2.00e-04</td><td>1.3e-02</td><td>165</td><td>0.8</td><td>8.4</td><td>1.08</td><td>0.45</td><td>600</td></tr><tr><td>MWC 480</td><td>0.16</td><td>5.8</td><td>1.41e-03</td><td>1.69e-04</td><td>6.1e-03</td><td>200</td><td>1.0</td><td>10.0</td><td>1.08</td><td>0.45</td><td>750</td></tr></tbody></table>

NOTE—(1) Source name; (2) total gas mass of the disk; (3) surface density of gas at $R_c^g$; (4) total mass of the large grain population; (5) total mass of the small grain population; (6) surface density of the small grain population at $R_c^g$; (7) characteristic radius of the gas and the small grain population; (8) surface density exponent in eq. (1); (9) scale height at 100 au; (10) power index of the radial dependence of scale height, eq. (3); (11) the inner edge of the disk model; (12) the outer edge of the disk model.

previous studies suggest the rate may be low in protoplanetary disks (Cleeves et al. 2015).

We let the chemistry evolve for 1 Myr for all models. Isotopologue fractionation is not considered in the chemical network, as fractionation is expected to be insignificant in massive disks (Miotello et al. 2014). We scale the output CO abundance by the local ISM ratios of CO/$^{13}$CO= 69, CO/C$^{18}$O= 557, CO/C$^{17}$O= 2005 to generate CO isotopologue abundances (Wilson 1999). We will discuss how this assumption may affect our results in Section 4.4.

### 3.1.5. Radiative Transfer modeling

For a given dust density structure and stellar spectrum, we calculate the dust temperature structure using the Monte Carlo radiative transfer code RADMC3D (Dullemond et al. 2012). These structures are then used to calculate synthetic SEDs and mm continuum images, using the ray-tracing module of RADMC3D. The disk inclination angle (i) and the major axis position angle (PA) are adopted from Table 1. Dust scattering is not included in our models of SED and continuum images, as it does not affect our results except for the inner 15 au region. The analysis of dust scattering effects of MAPS continuum images is explored in Sierra et al. (2021).

RADMC3D is also used to generate line image cubes. The CO line transitions are assumed to be in local thermal equilibrium (LTE). We generate model image cubes at 2-3 au resolution and 0.04 km s$^{-1}$ spectral resolution. The model channels are rebinned to the spectral resolution of MAPS observations and then are convolved with a Gaussian that has the same size and position angle of the beam of CO line observations. Then azimuthally-averaged radial intensity profiles of the model images

are generated using the radial_profile function in the Python package GoFish (Teague 2019), following the same procedure as those generated from the observations (Law et al. 2021a).

### 3.1.6. Constraining Disk Parameters

The disk structure model has 6 free parameters: 4 parameters to characterize the surface density profiles of gas and small dust, $\{\Sigma_c^g, \Sigma_c^{\mu m}, R_c, \gamma\}$, and 2 parameters to describe the scale height distribution $\{H_{100}, \psi\}$. Moreover, the surface density profile of the large grain population is determined by an iterative process to fit the mm continuum images (more details below).

The initial values of parameters are taken from literature for each disk (Cleeves et al. 2016; Macías et al. 2018; Fedele et al. 2018; Zhang et al. 2019; Liu et al. 2019). We then update the parameters by comparing model outputs with three types of observations: SEDs, mm-continuum images, and the CO emission surfaces (see detailed explanation below). The general modeling processes are outlined in Figure 3.

The mid and far-IR profile of the SED is mostly determined by the $\psi$ and $H_{100}$ parameters. Therefore, we run a small grid of $\psi$ and $H_{100}$ values to search for the best fitting value. The scale height of the large grain population is fixed to be 20% of the gas scale height (Dubrulle et al. 1995; Andrews et al. 2011). The SEDs from our best-fitting models are shown in Figure 4.

The second type of observational constraints is the spatially resolved mm-continuum images of the five disks. For the large dust grain population, we iterate an initial surface density profile until the azimuthally averaged surface brightness distribution matches with the high-resolution ALMA observations. The continuum
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Best-fit SED models of MAPS disks. The models are in blue and photometric data are in black. The photometric data of IM Lup, AS 209, and HD 163296 are adopted from Andrews et al. (2018), GM Aur from Macías et al. (2018), and MWC 480 from Anderson et al. (2013).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Surface density profiles of the gas, small, and large dust grain populations in our best-fit models of the five MAPS disks. The distributions of the gas and small dust grains are characterized by eq. (1), while the distributions of the large dust population are customized to match high-resolution continuum observations (see more details in Section 3.1.6).</figcaption></figure>

images of IM Lup, AS 209, and HD 163296 are from the DSHARP program (Andrews et al. 2018; Huang et al. 2018), MWC 480 image from Long et al. (2018) and GM Aur image from Huang et al. (2020).

Because our model is axisymmetric, we only compare azimuthally averaged radial profiles with observations. The ratio between the model image and ALMA observations is then used to update the input surface density

profile and generate the next model. After 3-4 iterations, the model image generally converges to a good match to observations. In some cases, the inner 10 au becomes optically thick, and thus increasing surface density alone cannot significantly increase the surface brightness. It suggests that our model underpredicts the mid-plane temperature inside 10 au, and the real dust scale height may deviate from the general power-law function we
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison between the observed 1.3 mm continuum images (top row) and simulated images from our thermo-chemical models (bottom row). The continuum images of the IM Lup, AS 209, and HD 163296 disks are from the DSHARP program (Andrews et al. 2018; Huang et al. 2018); the MWC 480 disk image from Long et al. (2018); Liu et al. (2019), and the GM Aur disk image from Huang et al. (2020). For each disk, the data and model panels use the same color scale. The model fitting processes are described in Section 3.1.6.</figcaption></figure>

used to match SEDs. Given our line observations have a spatial resolution of 15-24 au, we do not perturb the dust scale height to match the brightness in the inner 10 au. The surface density profiles and continuum images of best-fitting models are shown in Figure 5 and Figure 6, respectively.

The third type of observational constraint is the emission surfaces of optically thick CO and ¹³CO (2-1) lines. These emission surfaces are the vertical locations where a line has most of its emission at each radius. See Figure 7 for examples of emission surfaces. For optically thick lines, emission surfaces are roughly their $\tau \sim 1$ layer. With the high spatial resolution of MAPS data, the vertical locations of these surfaces can be directly measured using the asymmetry of the emission relative to the major axis of the disk. Please see Law et al. (2021b) for the exact procedure used for MAPS measurements.

Here, we use CO and ¹³CO (2-1) emission surfaces to constrain the vertical locations of the warm molecular layer of CO gas in our thermo-chemical models. Because the CO (2-1) line is highly optically thick, its emission surface is expected to be only slightly below the CO photodissociation layer (e.g., Pinte et al. 2018a). The ¹³CO (2-1) emission surface is expected to be slightly higher than the CO freeze-out surface. Therefore, these two emission surfaces provide rough upper and lower boundaries of the warm molecular layer of CO gas.

Comparing these emission surfaces to the CO abundance structures in thermo-chemical models provides a zeroth-order check on whether the models reasonably match with the observations. For example, the CO freeze-out surface in models should be always below the observed ¹³CO emission surface.

For a given stellar luminosity, the warm molecular layer of CO is mainly regulated by the scale height distribution and the amount of small grains in the disk. The amount of small grains affects how deep the stellar light penetrates and thus the $T_{gas} \simeq 20$ K layer where CO freezes out. We start with $H_{100}$ and $\psi$ parameters that match the SEDs. The large dust mass distribution is then constrained by mm-continuum images. For the third step, we run five models for each disk in which the mass fraction of the small grains takes 1%, 5%, 10%, 15%, and 20% of the total dust mass. We then compare the CO abundance structures of the five models with the observed CO emission surfaces and choose the one that best matches the observed surfaces by eye.

For the IM Lup, GM Aur, and HD 163296 disks, we find solutions that reasonably match both the SED and CO emission surfaces. For the AS 209 disk, we find that the thermo-chemical models with $\psi$ and $H_{100}$ reproducing the observed SED have a CO-freezeout layer much higher than the observed ¹³CO (2-1) emission surface. Changing the mass of small grains alone cannot solve the problem. We therefore choose to use scale height param-
eter values that match in CO emission surfaces instead of the SED, because our analysis is mainly based on the CO images. For the MWC 480 disk, our model matches the ¹³CO (2-1) emission surface but over predict the CO (2-1) emission surface. The best-fitting model parameters are listed in Table 2. Figure 7 shows the gas temperature and CO abundance structure in our best-fitting models.

### 3.1.7. Deriving CO column density

For a given disk, we derive the CO column density profiles ($N_{CO}$) from each CO line. Following the approach used in Zhang et al. (2019), we generate a grid of synthetic line images, by scaling the CO gas abundance in each model with a constant scaling factor throughout the disk. The model grid spans a wide range of scaling factor between 0.001 to 50, and the grid increases with a step-size of 1.1 in a logarithmic scale. For a given CO transition and a disk, we compare radial profiles from the model grid with the radial profile from observations. An example of the model grid is shown in Figure 8. At each radius, we adopt the local scaling factor from the best-fitting model and then create a composite 1D depletion profile of the CO gas. The depletion profile is then applied to the CO gas column density from our thermo-chemical model. A robustness test of this method is presented in the Appendix A. For a consistency check, we also apply the best-fitting 1D depletion profiles to CO abundance structures in thermo-chemical models and then generate simulated CO line radial profiles. Figure 26 in the Appendix shows that our best-fitting models match well with the observed profiles.

### 3.1.8. CO column density from thermo-chemical models

In Figure 9, we show profiles of CO gas column density ($N_{CO}$) for all five MAPS disks, based on four CO isotopologue transitions, i.e., (2-1) and (1-0) of C¹⁸O and ¹³CO. For the HD 163296 and MWC 480 disks, we also present $N_{CO}$ profiles derived from the C¹⁷O (1-0) line, as this line is sufficiently strong in these two disks.

In general, the $N_{CO}$ profiles are consistent among CO isotopologue lines, within a factor of 2. We also note that the spatial resolution of the (1-0) lines is a factor of two lower than the (2-1) lines, and therefore fewer substructures are seen in the $N_{CO}$ profiles from (1-0) lines. For the HD 163296 disk, the $N_{CO}$ derived from C¹⁷O and C¹⁸O lines are consistent outside 100 au but differ inside 100 au by up to a factor of five. We discuss this discrepancy in detail in Section 3.4. The derived $N_{CO}$ profiles show some small scale wiggles (10-20% variation on the scale of 10 au), which are numerical noise from our thermal chemical models. Here, we focus on larger scale features.

The $N_{CO}$ profiles from different transitions are in good agreement. Therefore we use the $N_{CO}$ profile from the C¹⁸O (2-1) line for general analysis since this line has a lower optical depth than ¹³CO lines, and has higher resolution and signal-to-noise ratio than the C¹⁸O (1-0) lines.

Among the five MAPS disks, their $N_{CO}$ distributions vary significantly in their general shapes and absolute values (see Figure 10). In the two disks around Herbig stars, the $N_{CO}$ profiles appear to be similar: the $N_{CO}$ shows a steep decrease with radius out to 100 au, which is followed by a slow decrease out to 400 au. In the three disks around T Tauri stars, however, the $N_{CO}$ profile varies significantly from one disk to another. For the IM Lup disk, the $N_{CO}$ shows a steep decrease with radius inside 20 au and a slow decrease between 50-500 au. The GM Aur disk shows a CO cavity inside 40 au. Beyond 40 au, its $N_{CO}$ profile is similar to that of the two disks around Herbig stars. For the AS 209 disk, $N_{CO}$ decreases with radius steeply inside 40 au, shows a broad shallow gap between 60-120 au, and a deeper gap around 240 au. The AS 209 disk has the smallest CO gas disk with a steep decrease between 150-300 au. Interestingly, the $N_{CO}$ of the remaining four disks show a very similar shallow slope between 150-400 au, see Figure 10. The slope can be characterized by a power-law function of $\Sigma(R) \propto R^{-2.4}$.

Comparing the absolute $N_{CO}$ value of the five disks, the largest differences occur inside 150 au, where the $N_{CO}$ in the two disks around Herbig stars are 1-3 orders of magnitude higher than that of the three disks around T Tauri stars. For individual disks, the absolute value of $N_{CO}$ varies by 4-5 orders of magnitude throughout a disk. The IM Lup disk shows the smallest range of $N_{CO}$ changes throughout its disk, while the MWC 480 disk shows the largest variation.

## 3.2. CO column density from empirical temperature structure

The $N_{CO}$ profiles derived above depend on gas temperature structures in thermo-chemical models. To test the robustness of these results, we employ a second method to measure $N_{CO}$ using the empirical gas temperature structure measured by Law et al. (2021b).

Gas temperature in the disk can be measured from the surface brightness of optically thick line emission provided the emission fills the beam (Weaver et al. 2018; Pinte et al. 2018a; Isella et al. 2018). Thanks to the high spatial resolution of MAPS data, the emission surfaces of lines can be directly measured (see our discussion in Section 3.1.6), which are a function of radius and vertical height $(r, z)$. For optically thick line emission, the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Gas temperature (top row) and CO abundance structure (bottom row) of the best-fit thermo-chemical models. The solid black lines are emission surfaces of ¹²CO and ¹³CO (2-1) lines measured from Law et al. (2021b). CO abundance structures agree with the observed emission surfaces relatively well, except for the MWC 480 disk.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** An example of a model grid used to derive the CO depletion profile. The example is for the ¹³CO (2-1) line of the HD 163296 disk. The grey dashed lines are models with different CO depletion factors, and the blue line is the radial profile from MAPS observations. At each radius, we adopt the local scaling factor of the best-fit model and then create a composite 1D depletion profile of the CO gas.</figcaption></figure>

gas temperature at these $(r, z)$ locations can be directly measured from their line surface brightness. Law et al. (2021b) measured gas temperatures at the emission surfaces of CO and ¹³CO (2-1) lines. These empirical temperature measurements are for discrete locations in a given disk. To construct the 2D gas temperature structure, Law et al. (2021b) used the two-layer model from

Dullemond et al. (2020) to fit the discrete temperature measurements.

In Figure 11, we compare the gas temperatures between our thermo-chemical models and empirical temperature structure from Law et al. (2021b). In general, thermo-chemical models match well with the empirical temperatures in the 20-40 K region. The temperature structures of IM Lup and AS 209 show the best match between these two temperature models. For the HD 163296 and MWC 480 disks, our thermo-chemical models underpredict the temperature at the CO (2-1) emission surfaces. The largest discrepancy is seen in the MWC 480 disk models: gas temperature in our thermo-chemical model is generally colder by 5-10 K than the empirical temperature structure. A caveat of our thermo-chemical models is that the gas temperature structure is calculated with the assumption of no-CO depletion. Because CO is an important coolant in the disk atmosphere, a strong depletion of CO can lead to a warmer atmosphere. Calahan et al. (2021b) and Schwarz et al. (2021) presented thermo-chemical models with spatially varying C/H elemental ratio for the HD 163296 and GM Aur disks, respectively. They found that the gas temperature becomes warmer at regions above Z/R >0.2 when accounting for CO depletion, but no significant changes were seen in the gas temperatures of the ¹³CO and C¹⁸O emission surfaces.

We replace the gas temperature distribution with the empirically-based temperature structure and then use the same approach described in Section 3.1.8 to derive
<figure><img src="image_10.png" /></figure>

**Figure 9.** CO column density profiles of five MAPS disks, based on $^{13}$CO, C$^{18}$O (2-1), (1-0) and C$^{17}$O (1-0) lines. The N$_{CO}$ are derived from thermo-chemical models described in section 3.1.8. More substructures can be seen in the N$_{CO}$ derived from (2-1) lines than that of (1-0) lines, because (2-1) line images have a factor of two higher resolutions (0.′′15 beams for (2-1) lines vs. 0.′′3 for (1-0) lines). The beam sizes are plotted at the upper right corner of each panel.

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

**Figure 10.** Left: Comparison of CO column density profiles among five MAPS disks. The N$_{CO}$ here is based on C$^{18}$O (2-1) line in thermo-chemical models. Right: The CO column density of four disks between 100 and 500 au. The dashed lines are power-law functions of $\Sigma(R) \propto R^{-2.4}$. It shows that the CO profiles of these four disks have a similar slope between 150 and 400 au.
the CO column density. The dust temperature and CO abundance structure are kept the same as the thermo-chemical models. Figure 26 in the Appendix shows our best-fitting models compared with the observed profiles.

Figure 12 shows the $N_{CO}$ profiles derived using the empirical temperature structures, compared with the profiles derived from our thermo-chemical models. The two $N_{CO}$ profiles of each disk generally have excellent agreement for most of the disk regions. For the IM Lup and AS 209 disks, the two $N_{CO}$ profiles are nearly identical. This is expected as these two disks have the best agreements in temperature structures between their thermo models and empirical temperature structures. For the GM Aur, AS 209, and HD 163296 disks, small differences are seen inside 100 au. The largest difference is in the MWC 480 models, as temperature from the RAC2D model generally underpredicts the gas temperature inside 150 au. The empirical temperature models also show a steep increase of $N_{CO}$ inside 100 au in the HD 163296 and MWC 480 disks, although the absolute $N_{CO}$ is smaller than the one based on the thermo-chemical models. These results indicate that exact $N_{CO}$ derived will depend on the temperature structure inside 100 au of these disks. Based on current best constraints of the gas temperature structures in the two disks, a steep increase of $N_{CO}$ is needed to reproduce the observed line intensity profiles.

### 3.3. CO column density from C$^{17}$O hyperfine line fitting

The third method we use to derive $N_{CO}$ is fitting C$^{17}$O hyperfine structure transitions. This method is independent of the temperature and CO abundance structures adopted in the thermo-chemical models.

The C$^{17}$O J=1-0 rotational transition is split into three components which can be fitted to obtain constraints on the C$^{17}$O column density, excitation temperature, and optical depth (Mangum & Shirley 2015). This method has been used recently in disk studies of C$_2$H, HCN, and CN (Bergner et al. 2019; Teague & Loomis 2020; Bergner et al. 2021; Guzmán et al. 2021). We present the first application of this method to C$^{17}$O in protoplanetary disks. For our sample, only the C$^{17}$O (1-0) images of the HD 163296 and MWC 480 disks have sufficient signal-to-noise ratio for this analysis.

The three components of C$^{17}$O (1-0) line have rest frequencies of 112.358777, 112.358982 and 112.360007 GHz (Klapper et al. 2003), as listed in Table 2 of Oberg et al. (2021). The 0.5 km s$^{-1}$ spectral resolution of the MAPS Band 3 data means that 2 of the components are always blended but the third component is ~2 km s$^{-1}$ offset

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Comparison between gas temperature structure of our thermo-chemical models (left) and the empirical temperature structure (right) from Law et al. (2021b). The circles are surface brightness temperature measurements from CO, $^{13}$CO, and C$^{18}$O (2-1) line observations. The emission surfaces of CO is the highest, $^{13}$CO is in the middle, and C$^{18}$O at the bottom, because of their relative abundances (Qi et al. 2011; Zhang et al. 2017). All panels use the same color scale.</figcaption></figure>

from the other two components, which is enough for us to model and fit the individual components.

Using the spectral shifting and stacking tool in Gofish, we generate spectra averaged over annuli across the disk. This technique removes most of the Keplerian broadening. This technique also increases the effective signal-to-noise of the data compared to a single-pixel spectrum. The HD 163296 and MWC 480 spectra were averaged over annuli 1/2 × the beam major axis resulting in radial bins of width 15 and 24 au, respectively.

Following the same method outlined in Bergner et al. (2021), we generate a model spectra for the C$^{17}$O (1-0)
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Comparison of CO gas column density derived using thermo-chemical models and CO depletion profiles and empirical temperature structures. The column density shown here is based on C¹⁸O (2-1) line.</figcaption></figure>

hyperfine transitions to fit the observed spectra. We assume LTE and that all three hyper-fine transitions have the same excitation temperature. The free parameters in the fit are: the total column density of C¹⁷O ($N_T$, cm⁻²), the excitation temperature ($T_{ex}$, K), the observed line width ($\sigma_{obs}$, km s⁻¹), and systematic velocity ($V_{lsrk}$, km s⁻¹). The $\sigma_{obs}$ parameter is larger than the local thermal+turbulence broadening, likely due to a combination of (1) emission contribution from the back and front sides of a flared disk, which have different projected velocities, (2) beam smearing, in which a wide range of Keplerian velocities are incorporated within the same beam, (3) Hanning smoothing of the data which effectively reduces the spectral resolution (Bergner et al. 2021).

The bounds for the column density and excitation temperature were set to $10^{10} - 10^{22}$ cm⁻² and 10–50 K. This encompasses the range of temperature for the C¹⁷O emitting region (see Section 3.2). $\sigma_{obs}$ is set to range from 0.01 to 8 km s⁻¹ and the $V_{lsrk}$ is between 0 and 10 km s⁻¹. The total line optical depth is the sum of the individual optical depths at the line center of each of the three transitions. This includes a correction factor to account for beam smearing that will artificially lower the inferred optical depth. For full details of the fitting procedure, see Bergner et al. (2021).

We use the MCMC package `emcee` (Foreman-Mackey et al. 2013) to sample the posterior distributions of the four fit parameters. The radially-averaged spectra from `Gofish` and the model fits for HD 163296 are shown in Figure 13. The same plot for the MWC 480 disk and the fits for N(C¹⁷O), $T_{ex}$, and $\tau$ for both disks are shown in the Appendix in Figure 27. The C¹⁷O line becomes optically thick ($\tau > 1$) within ≈55 au in the HD 163296 disk and ≈100 au in the MWC 480 disk. The $T_{ex}$ is poorly constrained as all of the transitions have the same upper-energy level but the derived column density is not sensitive to the assumed $T_{ex}$. We tested fixing the $T_{ex}$ to 20, 25, 30, 35 and 40 K for the HD 163296 disk and

in the region where the line is optically thin, the column densities were well within the fit uncertainties.

The N(C¹⁷O) measured from the hyperfine-line fitting method has an inner and outer limits of the radius. Inside the inner 1.5× beams, the fit is poor due to the complete blending of the lines and beam smearing. Beyond ~250 au the spectra are too noisy to obtain robust fits in the HD 163296 and MWC 480 disks.

The resulting N(C¹⁷O) radial profiles are then converted to total N(CO) radial profiles under the assumption that ¹⁶O/¹⁸O= 557 and ¹⁸O/¹⁷O= 3.6 (Wilson 1999).

Figure 14 presents the $N_{CO}$ profiles derived from the C¹⁷O fitting compared with the profiles derived in Section 3.1 for HD 163296 and MWC 480 from the C¹⁸O (2-1) and C¹⁷O (1-0) lines. In general, the profiles are consistent with each other, i.e., within a factor of a 2-3. For the HD 163296 disk, the $N_{CO}$ derived from C¹⁷O (1-0) using both methods suggests a higher CO column density between 30-150 au than that of C¹⁸O (2-1) line. See detailed discussions on the difference in the following section.

Overall, we have shown that with sufficiently high signal-to-noise data, fits to the C¹⁷O (1-0) hyperfine transitions can provide good constraints on the CO column density.

### 3.4. CO column density inside 100 au of HD 163296 and MWC 480

Inside 100 au, the HD 163296 and MWC 480 disks have a CO column density that is 1-2 orders of magnitude higher than that of the other three disks around T Tauri stars. The C¹⁸O (2-1) line emission interior to ~100 au in these two disks is optically thick. In these regions, the $N_{CO}$ derived from C¹⁸O thus relies on the robustness of the vertical temperature and CO abundance in models. The $N_{CO}$ derived from the C¹⁷O (1-0) line of HD 163296 shows a 2-6 times higher column density inside ~100 au, while the $N_{CO}$ inferred from the C¹⁷O and
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Radially averaged C<sup>17</sup>O (1-0) spectra (black) and the results from the hyperfine line fitting (blue) for the HD 163296 disk. The center radius of each bin is shown and the spectra are offset along the axis for clarity.</figcaption></figure>

C<sup>18</sup>O lines match well in the MWC 480 disk. Both disks have previous detections of rarer CO isotopologue lines, including <sup>13</sup>C<sup>18</sup>O and <sup>13</sup>C<sup>17</sup>O lines (Booth et al. 2019; Zhang et al. 2020a; Loomis et al. 2020), which provides constraints closer to the midplane in the two disks. Here we use the disk-integrated <sup>13</sup>C<sup>18</sup>O line spectra to test the robustness of the N<sub>CO</sub> results inside 100 au for the HD 163296 and MWC 480 disks. We apply the CO depletion profiles from the C<sup>18</sup>O and C<sup>17</sup>O lines to the CO abundance from our thermo-chemical models and generate simulated spectra for <sup>13</sup>C<sup>18</sup>O (2-1) or (3-2). The results are shown in Figure 15.

For HD 163296, the N<sub>CO</sub> from C<sup>17</sup>O (1-0) line has a better match to the <sup>13</sup>C<sup>18</sup>O (2-1) line spectrum. The high N<sub>CO</sub> is consistent with previous results from Zhang et al. (2020a), who found that the CO abundance inside 70 au in the HD 163296 disk is a factor of few larger than that in the outer disk region. As the <sup>13</sup>C<sup>18</sup>O and C<sup>18</sup>O (2-1) lines are at similar frequencies, the difference of N<sub>CO</sub> cannot be explained by optically thick dust. One possible explanation is that the CO-to-H<sub>2</sub> gas abundance ratio close to the mid-plane may be a factor of few higher than that of the disk atmosphere inside 70 au. This enhanced CO abundance at the mid-plane is possible if a large amount of icy pebbles evaporate their CO ice mantle at the mid-plane (Krijt et al. 2018, 2020). For the MWC 480 disk, the N<sub>CO</sub> discrepancy inside 100 au is smaller and both results match the <sup>13</sup>C<sup>18</sup>O (3-2) line within uncertainties. Higher spatial resolution observations of <sup>13</sup>C<sup>18</sup>O are necessary to better constrain the N<sub>CO</sub> inside 100 au of these two disks.

# 4. DISCUSSION

## 4.1. Minimum gas disk masses

Despite the large uncertainty of CO abundance in disks, the amount of CO gas in disks provides a lower limit of total gas mass. Here we estimate minimum gas disk masses in two ways: (1) we calculate the total CO gas mass for each disk, and then scale it with a constant CO-to-H<sub>2</sub> abundance ratio of 10<sup>-4</sup> to get a total disk gas mass, M<sub>gas1</sub>. This is a lower limit of disk gas mass, as it does not include gas mass in the CO freeze-out and photodissociation regions. (2) We scale the gas surface density in our disk models with the CO depletion profile, and calculate a total disk gas mass, M<sub>gas2</sub>. This approach assumes that the difference between thermo-chemical models and observations is caused by gas depletion and thus our CO depletion factor is a gas depletion factor. In this way, the gas masses are corrected for the CO freeze-out and photodissociation effects, but the result depends on the temperature and CO abundance structure in models. For the mass estimation, we use the
<figure><img src="image_15.png" /></figure>

**Figure 14.** Comparison of CO column densities derived from hyperfine line models with that from thermo-chemical models of the HD 163296 and MWC 480 disks. The blue and purple lines are $N_{\text{CO}}$ profiles from the C$^{18}$O (2-1) and C$^{17}$O (1-0) lines, respectively, in thermo-chemical models (Section 3.1.8). The purple colored markers are the best-fitting results for the hyperfine lines of C$^{17}$O (1-0) (Section 3.3). The x-scale error bars denote the size of the radial bin that the spectra are averaged over. In the results from hyperfine line fitting, we do not show points from the inner 1.5× the beam major axis, where the lines are too blended, or points from the outer disk, where the signal-to-noise of the spectra is too low.

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

**Figure 15.** Comparison of $^{13}$C$^{18}$O spectra (in grey) with models based on $N_{\text{CO}}$ from C$^{18}$O (2-1) and C$^{17}$O (1-0) lines in thermo-chemical models (Section 3.1).

CO depletion profiles based on the C$^{18}$O (2-1) line from empirical temperature models (Section 3.2). The only exception is that for the inner 150 au of the HD 163296 disk, we adopt the $N_{\text{CO}}$ from C$^{17}$O (1-0) line, due to the discrepancy discussed in Section 3.4. We note that our initial chemical condition in thermo-chemical models has a CO-to-H$_2$ abundance of $2.8 \times 10^{-4}$ (Lacy et al. 1994), and thus the undepleted CO abundance in warm molecular layer is between $1-2.8 \times 10^{-4}$, slightly higher than the canonic assumption of $10^{-4}$. Due to the differ-

ence in assumed CO-to-H$_2$ abundance ratios, even $M_{\text{gas2}}$ includes corrections of freeze-out and photodissociation, in the IM Lup and MWC 480 disks, the $M_{\text{gas2}}$ are still slightly smaller or the same as $M_{\text{gas1}}$.

The estimations of minimum gas disk masses are summarized in Table 3. The total dust masses are adopted from our best-fit thermo-chemical models (see Table 2). These five disks have a lower limit of gas-to-dust ratio between 1 and 24, demonstrating that even if no CO depletion is assumed, these disks still have a gas-to-dust ra-
tio ≥1. Based on the CO derived minimum gas masses, the two disks around Herbig stars have a higher gas-to-dust mass ratio than those around T Tauri stars. This is consistent with the expectation that disks around Herbig stars would have less depletion of CO, due to their warmer disk conditions (Bosman et al. 2018; Schwarz et al. 2018, 2019). Interestingly, the IM Lup disk has the lowest CO gas-to-dust mass ratio, despite the fact that it is the youngest disk in the sample.

In Figure 16, we compare the radial distributions of CO derived gas mass (method 2 above) with dust mass in our models. It shows that if the CO intensity variation is due to gas depletion, the IM Lup disk would have a gas-to-dust ratio ≤1 inside 200 au. Analyses of other molecules in the IM Lup disk showed it is unlikely the case, based on the observations of N-carriers (Cleeves et al. 2018). Therefore, the IM Lup disk likely has a low CO abundance rather than a gas-to-dust ratio of 1 in the disk.

## 4.2. Radial dependence of CO depletion in MAPS disks

In Figure 17, we show the CO depletion profiles derived from thermo-chemical models. Compared to our thermo-chemical models, observations require 1-2 orders of magnitude lower CO abundance in most of the disk region. This is consistent with previous findings that the disk-averaged CO gas may be heavily depleted in some gas-rich Class II disks (Schwarz et al. 2016; McClure et al. 2016; Zhang et al. 2017, 2019).

Various chemical and physical processes can reduce CO gas abundance in the disk atmosphere. Chemical processes can reprocess CO into other carbon-carriers (Aikawa et al. 1999; Bergin et al. 2014; Eistrup et al. 2016; Yu et al. 2016; Schwarz et al. 2018; Bosman et al. 2018). Additionally, dust growth and settling can sequester CO into the mid-plane and the radial drift of ice pebbles can bring CO into the inner disk region (Xu et al. 2017; Krijt et al. 2018; Booth & Ilee 2019; Krijt et al. 2020). These processes do not occur at the same rate at each radius, and thus the radial and vertical variation of CO abundance in disks can provide useful constraints on the CO depletion processes in disks.

Zhang et al. (2019) analyzed the radial variation of CO depletion in five disks using C$^{18}$O (2-1) or (3-2) line observations. They found that the CO gas abundance is heavily depleted in the region just beyond the mid-plane CO snowline and the CO abundance tends to be higher at the outermost region of the gas disk. They noted that the CO abundance inside the CO snowline of the HD 163296 disk is one order of magnitude higher than the region outside, which is further confirmed by a follow-up study of $^{13}$C$^{18}$O (2-1) (Zhang et al. 2020a).

The work here provides higher resolution constraints on the CO depletion profiles than previous observational studies. In Figure 17, we show that the two disks around Herbig stars (HD 163296 and MWC 480) exhibit a very similar depletion pattern: outside the mid-plane CO snowline, CO is depleted by a factor of ten or more; inside the CO snowline, CO gas abundance increases rapidly. This pattern is consistent with the dust evolution picture: CO is sequestered to the mid-plane as icy dust grains settle, and icy pebbles drift inwards into the region inside the mid-plane CO snowline, where CO ice evaporates and enhances the CO abundance at the inner disk regions (Booth et al. 2017; Booth & Ilee 2019; Krijt et al. 2018, 2020).

The three disks around T Tauri stars, however, do not show a clear trend. For the IM Lup disk, its CO abundance increases monotonically with radius. For the GM Aur disk, beyond its 40 au cavity, its CO abundance has a sharp decrease out to 80 au and then a slow decrease out to the edge of the disk. This profile is similar to that of the two disks around Herbig stars, but 80 au is much further out than the mid-plane CO snowline of 30 au in our chemical model of the GM Aur disk. The AS 209 disk has two local peaks around 30 and 150 au. In our models, the mid-plane CO snowlines of these T Tauri sources are comparable to or slightly smaller than our spatial resolution. But we do not see signals of CO enrichment like that in the two disks around Herbig stars.

## 4.3. CO column density distribution vs. mm-sized dust

### 4.3.1. CO gap properties

In Figure 18, we show the CO gas profiles and mm-sized dust distributions in our best-fit models. At our spatial resolution (10-24 au), a few substructures in $N_{CO}$ are noticeable. The most notable features are the AS 209 disk has a deep gap (~97% depletion) around 240 au, and GM Aur has a deep inner cavity inside 40 au
<fnref n="2" /> Except these features, other $N_{CO}$ substructures are more subtle compared to substructures seen in the dust distribution. To better characterize these substructures, we subtract a smooth function from the $N_{CO}$ profiles to identify significant features in the residuals. We then fit Gaussian functions to the residuals to estimate the widths and locations of gaps.

<footnote n="2">² We note that the $N_{CO}$ of the IM Lup disk appears to have a gap between 20-80 au. But it is possibly just a result of the high dust optical depth inside 20 au in our model (see Figure 16). For the IM Lup disk, all CO profiles show a depression inside 80 au, and its 1.3 mm continuum emission shows a steep increase inside 20 au. Given the large uncertainty of optical depth inside 20 au, we do not consider it as a real gas gap here.</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Radial profiles of gas and dust in MAPS disks. The gas surface density is derived by scaling gas surface density of thermo-chemical models with CO depletion profiles (see details in method 2 in Section 4.1). The dust surface density is from the best-fit thermo-chemical models. The red shaded region indicates the mid-plane CO snowline locations in our thermo-chemical models.</figcaption></figure>

<table><caption>**Table 3.** Minimum disk mass and snowline location</caption><thead><tr><th>Source</th><th>Mgas1<br>(10<sup>-3</sup> M<sub>☉</sub>)</th><th>Mgas2<br>(10<sup>-3</sup> M<sub>☉</sub>)</th><th>g2d1</th><th>g2d2</th><th>r<sub>CO</sub><br>(au)</th><th>r<sub>N<sub>2</sub></sub><br>(au)</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th></tr></thead><tbody><tr><td>IM Lup</td><td>1.9</td><td>1.9</td><td>1</td><td>1</td><td>15±5</td><td>25±5</td></tr><tr><td>GM Aur</td><td>8.0</td><td>19.9</td><td>12</td><td>29</td><td>30±5</td><td>70±5</td></tr><tr><td>AS 209</td><td>0.4</td><td>1.0</td><td>1</td><td>2</td><td>12±5</td><td>20±5</td></tr><tr><td>HD 163296</td><td>31.3</td><td>42.7</td><td>14</td><td>19</td><td>65±5</td><td>90±5</td></tr><tr><td>MWC 480</td><td>35.9</td><td>28.6</td><td>25</td><td>20</td><td>100±5</td><td>150±5</td></tr></tbody></table>

NOTE—(1) Source name (2) $M_{gas1}$: the minimum disk gas masses by scaling the total CO gas mass with $n_{CO}/n_{H_2} = 10^{-4}$. (3) $M_{gas2}$: the gas disk mass by scaling the gas surface density at each radius by CO depletion profiles derived in Section 3.2. (4) Minimum gas-to-dust mass ratio based on $M_{gas1}/M_{dust}$. (5) $M_{gas2}/M_{dust}$. (6) mid-plane CO snowline location in the our thermo-chemical model. (7) mid-plane $N_2$ snowline location.
<figure><img src="image_18.png" /><figcaption>**Figure 17.** CO depletion profiles derived from C<sup>18</sup>O (2-1) and C<sup>17</sup>O (1-0) line observations. We show both depletion profiles derived from thermo-chemical models (solid line, Section 3.1). The red shaded region indicates the mid-plane CO snowline location in our thermo-chemical models. The depletion profiles show that the CO abundance is depleted by a factor of 10-100 in the three disks around T Tauri stars. For the two disks around Herbig stars, the CO abundances are depleted by a factor of ~10 outside 150 au, but enhanced to the ISM level inside their mid-plane CO snowlines.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

**Figure 18.** Comparison of CO gas and dust gaps. Top plot: row1, CO column density distribution (solid line). The regions filled in light blue color are CO gaps identified. row 2: mm-dust distributions. The regions filled in blue are dust gaps identified. row 3: Comparison between gas and dust gaps in logarithmic scale. The bottom plot: Cross-comparison of the locations and depths of the CO gas and dust gaps among the five MAPS disks (depth is shown in linear scale).

Table 4 lists the CO gap properties. We note that the exact properties will likely vary with the particular choice of smooth functions. In general, when CO gaps are present in the pebble disk region, their gap locations roughly coincide with the gap locations of mm-sized dust, suggesting a correlation between these pairs. However, this is not a one-to-one correlation, as some of the deepest dust gaps do not have CO gap pairs.

The AS 209 disk shows a broad CO gap between 40-120 au, spanning the range of three deep dust gaps. Our fitting requires two overlapping Gaussian functions at 61 and 94 au to match this broad gap. The peaks of the two Gaussians are consistent with the locations of two deep dust gaps at 61 and 90 au. The AS 209 disk also has a deep CO gap around 236 au, which is beyond the pebble disk region. The GM Aur disk does not show clear
<figure><img src="image_20.png" /><figcaption>**Figure 19.** CO gaps in the HD 163296 disk. Left: Teague et al. (2018a) gas distribution model that reproduces the velocity deviations from a Keplerian rotation field in C$^{18}$O (2-1) line observations. The dashed line is the smooth function used to identify gap properties. Right: Comparison of gas gaps in Teague et al. (2018a) with the CO gaps found in this work.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Comparison of observed CO gap properties (depth and width) with predictions of planet-disk interactions. The top panels show the depths and widths of CO gaps (red diamonds) and predictions based on mm-continuum gaps (black circles). See detailed derivations in Section 4.3, and numbers are listed in the Appendix Table 5. The bottom left panel: the ratios of CO gap depth over the gas gap depth predicted. It shows that the CO gap depths are generally consistent with predictions based on mm-continuum gaps, but there are also cases where CO gaps are 5-10 times shallower compared to the predictions. Bottom right panel: the same as the bottom left panel, but for the ratio of gap width. It shows that the observed widths are slightly larger or comparable to the expectations.</figcaption></figure>
substructures in its CO distribution beyond its central cavity, even though it has two deep dust gaps at 70 and 120 au. The HD 163296 disk has three CO gaps at 43, 94, and 148 au, and these CO gaps are close to three dust gaps at 48, 86, and 145 au. The MWC 480 disk has one broad gap at 63 au, offset from its dust gap at 73 au. In terms of the gap depth, CO gaps inside the pebble disk region only show a modest depression, where the N<sub>CO</sub> is 38-68% of that at nearby regions. The deepest CO gaps in the sample is the 236 au gap of the AS 209 disk, which is outside of the pebble disk. The density inside this gap is only ~3% of the background density.

<table><caption>**Table 4.** Properties of CO gaps</caption><thead><tr><th>Source</th><th>r<sub>gap</sub><br>(au)</th><th>σ<br>(au)</th><th>depth</th></tr><tr><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th></tr></thead><tbody><tr><td>GM Aur</td><td>14±1</td><td>14±1</td><td>0.95±0.03</td></tr><tr><td rowspan="3">AS 209</td><td>59±9</td><td>13±6</td><td>0.47±0.27</td></tr><tr><td>92±12</td><td>18±8</td><td>0.55±0.12</td></tr><tr><td>238±1</td><td>13±1</td><td>0.97±0.08</td></tr><tr><td rowspan="3">HD 163296</td><td>44±1</td><td>7±1</td><td>0.49±0.04</td></tr><tr><td>93±1</td><td>16±1</td><td>0.42±0.03</td></tr><tr><td>148±1</td><td>4±1</td><td>0.32±0.05</td></tr><tr><td>MWC 480</td><td>63±1</td><td>17±1</td><td>0.62±0.03</td></tr></tbody></table>

NOTE—(1) Source name (2) CO gap location (3) gap width derived by fitting a Gaussian. σ is the standard deviation of the best-fitting Gaussian function. (4) gap depth, defined as $(\Sigma_{\text{base}} - \Sigma_{\text{gap}})/\Sigma_{\text{base}}$.

### 4.3.2. Testing the planet scenario: can CO and dust gaps be explained by the same planet masses?

It is still unclear what mechanisms cause these substructures in the CO gas distribution. Here we briefly discuss planet-disk interactions as one possible cause.

In our sample, the HD 163296 disk shows the strongest correlation between CO and dust gap locations. Various previous studies have suggested that the HD 163296 disk hosts multiple giant planets. One evidence is based on substructures in its 1.25 mm continuum image, for which three giant planets between 0.1-4 M<sub>J</sub> were needed to explain the dust gaps at 10, 48, and 86 au (Zhang et al. 2018). Another evidence is from local velocity deviations from Keplerian rotation (Teague et al. 2018a, 2019). Hydro-dynamical models showed that three giant planets at 54, 83, 136 au provide the best match to the

velocity perturbations<fnref n="3" />. In Figure 19, we compare our N<sub>CO</sub> gap locations and depths with gas surface density derived from velocity deviation of Teague et al. (2018a). The widths and depths of CO gaps are comparable to the gas gaps in Teague et al. (2018a) model, but the peaks differ by up to 10 au and the CO gap at 148 au is 2-3 times narrower than the one in their model. In addition to these putative planets within the pebble disk, Pinte et al. (2018b) suggested a 2 M<sub>J</sub> planet at 260 au, using local deviation from Keplerian velocity. This velocity deviation has been further confirmed in the CO observations of the DSHARP and MAPS programs (Pinte et al. 2020; Teague et al. 2021). But we do not find any clear CO gaps around 260 au.

For the AS 209 disk, we identify two CO gaps. The inner gap between 50-120 au coincides with two mm-continuum gaps at 61 and 94 au. Zhang et al. (2018) attributed these two dust gaps to a single planet (0.2-1.3 M<sub>J</sub>) at 99 au. Favre et al. (2019) showed that a planet of 0.2-0.3 M<sub>J</sub> at 100 au is roughly consistent with their C<sup>18</sup>O and <sup>13</sup>CO (2-1) observations. However, using higher spatial resolution CO and C<sub>2</sub>H observations, Alarcón et al. (2021) showed that the CO gap is at least partially caused by a local depletion of CO abundance and thus ruled out the presence of > 0.2 M<sub>J</sub> planet around 100 au. The AS 209 disk also has a deep CO gap around 236 au, which coincides with a local pressure minimum identified through velocity deviation (Teague et al. 2018b). This gap location also potentially coincides with a ring in scattered light image of the AS 209 disk (Avenhaus et al. 2018). The origin of the 236 au gap requires further investigation.

Figure 18 shows that CO gaps are much shallower than the mm-continuum gaps. This is qualitatively consistent with the predictions of planet-disk interactions that gas gaps are shallower than gaps of pebbles (e.g., Zhang et al. 2018). It is of particular interest to see whether the CO gas continuum gaps are quantitatively consistent with the expectation of planet-disk interactions.

Here we do a quick check on whether the CO and continuum gaps can be explained by the same planet masses. Following the method developed by Zhang et al. (2018), we first measure continuum gap widths, and then combine these widths with scale height ratios and gas surface densities in our best-fit disk parameters to derive planet masses (assuming a viscosity parameter α=10<sup>-3</sup>). We then use the inferred planet masses and

<footnote n="3">³ The original planet locations reported in Teague et al. (2018a) were based on a distance of 122 pc for HD 163296. Here we corrected the planet locations using the latest GAIA distance of 101 pc.</footnote>
eqs. (5-6) in Kanagawa et al. (2015) to derive the ex-
pected gas gap depths<fnref n="4" />. We use eq. (22) in Zhang et al.
(2018) to derive the expected gas gap width. The contin-
uum gap widths, disk parameters, and estimated planet
masses are listed in the Appendix in Table 5. We note
that these values are just used for a consistency check
of the inferred planets by CO and continuum, and the
results should not be considered rigorous derivations of
planet masses.

Figure 20 shows the comparison of CO gas properties with expected properties of gas gaps. For gap depths, the CO gas depths are generally consistent with the expectations within a factor of 2. But there are two cases where the CO gaps are 5-10 times shallower than the predictions: the 48 au gap of HD 163296, and the 67 au gap of GM Aur. For gap widths, the observed CO gap widths are consistent with the expected values within a factor of 3, but tend to be larger. This is partially due to our beam size of 15-24 au, but the beam size difference is not sufficient to explain the two shallow gaps mentioned above. Nevertheless, our simple estimation of gap depths does not consider possible complications, such as CO abundance variations across a gap, the strength of viscosity, and impact of planet migration. Therefore, the discrepancy between the observed CO gap depths and the predictions in some disks does not necessarily rule out the possibility of planets at these places. To further test the planet scenario, it will be beneficial to combine constraints from different methods, e.g., velocity deviations from Kepleration rotation (e.g., Teague et al. 2018b; Rab et al. 2020; Teague et al. 2021; Alarcón et al. 2021), and other chemical tracers like HCO⁺ (e.g., Smirnov-Pinchukov et al. 2020).

4.4. *Possible effects of isotope-selective photodissociation*

Our thermo-chemical models do not include isotope-
selective photodissociation effects. Here we discuss how
it might affect our results of the CO column densities.

CO is one of the few molecules that can self-shield
itself from UV radiation, as its photodissociation oc-
curs by absorbing UV photons at discrete wavelengths
and being subsequently excited to predissociated bound
states (e.g., van Dishoeck & Black 1988; Visser et al.
2009). Due to the differences in abundances, there is
a transition zone where the ratios of ¹²CO/C¹²O and
¹³CO/C¹⁸O are higher than the elemental isotopic ra-

tios. If the transition zone takes a significant fraction
of the total CO gas column, the NCO derived from
C¹⁸O and an ISM element isotopic ratio would be lower
than the true NCO (Miotello et al. 2014). Besides CO
itself, H₂ and small dust particles in the disk provide
additional attenuation to UV radiation. In general, the
more massive a disk is, the smaller the isotope-selective
photodissociation effect for the total CO column density.
For example, Miotello et al. (2014) showed that for an
ISM ratio of CO-to-H₂ abundance the isotope-selective
photodissociation is unimportant for disks >10⁻² M☉.

For the five MAPS disks studied here, we expect the
effects of isotope-selective photodissociation can only
change our results up to a factor of 2-3. On the global
disk scale, the CO column densities derived from ¹³CO
and C¹⁸O lines are consistent with a constant ratio of
¹³CO/C¹⁸O=557/69, as far out in the disks that C¹⁸O is
detected. In contrast, models with significant isotope-
selective photodissociation predict that ¹³CO/C¹⁸O ra-
tio is 3-10 times higher (e.g., Miotello et al. 2014). The
isotope-selective effect could also be important inside
dust gaps, if the gas and dust densities are reduced sig-
nificantly. However, as shown in Figure 18, the C¹⁸O de-
pletions inside gaps are only modest or even not detected
inside the mm-continuum gap locations. It is unlikely
that isotope-selective photodissociation is much stronger
inside these gaps compared to the nearby regions. Fur-
ther investigation of the UV intensity inside these gaps
can be done with modeling of small dust particles (based
on scattered light images), but it is beyond the scope of
this work. In short, isotope-selective photodissociation
effects would not significantly change our results.

4.5. *Gas distribution beyond 150 au: evidence of viscous disks*

The global mass distribution of a disk profoundly af-
fects planet formation and migration (e.g., Morbidelli &
Raymond 2016). Turbulent viscosity and disk winds are
currently two leading mechanisms proposed to explain
global disk mass transportation (Shakura & Sunyaev
1973; Balbus & Hawley 1991; Bai & Stone 2013). These
two mechanisms can result in very different gas mass
distributions, especially at the <10 au and >100 au re-
gions of the disk. For example, current disk wind models
generally predict steeper density profile beyond 100 au
compared to a pure viscous disk (Bai 2016; Suzuki et al.
2016). The difference in the outer disk is mainly because
a viscous disk spreads with time as angular momentum
moves to larger radii, while in the disk wind case the
disk size needs not grow with time. Global disk prop-
erties, such as accretion rates and disk sizes have been
tried to constrain viscous disk evolution models but ob-

<footnote n="4">⁴ We note that we did not use Zhang et al. (2018)'s eq. (24) to link planet masses with gas gap depth, because their definition of gas depth does not work well for the steep N<sub>CO</sub> profiles we derived in this work.</footnote>
servations showed relatively large scatter (Mulders et al. 2017; Najita & Bergin 2018; Hendler et al. 2020). The high resolution MAPS measurements of CO distribution now provide detailed constraints on the gas mass distribution in the outer disk region. Here we compare the N<sub>CO</sub> distributions with theoretical predictions of gas mass distribution at the outer disk region.

As shown in Figure 10, four out of five MAPS disks show similar N<sub>CO</sub> distributions between 150-400 au. Except for the AS 209 disk, the remaining four disks all show a smooth long tail, which can be well characterized by a power-law function of N<sub>CO</sub> ∝ r<sup>-2.4±0.2</sup>. This N<sub>CO</sub> tail is a very shallow profile compared to current predictions of disk wind models (Bai 2016; Suzuki et al. 2016) but comparable to that of viscously evolved disk models.

For illustrative purposes, we compare our N<sub>CO</sub> profiles with that from thermo-chemical models of viscous disks of Trapman et al. (2020). We adopt N<sub>CO</sub> profiles from three viscous disk models with α = 10<sup>-3</sup> around a 1 M<sub>☉</sub> star after 1 Myr of evolution. These models have an initial disk mass of 0.06 M<sub>☉</sub> and their initial disk size, R<sub>init</sub>, is 30, 50, and 100 au, respectively. These models predict that the slope of CO column density distribution is similar to the gas between 100-500 au region (see Figure 28 in the Appendix). We scale their absolute CO column density by a factor of 0.1 to match our N<sub>CO</sub> measurements. The differences are likely because these models do not consider CO depletion. Figure 21 shows the N<sub>CO</sub> measured in this work with that of viscous models. Between 150-400 au, the measured N<sub>CO</sub> profiles of the IM Lup, HD 163296, and MWC 480 disks match well to the R<sub>init</sub>=100 au model, and the GM Aur disk to the R<sub>init</sub>=50 au model. The AS 209 disk is similar to the R<sub>init</sub>=30 au model but likely requires an even smaller R<sub>init</sub> for better match. Beyond 400 au, the measured N<sub>CO</sub> profiles decrease faster compared to the models. The faster decreases might be caused by photodissociation of CO or real gas surface density changes due to photoevaporation.

In short, the N<sub>CO</sub> distributions in 150-400 au of the five MAPS disk match well with predictions of viscously evolving disk models. The initial evidence shows the importance of spatially resolved N<sub>CO</sub> profiles to constrain current disk evolution theories. A systematic study of N<sub>CO</sub> profiles in a larger disk sample and comparison with both viscously evolving and wind-driven models are needed for further investigation.

## 4.6. CO snowlines and continuum substructures

Snowlines of major volatile species have been proposed as a possible cause of substructures seen in (sub)mm

continuum of disks (Zhang et al. 2015; Okuzumi et al. 2016). The expectation is that the dust size distribution and mass surface density have a significant discontinuity at the snowline locations, as predicted by current models of icy dust growth (e.g., Ros & Johansen 2013; Banzatti et al. 2015; Okuzumi et al. 2016; Pinilla et al. 2017; Ros et al. 2019). Here we compare the locations of CO snowline in our models with continuum substructures.

In our thermo-chemical models, the mid-plane CO snowline is defined as the radius where the CO gas and ice abundance become equal at the disk mid-plane. The mid-plane snowline locations of IM Lup and AS 209 disks are slightly warmer than the CO condensation temperature (~20 K), because in our models a significant fraction of CO ices at the mid-plane has been processed into other species after 1 Myr). The snowline locations in our models are listed in Table 3. These snowline locations are consistent with the observations of C<sup>17</sup>O (1-0) line in the HD 163296 and MWC 480 disks, in which the C<sup>17</sup>O line becomes optically thick in the inner disks of at ~60 au and 80 au, respectively. Interestingly, a bright continuum ring is seen at the CO snowline location of both disks, as shown in Figure 16. This bright continuum ring is consistent with a pile-up of dust grains due to the sintering effect at snowlines (Okuzumi et al. 2016; Banzatti et al. 2015; Pinilla et al. 2017).

## 5. CONCLUSIONS

In this work, we present high resolution (15-24 au) observations and analysis of CO isotopologue lines from the MAPS ALMA Large Program. Our analysis employs C<sup>18</sup>O and <sup>13</sup>CO J=(2-1), (1-0), and C<sup>17</sup>O (1-0) lines of five protoplanetary disks around IM Lup, GM Aur, AS 209, HD 163296, and MWC 480. Our findings are summarized as follows:

* We retrieve CO gas density distributions, using three independent assumptions for the underlying gas temperature: (1) a thermo-chemical modeling framework based on the CO data, the broadband spectral energy distribution, and the mm continuum emission; (2) an empirical distribution based on optically thick CO emission lines; and (3) a direct estimate based on fits to the C<sup>17</sup>O hyperfine lines. Results from all three methods generally show excellent agreement.

* By comparing model temperature structures with empirical measurements from optically thick lines, we show that current thermo-chemical models can provide reliable temperature structures within the 15-40 K region.
<figure><img src="image_22.png" /><figcaption>**Figure 21.** Comparison of CO gas column density profiles of MAPS disks with that of viscously evolving disk models (Trapman et al. 2020). The viscous disk models are for $\alpha = 10^{-3}$ around a $1\,\text{M}_\odot$ star after 1 Myr evolution. The initial disk mass in the models is $0.06\,\text{M}_\odot$ and the initial radii, $R_{\text{init}}$, is 30, 50, and 100 au, respectively. Trapman et al. (2020) did not include CO depletion, and therefore we scaled their CO column density by a factor of 0.1 to better match our results. The CO distributions in the IM Lup, GM Aur, HD 163296, and MWC 480 disks match well with the viscous models of $R_{\text{init}}=50, 100$ au between the 150-400 au region. The CO profile of the AS 209 disk is more consistent with the $R_{\text{init}}=30$ au model.</figcaption></figure>

* We show that fitting $C^{17}\text{O}$ hyperfine structure components provides reliable constraints on the CO gas column density distribution in disks. The results are comparable to those derived from thermo-chemical models. The hyperfine analysis is independent of the disk model, or any previous knowledge about the temperature structure.

* Assuming a gas-to-dust mass ratio of 100, we find that all five disks have a CO-to-H$_2$ abundance of $10-100$ times lower than the ISM ratio of $10^{-4}$. The radial distribution of CO gas column density show significant variations in terms of the absolute value and the radial profile among the five disks.

* The MWC 480 and HD 163296 disks show a steep increase of CO gas column density inside their CO snowlines compared to predictions of thermo-chemical models, suggesting their CO abundance inside the CO snowline is significantly higher than the warm molecular layer in the outer disk region. This enhancement of CO abundance inside the CO

snowline is consistent with predictions of a large amount of icy pebble drift into the inner disk regions.

* Four of the five disks have remarkably similar surface density profiles between 150-400 au, which can be well characterized as $\Sigma_{\text{CO}} \propto R^{-2.4}$. We show that the $N_{\text{CO}}$ profiles of all five disks between 150-400 au are consistent with predictions of viscously evolving disks with an initial disk size between 30-100 au.

* We find that when CO gaps are present in the pebble disk region, their locations are correlated with locations of dust gaps. But some of the deepest dust gaps do not have a corresponding CO gap. These CO gas gaps are depleted by a factor of 30-62% of compared to nearby regions, a much shallower depression than the 1-2 orders of depletion seen in dust gaps. The relative depths of CO gaps and mm-continuum gaps are generally consistent with predictions of planet-disk interactions. But
there are also cases where CO gaps are 5-10 times
shallower than predictions based on continuum ob-
servations.

In summary, the MAPS observations show that high
spatial resolution observations of CO isotopologue lines
can provide direct constraints on fundamental disk prop-
erties, such as gas thermal structure and gas mass dis-
tributions. These constraints are key tests of current
theories of planet formation, including dust evolution,
planet-disk interaction, and global evolution of disk
structures. Similar CO observations for a larger sample
of protoplanetary disks will be essential for our under-
standing of typical environment and processes of planet
formation.

ACKNOWLEDGMENTS

This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2018.1.01055.L. ALMA is a part-
nership of ESO (representing its member states),
NSF (USA) and NINS (Japan), together with NRC
(Canada), MOST and ASIAA (Taiwan), and KASI (Re-
public of Korea), in cooperation with the Republic of
Chile. The Joint ALMA Observatory is operated by
ESO, AUI/NRAO and NAOJ. The National Radio As-
tronomy Observatory is a facility of the National Sci-
ence Foundation operated under cooperative agreement
by Associated Universities, Inc.

K.Z. acknowledges the support of the Office of the
Vice Chancellor for Research and Graduate Education
at the University of Wisconsin - Madison with fund-
ing from the Wisconsin Alumni Research Foundation.
K.Z., K.R.S., J.H., J.B., J.B.B., and I.C. acknowl-
edge the support of NASA through Hubble Fellow-
ship grant HST-HF2-51401.001, HST-HF2-51419.001,
HST-HF2-51460.001-A, HST-HF2-51427.001-A, HST-
HF2-51429.001-A, and HST-HF2-51405.001-A awarded
by the Space Telescope Science Institute, which is op-
erated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-
26555. C.J.L. acknowledges funding from the National
Science Foundation Graduate Research Fellowship un-
der Grant DGE1745303. A.D.B., E.A.B., and A.F. ac-
knowledges support from NSF AAG Grant #1907653.
V.V.G. acknowledges support from FONDECYT Ini-
ciación 11180904 and ANID project Basal AFB-170002.
Y.L. acknowledges the financial support by the Natu-
ral Science Foundation of China (Grant No. 11973090)
K.I.Ö acknowledges support from the Simons Foun-
dation (SCOL #321183) and an NSF AAG Grant
(#1907653). S.M.A. and J.H. acknowledge funding sup-
port from the National Aeronautics and Space Admin-

istration under Grant No. 17-XRP17 2-0012 issued
through the Exoplanets Research Program. J.D.I. ac-
knowledges support from the Science and Technology
Facilities Council of the United Kingdom (STFC) un-
der ST/T000287/1. C.W. acknowledges financial sup-
port from the University of Leeds and from the Sci-
ence and Technology Facilities Council (grant numbers
ST/R000549/1 and ST/T000287/1). R.T. and F.L. ac-
knowledges support from the Smithsonian Institution
as a Submillimeter Array (SMA) Fellow. M.L.R.H. ac-
knowledges support from the Michigan Society of Fel-
lows. Y. A. and G. C. acknowledge support by NAOJ
ALMA Scientific Research grant code 2019-13B. Y.A.
acknowledges support by Grant-in-Aid for Scientific Re-
search (S) 18H05222, and Grant-in-Aid for Transfor-
mative Research Areas (A) 20H05844 and 20H05847.
J.K.C. acknowledges support from the National Sci-
ence Foundation Graduate Research Fellowship under
Grant No. DGE 1256260 and the National Aeronau-
tics and Space Administration FINESST grant, under
Grant no. 80NSSC19K1534. A.S. acknowledges support
from ANID/CONICYT Programa de Astronomía Fondo
ALMA-CONICYT 2018 31180052. F.M. acknowledges
support from ANR of France under contract ANR-16-
CE31-0013 (Planet-Forming-Disks) and ANR-15-IDEX-
02 (through CDP "Origins of Life"). R.L.G. acknowl-
edges support from a CNES fellowship grant. H.N.
acknowledges support by NAOJ ALMA Scientific Re-
search grant code 2018-10B and Grant-in-Aid for Scien-
tific Research 18H05441. Y.Y. is supported by IGPEES,
WINGS Program, the University of Tokyo. L.M.P. ac-
knowledges support from ANID project Basal AFB-
170002 and from ANID FONDECYT Iniciación project
#11181068. T.T. is supported by JSPS KAKENHI
Grant Numbers JP17K14244 and JP20K04017.
*Facilities:* ALMA

*Software:* Astropy (Astropy Collaboration et al. 2013), CASA (McMullin et al. 2007), RAC2D (Du &

Bergin 2014), RADMC3D (Dullemond et al. 2012), Gofish (Teague 2019), emcee (Foreman-Mackey et al. 2013), dsharp-opac (Birnstiel et al. 2018).

APPENDIX

<figure><img src="image_23.png" /><figcaption>**Figure 22.** Matched filter responses of the GM Aur, IM Lup and AS 209 disks using a 100 au Keplerian mask. Horizontal dashed lines mark the 3 and 5 σ levels, and the vertical lines mark the frequencies of the C$^{17}$O (1-0) hyperfine structure transitions after correction for the velocities of the sources.</figcaption></figure>

A. ROBUSTNESS TEST OF THE CO RETRIEVAL METHOD

In Section 3.1, we compare a grid of CO models with the observed CO radial intensity profiles to retrieve CO column density distribution. Here we test the robustness of this approach, using the HD 163296 model as an example. We modify the initial CO abundance structure with three radial depletion profiles, a step function, an exponentially decreasing function, and a sinusoidal function (see Figure 25). Using the modified CO abundance structures, we then generate simulated CO observations and compare simulated radial profiles with our grid of CO models to retrieve radial CO depletion profiles. In Figure 25, we show that our method robustly recovered all three input depletion profiles, with an uncertainty less than 20%, consistent with the step size of our logarithmic grid.

REFERENCES

Aikawa, Y., Umebayashi, T., Nakano, T., & Miyama, S. 1998, Faraday Discussions, 109, 281, doi: 10.1039/a800258d

Aikawa, Y., Umebayashi, T., Nakano, T., & Miyama, S. M. 1999, ApJ, 519, 705, doi: 10.1086/307400

Alarcón, F., Teague, R., Zhang, K., Bergin, E. A., & Barraza-Alfaro, M. 2020, ApJ, 905, 68, doi: 10.3847/1538-4357/abc1d6

Alarcón, F., Bosman, A., Bergin, E., et al. 2021, arXiv e-prints, arXiv:2109.06263. https://arxiv.org/abs/2109.06263

Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20, doi: 10.1051/0004-6361/201629929

Anderson, D. E., Bergin, E. A., Maret, S., & Wakelam, V. 2013, ApJ, 779, 141, doi: 10.1088/0004-637X/779/2/141

Andrews, S. M., Wilner, D. J., Espaillat, C., et al. 2011, ApJ, 732, 42, doi: 10.1088/0004-637X/732/1/42

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741

Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46, doi: 10.3847/0004-637X/828/1/46

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068

Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44, doi: 10.3847/1538-4357/aab846
<figure><img src="image_24.png" /><figcaption>**Figure 23.** Deprojected radial intensity profiles of CO lines, for the MAPS sample, ordered by increasing stellar mass from left to right. The same as Figure 2, but on a logarithmic scale to highlight the outer disk emission. The beam sizes are plotted at the upper right corners.</figcaption></figure>

Bai, X.-N. 2016, ApJ, 821, 80,
doi: 10.3847/0004-637X/821/2/80

Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76,
doi: 10.1088/0004-637X/769/1/76

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214,
doi: 10.1086/170270

Banzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJ, 815,
L15, doi: 10.1088/2041-8205/815/1/L15

Bergin, E. A., Aikawa, Y., Blake, G. A., & van Dishoeck,
E. F. 2007, in Protostars and Planets V, ed. B. Reipurth,
D. Jewitt, & K. Keil, 751.
https://arxiv.org/abs/astro-ph/0603358

Bergin, E. A., Cleeves, L. I., Crockett, N., & Blake, G. A. 2014, Faraday Discussions, 168, 61,
doi: 10.1039/C4FD00003J

Bergin, E. A., & Williams, J. P. 2017, The Determination of Protoplanetary Disk Masses, Vol. 445 (Astrophysics and Space Science Library), 1,
doi: 10.1007/978-3-319-60609-5_1

Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature,
493, 644, doi: 10.1038/nature11805

Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, ApJ,
876, 25, doi: 10.3847/1538-4357/ab141e
—. 2020, ApJ, 898, 97, doi: 10.3847/1538-4357/ab9e71

Bergner, J. B., Oberg, K. I., Guzman, V. V., et al. 2021, arXiv e-prints, arXiv:2109.06694.
https://arxiv.org/abs/2109.06694

Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJL,
869, L45, doi: 10.3847/2041-8213/aaf743

Booth, A. S., Walsh, C., Ilee, J. D., et al. 2019, ApJL, 882, L31, doi: 10.3847/2041-8213/ab3645

Booth, R. A., Clarke, C. J., Madhusudhan, N., & Ilee, J. D. 2017, MNRAS, 469, 3994, doi: 10.1093/mnras/stx1103

Booth, R. A., & Ilee, J. D. 2019, MNRAS, 487, 3998,
doi: 10.1093/mnras/stz1488

Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182, doi: 10.1051/0004-6361/201833497

Calahan, J. K., Bergin, E., Zhang, K., et al. 2021a, ApJ,
908, 8, doi: 10.3847/1538-4357/abd255
<figure><img src="image_25.png" /><figcaption>**Figure 24.** Stellar spectra (top row) and dust absorption and scattering coefficients (bottom row) used in thermo-chemical models</figcaption></figure>

<table><caption>**Table 5.** Comparison of the CO gap properties with properties of gas gaps caused by planet-disk interactions</caption><thead><tr><th>Source</th><th>M<sub>*</sub><br>(M<sub>☉</sub>)</th><th>r<sub>CO</sub><br>(au)</th><th>r<sub>mm</sub><br>(au)</th><th>h/r<br>(g cm<sup>-2</sup>)</th><th>Σ<sub>g</sub><br>(au)</th><th>Δ<sub>mm</sub><br>(au)</th><th>δ<sub>CO</sub><br>(au)</th><th>CO<sub>fwhm</sub><br>(m<sub>Jup</sub>)</th><th>m<sub>p4</sub><br>(m<sub>Jup</sub>)</th><th>m<sub>p3</sub><br>(m<sub>Jup</sub>)</th><th>m<sub>p2</sub><br>(m<sub>Jup</sub>)</th><th>depth<br>(au)</th><th>width<br>(au)</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th><th>(10)</th><th>(11)</th><th>(12)</th><th>(13)</th><th>(14)</th></tr></thead><tbody><tr><td>AS 209</td><td>1.2</td><td>94</td><td>99</td><td>0.06</td><td>0.24</td><td>0.31</td><td>0.45</td><td>42.4</td><td>0.03</td><td>0.07</td><td>0.14</td><td>0.86</td><td>16.7</td></tr><tr><td>HD 163296</td><td>2.0</td><td>43</td><td>48</td><td>0.08</td><td>19.17</td><td>0.34</td><td>0.51</td><td>16.5</td><td>1.11</td><td>2.26</td><td>4.61</td><td>0.06</td><td>18.2</td></tr><tr><td>HD 163296</td><td>2.0</td><td>98</td><td>86</td><td>0.08</td><td>9.30</td><td>0.17</td><td>0.58</td><td>37.7</td><td>0.10</td><td>0.20</td><td>0.40</td><td>0.92</td><td>20.0</td></tr><tr><td>HD 163296</td><td>2.0</td><td>145</td><td>143</td><td>0.09</td><td>4.23</td><td>0.09</td><td>0.68</td><td>9.4</td><td>0.005</td><td>0.01</td><td>0.02</td><td>1.00</td><td>12.6</td></tr><tr><td>GM Aur</td><td>1.1</td><td>67</td><td>68</td><td>0.07</td><td>17.06</td><td>0.24</td><td>>0.90</td><td>...</td><td>0.17</td><td>0.34</td><td>0.70</td><td>0.26</td><td>19.3</td></tr><tr><td>MWC 480</td><td>2.1</td><td>63</td><td>73</td><td>0.10</td><td>10.98</td><td>0.29</td><td>0.38</td><td>40.0</td><td>0.66</td><td>1.35</td><td>2.75</td><td>0.37</td><td>22.1</td></tr></tbody></table>

NOTE— (1) Source name; (2) Stellar mass; (3) CO gas gap location; (4) 1.3 mm continuum gap location; (5) Scale height-to-radius ratio at r<sub>mm</sub>, based on model parameters in Table 2; (6) Gas surface density at r<sub>mm</sub>, based on model parameters in Table 2; (7) Δ<sub>mm</sub>: a width parameter measured from mm continuum radial profiles, defined as (r<sub>out</sub>−r<sub>in</sub>)/r<sub>out</sub>. See details in eq. (21) of Zhang et al. (2018); (8) CO gap depth measured in this work, see Table 4. Here δ<sub>CO</sub> = Σ<sub>gap</sub>/Σ<sub>0</sub>; (9) FWHM of CO gap measured in this work, adopted from Table 4; (10) Planet masses derived from {Δ<sub>mm</sub>, h/r, Σ<sub>g</sub>, and M<sub>*</sub>}, using Table 1 of Zhang et al. (2018) and α = 10<sup>−4</sup>; (11) the same as (10), but for α = 10<sup>−3</sup>; (12) the same as (10), but for α = 10<sup>−2</sup>; (13) expected gas gap depth, based on m<sub>p3</sub> and α = 10<sup>−3</sup>, using Kanagawa et al. (2015) eq. (5); (14) expected gas gap width based on m<sub>p3</sub> and α = 10<sup>−3</sup>, using Table 2 of Zhang et al. (2018).
<figure><img src="image_26.png" /><figcaption>**Figure 25.** Robustness test results of the CO abundance retrieval method. Row 1 shows the profiles of CO depletion used in the test. Row 2 shows the 2-dimensional CO abundance after scaling by the profiles. Row 3 and 4 show the result moment zero maps of the ¹³CO (2-1) line images and deprojected radial profiles. Row 5 shows the comparison between the input depletion profile and the retrieved CO depletion profile based on a grid of CO models with constant depletion throughout the whole disk.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Figure 26.** Comparison of radial profiles of the CO isotopologue lines and best-fit results from thermo-chemical models (red solid line, and see model details in Section 3.1.8) and empirical temperature models (blue dash line, and see model details in Section 3.2).</figcaption></figure>

Calahan, J. K., Bergin, E. A., Zhang, K., et al. 2021b, arXiv e-prints, arXiv:2109.06202.
https://arxiv.org/abs/2109.06202

Cieza, L. A., González-Ruilova, C., Hales, A. S., et al. 2021, MNRAS, 501, 2934, doi: 10.1093/mnras/staa3787

Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204, doi: 10.1088/0004-637X/799/2/204

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110, doi: 10.3847/0004-637X/832/2/110
—. 2018, ApJ, 865, 155, doi: 10.3847/1538-4357/aade96

Czekala, I., Loomis, R. A., Teague, R., et al. 2021, arXiv e-prints, arXiv:2109.06188.
https://arxiv.org/abs/2109.06188

Dionatos, O., Woitke, P., Güdel, M., et al. 2019, A&A, 625, A66, doi: 10.1051/0004-6361/201832860

Dong, R., Zhu, Z., & Whitney, B. 2015, ApJ, 809, 93, doi: 10.1088/0004-637X/809/1/93

Du, F., & Bergin, E. A. 2014, ApJ, 792, 2, doi: 10.1088/0004-637X/792/1/2

Du, F., Bergin, E. A., Hogerheijde, M., et al. 2017, ApJ, 842, 98, doi: 10.3847/1538-4357/aa70ee

Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237, doi: 10.1006/icar.1995.1058

Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075, doi: 10.1051/0004-6361:20040284

Dullemond, C. P., Isella, A., Andrews, S. M., Skobleva, I., & Dzyurkevich, N. 2020, A&A, 633, A137, doi: 10.1051/0004-6361/201936438

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library.
http://ascl.net/1202.015

Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A&A, 595, A83, doi: 10.1051/0004-6361/201628509

Espaillat, C., D'Alessio, P., Hernández, J., et al. 2010, ApJ, 717, 441, doi: 10.1088/0004-637X/717/1/441
<figure><img src="image_28.png" /><figcaption>**Figure 27.** Left: radially averaged C¹⁷O spectra (black) and the results from the hyperfine line fitting (blue) for the MWC 480 disk. The center radius of each bin is shown and the spectra are offset along the axis for clarity. Right: results from MCMC fit to the C¹⁷O spectra. From left to right are the C¹⁷O column density, excitation temperature (*T*<sub>ex</sub>) and optical depth (*τ*).</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 28.** Gas and CO column density profiles of viscously evolving disk models from Trapman et al. (2020). These models are for $\alpha = 10^{-3}$ around a 1 M$_{\odot}$ star after 1 Myr of evolution. The initial disk mass in the models is 0.06 M$_{\odot}$ and the initial radii, R$_{init}$, is 30, 50, and 100 au, respectively. These models show that between 150-400 au, the CO column density profile is similar to the gas surface density profile.</figcaption></figure>
Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104, doi: 10.1051/0004-6361/201731390

Fairlamb, J. R., Oudmaijer, R. D., Mendigutía, I., Ilee, J. D., & van den Ancker, M. E. 2015, MNRAS, 453, 976, doi: 10.1093/mnras/stv1576

Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJL, 776, L38, doi: 10.1088/2041-8205/776/2/L38

Favre, C., Fedele, D., Maud, L., et al. 2019, ApJ, 871, 107, doi: 10.3847/1538-4357/aaf80c

Fedele, D., Carney, M., Hogerheijde, M. R., et al. 2017, A&A, 600, A72, doi: 10.1051/0004-6361/201629860

Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24, doi: 10.1051/0004-6361/201731978

Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109, doi: 10.3847/1538-4357/ab8cc5

Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99, doi: 10.1088/0004-637X/813/2/99

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150, doi: 10.3847/1538-4357/aa79f9

Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68, doi: 10.1051/0004-6361/201424693

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051

Guzmán, V. V., Bergner, J. B., Law, C. J., et al. 2021, arXiv e-prints, arXiv:2109.06391. https://arxiv.org/abs/2109.06391

Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135, doi: 10.1146/annurev-astro-081915-023347

Hauschildt, P. H., Allard, F., Ferguson, J., Baron, E., & Alexander, D. R. 1999, ApJ, 525, 871, doi: 10.1086/307954

Hendler, N., Pascucci, I., Pinilla, P., et al. 2020, ApJ, 895, 126, doi: 10.3847/1538-4357/ab70ba

Herczeg, G. J., Linsky, J. L., Valenti, J. A., Johns-Krull, C. M., & Wood, B. E. 2002, ApJ, 572, 310, doi: 10.1086/339731

Herczeg, G. J., Wood, B. E., Linsky, J. L., Valenti, J. A., & Johns-Krull, C. M. 2004, ApJ, 607, 369, doi: 10.1086/383340

Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338, doi: 10.1126/science.1208931

Huang, J., Öberg, K. I., Qi, C., et al. 2017, ApJ, 835, 231, doi: 10.3847/1538-4357/835/2/231

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42, doi: 10.3847/2041-8213/aaf740
—. 2020, ApJ, 891, 48, doi: 10.3847/1538-4357/ab711e

Ingleby, L., Espaillat, C., Calvet, N., et al. 2015, ApJ, 805, 149, doi: 10.1088/0004-637X/805/2/149

Isella, A., Guidi, G., Testi, L., et al. 2016, Physical Review Letters, 117, 251101, doi: 10.1103/PhysRevLett.117.251101

Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJL, 869, L49, doi: 10.3847/2041-8213/aaf747

Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83, doi: 10.1051/0004-6361/201526991

Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2015, ApJL, 806, L15, doi: 10.1088/2041-8205/806/1/L15

Klapper, G., Surin, L., Lewen, F., et al. 2003, The Astrophysical Journal, 582, 262, doi: 10.1086/344615

Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134, doi: 10.3847/1538-4357/aba75d

Krijt, S., Ciesla, F. J., & Bergin, E. A. 2016, ApJ, 833, 285, doi: 10.3847/1538-4357/833/2/285

Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78, doi: 10.3847/1538-4357/aad69b

Lacy, J. H., Knacke, R., Geballe, T. R., & Tokunaga, A. T. 1994, ApJL, 428, L69, doi: 10.1086/187395

Law, C. J., Loomis, R. A., Teague, R., et al. 2021a, arXiv e-prints, arXiv:2109.06210. https://arxiv.org/abs/2109.06210

Law, C. J., Teague, R., Loomis, R. A., et al. 2021b, arXiv e-prints, arXiv:2109.06217. https://arxiv.org/abs/2109.06217

Liu, Y., Dipierro, G., Ragusa, E., et al. 2019, A&A, 622, A75, doi: 10.1051/0004-6361/201834157

Long, F., Herczeg, G. J., Pascucci, I., et al. 2017, ApJ, 844, 99, doi: 10.3847/1538-4357/aa78fc

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17, doi: 10.3847/1538-4357/aae8e1

Loomis, R. A., Öberg, K. I., Andrews, S. M., et al. 2018, AJ, 155, 182, doi: 10.3847/1538-3881/aab604
—. 2020, ApJ, 893, 101, doi: 10.3847/1538-4357/ab7cc8

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603, doi: 10.1093/mnras/168.3.603

Macías, E., Espaillat, C. C., Ribas, Á., et al. 2018, ApJ, 865, 37, doi: 10.3847/1538-4357/aad811

Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3, doi: 10.1051/0004-6361/201628549

Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266, doi: 10.1086/680323

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591

McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167, doi: 10.3847/0004-637X/831/2/167
McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127

Mendigutía, I., Brittain, S., Eiroa, C., et al. 2013, ApJ, 776, 44, doi: 10.1088/0004-637X/776/1/44

Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96, doi: 10.1051/0004-6361/201424712

Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2017, A&A, 599, A113, doi: 10.1051/0004-6361/201629556

Molyarova, T., Akimkin, V., Semenov, D., et al. 2017, ApJ, 849, 130, doi: 10.3847/1538-4357/aa9227

Montesinos, B., Eiroa, C., Mora, A., & Merín, B. 2009, A&A, 495, 901, doi: 10.1051/0004-6361:200810623

Morbidelli, A., & Raymond, S. N. 2016, Journal of Geophysical Research (Planets), 121, 1962, doi: 10.1002/2016JE005088

Mulders, G. D., Pascucci, I., Manara, C. F., et al. 2017, ApJ, 847, 31, doi: 10.3847/1538-4357/aa8906

Najita, J. R., & Bergin, E. A. 2018, ApJ, 864, 168, doi: 10.3847/1538-4357/aad80c

Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375, doi: 10.1016/0019-1035(86)90121-1

Nomura, H., Tsukagoshi, T., Kawabe, R., et al. 2016, ApJL, 819, L7, doi: 10.3847/2041-8205/819/1/L7

Oberg, K. I., Guzman, V. V., Walsh, C., et al. 2021, arXiv e-prints, arXiv:2109.06268. https://arxiv.org/abs/2109.06268

Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82, doi: 10.3847/0004-637X/821/2/82

Piétu, V., Dutrey, A., & Guilloteau, S. 2007, A&A, 467, 163, doi: 10.1051/0004-6361:20066537

Pinilla, P., Pohl, A., Stammler, S. M., & Birnstiel, T. 2017, ApJ, 845, 68, doi: 10.3847/1538-4357/aa7edb

Pinte, C., Ménard, F., Duchêne, G., et al. 2018a, A&A, 609, A47, doi: 10.1051/0004-6361/201731377

Pinte, C., Price, D. J., Ménard, F., et al. 2018b, ApJ, 860, L13, doi: 10.3847/2041-8213/aac6dc
—. 2020, ApJL, 890, L9, doi: 10.3847/2041-8213/ab6dda

Qi, C., D'Alessio, P., Oberg, K. I., et al. 2011, ApJ, 740, 84

Qi, C., Wilner, D. J., Calvet, N., et al. 2006, ApJL, 636, L157, doi: 10.1086/500241

Rab, C., Kamp, I., Dominik, C., et al. 2020, A&A, 642, A165, doi: 10.1051/0004-6361/202038712

Riols, A., & Lesur, G. 2019, A&A, 625, A108, doi: 10.1051/0004-6361/201834813

Ros, K., & Johansen, A. 2013, A&A, 552, A137, doi: 10.1051/0004-6361/201220536

Ros, K., Johansen, A., Riipinen, I., & Schlesinger, D. 2019, A&A, 629, A65, doi: 10.1051/0004-6361/201834331

Rosotti, G. P., Ilee, J. D., Facchini, S., et al. 2021, MNRAS, 501, 3427, doi: 10.1093/mnras/staa3869

Rosotti, G. P., Juhasz, A., Booth, R. A., & Clarke, C. J. 2016, MNRAS, 459, 2790, doi: 10.1093/mnras/stw691

Salyk, C., Herczeg, G. J., Brown, J. M., et al. 2013, ApJ, 769, 21, doi: 10.1088/0004-637X/769/1/21

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91, doi: 10.3847/0004-637X/823/2/91
—. 2018, ApJ, 856, 85, doi: 10.3847/1538-4357/aaae08

—. 2019, ApJ, 877, 131, doi: 10.3847/1538-4357/ab1c5e

Schwarz, K. R., Calahan, J. K., Zhang, K., et al. 2021, arXiv e-prints, arXiv:2109.06228. https://arxiv.org/abs/2109.06228

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, arXiv e-prints, arXiv:2109.06433. https://arxiv.org/abs/2109.06433

Simon, M., Guilloteau, S., Beck, T. L., et al. 2019, ApJ, 884, 42, doi: 10.3847/1538-4357/ab3e3b

Smirnov-Pinchukov, G. V., Semenov, D. A., Akimkin, V. V., & Henning, T. 2020, A&A, 644, A4, doi: 10.1051/0004-6361/202038572

Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2018, MNRAS, 477, 1239, doi: 10.1093/mnras/sty717

Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74, doi: 10.1051/0004-6361/201628955

Teague, R. 2019, The Journal of Open Source Software, 4, 1632, doi: 10.21105/joss.01632

Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378, doi: 10.1038/s41586-019-1642-0

Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018a, ApJ, 860, L12, doi: 10.3847/2041-8213/aac6d7

Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018b, ApJ, 868, 113, doi: 10.3847/1538-4357/aae836

Teague, R., & Loomis, R. 2020, ApJ, 899, 157, doi: 10.3847/1538-4357/aba956

Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49, doi: 10.1051/0004-6361/201628550

Teague, R., Bae, J., Aikawa, Y., et al. 2021, arXiv e-prints, arXiv:2109.06218. https://arxiv.org/abs/2109.06218

Trapman, L., Miotello, A., Kama, M., van Dishoeck, E. F., & Bruderer, S. 2017, A&A, 605, A69, doi: 10.1051/0004-6361/201630308
Trapman, L., Rosotti, G., Bosman, A. D., Hogerheijde, M. R., & van Dishoeck, E. F. 2020, A&A, 640, A5, doi: 10.1051/0004-6361/202037673

Umebayashi, T., & Nakano, T. 1981, Publications of the Astronomical Society of Japan, 33, 617

van der Marel, N., van Dishoeck, E. F., Bruderer, S., Pérez, L., & Isella, A. 2015, A&A, 579, A106, doi: 10.1051/0004-6361/201525658

van der Marel, N., Williams, J. P., & Bruderer, S. 2018, ApJL, 867, L14, doi: 10.3847/2041-8213/aae88e

van Dishoeck, E. F., & Black, J. H. 1988, ApJ, 334, 771, doi: 10.1086/166877

Visser, R., van Dishoeck, E. F., Doty, S. D., & Dullemond, C. P. 2009, A&A, 495, 881, doi: 10.1051/0004-6361/200810846

Walsh, C., Millar, T. J., Nomura, H., et al. 2014, A&A, 563, A33, doi: 10.1051/0004-6361/201322446

Weaver, E., Isella, A., & Boehler, Y. 2018, ApJ, 853, 113, doi: 10.3847/1538-4357/aaa481

Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59, doi: 10.1088/0004-637X/788/1/59

Wilson, T. L. 1999, Reports on Progress in Physics, 62, 143, doi: 10.1088/0034-4885/62/2/002

Xu, R., Bai, X.-N., & Öberg, K. 2017, ApJ, 835, 162, doi: 10.3847/1538-4357/835/2/162

Yu, M., Willacy, K., Dodson-Robinson, S. E., Turner, N. J., & Evans, II, N. J. 2016, ApJ, 822, 53, doi: 10.3847/0004-637X/822/1/53

Zhang, K., Bergin, E. A., Blake, G. A., et al. 2016, ApJL, 818, L16, doi: 10.3847/2041-8205/818/1/L16

Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, Nature Astronomy, 1, 0130, doi: 10.1038/s41550-017-0130

Zhang, K., Bergin, E. A., Schwarz, K., Krijt, S., & Ciesla, F. 2019, ApJ, 883, 98, doi: 10.3847/1538-4357/ab38b9

Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJL, 806, L7, doi: 10.1088/2041-8205/806/1/L7

Zhang, K., Bosman, A. D., & Bergin, E. A. 2020a, ApJL, 891, L16, doi: 10.3847/2041-8213/ab77ca

Zhang, K., Schwarz, K. R., & Bergin, E. A. 2020b, ApJL, 891, L17, doi: 10.3847/2041-8213/ab7823

Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47, doi: 10.3847/2041-8213/aaf744