© The Authors 2026

1 Introduction

Thousands of exoplanets have been discovered in the last 20 years¹, but many fundamental questions remain about how they were born. Our picture of planet formation is incomplete, mainly because we still do not understand how the protoplanetary disks that orbit young stars evolve and eventually form planets (Morbidelli & Raymond 2016). Surveys with the Atacama Large Millimeter/submillimeter Array (ALMA) have greatly improved our knowledge regarding protoplanetary disk evolution, enabling population studies to investigate how gas disks evolve in different star-forming regions (SFRs; e.g., Ansdell et al. 2018; Sanchis et al. 2021; Long et al. 2022; Zhang et al. 2025).

The typical observables for gas-disk evolution studies are line fluxes and radii. In particular, the latter can be used to constrain the dominant mechanisms controlling the evolution of protoplanetary disks (for a review, see Manara et al. 2023). For example, Tabone et al. (2025) recently used population-synthesis models to show that magnetohydrodynamic (MHD) wind-driven evolution best reproduces the typical masses and sizes of protoplanetary disks. The biggest challenge for measuring precise bulk gas-disk properties (such as disk radii and line-emission surfaces) for large samples is that most surveys are observed with moderate angular resolution and modest signal-to-noise ratios (S/N); this is a necessary trade-off to observe a large number of disks. For this reason, ALMA shallow surveys are often used as precursors for dedicated high-resolution observations of particularly bright and extended disks (see, e.g., Öberg et al. 2021; Teague et al. 2025).

One possibility for deriving bulk disk properties from moderate-angular-resolution data is offered by modeling spectral-line visibilities. So far, this technique has functioned via radiative-transfer codes (e.g., Guilloteau & Dutrey 1998; Simon et al. 2000; Dartois et al. 2003; Piétu et al. 2007; Rosenfeld et al. 2012; Flaherty et al. 2015) or parametric emission models (Kurtovic & Pinilla 2024) on a handful of resolved protoplanetary disks, but it has never been applied on a large survey: in this work, we filled this gap.

We implemented two new simple gas-disk models and used the software csalt² (Andrews et al., in prep.) to infer important quantities of protoplanetary disks, such as radius, line-emission layer, and integrated flux. We conducted our analysis on the ¹²CO J = 3–2 data presented by Carpenter et al. (2025) for the Upper Scorpius old SFR (4–14 Myr, Ratzenböck et al. 2023a; Ratzenböck et al. 2023b).

This paper is structured as follows. In Sect. 2, we present the sample over which we conducted our analysis. In Sect. 3, we show the parametric models adopted in this work and discuss the details of our modeling technique. In Sect. 4, we show an overview of our results and provide a comparison with literature values. In Sect. 5, we discuss our results and highlight confirmed and candidate multiple systems identified in our sample. Finally, in Sect. 6 we draw our conclusions.

In this paper, we do not comment on the stellar masses we derived, even though we report their values for completeness. This specific topic will be covered by Zallio et al. (submitted to A&A).

2 Sample selection

The sample was assembled from the work of Carpenter et al. (2025), which presented ALMA continuum and ¹²CO observations (project codes 2011.0.00526.S, 2012.1.00688.S, 2013.1.00395.S, and 2018.1.00564.S) for a total of 284 Class II objects. The observations were taken with the 12-m array using between 43 and 50 antennas, providing an angular resolution from 0.1″ − 0.3″, with a typical on-source integration time of 2.5 minutes. Among this sample, Carpenter et al. (2025) identified 202 sources associated with the Upper Scorpius SFR. The observations were obtained in Band 7, with spectral windows centered at 334.2, 336.1, 346.2, and 348.1 GHz. The window centered at 346.2 GHz included the ¹²CO J = 3–2 transition, with a channel width of 0.488 MHz (0.429 kms⁻¹ for ¹²CO J = 3–2). The data were calibrated with the ALMA calibration pipeline. Among the 202 sources, Carpenter et al. (2025) identified 83 ¹²CO J = 3–2 detections with S/N ≥ 3; these are the sources we analyzed.

In addition to the ALMA data of protoplanetary disks, we used the Very Large Telescope/X-Shooter spectra presented in Manara et al. (2020) (project codes 097.C-0378, 0101.C-0866) and new data that will be published in Empey et al. (in prep.) (project codes 105.2082.003, 113.26NN.001, 113.26NN.003, and 115.27XL.001) to derive the stellar properties of the disks that overlap with our ALMA sample.

3 Analysis

In this section, we present our analysis technique to constrain protoplanetary-disk parameters. Choosing a proper parametric model, it is possible to reproduce the channel maps (or their integrated maps) with high fidelity, and to interpret the best-fit parameters that maximize the specified likelihood (Teague et al. 2018; Pinte et al. 2018; Casassus & Pérez 2019; Izquierdo et al. 2021).

We performed our analysis using the software csalt (Andrews et al., in prep.), which allows robust inference on model parameters using visibility fitting. In this work, we present the first large-scale application of csalt to fit and extract information from moderate-angular-resolution ALMA observations. For each model calculation, we generated a set of model visibility spectra in the kinematic Local Standard of Rest (LSRK) frequencies at each timestamp, convolved it with the spectral response function (SRF) of the ALMA correlator, and then calculated a likelihood that explicitly treats the spectral covariance following Loomis et al. (2018). Finally, for user-specified priors, we sampled the posterior space using emcee (Foreman-Mackey et al. 2013).

3.1 Model description

In our sample, we identified two main groups of disks we could fit using csalt: “resolved” and “marginally resolved”. In the next sections, we present two different parametric prescriptions used to fit these two groups.

3.1.1 Parametric model A – detailed disk description

The first parametric model used in this work was designed to be simple and to mimic a zeroth-order radiative-transfer model. In particular, this parametric model is conceptually similar to the one implemented in DISCMINER (Izquierdo et al. 2021), but it uses a different parametrization for the line width. With this prescription, we fit the 33 “resolved” disks of our sample. We divided the model parameters into the four groups listed below.

1. Geometrical parameters. We fit for the protoplanetary disk inclination, PA, and offsets dRA and dDec. The PA notation was used to represent the angle to the redshifted major axis of the emission at the midplane when moving E of N in the plane of the sky.

2. Vertical height. We assumed that all the disk emission comes from two thin vertical layers, with the back surface beyond the front one from the perspective of the observer. We fit the height of the layer (vertical height) through an exponentially tapered power law. For simplicity, as the sources of our sample are observed at moderate-to-low angular resolution, and therefore we are not able to distinguish between the front and back surface layers, we set the back emission surface z_b equal to the front surface $z_f:z_b=-z_f^3$³. The vertical height equation is $$z(r)=z_1{\left(\frac{r}{1^{″}}\right)}^{{\psi }_z}\exp \left(-{\left(\frac{r}{r_z}\right)}^{{\varphi }_z}\right).$$(1)

3. Velocity. The model takes into account pure Keplerian rotation v_kep, which is projected into the observer frame. Then, the systemic velocity v_sys is added to obtain the velocity on the plane of the sky v_sky. The equations are $$v_{kep}(r,z)=\sqrt{\frac{G{M}_{\ast }r}{{\left(r^2+z^2\right)}^{3/2}}},$$(2) $$v_{\text{sky}}(r,z)=v_{\text{kep}}\cos (\theta )\sin (|inc|)+v_{\text{sys}},$$(3)

   where G is the gravitational constant, M_∗ is the stellar mass, θ is the azimuthal angle in cylindrical coordinates, and inc is the disk inclination.

4. Line parameters. The model uses a power-law prescription for the brightness temperature, similarly to other studies (see Law et al. 2021). We then adopted an exponentially tapered power law to spatially characterize the optical depth. This choice has the advantage of enabling a more accurate treatment of optically thin lines and line wings. Both temperature and optical depth profiles are assumed to be equal for the front and back surfaces. The equation for temperature reads $$T(r)=T_{10}{\left(\frac{r}{10\text{\hspace{0.17em}au}}\right)}^q,$$(4)

   while the equation describing the optical depth is $$\tau (r)={\tau }_{10}{\left(\frac{r}{10au}\right)}^{{\psi }_{\tau }}\exp \left(-{\left(\frac{r}{r_{\tau }}\right)}^{{\varphi }_{\tau }}\right),$$(5)

   where we considered log₁₀(τ₁₀) when sampling the posteriors. The line-profile shape is characterized by a Gaussian, such that $${\tau }_{\nu }(r,z)=\tau (r)\exp \left(-{\left(\frac{v-v_{sky}(r,z)}{\text{Δ}v}\right)}^2\right),$$(6)

   where the line width, ∆v, is purely thermal and described through $$\text{Δ}v(r)=\sqrt{\frac{2k_{B}T(r)}{\mu m_H}}.$$(7)

Here, k_B is the Boltzmann constant, µ is the molecular mass of the emitting molecule, and m_H is the hydrogen mass.

Finally, we computed the surface-brightness profile for the front surface as $$I_{v,f}(T,r,z)=B_v(T)\left(1-\exp \left(-{\tau }_v(r,z)\right)\right),$$(8) while the one for the back surface is $$I_{v,b}(T,r,z)=B_v(T)\left(1-\exp \left(-{\tau }_v(r,z)\right)\right)\exp \left(-{\tau }_v(r,z)\right),$$(9) where B_ν is the Planck function. Finally, the total surface brightness profile is computed as $$I_v(T,r,z)=I_{v,f}(T,r,z)+I_{v,b}(T,r,z).$$(10)

We summarize the equations of prescription A in Table 1, dividing them by the physical disk attribute they describe.

We stress that even though model A provides a full description of a protoplanetary disk relying on the assumptions discussed above, this is still very simplistic in comparison with the real complexity of the disk. However, given the limited quality of the data used in our study, this approximation is adequate.

3.1.2 Parametric model B – minimal disk description

To extract information from the four marginally resolved disks we have in the sample, we built a secondary disk prescription, which is more suited to describing their faint and compact line emission. The model B case assumes the emission is generated in the midplane (i.e., there is effectively no back surface, and the emission surface has z = 0 at all radii) and has a constant optical depth τ₀ out to radius R_out. We made this choice because the poor angular resolution did not allow us to resolve the optically thin outer disk, so we could not put any constraints on the disk’s vertical structure. In Table 2, we summarize the resulting equations for parametrization B.

3.2 Best-fit models, residuals, and gas-disk size extraction

In this section, we discuss the fitting procedure in detail, using the J16213469-2612269 disk as an example. This disk is compact, but still resolved in the image plane, and thus demonstrates the advantages of our fitting approach. The disk parametrization we adopted to fit this disk is prescription A (see Table A.1 for the collection of best-fit parameters).

3.2.1 Fitting procedure and residual plots

We used Markov chain Monte Carlo (MCMC) algorithms to model the posterior distribution of our parametric models, and we sampled it using emcee (Foreman-Mackey et al. 2013), which is implemented inside csalt. Inspecting the channel maps, we chose priors for inc, PA, v_sys, M_∗, r_τ and r_z for each disk. We chose a conservative number of burn-in steps (250), we let 90 walkers evolve for 1000 steps afterward, and we then checked the results. If the relative change in the median value of each parameter was below 5% over the last 500 steps, the fit was stopped; otherwise, the walkers were allowed to evolve for additional steps until this condition was met. The results that we report in Appendix A (Tables A.1, A.2) are the median as the best fit and the uncertainties as its ±68% confidence intervals. In Table 3, we summarize the priors we chose for the fit on J16213469-2612269 performed with model A. The setup for the priors of each fit on the full sample for both parametrizations can be found on GitHub⁴.

After convergence was met, we imaged the best-fit model and residual visibilities in the same way as the data: in Fig. 1 we show the residual plot for J16213469-2612269 (the full collection of residual plots can be found in Appendix B, while their corresponding corner plots can be found on GitHub⁴.). Here, in the first row we show the channel maps of the data. To produce these channel maps, we used the tclean routine in Common Astronomy Software Applications (CASA, McMullin et al. 2007) version 6.6.0.20. We chose a different velocity range for each disk so that it was possible to visualize the full-disk emission in ten channels.

For cleaning, we used an image size of 1500 pix, with a cell of 0.01″. We set four different scales (0, 10, 30, 60), with a Briggs weighting of robust=1, a threshold of 20 mJy per beam. We also specified a Keplerian mask based on the best-fit parameters for each specific disk. In the second row, we show the channel maps of the best-fit results: to clean the best-fit visibilities we used the same cleaning parameters as above. In the third row, we show the channel maps of residuals, which were obtained by cleaning the residual visibilities (data – model). These channel maps are particularly useful, because if any non-Keplerian motion is found in the data, it manifests as either a red or blue blob in the residual maps. In the case of J16213469-2612269, we only see residual maps of noise.

Fig. 1 Residual plot for J16213469-2612269. In the first row, we show the channel maps of the data; in the second row, we show the channel maps of the best-fit model; in the third row, we show the channel maps of the residuals (data – model). The residual plots for all the disks analyzed in this work can be found in Appendix B and on GitHub4.

3.2.2 Gas-disk radii extraction and comparison with image-plane analysis

Once we derived the best-fit model, we computed the gas-disk radii. This is usually done in the literature by first computing velocity-integrated (moment 0) maps, extracting the midplane disk coordinates from the sky projection, and azimuthally averaging to obtain the radial profile, which is then integrated to extract the radius that encompasses 68% and 90% of the disk flux (R_68% and R_90%, respectively; see Ansdell et al. 2018; Sanchis et al. 2021; Galloway-Sprietsma et al. 2025; Trapman et al. 2025). However, this methodology has an intrinsic problem, which is that the radius of the disk depends on the inclination and on the projected velocity field of the disk even after the disk was deprojected from the plane of the sky, as we discuss in Appendix C. To avoid this intrinsic effect, to determine the gas-disk radii, we first computed the peak surface-brightness profile of the disk front surface at the ¹²CO J = 3–2 rest frequency (ν = 345.7959 GHz) using Eq. (8), rather than calculating the moment 0 map. Secondly, we integrated the peak surface-brightness profile in cylindrical coordinates and extracted the radius that encompasses 68% and 90% of the line emission.

To determine the uncertainties of the gas-disk radii, we selected all random draws from the MCMC chain parameters needed in Eq. (8) and computed R_68% and R_90% as above. We then took the median value of the measured radii, R_90%, and the 16^th and 84^th percentiles to derive their upper and lower errors, respectively.

We are aware that the radius definition adopted in this work intrinsically differs from the one usually adopted in literature, as we integrated the peak line flux rather than the spectrally integrated line flux. For this reason, we stress that our results are not directly comparable with what has been done in the literature so far. To provide a comparison with what is usually done in the literature, in Fig. 2 we show the complementary cumulative distribution function (CCDF⁵) of R_90% of our sample and compare it with that of Zagaria et al., in prep., which analyzed the same sources of our sample, but using an image-plane analysis as described above. After running the Anderson-Darling two-sample test (AD; as implemented in scipy⁶; see Virtanen et al. 2020), we obtained an AD p-value p = 0.8 in the range of [0.6, 0.9], which means that the two distributions cannot be distinguished within uncertainties. This check shows that the methodology for extracting gas-disk radii presented in this paper provides results that are comparable to what is usually done in the literature when working on the same sample.

Fig. 2 CCDF of 12CO R90% computed for this work (in blue) and of the radii evaluated using image-plane analysis from Zagaria et al., in prep. (in orange).

3.2.3 Fitting small disks with prescription A

Using csalt, we can also retrieve bulk properties for small protoplanetary disks. In Fig. 3, we show the residual plot of J16062861-2121297, a compact disk (gas-disk size R_90% = 56 ± 1 au, with a beam of ∼25 au). Looking at Fig. 3, we see how csalt reproduces the disk morphology well using the parametric prescription A. We are also able to infer the vertical structure of the thick ¹²CO emitting layer from this fit. The channel maps of the residuals are almost completely clean, showing that an optimal fit was found for this source.

3.2.4 Dealing with cloud-absorbed channels

Some disks in our sample showed signs of foreground cloud absorption: the disk emission is either obscured in some channels or completely embedded in more diffuse cloud emission. To deal with such problematic features, we first checked the channel maps and the spectra of each of the disks in our sample and identified the frequency ranges that showed clear cloud absorption. Then, we used an option implemented in csalt that allows one to exclude specific frequency ranges from the MCMC likelihood evaluation; this allowed us to consider only clean disk emission in our fits. In Fig. 4, we show the residual plots we obtained for the cloud-contaminated disk J16123916-1859284. This plot shows clear 3 and 5σ residuals in the channel maps at v = 4.6, 5.5 km s⁻¹, which are the channels affected by cloud absorption. Moreover, in the channel map at v = 3.7 km s⁻¹ we show how csalt is able to separate the diffuse cloud emission, which results in red residuals, and Keplerian disk emission.

4 Results

In this section, we present our results and provide median values of the relevant physical quantities of the disks in Upper Scorpius. We were able to fit 37 sources out of a sample of 83, 33 with prescription A and four with prescription B. We present the best-fit results in Tables A.1 and A.2, respectively.

We attempted to perform our analysis on all 83 ¹²CO J = 3–2 detections, but we were unable to achieve any acceptable convergence for our parametric models on the remaining 46 disks for a variety of reasons. The first is the poor spectral sampling that we have for these sources: if ¹²CO is only weakly detected in a few channels (a bandwidth of ≲2 km s⁻¹), csalt cannot fit a Keplerian rotation spectrum due to lack of information from the data. The second is very low S/N; by looking at the integrated flux reported in Carpenter et al. (2025), most of these sources have integrated S/Ns of ≲10. We therefore empirically report that visibility fitting becomes difficult if the S/N of the integrated flux is ≲10. The third reason is poor angular resolution (the disk size is comparable with the size of the beam): if a disk is both spectrally and spatially unresolved, disk fitting fails even if visibilities are used. We report that a combination of these three effects is causing the fits to fail for 44 of these remaining objects. Finally, two peculiar objects that we were unable to fit are J16113134-1838259, which is a potential fly-by, and J16193570-1950426, which may hide at least one companion, as further detailed in Sect. 5.

4.1 Gas-disk radii

We find a median gas-disk size of R_90% ∼ 82 au, with the largest disk being J16042165-2130284, a well-known face-on disk (see, e.g., Stadler et al. 2023 and references therein) for which we derive R_90% = 247 ± 1 au. The most compact disk is J16095933-1800090, for which we derive R_90% = 22 ± 3 au. However, as said before, our results are sensitivity-limited. Since there is a tight correlation between the gas-disk radius and the line flux (Long et al. 2022; Zagaria et al. 2023, and discussed later on in this section), we conclude that the 44 disks that could not be fit (we excluded J16113134-1838259 and J16193570-1950426 as discussed before) are smaller than or of comparable size to J16095933-1800090. As a consequence, the real median gas-disk size for the ¹²CO detections should be, at the zeroth order, R_90% ≲ 22 au. A detailed discussion on this topic will be presented in Zagaria et al., in prep.

In Fig. 5, we report the CCDF of R_90% of our sample and compare it with that of Trapman et al. (2025) for the AGE-PRO Upper Scorpius sample. Our median gas-disk size is somewhat lower than the value R_CO ∼ 109.9 au reported in Agurto-Gangas et al. (2025) and Zhang et al. (2025). To determine if our results are consistent, we ran the AD two-sample test, and we obtained a p-value p = 0.7 in the range of [0.5, 0.8]. This means that the two distributions of R_90% cannot be distinguished within the uncertainties. We again emphasize the different techniques used in Trapman et al. (2025) and in our work to determine the gas-disk radii, and we note that the two studies involve different sample sizes and ALMA observations (with AGE-PRO observations conducted using ALMA ¹²CO J = 2–1 data with an integration time ∼1 hour on source).

4.2 Integrated fluxes

In Appendix D, we explain how we extracted the fluxes both from data and best-fit models, while in Table 4 we report them together with their associated uncertainties. We also include a column that shows the velocity ranges per source for which we identified foreground cloud absorption, and a column that reports whether after the fit any 3 or 5σ residuals are visible in the residual plots.

¹²CO flux – size correlation

In Fig. 6, we show the correlation between the ¹²CO sizes and the ¹²CO disk luminosity L_CO = F_CO × (d [pc] / 140 pc)² Jy km s⁻¹, where we highlight, with different markers and colors, the disks we fit with different parametrizations and those that showed clear absorption. We fit the correlation using the linmix⁷ package and find the best-fit correlation: log₁₀(R_90%) = (1.77 ± 0.03) + (0.49 ± 0.05) log₁₀(L_CO). The slope we derived is comparable to that found in Long et al. (2022), which reported the relationship log₁₀(R_CO) = (0.52 ± 0.05) log₁₀ (L_CO) + (2.07 ± 0.03) using¹²CO J = 2–1 data. The slope we measured is in agreement with the fact that the ¹²CO layer is optically thick⁸ and therefore can be used as a proxy for the gas temperature of the emitting layer (e.g., Sanchis et al. 2021; Long et al. 2022), as further discussed in Sect. 5. However, the intercept found in this work is different from the one reported by Long et al. (2022), with a difference of ∼0.3 in the intercept and a significance of ∼7σ; this is likely the result of using different techniques to extract gas-disk radii, combined with having different data qualities in different samples.

Fig. 3 Residual plot for J16062861-2121297, shown in the same way as in Fig. 1. The residual plots for all the disks analyzed in this work can be found in Appendix B and on GitHub4.

Fig. 4 Residual plot for J16123916-1859284, shown in the same way as in Fig. 1. We highlight the residuals in the channel maps at v = 4.6, 5.5 km s−1, which are those affected by cloud absorption. Moreover, in the channel map at v = 3.7 km s−1, we show how csalt is able to separate the diffuse cloud emission (the localized red emission) and the Keplerian disk emission. The residual plots for all the disks analyzed in this work can be found in Appendix B and on GitHub4.

Fig. 5 CCDF of 12CO R90% for our sample (in blue) and of the sources in Upper Scorpius reported from Trapman et al. (2025) (in orange). The gas-disk sizes reported in this figure are reported in Tables A.1 and A.2.

Fig. 6 Correlation between 12CO gas-disk size and luminosity. The values for R90% and the disk fluxes, together with the line absorption ranges, can be found in Tables A.1, A.2, and 4, respectively.

4.3 Vertical heights

Our first disk model includes information on the vertical location of the line-emitting layer. We derived a collection of ¹²CO J = 3–2 line-emission surfaces of 33 compact disks. We created the profiles according to Eq. (1) using random draws from the MCMC chains of z₁, ψ_z, r_z, and ϕ_z, evaluating the associated z(r) profiles, and taking their median value and the 16^th and 84^th percentiles for the associated uncertainties. In Fig. E.1 of Appendix E, we show the vertical-height profiles as a function of radius.

Through vertical-height profiles, we computed the values of the aspect ratio ⟨z/r⟩ in the same way as Law et al. (2021) and Paneque-Carreño et al. (2023) to homogenize our analysis. This was first done by defining r_cut as the value of the radius for which the vertical-height profile stops increasing. At this point, we defined r_taper = 0.8 × r_cut, and we evaluated ⟨z/r⟩, with 0 ≤ r ≤ r_taper.

We find a median value of ⟨z/r⟩ ∼ 0.16, which is lower than the one reported in Galloway-Sprietsma et al. (2025) (⟨z/r⟩ ∼ 0.28). This is likely because we analyzed a very different sample. The exoALMA sample (Teague et al. 2025) is composed of bright and extended protoplanetary disks; the same applies for the disks considered in the works of Law et al. (2022, 2023, 2024). Instead, we considered the disk population of Upper Scorpius, which is composed of both compact and extended disks. In Fig. 7, we show the Kaplan-Meier estimator (as implemented in lifelines⁹; see Davidson-Pilon et al. 2019) of the radially averaged ⟨z/r⟩ values for our sample and compare it with those of Galloway-Sprietsma et al. (2025), and the combined works of Law et al. (2022, 2023, 2024). We ran the log-rank test (again, as implemented in lifelines), which accounts for upper limits via the Kaplan-Meier estimator, finding that the distribution of ⟨z/r⟩ derived for this work differs significantly from the one of Galloway-Sprietsma et al. (2025) (p = 0.05 in a range [0.02, 0.12]) and the one of Law et al. (2022, 2023, 2024) (p = 0.03 in a range [0.01, 0.08]). A caveat of our analysis is that the ⟨z/r⟩ values derived for our disk sample are based on observations with significantly lower angular resolution compared to those used in Galloway-Sprietsma et al. (2025) and Law et al. (2022, 2023, 2024). Therefore, the differences observed in the ⟨z/r⟩ distributions cannot be solely attributed to the compactness of the disks considered in our analysis.

In Table 5, we report all the values of ⟨z/r⟩ used in Fig. 7. In Sect. 5, we show and discuss the correlation between ⟨z/r⟩ and R_90%.

Fig. 7 Kaplan–Meier estimator of radially averaged 12CO ⟨z/r⟩ values for our sample (in blue), compared with those of Galloway-Sprietsma et al. (2025) (in orange), and the combined works of Law et al. (2022, 2023, 2024) (in red).

5 Discussion

5.1 Brightness temperature – stellar luminosity correlation

We report a correlation between the brightness temperature of the ¹²CO layer at a given radius and the stellar luminosity. In the following, we chose a distance of 10 au from the central star to illustrate the results, as shown in Fig. 8. Since we are not directly constraining the ¹²CO temperature at 10 au, but rather inferring its value from a power-law parametric fit, we stress that the correlation should be interpreted as tracing how stellar luminosity affects the CO temperature structure in the outer disk, rather than a direct measurement of the temperature at 10 au. While here we choose to illustrate the results at 10 au because it is the quantity directly fit by the model, we verified that the correlation remains present when evaluating the temperature at larger radii (50–100 au). The correlation is log(T_b,10) = (2.10 ± 0.05) + (0.24 ± 0.07) log(L_∗), with a Pearson correlation coefficient of 0.6 in a range [0.4, 0.7]. This correlation is expected because more luminous stars irradiate the nearby circumstellar disk more strongly and therefore make it hotter. However, its slope and intercept have not been reported in the literature so far due to the lack of spatially resolved brightness temperature estimates on a large enough sample. We also compared these results to the total continuum flux analyzed in Carpenter et al. (2025), but found no obvious trend.

Additionally, we checked the presence of a correlation between T_b,10 and R_90% and found a positive correlation of T_b,10 = (70 ± 12) + (0.2 ± 0.1) R_90% at a 3σ confidence level, and a Pearson coefficient of r = 0.4. This hints toward a dependence of the CO brightness temperature on disk size, possibly reflecting underlying variations in disk mass and offering an alternative explanation in which the T_b,10 − L_∗ correlation reflects the correlation between the stellar luminosity and the stellar mass, which is in turn correlated to the disk mass. However, we found no empirical evidence in the literature for the presence of a correlation between gas disk size and gas disk mass. Pending a more detailed investigation on a larger sample that can disentangle the two explanations, in the following we proceed with the former since it provides a physical reason for the correlation we report.

The temperature values of the ¹²CO optically thick layer at a distance of 10 au from the central star were then compared to the temperature values of different dust grains at the same distance. First, as a reference, the red line shown in Fig. 8 represents the expected theoretical correlation of the equilibrium temperature, T_eq, at 10 au against the stellar luminosity. The equilibrium temperature is the temperature of a perfect blackbody dust grain irradiated by the central star: in this case, the slope of the correlation is 0.25. However, it is well-known that dust grains are not perfect blackbodies and that they absorb more efficiently than they re-emit (e.g., Draine 2003). For this reason, they undergo so-called superheating (i.e., not perfectly gray grains; see Chiang & Goldreich 1997). Second, we evaluated the expected temperature of the optically thin-to-stellar irradiation-superheated dust layer reported in Eq. (7) of Dullemond et al. (2001). We adopted a pure silicate dust composition for the Planck mean opacities¹⁰, evaluated using the dsharp_opac¹¹ package presented in Birnstiel et al. (2018): the result is shown as a dashed black line in Fig. 8.

There are two considerations we can make. The first is that the dashed black line could easily be interpreted as one of the gray best-fit lines in the background; this means that the dust at the optically thin-to-stellar-irradiation layer at 10 au has a similar temperature with respect to the optically thick ¹²CO layer that we extracted. The second is that the dust temperatures evaluated as a function of stellar luminosity are, on average, higher with respect to the best-fit correlation (solid black line). Since proto-planetary disks have a stratified temperature structure (Chiang & Goldreich 1997; Öberg et al. 2021), our results empirically imply that, on average, the optically thin-to-stellar-irradiation dust layer lies higher up than the ¹²CO optically thick emission layer. A similar result was reported by Law et al. (2024) and references therein, where a comparison between scattered-light images and the ¹²CO and ¹³CO layers is presented. We emphasize that, a priori, the stellar-irradiated layer and the CO emitting height could be totally different, given that they have vastly different opacities and the former is set by the optical depth from the star (e.g., D’Alessio et al. 1998, Bolchini et al., in prep.), while the latter comes from the observer (e.g., Rosotti et al. 2025). The fact that we find they are similar (albeit with a slightly lower ¹²CO emitting layer) is an empirical finding. Anecdotally, this is often found in radiative transfer calculations, but we are not aware of theoretical works that explicitly looked at the relative heights of these layers.

Finally, the slope of the observed correlation can be used to empirically measure how the Planck mean opacities scale as a function of temperature. If we suppose that the Planck mean opacities scale as a linear function of the temperature (κ_P ∝ αT^β), after some simple algebra, Eq. (7) of Dullemond et al. (2001) reads $$T_S^{4+\beta }=T_{*}^{\beta }{L}_{*}\frac{1}{16\pi \sigma R^2},$$(11) where T_s is the temperature of the optically thin-to-stellar-irradiation layer, T_∗ is the stellar surface temperature, σ is the Stefan-Boltzmann constant, and R is the distance from the stellar surface. From Eq. (11), we evaluated the exponent, β, which allowed us to retrieve the slope we evaluated for the correlation shown in Fig. 8; from this we found β ∼ 0.2. As shown in Fig. 8 of Birnstiel et al. (2018), the Planck mean opacities are expected to increase as a function of temperature: therefore, a slope β ∼ 0.2 is allowed. Moreover, in the temperature range of 10–200 K, Fig. 8 of Birnstiel et al. (2018) shows that the mean dependence of the Planck mean opacities from the temperature is κ_P ∝ T^0.45, the exponent of which is comparable to our β ∼ 0.2, even though it is not exactly equal. This result confirms the fact that the ¹²CO optically thick emission layer is more deeply embedded in the protoplanetary disk structure with respect to the optically thin to stellar irradiation dust layer. Finally, considering Fig. 8, we report that the median slope of the black dashed line is comparable with that of the solid black line (slope ∼0.22). Although not particularly surprising, the existence of this correlation offers empirical evidence that the basic dust models we used are reasonable.

Fig. 8 Correlation between 12CO brightness temperature at a distance of 10 au from the central star and stellar luminosity. The values of Tb,10 and log10(L*) were taken from Tables A.1 and A.2 and Empey et al., in prep., respectively. The values of Tb,10 were converted from cylindrical coordinates to spherical coordinates using the vertical-height-profile values reported in Table A.1. The values of Fcont were taken from Carpenter et al. (2025).

5.2 ⟨z/r⟩ – disk size correlation

Law et al. (2021) reported a correlation between ⟨z/r⟩ and the gas-disk size R_90%. We find that this still holds for our sample of compact disks orbiting T Tauri stars, as shown in Fig. 9. Here, we show ⟨z/r⟩ of our sample of disks as a function of R_90%.

The solid black curve represents the best fit, and in blue we over-plot the best-fit curve found by Galloway-Sprietsma et al. (2025) for comparison. We also color-code the data points by their integrated fluxes to show how the most vertically extended disks are, on average, also the brightest. We removed J16042165-2130284 from the sample, as it is found to be a face-on disk (inc ∼6^◦), for which the vertical-height extraction is very uncertain. We used J16142091-1906051, which is a known binary system, and J15583620-1946135 and J16132190-2136136, for which the resolution of the disk limits the derivation of reliable ⟨z/r⟩ values, as upper limits. We stress that the correlation would still be significant if we had not made such arrangements. The correlation is ⟨z/r⟩ = (0.05 ± 0.03) + (0.0009 ± 0.0003) R_90%, and it is found with a Pearson correlation coefficient of 0.6. Moreover, the fraction of posterior samples with an intercept >0 is ∼0.999, indicating a significance of ≳ 3σ. This result is in agreement within the 68th percentile with the results shown in Galloway-Sprietsma et al. (2025), and it shows that this correlation, to date found only for extended, massive, and bright protoplanetary disks, holds for fainter and more compact disks. Therefore, we have reason to think that the correlation between ⟨z/r⟩ and R_90% suggests something fundamental about all protoplanetary disks. However, it is not yet clear how this correlation is intimately related to the disk surface density. Paneque-Carreño et al. (2025) showed how the vertical structure of the disk as traced by CO is very sensitive to variations in disk mass by running thermochemical DALI models (Bruderer 2013) and accounting for hydrostatic equilibrium. More massive disks would therefore appear to be more extended both radially and vertically, and this would explain the observed correlation. However, CO depletion from the canonical ISM value (∼10⁻⁴) can also affect the vertical extent of the emitting layer (Paneque-Carreño et al. 2025; Rosotti et al. 2025), and the radial cutoff may be affected by the surface-density profile (Trapman et al. 2020).

Fig. 9 Correlation between 12CO ⟨z/r⟩ and R90%. The values of ⟨z/r⟩ were evaluated as discussed in Sect. 3, while the values of R90% were taken from Table A.1, and the integrated fluxes from Table 4.

5.3 Multiple systems

Among the sources we selected, we identified a group of multiple systems that show different characteristics. J16142091-1906051 shows a faint companion marginally visible in the continuum emission (see Carpenter et al. 2025 for the continuum gallery), while in the gas a spiral feature is clearly visible (see Fig. 10). However, the keplerian velocity pattern for the main source is perfectly visible and still well behaved: therefore, we were able to successfully run our analysis for it, and we list the results in Table A.1, and show the residuals in Fig. B.2f.

J16193570-1950426, shows evidence for what is likely to be an emission inflow/outflow (v = 1.5 km s⁻¹), which made it impossible to fit this source with the prescriptions adopted in this work (see Fig. 11). For this reason, we report the possibility that this source hides at least one hidden wide separated companion, which might be connected to the main star through the diffuse material. This source was not listed as a binary in Carpenter et al. (2025) or anywhere else, but was reported as a continuum burst-accreting source from Cody et al. (2017).

J16113134-1838259 (also known as AS 205, see e.g., Andrews et al. 2018) shows what is very likely to be a fly-by (see Fig. 12). For the latter, we tried to fit the two sources separately, but due to the lack of spectral resolution, we were unable to achieve a good fit.

In Fig. 13, we show the channel maps of J16185382-2053182. This system is marginally resolved, which makes it suitable for the analysis that we performed with csalt, and we show the residuals in Fig. B.3d. The results of the fit are listed in Table A.1. This interesting source with a spectral type of M0 was originally removed from the sample of Carpenter et al. (2025), as the disk classification of the star may be contaminated by another source in close proximity. We report the fact that we derived a stellar mass of ∼2.5 M_⊙ for this object, a result which is very discrepant with its associated spectral type and that could hint toward the presence of hidden companions.

Fig. 10 Channel maps of the binary system J16142091-1906051. The binary companion is placed at a RA offset of ∼+1″.

Fig. 11 Channel maps of J16193570-1950426. A bridge of extended emission is visible in the channel maps at v = 1.5 km s−1, and a clear spiral pattern is visible from the channel map at v = −0.4 km s−1 to the one at v = 6.1 km s−1.

Fig. 12 Channel maps of the binary system J16113134-1838259. A big tail of material is visible in the channel maps from v = 3.5–5.5 km s−1, showing that the two protoplanetary disks are actively interacting.

Fig. 13 Channel maps of J16185382-2053182. Blobs of extended emission are visible in the channel maps at v = 2.3, 4, 5.7, 12.5, and 14.2 km s−1.

6 Conclusions

In this work, we fit the visibilities of the uniformly surveyed Upper Scorpius SFR following the sample of Carpenter et al. (2025) and using two new parametric prescriptions to extract the disk bulk properties. For this scope, we used the software csalt (Andrews et al., in prep.) to model visibility spectra using two new parametric prescriptions. From our analysis, we derive the following conclusions:

• (i)

  We fit a total of 37 protoplanetary disks out of a sample of 83 ¹²CO J = 3–2 detection, and we report that visibility fitting becomes difficult if the S/N of the integrated flux is ≲10;

• (ii)

  We report the correlation between the ¹²CO brightness temperature at a given radius and the ¹²CO R_90%. We used this correlation to understand that the ¹²CO optically thick layer is more deeply embedded in the protoplanetary disk structure than the superheated-dust optically thin-to-stellar-irradiation layer;

• (iii)

  We extracted the ¹²CO J = 3–2-emitting layers and ⟨z/r⟩ values for “common” protoplanetary disks, deriving a median ⟨z/r⟩ ∼ 0.16 for our sample, which is smaller than what was reported for larger and brighter disks (Galloway-Sprietsma et al. 2025). Moreover, we find that the correlation between ⟨z/r⟩ and R_90% found by Law et al. (2021) for extended and bright disks still holds for the compact ones.

Acknowledgements

We thank the anonymous referee for their important suggestions which improved the quality and clarity of the paper. L.Z. is grateful to S. Facchini for some important informal suggestions for this work. L.Z. and G.R. acknowledge support from the European Union (ERC Starting Grant DiscEvol, project number 101039651), and from Fondazione Cariplo, grant No. 2022-1217. C.F.M. is funded by the European Union (ERC, WANDA, 101039452). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. T.P.C. was supported by the Heising-Simons Foundation through a 51 Pegasi b Fellowship. Support for C.J.L. was provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51535.001-A. Computational resources were provided by INDACO Platform, which is a project of High Performance Computing at the University of MILAN, and CINECA (http://www.unimi.it, https://www.cineca.it). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00526.S, ADS/JAO.ALMA#2012.1.00688.S, ADS/JAO.ALMA#2013.1.00395.S, ADS/-JAO.ALMA#2018.1.00564.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC 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. This paper also makes use of the following VLT/X-Shooter spectra: #097.C-0378, #0101.C-0866, #105.2082.003, #113.26NN.001, #113.26NN.003, #115.27XL.001.

References

Appendix A Tables of best-fit parameters

In Tables A.1, A.2, we report the best-fit parameters derived from our analysis.

Appendix B Residual plots

In Figs. B.1, B.2, B.3, B.4, B.5, B.6, we show the residual plots we created from the best-fit models following the procedure reported in Sect. 3.

Fig. B.1 Residual plots of J15583620-1946135 (a), J15583692-2257153 (b), J16035793-1942108 (c), J16052157-1821412 (d), J16062383-1807183 (e), J16064794-1841437(f). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be found on GitHub4.

Fig. B.2 Residual plots of J16095933-1800090 (a), J16101264-2104446 (b), J16120668-3010270 (c), J16132190-2136136 (d), J16140792-1938292 (e), J16142091-1906051(f). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be found on GitHub4.

Fig. B.3 Residual plots of J16145024-2100599 (a), J16152752-1847097 (b), J16181445-2319251 (c), J16185382-2053182 (d), J16202291-2227041 (e), J16202863-2442087(f). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be found on GitHub4.

Fig. B.4 Residual plots of J16203960-2634284 (a), J16215472-2752053 (b), J16221532-2511349 (c), J16222982-2002472 (d), J16230761-2516339 (e), J16253798-1943162(f). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be foundon GitHub4.

Fig. B.5 Residual plots of J16271273-2504017 (a), J16274905-2602437 (b), J16293267-2543291 (c), J16395577-2347355 (d), J16142029-1906481 (e), J16042165-2130284(f). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be found on GitHub4.

Fig. B.6 Residual plots of J16012268-2408003 (a), J16145026-2332397 (b), J16020757-2257467 (c), J16121242-1907191 (d). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be found on GitHub4.

Appendix C The intrinsic dependence of the gas-disk size from its inclination and from the projected velocity field

The disk geometrical framework of csalt is designed for inference. The disk model is created in the disk reference frame in cylindrical coordinates and is then projected into the Carte-sian plane of the sky for comparison with data. As a result, the parameters that describe the disk are generated and fit to the disk reference frame. Therefore, as described in Sec. 3, the most convenient way to evaluate quantities such as R_90% is to directly integrate the best-fit peak intensity profile in cylindrical coordinates. This approach allows one to get rid of the ‘deprojection problem’: as explained in Sec. 3, the usual way to extract the gas-disk radii is based on velocity-integrated (moment 0) maps. However, the moment 0 integrated intensities use the projected velocity field as the weighting function: therefore, same disk properties will result in different integrands depending on the value of inc and M_∗. Moreover, protoplanetary disks are flared: as a result of this, the disk is seen to be larger than it actually is when considering high inclinations, because the back surface becomes particularly relevant at large radii. While the latter effect can be mitigated from high resolution observations where the height of the emission surface can be well constrained, this is in general not possible for surveys and for this reason existing surveys have not corrected for this factor (e.g., Sanchis et al. 2021; Long et al. 2022; Trapman et al. 2025). In any case, the former effect remains unaccounted for.

We present an illustrative example in Fig. C.1. For this example, we generated five model cubes in the velocity range [−10; 20] km s⁻¹ with a spectral resolution of 0.429 km s⁻¹, to be consistent with the ALMA observations used in this work, using the dedicated framework of csalt. The model parameters were chosen following the best-fit model parameters of J16213469-2612269 (for which we derived ⟨z/r⟩ ∼ 0.14), which we showed in Sec. 3, and were kept fixed, while we set the inclination to five different values: 8°, 20°, 40°, 55°, 70°. Then, we created the integrated moment maps of our models and convolved them with the beam of our observations, which in the case of J16213469-2612269 was (0.4″, 0.3″), with PA = −81.5°. Then, using gofish, we extracted the radial profile of our convolved models and integrated it to derive R_90%, the radius that contains 90% of the total disk flux. Fig. C.1 shows the effect of extracting R_90% on velocity-integrated channel maps without keeping into consideration the flared nature of protoplanetary disks: the radius of a very inclined disk (70°) with an aspect ratio of 0.14 can be ∼40% greater than the true value.

Fig. C.1 Illustrative example of the intrinsic dependence between the gas-disk size and its inclination when extracting gas-disk radii from velocity-integrated maps.

Appendix D Data and model-integrated fluxes

To extract the fluxes from the data, we first generated moment 0 maps using a bandwidth of ±15 km s⁻¹ centered around the systemic velocity: this bandwidth allows us to include the entire disk emission for all sources. For this step, we used bettermoments (Teague & Foreman-Mackey 2018). Secondly, we evaluated geometric masks using the best-fit results reported in Tables A.1, A.2, and summed the emission within the mask for each moment map. The mask size was set to match R_90% + 0.2″. For the flux uncertainties, we moved the mask around the moment 0 maps in 8 different locations that were not enclosing any disk emission and not too far away from the image center, determined the flux in each of these locations and finally computed the standard deviation of the fluxes to obtain the integrated disk flux uncertainty. The flux from the models was computed in a similar way: we first generated datacubes using the best-fit parameters reported in Tables A.1, A.2 on the ±15 km s⁻¹ bandwidth, with the same spectral resolution as the observations. We then computed the disk fluxes using the same mask as for observations, and bootsrapped an uncertainty on the fluxes using the posterior distribution of the generated fluxes.

Appendix E Vertical heights

Fig. E.1 shows the vertical height profiles we derived in this work using model A (see Sect. 3) as a function of radius. We truncated each profile at a distance R_90% + 0.2″ from the central star.

Fig. E.1 Vertical heights of the 33 disks fit with the prescription A. We created the profiles according to Eq. (1) using all the random draws of the MCMC chains of z1, ψz, rz, and ϕz, evaluating the associated z(r) profiles, and taking their median value and the 16th and 84th percentiles for the associated uncertainties.

All Tables

All Figures

	Fig. 1 Residual plot for J16213469-2612269. In the first row, we show the channel maps of the data; in the second row, we show the channel maps of the best-fit model; in the third row, we show the channel maps of the residuals (data – model). The residual plots for all the disks analyzed in this work can be found in Appendix B and on GitHub4.
In the text

	Fig. 2 CCDF of 12CO R90% computed for this work (in blue) and of the radii evaluated using image-plane analysis from Zagaria et al., in prep. (in orange).
In the text

	Fig. 3 Residual plot for J16062861-2121297, shown in the same way as in Fig. 1. The residual plots for all the disks analyzed in this work can be found in Appendix B and on GitHub4.
In the text

Fig. 4 Residual plot for J16123916-1859284, shown in the same way as in Fig. 1. We highlight the residuals in the channel maps at v = 4.6, 5.5 km s−1, which are those affected by cloud absorption. Moreover, in the channel map at v = 3.7 km s−1, we show how csalt is able to separate the diffuse cloud emission (the localized red emission) and the Keplerian disk emission. The residual plots for all the disks analyzed in this work can be found in Appendix B and on GitHub4.

In the text

	Fig. 5 CCDF of 12CO R90% for our sample (in blue) and of the sources in Upper Scorpius reported from Trapman et al. (2025) (in orange). The gas-disk sizes reported in this figure are reported in Tables A.1 and A.2.
In the text

	Fig. 6 Correlation between 12CO gas-disk size and luminosity. The values for R90% and the disk fluxes, together with the line absorption ranges, can be found in Tables A.1, A.2, and 4, respectively.
In the text

	Fig. 7 Kaplan–Meier estimator of radially averaged 12CO ⟨z/r⟩ values for our sample (in blue), compared with those of Galloway-Sprietsma et al. (2025) (in orange), and the combined works of Law et al. (2022, 2023, 2024) (in red).
In the text

Fig. 8 Correlation between 12CO brightness temperature at a distance of 10 au from the central star and stellar luminosity. The values of Tb,10 and log10(L*) were taken from Tables A.1 and A.2 and Empey et al., in prep., respectively. The values of Tb,10 were converted from cylindrical coordinates to spherical coordinates using the vertical-height-profile values reported in Table A.1. The values of Fcont were taken from Carpenter et al. (2025).

In the text

	Fig. 9 Correlation between 12CO ⟨z/r⟩ and R90%. The values of ⟨z/r⟩ were evaluated as discussed in Sect. 3, while the values of R90% were taken from Table A.1, and the integrated fluxes from Table 4.
In the text

	Fig. 10 Channel maps of the binary system J16142091-1906051. The binary companion is placed at a RA offset of ∼+1″.
In the text

	Fig. 11 Channel maps of J16193570-1950426. A bridge of extended emission is visible in the channel maps at v = 1.5 km s−1, and a clear spiral pattern is visible from the channel map at v = −0.4 km s−1 to the one at v = 6.1 km s−1.
In the text

	Fig. 12 Channel maps of the binary system J16113134-1838259. A big tail of material is visible in the channel maps from v = 3.5–5.5 km s−1, showing that the two protoplanetary disks are actively interacting.
In the text

	Fig. 13 Channel maps of J16185382-2053182. Blobs of extended emission are visible in the channel maps at v = 2.3, 4, 5.7, 12.5, and 14.2 km s−1.
In the text

	Fig. B.1 Residual plots of J15583620-1946135 (a), J15583692-2257153 (b), J16035793-1942108 (c), J16052157-1821412 (d), J16062383-1807183 (e), J16064794-1841437(f). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be found on GitHub4.
In the text

	Fig. B.2 Residual plots of J16095933-1800090 (a), J16101264-2104446 (b), J16120668-3010270 (c), J16132190-2136136 (d), J16140792-1938292 (e), J16142091-1906051(f). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be found on GitHub4.
In the text

	Fig. B.3 Residual plots of J16145024-2100599 (a), J16152752-1847097 (b), J16181445-2319251 (c), J16185382-2053182 (d), J16202291-2227041 (e), J16202863-2442087(f). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be found on GitHub4.
In the text

	Fig. B.4 Residual plots of J16203960-2634284 (a), J16215472-2752053 (b), J16221532-2511349 (c), J16222982-2002472 (d), J16230761-2516339 (e), J16253798-1943162(f). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be foundon GitHub4.
In the text

	Fig. B.5 Residual plots of J16271273-2504017 (a), J16274905-2602437 (b), J16293267-2543291 (c), J16395577-2347355 (d), J16142029-1906481 (e), J16042165-2130284(f). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be found on GitHub4.
In the text

	Fig. B.6 Residual plots of J16012268-2408003 (a), J16145026-2332397 (b), J16020757-2257467 (c), J16121242-1907191 (d). The residual plots for all the disks analyzed in this work and their corresponding corner plots can be found on GitHub4.
In the text

	Fig. C.1 Illustrative example of the intrinsic dependence between the gas-disk size and its inclination when extracting gas-disk radii from velocity-integrated maps.
In the text

	Fig. E.1 Vertical heights of the 33 disks fit with the prescription A. We created the profiles according to Eq. (1) using all the random draws of the MCMC chains of z1, ψz, rz, and ϕz, evaluating the associated z(r) profiles, and taking their median value and the 16th and 84th percentiles for the associated uncertainties.
In the text