Testing protoplanetary disc evolution with CO ﬂuxes A proof of concept in Lupus and Upper Sco

,


Introduction
Over the last decades two disc evolution models, viscous theory (Lynden-Bell & Pringle 1974;Shakura & Sunyaev 1973) and the magnetohydrodynamic wind (MHD-wind) scenario (Blandford & Payne 1982), have been proposed (Manara et al. 2022).According to the viscous evolution model, the disc angular momentum is conserved and redistributed by turbulence: while a small fraction of the disc mass moves to larger sizes, the bulk is accreted.Instead, in the MHD-wind scenario, powerful magnetothermal winds are launched from the disc, allowing accretion to efficiently remove angular momentum.In addition to these mechanisms, thermal winds, not instrumental in driving the accretion process, are thought to play a key role in the disc dispersal phase (Pascucci et al. 2022), complicating the picture.Discriminating between these two scenarios requires large surveys targeting populations of discs of different ages in order to compare models and data in a statistical sense.In recent years, the Atacama Large Millimeter/submillimeter Array (ALMA) observed several nearby star-forming regions (SFRs) (e.g.Ansdell et al. 2016Ansdell et al. , 2018;;Pascucci et al. 2016;Barenfeld et al. 2016;Cieza et al. 2019;Cazzoletti et al. 2019) at moderate resolution (0.25 to 0.50 arcsec) and sensitivity (0.1 to 0.4 M C ), measuring fluxes and sizes for tens of discs (Manara et al. 2022;Miotello et al. 2022) from dust and CO rotational transitions.
Disc sizes have been particularly useful to study disc evolution because of the different trends predicted by models: while viscous discs are expected to get larger with time, in the MHDwind scenario discs either remain the same or shrink (Manara et al. 2022).In the case of dust, Rosotti et al. (2019) predicted that the disc radius (enclosing 95% of the total dust flux) expands with time in viscous models.However, if present, this behaviour can only be detected in very deep surveys, with a sensitivity that is fifty times better than in the available data.This sensitivity can be reached with roughly five hours on-source at an intermediate resolution (0.6 to 0.7 arcsec), which would be prohibitive for any future ALMA survey targeting hundreds of discs.Zagaria et al. (2022) extended the work of Rosotti et al. (2019), showing that this same factor of fifty is needed to distinguish between viscous and MHD-wind evolution.Furthermore, a direct comparison between models and data is made more difficult by the presence of substructures (Toci et al. 2021;Zormpas et al. 2022;Zagaria et al. 2022) since the observed sizes may trace the effects of discplanet interactions rather than disc evolution.
In the case of gas, following up on the early work of Najita & Bergin (2018), Trapman et al. (2020) used complex thermochemical models to show that small discs with low viscosities can explain most of the observationally inferred disc sizes in Lupus, but they spread too much to reproduce more compact discs in Upper Sco.MHD-wind models, instead, are broadly consistent with the gas disc sizes measured in both SFRs (Trapman et al. 2022).However, this comparison is affected by two main uncertainties: the small samples, particularly at the age of Upper Sco (Barenfeld et al. 2017), and the amount of carbon deple-tion.When the carbon abundance falls below x CO « 10 ´6, Trapman et al. (2022) showed that discs observed with low sensitivity could look up to 70% smaller or be unresolved.To mitigate this problem, integration times of one hour per source would be needed, which is challenging for large surveys.
However, targeting disc sizes is not the only possible strategy to study disc evolution.Here we introduce an alternative method based on 12 CO fluxes.Assuming that 12 CO emission is optically thick, CO fluxes scale as the disc surface area (i.e. the radius squared), suggesting that modelling fluxes is an indirect way of studying sizes since they would trace the same physicochemical processes in the disc.This assumption is supported by both models (Trapman et al. 2019(Trapman et al. , 2020(Trapman et al. , 2022;;Miotello et al. 2021) and by the data.For example, Long et al. (2022) showed that the observationally inferred CO fluxes and sizes correlate well, with R CO 9F 0.52˘0.05CO (see also Sanchis et al. 2021).Observing fluxes instead of sizes is less time consuming, firstly, because one would aim to detect, but not necessarily resolve, a target, and secondly, because there would be no need for very deep surveys targeting the faint outer disc regions that contribute marginally to the disc brightness.In this Letter we introduce a simple semi-analytical prescription to compute 12 CO disc fluxes under the optically thick assumption.We benchmark this prescription against a grid of full radiative transfer simulations and show that they agree, on average, within a factor of three.Then, as a proof of concept, we compare these models with Lupus and Upper Sco data, highlighting the main limitations of the available datasets and the foreseen improvements with future dedicated surveys.
This Letter is organised as follows.In Sect. 2 we introduce our semi-analytical method.In Sect. 3 we run a disc population synthesis model and compare viscous and MHD-wind predictions with Lupus and Upper Sco data.Our results are discussed in Sect.4, and in Sect. 5 we draw our conclusions.The code developed for this work is publicly available on github.

Methods
Here we summarise our assumptions and final equation to compute CO fluxes (see Sect.A for the full derivation).
We considered 12 CO emission to be optically thick and in local thermodynamical equilibrium.Under these assumptions where R is the cylindrical disc radius, i the disc inclination, and d its distance from the observer.The brightness profile in Eq. 1 can be written as where B ν 0 is the black-body emission at temperature T and frequency ν 0 , and the exponential term gives the thermal broadening of the line (Rybicki & Lightman 1986).Here m CO is the 12 CO molecular mass, c the speed of light, and k B the Boltzmann constant.We adopted a power-law temperature profile with exponent ´0.5 and normalisation 87.5 K at 20 au, in agreement with the inferences of Law et al. (2021Law et al. ( , 2022a,b),b).Our disc inclination was fixed to the sky-averaged value of cos i " π{4.We adopted R in " 10 ´2 au and R out " R CO , the radius where the gas surface density equals the column density, N CO " 5ˆ10 15 cm ´2 (a density slightly larger than the standard result of van Dishoeck & Black 1988), where 12 CO is not efficiently selfshielded against photodissociation and is quickly removed from the gas phase.To compute the gas surface density corresponding to N CO , we assumed the same carbon abundance of the diffuse ISM, x CO " 10 ´4.Although it is rather crude, this method is supported by the work of Trapman et al. (2019), who showed that R CO encloses all of the CO emission and is in good agreement with the results of complex thermochemical models.To test our method, we benchmarked our sizes and fluxes against the results of the thermochemical models of Miotello et al. (2016) and Trapman et al. (2020Trapman et al. ( , 2022)), run using the code DALI (Bruderer et al. 2012;Bruderer 2013).The results of this exercise are extensively discussed in Sect.B, where we show that our face-on fluxes underestimate DALI ones by a factor of three.

Population synthesis
We give a proof of concept of this new method comparing our semi-analytical predictions with the available Lupus (age À 3 Myr, Luhman & Esplin 2020) and Upper Sco (age 5 to 10 Myr, Luhman 2020) data.A quick description of the datasets can be found in Sect.C.Here we note that even with the limited data available, working with fluxes instead of sizes increases the samples by a factor of 1.3 (48 instead of 36 sources) in Lupus and by a substantial factor of 3.6 (32 instead of 9 sources) in Upper Sco.For this comparison we relied on a disc population synthesis approach: we prescribed a set of initial conditions and evolved our models in the viscous or MHD-wind case under the assumption that these two SFRs can be regarded as subsequent evolutionary stages of the same population (i.e. they have the same initial conditions).Unfortunately, these initial conditions are either unknown or very uncertain (e.g. they were inferred neglecting any contribution of dust to disc evolution; Lodato et al. 2017;Tabone et al. 2022b).Future more accurate distributions will allow more reliable comparisons between evolutionary models and data.

Viscous case
We used the Lynden-Bell & Pringle (1974) analytical solution to compute the CO radius.In this case, the surface density at a given time is a function of the viscous timescale (t acc ), initial disc mass (M 0 ), and scale radius (R 0 ).We assumed the viscous timescale to be distributed as logpt acc {yrq " Np5.8, 1.0q, where the notation Npµ, σq stands for a Gaussian distribution with mean µ and variance σ 2 .This distribution was inferred by Lodato et al. (2017) fitting the Lupus data in the 9 M acc ´Mdisc plane, under the assumption that viscosity is an increasing function of the disc radius with exponent γ " 1.5.For the initial disc mass distribution we considered logpM 0 {M d q " Np´2.7, 0.7q, similarly to Lodato et al. (2017).Even though young discs are known to be small (Maury et al. 2019;Maret et al. 2020;Tobin et al. 2020), their initial disc size distribution is not well constrained.To take into account possible envelope contributions, we adopted the best fit R disc « R 0 distribution of 25 Class 0 objects in Orion (VANDAM, Tobin et al. 2020) based on the radiative transfer models of Sheehan et al. (2022), under the assumption that gas and dust are co-located at such young ages: logpR 0 {auq " Np1.55, 0.4q.Finally, we assumed an age of logpt{yrq " Np5.9, 0.3q for Lupus (corresponding to our choice of t acc ; see Lodato et al. 2017) and 7  We checked the α SS (Shakura & Sunyaev 1973) distribution associated with our initial conditions (for a M d star): Here we considered h to be a power law with exponent 0.25 and normalisation h 0 " 0.1 at 10 au, in line with the results of Zhang et al. (2021).Our choices of t acc and R 0 give a distribution of log α SS « Np´2.42, 1.06q.
Our results are displayed in Fig. 1.Measured fluxes are shown as purple and orange patches for Lupus and Upper Sco; the survival functions and their 1σ spread were computed using the Kaplan-Meier estimator for left-censored datasets (see Sect.C).The survival functions for N discs " 3000 models are plotted as solid lines of the same colours.Our results under standard assumptions (see Sect. 2) are presented in the left panel.To get a better insight into these flux distributions, we follow the evolution of the median disc (i.e. the disc whose initial conditions are the median of our assumed distributions).This disc spreads viscously, getting bigger and brighter, until an inversion time t inv (Eq.12 of Toci et al. 2023).Then, the part of the disc that is viscously expanding falls below the CO photodissociation threshold, making the disc smaller and fainter.Our Upper Sco models are fainter than Lupus models because more discs (particularly those with larger R 0 , lower M 0 , and shorter t acc ) went past their inversion time.Nevertheless, our models are Á 10 times brighter than the data.To reconcile models and observations we introduced a column density fudge factor ξ, that makes photodissociation more efficient: N CO Ñ N CO {ξ.A discussion on the possible physico-chemical interpretation of such a factor can be found in Sect. 4. Our results are displayed in the right panel of Fig. 1, for ξ visc Lup " 2.5 ˆ10 ´2 and ξ visc Sco " 10 ´3: these fudge factors are able to reconcile models and observations both at the age of Lupus and Upper Sco.However, for the faintest discs in Lupus and the brightest in Upper Sco a smaller (larger) correction factor would be required.This effect can also be due to warmer (colder) discs than our average temperature profile.We note that these fudge factors are not an artefact of our initial conditions; in other words, we found no sensible combination of the initial parameters able to viscously reproduce both Lupus and Upper Sco observations with ξ visc Lup " ξ visc Sco " 1.

MHD-wind case
We used the Tabone et al. (2022a) analytical solution with constant magnetic field strength (ω " 1) to compute the CO radius.This solution can reproduce both the disc fraction decay with time and the Lupus data in the 9 M acc ´Mdisc plane (Tabone et al. 2022b).In this case the surface density at a given time is a function of the accretion timescale (t acc ), the initial disc mass (M 0 ), the initial scale radius (R 0 ), and the lever-arm parameter (λ).Because the wind-driven prescription accounts for disc dispersal after a finite time, knowledge of the disc fraction distribution can be used to infer a distribution of t acc (see Eq. A.2 of Tabone et al. 2022b).Following Tabone et al. (2022b), our initial disc mass distribution is log-normal with 1.0 dex spread and centred on M 0 " 2 ˆ10 ´3 M d , with corresponding mass ejectionto-accretion ratio f M " 0.6.We chose the same initial radius distribution of the viscous case.Then, the lever-arm parameter distribution can be computed from R 0 and f M (see Tabone et al. 2022b, where we fixed the innermost wind launching radius to 1 au).The parameter λ is distributed with µ « 4.8 and σ « 0.95.Finally, we assumed an age of 2 Myr for Lupus (corresponding to our choice of M 0 ; see Tabone et al. 2022b) and 7.5 Myr for Upper Sco.
Our results are displayed in the left panel of Fig. 2, using the same symbols as Fig. 1.As noted above, these MHD-wind models have a finite lifetime which was chosen to reproduce the observed age dependence of the disc fraction in SFRs.Since Upper Sco is older than Lupus and their samples are of similar sizes (Sect.C), a larger number of initial discs is needed to reproduce the number of sources observed in the former region: N 0,Lup " 196 and N 0,Sco " 2490 (in Fig. 2 20 times more models are shown to better explore the initial conditions).In the MHD-wind case the median disc evolves slowly and its radius is constant with time, until t{t depl Á 95%, when fast dispersal takes place and the disc gets abruptly smaller and dimmer.Consequently, we expect our CO flux distributions to be mostly dependent on the initial disc size.A clear difference with the viscous case is that the brightest MHD-wind models are almost as luminous as the brightest Lupus and Upper Sco data, while for fainter and fainter discs the discrepancy between models and data progressively increases.As a consequence, a constant fudge factor cannot reconcile models and observations; a disc mass dependent correction would need to be invoked.We obtained similar results for ω " 0.5 and the initial parameters explored by Tabone et al. (2022b).

Discussion
In this Letter we assumed that the Lupus and Upper Sco disc populations could be considered a respectively younger and older evolutionary stage of the same initial disc population.Under this hypothesis, we ran disc population synthesis models from sensible initial conditions and compared our CO flux distributions with the data.
To match models and data, in the viscous case we introduced a column density fudge factor that increases the CO photodissociation efficiency.This factor can be interpreted as the outcome of some processes depleting CO in protoplanetary discs (Miotello et al. 2022).The low gas masses estimated from CO in Lupus and Chamaeleon I discs (Ansdell et al. 2016(Ansdell et al. , 2018;;Miotello et al. 2017;Long et al. 2017) indicate that protoplanetary discs are fainter than expected in CO.This is supported by the few direct measurements of disc masses based on HD rotational line transitions, that require CO to be depleted by factors between 5 and 200 (Bergin et al. 2013;Favre et al. 2013;Mc-Clure et al. 2016;Schwarz et al. 2016).Two main processes have been proposed to explain CO underabundance: (i) because of CO chemical (gas-or ice-phase) conversion or evolution, carbon would be sequestered into more complex species, like CO 2 or hydrocarbons, that can freeze-out onto grains at higher temperatures than CO (e.g.Bosman et al. 2018;Schwarz et al. 2018); (ii) CO freeze-out on dust and subsequent grain growth into larger bodies that no longer participate in gas-phase chemistry would lock carbon up in the disc midplane and transport it radially (e.g.Krijt et al. 2016Krijt et al. , 2018;;Powell et al. 2022).A combination of the two processes is most likely to take place (e.g.Booth et al. 2017;Krijt et al. 2020) and is needed to explain depletion factors of 100 on a timescale of 1 to 3 Myr inferred from the comparison between Class I and Class II discs (Zhang et al. 2020).The carbon depletion scenario is supported by the detection of C 2 H in several protoplanetary discs, which can be explained by carbon and oxygen depletion with C{O Á 1 (Bergin et al. 2016;Cleeves et al. 2018;Miotello et al. 2019;Bosman et al. 2021).This result is consistent with the carbon-to-oxygen ratio inferred from the available NpCSq{NpSOq upper limits (Semenov et al. 2018;Facchini et al. 2021;Le Gal et al. 2021).
The fudge factor we introduced in the viscous case to explain Lupus data is in line with the observationally inferred CO depletion factors of two orders of magnitude (Miotello et al. 2022).Notably, in addition to CO fluxes, these carbon-depleted models can also reproduce very well the 9 M acc ´Mdisc correlation (by construction) and the CO size-luminosity correlation (see Sect.D).
Here we caution that our ξ visc Lup should be interpreted as a population average rather than an absolute depletion factor; other factors (e.g. a spread in the stellar mass and luminosity, dust evolution, and different source inclinations) can effect CO depletion, either by changing disc chemistry or the observed luminosities.The larger correction needed in the case of Upper Sco, instead, could be explained in three different ways.Firstly, CO is more depleted in Upper Sco than in Lupus.This hypothesis is supported by several evolutionary models, where the carbon depletion factor increases with time (Krijt et al. 2020;Powell et al. 2022).Furthermore, Anderson et al. (2019) showed that CO abundances ď 10 ´6 are needed to explain the observed N 2 H ànd CO line fluxes of two Upper Sco discs.Secondly, other processes are affecting disc evolution in Upper Sco.Removing the less bound material from the disc outer regions, external photoevaporation halts viscous spreading, making discs smaller and fainter (e.g.Clarke 2007).While in Lupus the irradiation levels are expected to be low (Cleeves et al. 2016) and photoevaporation to be negligible (with the possible exception of large discs, Haworth et al. 2017), Trapman et al. (2020) argued that the level of irradiation in Upper Sco can be « 100 times higher and photoevaporation more efficient.Thirdly, the initial conditions are different.In this case, Upper Sco cannot be regarded as the subsequent evolutionary stage of Lupus and the two regions must be modelled separately (e.g. if t acc is shorter in Upper Sco than in Lupus, its disc fluxes will decrease faster).
In the MHD-wind case, models and data have different shapes; no fudge factors are needed to explain the brightest sources, but fainter models require some corrections.Even though Trapman et al. (2022) invoked carbon depletion to reconcile MHD-wind models and the observationally inferred disc sizes in Upper Sco, we note that some of the brightest Upper Sco discs in our dataset were not included in their work because they were not well resolved (see black-contour dots in Fig. C.1).In any case, considering our previous arguments on carbon depletion, our conclusion that MHD-wind models do not need lower CO column densities to match the brightest sources is puzzling.A possible explanation is that our initial disc size distribution is not suitable; in the MHD-wind models of Tabone et al. (2022a,b), R 0 is constant with time, and thus it must match the observed disc sizes in Lupus and Upper Sco.For Lupus, assuming as initial disc size distribution logpR 0 {auq " Np1.75, 0.4q, whose trailing edge agrees with the observationally inferred R CO,68 distribution (Sanchis et al. 2021), MHD-wind models require a fudge factor of ξ wind Lup " 10 ´2, similar to the viscous one, to match the data.This is shown in the right panel of Fig. 2, where models with a reduced CO column density are plotted with a dotted line.As in the viscous case, these models can also reproduce very well the 9 M acc ´Mdisc correlation (by construction) and the CO size-luminosity correlation (see Sect.D).Instead, the very few constraints on the disc size distribution in Upper Sco are not in contrast with our assumption for the initial disc distribution in Sect.3. A possible explanation would be that Upper Sco discs were born more compact than Lupus discs (Barenfeld et al. 2016(Barenfeld et al. , 2017;;Miotello et al. 2021) or became smaller due to environmental effects such as photoevaporation (Trapman et al. 2020).Better data are needed to draw robust conclusions.

Conclusions
In this Letter we introduced a new method to study protoplanetary disc evolution based on 12 CO fluxes.Assuming optically thick emission, we built a semi-analytical model to compute disc fluxes; the results agree well (within a factor of three) with those of more time-expensive thermochemical models (DALI).Then, we simulated families of discs, evolving from the same initial conditions either under the effect of viscosity or MHD winds, and compared their fluxes with Lupus and Upper Sco data.Using fluxes instead of sizes increases our observational samples by a factor of 1.3 in Lupus and 3.6 in Upper Sco, allowing for a more robust comparisons between models and data.
In the viscous case, our models were brighter than the data.To match the observations, we introduced different column density fudge factors (ξ visc Lup " 2.5 ˆ10 ´2 and ξ visc Sco " 10 ´3) that can be explained by carbon depletion (Lupus), and the effects of thermal winds or different initial conditions (Upper Sco).In the MHD-wind case our models matched the brightest discs in Lupus and Upper Sco, but mass-dependent factors were needed to reproduce the fainter sources.In the case of Lupus, when larger initial disc sizes (compatible with the observed distribution) were prescribed, a constant factor (ξ wind Lup " 10 ´2, comparable with ξ visc Lup ) was needed to reproduce the observed fluxes.Unfortunately, our interpretation of the data is limited by the uncertainties on the initial conditions and the amount of carbon depletion.Nevertheless, our proof of concept shows the usefulness of CO fluxes to study disc evolution.Measuring fluxes instead of sizes is less time-consuming.Additionally, fluxes could be the only accessible observable in farther SFRs.Thanks to forthcoming surveys that will target tens of discs at limited resolution and with the potential to inform us about carbon depletion, we will be able to obtain the most knowledge about disc evolution from this new flux-oriented approach.
.5 Myr for Upper Sco.

Fig. 1 .
Fig. 1.Comparison of the data (patches) and viscous model (solid lines) survival functions at the age of Lupus (purple) and Upper Sco (orange).Left panel: Standard assumptions.Right panel: Reduced gas column density.Fudge factors (ξ visc , top right corner) are needed to match the data.

Fig. 2 .
Fig. 2. Comparison of the data (patches) and MHD-wind model (solid lines) survival functions at the age of Lupus (purple) and Upper Sco (orange).Left panel: Standard assumptions.A mass-dependent depletion factor needs to be invoked to match the data.Right panel: Larger initial disc size distribution in Lupus.A constant fudge factor can reproduce the data (dotted line for ξ wind Lup " 10 ´2).