Systematic Variations of CO Gas Abundance with Radius in Gas-rich Protoplanetary Disks

We employ a global surface density from the self-similar solution of a viscously evolving disk (Lynden-Bell & Pringle 1974), see Equation (1). This profile has been widely used in modeling of protoplanetary disks (e.g., Andrews et al. 2011; Zhang et al. 2014; Cleeves et al. 2016).

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

where Σ_c is the surface density at the characteristic radius R_c, and γ is the gas surface density exponent.

The disk model is composed of three mass components—gas, a small dust population, and a large dust grain population. The separation of different dust populations is used to mimic the effects of dust evolution in protoplanetary disks, including vertical settling and radial drifting (Birnstiel et al. 2012; Krijt et al. 2016). The mass distribution of all three components are assumed to follow the global surface density profile as Equation (1), but each component can have a different set of parameter values. The gas and the small grain population are assumed to be spatially coupled, while the large grain population has a different spatial distribution.

For each mass component, the vertical distribution is characterized by a scale height h(R) and the density ρ is computed through Equation (2)

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

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

where f_i is the mass fraction of each mass component, h₀ is the scale height at 100 au, ψ 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. This is modeled by the settling parameter χ_i, where χ = 1 of the gas and the small grain population, and χ < 1 for the large grain population.

Both dust populations are assumed to follow an MRN size distribution $n(a)\propto {a}^{-3.5}$ with a minimum grain size a_min = 0.005 μm, and a_max = 1 μm for the small grain population and a_max = 1 mm for the large grain population. We use Mie theory to compute the wavelength dependence of dust opacity.

3.2.1. Disk Parameters

All four disks are well-studied disks and their disk properties have been under intensive study in the literature. For our purpose of studying CO abundance structure, the most important parameters are the gas surface density profile and the disk flaring level. The former sets the baseline of CO gas distribution, and the latter largely determines the temperature structure in a disk. When gas/small dust distribution constraints are available from scattered light images, we adopt these values from previous studies. If not, we adopt γ_g = 1 for gas and small dust grains. The gas disk sizes are adopted from scattered light images or previous studies of multiple CO isotopologue lines. We vary the γ index of gas surface density profiles to test the uncertainties in the derived CO abundance structures in Section 4.

For a given surface density profile, the flaring level of the disk can be constrained by fitting its broadband spectral energy distribution (SED; Dullemond & Dominik 2004). The mass ratio of the small-to-large grains and the settling of large grains can affect the thermal structure and the UV penetration in a disk. Schwarz et al. (2018, 2019) had carried out systematic tests on the effect of small-to-large dust mass ratio on the CO processing chemistry. They showed that the CO depletion happened faster in models with higher small-to-large dust mass ratio. Here we adopt an intermediate level of 80% of the mass in large grains for modeling. The only exception is the IM Lup disk for which we adopt 99% of the mass in large grains as this value was used on the comprehensive model of the IM Lup disk of Cleeves et al. (2016). Here we use a settling parameter of χ = 0.2 for the large grain population. The detailed disk structures inside 20 au is not critical for our analysis as we focus on CO gas distribution on a scale of tens of astronomical units. We also ignore small-scale substructures in these disks (Andrews et al. 2016, 2018; Zhang et al. 2016), as this work focuses on the large-scale radial variations. Please see Section 5.2 for a more detailed discussion on possible effects of substructures on CO abundance distribution.

Two of these four disks (TW Hya and DM Tau) have HD (1-0) line observations and previous models suggest that their gas-to-dust mass ratios are close to 100 (Bergin et al. 2013; McClure et al. 2016). For simplicity, we use a gas-to-dust (small+large) mass ratio of 100 for all four disks. The uncertainty in the total gas mass is not crucial for the study here, as we focus on the radial variation patterns.

Table 3.  Source Information

Source	M⋆	R⋆	Teff	D	Md	Mgas	Vlsr	incl	PA	LX	References
	(M⊙)	(R⊙)	(K)	(pc)	(M⊙)	(M⊙)	(km s−1)	(deg)	(deg)	(erg s−1 cm−2)
DM Tau	0.53	1.25	3705	140	5.00E-04	5.00E-02	5.95	35	335	3.00E+29	1, 2, 3, 4, 5
TW Hya	0.8	1.04	4110	60	5.00E-04	5.00E-02	2.84	7	335	1.60E+30	2, 6, 7, 8
HD 163296	2.47	1.69	9250	101	1.40E-03	1.50E-01	5.74	49	312	4.00E+29	2, 9, 10, 11, 12, 16
IM Lup	1	2.5	3900	160	1.70E-03	1.70E-01	4.5	48	144	4.30E+30	2, 13, 14, 15

Note. (1) Kenyon & Hartmann (1995); (2) Gaia Collaboration et al. (2018); (3) Andrews et al. (2011); (4) Bergin et al. (2016); (5) Henning et al. (2010); (6) Andrews et al. (2012); (7) Huang et al. (2018a); (8) Brickhouse et al. (2010); (9) Natta et al. (2001); (10) Isella et al. (2016); (11) Rosenfeld et al. (2013a); (12) Swartz et al. (2005); (13) Pinte et al. (2008); (14) Cleeves et al. (2016); (15) Günther et al. (2010); (16) de Gregorio-Monsalvo et al. (2013).

Download table as:  ASCIITypeset image

Table 4 lists all the parameters adopted for our disk models. Here we briefly describe the adoption of parameters for each disk:

• 1.  
  DM Tau: we employ a γ_g = 1 for gas and small grains, as no scattered light image studies are available. The inner edge of gas is set to be 4 au (Calvet et al. 2005). The characteristic radius of gas, ${R}_{c}^{{\rm{g}}}$, is set to 270 au based on the extent of C¹⁸O (3–2) emission in channels close to the stellar velocity. For the large grain population, we use γ_mm = 1, ${R}_{c}^{\mathrm{mm}}=135$ au, and an inner radius of 20 au based on the model of Andrews et al. (2011), which matches the 880 μm continuum observations of DM Tau.
• 2.  
  TW Hya: we adopt the general gas surface density profile of TW Hya derived from scattered light images from van Boekel et al. (2017). This distribution has a global γ_g = 0.7 and a exponential taper beyond 104 au. The distribution of millimeter-grain is truncated at 70 au, with a ${R}_{c}^{\mathrm{mm}}$ = 60 au (Andrews et al. 2012).
• 3.  
  HD 163296: we employ the basic structure of Isella et al. (2016), which fitted the 1.3 mm continuum emission and the general sizes of CO/¹³CO/C¹⁸O (2–1) isotopologue emissions. The model uses ${R}_{c}^{{\rm{g}}}=165$ au and γ_g = 0.8 for gas and small grains, and ${R}_{c}^{\mathrm{mm}}=90$ au, γ_mm = 0.1 for the large grain population.
• 4.  
  IM Lup: we adopt ${\gamma }_{{\rm{g}}}=1$ for gas, and other surface density profile parameters from Cleeves et al. (2016), which carried out a comprehensive model of 875 μm continuum and multiple CO (2–1) isotopologue emissions.

Table 4.  Disk Model Parameters

	DM Tau	TW Hya	HD 163296	IM Lup	Units	Definition
γg	1	0.75	0.8	1		surface density index of gas
γmm	1	1	0.1	0.3		surface density index of large grain
${R}_{c}^{g}$	270	104	165	100	(au)	characteristic radius of gas
${R}_{c}^{\mathrm{mm}}$	135	60	90	100	(au)	characteristic radius of large grain
${M}_{\mathrm{mm}}$	4.00E-04	4.00E-04	1.20E-03	1.70E-03	(M⊙)	mass of large grain
Msmall	1.00E-04	1.00E-04	3.00E-04	1.70E-05	(M⊙)	mass of small grain
χmm	0.2	0.2	0.2	0.2		scale fraction of large grain
H100 (au)	6	6	8.5	8	(au)	gas scale height at 100 au
ψ	1.25	1.25	1.08	1.2		flaring parameter
${R}_{{in}}^{\mathrm{mm}}$	19	1	0.5	0.2	(au)	inner radius of large grain
R${}_{\mathrm{in}}^{{\rm{g}}}$	3	4	0.5	0.2	(au)	inner radius of gas/small grain
R${}_{\mathrm{out}}^{\mathrm{mm}}$	300	70	600	313	(au)	outer radius of large grain
R${}_{\mathrm{out}}^{{\rm{g}}}$	500	200	600	1200	(au)	outer radius of gas/small grain

Download table as:  ASCIITypeset image

3.2.2. Stellar Parameters

The general modeling processes are outlined in Figure 2.

Zoom In Zoom Out Reset image size

We first use the radiative transfer code RADMC3D to calculate the dust thermal structure in a disk (Dullemond et al. 2012) and then generate model SEDs to compare with observations. To constrain the vertical scale height of dust, we first run a grid of h₀ and ψ to search for best-fitting SED models of each disk. The h₀ range is from 0.02 to 0.2 with a step size of 0.02, and ϕ from 1 to 1.3 with a step size of 0.05. The best-fitting SED models are shown in Figure 3.

Zoom In Zoom Out Reset image size

Then we adopt the best-fitting values of density structure into a full 2D time-dependent thermo-chemical code RAC2D. This code computes the disk chemistry and the heating-cooling balance for both the gas and the dust self-consistently. We briefly describe the code here and interested readers can find a full description in Du & Bergin (2014). Given a static density structure of gas and dust, the code first solves the dust temperature and radiation field (from X-ray to centimeter wavelengths) in the disk using a Monte Carlo approach. The cosmic-ray ionization inside the disk is simulated with an attenuation length of 96 g cm⁻² (Umebayashi & Nakano 1981; Bergin et al. 2007). For the given radiation field, the code then simultaneously solves the chemical evolution and gas thermal structure in the disk. The chemical network has 467 species and 4801 reactions, including the full gas-phase network from the UMIST 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 chemical models are initialized with a composition listed in Table 1 of Du & Bergin (2014). In particular, all carbon is initially in CO gas with an abundance of 1.4 × 10⁻⁴ relative to H atoms and a binding energy of 855 K for CO ice (Öberg et al. 2005). We let the chemistry run for 1 Myr for all models. No isotopologue fractionation is considered in the chemical network, as the fractionation is expected to be insignificant in massive disks like this sample (Miotello et al. 2014). We rescale the output CO abundance by the local ISM ratios of CO/¹³CO = 69, and CO/C¹⁸O = 570 to generate CO isotopologue abundances (Wilson 1999).

Using the gas temperature and CO abundance structures, we compute line images using the ray-tracing module of RAC2D, assuming an LTE condition. We initially simulate the images at a 2 au spatial resolution and 10× spectral resolution as the observations, and then convolve model images to the same spatial and spectral resolution as observations.

We test model uncertainties of two important parameters: the cosmic-ray ionization rate (ζ_CR) and the γ index for surface density profiles.

Previous chemical studies suggested that an ISM level of cosmic-ray (CR) ionization rate (ζ_CR ≥ 1.36 × 10⁻¹⁷ s⁻¹) is critical for sufficient CO processing to other molecules (Bosman et al. 2018; Schwarz et al. 2018). This level of cosmic-ray rate is typical in ISM but it is unclear whether this is also present in disks, as disk wind or magnetic fields in actively accreting T Tauri disk systems may attenuate the cosmic-ray rate in disks (Cleeves et al. 2013, 2015). To test the impact of CR rate on our chemical models, we run two sets of models: a set of standard models with a ζ_CR = 1.36 × 10⁻¹⁸ s⁻¹ and the other set with a high CR rate of 1.36 × 10⁻¹⁷ s⁻¹.

We also test how uncertainties on the slope of surface density profiles would impact the constraints on CO depletion profile. Therefore, for each disk, we run two additional models with γ varies by ±0.2. The characteristic radii and total gas masses are kept the same.

In Figure 4, we show the CO gas abundance and gas temperature structures in the standard models (ζ_CR = 1.36 × 10⁻¹⁸ s⁻¹) of these four disks. All four disks show a molecular layer in which the CO gas abundance is close to 10⁻⁴. This layered structure has been seen in channel maps of CO line observations of HD 163296 and IM Lup (Rosenfeld et al. 2013b; Pinte et al. 2018a). The last column in Figure 4 shows the column density distributions CO gas and ice in the four disk models. The CO (gas+ice) column density distribution is close to 10⁻⁴ at all radii, indicating very little CO has been converted to other carbon carriers after 1 Myr in these low CR rate models.

Zoom In Zoom Out Reset image size

When no significant conversion of CO into other carriers occurs, the CO snowline location is predominantly determined by the thermal structure of the disk. DM Tau is the coldest disk in our sample and its CO midplane snowline is around 10–15 au in the standard model. In the TW Hya model, the CO midplane snowline is around 20 au, consistent with 21 ± 1.3 au derived from ¹³C¹⁸O (3–2) line observations (Zhang et al. 2017) and the N₂H⁺ line image (Qi et al. 2013). For the HD 163296 model, its CO midplane snowline is between 50 and 70 au, consistent with the previous estimation from C¹⁸O and N₂H⁺ line observations, after correcting to the Gaia2 distance (Qi et al. 2015). For the IM Lup model, its CO ice starts to appear around 30 au, but its CO ice column density increases slowly with radius compared with other disks, because the IM Lup disk is highly flaring and has a thicker CO gas layer than others.

In short, the four standard models (ζ_CR = 1.36 × 10⁻¹⁸ s⁻¹) all show a layered CO gas structure and the midplane CO snowline locations are consistent with existing observational constraints. In these models, less than 10% of total CO (gas+ice) have been converted to other carbon carriers after 1 Myr.

Figure 5 shows the comparison of CO gas column density distribution in high and standard models (ζ_CR = 1.36 × 10⁻¹⁷ s⁻¹ versus 1.36 × 10⁻¹⁸ s⁻¹). Models of DM Tau and TW Hya show substantial CO gas depletions between 10 and 100 au, by a factor of 5–20 after 1 Myr. The HD 163296 model shows an intermediate level of depletion by up to a factor of 10 between 30 and 200 au. The IM Lup model shows the smallest variation between high and low CR rate models with a depletion factor of less than 4 throughout the whole disk.

Zoom In Zoom Out Reset image size

In general, the chemical processing rates of CO gas depend on the thermal and ionization structures of disks, which in turn depends on many physical properties, including the stellar luminosity, UV luminosity, disk flaring structure, and dust size distribution. For a given cosmic-ray rate, the available surface areas on grain surface have a significant impact on the CO processing rate. Schwarz et al. (2018) and Bosman et al. (2018) have carried out systematic studies on how disk properties would impact the CO depletion level, including small/large grain mass ratio and temperature. In our models, there are some very noticeable effects. The effect of grain surface area is most obvious in the IM Lup disk models that have only 1% of the solid mass in small grains. Limited by the available grain surface areas, the chemical processing of CO in the IM Lup models is order(s) of magnitude smaller than that of other disks. The effect of thermal structure can be seen in HD 163296 models. As a Herbig disk, its gas temperature is much warmer than T Tauri disks and thus has less CO freeze-out regions for CO to be processed in the grain surface path.

Despite significant decreases in the CO gas column density, Figure 5 shows that the surface brightness profiles of C¹⁸O line emission are not very sensitive to the CO depletion caused by high CR ionization. This is because most of the depletion occurs in the deep region of these disks and the CO emission above the depletion region is already optically thick. These results indicate that the CO depletion by chemical processing alone, even by high CR rate, is unlikely to account for the low surface brightness observed in these four disks if their gas-to-dust ratio is still close to 100.

Figure 6 shows how the profiles of total gas surface column density (Σ_g), CO gas column density (Σ_CO), and the surface brightness of CO isotopologue line emission (I_CO) vary in models deviating from the γ of standard models by ±0.2. These results show that Σ_CO represents well the gas column density change within the midplane CO snowline, but it becomes less sensitive to the gas surface density changes beyond the CO snowline as the CO gas layer only accounts for a small fraction of the total gas column. The I_CO profiles only change less than 30% for uncertainty of 0.2 in the γ index.

Zoom In Zoom Out Reset image size

Here we compare the surface brightness distributions of CO lines (I_CO) between observations and our standard models to estimate the levels of variations in the radial CO gas abundance. For each disk, we generate a set of synthetic line images by rescaling the model CO gas abundance throughout the disk with a constant factor. The grid spans from a depletion factor of 0.1–400, increasing by a factor of 1.2 at each step. Figure 7 shows the line surface brightness profiles of models with exemplary depletion scales, comparing with the observed profiles in the four disks. Using the grid of models, we record which depletion factor matches the best with models at each radius.

Zoom In Zoom Out Reset image size

Figure 8 shows the required CO depletion distributions to match observations. On the disk-averaged scale, we find that significant depletion of CO gas abundance is needed to explain observations. This is consistent with the previous results of disk-averaged CO abundance studies (Favre et al. 2013; Flaherty et al. 2015; Cleeves et al. 2016; McClure et al. 2016). For the radial variation in individual disks, we find that one order of magnitude radial variation is usually needed to account for observations. Comparing to the changes of I_CO due to γ uncertainty, it is unlikely these large deviations can be explained by surface density profiles alone.

Zoom In Zoom Out Reset image size

For individual disks, the DM Tau disk shows a nearly constant level of CO depletion inside 100 au (by a factor of 12) and its CO abundance gradually increases with radius to an interstellar ratio at 280 au. TW Hya shows a more complicated pattern, with a sharp decrease of CO abundance between 10 and 30 au, and then increases with the radius until 60 au, and outside 100 au the CO abundance is extremely depleted (> a factor of 200). IM Lup shows a monotonic trend of CO abundance depletion, with an outer disk that has less depletion. HD 163296 shows a modest CO abundance enhancement inside 70 au (by a factor of 2–5), its 100–200 au shows the largest depletion and then a shallow increase of CO abundance with radius.

Except for the oldest disks TW Hya, the other three disks show a similar pattern in their CO abundance structures in regions beyond their CO snowline—the CO gas abundance is low in the intermediate disk region outside the midplane CO snowlines, but gradually increases with radius to a higher value in the outermost region. Interestingly, this pattern is consistent with the general trends of carbon elemental abundance suggested by other carbon carriers. Bergner et al. (2019) analyzed C₂H radial brightness distributions in the same three disks and their models required a high depletion of elemental carbon abundance (by a factor of 10–100) between 100–200 au and only a modest depletion beyond 200 au (by a factor of 1–10). Therefore, the CO depletion profiles likely represent a general depletion of carbon elemental abundance in the atmospheres of these disks.

Summary of model results: (1) In detailed thermo-chemical models for the four disks, a high CR ionization rate can significantly reduce the CO gas column between 50-100 au up to a factor of 20 after 1 Myr. However, the depletions by chemical processing mostly happen at the deep interior of these model disks, and the resulting CO abundance structures cannot reproduce the weak CO isotopologue lines observed. (2) Comparing these thermo-chemical models to spatially resolved observations, we find that the CO gas abundance likely varies significantly with radius. An interesting pattern in the radial profiles is that the CO gas abundance is low in the intermediate disk region outside the midplane CO snowlines, but gradually increases with radius to a higher value in the outermost region.

Current surveys of nearby star formation regions suggest that gas mass derived from CO observations is 1–2 orders of magnitude too low compared to the canonic CO ratio in ISM (Ansdell et al. 2016; Long et al. 2017; Miotello et al. 2017). In a broad sense, two types of depletion mechanisms have been proposed: (1) the first one is chemical processes that convert CO into other less volatile molecules (e.g., CO₂, CH₃OH, CH₄; Aikawa et al. 1999; Bergin et al. 2014; Reboussin et al. 2015; Eistrup et al. 2016; Yu et al. 2016; Bosman et al. 2018; Schwarz et al. 2018); (2) the second type is physical processes that CO gas abundance is altered as a result of the CO freeze-out onto dust grains and the subsequent dust growth, settling, and radial drifting in protoplanetary disks (Kama et al. 2016; Xu et al. 2017; Krijt et al. 2018). In other words, CO becomes trapped as ice in large grains that reside in the midplane or perhaps even in planetesimals. Previous studies of CO depletion mechanisms were mostly theoretical works. Here we discuss our results of radial CO variations in the context of these mechanisms proposed.

5.1.1. Chemical Depletion Mechanisms

5.1.2. Physical Depletion Mechanisms

5.1.3. Ways to Increase CO Gas Depletion in Disks

Recent ALMA long-baseline observations have revealed that many disks show 1–20 au scale substructures in their dust continuum emission (Andrews et al. 2018). Prominent substructures in the gas and dust mass distributions may alter the local CO gas emission. For example, Facchini et al. (2018) found that the gas inside gaps opened by 1–5 M_J planet becomes colder than the surrounding environment in their thermal-chemical models. The depletion of gas and temperature variation across gaps can affect a radial profile of the CO abundance derived from a smooth surface density structure such as ones used here.

Although the origins of these dust substructures are still under debate, proposed mechanisms generally predict that substructures in CO line emission would be much narrower and shallower than those of the dust emission (Facchini et al. 2018). Three of the four disks (TW Hya, DM Tau, and IM Lup) only show modest substructures in continuum emission beyond 20 au, with up to a factor of 3 depletion inside gaps (Huang et al. 2018b; Kudo et al. 2018). HD 163296 has the most significant substructures among the four disks, showing a factor of 33 depletion inside one of its gaps in the 1.3 mm continuum (Huang et al. 2018b). But even for the HD 163296 disk, its ¹³CO and C¹⁸O (2–1) line observations at 02 resolution only show subtle breaks in the slopes rather than clear gaps (Isella et al. 2016). In short, based on the depths and widths of known dust substructures in the four disks and current theoretical predictions, it is unlikely that the order of magnitude variations in CO abundance profiles are due to substructures alone.

On the other hand, it is important to study how the dust substructures may affect the local chemical substructures and as these local effects may significantly impact the compositions of forming planets. A Large ALMA Cycle 6 program is currently ongoing to study the detailed interplay between dust substructures and local chemical variations.

In our small sample, the radial profiles of CO depletion are diverse. Here we discuss how disk properties might play a role in regulating the radial CO depletion profile.

HD 163296 is the only Herbig disk in our sample and it shows a very similar depletion pattern to the dust evolution models of Krijt et al. (2018). Assuming a canonical gas-to-dust mass ratio of 100, it has a modest CO gas depletion in this sample (by a factor of 5 on disk average). Current chemical models suggest that the chemical processing of CO gas is most efficiently happening at the cold interior of a disk. Herbig disks have a larger fraction of the disk mass above the CO condensation temperature than that of T Tauri disks, and therefore they are expected to have less CO gas processing through the chemical pathways (Bosman et al. 2018). Besides HD 163296, another example of Herbig disk with modest carbon depletion is HD 100546. Kama et al. (2016) analyzed major carbon carriers (CO, C i, C ii) for the HD 100546 disk, and concluded that the elemental carbon abundance in its warm molecular layer is close to the interstellar ratio or only modestly depleted, with C/H = (0.1–1.5) × 10⁻⁴. Considering the old ages of the two Herbig disks (HD 163296 ≃12 Myr, HD 100546 between 4.1 and 16 Myr, Andrews et al. 2018; van der Marel et al. 2019), the chemical and physical processes should already have sufficient time to occur. The fact that the CO depletion is modest in these two warm disks suggests that temperature may play an important role in CO depletion processes. A larger sample of Herbig disks can easily confirm or rule out the effect of warm temperature in the CO depletion observed in protoplanetary disks.

Among the three T Tauri disks studied here, IM Lup is the most massive disk. Although it appears to be relatively young (∼1 Myr), its disk-averaged CO depletion is high, by a factor of 20 (Cleeves et al. 2016). Our model predicts that even with high cosmic-ray ionization rate, the expected CO depletion is a factor of 2–4. On the other hand, the high CO depletion in IM Lup is consistent with the expectation that dust evolution is faster in a massive disk.

DM Tau is the coldest disk in our sample but its disk-averaged CO depletion is just modest. This appears to be inconsistent with chemical processing of CO happening most efficiently at the cold region of a disk. On the other hand, the DM Tau disk is among the largest disks in both gas lines and (sub)millimeter continuum observations, and its column density is generally lower than those of the other three disks (see Figure 6). The low-density environment can lead to slower dust growth and less CO depletion.

The TW Hya disk shows the largest CO gas depletion on disk average. In particular, its outer disk region (>100 au) appears to be extremely depleted of CO (by a factor of 100–300). This extreme CO depletion in the outermost region is puzzling, as it can be explained by neither chemical processing nor dust evolution. In particular, molecular ion observations toward the TW Hya disks indicate its ionization level due to CR is quite low (ζ_CR < 10⁻¹⁹ s⁻¹) (Cleeves et al. 2015). On the dust evolution side, the highest depletion occurs in the outermost region, which is inconsistent with the predictions of dust evolution simulations.

Given these arguments, we suggest that the weak CO emission beyond 100 au of the TW Hya disk is more likely to be a result of gas dissipation beyond 100 au. TW Hya is one of the oldest T Tauri disks (5–10 Myr old) and thus is likely more evolved than other systems. Although its HD (1-0) line flux suggests its total gas disk mass is still high, the HD line mostly traces the relatively warm region inside 50 au (Bergin et al. 2013; Trapman et al. 2017). A lower limit of the gas mass beyond 100 au in the TW Hya disk has been constrained by ratios of resolved CS lines (Teague et al. 2018b). Comparing the difference between the minimum mass and the one used in our model, the CO abundance can be as high as only modestly depleted by a factor of 2–3. In summary, a gas-poor disk beyond 100 au is consistent with current observations of the TW Hya disk.