© The Authors 2024

1 Introduction

The total gas mass and gas density and temperature distributions are the fundamental parameters with which to understand the evolution of protoplanetary disks and the planet formation process (e.g., Mordasini et al. 2012). Probing the disk gas structure is a challenging task, as most of the mass is carried by H₂ molecules, which do not emit under the physical conditions of the outer disk regions probed by the far-infrared and radio observations. Therefore, the gas disk thermal and density distributions are studied through the emission of other, much less abundant molecules, in combination with theoretical models. The most common gas tracer is the CO molecule, as it is relatively abundant and emits at cold temperatures (~ 10−50 K), typical of the gas in the outer disk regions (e.g., Aikawa et al. 2002; Williams & Best 2014; Zhang et al. 2021). A variety of the CO isotopologues with abundances that can differ by factors of ≳70−150000, and, consequently, with emission lines of very different optical depths, enable one to trace the disk gas through a broad range of vertical heights in the molecular layer. CO molecules are the key carriers of the elemental carbon and oxygen, have low freeze-out temperatures, and are crucial to the chemical processes of gas-phase carbon and oxygen in the outer disk leading to the formation of more complex organic molecules (Walsh et al. 2014; Favre et al. 2018; Lee et al. 2019).

CO is a chemically stable molecule with a simple and well-known chemistry. The CO rotational transitions come from the so-called disk molecular layer, which extends above the disk cold midplane where it is removed from the gas phase by freeze-out, until the photodissociation region where CO is destroyed by the external or stellar radiation field (Aikawa et al. 2002; Dartois et al. 2003; Molyarova et al. 2017). Consequently, the temperature of the CO-emitting gas in disks should be higher than about 20 K; otherwise, CO molecules would be removed from the gas phase by freeze-out (Harsono et al. 2015; Schwarz et al. 2016; Pinte et al. 2018).

However, some recent observational studies have questioned the reliability of CO as a gas tracer in disks due to stronger CO depletion being predicted by contemporary models, leading to a much lower CO/H₂ ratio than the canonical interstellar medium (ISM) value of ~10⁻⁴ (e.g., Aikawa et al. 1997; Favre et al. 2013; Miotello et al. 2018; Schwarz et al. 2021). Recent surveys in the Chamaeleon and Lupus star formation regions found that weak CO emission could be common in protoplanetary disks (Ansdell et al. 2016; Long et al. 2017; Miotello et al. 2017). These studies indicate that the CO abundance can be lower by up to two orders of magnitude than the ISM value of ~10⁻⁴.

There are several processes that can affect the CO abundance in the gas phase, making the conversion from CO abundance to H₂ abundance uncertain. First, CO molecules can be frozen out in the cold disk midplane, or photodissociated in the upper disk layers by the far-ultraviolet (FUV) photons from the central star or by external FUV radiation (Miotello et al. 2016). Second, chemical processes can alter the CO abundance, for example by transforming it into other molecules such as CO₂ and complex organics (e.g., Bruderer et al. 2012; Bruderer 2013; Schwarz et al. 2018). The photodissociation or freeze-out results in overall lower C/H and O/H abundances in the gas, which are important parameters for disk chemistry and which can hence affect the abundances of other molecular disk tracers. Constraining the total CO depletion factor from the line observations is a challenging task, as it requires a good understanding of the underlying chemistry as well as dust and gas dynamics and properties (Krijt et al. 2018, 2020), but also optical depth effects and line excitation conditions, since even rare CO isotopologues can be optically thick in the dense inner disk within the CO snowline region.

Isotopologue-selective processes need to also be considered, such as self-shielding against photodissociation or low-temperature ¹²C/¹³C-fractionation (e.g., Woods & Willacy 2009; Visser et al. 2009). Disk gas mass estimates based on the CO isotopologue lines could be underestimated by up to two orders of magnitude if these processes are not fully taken into account (Miotello et al. 2023). These processes are accounted for in many contemporary physical-chemical models, but not all (e.g., Miotello et al. 2014, 2016; Ruaud et al. 2022). Other recent studies compared Class I and Class II disks and found that the CO gas-phase abundances decrease rapidly in older disks (e.g., Zhang et al. 2020). Thus, to use the CO emission as a reliable proxy of the disk gas masses, we need to employ more feasible disk physical-chemical models, as well as parameter space sampling algorithms.

In this paper, we model the millimeter emission of the optically thick and thin CO (2−1) isotopologue lines observed in several Class II T Tauri disks in the Taurus star formation region as part of the MPG-IRAM large guaranteed time project PRODIGE (PROtostars to Disks: Global Evolution; PIs: P. Caselli & Th. Henning); see Sect. 2. Previous works using the PRODIGE data are the studies of the Class 0/I protostars and streamers in the Perseus region (Hsieh et al. 2023; Valdivia-Mena et al. 2022). An overview of the CO observations for all targeted Class II sources and the first analysis and modeling of these data are presented in Semenov et al. (2024). To fit the CO NOEMA data, we developed a novel disk fitting tool, DiskCheF¹ (Disk Chemical Fitter), which uses a parametric disk physical model, ML-accelerated gas-grain chemistry, and line radiative transfer based on the RADMC-3D code (Sect. 3). Applying DiskCheF to our CO data, we infer the best-fit disk physical structures and the corresponding gaseous masses (Sect. 4). We discuss the results and limitations of our DiskCheF fitting and compare it with previous similar studies of disk masses in Sect. 5.

Recently, an ALMA study of molecular line emission in five protoplanetary disks was carried out by the MAPS collaboration ((PI: K. Oeberg, Öberg et al. 2021, see the following papers), with a different science goal than the one presented in this work. While the five MAPS sources comprise both T Tauri and Herbig Ae stars from different star-forming regions, the disks in the PRODIGE survey are all T Tauri systems and located in the same Taurus star-forming region.

2 Observations

We chose to perform the CO line fitting on six of the eight Class II disks observed in our sample; namely, CI, CY, DL, DM, DN, and IQ Tau (see Table 1). The selected systems are isolated T Tauri stars (M_⋆ < 1 M_⊙) with ages between 1 and 4 Myr (see Semenov et al. 2024 and references within), surrounded by extended disks (R_co ≳ 200 au) without strong substructures visible in the ALMA continuum observations (Long et al. 2018, 2019). The two systems removed from the modeling are DG Tau and UZ Tau E. The DG Tau system is young and still embedded in an extended envelope with a jet and an outflow, and would require more dedicated efforts and a detailed disk-envelope-outflow model to reproduce complex CO data. The UZ Tau E system is excluded because it is a part of a binary system, and shows a large-scale gaseous arc-like structure connecting the two circumstellar disks, likely caused by their gravitational interactions.

In this work, we focus on the CO, ¹³CO, and C¹⁸O J = 2−1 emission lines detected in the framework of the PRODIGE program. The spectral setup, observing details, data reduction, self-calibration, and imaging are described in the Class II PRODIGE paper by Semenov et al. (2024). Briefly, NOEMA/PolyFiX observations with ten antennas in C configuration (0.9″ angular resolution) were taken in Band 3 at 1.3 mm in 2020, using 62.5 kHz spectral chunks to spectrally resolve specific emission lines. We targeted the lines of ¹²CO, ¹³CO, C¹⁸O, para-H₂CO, DCO⁺, DCN, DNC, ¹³CN, cyclic C₃H₂, C₂D, HC₃N, N₂D⁺, and other rare species. Each source was observed independently (no track-sharing), using one to three tracks. The average on-source integration time was 5.5 h per disk, and the resulting 1σ rms noise is 3.7−9.2 mJy at 0.3 km s⁻¹ resolution.

A preliminary analysis of the combined CO (2−1) data in Semenov et al. (2024) shows that a foreground cloud partially obscures the optically thick ¹²CO (2−1) emission in the 4.56 km⁻¹ velocity range in the five disks, except for DM Tau, DN Tau, and IQ Tau. The ¹²CO (2−1) peak brightness temperatures range between 6.5 and 16.0 K, while for more optically thin ¹³CO (2−1) and C¹⁸O (2−1) these peak brightness temperatures are 1.2−7.5 K and 0.2−1.8 K, respectively. The moderate spatial resolution of our NOEMA observations (100−150 au) leads to smooth CO radial intensity profiles as reconstructed by the Kepler deprojection and azimuthal averaging, with no evidence of substructures. The apparent inner emission “holes” in the CO intensities at r ≲ 100−150 au are due to beam dilution. In Fig. 1, we show the moment zero emission maps of our selected sources.

Since our disks have many similar properties (e.g., low stellar luminosities, ages, etc.) and have been observed and imaged in the same way, our goal is to present a general model that can be homogeneously applied to the fitting of the selected disks in our sample. This homogenous modeling approach should allow us to better understand the similarities and differences within this class of disks, and how the differences in their structures affect their chemical composition. In Table 1 we present an overview of the stellar and disk properties as adopted from Semenov et al. (2024).

3 Disk model

To model the disk structure and chemistry, we developed the DiskCheF framework, of which this work is the first application (Smirnov-Pinchukov et al. 2022). This tool was developed to simulate molecular line data using a modular disk physical-chemical framework and can be applied either to fit disk interferometric data or as a standalone disk forward modeling tool. The code philosophy is to establish a series of customizable building blocks so that each code block can be modified without affecting the functionality of the other code structures. A resulting disk model can be built by choosing any combination of the individual disk’s physical, chemical, and radiative transfer models. The preset models can be easily subclassed to add custom functionalities when needed. The ^*Base objects described below are templates showing the information needed by each step of a DiskCheF modeling (namely, the disk physical model, chemical model, radiative transfer model, and fitter). These objects cannot be directly used to create a disk model, so the user has to either use them as a template for the desired theoretical disk model or choose the corresponding default models already implemented in DiskCheF.

The starting point of DiskCheF is to create a disk physical model (physics.PhysicsBase), with information about the gas and dust density and temperature distributions. This model is then used as input for the chemistry model, (chemistry.ChemistryBase), which computes the abundance of the chemical species using the disk physical structure provided by the previous disk physics module. DiskCheF also uses the radiative transfer modeling tool maps.RadMCBase to produce channel and moment maps (maps.RadMCRT) and dust continuum emission maps (maps.RadMCRTImage), and to calculate the disk temperature (maps.RadMCTherm). The workflow of the DiskCheF modeling is shown in Fig. 2. In the next sections, we describe in more detail the individual models implemented within the DiskCheF framework to fit the CO data.

3.1 Disk physical model

To model the disk’s physical and thermal structure, we adopted the Lynden-Bell & Pringle (1974) prescription. This model is implemented in the physics.WB100auWithSmoothInnerGap module. In this model, the radial surface density distribution of the gas is given by the Lynden-Bell & Pringle self-similar profile with a smooth inner gap: $$\text{Σ}(r)={\text{Σ}}_0{\left(\frac{r}{r_{\text{c}}}\right)}^{-\gamma }\text{ exp }\left[-{\left(\frac{r}{r_{\text{c}}}\right)}^{2-\gamma }\right]\text{ exp }\left[-{\left(\frac{r}{r_{in}}\right)}^{\gamma -2}\right]\text{, }$$(1) $${\text{Σ}}_0=(2-\gamma )\frac{M_{\text{disk }}}{2\pi r_{\text{c}}^2}\text{ exp }{\left(\frac{r_{in}}{r_{\text{c}}}\right)}^{2-\gamma },$$(2)

where γ is the tapering factor, r_c the tapering radius, r_in the inner disk radius, and M_disk the total disk mass. The vertical profile is given by integrating the vertical hydrostatic equilibrium equation: $$\frac{\partial \text{ ln }\rho }{\partial ɀ}=-\left[\left(\frac{G{M}_{\text{star }}ɀ}{{\left(r^2+z^2\right)}^{3/2}}\right)\left(\frac{\mu m_{\text{H}}}{kT}\right)+\frac{\partial \text{ ln }T}{\partial ɀ}\right].$$(3)

The midplane and atmosphere temperature of a disk is given by the radial power law distribution: $$T_{\text{mid/atm }}=T_{\text{mid/atm, }100}{\left(\frac{r}{100\text{\hspace{0.17em}au}}\right)}^{-q_{\text{mid/alm }}},$$(4)

with T_mid/atm,₁₀₀ being the temperature at a radius of 100 au and q_mid/atm the exponent of the distribution. The value of the temperature exponent has been measured from the CO data for disks in the Taurus region, and it is usually assumed that a slope of 0.55 for both the midplane and atmosphere temperature profile provides a good representation of the outer disk temperature distribution (Williams & Best 2014). The transition between the atmosphere and midplane temperature is parameterized with a sine function $$T_{\text{WB}}(r,z)=\{\begin{array}{ll}{T}_{\text{mid }}+\left(T_{\text{atm }}-T_{\text{mid }}\right){\left[\text{\hspace{0.17em}sin\hspace{0.17em}}\left(\frac{\pi z}{2z_q}\right)\right]}^{2\delta }\hfill & \text{ if }z<z_{q,}\hfill \\ T_{\text{atm }}\hfill & \text{ if }z\ge z_q\hfill \end{array}$$(5)

where δ describes the steepness of the profile and ɀ_q is the height over the disk atmosphere where the disk reaches the atmospheric temperature. Williams & Best (2014) explored the impact of these parameters and found that these parameters do not significantly affect the line luminosity. In the literature it is usually assumed that δ = 2 and ɀ_q = 4 H_p, where H_p is the pressure scale height: $$H_p=\sqrt{\frac{\kappa T_{\text{mid }}{r}^3}{G{M}_{\text{star }}\mu m_{\text{H}}}},$$(6)

where κ is the Boltzmann constant, G the gravitational constant, µ = 2.3 the mean molecular weight of the gas, and m_H the mass of atomic hydrogen. Moreover, the typical interstellar radiation field in star-forming regions prevents the disk from reaching unrealistically low temperatures. Following the prescription in Tazzari et al. (2021), we set a threshold of T_floor = 7 K, and we used an effective temperature of $$T^4=T_{\text{WB}}^4+T_{\text{floor }}^4.$$(7)

Fig. 1 Moment zero maps of the CO, 13CO and C18O (2−1) emission for the disks chosen for the DiskCheF fitting.

Fig. 2 Outline of the steps in a DiskCheF model. For each step, we indicate which template can be used to create a new custom model.

3.2 Chemical model

The main goal of DiskCheF is to fit disk models to the multi-line and/or multi-species interferometric data. When the data are gas emission lines, we need a chemical model that is fast enough yet feasible enough to be part of a fitting routine. A typical simplified network for the chemical evolution of a protoplanetary disk runs for a few seconds on a single CPU. However, we discussed in Sect. 1 that we have to account for all of the processes affecting CO abundances, such as self-shielding or photodissociation. The neglect of these processes can underestimate the disk mass by one or two orders of magnitude. Therefore, to fit the CO line emission, we need a more robust chemical network. The computational time for such networks ranges from tens of minutes to days per disk model. Since our fitting routine computes on the order of 10⁵ disk models to converge to a best-fit model, it adds up to unfeasible computational times.

One solution is to apply machine learning techniques to the results of time-dependent chemical kinetics models. These techniques find the correlation between the input and the output of the networks and the smallest set of disk parameters needed to predict chemical abundances without running the network for each disk model.Our chemical predictions are based on the results of the ANDES astrochemical model of a 2D azimuthally symmetric disk (Semenov & Wiebe 2011), applying the ALCHEMIC chemical code (Semenov et al. 2010). This analysis was performed using the Python machine learning library Scikit-learn (Pedregosa et al. 2011). We employed the KNeighborsRegressor estimator, which identifies the k nearest data points in the input feature space and performs the interpolation of the output feature values among them. The full derivation of the correlation between the disk physical parameters and the ANDES CO abundances, as well as the input parameters of the regressor, are described in detail in Smirnov-Pinchukov et al. (2022). Isotopologue selective processes are not included in the chemical network; therefore, we inferred the iso-topologue abundances using the standard ISM abundance ratios, CO/¹³CO=69 and C¹⁸O=557 (Wilson 1999).

The full gas-grain chemical network includes about 650 species and 7 000 reactions, including gas-phase and surface two-body reactions, adsorption and desorption, photoreactions and ionization/dissociation by X-ray, cosmic rays, short-lived radioactive nuclides, and reactive desorption. Following Eistrup et al. (2016), we adopted an icy molecular initial composition based on the abundance of ice in prestellar cores (Öberg et al. 2011). We ran the time-dependent chemical evolution till the age of 1 Myr. The disk physical structure was set through a stellar mass, M_⋆, a disk mass, M_disk, and a disk tapering radius, r_c. These parameters define the distribution of densities, temperatures, and the high-energy radiation field. The stellar mass also governs the stellar temperature and luminosity, which were calculated for the age of 1 Myr using the evolutionary model by Yorke & Bodenheimer (2008). The ionizing radiation field was computed using the Bruderer et al. (2009) X-ray prescription and the Padovani et al. (2018) cosmic ray prescription. The midplane temperature was computed by integrating stellar and accretion luminosities using a parametric approach, while the atmospheric temperature was computed by solving the vertical UV radiation transfer equation. These two temperatures were then connected using Eq. (5), as in Williams & Best (2014). Iterative calculations of density and temperature distributions were then conducted to achieve a self-consistent solution, as has been outlined by Molyarova et al. (2017).

We created a population of synthetic disks with different M_⋆, M_disk, r_c, and stellar X-ray luminosity, L_x, to cover a wide range of physical conditions typical for protoplanetary disks. Relying on the observational constraints on disk structure, these parameters were varied: M_⋆ = [0.3,2.5] M_⊙, M_disk = [10⁻³, 10⁻¹] × M_⋆, r_c = [20,100] au, and L_X = [10²⁹, 10³¹] erg s⁻¹. These parameters were used to compute 540 synthetic disk models, where each model is composed of 4000 physical cells. This amounts to 2 160000 data points, which took about one year of single CPU run time.

The performance of the machine learning chemical predictions is shown in Fig. 3 for the gas phase CO abundance, adapted from Smirnov-Pinchukov et al. (2022). The fitting procedure successfully reproduces the chemical model predictions, demonstrating minimal systematic errors in panel b and scatter in panel c. A noticeable increase in scatter is observed within a specific region of the parameter space, corresponding to the radiation-sensitive transition zone between the atmosphere and the other disk regions in the low-density area. This region constitutes only a marginal portion of the total CO inventory, and our conclusions remain unaffected by the scatter in CO abundance predictions within this particular disk region.

The main result of the machine learning analysis of this precomputed disk physical-chemical grid is that the disk CO abundances are mainly determined by the local density and temperature of the gas. Further parameters, such as the ionizing and UV radiation, increase the quality of the predictions but are a minor correction to the abundances predicted by the gas density and temperature alone. While the impinging FUV radiation does play a role in the chemical evolution, we found that the effect of the local UV field is correlated with the effect of the gas density and temperature. For instance, UV radiation is stronger in the atmosphere, where the gas density is lower and its temperature higher; hence, the CO abundance in this UV-irradiated region can be correctly predicted using the local low gas density and high temperature. In the DiskCheF framework, we used the same prescription for the X-ray and cosmic ray radiation part of the physics.PhysicsBase object and the CO abundances were predicted from the gas density, gas temperature, and ionizing radiation by the chemistry.SciKitChemistry object.

Fig. 3 Performance of ML-accelerated chemistry predictions for CO, adapted from Smirnov-Pinchukov et al. (2022). (a) Mean log10 predicted relative abundance as a function of local temperature, gas density, and ionization rate. Darker areas correspond to larger relative (to H atoms) abundances. (b) Median of the difference between the predicted values and test set data (bias, dex), in dex, as a function of temperature and density. Gray areas correspond to an unbiased fit. (c) Standard deviation between the predicted values and test set data, in dex (std, dex). (d) Relative density (histogram) of CO within the data points with contours, also shown in the other panels. Various regions of the protoplanetary disk are described in panel a.

3.3 Fitting the interferometric line data

Once we had a physical and chemical disk model, we used RADMC3D to calculate the line radiative transfer of the observed CO isotopologue lines and produce channel maps of these lines with the same spatial and spectral resolution as in the PRODIGE data. However, an interferometer like NOEMA does not produce an image of the source emission, but measures its visibilities, the complex values of its Fourier transform. There are two ways to compare the model channel maps to the data, either by working in the image plane by converting the observed visibilities to an image or by computing the synthetic visibilities of the model by knowing the antenna configuration of the interferometer. While more intuitive, working on the image plane requires several assumptions about the source, and there is no unique way to derive an image from visibility data. This effect is more severe with less sampling of the source visibilities. While NOEMA has a great spectral resolution ideal for deep molecular line surveys, it is more limited in the visibility sampling than ALMA with its >50 antennas.

The second approach is to directly operate in the uu-visibility plane, which represents the real instrument data, by taking the spectral resolution, sensitivity, and antenna configuration into account. This provides a more robust comparison with the data, but if used in a fitting algorithm requires the computation of synthetic visibilities at each likelihood evaluation. For this purpose, we used the GALARIO Python library (Tazzari et al. 2018), a computationally efficient tool for the generation of synthetic visibilities. In DiskCheF, the object uv.UVFits can be used to read, analyze, and visualize visibility data, including computing the χ² of visibility data to a channel map using GALARIO routines (uv.UVFits.chi2_with()).

With this likelihood evaluation, we could now fit the model to the observations. We adopted the nested sampling Monte Carlo algorithm MLFriends (Buchner 2016, 2019), using the UltraNest² package (Buchner 2021). This package provides computationally efficient and optimized MPI-cluster tools to find and analyze the posterior probability distribution of the model parameters. This allowed us to find more reliable results than traditional Markov chain Monte Carlo methods with a reasonable amount of computational resources. To fit the PRODIGE data, we deployed 2560 logical CPUs, corresponding to about 5 × 10⁵ likelihood evaluations within 24 h of computational time. The disk model that we employed in this first application of DiskCheF has only five free parameters: the total disk mass, M_disk, and the tapering radius, r_c, from Eq. (1), and the midplane and atmosphere temperatures, T_mid/atm,100, from Eq. (4). The disk model is not that strongly dependent on the other parameters mentioned in Sect. 3, and we kept them fixed to the values summarized in Semenov et al. (2024).

4 Results

The mass distribution is the most important result of our analysis for planet formation studies. The total disk gas mass can only be traced by optically thin emission. While ¹²CO (2−1) emission is optically thick, and hence can only provide information about the gas temperature at its emission surface, the more optically thin ¹³CO and C¹⁸O (2−1) lines can be used to constrain the disk gas mass (Dartois et al. 2003). To analyze the result of our fitting, we first present the best-fit model for the DM Tau disk. The data quality of this disk is the highest in our sample since it is the largest and brightest one, and we can use this disk to assess the feasibility and goodness of our best-fit model. In Fig. 4, we show the best-fit physical and chemical model for DM Tau. This figure shows how the ML-based chemical model correctly captures the CO freeze-out from the gas phase at low temperatures in the outer midplane, as well as the CO photodissociation in the disk atmosphere. In Fig. 5, the surface density and midplane temperature profiles of our best-fit models are shown.

HD- or CO-based mass estimates can be poorly constrained due to the uncertainties in the underlying gas thermal structure. However, our approach of fitting at the same time optically thin and optically thick CO isotopologue lines allows us to constrain both the disk mass and temperature distribution, leading to more reliable mass estimates. Law et al. (2021) found that the temperature profiles of the disks in the MAPS sample follow the power law in Eq. (4) outside the 100−150 au radius, but have a flatter profile in the inner disk. Thus, two power law profiles would provide a better fit for the high-resolution ALMA MAPS data.

However, the spatial resolution in our sample is about the size of the inner disk, and we cannot properly constrain the temperature profile in this region. Therefore, we assume a single power law distribution for the whole disk structure. We demonstrate the effect of this assumption on the DM Tau disk. In Fig. 6, we show the comparison of the best-fit model and the data in the CO (21) isotopologue channel maps, while in Fig. 7 the model-data residuals are shown. The residuals are larger in the inner disk, where the model underestimates the gas brightness temperature by about 2 K, and where the CO emission becomes optically thick, and hence more sensitive to the gas temperature than the gas density. However, the mass distribution is not affected by this assumption, since it is more strongly constrained by the outer disk emission, and thus we still get a reliable gas mass estimate.

The CO (2−1) line intensity of the output images is strongly affected by the disk temperature profile. Moreover, changes in the disk mass do affect the shape of the CO (2−1) emission lines as well, since the higher disk mass changes the Keple-rian velocity profile of the CO gas. Indeed, in more massive disks there is more mass at higher heights over the midplane, and the resulting projection of the Keplerian velocity profile is noticeably different. The ¹³CO (2−1) data, on the other hand, are more sensitive to the mass distribution than the gas temperature, as we would expect from this (partly) optically thin emission. Therefore, ¹³CO (2−1) data provide a more reliable disk gas mass estimate when fitted together with the ¹²CO (2−1) emission probing mainly the disk thermal structure. The C¹⁸O (2−1) images follow the same trend as the ¹³CO images but, being much fainter, this line is not always detected with a good-enough signal-to-noise ratio to provide additional good constraints on the disk gas mass.

Another parameter affecting the results of the DiskCheF fitting is the depletion of CO. Our predictions are based on the results of the gas-grain ANDES2 chemical model, which does not include isotopologue-selective CO photodissociation processes (e.g., Visser et al. 2009). This may lead to an underestimation of the disk masses probed via the minor CO isotopologue lines, and our estimates could be a lower limit (Miotello et al. 2014, 2016). This is a feature of the machine learning method adopted to speed up chemical modeling, as the ML-chemistry results inherit the same limitations that the chemical model used to produce the training data. The machine learning-accelerated chemistry based on the ANDES2 model made it possible for us to add a complex gas-grain chemical kinetics model into our fitting algorithm, instead of adopting a parameterized chemical model (Williams & Best 2014) or a reduced CO network (Deng et al. 2023). One of the main results of this paper is the application of the fast machine learning method to accelerate calculations of molecular abundances, which can be easily implemented within the flexible DiskCheF framework, using other, disk chemical kinetics models such as DALI or ProDiMo.

We expect the CO (2−1) emission to be in the local thermo-dynamical equilibrium (LTE) with the gas at the disk densities and temperatures from which it is emitted. The best-fit DM Tau model has CO (2−1) brightness temperatures matching the thermal profile of the underlying best-fit model, with a peak brightness temperature over the beam-emitting area of 17 K. In the other disks, the peak brightness temperature over the beam area ranges from 11 K in IQ Tau to 20 K in CI Tau. On the other hand, the ¹³CO and C¹⁸O (2−1) emissions are optically thin, with a ratio of ¹³CO-to-C¹⁸O emission equal to the isotopologue ratio (about 7.3). However, in the inner disk of the more massive disks, the emission ratio falls below the isotopologue ratio (about five in the most massive disk, DM Tau). This is an indication that the ¹³CO emission becomes partially optically thick in the inner disk as well (smaller than our beam size, about 100 au). The C¹⁸O emission remains optically thin at our spatial resolution, and our mass estimate of 0.04 M_⊙ is not affected by this opacity effect.

Fig. 4 Best-fit physical and chemical model for DM Tau: (top left) gas volume density distribution, (top right) CO to H2 number density ratio, (bottom left) gas temperature distribution, and (bottom right) CO isotopologues volume density distribution.

Fig. 5 Best-fit column density and midplane temperature profiles for the disks in the PRODIGE sample.

4.1 Foreground cloud absorption

The optically thick CO emission observed in CI, CY, DG, and DL Tau is affected by foreground absorption within the velocity channel range of approximately 4 to 6 km s⁻¹. Given that these disks exhibit a local standard of rest (LSR) velocity of around 6 km s⁻¹, the foreground cloud impacts the gas signature in the outer disk as it obscures the red-shifted emission from the Kep-lerian disk. This obstruction results in the asymmetrical features evident in the moment zero maps presented in Fig. 1.

To find the best fit in this case, we excluded these channels from our likelihood estimate. In doing so, we emphasize the inner disk emission in determining the best-fit model. This choice is made with consideration of the discussion in Sect. 4.2, where we acknowledge that the inner and outer disk structures can significantly differ, particularly in cases involving a flared outer disk structure. This discrepancy introduces uncertainty into the mass estimation for these partially obscured disks when using optically thin isotopologues. This is because the temperature in the outer disk is deduced from the inner disk structure, and the mass estimate from optically thin CO isotopologues is temperature-dependent.

4.2 Disk vertical structure

The disk physical model used to fit the data, described in Sect. 3.1, reconstructs the gas vertical distribution assuming hydrostatic equilibrium. However, recent studies have shown how the gas emission layer may differ from the one determined by the hydrostatic equilibrium. Law et al. (2021) found that the flaring indices of the CO emission surfaces for the disk in the MAPS sample can be larger than the one coming from hydrostatic equilibrium by a factor of a few. This effect was previously studied in the case of IM Lup, a particularly large and flared disk (Pinte et al. 2018; Franceschi et al. 2023). Similar results have been found for the largest disk in our sample, DM Tau (Dartois et al. 2003). The shape of the emitting region of CO isotopo-logues is determined by the stellar irradiation, rather than the disk scale height, tracing the surface of optical depth $\tau ~\frac{2}{3}-1$ with respect to the stellar irradiation, and this can affect our temperature estimate. Our best-fit model in hydrostatic equilibrium effectively reproduces the observational data for all the disks in our survey. However, a more sophisticated physical model may be necessary to study the disk structure in more detail, such as trying to model other molecules or by including dust continuum emission data into the fitting.

4.3 Temperature profiles

In our initial modeling attempts, we introduced two additional fitting parameters: the exponent governing the midplane and atmospheric temperature distributions, denoted as q_mid/atm in Eq. (4). However, our fitting algorithm revealed a complete degeneracy between T_mid and q_atm, as well as between T_atm and q_mid. This can be attributed to the fact that the observed flux arises from the integrated emission along the line of sight, rendering no formal distinction between the emission originating from a cold midplane and a warm atmosphere, or a warmer mid-plane and a colder atmosphere in the model and, hence, in our channel map data. This degeneracy could potentially be resolved through either higher-angular-resolution data, enabling precise measurements of the height of the emitting surface for each velocity channel, as has been suggested by previous studies (e.g., Dullemond et al. 2020; Law et al. 2021), or by adopting a more realistic underlying disk physical model.

We have adopted a commonly employed assumption in the literature and set the temperature exponent, q_mid/atm, to a fixed value of 0.55. Protoplanetary disks, such as those found in the PRODIGE sample and characterized by the absence of severe substructures or transient accretion outburst phenomena, generally exhibit temperature slopes in the range of 0.4−0.6 (Dutrey et al. 2014; Law et al. 2021; Zhang et al. 2021; Miotello et al. 2023). Given the relatively limited variation in temperature slope and our moderate angular resolution, minor adjustments to the temperature profile are unlikely to significantly impact our best-fit results.

Fig. 6 Comparison between the best-fit model and the observations of DM Tau for the CO, 13CO, and C18O isotopologue emission.

Fig. 7 Residuals of the DM Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.

4.4 Disk masses

The emission of CO isotopologues combined with the results of the ALCHEMIC chemical code allows us to derive the total disk masses. The estimates of the disk gas masses are summarized in Table 2, along with previous dust-based estimates. Generally, the masses obtained from our methodology align with the dust masses derived from continuum emissions, assuming a dust-togas ratio of 100, as was measured in Guilloteau et al. (2011, 2016); Güdel et al. (2018). Notably, exceptions are observed in the cases of DM Tau and DN Tau, where our mass estimates are about twice the values derived from dust emissions.

Several disks in our sample have mass estimates derived from alternative methodologies. For instance, McClure et al. (2016) computed a mass range of [1.0−4.7] × 10⁻² M_⊙ for the DM Tau disk based on HD emission, in accord with our estimate of 4.0 × 10⁻² M_⊙. Another independent measurement for the CY Tau disk mass was provided by Powell et al. (2019), which relates the size of the large grain-emitting region to the total gas mass of the disk. Their estimation of a substantial disk mass, 0.1 M_⊙, approximately 20% of the host star mass, contrasts with our estimate, which is approximately five times lower. This discrepancy can be explained by Franceschi et al. (2022): the method used by Powell et al. (2019) tends to overestimate the disk mass by up to one order of magnitude, justifying the divergence from our results.

The disks in our sample (except for DN Tau) have been previously characterized for their CO isotopologues emissions using a parametric model for the CO abundance by Williams & Best (2014). These mass estimates are also reported in Table 2. These values are approximately one order of magnitude lower than our measurements, indicating a gas-to-dust ratio close to ten. One plausible explanation for this disparity is the assumption made by Williams & Best (2014) that the CO abundance in the warm molecular layer is equivalent to that in the ISM (x_CO = 10⁻⁴). However, processes such as the conversion of CO to CO₂, ice chemistry driven by ultraviolet (UV) radiation, dust evolution and vertical transport, and isotope-selective photodissociation can lower the gas-phase CO abundances in disks considerably (e.g., Trapman et al. 2021; Zhang et al. 2021; Furuya et al. 2022a; Terwisscha van Scheltinga et al. 2022; Van Clepper et al. 2022), leading to a higher total mass estimate from the same observed CO isotopologue emissions.

In the case of the DM Tau disk, Williams & Best (2014) infer a higher gas-to-dust ratio of 60 and a disk gas mass of 0.9 × 10⁻² M_⊙, aligning their measurement more closely (approximately 25%) with our estimate and consistent with the lower edge of the mass range inferred from HD emission (McClure et al. 2016). This observation suggests that CO destruction processes are less pronounced in the bulk of the DM Tau structure. Indeed, our model indicates that this disk is warm, reducing the amount of CO depleted from the gas phase by freeze-out, and massive, increasing the amount of CO molecules surviving photodissociation through self-shielding. These effects converge to bring the CO abundance closer to the ISM value, thereby mitigating the disparity between our results and the mass estimate by Williams & Best (2014).

Indeed, Trapman et al. (2022) used thermochemical models to contain the CO abundance in the gas phase in the DM Tau, finding that x_CO = [0.23−1.3] × 10⁻⁴, supporting our argument that CO abundance in the gas phase is not much lower than its ISM value. Using this x_CO estimate, Trapman et al. (2022) find that the total mass of the DM Tau falls between [3.19.6] × 10⁻² M_⊙, in agreement with our estimate of 4.0 × 10⁻² M_⊙. Our results for DM Tau then align with both HD-based mass measurements and CO depletion studies based on thermochem-ical models, proving the robustness of our methodology for this particular source. This suggests that our findings are likely reliable for other sources as well, although additional independent mass measurements for the remaining sources in our sample would further demonstrate the robustness of our method.

5 Discussion

We employed our DiskCheF fitting model to analyze the disks within our sample, systematically evaluating the fidelity of our best-fit models in replicating the observed data. A comprehensive model-data comparison for all sources is presented in the appendix for brevity, while an overview of the derived best-fit disk physical parameters is presented in Table 2. It is worth noting that the mass estimates of our best-fit models are contingent upon both the luminosity of the CO emission and its kinematic properties, which collectively define the morphology of the emission across each velocity channel. The fitting of individual velocity channels consequently provides robust constraints on the disk masses.

A limitation shared by all our best-fit models pertains to the brightness of the CO isotopologues within the initial few tens of au in the inner disk. This luminosity is consistently brighter than what is observed in the data. This discrepancy can be attributed to the assumed power-law prescription for the temperature distribution. As is discussed in Sect. 4.3, this particular temperature distribution tends to overestimate temperatures within the inner disk, thereby leading to elevated residuals in this region. Nevertheless, most of the disk mass is found in the outer disk, and this discrepancy does not significantly affect our mass estimation.

Our estimates of disk mass align well with previously reported values found in the literature, summarized in Table 1, which are derived from dust emission analyses. An outlier is the DM Tau disk, which is the most massive source in the sample and much more radially extended even when compared to the DL Tau disk, which shares a similar mass. Disks with these characteristics are subjected to more rapid dust evolution processes. Consequently, they tend to exhibit higher dust-to-gas ratios than the commonly assumed value of 100, which is employed to estimate disk masses in Table 1 (e.g., Franceschi et al. 2023). This heightened dust-to-gas ratio justifies why our model mass estimate is double the value of the dust-based mass estimate from Table 1.

Our best-fit models effectively reproduce the observed emission within each channel map, reproducing both the CO line brightness and spatial distribution. The CI Tau, CY Tau, and DL Tau disks are, however, a notable exception. In these particular sources, our models reproduce well the observed ¹²CO (2−1) emission, not the observed optically thin minor CO isotopo-logue emission lines. This discrepancy does not originate from a wrong disk mass estimate, as modifying the mass would change the kinematic gas distribution, which, in contrast, is well matched by our models. This deviation is also not a consequence of temperature-related effects. Raising the atmospheric temperature, for instance, would result in an overly intense CO emission, while increasing the midplane temperature would lead to a similar outcome for the C¹⁸O emission. The only plausible explanation for this behavior, where optimizing the fit for one isotopologue negatively impacts the fit for another, is the assumption of an incorrect isotopologue ratio.

Our fitting routine, based on assumed CO isotopologue ratios for these sources, boosts the ¹³CO line by increasing its temperature, and since this isotopologue is more optically thin than CO, it is more sensitive to the midplane temperature, which is increased to values similar to the disk atmospheric temperature. The disk atmosphere is usually warmer by 10−20 K in the outer, ≳100−200 au T Tau disks, and the anomaly in the temperature profile can be used as a signature of a wrongly assumed¹² CO/¹³CO isotopologue ratio. Indeed, in the instance of the CI Tau disk, we refined our model fit to the observational data by increasing the¹³ CO isotopologue ratio by a factor of five compared to the ¹³C/¹²C ratio. Simultaneously, we adjusted the C¹⁸O abundance, reducing it by a factor of three in relation to the ¹⁸O/¹⁶O ratio. In the case of the CY Tau and DL Tau disks, we find that optimizing the ¹³CO isotopologue ratio by doubling it and adjusting the C¹⁸O isotopologue ratio by halving it yields a superior alignment with the data.

It is crucial to highlight that our approach does not rely solely on the isotopologue brightness profile to constrain the disk mass. Instead, it also considers their spatial distribution across each velocity channel, determined by the gas kinematic structure. The gas kinematics depend upon the disk mass, influencing both the brightness of the optically thin isotopologues and their spatial distribution.

According to our model, a correct prediction of the gas kinematic structure with an inaccurate optically thin isotopologue emission with respect to the observations is more plausible than a model with the wrong kinematic structure and accurate optically thin isotopologue emission. Therefore, the isotopologue ratios primarily impact the gas temperature estimate rather than the already constrained disk mass determined by gas kinematics. Consequently, adjustments to the isotopologue abundances do not significantly impact our overall mass estimates.

We stress that these estimates of the CO isotopologue ratios do not stem from our fitting algorithm, as the depletion factor is a degenerate parameter with the gas density and temperature. To obtain more feasible fits to the CO disk data with DiskCheF model, it is imperative to train the ML algorithm using a chemical network that accounts for ¹³C and ¹⁸O-fractionation and isotope-selective photodissociation, following the methodology described in this work.

It is interesting to note that previous modeling studies focused on fitting the observed CO emission in disks have also inferred lower ¹²CO/¹³CO ratios of ~ 20–40 than the local ISM ¹²C/¹³C isotopic ratio of about 69 (for a recent comprehensive review, see Nomura et al. 2023). Using the Plateau de Bure multi-line, CO isotopologue observations, Piétu et al. (2007) have derived the ¹²CO/¹³CO ratio of ~20 in the cold outer DM Tau T Tauri disk. Recently, Yoshida et al. (2022) have also derived a similar ¹²CO/¹³CO ratio of ~20 in the outer T Tauri disk around TW Hya, while in the inner, ~20 au region of the TW Hya disk this ¹²CO/¹³CO ratio seems to be higher, ~40 (Zhang et al. 2017). Furthermore, the H¹³CO⁺/HC¹⁸O⁺ ratio of ~8 has been measured in the TW Hya disk by Furuya et al. (2022b), which coincides with the local ISM value within the uncertainties.

The reason why the ¹²CO/¹³CO ratio measured in the T Tauri disks could be lower than the local ISM value is carbon fraction-ation in the cold, ~20−25 K gas, typical for the outer regions of these disks probed with CO rotational lines (e.g., Woods & Willacy 2009; Nomura et al. 2023). The key fractionation reaction is ¹³C⁺ + ¹²CO →¹³CO + ¹²C⁺ + ∆E, where ionized atomic carbon is produced by the cosmic ray ionization and the energy barrier, ∆E, is about 35 K. The forward pathway of this reaction is more efficient than the backward route in the gas with T ≲ 20–30 K, which corresponds to the outer, ≳50−100 au regions of T Tauri disks (the densest, bottom part of the outer, r ≳ 100 au molecular layer). Indeed, in the theoretical modeling of Woods & Willacy (2009), it has been found that the ¹²CO/¹³CO ratio could become as low as about 25 in a T Tauri disk due to this process. More detailed studies of carbon fractionation in a statistically significant sample of disks are needed to further understand the efficiency of this fractionation and the magnitude of its impact on the observed ¹²CO/¹³CO ratios and, hence, disk gas mass estimates using the CO disk emission.

6 Conclusions

The gas mass of a protoplanetary disk is an important but challenging quantity to measure. The CO emission remains the most reliable gas tracer, but CO-based mass estimates rely on the assumed gas temperature distribution and disk chemistry. A combination of optically thin and optically thick emission lines, however, can solve the degeneracy between the temperature and the abundance of CO isotopologues emission, with the help of chemical models. In this work, we propose a machine learning-based chemical model to predict a disk chemical composition from its physical structure without the need to run a chemical network. This makes it possible to reproduce observations of line emission through theoretical models using a fitting algorithm, which would otherwise require unreasonable computational resources. We tested the applicability of this method using the CO isotopologue emission data from the PRODIGE program on the NOEMA instrument. This data consists of CO, ¹³CO, and C¹⁸O isotopologue emission from six Class II pro-toplanetary disks from the Taurus region. Our findings can be summarized as:

• The combination of optically thin and optically thick CO iso-topologue emission, combined with the kinematic gas distribution in each velocity channel, can properly constrain the disk physical structure and is in good agreement with previous dust-based mass estimates. This approach gives tighter constraints on the disk’s thermal structure and provides a more reliable mass estimate;

• There is evidence that the emission surface of optically thick lines does not trace the gas pressure scale height, especially for a large disk with evidence of a flared structure, such as DM Tau. The vertical gas structure could not be in hydrostatic equilibrium in the inner disk, making it challenging to reproduce the inner and outer disk emission with a single temperature profile. However, most of the disk mass is found in the outer disk, and this effect does not significantly affect our mass estimates;

• The cloud absorption obscures the optically thick emission in the channels close to the rest frequency for a few disks in our sample. These channels show the emission coming from the outer regions of the disk and give us information on the temperature structure. For these obscured disks, the temperature structure is constrained only by the inner disk emission, which could have a different profile than the outer disk. In a particularly flared structure, the difference between the inner and outer disk structures adds to the uncertainty of the mass estimate;

• In some sources, the CO isotopologue abundances are determined by isotopologue-selective processes. In these sources, the mass is still well constrained by the gas kinematic structure, but the temperature profile is offset to boost the emission of under-abundant isotopologues, and we get a worse fit to the data. By changing the isotopologues ratios, we get a better fit to the data, allowing us to estimate the isotopologue abundances in these sources.

Acknowledgements

This work is based on observations carried out under project numbers L19ME and S19AW (for DNTau and CI Tau) with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). D.S., Th.H., G.S.-P., R.F, K.S., and S.vT, acknowledge support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant No. 832428-Origins. This work was partly supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. This research made use of NASA’s Astrophysics Data System.

Appendix A Model-data comparison and residual channel map

Fig. A.1 Comparison between the best-fit model and the observation of CI Tau for the CO, 13CO, and C18O isotopologue emission.

Fig. A.2 Residuals of CI Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.

Fig. A.3 Comparison between the best-fit model and the observation of CY Tau for the CO, 13CO, and C18O isotopologue emission.

Fig. A.4 Residuals of CY Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.

Fig. A.5 Comparison between the best-fit model and the observation of DL Tau for the CO, 13CO, and C18O isotopologue emission.

Fig. A.6 Residuals of DL Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.

Fig. A.7 Comparison between the best-fit model and the observation of DM Tau for the CO, 13CO, and C18O isotopologue emission.

Fig. A.8 Residuals of DM Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.

Fig. A.9 Comparison between the best-fit model and the observation of DN Tau for the CO, 13CO, and C18O isotopologue emission.

Fig. A.10 Residuals of DN Tau best-fit model for the CO,13CO, and C18O isotopologue emission.

Fig. A.11 Comparison between the best-fit model and the observation of DM Tau for the CO, 13CO, and C18O isotopologue emission.

Fig. A.12 Residuals of IQ Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.

References

All Tables

All Figures

	Fig. 1 Moment zero maps of the CO, 13CO and C18O (2−1) emission for the disks chosen for the DiskCheF fitting.
In the text

	Fig. 2 Outline of the steps in a DiskCheF model. For each step, we indicate which template can be used to create a new custom model.
In the text

Fig. 3 Performance of ML-accelerated chemistry predictions for CO, adapted from Smirnov-Pinchukov et al. (2022). (a) Mean log10 predicted relative abundance as a function of local temperature, gas density, and ionization rate. Darker areas correspond to larger relative (to H atoms) abundances. (b) Median of the difference between the predicted values and test set data (bias, dex), in dex, as a function of temperature and density. Gray areas correspond to an unbiased fit. (c) Standard deviation between the predicted values and test set data, in dex (std, dex). (d) Relative density (histogram) of CO within the data points with contours, also shown in the other panels. Various regions of the protoplanetary disk are described in panel a.

In the text

	Fig. 4 Best-fit physical and chemical model for DM Tau: (top left) gas volume density distribution, (top right) CO to H2 number density ratio, (bottom left) gas temperature distribution, and (bottom right) CO isotopologues volume density distribution.
In the text

	Fig. 5 Best-fit column density and midplane temperature profiles for the disks in the PRODIGE sample.
In the text

	Fig. 6 Comparison between the best-fit model and the observations of DM Tau for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. 7 Residuals of the DM Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.1 Comparison between the best-fit model and the observation of CI Tau for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.2 Residuals of CI Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.3 Comparison between the best-fit model and the observation of CY Tau for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.4 Residuals of CY Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.5 Comparison between the best-fit model and the observation of DL Tau for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.6 Residuals of DL Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.7 Comparison between the best-fit model and the observation of DM Tau for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.8 Residuals of DM Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.9 Comparison between the best-fit model and the observation of DN Tau for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.10 Residuals of DN Tau best-fit model for the CO,13CO, and C18O isotopologue emission.
In the text

	Fig. A.11 Comparison between the best-fit model and the observation of DM Tau for the CO, 13CO, and C18O isotopologue emission.
In the text

	Fig. A.12 Residuals of IQ Tau best-fit model for the CO, 13CO, and C18O isotopologue emission.
In the text