Gas Disk Sizes from CO Line Observations: A Test of Angular Momentum Evolution

The size of a disk encodes important information about its evolution. Combining new Submillimeter Array (SMA) observations with archival Atacama Large Millimeter Array (ALMA) data, we analyze mm continuum and CO emission line sizes for a sample of 44 protoplanetary disks around stars with masses of 0.15--2\,$M_{\odot}$ in several nearby star-forming regions. Sizes measured from $^{12}$CO line emission span from 50 to 1000\,au. This range could be explained by viscous evolution models with different $\alpha$ values (mostly of $10^{-4}-10^{-3}$) and/or a spread of initial conditions. The CO sizes for most disks are also consistent with MHD wind models that directly remove disk angular momentum, but very large initial disk sizes would be required to account for the very extended CO disks in the sample. As no CO size evolution is observed across stellar ages of 0.5--20\,Myr in this sample, determining the dominant mechanism of disk evolution will require a more complete sample for both younger and more evolved systems. We find that the CO emission is universally more extended than the continuum emission by an average factor of $2.9\pm1.2$. The ratio of the CO to continuum sizes does not show any trend with stellar mass, mm continuum luminosity, or the properties of substructures. The GO Tau disk has the most extended CO emission in this sample, with an extreme CO to continuum size ratio of 7.6. Seven additional disks in the sample show high size ratios ($\gtrsim4$) that we interpret as clear signs of substantial radial drift.


INTRODUCTION
Circumstellar disks are the cradles of young planets. Their physical structures and chemical conditions largely determine the properties of the resulting planetary systems. Detailed characterization of the gas and dust distributions in disks is therefore fundamental for developing theoretical models of planet formation (e.g., Benz et al. 2014, Morbidelli & Raymond 2016. These disk structures are closely linked to how disks evolve. However, the physical mechanisms that drive one key component of disk evolution -the angular momentum transport -are still debated (e.g., Turner et al. 2014).
In the classical viscous evolution model, mass accretion onto the central star requires the outward transport of angular momentum (Lynden-Bell & Pringle 1974;Hartmann et al. 1998). Consequently, the gas density distribution becomes more radially extended over time.
In the alternative magneto-hydrodynamical (MHD) disk wind model, angular momentum is transported in the wind, moving vertically away from the disk (Blandford & Payne 1982;Ferreira 1997;Gressel et al. 2015;Bai et al. 2016). That behavior can result in a more radially compact gas density distribution (e.g., Lesur 2021;Yang & Bai 2021). Comparing disk sizes in a large sam-ple that spans a range of evolutionary stages could help differentiate between the two scenarios.
Measurements of disk sizes are available from observations with (sub-)millimeter interferometers, which can spatially resolve disks in nearby ( 200 pc) star-forming regions. The (sub-)mm continuum emission from dust has been resolved in ∼200 disks (e.g., Tripathi et al. 2017;Andrews et al. 2018a;Hendler et al. 2020). However, these continuum observations primarily trace the mm-sized particles that have aerodynamically decoupled from the gas and migrated inwards toward any pressure maxima (Weidenschilling 1977). The extent of the continuum emission is therefore a diagnostic of the evolution of disk solids (e.g., Testi et al. 2014;Birnstiel & Andrews 2014). Observations of the molecular gas reservoir in disks are the key probes of the disk evolution tied to angular momentum transport. The most common tracer of the gas disk is the bright CO line emission. However, as the gas emission usually emerges in a narrow velocity range, detecting these lines is more challenging, and systematic CO disk surveys are still limited (e.g., Ansdell et al. 2016;Long et al. 2017;Boyden & Eisner 2020).
The size of the bright 12 CO line emission is found to be larger than that of the mm continuum emission in disks (e.g., Panić et al. 2009;Andrews et al. 2012;Ansdell et al. 2018). This size discrepancy has been widely interpreted as the consequence of dust evolution, associated with the growth and inward migration of solids (Birnstiel & Andrews 2014;Trapman et al. 2019). But some of the discrepancy can be attributed to the different optical depths of the size tracers (Hughes et al. 2008;Facchini et al. 2017). How well the CO-to-continuum size ratio probes the relative spatial distribution of the gas and solid reservoir masses, how it varies with stellar and disk properties, and how the inferred behaviors would relate to fundamental physical processes, are still open questions.
In this work, we combine new SMA CO line observations with archival ALMA data to build a collective view of CO and mm continuum disk sizes. With a sample of 44 targets across a wide range of stellar and disk properties, we explore how the CO and mm continuum sizes can help improve our understanding of disk evolution processes. We introduce the sample and the data in Section 2, and the method for size measurement in Section 3. The sizes and their relationships with stellar/disk properties are summarized in Section 4. In Section 5, we discuss the role of the CO emission sizes in differentiating between turbulent viscosity and MHD disk wind models and the CO-to-continuum size ratio in probing the dust disk evolution. Finally, we summarize the main findings in Section 6.

Three Disks with New SMA Observations
We observed DL Tau, GO Tau, and UZ Tau 1 with eight 6 m antennas of the SMA (Ho et al. 2004) in both extended and subcompact configurations. These targets were observed in two shared tracks in the extended configuration on 2019 November 10 and 11, and one shared track in the subcompact configuration on 2019 November 19. The dual-sideband receivers were both tuned to the same local oscillator (LO) frequency of 225.5 GHz to maximize the sensitivity in the 12 CO J = 2 − 1 line. The SWARM correlator was configured with 16384 channels in each of the four spectral chunks per sideband, with a spectral resolution of 140 kHz (or 0.18 km s −1 velocity resolution at 230.538 GHz) and covering the spectral ranges of 213.  GHz. The science targets were observed in an alternating sequence with the gain calibrators 3C111 and 0510+180. The bandpass calibrator 3C454.3 and flux calibrators MWC349A and Titan were observed at the beginning of each track.
The raw visibility data were calibrated using the IDLbased MIR software package 2 , following the standard SMA data reduction procedures. After the initial flagging of data with abnormal phase and amplitude variations, the bandpass response was determined with the observation of the bright quasar 3C454.3. The amplitude scale was set based on the frequent monitoring of MWC349A and Titan, with a typical systematic uncertainty of ∼ 10%. The repeated observations of 3C111 and 0510+180 provided the complex gain response of the system, which were then applied to the science targets. The calibrated visibility data 3 for individual spectral chunks were imported into CASA for further imaging. Using CASA version 5.6.0 (McMullin et al. 2007), the continuum baseline was first subtracted with the uvcontsub task in the spectral chunks covering the 12 CO J = 2 − 1 line. The continuum-subtracted visibility data from both receivers and antenna configurations of individual disks were then combined using the concat task. The 12 CO image cubes with channel widths of 0.25 km s −1 were produced using the tclean task with natural weighting to maximize the sensitivity, resulting in a typical beam size of 1. 3 × 0. 9 with an rms noise level of ∼80 mJy beam −1 per channel (∼5 K). Keplerian masks were applied in the CLEAN process, which were generated based on stellar and disk parameters from previous observations for each disk (Long et al. 2018) and designed to encompass the observed emission in individual channels.

Archival Data
Besides the three disks with new SMA observations, we collated a sample of disks for which spatially resolved 12 CO observations are available from the ALMA archive 4 for a feasible measurement of the CO disk size. The selected sample includes most well-studied systems and spans a wide range in stellar and disk properties.
The ALMA Large Project DSHARP (Disk Substructures at High Angular Resolution Project, Andrews et al. 2018b) observed twenty disks in the 1.3 mm continuum emission and 12 CO J = 2 − 1 lines at ∼ 0. 03 resolution. We selected the twelve disks with considerably low cloud contamination in the 12 CO emission. These disks are mostly located in the Lupus, Oph, and Upper Sco star-forming regions. The image products from the project website 5 were used in this study.
The ALMA Lupus survey reported gas disk sizes for 22 disks based on 12 CO J = 2−1 observations at ∼ 0. 25 resolution (Ansdell et al. 2018). As five disks (GW Lup, IM Lup, MY Lup, Sz 129, and HT Lup 6 ) have higher quality data from DSHARP, we included the remaining 17 Lupus disks, and directly adopted the size measurements in Ansdell et al. (2018) for the full sample analysis in this work.
Six disks around low mass stars (0.1-0.2 M ) in the Taurus star-forming region were included based on the recent work of Kurtovic et al. (2021), which targeted the 890 µm continuum emission and 12 CO J = 3 − 2 line at ∼ 0. 1 resolution. DM Tau is historically known to host an extended CO gas disk (Guilloteau & Dutrey 1998;Piétu et al. 2007). The ALMA data with ∼ 0. 3 resolution by Flaherty et al. (2020) were used for DM Tau in our analysis. The 12 CO J = 2 − 1 emission from the FP Tau disk was recently observed in an ALMA Chemistry program with ∼ 0. 2 resolution (Pegues et al. 2021). We also included the measurements for CX Tau, which was reported as the first example with very high gas-to-dust disk size ratios (R CO /R mm ∼ 4, Facchini et al. 2019).
Three additional sources from other star-forming regions with available 12 CO data were also included: TW Hya (Huang et al. 2018a), V4046 Sgr (Rosenfeld et al. 2012), and J11004022-7619280 (Pegues et al. 2021). Data images for the three sources and the above Taurus disks were kindly shared by authors of these references.

Host Star Properties
The full sample includes 12 disks from Taurus, 21 disks from Lupus, 5 disks from Oph, and 6 disks from an assortment of five other regions. Two sources (UZ Tau and V4046 Sgr) are spectroscopic binaries and host circumbinary disks. In addition, UZ Tau and Sz 123 A have stellar companions with moderate separations (3. 6 and 1. 7, respectively, Kraus & Hillenbrand 2009;Ghez et al. 1997). No additional stellar companions were reported for the remaining sample to our best knowledge (e.g., Kraus et al. 2012;Zurlo et al. 2020). Therefore, the effect of tidal truncation is minimal for this sample.
Target distances (d) were estimated based on Gaia DR2 parallax measurements (Gaia Collaboration et al. 2018). The stellar properties for individual disks were adopted from previous studies (accounted for updated Gaia distances) and summarized in Table 1. Though different evolutionary models were used to derive stellar parameters, no significant differences were observed. Specifically, for Taurus targets, effective temperature (T eff ) and luminosity (L * scaled by d 2 ) taken from the optical spectral survey of Herczeg & Hillenbrand (2014) were used to calculate stellar mass and age using Baraffe et al. (2015) and the non-magnetic Feiden (2016) models of pre-main sequence stellar evolution (Long et al. 2019). Stellar masses and ages for most Lupus disks and the DSHARP sample were estimated based on the MIST models using literature values of T eff and L * (Andrews et al. 2018a,b). For the few remaining targets, stellar properties were adopted from the corresponding references summarized in Table 1. The listed stellar masses for UZ Tau and V4046 Sgr are the total mass of the two stellar components, as derived from the gas disk rotation (Rosenfeld et al. 2012;Czekala et al. 2019). We also updated the stellar masses for a few systems when dynamical stellar mass measurements are available (see Table 1).
This sample covers a wide range of young star properties, with stellar masses of 0.15-2 M (spectral types from M5 to A1) and spanning more than two orders of magnitude in stellar luminosities. The stellar ages range from 0.5 to 20 Myr, while the average age is about 2-3 Myr as most of the disks are located in the Taurus Though the stellar mass distribution is more uniform, this sample lacks low mass host stars at younger and more evolved stages. A typical errorbar in log(M * ) and log(t * ) is shown in the lower right corner and corresponds to ±0.1 and ±0.3 dex, respectively. and Lupus regions. As shown in Figure 1, this sample is particularly incomplete for lower mass host stars towards younger and more evolved systems.

SIZE MEASUREMENTS
The most general definition of the size of an object is taken as the distance from its center to the outer edge. For astrophysical objects like disks, which fade around the edges and merge into the interstellar environment, the above definition is not applicable due to observational limitations. For practical purposes, Tripathi et al. (2017) introduced the disk effective radius (R eff ) concept, the radius that encompasses a fraction (x) of its total flux. For a constructed cumulative intensity profile, f ν (r) = 2π r 0 I ν (r )r dr , R eff is taken as f ν (R eff ) = xF ν , where F ν = f ν (∞) is the total flux. This effective radius definition has been widely applied in disk studies using the intensity profiles inferred from modeling the continuum visibilities Long et al. 2019;Hendler et al. 2020). Following a similar concept, some other works use increasing elliptical apertures in the image plane to construct a cumulative flux curve and determine the effective radius (Ansdell et al. 2018;Huang et al. 2018b).
In this paper, we define the size for a given disk tracer as the 90% effective radius based on the derived azimuthally averaged radial profiles from deprojected data images. The choice of x = 0.9 ensures the inclusion of most of the disk area. Meanwhile, this choice is consistent with the measurements of the 17 Lupus disks from the ALMA Lupus Survey (Ansdell et al. 2018), which enables a joint analysis by directly adopting their results. We also tabulate the measurements with x = 0.68, 0.95 in Appendix Table B2.
For the mm continuum size, we first deproject the continuum image using the known disk inclination and position angles from previous studies, and then extract the azimuthally averaged radial intensity profile. The size (R mm ) is then measured from the cumulative intensity profile, which is constructed by integrating the radial intensity profile. The statistical uncertainty of the size is estimated via bootstrapping over 1000 random sample profiles. Each sample profile is drawn by perturbing the intensity at individual radial bins assuming a Gaussian noise distribution, where the noise level at each radial bin is calculated as the 1σ scatter divided by the square root of the number of beams spanning the azimuthal angles over which the intensity is measured 7 . The adopted disk geometry parameters, continuum sizes, as well as some continuum observation details are summarized in Table 2.
The CO emission size is measured in the same way, using the radial emission profile extracted from the CO moment-zero map, obtained by integrating the data cube along the velocity axis. To boost SNR for the faint line emission, especially in the outer disk, we apply a Keplerian mask customized for each disk when creating the moment-zero map. This ensures that only regions associated with disk emission are included. This method has been used in previous work to enable weak line detection and increase SNR for faint emission (e.g., Salinas et al. 2017;Matrà et al. 2017;Yen et al. 2018;Teague et al. 2018;Loomis et al. 2018). The basic emission morphology at each velocity channel can be well predicted given the following parameters: stellar mass and source distance, disk inclination and position angle, phase center offsets in R.A. and Dec, the systemic velocity, and a mask size 8 . In some cases where both the disk front and back surfaces are well detected (e.g., HD 163296, IM Lup), we employ the aspect ratio of the emission surface (z/r) to describe the 3D geometry of the disk in the generation of the Keplerian mask. The parame-7 The final adopted uncertainty of the size also adds the beam radius in quadrature to represent possible systematic errors in the measurements. 8 The code to make Keplerian masks can be found at https:// github.com/richteague/keplerian mask; see also a detailed technical procedure for Keplarian mask making in the Appendix of Trapman et al. (2020a). For each disk, we create a series of masks with different sizes and determine the optimal mask as the one with which the measured gas radii no longer increase significantly (well within 1σ statistical uncertainty).   Figure 2). Each profile is normalized to its peak emission value. Each row has the same radial extent as labeled out in the x-axis labels of the left panels. The synthesized beam sizes are indicated by the horizontal bars. Rmm and RCO are marked as vertical dashed lines.
ters used for creating the mask for individual disks are summarized in Appendix Table B1. Figure 2 shows an example of the Keplerian mask applied to the GO Tau disk. With the extracted radial profile from the moment-zero map generated with the Keplerian mask, we then calculate the CO disk size (R CO ) and its uncertainty as described above for the dust disk size. For sources like GO Tau where cloud contamination dominates in some velocity channels, an additional mask (as shown in the bottom left panel of Figure 2) is applied to avoid using that portion of azimuthal angles in extracting the radial profile 9 . We usually use the half-side of the disk associated with less cloud con-ally exclude emission from central channels to estimate disk size for this demonstration).
The extraction of a radial emission profile from the moment-zero map usually assumes a flat disk geometry, since the exact flaring structures for most disks are unknown. For the well-studied HD 163296 disk, we have performed radial profile extraction considering the 12 CO emission surface ) 10 , and no significant differences in the derived CO sizes are observed when compared with the result based on a flat geometry. For disks with low to intermediate inclination angles, the flat disk geometry is a reasonable assumption. We have also tested with mock CO data with varying inclination angles. The comparison of high (70 • ) and low (10 • ) inclination models reveals an uncertainty of ∼10%. The derived radial profiles for both 12 CO and continuum emission for individual disks are shown in Figure 3, with the calculated disk sizes summarized in Table 2.
This effective radius definition depends on the distribution of emission within the disk. If the continuum emission is optically thin, R mm traces the outer boundary of the dust density distribution in the disk midplane (but see Rosotti et al. 2019). The 12 CO emission is usually optically thick across the disk and probes the gas temperature distribution in the elevated disk layers assuming local thermodynamic equilibrium. The momentzero maps could not directly reflect the mass distribution due to the high optical depth of CO line emission and the incorporation of the velocity information, for which line emission in the inner disk is associated with a broader velocity range thus higher weights. Nevertheless, the adopted R mm and R CO are reasonable metrics for the spatial extent of the disk dust and gas distributions, though we should keep in mind that the two disk sizes have different underlying physical meanings.

Overview of Disk Sizes
Following the procedure described in Section 3, we calculated R CO and R mm for 26 disks. build up a sample of 44 disks with both R CO and R mm measured in a similar manner from the images. The measurements are summarized in Table 2 and Figure 4. Recently, Sanchis et al. (2021) revisited the Lupus sample with new disk size measurements. They measured R CO from elliptical Gaussian or Nuker function fitting of the CO moment-zero maps, and R mm from modeling the continuum visibilities. Their R CO values are broadly consistent with those presented in Ansdell et al. (2018) (and adopted in this study), but the R mm derived from visibility modeling is generally 20-30% smaller than those derived in the image plane due to limited spatial resolution (see more detailed comparison in Appendix C). This is similar to the six very low mass Taurus stars, for which Kurtovic et al. (2021) measured R mm from uv-plane modeling profiles. For the other four Lupus disks where observations from DSHARP are available and adopted in this analysis, we find that our size measurements generally agree with Ansdell et al. (2018), except for the R CO of the IM Lup (Sz 82) disk. The clear detection of the diffuse outer gas emission in the deep DSHARP data leads to a larger R CO that is twice the size measured in Ansdell et al. (2018). Figure 4 compares R CO and R mm of this sample, which span ranges from ∼50-1000 au and ∼15-250 au, respectively. Disks with severe visible cloud absorption around the central velocity channels are marked as triangles, for which the R CO measurements are lower limits. The CO disks are all spatially resolved, and the majority of them are resolved in at least three indepen-    dent resolution elements (see Table B2 for the comparison of R CO and beam size). Though all the mm disks are also resolved, 9 of them are resolved in less than two independent resolution elements. R CO and R mm are positively correlated. Employing the Bayesian linear regression method of Kelly (2007) with its python implementation Linmix 12 , we find a best-fit relation of log R CO = (1.0±0.1)log R mm + (0.4±0.2), with 0.2 dex scatter of the correlation (as the standard deviation σ of an assumed Gaussian distribution around the mean relation). We also derive a best-fit relation in the linear space as R CO = (2.9±0.4) R mm -(7.5±39.5), demonstrating the average ratio of the two sizes. Though the two disk sizes are correlated, the derived scaling relations also suggest substantial dispersion. The CO disk of GO Tau is the most extended case, followed by the DM Tau disk. Both systems host large mm continuum disks, but are significantly offset upward from the bestfit scaling relation. In contrast, V1094 Sco has a CO disk (measured from the shallow Lupus survey, but the CO emission is spatially resolved), that is much smaller than suggested from its mm disk following the size relation.

Disk Size -Luminosity Plane
The disk size -luminosity scaling relation of R mm ∝ L 0.5 mm has been established in nearby star-forming regions (Tripathi et al. 2017;Andrews et al. 2018a). This correlation is also suggested to slightly flatten with evolution (Hendler et al. 2020). Though the origins of this relation are unclear, it provides a straightforward view of the typical disk characteristics. Figure 5a shows this sample in the L mm − R mm plane, which is well consistent with the established scaling relation 13 . The mm luminosities are normalized to a common distance of 140 pc and frequency of 225 GHz (1.3 mm), assuming a uniform diskintegrated spectral index of 2.3 (Andrews 2020;Tazzari et al. 2021). Data for a few disks were taken at 345 GHz and we assume R mm are similar at close frequencies of 225-345 GHz. The sample disks with R mm of ∼15-250 au almost cover the full range of dust disk sizes revealed from previous surveys (Andrews 2020). However, L mm of this sample spans from a few mJy to 300 mJy, mostly sampling the brighter end of the whole disk population, where more than half of known disks are fainter than 10 mJy at this wavelength.
The CO sizes and line luminosities are tightly correlated (Figure 5b). The line flux is calculated by integrating over the deprojected radial profile from the moment-zero map. The line fluxes at J = 3 − 2 transition in a few disks are converted into J = 2 − 1 fluxes using the ratio of the square of the line rest frequencies (assuming optically thick line emission), and assuming similar emitting radii for both transitions. We derive a best-fit relation of log R CO = (0.52±0.05) log L CO2−1 + (2.07±0.03), with a scatter of 0.15 dex. The slope is consistent with the fact that 12 CO emission is optically thick in disks and serves as a good proxy for the gas temperature of the emitting layer (see also Sanchis et al. 2021).

Disk Size -Host Star Properties
We consider here a search for any relationships between the disk sizes and host star properties, including mass (M * ), luminosity (L * ), and age (t * ). We find marginal relationships between disk sizes and M * . Regression analyses suggest best-fit relations of log R mm = (0.30±0.12) log M * + (1.96±0.05) with a scatter of 0.26 dex, and log R CO = (0.35±0.13) log M * + (2.40±0.06) with a scatter of 0.29 dex. These relationships are plotted in Figure 5c and (as also demonstrated in Andrews 2020), which is much steeper than what has been demonstrated by this sample. As we focus on bright disks, this inconsistency can be explained by the exclusion of a large number of faint and small disks around lower mass stars. Given the strong correlation between R CO and R mm , such bias in sample selection could also account for the derived relation between R CO and M * , and a steeper slope is expected when a more complete sample is considered.
Our linear regression analysis also suggests at most a weak relation between R CO and L * (R CO ∝ L 0.15±0.08 * , with a dispersion of 0.3 dex, see Figure 5e). This implies that the gas temperature in the disk surface layer may only weakly depend on stellar properties.
Although disks in this sample have stellar ages spanning 0.5 -20 Myr, no correlation between the age and the CO or continuum disk sizes is found (Figure 5f). This is probably because most of these disks are located in Taurus and Lupus and have ages of 1-3 Myr, accompanied with large uncertainties in individual age mea- surements. As the young disks (t <1 Myr) are all around K-type stars (see Figure 1), we also compare the CO disk sizes in different age ranges considering only earlier type stars (M * > 0.6 M ). By splitting the sample into t * < 1 Myr, 1 Myr < t * < 3 Myr, and t * > 3 Myr, We find comparable average CO sizes (also similar 1σ scatter of ∼0.2 dex) within the three bins.

CO-to-Continuum Disk Size Ratio
The sample presented in this work has universally larger R CO than R mm , with the R CO /R mm value spanning from 1.3 to the extreme case of 7.6. Most disks have R CO /R mm of 2-4, providing an average ratio of 2.9±1.2 ( Figure 6). Note that the measurements for 40% of the sample only provide a lower limit on R CO /R mm , because R CO is underestimated when cloud absorption is ob-served around the systemic velocity and/or R mm could be overestimated when the continuum emission is not well resolved (those disks are marked as grey triangles in Figure 6). Previous observations of individual bright disks and the ALMA-Lupus survey have demonstrated a general feature that the CO sizes are 2-3 time more extended than the mm continuum sizes (Isella et al. 2007;Panić et al. 2009;Andrews et al. 2012;Ansdell et al. 2018). Our analysis confirms and extends this finding in a sample with wider ranges of stellar and disk properties. The CX Tau disk, with R CO /R mm of 3.9, was considered as one of the most extreme cases ). In the sample presented here, we have six additional disks with higher R CO /R mm (another disk has ratio comparable to CX Tau). These disks with high R CO /R mm cover the full range of stellar mass and disk mm luminosity (see Figure 6).
We found no relationship between R CO /R mm with either stellar mass or disk mm luminosity. Theoretical calculation of dust evolution predicts higher drift velocity in disks around lower mass stars (Pinilla et al. 2013), thus larger R CO /R mm is expected for these disks. The lack of correlation suggests that pressure bumps may form in individual disks at certain radii or time that regulate the inward drift and break any potential relationships.
We explore here any relationship between the size ratio and the presence of substructures in disks. Within this sample, dust substructures have been identified in 31 disks, in the forms of inner cavities, gaps and rings, and/or spiral patterns (see corresponding references in Table 2). In the remaining 13 disks 14 , any substructure would be difficult to identify with the current data; four disks have large inclination angles (> 60 • ) and all but CX Tau in this subsample have R mm that is less than twice the data resolution. Figure 7 shows the histograms of R mm , R CO , and R CO /R mm by dividing the sample based on the observed presence of dust substructures. In the 31 disks of this sample where dust substructures have been identified, we found an average R CO /R mm of 3.0±1.3. Though the remaining 13 disks are in general more compact from both the dust and gas components, their average R CO /R mm (2.6±0.9) is comparable within uncertainties to that of the substructure group (see Figure 7). Those with high R CO /R mm include both extended substructure disks (GO Tau,  . The distribution of Rmm, RCO, RCO/Rmm for disks with observed dust substructures (grey) and disks that are featureless at the resolution of currently available data (blue). The latter group of disks are more compact, but their CO-tocontinuum size ratios are comparable to disks with substructures.
without large-scale substructures (CX Tau and CIDA 7). More specifically, the two largest CO disks, GO Tau and DM Tau, both surround M3 stars and host extended continuum disks with multiple gaps and rings detected (Long et al. 2018;Hashimoto et al. 2021). Around the same type of host star, the continuum disk of CX Tau only extends to 30 au and remains featureless with a high resolution of 5 au ). Our findings suggest that R CO /R mm may not strongly depend on the detailed dust distribution in disks (or every disk is highly structured and R mm is largely determined by where the last pressure bump can be built in the disk, so does R CO /R mm ). The evolution of R CO may also differ in different systems so that R CO /R mm is not a simple reflection of the dust radial drift efficiency.

DISCUSSION
In this section, we first discuss if CO disk sizes can reveal the dominant disk evolution mechanism, assuming R CO well traces the disk gas component. We then explore how CO-to-continuum size ratios probe the evolution and distribution of disk solids. Finally, we employ the GO Tau disk to demonstrate the influence of a large gas disk on dust evolution.

Gas Disk Evolution
In the classical theory of disk viscous evolution, turbulence transports angular momentum outward and enables disk material to be accreted onto the star (e.g., Hartmann et al. 2016). This re-distribution of angular momentum leads to the growth of the disk size. The expansion rate largely depends on the turbulent viscosity. Trapman et al. (2020b) recently explored how the disk sizes (using the same definition as 90% fractional radius) measured from 12 CO emission vary with time in viscous evolution models using a simplified α−prescription for the viscosity (Shakura & Sunyaev 1973;Lynden-Bell & Pringle 1974). This model involves three key parameters: viscosity coefficient α, initial disk mass M d,0 , and initial characteristic disk radius R c,0 . For a viscous disk, M d,0 is related to the stellar accretion rate (Hartmann et al. 1998), which scales with the central stellar mass. Therefore, we considered two sets of models with stellar masses of 0.32 and 1 M , and adopted the average stellar accretion rate among the corresponding host stars.
A comparison of the CO disk sizes derived here with model grids from Trapman et al. (2020b) suggests that most cases can be explained by viscous evolution models with low α in the range of 10 −4 -10 −3 and small R c,0 =10 au (Figure 8). The three most extended gas disks in this sample (GO Tau, DM Tau, and IM Lup) are better described by models with high α (∼0.01), though the large R CO of IM Lup was suggested to result from external photoevaporation in a weak radiation field (Haworth et al. 2017). Measuring the spectral line broadening provides observational constraints on disk turbulence, which suggests that weak turbulence (α < 10 −3 ) might be common in disks (Guilloteau et al. 2012;Flaherty et al. 2015Flaherty et al. , 2018Teague et al. 2018). DM Tau is the only case with measurable turbulence in the outer disk (Flaherty et al. 2020); such measurements are not yet available for GO Tau and IM Lup. An alternative solution is to start with larger initial sizes. Disk models with R c,0 = 50 au can spread to R CO = 500 au with a lower α of 10 −3 by 1 Myr (Trapman et al. 2020b). Future constraints on turbulence from nonthermal line broadening in these extended disks will potentially distinguish the two scenarios.
The magneto-rotational instability (MRI, Balbus & Hawley 1998) is taken as the leading mechanism in generating the needed turbulence for the viscous transport of disk angular momentum. However, the full operation of MRI throughout the disk is often questioned due to the weak ionization conditions in large disk areas (e.g., Turner et al. 2014). Numerical simulations with proper treatments of non-ideal MHD effects demonstrate that MRI turbulence is largely suppressed in the cold, low-  ionization disk regions (e.g., Bai & Stone 2013;Gressel et al. 2015). The magnetized disk wind concept (Blandford & Payne 1982) has thus re-emerged as a promising alternative (e.g., Bai et al. 2016;Lesur 2021). In this wind-driven accretion scenario, angular momentum is not transported within the disk but directly removed through MHD disk winds. Therefore the sizes of gas disks need not expand to enable mass accretion. Following the α-framework for viscous evolution, a simple parameterized description of disk evolution for a MHD disk wind was recently provided by Tabone et al. (2021a). This introduced a similar dimensionless param-eter α DW that describes the wind torque and relates to the local accretion rate driven by the wind. This framework therefore controls disk evolution with three key parameters: α DW , initial disk mass M d,0 , and initial characteristic disk radii R c,0 (the accretion timescale t acc,0 is related to α DW and R c,0 , which will be used in the following discussion). Based on the analytical solution from Tabone et al. (2021a then examined how R CO evolves in the pure disk wind scenario (ignoring the viscous term). When considering a constant α DW , R CO gradually declines with time (between 0.1 to 10 Myr, as shown in the left panel of Figure 9), which is in direct contrast to the case of viscous evolution. The declining rate of R CO depends on the choice of M d,0 and R c,0 , but in general models with higher M d,0 and larger R c,0 result in larger R CO at specific time steps. Therefore, disk wind models with different M d,0 (from 0.0001 to 0.1 M ) and R c,0 (from 20 to 90 au) can explain most of our disk sample. As R c,0 =90 au is the most extreme case considered in the models of Trapman et al. (2021) (guided by the Lupus survey, Ansdell et al. 2018), the few very extended disks and the disks around older systems in the sample are difficult to account for without even larger initial disk sizes.
The time evolution of R CO in this disk wind scenario largely depends on the rate at which the disk surface density decreases. If we consider a time-dependent α DW that scales with Σ c (t) −1 , the decrease of R CO becomes much shallower until a sharp drop occurs when the disk starts to dissipate (at 2t acc,0 , Tabone et al. 2021a). This difference is demonstrated in the two panels of Figure 9. In this time-dependent case, the evolution of R CO is mostly determined by the ratio between t acc,0 and t * . Therefore, old disks with large R CO of 200-300 au may be initialized with large t acc,0 and evolve with varying α DW . Considering a distribution of t acc,0 in any individual cluster that follows the evolution of disk fraction with cluster age (e.g., Fedele et al. 2010), Tabone et al. (2021b) show that the correlation between mass accretion rate and disk mass, as well as its large scatter, can be well reproduced by MHD disk wind models with time-dependent α DW . It is likely that these old and large disks represent the long-lived outliers (those with long initial accretion timescale) in individual starforming regions (e.g, TW Hya).
The comparisons above suggest that the measured CO sizes for this disk sample can be explained by either viscous models or MHD disk wind models, assuming varying initial conditions for individual disks. As shown in Figure 8 and 9, distinguishing between the two scenarios requires a more complete sample of older systems (see also Figure 1). Considering the typical disk dispersal timescale, the characterization of such disks will need deep observations (the ongoing ALMA Cycle 8 Large Program AGE-PRO will hopefully provide such measurements). The inclusion of a large sample is also necessary to account for any possible biases in the distribution of initial conditions.
With uniformly determined CO sizes, the sample of Class II disks presented here do not show any trend of R CO with evolution (with stellar ages from 0.5 to 20 Myr, Figure 5f). Our small number statistics (especially for older disks) limit the evidence for disk vis-cous evolution during the Class II phase. However, in a comparison of gas disk sizes between the younger Class 0/I and more evolved Class II disks, Najita & Bergin (2018) found that gas disk sizes of Class II disks are systematically larger and explained this as a result of viscous spreading. Given the lack of large disks with R CO > 800 au at t * < 0.5 Myr (R CO,Class 0/I ∼ 50-300 au, e.g., Harsono et al. 2014) 15 , MHD disk wind models alone would be challenging to explain the extended gas distribution in some Class II objects. We speculate that disks evolve in a hybrid mode, where substantial growth of sizes may occur early on for some systems to account for the overall size differences at different evolutionary stages. For individual disks, different mechanisms may dominate the evolution at different times, depending on the physical conditions in the disk and of the large-scale environment. It also remains unclear how the mass and angular momentum transfer from the envelope impacts the disk size evolution, as the envelope is not treated in the current model framework. One important caveat of the comparison in Najita & Bergin (2018) lies in the usage of various gas tracers and methods in calculating the disk size with literature results directly adopted. In addition, the difficulties in disentangling the disk and envelop emission result in only a small number of Class I sources with well-determined (Keplerian-rotating) disk radii. High spatial and spectral resolution ALMA observations towards these younger systems in the coming years will certainly enlarge the sample. A more robust test of viscous spreading would then be possible with sufficient lever arm in stellar age as well as comparable stellar mass distributions across the age spectrum.

R CO /R mm Tracing Dust Evolution
Our analysis shows that R CO /R mm can vary from 1.3 to 7.6, with an average ratio of 2.9±1.2 (see Sec 4.4). Such size differences have been widely attributed to two effects: emission optical depth and dust evolution. The 12 CO line optical depth is always much higher than the nearby continuum emission (e.g., Beckwith & Sargent 1993). The observed size discrepancy can simply be a manifestation of this optical depth difference, since the thicker line emission is easier to detect at low densities in the outer disk .
From the dust evolution perspective, dust particles in disks will grow to larger sizes, and the shorter grain growth timescale at smaller radial distances results in a concentration of larger grains in the inner disk (e.g., Testi et al. 2014). Meanwhile, once grains in the outer disk reach some critical sizes, the gas drag imposed by gas-dust rotational velocity difference will push the large grains to move inward toward higher pressure regions (Weidenschilling 1977). The combination of the two processes means that large grains preferentially accumulate in the inner disk. The millimeter continuum emission tracing these particles will therefore appear much more compact than the gas disk. In addition, Rosotti et al. (2019) have suggested that small grains (with sizes smaller than 100 µm), which are still well-coupled with the gas and present at large radial distances, could not produce sufficient mm continuum emission due to a sharp dust opacity drop.
The two effects of optical depth and dust evolution act in concert. However, their relative importance can not be easily retrieved without detailed knowledge of the gas surface density and other disk conditions that affect the grain growth and drift efficiencies (e.g., viscosity). Thermochemical models for the HD163296 disk have demonstrated that its R CO and R mm difference can be largely explained by the optical depth effect . Using a grid of thermochemical models with varying disk conditions, Trapman et al. (2019) suggest that R CO /R mm > 4 is a clear sign of dust evolution in action, regardless of the detailed stellar/disk properties. About 18% (8/44) of the disks in our sample have size ratio around or above 4. Presumably these disks have experienced substantial grain growth and radial drift. We recall that these high ratio systems do not occupy any preferred parameter space of stellar mass and disk millimeter brightness. Therefore, it remains unclear how R CO /R mm can be better employed as a generic diagnostic of dust evolution. Through quantitative comparisons to thermochemical models without dust radial drift for a few individual systems, Trapman et al. (2020a) also revealed that dust evolution is required in five Lupus disks to account for the observed disk size differences with R CO /R mm of 1.8-3; all of them are included in this sample.
Based on theoretical models, radial drift is expected to proceed in a short timescale (e.g., Takeuchi & Lin 2002, Brauer et al. 2008), which will lead to a very rapid shrinking of R mm . Toci et al. (2021) recently explored the secular evolution of dust and gas disk sizes with a set of model grids, considering grain growth, radial drift, and disk viscous evolution under smooth gas distributions. In model disks with R c =10 au surrounding solar-type stars, R CO reaches 50-300 au after 1-2 Myr, while R mm shrinks to 10-30 au, leading to universally CO line model data Figure 10. The comparison of GO Tau data radial profiles with model profiles convolved with our data resolution (blue for the CO emission and orange for the mm continuum emission). This fiducial model with a smooth gas distribution is directly adopted from Toci et al. (2021), selected to roughly match the large CO emission size of GO Tau. The grey dashed line shows the revised model with two pressure bumps added on top of the fiducial model aiming at reproducing the two observed dust rings (the corresponding CO profile is almost identical to the fiducial model).
high R CO /R mm (> 5). Though R CO in these models match with the peak distribution of the disk sample presented in this study (see Figure 7), R mm are too small compared to observations. The presence of local gas pressure maxima, due to either the influence of planets or fluid instabilities, will stop or slow down the dust inward drift (e.g., Pinilla et al. 2012). About two thirds of the disks in this sample show visible dust substructures, which may trace pressure bumps that help sustain large grains at large radial distances. The remaining sample, though with highly comparable R CO to models, show a wide distribution of R mm (15 ∼ 150 au) and low size ratios mainly clustering around 2 (Figure 7). None of the observations toward these disks are capable of detecting substructures with sizes comparable to the gas pressure scale height in their outer disk. Substructures may hide in current data for these disks given their R mm and R CO /R mm .

The Most Extended CO disk of GO Tau
The GO Tau disk has the largest R CO within this sample, approaching ∼1000 au. Viscous spreading likely plays a leading role in the past evolution for such extended disks (including e.g., DM Tau). Future measurements of the line broadening will help assess the disk turbulence level and therefore constrain some initial disk conditions. We also note that both GO Tau and DM Tau have comparably faint optical [O I] 6300Å line emission and without high velocity line component, indicating the absence of strong disk winds (Simon et al. 2016).
Such a large gas disk will affect the dust evolution in GO Tau. As a point of comparison, we selected a model from the Toci et al. (2021) grids that has a similar R CO to the GO Tau disk. We also re-ran the simulation with a perturbed version of this model, with two pressure bumps added to mimic the morphology of the mm continuum emission rings at 73 and 110 au, respectively. As seen in Figure 10, the radial profile from the pressure bump model matches well with the data. However, even without gas pressure bumps, the simulated mm continuum emission extends beyond the first prominent dust ring, preserving a large R mm . Therefore, a large R mm is not necessarily set by the presence of pressure bumps. It could also result from the extended gas distribution (or a large initial disk size), where grain growth takes longer and the drifting particles have just reached the current mm disk outer edge. The models presented here are not fine-tuned to match the details of the GO Tau disk, but are used for a qualitative demonstration. While a low R CO /R mm ratio suggests the presence of pressure bumps within the dust disk, a high ratio, like the case of GO Tau, only signifies prominent dust evolution outside the dust disk. This latter scenario could also be applied to the disks of CX Tau and CIDA 7, both of which have high R CO /R mm ratio. Their lack of dust substructures can be explained if the growth and subsequent radial migration of dust grains in these systems take quite a long time, though substructures may be present but harder to identify in these compact disks.

SUMMARY
We presented a joint analysis of disk sizes for a sample of 44 protoplanetary disks around central stars with masses of 0.15-2 M . The gas and dust disk sizes were calculated as the 90% flux fractional radii from the 12 CO line and millimeter continuum radial profiles. This sample covers a wide range in R mm of 15-250 au and R CO of 50-1000 au. Based on the distribution of the sample in the R mm -L mm plane, it is more representative of the brighter, larger disk population. We then explored how R CO and R CO /R mm are related to various stellar and disk properties, and considered their behavior in the context of disk evolution mechanisms. Our main results are summarized as follows.
1. R CO in this sample shows no evolution with stellar age within 0.5-20 Myr. Both viscous and MHD disk wind models can explain the sizes of individual disks in this sample, when considering varying initial conditions (e.g., initial disk characteristic radius, disk mass, accretion timescale). Disks with very large R CO are more readily explained in the context of the viscous models. Though the lack of apparent R CO evolution in this sample could be more consistent with the wind-driven accretion scenario, a more robust test would require a larger population study, especially towards both younger (Class I) objects and more evolved systems.
2. The disks in this sample have universally larger R CO than R mm . The measured R CO /R mm vary from 1.3 to 7.6, with an average value of 2.9±1.2. In 8/44 disks, high R CO /R mm values (>∼4) suggest that substantial grain growth and radial drift have already occurred. However, this sub-sample does not show any preference in stellar and disk properties. A significant fraction of the sample disks (31/44, 70%) exhibit dust substructures, which could mitigate dust radial drift and naturally explain their low R CO /R mm . The low disk size ratios in many currently featureless disks may then imply the presence of yet-unseen disk substructures.
3. Among this sample, the GO Tau disk stands out as an extreme outlier, with R CO extending to 1000 au and R CO /R mm of 7.6. There are two possible explanations for such an extended gas distribution: 1) a disk with an initial characteristic radius R c,0 =10 au experienced significant viscous spreading with α ∼ 0.01; 2) the disk was born large (R c,0 50 au). Measuring the non-thermal spectral line broadening that constrains the disk turbulence would be a promising approach to discern the two possibilities. Though GO Tau has two dust rings near the outer continuum edge, these substructures are not the prerequisite to sustain its large dust disk (∼130 au) due to the extended gas emission. From the dust evolution perspective, smooth disks with a high R CO /R mm are not prohibited if the gas disks were initially large enough.
We are still at the early stage in exploring the gas content of planet-forming disks and their roles in tracing disk evolution, especially considering the small number of disks that have high quality molecular line data. Improved investments in ALMA gas observations are fundamental in addressing many key questions of disk evolution and planet formation. The channel maps of the 12 CO J = 2 − 1 emission for DL Tau and UZ Tau from the SMA observations are shown in Figure A1 and Figure A2, respectively. Both disks suffer from severe cloud contamination near the systemic velocity (see Figure A3 for the extracted line profile).  Figure A2. Channel maps of 12 CO emission for UZ Tau stellar system, with Keplerian masks applied in each channel. Contours are starting from 5σ for the 1.3 mm continuum image. The close binary to the west (UZ Tau Wab) are also shown in continuum contours in the field, for which the 12 CO emission around this binary is not detected in our SMA observation. Cloud absorption is visible around central channels. Table B1 summarizes the mask size (R mask ) and the CO emission height (z/r) for individual disks when creating the Keplerian mask, which will ensure the inclusion of the full CO emission region. We also specify in Table B1 the status of cloud contamination and the azimuthal angle adopted for radial profile extraction in each disk.

B. NOTES ON KEPLERIAN MASKS AND CO SIZE MEASURMENTS FOR INDIVIDUAL DISKS
In Figure B4, we use HD 163296 to demonstrate how cloud absorption in central velocity channels affect the gas disk size measurement. Considering a typical cloud emission line width of 1-2 km s −1 (see Table B1), we calculate R CO by excluding a similar velocity range in the central channels. Our experiment shows that the disk sizes can be underestimated by 10-30%.

C. COMPARISON OF LUPUS DISKS
We adopted the mm continuum and CO disk size measurements from Ansdell et al. (2018) for 17 Lupus disks. Recently, Sanchis et al. (2021) provided new size measurements for this sample, but employing different methods. They calculated R CO through elliptical Gaussian or Nuker profile fitting on moment-zero maps, and R mm through Nuker profile fitting in the uv-plane. Figure C5 shows the comparison of disk sizes with their ratios from both studies. As Sanchis et al. (2021) only reported R CO,68% , we calculated the R CO,90% based on the model profile parameters provided in the paper to enable a fair comparison. The higher CO-to-mm size ratios in Sanchis et al. (2021) are largely attributed to the smaller R mm as estimated by visibility fitting, especially for sources that are only resolved in 2-3 beams in the mm images.  Figure D6 shows the distribution of R CO /R mm with R mm and stellar age. Combined with Figure 6, We do not find any correlation between the disk size ratios and stellar/disk properties.     (4) and (5) are the mask radius and the assumed CO emission height, respectively, for the generation of Keplerian mask. In IM Lup, AS 209, and DoAr 25, the mask is made with a 20% higher stellar mass, which matches better the line emission at high velocity channels. The parenthesis in column (6) denotes the cloud velocity range. Ansdell+2018 Sanchis+2021 Figure C5. Left: The comparison of disk CO-to-mm continuum size ratio for Lupus disks from Ansdell et al. (2018) and Sanchis et al. (2021). Right: The comparison of disk CO-to-mm continuum size ratio for Lupus disks with the mm continuum size. The vertical dashed line marks the typical beam size of the ALMA Lupus survey data.  Figure D6. The distribution of RCO/Rmm with Rmm (left) and stellar age (right). Grey triangles represent disks that either have cloud contamination around systemic velocities or the dust disk is only resolved in less than 2×beam sizes, in both cases RCO/Rmm is likely underestimated.