© The Authors 2024

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 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 Ata-cama 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_⊙ 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 Δx_max ~ 1.2 au. The starting point is a 1000 M_⊙ clump of temperature 10 K, initial radius ~0.38 pc, and uniform density 3 × 10⁻¹⁹g cm⁻³, initialised with supersonic turbulent velocities at 𝓜 = 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 µ = 10, which corresponds to a magnetic field of ~9.4 × 10⁻⁵ G. Sink particles (Bleuler & Teyssier 2014), which form when the density reaches n_thre = 10¹³ cm⁻³, 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_⋆, and radius, R_⋆, 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{\mathcal{G}{M}_{\star }{\dot{M}}_{\star }}{R_{\star }}$$(1)

where ƒ_acc is the fraction of the accretion gravitational energy radiated away and 𝓖 is the gravitational constant. We adopt ƒ_acc = 01 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_⋆ and its internal luminosity L_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_sink = 150 M_⊙. 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_sink = 100 M_⊙, 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.

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.

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_tot = L_acc + L_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_min = 1 nm and the maximum size a_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 ≳ 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 λ = 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⁸ 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.

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, Rsim, and the modelled disk radii, R90%, (detailed in Sect. 6) are noted for each configuration and projection.

4 Mock interferometric observations

The RADMC-3D maps provide us with the ideal images of different regions of the sky surrounding each individual proto-star. From these maps, by introducing proper instrumentational effects, we can produce realistic observations of the simulated sky as would otherwise be obtained with interferometers. To that end, we use thde 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 θ_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 }_{\text{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 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.

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 σdisk is surrounded by a circular envelope following a Plummer-like density profile with rotation-to-infall transitional radius, Ri, and truncated outer radius, Rout.

5.1 Two-component model for the disk and the envelope

In our synthetic Class 0/1 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₀ 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, σ_disk, as: $$I_{\nu ,\text{ disk }}(\overline{r})=f{I}_0\exp \left(-\frac{{\overline{r}}^2}{2{\sigma }_{\text{disk }}^2}\right),$$(2)

of ρ(r̄)∝r̄^−m and temperature of T(r̄)∝r̄⁻ⁿ (Adams 1991) as: $$I_{v,\text{ env }}(\overline{r})=\left\{\begin{array}{ll}\frac{(1-f)I_0}{{\left(1+{\left(\overline{r}/R_i\right)}^2\right)}^{-(m+n-1)/2}}\hfill & \text{ for }0\le \overline{r}\le R_{\text{out }}\hfill \\ 0\hfill & \text{ for }\overline{r}>R_{\text{out }}\hfill \end{array}\right\}.$$(3)

Here, 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 χ² 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σ 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 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.

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.

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.

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: υ_ϕ > 2υ_r,

• Rotating faster than falling vertically: υ_ϕ > 2υ_z,

• Composed of dense material of gas number density n > n_thre = 10⁹ cm⁻³.

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, σ_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 σ_disk as: $$\begin{array}{l}{R}_{68\%}={\sigma }_{\text{disk }}\times \sqrt{-2\times \ln (1-0.68)}\simeq 1.51\times {\sigma }_{\text{disk }},\\ R_{90\%}\simeq 1.42\times R_{68\%},\\ R_{95\%}\simeq 1.62\times R_{68\%},\end{array}$$(4)

In this two component model, since R_i and σ_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 }_{\text{res }}/(2\sqrt{2\ln 2})<{\sigma }_{\text{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 }_{\text{disk }}<{\theta }_{\text{res }}/(2\sqrt{2\ln 2})<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, σ_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 }_{\text{disk }},R_i<{\theta }_{\text{res }}/(2\sqrt{2\ln 2})$.

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, ƒI₀). 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₀ 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 R90% 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 σ_disk and to each other by a certain factor, the standard deviation σ_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 by a horizontal shift equal to the shifts of the median values indicated 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 σ_ratio ~ 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_68%/R_sim and R_90%/R_sim 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_68% and R_90% distributions have medians equally close to 1 in general, with R_68% slightly underestimates and R_90% slightly overestimates the disk radii. To include as much as possible the disk emission, we use R_90% for the histograms of the ratio of the disk radii, as the values for R_68% and R_95% can be related to it by a certain factor.

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 R90%/Rsim ratio for the three projections, shown at the bottom. Black dotted line marks R90%/Rsim = 1, where the two are in full agreement. Vertical dashed lines show the medians of the ratio distributions of three radii R68% (red), R90% (green), and R95% (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.

6.1.2 Effect of resolution

Here, we focus on how different resolutions of the observations could affect the disk sizes obtained with our visibility modelling. For that purpose, we used the same model and fitting method on the data obtained with C43-3 and C43-6, with angular resolutions of θ_res = 0.41″, 0.0887″, respectively. Figure 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 projections into one single sample for each configuration to maximise the sample size. From C43-3 to C43-4, as the angular resolution of the observation increases, the mean and median of the log(R_90%/R_sim) 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, θ_res = 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_90%/R_sim 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 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.

Fig. 7 Comparison of the disk sizes measured from the observations with three different angular resolutions (0.41″, 0.266″, 0.887″). R90% 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 R90%/Rsim ratio for the reported configurations, shown at the bottom.

6.1.3 Combined antenna arrays

With these observations of different resolutions in hand, it is possible to alternatively combine the data during the fitting to have a better coverage of the uv space and better separate the compact 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_90%, and the radius of the approximately uniform density region of the envelope, R_i, 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_90% and R_i 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_90% and R_i returned by the fits are much smaller than the angular resolution of the observations (marked by the dotted red lines). We consider those, namely, ${\sigma }_{\text{disk }}<{\theta }_{\text{res }}/(2\sqrt{2\ln 2})$ or $R_i<{\theta }_{\text{res}}/(2\sqrt{2\ln 2})$, to be cases where the fitting failed to detect the disk (in the case of σ_disk) or the transition between rotation and infall motions (in the case of R_i). 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_i for the disk size (as described in the sub-sections above) since the values of R_i for these disks are above the observation’s resolution, 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 observations. Consequently, this could help the Gaussian in detecting the correct disk structure, even when the statistics for the envelope’s R_i does not gain any enhancement. As a matter of fact, for C43-3, most of the results for R_i fall below the resolution of the configuration. Thus, the distributions for C43-3 consist of very small samples which do not give meaningful statistical 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).

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).

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}{\displaystyle {\int }_0^{R_{\text{obs}}}{B}_v}\left(T_{\text{dust }}\right)2\pi r\text{d}r.$$(5)

Here, d ~ 140 pc is the distance to the source and B_ν is the Planck function. The dust temperature radial profile T_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₀)(r/R₀)^–v to the radially averaged temperature shown in Fig. A.1. Using the 16th, 50th, and 84th percentiles of the fitted parameters, we found $\overline{T}\left(R_0\right)={181}_{-35}^{+74} \text{K},{\overline{R}}_0={35}_{-12}^{+36}\text{ au and }\overline{v}={0.52}_{-0.05}^{+0.03}$. To derive the dependency of the temperature on the stellar luminosity, T(R0) for the disks in the sample is then fitted as the power-law function of their central stars’ luminosities, L_⋆ = L_int + L_acc, such that T(L_⋆,R₀) = T₀(L_⋆/L_⋆)^u. This yields the following temperature description: $$T_{\text{dust }}\left(L_{\star },r\right)\simeq 71 \text{K}{\left(\frac{L_{\star }}{L_{\odot }}\right)}^{0.25}{\left(\frac{r}{35\text{\hspace{0.17em}au}}\right)}^{-0.52}.$$(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\left(T_{\text{dust }}\right)}\text{, }$$(7)

where κ_ν 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, κ_ν = κ₀(ν/ν₀)^β, with β being the dust emissivity index. Here, we adopted the commonly used value κ₀ = 10 cm² g^–1 for the normalisation coefficient and β = 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, λ = 0.89 mm, we get roughly the same value, namely, κ_ν ~ 3.5 cm² g⁻¹.

Another factor of uncertainty in the conversion comes from the average dust disk temperature (T_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_⊙, a 50 au circumstellar disk has an average dust temperature of 43 K; in addition, the general average temperature scales with the luminosity as ∝L^0.25, as well as with the disk radius, as ∝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_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 ${\overline{T}}_{\text{dust }}=122\hspace{0.17em}\text{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_obs/M_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₀, 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 thelinmix¹ 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 < α < 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_⊙). 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_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_dust and $T_{\text{dust }}\propto L_{\star }^u{R}_{\text{disk }}^{-v}$ 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_dust can shift the masses of the sample upward or downward. For instance, a much lower temperature such as T_dust = 30 K used by Tychoniec et al. (2020) or T_dust = 43 K for 50 au disks around 1 L_⊙ protostars in the VLA/ALMA Nascent Disk and Multiplicity (VANDAM) survey by Tobin et al. (2020) would increase M_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.

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 estimates 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).

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.

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).

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 Tdust(Rdisk) = 1122.76(Rdisk/1 au)−0.52 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 Mobs/Msim ratios in logarithmic scale. Vertical dashed blue and purple lines indicate the medians of the ratio distributions (shown on the right).

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 projection 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 measured 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 correctly 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 elongated structure with some inclination is needed to account for the asymmetry and irregular shapes. In other words, the emission profile of an asymmetric, but face-on disk in the simulation, when it is modelled as a 2D elliptical Gaussian, may be erroneously interpreted in our analysis as a circularly symmetric, but 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 spirals or streamers, in some cases also interferes significantly with respect to how the disks’ elongation is identified (a typical example 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 prescriptions. Conversely, the edge-on disks are obtained with less inclination with respect to the face-on position. This is somewhat 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, especially 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 temperature (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 survey 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 particular, 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 peculiar 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.

Fig. 12 Modelled inclinations of the Gaussian disk versus the inclinations of the disks in simulation based on the angular momentum vector orientation.

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 reassuring that with this method, we can still recover the disk sizes within acceptable uncertainty relatively well. Notably, the compact 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 determine 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 angular 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 accuracy 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 understanding 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 conversion 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 optically 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 former. 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_68%, R_90%, and R_95% above 1, with an excess of very large R_90%/R_sim ~ 10 values on the far right. Combined with the fact that the 1σ confidence interval of log(R_90%) 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 envelope are equally important, we caution that the disk sizes may not be accurately measured. Hence, more complicated properties such as the masses (provided that the optical thickness allows for a better conversion from the millimeter fluxes) would consequently be wrongly determined.

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).

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 obser-vationally 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 wold 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, ƒ_acc, taken quite arbitrarily as 0.1 in this study or additionally 0.5 in Lebreuilly et al. (2024). The value of ƒ_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 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_⊙ 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 (λ = 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 (θ_res = 0.41″–57.4au), C43-4 (θ_res = 0.266″37.2au), and C43-6 (θ_res = 0.0887"–12.4au) 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 Sklodowska-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 C52I13000140001). 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 50OO2206). 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.

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⁸ 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 (τ ~ 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 diving the total 10⁸ photon packages into N non-parallelised runs, each with 10⁸/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}{\displaystyle \sum _{i=1}^N{T}_i^4},$$(A.1)

where T_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 difussion 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_sink = 4Δ_max ~ 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 τ = 100 and high stellar temperature T_⋆ = 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 ≳ 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.

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 ≲ 4 au). Dashed, dotted, and dash-dotted lines shows 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.

Fig. A.2 Scatter plots for R90% 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.

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.

Fig. B.1 Same as Fig. 5, but for sink 17.

Fig. B.2 Same as Fig. 5, but for sink 20.

Fig. B.3 Same as Fig. 5, but for sink 30.

Fig. B.4 Same as Fig. 5, but for sink 42.

Fig. B.5 Same as Fig. 5, but for sink 50.

Fig. B.6 Same as Fig. 5, but for sink 52.

Fig. B.7 Same as Fig. 5, but for sink 53.

Fig. B.8 Same as Fig. 5, but for sink 65.

Fig. B.9 Same as Fig. 5, but for sink 67.

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.

Fig. B.11 Same as Fig. 5, but for sink 78.

Fig. B.12 Same as Fig. 5, but for sink 79.

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).

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 Mobs/Msim ratios in logarithmic scale. Vertical dashed lines indicate the medians of the ratio distributions, shown at the bottom.

References

All Figures

	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.
In the text

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, Rsim, and the modelled disk radii, R90%, (detailed in Sect. 6) are noted for each configuration and projection.

In the text

	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 σdisk is surrounded by a circular envelope following a Plummer-like density profile with rotation-to-infall transitional radius, Ri, and truncated outer radius, Rout.
In the text

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.

In the text

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.

In the text

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 R90%/Rsim ratio for the three projections, shown at the bottom. Black dotted line marks R90%/Rsim = 1, where the two are in full agreement. Vertical dashed lines show the medians of the ratio distributions of three radii R68% (red), R90% (green), and R95% (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.

In the text

	Fig. 7 Comparison of the disk sizes measured from the observations with three different angular resolutions (0.41″, 0.266″, 0.887″). R90% 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 R90%/Rsim ratio for the reported configurations, shown at the bottom.
In the text

	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).
In the text

	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.
In the text

	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).
In the text

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 Tdust(Rdisk) = 1122.76(Rdisk/1 au)−0.52 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 Mobs/Msim ratios in logarithmic scale. Vertical dashed blue and purple lines indicate the medians of the ratio distributions (shown on the right).

In the text

	Fig. 12 Modelled inclinations of the Gaussian disk versus the inclinations of the disks in simulation based on the angular momentum vector orientation.
In the text

	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).
In the text

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 ≲ 4 au). Dashed, dotted, and dash-dotted lines shows 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.

In the text

	Fig. A.2 Scatter plots for R90% 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.
In the text

	Fig. B.1 Same as Fig. 5, but for sink 17.
In the text

	Fig. B.2 Same as Fig. 5, but for sink 20.
In the text

	Fig. B.3 Same as Fig. 5, but for sink 30.
In the text

	Fig. B.4 Same as Fig. 5, but for sink 42.
In the text

	Fig. B.5 Same as Fig. 5, but for sink 50.
In the text

	Fig. B.6 Same as Fig. 5, but for sink 52.
In the text

	Fig. B.7 Same as Fig. 5, but for sink 53.
In the text

	Fig. B.8 Same as Fig. 5, but for sink 65.
In the text

	Fig. B.9 Same as Fig. 5, but for sink 67.
In the text

	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.
In the text

	Fig. B.11 Same as Fig. 5, but for sink 78.
In the text

	Fig. B.12 Same as Fig. 5, but for sink 79.
In the text

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 Mobs/Msim ratios in logarithmic scale. Vertical dashed lines indicate the medians of the ratio distributions, shown at the bottom.

In the text