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; Chunhua Qi; Merel L. R. van ’t Hoff; Feng Long; David J. Wilner; Jane Huang; Ian Czekala; John D. Ilee; Gianni Cataldi; Jennifer B. Bergner; Yuri Aikawa; Richard Teague; Jaehan Bae; Ryan A. Loomis; Jenny K. Calahan; Felipe Alarcón; François Ménard; Romane Le Gal; Anibal Sierra; Yoshihide Yamato; Hideko Nomura; Takashi Tsukagoshi; Laura M. Pérez; Leon Trapman; Yao Liu; Kenji Furuya

doi:10.3847/1538-4365/ac1580

1. Introduction

In protoplanetary disks, H₂ is the primary mass carrier, but it does not have strong transitions in the temperature range in which 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 2018). CO is also chemically stable (van Dishoeck & Black 1988; Visser et al. 2009) and stays in the gas phase at temperatures 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)millimeter CO (J ≤ 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; Long et al. 2017; Miotello 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, 2011; Piétu et al. 2007; Kama et al. 2016; Schwarz et al. 2016; Pinte et al. 2018a), kinematics in disks (Pinte et al. 2018b; Teague et al. 2018a, 2019), and turbulence levels (Flaherty et al. 2015, 2017, 2020; Teague et al. 2016). 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 has been 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 interstellar medium (ISM) ratio of 10⁻⁴ in the warm molecular layer throughout the entire disk lifetime (Bergin & Williams 2018). 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 ≤10 in the majority of disks (Ansdell et al. 2016; Long et al. 2017; Miotello et al. 2017). In contrast, the high gas accretion rates onto the central star suggest that 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. 2021).

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; Bosman et al. 2018; Schwarz 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 toward the midplane (Krijt et al. 2016, 2018, 2020; Xu et al. 2017). 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 ∼1 Myr (Bergner et al. 2020; Zhang et al. 2020b). 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)millimeter continuum emission of protoplanetary disks (e.g., Isella et al. 2016; Zhang 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 property 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 millimeter-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 previous studies of the CO gas mass distributions in Class II disks have relatively low spatial resolutions, between 25 and 100 au (e.g., Isella et al. 2016; Nomura et al. 2016; Schwarz et al. 2016; Fedele et al. 2017; Favre et al. 2019; Zhang et al. 2019; Rosotti et al. 2021). Although local coincidences in CO and dust gap locations have been reported in a few disks, previous observations 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 radius? (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 density 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 transitions (Öberg et al. 2021). The analysis presented here is based on the MAPS fiducial images, which have 0 farcs 15 and 0 farcs 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 observational logs, line frequencies, and image statistics are listed in Öberg et al. (2021). We use the integrated intensity 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.

**Figure 1.** Integrated intensity maps of CO,¹³CO, C¹⁸O, and C¹⁷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 in the lower right corner of each panel represents 50 au scale.
Download figure:
Standard image High-resolution image

**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 in the upper right corner of each panel. The beam sizes are 015 × 015 for all (2–1) lines and 03 × 03 for (1–0) lines. A plot of the intensity profiles on a logarithmic scale is shown in Appendix B in Figure 24.
Download figure:
Standard image High-resolution image

farcs — **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 in the upper right corner of each panel. The beam sizes are 015 × 015 for all (2–1) lines and 03 × 03 for (1–0) lines. A plot of the intensity profiles on a logarithmic scale is shown in Appendix B in Figure 24.
Download figure:
Standard image High-resolution image

We also present the C¹⁷O (1–0) data that are not covered by Law et al. (2021a). The C¹⁷O (1–0) line has three hyperfine components that are blended in the integrated intensity map. As shown in Figures 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 Appendix B in Figure 23. 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 column density varies with radial distance from the central star in the five MAPS disks (N_CO(R)), using spatially resolved CO isotopologue line images.

To test the robustness of the results, we employ three independent methods to constrain N_CO(R) profiles:

Method 1: we build thermochemical models for individual disks to constrain the temperature structures of gas and dust in these disks. These models are based on fitting broadband spectral energy distribution (SED), (sub)millimeter continuum images, and vertical locations of CO emission surfaces. We then use the thermal structure to retrieve the N_CO 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.

**Figure 3.** An outline of the modeling processes. Codes used are in light blue, input data are highlighted in filled blue boxes, and the output of each step is highlighted in red.
Download figure:
Standard image High-resolution image

Method 2: we replace the gas temperature in thermochemical models with empirically constructed temperature structures based on optically thick CO lines, and then we redo the N_CO 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 optical depth and excitation temperature. This analysis is described in Section 3.3.

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 at www.alma-maps.info.

3.1. Themochemical 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):

$\begin{eqnarray}&&{\rm{\Sigma }}(R)={{\rm{\Sigma }}}_{c}{\left(\frac{R}{{R}_{c}}\right)}^{-\gamma }\exp \left[-{\left(\frac{R}{{R}_{c}}\right)}^{2-\gamma }\right],\end{eqnarray} \tag{ 1 }$

where R_c is the characteristic scaling radius and γ is the gas surface density exponent.

For the large dust population (pebbles), its surface density is set by matching millimeter 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:

$\begin{eqnarray}&&{\rho }_{i}(R,Z)={f}_{i}\displaystyle \frac{{\rm{\Sigma }}(R)}{\sqrt{2\pi }{H}_{i}(R)}\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{Z}{{H}_{i}(R)}\right)}^{2}\right],\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&{H}_{i}(R)={\chi }_{i}{H}_{100}{\left(R/100\,\mathrm{au}\right)}^{\psi },\end{eqnarray} \tag{ 3 }$

where f_i is the mass fraction of each mass component, H₁₀₀ is the scale height at 100 au, and ψ is a parameter that characterizes the radial dependence of the scale height. The large-grain population is expected to be more settled compared to the gas and small grains (Nakagawa et al. 1986; Dullemond & Dominik 2004). This settling effect is characterized with the parameter χ_i. Here we fix χ = 1 for the gas and the small-grain population and χ = 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) ∝ 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 multiwavelength 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 and 3 mm beyond 25 au.

For the composition of the large dust population, we adopt the recommendation of the DSHARP collaboration (Birnstiel et al. 2018), which is a mixture of water ice, astronomical silicates, troilite, and refractory organics. 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 regions (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 wavelength dependence of dust opacity. We show the dust opacity profiles in Appendix B in Figure 25.

3.1.3. Stellar Spectrum

Young stars accrete disk materials actively and produce strong X-ray and UV radiation (Hartmann et al. 2016). This high-energy radiation is particularly important for chemistry in disks (Bergin et al. 2007). Here we make a composite spectrum for each source by combining a model stellar photosphere spectrum with UV and X-ray spectra from observations.

The stellar photospheric spectra are computed by interpolating the standard Nextgen spectral model library in the (log₁₀ T_eff, log₁₀ g) space, assuming a solar metallicity (Hauschildt et al. 1999). The UV and X-ray spectra 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 circumstellar materials and therefore do not represent the UV spectra received by the disk. We apply the extinction at UV wavelength based on A_V 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 Appendix B in Figure 25.

3.1.4. Thermochemical Modeling

We use the thermochemical 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 the 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 simultaneously solves the chemical evolution and gas thermal structure in the disk. The chemical network has 484 species and 5830 reactions, including the full gas-phase network from the UMIST 2006 database,²⁶ dissociation 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⁻¹ at the disk surface, because previous studies suggest that 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/¹³CO = 69, CO/C¹⁸O = 557, CO/C¹⁷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 millimeter 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).

Table 1. Stellar and Disk Properties

Source	SpT	d	Incl.	PA	T_eff	L_*	M_*	log₁₀( $\dot{M}$ )	v_sys	References
		(pc)	(deg)	(deg)	(K)	(L_⊙)	(M_⊙)	(M_⊙ yr⁻¹)	(km s⁻¹)
IM Lup	K5	158	47.5	144.5	4266	2.57	1.1	−7.9	4.5	1, 2, 3, 4, 5
GM Aur	K6	159	53.2	57.2	4350	1.2	1.1	−8.1	5.6	1, 6, 7, 8, 9
AS 209	K5	121	35.0	85.8	4266	1.41	1.2	−7.3	4.6	1, 2, 5, 10, 11
HD 163296	A1	101	46.7	133.3	9332	17.0	2.0	−7.4	5.8	1, 2, 5, 12, 13
MWC 480	A5	162	37.0	148.0	8250	21.9	2.1	−6.9	5.1	1, 14, 15, 16, 17, 18

Note. Reproduced from Öberg et al. (2021), where further details about the source characteristics are available.

References. (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.

Download table as: ASCII Typeset image

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⁻¹ spectral resolution. The model channels are rebinned to the spectral resolution of MAPS observations (0.2 km s⁻¹ and 0.5 km s⁻¹ for the (2-1) and (1-0) lines, respectively) 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 six free parameters: four parameters to characterize the surface density profiles of gas and small dust, $\left\{{{\rm{\Sigma }}}_{c}^{g},{{\rm{\Sigma }}}_{c}^{\mu {\rm{m}}},{R}_{c},\gamma \right\}$ , and two parameters to describe the scale height distribution $\left\{{H}_{100},\psi \right\}$ . Moreover, the surface density profile of the large-grain population is determined by an iterative process to fit the millimeter continuum images (more details below).

The initial values of parameters are taken from literature for each disk (Cleeves et al. 2016; Fedele et al. 2018; Macías et al. 2018; Liu et al. 2019; Zhang et al. 2019). We then update the parameters by comparing model outputs with three types of observations: SEDs, millimeter 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 ψ and H₁₀₀ parameters. Therefore, we run a small grid of ψ and H₁₀₀ 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.

**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), those of GM Aur from Macías et al. (2018), and those of MWC 480 from Anderson et al. (2013).
Download figure:
Standard image High-resolution image

The second type of observational constraints are the spatially resolved millimeter 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 images of IM Lup, AS 209, and HD 163296 are from the DSHARP program (Andrews et al. 2018; Huang et al. 2018), the MWC 480 image is from Long et al. (2018) and the GM Aur image is 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 three to four 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. This suggests that our model underpredicts the midplane temperature inside 10 au, and the real dust scale height may deviate from the general power-law function we used to match SEDs. Given that 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 Figures 5 and 6, respectively.

**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 Equation (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).
Download figure:
Standard image High-resolution image

**Figure 6.** Comparison between the observed 1.3 mm continuum images (top row) and simulated images from our thermochemical 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 is from Long et al. (2018) and Liu et al. (2019); and the GM Aur disk image is 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.
Download figure:
Standard image High-resolution image

The third type of observational constraint are 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 τ ∼ 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.

**Figure 7.** Gas temperature (top row) and CO abundance structure (bottom row) of the best-fit thermochemical 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.
Download figure:
Standard image High-resolution image

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 thermochemical 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 freezeout 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 thermochemical models provides a zeroth-order check on whether the models reasonably match with the observations. For example, the CO freezeout 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 ≃ 20 K layer where CO freezes out. We start with H₁₀₀ and ψ parameters that match the SEDs. The large dust mass distribution is then constrained by millimeter 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 thermochemical models with ψ and H₁₀₀ 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 parameter 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 overpredicts 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.

Table 2. Best-fit Disk Model Parameters

Source	M_g	${{\rm{\Sigma }}}_{c}^{g}$	M_mm	M_μm	${{\rm{\Sigma }}}_{c}^{\mu {\rm{m}}}$	${R}_{c}^{g}$	γ_g	H₁₀₀	ψ	r_in	r_out
	(M_⊙)	(g cm⁻²)	(M_⊙)	(M_⊙)	(g cm⁻²)	(au)		(au)		(au)	(au)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
IM Lup	0.2	28.4	1.97e−03	2.02e−05	2.9e−03	100	1.0	10.0	1.17	0.20	1200
GM Aur	0.2	9.4	5.94e−04	1.03e−04	6.0e−03	176	1.0	7.5	1.35	1.00	650
AS 209	0.0045	1.0	4.50e−04	5.23e−05	1.2e−02	80	1.0	6.0	1.25	0.50	300
HD 163296	0.14	8.8	2.31e−03	2.00e−04	1.3e−02	165	0.8	8.4	1.08	0.45	600
MWC 480	0.16	5.8	1.41e−03	1.69e−04	6.1e−03	200	1.0	10.0	1.08	0.45	750

Note. Column (1): source name. Column (2): total gas mass of the disk. Column (3): surface density of gas at ${R}_{c}^{g}$ . Column (4): total mass of the large-grain population. Column (5): total mass of the small-grain population. Column (6): surface density of the small-grain population at ${R}_{c}^{g}$ . Column (7): characteristic radius of the gas and the small-grain population. Column (8): surface density exponent in Equation (1). Column (9): scale height at 100 au. Column (10): power index of the radial dependence of scale height, Equation (3). Column (11): the inner edge of the disk model. Column (12): the outer edge of the disk model.

Download table as: ASCII Typeset image

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 and 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 thermochemical model. A robustness test of this method is presented in Appendix A. For a consistency check, we also apply the best-fitting 1D depletion profiles to CO abundance structures in thermochemical models and then generate simulated CO line radial profiles. Figure 26 in Appendix B shows that our best-fitting models match well with the observed profiles.

**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 gray 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.
Download figure:
Standard image High-resolution image

3.1.8. CO Column Density from Thermochemical 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 5. 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 and 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 and 120 au, and shows a deeper gap around 240 au. The AS 209 disk has the smallest CO gas disk with a steep decrease between 150 and 300 au. Interestingly, the N_CO of the remaining four disks show a very similar shallow slope between 150 and 400 au; see Figure 10. The slope can be characterized by a power-law function of Σ(R)∝ 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 those 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 thermochemical 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 that the emission fills the beam (Isella et al. 2018; Pinte et al. 2018a; Weaver 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 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 thermochemical models and empirical temperature structure from Law et al. (2021b). In general, thermochemical 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 thermochemical 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 thermochemical model is generally colder by 5–10 K than the empirical temperature structure. A caveat of our thermochemical 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. (2021) and Schwarz et al. (2021) presented thermochemical 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 the CO column density. The dust temperature and CO abundance structure are kept the same as the thermochemical models. Figure 26 in Appendix B 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 thermochemical 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 thermochemical 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¹⁷O Hyperfine Line Fitting

The third method we use to derive N_CO is fitting C¹⁷O hyperfine structure transitions. This method is independent of the temperature and CO abundance structures adopted in the thermochemical models.

The C¹⁷O J = 1–0 rotational transition is split into three components that can be fitted to obtain constraints on the C¹⁷O column density, excitation temperature, and optical depth (Mangum & Shirley 2015). This method has been used recently in disk studies of C₂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¹⁷O in protoplanetary disks. For our sample, only the C¹⁷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 the C¹⁷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 Öberg et al. (2021). The 0.5 km s⁻¹ spectral resolution of the MAPS Band 3 data means that two of the components are always blended but the third component is ∼2 km s⁻¹ offset 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 ratio 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 model spectra for the C¹⁷O (1–0) hyperfine transitions to fit the observed spectra. We assume LTE and that all three hyperfine 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 (σ_obs, km s⁻¹), and systematic velocity (V_lsrk, km s⁻¹). The σ_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; and (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²² cm⁻² and 10–50 K. This encompasses the range of temperature for the C¹⁷O-emitting region (see Section 3.2). σ_obs is set to range from 0.01 to 8 km s⁻¹, and 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 Markov Chain Monte Carlo (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 τ for both disks are shown in Appendix B in Figure 27. The C¹⁷O line becomes optically thick (τ > 1) within ≈55 au in the HD 163296 disk and within ≈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.

**Figure 13.** Radially averaged C¹⁷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.
Download figure:
Standard image High-resolution image

The N(C¹⁷O) measured from the hyperfine line fitting method has inner and outer limits of the radius. Inside the inner 1.5× beams, the fit is poor owing 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 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 and 150 au than that of the 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 ratio 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 C¹⁸O lines match well in the MWC 480 disk. Both disks have previous detections of rarer CO isotopologue lines, including ¹³C¹⁸O and ¹³C¹⁷O lines (Booth et al. 2019; Loomis et al. 2020; Zhang et al. 2020a), which provides constraints closer to the midplane in the two disks. Here we use the disk-integrated ¹³C¹⁸O line spectra to test the robustness of the N_CO results inside 100 au for the HD 163296 and MWC 480 disks. We apply the CO depletion profiles from the C¹⁸O and C¹⁷O lines to the CO abundance from our thermochemical models and generate simulated spectra for ¹³C¹⁸O (2–1) or (3–2). The results are shown in Figure 15.

**Figure 15.** Comparison of ¹³C¹⁸O spectra (in gray) with models based on N_CO from C¹⁸O(2–1) and C¹⁷O(1–0) lines in thermochemical models (Section 3.1).
Download figure:
Standard image High-resolution image

For HD 163296, the N_CO from the C¹⁷O (1–0) line has a better match to the ¹³C¹⁸O (2–1) line spectrum. The high N_CO 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 a few larger than that in the outer disk region. As the ¹³C¹⁸O and C¹⁸O (2–1) lines are at similar frequencies, the difference of N_CO cannot be explained by optically thick dust. One possible explanation is that the CO-to-H₂ gas abundance ratio close to the midplane may be a factor of a few higher than that of the disk atmosphere inside 70 au. This enhanced CO abundance at the midplane is possible if a large amount of icy pebbles evaporate their CO ice mantle at the midplane (Krijt et al. 2018, 2020). For the MWC 480 disk, the N_CO discrepancy inside 100 au is smaller, and both results match the ¹³C¹⁸O (3–2) line within uncertainties. Higher spatial resolution observations of ¹³C¹⁸O are necessary to better constrain the N_CO 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₂ abundance ratio of 10⁻⁴ to get a total disk gas mass, M_gas1. This is a lower limit of disk gas mass, as it does not include gas mass in the CO freezeout 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_gas2. This approach assumes that the difference between thermochemical 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 freezeout and photodissociation effects, but the result depends on the temperature and CO abundance structure in models. For the mass estimation, we use the CO depletion profiles based on the C¹⁸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_CO from the C¹⁷O (1–0) line, due to the discrepancy discussed in Section 3.4. We note that our initial chemical condition in thermochemical models has a CO-to-H₂ abundance of 2.8 × 10⁻⁴ (Lacy et al. 1994), and thus the undepleted CO abundance in the warm molecular layer is between 1 and 2.8 × 10⁻⁴, slightly higher than the canonic assumption of 10⁻⁴. Due to the difference in assumed CO-to-H₂ abundance ratios, even M_gas2 includes corrections of freezeout and photodissociation; in the IM Lup and MWC 480 disks, the M_gas2 are still slightly smaller than or the same as M_gas1.

The estimations of minimum gas disk masses are summarized in Table 3. The total dust masses are adopted from our best-fit thermochemical 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 ratio ≥ 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.

Table 3. Minimum Disk Mass and Snowline Location

Source	M_gas1	M_gas2	g2d1	g2d2	r_CO	${r}_{{{\rm{N}}}_{2}}$
	(10⁻³ M_⊙)	(10⁻³ M_⊙)			(au)	(au)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
IM Lup	1.9	1.9	1	1	15 ± 5	25 ± 5
GM Aur	8.0	19.9	12	29	30 ± 5	70 ± 5
AS 209	0.4	1.0	1	2	12 ± 5	20 ± 5
HD 163296	31.3	42.7	14	19	65 ± 5	90 ± 5
MWC 480	35.9	28.6	25	20	100 ± 5	150 ± 5

Note. Column (1): source name. Column (2): M_gas1, the minimum disk gas masses by scaling the total CO gas mass with n_CO/ ${n}_{{{\rm{H}}}_{2}}$ = 10⁻⁴. Column (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. Column (4): minimum gas-to-dust mass ratio based on M_gas1/M_dust. Column (5): M_gas2/M_dust. Column (6): midplane CO snowline location in our thermochemical model. Column (7): midplane N₂ snowline location.

Download table as: ASCII Typeset image

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 that this 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 thermochemical models. Compared to our thermochemical 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 (McClure et al. 2016; Schwarz 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; Bosman et al. 2018; Schwarz et al. 2018). Additionally, dust growth and settling can sequester CO into the midplane, 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¹⁸O (2–1) or (3–2) line observations. They found that the CO gas abundance is heavily depleted in the region just beyond the midplane 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 ¹³C¹⁸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 midplane CO snowline, CO is depleted by a factor of 10 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 midplane as icy dust grains settle, and icy pebbles drift inward into the region inside the midplane CO snowline, where CO ice evaporates and enhances the CO abundance at the inner disk regions (Booth et al. 2017; Krijt et al. 2018, 2020; Booth & Ilee 2019).

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 farther out than the midplane 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 midplane 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 versus Millimeter-sized Dust

4.3.1. CO Gap Properties

In Figure 18, we show the CO gas profiles and millimeter-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 that the AS 209 disk has a deep gap (∼97% depletion) around 240 au and GM Aur has a deep inner cavity inside 40 au.²⁷ Except for these features, other N_CO substructures are subtler 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.

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

Table 4. Properties of CO Gaps

Source	r_gap	σ	Depth
	(au)	(au)
(1)	(2)	(3)	(4)
GM Aur	14 ± 1	14 ± 1	0.95 ± 0.03
AS 209	59 ± 9	13 ± 6	0.47 ± 0.27
	92 ± 12	18 ± 8	0.55 ± 0.12
	238 ± 1	13 ± 1	0.97 ± 0.08
HD 163296	44 ± 1	7 ± 1	0.49 ± 0.04
	93 ± 1	16 ± 1	0.42 ± 0.03
	148 ± 1	4 ± 1	0.32 ± 0.05
MWC 480	63 ± 1	17 ± 1	0.62 ± 0.03

Note. Column (1): source name. Column (2): CO gap location. Column (3): gap width derived by fitting a Gaussian. σ is the standard deviation of the best-fitting Gaussian function. Column (4): gap depth, defined as (Σ_base − Σ_gap)/Σ_base.

Download table as: ASCII Typeset image

The AS 209 disk shows a broad CO gap between 40 and 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 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_CO is 38%–68% of that at nearby regions. The deepest CO gap 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.

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 and 4 M_J 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). Hydrodynamical models showed that three giant planets at 54, 83, and 136 au provide the best match to the velocity perturbations.²⁸ In Figure 19, we compare our N_CO 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 the 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_J 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.

**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¹⁸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.
Download figure:
Standard image High-resolution image

For the AS 209 disk, we identify two CO gaps. The inner gap between 50 and 120 au coincides with two millimeter 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_J) at 99 au. Favre et al. (2019) showed that a planet of 0.2–0.3 M_J at 100 au is roughly consistent with their C¹⁸O and ¹³CO (2–1) observations. However, using higher spatial resolution CO and C₂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 a > 0.2 M_J 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 the 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 millimeter 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⁻³). We then use the inferred planet masses and Equations (5)–(6) in Kanagawa et al. (2015) to derive the expected gas gap depths.²⁹ We use Equation (22) in Zhang et al. (2018) to derive the expected gas gap width. The continuum gap widths, disk parameters, and estimated planet masses are listed 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.

Table 5. Comparison of the CO Gap Properties with Properties of Gas Gaps Caused by Planet–Disk Interactions

Source	M_⋆	r_CO	r_mm	h/r	Σ_g	Δ_mm	δ_CO	CO_fwhm	m_p4	m_p3	m_p2	Depth	Width
	(M_⊙)	(au)	(au)		(g cm⁻²)			(au)	(m_Jup)	(m_Jup)	(m_Jup)		(au)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
AS 209	1.2	94	99	0.06	0.24	0.31	0.45	42.4	0.03	0.07	0.14	0.86	16.7
HD 163296	2.0	43	48	0.08	19.17	0.34	0.51	16.5	1.11	2.26	4.61	0.06	18.2
HD 163296	2.0	98	86	0.08	9.30	0.17	0.58	37.7	0.10	0.20	0.40	0.92	20.0
HD 163296	2.0	145	143	0.09	4.23	0.09	0.68	9.4	0.005	0.01	0.02	1.00	12.6
GM Aur	1.1	67	68	0.07	17.06	0.24	>0.90	⋯	0.17	0.34	0.70	0.26	19.3
MWC 480	2.1	63	73	0.10	10.98	0.29	0.38	40.0	0.66	1.35	2.75	0.37	22.1

Note. Column (1): source name. Column (2): stellar mass. Column (3): CO gas gap location. Column (4): 1.3 mm continuum gap location. Column (5): scale height-to-radius ratio at r_mm, based on model parameters in Table 2. Column (6): gas surface density at r_mm, based on model parameters in Table 2. Column (7): Δ_mm, a width parameter measured from millimeter continuum radial profiles, defined as (r_out–r_in)/r_out. See details in Equation (21) of Zhang et al. (2018). Column (8): CO gap depth measured in this work; see Table 4. Here δ_CO = Σ_gap/Σ₀. Column (9): FWHM of CO gap measured in this work, adopted from Table 4. Column (10): planet masses derived from {Δ_mm, h/r, Σ_g, and M_⋆}, using Table 1 of Zhang et al. (2018) and α = 10⁻⁴. Column (11): the same as Column (10), but for α = 10⁻³. Column (12): the same as Column (10), but for α = 10⁻². Column (13): expected gas gap depth, based on m_p3 and α = 10⁻³, using Kanagawa et al. (2015), Equation (5). Column (14): expected gas gap width based on m_p3 and α = 10⁻³, using Table 2 of Zhang et al. (2018).

Download table as: ASCII Typeset image

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; Alarcón et al. 2021; Teague 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 thermochemical 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 occurs 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 ratios. If the transition zone takes a significant fraction of the total CO gas column, the N_CO derived from C¹⁸O and an ISM element isotopic ratio would be lower than the true N_CO (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 that 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 the ¹³CO/C¹⁸O ratio 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 significantly. However, as shown in Figure 18, the C¹⁸O depletions inside gaps are only modest or even not detected inside the millimeter continuum gap locations. It is unlikely that isotope-selective photodissociation is much stronger inside these gaps compared to the nearby regions. Further 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 affects 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 regions 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 need not grow with time. Global disk properties, such as accretion rates and disk sizes, have been tried to constrain viscous disk evolution models, but observations 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_CO 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_CO distributions between 150 and 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_CO ∝ r^−2.4±0.2. This N_CO 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_CO profiles with those from thermochemical models of viscous disks of Trapman et al. (2020). We adopt N_CO profiles from three viscous disk models with α = 10⁻³ around a 1 M_⊙ star after 1 Myr of evolution. These models have an initial disk mass of 0.06 M_⊙, and their initial disk size, R_init, is 30, 50, and 100 au, respectively. These models predict that the slope of CO column density distribution is similar to the gas in the region between 100 and 500 au (see Figure 28 in Appendix B). We scale their absolute CO column density by a factor of 0.1 to match our N_CO measurements. The differences are likely because these models do not consider CO depletion. Figure 21 shows the comparison of N_CO measured in this work with that of viscous models. Between 150 and 400 au, the measured N_CO profiles of the IM Lup, HD 163296, and MWC 480 disks match well to the R_init = 100 au model, and the GM Aur disk to the R_init = 50 au model. The AS 209 disk is similar to the R_init = 30 au model but likely requires an even smaller R_init for better match. Beyond 400 au, the measured N_CO 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.

**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 α = 10⁻³ around a 1 M_⊙ star after 1 Myr evolution. The initial disk mass in the models is 0.06 M_⊙, and the initial radii, R_init, are 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_init = 50 and 100 au in the region between 150 and 400 au. The CO profile of the AS 209 disk is more consistent with the R_init = 30 au model.
Download figure:
Standard image High-resolution image

In short, the N_CO distributions in 150–400 au of the five MAPS disks match well with predictions of viscously evolving disk models. The initial evidence shows the importance of spatially resolved N_CO profiles to constrain current disk evolution theories. A systematic study of N_CO 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)millimeter 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 snowlines in our models with continuum substructures.

In our thermochemical models, the midplane CO snowline is defined as the radius where the CO gas and ice abundances become equal at the disk midplane. The midplane 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 midplane have 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 the C¹⁷O (1–0) line in the HD 163296 and MWC 480 disks, in which the C¹⁷O line becomes optically thick in the inner disks of at ∼60 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 pileup of dust grains due to the sintering effect at snowlines (Banzatti et al. 2015; Okuzumi et al. 2016; 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¹⁸O and ¹³CO J = (2–1), (1–0), and C¹⁷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:

1.
We retrieve CO gas density distributions, using three independent assumptions for the underlying gas temperature: (1) a thermochemical modeling framework based on the CO data, the broadband SED, and the millimeter continuum emission; (2) an empirical distribution based on optically thick CO emission lines; and (3) a direct estimate based on fits to the C¹⁷O hyperfine lines. Results from all three methods generally show excellent agreement.
2.
By comparing model temperature structures with empirical measurements from optically thick lines, we show that current thermochemical models can provide reliable temperature structures within the 15–40 K region.
3.
We show that fitting C¹⁷O hyperfine structure components provides reliable constraints on the CO gas column density distribution in disks. The results are comparable to those derived from thermochemical models. The hyperfine analysis is independent of the disk model, or any previous knowledge about the temperature structure.
4.
Assuming a gas-to-dust mass ratio of 100, we find that all five disks have a CO-to-H₂ abundance of 10–100 times lower than the ISM ratio of 10⁻⁴. The radial distribution of CO gas column density shows significant variations in terms of the absolute value and the radial profile among the five disks.
5.
The MWC 480 and HD 163296 disks show a steep increase of CO gas column density inside their CO snowlines compared to predictions of thermochemical models, suggesting that 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.
6.
Four of the five disks have remarkably similar surface density profiles between 150 and 400 au, which can be well characterized as Σ_CO ∝ R^−2.4. We show that the N_CO profiles of all five disks between 150 and 400 au are consistent with predictions of viscously evolving disks with an initial disk size between 30 and 100 au.
7.
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% 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 millimeter 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 observations.

In summary, the MAPS observations show that high spatial resolution observations of CO isotopologue lines can provide direct constraints on fundamental disk properties, such as gas thermal structure and gas mass distributions. 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 understanding of typical environments and processes of planet formation.

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

K.Z. acknowledges the support of the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin–Madison with funding from the Wisconsin Alumni Research Foundation. K.Z., K.R.S., J.H., J.B., J.B.B., and I.C. acknowledge the support of NASA through Hubble Fellowship grants HST-HF2-51401.001, HST-HF2-51419.001, HST-HF2-51460.001-A, HST-HF2-51427.001-A, HST-HF2-51429.001-A, and HST-HF2-51405.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. C.J.L. acknowledges funding from the National Science Foundation Graduate Research Fellowship under grant DGE1745303. A.D.B, E.A.B., and A.F. acknowledge support from NSF AAG grant No. 1907653. V.V.G. acknowledges support from FONDECYT Iniciación 11180904 and ANID project Basal AFB-170002. Y.L. acknowledges the financial support by the Natural Science Foundation of China (grant No. 11973090). K.I.Ö acknowledges support from the Simons Foundation (SCOL No. 321183) and an NSF AAG grant (No. 1907653). S.M.A. and J.H. acknowledge funding support from the National Aeronautics and Space Administration under grant No. 17-XRP17 2-0012 issued through the Exoplanets Research Program. J.D.I. acknowledges support from the Science and Technology Facilities Council of the United Kingdom (STFC) under ST/T000287/1. C.W. acknowledges financial support from the University of Leeds and from the Science and Technology Facilities Council (grant Nos. ST/R000549/1, ST/T000287/1, and MR/T040726/1). R.T. and F.L. acknowledges support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow. M.L.R.H. acknowledges support from the Michigan Society of Fellows. 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 Research (S) 18H05222, and Grant-in-Aid for Transformative Research Areas (A) 20H05844 and 20H05847. J.K.C. acknowledges support from the National Science Foundation Graduate Research Fellowship under grant No. DGE 1256260 and the National Aeronautics and Space Administration FINESST grant, under grant No. 80NSSC19K1534. 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. acknowledges support from a CNES fellowship grant. H.N. acknowledges support by NAOJ ALMA Scientific Research grant code 2018-10B and Grant-in-Aid for Scientific Research No. 18H05441. Y.Y. is supported by IGPEES, WINGS Program, the University of Tokyo. L.M.P. acknowledges support from ANID project Basal AFB-170002 and from ANID FONDECYT Iniciación project No. 11181068. T.T. is supported by JSPS KAKENHI grant Nos. JP17K14244 and JP20K04017.

Facility: 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 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 22). 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 22, 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.

**Figure 22.** 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 two-dimensional CO abundance after scaling by the profiles. Rows 3 and 4 show the resulting 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.
Download figure:
Standard image High-resolution image

**Figure 23.** 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¹⁷O (1–0) hyperfine structure transitions after correction for the velocities of the sources.
Download figure:
Standard image High-resolution image

**Figure 24.** 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 in the upper right corners.
Download figure:
Standard image High-resolution image

**Figure 25.** Stellar spectra (top row) and dust absorption and scattering coefficients (bottom row) used in thermochemical models.
Download figure:
Standard image High-resolution image

**Figure 26.** Comparison of radial profiles of the CO isotopologue lines and best-fit results from thermochemical models (red solid line; see model details in Section 3.1.8) and empirical temperature models (blue dashed line; see model details in Section 3.2).
Download figure:
Standard image High-resolution image

**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_ex), and optical depth (τ).
Download figure:
Standard image High-resolution image

**Figure 28.** Gas and CO column density profiles of viscously evolving disk models from Trapman et al. (2020). These models are for α = 10⁻³ around a 1 M_⊙ star after 1 Myr of evolution. The initial disk mass in the models is 0.06 M_⊙, and the initial radii, R_init, are 30, 50, and 100 au, respectively. These models show that between 150 and 400 au the CO column density profile is similar to the gas surface density profile.
Download figure:
Standard image High-resolution image

Appendix B: Additional Figures

Here we show additional figures mentioned in the main text. Figure 23 shows the matched filter responses for C¹⁷O (1–0) line in the GM Aur, IM Lup, and AS 209 disks. Figure 24 shows deprojected CO radial profiles in a logarithmic scale. Figure 25 shows the stellar spectra and dust-opacity profiles used for thermochemical models. Figure 26 compares the radial intensity profiles from models and observations. Figure 27 shows the model results from hyperfine-line fitting of the C¹⁷O (1–0) line. Figure 28 shows a comparison between CO column density profiles and gas-mass distributions in the viscous-evolving disk models of Trapman et al. (2020).

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

Article metrics

Permissions

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. Observations