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

HAL
open science

# The ALMA Survey of Gas Evolution of PROtoplanetary Disks (AGE-PRO). V. Protoplanetary Gas Disk Masses

Leon Trapman, Ke Zhang, Giovanni Rosotti, Paola Pinilla, Benoît Tabone, Ilaria Pascucci, Carolina Agurto-Gangas, Rossella Anania, John Carpenter, Lucas Cieza, et al.

► To cite this version:

Leon Trapman, Ke Zhang, Giovanni Rosotti, Paola Pinilla, Benoît Tabone, et al. The ALMA Survey of Gas Evolution of PROtoplanetary Disks (AGE-PRO). V. Protoplanetary Gas Disk Masses. The Astrophysical Journal, 2025, 989 (1), pp.5. 10.3847/1538-4357/adcd6e. hal-05411947

HAL Id: hal-05411947

https://hal.science/hal-05411947v1

Submitted on 17 Dec 2025

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

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

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

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

# The ALMA Survey of Gas Evolution of PROtoplanetary Disks (AGE-PRO). V. Protoplanetary Gas Disk Masses

Leon Trapman¹, Ke Zhang¹, Giovanni P. Rosotti², Paola Pinilla³, Benoît Tabone⁴, Ilaria Pascucci⁵,
Carolina Agurto-Gangas⁶, Rossella Anania², John Carpenter⁷, Lucas A. Cieza⁸, Dingshan Deng⁵, Camilo González-Ruilova⁸,9,10, Michiel R. Hogerheijde¹¹,12, Nicolás T. Kurtovic¹³,14, Aleksandra Kuznetsova¹⁵,
James Miley⁹,10,16, Laura M. Pérez⁶, Dary A. Ruiz-Rodriguez¹⁷, Kamber Schwarz¹⁴, Anibal Sierra³,6, Estephani TorresVillanueva¹, and Miguel Vioque¹⁸,19

¹Department of Astronomy, University of Wisconsin-Madison, 475 N Charter Street, Madison, WI 53706, USA; ltrapman@wisc.edu

²Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, I-20133 Milano, Italy

³Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK

⁴Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France

⁵Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA

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

⁷Joint ALMA Observatory, Avenida Alonso de Córdova 3107, Vitacura, Santiago, Chile

⁸Instituto de Estudios Astrofísicos, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile

⁹Millennium Nucleus on Young Exoplanets and their Moons (YEMS), Chile

¹⁰Center for Interdisciplinary Research in Astrophysics and Space Exploration (CIRAS), Universidad de Santiago, Chile

¹¹Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands

¹²Anton Pannekoek Institute for Astronomy, University of Amsterdam, The Netherlands

¹³Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse 1, D-85748 Garching, Germany

¹⁴Max-Planck-Institut fur Astronomie (MPIA), Königstuhl 17, 69117 Heidelberg, Germany

¹⁵Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10025, USA

¹⁶Departamento de Física, Universidad de Santiago de Chile, Av. Victor Jara 3659, Santiago, Chile

¹⁷National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

¹⁸European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei München, Germany

¹⁹Joint ALMA Observatory, Alonso de Córdova 3107, Vitacura, Santiago 763-0355, Chile

Received 2024 August 30; revised 2025 March 31; accepted 2025 April 14; published 2025 July 31

## Abstract

The evolution of the gas mass of planet-forming disks around young stars is crucial for our understanding of planet formation, yet it has proven hard to constrain observationally, due both to the difficulties of measuring gas masses and the lack of a homogeneous sample. Here we present a large grid of thermochemical models that we use to measure protoplanetary gas disk masses of AGE-PRO, the Atacama Large Millimeter/submillimeter Array survey of Gas Evolution in PROtoplanetary disks. AGE-PRO covers a sample of 30 disks around similar spectral type (M3-K6) stars with ages between 0.1 and 10 Myr. Our approach is to simultaneously fit observations of CO isotopologues and N₂H⁺, a complementary molecule produced when CO freezes out. We find that the median gas mass of the three regions decreases over time, from 7.0±4.4 to 10⁻³ M☉ in Ophiuchus (≤1 Myr) to 9.4±5.4 to 10⁻⁴ M☉ for Lupus (~1–3 Myr) and 6.8±5.1 to 10⁻⁴ M☉ for Upper Sco (~2–6 Myr), with ~1 dex scatter in gas mass in each region. We note that the gas mass distributions for Lupus and Upper Sco look very similar, which could be due to survivorship bias for the latter. The median bulk CO abundance in the CO emitting layer is found to be a factor ~10 lower than the interstellar medium value but does not significantly change between Lupus and Upper Sco. From Lupus to Upper Sco, the median gas-to-dust mass ratio increases by a factor ~3 from ~40 to ~120, suggesting efficient inward pebble drift and/or the formation of planetesimals.

*Unified Astronomy Thesaurus concepts: Astrochemistry (75); Protoplanetary disks (1300); Radiative transfer (1335)*

*Materials only available in the online version of record: figure set*

## 1. Introduction

The gas mass of planet-forming disks around young stars is a vital yet elusive part of planet formation theories and very much needed for our understanding of observed planetary systems (e.g., A. Morbidelli & S. N. Raymond 2016; J. Dražkowska et al. 2023). In its simplest form, the disk gas mass represents the total reservoir of material available for the formation of gas giants. Knowing the evolution of protoplanetary gas disk masses

is also highly desirable, as it sets how much time giant planets have to accrete their atmosphere. In addition, these giant planets can also influence the formation and properties of terrestrial planets (e.g., B. Bitsch & A. Izidoro 2023). But gas is also important for all steps of planet formation. The gas density and local gas-to-dust mass ratio regulate the dynamics of dust grains and larger bodies that could eventually grow to form rocky planets and giant planet cores. Dust growth, settling rates, and the speed at which larger grains radially drift inward all depend on the gas content, as well as the migration of young planets that have already formed (see, e.g., recent reviews by J. Dražkowska et al. 2023; S. Paardekooper et al. 2023; T. Birnstiel 2024).

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Measuring gas disk masses has proven difficult (e.g., E. A. Bergin & J. P. Williams 2017), primarily because the main constituent of the gas, molecular hydrogen ($H_2$), is both difficult to observe in disks and the few available observations do not trace the bulk gas in protoplanetary disks. $H_2$ is a light symmetric molecule that lacks a permanent dipole, limiting its rotational lines to quadrupole transitions ($\Delta J = 2$). Even the lowest transition of $H_2$, $S(0)$ ($J = 2-0$) at 28 µm, requires exciting the $J = 2$ level at 549.2 K, which only happens to any appreciable degree for gas temperatures above 100 K (e.g., W. F. Thi et al. 2001; A. Carmona et al. 2011). However, the bulk of the gas in protoplanetary disks resides at temperatures below ~30 K (e.g., C. Pinte et al. 2018; K. R. Schwarz et al. 2018; C. J. Law et al. 2021, 2022; T. Paneque-Carreño et al. 2023), meaning that $H_2$ emission does not trace the bulk of the mass (e.g., I. Pascucci et al. 2013).

The gas mass therefore has to be measured using indirect tracers. By far the most numerous have been disk masses derived from observations of the (sub)millimeter continuum emission of dust grains in protoplanetary disks (e.g., S. V. W. Beckwith et al. 1990; A. Dutrey et al. 1996; J. P. Williams et al. 2005; S. M. Andrews et al. 2013; M. Ansdell et al. 2016; S. A. Barenfeld et al. 2016; I. Pascucci et al. 2016; L. A. Cieza et al. 2019; J. J. Tobin et al. 2020; S. E. van Terwisga et al. 2022). Measuring a disk mass from a millimeter continuum flux relies on assumptions of the dust temperature, optical depth, dust opacity, and grain size distribution as well as the inclusion of scattering effects, each with their uncertainties (e.g., J. P. Williams & L. A. Cieza 2011; Z. Zhu et al. 2019; A. Miotello et al. 2023). But most importantly, to determine the gas mass, a gas-to-dust mass ratio has to be assumed. It is common to use the interstellar medium (ISM) value of 100, as this is likely the ratio at the formation of the disk, but dust evolution models predict that the gas-to-dust mass ratio would vary across the disk as it evolves due to processes such as photoevaporation, mass loss due to disk winds, grain growth, dust settling, and radial drift (e.g., T. Takeuchi et al. 2005; T. Birnstiel et al. 2012; P. Pinilla et al. 2012; L. Testi et al. 2014; G. P. Rosotti et al. 2019). Moreover, the presence of substructures in millimeter continuum observations of most disks strongly suggests that the gas-to-dust mass ratio at least varies radially in protoplanetary disks (e.g., S. M. Andrews et al. 2018; J. Huang et al. 2018; F. Long et al. 2019; L. A. Cieza et al. 2021; see S. M. Andrews 2020; J. Bae et al. 2023 for recent reviews).

Among the gaseous species, hydrogen deuteride (HD) has been shown as a promising indirect tracer of the gas mass. Due to their chemical similarities, HD closely follows $H_2$ (see L. Trapman et al. 2017); but contrary to $H_2$, HD has a small dipole moment and requires much lower temperatures (~30–50 K; e.g., E. A. Bergin et al. 2013; L. Trapman et al. 2017; J. K. Calahan et al. 2021) to excite and produce significant emission. HD $J = 1-0$, emitted at 112.07 µm, was detected with Herschel toward three protoplanetary disks, TW Hya, DM Tau, and GM Aur (E. A. Bergin et al. 2013; M. K. McClure et al. 2016). By using these detections, gas masses were derived for each disk (e.g., E. A. Bergin et al. 2013; M. K. McClure et al. 2016; L. Trapman et al. 2017, 2022b; J. K. Calahan et al. 2021; K. R. Schwarz et al. 2021). However, with the end of the Herschel mission,

there is currently no far-infrared observatory available to observe HD in more disks.

Gas masses are therefore most commonly measured from the rotational emission of carbon monoxide (CO), specifically from its more optically thin isotopologues $^{13}$CO and $C^{18}$O (e.g., J. P. Williams & W. M. J. Best 2014; F. Long et al. 2017; A. Miotello et al. 2017). CO is the second most abundant molecule in the gas of protoplanetary disks, and its emission is bright at millimeter wavelengths. Using CO as a gas mass tracer requires knowledge of the CO abundance with respect to $H_2$ ($x_{CO}$). CO is a relatively chemically stable molecule that is expected to have an abundance of $x_{CO} \approx 10^{-4}$ in the warm molecular layer of protoplanetary disks (e.g., Y. Aikawa et al. 2002). Two processes reduce $x_{CO}$, namely photodissocation in the surface layers of the disk and freeze-out close to the midplane. Both of these processes have been extensively studied in the context of disks (e.g., Y. Aikawa et al. 2002; R. Visser et al. 2009; A. Miotello et al. 2014; M. Ruaud & U. Gorti 2019), and by including them in thermochemical disk models, we can account for these processes and derive a relation between $^{13}$CO and $C^{18}$O fluxes and disk mass (e.g., J. P. Williams & W. M. J. Best 2014; A. Miotello et al. 2016; D. Deng et al. 2023).

Applying this method to derive gas masses for a large number of disks, it has become increasingly clear that are other processes affecting CO in protoplanetary disks that remain unaccounted for (e.g., C. Favre et al. 2013; F. Long et al. 2017; A. Miotello et al. 2017). Even after accounting for (isotope-selective) photodissociation and freeze-out, gas masses derived from CO are found to be ~5–100× lower than those derived independently from HD (e.g., C. Favre et al. 2013; M. Kama et al. 2016; M. K. McClure et al. 2016; K. R. Schwarz et al. 2016; L. Trapman et al. 2017; J. K. Calahan et al. 2021). Similarly, gas-to-dust mass ratios from CO-based gas masses are found to be low, ~1–10, compared to the expected value of ~100 the disk should have inherited from the ISM (e.g., M. Ansdell et al. 2016; F. Long et al. 2017; A. Miotello et al. 2017). Suggestions for these unaccounted for processes that reduce $x_{CO}$ include the chemical conversion of CO into more complex, less volatile species (e.g., Y. Aikawa et al. 1997; K. Furuya & Y. Aikawa 2014; M. Yu et al. 2016, 2017; A. D. Bosman et al. 2018; K. R. Schwarz et al. 2018; K. Furuya et al. 2022; M. Ruaud et al. 2022) and locking up CO into larger dust bodies that settle toward the midplane (e.g., E. Bergin et al. 2010; E. A. Bergin et al. 2016; M. Kama et al. 2016; S. Krijt et al. 2018, 2020; D. Powell et al. 2022; E. Van Clepper et al. 2022).

Models using a self-consistent disk structure, and accounting for the photodissociation, freeze-out, and CO conversion into $CO_2$ ice can reproduce the observed line fluxes of CO isotopologues as well as the main dissociation product of CO, i.e., atomic carbon, without significantly reducing the gas-to-dust mass ratio or the bulk carbon abundance (M. Ruaud et al. 2022; D. Deng et al. 2023; I. Pascucci et al. 2023).

Rather than focusing on one or more of these processes, each with their own uncertainties and assumptions, recent work has proposed a different strategy, namely constraining the bulk $x_{CO}$ of disk using observations of diazenylium ($N_2H^+$; D. E. Anderson et al. 2019, 2022; L. Trapman et al. 2022b). $N_2H^+$ is formed through proton transfer of $H_3^-$ to $N_2$. If CO is present, it competes with $N_2$ for the available $H_3^-$ to form $HCO^+$. Additionally, the main destruction pathway of $N_2H^+$ is
through reacting with CO. As a result of both of these factors, N₂H⁺ is a molecule that is particularly sensitive to the gas-phase CO. L. Trapman et al. (2022b) compared the gas masses derived from the CO isotopologue and N₂H⁺ method with the masses from HD (1-0) line fluxes in three protoplanetary disks. The results demonstrated that, even including the two orders-of-magnitude uncertainty on the ionization rate, the two methods give consistent gas mass estimates. Similarly, D. E. Anderson et al. (2019, 2022) demonstrated (using an independent thermochemical code) that the N₂H⁺-to-CO isotopologue line flux ratio is highly sensitive to the gas-phase xCO and the inclusion of N₂H⁺ significantly improves the gas mass accuracy.

Based on previous successful benchmark tests, in this work we apply this CO and N₂H⁺ method on new observations of the AGE-PRO Atacama Large Millimeter/submillimeter Array (ALMA) large program (K. Zhang et al. 2025) to measure gas masses, gas-to-dust mass ratios and bulk CO abundances for the 20 Lupus and Upper Sco Class II disks in the AGE-PRO sample. The 10 Ophiuchus Class I disks in the AGE-PRO sample are young enough that the aforementioned processes affecting CO have not had the time to run their course (e.g., K. Zhang et al. 2020). When measuring their gas masses, we can therefore limit ourselves to photodissociation of CO in the disk surface and freeze-out in the midplane when interpreting the CO isotopologue observations of Ophiuchus, without having to include N₂H⁺ to constrain their bulk CO abundances.

We run a large grid of thermochemical models, similar in spirit to the grids in J. P. Williams & W. M. J. Best (2014) and A. Miotello et al. (2016) but now also covering a range of CO abundances. The AGE-PRO program observed ¹²CO, ¹³CO, C¹⁸O, and C¹⁷O for its sources in Ophiuchus (D. A. Ruiz-Rodriguez et al. 2025) and ¹²CO, ¹³CO, C¹⁸O, and N₂H⁺ for its sources in Lupus (D. Deng et al. 2025) and Upper Sco (C. Agurto-Gangas et al. 2025). We extract synthetic observations of these lines from our models and compare them to the AGE-PRO observations to derive gas masses. Combining this with dust masses derived from the millimeter continuum, we can provide gas-to-dust mass ratios for the whole AGE-PRO sample. These newly derived gas masses are an essential input for the disk evolution models of B. Tabone et al. (2025).

The structure of this paper is as follows: In Section 2 we describe the model grid. In Section 3 we first examine the simulated observables of the model after which we compare them to the observations using a Markov Chain Monte Carlo (MCMC) method to derive gas masses and their associated uncertainties. In Section 4 we compare our gas masses to stellar accretion rates and discuss the resulting disk lifetimes. We also look for and discuss correlations between gas mass and other disk and stellar properties. Finally, we discuss caveats that could affect the derived gas masses, and we summarize our conclusions in Section 5.

## 2. DALI Models

In this work, we used the thermochemical code DALI (S. Bruderer et al. 2012; S. Bruderer 2013). For a given density structure and stellar spectrum, DALI self-consistently calculates the chemical and thermal structure of the disk. First, the dust temperature and internal radiation field are calculated using a 2D Monte Carlo method to solve the radiative transfer

equation. Next, DALI computes the time-dependent chemistry, calculates the atomic and molecular excitation levels, and balances the heating and cooling processes to determine the gas temperature. Due to the inter-dependencies of these processes, they are calculated iteratively until a self-consistent solution is found. The resulting temperature structure is then used to run the larger chemical network presented in A. Miotello et al. (2014), which includes CO isotope-selective chemistry and photodissociation. We note that in this final step, there is no further iteration on the disk temperature. For a more detailed description of DALI, see Appendix A in S. Bruderer et al. (2012).

Both chemical networks mentioned previously do not include the chemical conversion of CO into more complex, less volatile species. Instead, we opt to mimic these processes in a post-processing step with the peak CO abundance (xCO) as a free parameter. Here we calculate the N₂H⁺ chemistry using the chemical network presented in M. L. R. van't Hoff et al. (2017) for a range of CO abundances and cosmic-ray ionization rates (ζcr). The N₂H⁺ abundance structure is then reinserted into the DALI model. Similarly, the existing CO abundance structure, and those of its isotopologues are scaled by factor xCO/10⁻⁴ and reinserted into the DALI model. Note that at this point, the model is not longer fully self-consistent, as CO abundance also affects the chemistry and more importantly the gas temperature. However, in the emitting regions of ¹³CO, C¹⁸O, C¹⁷O, and N₂H⁺, i.e., the lines of interest for this work, the gas and dust are still well coupled, meaning that the gas temperature is set by dust temperature. Finally, the excitation is recalculated, and the model is ray-traced to produce synthetic line emission.

For the input surface density, we use a tapered power law, the self-similar solution of a viscously evolving disk (e.g., D. Lynden-Bell & J. E. Pringle 1974)

$$ \Sigma_{\text{gas}}(R) = \frac{(2 - \gamma)M_{\text{gas}}}{2\pi R_c^2} \left(\frac{R}{R_c}\right)^{-\gamma} \exp\left[-\left(\frac{R}{R_c}\right)^{2-\gamma}\right]. \quad (1) $$

Here Mgas is the disk mass, Rc is the characteristic radius, and γ is the slope of the surface density, which, for a fully viscous disk coincides with the slope of the kinetic viscosity. Note that for this work, we set γ = 1 for all disks. For a typical temperature profile, this corresponds to the case where the dimensionless α viscosity parameter stays constant with radius.

The vertical density structure is given by a Gaussian with a width that increases with radius as a power law to simulate the effect that disks are flared (e.g., C. P. Dullemond & C. Dominik 2005; H. Avenhaus et al. 2018; C. J. Law et al. 2021, 2022; T. Paneque-Carreño et al. 2023):

$$ H(R) = R \ h_{100} \left( \frac{R}{100 \text{ au}} \right)^{\psi}. \quad (2) $$

Here, h₁₀₀ is the opening angle at 100 au, and ψ is the flaring angle.

Dust in the model is split into two populations (following, e.g., S. M. Andrews et al. 2011). Both grain populations follow a Mathis-Rumpl-Nordsieck size distribution with a slope of -3.5 (e.g., J. S. Mathis et al. 1977). The largest fraction (flarge) of the total dust mass is in large dust grains [1 µm–1 mm]. To simulate the fact that these grains are settled toward the midplane, their scale height is multiplied by a factor χ with
respect to that of the gas. Small grains [0.005 µm–1 µm] make up the remaining fraction (1 − f<sub>large</sub>) of the dust mass. As these grains are still well coupled to the gas, their scale height is set equal to the gas scale height.

Finally, the stellar spectrum is assumed to be a blackbody with an effective temperature $T_{\text{eff}}$ and a stellar luminosity $L_*$. Added to spectrum is a 10,000 K blackbody representing excess ultraviolet (UV) radiation released by stellar accretion. This UV radiation is of particular importance for CO as it sets its photodissociation rate. The luminosity of this component is determined by stellar mass accretion rate, under the assumption that half of the gravitational potential energy is released as radiation (e.g., M. Kama et al. 2015).

## 2.1. The AGE-PRO Model Grid

While the CO isotopologue lines predominantly trace the gas disk mass, their brightness is also affected by the physical structure of the disk as well as the characteristics of the disk-hosting star. When trying to measure the gas disk mass from $^{13}$CO and C$^{18}$O line fluxes, the disk and stellar properties effectively become nuisance parameters, reducing the accuracy with which disk mass can be obtained. For individual sources, it is possible to constrain the disk and stellar properties using additional observations, for example, by fitting the observed spectral energy distribution (e.g., E. I. Chiang & P. Goldreich 1997; P. D'Alessio et al. 1998; L. Hartmann et al. 1998; S. M. Andrews et al. 2011; M. Kama et al. 2016; K. Zhang et al. 2019, 2021; D. Deng et al. 2023). This is, in particular, the case for well-known sources, for which a large suite of ancillary observations is available. Alternatively, if these observations are not available or if the number of sources involved makes modeling them individually too time expensive, one can instead create a grid of models that covers a range of values for both the disk mass and the other parameters (e.g., J. P. Williams & W. M. J. Best 2014; A. Miotello et al. 2016, 2021; M. Ruaud et al. 2022; I. Pascucci et al. 2023). In this case, no model is expected to exactly match the observed source, but the grid as a whole can be used to estimate, within the assumed modeling framework, the uncertainty on the obtained gas mass. Given that the AGE-PRO sample includes both well-known sources, e.g., GW Lup, and sources with limited previous observations, e.g., Sz 95, we opt here for the model grid approach, leaving the modeling of individual sources for future studies (e.g., A. Sierra et al. 2024; E. TorresVillaneuva & K. Zhang 2025, in preparation). Below, we describe the modeling procedure and discuss the range of each of the parameters varied in our grid.

### 2.1.1. Modeling Procedure

The modeling begins by running DALI on a protoplanetary disk structure defined by initial parameters: gas and dust density distributions, stellar spectrum, and elemental abundances for carbon and oxygen at ISM levels (~$10^{-4}$). From this, DALI calculates fundamental disk properties such as dust temperature, gas temperature, and the intensity of ultraviolet radiation throughout the disk's radial and vertical extent. It also determines the baseline CO and $N_2$ abundance profiles, reflecting how molecules distribute and evolve under the influence of stellar irradiation, density gradients, and temperature variations.

In post-processing, the CO abundance is systematically and globally scaled down to mimic various degrees of depletion.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Flowchart of processes to generate the AGE-PRO model grid. The square boxes show actions and ellipses show outputs after each action step.</figcaption></figure>

We use the simplified chemical network of M. L. R. van't Hoff et al. (2017), paired with different cosmic-ray ionization rates, to calculate the corresponding $N_2H^+$ abundance. This updated CO and $N_2H^+$ structure is reinserted into the initial DALI model, and then we use ray-tracing function of DALI to produce simulated line fluxes for $^{12}$CO, $^{13}$CO, C$^{18}$O (2–1), $N_2H^+$ (3–2), and the 1.3 mm dust continuum, as well as the 90% flux radius of $^{12}$CO (2–1) line emission. Figure 1 shows a flowchart of the modeling processes.

### 2.1.2. Model Parameter Ranges

**Gas disk mass.** The gas disk mass ($M_{\text{disk}}$) is covered by seven values from $M_{\text{disk}} = 10^{-6} M_{\odot}$ up to $0.5 M_{\odot}$. The range of masses in the grid grew dynamically over the course of this work in order to cover all AGE-PRO sources.

**Disk size.** Similarly, the range for disk size ($R_c$), seven values between 1 and 300 au, was chosen to match the range of observed sizes within the AGE-PRO sample. We note here that
for disks with an exponential taper, the observed gas disk size, e.g., the radius that encloses 90% of the $^{12}$CO, can be up to 10 times larger than $R_c$, depending on the mass of the disk and its radial distribution (e.g., L. Trapman et al. 2023).

**Gas-to-dust mass ratio and dust properties.** For the gas-to-dust mass ratio ($\Delta_{gd}$), we assume three values (10, 100, 1000) representing a 1 dex spread around $\Delta_{gd} = 100$, which is the ratio that disks are expected to inherit from their parent cloud (e.g., P. F. Goldsmith et al. 1997). We also include two combinations of the mass fraction ($f_{large}$) and degree of settling ($\chi$) of the large grains: a “young” population ($f_{large} = 0.8, \chi = 0.6$) where large grains have started to form but have not fully settled to the midplane yet (e.g., Z.-Y. D. Lin et al. 2023; M. Villenave et al. 2023) and an “evolved” population ($f_{large} = 0.9, \chi = 0.2$) where the grains have finished settling. In the dust evolution models from N. T. Kurtovic et al. (2025), these are the approximate global, i.e., dust mass weighted and radially averaged, values of the dust population at 0.1 Myr and 1 Myr, respectively.

**Vertical structure.** The thermal structure of a disk is, to a large degree, set by how much stellar light the disk intercepts, which is determined by its vertical structure. In our grid, we cover two cases, a puffed-up disk with a relatively flared surface ($h_{100} = 0.1, \psi = 0.25$), and a much thinner disk, whose surface is only weakly flared ($h_{100} = 0.05, \psi = 0.1$). We will refer to these two cases as “flared” and “flat,” respectively. These two cases represent the range of temperature structures we can expect for protoplanetary disks (e.g., E. I. Chiang & P. Goldreich 1997; C. P. Dullemond et al. 2002; C. J. Law et al. 2021; K. Zhang et al. 2021).

**Stellar properties.** The AGE-PRO sample was selected to have a narrow range in stellar spectral types (M3-K6). Within the sample, the stellar luminosities range from $L_* = 0.08L_\odot$ to $L_* = 1.15L_\odot$. Because disk properties are known to be stellar-mass dependent (e.g., I. Pascucci et al. 2009), and the disk temperature structure depends heavily on the stellar luminosity as the main heating source, we cover this stellar luminosity range with four stellar spectra that have $L_* = [0.1, 0.25, 0.5, 1.0]L_\odot$, which we will refer to as “pale,” “faint,” “intermediate,” and “bright,” respectively. For each spectrum, we use the observations in J. Alcalá et al. (2017) to determine a typical stellar mass accretion rate for sources with similar stellar luminosities as the spectrum. Using this accretion rate, we add a far-UV (FUV) component to the spectrum, as described at the end of the previous section.

**CO abundance.** Our model does not include the chemical conversion of CO into $CO_2$ ice or the locking up of CO into larger icy bodies that have been proposed to explain the observed low CO isotopologue fluxes. Rather than including these processes, as was done in, e.g., S. Krijt et al. (2020), L. Trapman et al. (2021), D. Powell et al. (2022), M. Ruaud et al. (2022), and E. Van Clepper et al. (2022), here we instead include the peak CO gas abundance as a free parameter. This allows us to cover a much larger parameter space in a reasonable amount of time, and it limits the number of additional parameters that would need to be included and explored to a single one. In our grid, we include six values of the peak CO gas abundance between $x_{CO} = 3 \times 10^{-7}$ and $10^{-4}$, spanning the range of previous estimates (e.g., E. A. Bergin et al. 2013; C. Favre et al. 2013; L. I. Cleeves et al. 2015; L. Trapman et al. 2017; J. K. Calahan et al. 2021). One caveat of the abundance scaling method is that a simple scaling may

overpredict the abundances of rare isotopologues like C$^{18}$O due to isotopologue-selective photodissociation, when the CO gas column becomes very low (A. Miotello et al. 2014). To evaluate the overprediction levels, in Appendix A, we run additional models with initial low CO abundance and low gas disk masses to predict line fluxes and compare with that of our scaling method. The results show that in the line flux range of the AGE-PRO sample, the line fluxes between the two methods are consistent within 30%. It should be kept in mind however that this is an approximation of more complex chemical and physical processes occurring in the disk. The CO abundances in this work should therefore be seen as the bulk gas CO abundance in the $^{13}$CO and C$^{18}$O emitting layers.

**Ionization rate.** The formation of $N_2H^+$ requires $H_3^+$ that are products of cosmic-ray ionization or X-ray ionization of molecular hydrogen (T. Oka 2006). These ionization rates thus directly affect the abundance of $N_2H^+$. It is hard for X-rays to penetrate beyond the upper atmosphere of the disk. In contrast, cosmic rays penetrate deeper in the disk interior where CO freezes out and $N_2H^+$ is more likely to form. Detailed disk ionization models by D. E. Anderson et al. (2022) showed $N_2H^+$ line flux is insensitive to X-ray luminosity below $10^{30}$ erg s$^{-1}$, and only increases by a factor of 2 when X-ray luminosity increases from $10^{30}$ to $10^{31}$ erg s$^{-1}$. The cosmic-ray ionization rate is expected to have a larger effect on the abundance of $N_2H^+$. The cosmic-ray ionization rate in disks is not well constrained. It is often assumed that the rate is similar to that found in molecular clouds ($\zeta_{CR} \approx 10^{-17}$ s$^{-1}$; e.g., F. F. S. van der Tak & E. F. van Dishoeck 2000), but analysis of individual disks shows that it can be much lower ($\zeta_{CR} \approx 10^{-19}-10^{-18}$ s$^{-1}$; e.g., L. I. Cleeves et al. 2015; Y. Aikawa et al. 2021). The sensitivity of $N_2H^+$ to the cosmic-ray ionization rate translates into an uncertainty on how well it can be used to constrain the CO abundance and measure the disk mass (e.g., M. L. R. van't Hoff et al. 2017; D. E. Anderson et al. 2022; L. Trapman et al. 2022b; J. A. Sturm et al. 2023). To properly include this uncertainty, we include three cosmic-ray ionization rates in our grid, $\zeta_{CR} \in [10^{-19}, 10^{-18}, 10^{-17}]$ s$^{-1}$.

It is worth pointing out that the ionization rate can also play a role in the chemical conversion of CO into $CO_2$ ice and/or hydrocarbons. In the gas phase, CO is destroyed through reaction with $He^+$ created by cosmic rays or X-rays. On the other hand, CO ice on the grain surface can be converted into $CO_2$ and/or hydrocarbon via grain surface reactions with OH and H. The amount of OH on the grain surface depends on photodissociation rate of H$_2$O ice by cosmic rays as well as UV radiation (see A. D. Bosman et al. 2018 for a detailed discussion; see also M. Ruaud & U. Gorti 2019; M. Ruaud et al. 2022). Previous theoretical studies of CO depletion in disks suggested the chemical conversion timescale of CO is particularly sensitive to the cosmic-ray ionization rate (e.g., L. Reboussin et al. 2015; A. D. Bosman et al. 2018; C. Eistrup et al. 2018; K. R. Schwarz et al. 2018). As such, the uncertainty on the ionization rate also translates into an uncertainty on the efficiency of the chemical conversion of CO.

**Model summary.** In total, we ran 2232 models,$^{20,21}$ each with 18 different $N_2H^+$ chemistries, i.e., different ($x_{CO}, \zeta_{CR}$),

<footnote>$^{20}$ Note that not all $M_{disk}-R_c$ combinations are available. For $M_{disk} = 0.5M_\odot$ and $0.1M_\odot$, the minimum sizes are $R_c = 15$ au and 5 au, respectively. For $M_{disk} = 10^{-4}M_\odot, 10^{-5}M_\odot, \text{and } 10^{-4}M_\odot$, the maximum sizes are $R_c = 180$ au, 120 au, and 120 au, respectively.</footnote>

<footnote>$^{21}$ Model grid and fitting code are available here: https://agepro.das.uchile.cl.</footnote>
<table><caption>**Table 1**
DALI AGE-PRO Model Grid Parameters</caption><thead><tr><th>Parameter</th><th>Range</th></tr></thead><tbody><tr><td>Chemistry</td><td></td></tr><tr><td>Chemical age</td><td>1 Myr</td></tr><tr><td>x<sub>CO</sub></td><td>[0.003, 0.01, 0.03, 0.1, 0.3, 1] × 10<sup>-4</sup></td></tr><tr><td>Physical structure</td><td></td></tr><tr><td>γ</td><td>1.0</td></tr><tr><td>R<sub>c</sub></td><td>[1, 5, 15, 30, 60, 120, 180, 300] au</td></tr><tr><td>M<sub>gas</sub></td><td>[10<sup>-6</sup>, 10<sup>-5</sup>, 10<sup>-4</sup>, 10<sup>-3</sup>, 10<sup>-2</sup>, 10<sup>-1</sup>, 0.5] M<sub>☉</sub></td></tr><tr><td>Vertical structure</td><td>"flat" (h<sub>100</sub> = 0.05, ψ = 0.1)<br>"flared" (h<sub>100</sub> = 0.1, ψ = 0.25)</td></tr><tr><td>Gas-to-dust ratio</td><td>[10, 100, 1000]</td></tr><tr><td>Dust properties</td><td></td></tr><tr><td>Dust population</td><td>"young" (χ = 0.6, f<sub>large</sub> = 0.8)<br>"evolved" (χ = 0.2, f<sub>large</sub> = 0.9)</td></tr><tr><td>Composition</td><td>standard ISM<sup>a</sup></td></tr><tr><td>Stellar spectrum</td><td></td></tr><tr><td>"pale" (T<sub>eff</sub> = 3500 K, L<sub>*</sub> = 0.1L<sub>☉</sub>,<br>M<sub>*</sub> = 0.1M<sub>☉</sub>,<br>M<sub>acc</sub> = 10<sup>-11</sup>M<sub>☉</sub> yr<sup>-1</sup>)</td><td></td></tr><tr><td>"faint" (T<sub>eff</sub> = 4000 K, L<sub>*</sub> = 0.25L<sub>☉</sub>,<br>M<sub>*</sub> = 0.2M<sub>☉</sub>,<br>M<sub>acc</sub> = 10<sup>-9</sup>M<sub>☉</sub> yr<sup>-1</sup>)</td><td></td></tr><tr><td>"intermediate" (T<sub>eff</sub> = 4000 K, L<sub>*</sub> = 0.5L<sub>☉</sub>,<br>M<sub>*</sub> = 0.4M<sub>☉</sub>,<br>M<sub>acc</sub> = 10<sup>-9</sup>M<sub>☉</sub> yr<sup>-1</sup>)</td><td></td></tr><tr><td>"bright" (T<sub>eff</sub> = 4000 K, L<sub>*</sub> = 1.0L<sub>☉</sub>,<br>M<sub>*</sub> = 1.0M<sub>☉</sub>,<br>M<sub>acc</sub> = 10<sup>-8</sup>M<sub>☉</sub> yr<sup>-1</sup>)</td><td></td></tr><tr><td>ζ<sub>cr</sub><sup>b</sup></td><td>[10<sup>-19</sup>, 10<sup>-18</sup>, 10<sup>-17</sup>] s<sup>-1</sup></td></tr><tr><td>Geometry</td><td></td></tr><tr><td>i</td><td>[0, 50, 80]<sup>c</sup></td></tr><tr><td>PA</td><td>0°</td></tr><tr><td>d</td><td>150 pc</td></tr></tbody></table>

**Notes.**

post-processing steps, resulting in 40,176 distinct points in our model grid. A summary of the grid parameters is given in Table 1.

## 3. Results

### 3.1. Effect of Individual Parameters on CO and N<sub>2</sub>H<sup>+</sup> Line Fluxes

Before comparing our models to the AGE-PRO observations, we first examine the models themselves. Figure 2 shows how disk parameters other than disk mass and CO abundance affect the <sup>13</sup>CO and C<sup>18</sup>O 2–1 line luminosities. Flared disks have an overall higher disk temperature, which results in having more of the gas above the CO freeze-out temperature and thus more gas-phase CO. This effect increases the C<sup>18</sup>O luminosity more than <sup>13</sup>CO, because the more optically thin C<sup>18</sup>O is a better tracer of the full CO reservoir. For the characteristic size, instead the models show that the size has a larger effect on the more optically thick <sup>13</sup>CO, as its integrated line luminosity more strongly scales with the emitting area (\*R<sub>c</sub><sup>2</sup>). For the stellar parameters, gas-to-dust mass ratios, and dust properties, the effects are more mixed. A more luminous star raises the disk temperature, but the FUV luminosity also

increases, which photodissociates more CO. Similarly, a disk with a higher gas-to-dust ratio, which, for a fixed gas, is the same as a lower dust mass, has fewer dust particles that are able to cool the disk but also fewer small grains that absorb FUV photons. Overall, our results are consistent with the similar model grids presented in A. Miotello et al. (2016) and A. Miotello et al. (2021).

Turning now to N<sub>2</sub>H<sup>+</sup>, shown in the bottom panels of Figure 2, it immediately stands out that the range of N<sub>2</sub>H<sup>+</sup> line luminosities is somewhat larger than the range in C<sup>18</sup>O line luminosities. Within the mass bin that is shown in the figure, the N<sub>2</sub>H<sup>+</sup> line luminosity has a ~4 dex spread, whereas the C<sup>18</sup>O line luminosities have a spread of ~2 dex. The middle panel at the bottom of Figure 2 shows that a large part of this spread is driven by the compact disk models (R<sub>c</sub> ≲ 15 au). Closer inspection of these disks reveals that their CO snowline is located beyond R<sub>c</sub>, in the exponential taper of the surface density profile. Due to the steep dependence of the surface density on radius, small changes in the exact location of the CO snowline have large effects on the column density of N<sub>2</sub>H<sup>+</sup>, which is abundant between the CO and N<sub>2</sub> snowlines (e.g., C. Qi et al. 2013; M. L. R. van’t Hoff et al. 2017). For models with R<sub>c</sub> ≳ 30 au, the CO snowline lies consistently inside R<sub>c</sub>, i.e., in the power-law part of the surface density profile. Here, the exact location of the CO snowline has a much smaller effect on N<sub>2</sub>H<sup>+</sup> column density, which is reflected in the relatively small spread (≈1 dex) of N<sub>2</sub>H<sup>+</sup> line luminosities for models with R<sub>c</sub> ≳ 30 au under a single value of cosmic-ray ionization rate. Finally, the bottom-right panel of Figure 2 shows the effect of the cosmic-ray ionization rate. Increasing or decreasing ζ<sub>CR</sub> by an order of magnitude increases or decreases the N<sub>2</sub>H<sup>+</sup> flux by a factor ~3, largely independent of the other model parameters (see also J. A. Sturm et al. 2023).

### 3.2. AGE-PRO <sup>12</sup>CO, <sup>13</sup>CO, and C<sup>18</sup>O Observations versus Models

The AGE-PRO large program observed multiple CO isotopologue lines for 30 disks distributed over three star-forming regions (K. Zhang et al. 2025). In this work, we focus on the following subset: the 225 GHz continuum (100%) and C<sup>17</sup>O 2–1 (80%) integrated line fluxes for disks in Ophiuchus (D. A. Ruiz-Rodriguez et al. 2025). Here, the percentages refer to the fraction of detections (\>3σ) for that observable. For disks in Lupus, we use the 234 GHz continuum (100%), <sup>13</sup>CO 2–1 (100%), C<sup>18</sup>O 2–1 (70%), and N<sub>2</sub>H<sup>+</sup> 3–2 (40%) integrated fluxes and R<sub>CO</sub>, 90%, the radius that encloses 90% of the <sup>12</sup>CO 2–1 integrated flux (100%; D. Deng et al. 2025). For disks in Upper Sco, we use the same set of observables, i.e., the 234 GHz continuum (100%), <sup>13</sup>CO 2–1 (80%), C<sup>18</sup>O 2–1 (70%), and N<sub>2</sub>H<sup>+</sup> 3–2 (40%) integrated fluxes and R<sub>CO</sub>, 90% (100%; C. Agurto-Gangas et al. 2025).

We can compare the CO isotopologue fluxes of our models to the AGE-PRO observations. Figure 3 shows the <sup>13</sup>CO and C<sup>18</sup>O line luminosities of the 20 AGE-PRO sources in Lupus and Upper Sco on top of our models (see D. Deng et al. 2025 and C. Agurto-Gangas et al. 2025 for details on the observations). Most of the Lupus disks are consistent with 10<sup>-5</sup> ≲ M<sub>gas</sub> ≲ 10<sup>-3</sup> M<sub>☉</sub> (when assuming an ISM CO abundance), which is in line with the previous gas estimates for these disks from M. Ansdell et al. (2016) and A. Miotello et al. (2017) that were based solely on CO isotopologue line fluxes.

<footnote><sup>a</sup> J. C. Weingartner & B. Draine (2001); see also S. Facchini et al. (2017).</footnote>

<footnote><sup>b</sup> Range of cosmic-ray ionization rates only used for computing the N<sub>2</sub>H<sup>+</sup> chemistry.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Top panels: $^{13}$CO and C$^{18}$O 2–1 line luminosities ($i = 0^{\circ}$) for models with a gas mass of $M_{\text{gas}} = 10^{-3} M_{\odot}$. The same models are shown in each panel, with colors highlighting the effect of different disk parameters on the line luminosity (see Table 1). Bottom panels: the same as above, but showing C$^{18}$O 2–1 and N$_{2}$H$^{+}$ 3–2 line luminosities ($i = 0^{\circ}$) for models with a gas mass of $M_{\text{gas}} = 10^{-3} M_{\odot}$, $x_{\text{CO}} = 10^{-5}$, and $\zeta_{\text{CR}} = 10^{-18} \text{ s}^{-1}$. The bottom-right panel instead shows models with $\zeta_{\text{CR}} \in [10^{-19}, 10^{-18}, 10^{-17}] \text{ s}^{-1}$.</figcaption></figure>

However, the observations are also consistent with much more
massive disks with small $R_c$, as was previously proposed by
A. Miotello et al. (2021). Fortunately, the AGE-PRO

observations also provide some constraints on $R_c$, which can
be used to better estimate disk masses (see Appendix I). We
will discuss this further in Section 3.4.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** $^{13}$CO and C$^{18}$O J = 2–1 line luminosities for the AGE-PRO sources in Lupus and Upper Sco compared to our grid of thermochemical disk models with $x_{CO} = 10^{-4}$. The rightmost panels zoom in on the part of the main figure that contains the bulk of the AGE-PRO observations. One source, Lupus 10 (V1094 Sco) lies outside this region. A diagonal shaded region shows the absolute flux uncertainty expected for ALMA Band 6. The color of the model points show either the mass (left panel and panel (a)) or the characteristic radius (panel (b)). In the top-right panel, the marker size is scaled based on the $R_c$ of the model.</figcaption></figure>

The disks in Upper Sco cover a similar mass range as the Lupus disks, but this includes three sources (Upper Sco 2, 4, and 5) where $^{13}$CO and C$^{18}$O are either not or only tentatively detected, meaning their masses can be much lower. Lupus 10 (V1094 Sco) is the clear outlier with respect to the rest of the sample, being an order of magnitude brighter in both $^{13}$CO and C$^{18}$O compared to the other Lupus and Upper Sco disks (see S. E. van Terwisga et al. 2018; I. Pascucci et al. 2023; D. Deng et al. 2025 for a more detailed discussion on this source). It is also worth highlighting that several sources, specifically Upper Sco 4, and 5 and Lupus 9, appear to lie outside of the model grid, but since $^{13}$CO and C$^{18}$O are not detected, within the observational uncertainties, they are therefore consistent with $M_{disk} \le 10^{-5} M_{\odot}$. Furthermore, if we compare the line luminosity of $^{12}$CO 2–1, which was detected for all three sources, to our models, it suggests gas masses between $10^{-6}$ and $10^{-5} M_{\odot}$ (see Figure 18 in Appendix F).

### 3.3. Measuring CO Abundance from C$^{18}$O and $N_2H^+$

Moving our attention to $N_2H^+$ and the question of the bulk CO abundance of the AGE-PRO disks, Figure 4 shows the C$^{18}$O 2–1 and $N_2H^+$ line luminosities for models within four of the gas mass bins, with colors showing the CO abundance ($x_{CO}$) of the models. On top of these models, we show the AGE-PRO observations of the 20 AGE-PRO disks in Lupus

and Upper Sco. Focusing on the models first, we see a clear stratification based on CO abundance. Models with a lower CO abundance have fainter C$^{18}$O line luminosities but a constant or increasing $N_2H^+$ line luminosity, which is consistent with previous work (e.g., D. E. Anderson et al. 2019, 2022; L. Trapman et al. 2022a; J. A. Sturm et al. 2023). The figure also shows a not insignificant amount of overlap between models with different CO abundances, which reduces how well the bulk CO abundance of a disk can be determined.

From the 20 AGE-PRO sources in the figure, eight have detections of both C$^{18}$O and $N_2H^+$. With the exception of Lupus 3, these sources overlap with the models that have $x_{CO} = 10^{-5.5}-10^{-4.5}$ (when assuming $\zeta_{CR} = 10^{-18} \text{ s}^{-1}$), suggesting that the CO abundance in their warm molecular layer has been reduced by approximately 1 order of magnitude with respect to the ISM value. Two other disks, Lupus 3 and Upper Sco 7, the latter with no C$^{18}$O detection, are consistent with a lower CO abundance, $x_{CO} \approx 10^{-6}$. For the remaining 11 disks, $N_2H^+$ is not detected at the 3 mJy km s$^{-1}$ level, and C$^{18}$O is only detected for half of them. For these sources, it is difficult to determine how much CO, if any, has been removed from their warm molecular layer. It is also worth pointing out that the models shown in Figure 4 only represent a single cosmic-ray ionization rate, which Figure 2 showed to have a clear effect on the $N_2H^+$ flux. Including this effect, and the
<figure><img src="image_8.png" /><figcaption>**Figure 4.** $N_2H^+$ 3-2 and $C^{18}O$ 2-1 line luminosities for our models and the AGE-PRO sources. Observations for Lupus and Upper Sco are shown in black and purple, respectively. Circles denote detections in both lines, and diamonds denote $1\sigma$ for nondetections in both lines. In the case that only one of the lines is detected, the observations are shown as triangles pointing in the direction of the $1\sigma$ upper limit. Each panel shows models with a single disk mass, where the color denotes the bulk CO abundance of the model. (Note the different ranges on x- and y-axes for each panel). All models shown here were run using a cosmic ionization rate $\zeta_{\text{CR}} = 10^{-18} \text{ s}^{-1}$.</figcaption></figure>

effect of other parameters discussed previously, on measuring the gas masses and CO abundances will therefore require a different approach.

### 3.4. Measuring Gas Masses and Their Uncertainties

In the previous two sections, we have compared the observations of the AGE-PRO disks to our models, which gives us a rough estimate of their mass, bulk CO abundance, and characteristic size. However, these comparisons also revealed an inherent problem with deriving gas masses. Depending on their disk parameters, which are often unknown or poorly constrained, models with substantially different disk masses and CO abundances can produce similar $^{13}$CO, $C^{18}$O, and $N_2H^+$ fluxes (see, e.g., panel (a) in Figure 3). Here we can turn to Bayesian statistics, which provides a framework for answering the question: what is the probability of $M_{\text{gas}}$, given the observations and our prior knowledge of the disk parameters, i.e., $P(M_{\text{gas}}|F_{^{13}\text{CO}}, \dots, R_c, \dots)$. This can be obtained constructing a joint posterior probability distribution $P(M_{\text{gas}}, R_c, \dots|F_{^{13}\text{CO}}, \dots)$ and marginalizing over all disk parameters but the gas mass.

To sample the probability distribution $P(M_{\text{gas}}, R_c, \dots|F_{^{13}\text{CO}}, \dots)$ using our models, we have to make one assumption, namely that our grid is fine enough (i.e., not too coarse) that a piecewise linear interpolation accurately captures the dependence of flux on each of our disk parameters. This is a reasonable assumption for some parameters, such as $M_{\text{disk}}$ and $R_c$, where the full domain is covered in multiple steps. However, for parameters with only two values, like flaring or dust properties, it is unclear how well a linear interpolation represents the effect of these parameters on the line fluxes over the full domain included in the fit. We will discuss this further in Appendix D.

For each disk, our process to measure the gas mass is the following: We first obtain integrated $^{13}\text{CO}$ 2–1, $C^{18}\text{O}$ 2–1, $N_2\text{H}^+$ 3–2, and 1.3 mm continuum fluxes and the radius that encloses 90% of the $^{12}\text{CO}$ 2–1 line flux ($R_{\text{CO}}$, 90%) from our models for the inclination that most closely matches the inclination of the source ($i_{\text{model}} = 0^\circ$ if $i_{\text{obs}} < 30^\circ$; $i_{\text{model}} = 50^\circ$ if $30^\circ \le i_{\text{obs}} < 65^\circ$; $i_{\text{model}} = 80^\circ$ if $i_{\text{obs}} \ge 65^\circ$).

Our analysis constrains the disk gas mass by leveraging multiple lines. The $^{13}\text{CO}$ line flux provides an initial gas mass range based on comparison with our model grid. If $^{13}\text{CO}$ is optically thick, its saturation leads to a very high upper limit on the possible mass. The $C^{18}\text{O}$ flux narrows the possible mass
range. In the case of a C¹⁸O nondetection, the nondetection sets an upper limit of the gas mass, as an optically thick C¹⁸O line would be detectable. Furthermore, fitting the emitting radius of the ¹²CO constrains the compactness of the disk and thus the emitting areas of ¹³CO and C¹⁸O, thereby refining the gas mass estimate. An important additional constraint comes from the flux of N₂H⁺, which is sensitive to the CO abundance. In Appendix C, we show tests on the effect of including N₂H⁺ on the gas mass constraints. We compare the resulting gas masses based on CO isotopologues alone with masses from the combined results of CO and N₂H⁺ line fluxes. The results show that including N₂H⁺ significantly improves both the accuracy and precision of the gas mass determination.

In the case of N₂H⁺ nondetection, the minimum gas mass is set by the lowest ¹³CO+C¹⁸O gas mass corresponding to the ISM CO abundance, while the maximum is determined by the highest gas mass allowed under the maximum CO depletion consistent with the N₂H⁺ nondetection.

To calculate the probability of observed line fluxes, we model each measurement as a Gaussian distribution with the observed flux as the mean and the corresponding uncertainty as the standard deviation, regardless of whether the flux exceeds a 3σ threshold. However, for cases where the observed flux is negative, we set the value to zero. This approach is motivated by two considerations. First, for disks where we report C¹⁸O upper limits, most measured fluxes fall between the 2σ and 3σ levels, making a Gaussian treatment more representative than a hard cutoff. Second, the N₂H⁺ detections exhibit a bimodal behavior: the line is either robustly detected (>4.5σ) or detected at a low level (around 1σ) or even appears negative. We compared two methods for constraining the gas mass: (1) adopting a hard 3σ upper limit, and (2) using the full Gaussian probability distribution defined by the flux and its uncertainty. Although both methods yield similar gas mass estimates, the hard upper limit method tends to predict slightly higher masses because it assumes that all models below the upper limit are equally probable, thereby allowing for higher average N₂H⁺ fluxes. Consequently, we adopt the Gaussian distribution approach for all cases, as it better captures the statistical nature of our flux measurements.

For R<sub>CO</sub>, 90%, we also have to consider the spatial resolution of the ¹²CO observations, especially since some sources are only marginally resolved. Across the AGE-PRO sample, the ¹²CO resolution is ~0.34 with some outliers at 0.16 and 0.4. For each model in our grid, we created synthetic ¹²CO 2–1 channel maps and convolved them with circular Gaussian beams of varying size (BMAJ ∈ [0.1, 0.2, 0.3, 0.4, 0.5]). We integrate these convolved channel maps over the velocity axis to create integrated intensity (moment zero) maps, from which we measure R<sub>CO</sub>, 90% in the same manner as was done for the observations (see C. Agurto-Gangas et al. 2025; D. Deng et al. 2025 for details). For each source, we then select the convolved model radii that were created with circular Gaussian beam whose size (BMAJ<sub>model</sub>) is closest but smaller than the ¹²CO resolution of that source (BMAJ<sub>obs</sub>). To correct for the small difference between BMAJ<sub>model</sub> and BMAJ<sub>obs</sub>, we increase the model R<sub>CO</sub>, 90% by 1/(2√2ln2)(BMAJ<sub>obs</sub>−BMAJ<sub>model</sub>)². Here, we have made use of the fact that a convolution with a Gaussian with variance σ<sup>2</sup><sub>wide</sub> can be written as consecutive convolutions with two Gaussians with smaller variances σ<sup>2</sup><sub>1</sub>, σ<sup>2</sup><sub>2</sub> where σ<sup>2</sup><sub>wide</sub> = σ<sup>2</sup><sub>1</sub> + σ<sup>2</sup><sub>2</sub>.

<figure><img src="image_9.png" /><figcaption>Figure 5. Example posterior distributions of M<sub>gas</sub> and x<sub>CO</sub> for Lupus 2 (Sz 71/GW Lup), using constraints from fluxes of the ¹³CO, C¹⁸O, N₂H⁺ lines, and 1.3 mm continuum. Shown at the top of each panel and by the vertical dashed lines are the 16th, 50th, and 84th quantiles of each distribution.</figcaption></figure>

Using MCMC as implemented in the python package `emcee` (D. Foreman-Mackey et al. 2013), we sample the posterior probability distribution P(M<sub>gas</sub>, R<sub>c</sub>, ..|F<sub>13CO</sub>, ..) by minimizing

$$ \ln P(Y_{\text{model}}|X_{\text{obs}}) = -\frac{1}{2} \sum_{x \in X_{\text{obs}}} \left[ \frac{(x_{\text{model}}(Y_{\text{model}}) - x_{\text{obs}})^2}{\sigma_{x, \text{obs}}^2} \right]. \quad (3) $$

Here, $X_{\text{obs}} = \{F_{13CO}, F_{C^{18}O}, F_{N_2H^+}, F_{1.3\text{mm}}, R_{\text{CO}}, 90\%\}$ is the set of observables for each disk, and $Y_{\text{model}} = \{M_{\text{gas}}, R_c, \Delta_{\text{gd}}, \psi, f_{\text{large}}, L_*, x_{\text{CO}}, \zeta_{\text{CR}}\}$ is the set of model parameters that are varied in the grid.

Note that we sample M<sub>gas</sub>, Δ<sub>gd</sub>, x<sub>CO</sub>, and ζ<sub>CR</sub> in logspace while running the MCMC. Also note that ψ, f<sub>large</sub>, and L* are used as numerical proxies for the vertical structure, the dust properties, and the stellar properties, respectively. For example, the flux of an “intermediately flared” (ψ = 0.175) model is calculated by linearly interpolating between fluxes of a “flat”(ψ = 0.1, h<sub>100</sub> = 0.05) and “flared”(ψ = 0.25, h<sub>100</sub> = 0.1) disk model (see Table 1).

We assume uniform priors for R<sub>c</sub>, ψ, f<sub>large</sub> and loguniform priors for M<sub>gas</sub>, Δ<sub>gd</sub>, x<sub>CO</sub>, ζ<sub>CR</sub> (see Table 5 for the ranges used for each source). For the stellar luminosity, we instead assume a Gaussian prior based on the observed stellar luminosity and assuming a 30% uncertainty (see, e.g., J. Alcalá et al. 2014, 2017; J. M. Alcalá et al. 2019; C. F. Manara et al. 2020; J. M. Carpenter et al. 2025). The MCMC is computed using 1152 walkers and 4000 steps, which, after discarding the first 250 steps and thinning by samples by 40, gives a set of 107,136 samples from which we compute the posterior probability distribution. Example posterior probability distributions of M<sub>gas</sub> and x<sub>CO</sub> for Lupus 2 (Sz 71/GW Lup) are shown in Figure 5. Figure 6 shows the M<sub>gas</sub> posterior distributions for all 20 AGE-PRO disks in Lupus and Upper Sco (see Figure 19 in Appendix G for joint posterior
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Posterior probability distribution of $M_{\text{gas}}$ obtained from fitting either the observed integrated $^{13}\text{CO}$ 2–1, $^{18}\text{O}$ 2–1, $N_2\text{H}^+$ 3–2, and 1300 µm continuum fluxes (Lupus and Upper Sco) or the $^{17}\text{O}$ 2–1 and 1.300 µm continuum fluxes (Ophiuchus) using the MCMC routine outline in Section 3.4. In the case of nondetections, upper limits are used for the MCMC constraints. Orange dashed lines show the 16th, 50th, and 84th quantiles of the distribution.</figcaption></figure>

probability distributions of $M_{\text{gas}}$ and $x_{\text{CO}}$ for these sources. Full corner plots can be found in Appendix K).

As a test of the reliability of the MCMC method, we used it to measure the gas masses for the three disks where HD was

detected and compared our gas masses estimates with the independent HD-based gas masses. We find that the MCMC-estimated gas masses agree with HD-based gas masses within their respective uncertainties, similar to what was found by
L. Trapman et al. (2022b). A full description of this test can be found in Appendix B.

### 3.5. Gas Masses for the Ophiuchus Disks

The young Ophiuchus disks in the AGE-PRO sample differ from the older disks in Lupus and Upper Sco in several ways. First, these younger disks are still partially embedded in the remnants of their envelope. From the observational side, the presence of this envelope makes it harder to separate the emission originating from the disk from the contribution of the envelope (see, e.g., M. L. R. van't Hoff et al. 2018 or D. A. Ruiz-Rodriguez et al. 2025; C. Agurto-Gangas et al. (2025) for details specific to the AGE-PRO Ophiuchus sources). Furthermore, the Ophiuchus star-forming region has a larger extinction than both Lupus and Upper Sco. As a result of this, cloud emission and absorption have a significant effect on not only $^{12}$CO, but also on the $^{13}$CO and $^{18}$O emission (see D. A. Ruiz-Rodriguez et al. 2025 for details).

From the disk structure side, accretion from the envelope onto the disk, higher stellar mass accretion, and backwarming from the envelope through micron-sized grains scattering photons down to the midplane, mean that young embedded disks have a higher overall disk temperature structure than older class II protoplanetary disks (see, e.g., P. D'Alessio et al. 1998; B. A. Whitney et al. 2013; M. L. R. van't Hoff et al. 2018; M. L. R. van't Hoff et al. 2020; A. Kuznetsova et al. 2022). Due to their young age, the Ophiuchus disks are also expected to have undergone less grain growth, dust settling, and chemical conversion of CO than those in Lupus and Upper Sco, meaning that their bulk CO abundance is expected to still be close to the ISM value (e.g., K. Zhang et al. 2020).

With these points in mind, we make a few changes to the gas mass measurement method described in Section 3.4. First, we only fit the C$^{17}$O 2–1 integrated line and 1.3 millimeter flux, as C$^{17}$O is the most optically thin CO line in the AGE-PRO data and the one least affected by cloud emission/absorption (see D. A. Ruiz-Rodriguez et al. 2025). To further focus on the disk emission, we use the C$^{17}$O 2–1 measured from a tight Keplerian mask around the disk to mask out any significant emission not originating from the disk, such as the infalling material seen in Oph 7 (e.g., C. Flores et al. 2023). In some cases, absorption can also be seen in channels close to the cloud velocity. If present, these channels are masked to prevent them from lowering the C$^{17}$O 2–1 flux. We do note that even when masked out in this manner, the presence of absorption means that the integrated C$^{17}$O and the gas mass derived from it will be underestimated. For Oph 2, 6, and 9, the C$^{18}$O emission also appears reasonably clear of cloud contamination. For these sources, we repeat the fit, adding the C$^{18}$O emission and find that the resulting gas masses are in good agreement with those obtained from only fitting C$^{17}$O.

From the side of the models, we fix the dust population to “young” (see Table 1) to represent the fact that the large dust grains have not fully settled yet. Similarly, we limit $R_c \le 60$ au, as the Ophiuchus sources appear to be relatively compact, and previous studies have found young Keplerian disks to be compact (see e.g., J. R. Najita & E. A. Bergin 2018; S. Maret et al. 2020). Finally, we limit the CO abundance to a narrow range of $3 \times 10^{-5} \le x_{CO} \le 10^{-4}$. This represents the fact that, as mentioned earlier, the bulk CO abundances in the Ophiuchus are expected to be approximately the ISM value,

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Comparison of the Ophiuchus gas disk masses derived using cold ($L_* = L_{*\text{obs}}$) and hot ($L_* = 10 \times L_{*\text{obs}}$) models. A comparison of the posterior distributions is presented in Figure 20.</figcaption></figure>

but allow for the possibility that some of the processes that lower this value have already started.

To account for the overall higher disk temperature of young disks, we run two sets of gas mass fits for each source. The first, which we will refer to as the “cold model,” uses only the models with a stellar luminosity closest to the observed stellar luminosity (see Table 1 in D. A. Ruiz-Rodriguez et al. 2025). In this model, the only increase in disk temperature comes from the overall higher stellar luminosity of the Ophiuchus sources compared to the other two regions. However, as mentioned previously, the envelope will also increase the disk temperature (e.g., P. D'Alessio et al. 1998; B. A. Whitney et al. 2013; A. Kuznetsova et al. 2022; F. Han et al. 2023. The Ophiuchus sources are all Class I or Flat-Spectrum type objects and are therefore likely in the process of losing their envelope but can still have enough of an envelope present to affect the disk temperature. We therefore run a second fit, which we will refer to as the “warm model,” using models with a stellar luminosity that is ~10× the observed stellar luminosity. Given that the disk temperature is expected to scale as $T \propto L_*^{0.25}$, the “warm models” have an average disk temperature that is ~75% higher compared to the “cold models.” From the modeling side, the “cold models” can be considered an upper limit on the gas mass and the “warm models” a lower limit. We should note here however that increasing disk temperature through increasing the stellar luminosity is not strictly equivalent to including heating by an envelope, as the latter is expected to mostly affect larger radii from the star.

Figure 7 shows the effect of this disk temperature increase on the derived gas disk masses for the Ophiuchus sources. As expected, gas masses are lower if the warmer disk models are used. Based on the bottom panel of the figure, an increase of ~75% in disk temperature leads to a ~2–3× decrease in gas disk mass. It should be noted here that for Ophiuchus 1 and 7, the “cold model” gas mass posterior peaks at the upper edge of our model grid, meaning that these gas masses should be viewed as upper limits (see, e.g., 20). For their “warm model,” this is no longer the case, which means that the gas mass ratio
<figure><img src="image_12.png" /><figcaption>**Figure 8.** From top to bottom, the left panels show the gas disk masses, gas-to-dust mass ratios, and bulk CO gas abundance for sources in Ophiuchus (gray pentagons), Lupus (brown circles), and Upper Sco (green diamonds), all set against their dust disk mass. The markers show the median of the posterior distribution, and the vertical lines denote its 16th and 84th quantiles (see Section 3.4). The horizontal gray dashed lines in the top panel denote the minimum mass solar nebula (MMSN; C. Hayashi 1981) and the mass of Jupiter. The diagonal dotted lines show a constant gas-to-dust mass ratio. The right panels show histograms of the gas disk masses, gas-to-dust mass ratios and bulk CO gas abundance, obtained by summing the posterior distributions of individual sources and normalizing the histogram so that the sum over it equals 10.</figcaption></figure>

seen in Figure 7 underestimates the true gas mass ratio for these sources.

While the comparison of the “cold” and “warm” models gives a quantitative handhold on the effect of disk temperature on the derived gas mass, it is unclear if the Ophiuchus disks in the AGE-PRO sample are warmer than expected from just their stellar luminosity. For the rest of this work, we therefore adopt the “cold model” gas masses with the explicit caveat that the gas masses could be lower if the disks turn out to be warmer.

### 3.6. AGE-PRO Gas Masses

The top panels of Figure 8 show the gas masses for the AGE-PRO disks in Ophiuchus, Lupus, and Upper Sco (see also Table 2). The gas mass scales with the dust mass, taken from C. F. Manara et al. (2023). Throughout this work, we will use these dust masses, as they facilitate easy comparisons with previous work. The small exception is Upper Sco, where dust masses are available only for half of the sources, and these are based on ~335 GHz observations rather than the ~225 GHz
<table><caption>**Table 2**
Derived Gas Masses and Bulk CO Abundances</caption><thead><tr><th rowspan="3">Name</th><th rowspan="3">AGEPRO ID</th><th rowspan="3">M<sub>dust</sub></th><th rowspan="3">lg<sub>10</sub> M<sub>gas</sub></th><th rowspan="3">lg<sub>10</sub> x<sub>CO</sub></th><th colspan="7">Observational Constraints</th></tr><tr><th colspan="6">(d = detection, u = upper limit (&lt;3σ))</th></tr><tr><th>1.3 mm</th><th><sup>13</sup>CO 2–1</th><th>C<sup>18</sup>O 2–1</th><th>C<sup>17</sup>O 2–1</th><th>N<sub>2</sub>H<sup>+</sup> 3–2</th><th>R<sub>CO</sub>, 90%</th></tr></thead><tbody><tr><td>J1626-2424</td><td>Oph 1</td><td>63.6</td><td>-0.63<sup>+0.23</sup><sub>-0.33</sub></td><td></td><td>d</td><td>...</td><td>...</td><td>d</td><td>...</td><td>...</td></tr><tr><td>J1627-2420</td><td>Oph 2</td><td>16.9</td><td>-2.12<sup>+0.53</sup><sub>-0.39</sub></td><td></td><td>d</td><td>...</td><td>...</td><td>d</td><td>...</td><td>...</td></tr><tr><td>J1627-2428</td><td>Oph 3</td><td>12.5</td><td>-2.43<sup>+0.45</sup><sub>-0.26</sub></td><td></td><td>d</td><td>...</td><td>...</td><td>d</td><td>...</td><td>...</td></tr><tr><td>J1627-2427</td><td>Oph 4</td><td>1.1</td><td>-3.39<sup>+0.26</sup><sub>-0.24</sub></td><td></td><td>d</td><td>...</td><td>...</td><td>u</td><td>...</td><td>...</td></tr><tr><td>J1627-2442</td><td>Oph 5</td><td>2.6</td><td>-3.46<sup>+0.29</sup><sub>-0.20</sub></td><td></td><td>d</td><td>...</td><td>...</td><td>u</td><td>...</td><td>...</td></tr><tr><td>J1627-2440</td><td>Oph 6</td><td>21.2</td><td>-0.96<sup>+0.37</sup><sub>-0.40</sub></td><td></td><td>d</td><td>...</td><td>...</td><td>d</td><td>...</td><td>...</td></tr><tr><td>J1631-2401</td><td>Oph 7</td><td>188.2</td><td>-0.39<sup>+0.06</sup><sub>-0.10</sub></td><td></td><td>d</td><td>...</td><td>...</td><td>d</td><td>...</td><td>...</td></tr><tr><td>J1623-2302</td><td>Oph 8</td><td>2.5</td><td>-3.21<sup>+0.64</sup><sub>-0.53</sub></td><td></td><td>d</td><td>...</td><td>...</td><td>d</td><td>...</td><td>...</td></tr><tr><td>J1626-2441</td><td>Oph 9</td><td>4.7</td><td>-1.99<sup>+0.20</sup><sub>-0.26</sub></td><td></td><td>d</td><td>...</td><td>...</td><td>d</td><td>...</td><td>...</td></tr><tr><td>J1627-2439</td><td>Oph 10</td><td>16.3</td><td>-2.24<sup>+0.43</sup><sub>-0.30</sub></td><td></td><td>d</td><td>...</td><td>...</td><td>d</td><td>...</td><td>...</td></tr><tr><td>IK Lup</td><td>Lupus 1</td><td>21.2</td><td>-2.72<sup>+0.20</sup><sub>-0.17</sub></td><td>-4.57<sup>+0.16</sup><sub>-0.16</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>d</td><td>d</td></tr><tr><td>GW Lup</td><td>Lupus 2</td><td>50.1</td><td>-2.04<sup>+0.25</sup><sub>-0.26</sub></td><td>-4.86<sup>+0.14</sup><sub>-0.12</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>d</td><td>d</td></tr><tr><td>J1612-3815</td><td>Lupus 3</td><td>9.0</td><td>-2.83<sup>+0.22</sup><sub>-0.24</sub></td><td>-5.79<sup>+0.14</sup><sub>-0.24</sub></td><td>d</td><td>d</td><td>u</td><td>...</td><td>d</td><td>d</td></tr><tr><td>HM Lup</td><td>Lupus 4</td><td>4.1</td><td>-3.72<sup>+0.42</sup><sub>-0.28</sub></td><td>-5.21<sup>+0.26</sup><sub>-0.25</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>u</td><td>d</td></tr><tr><td>Sz 77</td><td>Lupus 5</td><td>1.4</td><td>-4.35<sup>+0.68</sup><sub>-0.35</sub></td><td>-4.80<sup>+0.35</sup><sub>-0.32</sub></td><td>d</td><td>d</td><td>u</td><td>...</td><td>u</td><td>d</td></tr><tr><td>J1608-3914</td><td>Lupus 6</td><td>6.4</td><td>-3.27<sup>+0.43</sup><sub>-0.29</sub></td><td>-4.76<sup>+0.29</sup><sub>-0.27</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>u</td><td>d</td></tr><tr><td>Sz 131</td><td>Lupus 7</td><td>2.6</td><td>-3.79<sup>+0.71</sup><sub>-0.38</sub></td><td>-4.22<sup>+0.15</sup><sub>-0.21</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>u</td><td>d</td></tr><tr><td>Sz 66</td><td>Lupus 8</td><td>4.7</td><td>-3.12<sup>+0.40</sup><sub>-0.33</sub></td><td>-4.93<sup>+0.23</sup><sub>-0.22</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>u</td><td>d</td></tr><tr><td>Sz 95</td><td>Lupus 9</td><td>1.4</td><td>-3.42<sup>+0.90</sup><sub>-0.71</sub></td><td>-5.62<sup>+0.51</sup><sub>-0.55</sub></td><td>d</td><td>d</td><td>u</td><td>...</td><td>u</td><td>d</td></tr><tr><td>V1094 Sco</td><td>Lupus 10</td><td>135.1</td><td>-1.19<sup>+0.11</sup><sub>-0.13</sub></td><td>-4.76<sup>+0.14</sup><sub>-0.12</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>d</td><td>d</td></tr><tr><td>J1612-3010</td><td>UppSco 1</td><td>8.5</td><td>-2.49<sup>+0.15</sup><sub>-0.11</sub></td><td>-4.82<sup>+0.11</sup><sub>-0.18</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>d</td><td>d</td></tr><tr><td>J1605-2023</td><td>UppSco 2</td><td>2.5</td><td>-3.90<sup>+0.38</sup><sub>-0.24</sub></td><td>-4.61<sup>+0.27</sup><sub>-0.27</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>u</td><td>d</td></tr><tr><td>J1602-2257</td><td>UppSco 3</td><td>1.0</td><td>-3.19<sup>+0.24</sup><sub>-0.43</sub></td><td>-4.04<sup>+0.03</sup><sub>-0.07</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>u</td><td>d</td></tr><tr><td>J1611-1918</td><td>UppSco 4</td><td>0.08</td><td>-3.69<sup>+0.92</sup><sub>-1.33</sub></td><td>-4.28<sup>+0.20</sup><sub>-0.26</sub></td><td>d</td><td>u</td><td>u</td><td>...</td><td>u</td><td>d</td></tr><tr><td>J1614-2332</td><td>UppSco 5</td><td>0.4</td><td>-5.36<sup>+0.27</sup><sub>-0.20</sub></td><td>-4.78<sup>+0.28</sup><sub>-0.36</sub></td><td>d</td><td>u</td><td>u</td><td>...</td><td>u</td><td>d</td></tr><tr><td>J1616-2521</td><td>UppSco 6</td><td>0.8</td><td>-4.16<sup>+0.30</sup><sub>-0.33</sub></td><td>-4.45<sup>+0.25</sup><sub>-0.23</sub></td><td>d</td><td>d</td><td>u</td><td>...</td><td>u</td><td>d</td></tr><tr><td>J1620-2442</td><td>UppSco 7</td><td>1.1</td><td>-3.38<sup>+0.29</sup><sub>-0.23</sub></td><td>-5.02<sup>+0.18</sup><sub>-0.20</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>d</td><td>d</td></tr><tr><td>J1622-2511</td><td>UppSco 8</td><td>7.4</td><td>-2.41<sup>+0.27</sup><sub>-0.35</sub></td><td>-4.22<sup>+0.15</sup><sub>-0.19</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>u</td><td>d</td></tr><tr><td>J1608-1930</td><td>UppSco 9</td><td>9.3</td><td>-1.28<sup>+0.20</sup><sub>-0.38</sub></td><td>-4.88<sup>+0.21</sup><sub>-0.24</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>d</td><td>d</td></tr><tr><td>J1609-1908</td><td>UppSco 10</td><td>12.9</td><td>-2.14<sup>+0.43</sup><sub>-0.45</sub></td><td>-5.52<sup>+0.22</sup><sub>-0.32</sub></td><td>d</td><td>d</td><td>d</td><td>...</td><td>d</td><td>d</td></tr></tbody></table>

Note. From left to right, the first five columns show the dust mass, obtained from converting the 1.3 mm flux under the assumptions of optically thin emission and a dust temperature of 20 K, the median gas mass, and the median bulk CO abundance, the latter two obtained from the MCMC fit (see details in Section 4.3). For Ophiuchus disks, the gas masses are derived by assuming ISM level CO abundance. Uncertainties on $M_{gas}$ (and $x_{CO}$ for disks in Lupus and Upper Sco) correspond to the 16th and 84th quartiles of their respective posterior distributions. The remaining five columns list the observational constraints used in the MCMC fit and whether it is detected at $\ge 3\sigma$ for a given source.

observations used for Ophiuchus and Lupus. To facilitate a consistent comparison, we compute the dust masses for our 10 Upper Sco sources using the 225 GHz AGE-PRO continuum fluxes from C. Agurto-Gangas et al. (2025) and the same equation and assumptions as were used in C. F. Manara et al. (2023; see also M. Ansdell et al. 2016). The dust masses are tabulated in Table 2. It should be kept in mind, however, that these dust masses are simple conversions of the millimeter flux, assuming a single temperature and dust opacity. Uncertainties in, e.g., the dust composition and grain size distribution, mean that the dust masses can only be assumed to be correct to within an order of magnitude (see, e.g., L. Ricci et al. 2010; N. P. Ballering & J. A. Eisner 2019; Á. Ribas et al. 2020; Ł. Tychoniec et al. 2020; E. Macías et al. 2021).

To measure the median gas mass in region while including the uncertainties on each individual gas mass, we use a Monte Carlo approach. Note that K. Zhang et al. (2025) calculated the

median gas mass using a Kaplar-Meier estimator instead but found similar values. For each disk in the region, we randomly draw a gas mass from its posterior distribution and compute the median for these gas masses. We repeat this process 1000 times for each star-forming region to produce a distribution of median gas masses, which are shown in Figure 9. The median gas disk mass is $7.0_{-2.6}^{+4.4} \times 10^{-3} M_{\odot}$ in Ophiuchus ($\sim 0.5$ Myr; e.g., J. I. Evans Neal et al. 2009; D. A. Ruiz-Rodriguez et al. 2025), $9.4_{-3.4}^{+5.4} \times 10^{-4} M_{\odot}$ for Lupus ($\sim 1–3$ Myr; F. Comerón 2008; P. A. B. Galli et al. 2020; D. Deng et al. 2025), and $6.8_{-2.8}^{+5.1} \times 10^{-4} M_{\odot}$ for Upper Sco ($\sim 2–6$ Myr; e.g., M. J. Pecaut et al. 2012; G. Briceño-Morales & J. Chanamé 2023; S. Ratzenböck et al. 2023; C. Agurto-Gangas et al. 2025; R. Anania et al. 2025), with approximately a 1 dex scatter in gas disk mass in all three regions. Note that the uncertainties on the medians only include the observational uncertainties on the individual gas masses and not any statistical uncertainties
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Distribution of median gas masses for the three regions, obtained using the Monte Carlo process described in Section 3.6. Circles above each distribution show its median, and the vertical line stretches from its 16th to its 84th quantiles. Open diamonds show the median values from K. Zhang et al. (2025). Note that the uncertainties shown here only include the observational uncertainties on the individual gas masses and not any statistical uncertainties due to the limited sample size of 10 sources.</figcaption></figure>

due to the limited sample size of 10 sources and systematic uncertainties due to the models (see Section 4.3).

Comparing the three regions suggests that median gas mass decreases over time, dropping by a factor of ~10 over 10 Myr. Most of this decrease happens between the ages of Ophiuchus and Lupus, as the median gas masses of Lupus and Upper Sco agree with each other within their uncertainties. In addition to that, Figure 8 shows that the gas mass distribution of disks in Lupus looks very similar to the distribution for disks in Upper Sco. At face value, this suggests there is very little evolution of the gas reservoir. However, such a direct comparison between the two regions does not include survivorship bias. During the evolution of a population, the low-mass disks could disperse first, meaning that our sample of Upper Sco disks likely represents the large massive disks that have survived up to this point (see B. Tabone et al. 2025 for a detailed discussion).

Continuing on the topic of age, Figure 10 shows the derived gas mass set against the derived stellar ages of the individual sources (see C. Agurto-Gangas et al. 2025; D. Deng et al. 2025; D. A. Ruiz-Rodriguez et al. 2025; K. Zhang et al. 2025 for details on the age derivation). Here, we see a clear evolution of the median gas mass with time from Ophiuchus (<1 Myr) to the older samples from the Lupus and Upper Sco regions. However, between Lupus and Upper Sco sources, there is only a slight decrease in the median gas mass (from 6.5 to $4.2 \times 10^{-4} M_{\odot}$) that is within the uncertainties of the median gas masses. Contributing to the overlap in the median gas masses is the large spread of individual disks, including a few relatively old (>4 Myr) Upper Sco disks with high gas masses. There is also a significant overlap in stellar age between the sources in Lupus and some of the sources in Upper. These younger Upper Sco sources tend to be more massive than the older ones. A partial age overlap can therefore help explain why the median disk masses of Lupus and Upper Sco are so similar. Considering the individual stellar ages have large

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Gas disk mass vs. stellar age for individual disks.</figcaption></figure>

uncertainties, especially some of the Upper Sco sources, constraining the gas mass evolution in disks older than 1 Myr will need a larger sample and particularly more disks older than 5 Myr.

Comparing the obtained gas masses to the dust masses, the middle panels of Figure 8 show the gas-to-dust mass ratios of the 30 disks. In Ophiuchus, we find gas-to-dust ratios from ~40 to ~1000, with a median value of $\Delta_{gd} \approx 250$. Notably, four Ophiuchus sources show very high gas-to-dust ratio, ~1000. However, we caution that dust masses for these sources may be underestimated, as all dust masses are derived under the assumption of optically thin emission and using the same dust opacity across the three regions for a homogeneous analysis. In the most-massive disks, the dust emission may be optically thick, and dust opacity in the younger disks may be different from older Class ones, as shown in N. T. Kurtovic et al. (2025), which compares dust evolution simulations with AGE-PRO continuum observations. Therefore, these high gas-to-dust ratios in the most-massive disks need to be looked at with caution.

Gas-to-dust mass ratios are lower in Lupus, ranging from ~13 to ~140 with a median value of $\Delta_{gd} = 40$. In Upper Sco, we find higher ratios again, between ~13 and 500, with a median gas-to-dust ratio $\Delta_{gd} = 116$. The increase in gas-to-dust mass ratio with age between Lupus and Upper Sco is in line with predictions from dust evolution models with no dust traps, where over time, radial drift removes dust from the disk (e.g., S. J. Weidenschilling 1977; T. Birnstiel et al. 2012; J. Dražkowska et al. 2023; T. Birnstiel 2024). It could also be indicative of the formation of planets, where dust traced by millimeter emission has been incorporated into planetesimals and planets, thereby removing it from the disk. For a more detailed analysis of dust evolution in the AGE-PRO sample, we refer the reader to N. T. Kurtovic et al. (2025).

As a whole, the oscillation of the median gas-to-dust mass ratio, i.e., a decrease from Ophiuchus to Lupus followed by a increase when moving to Upper Sco, is a somewhat puzzling. It suggests that there are different timescales for the evolution of the gas and the dust. While the dust disk mass continuously decreases across different evolutionary stages, most of the gas disk mass evolution happens between Ophiuchus and Lupus, and at a faster rate than the dust mass evolution. At later stages, there is no significant change between Lupus and
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Top-left panel: the correlation between $M_{\text{gas}}$ and $R_{\text{CO}}$, 90%, the latter being obtained from the best-fit Nuker profile of the $^{12}\text{CO}$ 2–1 emission (see L. Trapman et al. 2025 for details). Black circles and purple diamonds represent sources in Lupus and Upper Sco, respectively. Blue hexagons show the three disks with HD-based gas masses, TW Hya, DM Tau, and GM Aur, denoted by “T,” “D,” and “G.” The orange dashed line shows a power-law fit to the observations. Shown in gray are analytically calculated $R_{\text{CO}}$, 90% for different combinations of $M_{\text{gas}}$ and $R_c$ using the analytical equation for $R_{\text{CO}}$, 90% from L. Trapman et al. (2023). Here we used $L_{*} = 0.3L_{\odot}$, $\gamma = 1$, and $\xi = 0$ (see L. Trapman et al. 2023 for details). Top-right panel: $R_{\text{dust}}$, 90% vs. $M_{\text{gas}}$ for the same sources. Note that the $R_{\text{dust}}$, 90% shown here were measured from the intensity profiles fitted to the visibilities (see M. Vioque et al. 2025 for details). Bottom panel: the correlation between $M_{\text{gas}}$ and $R_c$ analytically derived from $M_{\text{gas}}$ and $R_{\text{CO}}$, 90% following the method from L. Trapman et al. (2023; see Appendix I for details.)</figcaption></figure>

Upper Sco. As a result, the gas-to-dust mass ratio changes from ISM values in Ophiuchus, to lower values in Lupus, to ISM values in Upper Sco.

The bottom panels of Figure 8 show bulk CO abundances in the $^{13}\text{CO}$ and C$^{18}\text{O}$ emitting layers. Most disks in Lupus and Upper Sco have CO abundances between 0.1× and 1× the ISM value (~$10^{-4}$; e.g., J. Lacy et al. 1994). Both regions have very similar median CO abundances, 1.4 × $10^{-5}$ and 2.0 × $10^{-5}$ for Lupus and Upper Sco, respectively. Given the scattering of the CO abundances in these two regions, the median values are statistically indistinguishable. As with the gas masses, the similarities between the CO abundances in both regions are surprising. This suggests that the removal and production of CO from the gas reaches a steady-state after ~1 Myr. For example, if the removal of CO is driven by cosmic rays (e.g., A. D. Bosman et al. 2018; K. R. Schwarz et al. 2018), the rate impinging on the disk would need to be high in the first Myr of the disk evolution, followed by an extended period of low cosmic-ray rates. Significant cosmic-ray production by protostellar jets (e.g., M. Padovani et al. 2016) during the class I stage would fit this description.

However, CO removal could also be driven by another process, such as UV processing of ices (e.g., M. Ruaud & U. Gorti 2019; K. Furuya et al. 2022) or sequestering of CO in large dust grains (e.g., S. Krijt et al. 2020). However, a full quantitative comparison between these processes and how well they describe the observations lies beyond the scope of this work.

# 4. Discussion

## 4.1. Correlations between Disk and Stellar Parameters

Having derived bulk disk properties ($M_{\text{gas}}$, $\Delta_{\text{gd}}$, $x_{\text{CO}}$), an interesting question to ask is how these correlate with other properties of the disk and of the host star. Here, we compare the gas masses derived for Lupus and Upper Sco to $T_{\text{eff}}$, $L_*$, $M_*$, $\dot{M}_{\text{acc}}$, $R_{\text{CO},90\%}$, and $R_{\text{dust},90\%}$ (J. Alcalá et al. 2014, 2017; J. M. Alcalá et al. 2019; C. F. Manara et al. 2020; J. M. Carpenter et al. 2025; see K. Zhang et al. 2025 for details). Note that here we take the $R_{\text{dust}}$, 90% from M. Vioque et al. (2025) and the $R_{\text{CO}}$, 90% from L. Trapman et al. (2025), who corrected for the convolution with the observational beam
<table><caption>**Table 3**
Correlations between Disk and Stellar Properties</caption><thead><tr><th>x</th><th>y</th><th>Spearman ρ</th><th>p-value</th><th>Correlation?</th></tr></thead><tbody><tr><td>M<sub>gas</sub></td><td>T<sub>eff</sub></td><td>0.35<sup>+0.19</sup><sub>-0.09</sub></td><td>0.07<sup>+0.11</sup><sub>-0.04</sub></td><td>positive (tentative)</td></tr><tr><td>M<sub>gas</sub></td><td>L<sub>*</sub></td><td>0.37<sup>+0.18</sup><sub>-0.09</sub></td><td>0.05<sup>+0.09</sup><sub>-0.03</sub></td><td>positive (tentative)</td></tr><tr><td>M<sub>gas</sub></td><td>M<sub>*</sub></td><td>0.38<sup>+0.16</sup><sub>-0.08</sub></td><td>0.05<sup>+0.08</sup><sub>-0.03</sub></td><td>positive (tentative)</td></tr><tr><td>M<sub>gas</sub></td><td>M&#770;<sub>acc</sub></td><td>0.16<sup>+0.22</sup><sub>-0.11</sub></td><td>0.28<sup>+0.26</sup><sub>-0.12</sub></td><td>...</td></tr><tr><td>M<sub>gas</sub></td><td>R<sub>CO</sub> 90%</td><td>0.68<sup>+0.17</sup><sub>-0.10</sub></td><td>0.0005<sup>+0.0032</sup><sub>-0.0004</sub></td><td>positive</td></tr><tr><td>M<sub>gas</sub></td><td>R<sub>dust</sub> 90%</td><td>0.49<sup>+0.15</sup><sub>-0.07</sub></td><td>0.01<sup>+0.03</sup><sub>-0.01</sub></td><td>positive</td></tr><tr><td>M<sub>gas</sub></td><td>R<sub>c</sub></td><td>0.46<sup>+0.20</sup><sub>-0.11</sub></td><td>0.02<sup>+0.06</sup><sub>-0.02</sub></td><td>positive (tentative)</td></tr><tr><td>Δ<sub>gd</sub></td><td>T<sub>eff</sub></td><td>0.32<sup>+0.24</sup><sub>-0.13</sub></td><td>0.08<sup>+0.18</sup><sub>-0.06</sub></td><td>positive (tentative)</td></tr><tr><td>Δ<sub>gd</sub></td><td>L<sub>*</sub></td><td>0.10<sup>+0.26</sup><sub>-0.13</sub></td><td>0.67<sup>+0.38</sup><sub>-0.21</sub></td><td>...</td></tr><tr><td>Δ<sub>gd</sub></td><td>M<sub>*</sub></td><td>0.38<sup>+0.23</sup><sub>-0.12</sub></td><td>0.05<sup>+0.11</sup><sub>-0.03</sub></td><td>positive (tentative)</td></tr><tr><td>Δ<sub>gd</sub></td><td>M&#770;<sub>acc</sub></td><td>-0.07<sup>+0.26</sup><sub>-0.13</sub></td><td>0.41<sup>+0.35</sup><sub>-0.17</sub></td><td>...</td></tr><tr><td>Δ<sub>gd</sub></td><td>R<sub>CO</sub> 90%</td><td>0.35<sup>+0.23</sup><sub>-0.12</sub></td><td>0.06<sup>+0.04</sup><sub>-0.04</sub></td><td>positive (tentative)</td></tr><tr><td>Δ<sub>gd</sub></td><td>R<sub>dust</sub> 90%</td><td>0.48<sup>+0.21</sup><sub>-0.11</sub></td><td>0.98<sup>+0.05</sup><sub>-0.04</sub></td><td>...</td></tr><tr><td>Δ<sub>gd</sub></td><td>R<sub>c</sub></td><td>0.22<sup>+0.28</sup><sub>-0.14</sub></td><td>0.19<sup>+0.31</sup><sub>-0.12</sub></td><td>...</td></tr><tr><td>x<sub>CO</sub></td><td>T<sub>eff</sub></td><td>-0.09<sup>+0.27</sup><sub>-0.13</sub></td><td>0.36<sup>+0.40</sup><sub>-0.18</sub></td><td>...</td></tr><tr><td>x<sub>CO</sub></td><td>L<sub>*</sub></td><td>-0.37<sup>+0.25</sup><sub>-0.12</sub></td><td>0.05<sup>+0.13</sup><sub>-0.04</sub></td><td>negative (tentative)</td></tr><tr><td>x<sub>CO</sub></td><td>M<sub>*</sub></td><td>-0.09<sup>+0.26</sup><sub>-0.13</sub></td><td>0.35<sup>+0.39</sup><sub>-0.18</sub></td><td>...</td></tr><tr><td>x<sub>CO</sub></td><td>M&#770;<sub>acc</sub></td><td>-0.36<sup>+0.28</sup><sub>-0.13</sub></td><td>0.09<sup>+0.18</sup><sub>-0.06</sub></td><td>negative (tentative)</td></tr><tr><td>x<sub>CO</sub></td><td>R<sub>CO</sub> 90%</td><td>0.03<sup>+0.25</sup><sub>-0.12</sub></td><td>0.54<sup>+0.40</sup><sub>-0.20</sub></td><td>...</td></tr><tr><td>x<sub>CO</sub></td><td>R<sub>dust</sub> 90%</td><td>-0.04<sup>+0.25</sup><sub>-0.12</sub></td><td>0.43<sup>+0.39</sup><sub>-0.18</sub></td><td>...</td></tr><tr><td>x<sub>CO</sub></td><td>R<sub>c</sub></td><td>-0.05<sup>+0.30</sup><sub>-0.15</sub></td><td>0.66<sup>+0.55</sup><sub>-0.31</sub></td><td>...</td></tr></tbody></table>

Note. Correlations with $p \le 0.01$ are assumed to be statistically significant. Those with $p - \sigma_p \le 0.05 \le p + \sigma_p$ are denoted as tentative.

by modeling the emission in the visibility and image plane, respectively.

Visual inspection reveals positive correlations between $M_{gas}$ and the two radii ($R_{CO}$, 90%, $R_{dust}$, 90%; see Figure 11), but for the other parameters, it is not clear if they correlate with $M_{gas}$ (see Figure 23 in Appendix J). To test this statistically, we carried out Spearman rank test for correlations between $x \in \{M_{gas}, \Delta_{gd}, x_{CO}\}$ and $y \in \{T_{eff}, M_{*}, L_{*}, R_{CO,90\%}, R_{dust,90\%}\}$. To include, in some cases substantial, uncertainties on $M_{gas}$, $\Delta_{gd}$ and $x_{CO}$, we used a Monte Carlo approach. Briefly, for each source, we draw a value for $x$ from their posterior distribution of the relevant parameter (see Section 3.4) and use these plus the $y$ for each source to compute the Spearman $\rho$ and the associated $p$-value. By repeating this exercise 10,000 times, we build up a distribution of Spearman $\rho$s from which we can obtain its median value and uncertainties, taken to be the 16th and 84th quantiles (see Figures 24, 26, and 28 in Appendix J).

Table 3 summarizes the results of the Spearman rank tests. In total, we find two strong ($\langle p \rangle \le 0.01$) positive correlations, between $M_{gas}$ and $R_{CO}$, 90%, and between $M_{gas}$ and $R_{dust}$, 90%. If we fit the relation between $M_{gas}$ and $R_{CO,90\%}$ with a power law, we obtain $M_{gas} \propto R_{CO,90\%}^{2.1}$. That more-massive disks have a larger observed gas disk size is to some extent expected, as the gas mass is one of the parameters that determines the observed gas disk size (see, e.g., L. Trapman et al. 2020; C. Toci et al. 2023; L. Trapman et al. 2023; F. Zagaria et al. 2023). Figure 11 shows the $M_{gas}-R_{CO}$, 90% dependence, calculated for different values of $R_c$ using the analytical formula for $R_{CO}$, 90% presented in L. Trapman et al. (2023). Based on the figure, the disks in the sample have $R_c$ between ~5 and ~60 au, with the exception of Lupus 10, which has a much larger $R_c$. In Appendix I we use the same

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Disk mass vs. stellar mass accretion rate for disks in Lupus (black) and Upper Sco (purple). The gray error bar shows typical uncertainty of the accretion rate (see J. Alcalá et al. 2017; C. F. Manara et al. 2020). Dotted gray lines show different disk lifetimes, defined as $M_{disk}/M_{acc}$.</figcaption></figure>

analytical formula to derive $R_c$ for each of the disks in our sample. Figure 11 shows that $M_{gas}$ is also positively correlated with $R_c$, suggesting that more-massive disks are also physically more extended, i.e., their mass reservoir is spread out over a larger area.

In the right panel of Figure 11, we also see a positive correlation between $M_{gas}$ and $R_{dust}$, 90%. N. T. Kurtovic et al. (2025) demonstrated that to explain the continuum observations of the AGE-PRO sample, substructures in the disks are needed. If substructures are present in the AGE-PRO sources (see also M. Vioque et al. 2025), $R_{dust}$, 90% traces the location of the farthest away dust trap. The positive correlation between $M_{gas}$ and $R_{dust}$, 90% is therefore indirect evidence that more-massive disks can develop substructures farther away from the star.

Several other parameter combinations show tentative positive or negative correlations (see Table 3). For example, there is some evidence for a positive correlation between $M_{gas}$ and stellar mass, but its median $p$-value is 0.06, meaning there is a reasonable chance that the correlation is not statistically significant. This is not unsurprising, as by design, the AGE-PRO sample has a narrow stellar mass and spectral type range and lacks the baseline with which previous correlations were found (e.g., I. Pascucci et al. 2016; L. Testi et al. 2022). The forthcoming ALMA cycle 9 large program DECO (2022.1.00875.L, PI: Ilse Cleeves), which covers a similar spectral set-up as AGE-PRO, will significantly expand the stellar mass and spectral type baselines, meaning a combined analysis of these two programs will significantly increase our understanding of the dependence of disk properties such as gas disk mass on the properties of their host star.

## 4.2. Estimating Disk Lifetimes

One of key parameters for planet formation is the lifetime of the protoplanetary disk. A rough estimate of the disk lifetime can be obtained from the ratio of the disk mass and the stellar accretion rate, $M_{gas}/M_{acc}$ (see, e.g., M. G. Jones et al. 2012; G. Lodato et al. 2017; G. P. Rosotti et al. 2017). Of the 30 disks in the AGE-PRO sample, all 10 disks in Lupus have measured accretion rates from near-UV/optical spectroscopy (J. Alcalá et al. 2014, 2017; J. M. Alcalá et al. 2019) and six of
the 10 disks in Upper Sco have reported values as well (C. F. Manara et al. 2020; M. Fang et al. 2023).

Figure 12 sets $M_{\text{gas}}$ from this work against the observed stellar mass accretion rate. Six of the 10 disks in Lupus occupy the part of the $\dot{M}_{\text{acc}}-M_{\text{disk}}$ plane with disk lifetimes of 1–10 Myr, which is consistent with the age of the star-forming region. Two sources, Lupus 7 and 8, lie between 0.1 and 1.0 Myr, and the remaining two sources, Lupus 4 and Lupus 5, have estimated disk lifetimes $\lesssim 0.1$ Myr, much shorter than the age of the region. Especially for the latter two sources, independent of the mechanism driving their evolution, these disks should have already been dispersed based on their current gas mass and stellar accretion rate, assuming that the accretion rate has been constant during the disk’s evolution. For Upper Sco, five of the six disks have disk lifetimes consistent with the age of Upper Sco, but one disk, Upper Sco 2, has a disk lifetime of $\lesssim 1$ Myr and should have already been dispersed.

This raises the question, how are the three disks with lifetimes shorter than $\sim 10 \times$ their star-forming region age still around? We will discuss three possible sets of explanations: (1) the gas masses are underestimated, (2) the stellar accretion rates are overestimated, and (3) the lifetime computed from the disk mass and stellar accretion disk mass does not accurately represent the true lifetime of the disk. S. L. Grant et al. (2023) recently showed that a similar disk lifetime problem exists for disks of Herbig Ae/Be stars, which they suggest could partially be explained by a bias in the classification of Herbig stars. A similar bias is less likely for T-Tauri stars, but it could be a contributing factor.

The first possibility is that the gas masses for these sources are underestimated. We have dedicated Section 4.3 to discussing different scenarios by which the gas masses presented in this work could be underestimated. Looking at Figure 12, requiring that the disk lifetime be equal to or longer than the star-forming region age would require increasing the gas mass by a factor $\sim 10–50$. It is interesting to note here that all of these disks are also compact ($R_{\text{CO}}$, 90% $\lesssim 100$ au; C. Agurto-Gangas et al. 2025; D. Deng et al. 2025). Whether this indicates a different evolutionary path or simply highlights a difficulty with measuring gas masses and/or stellar mass accretion rates for compact disks remains unclear.

Alternatively, the stellar mass accretion rates could be overestimated. The accretion rate is based on the accretion luminosity ($L_{\text{acc}}$), which was obtained by fitting the observed UV spectrum of each source with a template spectrum of a nonaccreting class III star plus a slab of hydrogen (see L. Hartmann et al. 2016; C. F. Manara et al. 2016; J. Alcalá et al. 2017 for details on how accretion columns are modeled). The mass accretion is then computed using $\dot{M}_{\text{acc}} \approx 1.25 L_{\text{acc}} R_{*}/GM_{*}$, where $R_{*}$ and $M_{*}$ are the stellar radius and mass, respectively (see, e.g., J. Alcalá et al. 2017; C. F. Manara et al. 2020).

The main sources of uncertainty for $\dot{M}_{\text{acc}}$ are the extinction, the choice of stellar template, and the uncertainties on the stellar parameters, the latter of which depend on the assumed stellar evolutionary track. In aggregate, J. Alcalá et al. (2014) estimated that this leads to a $\sim 0.35$ dex, or a factor of $\sim 2.3$, total uncertainty on $\dot{M}_{\text{acc}}$ (see also J. Alcalá et al. 2017; C. F. Manara et al. 2023). Also, we note that C. V. Pittman et al. (2022) used near-UV data to show that accretion rates derived with the before discussed methodology can be underestimated by a factor of 3–4, which would lower the

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Lupus and Upper Sco gas masses derived using either $\zeta_{\text{CR}} \in [10^{-19}, 10^{-17}]$ s$^{-1}$ (clear markers) or $\zeta_{\text{CR}} \in [10^{-20}, 10^{-19}]$ s$^{-1}$ (solid markers) as a prior for the cosmic-ray ionization rate (see Section 3.4 for details). Note that in most cases, decreasing $\zeta_{\text{CR}}$ increases $M_{\text{gas}}$.</figcaption></figure>

disk lifetime even more. This suggests that the uncertainty on the stellar mass accretion rate alone cannot explain the short disk lifetimes for five of our sources.

Finally, the fault could be in our interpretation of $M_{\text{disk}}/\dot{M}_{\text{acc}}$ as the disk lifetime. Stellar accretion is known to be a stochastic rather than a smooth process (see W. J. Fischer et al. 2023 for a recent review), meaning the currently observed stellar accretion rate need not be the same as the disk-lifetime integrated accretion rate. In addition, accretion bursts of a factor $\sim 10$, the approximate factor the lifetime-averaged $\dot{M}_{\text{acc}}$ would need to be lower than the current value to provide a disk lifetime of $\sim 1$ Myr, are within the expected range of variability (e.g., W. J. Fischer et al. 2023). However, we do not expect accretion bursts to be the general explanation, as they are uncommon, and our samples were not selected to favor strong accretors.

### 4.3. Could the Gas Masses Be Underestimated?

As discussed previously in Section 4.2, disk lifetimes (based on estimated disk masses and accretion rates) for some disks are lower than expected and apparently inconsistent with the age of the star-forming region in which these disks reside. It is therefore worthwhile to discuss whether the disks masses measured in this work are underestimated. In this section, we will examine a few possible reasons.

Low cosmic ionization rate ($\zeta_{\text{CR}}$). as an ion, the abundance of $N_2H^+$ is directly affected by the amount of ionization in the disk. For $N_2H^+$, which is abundant close to the disk midplane, the ionization is dominated by cosmic rays. This creates a link between $\zeta_{\text{CR}}$ and $M_{\text{gas}}$: a disk with a lower $\zeta_{\text{CR}}$ requires a larger $M_{\text{gas}}$ to reproduce the observed $N_2H^+$ emission. Or in the case of an upper limit, a lower $\zeta_{\text{CR}}$ would mean that a higher gas mass is still consistent with the observed upper limit. In our analysis, we have assumed a wide range of cosmic-ray ionization rates, $\zeta_{\text{CR}} \in [10^{-19}, 10^{-17}]$ s$^{-1}$. However, work by L. I. Cleeves et al. (2015) on the disk around TW Hya suggests that the cosmic-ray ionization rate in disks could be even lower, $\zeta_{\text{CR}} \sim 10^{-20}$ s$^{-1}$. If this is also the case, in some or all of the AGE-PRO disks, our gas masses could be underestimated.

In order to quantify this effect, we compute $N_2H^+$ fluxes for all models assuming $\zeta_{\text{CR}} = 10^{-20}$ s$^{-1}$ and derive the gas masses
using a low $\zeta_{CR}$ prior, $\zeta_{CR} \in [10^{-20}, 10^{-19}]$ s$^{-1}$, consistent with the range found by L. I. Cleeves et al. (2015). These new gas masses are compared to our canonical gas masses in Figure 13. In all cases, the low-$\zeta_{CR}$-gas mass is indeed larger the canonical gas mass, usually by a factor $\sim$2–4. This is in particular the case for disks where $N_2H^+$ was detected, as a detection provides a harder constraint while fitting the observations. As discussed in Section 3.4, when $N_2H^+$ is not detected, the gas mass upper limit is set by the maximum CO depletion that is consistent with the $N_2H^+$ nondetection. With a lower cosmic-ray ionization rate, the overall $N_2H^+$ abundance is lower, meaning that more-massive disks with a lower CO abundance still produce an $N_2H^+$ flux that is consistent with the observed upper limit. To summarize, assuming a lower cosmic-ray ionization rate does increase $M_{gas}$, but it cannot fully explain the short disk lifetimes that were found in Section 4.2.

**Low $N_2$ abundance.** One of the underlying assumptions in our $N_2H^+$ chemical network is that $N_2$ is the main carrier of nitrogen in protoplanetary disks. If this is not the case or if $N_2$ is underabundant in a manner similar to CO, then our models would be overestimating the emission of $N_2H^+$, meaning that the gas masses derived here are underestimated. So far, studies seem to suggest that disks are not underabundant in nitrogen (e.g., L. I. Cleeves et al. 2018) and that other species like ammonia are not a dominant nitrogen carrier (e.g., V. N. Salinas et al. 2016). These results are consistent with recent modeling results, which show that degree of nitrogen depletion from the warm molecular layer is much lower than that of carbon and oxygen (K. Furuya et al. 2022). However, the sample of sources with a well-studied nitrogen reservoir is still small, and there is some evidence from solar system bodies that some fraction of the available nitrogen could be locked up in ammonia salts (e.g., A. Boogert et al. 2015; K. Altwegg et al. 2020).

Lowering the $N_2$ abundance is similar to lowering $\zeta_{CR}$, in that both reduce the amount of $N_2H^+$, which in turn means that a larger gas mass is needed to match the observations. Therefore, if the $N_2$ abundance is indeed lower than what is assumed in our models, we would expect that the effect on the gas mass is similar to the low $\zeta_{CR}$ experiment discussed above. For disks with $N_2H^+$ detections, the gas mass will increase by a factor inversely proportional to the amount of $N_2$ removed, but for disks with $N_2H^+$ upper limits, the effects are likely minimal. Since the disks with the lowest disk lifetimes also do not have detected $N_2H^+$ emission, neither lowering the cosmic-ray ionization rate nor reducing the bulk $N_2$ abundance is expected to solve the discrepancy between the derived disk lifetimes and the age of the star formation region.

**Model assumptions.** The disk models used in this work to derive gas masses from the observations depend on a large number of parameters and assumptions. A number of these, for example, the disk size, were varied within the model grid, and their effect was therefore included in the gas mass measurement, but many others were not. For example, all models have an exponentially tapered gas surface density, whereas we know at least Upper Sco 1 to have a gas gap (A. Sierra et al. 2024). This gas surface density profile is expected for turbulently evolving disks (e.g., D. Lynden-Bell & J. E. Principle 1974; L. Hartmann 1998; B. Tabone et al. 2022), but observational constraints on the gas surface density in the outer disk are limited (e.g., C. P. Dullemond et al. 2020; K. Zhang et al. 2021). Moreover, this theoretically predicted gas surface density profile neglects the effects of external photoevaporation, an assumption that likely cannot be made

for the disks in Upper Sco (see R. Anania et al. 2025 for a detailed discussion on the effect on external photoevaporation on the AGE-PRO Upper Sco disks.)

As another example, the models use a Gaussian vertical density distribution expected for vertically isothermal disks, but it is well known that disks are not vertically isothermal. M. Ruaud et al. (2022) recently showed that this assumption can have a significant effect on the CO isotopologue line fluxes, although, as A. D. Bosman et al. (2022) argued, their analysis is also not free of its own underlying assumptions. With these points in mind, it is therefore interesting to note that for the Lupus disks in AGE-PRO, D. Deng et al. (2025) showed that gas masses obtained using the M. Ruaud et al. (2022) models are within a factor $\sim$3–5 from the values reported in this work (see D. Deng et al. 2025 for details).

Added up, these and the other unmentioned model assumptions will have an effect on the measured gas masses. While it is hard to determine how much they would affect the gas mass estimates or in what direction they would change them, assuming an additional uncertainty of a factor $\sim$2–3 is likely advisable based on previous benchmarking of similar codes (e.g., M. Röllig et al. 2007).

**Compact disks:** A unifying characteristic of the disks with short disk lifetimes is that all of them are compact ($R_{CO}$, 90% $\lesssim$ 100 au, a measurement that still includes convolution with a beam that is $\sim$56 au in diameter). A. Miotello et al. (2021) proposed that compact, massive disks provided an alternative explanation for the faint $^{12}$CO and $^{13}$CO line emission observed in Lupus (see M. Ansdell et al. 2016, 2018). Fitting the observations with our model grid, which includes disks with $R_c = 1$ au and $M_{gas}$ up to 0.1 $M_\odot$, shows that this cannot be the sole explanation. Examining the $R_c$-posteriors of the individual disks reveals that the fit does favor compact disks ($R_c \approx$ 1–25 au; see Appendix K), but compact massive disks can be ruled out either by the need to reproduce $R_{CO}$, 90% (which sets a minimum on $R_c$) and/or the requirement for $N_2H^+$ to remain undetected (which sets a maximum on the gas mass). The latter can also be seen in the fact that the $\zeta_{CR}$-posterior highly favors a low cosmic-ray ionization rate for these disks.

That being said, however, these compact disks have also been considerably less well studied than their larger, brighter counterparts (e.g., N. T. Kurtovic et al. 2021; J. Pegues et al. 2021). Most of our current understanding about protoplanetary disks is therefore based on the properties of large disks and the models that try to reproduce those. For example, the $N_2H^+$ + $C^{18}O$ gas mass measurement was tested against the three disks with HD detections, all of which are large disks (L. Trapman et al. 2022b). If compact disks are not simply scaled-down large disks and instead have different physical structure and/or chemical composition, their gas masses estimates would be inaccurate and could be underestimated. This suggests the need for future continuum and gas observations focused on more compact disks ($R_c \lesssim 15$ au). Such compact disks are much more common than previously thought and crucial to our understanding of planet formation in nearby star-forming regions.

## 5. Conclusions

In this work, we examined line observations of $^{12}$CO, $^{13}$CO, $C^{18}O$, $C^{17}O$ 2–1, and $N_2H^+$ 3–2 of the AGE-PRO large program and compared them to a large grid of protoplanetary disk models. We combine a piecewise linear interpolation of
the model ¹³CO, C¹⁸O, and N₂H⁺ line fluxes and ¹²CO gas
disk size with an MCMC routine to derive gas masses and their
associated uncertainties for the 20 AGE-PRO disks in the
Lupus and Upper Sco star-forming regions. For the 10 AGE-
PRO disks in the Ophiuchus star-forming region, we compare
the observed C¹⁷O 2–1 line fluxes to a wide range of models to
estimate their gas masses. Our conclusions are as follows:

1. Using our MCMC routine, we measure the gas masses of class II disks with a typical precision of ~0.35 dex. Uncertainties for the younger class I sources are much larger due to uncertainties in their temperature structure and difficulties due to cloud contamination and envelope contribution to the line fluxes. The accuracy of the gas masses is harder to quantify, but analysis in D. Deng et al. (2025) suggests an overall good agreement, i.e., within a factor of ~3–5, between gas masses derived using different models and methods.

2. The median gas mass of the three regions decreases over time, from $7.0_{-2.6}^{+4.4} \times 10^{-3} M_{\odot}$ in Ophiuchus ($\lesssim 1$ Myr) to $9.4_{-3.4}^{+5.4} \times 10^{-4} M_{\odot}$ for Lupus ($\sim 1-3$ Myr) and $6.8_{-2.8}^{+5.1} \times 10^{-4} M_{\odot}$ for Upper Sco ($\sim 2-6$ Myr), with $\sim 1$ dex scatter in gas mass in each region. Almost all of the gas mass evolution occurs between the ages of Ophiuchus and Lupus, given that the median gas masses of Lupus and Upper Sco agree within their uncertainties and the gas mass distributions of Lupus and Upper Sco appear very similar. This suggests limited evolution of the gas mass distribution at later times, but it might be due to survivorship bias, i.e., disk dispersal selects for the most-massive disks.

3. The median gas-to-dust mass ratio of disks in Upper Sco ($\langle\Delta_{gd}\rangle \approx 120$) is a factor of $\sim 3$ higher than that of disks in Lupus ($\langle\Delta_{gd}\rangle \approx 40$). Given that the gas masses of Lupus and Upper Sco are similar, this evolution in the gas-to-dust mass ratio is driven by the evolution of the dust, likely efficient inward drift of pebbles and/or the formation of planetesimals and planets.

4. Median bulk CO abundances in the warm molecular layer are about an order of magnitude below the ISM value for both disks in Lupus and Upper Sco. This suggests that CO is efficiently removed from the warm molecular layer in the first $\sim 1$ Myr, followed by a steady state that lasts until the end of the gas disk lifetime.

5. We find strong correlations between $M_{\text{gas}}$ and the gas disk size $R_{\text{CO}}$, 90% and between $M_{\text{gas}}$ and the dust disk size $R_{\text{dust}}$, 90%. The former is in some ways expected, as $M_{\text{gas}}$ plays an important role in setting $R_{\text{CO}}$, 90%, but further analysis shows that there is also a correlation between $M_{\text{gas}}$ and the characteristic radius $R_c$, suggesting that more-massive disks are physically more extended than lower-mass disks. The correlation between $M_{\text{gas}}$ and $R_{\text{dust}}$, 90% could be indirect evidence that more-massive disks can develop substructures farther away from the star.

6. Most AGE-PRO disks have disk lifetimes, computed as $M_{\text{gas}}/M_{\text{acc}}$, that are similar or longer than the age of the star-forming region in which they reside. However, a number of disks with measured $\dot{M}_{\text{acc}}$ have a disk lifetime that is $\gtrsim 10 \times$ shorter than their star-forming region age. This suggests that they have not maintained their current accretion rate throughout their lifetime, that their gas mass is underestimated, or that they will disperse in the near future.

Acknowledgments

This paper makes use of the following ALMA data: ADS/
JAO.ALMA#2021.1.00128.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 Asso-
ciated Universities, Inc.

We thank the anonymous reviewer for providing helpful
comments. L.T. and K.Z. acknowledge the support of the NSF
AAG grant No. 2205617. G.R. acknowledges funding from the
Fondazione Cariplo, grant No. 2022-1217, and the European
Research Council (ERC) under the European Union's Horizon
Europe Research & Innovation Programme under grant
agreement No. 101039651 (DiscEvol). Views and opinions
expressed are however those of the author(s) only, and do not
necessarily reflect those of the European Union or the European
Research Council Executive Agency. Neither the European
Union nor the granting authority can be held responsible for
them. P.P. and A.S. acknowledge the support from the UK
Research and Innovation (UKRI) under the UK government's
Horizon Europe funding guarantee from ERC (under grant
agreement No. 101076489). A.S. also acknowledges support
from FONDECYT de Postdoctorado 2022 No. 3220495. B.T.
acknowledges support from the Programme National "Physique
et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU
with INC/INP and co-funded by CNES. I.P. and D.D.
acknowledge support from Collaborative NSF Astronomy &
Astrophysics Research grant (ID: 2205870). C.A.G. and L.P.
acknowledge support from FONDECYT de Postdoctorado 2021
No. 3210520. L.P. also acknowledges support from ANID
BASAL project FB210003. L.A.C and C.G.R. acknowledge
support from the Millennium Nucleus on Young Exoplanets and
their Moons (YEMS), ANID-Center Code NCN2024_001, and
L.A.C. also acknowledges support from the FONDECYT grant
No. 1241056. N.T.K. acknowledges support provided by the
Alexander von Humboldt Foundation in the framework of the
Sofja Kovalevskaja Award endowed by the Federal Ministry of
Education and Research. K.S. acknowledges support from the
European Research Council under the Horizon 2020 Framework
Program via the ERC Advanced Grant Origins 83 24 28. The
National Radio Astronomy Observatory and Green Bank
Observatory are facilities of the U.S. National Science
Foundation operated under cooperative agreement by Asso-
ciated Universities, Inc. All figures were generated with the
PYTHON-based package MATPLOTLIB (J. D. Hunter 2007).
This research made use of Astropy,<fnref n="22" /> a community-developed
core Python package for Astronomy (Astropy Collaboration
et al. 2013, 2018), and Scipy (P. Virtanen et al. 2020).

**Appendix A**

Line Flux Comparisons

In Section 2.1, we first compute gas temperature and CO
isotopologue abundances from an ISM level of C/H and O/H
elemental ratio. We then mimic the CO depletion by scaling
down the peak CO abundance and calculating line fluxes. This

<footnote n="22">²² http://www.astropy.org</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Top row: line fluxes of models with low C/H and O/H elemental ratios. The gray dashed line indicates the 3σ detection level of AGE-PRO observations. Bottom row: line flux ratios between the models with initial low elemental ratios with models with scaled CO abundance structures.</figcaption></figure>

cost-effective approach allows for efficient exploration of a broad range of parameter space. Although isotopologue-selective photodissociation is included in the initial models, when the CO gas column becomes very low (A. Miotello et al. 2014), simply scaling the abundances may overpredict the abundances of rare isotopologues like C¹⁸O, because these rare isotopologues experience less self-shielding than the more abundant ones. To assess this overestimation, we run additional models with initial low CO abundance and low gas disk masses and compare the line fluxes with those from the scaling method. Specifically, we adopt C/H ratios of 1.4 × 10⁻⁵ and 1.4 × 10⁻⁶, representing CO depletion by a factor of 10 and 100 times relative to the ISM levels. The O/H ratio in these models is also scaled similarly. We test disk gas masses of 10⁻⁴–10⁻³ M☉, characteristic radii R_c of 5–60 au, an L_★ = 0.5 L☉, a gas-to-dust mass ratio of 100, a flat disk structure, and an evolved dust population, which spans the bulk parameter range of the AGE-PRO Lupus and Upper Sco sample. These models use a cosmic-ray ionization rate of 10⁻¹⁷ s⁻¹ and run for 1 Myr.

Figure 14 shows the line fluxes from the additional models, and the line flux ratios of ¹³CO and C¹⁸O (2–1), calculated by dividing fluxes from the initial low abundance models by those from the scaling models. For the $x_{CO} = 10^{-5}$ case, both ¹³CO and C¹⁸O flux ratios are consistent with 30% for gas disk masses of 10⁻⁴ and 10⁻³ M☉. For the extremely depleted

$x_{CO} = 10^{-6}$, the ¹³CO flux ratios stay within 30%, except for the most compact disks ($R_c = 5$ au). For C¹⁸O in $x_{CO} = 10^{-6}$ models, the line ratios vary significantly, ranging between 1 and 4. However, these cases have C¹⁸O line fluxes much lower than the 3σ flux detection limits of our AGE-PRO measurements. In these cases, their gas disk masses are mainly constrained by their ¹³CO flux. Therefore, for the line flux ranges probed by the AGE-PRO sample, we consider the line flux predictions from abundance scaling method valid.

## Appendix B

### Disk Mass Retrieval Comparison with HD-based Gas Mass Measurements

In this work, we demonstrated an approach for deriving gas disk masses using an MCMC fit of the ¹³CO 2–1, C¹⁸O 2–1, N₂H⁺ 3–2 line fluxes and the 1.3 mm continuum flux. Since this is the first application of this approach, it is highly worth testing it against an independent method of measuring gas masses, such as using hydrogen deuteride. The HD 1-0 line at 112 µm has been detected in three protoplanetary disks with Herschel and has proven to be a robust measurement of the gas mass (e.g., E. A. Bergin et al. 2013; M. K. McClure et al. 2016; K. R. Schwarz et al. 2016; L. Trapman et al. 2017; J. K. Calahan et al. 2021; K. R. Schwarz et al. 2021). A comparison between HD-based gas masses and gas masses derived from the combination of C¹⁸O and N₂H⁺ was recently presented
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Posterior distribution of $M_{gas}$, $R_c$, $x_{CO}$, $\zeta_{cr}$, and $L_*$ for TW Hya (top), DM Tau (middle), and GM Aur (bottom). Black vertical dashed lines show the 16th, 50th, and 84th quantiles of the distribution. Shown in orange are the HD-based gas mass and the corresponding $x_{CO}$ from L. Trapman et al. (2022b), with the dashed and dotted lines showing the measured value and its uncertainties, respectively. The orange dashed lines in the $R_c$ and $L_*$ panels show the values used in the source-specific models used to derive the HD-based gas mass (see M. Kama et al. 2016; K. Zhang et al. 2019, 2021; K. R. Schwarz et al. 2021).</figcaption></figure>

<table><caption>**Table 4**
Source Information</caption><thead><tr><th>Source</th><th>L<sub>*</sub><br>(L<sub>☉</sub>)</th><th>dist<br>(pc)</th><th>1.3 mm<br>(mJy)</th><th><sup>13</sup>CO 2–1<br>(Jy km s<sup>-1</sup>)</th><th>C<sup>18</sup>O 2–1<br>(Jy km s<sup>-1</sup>)</th><th>N<sub>2</sub>H<sup>+</sup> 3–2<br>(Jy km s<sup>-1</sup>)</th><th>HD-based M<sub>gas</sub><br>(M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>TW Hya</td><td>0.28</td><td>59</td><td>...</td><td>1.89 ± 0.18</td><td>0.57 ± 0.06</td><td>2.0 ± 0.21</td><td>1.15 × 10<sup>-2</sup> (7.5 × 10<sup>-3</sup>, 1.9 × 10<sup>-2</sup>)</td></tr><tr><td>DM Tau</td><td>0.24</td><td>145</td><td>...</td><td>4.2 ± 0.076</td><td>1.11 ± 0.11</td><td>1.286 ± 0.17</td><td>1.3 × 10<sup>-1</sup> (4.5 × 10<sup>-2</sup>, 2 × 10<sup>-1</sup>)</td></tr><tr><td>GM Aur</td><td>1.2</td><td>159</td><td>...</td><td>5.028 ± 0.048</td><td>1.092 ± 0.039</td><td>1.487 ± 0.18</td><td>1.15 × 10<sup>-1</sup> (4 × 10<sup>-2</sup>, 2.15 × 10<sup>-1</sup>)</td></tr><tr><td>References</td><td>(1, 2, 3)</td><td></td><td></td><td>(4, 5)</td><td>(6, 7, 5)</td><td>(8, 9)</td><td>(10)</td></tr></tbody></table>

**Note.** ALMA line fluxes include a 10% systematic flux uncertainty.

**References.** (1) S. M. Andrews et al. (2011), (2) S. J. Kenyon & L. Hartmann (1987), (3) Gaia Collaboration et al. (2018), (4) C. Favre et al. (2013), (5) K. I. Öberg et al. (2021), (6) J. K. Calahan et al. (2021), (7) J. B. Bergner et al. (2019), (8) C. Qi et al. 2025, in preparation, (9) C. Qi et al. (2019), (10) L. Trapman et al. (2022b).

by L. Trapman et al. (2022b). The test presented here is similar, but instead of using the source-specific models used in the aforementioned work, we instead derive gas masses using the MCMC method outlined in this work and compare them to the HD-based gas masses from L. Trapman et al. (2022b).

We ran the MCMC method as outlined in Section 3.4 for TW Hya, DM Tau, and GM Aur using the integrated fluxes collated in Table 4. For GM Aur, a small modification was made. Its stellar luminosity lies outside of our grid of models (the “bright” star had $L_* = 1.0L_☉$; see Table 1), so we limit the

models to only “bright” stars. Figure 15 shows the resulting posterior distributions for $M_{gas}$, $R_c$, $x_{CO}$, $\zeta_{cr}$, and $L_*$. For DM Tau and GM Aur, we find good agreement between the two gas mass estimates, with the median gas mass obtained from the MCMC lying with the uncertainties on the HD-based gas mass. For TW Hya, the gas mass posterior distribution is bimodal, with a lower-mass peak (around $\sim 10^{-2}M_☉$) coinciding with the HD-based gas mass and a high-mass peak that lies outside of our prior ($M_{gas} > 0.1M_☉$). Closer inspection reveals that this bimodality is correlated with the cosmic-ray
ionization rate. The lower gas mass corresponds to $\zeta \approx 10^{-17} s^{-1}$, whereas the higher gas mass corresponds to a lower rate of $\zeta_{CR} \lesssim 10^{-18} s^{-1}$.

These comparisons show that the grid model approach recovers mass well but does not sufficiently constrain the $R_c$. The previous source-specific models had various assumptions, and the $R_c$ estimations do not have any uncertainties, which suggest that the $R_c$ solution may not be unique in these models. Therefore, surface mass distribution needs further investigation. E. Torres Villanueva et al. (2025, in preparation) sourced specific models to match the observed radial profiles of $^{13}$CO and C$^{18}$O line emission, which will provide more insights on the mass distribution.

## Appendix C

### Testing Gas Mass Constraints with and without $N_2H^+$

To test the effect of including $N_2H^+$ on gas mass constraints, we compare MCMC results of gas masses derived from CO isotopologues alone with those derived from both CO isotopologues and $N_2H^+$ line fluxes. Figure 16 shows the resulting gas mass constraints without (left) and with (right) $N_2H^+$ for the Lupus 3 source as an example. After marginalizing over all other disk properties—including the uncertainty of cosmic-ray ionization rate—we find that including $N_2H^+$ significantly increases both the accuracy and precision of the gas mass estimate.

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Comparison of the accuracy of gas mass constraints without (left) and with (right) $N_2H^+$ for the Lupus 3 source.</figcaption></figure>
## Appendix D

### Testing the (Piecewise) Linear Interpolation Assumption

One of the main assumptions we have to make in order to fit the AGE-PRO observations to our models using MCMC is that we can linearly interpolate fluxes between models with different parameters. For most of the parameters, the model grid covers the parameter space in several steps (for example, there are three gas-to-dust mass ratios), but for two parameters, the vertical structure and the dust properties, they are only covered by two steps, one at each extreme of the parameters range. Without any constraint in between, there is a possibility that our linear interpolation is misrepresenting the actual dependence of, for example, the ¹³CO 2–1 flux on the vertical structure of the disk.

To investigate this possibility, we ran a series of models with vertical structures in between “flat” ($h_{100}$ = 0.05, $\psi$ = 0.1) and “flared” ($h_{100}$ = 0.1, $\psi$ = 0.25) and with dust populations between “young” ($\chi$ = 0.6, $f_{large}$ = 0.8) and “evolved” ($\chi$ = 0.2, $f_{large}$ = 0.9) as defined in Section 2. Using

the vertical structure as an example, we scale the relevant parameters as $h_{100,x\%} = x \times (h_{100,flared} - h_{100,flat}) + h_{100,flat}$ and $\psi_{x\%} = x \times (\psi_{flared} - \psi_{flat}) + \psi_{flat}$. These models were run for $M_{gas} \in \{10^{-5}, 10^{-4}, 10^{-3}, 10^{-2}\} M_{\odot}$, $R_c \in \{15, 60\}$ au, an “intermediate” star and $\Delta_{gd} = 100$.

The left side of Figure 17 shows the ¹³CO 2–1, C¹⁸O 2–1, and N₂H⁺ 3–2 line fluxes as a function of vertical structure for a disk with $M_{gas} = 10^{-3} M_{\odot}$ and $R_c = 60$ au. The line fluxes follow an approximately linear trend with increasing $\psi$ and $h_{100}$, which supports our use of a linear interpolation between “flat” and “flared” disks, but it is also clear that the linear interpolation consistently underestimates the ¹³CO and C¹⁸O 2–1 fluxes. To quantify this, we calculated the difference between the flux obtained from the linear interpolation, $F_{interp}$, and true flux $F_{true}$ of a model with that specific vertical structure. The right panels of Figure 17 show this difference, expressed as a fraction of $F_{interp}$, as a function of gas disk mass. We find that the linear interpolation consistently underestimates the line fluxes by ~5%–20%. This suggests that we are underestimating line fluxes in our MCMC fits,

<figure><img src="image_21.png" /><figcaption>**Figure 17.** Left panels: ¹³CO 2–1, C¹⁸O 2–1, and N₂H⁺ 3–2 line fluxes as a function of “flaring”-ness (i.e., ($\psi_{x\%}$, $h_{100,x\%}$)) for an example disk with $M_{gas} = 10^{-3} M_{\odot}$ and $R_c = 60$ au. Black crosses denote “flat” and “flared” models that are in the model grid, with the dashed gray line showing the linear interpolation used in the MCMC. The blue circles show models with intermediate vertical structures. Right panels: the difference between the line flux obtained from linearly interpolating between “flat” and “flared” models and the true line flux for that vertical structure, expressed as the fraction of the interpolated flux. The gas mass of each model is shown on the x-axis, with each gas mass bin containing six models: two disk sizes and three intermediate vertical structures.</figcaption></figure>
which could potentially result in a small overestimation of the disk mass, as a slightly lower disk mass with the same vertical structure can now also reproduce the observations. However, this offset is much smaller than the effect that changes of $M_{\text{gas}}$ and $x_{\text{CO}}$ will have on the line fluxes; e.g., an increase in gas mass by a factor of 3 will increase the C$^{18}$O 2–1 line by a similar factor.

## Appendix E
### MCMC Gas Mass Fitting Priors

This appendix contains Table 5, which provides the details on the priors used in Section 3.4. to measure the gas disk masses of the 20 disks in the Lupus and Upper Sco star-forming regions.

<table><caption>**Table 5**
MCMC Priors</caption><thead><tr><th>AGEPRO ID</th><th>lg<sub>10</sub> M<sub>gas</sub></th><th>R<sub>c</sub></th><th>lg<sub>10</sub> Δ<sub>gd</sub></th><th>ψ</th><th>f<sub>large</sub></th><th>L<sub>*</sub><sup>a</sup></th><th>lg<sub>10</sub> x<sub>CO</sub></th><th>lg<sub>10</sub> ζ<sub>CR</sub></th></tr></thead><tbody><tr><td>Lupus 1</td><td>-5 to -1</td><td>15-60</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.87, 0.26)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>Lupus 2</td><td>-5 to -1</td><td>15-120</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.33, 0.10)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>Lupus 3</td><td>-4 to -1</td><td>15-180</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.39, 0.12)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>Lupus 4</td><td>-6 to -2</td><td>1-60</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.27, 0.08)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>Lupus 5</td><td>-6 to -2</td><td>1-60</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.59, 0.18)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>Lupus 6</td><td>-6 to -2</td><td>1-60</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.20, 0.06)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>Lupus 7</td><td>-6 to -2</td><td>1-60</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.15, 0.04)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>Lupus 8</td><td>-6 to -2</td><td>1-60</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.22, 0.07)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>Lupus 9</td><td>-5 to -2</td><td>1-120</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.27, 0.08)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>Lupus 10</td><td>-3 to -0.3</td><td>60-300</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>1.15<sup>b</sup></td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>UppSco 1</td><td>-5 to -1</td><td>15-120</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.25, 0.07)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>UppSco 2</td><td>-6 to -2</td><td>1-60</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.09, 0.03)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>UppSco 3</td><td>-5 to -2</td><td>1-60</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.08, 0.02)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>UppSco 4</td><td>-6 to -2</td><td>1-60</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.32, 0.10)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>UppSco 5</td><td>-6 to -2</td><td>1-60</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.09, 0.03)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>UppSco 6</td><td>-6 to -2</td><td>1-120</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.17, 0.05)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>UppSco 7</td><td>-4 to -1</td><td>15-180</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.18, 0.05)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>UppSco 8</td><td>-5 to -1</td><td>15-120</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.10, 0.03)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>UppSco 9</td><td>-5 to -1</td><td>5-120</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.32, 0.10)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr><tr><td>UppSco 10</td><td>-5 to -1</td><td>5-120</td><td>1-3</td><td>0.1-0.25</td><td>0.8-0.9</td><td>(0.32, 0.10)</td><td>-6.5 to -4.0</td><td>-19 to -17</td></tr></tbody></table>

**Notes.**

<footnote><sup>a</sup> Values for the Gaussian prior for stellar luminosity are the mean and standard deviation, respectively.</footnote>

<footnote><sup>b</sup> For Lupus 10, no Gaussian prior is used, but instead, only models with $L_{*} = 1.0L_{⊙}$ were used in the fit.</footnote>
Appendix F

Comparison of the AGE-PRO <sup>12</sup>CO Observations and the
Models

This appendix contains Figure 18, which shows a compar-
ison between the observed <sup>12</sup>CO 2–1 and <sup>13</sup>CO 2–1 integrated
fluxes from the 20 AGE-PRO disks in Lupus and Upper Sco
(C. Agurto-Gangas et al. 2025; D. Deng et al. 2025) and the
fluxes from the models presented in this work.

<figure><img src="image_22.png" /><figcaption>**Figure 18.** The same as Figure 3, but now showing the <sup>12</sup>CO and <sup>13</sup>CO *J* = 2–1 line luminosities. Note that due to cloud absorption seen toward most sources in Lupus, their observed <sup>12</sup>CO line luminosities should be considered as lower limits (see D. Deng et al. 2025 for details). Accounting for this effect would move the points to the right of the figure. The zoomed-in portions shown in panels (a) and (b) now focus on the sources in Lupus and Upper Sco where <sup>13</sup>CO and C<sup>18</sup>O are not detected.</figcaption></figure>
## Appendix G
### Joint Posterior $M_{\text{gas}}$ and $x_{\text{CO}}$ Probability Distributions

This appendix contains Figure 19, which shows the joint posterior probability distributions of $M_{\text{gas}}$ and $x_{\text{CO}}$ from the MCMC fits discussed in Section 3.6.

<figure><img src="image_23.png" /><figcaption>**Figure 19.** Joint posterior probability distribution of $M_{\text{gas}}$ and $x_{\text{CO}}$ obtained from fitting the observed integrated $^{13}$CO 2–1, C$^{18}$O 2–1, N$_{2}$H$^{+}$ 3–2, and 1300 µm continuum fluxes using MCMC. Orange dashed lines show the best-fit $M_{\text{gas}}$ and $x_{\text{CO}}$ for each disk.</figcaption></figure>
## Appendix H
### Ophiuchus Gas Mass Posteriors from Cold and Warm Disk Models

This appendix contains Figure 20, which shows the $M_{\text{gas}}$ posterior probability distributions of the MCMC fits of the Ophiuchus disks discussed in Section 3.5.

<figure><img src="image_24.png" /><figcaption>**Figure 20.** Comparison of the gas mass posterior distribution obtained with fitting the Ophiuchus observations using cold ($L_* = L_{*\text{,obs}}$; in blue) and warm ($L_* = 10 \times L_{*\text{,obs}}$; in red) disk models. See Section 3.5 for details.</figcaption></figure>

## Appendix I
### Derivation of the Characteristic Radius

Recently, C. Toci et al. (2023) and L. Trapman et al. (2023) demonstrated that the observed gas disk size $R_{\text{CO}}$, 90% is directly related to a critical surface density and can be written down analytically as a function of $M_{\text{gas}}$, $R_c$, ($L_*$), and the slope of the surface density $\gamma$. L. Trapman et al. (2023) used this expression and an estimate of the gas mass based on the dust mass to derive $R_c$ from the observed $R_{\text{CO}}$, 90% of disks in Lupus, Taurus, Upper Sco, and the DSHARP sample. Here, we repeat that exercise for the AGE-PRO disks in Lupus and Upper Sco using the gas masses derived in this work. Briefly, for each source, we compute $N$ estimates for $R_{\text{CO},90\%\text{,est}}$ using

Equation (9) from L. Trapman et al. (2023) using the gas mass of the source and a range of $NR_c$ values. Here, we also include the stellar luminosity dependence of critical surface density (Equation (7) in L. Trapman et al. 2023), taking the stellar luminosities of each source from the literature (J. M. Alcalá et al. 2019; C. F. Manara et al. 2020; J. M. Carpenter et al. 2025; see K. Zhang et al. 2025 for details.). The best-fitting $R_c$ is then selected by minimizing $|R_{\text{CO},90\%\text{,est}} - R_{\text{CO},90\%\text{,obs}}^{\text{fit}}|$, where $R_{\text{CO},90\%\text{,obs}}^{\text{fit}}$ is the gas disk size from L. Trapman et al. (2025), who fitted the ¹²CO moment zero map with a beam-convolved Nuker profile to obtain an $R_{\text{CO},90\%}$ corrected for the resolution of the observations (see L. Trapman et al. 2025 for details).

To include the uncertainties on $M_{\text{gas}}$ and $R_{\text{CO}}$, 90%, obs, this procedure is repeated $M$ times in a Monte Carlo fashion, each

<figure><img src="image_25.png" /><figcaption>**Figure 21.** Left panel: distribution of derived $R_c$ for Lupus 3 based on its $R_{\text{CO}}$, 90%, $M_{\text{gas}}$, and $L_*$. Middle panel: difference between the best-fitting disk size ($R_{\text{CO},90\%\text{,est}}$) and the observed value ($R_{\text{CO},90\%\text{,obs}}$). Right panel: comparison of the full posterior distribution of $M_{\text{gas}}$ and the portion sampled during the Monte Carlo process used to build up the distribution of $R_c$ in the left panel.</figcaption></figure>
<table><caption>**Table 6**
Characteristic Radii for Lupus and Upper Sco Disks</caption><thead><tr><th>Name</th><th>M<sub>gas</sub></th><th>L<sub>*</sub></th><th>R<sub>CO</sub>, 90%</th><th>R<sub>c</sub><br>(Analytical)</th><th>R<sub>c</sub><br>(MCMC)</th></tr><tr><th>Comments</th><th>lg (M<sub>CO</sub>)</th><th>(L<sub>CO</sub>)</th><th>(au)</th><th>(au)</th><th>(au)</th></tr></thead><tbody><tr><td>Lupus 1</td><td>-2.72<sup>+0.20</sup><sub>-0.17</sub></td><td>0.87</td><td>169.8<sup>+1.4</sup><sub>-2.1</sub></td><td>32<sup>+2</sup><sub>-2</sub></td><td>31<sup>+4</sup><sub>-4</sub></td></tr><tr><td>Lupus 2</td><td>-2.04<sup>+0.25</sup><sub>-0.26</sub></td><td>0.33</td><td>312.7<sup>+0.8</sup><sub>-0.6</sub></td><td>65<sup>+8</sup><sub>-6</sub></td><td>71<sup>+10</sup><sub>-9</sub></td></tr><tr><td>Lupus 3</td><td>-2.83<sup>+0.22</sup><sub>-0.24</sub></td><td>0.39</td><td>106.6<sup>+3.7</sup><sub>-9.1</sub></td><td>18<sup>+2</sup><sub>-2</sub></td><td>59<sup>+14</sup><sub>-11</sub></td></tr><tr><td>Lupus 4</td><td>-3.72<sup>+0.42</sup><sub>-0.28</sub></td><td>0.27</td><td>43.1<sup>+29.8</sup><sub>-39.2</sub></td><td>8<sup>+8</sup><sub>-5</sub></td><td>17<sup>+7</sup><sub>-5</sub></td></tr><tr><td>Lupus 5</td><td>-4.35<sup>+0.68</sup><sub>-0.35</sub></td><td>0.59</td><td>53.7<sup>+17.5</sup><sub>-9.9</sub></td><td>10<sup>+5</sup><sub>-4</sub></td><td>11<sup>+12</sup><sub>-6</sub></td></tr><tr><td>Lupus 6</td><td>-3.27<sup>+0.43</sup><sub>-0.29</sub></td><td>0.20</td><td>56.6<sup>+1.8</sup><sub>-1.6</sub></td><td>9<sup>+1</sup><sub>-1</sub></td><td>19<sup>+7</sup><sub>-5</sub></td></tr><tr><td>Lupus 7</td><td>-3.79<sup>+0.71</sup><sub>-0.38</sub></td><td>0.15</td><td>88.0<sup>+2.1</sup><sub>-1.8</sub></td><td>20<sup>+4</sup><sub>-4</sub></td><td>12<sup>+6</sup><sub>-5</sub></td></tr><tr><td>Lupus 8</td><td>-3.12<sup>+0.40</sup><sub>-0.33</sub></td><td>0.22</td><td>58.6<sup>+3.0</sup><sub>-1.6</sub></td><td>9<sup>+1</sup><sub>-1</sub></td><td>14<sup>+3</sup><sub>-4</sub></td></tr><tr><td>Lupus 9</td><td>-3.42<sup>+0.90</sup><sub>-0.71</sub></td><td>0.27</td><td>7.2<sup>+7.2</sup><sub>-7.2</sub></td><td>2<sup>+1</sup><sub>-1</sub></td><td>5<sup>+7</sup><sub>-2</sub></td></tr><tr><td>Lupus 10</td><td>-1.19<sup>+0.11</sup><sub>-0.13</sub></td><td>1.15</td><td>838.5<sup>+55.3</sup><sub>-58.8</sub></td><td>184<sup>+24</sup><sub>-21</sub></td><td>278<sup>+15</sup><sub>-21</sub></td></tr><tr><td>UppSco 1</td><td>-2.49<sup>+0.15</sup><sub>-0.11</sub></td><td>0.25</td><td>167.5<sup>+7.8</sup><sub>-2.2</sub></td><td>31<sup>+2</sup><sub>-2</sub></td><td>49<sup>+7</sup><sub>-5</sub></td></tr><tr><td>UppSco 2</td><td>-3.90<sup>+0.38</sup><sub>-0.24</sub></td><td>0.09</td><td>25.7<sup>+30.8</sup><sub>-19.4</sub></td><td>4<sup>+5</sup><sub>-3</sub></td><td>14<sup>+6</sup><sub>-5</sub></td></tr><tr><td>UppSco 3</td><td>-3.19<sup>+0.24</sup><sub>-0.43</sub></td><td>0.08</td><td>34.9<sup>+1.8</sup><sub>-4.5</sub></td><td>5<sup>+1</sup><sub>-1</sub></td><td>12<sup>+3</sup><sub>-3</sub></td></tr><tr><td>UppSco 4</td><td>-3.69<sup>+0.92</sup><sub>-1.33</sub></td><td>0.32</td><td>46.3<sup>+1.1</sup><sub>-1.5</sub></td><td>8<sup>+5</sup><sub>-2</sub></td><td>3<sup>+6</sup><sub>-1</sub></td></tr><tr><td>UppSco 5</td><td>-5.36<sup>+0.27</sup><sub>-0.20</sub></td><td>0.09</td><td>30.4<sup>+1.4</sup><sub>-1.0</sub></td><td>7<sup>+1</sup><sub>-1</sub></td><td>19<sup>+7</sup><sub>-6</sub></td></tr><tr><td>UppSco 6</td><td>-4.16<sup>+0.30</sup><sub>-0.33</sub></td><td>0.17</td><td>155.9<sup>+16.7</sup><sub>-13.9</sub></td><td>69<sup>+28</sup><sub>-19</sub></td><td>103<sup>+12</sup><sub>-18</sub></td></tr><tr><td>UppSco 7</td><td>-3.38<sup>+0.29</sup><sub>-0.23</sub></td><td>0.18</td><td>160.5<sup>+1.4</sup><sub>-1.5</sub></td><td>44<sup>+7</sup><sub>-6</sub></td><td>97<sup>+27</sup><sub>-23</sub></td></tr><tr><td>UppSco 8</td><td>-2.41<sup>+0.27</sup><sub>-0.35</sub></td><td>0.10</td><td>144.1<sup>+1.7</sup><sub>-2.5</sub></td><td>26<sup>+3</sup><sub>-2</sub></td><td>36<sup>+5</sup><sub>-5</sub></td></tr><tr><td>UppSco 9</td><td>-1.28<sup>+0.20</sup><sub>-0.38</sub></td><td>0.32</td><td>181.8<sup>+5.2</sup><sub>-8.5</sub></td><td>25<sup>+3</sup><sub>-2</sub></td><td>19<sup>+7</sup><sub>-4</sub></td></tr><tr><td>UppSco 10</td><td>-2.14<sup>+0.43</sup><sub>-0.45</sub></td><td>0.32</td><td>75.0<sup>+1.8</sup><sub>-4.0</sub></td><td>9<sup>+1</sup><sub>-1</sub></td><td>19<sup>+6</sup><sub>-4</sub></td></tr></tbody></table>

time drawing a random $M_{gas,i}$ from its posterior probability distribution (see Section 3.4) and a random $R_{CO,90\%\_obs}^{\text{fit},i}$ from within its uncertainties (see Table 2 in L. Trapman et al. 2025). Figure 21 shows an example of the resulting distribution of $R_c$ for Lupus 3, and estimates of $R_c$ for all 20 sources can be found in Table 6.

Another estimate for $R_c$ comes from the MCMC fit to the CO isotopologue line fluxes and $R_{CO}$. 90% discussed in Section 3.4. The characteristic radius is one of the model parameters we marginalize over to obtain the gas mass posterior distribution, but we can also marginalize over all parameters except $R_c$ to obtain its posterior distribution (see Table 6).

Before comparing the two characteristic radii, it is useful to first briefly discuss the main difference between the analytical and the MCMC-derived $R_c$. Both are, in their own way, linked to the observed gas disk size $R_{CO}$, 90%. In the case of the analytical $R_c$, $R_{CO}$, 90% is obtained from a best-fit model intensity profile fitted to the $^{12}\text{CO}$ 2–1 moment zero map. For the fit, the model is convolved with the clean beam of the observations. $R_{CO}$, 90% is measured from the unconvolved model to correct for the effect of the resolution of the observations on the measurement of $R_{CO}$, 90% (see L. Trapman et al. 2025 for details). For the MCMC-$R_c$, the approach is more forward-modeling, where the synthetic observations of the DALI models were convolved with the clean beam of the observations before being compared to the observed $R_{CO}$, 90% (see Section 3.4). While these two approaches are different, to first order we expect them to provide a similar correction of $R_{CO}$, 90% for the resolution of the observations.

The largest difference lies with the fact that the MCMC-$R_c$ is compared to the $^{13}\text{CO}$ and $C^{18}\text{O}$ line fluxes. As can be seen

in panel (b) in Figure 3, in particular the $^{13}\text{CO}/C^{18}\text{O}$ 2–1 line ratio is sensitive to $R_c$ for small $R_c$, with the line ratio decreasing toward 1 as $R_c$ decreases. The main cause for this is optical depth. For these disks, the $^{13}\text{CO}$ and $C^{18}\text{O}$ integrated line fluxes are dominated by optically thick emission, meaning that the line ratio scales as the ratio of their emitting areas. Due to its lower abundance of $C^{18}\text{O}$, its emission will become optically thin at a higher column density, and therefore a smaller radius, than the more abundant $^{13}\text{CO}$, meaning that the optically thick emitting area of $C^{18}\text{O}$ is always smaller than that of $^{13}\text{CO}$. However, if the slope of the surface density profile is steep, such as in the exponential taper of the surface density profile, the radii at which $C^{18}\text{O}$ and $^{13}\text{CO}$ become optically thick will lie close together, and the line ratio will be close to 1. Decreasing $R_c$ moves these two radii farther out into the exponential taper and thereby closer together, resulting in the line ratio decreasing toward 1.

Figure 22 compares the MCMC-derived $R_c$ to the analytically computed $R_c$ discussed previously. In almost all cases, the two radii agree to within a factor of 2, although for the more compact disks ($R_c \lesssim 10$ au), the two estimates start to diverge. Among the larger disks, Lupus 3 is a clear outlier with MCMC-$R_c$ that is $\sim 3 \times$ larger than its analytical $R_c$. This source has a $^{13}\text{COC}^{18}\text{O}$ 2–1 line flux ratio of $\sim 16$ (D. Deng et al. 2025), which is too large to be consistent with the smaller $R_c$ obtained from the analytical calculation. The discrepancy between these two characteristic radii could be an indication that the surface density of this disk does not follow an exponential taper but has a more shallow decline with radius.
<figure><img src="image_26.png" /><figcaption>**Figure 22.** Comparison of the characteristic radii ($R_c$) obtained from the MCMC fit described in Section 3.4 to the analytically computed $R_c$ discussed in Appendix I.</figcaption></figure>
## Appendix J
### Extra Figures of the Correlations between Disk and Stellar Properties

This appendix contains Figures 23–28, which show the correlations between disk and stellar parameters discussed in Section 4.1.

<figure><img src="image_27.png" /><figcaption>**Figure 23.** Gas mass vs. effective temperature ($T_{\text{eff}}$), stellar mass ($M_{*}$), stellar luminosity ($L_{*}$), gas disk outer radius ($R_{\text{CO,90\%}}$), dust disk outer radius ($R_{\text{dust,90\%}}$), and the characteristic radius $R_{c}$ for the 20 disks in Lupus (black) and Upper Sco (purple). See Appendix I for details on how $R_{c}$ was derived for each source.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>**Figure 24.** Results of the Monte Carlo Spearman rank tests, quantified by $\rho$, between $M_{\text{gas}}$ and various stellar and disk parameters. The left panels show the distribution of $\rho$, and the right panels show the associated $p$-value. Orange vertical lines in each panel show the 16th, 50th, and 84th quantiles of the distribution. The red vertical line denotes $p = 0.05$ to provide a guide for how statistically significant the correlation is.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Figure 25.** The same as Figure 23, but comparing $\Delta_{gd}$ to other disk and stellar parameters.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 26. The same as Figure 24, but for correlations with $\Delta_{gd}$.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Figure 27.** The same as Figure 23, but comparing $x_{CO}$ to other disk and stellar parameters.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 28. The same as Figure 24, but for correlations with $x_{CO}$.</figcaption></figure>

## Appendix K
### Corner Plots of the MCMC

This appendix contains Figures 29.1–20, which show corner plots of the posterior probability distributions of the model

parameters obtained from fitting observed integrated $^{13}$CO 2–1, C$^{18}$O 2–1, N$_{2}$H$^{+}$ 3–2, and 1300 µm continuum fluxes and the gas disk size ($R_{CO}$, 90%) using MCMC for all 20 disks in the Lupus and Upper Sco star-forming regions (see Section 3.4.).
<figure><img src="image_33.png" /><figcaption>**Figure 29.** Posterior probability distributions of the model parameters obtained from fitting the observed integrated $^{13}$CO 2–1, C$^{18}$O 2–1, N$_{2}$H$^{+}$ 3–2, and 1300 µm continuum fluxes and the gas disk size ($R_{\text{CO}}$, 90%) of Lupus 1 (Sz 65) using MCMC. Vertical dashed lines show the 16th, 50th, and 84th quantiles of the distribution.</figcaption></figure>

(The complete figure set (20 images) is available in the online article.)

ORCID iDs

Leon Trapman https://orcid.org/0000-0002-8623-9703

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

Giovanni P. Rosotti https://orcid.org/0000-0003-4853-5736

Paola Pinilla https://orcid.org/0000-0001-8764-1780

Benoît Tabone https://orcid.org/0000-0002-1103-3225

Ilaria Pascucci https://orcid.org/0000-0001-7962-1683

Carolina Agurto-Gangas https://orcid.org/0000-0002-7238-2306

Rossella Anania https://orcid.org/0009-0004-8091-5055

John Carpenter https://orcid.org/0000-0003-2251-0602

Lucas A. Cieza https://orcid.org/0000-0002-2828-1153

Dingshan Deng https://orcid.org/0000-0003-0777-7392

Camilo González-Ruilova https://orcid.org/0000-0003-4907-189X

Michiel R. Hogerheijde https://orcid.org/0000-0001-5217-537X

Nicolás T. Kurtovic https://orcid.org/0000-0002-2358-4796

Aleksandra Kuznetsova https://orcid.org/0000-0002-6946-6787

James Miley https://orcid.org/0000-0002-1575-680X

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

Dary A. Ruiz-Rodriguez https://orcid.org/0000-0003-3573-8163

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

Anibal Sierra https://orcid.org/0000-0002-5991-8073
Estephani TorresVillanueva https://orcid.org/0000-0001-9961-8203

Miguel Vioque https://orcid.org/0000-0002-4147-3846

## References

Agurto-Gangas, C., Pérez, L. M., Sierra, A., et al. 2025, ApJ, 989, 4

Aikawa, Y., Cataldi, G., Yamato, Y., et al. 2021, ApJ, 257, 10

Aikawa, Y., Umebayashi, T., Nakano, T., & Miyama, S. M. 1997, ApJL, 486, L51

Aikawa, Y., van Zadelhoff, G., van Dishoeck, E. F., & Herbst, E. 2002, A&A, 386, 622

Alcalá, J., Manara, C., Natta, A., et al. 2017, A&A, 600, A20

Alcalá, J., Natta, A., Manara, C., et al. 2014, A&A, 561, A2

Alcalá, J. M., Manara, C. F., France, K., et al. 2019, A&A, 629, A108

Altwegg, K., Balsiger, H., Hänni, N., et al. 2020, NatAs, 4, 533

Anania, R., Rosotti, G. P., Gárate, M., et al. 2025, ApJ, 989, 8

Anderson, D. E., Blake, G. A., Bergin, E. A., et al. 2019, ApJ, 881, 127

Anderson, D. E., Cleeves, L. I., Blake, G. A., et al. 2022, ApJ, 927 229

Andrews, S. M. 2020, ARA&A, 58, 483

Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129

Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018, ApJ, 865, 157

Andrews, S. M., Wilner, D. J., Espaillat, C., et al. 2011, ApJ, 732, 42

Andrews, S. M., Wilner, D. J., Hughes, A., et al. 2011, ApJ, 744, 162

Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21

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

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

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

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

Bae, J., Isella, A., Zhu, Z., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 423

Ballering, N. P., & Eisner, J. A. 2019, AJ, 157, 144

Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142

Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924

Bergin, E., Hogerheijde, M., Brinch, C., et al. 2010, A&A, 521, L33

Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 10

Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Natur, 493, 644

Bergin, E. A., & Williams, J. P. 2017, in Astrophysics and Space Science Library, Vol. 445, Formation, Evolution, and Dynamics of Young Solar Systems, ed. M. Pessah & O. Gressel (Cham: Springer), 1

Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, ApJ, 876, 25

Birnstiel, T. 2024, ARA&A, 62, 157

Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148

Bitsch, B., & Izidoro, A. 2023, A&A, 674, A178

Boogert, A., Gerakines, P., & Whittet, D. 2015, ARA&A, 53, 541

Bosman, A. D., Trapman, L., Sturm, A., et al. 2022, RNAAS, 6, 176

Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182

Briceño-Morales, G., & Chanamé, J. 2023, MNRAS, 522, 1288

Bruderer, S. 2013, A&A, 559, A46

Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91

Calahan, J. K., Bergin, E., Zhang, K., et al. 2021, ApJ, 908, 8

Carmona, A., van der Plas, G., van den Ancker, M. E., et al. 2011, A&A, 533, A39

Carpenter, J. M., Esplin, T. L., Luhman, K. L., Mamajek, E. E., & Andrews, S. M. 2025, ApJ, 978, 117

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

Cieza, L. A., González-Ruilova, C., Hales, A. S., et al. 2021, MNRAS, 501, 2934

Cieza, L. A., Ruiz-Rodriguez, D., Hales, A., et al. 2019, MNRAS, 482, 698

Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2018, ApJ, 865, 155

Comerón, F. 2008, in Handbook of Star Forming Regions, Volume II, The Southern Sky, ed. Bo Reipurth (San Francisco, CA: ASP), 295

D'Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411

Deng, D., Vioque, M., & Pascucci, I. 2025, ApJ, 989, 3

Deng, D., Ruaud, M., Gorti, U., & Pascucci, I. 2023, ApJ, 954, 165

Drażkowska, J., Bitsch, B., Lambrechts, M., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII Vol. 534 ed. S. Inutsuka et al. (San Francisco, CA: ASP), 717

Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971

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

Dullemond, C. P., van Zadelhoff, G. J., & Natta, A. 2002, A&A, 389, 464

Dutrey, A., Guilloteau, S., Duvert, G., et al. 1996, A&A, 309, 493

Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 613, A14

Evans Neal, J. I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321
Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 605, A16

Fang, M., Pascucci, I., Edwards, S., et al. 2023, ApJ, 945, 112

Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJL, 776, L38

Fischer, W. J., Hillenbrand, L. A., Herczeg, G. J., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 355

Flores, C., Ohashi, N., Tobin, J. J., et al. 2023, ApJ, 958, 98

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

Furuya, K., & Aikawa, Y. 2014, ApJ, 790, 97

Furuya, K., Lee, S., & Nomura, H. 2022, ApJ, 938, 29

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Galli, P. A. B., Bouy, H., Olivares, J., et al. 2020, A&A, 643, A148

Goldsmith, P. F., Bergin, E. A., & Lis, D. C. 1997, ApJ, 491, 615

Grant, S. L., Stapper, L. M., Hogerheijde, M. R., et al. 2023, AJ, 166, 147

Han, F., Hartmann, L., Calvet, N., & Franco-Hernández, R. 2023, MNRAS, 526, 2566

Hartmann, L. 1998, Accretion processes in Star formation, vol. 32 (Cambridge: Cambridge Univ. Press)

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135

Hayashi, C. 1981, PThPS, 70, 35

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJL, 869, L42

Hunter, J. D. 2007, CSE, 9, 90

Jones, M. G., Pringle, J. E., & Alexander, R. D. 2012, MNRAS, 419, 925

Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83

Kama, M., Folsom, C. P., & Pinilla, P. 2015, A&A, 582, L10

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

Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134

Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78

Kurtovic, N. T., Gárate, M., & Pinilla, P. 2025, ApJ, 989, 6

Kurtovic, N. T., Pinilla, P., Long, F., et al. 2021, A&A, 645, A139

Kuznetsova, A., Bae, J., Hartmann, L., & Mac Low, M.-M. 2022, ApJ, 928, 92

Lacy, J., Knacke, R., Geballe, T., & Tokunaga, A. 1994, ApJL, 428, L69

Law, C. J., Crystian, S., Teague, R., et al. 2022, ApJ, 932, 114

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

Lin, Z.-Y. D., Li, Z.-Y., Tobin, J. J., et al. 2023, ApJ, 951, 9

Lodato, G., Scardoni, C. E., Manara, C. F., & Testi, L. 2017, MNRAS, 472, 4700

Long, F., Herczeg, G. J., Harsono, D., et al. 2019, ApJ, 882, 49

Long, F., Herczeg, G. J., Pascucci, I., et al. 2017, ApJ, 844, 99

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A, 648, A33

Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 539

Manara, C. F., Fedele, D., Herczeg, G. J., & Teixeira, P. S. 2016, A&A, 585, A136

Manara, C. F., Natta, A., Rosotti, G. P., et al. 2020, A&A, 639, A58

Maret, S., Maury, A. J., Belloche, A., et al. 2020, A&A, 635, A15

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, AJ, 831, 167

Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96

Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 501

Miotello, A., Rosotti, G., Ansdell, M., et al. 2021, A&A, 651, A48

Miotello, A., van Dishoeck, E., Williams, J., et al. 2017, A&A, 599, A113

Miotello, A., van Dishoeck, E. F., Kama, M., & Bruderer, S. 2016, A&A, 594, A85

Morbidelli, A., & Raymond, S. N. 2016, JGRE, 121, 1962

Najita, J. R., & Bergin, E. A. 2018, ApJ, 864, 168

Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Oka, T. 2006, PNAS, 103, 12235

Paardekooper, S., Dong, R., Duffell, P., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 685
Padovani, M., Marcowith, A., Hennebelle, P., & Ferrière, K. 2016, A&A, 590, A8

Paneque-Carreño, T., Miotello, A., van Dishoeck, E. F., et al. 2023, A&A, 669, A126

Pascucci, I., Apai, D., Luhman, K., et al. 2009, ApJ, 696, 143

Pascucci, I., Herczeg, G., Carr, J., & Bruderer, S. 2013, ApJ, 779, 178

Pascucci, I., Skinner, B. N., Deng, D., et al. 2023, ApJ, 953, 183

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

Pecaut, M. J., Mamajek, E. E., & Bubar, E. J. 2012, ApJ, 746, 154

Pegues, J., Öberg, K. I., Bergner, J. B., et al. 2021, ApJ, 911, 150

Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114

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

Pittman, C. V., Espaillat, C. C., Robinson, C. E., et al. 2022, AJ, 164, 201

Powell, D., Gao, P., Murray-Clay, R., & Zhang, X. 2022, NatAs, 6, 1147

Qi, C., Öberg, K. I., Espaillat, C. C., et al. 2019, ApJ, 882, 160

Qi, C., Öberg, K. I., Wilner, D. J., et al. 2013, Sci, 341, 630

Ratzenböck, S., Großschedl, J. E., Alves, J., et al. 2023, A&A, 678, A71

Reboussin, L., Wakelam, V., Guilloteau, S., Hersant, F., & Dutrey, A. 2015, A&A, 579, A82

Ribas, Á., Espaillat, C. C., Macías, E., & Sarro, L. M. 2020, A&A, 642, A171

Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15

Röllig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187

Rosotti, G. P., Clarke, C. J., Manara, C. F., & Facchini, S. 2017, MNRAS, 468, 1631

Rosotti, G. P., Tazzari, M., Booth, R. A., et al. 2019, MNRAS, 486, 4829

Ruaud, M., & Gorti, U. 2019, ApJ, 885, 146

Ruaud, M., Gorti, U., & Hollenbach, D. J. 2022, ApJ, 925, 49

Ruiz-Rodríguez, D. A., González-Ruilova, C., Cieza, L. A., et al. 2025, ApJ, 989, 2

Salinas, V. N., Hogerheijde, M. R., Bergin, E. A., et al. 2016, A&A, 591, A122

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2018, ApJ, 856, 85

Schwarz, K. R., Calahan, J. K., Zhang, K., et al. 2021, ApJS, 257, 20

Sierra, A., Pérez, L. M., Agurto-Gangas, C., et al. 2024, ApJ, 974, 102

Sturm, J. A., Booth, A. S., McClure, M. K., Leemker, M., & van Dishoeck, E. F. 2023, A&A, 670, A12

Tabone, B., Rosotti, G. P., Cridland, A. J., Armitage, P. J., & Lodato, G. 2022, MNRAS, 512, 2290

Tabone, B., Rosotti, G. P., Trapman, L., et al. 2025, ApJ, 989, 7

Takeuchi, T., Clarke, C. J., & Lin, D. N. C. 2005, ApJ, 627, 286

Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 339

Testi, L., Natta, A., Manara, C. F., et al. 2022, A&A, 663, A98

Thi, W. F., van Dishoeck, E. F., Blake, G. A., et al. 2001, ApJ, 561, 1074

Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130

Toci, C., Lodato, G., Livio, F. G., Rosotti, G., & Trapman, L. 2023, MNRAS, 518, L69

Trapman, L., Bosman, A. D., Rosotti, G., Hogerheijde, M. R., & van Dishoeck, E. F. 2021, A&A, 649, A95

Trapman, L., Miotello, A., Kama, M., van Dishoeck, E. F., & Bruderer, S. 2017, A&A, 605, A69

Trapman, L., Rosotti, G., Bosman, A. D., Hogerheijde, M. R., & van Dishoeck, E. F. 2020, A&A, 640, A5

Trapman, L., Rosotti, G., Zhang, K., & Tabone, B. 2023, ApJ, 954, 41

Trapman, L., Tabone, B., Rosotti, G., & Zhang, K. 2022a, ApJ, 926, 61

Trapman, L., Vioque, M., Kurtovic, N. T., et al. 2025, ApJ, 989, 10

Trapman, L., Zhang, K., van't Hoff, M. L. R., Hogerheijde, M. R., & Bergin, E. A. 2022b, ApJL, 926, L2

Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19

Van Clepper, E., Bergner, J. B., Bosman, A. D., Bergin, E., & Ciesla, F. J. 2022, ApJ, 927, 206

van der Tak, F. F. S., & van Dishoeck, E. F. 2000, A&A, 358, L79

van't Hoff, M. L. R., Tobin, J. J., Harsono, D., & van Dishoeck, E. F. 2018, A&A, 615, A83

van't Hoff, M. L. R., Walsh, C., Kama, M., Facchini, S., & van Dishoeck, E. F. 2017, A&A, 599, A101

van Terwisga, S. E., Hacar, A., van Dishoeck, E. F., Oonk, R., & Portegies Zwart, S. 2022, A&A, 661, A53

van Terwisga, S. E., van Dishoeck, E. F., Ansdell, M., et al. 2018, A&A, 616, A88

van't Hoff, M. L. R., Harsono, D., Tobin, J. J., et al. 2020, ApJ, 901, 166

Villenave, M., Podio, L., Duchêne, G., et al. 2023, ApJ, 946, 70

Vioque, M., Kurtovic, N. T., Trapman, L., et al. 2025, ApJ, 989, 9

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

Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Weingartner, J. C., & Draine, B. 2001, ApJ, 548, 296

Whitney, B. A., Robitaille, T. P., Bjorkman, J. E., et al. 2013, ApJS, 207, 30

Williams, J. P., Andrews, S. M., & Wilner, D. J. 2005, ApJ, 634, 495

Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59

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

Yu, M., Evans Neal, J. I., Dodson-Robinson, S. E., Willacy, K., & Turner, N. J. 2017, ApJ, 841, 39

Yu, M., Willacy, K., Dodson-Robinson, S. E., Turner, N. J., & Evans Neal, J. I. 2016, ApJ, 822, 53

Zagaria, F., Facchini, S., Miotello, A., et al. 2023, A&A, 672, L15

Zhang, K., Bergin, E. A., Schwarz, K., Krijt, S., & Ciesla, F. 2019, ApJ, 883, 98

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

Zhang, K., Pérez, L. M., Pascucci, I., et al. 2025, ApJ, 989, 1

Zhang, K., Schwarz, K. R., & Bergin, E. A. 2020, ApJL, 891, L17

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