© The Authors 2025

1 Introduction

Star formation, particularly in high-mass star-forming regions, is a multi-scale process (e.g., Vázquez-Semadeni et al. 2019; Liu et al. 2023; Yang et al. 2023). It develops within giant molecular clouds and filaments, where the assembly of parsec-scale clumps takes place (e.g., Molinari et al. 2010), and then progresses through a hierarchy of nested structures on scales as small as ∼1000 AU (e.g., Sanhueza et al. 2019; Svoboda et al. 2019; Motte et al. 2022; Traficante et al. 2023; Morii et al. 2024; Liu et al. 2024; Xu et al. 2024; Coletta et al. 2025).

A crucial aspect of this process is represented by the physical mechanisms that regulate the collapse of clumps of mass of the order of ∼10³ M_⊙, ultimately leading to the formation of clusters of new stars, including future high-mass stars (e.g., Lada & Lada 2003; Henning et al. 2010; Beuther et al. 2018; Pouteau et al. 2022). Observational studies and numerical simulations have been employed to quantify the impact of gravity, turbulence, magnetic fields, and feedback mechanisms – as well as that of the potential interplay of these forces – on clump fragmentation properties. However, results from these works have not been conclusive, encompassing a range of viable star-formation models beyond the long-established dichotomy between the core-fed (e.g., McKee & Tan 2003; Tan et al. 2014) and the clump-fed scenarios (e.g., Bonnell et al. 1997, 2001; Klessen & Burkert 2000; Bonnell & Bate 2006; Girichidis et al. 2011, 2012; Hennemann et al. 2012; Vázquez-Semadeni et al. 2019; Peretto et al. 2020; Vázquez-Semadeni et al. 2024b).

The most recent studies involving high-mass star-forming regions exhibit shared evidence favoring a scenario in which fragments are not isolated from the environment. While some works report thermal Jeans fragmentation as the leading starformation mechanism (e.g., Palau et al. 2015; Sanhueza et al. 2019, Beuther et al. 2021; Morii et al. 2024, Schisano et al. (submitted)), others indicate that additional support from turbulence, magnetic pressure, and feedback mechanisms or a continuous mass accretion onto the fragments must be in place during clump collapse (e.g., Ballesteros-Paredes et al. 2007; Zhang et al. 2009; Hennebelle et al. 2011; Pillai et al. 2011; Wang et al. 2011; Fontani et al. 2018; Traficante et al. 2023; Xu et al. 2023; Vázquez-Semadeni et al. 2024a). Magnetic regulation, in particular, seems to be crucial in the case of massive clumps that exhibit a few fragments and/or the presence of super-Jeans fragments (e.g., Zhang et al. 2015; Pillai et al. 2016; Hennebelle et al. 2022; Pattle et al. 2023); yet some works have indicated only partial agreement or conflicting evidences (e.g., Palau et al. 2021; Beuther et al. 2024).

Although conflicting pieces of evidence are often due to the diversity of the explored environments, they may also be due to the spatial scales probed in the various works (e.g., Louvet et al. 2021). Indeed, the relative importance of gravitational energy, magnetic field strength, turbulence, and radiative pressure might vary across scales, leading to scale-dependent fragmentation patterns regulated by the dominant force at the scale under analysis, and across environments (e.g., Fontani et al. 2018; Tang et al. 2019; Motte et al. 2022; Wang et al. 2024). While defining a universal model for clump fragmentation can be particularly challenging under these circumstances, the comparison between observations and simulations has proven to be a useful tool for constraining star-formation mechanisms, especially when involving synthetic observations (e.g., Krumholz et al. 2007b; Goodman 2011; Offner & Arce 2014; Koepferl et al. 2017; Svoboda et al. 2019; Jáquez-Domínguez et al. 2023).

With a specific focus on viable scenarios of clump fragmentation, Padoan et al. (2020) used synthetic observations to test the core collapse and the competitive accretion models against Herschel and Atacama Large Millimeter/submillimeter Array (ALMA) observations. They proposed the inertial-inflow model to describe turbulent fragmentation and its potential connection to the origin of massive stars. Fontani et al. (2018) carried out analogous studies to analyze the impact of the interplay of turbulence and magnetic fields within simulated clumps. Comparing their results with the fragmentation properties of massive clumps observed with ALMA, they suggest that higher levels of turbulence produce more fragments and that higher magnetization favors the formation of filamentary structures. At the same time, works by Malinen et al. (2011), Harries et al. (2017), and Padoan et al. (2023), for example, have proven the utility of synthetic observations in evaluating technical aspects concerning real observations. They used synthetic data to quantify the reliability of clump, core, and protostar mass estimates derived from submillimeter and infrared dust emission. Mairs et al. (2014), Beuther et al. (2018), and Bonanomi et al. (2024) used synthetic observations to further address the difference in sensitivity to the more extended emission related to single-dish and inter-ferometeric observations, including the ability of detecting the smallest structures and the non-negligible impact of large-scale interferometric filtering.

Despite the advent of synthetic observations have led to significant advancements in how the star-formation process is understood (e.g., Haworth et al. 2018; Rosen et al. 2020), the majority of works involving synthetic data are still subject to some limitations. In particular, the theoretical setups are typically tuned to mimic a specific environment and/or target source, hindering the chance to perform a broad and self-consistent exploration of the parameter space and therefore failing to capture the variety of real star-forming systems (e.g., McKee & Ostriker 2007; Hennebelle & Falgarone 2012; Klessen & Glover 2016). The comparison between observations and synthetic data itself represents an issue. Observational features caused by the resolution limit and by spatial filtering, as well as by the presence of background and interferometric noise, are often not well reproduced in the synthetic observations, precluding a one-to-one comparison with available real data. To effectively address the physics governing star-forming sites at the scales probed by the observations, it is therefore necessary to apply a rigorous post-processing routine to the simulations and ensure that these also mimic the observational features (instrumental effects and observational biases) of the dataset they will be compared against. While several studies in the literature employ post-processing routines that incorporate customized radiative transfer models and observational feature integration, they primarily address either the physics of isolated cores and the disk scales (e.g., Krumholz et al. 2007a; Commerçon et al. 2012; Dipierro et al. 2015; Szulágyi et al. 2019; Curone et al. 2022; Maury et al. 2022; Redaelli et al. 2024; Tung et al. 2024) or the large-scale interstellar medium (e.g., Smith et al. 2014; Duarte-Cabral et al. 2015; Sormani et al. 2018; Juvela et al. 2019). A self-consistent and systematic comparison between actual observations and dust continuum emission synthetic observations, which could provide a quantitative characterization of how the initial conditions of clumps and environment affect the observed clump fragmentation properties, is currently missing in this panorama. To address this problem, we developed the Rosetta Stone project¹: an end-to-end (simulation ⇔ observations) framework designed to systematically produce realistic synthetic observations of high-mass clump fragmentation and compare them with real interferometric data.

In the present work, we provide a detailed description of the Rosetta Stone framework. In Section 2, we present our method, and we focus on the strategy designed for the production of synthetic observations. In Section 3, we report on the physical properties of fragments identified in the synthetic fields. Section 4 is dedicated to the analysis of the fragmentation properties (fragmentation level, mass accretion onto the fragments, relative distances between the fragments), and the impact of the initial conditions used in the models. In Section 5, we compare evidence from synthetic observations with both models and actual observations. Finally, in Section 6 we draw the most relevant conclusions, and we provide insights into future research directions.

2 The Rosetta Stone project

Below, we introduce the Rosetta Stone post-processing pipeline aimed at producing realistic synthetic observations (see Section 2.2) starting from the output of radiative magnetohydro-dynamical (RMHD) simulations (see Section 2.1). The pipeline is designed to be flexible and can already be adapted to mimic any continuum survey performed with ALMA, the Very Large Array (VLA), and the Northern Extended Millimetre Array (NOEMA) interferometers. In Fig. 1, we provide the flowchart of the Rosetta Stone end-to-end framework and its main building blocks, which are, respectively:

1. The reference observations – our first case study takes as reference ALMA 1.3 mm continuum dust emission observations from the Star formation in QUiescent And Luminous Objects (SQUALO) survey by Traficante et al. (2023). The SQUALO sample consists of 13 high-mass star-forming regions of the Galaxy at different evolutionary stages, and all have evidence of infall motions at the clump scales.

2. The suite of RMHD simulations – we rely on a new set of 24 radiative magnetohydrodynamical simulations of high-mass clumps fragmentation produced by means of the RAMSES code by Lebreuilly et al. 2025 (hereafter Paper I).

3. The radiative transfer – our simulations are post-processed at Spitzer 24 µm, Herschel and ALMA Band 6 wavelengths by means of the RADMC-3D code by Tung et al. 2025 (hereafter Paper II).

4. The synthetic observations – after the radiative transfer, our simulations are post-processed by means of the CASA software (McMullin et al. 2007) to mimic the SQUALO survey observing strategy (ACA + 12 m array) and data reduction. In this paper we focus on the systematic production of synthetic observations of continuum dust emission, coherently with the reference observations.

5. The source extraction and photometry – we identify compact sources and we compute the photometry using the Hyper software (Traficante et al. 2015).

6. The common method of analysis and the comparison of the real and synthetic datasets – we compare the fragmentation properties recovered from the analysis of the synthetic sources with those obtained from the analysis of ALMA observations in Traficante et al. (2023).

7. The refinement of models and follow-up observations – the comparison of the two datasets allows us to better tune the initial conditions of the models and optimize the observational setup and strategy for follow-up observations.

Fig. 1 Flowchart of the Rosetta Stone end-to-end framework: scheme of the comparison between observations and simulations by means of the production of realistic synthetic observations. Real data and inherent information are highlighted in green, synthetic data (RMHD simulations, radiative transfer, synthetic observations) in purple, and the steps of the analysis performed on both real and synthetic maps in orange. In the lighter-shaded boxes we specify how the building blocks of the Rosetta Stone framework have been adapted in this work for the SQUALO project (Traficante et al. 2023) science case.

2.1 Setup: Simulations and radiative transfer

The RMHD simulations of high-mass clump fragmentation have been newly developed within the framework of the Rosetta Stone project and are described in detail in Paper I. Specifically, we rely on a first suite of simulations labeled RS1.0, together with a supplementary collection of models used to explore the effect of different numerical prescriptions and to tune the fiducial setup. The RS1.0 suite is part of a planned series of simulations specifically tailored to investigate the formation of high-mass protoclusters at the sub-parsec scales.

Our simulations solve the three-dimensional equations for the conservation of gas mass, momentum, total energy, and radiation total energy, as well as the evolution of the magnetic field under the assumption of ideal MHD. In these models, which are computed by means of the adaptive mesh refinement finite-volume code RAMSES (Teyssier 2002; Fromang et al. 2006), individual forming stars are represented by sink particles (Bleuler & Teyssier 2014) that form self-consistently above a threshold of n_thre = 10⁹cm⁻³.

Since our ultimate goal is to mimic the observations of real star-forming high-mass clumps in the early phases of their evolution, we simulate the collapse of uniform clouds initialized with values of mass, radius, and temperature chosen in agreement with observations. These values correspond to clump masses of 500 and 1000 M_⊙, a radius of ∼0.4 pc (e.g., Urquhart et al. 2014; Traficante et al. 2017; Elia et al. 2021), and a temperature of 10 K, as suggested by the peak temperature derived for prestellar clumps in the Galaxy (e.g., Elia et al. 2021). To also reproduce a realistic clump environment, we initialize our models including self-gravity, stellar radiation due to mass accretion onto the sink particles, turbulence, and magnetic fields. In particular, we explore two different initial levels of turbulence quantified by the Mach number: $$\mathcal{M}=\frac{v}{c_s},$$(1)

where v is the gas speed, and c_s is the speed of sound in the medium. According to the results from Traficante et al. (2017), we set M to 7 and 10, respectively. These values are also supported by the analysis of turbulence in clumps performed by e.g., Rygl et al. (2010), and Colman et al. (2024).

Given our limited understanding of the magnetic field impact on fragmentation (from observations alone), we explore three potential scenarios: a quasi-hydrodynamical one, a moderately magnetized one, and a more magnetized one. The clump magnetization of each scenario is quantified by the normalized parameter mass-over-flux to critical-mass-over-flux ratio (also known as normalized mass-to-flux ratio): $$\mu =(\frac{M_0}{\varphi })/(\frac{M}{\varphi }{)}_c,$$(2)

where M₀ corresponds to the mass of the clump, and $(M/\varphi {)}_c=(0.53/\pi )\sqrt{5/G}$ corresponds to the critical mass-to-flux ratio at which the clump is magnetically stabilized against gravitational collapse (Mouschovias & Spitzer 1976). The stronger the magnetic field, the lower is µ for a clump of fixed mass. We set µ to ∼3, 10 and 100, respectively (Paper I). As initial condition, the magnetic field is uniform and oriented along the z-axis.

As anticipated above, this first suite of simulations has been developed including feedback in the form of radiation from luminosity accretion. This approach enables us to isolate the individual contributions and focus on the impact of turbulence and magnetic fields on clump fragmentation properties. Moreover, the computational cost savings achieved through this choice of feedback mechanisms allow us to explore a significantly broader range of parameters (Paper I). Outflows and Hii regions are planned to be incorporated in a future set of Rosetta Stone models.

With fixed physical initial conditions, we let each model evolve in time with two different turbulent seeds, 1 and 2 respectively, to account for the stochastic fluctuations in the observed fragmentation process. We also test different time-step samplings across the various realizations, resulting in a varying number of saved time steps per realization. A larger collection of time steps is available for the realizations initialized with Seed 2. Each computational time step is associated with a specific time in years and a sink formation efficiency (SFE) computed as: $$\text{SFE}=\frac{M_{sinks}}{(M_{gas}+M_{sinks})},$$(3)

where $M_{sinks}=\sum _{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{k}\mathrm{s}}{M}_{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{k}}$ and (M_gas + M_sinks) is the initial mass in the simulated box. As showed in Paper II, the SFE is the most reliable parameter for defining the evolutionary stage of the clumps, providing a direct comparison with observational parameters such as L/M. To sample an evolution range of ∼10⁵ years from the formation of the first sink particle, we stop the calculations at a SFE of ∼15% and ∼30% for the M500 ² and M1000 realizations, respectively (Paper I).

For each time step, three orthogonal projections are analyzed. These are fixed for all realizations, and correspond to actual projections in the sky along the x, y, and z axes. Since our observational vantage point is fixed in real observations, such a study is of crucial relevance to quantify the impact of the line-of-sight projection effects on the analysis of the clump fragmentation properties. In total, we rely on a set of 24 available combinations of initial conditions (see Table 1), corresponding to a statistically significant sample of 244 data cubes and 732 2D fields. For each data cube, we have access to the density, velocity, and B fields, and to the gas temperature. For this first work, we focus on the 2D information relative to the projected density fields and projected temperatures. The column density maps also contain information about the projected positions of the sink particles as shown in Fig. 2a.

Following the philosophy of the Rosetta Stone project, performing a comparison between this set of simulations and our reference observations requires two additional steps from our end-to-end pipeline. These consist in post-processing physical quantities from the simulations into observable radiation and incorporating observational features into the simulated data to generate realistic synthetic observations, respectively. At this stage, we have to choose which set of observations we want to compare our models to. This is fundamental to replicate the emission recovered at a given wavelength, for a simulated field placed at a given heliocentric distance, as if it was observed with a specific facility.

In this paper, we employ as reference the ALMA 1.3 mm continuum dust emission observations from the SQUALO project by Traficante et al. (2023). Clump and fragmentation properties from the SQUALO survey that are useful for the analysis carried out in this paper are reported in Appendix A. In particular, we reproduce the observational strategy used to observe the target source HIGALBM24.0116+0.4897 (HG24), on which forthcoming studies from the SQUALO survey are focusing (Jones et al. in prep.). The source is located at a heliocentric distance of 5.2 kpc (Mège et al. 2021). At the SQUALO angular resolution of ~1.4″, we are therefore able to probe the fragmentation process at ~7000 AU scales. We emphasize that our post-processing pipeline is highly flexible and can be adapted to different surveys. In a forthcoming paper we will adapt our pipeline to reproduce the observing strategy of the ALMAGAL survey (Molinari et al. 2025) and discuss clump fragmentation at ~2000 AU scales.

As a first step in the post-processing, we run the previously saved snapshots of our simulations through the radiative transfer code RADMC-3D (Dullemond et al. 2012). This includes the physics of photon emission and absorption, coupled with a raytracing approach: as UV photon packages move through the grid they get absorbed by dust and re-emitted in a new direction, with a new wavelength, both derived according to the prescription by Bjorkman & Wood (2001). To initialize the radiative transfer, we make use of the distribution of sink particles in the 3D-space (volume density cube) and the gas temperature cube ought to derive the dust temperature using a thermal Monte Carlo (MC) computation (Paper II). Once the sources of radiation are identified, the total luminosity is divided into photon packages. We refer to Paper II for a detailed description of the radiative transfer method.

We focus on the radiative post-processed dust continuum emission maps at 1.3 mm as the one displayed in Fig. 2b. The flux is here computed with the same opacity assumption in use for the reference observation, i.e. k_1.3 = 0.005 g cm⁻² (Preibisch et al. 1993), that already includes a gas-to-dust mass ratio of 100. In Paper II, we also post-process the available simulations at Spitzer and Herschel wavelengths (24 µm and from 70 to 500 µm, respectively) to derive an estimate of the clump-averaged bolometric luminosity over mass ratio L/M, that Molinari et al. (2016) demonstrated to be a reliable evolutionary parameter. Specifically, Paper II establishes the quantitative relationship between the timing in SFE from the simulations in Paper I and the observational L/M parameter: $$\log (L/M)={1.20}_{-0.02}^{+0.02}\log (\text{SFE})+{3.28}_{-0.03}^{+0.03}.$$(4)

We use this result in Section 4 and in Section 5 to investigate if and how the global evolution of our clumps – as seen through the eyes of an observer – correlates with fragmentation properties. The conversion from a numerical to an observationally-defined metric of clump evolution allows us to compare the fragmentation properties recovered in Traficante et al. (2023) with the ones from synthetic clumps.

Fig. 2 Example of the post-processing routine on a snapshot of clump collapse. This snapshot, belonging to the Μ1000_μ 10_M7_S2 realization, corresponds to an intermediate level of SFE (~15%). The displayed projection corresponds to the x-axis, i.e., the y - z plane. Left to right: (a) column density map obtained projecting the RAMSES volume density cube. The blue stars mark the projected positions of the sink particles; (b) intensity map after the RADMC-3D computation of the radiative transfer at 1.3 mm; (c) synthetic intensity map after the post processing performed with the CASA software; (d) source extraction and photometry performed by means of the Hyper software. The peaks identification is performed on the synthetic map in (c). Then, source photometry is computed on the PB corrected map, which takes into account the synthetic beam obtained by the combination of synthetic observations with both the 7 m and 12 m arrays. The black crosses and ellipses mark the centroids and the contours of the 8 identified fragments. Sink particles with and without a fragment counterpart are shown as purple and magenta stars, respectively. In each map we show the linear scale corresponding to 0.1 pc, and in the steps involving CASA software we also show the 1.29″×1.51″ synthetic beam footprint.

2.2 Synthetic observations

Actual observations result from the combination of real sky emission with the response of the imaging system, a process that cannot be replicated by the raw simulations. Effects as the smoothing caused by a finite beam size and the filtering of small or large-scale structures caused by interferometric observations, are examples of possible biases that the radiative transfer computation alone does not address. To mimic the features of the SQUALO project ALMA observations, we post-process our data using the CASA software (McMullin et al. 2007), version 6.5.4.

2.2.1 Creating synthetic observations with the CASA software: Main challenges

The CASA software allows the creation of ALMA-like datasets starting from a sky model and incorporate both atmospheric and instrumental effects. For a given simulated source, visibilities³ can be computed using the simalma⁴ task for ideal ALMA interferometer arrays setups. However, during actual observing sessions the configuration of antennas may not match the ideal setup: some antennas might be out of use or poorly calibrated. Together with the elevation of the source, which often does neither reach the zenith, nor its culmination point, this causes the telescope beam to be more elongated than in the ideal case. Additionally, weather conditions might affect the data quality within the required integration time. These factors introduce layers of complexity that must be considered to produce accurate synthetic observations, and can be included by tuning the CASA task simobserve⁵.

2.2.2 Rosetta Stone (SQUALO) pipeline

To address the aforementioned aspects and reproduce the observing conditions in the most accurate way possible, we make use of actual measurement sets and observations reports available for the SQUALO observations. Specifically, the archival data for the 2018.1.00443.S project (ALMA Cycle 6) include a single observing block on the 12 m array (C43-1 configuration), and five observing blocks on the 7 m array for the HG24 reference source. We used the observations reports to recover the following:

• The frequency of the central channel: 231.1 GHz, that is the 1.3 mm wavelength used for the radiative transfer;

• The fiducial antenna configurations, made up of 10 to 12 antennas for the execution blocks on the 7 m array, and 47 antennas for the 12 m array;

• The time on source, with ~470 s of integration for each block on the 7 m array, and ~360 s for the 12 m array;

• The time off source due to calibration necessity, which further splits each observing block into four intervals;

• The observing conditions, with a mean precipitable water vapor (PWV) of ~1.8 and ~2.9 mm computed over the observing blocks in the two different configurations;

• The (average) source elevation, covering a range from ~45° to ~72° in the different observing blocks.

Taking into account the details of each execution block, we used the CASA task simobserve to create synthetic measurement sets of the 1.3 mm sky models produced in Paper II, including the corresponding uv coverage and visibilities. Then, the joint deconvolution with cleaning and masking was performed interactively with the task tclean⁶ with the same pipeline settings that were applied in the imaging of the reference SQUALO source, namely:

• Spectral definition mode: Mfs, to get continuum imaging (Stokes I) with only one output image channel.

• Gridding options: mosaic, to resample visibilities onto a regular uv-grid with azimuthally symmetric beams without side lobes.

• Minor cycle algorithm (deconvolver): multiscale with scales = [0, 7, 21], to look for point-like sources and multiscale (multiples of the beam in pixels, corresponding to the beam and three times the beam) components left in the residual image to add them to the model image.

• Weighting scheme: briggs with robustness parameter of 0.5. We use a flexible weighting scheme based on the signal-to-noise ratio of the MeasurementSets and on the noise threshold. In Traficante et al. (2023) it was possible to independently tune the cleaning noise threshold for the 13 sources, here the threshold is automatically set according to the source evolutionary stage.

• Type of mask for deconvolution: auto-multithresh, to automatically use multiple thresholds for deconvolution.

• Cell size: 0.2″, to sample our 1.29″×1.51″ synthetic beam with ~7 pixels. With this setup, we reach a linear resolution of ~7000 AU (or ~0.036 pc) at the reference heliocentric distance of 5.2 kpc.

Cleaned images with and without primary beam (PB) correction, as the ones displayed in Figs. 2c and d respectively, are produced along with the residual maps and the sources masks. In the PB corrected maps, we recover from 80% to 95% of the flux present in the RADMC-3D projected maps as a function of clump evolution. Specifically, the interferometric filtering of the extended emission has a larger impact toward earlier stages of clump collapse. On the other hand, the flux retrieval is more effective toward advanced stages of clump evolution where most of the emission comes from the compact structures. In Appendix B, we report on further consistency checks performed on the recovered r.m.s. In Appendix C, we address the differences between the results obtained with the Rosetta Stone pipeline and those obtained with simalma.

2.3 Source extraction and photometry

As for the previous post-processing steps, we want the source extraction and photometry to match that performed on the reference observations as much as possible. This approach further reduces the chance of introducing additional biases when comparing observations to simulations. Therefore, the source extraction is performed by using the extraction code Hyper⁷ (Traficante et al. 2015), which went through a fine-tuning procedure by Traficante et al. (2023). Hyper has been originally designed to account for a multi-wavelength source identification and photometry of point-like and compact (elliptical) sources in regions (including crowded regions) where there is significant and variable background emission, and it is thus an optimal tool to use on interferometric observations of clusters formation.

Even though our strategy involves performing the source extraction as on the reference observations, we have also conducted an independent tuning of the software parameters. We have accounted for different detection thresholds, polynomials for background modeling, source sizes and shapes (defined by centroid position, position angle, pixel weighting). In agreement with the settings used in the SQUALO project, we considered as compact sources all identified objects with semi-axes⁸ between one and two times the telescope beam size, having a maximum aspect ratio of 1.5. We modeled the background with polynomials up to the second order, and we performed the peak identification on the PB non-corrected images at three times above the locally estimated r.m.s., to then compute the photometry on the relative PB corrected images, as shown in Fig. 2d.

For each map in our sample of 732, we used Hyper to estimate the fragment and field parameters. Specifically, we derived the number of identified fragments and the coordinates of their centroids. We also determined the semi-minor and semi-major axes, along with the position angle, to characterize the region where the aperture photometry is performed. Within this same region, we measured the peak of emission, and the integrated flux. Additionally, we assessed the local background level for each fragment and the r.m.s. level characteristic of each field.

In Sections 3, 4, and 5.1, we report the results from the realizations with turbulent seed 2, for which we sampled the clump evolution at a higher cadence obtaining a sample of 393 maps. Results from seed 1 are reported in Appendix D. In Section 5.2, we use all the 732 maps to reconstruct the parameter space covered by the whole RS1.0 collection of synthetic observations.

3 Results

Based on the criteria provided in Section 2.3, we are able to identify (and visually validate) 2749 fragments across the 393 maps under analysis. In the following, we report on the physical properties of the identified sample, with specific focus on the mass distribution and the associated uncertainties.

3.1 Physical properties of the fragments

The sample of validated fragments is characterized by radii in the range 0.036 ≤ R_f ≤ 0.071 pc, computed as the geometrical mean of the two FWHMs estimated by Hyper. About 65% of these fragments are symmetrical, of the remaining ∼35%, ∼20% have an aspect ratio between 1 and 1.25, while ~15% have an aspect ratio between 1.25 and 1.5 (the maximum aspect ratio allowed).

The distribution of mass of the fragments was obtained by converting fluxes into masses as follows (Hildebrand 1983): $$M_f=\frac{D^2{S}_{1.3}}{B_{1.3}(T_f){\kappa }_{1.3}},$$(5)

where D is the distance of the reference source, S_1.3 the integrated flux at 1.3 mm estimated by Hyper for each fragment, B_1.3(T_f) is the Planck function computed at 1.3 mm at the chosen temperature T_f for the dust envelope of each fragment, and κ_1.3 the dust absorption coefficient at 1.3 mm.

As in Traficante et al. (2023), we adopt the single temperature model to derive the mass of the fragments associated with the same parent clump. We first assign a temperature to each clump based on the L/M parameter, and then we compute the mass of the fragments therein. Thanks to the Rosetta Stone framework, we are now able to assign a value of L/M to each individual time step on the basis of the associated SFE (see Section 2 and Paper II). The relative mass uncertainties due to the possible error on the temperature estimate is computed via 1000 MC simulations of fragments temperature, which let the temperature vary within fixed ranges. The L/M ranges, the reference temperatures, and the relative temperature ranges are provided in Table 2.

We do not expect the real dust temperature to significantly vary from the adopted one in the initial evolutionary stages. However, according to recent studies (e.g., Motte et al. 2025; Coletta et al. 2025, Jones et al. in prep.), a larger mismatch between the real temperatures and the ranges indicated in Table 2 is expected in the more evolved stages. In this case, we face an additional limitation: the SQUALO clumps reached a maximum L/M of 107 L_⊙/M_⊙, and therefore the L/M range from 10 L_⊙/M_⊙ onward did not require further division into smaller sub-ranges. As a consequence, the reference temperature of 40 K was not conceived to represent clumps with L/M > 100 L_⊙/M_⊙.

Before assessing the clump fragmentation properties, we evaluate how effectively we can detect the faintest objects in each synthetic observation depending on the r.m.s. of the map. As in Traficante et al. (2023), we consider the 1σ value of the local r.m.s. as a lower limit to the peak emission that can be recovered in each map. Assuming that the 1σ value corresponds to the peak flux emitted by a point-like source, we compute the equivalent mass with Eq. (5). From this procedure, we estimate a pointlike source mass sensitivity in the range between ~0.02 and ~5.4M_⊙ across the whole sample. These two-order magnitude span depends on the temperature associated with each clump in a given evolutionary stage (see Table 2) and on the r.m.s. level characteristic of each map, which variations are characterized in Appendix B.

Using the single temperature assumption and Eq. (5), we obtain a mass estimate for each of the 2749 identified fragments. The overall sample spans a range between ~0.02 and ∼1010 M_⊙, as displayed in Fig. 3. For the M500 realizations, the range of mass covers structures from ~0.02 to ~314M_⊙. For the M1000 realizations the range in mass is wider, from ~0.02 to ~1010 M_⊙. We address this unrealistically high value in Section 3.2.1 attributing it to the uncertainties on the clump temperatures and in Section 3.2.2 where we also compare synthetically observed masses with the values obtained from the RAMSES column density maps.

For each mass distribution we compute the median value, which corresponds to ∼4 and ∼11 M_⊙, respectively. Taking into account the span in mass, and the median value of mass of each sample, we find that the initial reservoir has a role in the mass distribution of the object we recover: the most massive fragments are preferentially observed in the most massive clumps. The SQUALO distribution, ranging from 0.4 and 309 M_⊙, falls exactly within the range spanned by both M500 and M1000, as expected. The minimum mass in the SQUALO distribution is one order of magnitude larger than the minimum of the synthetic distribution. This can be a consequence of the range in mass sensitivity associated with the SQUALO sample, from 0.11 to 4.74 M_⊙. While the synthetic clumps are all located at the same distance from the observer, SQUALO clumps are distributed across 13 different distances, resulting in observations at varying resolutions. We acknowledge that variations in clump distances are expected to introduce a slight bias on mass recovery.

We test the robustness of the results on the mass distribution using 1000 MC runs in which a random temperature in the ranges presented in Table 2 is assigned to each fragment. Following this strategy, the ranges of mass are updated: between ~0.01 and ~724 M_⊙, and between ~0.01 and ~2338 M_⊙, for the M500 and M1000 realizations, respectively. Despite fluctuations in the mass range, the median values of the distributions remain stable. The current analysis demonstrates the importance of providing an accurate temperature estimate for our fragments in order to derive reliable mass measurements. For a more detailed discussion on the impact of temperature uncertainties on the mass estimation of our fragments, we refer to the following section.

Fig. 3 Mass distribution. The distributions relative to the M500 and M1000 samples are displayed in purple and orange, respectively. The dashed vertical lines represent the median values of 4 and 11 M⊙. The mass distribution relative to the SQUALO sample is overlaid in blue.

3.2 Temperature uncertainties

Uncertainties on the temperatures of the dust envelope of the fragments – necessary to derive masses from both reference and synthetic observations – contribute significantly to the uncertainties associated with fragment mass estimates. The MC runs allow us to test the assumption of multiple temperatures for multiple fragments within the same parent clump. However, MC simulations are conducted on individual fields, without considering that we are observing subsequent time steps in the evolution of the same clumps. This approach hampers our ability to track potential temperature trends throughout the clump evolution. In the case of synthetic observations, we can address these limitations using the projected temperature maps from the simulations.

3.2.1 Temperature consistency evaluation

Using the known position of the peaks and the area covered by each fragment, we obtain the footprints of the fragments on the RAMSES projected temperature maps. A mean fragment temperature can be then evaluated within these footprints. This kind of analysis, which cannot be reproduced on real observations, is useful to quantify the uncertainties on the assumed temperature values and ranges, along with its impact on mass estimates.

In Fig. 4, we compare the fragment temperatures and ranges adopted in Traficante et al. (2023) with the fragment temperatures derived for the Μ1000_μ 10_M7_S2 realization (the projection in use corresponds to the x-axis, i.e. the y - z plane). The simulations-based recovered temperatures gradually increase with the evolution of the clump, as predicted for the observations-based T_f. However, in absolute values, the latter does not always prove to have a good match with the temperatures computed from the simulations. By fitting the RAMSES mean temperature, computed across ~7000 AU fragments in the same synthetic field, we derive a power-law of: $$T_{f,mean}\text{(K)}\sim 10.4+0.7\times (L/M{)}^{0.75}.$$

This trend suggests that temperature at the fragment scales increases as a function of clump evolution faster than what observed at the clump scales, e.g., T_clump (K) ∝ (L/M)^0.22 for L/M ≥10L_⊙/M_⊙ (Molinari et al. 2016).

As observed in Fig. 4, toward the evolutionary stages below ~30 L_⊙/M_⊙, the fragment temperatures populate only the lower end of the SQUALO temperature ranges. This could be due to the absence of feedback mechanisms, as well as a general overestimation of the fragment temperatures in the SQUALO sample. Clump temperatures in the early stages of collapse may also be strongly influenced by the external radiation field. Between L/M ~30L_⊙/M_⊙ and ~300L_⊙/M_⊙, the prescription well represents both the mean temperature and the temperature range associated with the fragments in the field. Above ~300 L_⊙/M_⊙, few fragments populate the upper end of the SQUALO temperature ranges. The others exceed the predicted maximum temperature (60 K), also reaching temperatures above 100 K. As aforementioned, the reference temperature of 40 K, as well as the 20–60 K span, have not been conceived to represent clumps with L/M > 100 L_⊙/M_⊙, which are outside the SQUALO L/M range. The mismatch between the SQUALO assumption and the temperature estimates derived from the numerical simulation causes a systematical overestimation of the mass of the observed fragments in the most evolved clumps.

The results obtained from the Rosetta Stone pipeline are already noteworthy, particularly in helping observers define an appropriate range of temperatures for their ∼7000 AU resolution surveys. A more detailed investigation of temperature and mass uncertainties arising from observational assumptions will be presented in future works, through the analysis of a new set of simulations that include jets and Hii regions. By incorporating the effects of these feedback mechanisms on local temperatures and, consequently, on the multiplicity and masses of the recovered fragments, we will be able to more realistically model the advanced stages of cluster formation.

Fig. 4 Evolution of fragments temperature (top) and mass (bottom) as a function of the clump L/M for the Μ1000_μ 10_M7_S2 realization. The projection corresponds to the x-axis, i.e. the y - z plane. In the top panel, dark blue stars mark SQUALO reference temperatures according to the evolutionary stage of the clump (see Table 2). Light blue stars mark the reference temperature in the range of L/M not covered by the SQUALO sample. The shaded area corresponds to the temperature ranges used in the 1000 MC simulations of fragments temperature. Small and big orange dots mark, at each given time step, the RAMSES individual temperature of the identified fragments and the RAMSES mean temperature computed across all fragments in the same synthetic field, respectively. The fit of the RAMSES mean temperature as a function of clump L/M is provided as a reference and corresponds to the orange solid line. In the bottom panel, stars and dots indicate mass estimates obtained using the temperatures in the top panel (same color code). Red crosses indicate the mass estimates obtained by integrating the column density maps from the RAMSES simulations within the areas corresponding to the footprints of the identified fragments.

3.2.2 Mass consistency evaluation

In this section we compare the variations in mass estimates caused by the different temperature assumptions with the masses obtained by integrating the column density maps from simulations. This analysis allows us to quantify to what extent the synthetically observed fluxes yield to outcomes that are consistent with the density distribution observed in the numerical models (Paper I). In particular, in Fig. 4 (bottom panel) we compare the total mass in fragments obtained in four different cases: (i) using the SQUALO reference temperatures and temperature ranges; (ii) using the mean temperature across fragments derived from the RAMSES temperature maps; (iii) using the mean temperature of individual fragments derived from the RAMSES temperature maps; (iv) integrating the column densities in the simulations (within the footprints corresponding to the identified fragments).

Where the SQUALO temperature assumptions are comparable to the average temperatures in the 2D temperature maps from simulations (see Fig. 4, top panel), we obtain consistent mass estimates across the four methods. However, this is not always the case, as observed in e.g., the first evolutionary step. In this case, the SQUALO temperature assumption, as the averaged RAMSES temperatures, do not reflect the actual RAMSES temperature along the line of sight causing the underestimation of the total mass in fragments. We observe a factor ∼10 of difference with respect to the mass derived from the column density. Analogous results are discussed in Paper II.

Additionally, the SQUALO-like mass estimates in the latest evolutionary stages exhibit a factor ∼2 of difference with respect to the results obtained using RAMSES temperatures or column density maps. This mismatch occurs because the SQUALO assumptions were designed to cover a narrower range of L/M, ending at ∼100 L_⊙/M_⊙. In the most extreme cases, because of the lower end of the assumed temperature range (see Fig. 4 and Table 2), the most evolved clumps are observed to have a total mass in fragments that exceeds the mass assigned as initial condition in the simulations (see Fig. 6). These results indicate that a finer prescription on the reference temperature (and associated temperature range) should be made to correctly characterize the observed evolved clumps, especially in the case of L/M > 300 L_⊙/M_⊙.

4 Analysis

In this section, we report on the analysis of the statistical fragmentation properties derived from the synthetic observations as a function of the L/M clump evolutionary parameter estimated in Paper II. Specifically, we investigate the level of fragmentation in each synthetic map, the accretion of mass onto the fragments (total mass recovered in fragments with respect to the clump mass and evolutionary stage), and the relative (minimum and maximum) distances between the fragments. Specific focus is on the relative contributions of turbulence and magnetic fields, which cannot be directly studied using the 1.3 mm continuum data from Traficante et al. (2023). We use the variations observed across the three available projections to provide an estimate of the uncertainties associated with the derived mean fragment properties at ∼7000 AU.

Determining the most accurate estimate of fragmentation properties from synthetic maps is a challenging task, as demonstrated by the initial assessment of how well the masses can be recovered using observationally driven assumptions (see Section 3.2.2). A more robust characterization would require identifying three-dimensional, fragment-like structures directly within the simulations smoothed at ∼7000 AU, which is currently beyond the scope of this work. Rather than ensuring 100% consistency between simulated and synthetic quantities recovered at 1.3 mm, our primary goal is to compare the fragmentation properties from the post-processed maps with the SQUALO results within a self-consistent framework (see Section 5.2). To this end, we rely on the observationally driven assumptions presented in Traficante et al. (2023), even if this approach introduces some known biases.

4.1 Fragmentation level

As the clumps collapse and sinks begin to form in the models, we can trace the emergence of fragments in the synthetic observations. In Fig. 5, we present the evolution of the number of identified fragments as a function of the L/M evolutionary parameter. This analysis traces the clump evolution for the available combinations of initial conditions, considering both the mean values and the spread observed across the three different projections. For clarity, the plot is divided into two panels to separate the M500 and M1000 realizations.

The fragmentation levels we observe within individual fields range from 2 to 14 fragments, with an average of ∼7 fragments. Across the 393 maps, we do not find cases in which fragmentation is suppressed. On the contrary, our analysis indicates that fragmentation takes place within the clumps at every evolutionary stage, and that the number of fragments within each field is only weakly dependent on the initial mass of the clump. This suggests that physical factors beyond the initial mass are more relevant in determining the fragmentation patterns.

When examining the distribution of fragment multiplicities for individual realizations, we do not observe a clear, consistent trend with evolution. Instead, we note a general decrease in the number of fragments toward the more advanced phases of clump collapse. This could be due to both pre-existing and newly formed fragments moving closer to the gravitational potential well(s) and merging as the clump evolves. At ∼7000 AU resolution, what appears to be a merger of fragments, is not necessarily a physical merger. It is likely the result of multiple structures that may get closer than our resolution limit (see Section 4.3). From this perspective, the global collapse of the clump appears to dominate the number of detectable fragments at ∼7000 AU resolution, especially in the latest evolutionary stages.

Comparing the level of fragmentation for the three available clump magnetization (i.e., μ = 3, 10, and 100) we observe that, for L/M > 20 L_⊙/M_⊙, a larger number of fragments are preferentially formed when the magnetic parameter μ is set to the highest of the three explored values. However, while magnetized and non-magnetized scenarios can be disentangled within the uncertainties, it is not possible to distinguish between the two magnetized realizations based solely on the number of identified fragments: multiplicity in the μ3 realizations is higher than that in the μ 10 realizations during the later stages of the M500 clump collapse, while it is comparable in the M1000 case.

In addition to the influence of magnetic fields, we investigate the role of turbulence in driving the fragmentation properties of our synthetic clumps. As shown in Fig. 5, the explored levels of turbulence do not appear to have a significant effect on the multiplicity of fragments. Realizations that differ only for the initial value of M, are indeed comparable. It is possible that the values chosen for the Mach number, 7 and 10, are too close to induce any significant variation in the properties of the structures identified at ∼7000 AU.

To further characterize potential relationships between the number of identified fragments and the initial conditions of clump collapse, we also use a suite of out-of-the-grid control models. Specifically, we investigate the impact of the initial volume density profile by setting a Bonnor-Ebert profile (realizations Μ1000_μ3_M7_S2 and Μ1000_μ 10_M 7_S2, respectively), the impact of a lower level of turbulence by setting the initial Mach number equal to 3 (realization Μ1000_μ 10_M3_S2), the impact of stochastic processes by including two additional seeds (realizations Μ1000_μ 10_M7_S3 and Μ1000_μ 10_M7_S4, respectively). Moreover, we evaluate the impact of the numerical prescription for the accretion luminosity coefficient f_acc, used to compute the efficiency of the accretion luminosity by testing different values of this parameter with respect to the reference value of 0.1 (for the Μ1000_μ 10_M7_S1 realization; Paper I). Results from the analysis of these control models are reported in Appendix E.

Fig. 5 Number of fragments identified in each clump as a function of the clump L/M. Left: realizations with initial clump mass of 500 M⊙. Right: realizations with initial clump mass of 1000 M⊙. Shaded areas illustrate the range in number of fragments recovered across the three projections, indicating the uncertainty in the distribution. The color code accounts for the different initial conditions of the realizations, with particular interest on the values assumed by the μ parameter: μ = 3 for the green dots, μ = 10 for the purple dots, and μ = 100 for the orange dots, respectively.

Fig. 6 Total mass accreted onto the fragments as a function of the clump L/M. Left: realizations with initial clump mass of 500 M⊙. Right: realizations with initial clump mass of 1000 M⊙. The temperature of each fragment used to convert the fluxes into total mass (scatter plot) has been assigned using MC simulations, following the ranges in Table 2. The span in mass marked by the shaded areas accounts for stochastic variations due to 1000 different MC runs. The color code accounts for the different initial conditions of the realizations, as described in Fig. 5. The horizontal dashed lines correspond to the initial condition on the clump mass: 500 and 1000 M⊙, respectively.

4.2 Accretion of mass onto the fragments

After presenting the characteristic mass distribution of the fragments across all the available evolutionary stages of the parent clumps (see Fig. 3), we now turn our attention to how the total mass recovered in fragments changes as a function of the L/M evolutionary parameter. This is done by calculating the mass of each fragment, assuming a single-temperature model (see Table 2), and then summing these values for each field. In Fig. 6, we present the mean total mass accreted onto the fragments, which is calculated by considering the three different projections available at each time step. In this case, we use the 1000 MC simulations (see Section 3.2) to indicate the uncertainty on the distribution.

The results reveal that uncertainties in the temperature of the dust envelope surrounding the fragments contribute significantly to the uncertainties in the total mass estimate. The variation in mass due to different temperature assumptions can lead to total mass estimates that are up to ∼5 times higher than the lowest value they can assume. As mentioned in Section 3.2.1 (see also Section 3.2.2), the lower ends of the temperature ranges are responsible for overestimating the fragment masses, particularly in the more evolved stages of clump collapse. As a consequence, in most of the MC runs, the most evolved clumps happen to have a total mass in fragments that exceeds the mass assigned as initial conditions.

By examining the evolution of the total mass in fragments at ∼7000 AU, we find that a continuous accretion process from the parent clumps to the fragments is taking place. The seeds of star formation initially have low masses, often less than 1 M_⊙, and their final masses are progressively built up over time. Even relatively large fragments appear to be accreting mass from the parent clump. This provides clear evidence that the fragments do not evolve in isolation, but are instead fed by material from the parent clump which act as a mass reservoir. It is significant that this trend holds regardless of the different environmental initial conditions explored in the simulations.

Fig. 7 Minimum (top) and maximum (bottom) relative distance between the fragments as a function of the clump L/M. Left: realizations with initial clump mass of 500 M⊙. Right: realizations with initial clump mass of 1000 M⊙. Shaded areas illustrate the range in relative distance recovered across the three projections, indicating the uncertainty in the distribution. The color code accounts for the different initial conditions of the realizations, as described in Fig. 5. The horizontal dashed line corresponds to the resolution limit of the synthetic observations, which is the minimum observable source separation.

4.3 Relative distance between the fragments

Among the temperature-independent quantities of interest, we investigate the relative distances between the fragments as a function of the L/M evolutionary parameter. Specifically, we focus on the minimum 2D distance (d_min,2D) between the fragments, already explored in Traficante et al. (2023), as well as on their maximum 2D distance (d_max,2D). Both quantities are derived from the distance matrix of the fragments, using the centroids of the fragments as edges to compute the relative distances. Together, they provide complementary information about the distribution of the fragments within the clumps at different evolutionary stages, potentially revealing important information about the mechanisms driving the approach or separation of the fragments during clump evolution.

In Fig. 7, we observe the oscillation in time of both d_min,2D and d_max,2D. The oscillation occurs regardless of the initial conditions of the simulations, suggesting that the observed behavior is either a general characteristic of the fragment dynamics at ∼7000 AU scales or the result of a systematic observational bias, rather than being driven by specific starting configurations. The resolved fragments appear to be distributed rather evenly across the parent clumps at each evolutionary stage. However, this result might be strongly biased by the linear resolution of ~0.036 pc, which is the minimum observable source separation determined by the combination of the synthetic beam and the reference distance of our sources. We acknowledge that, in the case of rapid collapse, we would be unable to detect new structures forming at scales below the available resolution, or to track already identified structures that get closer than this minimum separation.

We find that slightly smaller mean values of d_max,2D – reflecting a smaller separation between the most distant fragments – are preferentially observed toward evolved, magnetized clumps. As in the case of the number of fragments (see Section 4.1), the evolutionary stage after which the realizations with different clump magnetization start to behave differently corresponds to L/M ∼ 20 L_⊙/M_⊙. Although we see a difference in the mean values of d_max,2D, this result is not always robust across the different projections. For the µ3 realizations, the fluctuations associated with the three lines of sight are of similar magnitude to the variations introduced by different initial magnetic field conditions, causing its distribution to overlap with that of the non-magnetized scenario at the upper limit.

The different values observed for d_min,2D and d_max,2D may possibly arise from geometrical effects related to the combination between the size of the parent clumps and the number of fragments formed within them. In our case, it is likely that stronger magnetic fields are shaping the interstellar gas leading to more tightly bound configurations of the parent clumps and, consequently, to smaller d_max,2D. However, the observed oscillations and the spread in the data caused by the line of sight effect limit our ability to definitively attribute the recovered fragment distribution to the clump magnetization alone.

5 Discussion

Building on the results and the analysis presented in Section 3 and Section 4, we compare the synthetic observations produced at ∼7000 AU resolution with the output of the simulations presented in Paper I. Furthermore, we compare results from our synthetic observations with actual observations, with particular focus on the fragmentation properties analyzed in the SQUALO project (Traficante et al. 2023, see Table A.1).

5.1 Comparison with numerical models

The systematic production of synthetic observations of clumps undergoing the fragmentation process is a yet to be established approach (e.g., Haworth et al. 2018; Rosen et al. 2020), and to the best our knowledge, there are no existing examples of systematic comparisons of fragmentation properties between the synthetic observations and the numerical simulations from which they are derived. Such comparison serves a dual purpose: first, to assess whether the synthetic observations are reliable proxies for the physical outcomes of the underlying models, and second, to ensure that the numerical computation – native at a numerical resolution of ∼40 AU (Paper I) – accurately captures the key physical processes, even on larger scales.

5.1.1 Fragments versus sinks multiplicity

A substantial discrepancy exists between the multiplicity observed in the simulations and in the corresponding synthetic observations. The number of sink particles identified in the models can reach up to the order of 100 in the later stages of clump evolution, with this value being highly sensitive to the initial conditions of the clump collapse, particularly the clump magnetization (up to ∼50 sinks of difference across different clump magnetization; Paper I). In contrast, the identification of fragments is strongly affected by the resolution limit covering a narrow range of observed fragment multiplicity (2–14), as discussed in Section 4.1, and supported by results from Fontani et al. (2018) and Louvet et al. (2021).

To address the discrepancy in sinks/fragments multiplicity, we make use of the 3D positional information of the sink particles within the simulated volume to then examine the correspondence between the 2D projected positions of the sinks and the locations of the identified fragments in our synthetic observations. This analysis allows us to evaluate the nature of the structures identified at ∼7000 AU and assess the extent to which the fragmentation properties at this scale reflect the underlying physical processes in the simulations. As illustrated in Fig. 2d, we distinguish different scenarios of fragment-sink(s) correspondence, with their respective statistics reported in Table 3. Focusing on the fragments:

• ∼75% of the identified fragments have a correspondence with sink particles. Only ∼17% of the total are associated with a single sink particle;

• ∼25% of the identified fragments do not correspond to any sink particle.

The existence of multiple sinks within individual fragments supports a hierarchical and/or multi-layered fragmentation process, suggesting that fragmentation occurs across multiple scales simultaneously. Larger clumps break into smaller fragments, and within these fragments, further fragmentation can occur, leading to the formation of one or more sinks. The differing sensitivity of fragment and sink multiplicity to the clump magnetization indicates that magnetic fields influence the fragmentation process at different scales, though with varying degrees of significance.

The case of identified fragments that do not correspond to any sink particle, and therefore are not actively forming stars at a given time, can be explained by the following scenarios: (i) they correspond to overdensities that have not yet reached the volume density threshold required to form a sink particle but will reach it in the subsequent steps of clump evolution. (ii) They correspond to overdensities that are transient and will not form a sink particle in the subsequent steps of clump evolution. Some of these will eventually be dispersed. (iii) They appear as overdensities due to projection along the line of sight and/or are the artificial result of the convolution of diffuse structures with the beam of the telescope.

Focusing on the sink particles, we find that ∼76% of them exhibit a corresponding fragment, while the remaining ∼24% do not have a correspondence with the identified fragments. Our analysis reveals that this lack of correspondence never involves the most massive sink particle in each field. It is also noteworthy that the absolute number of sink particles without a matching clump varies with the chosen value of µ. For the µ100, µ10, and µ3 realizations, this goes up to ∼20, 13, and 8 sinks per field, respectively. This result confirms our ability to best recover compact structures within magnetized clumps. Magnetic fields cause initially homogeneous clumps to collapse in dense and compact sheets and/or filaments (e.g., Paper I, Hennebelle & Inutsuka 2019), reducing the presence of diffused matter more efficiently than in the quasi-hydrodynamical case. As a potential consequence, compact fragments that exceed the 3σ threshold might be more effectively formed in the magnetized scenarios.

We will be able to address these specific aspects more effectively by closely tracking the evolution of individual fragments through the use of tracer particles, which we plan to include in a forthcoming suite of simulations. As of now, the number of fragments does not correspond to the number of sink particles due to limitations in both spatial resolution and sensitivity (some sinks may not be bright enough to be detected as fragments arising above the 3σ level). As a result, the increasing number of sink particles⁹ with the progressive evolution of individual clumps and the impact of clump magnetization – discussed in Paper I, and in line with the models by e.g., Peters et al. (2011), Fontani et al. (2018), Hennebelle et al. (2022) – are not equally well observed in the number of fragments at ∼7000 AU. Beside the observational biases, this evidence might suggest that the physical mechanisms explored in the simulations could have different significance at the different physical scales.

Fig. 8 Fragment formation efficiency vs. sink formation efficiency. The color code accounts for the different initial conditions of the realizations, as described in Fig. 5. Specifically, the crosses and the points mark the M500 and M1000 realizations, respectively. Shaded areas illustrate the range in FFE recovered across the three projections, indicating the uncertainty in the distribution. The diagonal dashed line marks the equality between fragment formation efficiency and sink formation efficiency. Data points are expected to lay above the line. Only if the sink particles were resolved, the fragment formation efficiency could align with the one measured for the sinks.

5.1.2 Comparison with sink formation efficiency from the simulations

As previously discussed, the SFE parameter, defined as the total mass recovered in sink particles with respect to the mass of the clump (see Eq. (3)), can be used in simulations as a proxy of the evolution. In a similar manner, we can define a corresponding quantity at the spatial scales probed by the synthetic observations: the fragment formation efficiency (FFE). The FFE is calculated by dividing the total mass recovered in fragments by the mass of their parent clump, offering a way to track the evolution of the clump at intermediate scales.

The comparison between the FFE and the SFE provides a valuable tool for interpreting the distribution of mass into compact structures at various physical scales. By examining how much mass is condensed into ∼7000 AU fragments and sink particles at different stages of clump evolution, we gain first insights into the hierarchical fragmentation process. As shown in Fig. 8, the FFE consistently exceeds the SFE at each evolutionary stage. Fragments at ∼7000 AU are more than the simple product of the convolution of one/more sink particles with the telescope beam. Rather, they are real, compact structures that are made up of gas which is more diffuse in comparison to the sinks. This evidence suggests that a significant fraction of the mass within the clump is condensed into fragments, and only a small fraction of it is ultimately converted into stars.

In particular, in the earliest stages of clump collapse, at SFE ∼ 10⁻³ , the majority of the gas locked into the fragments (FFE ∼ 10⁻²) has not yet been converted into sinks. The difference in these initial efficiencies is further evidence in favor of a multi-layered fragmentation process, where small-scale sinks form within intermediate-scale fragments. With the clump evolution, the mass reservoir available to the fragments diminishes. On the contrary, the sinks accrete with high efficiency from the fragments, bringing the FFE vs. SFE relation closer to the one-to-one line toward the most advanced evolutionary stages. Thanks to the forthcoming RS2.0 suite of simulations, we will explore the same relation in the case of clumps that are formed within filaments. Such analysis will allow us to assess if a depletion in the clump reservoir is observed also when the parent clump can potentially accrete from a larger filament (in opposition to the isolated box).

In actual observations, the trend in accretion efficiency onto compact structures across scales shows a marked decrease, from values up to ≥40% when estimated from clumps down to ∼5000 AU structures, to ≤10% when considering the efficiency from clumps down to the smallest fragments and cores¹⁰, as outlined by Traficante et al., in prep.. We will also test this result thanks to the production of high-resolution (∼2000 AU) synthetic observations.

5.2 Comparison with observations

Synthetic observations serve as a bridge between theory and actual observations, allowing for a self-consistent comparison in terms of common metrics and scales under analysis. The following comparison with the reference observations further proves the predictive power of synthetic observations in the context of star formation. It also represents an important step in the scientific validation of both the RS1.0 simulations and their derived synthetic counterparts.

Previous studies have certainly benefited from direct comparisons between observations and numerical models (e.g., Federrath & Klessen 2012; Peretto et al. 2013; Sanhueza et al. 2019; Añez-López et al. 2020; Avison et al. 2021; Xu et al. 2024), but additional layers of complexity must be taken into account when interpreting actual observations on the basis of numerical outputs. Specifically, observers must first understand what kinds of structures and fragmentation properties the models generate at the scales they are interested in, since these are likely to be different from the observed sinks properties (see Section 5.1).

5.2.1 Fragmentation level

Based on the results from synthetic observations alone, the level of fragmentation observed in clumps at a spatial resolution of ∼7000 AU does not allow us to resolve the degeneracy concerning the physical mechanisms regulating the fragmentation properties. Extreme cases with more than ∼11 fragments appear to be an exception (see Fig. 5), as they are only associated with the μ 100 realizations in our synthetic observations. According to these findings, clumps with significantly high levels of fragmentation, likely originate from clumps in which the magnetic field is sub-dominant and vice-versa, in line with the observations by e.g., Pillai et al. (2015) and Añez-López et al. (2020). Given this result, we hypothesize that the SQUALO clumps, hosting from 1 to 9 fragments, are more likely to correspond to magnetized clumps.

The comparison between fragment multiplicity in the SQUALO clumps and the parameter space covered by the synthetic observations (see Fig. 9) allows us to evaluate how well the combinations of initial conditions explored in the simulations (see Table 1), together with the post-processing strategy, reproduce the fragmentation properties observed in actual interferometric data. We observe that 7 out of the 13 SQUALO clumps are well represented in terms of fragmentation level by the overall distributions. Among the outliers, two clumps, exhibiting 3 and 4 fragments, respectively, are very close to the distribution. Four other clumps, exhibiting 2 fragments or no fragmentation, are poorly represented by the synthetic sample.

Unlike actual observations, our synthetic observations do not reproduce clumps with a very low number of fragments: only a few maps are characterized by the presence of two fragments, and no cases with a single fragment. Since in our data there is an indication that lower values of μ may lead to fewer fragments at a resolution of ∼7000 AU, observations with one or two fragments may require stronger magnetic fields than those considered in this study. Alternatively, they could suggest the involvement of additional physical mechanisms to explain the observed fragmentation. Following the philosophy of the Rosetta Stone approach, we begin investigating this hypothesis by utilizing the aforementioned collection of control models, which lie outside the parameter grid presented in Table 1 (see Appendix E).

By looking at the overall coverage of the parameter space, we observe that higher levels of fragmentation are generally associated with more advanced evolutionary stages. This pattern also holds for the SQUALO sample, where high fragmentation is not observed in clumps with low L/M. Similar findings, supporting the link between high fragmentation and evolved clumps, have been reported in studies based on ∼1500 AU resolution observations, such as those by Coletta et al. (2025) and Elia et al. (submitted). Despite this association, a strict correspondence between the number of fragments and the clump evolutionary stage is not found, either in observations or in synthetic observations. In fact, there are many cases where evolved clumps exhibit low fragmentation levels. This suggests that peculiar interplay of the physical forces and/or the impact of larger-scales dynamics concur to the observed variability in fragmentation levels.

While our synthetic observations predict the formation of up to 14 fragments, higher-resolution observations reveal up to more than 20 fragments (e.g., Beuther et al. 2018; Pouteau et al. 2022; Coletta et al. 2025). This discrepancy suggests that further fragmentation, occurring at smaller scales that our current observations cannot resolve, is expected. As an example, the reference source HG24 exhibits at ∼2000 AU (Coletta et al. 2025) one more fragment with respect to the 2 fragments resolved at ∼7000 AU (Traficante et al. 2023).

Fig. 9 Number of fragments identified in each clump as a function of the clump L/M and comparison with the SQUALO sample. All the combinations of parameters presented in Table 1 are included in the gray area. The stars mark the fragmentation level of the 13 clumps of the SQUALO sample.

5.2.2 Comparison with mass accretion and fragment formation efficiency from the SQUALO sample

Building on prior findings, we now compare the total mass in fragments and the FFE of the SQUALO clumps with the parameter space covered by the synthetic observations. As observed in Section 4.2, the fragments do not evolve in isolation, but remain closely linked to the parent clump throughout evolution, supporting the clump-fed scenario for star formation (e.g., Bonnell et al. 1997, 2001; Klessen & Burkert 2000; Bonnell & Bate 2006; Girichidis et al. 2011, 2012; Hennemann et al. 2012; Vázquez-Semadeni et al. 2019; Peretto et al. 2020; Traficante et al. 2023).

In Fig. 10 (top panel), which shows the total mass accreted onto the fragments as a function of the L/M ratio, only a few outliers (and not too far from the overall covered parameter space) were observed among the 13 SQUALO clumps. The discrepancy likely arises from differences in the initial clump mass between the SQUALO clumps (ranging from 177 to 2689 M_⊙) and the synthetic models, which use two fixed initial mass values. In this case, the normalization by the parent clump mass offers a more direct and meaningful way to compare the two datasets. Indeed, in Fig. 10 (bottom panel), the FFE values from the synthetic models align well with the entire SQUALO sample, suggesting consistent trends across datasets. We acknowledge that the different distances of the SQUALO sources introduce an observational bias in the process of fragment mass estimation (see Eq. (5)), though this bias should not be significant (except for sensitivity issues) when normalizing the mass in fragments by the mass of the parent clump (which is computed as a function of the same distance). Beside the distance correction, the agreement between datasets implies that the parameter space explored in the RS1.0 suite of simulations successfully reproduces the SQUALO sample in terms of FFE, and further supports the agreement with the clump-fed scenario of star formation already proposed for the SQUALO sources by Traficante et al. (2023).

5.2.3 Comparison with relative distances from the SQUALO sample

When comparing the relative distances between fragments from the SQUALO sample with the parameter space covered by the synthetic observations, we observe some discrepancies across datasets regarding both d_min,2D and d_max,2D. The associated trends are shown in Fig. 11. In particular, the d_min,2D distribution in the SQUALO sample shows a decreasing trend over time, which contrasts with our findings (see also the comparison between the SQUALO and CORE surveys presented by Traficante et al. 2023).

In our synthetic observations, we do not observe a decrease in d_min,2D. The large initial spread is due to very high uncertainties (across projections) in the early evolutionary stages of a single realization (see Fig. D.1 for details). These uncertainties cover a shaded area that may give the false impression of a general decreasing trend, even though this behavior does not have a match across all the other realizations. The distribution associated with the synthetic observations cannot fully cover all the SQUALO data due to a known resolution issue. Specifically, we are replicating the observations of the farthest source within the SQUALO sample, which corresponds to the lowest linear resolution in that dataset. As a result, in the comparison, we observe SQUALO points falling below our resolution threshold, reflecting the limitations in resolving smaller-scale structures (which we would resolve by mimicking the high-resolution end of the sample). The resolution issue limits our ability to accurately capture the relative distances between fragments, and to compare these values with those obtained from observations at different resolutions.

The d_max,2D distribution should, in principle, not be directly affected by resolution limits. However, the synthetic observations do not accurately reproduce the spread observed in the SQUALO sample, especially for the most compact sources. In this case, d_max,2D is likely not recovered properly due to the absence of large-scale dynamics in our simulations, as we are simulating clumps within isolated boxes. The larger-scale interactions, such as the global collapse of a filament or the interactions with neighboring structures, can significantly influence the morphology of the clumps and the clustering of fragments therein. Without these effects, the fragment distribution in our synthetic clumps seems to be on average less compact than in the SQUALO sample. The presence of feedback mechanisms might as well affect the observed clustering properties. We plan to address these limitations in the future RS2.0 suite of models by simulating the collapse of non-isolated clumps incorporating feedback mechanisms.

Fig. 10 Top: total mass accreted onto the fragments as a function of the clump L/M and comparison with the SQUALO sample. Bottom: FFE as a function of the clump L/M and comparison with the SQUALO sample. All the combinations of parameters presented in Table 1 are included in the gray area. The stars mark the total mass accreted on fragments for the 13 clumps of the SQUALO sample, and the relative FFE. In the top panel, the total mass in fragments is size-coded according to the SQUALO clump mass. We also provide a reference of the clump mass–marker size relation for the most and least massive SQUALO clumps, of 2689 and 177 M⊙, respectively.

Fig. 11 Top: minimum distance between fragments as a function of the clump L/M and comparison with the SQUALO sample. Bottom: maximum distance between fragments as a function of the clump L/M and comparison with the SQUALO sample. All the combinations of parameters presented in Table 1 are included in the gray area. The stars mark the minimum and maximum distances between fragments for the 13 clumps of the SQUALO sample.

6 Summary and conclusions

In this paper, we have presented the general workflow of the Rosetta Stone project: an end-to-end framework that aims to probe the fragmentation process of massive clumps by comparing observations and dedicated numerical simulations using realistic synthetic observations (refer to Paper I + Paper II + Paper III for the complete RS1.0 framework). Specifically, we have focused on the comparison of the post-processed simulations produced for this work with real ALMA 1.3 mm observations conducted as part of the SQUALO survey by Traficante et al. (2023).

We have demonstrated that, to effectively address the physics governing star-forming sites at the scales probed by actual observations, it is not only necessary to produce readable synthetic observations, but also to ensure that these mimic the observational features (instrumental effects and observational biases) of the dataset they are compared against. This approach has allowed us to perform a self-consistent and systematic comparison between actual observations and synthetic observations to provide a quantitative characterization of how the initial conditions of clumps and environment affect the observed clump fragmentation properties. This comparison was made possible by the estimation of the L/M evolutionary parameter, which was computed for the synthetic clumps in Paper II.

By the analysis of the fragmentation properties of synthetic SQUALO-like clumps performed at ∼7000 AU resolution, we found the following:

• Fragmentation occurs within the clumps at each evolutionary stage. The fragmentation level observed within individual fields ranges from 2 to 14 fragments, with an average of ∼7 fragments, and is comparable to the fragment multiplicity recovered in SQUALO (1–9). As of now, our suite of simulations does not systematically reproduce cases with 1 (zero occurrence) or 2 fragments, which are instead observed in the SQUALO sample;

• At a fixed linear resolution, the number of identified fragments mainly depends on the magnetization of the clump, with high levels of fragmentation preferentially associated with clumps in which the magnetic field is subdominant, and vice-versa. This effect is particularly evident in advanced stages of clump collapse corresponding to L/M ≳ 20 L_⊙/M_⊙. In some instances, stochastic processes limit our ability to discern the impact of clump magnetization on the fragments multiplicity;

• It is possible that magnetic fields are modulating the fragmentation process also in terms of fragmentation patterns, leading to more tightly bound clump configurations and fragment distributions, i.e. to smaller d_max,2D between fragments, in the evolved magnetized clumps;

• The d_min,2D does not exhibit any trend with evolution. However, this result might be strongly biased by the available linear resolution of ~0.036 pc, which is the minimum observable source separation;

• Among identified fragments, ∼75% have a match with one or multiple sink particles. This evidence supports a hierarchical and/or multi-layered fragmentation process;

• Among the identified fragments that do not have a counterpart in sink particles (∼25%): some will develop one (or more) sink(s) with the evolution, some are transients, and some others originate from the convolution of the beam with diffuse gas along the line of sight. In the last scenario, identified fragments do not correspond to real compact structures;

• Our synthetic observations support a clump-fed starformation scenario in which fragments are not isolated from the environment. The accretion of mass onto the fragments continues throughout the whole clump evolution. Based on the comparison with the SFE, the efficiency of the mass accretion process appears to be scale dependent. This represents further evidence of the hierarchical nature of the fragmentation process;

• Overall, the parameter space of fragmentation properties explored with the RS1.0 suite of simulations is in good agreement with the SQUALO observations. The partial mismatch concerning the number of identified fragments suggests that other physical processes with respect to those explored in our grid may be at play in these particular clumps, and has been the primary motivation for developing a set of control models (see Appendix E) outside the grid presented in Table 1.

As anticipated, the upcoming Rosetta Stone catalogs will feature simulations with a broader parameter space, incorporating more physical effects (such as jets and Hii regions) and a more extensive range of scales (spanning from kpc scales to protostellar disks scales). Thanks to its flexibility, the postprocessing pipeline will be also tuned accordingly. We hope that this end-to-end framework will inspire the community to more systematically compare models and observations by means of synthetic observations to advance our understanding of star formation across all scales.

Acknowledgements

We thank the anonymous referee for the useful feedback, which has helped to improve the quality of the paper. This project was funded by the European Research Council via the ERC Synergy Grant “ECOGAL” (project ID 855130). We thank the whole consortium for the stimulating discussions which helped us tremendously through this process. AT gratefully acknowledges support from a mini-grant funded by INAF. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.00443.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. We acknowledge PRACE for awarding us access to the JUWELS supercomputer. This work was also granted access to HPC resources of CINES and CCRT under the allocation A0130407023 made by GENCI (Grand Equipement National de Calcul Intensif). RSK also acknowledges financial support from the German Excellence Strategy via the Heidelberg Cluster “STRUCTURES” (EXC 2181 – 390900948) and from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 50OO2206). RSK thanks the 2024/25 Class of Radcliffe Fellows for highly interesting and stimulating discussions. G.A.F gratefully acknowledges the Deutsche Forschungsgemeinschaft (DFG) for funding through SFB 1601 “Habitats of massive stars across cosmic time” (sub-project B1) and support from the University of Cologne and its Global Faculty programme. D.W. acknowledges support from the Science and Technology Facilities Council (STFC) [grant number ST/Y004108/1].

Appendix A The SQUALO sample

In this Appendix we report the table (see Table A.1) with all the properties of the 13 clumps from the SQUALO survey (Traficante et al. 2023). Specifically, the clump properties have been taken from the latest release of the Hi-GAL clump catalog (Elia et al. 2021) and fragment-scale properties have been taken from Traficante et al. (2023) or derived from the analysis therein.

Appendix B Noise consistency analysis

In order to ensure data consistency between the SQUALO reference observations and the synthetic data, we compare the r.m.s. properties obtained for the two datasets (see Fig. B.1). The estimated r.m.s. for the set of synthetic observations, calculated using the Hyper software directly on the cleaned (primary beam non-corrected) maps, ranges from ~0.02 to ~15 mJy/beam. While its values varies significantly across the sample, it is comparable to the r.m.s. range of the cleaned reference observations, which spans from ~0.8 to ~15 mJy/beam. The mean r.m.s. of 2.9 mJy/beam across the sample is consistent with the 1.42 mJy/beam r.m.s. of the HG24 reference source.

Notably, the r.m.s. values for the synthetic sample increases with evolutionary stage. This effect may be possibly linked to the variation in clump brightness and morphology over time. The integrated signal reflects both the increased flux from compact sources and the ongoing collapse of filamentary structures; indeed, these structures become progressively brighter and less diffuse as the clumps evolve. Despite the ∼3-order-of-magnitude span in r.m.s. values, the corresponding dynamic range, computed as the ratio between the peak flux of the brightest source and the r.m.s. of each field, varies by less than ∼2 orders of magnitude. The dynamic range estimated across all realizations and all evolutionary steps goes from ∼2 to ∼115. This relatively narrow range suggests that the increase in peak source brightness and in the background fluctuations are correlated.

Fig. B.1 Root mean square level obtained from Hyper as a function of the clump L/M and comparison with the SQUALO sample. All the combinations of parameters presented in Table 1 are included in the gray area. The stars mark the r.m.s. level of the 13 clumps of the SQUALO sample computed by Traficante et al. (2023).

Appendix C Pipelines comparison

In Fig. C.1, we present a visual comparison that illustrates the impact of different post-processing pipelines on the resulting synthetic images. Specifically, we compare the results obtained using the Rosetta Stone pipeline - which includes a tailored computation of visibilities and the tclean parameters of the SQUALO project, optimized for clump fragmentation observations - with those from the simalma pipeline implemented in the CASA software. This comparison is based on two distinct evolutionary stages of clump collapse, corresponding to L/M ~ 10 and 250 L_⊙/M_⊙, respectively.

As mentioned in Section 2.2, if not tuned to reproduce a specific observation, the simalma pipeline operates under default conditions such as ideal instrumental and weather parameters, and continuous integration time. The default settings also include general-purpose data reduction and cleaning steps, which are well suited for a wide range of applications. However, for our specific science goals, a more customized approach was necessary to fully reproduce the complexity of our reference observations. Taking as a reference the outputs from the Rosetta Stone (SQUALO) pipeline, the simalma setup causes a loss of flux from the large-scale emission up to ∼80% (∼40%) toward the overall clump structure and up to ∼50% (∼20%) toward the individual fragments in the earliest (latest) evolutionary stage. This second estimation is possible only for the cases in which the same fragment is detected in both maps. It is likely that a general-purpose algorithm fails to recover most large-scale emission from the 7 m visibilities causing, for example, the presence of dips surrounding the main structure in the most evolved clumps (typically indicative of missing short-spacings).

Two main factors contribute to the observed flux loss. First, the synthetic acquisition of visibilities in the simalma pipeline employs a single-track observation strategy, leading to a different UV-plane sampling with respect to the coverage achieved with multiple execution blocks in real observations. The impact of the different sampling also affects the size of the synthetic beam: in the simalma case the beam is ~1.2 times the SQUALO one, while the one obtained with the Rosetta Stone pipeline is comparable to our reference. Second, the cleaning procedure in simalma lacks tools such as auto-masking and multi-scale cleaning, which are typically used to recover large-scale emission components in the final cleaned image when properly tuned. The use of different execution blocks to sample the UV-plane, the incorporation of realistic instrumental setups, the modeling of weather conditions, and the use of joint deconvolution with cleaning and masking techniques implemented by Traficante et al. (2023), make the Rosetta Stone pipeline the optimal pipeline for interpreting the reference observations using synthetic observations as a tool.

Following the production of realistic synthetic maps, the next step in our methodology involves extracting candidate sources from the images and assessing their photometric properties, which allow us to explore the fragmentation properties of the sample (see Section 4). An example of source extraction for the fields in Fig. C.1 is provided to demonstrate how the choice of pipeline also affects source identification. In particular, the different UV-plane coverage and deconvolution strategies significantly affect the quality of the cleaned images, and consequently, the accuracy of source identification. This is particularly evident in the early stages of clump collapse. In the simalma cleaned map at L/M ~ 10 L_⊙/M_⊙, the unrealistic UV coverage leads to the production of false positives and negatives (see Fig. C.1), which make source identification and photometry unreliable unless the source extraction algorithm is carefully tuned. The L/M ~ 10 L_⊙/M_⊙ simalma cleaned image also exhibits a candidate fragment (peaks exceeding three times the local r.m.s.) outside the main emitting structure (gray contours, identified in the model). Two other false candidate fragments are identified outside of this map cut. With the results from this comparison, we further validate the data reduction and the production of final images with the fine-tuned Rosetta Stone pipeline.

Fig. C.1 Post-processing pipeline comparison chart. Top to bottom: two evolutionary phases of the realization with Μ1000_μ 10_M7_S2. The snapshots correspond to ∼1 and 15% SFE, respectively. Left to right: a) synthetic surface brightness map after the radiative transfer; b) synthetic surface brightness map after the post-processing with the simalma pipeline of the CASA software; c) synthetic surface brightness map after the post-processing with the carefully tuned CASA-based Rosetta Stone pipeline, including the SQUALO tclean pipeline. The blue stars mark the positions of the sink particles. The gray contours mark the main emitting area in the RADMC-3D map and are overlaid on the post-processed maps for comparison. The black crosses and ellipses mark the centroids and the contours of the identified fragments. Source identification is performed as described in Section 2.3. In the RADMC-3D map we show the linear scale corresponding to 0.1 pc, while in the following maps we show the relative synthetic beam footprint.

Appendix D Fragmentation properties from S1 realizations

In Fig. D.1, we provide the distribution of the number of fragments identified in each clump, the total mass in fragments and the relative distance between the fragments as a function of the clump L/M for the S1 realizations. In contrast to S2 results, the impact of the clump magnetization on the number of fragments is less clear, especially across the M500 realizations. This discrepancy is mainly due to the peculiar stochastic processes at play when the same models are initialized with a different seed. In the S1 case, the stochasticity has larger impact on the fragmentation, and its effects seem to counteract the magnetic regulation. However, we cannot rule out the possibility that the combination of the stochastic seed and our fixed orthogonal projections may be allowing us to better address the impact of the physical initial conditions of the simulations for the S2 realizations. General considerations on the fragmentation properties reported in Section 4 hold true also for the S1 realizations.

Fig. D.1 Top to bottom: Number of fragments identified in each clump, total mass accreted onto the fragments, minimum and maximum distances between the fragments as a function of the clump L/M. Results from S1 realizations. Left: Realizations with initial clump mass of 500 M⊙. Right: Realizations with initial clump mass of 1000 M⊙. Shaded areas illustrate the span of the analyzed quantities across the three projections, indicating the uncertainty in the respective distributions. Only in the case of the total mass in fragments, the shaded areas represent the stochastic variations of fragments temperature due to 1000 MC runs. The color code accounts for the different initial conditions of the realizations, as described in Fig. 5.

Appendix E The necessity of having an end-to-end framework

Following the philosophy of the Rosetta Stone project, we investigate how to recover clumps with low levels of fragmentation by utilizing a collection of control models that lie outside the parameter grid (see Table 1) used to derive the results presented in the main body of this paper. These different models, shown in Fig. E.1, include two realizations with an initial Bonnor-Ebert (BE) profile, a realization with M = 3, and two realizations in which we explore the impact of different initial seeds, S3 and S4, respectively.

None of these control models allows us to systematically reproduce the cases of low or no fragmentation. Only in the case of the Bonnor-Ebert sphere with initial µ = 3 (top panel), and in the case with M = 3 (central panel), we cover two small portions of the parameter space which are not populated by the RS1.0 models.

The case of the comparison of the different seeds is of particular interest (see Fig. E.1, bottom panel). Variations in the initial seed, which influence the morphology of the evolving clumps, typically induce fluctuations of similar magnitude to those associated with projection effects (i.e., variations of a few fragments). However, in some instances the stochastic processes associated with the choice of the initial seed can completely dominate the fragmentation (as in e.g., Fig. E.1, realization M1000_µ10_M7_S4). Also in less extreme examples, their impact might hinder the chance to disentangle the role of the other physical initial conditions, as assessed in the case of M500_S1 realizations (see Appendix D).

Fig. E.1 Parameter space exploration of the control models. Number of fragments as function of L/M: comparison with the SQUALO sample. All the combinations of parameters presented in Table 1 are included in the gray area. In the top panel, the blue and the purple shaded areas identify two models initialized with a Bonnor-Ebert density profile: M1000_µ3_M7_S 2_BE and M1000_µ10_M7_S 2_BE, respectively. In the central panel, the yellow shaded area identifies the M1000_µ10_M3_S2 realization produced to test the effect on the fragmentation of a lower value of turbulence. It is compared with two models from the RS1.0 suite, i.e., M1000_µ10_M7_S2 and M1000_µ10_M10_S2, in purple and blue, respectively. In the bottom panel, the yellow and the cyan shaded areas identify the M1000_µ10_M7_S3 and M1000_µ10_M7_S4 realizations, respectively, produced to test the effect of of two additional seeds on the fragmentation. They are compared with two models from the RS1.0 suite, i.e., M1000_µ10_M7_S1 and M1000_µ10_M7_S2, in green and purple, respectively.

Beside changing the physical initial conditions, we test the potential impact on the number of observed fragments due to a variation of the numerical prescriptions for the luminosity accretion factor ( f_acc). The efficiency of the accretion affects both the thermal feedback from accretion processes and the total luminosity of each sink. The test is conducted at an intermediate evolutionary step of the M1000_µ100_M7_1 realization. In our suite of simulations RS1.0, we compute the accretion luminosity using a fixed factor of 0.1. To further investigate the impact of this parameter on the number of detected fragments, we explore four additional values of f_acc: 0.01, 0.05, 0.5, and 1, respectively. These values represent a broad range of possible luminosity accretion conditions. As a first test, we change f_acc in the radiative transfer calculations, using the POLARIS¹¹ code (Reissl et al. 2016, 2018) to calculate the temperature and the resulting intensity maps. As the variation of the f_acc parameter is introduced only in the radiative transfer phase rather than in the underlying RMHD simulation, it does not directly change the fragmentation of the clumps and the number of sinks. Rather, it modifies our ability to detect the fragments at a resolution of ∼7000 AU that may correspond to those sinks.

As shown in Fig. E.2, we observe a slight change in the number of fragments as the f_acc parameter is varied. Specifically, higher values of f_acc (i.e., 0.5 and 1) lead to slightly more fragments being detected compared to the lower and intermediate values. This suggests that increased luminosity from accretion luminosity contributes to the identification of ∼7000 AU fragments. We also find that lowering f_acc primarily affects the total luminosity of the low-mass sinks, which derive most of their luminosity from accretion luminosity (in contrast to high-mass sinks which are dominated by internal luminosity). As a consequence, in the CASA post-processed images we still detect fragments associated with high-mass stars, while it becomes more challenging to recover the population of the faintest sources linked to lower-mass sinks. Despite the observed variation in the number of fragments, the accretion luminosity factor does not appear to be the primary factor affecting the recovered fragments multiplicity in this context.

Although the low fragmentation cases are not reproduced within the current simulated parameter space, the outputs from both the control models and the RS1.0 grid help refine the initial conditions and the numerical assumptions for our simulations, guiding future investigations.

Fig. E.2 Number of fragments identified in each clump as a function of facc. Shaded areas illustrate the range in number of fragments recovered across the three projections, indicating the uncertainty in the distribution.

References

All Tables

All Figures

Fig. 1 Flowchart of the Rosetta Stone end-to-end framework: scheme of the comparison between observations and simulations by means of the production of realistic synthetic observations. Real data and inherent information are highlighted in green, synthetic data (RMHD simulations, radiative transfer, synthetic observations) in purple, and the steps of the analysis performed on both real and synthetic maps in orange. In the lighter-shaded boxes we specify how the building blocks of the Rosetta Stone framework have been adapted in this work for the SQUALO project (Traficante et al. 2023) science case.

In the text

Fig. 2 Example of the post-processing routine on a snapshot of clump collapse. This snapshot, belonging to the Μ1000_μ 10_M7_S2 realization, corresponds to an intermediate level of SFE (~15%). The displayed projection corresponds to the x-axis, i.e., the y - z plane. Left to right: (a) column density map obtained projecting the RAMSES volume density cube. The blue stars mark the projected positions of the sink particles; (b) intensity map after the RADMC-3D computation of the radiative transfer at 1.3 mm; (c) synthetic intensity map after the post processing performed with the CASA software; (d) source extraction and photometry performed by means of the Hyper software. The peaks identification is performed on the synthetic map in (c). Then, source photometry is computed on the PB corrected map, which takes into account the synthetic beam obtained by the combination of synthetic observations with both the 7 m and 12 m arrays. The black crosses and ellipses mark the centroids and the contours of the 8 identified fragments. Sink particles with and without a fragment counterpart are shown as purple and magenta stars, respectively. In each map we show the linear scale corresponding to 0.1 pc, and in the steps involving CASA software we also show the 1.29″×1.51″ synthetic beam footprint.

In the text

	Fig. 3 Mass distribution. The distributions relative to the M500 and M1000 samples are displayed in purple and orange, respectively. The dashed vertical lines represent the median values of 4 and 11 M⊙. The mass distribution relative to the SQUALO sample is overlaid in blue.
In the text

Fig. 4 Evolution of fragments temperature (top) and mass (bottom) as a function of the clump L/M for the Μ1000_μ 10_M7_S2 realization. The projection corresponds to the x-axis, i.e. the y - z plane. In the top panel, dark blue stars mark SQUALO reference temperatures according to the evolutionary stage of the clump (see Table 2). Light blue stars mark the reference temperature in the range of L/M not covered by the SQUALO sample. The shaded area corresponds to the temperature ranges used in the 1000 MC simulations of fragments temperature. Small and big orange dots mark, at each given time step, the RAMSES individual temperature of the identified fragments and the RAMSES mean temperature computed across all fragments in the same synthetic field, respectively. The fit of the RAMSES mean temperature as a function of clump L/M is provided as a reference and corresponds to the orange solid line. In the bottom panel, stars and dots indicate mass estimates obtained using the temperatures in the top panel (same color code). Red crosses indicate the mass estimates obtained by integrating the column density maps from the RAMSES simulations within the areas corresponding to the footprints of the identified fragments.

In the text

Fig. 5 Number of fragments identified in each clump as a function of the clump L/M. Left: realizations with initial clump mass of 500 M⊙. Right: realizations with initial clump mass of 1000 M⊙. Shaded areas illustrate the range in number of fragments recovered across the three projections, indicating the uncertainty in the distribution. The color code accounts for the different initial conditions of the realizations, with particular interest on the values assumed by the μ parameter: μ = 3 for the green dots, μ = 10 for the purple dots, and μ = 100 for the orange dots, respectively.

In the text

Fig. 6 Total mass accreted onto the fragments as a function of the clump L/M. Left: realizations with initial clump mass of 500 M⊙. Right: realizations with initial clump mass of 1000 M⊙. The temperature of each fragment used to convert the fluxes into total mass (scatter plot) has been assigned using MC simulations, following the ranges in Table 2. The span in mass marked by the shaded areas accounts for stochastic variations due to 1000 different MC runs. The color code accounts for the different initial conditions of the realizations, as described in Fig. 5. The horizontal dashed lines correspond to the initial condition on the clump mass: 500 and 1000 M⊙, respectively.

In the text

Fig. 7 Minimum (top) and maximum (bottom) relative distance between the fragments as a function of the clump L/M. Left: realizations with initial clump mass of 500 M⊙. Right: realizations with initial clump mass of 1000 M⊙. Shaded areas illustrate the range in relative distance recovered across the three projections, indicating the uncertainty in the distribution. The color code accounts for the different initial conditions of the realizations, as described in Fig. 5. The horizontal dashed line corresponds to the resolution limit of the synthetic observations, which is the minimum observable source separation.

In the text

Fig. 8 Fragment formation efficiency vs. sink formation efficiency. The color code accounts for the different initial conditions of the realizations, as described in Fig. 5. Specifically, the crosses and the points mark the M500 and M1000 realizations, respectively. Shaded areas illustrate the range in FFE recovered across the three projections, indicating the uncertainty in the distribution. The diagonal dashed line marks the equality between fragment formation efficiency and sink formation efficiency. Data points are expected to lay above the line. Only if the sink particles were resolved, the fragment formation efficiency could align with the one measured for the sinks.

In the text

	Fig. 9 Number of fragments identified in each clump as a function of the clump L/M and comparison with the SQUALO sample. All the combinations of parameters presented in Table 1 are included in the gray area. The stars mark the fragmentation level of the 13 clumps of the SQUALO sample.
In the text

Fig. 10 Top: total mass accreted onto the fragments as a function of the clump L/M and comparison with the SQUALO sample. Bottom: FFE as a function of the clump L/M and comparison with the SQUALO sample. All the combinations of parameters presented in Table 1 are included in the gray area. The stars mark the total mass accreted on fragments for the 13 clumps of the SQUALO sample, and the relative FFE. In the top panel, the total mass in fragments is size-coded according to the SQUALO clump mass. We also provide a reference of the clump mass–marker size relation for the most and least massive SQUALO clumps, of 2689 and 177 M⊙, respectively.

In the text

Fig. 11 Top: minimum distance between fragments as a function of the clump L/M and comparison with the SQUALO sample. Bottom: maximum distance between fragments as a function of the clump L/M and comparison with the SQUALO sample. All the combinations of parameters presented in Table 1 are included in the gray area. The stars mark the minimum and maximum distances between fragments for the 13 clumps of the SQUALO sample.

In the text

	Fig. B.1 Root mean square level obtained from Hyper as a function of the clump L/M and comparison with the SQUALO sample. All the combinations of parameters presented in Table 1 are included in the gray area. The stars mark the r.m.s. level of the 13 clumps of the SQUALO sample computed by Traficante et al. (2023).
In the text

Fig. C.1 Post-processing pipeline comparison chart. Top to bottom: two evolutionary phases of the realization with Μ1000_μ 10_M7_S2. The snapshots correspond to ∼1 and 15% SFE, respectively. Left to right: a) synthetic surface brightness map after the radiative transfer; b) synthetic surface brightness map after the post-processing with the simalma pipeline of the CASA software; c) synthetic surface brightness map after the post-processing with the carefully tuned CASA-based Rosetta Stone pipeline, including the SQUALO tclean pipeline. The blue stars mark the positions of the sink particles. The gray contours mark the main emitting area in the RADMC-3D map and are overlaid on the post-processed maps for comparison. The black crosses and ellipses mark the centroids and the contours of the identified fragments. Source identification is performed as described in Section 2.3. In the RADMC-3D map we show the linear scale corresponding to 0.1 pc, while in the following maps we show the relative synthetic beam footprint.

In the text

Fig. D.1 Top to bottom: Number of fragments identified in each clump, total mass accreted onto the fragments, minimum and maximum distances between the fragments as a function of the clump L/M. Results from S1 realizations. Left: Realizations with initial clump mass of 500 M⊙. Right: Realizations with initial clump mass of 1000 M⊙. Shaded areas illustrate the span of the analyzed quantities across the three projections, indicating the uncertainty in the respective distributions. Only in the case of the total mass in fragments, the shaded areas represent the stochastic variations of fragments temperature due to 1000 MC runs. The color code accounts for the different initial conditions of the realizations, as described in Fig. 5.

In the text

Fig. E.1 Parameter space exploration of the control models. Number of fragments as function of L/M: comparison with the SQUALO sample. All the combinations of parameters presented in Table 1 are included in the gray area. In the top panel, the blue and the purple shaded areas identify two models initialized with a Bonnor-Ebert density profile: M1000_µ3_M7_S 2_BE and M1000_µ10_M7_S 2_BE, respectively. In the central panel, the yellow shaded area identifies the M1000_µ10_M3_S2 realization produced to test the effect on the fragmentation of a lower value of turbulence. It is compared with two models from the RS1.0 suite, i.e., M1000_µ10_M7_S2 and M1000_µ10_M10_S2, in purple and blue, respectively. In the bottom panel, the yellow and the cyan shaded areas identify the M1000_µ10_M7_S3 and M1000_µ10_M7_S4 realizations, respectively, produced to test the effect of of two additional seeds on the fragmentation. They are compared with two models from the RS1.0 suite, i.e., M1000_µ10_M7_S1 and M1000_µ10_M7_S2, in green and purple, respectively.

In the text

	Fig. E.2 Number of fragments identified in each clump as a function of facc. Shaded areas illustrate the range in number of fragments recovered across the three projections, indicating the uncertainty in the distribution.
In the text