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

# The ALMA Survey of Gas Evolution of PROtoplanetary Disks (AGE-PRO). VII. Testing Accretion Mechanisms from Disk Population Synthesis

Benoît Tabone¹ ID, Giovanni P. Rosotti² ID, Leon Trapman³ ID, Paola Pinilla⁴ ID, Ilaria Pascucci⁵ ID, Alice Somigliana⁶ ID,
Richard Alexander⁷ ID, Miguel Vioque⁶,⁸ ID, Rossella Anania² ID, Aleksandra Kuznetsova⁹ ID, Ke Zhang³ ID, Laura M. Pérez¹⁰ ID,
Lucas A. Cieza¹¹ ID, John Carpenter⁸ ID, Dingshan Deng⁵ ID, Carolina Agurto-Gangas¹⁰ ID, Dary A. Ruiz-Rodriguez¹² ID,
Anibal Sierra¹⁰ ID, Nicolás T. Kurtovic¹³,¹⁴ ID, James Miley¹⁵,¹⁶,¹⁷ ID, Camilo González-Ruilova¹¹,¹⁶,¹⁷ ID,
Estephani Torres-Villanueva³ ID, Michiel R. Hogerheijde¹⁸,¹⁹ ID, Kamber Schwarz¹⁴ ID, Claudia Toci⁶,²⁰ ID,
Leonardo Testi²¹ ID, and Giuseppe Lodato² ID

¹ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France; benoit.tabone@cnrs.fr

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

³ Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St., Madison, WI 53706, USA

⁴ Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey RH5 6NT, UK
⁵ Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA

⁶ European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching bei München, Germany

⁷ School of Physics & Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK

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

⁹ Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Ave., New York, NY 10025, USA

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

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

¹² National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA

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

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

¹⁵ Departamento de Física, Universidad de Santiago de Chile, Av. Victor Jara 3659, 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 de Chile, Chile

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

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

²⁰ INAF, Osservatorio Astrofísico di Arcetri, 50125, Firenze, Italy

²¹ Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, I-40122 Bologna, Italy

Received 2024 September 17; revised 2025 February 17; accepted 2025 February 18; published 2025 July 31

## Abstract

The architecture of planetary systems depends on the evolution of the disks in which they form. In this work, we develop a population synthesis approach to interpret the Atacama Large Millimeter/submillimeter Array survey of Gas Evolution of PROtoplanetary Disks (AGE-PRO) measurements of disk gas mass and size considering two scenarios: turbulence-driven evolution with photoevaporative winds and MHD wind-driven evolution. A systematic method is proposed to constrain the distribution of disk parameters from the disk fractions, accretion rates, disk gas masses, and CO gas sizes. We find that turbulence-driven accretion with initially compact disks ($R_0 \simeq 5-20$ au), low mass-loss rates, and relatively long viscous timescales ($t_{v,0} \simeq 0.4-3$ Myr or $\alpha_{SS} \simeq 2-4 \times 10^{-4}$) can reproduce the disk fractions and gas sizes. However, the distribution of apparent disk lifetimes defined as the $M_D/M_*$ ratio is severely overestimated by turbulence-driven models. On the other hand, MHD wind-driven accretion can reproduce the bulk properties of disk populations from Ophiuchus to Upper Scorpius assuming compact disks with an initial magnetization of about $\beta \simeq 10^5$ ($\alpha_{DW} \simeq 0.5-1 \times 10^{-3}$) and a magnetic field that declines with time. More studies are needed to confirm the low masses found by AGE-PRO, notably for compact disks that question turbulence-driven accretion. The constrained synthetic disk populations can now be used for realistic planet population models to interpret the properties of planetary systems on a statistical basis.

*Unified Astronomy Thesaurus concepts: Exoplanet formation (492); Protoplanetary disks (1300)*

## 1. Introduction

With more than 5000 exoplanets detected so far,<fnref n="22" /> planet formation appears to be a widespread phenomenon in disks orbiting young stars (R. B. Fernandes et al. 2019). The more

striking results obtained by the exoplanet community are the large diversity in the properties of exoplanetary systems and the occurrence of specific kinds of systems. The next challenge is to understand how different evolutionary pathways of planet-forming disks lead to the diversity and occurrence of planetary systems (A. Morbidelli & S. N. Raymond 2016).

Recent models have already attempted to include many steps of planet formation in disk evolution models to predict the properties of exoplanets at the population level (J. Drazkowska et al. 2023). Yet the evolution of gaseous disks remains vastly unknown despite its paramount importance in every step of planet formation (A. Miotello et al. 2023). One of

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

<footnote n="22">²² As of 2024 April: https://exoplanets.nasa.gov/.</footnote>
the main limiting factors is the need for unbiased surveys of
young stars of different ages since an individual disk cannot be
followed up across its lifetime. Pioneering surveys have
already unveiled some of the basic properties of disks. IR
excess and accretion signatures in clusters of different ages
indicate that disks disperse quickly after a typical timescale of
2–5 Myr (D. Fedele et al. 2010; Á. Ribas et al. 2014),
depending on the mass of the host star and the vicinity of
massive stars (Á. Ribas et al. 2015; A. J. Winter &
T. J. Haworth 2022). Surveys from the UV to the near-IR
have also shown that almost all disks exhibit signs of accretion
with accretion rates scaling with stellar mass and declining
with time (L. Hartmann et al. 2016; J. M. Alcalá et al. 2017;
C. F. Manara et al. 2017; M. Fang et al. 2023).

Over the past decades, disk dispersal and disk accretion
have been explained by two distinct processes. In the classical
picture of viscous disks, popularized by N. I. Shakura &
R. A. Sunyaev (1973), disks accrete via the radial transport
of angular momentum mediated by turbulence. The origin of
turbulence can be the magnetorotational instability (MRI) or
hydrodynamical instabilities such as the vertical shear
instability or the gravitational instability (K. Kratter &
G. Lodato 2016; G. R. J. Lesur et al. 2023). On the other
hand, disk dispersal is explained by photoevaporative winds
launched from the disk atmosphere heated by energetic
radiation from the star (D. Hollenbach et al. 1994;
R. Alexander et al. 2014). The relevant energy of the photons
remains debated. Hydrogen ionizing UV photons (EUV)
would have a limited impact due to efficient screening of the
disk layers but either X-ray or far-ultraviolet (FUV) photons
can penetrate the disk's upper layers and deposit enough
energy to drive a thermal wind outward of 1–10 au (U. Gorti &
D. Hollenbach 2009; J. E. Owen et al. 2010; R. Nakatani et al.
2018). Early models combining turbulence-driven accretion
and photoevaporation were able to account for the observed
accretion rates and disk dispersal times (e.g., R. D. Alexander
et al. 2006; J. E. Owen et al. 2012).

However, theoretical works including the weak coupling
between the magnetic field and the gas show that MRI cannot
be triggered in large portions of disks called dead zones
(C. F. Gammie 1996; X.-N. Bai & J. M. Stone 2011; T. N. Delage
et al. 2022). Alternatively, MHD disk winds can efficiently drive
accretion across both the dead zones and MRI active regions by
extracting angular momentum vertically (R. D. Blandford &
D. G. Payne 1982; J. Ferreira 1997; X.-N. Bai & J. M. Stone
2013). MHD disk winds would naturally account for the
observations of winds revealed by a wide range of tracers from
the millimeter to the optical (e.g., B. Tabone et al. 2017; F. Louvet
et al. 2018; I. Pascucci et al. 2025). The only requirement is that
the disk is embedded in a magnetic field with a nonvanishing
magnetic flux. It should be noted that while early MHD disk-wind
models used strongly magnetized disks, resulting in very low disk
masses, recent theoretical work demonstrated the ability of MHD
disk winds to operate even with a midplane magnetization as low
as $\beta = 10^5$ where the $\beta$ parameter is the ratio between thermal and
magnetic pressure (W. Béthune et al. 2017). Disk dispersal in
MHD wind-driven accretion remains an open question due to the
lack of numerical simulations running on long timescales. Still,
simple disk evolution models show that if the magnetic field
strength does not decline too fast, disks can be dispersed in a short
dispersal phase (P. J. Armitage et al. 2013; X.-N. Bai 2016;
B. Tabone et al. 2022a).

Today, our lack of knowledge of which mechanism drives
disk accretion and dispersal constitutes a fundamental hurdle
in our understanding of planet formation and migration. If
accretion is driven by MHD disk winds, every step of planet
formation will be impacted. To name a few examples, MHD
disk winds can create pressure bumps where dust grows
(S. S. Suriano et al. 2018; A. Riols et al. 2020), pebble
accretion is expected to be amplified in low-turbulence disks
(T. E. Yap & K. Batygin 2024), and planet migration patterns
can be profoundly impacted by different local gas flows
(E. Lega et al. 2022; G. Wafflard-Fernandez & G. Lesur 2023)
with possibly outward Type II planet migration for highly
magnetized disks (C. N. Kimmig et al. 2020). 1D disk
evolution models, under specific prescriptions, also predict the
formation of a gas depletion in the inner disk, altering Type I
migration (M. Ogihara et al. 2015, 2018), a proposal that calls
for global numerical simulations running on secular time-
scales. In other words, if accretion is effectively driven by
MHD winds, planet-formation scenarios and planet population
synthesis models need to be revised.

However, IR excess or accretion signatures are indirect
probes of the evolution of the bulk disk. It is only with the
Atacama Large Millimeter/submillimeter Array (ALMA) that
fundamental properties like disk masses and sizes can now be
probed in statistically representative samples. The first surveys
focused on millimeter dust emission, giving access to the
continuum flux and continuum size in close by star-forming
regions (e.g., S. M. Andrews & J. P. Williams 2005; L. Ricci
et al. 2010; S. M. Andrews et al. 2013; S. A. Barenfeld et al.
2016; M. Ansdell et al. 2017; P. Cazzoletti et al. 2019;
S. E. van Terwisga et al. 2022). It is unclear though to what
extent dust millimeter flux can inform us about the gas.
Because of dust radial drift, one would naturally expect the
continuum disk size and the continuum luminosity to be poor
proxies for the disk gas size and mass (L. Testi et al. 2014;
A. D. Sellek et al. 2020; C. Toci et al. 2021; F. Zagaria et al.
2022). Yet, using millimeter continuum luminosities,
C. F. Manara et al. (2016) found a correlation between disk
mass (assumed to be 100 times the dust disk mass) and
accretion rate with a typical ratio of $M_D/\dot{M}_* \simeq 1-3$ Myr.

Surveys at the population level offer a unique opportunity to
challenge disk evolution scenarios and recent ALMA observa-
tions have triggered a renewed interest in disk evolution
models. G. Lodato et al. (2017) and G. D. Mulders et al. (2017)
showed that the $M_D - \dot{M}_*$ correlation can be naturally
reproduced by viscous models with a significant spread in
the viscous timescale and disk age. Yet, these models did not
consider disk dispersal, which can affect the $M_D - \dot{M}_*$
correlation (A. Somigliana et al. 2020). When photoevapora-
tion is included, viscous models can still reproduce the overall
$M_D - \dot{M}_*$ relation but remain dependent on dust evolution
models (C. F. Manara et al. 2019; A. D. Sellek et al. 2020;
A. Emsenhuber et al. 2023). The development of 1D disk
evolution models including MHD disk winds (X.-N. Bai 2016;
T. K. Suzuki et al. 2016; B. Tabone et al. 2022a) now permits
confronting the predictions of MHD wind-driven accretion
with disk demographics. In this context, B. Tabone et al.
(2022b) showed that simple MHD wind-driven evolution
models can naturally reproduce the $M_D - \dot{M}_*$ correlation and
the rapidity of disk dispersal.

Today, the ALMA survey of Gas Evolution of PROto-
planetary Disks (AGE-PRO; K. Zhang et al. 2025), gives
access to gas evolution thanks to deep observations of molecular lines over a sample of 30 disks located in three star-forming regions of different ages: Ophiuchus (Oph; D. A. Ruiz-Rodriguez et al. 2025), Lupus (D. Deng et al. 2025), and Upper Scorpius (Upper Sco; C. Agurto-Gangas et al. 2025). The outcome of AGE-PRO is twofold. First, the inclusion of $N_2H^+$ and CO isotopolgues in the spectral setup permits estimates of disk gas masses (D. E. Anderson et al. 2019; L. Trapman et al. 2022) using extensive grids of thermochemical models (L. Trapman et al. 2025a). The disk mass estimates are therefore not hampered by assumptions of the dust's optical properties and gas-to-dust ratio. Second, the CO disk sizes, a key diagnostic for disk physical size (C. Toci et al. 2023; L. Trapman et al. 2023), are now measured taking into account beam convolution effects (L. Trapman et al. 2025b).

In this paper, we develop a disk population synthesis approach to interpret the results of the ALMA AGE-PRO large program. A systematic approach is proposed to constrain the distribution of disk parameters from not only the disk mass and size inferred by AGE-PRO but also accretion properties and disk fraction. In Section 2, the disk evolution models are presented along with observational constraints and the population synthesis approach. The fitting procedure is then applied to each scenario in Section 3. The implications of the present work and its limitations are discussed in Section 4 and the findings are summarized in Section 5.

## 2. Model

### 2.1. Individual Disks

The evolution of an individual disk is computed using a 1D formalism, where the MHD equations are integrated in the vertical direction. The radial transport of angular momentum is described following the N. I. Shakura & R. A. Sunyaev (1973) parameterization with the generalization proposed by B. Tabone et al. (2022a). Our choice of the parameterization of disk evolution has the major advantage of facilitating the comparison between turbulence-driven and MHD wind-driven evolution. The values of the parameters can also be easily compared with the results of numerical simulations. The parameters of the model are summarized in Table 1. In both scenarios, the evolution of an individual disk is set by four parameters: initial disk mass $M_0$ and radius $R_0$, the generalized $\alpha$ parameter (or equivalently the accretion timescale $t_{acc,0}$ or the viscous timescale $t_{v,0}$), and a parameter quantifying the wind mass-loss rate. For the MHD disk-wind model, an additional parameter, $\omega$, controls the secular evolution of the magnetic field. We only consider isolated disks that are not fed by any envelope or streamer. The time $t=0$ can therefore be considered as the end of the Class I phase, and as several of the AGE-PRO Oph targets are Class I, we assume that the Oph disks represent the time zero in our models. Neglecting the impact of early accretion is further justified by the mass of the residual envelope around the Oph sources being smaller than the disk masses (D. A. Ruiz-Rodriguez et al. 2025). We further neglect the impact of mild external UV irradiation, which is studied in the case of turbulence-driven accretion by R. Anania et al. (2025) for the Upper Sco AGE-PRO sample.

<table><caption>**Table 1**
Summary of the Parameters of the Population Models with the Free Parameters Outlined in Bold Font</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Range of median value</th><th>Adopted spread<sup>a</sup></th></tr></thead><tbody><tr><td colspan="4">Initial conditions</td></tr><tr><td>Initial disk size</td><td>R<sub>0</sub></td><td>5–20 au</td><td>0.3</td></tr><tr><td>Initial disk mass</td><td>M<sub>0</sub></td><td>10<sup>-4</sup>–10<sup>-1</sup> M<sub>☉</sub></td><td>0.6</td></tr><tr><td colspan="4">MHD wind model</td></tr><tr><td>α parameter</td><td>α<sub>DW</sub></td><td>10<sup>-4</sup>–10<sup>-2</sup></td><td>0.2</td></tr><tr><td>Accretion timescale</td><td>t<sub>acc,0</sub></td><td>0.1–2 Myr</td><td><sup>b</sup></td></tr><tr><td>ω parameter</td><td>ω</td><td>0.25–0.5</td><td>0</td></tr><tr><td>Magnetic lever arm</td><td>λ</td><td>2–20</td><td>0</td></tr><tr><td>Ejection-to-accretion ratio</td><td>f<sub>M</sub></td><td>0.1–10</td><td><sup>b</sup></td></tr><tr><td colspan="4">Turbulent model</td></tr><tr><td>α parameter</td><td>α<sub>SS</sub></td><td>10<sup>-4</sup>–10<sup>-2</sup></td><td>0.2</td></tr><tr><td>Viscous timescale</td><td>t<sub>v,0</sub></td><td>0.1–2 Myr</td><td><sup>b</sup></td></tr><tr><td>Wind mass-loss rate</td><td>M&#775;<sub>PEW</sub></td><td>10<sup>-11</sup>–10<sup>-6</sup> M<sub>☉</sub> yr<sup>-1</sup></td><td>0.3</td></tr></tbody></table>

**Notes.** The free parameters are outlined in bold font.

<sup>a</sup> This corresponds to the spread of the lognormal distributions used for synthetic populations and measured in units of dex. Zero range corresponds to parameters constant across a disk population.

<sup>b</sup> Corresponds to parameters which are derived from other parameters.

#### 2.1.1. Turbulence-driven Evolution

In the turbulence-driven scenario, the evolution of the surface density profile $\Sigma(r, t)$ is governed by

$$ \frac{\partial \Sigma}{\partial t} = \frac{3}{R} \frac{\partial}{\partial R} \left\{ \frac{1}{R \Omega} \frac{\partial}{\partial R} (R^2 \alpha_{SS} \Sigma c_s^2) \right\} - \dot{\Sigma}_{\text{PEW}}(R), \quad (1) $$

where $\alpha_{SS}$ is the Shakura–Sunyaev parameter and $\dot{\Sigma}_{\text{PEW}}(R)$ is the photoevaporative wind mass-loss rate per unit surface area. In this work, we adopt the radial profile of $\dot{\Sigma}_{\text{PEW}}(R)$ computed by G. Picogna et al. (2021) for a stellar mass of $M_* = 0.5 M_☉$. The latter work utilizes hydrodynamical simulations including heating by stellar X-rays and neglecting molecule formation or shielding by a putative inner disk wind ($\lesssim 1$ au). Instead of relying on the values of mass-loss rates computed by numerical simulations, we take the total mass-loss rate as a free parameter to scale the mass-loss rate profile. This allows us to mitigate the discrepancies in terms of mass-loss rates found in different numerical simulations (see the review by I. Pascucci et al. 2023). The importance of the radial profile of $\dot{\Sigma}_{\text{PEW}}(R)$ and the comparison with the total mass-loss rate provided in the literature is discussed in Section 4. We further assume that the mass-loss rate is constant over time. This implies that we neglect the dependency of the photoevaporation rate on the disk structure, in line with the findings of J. E. Owen et al. (2012); the effect of the attenuation of the stellar radiation field by an evolving inner wind; as well as an evolving stellar radiation field and flares. We stress that the photoevaporative mass-loss rates $\dot{M}_{\text{PEW}}$ reported in this study are effective mass-loss rates formally calculated as $\dot{M}_{\text{PEW}} = \int_{0}^{+\infty} \dot{\Sigma}_{\text{PEW}}(R) 2\pi R dR$ and used to normalize the $\dot{\Sigma}_{\text{PEW}}(R)$ profile. When the local surface density drops to zero, no mass-loss rate is assumed. Therefore, for compact disks for which the surface density rapidly drops to zero in the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Example of turbulence-driven evolution for a set of parameters illustrating two different pathways of dispersal. Left: surface density profiles from the initial time (red) to dispersal (green). Right: evolution of the disk mass and accretion rate. The first solution is an example of the popular inside-out dispersal pathway with gap opening occurring for large disks and leading to relic disks. The second solution exhibits an outside-in dispersal route due to the wind extracting the mass residing in the outer disk. The two solutions share the same viscous timescale, initial mass, and wind mass-loss rate but differ in their initial disk size as indicated in the left panels. The evolution of the disk mass and accretion rate in the absence of photoevaporation is depicted in dotted lines.</figcaption></figure>

outer region where $\dot{\Sigma}_{PEW}(R)$ is nonvanishing, or for disks with gaps, the true mass-loss rate is lower than the reported $\dot{M}_{PEW}$ value.

Numerical simulations of MRI turbulence predict a weak correlation between the value of $\alpha_{SS}$ and the disk magnetization as quantified by $\beta$, the ratio between the gas and magnetic pressure (G. R. J. Lesur et al. 2023). As a consequence, $\alpha_{SS}$ could vary in time as the disk evolves but this also depends on the poorly known evolution of the magnetic field strength. For simplicity, we assume $\alpha_{SS}$ to be constant in time.

The radial profile of $\alpha_{SS}$ is taken to be constant, as often assumed in disk evolution models (e.g., G. D. Mulders et al. 2017; G. Lodato et al. 2017). Theoretical works predict steep variations in the $\alpha_{SS}$ profile depending on the ability of the MRI to be triggered. Notably within 0.1–0.5 au, the ionization of alkali metals ensures efficient MRI, possibly coupled to MHD disk winds. The bulk part of the disk and the stellar accretion rate are not expected to be affected by such a high $\alpha_{SS}$ region as gas that slowly advected from the outer to the inner regions will rapidly flow across the inner highly turbulent region down to the magnetosphere of the star. Our value of $\alpha_{SS}$ is therefore representative of the bulk part of the disk and we adopt an inner disk radius of 0.1 au.

Initially, the surface density profile follows the self-similar solution of D. Lynden-Bell & J. E. Pringle (1974) with

$$ \Sigma(R, t) = \frac{M_0}{2\pi R_0^2} \left(\frac{R}{R_0}\right)^{-1} e^{-R/R_0}, \qquad (2) $$

where $M_0$ is the initial disk mass and $R_0$ is the initial disk characteristic size. We further define the initial viscous timescale as

$$ t_{\nu,0} = \frac{R_0}{3\epsilon_0 c_{S,0} \alpha_{SS}}, \qquad (3) $$

where $c_{S,0}$ and $\epsilon_0$ are the sound speed and the disk aspect ratio at $R_0$, respectively. In this work, we assume that $T \propto R^{-1/2}$ and

$\epsilon(R = 1 \text{ au}) = 3.33 \times 10^{-2}(M_*/1M_{\odot})^{-1/2}$. This leads to the following relation:

$$ t_{\nu,0} = 0.34 \text{ Myr} \left(\frac{\alpha_{SS}}{10^{-3}}\right)^{-1} \left(\frac{R_0}{10 \text{ au}}\right) \left(\frac{M_*}{0.5 M_{\odot}}\right)^{1/2} \qquad (4) $$

The different evolution pathways under the effect of turbulence and photoevaporative wind are illustrated in Figure 1 with two solutions. The solutions share the same viscous timescale, initial disk mass, and total mass-loss rate but differ in their initial disk size. Both lead to a disk dispersal time, defined as the time when accretion stops, of about 3–4 Myr. The first example (top left panel) has already been extensively described in the literature and corresponds to a relatively large disk of $R_0 = 20$ au. As the disk evolves, it spreads radially until photoevaporation opens a gap around 2–5 au where the mass-loss rate profile $\dot{\Sigma}_{PEW}$ peaks (see Figure 7 in G. Picogna et al. 2021). Following the gap opening, the inner disk is rapidly drained and a relic disk remains. This behavior is reflected in the evolution of the accretion rate and disk mass (right panels, blue lines): after a decline of a typical timescale of $t_{\nu,0}$, the accretion rate quickly vanishes whereas the disk mass remains relatively constant, even after the accretion stops. We however stress that this so-called “relic disk problem” is only the result of our simple prescription of the photoevaporation mass-loss rate profile, which does not account for the direct irradiation of the cavity following gap opening. In the following, a disk is considered to be dispersed when the stellar accretion rate drops below a threshold of $10^{-12} M_{\odot} \text{ yr}^{-1}$.

In the second solution, the initially more compact disk ($R_0 = 10$ au) does not spread despite the radial transport of angular momentum. This surprising feature is because the bulk of the disk and its outermost region reside where the wind mass-loss rate is significant. The consequence is that the wind does not open a gap but truncates the disk. This results in a drastically different evolution: the outer disk is eroded by the wind and the inner disk is drained due to turbulent accretion.
Therefore, the disk mass, accretion rate, and size jointly decrease before full disk dispersal. This evolution is somewhat similar to that of a large disk exposed to external UV radiation since external photoevaporation removes mass from the outer disk (e.g., C. J. Clarke 2007).

### 2.1.2. Magnetohydrodynamic Wind-driven Evolution

Following the formalism of B. Tabone et al. (2022a) the evolution of the surface density under the effect of an MHD disk wind is described as

$$ \frac{\partial \Sigma}{\partial t} = \frac{3}{2R} \frac{\partial}{\partial R} \left[ \frac{\alpha_{\text{DW}} \Sigma c_s^2}{\Omega} \right] + \frac{3}{R} \frac{\partial}{\partial R} \left\{ \frac{1}{r\Omega} \frac{\partial}{\partial R} (r^2 \alpha_{\text{SS}} \Sigma c_s^2) \right\} - \dot{\Sigma}_{\text{MHD-DW}}, \quad (5) $$

where $\alpha_{\text{DW}}$ parameterizes the wind torque and $\dot{\Sigma}_{\text{MHD-DW}}$ is the mass-loss rate of the wind. In contrast with photoevaporative winds, the mass-loss rate of the MHD disk wind depends on the extraction of angular momentum as

$$ \dot{\Sigma}_{\text{MHD-DW}} = \frac{3\alpha_{\text{DW}}\Sigma c_s^2}{4(\lambda - 1)R^2\Omega}, \quad (6) $$

where $\lambda$, the magnetic lever arm parameter (R. D. Blandford & D. G. Payne 1982), quantifies the efficiency of the wind to carry out angular momentum. The lower the $\lambda$, the more mass is required to be launched to drive disk accretion. A considerable advantage of this parameter is that $\lambda$ can be estimated from wind observations (J. Ferreira et al. 2006; B. Tabone et al. 2020). Observations typically find $\lambda \simeq 2-6$ (B. Tabone et al. 2017; A. de Valon et al. 2020; A. S. Booth et al. 2021; P. Nazari et al. 2024) and in this work we explore a wide range of $\lambda$ values (Table 1). In Equation (5), $\alpha_{\text{SS}}$ quantifies radial transport by either turbulence or by a laminar torque emerging from the radial component of the magnetic field, which is usually negligible compared to vertical torque.

In this work, we adopt the analytical solutions of B. Tabone et al. (2022a) in the absence of radial transport of angular momentum ($\alpha_{\text{SS}} = 0$). The analytical solutions further assume a constant $\lambda$ value in space and time and a constant $\alpha_{\text{DW}}$ across the disk. The surface density profile is given by

$$ \Sigma(R, t) = \Sigma_c(t) (R/R_0)^{-1+\xi} e^{-R/R_0}, \quad (7) $$

where $\Sigma_c(t)$ is the characteristic surface density, $\xi = 1/[2(\lambda - 1)]$ is the mass ejection index, and $R_0$ is the characteristic disk radius. A considerable source of uncertainty in MHD wind-driven disk evolution is the secular evolution of the magnetic field strength since $\alpha_{\text{DW}} \propto B_z B_\phi / \Sigma$. Here, we use the “$\Sigma_c$ dependent wind torque” solutions of B. Tabone et al. (2022a) where the evolution of the wind torque is parameterized as $\alpha_{\text{DW}} = \alpha_{\text{DW}}(0)(\Sigma_c(t)/\Sigma_c(0))^{-\omega}$, where $\alpha_{\text{DW}}(0)$ is the initial $\alpha_{\text{DW}}$ and $\omega$ is a free parameter. The case $\omega = 1$ describes a constant magnetic field strength (assuming a constant $B_z/B_\phi$ ratio). This leads to the following evolution of the disk mass

and accretion rate:

$$ \begin{gathered} M_D(t) = M_0 \left( 1 - \frac{\omega}{2t_{\text{acc},0}} t \right)^{1/\omega}, \\ \dot{M}_*(t) = \frac{M_0}{2t_{\text{acc},0}(1+f_M)} \left( 1 - \frac{\omega}{2t_{\text{acc},0}} t \right)^{-1+1/\omega}, \quad (8) \end{gathered} $$

where

$$ f_M \equiv \dot{M}_W / \dot{M}_* = (R_0 / R_{\text{in}})^\xi - 1, \quad (9) $$

is the mass ejection-to-accretion ratio. $t_{\text{acc},0}$ is a generalization of the viscous timescale (see Equation (3)). Following the same assumption about the disk thermal structure as for the turbulent case, the accretion timescale is

$$ t_{\text{acc},0} = 0.34 \text{Myr} \left( \frac{\alpha_{\text{DW}}}{10^{-3}} \right)^{-1} \left( \frac{R_0}{10 \text{ au}} \right) \left( \frac{M_*}{0.5 M_{\odot}} \right)^{1/2} . \quad (10) $$

The essential features of the solution are illustrated in Figure 2 with two values of $\omega$. In both cases, the surface density profile preserves its shape while decreasing in absolute value. The angular momentum is indeed transported vertically without the need for radial expansion. We however stress that this is a property of the solution related to its initial conditions. For sharp cuts in the outer disk, the disk can shrink since wind-driven accretion acts as an advection of gas. The latter situation is however not realistic since any turbulent or laminar radial stress in the early phase leads to a smoothly declining outer surface density profile. We note that the alternative model of T. K. Suzuki et al. (2016) predicts the formation of a cavity, reminiscent of transition disks. Whereas cavities sustained by MHD disk winds are a stable configuration (É. Martel & G. Lesur 2022), gap opening has never been characterized in global numerical simulations. Pending simulations on secular timescales, we adopt here disk evolution solutions that do not introduce artificial features in the surface density profile. The second property of MHD wind-driven evolution is the dispersal of the disk after a finite time due to rapid disk draining. Therefore, our MHD wind-driven solution accounts for both accretion and dispersal. More specifically, the disk dispersal time depends on $t_{\text{acc},0}$ and $\omega$ as

$$ t_{\text{disp}} = 2t_{\text{acc},0}/\omega. \quad (11) $$

In both solutions, we choose $t_{\text{acc},0}$ so that the two cases share the same dispersal time (see Figure 2, right panels). This highlights the essential effect of the $\omega$ parameter, which describes the evolution of the magnetic field strength. Higher values of $\omega$ lead to a rather constant disk mass and a high stellar accretion rate before dispersal. Therefore, for a given dispersal time, a disk with a higher $\omega$ value spends more time in a high mass and high accretion rate phase. This property is essential for the evolution of synthetic populations. Finally, we note that from Equation (8) the value of $\lambda$ does not affect the disk mass but the accretion rate by a factor $f_M + 1 = (R_0/R_{\text{in}})^{1/2}[\lambda-1]$.

### 2.2. Population Synthesis with Diskpop

Observations provide us with the properties of disk populations of different ages and not with the evolution of individual disks. To simulate the evolution of disk populations,
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Example of MHD wind-driven evolution for a set of parameters illustrating the impact of time evolution of the disk magnetization as parameterized by $\omega$. The left panels show the surface density profile from the initial time (red) to dispersal (green). The right panels show the evolution of the disk mass and accretion rate. The two solutions share the same value of $\lambda$ and initial mass and size but differ by the value of $\omega$. The accretion timescale is chosen to get the same dispersal time according to Equation (11). In practice, the disk is dispersed slightly earlier because either the accretion rate or the disk mass dropped to their threshold disk dispersal values.</figcaption></figure>

we use the *Diskpop* Python package (A. Somigliana et al. 2024). As a first step, the code randomly picks the values of the parameters for a sample of disks following a user-defined probability distribution and possible correlations between these parameters. The gas evolution of each disk is then computed by solving numerically Equation (1) or (5), or using analytical solutions with the option of including dust evolution. The output of *Diskpop* is then analyzed using the *Popcorn* package (A. Somigliana et al. 2024).

Throughout this paper, the synthetic populations are made of 400 disks by default and we consider only gas evolution. When building a synthetic disk population, the four free disk parameters outlined in bold font in Table 1 and detailed above translate into four probability distributions. We assume that the distributions of the four disk parameters are not correlated and follow a lognormal distribution. Therefore, our population model is controlled by eight free population parameters: a median and a spread for each disk parameter. The MHD disk-wind model has a fifth free parameter, $\omega$, for which only two values are explored. Following B. Tabone et al. (2022b), a too low value of $\omega$ produces disks with no detectable accretion signatures but significant mass, which would lead to a strong difference between the dispersal time inferred from accretion signatures and IR excess. Values of $\omega > 1$ are also unrealistic since they would correspond to disks with increasing magnetic field strength (B. Tabone et al. 2022a).

Fitting the properties of disk populations using a population synthesis approach can appear as an underconstrained problem. Here, we propose a rationalized approach to map the parameter space and alleviate the degeneracies in a step-by-step process.

### 2.2.1. Observational Data

The main observational data to be confronted with population synthesis are shown in Figure 3. Our synthetic populations are designed to describe the evolution of isolated

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Summary of the observational data to be reproduced by our population synthesis model. Top left: disk fraction probed by the IR excess in several star-forming regions. The disk fractions compiled by D. Fedele et al. (2010) are rescaled by a factor of 1.2 to account for the short-lived disks around binaries. Top right: stellar accretion rates from the nearly complete surveys of the Lupus and Upper Sco star-forming regions (J. M. Alcalá et al. 2017; M. Fang et al. 2023). The boxes represent the median accretion rates and the dots and triangles are the measurement and upper limits of the individual sources. Bottom: AGE-PRO estimates of the median CO gas size (L. Trapman et al. 2025b) and gas mass (L. Trapman et al. 2025a), along with the measurements for individual sources. We note that the AGE-PRO Upper Sco sample is somewhat younger than the bulk region over which the accretion rates are measured.</figcaption></figure>

disks. Therefore, we consider the $t=0$ age as the end of the Class I phase and take the Oph AGE-PRO sample as representative of this initial condition. As discussed in K. Zhang et al. (2025), there is a significant spread in the isochronal age of each source for the Lupus and Upper Sco AGE-PRO sample with overlap between the two regions. Considering the considerable systematic uncertainties in the
age determination, we adopt a median age of 2 Myr for the Lupus and 4 Myr for the Upper Sco AGE-PRO sample with a typical range of 2 and 4 Myr, respectively. We note that the sources selected by AGE-PRO in Upper Sco are part of the youngest sources of this star-forming region. Therefore, our adopted median age of the AGE-PRO Upper Sco sample is shorter than the median age of the entire cluster, which is taken to be 7 Myr. This results in an age difference between the median disk mass and size, estimated from the AGE-PRO subsample (Figure 3, bottom panels), and the median accretion rate evaluated from a less biased sample (top right panel).

The disk fraction (top left panel) stems from the compilation of D. Fedele et al. (2010) of IR excess, which is consistent with other works (Á. Ribas et al. 2014) and representative of the stellar mass targeted by AGE-PRO (Á. Ribas et al. 2015). In our models, we neglect the effect of multiplicity, which significantly reduces the disk lifetime for close-in binaries (≤40 au; A. L. Kraus et al. 2012). We corrected the disk fraction of D. Fedele et al. (2010) by applying a multiplicative factor of 1.2 to match the disk fraction around single or wide binaries determined by A. L. Kraus et al. (2012). The disk fraction is interpreted as the cumulative probability distribution of the disk lifetime. By assuming that disk lifetimes follow a lognormal distribution, we fit the disk fraction by an error function (see black line Figure 3, top left). We estimate a spread of about 0.3 dex and a median value of the disk lifetime of 3 Myr.

The accretion rates shown in Figure 3 stem from the compilation of J. M. Alcalá et al. (2014, 2017) and C. F. Manara et al. (2023) for Lupus and M. Fang et al. (2023) for the Upper Sco star-forming regions. We selected the stars that are within the bin of AGE-PRO stellar masses ($M_{*} = 0.3-0.7 M_{\odot}$). We stress that the sample of M. Fang et al. (2023) reports upper limits for 50% of the sources. Therefore, the median value of the accretion rate should be considered as an upper limit. Compared to the pioneering survey of C. F. Manara et al. (2020) the sample of M. Fang et al. (2023) is larger (12 versus 29 sources within the AGE-PRO stellar mass range) thought using only the Hα line flux as a proxy of the accretion rate. Still, the two studies lead to similar median accretion rates for Upper Sco. Overall, we find that the median accretion rate decreases by a factor of at least 3 from Lupus to Upper Sco with a large spread of 0.8 dex for each region.

The distribution of disk size and mass obtained by the AGE-PRO program is also represented in Figure 3 as derived in L. Trapman et al. (2025b) and L. Trapman et al. (2025a), respectively. The observed CO gas sizes, denoted as $R_{CO,90\%}$ or simply $R_{CO}$, correspond to the radius that encompasses 90% of the total CO flux and are evaluated by correcting for the beam of the observations. A major caveat when analyzing disk size is the conversion between the actual disk size, i.e., the radius encompassing the majority of disk mass, and the observed CO disk size (L. Trapman et al. 2023). A possible strategy to extract the predicted value of $R_{CO,90\%}$ in the population model is to run radiative transfer calculations using the surface density calculated by the 1D disk evolution model and assuming either standard CO abundances or CO abundances computed with thermochemical models. This method is however computationally prohibitive and in this work, we use a simple analytical formula derived from L. Trapman et al. (2023), who show that the $R_{CO,90\%}$ size is given by the radius where the total gas column density drops

below a given critical column density of hydrogen denoted as $N_{gas}(R_{CO,90\%})$. However, L. Trapman et al. (2023) do not consider the depletion of gas-phase carbon, which can be substantial in Class II disks (C. Favre et al. 2013). Following C. Toci et al. (2023), we increased the threshold column density by the carbon depletion fraction as

$$N_{\text{gas}}(R_{\text{CO},90\%}) = 3.7 \times 10^{21} \delta_{\text{C}}^{-1} (M_{\text{D}}/M_{\odot})^{0.34} \text{cm}^{-2}. \quad (12)$$

Throughout this work we adopt $\delta_{\text{C}} = 0.2$, in line with the underabundance of CO found in L. Trapman et al. (2025a) for the Lupus and Upper Sco AGE-PRO samples. Therefore, by finding the radius where the column density reaches $N_{gas}(R_{CO,90\%})$, we determine for each simulated disk the $R_{CO,90\%}$ radius.

The disk masses in Lupus and Upper Sco are calculated by matching the $N_{2}H^{+}$ and CO line fluxes, and CO gas size (L. Trapman et al. 2025a). For the younger disk population of Oph, only CO isotopologue lines are used. Overall, we find a drop in disk mass from Oph to Lupus by about an order of magnitude. From Lupus to Upper Sco, the median CO size and mass do not change significantly.

### 2.2.2. Fitting Procedure

Our goal is to reproduce the disk fraction and the median values of the accretion rate, disk mass, and disk size. The spreads of the parameters are therefore of less importance and we keep the four spreads fixed to meaningful values provided in Table 1. As explained below (see Section 3 and 4), these values adequately fit the dropping rate of the disk fraction, and the spread in the observed disk masses and disk sizes. The main challenge is then to constrain the four free population parameters. These four free parameters are constrained by sequentially reproducing (1) the disk fraction (through the disk dispersal timescale), (2) the median accretion rate in Lupus, (3) the median CO gas size inferred in Lupus, and (4) the median gas mass estimated in Lupus. The workflow is summarized in Figure 4 for the two scenarios. The middle-aged Lupus population is used to constrain the parameters before a detailed comparison with the old Upper Sco and the young Oph populations. The choice of the Lupus population as a reference point is justified by the systematic measurement of accretion rates for all the Lupus AGE-PRO sources, the relatively narrow range in stellar ages compared with Upper Sco, and the measurement of the CO disk size, which is challenging in Oph due to the presence of envelopes (D. A. Ruiz-Rodriguez et al. 2025).

Because disk dispersal is not controlled by the same parameters (see Section 2.1), the parameter adjusted for each stage differs in the two scenarios. In the first stage (fit of the disk fraction), it is the median mass-loss rate $\langle M_{\text{PEW}} \rangle$ that is adjusted for the turbulence-driven scenario, while for the MHD disk-wind scenario it is the accretion timescale $\langle t_{\text{acc},0} \rangle$. For the second stage (fit of the median accretion rate in Lupus), the initial disk mass is adjusted in both scenarios, since the accretion rates are proportional to $M_{0}$. The third stage (fit of the median CO gas size in Lupus) constrains the initial disk size $R_{0}$, and the last stage constrains the remaining parameters, namely the viscous timescale for the turbulence-driven case and the wind-to-accretion mass ratio in the MHD disk-wind case. These four stages allow us to confront the predictions of the population synthesis models to the full
<figure><img src="image_6.png" /><figcaption>Figure 4. Schematic view of the fitting stages for the turbulent (left) and the MHD disk-wind (right) models. This step-by-step approach allows us to build synthetic populations consistent with the observations in a rationalized approach.</figcaption></figure>

AGE-PRO data with very limited degeneracy in the parameters.

## 3. Results

### 3.1. Magnetohydrodynamic Wind-driven Evolution

#### 3.1.1. Disk Dispersal

In the analytical solution of B. Tabone et al. (2022a), the theoretical disk dispersal time, as defined by the time when disk mass drops to zero, depends only on the accretion timescale $t_{acc,0}$ and the $\omega$ parameter, which quantifies the evolution of the magnetic field strength. For $\omega > 0$ and a stellar mass of $M_* = 0.5 M_\odot$, an individual disk disperses at a time

$$t_{\text{disp}} = 2t_{\text{acc,0}}/\omega = 0.68 \text{ Myr} \left(\frac{\alpha_{\text{DW}}}{10^{-3}}\right)^{-1} \left(\frac{R_0}{10 \text{ au}}\right) \omega^{-1}. \quad (13)$$

In this work, our definition of disk dispersal time as the time when the stellar accretion rate drops below $10^{-12} M_\odot \text{ yr}^{-1}$, leads to very similar disk dispersal time for $\omega \gtrsim 0.2$.

Therefore, one can associate for each individual disk, a dispersal time $t_{disp}$. Because $\alpha_{DW}$ and $R_0$ follow lognormal probability distributions, the disk dispersal time will also follow a lognormal distribution with a median value and a spread of

$$\langle t_{\text{disp}} \rangle = 0.68 \text{ Myr} \left( \frac{\langle \alpha_{\text{DW}} \rangle}{10^{-3}} \right)^{-1} \left( \frac{\langle R_0 \rangle}{10 \text{ au}} \right) \omega^{-1}, \\ \sigma_{t_{\text{disp}}}^2 = \sigma_{\alpha_{\text{DW}}}^2 + \sigma_{R_0}^2, \quad (14)$$

where $\sigma_X$ is the spread of the quantity X measured in units of dex.

In the first stage of the fitting process, we consider $R_0$ as a free parameter, which will be constrained from the observed CO size in stage 3 (see Figure 4). Therefore, to fit a median disk dispersal time of 3 Myr, one needs to set the median value of $\alpha_{DW}$ for each value of $R_0$ to

$$\langle \alpha_{\text{DW}} \rangle = 2.3 \times 10^{-4} \left( \frac{\langle R_0 \rangle}{10 \text{ au}} \right) \omega^{-1}. \quad (15)$$

Equation (14) also justifies the adopted spread in $R_0$ and $\alpha_{DW}$ (see values in Table 1). We indeed choose the spreads such that $\sqrt{\sigma_{\alpha_{DW}}^2 + \sigma_{t_{R0}}^2} = 0.36$, a value close to the spread in the disk lifetimes inferred from the decline of the disk fraction.

#### 3.1.2. Accretion Rates

As discussed in B. Tabone et al. (2022b), once the distribution of $t_{acc,0}$ is set, some essential properties of accretion are set since $t_{acc,0}$ also controls the ratio between the accretion rate and the disk mass. In particular, Equation (8) shows that the absolute value of the accretion rate scales with the initial disk mass. Therefore, in the second stage of our fitting procedure, we adjust the median disk mass $\langle M_0 \rangle$ for each value of $\omega$, $f_M$, and $\langle R_0 \rangle$ to match a median accretion rate of $\dot{M}_* = 10^{-9} M_\odot \text{ yr}^{-1}$ measured in Lupus. Practically, we start with a first guess of the initial disk mass to compute the evolution of a population and measure the median accretion rate at an age of 2 Myr. The initial median disk mass is subsequently increased (decreased) if the median accretion is higher (lower) than $10^{-9} M_\odot \text{ yr}^{-1}$. A tolerance of 10% on the median accretion rate is adopted.

The resulting values of $\langle M_0 \rangle$ as a function of the remaining free parameters ($\langle R_0 \rangle$, $\omega$, and $\langle f_M \rangle$) are shown in Figure 5. When plotted as a function of the ejection-to-accretion mass ratio $f_M = (R_0/R_{in})^{1/[2(\lambda-1)]} - 1$, the results do not depend on the assumed initial disk sizes. This is because the time evolution of the accretion rate is primarily driven by the accretion timescale $t_{acc,0}$ and not directly the disk size.

From Equation (8), increasing the ejection-to-accretion mass ratio $f_M$ reduces all the accretion rates by a factor of $1/(1+f_M)$. Therefore, the value of $\langle M_0 \rangle$ required to reproduce the accretion rate in Lupus increases with $f_M$ as $1+f_M$ (see the black curve in Figure 5), and the impact of the mass-loss rate is significant only for $f_M \gtrsim 1$. Finally, lower values of $\omega$ require higher disk masses to fit the accretion rates. The effect of $\omega$ is more subtle to grasp but can be understood from an inspection of Figure 2. For the same disk dispersal time, a lower value of $\omega$ leads to an accretion rate with a steeper drop before dispersal. This implies that lower values of $\omega$ lead to lower accretion rates at 2 Myr. As a consequence, the initial disk
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Constraints on the initial median disk mass $\langle M_0 \rangle$ obtained by reproducing the median accretion rate of Lupus ($\dot{M}_* \simeq 10^{-9} M_\odot \text{ yr}^{-1}$) for the MHD disk-wind case. The gray area indicates the median disk mass derived by AGE-PRO in Oph. Each dot corresponds to a synthetic population with a medial $\alpha_{DW}$ given by Equation (15) matching the observed disk dispersal.</figcaption></figure>

mass $\langle M_0 \rangle$ required to fit the observed accretion rate in Lupus needs to be increased when decreasing the $\omega$ value.

Overall, we find that initial median disk masses of at least $\langle M_0 \rangle \simeq 8 \times 10^{-3} M_\odot$ are required to reproduce the median accretion rate measured in Lupus. The median gas mass of $\langle M_D \rangle \simeq 7_{-2}^{+4} \times 10^{-3} M_\odot$ obtained by AGE-PRO in Oph (D. A. Ruiz-Rodriguez et al. 2025) favors low ejection-to-accretion ratios $f_M \lesssim 1$ and relatively high values of $\omega$. Interestingly, large mass ejection ratios $f_M$ are excluded since this would require disks that are initially unstable against gravitational instability.

### 3.1.3. Constraints from the AGE-PRO Results for Lupus

At this stage, we are in a position to compare the results of AGE-PRO to synthetic populations which are consistent with disk dispersal and the accretion rate in Lupus (stage 3 and 4 in Figure 4). The initial disk mass and accretion timescale are already constrained and only the initial disk size, the $f_M$ parameter, and $\omega$ remain free parameters. In Figure 6, we show the predicted median CO disk size and the median disk mass after 2 Myr, corresponding to the age of Lupus. As expected, a population of initially larger disks has a larger $\langle R_{CO} \rangle$ size after 2 Myr (Figure 6 top panel). It is immediately apparent that only initially compact disks can fit the observed CO disk sizes in Lupus with a median disk radius of $\langle R_0 \rangle \simeq 10$ au. We stress that this does not imply that the characteristic size at the age of Lupus is that small since compact disks have shorter accretion timescales and therefore are dispersed first (see Section 4).

The median disk mass $\langle M_D \rangle$ at 2 Myr increases with the mass ejection-to-accretion ratio $f_M$. This is simply because higher $f_M$ values correspond to a higher initial mass $\langle M_0 \rangle$—a requirement to fit the median accretion rate in Lupus (see Figure 5 and Section 3.1.2). Interestingly, we find that lower values of $\omega$ produce lower disk masses even though the initial disk mass is larger. This is an essential effect of $\omega$: for a single disk, the mass decreases more steeply before disk dispersal for low $\omega$ values (see example Figure 2). Therefore, a lower value of $\omega$ results in a lower disk mass at the age of Lupus.

When compared with the median disk mass estimated by AGE-PRO in Lupus (Figure 6, bottom panel), half of the synthetic populations that fit the disk dispersal, accretion rate,

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Predicted median CO gas size and disk mass at the age of Lupus for MHD wind-driven synthetic populations matching the disk fraction (see Equation (15)) and the Lupus median accretion rate (see Figure 5). The gray area shows the median CO disk size and disk mass as measured by AGE-PRO.</figcaption></figure>

and disk size are excluded. Most of the predicted median disk masses lie indeed above the median disk mass of Lupus. For $\omega = 0.25$, values of $f_M$ smaller than one reproduce well the low disk masses measured in Lupus even though $\omega = 0.5$, and $f_M \lesssim 1$ remains a good fit considering the large uncertainty on the estimated median disk mass.

### 3.1.4. Best-fit Magnetohydrodynamic Disk-wind Model

Putting all the constraints stemming from the properties of the Lupus population and disk fraction together, we find best-fit models with an $\omega = 0.25-0.5$, low mass ejection-to-accretion ratios $f_M \lesssim 1$, initially compact disks $\langle R_0 \rangle = 10$ au, and $\alpha_{DW} \simeq 5 \times 10^{-4}-10^{-3}$. Reproducing the essential features of the Lupus population constitutes the first success of the MHD disk-wind models. This also considerably reduces the parameter space before comparing the time evolution of the median accretion rate, mass, and CO gas size. Here, we retained two best-fit models from stage 4, which correspond to $\omega = 0.25$ (steep decline) and $\omega = 0.5$ (shallow decline).

In Figure 7, we present the time evolution of the median disk mass $\langle M_D \rangle$, accretion rate $\langle \dot{M}_* \rangle$, and CO gas size $\langle R_{CO} \rangle$ and compare them to the available observational constraints. By construction, the two synthetic populations reproduce the Lupus population at 2 Myr. We recall that higher $\omega$ values (right panels) correspond to more shallow evolution of the medians. Of the two models, the synthetic population with a higher value of $\omega = 0.5$ reproduces better the evolution from Oph to Upper Sco. Overall, the match between the observations and the best-fit model is excellent. This is particularly clear with the disk mass. AGE-PRO finds a decrease in the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Best-fit MHD disk-wind model for $\omega = 0.25$ (left) and $\omega = 0.5$ (right) compared with the observations. The top panels show the total disk gas mass, the middle panels the accretion rates, and the bottom panels the CO disk size. Individual sources are represented by dots and the corresponding median with uncertainties by colored boxes. AGE-PRO measurements are from L. Trapman et al. (2025a, top panels) and L. Trapman et al. (2025b, bottom panels) whereas the accretion rates were compiled by C. F. Manara et al. (2023) for Lupus and estimated by M. Fang et al. (2023) for Upper Sco. The gray shaded areas depict the first and third quartiles of synthetic subpopulations made of 10 disks. This allows us to depict the uncertainties related to the limited size of the samples. The synthetic populations are made of 2000 disks. We stress that the medians include the survivor bias since they are calculated over the disk-bearing sources of the synthetic population.</figcaption></figure>

median disk mass of about 10 from Oph to Lupus, which is well reproduced by $\omega = 0.5$ and overestimated by $\omega = 0.25$. The importance of having estimates for a young star-forming region was also highlighted in Figure 5 where we present the initial disk mass required to fit disk fraction and accretion rate: only $\omega = 0.5$ and low $f_M$ values led to initial disk masses consistent with Oph. An $\omega = 0.5$ value is also consistent with the relatively constant median disk mass found from Lupus to Upper Sco (Figure 7).

Both models also reproduce the shallow evolution of the CO disk size with a small increase in $\langle R_{CO} \rangle$ from Lupus to Upper Sco. In the models, the increase in $\langle R_{CO} \rangle$ is exclusively a survivorship bias. In fact, for an individual disk, $R_{CO}$ is rather slowly decreasing as an effect of the decrease in the disk mass. However, small disks are dispersed first due to their smaller $t_{acc,0}$, driving an increase in the median disk size. The result on the median accretion rate is similar to the disk mass: an $\omega = 0.5$ value does a slightly better job at reproducing the shallow evolution from Lupus to Upper Sco.

## 3.2. Turbulence-driven Evolution with Internal Photoevaporation

### 3.2.1. Disk Dispersal

In the case of turbulent disks, the dispersal time for an individual disk depends not only on the accretion timescale (often also referred to as the viscous timescale), but also on the

three other disk parameters: the wind mass-loss rate $\dot{M}_{PEW}$, the initial disk mass $M_0$ and, to a lesser extent, the disk size $R_0$. This means that fitting disk dispersal does not provide a constraint on only one parameter, as in the case of MHD wind-driven accretion, but a joint constraint on the four parameters. In the first stage (see workflow in Figure 4), we determine the median photoevaporative wind mass-loss rate $\langle \dot{M}_{PEW} \rangle$ required to fit disk fraction for all possible values of the viscous timescale $\langle t_{v,0} \rangle$ (or, equivalently $\langle \alpha_{SS} \rangle$), initial disk size $\langle R_0 \rangle$, and mass $\langle M_0 \rangle$.

Quantitatively, we find that for $t_{disp} > t_{v,0}$ the dispersal time of an individual disk is roughly a power law of the disk parameters:

$$t_{disp} \approx 2.9 \text{ Myr} \left( \frac{t_{v,0}}{1 \text{ Myr}} \right)^{1/3} \left( \frac{M_0}{10^{-1} M_{\odot}} \right)^{2/3} \times \left( \frac{\dot{M}_W}{1 \times 10^{-8} M_{\odot} \text{yr}^{-1}} \right)^{-2/3} \quad (16)$$

This equation is motivated by the idea that a disk disperses when its accretion rate as predicted by the pure viscous evolution drops below the wind mass-loss rate (C. J. Clarke et al. 2001; R. Alexander 2008). This reasoning provides the power-law indices in Equation (16). We also note that this simple equation predicts that to first order, the dispersal time is sensitive to the viscous timescale, which is proportional to the ratio $R_0/\alpha_{SS}$ (Equation (3)).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Illustration of the method used to constrain the median mass-loss rate $\langle \dot{M}_{\text{PEW}} \rangle$ from the disk fraction, assuming fixed values of $R_0$, $\langle M_0 \rangle$, and $\langle t_{\nu,0} \rangle$. The predicted disk fraction as defined by nonaccreting disks is shown as a solid line for various values of the median mass-loss rate. The other population parameters are indicated in the top left. The fitted values of $\langle \dot{M}_{\text{PEW}} \rangle$ for $R_0 = 10$ au and wide ranges of $\langle M_0 \rangle$ and $\langle t_{\nu,0} \rangle$ are provided in Figure 9. The data compiled by D. Fedele et al. (2010) are in blue.</figcaption></figure>

The idea behind the fitting procedure is illustrated in Figure 8: for given values of the median $\langle R_0 \rangle$, $\langle M_0 \rangle$, and $\langle t_{\nu,0} \rangle$, the decline of disk fraction with the age of the population depends on the value of $\langle \dot{M}_{\text{PEW}} \rangle$. We therefore assume fixed values of $\langle R_0 \rangle$, $\langle M_0 \rangle$, and $\langle t_{\nu,0} \rangle$, and run a synthetic disk population using a guessed value for $\langle \dot{M}_{\text{PEW}} \rangle$. From Equation (16), the guessed value to obtain a disk dispersal time of $t_{\text{disp}} = 3$ Myr is

$$ \begin{aligned} \langle \dot{M}_{\text{PEW}} \rangle &\approx 1 \times 10^{-8} M_{\odot} \text{yr}^{-1} \left( \frac{\langle t_{\nu,0} \rangle}{1 \text{Myr}} \right)^{1/2} \left( \frac{\langle M_0 \rangle}{0.1 M_{\odot}} \right) \\ &\approx 5.8 \times 10^{-9} M_{\odot} \text{yr}^{-1} \left( \frac{\alpha_{\text{SS}}}{10^{-3}} \right)^{-1/2} \left( \frac{R_0}{10 \text{au}} \right)^{1/2} \\ &\quad \times \left( \frac{\langle M_0 \rangle}{0.1 M_{\odot}} \right). \end{aligned} \tag{17} $$

If the disk fraction at the age of 3 Myr is larger (smaller) than 50%, the mass-loss rate is reduced (increased) in proportion, and a new population is calculated. Convergence is typically reached after 2–4 iterations with a tolerance on the disk fraction of 6%.

The inferred median mass-loss rate $\langle \dot{M}_{\text{PEW}} \rangle$ as a function of the assumed $\langle t_{\nu,0} \rangle$ and $\langle M_0 \rangle$ is shown in Figure 9(a) for $\langle R_0 \rangle = 10$ au. The global evolution of $\langle \dot{M}_{\text{PEW}} \rangle$ with $\langle t_{\nu,0} \rangle$ is similar for different values of $\langle R_0 \rangle$ as shown in the Appendix (see Figure 15). As expected from Equation (17), the inferred mass-loss rate $\langle \dot{M}_{\text{PEW}} \rangle$ is larger for a higher disk mass $\langle M_0 \rangle$, and a longer viscous timescale $\langle t_{\nu,0} \rangle$. The mass-loss rate is also relatively insensitive to the initial disk size $\langle R_0 \rangle$ solely but on the viscous timescale $\langle t_{\nu,0} \rangle$, which is the ratio between the initial characteristic size $R_0$ and $\alpha_{\text{SS}}$. Interestingly, the increase in the mass-loss rate $\langle \dot{M}_{\text{PEW}} \rangle$ with $\langle t_{\nu,0} \rangle$ is steeper above $\langle t_{\nu,0} \rangle \approx 2$ Myr. This is because, for very low values of $\alpha_{\text{SS}}$ (i.e., long $t_{\nu,0}$), the timescale to drain the innermost regions of the disk becomes of the same order of magnitude as the dispersal timescale. Therefore, the outer disk, where the photoevaporative wind is emitted from, needs to be very quickly dispersed, explaining the steep increase in the wind mass-loss rate$^{23}$ with $t_{\nu,0}$. Even longer viscous timescales are

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Constraints on $\langle \dot{M}_{\text{PEW}} \rangle$ obtained by fitting a median disk dispersal time of $t_{\text{disp}} = 3$ Myr and the resulting median accretion rate at the initial time and after 2 Myr for a median initial disk size of $\langle R_0 \rangle = 10$ au. The method used to find the value of $\langle \dot{M}_{\text{PEW}} \rangle$ that fits disk fraction for each value of $\langle M_0 \rangle$, $\langle t_{\nu,0} \rangle$, and $\langle R_0 \rangle$ is illustrated in Figure 8. For $\langle R_0 \rangle = 10$ au, viscous timescales longer than 3 Myr are not able to reproduce the disk fraction and median accretion rates simultaneously at 2 Myr. Results of the fit for other values of the initial median disk size $\langle R_0 \rangle = 5, 20$, and 50 au are provided in the Appendix (see Figure 15). The solutions with unrealistically large mass-loss rates of $\langle \dot{M}_{\text{PEW}} \rangle \gtrsim 10^{-6} M_{\odot} \text{yr}^{-1}$ are not displayed.</figcaption></figure>

excluded since they lead to too long dispersal times due to the long lifetimes of the inner regions that are not affected by photoevaporation.

### 3.2.2. Accretion Rates

The constraints on the median mass-loss rate obtained from the disk fraction are converted into median initial accretion rates in Figure 9(b) for a reference initial disk size of $\langle R_0 \rangle = 10$ au. The initial $\langle \dot{M}_* \rangle$ scales as $\propto t_{\nu,0}^{-1}$ since $\langle \dot{M}_* \rangle = \langle M_0/2t_{\nu,0} \rangle$ and since the mass-loss rate does not influence the initial accretion rate. In Figure 9(c), we further plot the median accretion rate of each synthetic population after 2 Myr. Due to the decline of the accretion rate of each disk, the median accretion rate is shown to drop between $t=0$ and 2 Myr. The drop is more pronounced for short viscous timescales due to the faster disk evolution. We also note that photoevaporation tends to disperse the lightest disks, which have higher accretion rates, limiting the drop in the median

<footnote>$^{23}$ We recall here that the parameter $\dot{M}_{\text{PEW}}$ is an effective mass-loss rate used to scale the $\dot{\Sigma}_{\text{PEW}}$ profile. When a disk is compact, which is the case for long $t_{\nu,0}$ and small $R_0$, the true mass-loss rate is smaller.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Summary of the constraints on turbulence-driven disk evolution obtained from the disk dispersal time and median accretion rate measured in Lupus (namely, after following Stage 2 of the analysis scheme). The fitted values of the $\langle M_0 \rangle$ (panel (a)) and $\langle \dot{M}_{PEW} \rangle$ (panel (b)) parameters are plotted as a function of $\langle t_{\nu,0} \rangle$ for various values of $\langle R_0 \rangle$. The calculated median $\langle R_0 \rangle$ and $\langle M_D \rangle$ at the estimated age of Lupus (2 Myr) are shown in panels (c) and (d). Median CO gas size and disk mass estimated by AGE-PRO are represented by a gray-shaded area. None of the synthetic populations can simultaneously reproduce the low disk mass and the disk size. Focusing on $\langle R_{CO} \rangle$, we selected three best-fit models outlined by black circles and boxes.</figcaption></figure>

accretion rates. Overall, this leads to the counterintuitive result that the viscous timescale has little impact on the median accretion rate at 2 Myr, except for very long viscous timescales.

Figure 9(c) illustrates our method to constrain the initial disk mass from the median accretion rate estimated for Lupus using the example of an initial disk size of $\langle R_0 \rangle = 10$ au. The median accretion rate of the Lupus disk population is depicted by a gray area. From this figure, it is evident that for $\langle R_0 \rangle = 10$ au and $\langle t_{\nu,0} \rangle \simeq 0.2–1.5$ Myr, an initial median disk mass of about $1 \times 10^{-2} M_\odot$ is required to match the Lupus median accretion rate. For long viscous timescales, the behaviors change. For $\langle t_{\nu,0} \rangle = 2$ Myr, slightly higher disk masses of $\simeq 3 \times 10^{-2} M_\odot$ are required. Longer viscous timescales are excluded since they produce too low accretion rates at 2 Myr even for gravitationally unstable disks ($\langle M_0 \rangle = 0.1 M_\odot$). As mentioned above, for even longer viscous timescales, the synthetic population cannot match the disk fraction.

To find the value of the median disk mass that matches the median accretion rate for every value of $\langle t_{\nu,0} \rangle$ and $\langle R_0 \rangle$, we interpolate the value of $\langle M_0 \rangle$ and $\langle \dot{M}_{PEW} \rangle$. For that, we start with a grid of $\langle R_0 \rangle$, $\langle t_{\nu,0} \rangle$, and $\langle M_0 \rangle$ over which the value of $\langle \dot{M}_{PEW} \rangle$ has been computed in stage 1. We then find, for each value of $\langle t_{\nu,0} \rangle$ and $\langle R_0 \rangle$, the two populations that predict the accretion rates at the age of Lupus bracketing the observed rate of $10^{-9} M_\odot$ yr$^{-1}$. We finally compute the values of $\langle M_0 \rangle$ and $\langle \dot{M}_{PEW} \rangle$ using a linear interpolation. We further check that the populations with the interpolated values of $\langle M_0 \rangle$ and $\langle \dot{M}_{PEW} \rangle$ does predict a disk dispersal time of about 3 Myr and a median accretion rate of $10^{-9} M_\odot$ yr$^{-1}$ at 2 Myr.

The result of the method is shown in Figure 10, top panels. It highlights that the two observational constraints used so far, namely the disk fraction and accretion rate in Lupus, can be framed as constrained mass-loss rates $\langle \dot{M}_{PEW} \rangle$ and initial disk

masses $\langle M_0 \rangle$ as a function of the two free parameters: the viscous timescale $\langle t_{\nu,0} \rangle$ and initial disk size $\langle R_0 \rangle$. In Figure 10, we recover the result of Figure 9(c) that the initial median disk mass required to fit disk dispersal and accretion rate in Lupus for $\langle R_0 \rangle = 10$ au is constant for short viscous timescales and increases at longer viscous timescales. The corresponding mass-loss rate also increases with the viscous timescale as suggested by Equation (17). Its value is typically above $10^{-9} M_\odot$ yr$^{-1}$ and can reach a few times $10^{-8} M_\odot$ yr$^{-1}$ for long viscous timescales and compact disks. Interestingly, compared to the value predicted by the hydrodynamical model of G. Picogna et al. (2021), we find that only very specific values of $\langle t_{\nu,0} \rangle$ and $\langle R_0 \rangle$ are in line with their theoretical results.

### 3.2.3. Constraints from the AGE-PRO Results for Lupus

At this stage, we are in a position to compare the results of AGE-PRO to synthetic populations which are consistent with disk dispersal and the median accretion rate in Lupus. The leftover free parameters are $\langle t_{\nu,0} \rangle$ (or, equivalently, $\langle \alpha_{SS} \rangle$) and $\langle R_0 \rangle$. In Figure 10, bottom panels, we show the Lupus median CO size and median gas mass for all the synthetic populations that fit the disk fractions and the median accretion rate in Lupus. As expected, a population of initially larger disks or shorter viscous timescales predicts a larger $\langle R_{CO} \rangle$ size after 2 Myr. The median CO gas size is particularly sensitive to $\langle t_{\nu,0} \rangle$ for initially compact disk populations with a drop in $\langle R_{CO} \rangle$ with $\langle t_{\nu,0} \rangle$ for $\langle R_0 \rangle \le 20$ au. This is not only due to limited viscous spreading, but also to the onset of disk dispersal by outside-in pathways occurring for disks that remain compact over their lifetime.

Compared with the AGE-PRO measurements of the CO gas size in Lupus, it is immediately apparent that only initially compact disks can fit the observations with a median disk radius of $\langle R_0 \rangle \simeq 5–20$ au (see models outlined by black circles and squares in Figure 10(c)). We also find a degeneracy
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Best-fit viscous model for $t_{v,0} = 0.4$ Myr (left), $t_{v,0} = 1.0$ Myr (middle), and $t_{v,0} = 3.0$ Myr (left) compared with the observations. The figure follows the same convention as Figure 7.</figcaption></figure>

between initially compact disks with short viscous timescales and larger disks with longer viscous timescales. This results in a relatively well-constrained $\alpha_{SS}$ parameter of $\alpha_{SS} = 2-4 \times 10^{-4}$. In the next section, we adopt the three best-fit models highlighted by black points in Figure 10.

The predicted median disk mass at 2 Myr appears to be the main tension with the observations. All the models that reproduce the median CO size consistently produce a median disk mass of about $5 \times 10^{-3} M_{\odot}$ (see models outlined by black circles in Figure 10(d)). This is a factor of 5 to 10 above the disk mass estimated for Lupus by AGE-PRO (see the gray shaded area in Figure 10, bottom right). The only synthetic populations that are closer to the median Lupus disk mass have compact disks and long viscous timescales. In this extreme regime, photoevaporation removes a significant fraction of the disk mass because most of the disk mass resides where the mass-loss rate profile peaks. However, these models strongly underpredict the $\langle R_{CO} \rangle$ size since the low disk masses are reproduced at the expense of a dramatic shrinking of the disk.

### 3.2.4. Best-fit Viscous Model with Photoevaporative Winds

Putting together all the constraints stemming from the properties of the Lupus population and disk fraction, we find that the disks should be initially compact ($R_0 \simeq 5-20$ au) with relatively long viscous timescales of $t_{v,0} \simeq 0.4-3$ Myr even though no models can simultaneously reproduce the disk accretion rate, disk mass, and size in Lupus. We show three best-fit models in Figure 11 encompassing the parameter space of the best fit and already outlined in Figure 10 (see black circles and squares). By construction, the models reproduce the median values of $\langle R_{CO} \rangle$, and $\langle M_* \rangle$ at the age of Lupus (2 Myr).

We recover that the essential tension between the models and the observations is the median disk mass. In both models, the median disk mass is remarkably constant even after 10 Myr. This is not only the result of a slow evolution of

individual disks since the disk mass is rather flat even after 10 times the viscous timescale but also the result of disk dispersal, which preferentially removes light disks. Such a survivorship bias on median disk mass or accretion rates can also be seen in former population models but, to our knowledge, has never been discussed (e.g., R. D. Alexander & P. J. Armitage 2009). The drop in median disk mass between Oph and Lupus cannot be reproduced by the turbulence-driven model.

The median accretion rate is predicted to drop by about a factor of 5–10 within the first 2 Myr before experiencing a relatively shallow decline at a longer time. Compared with the evolution of the median accretion rates between Lupus and Upper Sco, the best-fit models tend to overestimate the accretion rate in Upper Sco.

In the three models, we predict an increase in the CO gas size with time but the Upper Sco AGE-PRO population is not old enough to test this increase. Therefore, the median CO disk size for Upper Sco is consistent with the models. An increase in gas size is expected from turbulence-driven models due to disk spreading (L. Trapman et al. 2020; C. Toci et al. 2023). However, disk dispersal also removes disks with small radii (i.e., short viscous timescales), enhancing the increase in disk size. We also refer to R. Anania et al. (2025) for a detailed discussion of the effect of external photoevaporation in Upper Sco.

# 4. Discussion

## 4.1. Constraints on Disk Evolution Mechanisms

### 4.1.1. Summary of the Results

Our study demonstrates that reproducing only the median values of the fundamental population parameters such as the accretion rate, disk fraction, disk size, and disk mass is a discriminant test for disk evolution models. We choose to
<table><caption>**Table 2**
Summary of the Best-fit Parameters for the Two Scenarios</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Best-fit Value</th></tr></thead><tbody><tr><td colspan="3">MHD wind model</td></tr><tr><td>Initial disk size</td><td>R<sub>0</sub></td><td>10 au</td></tr><tr><td>Initial disk mass</td><td>M<sub>0</sub></td><td>7.5 × 10<sup>-3</sup> M<sub>☉</sub></td></tr><tr><td>α parameter</td><td>α<sub>DW</sub></td><td>4.6 × 10<sup>-4</sup></td></tr><tr><td>Accretion timescale</td><td>t<sub>acc,0</sub></td><td>0.75 Myr</td></tr><tr><td>ω parameter</td><td>ω</td><td>0.5</td></tr><tr><td>Magnetic lever arm</td><td>λ</td><td>>8<sup>b</sup></td></tr><tr><td>Ejection-to-accretion ratio</td><td>f<sub>M</sub></td><td>&lt;0.5</td></tr><tr><td colspan="3">Turbulent model<sup>a</sup></td></tr><tr><td>Initial disk size</td><td>R<sub>0</sub></td><td>10 au</td></tr><tr><td>Initial disk mass</td><td>M<sub>0</sub></td><td>1.2 × 10<sup>-2</sup> M<sub>☉</sub></td></tr><tr><td>α parameter</td><td>α<sub>SS</sub></td><td>3.4 × 10<sup>-4</sup></td></tr><tr><td>Viscous timescale</td><td>t<sub>ν,0</sub></td><td>1.0 Myr</td></tr><tr><td>Wind mass-loss rate</td><td>M&#775;<sub>PEW</sub></td><td>4.4 × 10<sup>-9</sup> M<sub>☉</sub> yr<sup>-1</sup></td></tr></tbody></table>

**Notes.**

<sup>a</sup> No set of parameters for the turbulence-driven model fits all the observational constraints. The parameters correspond to a model that fits the accretion rate but not the disk mass. There is also a degeneracy between R<sub>0</sub> and t<sub>ν,0</sub> outlined in Figure 10.

<sup>b</sup> Above this value the accretion-to-ejection ratio is small and λ has a negligible impact on the properties of the synthetic population.

follow the step-by-step approach summarized in Figure 4, adopting simple disk evolution models and assuming that Oph, Lupus, and Upper Sco are representative of the same population at different ages. The best-fit parameters obtained in the previous section are summarized in Table 2.

We find that the adopted MHD wind-driven model can reproduce very well the main properties of the three populations with the gas mass and gas size estimates provided by the AGE-PRO ALMA large program (see Figure 7). Notably, the best-fit accretion timescale of t<sub>acc,0</sub> ≈ 0.75 Myr along with an initial disk size of ⟨R<sub>0</sub>⟩ = 10 au give an initial value of α<sub>DW</sub> = 4.6 × 10<sup>-4</sup>. Following B. Tabone et al. (2022a), α<sub>DW</sub> is related to the disk magnetization often parameterized by the thermal-to-magnetic pressure ratio denoted as β. Our constraints translate to an initial value of β<sub>0</sub> ≈ 10<sup>5</sup>, depending on the unknown B<sub>☉</sub>/B<sub>z</sub> in the disk atmosphere (see Section 4.3 in B. Tabone et al. 2022a, for a discussion).

In contrast, following our fitting procedure, we find that turbulence-driven models with constant α<sub>SS</sub> and the photo-evaporative mass-loss rate profile of G. Picogna et al. (2021) can match successively the disk fraction, disk accretion rate, and disk size (Figure 11). This requires relatively long viscous timescales t<sub>ν,0</sub> = 0.4–3 Myr, compact disks of ⟨R<sub>0</sub>⟩ ≈ 5–20 au, and low wind mass-loss rates of ⟨M&#775;<sub>PEW</sub>⟩ ≲ 10<sup>-8</sup> M<sub>☉</sub> yr<sup>-1</sup>. However, when the three previously mentioned observational constraints are matched, the median disk mass at the age of Lupus and Upper Sco is overestimated by a factor of 3–10.

### 4.1.2. Origin of the Discrepancy with Turbulent Models

Compared with e.g., Markov Chain Monte Carlo fits, our step-by-step fitting method summarized in Figure 4 has the advantage of being intuitive and computationally efficient at finding best-fit solutions. However, in the absence of a set of parameters that fits all the constraints, the order in which the observations are fitted matters. In the case of turbulent models,

one can naively conclude that the disk mass is the experimentum crucis. However, the observed Lupus median disk mass can also be matched in stage 2 instead of adjusting for the observed accretion rate. By doing so, we identified the two fundamental origins of the discrepancy between the turbulence-driven scenario and the observations.

First, matching median disk mass derived by AGE-PRO in stage 2, requires low initial disk masses in the range of 1–2 × 10<sup>-3</sup> M<sub>☉</sub>. This is a factor of 2–5 smaller than the median disk mass inferred in the young Oph sample. In other words, the turbulence-driven model fails to reproduce the significant drop in median disk mass between Oph and Lupus. This is primarily due to the slow evolution of disk mass for viscous disks: ignoring dispersal, the mass of individual disks drops by a factor of (1 + t/t<sub>ν,0</sub>)<sup>1/2</sup>. Even for t<sub>ν,0</sub> = 0.2 Myr, this only gives a drop of a factor of 3 from Oph to Lupus. Disk dispersal worsens this trend by preferentially removing light disks and reducing the dropping rate of the median disk mass, in particular at the age of Upper Sco. In contrast, MHD wind-driven accretion can produce a quick drop in disk mass and match the decline of the median disk mass from Oph to Lupus. However, we acknowledge that the disk masses in Oph versus Lupus and Upper Sco are estimated following slightly different approaches (L. Trapman et al. 2025a) such that the dropping rate of disk mass remains a loose constraint.

The major tension between the turbulence-driven model and the observations lies in the ratio between the disk mass and the accretion rate, also called the apparent disk lifetime t<sub>lt</sub> = M<sub>D</sub>/M&#775;<sub>*</sub>. Specifically, the fact that the apparent disk lifetime estimated by AGE-PRO is on average shorter than the age of the disk excludes turbulence-driven accretion with a constant α<sub>SS</sub> profile. We show in Figure 12 the AGE-PRO sources in the M<sub>D</sub> − M&#775;<sub>*</sub> plane, focusing on stellar ages between 1 and 3 Myr. This includes the 10 Lupus and two Upper Sco sources. Figure 12 also depicts the predictions for the best-fit MHD disk-wind model and turbulence-driven models as probability density maps. Because the uncertainties on the disk mass and the variability in the accretion rates can strongly affect the correlations in the M<sub>D</sub> − M&#775;<sub>*</sub> plane, we take into account a lognormal error on the disk mass of σ = 0.7 dex (L. Trapman et al. 2025a) and on the accretion rate of σ = 0.3 dex (L. Venuti et al. 2014).

We recover that the turbulence-driven model predicts a correlation between disk mass and accretion rate (G. P. Rosotti et al. 2017). This is reminiscent of the population synthesis models of G. Lodato et al. (2017) and G. D. Mulders et al. (2017) with the addition that our population model also includes disk dispersal by photoevaporation, which alters the simple linear relationship. In particular, the predicted correlation flattens toward low-mass disks producing a subpopulation of low-mass disks (M<sub>D</sub> ≲ 10<sup>-4</sup> M<sub>☉</sub>) with a short apparent disk lifetime of M<sub>D</sub>/M&#775;<sub>*</sub> ≈ 0.1 – 1 Myr. This behavior contrasts with the results of A. Somigliana et al. (2020) who find that photoevaporation tends to produce lower accretors in the low-mass part of the M<sub>D</sub> − M&#775;<sub>*</sub> plane (bottom left region in Figure 12). This difference is due to the initial disk sizes, which are now constrained by the AGE-PRO data. In our best-fit model, the disks, which are on average compact, tend to experience an outside-in dispersal pathway outlined in Figure 1, leading to short M<sub>D</sub>/M&#775;<sub>*</sub> before dispersal.

Our AGE-PRO sample does not recover the correlation between disk masses and accretion rates found by
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Comparison in the disk mass-accretion rate plane between the middle-aged AGE-PRO sources and our best-fit population models for the turbulence-driven and MHD wind-driven case. The AGE-PRO sources with accretion rates from C. F. Manara et al. (2020) are shown in red. The simulated disks are presented in black points and correspond to a population of 2000 disks. The probability density is plotted in blue shade and the 2D cumulative density distribution is in gray contours (contours encompassing 50% and 90% of the disks). To take into account the spreads due to the variability of the accretion rate, and the uncertainties on the measured accretion rate and disk mass, we assume that each modeled disk generates a lognormal probability distribution in the $\dot{M}_{*} - M_{D}$ plane. The best-fit MHD disk-wind model corresponds to $\omega = 0.5$ and the best-fit turbulence-driven model to $\langle t_{v,0} \rangle = 1.0$ Myr and $\langle R_{0} \rangle = 10$ au. The modeled distributions are computed assuming a detection threshold of $M_{D} = 10^{-5} M_{\odot}$.</figcaption></figure>

C. F. Manara et al. (2016) for larger samples and using dust emission. This is likely due to too low statistics and large uncertainties in the retrieved gas mass. As detailed in Section 3 and further highlighted in Figure 13 (top and middle panels) the turbulent model fits well the accretion rates but over-predicts the disk masses. However, we recall that the model can alternatively fit the disk masses but would then under-predict the accretion rates. The true discrepancy between the turbulent model and the observations is the short apparent disk lifetime $M_{D}/\dot{M}_{*}$ found in AGE-PRO. As seen in the $M_{D}-\dot{M}_{*}$ plane (Figure 12) and in the cumulative distribution of $M_{D}/\dot{M}_{*}$ (Figure 13, bottom panel), half of the AGE-PRO sample is below 1 Myr with two disks below 0.1 Myr. As demonstrated in M. G. Jones et al. (2012) and G. P. Rosotti et al. (2017), turbulent-driven evolution leads to $M_{D}/\dot{M}_{*}$ longer or about the disk age, regardless of the details of the turbulence-driven model. This result is fully recovered in our synthetic population in Figure 12 (right panel) and Figure 13 (bottom panel) with the difference that the large uncertainty in disk mass makes very short $M_{D}/\dot{M}_{*}$ a possible but unlikely outcome of retrieval. In addition, outside-in disk dispersal, discussed by G. P. Rosotti et al. (2017) in the context of external photoevaporation, produces low-mass disks with short $M_{D}/\dot{M}_{*}$ ratios. Still, because these disks are short-lived, they should represent a minor fraction of sources. Moreover, they cluster in the very low disk mass region. Therefore outside-in disk dispersal cannot account for the short disk lifetimes found by AGE-PRO. We conclude that the tension with turbulence-driven models arises from the combination of accretion rate and disk mass with a disk-mass-to-accretion-rate ratio too short to be matched by population synthesis models.

In contrast, the MHD wind-driven model reproduces very well the distributions of disk gas masses, accretion rates, and disk lifetimes (Figure 12, left panel, and Figure 13, red curve). The simulated population exhibits a large scatter, even before taking into account the uncertainty on the derived disk mass and the variability in the accretion rate (see black points versus probability distribution in Figure 12). Overall, the 2D probability distribution encompasses very well the AGE-

PRO measurements. The short disk lifetimes reproduced by the model are due to the evolutionary tracks of individual disks (see B. Tabone et al. 2022a, for a detailed discussion): before disk dispersal, a disk moves horizontally in the $M_{D}-\dot{M}_{*}$ plane, producing a population of low-mass disks with high accretion rates. It is this population of disks about to be dispersed that matches the disks with a very short disk lifetime. Interestingly, disks with low masses and short disk lifetimes are also compact. This is consistent with the MHD wind dispersal process since the disk dispersal time is set by the ratio between disk size and $\alpha_{DW}$.

## 4.2. Perspectives on the Observational Data

### 4.2.1. Representativeness of the AGE-PRO Star-forming Regions

One of the fundamental assumptions behind our analysis is that the Oph, Lupus, and Upper Sco populations are representative of the same population at different ages. Based on the analysis presented in the Oph, Lupus, and Upper Sco AGE-PRO papers (C. Agurto-Gangas et al. 2025; D. Deng et al. 2025; D. A. Ruiz-Rodriguez et al. 2025) the AGE-PRO sample is representative of the whole population in each region, based on a comparison with the available accretion rates (not applicable to Oph) and millimeter luminosities for similar spectral types. Surveys of young star-forming regions of a similar age as Lupus (1–3 Myr) show little variation in their distribution of continuum millimeter luminosity and accretion rates except for Corona Australis, which shows weaker millimeter continuum fluxes (P. Cazzoletti et al. 2019). Our knowledge of old star-forming regions similar to Upper Sco is more limited. It remains unclear whether the median accretion rate declines with cluster age, as the measurements of M. Fang et al. (2023) used in the present work suggest, or if they remain comparable to young star-forming regions even when corrected for the stellar mass dependency (see, e.g., L. Testi et al. 2022, Figure 11). Surveys of old star-forming regions in the millimeter are also lacking, making Upper Sco our unique reference point with a complex star formation history (N. Miret-Roig et al. 2022;
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Cumulative distribution of the accretion rate, the disk mass, and the disk lifetime for the best-fit MHD disk-wind (red) and turbulence-driven models (blue) vs. the AGE-PRO sample focusing on 1–3 Myr disks (black). The MHD wind model reproduces well the distributions. The salient discrepancy between the turbulence-driven model and the observations lies in the $M_D/\dot{M}_*$ ratio. The modeled cumulative distributions are computed assuming a detection threshold of $M_D = 10^{-5} M_\odot$.</figcaption></figure>

S. Ratzenböck et al. 2023). Therefore, variations in terms of gas content between different star-forming regions remains to be probed, targeting the well-known middle-aged but also less-known old star-forming regions.

Another important caveat for the Upper Sco region is the importance of external irradiation in the disk evolution. Throughout this work, we neglected external photoevaporation, which can affect disk evolution (G. A. L. Coleman & T. J. Haworth 2022). Only strong external FUV radiation of $G_0 > 10^3$ is expected to significantly affect the dispersal time and all the properties of the disks (M. Ansdell et al. 2017; K. Maucó et al. 2023; M. Gárate et al. 2024), but the modest FUV strength experienced by the AGE-PRO Upper Sco sample can slightly affect the disk gas mass but largely decrease disk gas size (R. Anania et al. 2025).

### 4.2.2. Disk Gas Mass

It is commonly suggested that the measurement of the disk size alone constitutes a decisive test for disk evolution models. This comes from the idea that turbulence-driven accretion

requires viscous expansion and MHD wind-driven accretion does not. Our work demonstrates that disk sizes constitute an essential constraint for disk evolution models but also shows that the latter idea is somewhat misleading. First, as demonstrated in Figure 1, turbulence-driven disks can shrink over time due to internal photoevaporation. The requirement is that the disk is compact enough such that the local wind mass-loss rate is large enough out to the disk edge. Since AGE-PRO shows that the disks are on average compact ($R_{CO} \simeq 70–110$ au), internal photoevaporation is expected to affect the evolution of the disk size in the case of turbulence-driven evolution. External photoevaporation would be even more efficient (R. Anania et al. 2025). On the other hand, the survivorship bias can greatly affect disk sizes: in MHD wind-driven models, we find that CO disk size slightly increases due to the removal of compact disks. The measurement of CO size alone, even in regions of different age, is not a decisive test for disk evolution scenarios.

However, combining disk gas size, disk gas mass, and stellar accretion rate measurements brings solid constraints on disk evolution models. This is demonstrated in Figure 10, which shows that turbulence-driven models cannot simultaneously reproduce disk mass, accretion rate, and CO size. Overall, we find that disk mass is a crucial parameter along with accretion rate, but is indubitably the most uncertain. It is the relatively low median disk gas masses found in Lupus and Upper Sco that question the turbulence-driven paradigm.

If disk masses in these regions were to be systematically underestimated by a factor of 5 to 10, turbulence-driven models would be compatible with the data with an apparent lifetime of $M_D/\dot{M}_* > 1$ Myr. MHD wind-driven models would also be able to fit the data with a higher ejection-to-accretion mass ratio $f_M$ and higher initial disk mass $\langle M_0 \rangle$. This can be seen in Figure 6: if the median disk mass at Lupus' age were to be an order of magnitude higher, the best-fit ejection-to-accretion mass ratio would be about $f_M = 1–10$. Such high values of the ejection rate are in fact not unrealistic compared with theoretical MHD wind models and numerical simulations (e.g., L. Wang et al. 2019). In this case, the initial disk mass would need to be increased (see Figure 5).

The estimates of the disk gas mass are subject to considerable challenges (A. Miotello et al. 2023). Within the AGE-PRO large program, the DALI thermochemical model is used (L. Trapman et al. 2025a) to retrieve disk gas masses from CO and $N_2H^+$ but the retrieval based solely on the CO line fluxes and using the M. Ruaud et al. (2022) models give similar results in Lupus (see D. Deng et al. 2025). Overall, our estimates depend on our understanding of the chemistry and physics of disks and one can expect considerable progress in the future. Benchmarks between mass estimates based on $N_2H^+$ (D. E. Anderson et al. 2019; L. Trapman et al. 2022), kinematics (P. Martire et al. 2024), and hydrodynamics (L. Trapman et al. 2017) show relatively consistent results for large and massive disks (see A. Miotello et al. 2023, for a review). However, we stress that the lowest disk masses found in AGE-PRO are for very compact disks. Deep observations designed to characterize the chemistry of compact disks ($\langle R_{CO} \rangle \lesssim 60$ au) should be performed. Our fitting strategy allows one to investigate how the best-fit model will be affected by a revision of the disk gas masses since the constraints from the gas masses are, on purpose, taken into account at the final stage of the fit (see Figure 4).
<figure><img src="image_16.png" /><figcaption>**Figure 14.** The survivorship bias illustrated by the initial median value of the disk mass (left) and size (right) calculated over the surviving disks. The blue and red lines correspond to the best-fit MHD wind-driven and turbulence-driven populations (see Table 2).</figcaption></figure>

A final caveat behind the disk mass estimates is that the molecular lines trace the outer disk where it is assumed that the bulk mass resides. The regions inside of ≈10 au could hide a large fraction of the mass, which would remain poorly traced by ALMA due to optical depth effects. This would require a significant jump in the surface density profile. Since disks accrete, a jump in the surface density profile implies a jump in $\alpha$ values (either radial or vertical torque). In the turbulence-driven scenario, a massive inner disk with very low $\alpha$ values is expected to have a long survival timescale since photoevaporation does not operate in that region (L. Hartmann et al. 2006; R. Morishima 2012; M. Gárate et al. 2021). One can then speculate that such a setup would lead to a population of extremely compact but accreting disks. Detailed models including radial variation in $\alpha$ values (e.g., S. Tong et al. 2024) are deferred to future work. In parallel, constraints on the absolute value of the column density in the inner disk can be obtained by line pressure broadening of CO (T. C. Yoshida et al. 2022).

### 4.3. Disk Dispersal and the Survivorship Bias

Because the evolution of a disk cannot be observationally monitored, the use of surveys at a population level is needed. However, in the presence of disk dispersal, the evolution of the properties of a population is not only determined by the evolution of individual disks but also by the properties of the survivors. Our population synthesis approach highlights the importance of this bias, aka survivorship bias, when interpreting the results of surveys. The evolution of the median disk mass in the turbulence-driven scenario is undoubtedly the most striking example: whereas the mass of individual disks declines, disk dispersal removes light disks resulting in a shallow evolution of the median disk mass. This selection process, identified by A. Somigliana et al. (2020), is evidenced by the dramatic increase in the initial mass of the surviving disks shown in Figure 14 (see left panel). MHD wind-driven evolution also introduces a survivorship bias: compact disks evolve and disperse faster since they have shorter accretion timescales $t_{acc,0}$. This results in an increase in the initial median disk size of the surviving disks highlighted in Figure 14 (right panel). This effect drives the small increase in CO gas size seen in Figure 7.

We stress that the survivorship bias is different from the selection bias induced by the definition of a sample. It is a bias that does not depend on the choice of the observer but arises because the defined population of disk loses members due to dispersal. One could instead follow the population of all the stars but in this case, the median disk mass drops to zero as soon as the disk fraction

drops below 50%. We call it bias because it is a selection process that prevents us from probing directly the evolution of the disk properties. The implication is that the survivorship bias is a fundamental bias that is difficult to correct in the observations. This requires knowing which disk in a young population survives at a later stage. For a given disk evolution model, the effect of the survivorship bias can still be qualitatively assessed since the population parameters that are the most affected are those that control disk dispersal. On the other hand, survivorship bias constitutes a promising way to constrain disk evolution mechanisms. As an example, X-ray photoevaporation is expected to disperse not only light disks but also disks around X-ray bright stars. This idea has been explored to test the relevance of X-ray photoevaporation (T. Preibisch et al. 2005; S. Flaischlen et al. 2021; S. Laos et al. 2022) even though the interpretations of the results are complicated by the possible quenching of X-ray emission by the accretion flow.

### 4.4. Spreads and Correlations between Disk Parameters

Throughout this work, we assumed that $\alpha$, $R_0$, $M_0$, and the parameter controlling the wind mass-loss rate ($\dot{M}_{PEW}$ or $f_M$) are independent and follow lognormal distributions with fixed spreads. We recall that the spreads have been chosen to reproduce the decline of the disk fraction but no further adjustment has been made to fit the spread in the observed disk masses and accretion rates. We however note that our best-fit models reproduce well the spread in the accretion rates (Figure 13, middle panel), especially if we consider the effect of variability. The spread in the distribution of disk masses is also well matched (see Figure 13, top panel) even though the spreads are dominated by the large statistical uncertainties in the retrieved gas masses.

Interestingly, we note that the spreads in the disk parameters impact the evolution of the medians. This is again due to the survivorship bias. The clearest example is the evolution of the median CO disk size in MHD wind-driven evolution. Disks are removed according to their $t_{acc,0}$ value, which is the ratio between $\alpha_{DW}$ and $R_0$. This drives the increase in $\langle R_{CO} \rangle$. However, in the absence of any spread in initial disk size, $\langle R_{CO} \rangle$ would remain constant in time.

The effect of the correlations between the disk parameters is an interesting avenue to explore. Due to survivorship bias, a parameter that is not related to disk dispersal can be greatly affected if it is correlated with a key parameter for dispersal. The most obvious correlation is probably between the initial disk mass and the initial disk size even though observations of Class I disks suggest a rather weak dependency between disk mass and disk size (J. J. Tobin et al. 2020). For example, the disks could be born with similar values of the surface density and differ in mass only because they differ in size. In the case of our MHD wind-driven model, small disks are dispersed first such that a positive correlation would produce a shallower evolution in disk mass. For our turbulence-driven model, light disks are dispersed first such that a positive correlation would steepen the evolution of the CO disk size.

All in all, our best-fit models constitute a baseline for a full exploration of the effect of the spreads and the correlations between the disk parameters. Such models will add new free parameters and additional observational constraints should be taken into account to limit the degeneracies.
## 4.5. Limits and Perspective on One-dimensional Disk Evolution Models

Our work aims to test and quantify disk evolution processes by confronting existing surveys and in particular AGE-PRO to population synthesis models. Our goal is not to propose new disk evolution models but to provide reference population models based on existing 1D evolution models. These models are meant to capture the essential features of either turbulence-driven evolution or MHD wind-driven evolution. Here, we outline the simplifications of these models and discuss future disk evolution models that remain to be built.

### 4.5.1. Turbulence-driven Case

For the turbulence-driven model, one of the uncertainties is the photoevaporative wind model (I. Pascucci et al. 2023). The wind properties are the result of the spectral distribution of high-energy photons (from ≈10 eV to ≈1 keV), radiative transfer, and gas thermal balance of the disk upper layers. Over the past decades, several studies have been conducted to establish mass-loss rate profiles along with absolute mass-loss rates with contradictory results (L. Wang & J. Goodman 2017; A. Komaki et al. 2021; G. Picogna et al. 2021). In our work, we kept the total mass-loss rate as a free parameter to avoid biases due to the use of a specific study and we adopted the mass-loss rate profile of G. Picogna et al. (2021) for a 0.5 $M_\odot$ mass star. Our best-fit models have a typical mass-loss rate of ≈2 × 10⁻⁹–10⁻⁸ $M_\odot$ yr⁻¹, which is smaller than the mass-loss rate of ≈2 × 10⁻⁸ $M_\odot$ yr⁻¹ predicted by G. Picogna et al. (2021) for a 0.5 $M_\odot$ stellar mass. Similar results were also found by analyzing the distribution of the accretion rates alone (R. Alexander et al. 2023). Our findings are also well in line with the recent work of A. D. Sellek et al. (2024) who showed that the inclusion of collisions of O with neutral H, which was not included in G. Picogna et al. (2021), and subsequent line cooling, decreases the mass-loss rate by a factor of 5–10.

Our study demonstrates that internal photoevaporative winds control not only the disk dispersal time but also the overall evolution of the disk properties, notably the disk size as illustrated in Figure 1. In this context, the exact mass-loss rate profile matters. If the local mass-loss rate were to be significant at large distances, as in the calculations of U. Gorti & D. Hollenbach (2009), disks could be initially more extended.

An important simplification of the model is the absence of attenuation of stellar radiation by an inner MHD disk wind. MHD disk winds launched from the inner disks (≈1 au) are likely ubiquitous (M. N. Simon et al. 2016; M. Fang et al. 2018; I. Pascucci et al. 2020) and attenuation could reduce the heating of the outer disk where photoevaporative winds are launched. Our estimate of the median wind mass-loss rate should therefore be considered as time averaged. Consistent models would require identifying which photon energy drives the wind and the composition of the inner MHD wind (dust content or metallicity). More detailed models including attenuation by an evolving inner MHD disk wind should be explored (see also J. Weder et al. 2023).

The assumption of a constant $\alpha_{SS}$ in space and time is also an important simplification of the turbulence-driven model used in this work. The presence of a dead zone with a limited radial extent can create a steep transition in the $\alpha_{SS}$ profile (T. N. Delage et al. 2022). A systematic study exploring sophisticated profiles of $\alpha_{SS}$ is beyond the scope of the paper. Preliminary tests show that a population of disks with a highly turbulent outer region ($r > 20$ au)

and low-$\alpha_{SS}$ inner disk tends to behave like a population of disks with constant $\alpha_{SS}$ values. In particular, an $\alpha_{SS}$ with a step function does not seem to solve the tension with the observations but more work remains to be done in that regard.

### 4.5.2. Magnetohydrodynamic Wind-driven Case

In contrast to turbulence-driven models, 1D disk evolution models including MHD disk winds are relatively recent (X.-N. Bai 2016; T. K. Suzuki et al. 2016; B. Tabone et al. 2022a). Throughout this work, we adopted the simple model of B. Tabone et al. (2022a), which has the advantage of relying on limited assumptions and being directly comparable to numerical simulations. The wind torque is parameterized by the $\alpha_{DW}$ parameter, which is proportional to the mid-lane magnetization, traditionally quantified by the thermal-to-magnetic pressure ratio $\beta$.

The secular evolution of the magnetic field strength, which controls the wind torque, is notably a major uncertainty in secular evolution models. As shown by P. J. Armitage et al. (2013) and X.-N. Bai (2016), if the disk preserves its magnetic field, the disk can be quickly dispersed. This is the situation explored in this work where the $\omega$ parameter controls the time evolution of $\alpha_{DW}$. In the framework of this model, we show that a slowly decreasing magnetic field producing an increased magnetization can naturally account for the observational data. Numerical simulations running on a secular timescale opened by GPU-accelerated codes (see, e.g., G. R. J. Lesur et al. 2023) are warranted to interpret these constraints.

The wind mass-loss rate and the impact of high-energy photons are also major uncertainties in MHD wind-driven models. In this work, we use the magnetic lever arm parameter $\lambda$ to parameterize the wind mass-loss rate (see Equation 6). In a realistic setup, the value of $\lambda$ is the result of the field morphology, field strength, and heating of the disk surface (see, e.g., X.-N. Bai 2016). The contrast between photoevaporative winds and MHD disk winds is commonly made in the literature. We recall that MHD wind winds emerge when a magnetic field with a net flux threatens the disk but thermal effects are fundamental in the wind launching process (F. Casse & J. Ferreira 2000; X.-N. Bai 2016). The transition between photoevaporative winds and MHD disk winds is often presented as the transition when the MHD wind mass-loss rate is about that predicted in the absence of any net magnetic flux. W. Béthune et al. (2017) found that the wind mass-loss rate decreases with disk magnetization as $\dot{M}_{MHD-DW} \propto \beta^{-1/2}$. P. J. Rodenkirch et al. (2020) further confirmed this trend showing that for a given disk setup, the wind mass-loss rate reaches the threshold value of photoevaporation for $\beta \gtrsim 10^7$. Within our modeling framework, it means that $\lambda$ is predicted to depend on $\beta$ or, equivalently on $\alpha_{DW}$. However, in the absence of systematic theoretical constraints on how $\lambda$ varies with the disk’s fundamental properties, we simply assumed $\lambda$ to be constant in space and time. We also note that our best-fit MHD disk-wind model predicts disks with increasing magnetizations but decreasing mass-loss rates. In the literature, 1D evolution models have attempted to couple MHD disk winds and photoevaporative winds by simply adding the MHD disk-wind mass-loss rate and the photoevaporative mass-loss rate (M. Kunitomo et al. 2020; J. Weder et al. 2023). A promising avenue to consistently describe the transition between photoevaporation and MHD disk wind would be to determine how $\lambda$ varies with at least $\alpha_{DW}$, and with the
impinging radiation field using a wide grid of numerical
simulations including detailed microphysics.

5. Conclusion

In this work, we developed a disk population synthesis approach to challenge existing disk evolution models in the framework of the AGE-PRO large program. The use of population synthesis is essential since the evolution of a population is not only set by the evolution of each individual disk but also by the properties of the survivors. To obtain robust constraints we used not only the results of AGE-PRO, namely the disk gas sizes and masses, but also the disk fractions and the accretion rates. Each free population parameter is constrained sequentially using the disk fractions and the properties of the middle-aged Lupus population. This allows us to considerably reduce the parameter space before confronting the two scenarios with the AGE-PRO results. Our results are the following.

1. The analysis of disk surveys requires adopting a disk population synthesis approach to reproduce all the available observational constraints simultaneously. This approach is key to include the fact that the evolution of a population is not only the result of the evolution of individual disks but also of disk dispersal, which selects disks with specific properties, a phenomenon that we call here survivorship bias.

2. MHD wind-driven models mimicking the evolution of the disk magnetization can reproduce the disk fractions and the overall properties of the Oph, Lupus, and Upper Sco populations. It requires initially compact ($R_0 \simeq 10$ au) disks with a moderate mass-loss rate $\dot{M}_W/\dot{M}_* \lesssim 1$, and an accretion timescale of about $t_{acc,0} \simeq 0.4-0.8$ Myr ($\alpha_{DW,0} = 5 \times 10^{-4}-10^{-3}$). This corresponds to disks with an initial magnetization of about $\beta \simeq 10^5$.

3. Turbulence-driven models with internal photoevaporation can reproduce disk dispersal and gas size but fail to reproduce the accretion rates and the disk masses simultaneously. The short apparent disk lifetimes of $M_D/\dot{M}_* \simeq 0.1-1$ Myr for half of the middle-aged disks are the major tension with turbulence-driven models. Still, detailed studies of compact disks that have the lowest masses and the shortest apparent disk lifetimes should be conducted to confirm their low masses.

This work tends to favor the emerging MHD wind-driven
disk accretion scenario, which has major implications for
planet formation and migration. The increasing number of disk
wind detections with ALMA and JWST will enable constrain-
ing some of the key input parameters of these models such as
the wind mass-loss rates and the lever arm. This work provides
the community with observationally constrained synthetic disk
populations that are required to study dust evolution and planet
formation. We, however, acknowledge that the disk evolution
models cannot capture the full complexity of the disk physics.
This first study warrants the development of more sophisti-
cated disk evolution models that are timely in the era of disk
surveys.

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 Associated Universities, Inc.

The authors thank the referee for a rigorous and constructive
report. B.T. thanks Min Fang for sharing data on Upper Sco and
Antonella Natta for useful discussions. This work was supported
by the Programme National PCMI of CNRS/INSU with INC/
INP cofunded by CEA and CNES. 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. R.A.
acknowledges funding from the Science & Technology Facilities
Council (STFC) through Consolidated Grant ST/W000857/1. L.
T. and K.Z. acknowledge the support of the NSF AAG grant
\#2205617. 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 \#3220495. 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
\#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 \#1241056.
N.T.K. acknowledges support provided by the Alexander von
Humboldt Foundation in the framework of the Sofja Kovalevs-
kaja 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. G.L. has
received funding from the European Union’s Horizon 2020
research and innovation program under the Marie-Sklodowska
Curie grant agreement No. 823823738 (DUSTBUSTERS) and
from PRIN-MUR 20228JPA3A. All figures were generated with
the PYTHON-based package MATPLOTLIB (J. D. Hunter 2007).
This research made use of Astropy,<fnref n="24" /> a community-developed
core Python package for Astronomy (Astropy Collaboration
et al. 2018).

Data Availability

The best-fit synthetic disk populations are available on the
AGE-PRO website at: http://agepro.das.uchile.cl/.

Facility: ALMA.

Software: astropy (Astropy Collaboration et al. 2018),
Diskpop (A. Somigliana et al. 2024).

<footnote n="24">²⁴ http://www.astropy.org</footnote>
# Appendix
## Constraints on $\langle \dot{M}_{\text{PEW}} \rangle$ and $\langle M_0 \rangle$ for Turbulence-driven Populations

In this appendix, we present the constraints on the median mass-loss rate, $\langle \dot{M}_{\text{PEW}} \rangle$, for various initial median disk sizes, $\langle R_0 \rangle$, obtained by matching the disk fractions in the turbulence-driven scenario. Figure 15 complements Figure 9 in the main text, with the latter corresponding to a single value of $\langle R_0 \rangle$. Following Figure 9, we also show the predicted median accretion rate at the initial time and at 2 Myr.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Summary of the constraints on turbulence-driven disk evolution obtained from the disk dispersal time. The fitted values of $\langle \dot{M}_{\text{PEW}} \rangle$ are plotted as a function of $\langle t_{\text{v,0}} \rangle$ for various values of $R_0$ (top row). The median accretion rates at $t=0$ and $t=2$ Myr are shown in the middle and bottom panels, respectively.</figcaption></figure>
ORCID iDs

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

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

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

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

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

Alice Somigliana https://orcid.org/0000-0003-2090-2928

Richard Alexander https://orcid.org/0000-0001-6410-2899

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

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

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

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

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

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

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

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

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

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

Anibal Sierra https://orcid.org/0000-0002-5991-8073

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

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

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

Estephani TorresVillanueva https://orcid.org/0000-0001-9961-8203

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

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

Claudia Toci https://orcid.org/0000-0003-2090-2928

Leonardo Testi https://orcid.org/0000-0003-1859-3070

Giuseppe Lodato https://orcid.org/0000-0002-2357-7692

References

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

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

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

Alexander, R. 2008, NewAR, 52, 60

Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VII, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 475

Alexander, R., Rosotti, G., Armitage, P. J., et al. 2023, MNRAS, 524, 3948

Alexander, R. D., & Armitage, P. J. 2009, ApJ, 704, 989

Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006, MNRAS, 369, 229

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

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

Andrews, S. M., & Williams, J. P. 2005, ApJ, 631, 1134

Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240

Armitage, P. J., Simon, J. B., & Martin, R. G. 2013, ApJL, 778, L14

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

Bai, X.-N. 2016, ApJ, 821, 80

Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144

Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76

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

Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Booth, A. S., Tabone, B., Ilee, J. D., et al. 2021, ApJS, 257, 16

Casse, F., & Ferreira, J. 2000, A&A, 361, 1178

Cazzoletti, P., Manara, C. F., Liu, H. B., et al. 2019, A&A, 626, A11

Clarke, C. J. 2007, MNRAS, 376, 1350

Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485

Coleman, G. A. L., & Haworth, T. J. 2022, MNRAS, 514, 2315

de Valon, A., Dougados, C., Cabrit, S., et al. 2020, A&A, 634, L12

Delage, T. N., Okuzumi, S., Flock, M., Pinilla, P., & Dzyurkevich, N. 2022, A&A, 658, A97

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

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

Emsenhuber, A., Burn, R., Weder, J., et al. 2023, A&A, 673, A78

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

Fang, M., Pascucci, I., Edwards, S., et al. 2018, ApJ, 868, 28

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

Fedele, D., van den Ancker, M. E., Henning, T., Jayawardhana, R., & Oliveira, J. M. 2010, A&A, 510, A72

Fernandes, R. B., Mulders, G. D., Pascucci, I., Mordasini, C., & Emsenhuber, A. 2019, ApJ, 874, 81

Ferreira, J. 1997, A&A, 319, 340

Ferreira, J., Dougados, C., & Cabrit, S. 2006, A&A, 453, 785

Flaischlen, S., Preibisch, T., Manara, C. F., & Ercolano, B. 2021, A&A, 648, A121

Gammie, C. F. 1996, ApJ, 457, 355

Gárate, M., Delage, T. N., Stadler, J., et al. 2021, A&A, 655, A18

Gárate, M., Pinilla, P., Haworth, T. J., & Facchini, S. 2024, A&A, 681, A84

Gorti, U., & Hollenbach, J. 2009, ApJ, 690, 1539

Hartmann, L., D'Alessio, P., Calvet, N., & Muzerolle, J. 2006, ApJ, 648, 484

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

Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654

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

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

Kimmig, C. N., Dullemond, C. P., & Kley, W. 2020, A&A, 633, A4

Komaki, A., Nakatani, R., & Yoshida, N. 2021, ApJ, 910, 51

Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271

Kraus, A. L., Ireland, M. J., Hillenbrand, L. A., & Martinache, F. 2012, ApJ, 745, 19

Kunitomo, M., Suzuki, T. K., & Inutsuka, S.-i. 2020, MNRAS, 492, 3849

Laos, S., Wisniewski, J. P., Kuchner, M. J., et al. 2022, ApJ, 935, 111

Lega, E., Morbidelli, A., Nelson, R. P., et al. 2022, A&A, 658, A32

Lesur, G. R. J., Flock, M., Ercolano, B., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 465

Lesur, G. R. J., Baghdadi, S., Wafflard-Fernandez, G., et al. 2023, A&A, 677, A9

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

Louvet, F., Dougados, C., Cabrit, S., et al. 2018, A&A, 618, A120

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

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

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

Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3

Manara, C. F., Tazzari, M., Long, F., et al. 2019, A&A, 628, A95

Manara, C. F., Testi, L., Herczeg, G. J., et al. 2017, A&A, 604, A127

Martel, É., & Lesur, G. 2022, A&A, 667, A17

Martire, P., Longarini, C., Lodato, G., et al. 2024, A&A, 686, A9

Maucó, K., Manara, C. F., Ansdell, M., et al. 2023, A&A, 679, A82

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. (Francisco, CA: ASP), 501

Miret-Roig, N., Galli, P. A. B., Olivares, J., et al. 2022, A&A, 667, A163

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

Morishima, R. 2012, MNRAS, 420, 2851

Mulders, G. D., Pascucci, I., Manara, C. F., et al. 2017, ApJ, 847, 31

Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018, ApJ, 857, 57

Nazari, P., Tabone, B., Ahmadi, A., et al. 2024, A&A, 686, A201

Ogihara, M., Kokubo, E., Suzuki, T. K., & Morbidelli, A. 2018, A&A, 615, A63

Ogihara, M., Morbidelli, A., & Guillot, T. 2015, A&A, 584, L1

Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880

Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415

Pascucci, I., Banzatti, A., Gorti, U., et al. 2020, ApJ, 903, 78

Pascucci, I., Beck, T. L., Cabrit, S., et al. 2025, NatAs, 9, 81
Pascucci, I., Cabrit, S., Edwards, S., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 567

Picogna, G., Ercolano, B., & Espaillat, C. C. 2021, MNRAS, 508, 3611

Preibisch, T., Kim, Y.-C., Favata, F., et al. 2005, ApJS, 160, 401

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

Ribas, Á., Bouy, H., & Merín, B. 2015, A&A, 576, A52

Ribas, Á., Merín, B., Bouy, H., & Maud, L. T. 2014, A&A, 561, A54

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

Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95

Rodenkirch, P. J., Klahr, H., Fendt, C., & Dullemond, C. P. 2020, A&A, 633, A21

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

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

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

Sellek, A. D., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 498, 2845

Sellek, A. D., Grassi, T., Picogna, G., et al. 2024, A&A, 690, A296

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33

Simon, M. N., Pascucci, I., Edwards, S., et al. 2016, ApJ, 831, 169

Somigliana, A., Toci, C., Lodato, G., Rosotti, G., & Manara, C. F. 2020, MNRAS, 492, 1120

Somigliana, A., Testi, L., Rosotti, G., et al. 2024, A&A, 689, A285

Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2018, MNRAS, 477, 1239

Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74

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

Tabone, B., Cabrit, S., Bianchi, E., et al. 2017, A&A, 607, L6

Tabone, B., Cabrit, S., Pineau des Forêts, G., et al. 2020, A&A, 640, A82

Tabone, B., Rosotti, G. P., Lodato, G., et al. 2022b, MNRAS, 512, L74

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

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

Toci, C., Rosotti, G., Lodato, G., Testi, L., & Trapman, L. 2021, MNRAS, 507, 818

Tong, S., Alexander, R., & Rosotti, G. 2024, MNRAS, 533, 1211

Trapman, L., Zhang, K., Rosotti, G. P., et al. 2025a, ApJ, 989, 5

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

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., Zhang, K., van't Hoff, M. L. R., Hogerheijde, M. R., & Bergin, E. A. 2022, ApJL, 926, L2

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

Venuti, L., Bouvier, J., Flaccomio, E., et al. 2014, A&A, 570, A82

Wafflard-Fernandez, G., & Lesur, G. 2023, A&A, 677, A70

Wang, L., Bai, X.-N., & Goodman, J. 2019, ApJ, 874, 90

Wang, L., & Goodman, J. 2017, ApJ, 847, 11

Weder, J., Mordasini, C., & Emsenhuber, A. 2023, A&A, 674, A165

Winter, A. J., & Haworth, T. J. 2022, EPJP, 137, 1132

Yap, T. E., & Batygin, K. 2024, Icar, 417, 116085

Yoshida, T. C., Nomura, H., Tsukagoshi, T., Furuya, K., & Ueda, T. 2022, ApJL, 937, L14

Zagaria, F., Rosotti, G. P., Clarke, C. J., & Tabone, B. 2022, MNRAS, 514, 1088

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