# Accuracy of ALMA estimates of young disk radii and masses

Ngo-Duy Tung, Leonardo Testi, Ugo Lebreuilly, Patrick Hennebelle, Anaëlle Maury, Ralf S. Klessen, Luca Cacciapuoti, Matthias González, Giovanni Rosotti, Sergio Molinari

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

## ► To cite this version:

Ngo-Duy Tung, Leonardo Testi, Ugo Lebreuilly, Patrick Hennebelle, Anaëlle Maury, et al.. Accuracy of ALMA estimates of young disk radii and masses: Predicted observations from numerical simulations. Astronomy and Astrophysics - A&A, 2024, 684, pp.A36. 10.1051/0004-6361/202348730. cea-04531552

HAL Id: cea-04531552

https://cea.hal.science/cea-04531552

Submitted on 3 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Accuracy of ALMA estimates of young disk radii and masses
## Predicted observations from numerical simulations

Ngo-Duy Tung<sup>1,2</sup> <sup>●</sup>, Leonardo Testi<sup>3,4</sup> <sup>●</sup>, Ugo Lebreuilly<sup>1</sup> <sup>●</sup>, Patrick Hennebelle<sup>1</sup> <sup>●</sup>, Anaëlle Maury<sup>1</sup>, Ralf S. Klessen<sup>5,6</sup> <sup>●</sup>, Luca Cacciapuoti<sup>7,4,8</sup> <sup>●</sup>, Matthias González<sup>2</sup> <sup>●</sup>, Giovanni Rosotti<sup>9,10</sup> <sup>●</sup>, and Sergio Molinari<sup>11</sup>

Received 25 November 2023 / Accepted 20 January 2024

## ABSTRACT

**Context.** Protoplanetary disks, which are the natural consequence of the gravitational collapse of the dense molecular cloud cores, host the formation of the known planetary systems in our universe. Substantial efforts have been dedicated to investigating the properties of these disks in the more mature Class II stage, either via numerical simulations of disk evolution from a limited range of initial conditions or observations of their dust continuum and line emission from specific molecular tracers. The results coming from these two standpoints have been used to draw comparisons. However, few studies have investigated the main limitations at work when measuring the embedded Class 0/I disk properties from observations, especially in a statistical fashion.

**Aims.** In this study, we provide a first attempt to compare the accuracy of some critical disk parameters in Class 0/I systems, as derived on real ALMA observational data, with the corresponding physical parameters that can be directly defined by theoreticians and modellers in numerical simulations. The approach we follow here is to provide full post-processing of the numerical simulations and apply it to the synthetic observations the same techniques used by observers to derive the physical parameters.

**Methods.** We performed 3D Monte Carlo radiative transfer and mock interferometric observations of the disk populations formed in a magnetohydrodynamic (MHD) simulation model of disk formation through the collapse of massive clumps with the tools RADMC-3D and CASA, respectively, to obtain their synthetic observations. With these observations, we re-employed the techniques commonly used in disk modelling from their continuum emissions to infer the properties that would most likely be obtained with real interferometers. We then demonstrated how these properties may vary with respect to the gas kinematics analyses and dust continuum modelling.

**Results.** Our modelling procedure, based on a two-component model for the disk and the envelope, shows that the disk sizes can be properly recovered from observations with sufficient angular resolutions, with an uncertainty of a factor ≈1.6–2.2, whereas their masses cannot be accurately measured. Overall, the masses are predominantly underestimated for larger, more massive disks by a median factor of ≈2.5, and even up to 10 in extreme cases, with the conversion from flux to dust mass under the optically thin assumption. We also find that the single Gaussian fittings are not a reliable modelling technique for young, embedded disks characterised by a strong presence of the envelopes. Thus, such an approach is to be used with caution.

**Conclusions.** The radiative transfer post-processing and synthetic observations of MHD simulations offer genuine help in linking important observable properties of young planet-forming disks to their intrinsic values in simulations. Further extended investigations that tackle the caveats of this study, such as the lack of variation in the dust composition and distribution, dust-to-gas ratio, and other shortcomings in the numerical models, would be essential for setting constraints on our understanding of disk and planet formations.

**Key words.** protoplanetary disks – circumstellar matter – submillimeter: general

## 1. Introduction

The formation of stars and planetary systems starts within a dense molecular cloud core: the high-density region of the cloud where the majority of star formation activities take place thanks

to the dominance of self-gravity over gas pressure. This dominance also gives rise to gravitational collapse of the core. In a simplistic picture, the angular momentum conservation of the core during the collapse causes most of the infalling matter to form a circumstellar disk around the protostar, instead of falling

<footnote n="1"><sup>1</sup> Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France
e-mail: *duy-tung.ngo@cea.fr*</footnote>

<footnote n="2"><sup>2</sup> Université Paris Cité, Université Paris-Saclay, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France
e-mail: *duy-tung.ngo@etu.u-paris.fr*</footnote>

<footnote n="3"><sup>3</sup> Alma Mater Studiorum, Università di Bologna, Dipartimento di Fisica e Astronomia (DIFA), Via Gobetti 93/2, 40129 Bologna, Italy</footnote>

<footnote n="4"><sup>4</sup> INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy</footnote>

<footnote n="5"><sup>5</sup> Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Üeberle-Str. 2, 69120 Heidelberg, Germany</footnote>

<footnote n="6"><sup>6</sup> Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany</footnote>

<footnote n="7"><sup>7</sup> European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="8"><sup>8</sup> Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany</footnote>

<footnote n="9"><sup>9</sup> Dipartimento di Fisica, Università degli Studi di Milano, Via Giovanni Celoria, 16, 20133 Milano, Italy</footnote>

<footnote n="10"><sup>10</sup> Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands</footnote>

<footnote n="11"><sup>11</sup> INAF-Istituto di Astrofisica e Planetologia Spaziali (INAF-IAPS), Via Fosso del Cavaliere 100, 00133 Roma, Italy</footnote>
# Accuracy of ALMA estimates of young disk radii and masses

## Predicted observations from numerical simulations

Ngo-Duy Tung<sup>1,2</sup>①, Leonardo Testi<sup>3,4</sup>①, Ugo Lebreuilly<sup>1</sup>①, Patrick Hennebelle<sup>1</sup>①, Anaëlle Maury<sup>1</sup>,
Ralf S. Klessen<sup>5,6</sup>①, Luca Cacciapuoti<sup>7,4,8</sup>①, Matthias González<sup>2</sup>①, Giovanni Rosotti<sup>9,10</sup>①, and Sergio Molinari<sup>11</sup>

Received 25 November 2023 / Accepted 20 January 2024

## ABSTRACT

**Context.** Protoplanetary disks, which are the natural consequence of the gravitational collapse of the dense molecular cloud cores, host the formation of the known planetary systems in our universe. Substantial efforts have been dedicated to investigating the properties of these disks in the more mature Class II stage, either via numerical simulations of disk evolution from a limited range of initial conditions or observations of their dust continuum and line emission from specific molecular tracers. The results coming from these two standpoints have been used to draw comparisons. However, few studies have investigated the main limitations at work when measuring the embedded Class 0/I disk properties from observations, especially in a statistical fashion.

**Aims.** In this study, we provide a first attempt to compare the accuracy of some critical disk parameters in Class 0/I systems, as derived on real ALMA observational data, with the corresponding physical parameters that can be directly defined by theoreticians and modellers in numerical simulations. The approach we follow here is to provide full post-processing of the numerical simulations and apply it to the synthetic observations the same techniques used by observers to derive the physical parameters.

**Methods.** We performed 3D Monte Carlo radiative transfer and mock interferometric observations of the disk populations formed in a magnetohydrodynamic (MHD) simulation model of disk formation through the collapse of massive clumps with the tools RADMC-3D and CASA, respectively, to obtain their synthetic observations. With these observations, we re-employed the techniques commonly used in disk modelling from their continuum emissions to infer the properties that would most likely be obtained with real interferometers. We then demonstrated how these properties may vary with respect to the gas kinematics analyses and dust continuum modelling.

**Results.** Our modelling procedure, based on a two-component model for the disk and the envelope, shows that the disk sizes can be properly recovered from observations with sufficient angular resolutions, with an uncertainty of a factor ≈1.6–2.2, whereas their masses cannot be accurately measured. Overall, the masses are predominantly underestimated for larger, more massive disks by a median factor of ≈2.5, and even up to 10 in extreme cases, with the conversion from flux to dust mass under the optically thin assumption. We also find that the single Gaussian fittings are not a reliable modelling technique for young, embedded disks characterised by a strong presence of the envelopes. Thus, such an approach is to be used with caution.

**Conclusions.** The radiative transfer post-processing and synthetic observations of MHD simulations offer genuine help in linking important observable properties of young planet-forming disks to their intrinsic values in simulations. Further extended investigations that tackle the caveats of this study, such as the lack of variation in the dust composition and distribution, dust-to-gas ratio, and other shortcomings in the numerical models, would be essential for setting constraints on our understanding of disk and planet formations.

**Key words.** protoplanetary disks – circumstellar matter – submillimeter: general

## 1. Introduction

The formation of stars and planetary systems starts within a dense molecular cloud core: the high-density region of the cloud where the majority of star formation activities take place thanks

to the dominance of self-gravity over gas pressure. This dominance also gives rise to gravitational collapse of the core. In a simplistic picture, the angular momentum conservation of the core during the collapse causes most of the infalling matter to form a circumstellar disk around the protostar, instead of falling

<footnote n="1"><sup>1</sup> Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France
e-mail: *duy-tung.ngo@cea.fr*</footnote>

<footnote n="2"><sup>2</sup> Université Paris Cité, Université Paris-Saclay, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France
e-mail: *duy-tung.ngo@etu.u-paris.fr*</footnote>

<footnote n="3"><sup>3</sup> Alma Mater Studiorum, Università di Bologna, Dipartimento di Fisica e Astronomia (DIFA), Via Gobetti 93/2, 40129 Bologna, Italy</footnote>

<footnote n="4"><sup>4</sup> INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy</footnote>

<footnote n="5"><sup>5</sup> Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-UEberle-Str. 2, 69120 Heidelberg, Germany</footnote>

<footnote n="6"><sup>6</sup> Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany</footnote>

<footnote n="7"><sup>7</sup> European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="8"><sup>8</sup> Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany</footnote>

<footnote n="9"><sup>9</sup> Dipartimento di Fisica, Università degli Studi di Milano, Via Giovanni Celoria, 16, 20133 Milano, Italy</footnote>

<footnote n="10"><sup>10</sup> Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands</footnote>

<footnote n="11"><sup>11</sup> INAF-Istituto di Astrofisica e Planetologia Spaziali (INAF-IAPS), Via Fosso del Cavaliere 100, 00133 Roma, Italy</footnote>
directly onto it (see Shu et al. 1987; Williams & Cieza 2011; Tsukamoto et al. 2023 for reviews). This disk of gas and dust is often referred to as the protostellar disk in the broad sense and also referred to as protoplanetary disk (hereafter PPD) in reference to its role as a (proto-)planet formation site. Throughout all the evolutionary stages of PPDs (and most evident in the later ones), dust grains play a central role in the evolution of PPDs as they grow into larger aggregates, then centimeter-sized pebbles, which eventually become planetesimals, followed by kilometer-sized bodies that make up the building blocks of planetary systems (Testi et al. 2014).

Based on their spectral energy distributions (SEDs), the evolution of such star-disk systems can be classified into different stages, ranging from the cold, dense pre-stellar core phase to Class 0-I protostellar phase, where the transition from cores to disks happens. This is followed by the pre-main sequence Class II-II disks, where most of the planet formation activities are widely believed to take place (Andre et al. 2000).

So far, much of our understanding of the evolution of PPDs and planet formation is based on theory and observations of isolated, more evolved Class II objects. However, our Sun, as most stars, was not initially isolated; it formed in a stellar cluster environment (Adams 2010; Pfalzner et al. 2015), where the radiation field (photoevaporation) and (to a much lesser extent) tidal interactions from other forming stars strongly influenced the formation of the Earth and other planets in the Solar System (Winter et al. 2018). Thus, the need for modelling star and planet formation in massive star-forming complexes becomes ineluctable. Added to this is the fact that isolated simulations of disk formation and early evolution of disks from single core collapses in different environments often yield different outcomes regarding the disk properties (Hennebelle et al. 2020; Lee et al. 2021). Therefore, it is of great importance to model these processes starting from the large-scale environment and compare the outputs of these simulations with observations of young stellar objects (YSOs), with an aim to answer the question of whether the physical ingredients in these simulations are sufficient to represent the reality of disks in different galactic environments. In this light, Lebreuilly et al. (2021, 2024) recently performed simulations of protoplanetary disk formation, starting from the collapse of massive star-forming clumps to the formation of a cluster of YSOs in the Class 0/early Class I stage. These simulations open up new horizons for the comparison between models and observations of young disk populations, as they not only investigate the influence of various initial conditions on the self-consistent disk population within their parent environments, they also provide large samples of disks that form the basis for insightful statistical analyses.

On the observational side, the disk properties inferred from dust continuum and molecular line emission often suffer from uncertainties arising from noise, resolution limitations, and instrumental effects, as well as the modelling techniques to derive the physical parameters (see, e.g. Koepferl & Robitaille 2017). Another complication in the modelling of Class 0/I disks comes from the fact that these small embedded objects are still surrounded by a circumstellar envelope, which makes it difficult to separate the emission of the two components in sub-millimetric observations (Tobin et al. 2015; Segura-Cox et al. 2018; Maury et al. 2019). As a result, the comparison between simulated disk populations and their synthetic observations is also expected to help to test how closely the disk properties that we obtain from these modelling processes resemble their true properties.

The aim of the work presented in this paper is the first step toward our overarching effort of the Synthetic Populations of Protoplanetary Disks project (see also Lebreuilly et al. 2024). This endeavour is aimed at connecting the disk population properties from simulations of cluster-forming clouds to observations of young disk populations in our Galaxy obtained with the Atacama Large Millimeter/submillimeter Array (ALMA). We centre our study on the comparison methodology as well as the issue of how to design observations and data analyses to derive the disk parameters in embedded protostars. For our purposes, we post-processed an MHD simulation of protoplanetary disk population through the collapse of a 1000 $M_{\odot}$ clump in Lebreuilly et al. (2024) by means of continuum radiative transfer using the Monte Carlo radiative transfer code RADMC-3D (Dullemond et al. 2012). We produced realistic interferometric observation simulations using the Common Astronomy Software Applications (CASA) tool (McMullin et al. 2007; CASA Team et al. 2022), which provides us with a simulator for ALMA. We then treated them as real observational data and employed the modelling methods from recent studies and surveys (Maury et al. 2019; Tazzari et al. 2021a,b) to extract the disks' parameters.

The paper is structured as follows. In Sect. 2, we present the simulation model and snapshot used for the demonstration of our methodology. Sections 3 and 4 provide the technical details for the radiative transfer and interferometric post-processing, respectively. The modelling of the synthetic observations is described in Sect. 5, followed by an in-depth analysis of the modelled data regarding the disk properties in Sect. 6. We then spend Sect. 7 to discuss various implications of our results on the modelling of Class 0 disks from real observations. A summary of our finding is given in Sect. 8.

## 2. Simulation models

We begin by briefly presenting the simulation model (NMHD-F01 of Lebreuilly et al. 2024) used for the study in this paper. More details can be found in the corresponding paper as well as in Lebreuilly et al. (2021). The simulation is performed with the adaptive mesh refinement (AMR; Berger & Oliger 1984) code RAMSES (Teyssier 2002; Fromang et al. 2006), which solves the MHD equations using the finite volume method and enables us to resolve up to the disk formation scale with a maximum resolution of $\Delta x_{\max} \sim 1.2$ au. The starting point is a 1000 $M_{\odot}$ clump of temperature 10 K, initial radius $\sim 0.38$ pc, and uniform density $3 \times 10^{-19}$ g cm$^{-3}$, initialised with supersonic turbulent velocities at $M = 7$ with random phases and a Kolmogorov power spectrum of $k^{-11/3}$. The initially vertical magnetic field is set according to the mass-to-flux to critical-mass-to-flux ratio such as $\mu = 10$, which corresponds to a magnetic field of $\sim 9.4 \times 10^{-5}$ G. Sink particles (Bleuler & Teyssier 2014), which form when the density reaches $n_{\text{thre}} = 10^{13}$ cm$^{-3}$, were used as sub-grid models to account for the presence of fully formed stars. These sinks were not allowed to merge.

There are two types of luminosity associated with these sinks. The first one, accretion luminosity, arises when a star of mass, $M_{\star}$, and radius, $R_{\star}$, accretes mass and the incoming kinetic energy of the gas is radiated away. This luminosity has a value of

$$L_{\text{acc}} = f_{\text{acc}} \frac{GM_{\star} \dot{M}_{\star}}{R_{\star}}, \qquad (1)$$
directly onto it (see Shu et al. 1987; Williams & Cieza 2011; Tsukamoto et al. 2023 for reviews). This disk of gas and dust is often referred to as the protostellar disk in the broad sense and also referred to as protoplanetary disk (hereafter PPD) in reference to its role as a (proto-)planet formation site. Throughout all the evolutionary stages of PPDs (and most evident in the later ones), dust grains play a central role in the evolution of PPDs as they grow into larger aggregates, then centimeter-sized pebbles, which eventually become planetesimals, followed by kilometer-sized bodies that make up the building blocks of planetary systems (Testi et al. 2014).

Based on their spectral energy distributions (SEDs), the evolution of such star-disk systems can be classified into different stages, ranging from the cold, dense pre-stellar core phase to Class 0-I protostellar phase, where the transition from cores to disks happens. This is followed by the pre-main sequence Class II-II disks, where most of the planet formation activities are widely believed to take place (Andre et al. 2000).

So far, much of our understanding of the evolution of PPDs and planet formation is based on theory and observations of isolated, more evolved Class II objects. However, our Sun, as most stars, was not initially isolated; it formed in a stellar cluster environment (Adams 2010; Pfalzner et al. 2015), where the radiation field (photoevaporation) and (to a much lesser extent) tidal interactions from other forming stars strongly influenced the formation of the Earth and other planets in the Solar System (Winter et al. 2018). Thus, the need for modelling star and planet formation in massive star-forming complexes becomes ineluctable. Added to this is the fact that isolated simulations of disk formation and early evolution of disks from single core collapses in different environments often yield different outcomes regarding the disk properties (Hennebelle et al. 2020; Lee et al. 2021). Therefore, it is of great importance to model these processes starting from the large-scale environment and compare the outputs of these simulations with observations of young stellar objects (YSOs), with an aim to answer the question of whether the physical ingredients in these simulations are sufficient to represent the reality of disks in different galactic environments. In this light, Lebreuilly et al. (2021, 2024) recently performed simulations of protoplanetary disk formation, starting from the collapse of massive star-forming clumps to the formation of a cluster of YSOs in the Class 0/early Class I stage. These simulations open up new horizons for the comparison between models and observations of young disk populations, as they not only investigate the influence of various initial conditions on the self-consistent disk population within their parent environments, they also provide large samples of disks that form the basis for insightful statistical analyses.

On the observational side, the disk properties inferred from dust continuum and molecular line emission often suffer from uncertainties arising from noise, resolution limitations, and instrumental effects, as well as the modelling techniques to derive the physical parameters (see, e.g. Koepferl & Robitaille 2017). Another complication in the modelling of Class 0/I disks comes from the fact that these small embedded objects are still surrounded by a circumstellar envelope, which makes it difficult to separate the emission of the two components in sub-millimetric observations (Tobin et al. 2015; Segura-Cox et al. 2018; Maury et al. 2019). As a result, the comparison between simulated disk populations and their synthetic observations is also expected to help to test how closely the disk properties that we obtain from these modelling processes resemble their true properties.

The aim of the work presented in this paper is the first step toward our overarching effort of the Synthetic Populations of Protoplanetary Disks project (see also Lebreuilly et al. 2024). This endeavour is aimed at connecting the disk population properties from simulations of cluster-forming clouds to observations of young disk populations in our Galaxy obtained with the Atacama Large Millimeter/submillimeter Array (ALMA). We centre our study on the comparison methodology as well as the issue of how to design observations and data analyses to derive the disk parameters in embedded protostars. For our purposes, we post-processed an MHD simulation of protoplanetary disk population through the collapse of a 1000 $M_{\odot}$ clump in Lebreuilly et al. (2024) by means of continuum radiative transfer using the Monte Carlo radiative transfer code RADMC-3D (Dullemond et al. 2012). We produced realistic interferometric observation simulations using the Common Astronomy Software Applications (CASA) tool (McMullin et al. 2007; CASA Team et al. 2022), which provides us with a simulator for ALMA. We then treated them as real observational data and employed the modelling methods from recent studies and surveys (Maury et al. 2019; Tazzari et al. 2021a,b) to extract the disks' parameters.

The paper is structured as follows. In Sect. 2, we present the simulation model and snapshot used for the demonstration of our methodology. Sections 3 and 4 provide the technical details for the radiative transfer and interferometric post-processing, respectively. The modelling of the synthetic observations is described in Sect. 5, followed by an in-depth analysis of the modelled data regarding the disk properties in Sect. 6. We then spend Sect. 7 to discuss various implications of our results on the modelling of Class 0 disks from real observations. A summary of our finding is given in Sect. 8.

## 2. Simulation models

We begin by briefly presenting the simulation model (NMHD-F01 of Lebreuilly et al. 2024) used for the study in this paper. More details can be found in the corresponding paper as well as in Lebreuilly et al. (2021). The simulation is performed with the adaptive mesh refinement (AMR; Berger & Oliger 1984) code RAMSES (Teyssier 2002; Fromang et al. 2006), which solves the MHD equations using the finite volume method and enables us to resolve up to the disk formation scale with a maximum resolution of $\Delta x_{max} \sim 1.2$ au. The starting point is a 1000 $M_{\odot}$ clump of temperature 10 K, initial radius $\sim 0.38$ pc, and uniform density $3 \times 10^{-19}$ g cm$^{-3}$, initialised with supersonic turbulent velocities at $M = 7$ with random phases and a Kolmogorov power spectrum of $k^{-11/3}$. The initially vertical magnetic field is set according to the mass-to-flux to critical-mass-to-flux ratio such as $\mu = 10$, which corresponds to a magnetic field of $\sim 9.4 \times 10^{-5}$ G. Sink particles (Bleuler & Teyssier 2014), which form when the density reaches $n_{thre} = 10^{13}$ cm$^{-3}$, were used as sub-grid models to account for the presence of fully formed stars. These sinks were not allowed to merge.

There are two types of luminosity associated with these sinks. The first one, accretion luminosity, arises when a star of mass, $M_{\star}$, and radius, $R_{\star}$, accretes mass and the incoming kinetic energy of the gas is radiated away. This luminosity has a value of

$$L_{\text{acc}} = f_{\text{acc}} \frac{GM_{\star} \dot{M}_{\star}}{R_{\star}}, \qquad (1)$$
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Column density map integrated along the z-direction of the NMHD-F01 model from Lebreuilly et al. (2024) at the simulation snapshot used in this study, seen at a full simulation box's length. The star symbols indicate the sink particle positions.</figcaption></figure>

where $f_{\text{acc}}$ is the fraction of the accretion gravitational energy radiated away and $G$ is the gravitational constant. We adopt $f_{\text{acc}} = 0.1$ in the model to investigate the influence of the radiative feedback. The second source of luminosity comes from the stellar radiation. The star radius $R_{\star}$ and its internal luminosity $L_{\text{int}}$ are inferred using the stellar models from Kuiper & Yorke (2013).

The simulation is evolved until the total mass in all the sinks of the simulations reaches $M_{\text{sink}} = 150 M_{\odot}$. For the calculations presented in the following sections, we use a particular output of the simulations where the star formation efficiency (SFE) is ≈10%, which corresponds to $M_{\text{sink}} = 100 M_{\odot}$, and whose physical time is 110.95 kyr. Figure 1 shows the column density map of the simulation box (~1.2 pc in length) integrated along the z-direction at the chosen output.

### 3. Radiative transfer post-processing

From the simulation output, we performed 3D Monte Carlo (MC) radiative transfer calculations on the native AMR grid using the RADMC-3D code (Dullemond et al. 2012) to ultimately simulate the sky to be observed and obtain the intensity maps of the disk population at a resolution equivalent to the maximum resolution the simulation.

One of the most important properties of the dust and gas components in the simulation grid is their temperature, which is critical for the later derivation of their emissions. On the one hand, RAMSES's extension to radiative transfer in the flux limited diffusion (FLD) approximation (Levermore & Pomraning 1981; Commerçon et al. 2011, 2014) provides us with the temperature profile of the gas and dust in thermal equilibrium. On the other hand, RADMC-3D enables us to re-calculate the dust temperature using a full Monte Carlo approach, without the FLD approximation limitation. Therefore, it is useful to compare the results of the two methods to assess the accuracy of the temperature computation.

For the RADMC-3D MC thermal calculation (as well as the imaging process through ray-tracing), a dust-to-gas ratio of 0.01 is assumed, similarly to what is done in the simulation, to infer the dust density profile for RADMC-3D from the gas density calculated with RAMSES. Then, multiple sources of photons are input to RADMC-3D from the sink properties. In particular, the sources' luminosity are taken to be the sum of their internal and accretion luminosities $L_{\text{tot}} = L_{\text{acc}} + L_{\text{int}}$. These luminosities, along with other stellar properties namely the effective temperature and the star's radius, are computed in the simulation run. Another important input for the radiative transfer modelling is the dust opacity. Since RAMSES only works with Planck and Rosseland mean opacity which are not compatible with the wavelength-dependent opacity required for RADMC-3D, here we adopted the DIANA standard dust model and computed the opacity file for RADMC-3D using the tool for dust particle opacities calculations, OPTOOL (Dominik et al. 2021). This model comprises a specific pyroxene (70% Mg) and carbon, in a mass ratio of 0.87/0.13, and with a porosity of 25% (Woitke et al. 2016). An MRN power-law sized distribution with the minimum grain size $a_{\text{min}} = 1$ nm and the maximum size $a_{\text{max}} = 10$ µm was assumed for this dust model. To ease the computation cost, we treated this size distribution as one single dust population in the opacity file.

We then used both temperature profiles to carry out the disk size analyses. We give a comparison of the results in Appendix A. In general, we find that the two temperature profiles do not differ significantly in the outer region of the disks (where $r \gtrsim 10$ au), which is expected if the temperature is mainly determined by the reprocessing of radiation, unless there is substantial dissipation via turbulence; namely, high turbulence could result in strong local heating. In a standard model of an accretion disk, most of the viscous heating takes place in the inner region of the disk and the outer disk is always dominated by passive heating. As a result, the disk parameters obtained with both profiles are quite similar. However, the temperature obtained with RADMC-3D sometimes introduces over-heated region surrounding the photon source which results in a spike in the central flux and causes difficulties in the subsequent modelling of the post-processed results, whereas in RAMSES, the inclusion of a smoothing kernel for the photon flux within the ≈5 au radius of the sink particle and dust sublimation maintains the consistency of the central intensity. Therefore, we chose to continue with the dust temperature from RAMSES.

Finally, we imaged all the disks present in the simulation snapshot by performing the ray-tracing through the grid to get their continuum emissions at ALMA's band 7 wavelength of $\lambda = 0.89$ mm. The imaging on the plane of the sky is done through 2D regularly-spaced cameras observing the 3D AMR grids and, thus, for practical computational reasons, we proceeded with the second-order ray-tracing by zooming on individual sinks particles whose disks are detected. We carried out the imaging procedure in sub-regions of the grid around the sink positions with an 1 au pixel resolution, which allowed us to emulate the maximum numerical resolution of ≈1.2 au in the most resolved regions of the disks. For each of these regions, $10^8$ scattering photon packages were used for the MC scattering calculations before the ray-tracing to get the intensity maps. We accounted for the variation of the viewing angle in observations by varying the camera positions to image the disks in three different planes of the internal simulation grid: xy, xz, and yz. We note that these projections do not necessarily correspond to the direction of the angular momentum vector of the disks, which are randomly oriented by nature.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Column density map integrated along the z-direction of the NMHD-F01 model from Lebreuilly et al. (2024) at the simulation snapshot used in this study, seen at a full simulation box's length. The star symbols indicate the sink particle positions.</figcaption></figure>

where $f_{\text{acc}}$ is the fraction of the accretion gravitational energy radiated away and $G$ is the gravitational constant. We adopt $f_{\text{acc}} = 0.1$ in the model to investigate the influence of the radiative feedback. The second source of luminosity comes from the stellar radiation. The star radius $R_{\star}$ and its internal luminosity $L_{\text{int}}$ are inferred using the stellar models from Kuiper & Yorke (2013).

The simulation is evolved until the total mass in all the sinks of the simulations reaches $M_{\text{sink}} = 150 M_{\odot}$. For the calculations presented in the following sections, we use a particular output of the simulations where the star formation efficiency (SFE) is ≈10%, which corresponds to $M_{\text{sink}} = 100 M_{\odot}$, and whose physical time is 110.95 kyr. Figure 1 shows the column density map of the simulation box (~1.2 pc in length) integrated along the z-direction at the chosen output.

### 3. Radiative transfer post-processing

From the simulation output, we performed 3D Monte Carlo (MC) radiative transfer calculations on the native AMR grid using the RADMC-3D code (Dullemond et al. 2012) to ultimately simulate the sky to be observed and obtain the intensity maps of the disk population at a resolution equivalent to the maximum resolution of the simulation.

One of the most important properties of the dust and gas components in the simulation grid is their temperature, which is critical for the later derivation of their emissions. On the one hand, RAMSES's extension to radiative transfer in the flux limited diffusion (FLD) approximation (Levermore & Pomraning 1981; Commerçon et al. 2011, 2014) provides us with the temperature profile of the gas and dust in thermal equilibrium. On the other hand, RADMC-3D enables us to re-calculate the dust temperature using a full Monte Carlo approach, without the FLD approximation limitation. Therefore, it is useful to compare the results of the two methods to assess the accuracy of the temperature computation.

For the RADMC-3D MC thermal calculation (as well as the imaging process through ray-tracing), a dust-to-gas ratio of 0.01 is assumed, similarly to what is done in the simulation, to infer the dust density profile for RADMC-3D from the gas density calculated with RAMSES. Then, multiple sources of photons are input to RADMC-3D from the sink properties. In particular, the sources' luminosity are taken to be the sum of their internal and accretion luminosities $L_{\text{tot}} = L_{\text{acc}} + L_{\text{int}}$. These luminosities, along with other stellar properties namely the effective temperature and the star's radius, are computed in the simulation run. Another important input for the radiative transfer modelling is the dust opacity. Since RAMSES only works with Planck and Rosseland mean opacity which are not compatible with the wavelength-dependent opacity required for RADMC-3D, here we adopted the DIANA standard dust model and computed the opacity file for RADMC-3D using the tool for dust particle opacities calculations, OPTOOL (Dominik et al. 2021). This model comprises a specific pyroxene (70% Mg) and carbon, in a mass ratio of 0.87/0.13, and with a porosity of 25% (Woitke et al. 2016). An MRN power-law sized distribution with the minimum grain size $a_{\text{min}} = 1$ nm and the maximum size $a_{\text{max}} = 10$ µm was assumed for this dust model. To ease the computation cost, we treated this size distribution as one single dust population in the opacity file.

We then used both temperature profiles to carry out the disk size analyses. We give a comparison of the results in Appendix A. In general, we find that the two temperature profiles do not differ significantly in the outer region of the disks (where $r \gtrsim 10$ au), which is expected if the temperature is mainly determined by the reprocessing of radiation, unless there is substantial dissipation via turbulence; namely, high turbulence could result in strong local heating. In a standard model of an accretion disk, most of the viscous heating takes place in the inner region of the disk and the outer disk is always dominated by passive heating. As a result, the disk parameters obtained with both profiles are quite similar. However, the temperature obtained with RADMC-3D sometimes introduces over-heated region surrounding the photon source which results in a spike in the central flux and causes difficulties in the subsequent modelling of the post-processed results, whereas in RAMSES, the inclusion of a smoothing kernel for the photon flux within the ≈5 au radius of the sink particle and dust sublimation maintains the consistency of the central intensity. Therefore, we chose to continue with the dust temperature from RAMSES.

Finally, we imaged all the disks present in the simulation snapshot by performing the ray-tracing through the grid to get their continuum emissions at ALMA's band 7 wavelength of $\lambda = 0.89$ mm. The imaging on the plane of the sky is done through 2D regularly-spaced cameras observing the 3D AMR grids and, thus, for practical computational reasons, we proceeded with the second-order ray-tracing by zooming on individual sinks particles whose disks are detected. We carried out the imaging procedure in sub-regions of the grid around the sink positions with an 1 au pixel resolution, which allowed us to emulate the maximum numerical resolution of ≈1.2 au in the most resolved regions of the disks. For each of these regions, $10^8$ scattering photon packages were used for the MC scattering calculations before the ray-tracing to get the intensity maps. We accounted for the variation of the viewing angle in observations by varying the camera positions to image the disks in three different planes of the internal simulation grid: xy, xz, and yz. We note that these projections do not necessarily correspond to the direction of the angular momentum vector of the disks, which are randomly oriented by nature.
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Images of the disk around sink 57 viewed in different projections and through different stages of the pipeline. From left to right: column density maps of the 1000 au region centred around the sink particle in 3 projections along the internal axes from the RAMSES simulation; continuum emission maps obtained with RADMC-3D of the same regions; and CASA images obtained with tclean from the C43-3, C43-4 and C43-6 visibilities. The disk sizes inferred from the gas kinematics, R<sub>sim</sub>, and the modelled disk radii, R<sub>90%</sub>, (detailed in Sect. 6) are noted for each configuration and projection.</figcaption></figure>

**4. Mock interferometric observations**

The RADMC-3D maps provide us with the ideal images of different regions of the sky surrounding each individual protostar. From these maps, by introducing proper instrumental effects, we can produce realistic observations of the simulated sky as would otherwise be obtained with interferometers. To that end, we use the ALMA simulator of the CASA software (McMullin et al. 2007; CASA Team et al. 2022) to 'observe' each of the stars and simulate a real observing session.

In particular, each emission map from RADMC-3D was used as the sky model for the simobserve task to get three simulated measurement sets for compact and extended ALMA antenna configurations, C43-3, C43-4, and C43-6, in Band 7 (at 0.89 mm, or 345 GHz); this corresponds to three angular resolutions of $\theta_{res} = 0.41, 0.266, 0.0887''$, respectively. For C43-3 and C43-4, which are relatively lower resolution, we integrated, in a "survey" fashion, for 300 s; this is comparable to the integration time of ALMA Band 7 observations in Evans et al. (2015) and Ohashi et al. (2022). For the high-resolution observations with C43-6, we increased the integration time as inversely proportional to $\theta_{res}^2$ to ~2700 s to achieve the same signal-to-noise ratio (S/N) per resolution element. A precipitable water vapor (PWV) value of 0.7 mm (as under good atmospheric conditions) was assumed. Thermal noise (default ATM model; Pardo 2019) was added to the simulated data. The measurement sets simulated by simobserve were then imaged and analyzed using the simanalyze task, which uses the tclean algorithm to perform an iterative cleaning process based on the 'clark' deconvolver (Clark 1980) and the 'briggs' weighting scheme (Briggs 1995), with a robust factor of 0.5.

Figure 2 illustrates (from left to right) the outcome at each step of the pipeline, namely: the column density integrated along one of the three coordinate axes x, y, or z through the entire grid, the dust continuum emission maps from RADMC-3D, and the visibility data achieved with C43-3, C43-4, and C43-6 imaged with the tclean algorithm integrated in CASA for one disk in the population. The structures are filtered out in the CASA images due to an incomplete coverage of the uv-plane; thus, recovering the disk properties at this stage is not trivial, just like recovering the 'real-life' disk properties in real observations.

Based on the images produced by tclean and their radial intensity profiles from the star positions, we selected the isolated disks for our subsequent analyses from the uv-data. Since the modelling of binaries and triple systems in the uv-space is not trivial without the use of a complex analytical model, such systems (which comprise of 10–11 out of the 26 disks formed in this simulation) were excluded in our current study. However, we note that we selected the single systems by looking at the observations imaged with tclean – and not at the images of the original disks; this was done to avoid prior knowledge of the real objects. As such, any imaged system in which there is more than one source to model within the primary beam is considered a binary or multiple, provided that the angular resolution of the observation is sufficient to produce a noticeable separation of the multiple sources.

**5. Synthetic modelling**

We now treat these measurement sets produced by CASA as if they were real observations and we applied similar techniques
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Images of the disk around sink 57 viewed in different projections and through different stages of the pipeline. From left to right: column density maps of the 1000 au region centred around the sink particle in 3 projections along the internal axes from the RAMSES simulation; continuum emission maps obtained with RADMC-3D of the same regions; and CASA images obtained with tclean from the C43-3, C43-4 and C43-6 visibilities. The disk sizes inferred from the gas kinematics, R<sub>sim</sub>, and the modelled disk radii, R<sub>90%</sub>, (detailed in Sect. 6) are noted for each configuration and projection.</figcaption></figure>

# 4. Mock interferometric observations

The RADMC-3D maps provide us with the ideal images of different regions of the sky surrounding each individual protostar. From these maps, by introducing proper instrumental effects, we can produce realistic observations of the simulated sky as would otherwise be obtained with interferometers. To that end, we use the ALMA simulator of the CASA software (McMullin et al. 2007; CASA Team et al. 2022) to 'observe' each of the stars and simulate a real observing session.

In particular, each emission map from RADMC-3D was used as the sky model for the simobserve task to get three simulated measurement sets for compact and extended ALMA antenna configurations, C43-3, C43-4, and C43-6, in Band 7 (at 0.89 mm, or 345 GHz); this corresponds to three angular resolutions of $\theta_{res} = 0.41, 0.266, 0.0887''$, respectively. For C43-3 and C43-4, which are relatively lower resolution, we integrated, in a "survey" fashion, for 300 s; this is comparable to the integration time of ALMA Band 7 observations in Evans et al. (2015) and Ohashi et al. (2022). For the high-resolution observations with C43-6, we increased the integration time as inversely proportional to $\theta_{res}^2$ to ~2700 s to achieve the same signal-to-noise ratio (S/N) per resolution element. A precipitable water vapor (PWV) value of 0.7 mm (as under good atmospheric conditions) was assumed. Thermal noise (default ATM model; Pardo 2019) was added to the simulated data. The measurement sets simulated by simobserve were then imaged and analyzed using the simanalyze task, which uses the tclean algorithm to perform an iterative cleaning process based on the 'clark' deconvolver (Clark 1980) and the 'briggs' weighting scheme (Briggs 1995), with a robust factor of 0.5.

Figure 2 illustrates (from left to right) the outcome at each step of the pipeline, namely: the column density integrated along one of the three coordinate axes x, y, or z through the entire grid, the dust continuum emission maps from RADMC-3D, and the visibility data achieved with C43-3, C43-4, and C43-6 imaged with the tclean algorithm integrated in CASA for one disk in the population. The structures are filtered out in the CASA images due to an incomplete coverage of the uv-plane; thus, recovering the disk properties at this stage is not trivial, just like recovering the 'real-life' disk properties in real observations.

Based on the images produced by tclean and their radial intensity profiles from the star positions, we selected the isolated disks for our subsequent analyses from the uv-data. Since the modelling of binaries and triple systems in the uv-space is not trivial without the use of a complex analytical model, such systems (which comprise of 10–11 out of the 26 disks formed in this simulation) were excluded in our current study. However, we note that we selected the single systems by looking at the observations imaged with tclean – and not at the images of the original disks; this was done to avoid prior knowledge of the real objects. As such, any imaged system in which there is more than one source to model within the primary beam is considered a binary or multiple, provided that the angular resolution of the observation is sufficient to produce a noticeable separation of the multiple sources.

# 5. Synthetic modelling

We now treat these measurement sets produced by CASA as if they were real observations and we applied similar techniques
<figure><img src="image_3.png" /><figcaption>Fig. 3. Schematic illustration of the disk and envelope model used to describe the emission of the YSOs formed in the simulation snapshot. A disk whose emission is described by an elliptic 2D Gaussian of semi-major axis $\sigma_{disk}$ is surrounded by a circular envelope following a Plummer-like density profile with rotation-to-infall transitional radius, $R_i$, and truncated outer radius, $R_{out}$.</figcaption></figure>

to model the “observed” disks to those in recent observational studies of YSOs (Maury et al. 2019; Tazzari et al. 2021a,b). This experiment allowed us to obtain their observable properties, namely: disk sizes, masses, and inclination angles.

## 5.1. Two-component model for the disk and the envelope

In our synthetic Class 0/I disk sample, due to the strong presence of the envelopes surrounding the young protostars, we opted to model each of these objects as a disk-envelope system, as illustrated in Fig. 3. Thus, in addition to the Gaussian profile for the disk, we adopted the empirical analytical model of a truncated Plummer envelope according to the method presented in Maury et al. (2019) to account for the more extended emission.

We let $I_0$ be the total intensity at $r = 0$, of which a fraction, $f$, is attributed to the disk. We can then express the surface brightness profile of the 2D disk of semi-major axis, $\sigma_{disk}$, as:

$$I_{v,disk}(\bar{r}) = f I_0 \exp\left(-\frac{\bar{r}^2}{2\sigma_{disk}^2}\right), \qquad (2)$$

and that of the circularly symmetric Plummer envelope truncated at outer radius, $R_{out}$, with a power-law distribution of the density

of $\rho(\bar{r}) \propto \bar{r}^{-m}$ and temperature of $T(\bar{r}) \propto \bar{r}^{-n}$ (Adams 1991) as:

$$I_{v,env}(\bar{r}) = \begin{cases} \frac{(1-f)I_0}{(1+(\bar{r}/R_i)^2)^{-(m+n-1)/2}} & \text{for } 0 \le \bar{r} \le R_{out} \\ 0 & \text{for } \bar{r} > R_{out} \end{cases}. \qquad (3)$$

Here, $\bar{r}$ is the projected radius in the sky and $R_i$ is the envelope’s radius at which the approximately uniform density of the inner region ends, which can be physically interpreted as the transition between envelope motions at $r > R_i$ and rotational motions within $r < R_i$.

## 5.2. Parameter fitting

Armed with these models, we could then use the galario (Tazzari et al. 2018) and emcee (Foreman-Mackey et al. 2013) libraries for Python to fit the visibilities obtained with CASA. As the system is comprised of an inclined disk and a spherical envelope, for each iteration, we used the radial profiles to construct two separate 2D images of the disk with an inclination angle, $i$, introduced, along with a flat circular envelope with the above-mentioned parameters; in addition an offset (dRA, dDec) was added from the centre and position angle (PA). Then, we superposed the two images and computed the synthetic visibility data for the model. The value of $\chi^2$ can be obtained by comparing these visibilities with the “observed” ones. Finally, we ran the affine invariant Markov chain Monte Carlo (MCMC) from Goodman & Weare (2010) with 80 walkers to explore the parameter space for each single disk present in the simulation output. A “burn-in” period of 4000 steps is initialised, followed by 20 000 main steps to guarantee convergence.

Figure 4 shows an example of the results obtained with this fitting method for one of the disks in the simulation snapshot. The fits for the visibility data returned by simobserve for three configurations C43-4, C43-4, and C43-6 with the disk-envelope model are provided in the upper panels, along with the triangular plots showing the correlation between the fitted parameters in the lower ones. The real parts of the uv data are aptly fitted with the models, whereas there is an offset at small uv distance in the imaginary parts, suggesting that the large scale components (i.e. the envelopes) are not perfectly modelled. Nevertheless, it might be more intuitive to check the goodness of the modelling by looking at the easier-to-visualise image plane. Thus, for illustrative and comparative purposes, an image of the model after interferometric filtering can be obtained by CASA’s tclean to compare with the synthetic image of the simulated disks by looking at the residual of the two images. An example of such images is shown in the lower panels of Fig. 5. The residual map is within the 3$\sigma$ noise level, showing relatively well modelled systems in the image plane. For the remainder of the disks in our simulation in the same xy projection, we refer to Appendix B for the model images and the comparison with their corresponding imaged observations from the three antenna configurations.

The results returned by galario have allowed us to obtain a simple analytical model of the disks. This kind of model has the advantage that it enables us to investigate the disk properties directly from the parameters in a self-consistent way. Moreover, it is also the method that is most widely used in disk modelling from observations (see, e.g. Maury et al. 2019; Sanchis et al. 2020); therefore, it is more comparable with the properties extracted from real observations than those directly inferred from simulations.

# 6. Analysis

The disk properties obtained from the modelling of observations often suffer from numerous uncertainties due to instrumental and
<figure><img src="image_5.png" /><figcaption>Fig. 3. Schematic illustration of the disk and envelope model used to describe the emission of the YSOs formed in the simulation snapshot. A disk whose emission is described by an elliptic 2D Gaussian of semi-major axis $\sigma_{disk}$ is surrounded by a circular envelope following a Plummer-like density profile with rotation-to-infall transitional radius, $R_i$, and truncated outer radius, $R_{out}$.</figcaption></figure>

to model the “observed” disks to those in recent observational studies of YSOs (Maury et al. 2019; Tazzari et al. 2021a,b). This experiment allowed us to obtain their observable properties, namely: disk sizes, masses, and inclination angles.

### 5.1. Two-component model for the disk and the envelope

In our synthetic Class 0/I disk sample, due to the strong presence of the envelopes surrounding the young protostars, we opted to model each of these objects as a disk-envelope system, as illustrated in Fig. 3. Thus, in addition to the Gaussian profile for the disk, we adopted the empirical analytical model of a truncated Plummer envelope according to the method presented in Maury et al. (2019) to account for the more extended emission.

We let $I_0$ be the total intensity at $r = 0$, of which a fraction, $f$, is attributed to the disk. We can then express the surface brightness profile of the 2D disk of semi-major axis, $\sigma_{disk}$, as:

$$I_{v,disk}(\bar{r}) = f I_0 \exp\left(-\frac{\bar{r}^2}{2\sigma_{disk}^2}\right), \qquad (2)$$

and that of the circularly symmetric Plummer envelope truncated at outer radius, $R_{out}$, with a power-law distribution of the density

of $\rho(\bar{r}) \propto \bar{r}^{-m}$ and temperature of $T(\bar{r}) \propto \bar{r}^{-n}$ (Adams 1991) as:

$$I_{v,env}(\bar{r}) = \begin{cases} \frac{(1-f)I_0}{(1+(\bar{r}/R_i)^2)^{-(m+n-1)/2}} & \text{for } 0 \le \bar{r} \le R_{out} \\ 0 & \text{for } \bar{r} > R_{out} \end{cases}. \qquad (3)$$

Here, $\bar{r}$ is the projected radius in the sky and $R_i$ is the envelope’s radius at which the approximately uniform density of the inner region ends, which can be physically interpreted as the transition between envelope motions at $r > R_i$ and rotational motions within $r < R_i$.

### 5.2. Parameter fitting

Armed with these models, we could then use the galario (Tazzari et al. 2018) and emcee (Foreman-Mackey et al. 2013) libraries for Python to fit the visibilities obtained with CASA. As the system is comprised of an inclined disk and a spherical envelope, for each iteration, we used the radial profiles to construct two separate 2D images of the disk with an inclination angle, $i$, introduced, along with a flat circular envelope with the above-mentioned parameters; in addition an offset (dRA, dDec) was added from the centre and position angle (PA). Then, we superposed the two images and computed the synthetic visibility data for the model. The value of $\chi^2$ can be obtained by comparing these visibilities with the “observed” ones. Finally, we ran the affine invariant Markov chain Monte Carlo (MCMC) from Goodman & Weare (2010) with 80 walkers to explore the parameter space for each single disk present in the simulation output. A “burn-in” period of 4000 steps is initialised, followed by 20 000 main steps to guarantee convergence.

Figure 4 shows an example of the results obtained with this fitting method for one of the disks in the simulation snapshot. The fits for the visibility data returned by simobserve for three configurations C43-4, C43-4, and C43-6 with the disk-envelope model are provided in the upper panels, along with the triangular plots showing the correlation between the fitted parameters in the lower ones. The real parts of the uv data are aptly fitted with the models, whereas there is an offset at small uv distance in the imaginary parts, suggesting that the large scale components (i.e. the envelopes) are not perfectly modelled. Nevertheless, it might be more intuitive to check the goodness of the modelling by looking at the easier-to-visualise image plane. Thus, for illustrative and comparative purposes, an image of the model after interferometric filtering can be obtained by CASA’s tclean to compare with the synthetic image of the simulated disks by looking at the residual of the two images. An example of such images is shown in the lower panels of Fig. 5. The residual map is within the 3$\sigma$ noise level, showing relatively well modelled systems in the image plane. For the remainder of the disks in our simulation in the same xy projection, we refer to Appendix B for the model images and the comparison with their corresponding imaged observations from the three antenna configurations.

The results returned by galario have allowed us to obtain a simple analytical model of the disks. This kind of model has the advantage that it enables us to investigate the disk properties directly from the parameters in a self-consistent way. Moreover, it is also the method that is most widely used in disk modelling from observations (see, e.g. Maury et al. 2019; Sanchis et al. 2020); therefore, it is more comparable with the properties extracted from real observations than those directly inferred from simulations.

## 6. Analysis

The disk properties obtained from the modelling of observations often suffer from numerous uncertainties due to instrumental and
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Fitting results from galario for the synthetic observations of the disk around sink 57, seen in the *xy* projection. Left: observed (black dots with error bars) and modelled (red curve) visibilities in real (top) and imaginary parts (lower panel) as a function of the baseline (in kλ) for three configurations C43-3 (upper), C43-4 (middle), and C43-6 (lower). Right: corner plots of the galario fittings of the visibility data with the disk and envelope model. The top sub-panels show the 1D histograms of the free parameters from the MCMC chains. The rest of the panels are the 2D histograms between each pair of parameters.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Fitting results from galario for the synthetic observations of the disk around sink 57, seen in the *xy* projection. Left: observed (black dots with error bars) and modelled (red curve) visibilities in real (top) and imaginary parts (lower panel) as a function of the baseline (in kλ) for three configurations C43-3 (upper), C43-4 (middle), and C43-6 (lower). Right: corner plots of the galario fittings of the visibility data with the disk and envelope model. The top sub-panels show the 1D histograms of the free parameters from the MCMC chains. The rest of the panels are the 2D histograms between each pair of parameters.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Modelling results imaged for the disk around sink 57. Column density plot from the RAMSES output with the circle enclosing the disk size (top panel) derived from simulation (left), and dust continuum image obtained with RADMC-3D (right) of the 1000 au region around the star. The CASA emission map imaged (bottom panel) with tclean from the simobserve visibilities with C43-4 (left), emission map of the best-fit model obtained with tclean (middle), and the corresponding residual (right). Contour lines are plotted at 3, 6, 12, 24, 48, and 96 multiples of the rms, σ. The disk properties inferred from the simulation and their modelled properties from observations are noted in the corresponding panels.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Modelling results imaged for the disk around sink 57. Column density plot from the RAMSES output with the circle enclosing the disk size (top panel) derived from simulation (left), and dust continuum image obtained with RADMC-3D (right) of the 1000 au region around the star. The CASA emission map imaged (bottom panel) with tclean from the simobserve visibilities with C43-4 (left), emission map of the best-fit model obtained with tclean (middle), and the corresponding residual (right). Contour lines are plotted at 3, 6, 12, 24, 48, and 96 multiples of the rms, σ. The disk properties inferred from the simulation and their modelled properties from observations are noted in the corresponding panels.</figcaption></figure>
interferometric effects. The information regarding their intrinsic properties that we are able to derive from our simulations of their formation allows us to quantify the accuracy of these observational measurements. Thus, in this section, we take the simulation-derived properties as the “ground truth” that forms the basis of the subsequent assessment.

## 6.1. Disk sizes

### 6.1.1. Radius measurements

The most accessible and best studied physical property of the disks from continuum observations is probably the disk size, whose measurements are essential to study the disk evolution driven either by the viscosity or MHD disk wind (Toci et al. 2021, 2023; Rosotti et al. 2019; Miotello et al. 2023). We therefore use our synthetic observations to assess the relationship between the radii derived by the simple analytical modelling of the observations with the theoretically defined disk radii as measured in numerical simulations.

On the theoretical side, accompanied with the simulations of the PPD formation, Lebreuilly et al. (2021) proposed a method known as “disk finder” to select the disk materials from the gas kinematics based on the criteria from Joos et al. (2012), which can be highlighted as follows. In cylindrical coordinates (R; φ; z), co-moving with the central star, the disk materials will be selected if they satisfied the following condition:

- Rotating faster than falling radially: $v_φ > 2v_r$,

- Rotating faster than falling vertically: $v_φ > 2v_z$,

- Composed of dense material of gas number density $n > n_{thre} = 10^9 \text{ cm}^{-3}$.

It is worth pointing out that while this result is closer to the true properties, this is still more of a statistical method for determining the parameters. As such, it does not aim to represent the physical reality of the object as disks in these simulations are often not well defined and separated from the collapsing flows. This is in line with what was found by Aso & Machida (2020) in a more extended investigation of the disk identification methods from both simulations and synthetic observations of an isolated protostellar disk formed in the single core-collapse scenario. Observationally, it is also the realistic picture confirmed by recent studies of Class 0/I systems in which the presence of streamers (see, e.g. Pineda et al. 2023; Bianchi et al. 2023; Cacciapuoti et al. 2023) or bridges (Sadavoy et al. 2018) were detected. With that in mind, we re-applied the same method to derive the disk radii, $R_{sim}$, from the simulations.

From the observations, the MCMC visibility probability fitting gives us two important parameters in the models: the semi-major axis of the Gaussian disk, $\sigma_{disk}$, and the radius of the inner region of the Plummer envelope, $R_i$. While $R_i$ could be a good estimation of the disk size as inferred from the envelope’s geometry, in parallel we used a more robust definition of the disk radius from the Gaussian component, taking its emission into account. This employs a characterisation of the disk size by the radius within which a certain fraction of the total emission is contained. The most commonly used fractions in the literature are 68% (Tripathi et al. 2017; Andrews et al. 2018; Facchini et al. 2019; Long et al. 2019; Manara et al. 2019), 90%, or 95% (Tazzari et al. 2017). The corresponding radii can be related to $\sigma_{disk}$ as:

$$R_{68\%} = \sigma_{\text{disk}} \times \sqrt{-2 \times \ln(1 - 0.68)} \approx 1.51 \times \sigma_{\text{disk}}, \quad (4)$$

$$R_{90\%} \approx 1.42 \times R_{68\%},$$

$$R_{95\%} \approx 1.62 \times R_{68\%},$$

In this two component model, since $R_i$ and $\sigma_{disk}$ are not explicitly related according to our analytical description, there could be four scenarios for the results of the fits regarding the disk sizes:

1. The Gaussian disk is within the rotation-to-infall transition marked by $R_i$ and both are above the resolution threshold of the observation, namely, $\theta_{res}/(2\sqrt{2\ln2}) < \sigma_{disk} < R_i$.

2. The Gaussian disk size is below the resolution of the observation and the rotation-to-infall transition is above, namely, $\sigma_{disk} < \theta_{res}/(2\sqrt{2\ln2}) < R_i$. This usually happens for large disks in the more resolved observations (mostly C43-6), whose low-surface-brightness structure is missed by the Gaussian, which instead ends up fitting the narrow central peak.

3. The Gaussian disk radius is beyond the rotation-to-infall transition and its peak intensity is lower than that of the envelope, namely, $\sigma_{disk} > R_i$ and $f < 0.5$. We consider the Gaussian to be tracing the envelope’s more extended emission in such case.

4. The Gaussian disk radius and the rotation-to-infall transition are both below the resolution threshold of the observation, namely, $\sigma_{disk}, R_i < \theta_{res}/(2\sqrt{2\ln2})$.

In the first scenario, we used the three Gaussian radii given by Eq. (4) to determine the observational disk radius, $R_{obs}$, and subsequently, the disk mass $M_{obs}$ (in combination with the peak intensity of the disk, $fI_0$). In the next two cases, we switched the roles of the two analytical components, taking the inner structure of the Plummer-like component as the disk and the Gaussian the envelope, namely: using $R_i$ and $(1-f)I_0$ as the radius and the peak intensity of the disk instead, so as to minimise the number of outliers. In the final case, we used the greater radius of the two and the corresponding properties. Our pre-testing found that this alternative approach to using the Gaussian’s radii mainly helps to improve the results for C43-6 high angular resolution observations. This is because it works to compensate for the weakness of the Gaussian in tracing the low-surface-brightness structure of the larger disks; whereas for the other two lower-resolution configurations, it has a negligible impact on the results.

We started with the C43-4 configuration, which allowed us to sufficiently resolve most of the disks in the sample with its angular resolution of 0.266" (≈42 au), comparable to that of recent surveys, such as Tobin et al. (2020) and Sheehan et al. (2022). Figure 6 compares the value we get at the two ends. We plot in the upper panels the disk sizes obtained from the synthetic modelling of the uv-data ‘observed’ with C43-4 and Eq. (4) versus their sizes from the disk finder method for the RAMSES simulations. The uncertainties (marked by the corresponding error bars) are determined by a lower and upper boundary at the 16th and 84th percentiles of the samples, respectively. The $R_{90\%}$ radii follow relatively well the one-to-one correlation shown by the black straight curve. So does $R_{68\%}$, although it does fall slightly below the one-to-one line for the most part. On the other hand, $R_{95\%}$ tends to overestimate disk size most of the time.

As a more direct comparison, we divided the three radii by the disk radius values from the simulation and plotted the histograms of the logarithms (base 10) of these ratios in the three projections: xy, yz, xz according to the simulation grid (shown in the lower panels). Since all these ratios can be related to $\sigma_{disk}$ and to each other by a certain factor, the standard deviation $\sigma_{ratio}$ calculated for $R_{68\%}$, $R_{90\%}$, and $R_{95\%}$ does not differ. Thus, the distributions of the three ratios in logarithmic scale can be related
interferometric effects. The information regarding their intrinsic properties that we are able to derive from our simulations of their formation allows us to quantify the accuracy of these observational measurements. Thus, in this section, we take the simulation-derived properties as the “ground truth” that forms the basis of the subsequent assessment.

## 6.1. Disk sizes

### 6.1.1. Radius measurements

The most accessible and best studied physical property of the disks from continuum observations is probably the disk size, whose measurements are essential to study the disk evolution driven either by the viscosity or MHD disk wind (Toci et al. 2021, 2023; Rosotti et al. 2019; Miotello et al. 2023). We therefore use our synthetic observations to assess the relationship between the radii derived by the simple analytical modelling of the observations with the theoretically defined disk radii as measured in numerical simulations.

On the theoretical side, accompanied with the simulations of the PPD formation, Lebreuilly et al. (2021) proposed a method known as “disk finder” to select the disk materials from the gas kinematics based on the criteria from Joos et al. (2012), which can be highlighted as follows. In cylindrical coordinates (R; φ; z), co-moving with the central star, the disk materials will be selected if they satisfied the following condition:

- Rotating faster than falling radially: $v_φ > 2v_r$,

- Rotating faster than falling vertically: $v_φ > 2v_z$,

- Composed of dense material of gas number density $n > n_{thre} = 10^9 \text{ cm}^{-3}$.

It is worth pointing out that while this result is closer to the true properties, this is still more of a statistical method for determining the parameters. As such, it does not aim to represent the physical reality of the object as disks in these simulations are often not well defined and separated from the collapsing flows. This is in line with what was found by Aso & Machida (2020) in a more extended investigation of the disk identification methods from both simulations and synthetic observations of an isolated protostellar disk formed in the single core-collapse scenario. Observationally, it is also the realistic picture confirmed by recent studies of Class 0/I systems in which the presence of streamers (see, e.g. Pineda et al. 2023; Bianchi et al. 2023; Cacciapuoti et al. 2023) or bridges (Sadavoy et al. 2018) were detected. With that in mind, we re-applied the same method to derive the disk radii, $R_{sim}$, from the simulations.

From the observations, the MCMC visibility probability fitting gives us two important parameters in the models: the semi-major axis of the Gaussian disk, $\sigma_{disk}$, and the radius of the inner region of the Plummer envelope, $R_i$. While $R_i$ could be a good estimation of the disk size as inferred from the envelope’s geometry, in parallel we used a more robust definition of the disk radius from the Gaussian component, taking its emission into account. This employs a characterisation of the disk size by the radius within which a certain fraction of the total emission is contained. The most commonly used fractions in the literature are 68% (Tripathi et al. 2017; Andrews et al. 2018; Facchini et al. 2019; Long et al. 2019; Manara et al. 2019), 90%, or 95% (Tazzari et al. 2017). The corresponding radii can be related to $\sigma_{disk}$ as:

$$R_{68\%} = \sigma_{\text{disk}} \times \sqrt{-2 \times \ln(1 - 0.68)} \approx 1.51 \times \sigma_{\text{disk}}, \quad (4)$$

$$R_{90\%} \approx 1.42 \times R_{68\%},$$

$$R_{95\%} \approx 1.62 \times R_{68\%},$$

In this two component model, since $R_i$ and $\sigma_{disk}$ are not explicitly related according to our analytical description, there could be four scenarios for the results of the fits regarding the disk sizes:

1. The Gaussian disk is within the rotation-to-infall transition marked by $R_i$ and both are above the resolution threshold of the observation, namely, $\theta_{res}/(2\sqrt{2\ln2}) < \sigma_{disk} < R_i$.

2. The Gaussian disk size is below the resolution of the observation and the rotation-to-infall transition is above, namely, $\sigma_{disk} < \theta_{res}/(2\sqrt{2\ln2}) < R_i$. This usually happens for large disks in the more resolved observations (mostly C43-6), whose low-surface-brightness structure is missed by the Gaussian, which instead ends up fitting the narrow central peak.

3. The Gaussian disk radius is beyond the rotation-to-infall transition and its peak intensity is lower than that of the envelope, namely, $\sigma_{disk} > R_i$ and $f < 0.5$. We consider the Gaussian to be tracing the envelope’s more extended emission in such case.

4. The Gaussian disk radius and the rotation-to-infall transition are both below the resolution threshold of the observation, namely, $\sigma_{disk}, R_i < \theta_{res}/(2\sqrt{2\ln2})$.

In the first scenario, we used the three Gaussian radii given by Eq. (4) to determine the observational disk radius, $R_{obs}$, and subsequently, the disk mass $M_{obs}$ (in combination with the peak intensity of the disk, $fI_0$). In the next two cases, we switched the roles of the two analytical components, taking the inner structure of the Plummer-like component as the disk and the Gaussian the envelope, namely: using $R_i$ and $(1-f)I_0$ as the radius and the peak intensity of the disk instead, so as to minimise the number of outliers. In the final case, we used the greater radius of the two and the corresponding properties. Our pre-testing found that this alternative approach to using the Gaussian’s radii mainly helps to improve the results for C43-6 high angular resolution observations. This is because it works to compensate for the weakness of the Gaussian in tracing the low-surface-brightness structure of the larger disks; whereas for the other two lower-resolution configurations, it has a negligible impact on the results.

We started with the C43-4 configuration, which allowed us to sufficiently resolve most of the disks in the sample with its angular resolution of 0.266" (≈42 au), comparable to that of recent surveys, such as Tobin et al. (2020) and Sheehan et al. (2022). Figure 6 compares the value we get at the two ends. We plot in the upper panels the disk sizes obtained from the synthetic modelling of the uv-data ‘observed’ with C43-4 and Eq. (4) versus their sizes from the disk finder method for the RAMSES simulations. The uncertainties (marked by the corresponding error bars) are determined by a lower and upper boundary at the 16th and 84th percentiles of the samples, respectively. The $R_{90\%}$ radii follow relatively well the one-to-one correlation shown by the black straight curve. So does $R_{68\%}$, although it does fall slightly below the one-to-one line for the most part. On the other hand, $R_{95\%}$ tends to overestimate disk size most of the time.

As a more direct comparison, we divided the three radii by the disk radius values from the simulation and plotted the histograms of the logarithms (base 10) of these ratios in the three projections: xy, yz, xz according to the simulation grid (shown in the lower panels). Since all these ratios can be related to $\sigma_{disk}$ and to each other by a certain factor, the standard deviation $\sigma_{ratio}$ calculated for $R_{68\%}$, $R_{90\%}$, and $R_{95\%}$ does not differ. Thus, the distributions of the three ratios in logarithmic scale can be related
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Comparison of the disk sizes measured from the observations with the C43-4 configuration and the theoretical radii. Disk sizes in arcsecs obtained from the modelling of CASA observations vs. the values derived from the gas kinematics in the simulation for all single star-disk systems present in the output, seen in three projections: xy (left), yz (middle), and xz (right), shown at the top. Black dashed line shows the one-to-one correlation. Histograms of the R<sub>90%</sub>/R<sub>sim</sub> ratio for the three projections, shown at the bottom. Black dotted line marks R<sub>90%</sub>/R<sub>sim</sub> = 1, where the two are in full agreement. Vertical dashed lines show the medians of the ratio distributions of three radii R<sub>68%</sub> (red), R<sub>90%</sub> (green), and R<sub>95%</sub> (blue). The values within the 16th to the 84th percentiles are covered by the shaded areas. The statistics (including the mean, median, and standard deviation) for the three radii are indicated above each panel.</figcaption></figure>

by a horizontal shift equal to the shifts of the median values indi-
cated by the vertical dashed lines; this is also the case for the
size of the 1σ ranges indicated by the shaded areas under the
histograms. The values of σ<sub>ratio</sub> ~ 0.22–0.23 indicate a factor of
~1.6–1.7 in the uncertainty of the estimates. For the *xy* and *yz*
projections, the distributions for *R*<sub>68%</sub>/*R*<sub>sim</sub> and *R*<sub>90%</sub>/*R*<sub>sim</sub> peaks
near unity, indicated by the vertical dotted black lines, as well
as the means and the medians of their distributions. Based on
their statistics, *R*<sub>68%</sub> and *R*<sub>90%</sub> distributions have medians equally
close to 1 in general, with *R*<sub>68%</sub> slightly underestimates and *R*<sub>90%</sub>
slightly overestimates the disk radii. To include as much as pos-
sible the disk emission, we use *R*<sub>90%</sub> for the histograms of the
ratio of the disk radii, as the values for *R*<sub>68%</sub> and *R*<sub>95%</sub> can be
related to it by a certain factor.

6.1.2. Effect of resolution

Here, we focus on how different resolutions of the observa-
tions could affect the disk sizes obtained with our visibility
modelling. For that purpose, we used the same model and fit-
ting method on the data obtained with C43-3 and C43-6, with
angular resolutions of θ<sub>res</sub> = 0.41'', 0.0887'', respectively. Fig-
ure 7 shows the histogram of the ratio between the modelled
radii over the simulations' radii for the data from C43-3 (left
panel) and C43-6 (right panel), along with the results for C43-
4 with intermediate resolution that we analyzed in the previous

sub-section. We merged the projected disks in three projec-
tions into one single sample for each configuration to maximise
the sample size. From C43-3 to C43-4, as the angular reso-
lution of the observation increases, the mean and median of
the log(R<sub>90%</sub>/R<sub>sim</sub>) distribution get closer to 0, which represents
unity in linear scale. The distribution for C43-4 also appears
more centred and peaked, and the factor of uncertainty decreases
from ~2.3 to ~1.7. From C43-4 to C43-6, the distribution keeps
getting slightly better mean (factor of 1.1 difference compared to
1.3), but the uncertainty now slightly rises up to ~2.2 as more
outliers appear.

As suggested by the dotted red lines in the upper panels of Fig. 7 (displaying the angular resolution of the corresponding observations), the resolution achievable by C43-4, namely, θ<sub>res</sub> = 0.266'', allows us to sufficiently observe all the disks in this sample, whereas with C43-3, a significant portion of it falls under the resolution limit of the observation. This explains the right-skewness and the offsets of the mean and median of the R<sub>90%</sub>/R<sub>sim</sub> ratio distribution for C43-3 and the smaller uncertainty in the case of C43-4. Furthermore, the fact that C43-4's angular resolution is very close to the smallest disk sizes in the sample means that a number of disks are only partially resolved. This is also the situation where the Gaussian profile describes best the disk emission due to the effect of the beam. C43-6, on the other hand, provides an angular resolution that is more than enough to resolve even the smallest disks in the sample and therefore suffers from larger uncertainty. This is because the well
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Comparison of the disk sizes measured from the observations with the C43-4 configuration and the theoretical radii. Disk sizes in arcsecs obtained from the modelling of CASA observations vs. the values derived from the gas kinematics in the simulation for all single star-disk systems present in the output, seen in three projections: xy (left), yz (middle), and xz (right), shown at the top. Black dashed line shows the one-to-one correlation. Histograms of the R<sub>90%</sub>/R<sub>sim</sub> ratio for the three projections, shown at the bottom. Black dotted line marks R<sub>90%</sub>/R<sub>sim</sub> = 1, where the two are in full agreement. Vertical dashed lines show the medians of the ratio distributions of three radii R<sub>68%</sub> (red), R<sub>90%</sub> (green), and R<sub>95%</sub> (blue). The values within the 16th to the 84th percentiles are covered by the shaded areas. The statistics (including the mean, median, and standard deviation) for the three radii are indicated above each panel.</figcaption></figure>

by a horizontal shift equal to the shifts of the median values indi-
cated by the vertical dashed lines; this is also the case for the
size of the 1σ ranges indicated by the shaded areas under the
histograms. The values of σ<sub>ratio</sub> ~ 0.22–0.23 indicate a factor of
~1.6–1.7 in the uncertainty of the estimates. For the *xy* and *yz*
projections, the distributions for *R*<sub>68%</sub>/*R*<sub>sim</sub> and *R*<sub>90%</sub>/*R*<sub>sim</sub> peaks
near unity, indicated by the vertical dotted black lines, as well
as the means and the medians of their distributions. Based on
their statistics, *R*<sub>68%</sub> and *R*<sub>90%</sub> distributions have medians equally
close to 1 in general, with *R*<sub>68%</sub> slightly underestimates and *R*<sub>90%</sub>
slightly overestimates the disk radii. To include as much as pos-
sible the disk emission, we use *R*<sub>90%</sub> for the histograms of the
ratio of the disk radii, as the values for *R*<sub>68%</sub> and *R*<sub>95%</sub> can be
related to it by a certain factor.

6.1.2. Effect of resolution

Here, we focus on how different resolutions of the observa-
tions could affect the disk sizes obtained with our visibility
modelling. For that purpose, we used the same model and fit-
ting method on the data obtained with C43-3 and C43-6, with
angular resolutions of θ<sub>res</sub> = 0.41'', 0.0887'', respectively. Fig-
ure 7 shows the histogram of the ratio between the modelled
radii over the simulations' radii for the data from C43-3 (left
panel) and C43-6 (right panel), along with the results for C43-
4 with intermediate resolution that we analyzed in the previous

sub-section. We merged the projected disks in three projec-
tions into one single sample for each configuration to maximise
the sample size. From C43-3 to C43-4, as the angular reso-
lution of the observation increases, the mean and median of
the log(R<sub>90%</sub>/R<sub>sim</sub>) distribution get closer to 0, which represents
unity in linear scale. The distribution for C43-4 also appears
more centred and peaked, and the factor of uncertainty decreases
from ~2.3 to ~1.7. From C43-4 to C43-6, the distribution keeps
getting slightly better mean (factor of 1.1 difference compared to
1.3), but the uncertainty now slightly rises up to ~2.2 as more
outliers appear.

As suggested by the dotted red lines in the upper panels of Fig. 7 (displaying the angular resolution of the corresponding observations), the resolution achievable by C43-4, namely, θ<sub>res</sub> = 0.266'', allows us to sufficiently observe all the disks in this sample, whereas with C43-3, a significant portion of it falls under the resolution limit of the observation. This explains the right-skewness and the offsets of the mean and median of the R<sub>90%</sub>/R<sub>sim</sub> ratio distribution for C43-3 and the smaller uncertainty in the case of C43-4. Furthermore, the fact that C43-4's angular resolution is very close to the smallest disk sizes in the sample means that a number of disks are only partially resolved. This is also the situation where the Gaussian profile describes best the disk emission due to the effect of the beam. C43-6, on the other hand, provides an angular resolution that is more than enough to resolve even the smallest disks in the sample and therefore suffers from larger uncertainty. This is because the well
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Comparison of the disk sizes measured from the observations with three different angular resolutions (0.41", 0.266", 0.887"). R<sub>90%</sub> radius of the observed disks versus their theoretical sizes for three configuration C43-3 (left), C43-4 (middle), and C43-6 (right), shown at the top. Histograms of the R<sub>90%</sub>/R<sub>sim</sub> ratio for the reported configurations, shown at the bottom.</figcaption></figure>

resolved disks would likely require a more complex description
than a Gaussian profile. This is visible in the scatter plot for
C43-6, where the few smaller disks that are not resolved (and
thus overestimated) by C43-4 are now slightly better modelled
with higher resolution observation; at the same time, a number
of larger disks are severely underestimated as the Gaussian is
trying to fit the more compact central peak rather than the more
extended disk structure that is filtered out by the observation.

6.1.3. Combined antenna arrays

With these observations of different resolutions in hand, it is pos-
sible to alternatively combine the data during the fitting to have
a better coverage of the uv space and better separate the com-
pact and extended emissions from the disk and the envelope,
respectively, as done by Cacciapuoti et al. (2023) to study the
dust properties within the envelope. Here, we adopt the typical
ALMA combination of the antenna arrays C43-3 and C43-6 for
our modelling. The χ² now is the sum of the values with the two
data sets.

A comparison of the results obtained with this strategy and
those with the separate visibilities from the two configurations
in the previous sub-section can be found in Fig. 8. Since the
disk radius, R<sub>90%</sub>, and the radius of the approximately uniform
density region of the envelope, R<sub>i</sub>, are two free parameters that
are completely unrelated to one another in the fittings (and both
more or less trace the outer disk radius) we opted to compare
them individually with the ‘true’ disk size. This was done to
evaluate the accuracy of the model in recovering and separating
the two components, contrary to what is detailed in Sect. 6.1.1.
Therefore, we show (in the upper panels) the values of R<sub>90%</sub>

and Rᵢ radii versus the disk size from simulations and (in the
lower panels) the histograms of the ratios of the two radii
modelled from observations over the simulation sizes, without
the fine-tuning criteria used in the two previous sub-sections. It
can be seen from the scatter plots that there are cases where the
values of R<sub>90%</sub> and Rᵢ returned by the fits are much smaller than
the angular resolution of the observations (marked by the dotted
red lines). We consider those, namely, σ<sub>disk</sub> < θ<sub>res</sub>/(2√2 ln 2) or
Rᵢ < θ<sub>res</sub>/(2√2 ln 2), to be cases where the fitting failed to detect
the disk (in the case of σ<sub>disk</sub>) or the transition between rotation
and infall motions (in the case of Rᵢ). Thus, we did not include
them in the samples for the lower panels' statistics.

The histograms suggest that the combined data genuinely
improve the measurements of the disk radii by eliminating both
the underestimated points returned by the C43-6 fittings and the
heavily overestimated ones by the C43-3 fittings due to the lack
of resolution. It is worth keeping in mind that we obtain three
non-detected Gaussian disks with the combined data; while that
is not the case with C43-6, these can still be substituted with Rᵢ
for the disk size (as described in the sub-sections above) since
the values of Rᵢ for these disks are above the observation's reso-
lution, which does not decrease nor increase the overall accuracy
of the disk radius extraction. This indicates that the extension
of the uv coverage could play a significant role in subtracting
the envelope's emission from that of the disk in the observa-
tions. Consequently, this could help the Gaussian in detecting
the correct disk structure, even when the statistics for the enve-
lope's Rᵢ does not gain any enhancement. As a matter of fact,
for C43-3, most of the results for Rᵢ fall below the resolution
of the configuration. Thus, the distributions for C43-3 consist
of very small samples which do not give meaningful statistical
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Comparison of the disk sizes measured from the observations with three different angular resolutions (0.41", 0.266", 0.887"). R<sub>90%</sub> radius of the observed disks versus their theoretical sizes for three configuration C43-3 (left), C43-4 (middle), and C43-6 (right), shown at the top. Histograms of the R<sub>90%</sub>/R<sub>sim</sub> ratio for the reported configurations, shown at the bottom.</figcaption></figure>

resolved disks would likely require a more complex description
than a Gaussian profile. This is visible in the scatter plot for
C43-6, where the few smaller disks that are not resolved (and
thus overestimated) by C43-4 are now slightly better modelled
with higher resolution observation; at the same time, a number
of larger disks are severely underestimated as the Gaussian is
trying to fit the more compact central peak rather than the more
extended disk structure that is filtered out by the observation.

6.1.3. Combined antenna arrays

With these observations of different resolutions in hand, it is pos-
sible to alternatively combine the data during the fitting to have
a better coverage of the uv space and better separate the com-
pact and extended emissions from the disk and the envelope,
respectively, as done by Cacciapuoti et al. (2023) to study the
dust properties within the envelope. Here, we adopt the typical
ALMA combination of the antenna arrays C43-3 and C43-6 for
our modelling. The χ² now is the sum of the values with the two
data sets.

A comparison of the results obtained with this strategy and
those with the separate visibilities from the two configurations
in the previous sub-section can be found in Fig. 8. Since the
disk radius, R<sub>90%</sub>, and the radius of the approximately uniform
density region of the envelope, R<sub>i</sub>, are two free parameters that
are completely unrelated to one another in the fittings (and both
more or less trace the outer disk radius) we opted to compare
them individually with the ‘true’ disk size. This was done to
evaluate the accuracy of the model in recovering and separating
the two components, contrary to what is detailed in Sect. 6.1.1.
Therefore, we show (in the upper panels) the values of R<sub>90%</sub>

and Rᵢ radii versus the disk size from simulations and (in the
lower panels) the histograms of the ratios of the two radii
modelled from observations over the simulation sizes, without
the fine-tuning criteria used in the two previous sub-sections. It
can be seen from the scatter plots that there are cases where the
values of R<sub>90%</sub> and Rᵢ returned by the fits are much smaller than
the angular resolution of the observations (marked by the dotted
red lines). We consider those, namely, σ<sub>disk</sub> < θ<sub>res</sub>/(2√2 ln 2) or
Rᵢ < θ<sub>res</sub>/(2√2 ln 2), to be cases where the fitting failed to detect
the disk (in the case of σ<sub>disk</sub>) or the transition between rotation
and infall motions (in the case of Rᵢ). Thus, we did not include
them in the samples for the lower panels' statistics.

The histograms suggest that the combined data genuinely
improve the measurements of the disk radii by eliminating both
the underestimated points returned by the C43-6 fittings and the
heavily overestimated ones by the C43-3 fittings due to the lack
of resolution. It is worth keeping in mind that we obtain three
non-detected Gaussian disks with the combined data; while that
is not the case with C43-6, these can still be substituted with Rᵢ
for the disk size (as described in the sub-sections above) since
the values of Rᵢ for these disks are above the observation's reso-
lution, which does not decrease nor increase the overall accuracy
of the disk radius extraction. This indicates that the extension
of the uv coverage could play a significant role in subtracting
the envelope's emission from that of the disk in the observa-
tions. Consequently, this could help the Gaussian in detecting
the correct disk structure, even when the statistics for the enve-
lope's Rᵢ does not gain any enhancement. As a matter of fact,
for C43-3, most of the results for Rᵢ fall below the resolution
of the configuration. Thus, the distributions for C43-3 consist
of very small samples which do not give meaningful statistical
<figure><img src="image_8.png" /><figcaption>Fig. 8. Comparison of different fitting strategies: with data from individual configurations (left and middle panels for C43-3 and C43-6, respectively) and from combined antenna arrays (right panels for the combined C43-3 and C43-6).</figcaption></figure>

indications. C43-6 appears to have the best measurements for $R_i$, although they slightly overestimate the disk radius in general. The combined $R_i$ distributions have a broader spread within the 16th-84th percentile range, which displays quite a flat shape, without any real centralised peak. This may be an indication that $R_i$ might not be a meaningful parameter to study the envelope's structure with the Plummer density profile, but is nonetheless useful when employed in combination with the Gaussian radii, to provide a more accurate measurement of the disk sizes when the Gaussian encounters difficulties locating the right disk structure, as can be seen in the results for C43-6 (lower-right panel of Fig. 7 and lower-middle panel of Fig. 8).

## 6.2. Disk masses

Another important property of the disks is their masses, as these values offer insights into the mass budget for the formation of planets in their later evolutionary stages.

From the central peak intensity that we modelled with the Gaussian profile in Eq. (3), we performed an integration to get the millimeter flux, $F_ν$, of the disks. These fluxes are shown in Fig. 9 as a function of the modelled disk radii, $R_{obs} = R_{90\%}$. As a reference, we plotted the fluxes derived analytically from $R_{obs}$ while assuming that the disks are fully optically thick:

$$F_{v,\text{th}} = \frac{1}{d^2} \int_{0}^{R_{\text{obs}}} B_{v}(T_{\text{dust}}) 2\pi r dr. \qquad (5)$$

Here, $d \sim 140$ pc is the distance to the source and $B_v$ is the Planck function. The dust temperature radial profile $T_{\text{dust}}(r)$ for each disk is inferred from the 3D profile in the simulation by fitting a power-law function of $T(r) = T(R_0)(r/R_0)^{-\nu}$ to the radially averaged temperature shown in Fig. A.1. Using the 16th,

50th, and 84th percentiles of the fitted parameters, we found $\bar{T}(R_0) = 181_{-35}^{+74}$ K, $\bar{R}_0 = 35_{-12}^{+36}$ au and $\bar{v} = 0.52_{-0.05}^{+0.03}$. To derive the dependency of the temperature on the stellar luminosity, $T(R_0)$ for the disks in the sample is then fitted as the power-law function of their central stars' luminosities, $L_\star = L_{\text{int}} + L_{\text{acc}}$, such that $T(L_\star, R_0) = T_0(L_\star/L_\star)^\nu$. This yields the following temperature description:

$$T_{\text{dust}}(L_\star, r) \approx 71 \, \text{K} \left(\frac{L_\star}{L_\odot}\right)^{0.25} \left(\frac{r}{35 \, \text{au}}\right)^{-0.52}. \qquad (6)$$

Next, we compute the dust mass in each disk with the optically thin emission assumption:

$$M_{\text{dust}} = \frac{F_{v}d^{2}}{\kappa_{v}B_{v}(T_{\text{dust}})}, \qquad (7)$$

where $\kappa_v$ is the dust opacity at the observations' frequency. This assumption likely does not hold for the young disks in consideration in our study. Indeed, Fig. 9 indicates that the fluxes of our modelled disks are quite close to the optically thick fluxes, with a difference that is within one order of magnitude. However, we aim to pinpoint how far the results from this simple approximation could be from the masses of the identified disk material. For the dust opacity, we used the frequency-dependent power law approximation, $\kappa_v = \kappa_0(\nu/v_0)^\beta$, with $\beta$ being the dust emissivity index. Here, we adopted the commonly used value $\kappa_0 = 10 \, \text{cm}^2 \, \text{g}^{-1}$ for the normalisation coefficient and $\beta = 1$ for the emissivity index in PPDs (Draine 2006; Andrews & Williams 2005) in accordance with our chosen value of $q = 3.5$ for the grain size distribution. Figure 10 compares the opacity from the power law approximation and the one from OPTOOL. We note that at the wavelength of interest in this study, $\lambda = 0.89$ mm, we get roughly the same value, namely, $\kappa_v \sim 3.5 \, \text{cm}^2 \, \text{g}^{-1}$.
<figure><img src="image_10.png" /><figcaption>Fig. 8. Comparison of different fitting strategies: with data from individual configurations (left and middle panels for C43-3 and C43-6, respectively) and from combined antenna arrays (right panels for the combined C43-3 and C43-6).</figcaption></figure>

indications. C43-6 appears to have the best measurements for $R_i$, although they slightly overestimate the disk radius in general. The combined $R_i$ distributions have a broader spread within the 16th-84th percentile range, which displays quite a flat shape, without any real centralised peak. This may be an indication that $R_i$ might not be a meaningful parameter to study the envelope's structure with the Plummer density profile, but is nonetheless useful when employed in combination with the Gaussian radii, to provide a more accurate measurement of the disk sizes when the Gaussian encounters difficulties locating the right disk structure, as can be seen in the results for C43-6 (lower-right panel of Fig. 7 and lower-middle panel of Fig. 8).

## 6.2. Disk masses

Another important property of the disks is their masses, as these values offer insights into the mass budget for the formation of planets in their later evolutionary stages.

From the central peak intensity that we modelled with the Gaussian profile in Eq. (3), we performed an integration to get the millimeter flux, $F_v$, of the disks. These fluxes are shown in Fig. 9 as a function of the modelled disk radii, $R_{obs} = R_{90\%}$. As a reference, we plotted the fluxes derived analytically from $R_{obs}$ while assuming that the disks are fully optically thick:

$$F_{v,\text{th}} = \frac{1}{d^2} \int_{0}^{R_{\text{obs}}} B_v(T_{\text{dust}}) 2\pi r dr. \quad (5)$$

Here, $d \sim 140$ pc is the distance to the source and $B_v$ is the Planck function. The dust temperature radial profile $T_{\text{dust}}(r)$ for each disk is inferred from the 3D profile in the simulation by fitting a power-law function of $T(r) = T(R_0)(r/R_0)^{-\nu}$ to the radially averaged temperature shown in Fig. A.1. Using the 16th,

50th, and 84th percentiles of the fitted parameters, we found $\bar{T}(R_0) = 181_{-35}^{+74}$ K, $\bar{R}_0 = 35_{-12}^{+36}$ au and $\bar{v} = 0.52_{-0.05}^{+0.03}$. To derive the dependency of the temperature on the stellar luminosity, $T(R_0)$ for the disks in the sample is then fitted as the power-law function of their central stars' luminosities, $L_* = L_{\text{int}} + L_{\text{acc}}$, such that $T(L_*, R_0) = T_0(L_*/L_*)^\nu$. This yields the following temperature description:

$$T_{\text{dust}}(L_*, r) \approx 71 \, \text{K} \left(\frac{L_*}{L_{\odot}}\right)^{0.25} \left(\frac{r}{35 \, \text{au}}\right)^{-0.52} . \quad (6)$$

Next, we compute the dust mass in each disk with the optically thin emission assumption:

$$M_{\text{dust}} = \frac{F_v d^2}{\kappa_v B_v(T_{\text{dust}})}, \quad (7)$$

where $\kappa_v$ is the dust opacity at the observations' frequency. This assumption likely does not hold for the young disks in consideration in our study. Indeed, Fig. 9 indicates that the fluxes of our modelled disks are quite close to the optically thick fluxes, with a difference that is within one order of magnitude. However, we aim to pinpoint how far the results from this simple approximation could be from the masses of the identified disk material. For the dust opacity, we used the frequency-dependent power law approximation, $\kappa_v = \kappa_0(\nu/v_0)^\beta$, with $\beta$ being the dust emissivity index. Here, we adopted the commonly used value $\kappa_0 = 10 \, \text{cm}^2 \, \text{g}^{-1}$ for the normalisation coefficient and $\beta = 1$ for the emissivity index in PPDs (Draine 2006; Andrews & Williams 2005) in accordance with our chosen value of $q = 3.5$ for the grain size distribution. Figure 10 compares the opacity from the power law approximation and the one from OPTOOL. We note that at the wavelength of interest in this study, $\lambda = 0.89$ mm, we get roughly the same value, namely, $\kappa_v \sim 3.5 \, \text{cm}^2 \, \text{g}^{-1}$.
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Integrated millimeter flux of the Gaussian disks as a function of the modelled disk radii from C43-4 data. Black points show the analytical fluxes in the optically thick regime.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Absorption opacity from OPTOOL (dashed black curve) used in Sect. 3 and from the frequency-dependent power law approximation with the powers of 1 and 1.5 (red and blue lines, respectively).</figcaption></figure>

Another factor of uncertainty in the conversion comes from the average dust disk temperature ($T_{\text{dust}}$). This temperature has been, for a long time, taken as a constant value for a specific evolutionary stage (see, e.g. Tychoniec et al. 2020; Andrews & Williams 2005; Pascucci et al. 2016; Ansdell et al. 2016, 2018; Sanchis et al. 2020). More recently, Tobin et al. (2020) used a grid of radiative transfer models with RADMC-3D to sample the parameter space of a flared disk model embedded within a rotating envelope and found that for a protostar of 1 $L_{\odot}$, a 50 au circumstellar disk has an average dust temperature of 43 K; in addition, the general average temperature scales with the luminosity as $\propto L^{0.25}$, as well as with the disk radius, as $\propto R^{-0.5}$. We note that this dependency on the disk radius is close to what we found with the radial temperature profile in our simulation, to the power of -0.52; however, is generally lower by a factor of ~1.3. In the following analyses, we chose to put some restriction on the uncertainty in the flux-to-mass conversion by assuming that we have a good knowledge of either the temperature scaling of our disk sample in relation to the disk radius and to the star's luminosity (e.g. from the literature), but not the radius itself; or, alternatively, of the median temperature of the sample. The latter is aimed at mimicking the use of a constant temperature for the

mass determination of a disk population. With that in mind, we use the disk radius, $R_{\text{obs}} = R_{90\%}$, from the synthetic modelling to estimate the disk temperature based on the power-law fits of the RAMSES radial temperature profiles shown in Fig. A.1. In particular, we employed the temperature scaling in Eq. (6) in parallel with the median temperature of $\bar{T}_{\text{dust}} = 122$ K for the entire sample, which was determined by the median values of the fitted parameters.

Since the disk component is better recovered with the intermediate-resolution single-configuration setup as discussed in the previous sub-section, here we show only the results for C43-4, which gives the best values for the disk sizes (see Fig. 7). We also include the plots for the other configurations in Appendix C. Figure 11 presents the masses derived from this approximation, in comparison with the values inferred from the material that belongs to the disks based on the gas kinematics in the simulation. Similarly to the way we obtained the disk radii, we calculated the ratios of the masses in the modelling over the masses inferred from the simulations, $r_M = M_{\text{obs}}/M_{\text{sim}}$. The left panel shows the values obtained from our modelling, with constant and scaled temperatures versus those derived from the gas kinematics from the simulation; while the right panel the histograms of the logarithms of their ratios $r_M$. The error bars come from the uncertainties in the peak flux, $I_0$, and the flux fraction of the disk, $f$. Unlike the disk radius, the mass cannot be retrieved with a good level of accuracy. We also provided a simple log-linear fit to the data points in the scatter plots using the `linmix`<fnref n="1" /> library (Kelly 2007). In general, we are underestimating the disk mass using the proposed temperature, as indicated by the shaded 16th-84th percentile uncertainty ranges in the histogram in the right panel. Notably, all of the log-linear fits return very shallow slopes (~0.37–0.39), compared to the $x = y$ line which represents the perfect match between the values on the two axes. This $M_{\text{obs}} \propto M_{\text{sim}}^{\alpha}$ where $0 < \alpha < 1$ relation between the observed and the ‘real’ masses would translate into larger negative offsets when we reach very high disk masses, suggesting a discrepancy of one order of magnitude toward the most massive disks (~1 $M_{\odot}$). We note that here we use $R_{90\%}$ for the disk radius to avoid the abuse of prior knowledge of real objects, which slightly overestimates their ‘actual’ sizes and, hence, underestimates their temperatures. Otherwise speaking, the dust temperatures currently used to derive the masses are lower than the ones obtained with the power-law fits given by the black solid lines in Fig. A.1. The ‘correct’ $T_{\text{dust}}$, namely, the one that aptly describes the temperature profile of the disks in the simulation, could in fact reduce a little further the modelled disk masses. Interestingly, the results with the constant $T_{\text{dust}}$ and $T_{\text{dust}} \propto L_{\star}^{\mu} R_{\text{disk}}^{-\nu}$ are not in too much of a disagreement, with the latter even giving a slightly larger spread.

Now, if we relax the assumption on the prior knowledge of the temperature, an alteration of $T_{\text{dust}}$ can shift the masses of the sample upward or downward. For instance, a much lower temperature such as $T_{\text{dust}} = 30$ K used by Tychoniec et al. (2020) or $T_{\text{dust}} = 43$ K for 50 au disks around 1 $L_{\odot}$ protostars in the VLA/ALMA Nascent Disk and Multiplicity (VANDAM) survey by Tobin et al. (2020) would increase $M_{\text{obs}}$ significantly; however, the change would be systematic throughout the population.

In short, the disk masses cannot be easily recovered from the millimeter fluxes. This can be attributed to two possible sources of uncertainty, the conversion formula in Eq. (7) and the optical thickness of the early Class 0/I disks in question.

<footnote n="1">¹ https://github.com/jmeyers314/linmix</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Integrated millimeter flux of the Gaussian disks as a function of the modelled disk radii from C43-4 data. Black points show the analytical fluxes in the optically thick regime.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Absorption opacity from OPTOOL (dashed black curve) used in Sect. 3 and from the frequency-dependent power law approximation with the powers of 1 and 1.5 (red and blue lines, respectively).</figcaption></figure>

Another factor of uncertainty in the conversion comes from the average dust disk temperature ($T_{\text{dust}}$). This temperature has been, for a long time, taken as a constant value for a specific evolutionary stage (see, e.g. Tychoniec et al. 2020; Andrews & Williams 2005; Pascucci et al. 2016; Ansdell et al. 2016, 2018; Sanchis et al. 2020). More recently, Tobin et al. (2020) used a grid of radiative transfer models with RADMC-3D to sample the parameter space of a flared disk model embedded within a rotating envelope and found that for a protostar of $1 L_{\odot}$, a 50 au circumstellar disk has an average dust temperature of 43 K; in addition, the general average temperature scales with the luminosity as $\propto L^{0.25}$, as well as with the disk radius, as $\propto R^{-0.5}$. We note that this dependency on the disk radius is close to what we found with the radial temperature profile in our simulation, to the power of -0.52; however, is generally lower by a factor of ~1.3. In the following analyses, we chose to put some restriction on the uncertainty in the flux-to-mass conversion by assuming that we have a good knowledge of either the temperature scaling of our disk sample in relation to the disk radius and to the star's luminosity (e.g. from the literature), but not the radius itself; or, alternatively, of the median temperature of the sample. The latter is aimed at mimicking the use of a constant temperature for the

mass determination of a disk population. With that in mind, we use the disk radius, $R_{\text{obs}} = R_{90\%}$, from the synthetic modelling to estimate the disk temperature based on the power-law fits of the RAMSES radial temperature profiles shown in Fig. A.1. In particular, we employed the temperature scaling in Eq. (6) in parallel with the median temperature of $\bar{T}_{\text{dust}} = 122$ K for the entire sample, which was determined by the median values of the fitted parameters.

Since the disk component is better recovered with the intermediate-resolution single-configuration setup as discussed in the previous sub-section, here we show only the results for C43-4, which gives the best values for the disk sizes (see Fig. 7). We also include the plots for the other configurations in Appendix C. Figure 11 presents the masses derived from this approximation, in comparison with the values inferred from the material that belongs to the disks based on the gas kinematics in the simulation. Similarly to the way we obtained the disk radii, we calculated the ratios of the masses in the modelling over the masses inferred from the simulations, $r_M = M_{\text{obs}}/M_{\text{sim}}$. The left panel shows the values obtained from our modelling, with constant and scaled temperatures versus those derived from the gas kinematics from the simulation; while the right panel the histograms of the logarithms of their ratios $r_M$. The error bars come from the uncertainties in the peak flux, $I_0$, and the flux fraction of the disk, $f$. Unlike the disk radius, the mass cannot be retrieved with a good level of accuracy. We also provided a simple log-linear fit to the data points in the scatter plots using the `linmix`<fnref n="1" /> library (Kelly 2007). In general, we are underestimating the disk mass using the proposed temperature, as indicated by the shaded 16th-84th percentile uncertainty ranges in the histogram in the right panel. Notably, all of the log-linear fits return very shallow slopes (~0.37–0.39), compared to the $x = y$ line which represents the perfect match between the values on the two axes. This $M_{\text{obs}} \propto M_{\text{sim}}^{\alpha}$ where $0 < \alpha < 1$ relation between the observed and the ‘real’ masses would translate into larger negative offsets when we reach very high disk masses, suggesting a discrepancy of one order of magnitude toward the most massive disks (~1 $M_{\odot}$). We note that here we use $R_{90\%}$ for the disk radius to avoid the abuse of prior knowledge of real objects, which slightly overestimates their ‘actual’ sizes and, hence, underestimates their temperatures. Otherwise speaking, the dust temperatures currently used to derive the masses are lower than the ones obtained with the power-law fits given by the black solid lines in Fig. A.1. The ‘correct’ $T_{\text{dust}}$, namely, the one that aptly describes the temperature profile of the disks in the simulation, could in fact reduce a little further the modelled disk masses. Interestingly, the results with the constant $T_{\text{dust}}$ and $T_{\text{dust}} \propto L_{\star}^{\mu} R_{\text{disk}}^{-\nu}$ are not in too much of a disagreement, with the latter even giving a slightly larger spread.

Now, if we relax the assumption on the prior knowledge of the temperature, an alteration of $T_{\text{dust}}$ can shift the masses of the sample upward or downward. For instance, a much lower temperature such as $T_{\text{dust}} = 30$ K used by Tychoniec et al. (2020) or $T_{\text{dust}} = 43$ K for 50 au disks around $1 L_{\odot}$ protostars in the VLA/ALMA Nascent Disk and Multiplicity (VANDAM) survey by Tobin et al. (2020) would increase $M_{\text{obs}}$ significantly; however, the change would be systematic throughout the population.

In short, the disk masses cannot be easily recovered from the millimeter fluxes. This can be attributed to two possible sources of uncertainty, the conversion formula in Eq. (7) and the optical thickness of the early Class 0/I disks in question.

<footnote n="1">¹ https://github.com/jmeyers314/linmix</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Scatter plots of the modelled disk masses using a constant dust temperature *T* = 123 K for all the disks (in blue) and a temperature scaling based on the disk radius *T*<sub>dust</sub>(*R*<sub>disk</sub>) = 1122.76(*R*<sub>disk</sub>/1 au)<sup>-0.52</sup> K (in purple) versus the masses inferred from the simulation (shown on the left). The log-linear fits to the data points with corresponding 1σ confidence intervals based on the 16th and 84 percentiles are plotted with the solid lines. Histograms of the *M*<sub>obs</sub>/<i>M</i><sub>sim</sub> ratios in logarithmic scale. Vertical dashed blue and purple lines indicate the medians of the ratio distributions (shown on the right).</figcaption></figure>

While the validity of the former is debatable, the latter can be
confirmed by the small difference (within less than one order
of magnitude) between the fully optically thick fluxes and the
fluxes measured from our synthetic observations. The shallow
slopes of the observed-versus-real-mass plots in Fig. 11 also
point to the effects of high optical depth in hindering the esti-
mates of the masses from the fluxes in the case of more massive
disks, which is a confirmed obstacle in the modelling of young
disks from millimeter observations (Maureira et al. 2022;
Sharma et al. 2023).

6.3. Disk inclinations

Another property that we have access to in our modelling is the
disk inclination, *i*, which can give insights into how the 2D pro-
jection of the disk on the plane of the sky affects the observed
fluxes. In practical terms, this can be directly inferred from the
semi-major axis, *a*, and semi-minor axis, *b*, of the projected
elliptic Gaussian disk as cos(*i*) = *a*/*b*. In the simulation, we mea-
sured this quantity by determining the orientation of the angular
momentum and calculated the angle between this vector and the
observation plane.

These results are depicted in Fig. 12. Though most of the
disks in our simulation are close to either edge-on (90°) or face-
on (0°) geometries when seen in the three native projections xy,
yz, and xz. This is due to the numerical tendency for the disks
to align with the Cartesian axes of the grid, they are not cor-
rectly reconstructed in our modelling from observations. The
face-on ones appears much more inclined, probably due to the
fact that the original disks are mostly not circular and an elon-
gated structure with some inclination is needed to account for
the asymmetry and irregular shapes. In other words, the emis-
sion profile of an asymmetric, but face-on disk in the simulation,
when it is modelled as a 2D elliptical Gaussian, may be erro-
neously interpreted in our analysis as a circularly symmetric, but

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Modelled inclinations of the Gaussian disk versus the inclina-
tions of the disks in simulation based on the angular momentum vector
orientation.</figcaption></figure>

inclined disk. This could potentially also be the case for an actual
observation or a 'real' disk modelled with the same method.
The presence of strong non-symmetric structures, such as spi-
rals or streamers, in some cases also interferes significantly with
respect to how the disks' elongation is identified (a typical exam-
ple could be seen in the nearly face-on disk in Fig. B.3, which
is modelled as a larger and much inclined disk from the C43-
4's observation). Thus, such realistically well grounded scenario
supported by observations (see, e.g., Pineda et al. 2020; Tobin
et al. 2016) may have an important impact on the continuum
modelling outcome of young disks using simple analytical pre-
scriptions. Conversely, the edge-on disks are obtained with less
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Scatter plots of the modelled disk masses using a constant dust temperature *T* = 123 K for all the disks (in blue) and a temperature scaling based on the disk radius *T*<sub>dust</sub>(*R*<sub>disk</sub>) = 1122.76(*R*<sub>disk</sub>/1 au)<sup>-0.52</sup> K (in purple) versus the masses inferred from the simulation (shown on the left). The log-linear fits to the data points with corresponding 1σ confidence intervals based on the 16th and 84 percentiles are plotted with the solid lines. Histograms of the *M*<sub>obs</sub>/<i>M</i><sub>sim</sub> ratios in logarithmic scale. Vertical dashed blue and purple lines indicate the medians of the ratio distributions (shown on the right).</figcaption></figure>

While the validity of the former is debatable, the latter can be
confirmed by the small difference (within less than one order
of magnitude) between the fully optically thick fluxes and the
fluxes measured from our synthetic observations. The shallow
slopes of the observed-versus-real-mass plots in Fig. 11 also
point to the effects of high optical depth in hindering the esti-
mates of the masses from the fluxes in the case of more massive
disks, which is a confirmed obstacle in the modelling of young
disks from millimeter observations (Maureira et al. 2022;
Sharma et al. 2023).

6.3. Disk inclinations

Another property that we have access to in our modelling is the
disk inclination, *i*, which can give insights into how the 2D pro-
jection of the disk on the plane of the sky affects the observed
fluxes. In practical terms, this can be directly inferred from the
semi-major axis, *a*, and semi-minor axis, *b*, of the projected
elliptic Gaussian disk as cos(*i*) = *a*/*b*. In the simulation, we mea-
sured this quantity by determining the orientation of the angular
momentum and calculated the angle between this vector and the
observation plane.

These results are depicted in Fig. 12. Though most of the
disks in our simulation are close to either edge-on (90°) or face-
on (0°) geometries when seen in the three native projections xy,
yz, and xz. This is due to the numerical tendency for the disks
to align with the Cartesian axes of the grid, they are not cor-
rectly reconstructed in our modelling from observations. The
face-on ones appears much more inclined, probably due to the
fact that the original disks are mostly not circular and an elon-
gated structure with some inclination is needed to account for
the asymmetry and irregular shapes. In other words, the emis-
sion profile of an asymmetric, but face-on disk in the simulation,
when it is modelled as a 2D elliptical Gaussian, may be erro-
neously interpreted in our analysis as a circularly symmetric, but

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Modelled inclinations of the Gaussian disk versus the inclina-
tions of the disks in simulation based on the angular momentum vector
orientation.</figcaption></figure>

inclined disk. This could potentially also be the case for an actual
observation or a 'real' disk modelled with the same method.
The presence of strong non-symmetric structures, such as spi-
rals or streamers, in some cases also interferes significantly with
respect to how the disks' elongation is identified (a typical exam-
ple could be seen in the nearly face-on disk in Fig. B.3, which
is modelled as a larger and much inclined disk from the C43-
4's observation). Thus, such realistically well grounded scenario
supported by observations (see, e.g., Pineda et al. 2020; Tobin
et al. 2016) may have an important impact on the continuum
modelling outcome of young disks using simple analytical pre-
scriptions. Conversely, the edge-on disks are obtained with less
inclination with respect to the face-on position. This is some-
what expected in our type of modelling due to the impossibility
to obtain an infinitely flat Gaussian emission required for the 90°
edge-on inclination.

Another possible explanation could be the increase in the
apparent height when these disks are viewed from the side, espe-
cially for the points that are too far off in the plot where the
90° inclination is reproduced with a model fit close to face-on
0–20° structures. One example is the increase caused by the
envelope back-warming that can heat the disk to higher temper-
ature (Butner et al. 1991, 1994; D’Alessio et al. 1997; see also
Agurto-Gangas et al. 2019) and potentially enhance its flux in the
perpendicular direction to blur its elongation. In that same vein,
Han et al. (2023) recently attempted to revisit the observational
study of a subset of protostars in Orion from the VANDAM sur-
vey by Karnath et al. (2020), which excluded protostellar disks as
a source of the observed irregular mm/sub-mm emission due to
the lack of elongation in the supposed edge-on emission. In par-
ticular, the authors investigated the effect of the back-warming
of the envelope on the gas pressure scale height of the disk.
They reported that in the presence of an envelope around the
disk, its irradiation onto the disk could significantly increase
its millimeter flux, but only raise the scale height moderately.
Here, we suggest an additional effect which acts on the disk
scale height when looking at the geometry of the emission, also
caused by the envelope back-warming, which could contribute to
the arguments in favor of the disk geometry regarding the pecu-
liar observed structures. It is worth pointing out that, contrary
to the face-on cases, the edge-on results have greater error bars
that indicate the possible range of values is very large, oftentimes
covering the ‘original’ inclination found in the simulation.

7. Discussion

7.1. Implications for disk modelling from ALMA observations

Our analyses with simple modelling of the synthetic ALMA
observations with a Gaussian disk and a surrounding spherical
envelope suggest that with angular resolution enough to resolve
the disk component, we could reproduce their size properties
with a good accuracy, within a factor of 1.6–2.2 in uncertainty.
This idealised model, despite being somewhat simplified, is still
being extensively used in the flux modelling of protoplanetary
disks and protostellar envelopes from continuum observations
(Maury et al. 2019; Cacciapuoti et al. 2023). It is now reassur-
ing that with this method, we can still recover the disk sizes
within acceptable uncertainty relatively well. Notably, the com-
pact configuration C43-4 with which we are simulating for the
close star-forming regions, whose physical resolution is ≈37 au,
gives the best results in term of uncertainty for the disks of radii
between 50 and 150 au, when those disks are modelled with a 2D
Gaussian. This can be scaled for more distant regions to deter-
mine the most suitable antenna configuration and observing time
for observations.

Our further analyses with varied angular resolutions and
combined antenna arrays showed that having not enough angu-
lar resolution, that is, a resolution lower than the ‘actual’ sizes
of the disks will result in over-estimating their radii. Very high
resolutions, on the other hand, do not really improve the accu-
racy of the idealised Gaussian disk models and could potentially
degrade the large-scale component (see the lower right panel
of Fig. 7). A more complex disk model, such as one with a
power-law intensity profile and a rotating circumstellar envelope,

potentially coupled with radiative transfer as in Sheehan et al.
(2022), would probably help in the case of more resolved disks.

The disk masses, on the other hand, are poorly modelled,
even when we have a reasonable estimation of the temperature,
and overall underestimated with the optically thin conversion
from the millimeter fluxes, which is a common method to this
date to derive the mass in the modelling of young disks from
continuum observations (Tychoniec et al. 2018, 2020; Tobin et al.
2020). As it stands, this would challenge our current under-
standing of the mass budget for planet formation. The poor
correlation between the masses of the Gaussian disks and the
masses estimated in the simulation suggests that either the con-
version between the millimeter flux to the dust mass is to be
revised, or the wavelength at which we performed our synthetic
observations needs to be extended, so that the disks are less opti-
cally thick to satisfy the assumption of Eq. (7). The fact that
the fluxes that we modelled for the disks are relatively close
to their fluxes in the optically thick regime (as suggested by
Fig. 9) strongly supports the latter, while not excluding the for-
mer. We believe that observations at longer wavelengths, for
example, those covered by ALMA Band 1 (6–8.6 mm), could
help us to probe better the disk mass, as Tazzari et al. (2021a)
showed that at 3.1 mm, the optically thick fraction of the fluxes
of Class II disks in the Lupus star-forming region is already
significantly reduced compared to at 0.9 mm.

7.2. More simplistically: Considering whether a single Gaussian be adequate

The simplest analytical prescription that often comes to mind when attempting to model the isolated disk emission is probably the single elliptic 2D Gaussian, thanks to its limited number of parameters. This limited set enables modellers to easily fit huge samples of sources in large surveys without massive computational trade-off, at the risk of missing more refined and complex structures. In the literature, the Gaussian has often been used for more evolved, isolated Class II disks (see, e.g. Tazzari et al. 2017; Manara et al. 2019; Sanchis et al. 2020), but sometimes also for embedded disks in younger class I or even class 0 stages (see, e.g. Tobin et al. 2020; Ohashi et al. 2023).

With our simulation and synthetic observation data, we can
substitute the two-component model presented in Sect. 5 by the
single Gaussian and re-employ the fitting methodology to test
the accuracy of this simplistic model. The results for the disk
radii measured with C43-4's data, which provided us with the
best estimations of the disk sizes with the disk-envelope model,
are plotted in Fig. 13. Owing to the limitations of the single
Gaussian in detecting the compact, embedded structures in the
presence of the extended emission of the envelope, most of the
disks (especially the partially resolved ones) are over-estimated
in size; this drives the means and medians of the ratios for all
three radii *R*<sub>68%</sub>, *R*<sub>90%</sub>, and *R*<sub>95%</sub> above 1, with an excess of very
large *R*<sub>90%</sub>/*R*<sub>sim</sub> ~ 10 values on the far right. Combined with the
fact that the 1σ confidence interval of log(*R*<sub>90%</sub>) lying entirely on
the positive side, this is an indication that the Gaussian tends to
bias toward larger disk radii by an average factor of ~2. Thus, this
overly simplistic model is inadequate for the type of embedded
disks formed in our simulations and studied in this paper.

In real applications of such methodology to modelling young
sources where the contribution of both the disk and the enve-
lope are equally important, we caution that the disk sizes may
not be accurately measured. Hence, more complicated prop-
erties such as the masses (provided that the optical thickness
inclination with respect to the face-on position. This is some-
what expected in our type of modelling due to the impossibility
to obtain an infinitely flat Gaussian emission required for the 90°
edge-on inclination.

Another possible explanation could be the increase in the
apparent height when these disks are viewed from the side, espe-
cially for the points that are too far off in the plot where the
90° inclination is reproduced with a model fit close to face-on
0–20° structures. One example is the increase caused by the
envelope back-warming that can heat the disk to higher temper-
ature (Butner et al. 1991, 1994; D’Alessio et al. 1997; see also
Agurto-Gangas et al. 2019) and potentially enhance its flux in the
perpendicular direction to blur its elongation. In that same vein,
Han et al. (2023) recently attempted to revisit the observational
study of a subset of protostars in Orion from the VANDAM sur-
vey by Karnath et al. (2020), which excluded protostellar disks as
a source of the observed irregular mm/sub-mm emission due to
the lack of elongation in the supposed edge-on emission. In par-
ticular, the authors investigated the effect of the back-warming
of the envelope on the gas pressure scale height of the disk.
They reported that in the presence of an envelope around the
disk, its irradiation onto the disk could significantly increase
its millimeter flux, but only raise the scale height moderately.
Here, we suggest an additional effect which acts on the disk
scale height when looking at the geometry of the emission, also
caused by the envelope back-warming, which could contribute to
the arguments in favor of the disk geometry regarding the pecu-
liar observed structures. It is worth pointing out that, contrary
to the face-on cases, the edge-on results have greater error bars
that indicate the possible range of values is very large, oftentimes
covering the ‘original’ inclination found in the simulation.

7. Discussion

7.1. Implications for disk modelling from ALMA observations

Our analyses with simple modelling of the synthetic ALMA
observations with a Gaussian disk and a surrounding spherical
envelope suggest that with angular resolution enough to resolve
the disk component, we could reproduce their size properties
with a good accuracy, within a factor of 1.6–2.2 in uncertainty.
This idealised model, despite being somewhat simplified, is still
being extensively used in the flux modelling of protoplanetary
disks and protostellar envelopes from continuum observations
(Maury et al. 2019; Cacciapuoti et al. 2023). It is now reassur-
ing that with this method, we can still recover the disk sizes
within acceptable uncertainty relatively well. Notably, the com-
pact configuration C43-4 with which we are simulating for the
close star-forming regions, whose physical resolution is ≈37 au,
gives the best results in term of uncertainty for the disks of radii
between 50 and 150 au, when those disks are modelled with a 2D
Gaussian. This can be scaled for more distant regions to deter-
mine the most suitable antenna configuration and observing time
for observations.

Our further analyses with varied angular resolutions and
combined antenna arrays showed that having not enough angu-
lar resolution, that is, a resolution lower than the ‘actual’ sizes
of the disks will result in over-estimating their radii. Very high
resolutions, on the other hand, do not really improve the accu-
racy of the idealised Gaussian disk models and could potentially
degrade the large-scale component (see the lower right panel
of Fig. 7). A more complex disk model, such as one with a
power-law intensity profile and a rotating circumstellar envelope,

potentially coupled with radiative transfer as in Sheehan et al.
(2022), would probably help in the case of more resolved disks.

The disk masses, on the other hand, are poorly modelled,
even when we have a reasonable estimation of the temperature,
and overall underestimated with the optically thin conversion
from the millimeter fluxes, which is a common method to this
date to derive the mass in the modelling of young disks from
continuum observations (Tychoniec et al. 2018, 2020; Tobin et al.
2020). As it stands, this would challenge our current under-
standing of the mass budget for planet formation. The poor
correlation between the masses of the Gaussian disks and the
masses estimated in the simulation suggests that either the con-
version between the millimeter flux to the dust mass is to be
revised, or the wavelength at which we performed our synthetic
observations needs to be extended, so that the disks are less opti-
cally thick to satisfy the assumption of Eq. (7). The fact that
the fluxes that we modelled for the disks are relatively close
to their fluxes in the optically thick regime (as suggested by
Fig. 9) strongly supports the latter, while not excluding the for-
mer. We believe that observations at longer wavelengths, for
example, those covered by ALMA Band 1 (6–8.6 mm), could
help us to probe better the disk mass, as Tazzari et al. (2021a)
showed that at 3.1 mm, the optically thick fraction of the fluxes
of Class II disks in the Lupus star-forming region is already
significantly reduced compared to at 0.9 mm.

7.2. More simplistically: Considering whether a single
Gaussian be adequate

The simplest analytical prescription that often comes to mind when attempting to model the isolated disk emission is probably the single elliptic 2D Gaussian, thanks to its limited number of parameters. This limited set enables modellers to easily fit huge samples of sources in large surveys without massive computational trade-off, at the risk of missing more refined and complex structures. In the literature, the Gaussian has often been used for more evolved, isolated Class II disks (see, e.g. Tazzari et al. 2017; Manara et al. 2019; Sanchis et al. 2020), but sometimes also for embedded disks in younger class I or even class 0 stages (see, e.g. Tobin et al. 2020; Ohashi et al. 2023).

With our simulation and synthetic observation data, we can
substitute the two-component model presented in Sect. 5 by the
single Gaussian and re-employ the fitting methodology to test
the accuracy of this simplistic model. The results for the disk
radii measured with C43-4’s data, which provided us with the
best estimations of the disk sizes with the disk-envelope model,
are plotted in Fig. 13. Owing to the limitations of the single
Gaussian in detecting the compact, embedded structures in the
presence of the extended emission of the envelope, most of the
disks (especially the partially resolved ones) are over-estimated
in size; this drives the means and medians of the ratios for all
three radii *R*<sub>68%</sub>, *R*<sub>90%</sub>, and *R*<sub>95%</sub> above 1, with an excess of very
large *R*<sub>90%</sub>/*R*<sub>sim</sub> ~ 10 values on the far right. Combined with the
fact that the 1σ confidence interval of log(*R*<sub>90%</sub>) lying entirely on
the positive side, this is an indication that the Gaussian tends to
bias toward larger disk radii by an average factor of ~2. Thus, this
overly simplistic model is inadequate for the type of embedded
disks formed in our simulations and studied in this paper.

In real applications of such methodology to modelling young
sources where the contribution of both the disk and the enve-
lope are equally important, we caution that the disk sizes may
not be accurately measured. Hence, more complicated prop-
erties such as the masses (provided that the optical thickness
<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Disk sizes in arcsecs obtained from the modelling of C43-4 observations with a single Gaussian model vs. their radii as inferred from the simulation (left panel) and the distribution of the ratio between the observed and theoretical radii (right panel).</figcaption></figure>

allows for a better conversion from the millimeter fluxes) would consequently be wrongly determined.

### 7.3. Considering the comparison between simulations and observations of planet-forming disks

The overarching goal of the post-processing and synthetic observation study presented in this paper is to provide a better connection between numerical simulations and interferometric observations of young disks in Class 0/I systems, with detailed physics of the transfer of radiation and various instrumental effects included. As a result, the outputs of our pipeline from the simulation results could then be used for the comparison of the disk properties with dust continuum observations. This could also help constraining the initial ingredients of the numerical simulations of disk formation.

Interestingly, there has been significant tension between simulations and observations regarding the theoretically and observationally constrained disk mass (Bate 2018; Tobin et al. 2020; Lebreuilly et al. 2024), although the disk radii show relatively good agreements (Lebreuilly et al. 2021, 2024). This difference by one to two orders of magnitude in the simulated and the observed disk mass adds another layer of complexity to our question, which specifically asks whether the mass budget of the disks is sufficient to form Solar-like planetary systems according to the Minimum Solar Mass Nebula (Hayashi 1981; see also, e.g. Desch 2007) and whether simulations and observations are correctly predicting these masses. Here, we have demonstrated that the masses modelled from the observed millimeter fluxes could be 2.5–10 times lower than their actual masses. By extending this study to different simulation models and observations of various disk populations from various clouds, it is possible to constrain the choice of initial parameters that would be able to reconcile this current mismatch.

Ultimately, our final goal would be to explore, with the standardised methodology presented in this paper, the larger parameter space to put more constraints on the physics in the simulations, and more broadly, theories of the evolution of PPDs. This includes the dust composition and distribution that

determine the dust opacity, which in turn controls the results of the radiative transfer, as well as more complex processes such as dust growth and fragmentation (Lebreuilly et al. 2023). The changes in the opacity could consequently lead to a difference in the way the disk masses would be obtained, although it would not be likely to optimally effect at the current 0.89 mm wavelength, which urges a multi-wavelength extension to other ALMA bands of this study. The disk sizes would not potentially be overly affected due to the constant dust-to-gas ratio assumption of 0.01 throughout the simulation grid, which is another strong caveat that ought to be addressed. In addition, we must not neglect the criteria of the disk extraction from the gas kinematics, which can alter the theoretical disk properties that we are currently taking blindly as the ‘ground truth’. The same goes for the techniques used to model the observed fluxes of the disks in observations, as suggested by the difference between disk sizes obtained with the one-component and two-component models. Different tracers and observation techniques, such as molecular lines, channel maps, and PV diagrams, have also been shown to produce uncertainties in the disk physical quantities. This was demonstrated in detail by Aso & Machida (2020) with regards to the disk size using single core collapse simulation and synthetic observations.

Lastly, we would need also to address other shortcomings of our simulation models, such as the accretion luminosity, $L_{acc}$, that is so far converted from the accretion gravitational energy with an efficiency factor, $f_{acc}$, taken quite arbitrarily as 0.1 in this study or additionally 0.5 in Lebreuilly et al. (2024). The value of $f_{acc}$, which eventually determines $L_{acc}$, has been shown to impact not only the radiative transfer photon energy, but also the disk temperature and fragmentation from a numerical standpoint, even though it has little impact on the disk size and mass. Sink particles, which are a sub-grid modelling of stars due to the lack of resolution to resolve the disk-star interaction, with their own set of fine-tuning parameters including the accretion threshold, $n_{thre}$, and sink radii, contribute further to the uncertainty of the disk properties modelled with our simulations. In particular, Hennebelle et al. (2020) showed that the choice of $n_{thre}$ could significantly alter the disk mass, which has been studied extensively in this paper. As a side note, we
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Disk sizes in arcsecs obtained from the modelling of C43-4 observations with a single Gaussian model vs. their radii as inferred from the simulation (left panel) and the distribution of the ratio between the observed and theoretical radii (right panel).</figcaption></figure>

allows for a better conversion from the millimeter fluxes) would consequently be wrongly determined.

### 7.3. Considering the comparison between simulations and observations of planet-forming disks

The overarching goal of the post-processing and synthetic observation study presented in this paper is to provide a better connection between numerical simulations and interferometric observations of young disks in Class 0/I systems, with detailed physics of the transfer of radiation and various instrumental effects included. As a result, the outputs of our pipeline from the simulation results could then be used for the comparison of the disk properties with dust continuum observations. This could also help constraining the initial ingredients of the numerical simulations of disk formation.

Interestingly, there has been significant tension between simulations and observations regarding the theoretically and observationally constrained disk mass (Bate 2018; Tobin et al. 2020; Lebreuilly et al. 2024), although the disk radii show relatively good agreements (Lebreuilly et al. 2021, 2024). This difference by one to two orders of magnitude in the simulated and the observed disk mass adds another layer of complexity to our question, which specifically asks whether the mass budget of the disks is sufficient to form Solar-like planetary systems according to the Minimum Solar Mass Nebula (Hayashi 1981; see also, e.g. Desch 2007) and whether simulations and observations are correctly predicting these masses. Here, we have demonstrated that the masses modelled from the observed millimeter fluxes could be 2.5–10 times lower than their actual masses. By extending this study to different simulation models and observations of various disk populations from various clouds, it is possible to constrain the choice of initial parameters that would be able to reconcile this current mismatch.

Ultimately, our final goal would be to explore, with the standardised methodology presented in this paper, the larger parameter space to put more constraints on the physics in the simulations, and more broadly, theories of the evolution of PPDs. This includes the dust composition and distribution that

determine the dust opacity, which in turn controls the results of the radiative transfer, as well as more complex processes such as dust growth and fragmentation (Lebreuilly et al. 2023). The changes in the opacity could consequently lead to a difference in the way the disk masses would be obtained, although it would not be likely to optimally effect at the current 0.89 mm wavelength, which urges a multi-wavelength extension to other ALMA bands of this study. The disk sizes would not potentially be overly affected due to the constant dust-to-gas ratio assumption of 0.01 throughout the simulation grid, which is another strong caveat that ought to be addressed. In addition, we must not neglect the criteria of the disk extraction from the gas kinematics, which can alter the theoretical disk properties that we are currently taking blindly as the ‘ground truth’. The same goes for the techniques used to model the observed fluxes of the disks in observations, as suggested by the difference between disk sizes obtained with the one-component and two-component models. Different tracers and observation techniques, such as molecular lines, channel maps, and PV diagrams, have also been shown to produce uncertainties in the disk physical quantities. This was demonstrated in detail by Aso & Machida (2020) with regards to the disk size using single core collapse simulation and synthetic observations.

Lastly, we would need also to address other shortcomings of our simulation models, such as the accretion luminosity, $L_{acc}$, that is so far converted from the accretion gravitational energy with an efficiency factor, $f_{acc}$, taken quite arbitrarily as 0.1 in this study or additionally 0.5 in Lebreuilly et al. (2024). The value of $f_{acc}$, which eventually determines $L_{acc}$, has been shown to impact not only the radiative transfer photon energy, but also the disk temperature and fragmentation from a numerical standpoint, even though it has little impact on the disk size and mass. Sink particles, which are a sub-grid modelling of stars due to the lack of resolution to resolve the disk-star interaction, with their own set of fine-tuning parameters including the accretion threshold, $n_{thre}$, and sink radii, contribute further to the uncertainty of the disk properties modelled with our simulations. In particular, Hennebelle et al. (2020) showed that the choice of $n_{thre}$ could significantly alter the disk mass, which has been studied extensively in this paper. As a side note, we
refer to Lebreuilly et al. (2024) for a complete discussion on the
caveats of the simulation models. Such detailed investigation
will be presented in our follow-up studies.

8. Summary

In this paper, we study two critical disk parameters in Class 0/I
systems, namely, the disk size and mass, by means of realis-
tic radiative post-processing and synthetic ALMA observations
of MHD simulations combined with observational modelling
techniques. Our principal results can be summarised as follows:

1. We developed a procedure to compare disks' observational properties with the properties in the simulations, including 3D Monte Carlo radiative transfer with the RADMC-3D code, mock interferometric observations with the CASA software, and parameter fitting with the galario library.

2. Our analysis of the disk properties from synthetic observations suggests that observations with enough angular resolution to resolve the disks should allow us to recover their sizes with an uncertainties of a factor of ~1.6–2.2.

3. Their masses, on the other hand, are poorly reproduced, with increasingly larger negative offsets towards very high disk masses that could reach an order of magnitude difference at the higher end (~1 $M_\odot$ disks). The potential problem might lie in the optically thin assumption for the conversion from the millimeter flux to the mass. This is due to the fact that at the wavelength we use to simulate the ALMA Band 7 observations ($\lambda$ = 0.89 mm), which corresponds to a frequency of 345 GHz, the dust disks are still optically thick.

4. We varied the angular resolutions of our observations by using three antenna configurations C43-3 ($\theta_{res}$ = 0.41"-57.4 au), C43-4 ($\theta_{res}$ = 0.266"-37.2 au), and C43-6 ($\theta_{res}$ = 0.0887"-12.4 au) and found that different resolutions would return different disk properties. Low resolution (C43-3) would result in overestimating the disk sizes, a sufficient angular resolution (C43-4) to resolve the disks would be the most suitable for the Gaussian modelling of the disk fluxes, whereas overly high-resolution (C43-6) observations might need a more complex description for the disk intensity profile.

5. Another experiment to use the visibility data from two antenna configurations C43-3 and C43-6 simultaneously for the fitting to extend the uv coverage shows improvements in the results obtained for the disk, which suggests that by extending the uv coverage of the observations, we can better separate the compact (disk) and extended (envelope) emission.

6. Compared to the two-component model, a single Gaussian model tends to bias the disk size toward larger radii and overall over-estimate the disk radius by a factor of ~2; therefore, it should be used with caution, especially for the modelling of young sources with significant contribution of the envelope to the total emission.

In conclusion, we have demonstrated, via synthetic observations
of self-consistent MHD simulations of disk population forma-
tion, how accurately the two important observable properties
of Class 0/I disks, their sizes and masses, can or cannot be
retrieved from ALMA observations, as well as the impact of
certain observation designs and data analyses on the modelling
results. Future studies making use of this methodology with
potentially more advanced observational modelling techniques
and greater parameter space would be essential to bridge the gap

between the theoretically motivated numerical simulations and
interferometric observations of young disks.

Acknowledgements. We are grateful to the referee for the very constructive and valuable report that improved immensely the overall quality of the paper, and Agnes Monod-Gayraud for her thorough proof-reading of the text. We thank Sylvie Cabrit for her suggestion on the inclusion of the single Gaussian fitting assessment, and Cornelis Dullemond for his tremendous help with the RADMC-3D code. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 823823 (DUSTBUSTERS) and from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130). This work was partly supported by the Italian Ministero dell'Istruzione, Università e Ricerca through the grant Progetti Premiali 2012 – iALMA (CUP C52113000140001). We acknowledge PRACE for awarding us access to the JUWELS supercomputer. M.G. acknowledges the support of the French Agence Nationale de la Recherche (ANR) through the project COSMHIC (ANR-20-CE31-0009). R.S.K. acknowledges financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 - 390900948) “STRUCTURES”, and from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 50002206). R.S.K. also thanks for computing resources provided by the Ministry of Science, Research and the Arts (MWK) of the State of Baden-Württemberg through bwHPC and the German Science Foundation (DFG) through grants INST 35/1134-1 FUGG and 35/1597-1 FUGG, and for data storage at SDS@hd funded through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG. G.R. acknowledges the support from the Fondazione Cariplo, grant no. 2022-1217, and the European Research Council (ERC) under the European Union's Horizon Europe Research & Innovation Programme under grant agreement no. 101039651 (DiscEvol). Views and opinions expressed are however those of the author(s) only, and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

References

Adams, F. C. 1991, ApJ, 382, 544

Adams, F. C. 2010, ARA&A, 48, 47

Agurto-Gangas, C., Pineda, J. E., SzHucs, L., et al. 2019, A&A, 623, A147

Andre, P., Ward-Thompson, D., & Barsony, M. 2000, in Protostars and Planets IV, eds. V. Mannings, A. P. Boss, & S. S. Russell (Tucson: University of Arizona Press), 59

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

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

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

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

Aso, Y., & Machida, M. N. 2020, ApJ, 905, 174

Bate, M. R. 2018, MNRAS, 475, 5618

Berger, M. J., & Oliger, J. 1984, J. Computat. Phys., 53, 484

Bianchi, E., López-Sepulcre, A., Ceccarelli, C., et al. 2023, Faraday Discuss., 245, 164

Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015

Briggs, D. S. 1995, PhD thesis, New Mexico Institute of Mining and Technology, Mexico

Butner, H. M., Evans, N. J., I., Lester, D. F., Levreault, R. M., & Strom, S. E. 1991, ApJ, 376, 636

Butner, H. M., Natta, A., & Evans, Neal J., I. 1994, ApJ, 420, 326

Cacciapuoti, L., Macias, E., Maury, A. J., et al. 2023, A&A, 676, A4

CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501

Clark, B. G. 1980, A&A, 89, 377

Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35

Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11

D'Alessio, P., Calvet, N., & Hartmann, L. 1997, ApJ, 474, 397

Desch, S. J. 2007, ApJ, 671, 878

Dominik, C., Min, M., & Tazaki, R. 2021, Astrophysics Source Code Library [record ascl:2104.010]

Draine, B. T. 2006, ApJ, 636, 1114

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, Astrophysics Source Code Library [record ascl:1202.015]

Evans, Neal J., I., Di Francesco, J., Lee, J.-E., et al. 2015, ApJ, 814, 22

Facchini, S., van Dishoeck, E. F., Manara, C. F., et al. 2019, A&A, 626, L2

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

Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371

Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Computat. Sci., 5, 65
refer to Lebreuilly et al. (2024) for a complete discussion on the caveats of the simulation models. Such detailed investigation will be presented in our follow-up studies.

## 8. Summary

In this paper, we study two critical disk parameters in Class 0/I systems, namely, the disk size and mass, by means of realistic radiative post-processing and synthetic ALMA observations of MHD simulations combined with observational modelling techniques. Our principal results can be summarised as follows:

1. We developed a procedure to compare disks' observational properties with the properties in the simulations, including 3D Monte Carlo radiative transfer with the RADMC-3D code, mock interferometric observations with the CASA software, and parameter fitting with the galario library.

2. Our analysis of the disk properties from synthetic observations suggests that observations with enough angular resolution to resolve the disks should allow us to recover their sizes with an uncertainties of a factor of ~1.6–2.2.

3. Their masses, on the other hand, are poorly reproduced, with increasingly larger negative offsets towards very high disk masses that could reach an order of magnitude difference at the higher end (~1 $M_\odot$ disks). The potential problem might lie in the optically thin assumption for the conversion from the millimeter flux to the mass. This is due to the fact that at the wavelength we use to simulate the ALMA Band 7 observations ($\lambda$ = 0.89 mm), which corresponds to a frequency of 345 GHz, the dust disks are still optically thick.

4. We varied the angular resolutions of our observations by using three antenna configurations C43-3 ($\theta_{res}$ = 0.41"–57.4 au), C43-4 ($\theta_{res}$ = 0.266"–37.2 au), and C43-6 ($\theta_{res}$ = 0.0887"–12.4 au) and found that different resolutions would return different disk properties. Low resolution (C43-3) would result in overestimating the disk sizes, a sufficient angular resolution (C43-4) to resolve the disks would be the most suitable for the Gaussian modelling of the disk fluxes, whereas overly high-resolution (C43-6) observations might need a more complex description for the disk intensity profile.

5. Another experiment to use the visibility data from two antenna configurations C43-3 and C43-6 simultaneously for the fitting to extend the uv coverage shows improvements in the results obtained for the disk, which suggests that by extending the uv coverage of the observations, we can better separate the compact (disk) and extended (envelope) emission.

6. Compared to the two-component model, a single Gaussian model tends to bias the disk size toward larger radii and overall over-estimate the disk radius by a factor of ~2; therefore, it should be used with caution, especially for the modelling of young sources with significant contribution of the envelope to the total emission.

In conclusion, we have demonstrated, via synthetic observations of self-consistent MHD simulations of disk population formation, how accurately the two important observable properties of Class 0/I disks, their sizes and masses, can or cannot be retrieved from ALMA observations, as well as the impact of certain observation designs and data analyses on the modelling results. Future studies making use of this methodology with potentially more advanced observational modelling techniques and greater parameter space would be essential to bridge the gap

between the theoretically motivated numerical simulations and interferometric observations of young disks.

**Acknowledgements.** We are grateful to the referee for the very constructive and valuable report that improved immensely the overall quality of the paper, and Agnes Monod-Gayraud for her thorough proof-reading of the text. We thank Sylvie Cabrit for her suggestion on the inclusion of the single Gaussian fitting assessment, and Cornelis Dullemond for his tremendous help with the RADMC-3D code. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 823823 (DUSTBUSTERS) and from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130). This work was partly supported by the Italian Ministero dell'Istruzione, Università e Ricerca through the grant Progetti Premiali 2012 – iALMA (CUP C52113000140001). We acknowledge PRACE for awarding us access to the JUWELS supercomputer. M.G. acknowledges the support of the French Agence Nationale de la Recherche (ANR) through the project COSMHIC (ANR-20-CE31-0009). R.S.K. acknowledges financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 - 390900948) "STRUCTURES", and from the German Ministry for Economic Affairs and Climate Action in project "MAINN" (funding ID 50002206). R.S.K. also thanks for computing resources provided by the Ministry of Science, Research and the Arts (MWK) of the State of Baden-Württemberg through bwHPC and the German Science Foundation (DFG) through grants INST 35/1134-1 FUGG and 35/1597-1 FUGG, and for data storage at SDS@hd funded through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG. G.R. acknowledges the support from the Fondazione Cariplo, grant no. 2022-1217, and the European Research Council (ERC) under the European Union's Horizon Europe Research & Innovation Programme under grant agreement no. 101039651 (DiscEvol). Views and opinions expressed are however those of the author(s) only, and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

## References

Adams, F. C. 1991, ApJ, 382, 544

Adams, F. C. 2010, ARA&A, 48, 47

Agurto-Gangas, C., Pineda, J. E., SzHucs, L., et al. 2019, A&A, 623, A147

Andre, P., Ward-Thompson, D., & Barsony, M. 2000, in Protostars and Planets IV, eds. V. Mannings, A. P. Boss, & S. S. Russell (Tucson: University of Arizona Press), 59

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

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

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

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

Aso, Y., & Machida, M. N. 2020, ApJ, 905, 174

Bate, M. R. 2018, MNRAS, 475, 5618

Berger, M. J., & Oliger, J. 1984, J. Computat. Phys., 53, 484

Bianchi, E., López-Sepulcre, A., Ceccarelli, C., et al. 2023, Faraday Discuss., 245, 164

Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015

Briggs, D. S. 1995, PhD thesis, New Mexico Institute of Mining and Technology, Mexico

Butner, H. M., Evans, N. J., I., Lester, D. F., Levreault, R. M., & Strom, S. E. 1991, ApJ, 376, 636

Butner, H. M., Natta, A., & Evans, Neal J., I. 1994, ApJ, 420, 326

Cacciapuoti, L., Macias, E., Maury, A. J., et al. 2023, A&A, 676, A4

CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501

Clark, B. G. 1980, A&A, 89, 377

Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35

Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11

D'Alessio, P., Calvet, N., & Hartmann, L. 1997, ApJ, 474, 397

Desch, S. J. 2007, ApJ, 671, 878

Dominik, C., Min, M., & Tazaki, R. 2021, Astrophysics Source Code Library [record ascl:2104.010]

Draine, B. T. 2006, ApJ, 636, 1114

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, Astrophysics Source Code Library [record ascl:1202.015]

Evans, Neal J., I., Di Francesco, J., Lee, J.-E., et al. 2015, ApJ, 814, 22

Facchini, S., van Dishoeck, E. F., Manara, C. F., et al. 2019, A&A, 626, L2

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

Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371

Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Computat. Sci., 5, 65
Han, F., Hartmann, L., Calvet, N., & Franco-Hernández, R. 2023, MNRAS, 526, 2566

Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67

Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128

Karnath, N., Megeath, S. T., Tobin, J. J., et al. 2020, ApJ, 890, 129

Kelly, B. C. 2007, ApJ, 665, 1489

Koeperfer, C. M., & Robitaille, T. P. 2017, ApJ, 849, 3

Kuiper, R., & Yorke, H. W. 2013, ApJ, 772, 61

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10

Lebreuilly, U., Vallucci-Goy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2024, A&A, 682, A30

Lee, Y.-N., Charnoz, S., & Hennebelle, P. 2021, A&A, 648, A101

Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321

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

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

Maureira, M. J., Gong, M., Pineda, J. E., et al. 2022, ApJ, 941, L23

Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127

Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2020, A&A, 635, A42

Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155

Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, ASP Conf. Ser., 534, 501

Ohashi, S., Kobayashi, H., Sai, J., & Sakai, N. 2022, ApJ, 933, 23

Ohashi, N., Tobin, J. J., Jørgensen, J. K., et al. 2023, ApJ, 951, 8

Pardo, J. R. 2019, in ALMA Development Workshop, 36

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

Pfalzner, S., Davies, M. B., Gounelle, M., et al. 2015, Phys., Scrip., 90, 068001

Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nat. Astron., 4, 1158

Pineda, J. E., Arzoumanian, D., Andre, P., et al. 2023, ASP Conf. Ser., 534, 233

Robitaille, T. P. 2010, A&A, 520, A70

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

Sadavoy, S. I., Myers, P. C., Stephens, I. W., et al. 2018, ApJ, 869, 115

Sanchis, E., Testi, L., Natta, A., et al. 2020, A&A, 633, A114

Segura-Cox, D. M., Looney, L. W., Tobin, J. J., et al. 2018, ApJ, 866, 161

Sharma, R., Jørgensen, J. K., Gavino, S., et al. 2023, ApJ, 954, 69

Sheehan, P. D., Tobin, J. J., Looney, L. W., & Megeath, S. T. 2022, ApJ, 929, 76

Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23

Tazzari, M., Testi, L., Natta, A., et al. 2017, A&A, 606, A88

Tazzari, M., Beaujean, F., & Testi, L. 2018, MNRAS, 476, 4527

Tazzari, M., Clarke, C. J., Testi, L., et al. 2021a, MNRAS, 506, 2804

Tazzari, M., Testi, L., Natta, A., et al. 2021b, MNRAS, 506, 5117

Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 339

Teyssier, R. 2002, A&A, 385, 337

Tobin, J. J., Looney, L. W., Wilner, D. J., et al. 2015, ApJ, 805, 125

Tobin, J. J., Kratter, K. M., Persson, M. V., et al. 2016, Nature, 538, 483

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

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

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

Tripathi, A., Andrews, S. M., Birnstiel, T., & Wilner, D. J. 2017, ApJ, 845, 44

Tsukamoto, Y., Maury, A., Commercon, B., et al. 2023, ASP Conf. Ser., 534, 317

Tychoniec, Ł., Tobin, J. J., Karska, A., et al. 2018, ApJS, 238, 19

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

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

Winter, A. J., Clarke, C. J., Rosotti, G., et al. 2018, MNRAS, 478, 2700

Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103
Han, F., Hartmann, L., Calvet, N., & Franco-Hernández, R. 2023, MNRAS, 526, 2566

Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67

Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128

Karnath, N., Megeath, S. T., Tobin, J. J., et al. 2020, ApJ, 890, 129

Kelly, B. C. 2007, ApJ, 665, 1489

Koeperfer, C. M., & Robitaille, T. P. 2017, ApJ, 849, 3

Kuiper, R., & Yorke, H. W. 2013, ApJ, 772, 61

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10

Lebreuilly, U., Vallucci-Goy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2024, A&A, 682, A30

Lee, Y.-N., Charnoz, S., & Hennebelle, P. 2021, A&A, 648, A101

Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321

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

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

Maureira, M. J., Gong, M., Pineda, J. E., et al. 2022, ApJ, 941, L23

Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127

Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2020, A&A, 635, A42

Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155

Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, ASP Conf. Ser., 534, 501

Ohashi, S., Kobayashi, H., Sai, J., & Sakai, N. 2022, ApJ, 933, 23

Ohashi, N., Tobin, J. J., Jørgensen, J. K., et al. 2023, ApJ, 951, 8

Pardo, J. R. 2019, in ALMA Development Workshop, 36

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

Pfalzner, S., Davies, M. B., Gounelle, M., et al. 2015, Phys., Scrip., 90, 068001

Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nat. Astron., 4, 1158

Pineda, J. E., Arzoumanian, D., Andre, P., et al. 2023, ASP Conf. Ser., 534, 233

Robitaille, T. P. 2010, A&A, 520, A70

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

Sadavoy, S. I., Myers, P. C., Stephens, I. W., et al. 2018, ApJ, 869, 115

Sanchis, E., Testi, L., Natta, A., et al. 2020, A&A, 633, A114

Segura-Cox, D. M., Looney, L. W., Tobin, J. J., et al. 2018, ApJ, 866, 161

Sharma, R., Jørgensen, J. K., Gavino, S., et al. 2023, ApJ, 954, 69

Sheehan, P. D., Tobin, J. J., Looney, L. W., & Megeath, S. T. 2022, ApJ, 929, 76

Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23

Tazzari, M., Testi, L., Natta, A., et al. 2017, A&A, 606, A88

Tazzari, M., Beaujean, F., & Testi, L. 2018, MNRAS, 476, 4527

Tazzari, M., Clarke, C. J., Testi, L., et al. 2021a, MNRAS, 506, 2804

Tazzari, M., Testi, L., Natta, A., et al. 2021b, MNRAS, 506, 5117

Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 339

Teyssier, R. 2002, A&A, 385, 337

Tobin, J. J., Looney, L. W., Wilner, D. J., et al. 2015, ApJ, 805, 125

Tobin, J. J., Kratter, K. M., Persson, M. V., et al. 2016, Nature, 538, 483

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

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

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

Tripathi, A., Andrews, S. M., Birnstiel, T., & Wilner, D. J. 2017, ApJ, 845, 44

Tsukamoto, Y., Maury, A., Commercon, B., et al. 2023, ASP Conf. Ser., 534, 317

Tychoniec, Ł., Tobin, J. J., Karska, A., et al. 2018, ApJS, 238, 19

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

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

Winter, A. J., Clarke, C. J., Rosotti, G., et al. 2018, MNRAS, 478, 2700

Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103
## Appendix A: RADMC-3D and RAMSES temperature profiles

To obtain the temperature profile with RADMC-3D, we set up the thermal Monte Carlo (MC) runs with $10^8$ photon packages for all the sources present in the specific RAMSES output that we want to post-process. Due to the fact that the very few cells around the sink particles are a region of overdense gas and hence high optical depth ($\tau \sim 10$), the photon packages from the position of the star (of point-source type in RADMC-3D) would have to reprocess multiple times before being able to get out of this region, with the risk of being stuck there for a large amount of time. This translates into a long runtime for the MC runs. Therefore, the modified random walk (MRW; Min et al. 2009; Robitaille 2010) mode is switched on to speed up the code. Additionally, we chose an embarrassingly parallel approach as an alternative to the parallelisation implemented in the code since its efficiency decreases when increasing the number of computer's threads in use. This is done by dividing the total $10^8$ photon packages into $N$ non-parallelised runs, each with $10^8/N$ photon packages and different MC seeds, to calculate the dust temperature. Then, we averaged the temperature profiles from these $N$ runs according to their radiant heat energy:

$$T_{\text{avg}}^4 = \frac{1}{N} \sum_{i=1}^{N} T_i^4, \qquad (A.1)$$

where $T_{\text{avg}}$ is the final average dust temperature for each cell in the grid and $T_i$ is the dust temperature in the same cell from the $i$-th non-parallelised run. We used $N = 100$ runs to achieve the temperature profile with low computational cost.

On the other hand, RAMSES, with its radiative transfer solver using the gray FLD method, provides us with the temperature profile of the gas and dust in thermal equilibrium within the simulation grid. This radiative transfer approach includes approximation of the diffusion of radiation, local production of heat, and gray approximation of the dust opacity. Unlike in RADMC-3D, the energy from the central source is smoothed by a Gaussian kernel centred around the star position and distributed within the radius of the sink particle ($r_{\text{sink}} = 4\Delta_{\text{max}} \sim 5$ au), and each cell emits like a black body. This is a possible explanation for the difference usually observed within the inner most 10 au of the disks, as suggested in the temperature profiles obtained with the two codes shown in Fig. A.1, which causes difficulties in the parameter fitting and modelling process as the Gaussian is mistakenly fitted to the over-fluxed region instead of the disk. The lack of dust sublimation in RADMC-3D, which is implemented in RAMSES could also account for the discrepancy. The same difference was reported in Mignon-Risse et al. (2020) for their model with high opacity $\tau = 100$ and high stellar temperature $T_{\star} = 15000$ K (see Fig. 4 in Mignon-Risse et al. 2020).

We continue and use both temperature profiles, which we found to be relatively similar in the outer disk (where $r \ge 10$ au) and we performed the ray-tracing, mock interferometric observations, and synthetic modelling to extract the disks in the fiducial case with C43-4. Our goal is to test whether the final results change significantly when changing the profiles as the input for the imaging processes. Figure A.2 compares the disk sizes obtained with the temperature from RADMC-3D and RAMSES. In most cases, the results are identical, except for the disks in which the central flux resulted from RADMC-3D temperature profile is mistaken by a narrow Gaussian, and very few other disks which are not properly modelled either with RADMC-3D or RAMSES profiles.

We note, however, that with C43-6, using RAMSES would only eliminate a good portion of the problematic cases, but not all. In other words, some disks still have the central flux mistaken by the Gaussian component. In those very few cases, we could alternatively use the $R_i$ parameter of the envelope, which describes the disk-envelope structure well, compared to the value obtained from the simulation (see, e.g. Fig. 8), as the radius of the disk.
## Appendix A: RADMC-3D and RAMSES temperature profiles

To obtain the temperature profile with RADMC-3D, we set up the thermal Monte Carlo (MC) runs with $10^8$ photon packages for all the sources present in the specific RAMSES output that we want to post-process. Due to the fact that the very few cells around the sink particles are a region of overdense gas and hence high optical depth ($\tau \sim 10$), the photon packages from the position of the star (of point-source type in RADMC-3D) would have to reprocess multiple times before being able to get out of this region, with the risk of being stuck there for a large amount of time. This translates into a long runtime for the MC runs. Therefore, the modified random walk (MRW; Min et al. 2009; Robitaille 2010) mode is switched on to speed up the code. Additionally, we chose an embarrassingly parallel approach as an alternative to the parallelisation implemented in the code since its efficiency decreases when increasing the number of computer's threads in use. This is done by dividing the total $10^8$ photon packages into $N$ non-parallelised runs, each with $10^8/N$ photon packages and different MC seeds, to calculate the dust temperature. Then, we averaged the temperature profiles from these $N$ runs according to their radiant heat energy:

$$T_{\text{avg}}^4 = \frac{1}{N} \sum_{i=1}^{N} T_i^4, \qquad (A.1)$$

where $T_{\text{avg}}$ is the final average dust temperature for each cell in the grid and $T_i$ is the dust temperature in the same cell from the $i$-th non-parallelised run. We used $N = 100$ runs to achieve the temperature profile with low computational cost.

On the other hand, RAMSES, with its radiative transfer solver using the gray FLD method, provides us with the temperature profile of the gas and dust in thermal equilibrium within the simulation grid. This radiative transfer approach includes approximation of the diffusion of radiation, local production of heat, and gray approximation of the dust opacity. Unlike in RADMC-3D, the energy from the central source is smoothed by a Gaussian kernel centred around the star position and distributed within the radius of the sink particle ($r_{\text{sink}} = 4\Delta_{\text{max}} \sim 5$ au), and each cell emits like a black body. This is a possible explanation for the difference usually observed within the inner most 10 au of the disks, as suggested in the temperature profiles obtained with the two codes shown in Fig. A.1, which causes difficulties in the parameter fitting and modelling process as the Gaussian is mistakenly fitted to the over-fluxed region instead of the disk. The lack of dust sublimation in RADMC-3D, which is implemented in RAMSES could also account for the discrepancy. The same difference was reported in Mignon-Risse et al. (2020) for their model with high opacity $\tau = 100$ and high stellar temperature $T_{\star} = 15000$ K (see Fig. 4 in Mignon-Risse et al. 2020).

We continue and use both temperature profiles, which we found to be relatively similar in the outer disk (where $r \ge 10$ au) and we performed the ray-tracing, mock interferometric observations, and synthetic modelling to extract the disks in the fiducial case with C43-4. Our goal is to test whether the final results change significantly when changing the profiles as the input for the imaging processes. Figure A.2 compares the disk sizes obtained with the temperature from RADMC-3D and RAMSES. In most cases, the results are identical, except for the disks in which the central flux resulted from RADMC-3D temperature profile is mistaken by a narrow Gaussian, and very few other disks which are not properly modelled either with RADMC-3D or RAMSES profiles.

We note, however, that with C43-6, using RAMSES would only eliminate a good portion of the problematic cases, but not all. In other words, some disks still have the central flux mistaken by the Gaussian component. In those very few cases, we could alternatively use the $R_i$ parameter of the envelope, which describes the disk-envelope structure well, compared to the value obtained from the simulation (see, e.g. Fig. 8), as the radius of the disk.
<figure><img src="image_14.png" /><figcaption>**Fig. A.1.** Radial temperature profiles from RADMC-3D (red curves) and RAMSES (blue curves). The difference is more pronounced within the sink particle radius ($r \le 4$ au). Dashed, dotted, and dash-dotted lines show the profiles smoothed by a Gaussian whose FWHM equal to that of the three antenna configurations in this study. The shaded areas marked the resolution limit corresponding to said configurations. The power-law fits to be used for the analytical calculations of the temperature are shown by the solid black line.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. A.2.** Scatter plots for $R_{90\%}$ extracted with the temperature profiles from RADMC-3D (green dots) and RAMSES (orange triangles) for all the single disks in one population versus the disk sizes extracted from the simulation. The results are similar in most cases, except notably for the disks where the central flux produced by RADMC-3D are modelled as the very narrow Gaussian.</figcaption></figure>

## Appendix B: Synthetic observations and modelling results for the disk population

Here, we present the results of our modelling from synthetic observations with CASA with three antenna configuration C43-3, C43-4, and C43-6 for all the disks around single stars seen in one projection *xy*, in the same fashion as Fig. 5. It is worth noting that the selection of single systems is based on the images of the observations obtained with tclean, not the column density plot from RAMSES or the sky model from RADMC-3D; that is, without any knowledge of the real objects.
<figure><img src="image_16.png" /><figcaption>**Fig. A.1.** Radial temperature profiles from RADMC-3D (red curves) and RAMSES (blue curves). The difference is more pronounced within the sink particle radius ($r \le 4$ au). Dashed, dotted, and dash-dotted lines show the profiles smoothed by a Gaussian whose FWHM equal to that of the three antenna configurations in this study. The shaded areas marked the resolution limit corresponding to said configurations. The power-law fits to be used for the analytical calculations of the temperature are shown by the solid black line.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. A.2.** Scatter plots for $R_{90\%}$ extracted with the temperature profiles from RADMC-3D (green dots) and RAMSES (orange triangles) for all the single disks in one population versus the disk sizes extracted from the simulation. The results are similar in most cases, except notably for the disks where the central flux produced by RADMC-3D are modelled as the very narrow Gaussian.</figcaption></figure>

## Appendix B: Synthetic observations and modelling results for the disk population

Here, we present the results of our modelling from synthetic observations with CASA with three antenna configuration C43-3, C43-4, and C43-6 for all the disks around single stars seen in one projection *xy*, in the same fashion as Fig. 5. It is worth noting that the selection of single systems is based on the images of the observations obtained with tclean, not the column density plot from RAMSES or the sky model from RADMC-3D; that is, without any knowledge of the real objects.
<figure><img src="image_16.png" /><figcaption>Fig. B.1. Same as Fig. 5, but for sink 17.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. B.1. Same as Fig. 5, but for sink 17.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. B.2. Same as Fig. 5, but for sink 20.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. B.2. Same as Fig. 5, but for sink 20.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. B.3.** Same as Fig. 5, but for sink 30.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. B.3. Same as Fig. 5, but for sink 30.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. B.4. Same as Fig. 5, but for sink 42.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. B.4. Same as Fig. 5, but for sink 42.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. B.5.** Same as Fig. 5, but for sink 50.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. B.5. Same as Fig. 5, but for sink 50.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. B.6. Same as Fig. 5, but for sink 52.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. B.6. Same as Fig. 5, but for sink 52.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. B.7. Same as Fig. 5, but for sink 53.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. B.7. Same as Fig. 5, but for sink 53.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. B.8. Same as Fig. 5, but for sink 65.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. B.8. Same as Fig. 5, but for sink 65.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. B.9. Same as Fig. 5, but for sink 67.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. B.9. Same as Fig. 5, but for sink 67.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. B.10.** Same as Fig. 5, but for sink 70. For configurations C43-3 and C43-4, the imaged disk appear as a single system; therefore, it is included in the statistic, whereas for C43-6 it is not the case.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Fig. B.10.** Same as Fig. 5, but for sink 70. For configurations C43-3 and C43-4, the imaged disk appear as a single system; therefore, it is included in the statistic, whereas for C43-6 it is not the case.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. B.11. Same as Fig. 5, but for sink 78.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. B.11. Same as Fig. 5, but for sink 78.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. B.12. Same as Fig. 5, but for sink 79.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Fig. B.12. Same as Fig. 5, but for sink 79.</figcaption></figure>
**Appendix C: Disk masses modelled from C43-3's
and C43-6's data**

Below, we show the disk masses we modelled from the obser-
vations with C43-3 and C43-6. Similarly to the C43-4 results,
the masses are not correctly measured, with shallow slope and
large spread in the log-linear fits. With a disk temperature loosely
derived from the simulation temperature profiles, we tend to
underestimate the mass by a factor of 2 – 3 (and in extreme cases:
10).

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

**Fig. C.1.** Scatter plots of the modelled disk masses versus the ones inferred from the simulation for C43-3 (left) and C43-6 (right) using the constant median (in blue) and variable disk temperature (in purple), shown at the top. Solid colored show the log-linear fits for the data points and the shaded areas indicate the corresponding confidence intervals. Histograms of the *M*<sub>obs</sub>/*M*<sub>sim</sub> ratios in logarithmic scale. Vertical dashed lines indicate the medians of the ratio distributions, shown at the bottom.
## Appendix C: Disk masses modelled from C43-3's and C43-6's data

Below, we show the disk masses we modelled from the observations with C43-3 and C43-6. Similarly to the C43-4 results, the masses are not correctly measured, with shallow slope and large spread in the log-linear fits. With a disk temperature loosely derived from the simulation temperature profiles, we tend to underestimate the mass by a factor of 2–3 (and in extreme cases: 10).

<figure><img src="image_30.png" /><figcaption>**Fig. C.1.** Scatter plots of the modelled disk masses versus the ones inferred from the simulation for C43-3 (left) and C43-6 (right) using the constant median (in blue) and variable disk temperature (in purple), shown at the top. Solid colored show the log-linear fits for the data points and the shaded areas indicate the corresponding confidence intervals. Histograms of the $M_{\text{obs}}/M_{\text{sim}}$ ratios in logarithmic scale. Vertical dashed lines indicate the medians of the ratio distributions, shown at the bottom.</figcaption></figure>