Accuracy of ALMA estimates of young disk radii and masses Predicted observations from numerical simulations

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 R ADMC -3 D and C ASA , 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.


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 kilometersized 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. (2021Lebreuilly et al. ( , 2023a) ) 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. 2023a).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 ⊙ clump in Lebreuilly et al. (2023a) 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;Tazzari et al. 2021b) 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.

Simulation models
We begin by briefly presenting the simulation model (NMHD-F01 of Lebreuilly et al. 2023a) 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 −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 µ = 10, which corresponds to a magnetic field of ∼ 9.4 × 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 ⋆ , and radius, R ⋆ , accretes mass and the incoming kinetic energy of the gas is radiated away.This luminosity has a value of where f acc is the fraction of the accretion gravitational energy radiated away and G is the gravitational constant.We adopt f 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 ⋆ 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 = 150M ⊙ .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 = 100M ⊙ , 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.

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;Commerçon et al. 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 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.

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 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 θ 2 res 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.

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;Tazzari et al. 2021b).This experiment allowed us to obtain their observable properties, namely: disk sizes, masses, and inclination angles.

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, σ disk , as: and that of the circularly symmetric Plummer envelope truncated at outer radius, R out , with a power-law distribution of the density of ρ(r) ∝ r−m and temperature of T (r) ∝ r−n (Adams 1991) as: 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 .

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 abovementioned 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 χ 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 20000 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 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.

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.

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;Toci et al. 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 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 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, R i , and truncated outer radius, R out .
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 semimajor 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: 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, θ res /(2 The Gaussian disk size is below the resolution of the observation and the rotation-to-infall transition is above, namely, 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, 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, f I 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 σ 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.

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.

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 χ 2 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, σ disk < θ res /(2 √ 2 ln 2) or R i < θ res /(2 √ 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).

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: 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 int + L acc , such that T (L ⋆ , R 0 ) = T 0 (L ⋆ /L ⋆ ) u .This yields the following temperature description: Next, we compute the dust mass in each disk with the optically thin emission assumption: 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, κ ν = κ 0 (ν/ν 0 ) β , with β being the dust emissivity index.Here, we adopted the commonly used value κ 0 = 10 cm 2 /g 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.[cm 2 /g] 0.89 mm OpTool = 0 ( / 0 ) , = 1 = 1.5 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).
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 2 /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, Ansdell et al. 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 powerlaw fits of the Ramses radial temperature profiles shown in Fig.
A. In particular, we employed the temperature scaling in Eq. ( 6) in parallel with the median temperature of Tdust = 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 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 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 loglinear fit to the data points in the scatter plots using the linmix1 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 obs ∝ M α sim 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 (∼ 1M ⊙ ).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. 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 dust ∝ L u ⋆ R −v disk 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 C43-4 (0.27 37.2 au) T = 122 K: log(y) = 0.39 +0.11 0.10 log(x) 0.97 +0.12 0.11 , 2 = 0.11 +0.04 0.03 T L 0.25 R 0.52 disk : log(y) = 0.37 +0.13 0.13 log(x) 0.99 +0.15  counts log(M obs /M sim ), T = 122 K: mean = -0.39,median = -0.41,std = 0.45 log(M obs /M sim ), T L 0.25 R 0.52 disk : mean = -0.41,median = -0.39,std = 0.53 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 dust (R disk ) = 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 M obs /M sim ratios in logarithmic scale.Vertical dashed blue and purple lines indicate the medians of the ratio distributions (shown on the right).
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).

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 faceon (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 erro- neously 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;Butner et al. 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.

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 starforming regions, whose physical resolution is ≈ 37 au, gives the best results in term of uncertainty for the disks of radii between 50 − 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 powerlaw 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;Tobin et al. 2020;Tychoniec 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.

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 proper-   ties such as the masses (provided that the optical thickness allows for a better conversion from the millimeter fluxes) would consequently be wrongly determined.

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. 2023a), although the disk radii show relatively good agreements (Lebreuilly et al. 2021;Lebreuilly et al. 2023a).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 simula-tions, 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. 2023b).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, f acc , taken quite arbitrarily as 0.1 in this study or additionally 0.5 in Lebreuilly et al. (2023a).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. (2023a) for a complete discussion on the caveats of the simulation models.Such detailed investigation will be presented in our follow-up studies.

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 (∼ 1M ⊙ 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.4 au), C43-4 (θ res = 0.266 ′′ − 37.2 au), and C43-6 (θ 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.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. 1.Column density map integrated along the z−direction of the NMHD-F01 model from Lebreuilly et al. (2023a) at the simulation snapshot used in this study, seen at a full simulation box's length.The star symbols indicate the sink particle positions.

FigFig. 2 .
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 sim , and the modelled disk radii, R 90% , (detailed in Sect.6) are noted for each configuration and projection.
Article number, page 6 of 32 Ngo-Duy Tung et al.: Synthetic Observations of Protoplanetary Disk Simulations

Fig. 4 .
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 .
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.

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

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).
) 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 0 )(r/R 0 ) −v to the radially averaged temperature shown in Fig. A. Using the 16th, 50th, and 84th percentiles of the fitted parameters, we found T (R 0 ) = 181 +74 −35 K, R0 = 35 +36 −12 au and v = 0.52 +0.03 −0.05 .

Fig. 9 .
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. 12 .
Fig. 12. Modelled inclinations of the Gaussian disk versus the inclinations of the disks in simulation based on the angular momentum vector orientation.

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