© The Authors 2025

1 Introduction

The majority of stars, including the Sun (Adams 2010; Gounelle & Meynet 2012; Pfalzner et al. 2015), form within crowded clusters hosting at least one high-mass star (M_⋆ ≥ 8M_⊙, see, e.g., Carpenter 2000; Lada & Lada 2003). Therefore, studies of the physical properties of high-mass star-forming regions (HMSFRs) and, in particular, the mechanisms regulating their evolution, provide crucial information on the birth processes of stars and, ultimately, of planetary systems.

Massive stars influence their surrounding environment in multiple ways (see Krumholz et al. 2014; Rosen & Krumholz 2020), for instance, through gravitational, mechanical (winds, outflows), and radiative (radiation pressure) interaction, and eventually through supernova explosions that enrich the interstellar medium (ISM) with heavy elements. Stellar feedback processes turn out to be relevant not only locally, but also on Galactic scales (Bolatto et al. 2013; Girichidis et al. 2016; Rathjen et al. 2021), and their effects vary as the physical and chemical properties of the involved regions change with evolution (e.g., Caselli 2005; Molinari et al. 2008; Klessen & Glover 2016; Molinari et al. 2019; Coletta et al. 2020).

The formation of high-mass stars (see, e.g., reviews by Stahler et al. 2000; Beuther et al. 2007; Zinnecker & Yorke 2007; Motte et al. 2018a) takes place in cold and dense regions within interstellar molecular clouds (MCs) called clumps, which typically have sizes of the order of 0.1–2 pc, temperatures of ~10–60 K, masses of ~10²–10⁴ M_⊙, and average densities of ~10³–10⁶cm⁻³ (e.g., Elia et al. 2017, 2021). In detail, massive stars form within the densest, relatively hot compact substructures called cores (n ≥ 10⁶ cm⁻³, R ≤ 0.05 pc, e.g., Garay & Lizano 1999; Beuther 2007; Zhang et al. 2009; Sánchez-Monge et al. 2013; Yamamoto 2017; Sanhueza et al. 2019; Sadaghiani et al. 2020; Pouteau et al. 2022), which can be heated by the forming protostar(s) to temperatures up to ~100 K or even more (e.g., Beuther et al. 2005; Kurtz et al. 2000; Cesaroni et al. 2007; Beltrán et al. 2011; Silva et al. 2017).

The formation of a (proto)cluster involves the breaking up of the molecular clump into cores via dynamical fragmentation during its gravitational collapse (Stahler & Palla 2004). The outcome of an individual collapse (e.g., size, mass, and spatial distribution of the formed compact fragments, i.e., cores) depends on the initial conditions of the parent clump, for example, in terms of morphology, mass, and density. Observationally, clumps present a variety of shapes and emission patterns (ranging from extended and diffuse, to filamentary-like, to compact structures) and different degrees of fragmentation into substructures (e.g., Sanhueza et al. 2019; Svoboda et al. 2019). Besides a pure gravitationally driven scenario, other effects have been also brought into play to explain the fragmentation mechanisms leading to the formation of high-mass stars, such as turbulent support (McKee & Tan 2003; Zhang et al. 2009), interstellar magnetic fields (Commerçon et al. 2011; Hennebelle et al. 2011; Zhang et al. 2014; Tang et al. 2019; Palau et al. 2021), radiative feedback (Krumholz et al. 2009; Peters et al. 2010; Hennebelle et al. 2020), and mass flow along filaments (Peretto et al. 2013; Smith et al. 2014, 2016; Lu et al. 2018; Wells et al. 2024).

The competitive accretion scenario (or clump-fed, e.g., Zinnecker 1982; Klessen & Burkert 2000, 2001; Bonnell et al. 2001, 2004; Bonnell & Bate 2006; Smith et al. 2009; Vázquez-Semadeni et al. 2019; Padoan et al. 2020; Traficante et al. 2023; Morii et al. 2024) depicts a dynamic, multi-scale hierarchical framework, where the material that accretes onto the cores (and ultimately onto one or more protostars) is infalling from the larger scale intra-clump medium. Multiple small, low-mass seeds resulting from the fragmentation of a massive clump compete to access and accrete the available reservoir. This implies that not all the fragments will manage to produce high-mass stars. Moreover, objects that succeed in consistently accreting mass will grow more efficiently with time (e.g., Bonnell et al. 2001; Wang et al. 2010; Battersby et al. 2017), thanks to their increasing gravitational attraction on larger portions of the surrounding material. This scenario is nowadays generally favored over the core (monolithic) accretion theory (or core-fed, e.g., McKee & Tan 2003; Tan & McKee 2003; Krumholz et al. 2005; Tan et al. 2014), where the reservoir for accretion is entirely initially available within a massive, isolated core. While in fact substantial proofs of existence of massive prestellar cores are missing, despite dedicated efforts (e.g., Zhang et al. 2009; Wang et al. 2014; Zhang et al. 2015; Sanhueza et al. 2017, 2019; Svoboda et al. 2019; Morii et al. 2023; Mai et al. 2024), multiple evidence in support of clump-fed mechanisms have been produced, both theoretically (e.g., Bonnell & Bate 2006; Bonnell et al. 2007; Smith et al. 2009; Grudić et al. 2022; Hennebelle et al. 2022) and observationally (e.g., Sanhueza et al. 2019; Anderson et al. 2021; Liu et al. 2023; Traficante et al. 2023; Morii et al. 2024; Wells et al. 2024; Xu et al. 2024). These and other works revealed clump-to-core infall motions and core mass growth with evolution.

Clumps and cores are observed at millimeter and submillimeter (mm and submm) wavelengths via the continuum emission of dust grains. The introduction of modern radio mm/submm interferometers, such as the Submillimeter Array (SMA, Ho et al. 2004), the NOrthern Extended Millimeter Array (NOEMA, Chenu et al. 2016), and, especially, the Atacama Large Millimeter/submillimeter Array (ALMA, Wootten & Thompson 2009), has recently made it possible to investigate statistically significant samples of high-mass star-forming regions at an unprecedented level of detail. High-resolution (<1″) observations are required to locate and resolve the compact substructures (cores) within massive clumps, which are typically located at large distances (>1 kpc). Various surveys have recently been conducted to study the fragmentation properties of candidate star-forming clumps at different evolutionary stages.

For example, InfraRed Dark Clouds (IRDCs) were targeted with the SMA (e.g., Zhang et al. 2009; Zhang & Wang 2011; Wang et al. 2011, 2014; Sanhueza et al. 2017). The CORE survey (Beuther et al. 2018) observed 20 protostellar regions at ≲ 1000 au spatial resolution with NOEMA. Among the surveys employing ALMA, ASHES (Sanhueza et al. 2019; Morii et al. 2023, 2024) performed mosaic observations of 39 70 µm-dark regions at ~2000 au resolution, ATOMS (Liu et al. 2020, 2022b,a) studied 146 hyper- and ultra-compact HII regions (HC/UCHIIs) at ~2000 au resolution, while SQUALO (Traficante et al. 2023) inspected a varied sample of 13 regions at different evolutionary stages (from 70 µm-dark to HII) at a ~2000 au resolution. The TEMPO survey (Avison et al. 2023) observed 38 massive star-forming regions with a wide range of evolutionary stages at a ~2000 au resolution. Most recently, ASSEMBLE (Xu et al. 2024) targeted 11 evolved massive clumps at a 2200 au maximum resolution, while the DIHCA survey (Ishihara et al. 2024) observed 30 clumps at a ~900 au resolution. These and other studies (e.g., Palau et al. 2014; Svoboda et al. 2019) have revealed a variety of clump fragmentation degrees, going from fields with only 1 compact fragment to more crowded ones (with up to 40 fragments). Wide ranges of core sizes (from ~600 to a few thousands of au) and masses (from ~0.1 to ~300 M_⊙) were also estimated. The core mass function (CMF) was recently examined in some detail by, for example, Sanhueza et al. (2019), Lu et al. (2020), Sadaghiani et al. (2020), and Pezzuto et al. (2023). In particular, the ALMA-IMF survey (Motte et al. 2022) inspected 15 massive protoclusters, discussing their inner and overall CMF and its variation with evolution (Pouteau et al. 2022, 2023; Louvet et al. 2024).

It must be noted that such analyses are always sensitive to the achieved spatial resolution (see, e.g., Sadaghiani et al. 2020; Louvet et al. 2021), mass sensitivity, and employed observational technique (e.g., single pointing or mosaic). Moreover, importantly, different algorithms are usually used to extract compact sources from continuum maps; for instance, CuTEx (Curvature Thresholding Extractor, Molinari et al. 2016b; Elia et al. 2017, 2021), astrodendro (e.g., Sanhueza et al. 2019; Svoboda et al. 2019; Anderson et al. 2021; Morii et al. 2023), getsf (e.g., Pouteau et al. 2022; Xu et al. 2024), and hyper (Traficante et al. 2023). Each approach inherently defines the properties of objects to be revealed (e.g., in terms of size and shape). This aspect must be considered when comparing results based on different assumptions.

However, a comprehensive understanding of the clump fragmentation process and the mechanisms of massive star formation is still lacking. Large, representative samples observed with high resolution and sensitivity are needed to investigate the physical properties of the cores, their relation with hosting clump properties, and how they evolve over time (e.g., Csengeri et al. 2017; Fontani et al. 2018; Sanhueza et al. 2019; Anderson et al. 2021; Traficante et al. 2023; Xu et al. 2024).

The ALMAGAL survey (see Molinari et al. 2025, Sánchez-Monge et al. 2025), an ALMA Cycle 7 Large Program (2019.1.00195.L, PIs: Sergio Molinari, Cara Battersby, Paul Ho, Peter Schilke), is now providing the opportunity to study all the different aspects of the star formation process in our Galaxy with an unprecedented level of statistical relevance, robustness, and detail. This is thanks to the size and variety of its target sample, featuring 1013 massive dense clumps spread across the Galactic plane and spanning the full evolutionary sequence from IRDCs to HII regions, and its high-resolution and high-sensitivity 1.38 mm observations, making it possible to assess fragmentation down to the typical compact cores scales (≲ 1000 au).

In this paper, we investigate the clump fragmentation process and its outcome by taking full advantage of the advanced capabilities provided by the ALMAGAL survey. The paper is structured as follows. In Sect. 2, we present the full ALMAGAL target sample and the observations carried out with ALMA, also characterizing the main properties of the dust continuum images. In Sect. 3, we describe the compact source extraction procedure, which was developed using the CuTEx algorithm (Molinari et al. 2011) and applied to the maps to reveal the dense substructures (cores) within the whole ALMAGAL clump sample. In Sect. 4, the obtained ALMAGAL catalog of compact sources is presented, along with an analysis of the revealed fragmentation properties and measured photometric parameters of the sources (e.g., fluxes and angular sizes). We derive and discuss the main physical properties of the cores (e.g., physical sizes, masses, and densities) in Sect. 5. The evolution of the CMF is investigated in Sect. 6. Further detailed analysis on the evolution of core masses and its relation with fragmentation is performed in Sect. 7. Lastly, in Sect. 8, we provide a summary of the analyses and the main results of the paper and draw suitable conclusions.

2 Observations

2.1 The ALMAGAL clump sample

The ALMAGAL target sample (see Molinari et al. 2025 for full description) includes 1013 dense clumps with declination δ ≤ 0°, which have been extracted as candidate hosts of high- mass star formation from the large infrared surveys of young massive protostellar regions: Herschel Infrared Galactic Plane Survey (Hi-GAL, Molinari et al. 2010, 2016b; Elia et al. 2017, 2021, 915 targets) and Red MSX Source survey (RMS, Hoare et al. 2005; Urquhart et al. 2007; Lumsden et al. 2013, 98 targets).

The sample covers a broad range of heliocentric distances (~2–8 kpc), masses (~10²–10⁴M_⊙), and surface densities (~0.1–10 g cm⁻², compatible with high-mass star formation, Kauffmann & Pillai 2010; Butler & Tan 2012; Krumholz et al. 2014; Tan et al. 2014). The ALMAGAL clumps probe different Galactic environments, as they are distributed across the disc of the Milky Way, from the central bar to the outskirts of the outer spiral arm, spread over 3–14 kpc in Galactocentric distance. Moreover, target clumps cover the full evolutionary path, representing different evolutionary stages across the star formation process, from IRDCs (e.g., Kauffmann & Pillai 2010; Sanhueza et al. 2012; Barnes et al. 2021) to HII Regions (e.g., Hoare et al. 2005); namely, from prestellar to protostellar phases (see Elia et al. 2017, 2021). This feature is affirmed by the wide range of luminosity over mass ratio covered (~0.05 < L/M < 450 L_⊙/M_⊙), that can be adopted as an evolutionary proxy (e.g., Molinari et al. 2008; Elia et al. 2010; Molinari et al. 2016a; Elia et al. 2017, 2021; Traficante et al. 2023).

The extensive and uniform coverage of all these physical parameters, which are key aspects of high-mass star formation environments, ensures our sample and results have a very high level of statistical significance and robustness that is certainly unprecedented in this field of study.

2.2 ALMA observations

Full details of the ALMAGAL observations can be found in Molinari et al. (2025) and Sánchez-Monge et al. (2025). We provide an overview in the following. Single-pointing observations of the 1013 ALMAGAL targets were conducted from October 2019 to March 2020, and from March 2021 to July 2022 with the ALMA interferometer (Wootten & Thompson 2009) in the spectral Band 6 at around 1.38 mm (~217 GHz, Kerr et al. 2014), using both ACA-7 m (Atacama Compact Array, Iguchi et al. 2009; Kamazaki et al. 2012) and 12m (Escoffier et al. 2007) arrays. The observing time for each target ranged from 1 to 4 minutes. Two distinct observing setups consisting of different configurations of the ALMA arrays were employed: i) near sample targets (d ≲ 4.7 kpc, 535 fields) were observed using the 7m array of ACA (also labeled as 7M throughout the paper) plus the C-2 and C-5 configurations of the 12m array; ii) far sample targets (d ≳ 4.7 kpc, 478 fields) were observed using the 7M array plus the C-3 and C-6 configurations of the 12m array. The C-2 and C-3 (hereafter also referred to as TM2) are the compact configurations of the 12m array (baselines up to 500 m). Instead, C-5 and C-6 (hereafter also referred to as TM1) represent extended configurations of the 12m array (i.e., featuring the longest baselines, up to 2.5 km). The maximum (i.e., highest) angular resolution achieved in Band 6 with the two setups, defined by C-5 and C-6 configurations, is ~0.3″ and ~0.15″, respectively. The maximum recoverable scale (MRS) is set by the minimum baseline length (9 m, given by the 7M array) and is equal to ~29″, corresponding to average physical scales of ~0.5 pc in the near sample and ~0.8 pc in the far sample fields. By combining such observational configurations, we can recover the range of spatial scales of interest for the clump fragmentation process, with the (sub)parsec clump-scale structured and diffuse emission revealed by the 7m array and more compact configurations, and the smaller scale structure (down to the core scales of ~1000 au) traced by more extended arrays (TM2 and especially TM1 configurations). Most importantly, the two setups allow us to obtain a nearly uniform resulting spatial resolution of ~800– 2000 au (median of ~1400 au) for all clumps (see Sect. 2.3), despite their wide spread in distance. This has made it possible to perform a consistent and detailed analysis of fragmentation statistics and physical properties over the whole sample. The main specifications of the observing setups used are summarized in Table 1.

2.3 Continuum maps properties

Data from the different observational configurations employed were combined to produce continuum maps (called 7M+TM2+TM1 images) using the Common Astronomy Software Applications (CASA, McMullin et al. 2007; CASA Team 2022) package. A fully detailed technical description of the ALMAGAL data acquisition and processing pipeline, and characterization of the obtained continuum maps can be found in Sánchez-Monge et al. (2025). Here, we discuss the properties that are most relevant to the interpretation of the analysis and results of this paper.

The left panel of Fig. 1 shows the distribution of the angular size (i.e., circularized full width at half maximum, FWHM) for the synthesized beam of the maps (θ_beam). As anticipated in Sect. 2, due to the different observing setups used, the distribution is clearly bimodal, peaking at ~0.2″ for far sample fields and at ~0.4″ for near sample ones. Overall, the maximum angular resolution, achieved through the addition of the most extended baselines configurations (TM1) to the combined data, is ~0.15″, which leads in our sample to a corresponding maximum spatial resolution ≲ 1000 au, more than sufficient to resolve the typical core scales (<0.05 pc or ~10000 au, e.g., Sánchez-Monge et al. 2013; Beuther et al. 2018; Sanhueza et al. 2019; Sadaghiani et al. 2020; Pouteau et al. 2022). We refer here to circularized sizes, derived as the geometric mean of the major and minor axes of the elliptical map beam (see Sánchez-Monge et al. 2025). Such form is more meaningful for this work since circularized angular sizes have been used to estimate the physical sizes of the extracted compact sources (see Sect. 5.1).

Using those two different setups is justified when converting beam angular sizes into corresponding linear sizes, shown in the right panel of Fig. 1. Although the near and far subsamples do not perfectly overlap, the median values of their distributions (about 1600 au and 1300 au, respectively) are comparable within ~20%, and combine for a fairly homogeneous overall distribution, thus ensuring consistent comparison among targets at considerably different distances. The resulting range of resolved physical scales is within ~800−2000 au in most cases (~90%), with a median of 1400 au. Fields showing very high (<500 au) or low (>3000 au) linear resolutions are outliers (<5%) mainly produced by the refinement of clump distances performed between the time of the proposal preparation and the present work (see Molinari et al. 2025 and Benedettini et al. in prep. for more details). That caused 81 targets (~8% of the sample) which were observed with the near sample setup to assume a new distance compatible with the far sample (thus degrading the corresponding linear resolution) or vice versa (leading to linear sizes well below 1000 au, see Sect. 2.3). As a consequence, in the latter case, the survey resolved even smaller structures than what was originally planned. However, such variations do not affect the statistical significance of the analysis performed in this work.

Figure 2 reports the beam linear size (i.e., circularized diameter) as a function of the target distance. The two increasing trends within near and far subsamples are due to the effect of distance, while the spread is due to the beam size distribution (Fig. 1, left panel). In addition, for visualization and future analysis purposes, we split the sample into two groups according to the achieved spatial resolution, using as discriminator 1500 au, which roughly corresponds to the mode of the beam size distribution (right panel of Fig. 1). Moreover, we distinguish between targets having distances below or above ~4.7 kpc (see Sect. 2.2). Due to the distance refinements described above, these two ranges do not correspond to the near/far sample observational setups for all ALMAGAL fields.

The root mean square (rms) noise level of the maps (σ_rms) was estimated as the standard deviation of the residual image once masked to exclude regions with significantly bright emission in the final intensity image, and thus corresponds to the AGSTDREM keyword stored in the header of the FITS images of ALMAGAL targets (see Sánchez-Monge et al. 2025). These estimates have been used for flux thresholding during the compact source extraction procedure (see Sect. 3). The distribution of rms noise levels across our sample is reported in Fig. 3. Values mostly range between ~0.05 and ~0.2 mJy/beam. The subdivision shown in the plot, identifying three different groups of targets based on their rms level, is introduced and used in Sect. 3.3. Flux sensitivity, intended as the minimum peak flux detected at 1σ level of the map rms noise, is then ~0.05 mJy/beam, satisfying (and actually improving) the requested sensitivity of ~0.1 mJy/beam (see Sánchez-Monge et al. 2025).

Fig. 1 Beam properties of the ALMAGAL continuum maps. Left panel: probability density distribution of the circularized angular sizes of the map beam, separately for near (red) and far (blue) observational samples, respectively (see Sánchez-Monge et al. 2025 for more details). Right panel: distribution of the corresponding linear sizes (i.e., achieved spatial resolution). The grey histogram represents the overall distribution, with the vertical dashed black line marking its median value (1400 au). Near and far subsamples are overplotted with the same color coding of the left panel. It can be noted that, due to the distance variation across our sample, the bimodal distribution in the left panel converts into a fairly uniform beam linear size overall distribution (see text for details).

Fig. 2 Linear size of the map beam as a function of the heliocentric distance of the ALMAGAL clumps. Targets observed with near or far sample configurations appear as red or blue points, respectively. The vertical dashed grey line marks the nominal ~4.7 kpc distance limit separating the near and far subsamples (with the few exceptions described in Sect. 2.3). The horizontal dashed grey line separates higher spatial resolution fields (i.e., beam linear size <1500 au) from lower resolution ones (>1500 au). Note: for better visualization purposes, distance and/or resolution outliers are not shown in this plot (~30 targets, ≲3%).

3 Compact source extraction

3.1 Source extraction strategy

The ALMAGAL survey poses a considerable challenge in terms of adopting the best technique to perform automatic extraction and flux estimation of compact sources across an extensive sample of star-forming clumps in the Galactic plane. The continuum images of the ALMAGAL clumps are characterized by various morphologies, spatial emission distributions, and background properties, the latter being particularly relevant due to the thermal emission of cold dust dominating these environments at mm wavelengths.

Moreover, choosing a particular algorithm inherently selects the nature and properties of objects that will be extracted, for example: size, shape, and compactness. For instance, the widely used code astrodendro (Rosolowsky et al. 2008; see, e.g., Sanhueza et al. 2019; Svoboda et al. 2019; Anderson et al. 2021; Morii et al. 2023; Ishihara et al. 2024) does not make any assumption on the shape of the sources, usually resulting in more irregular structures in which also more extended, diffuse emission might be included. These features make it more suitable for studying the hierarchical structuring of star-forming regions (see Wallace et al., in prep.), rather than performing a physical characterization of the population of compact, coherent cores. Furthermore, astrodendro performs a plain sum of the emission within the structure contour starting from the zero-level flux and does not include by default any background emission estimation and subtraction (Rosolowsky et al. 2008); thus, source fluxes tend to be frequently overestimated. Instead, codes such as getsf (Men’shchikov 2021; see, e.g., Pouteau et al. 2022; Xu et al. 2024), hyper (Traficante et al. 2015; see Traficante et al. 2023), and the here used CuTEx (Molinari et al. 2011, 2017; see, e.g., Molinari et al. 2016b; Elia et al. 2017, 2021, and below), although differing in the flux estimation method, focus on more compact, roundish (elliptical) sources. This latter approach is particularly suited for ALMA continuum images, since they are obtained through a cleaning procedure in which the model image is convolved with an elliptical Gaussian beam (see, e.g., Högbom 1974; Cornwell 2008), so that ultimately small-scale substructures tend to resemble, to some extent, the size and shape of the beam. Such different approaches in extracting substructures make it essential to specify which kind of “compact” sources the detection is aimed at and on what physical scales (given the resolution available in the observations). Namely, it is essential to first clarify what is intended as a core. These aspects have to be properly taken into account when comparing the properties of objects detected under different assumptions. Here, we define as cores the dense, compact sources of roundish (or slightly elliptical) shape with a centrally peaked brightness profile, detected above a certain flux threshold at ~1400 au average scales (see Fig. 1). We fit the sources with a 2D Gaussian, and, to properly recover the true source flux, we simultaneously perform an estimation and subtraction of the background emission (see Sect. 3.2 for full description).

Moreover, to test two different approaches, we injected synthetic sources into selected ALMAGAL maps and conducted source extraction runs with CuTEx and astrodendro. After comparing their performance, CuTEx proved to be more efficient in detecting compact sources and accurately measuring their flux. Additionally, we performed analogous tests with hyper, which is based on aperture photometry. Full details of these tests are given in Appendix A.

Fig. 3 Distribution of the map rms noise level σrms across the ALMA- GAL maps sample (see Sánchez-Monge et al. 2025 for more details). The vertical dashed lines mark the subdivision into the three reported groups of fields that are employed in Sect. 3.3.

Fig. 4 Diagram of the main phases of the source detection and photometry pipeline developed and used in this work. Underlined black steps were already present in the original CuTEx version (but were still further implemented in this work), while those in red have been newly introduced in the context of the present work.

3.2 Source detection and photometry pipeline

The identification and measurement of compact sources were performed on the processed continuum images with the CuTEx algorithm (Molinari et al. 2011, 2016b, 2017). A detailed characterization and operational description of the original version of CuTEx is reported in Molinari et al. (2011). The code has been specifically further developed and improved during this work to best adapt to the properties of the ALMAGAL maps, which partly differ from IR/submm images, like the Hi-GAL maps for which it was originally designed. The key difference is that the beam is sampled with ~6 pixels in ALMAGAL maps (rather than ~3 pixels in Hi-GAL maps). As a consequence, the derivative filter which is applied to the signal image to obtain the “curvature” image (see Molinari et al. 2016b for details) has been modified in this CuTEx version to take into account the beam sampling while preserving the ability to make ~beam-size compact sources stand out, by selectively filtering out the larger spatial scales.

Figure 4 reports the workflow of the final source extraction pipeline that was applied to the ALMAGAL 7M+TM2+TM1 continuum images. The detection procedure uses the nonprimary beam-corrected intensity images, ensuring a homogeneous noise level across the entire field. In more detail, the detection (see Fig. 5, left panel) is carried out on the curvature images obtained with a double differentiation of the intensity image in four different directions (x, y, and the two diagonals). This approach, peculiar with respect to the most commonly used algorithms based on plain thresholding on direct intensity image, allows for an easier identification (using a curvature threshold) and size estimation of compact sources, which stand out in the curvature image, whereas slower varying diffuse and background emission is strongly dampened. Adjacent pixels exceeding a curvature threshold (in our case three times the rms of the curvature image, ${\sigma }_{{\partial }^2}$) are grouped and taken as candidate source(s). The threshold must be exceeded over all four differentiation directions to ensure that possible filamentary-like structures are not extracted as compact objects. The position of the compact source(s) contained within the cluster of pixels is placed at significant local maxima (curvature value at least 1 ${\sigma }_{{\partial }^2}$ above one of all surrounding pixels) on the average of the four curvature images or, if none, at the average coordinates of all pixels in the cluster. The list of sources produced by the detection phase undergoes a first stage of filtering requiring that their peak flux (F_PEAK) is higher than four times the intensity map rms σ_rms. This step provides the advantage of running detection with a relatively relaxed curvature threshold, but then obtaining a more robust catalog with only significant sources, while excluding possible false positives. Moreover, the local value of the curvature is linked to both the source intensity and compactness (size), so that the group of pixels above the threshold could be either sufficiently bright and/or compact. Applying the first stage filtering allows us to keep, among the initially selected cluster of pixels, only the ones that are sufficiently bright to be cleaned during the imaging. In this way, spurious sources produced, for example, by correlated noise and small blobs can be excluded.

Source photometry (see Fig. 5, right panel) is based on the assumption that the source brightness profile can be represented by a 2D elliptical Gaussian with variable parameters, including size and orientation. As anticipated in Sect. 3.1, this approach is suited to detect and measure compact sources in ALMA images, given that the cleaning procedure involves a Gaussian beam. A relatively accurate initial guess for source size parameters (i.e., source FWHM values and position angle, λ_PA, of the Gaussian ellipse) can be obtained using information from the curvature image, namely, by locating the points where the second derivatives reach a minimum that corresponds to a change in the brightness profile curvature, indicating that the source is starting to merge with the underlying background. The eight points resulting from the four curvature images are used to fit an ellipse with variable semi-major and semi-minor axes and orientation, which provides the photometry with the initial guesses for the source size (in terms of the standard deviations σ_a/b) and corresponding orientation, λ_PA, of the Gaussian fit. Additionally, to account for the fact that detected sources may be lying on intense and spatially varying backgrounds, each source is simultaneously fit with a plane of variable inclination degree and direction, to ensure an accurate estimation of the background emission to be then subtracted from the measured brightness of the source. Furthermore, to account for possible partially blended or nearby sources that may contaminate one another, a detected source is simultaneously fit together with its neighbors (through a multiple Gaussians fit) within a given radius, set to 0.75 times the guessed source size. The fit is performed over an image area that is three times the beam width around the source position, thus providing sufficient space to reliably estimate the local background as the median of the pixel fluxes outside the source(s) mask. In particular, the fit is executed on the primary beam- corrected image using a total of nine parameters: six needed to describe the assumed Gaussian profile of each source simultaneously fitted (peak flux F_PEAK, central coordinates x₀ and y₀, standard deviations σ_a and σ_b, and orientation λ_PA) and three for the planar background plateau. To permit its refinement on the signal image, the source size estimated by the fit is left free to vary within 0.5–1.5 times the initial guess, assuming a defined lower limit of 0.95 times the map beam size (in terms of the minimum axis of the Gaussian ellipse, ${\theta }_{\text{beam}}^{\min }$).

An extensive and varied series of tests were performed on both observed and synthetic maps, generated by injecting Gaussian sources at different flux levels into ALMAGAL images to reproduce in a controlled way the conditions addressed in the “real” source extraction run (see Appendices A and B for details). These testing and tuning runs were performed to maximize source detection efficiency and photometry accuracy by finding the most suitable set of extraction parameters, able to account for the diverse physical properties of the targets (e.g., source emission and background structure) and the technical features of the maps (e.g., noise and dynamic range properties). For each of the test runs, the outcome of the source extraction procedure across the whole clump sample was evaluated through a thorough visual inspection. The integrated flux F_INT of the detected source is computed as (Molinari et al. 2011) $$F_{\text{INT}}=2\pi {\sigma }_{\text{a}}{\sigma }_{\text{b}}{F}_{\text{PEAK}}/{\text{Ω}}_{\text{beam}},$$(1)

where F_PEAK is the source peak flux subtracted by the background at the source location (expressed in Jy/beam), σ_a/b are the standard deviations of the fitted source profile, and Ω_beam is the beam solid angle: $${\text{Ω}}_{\text{beam\hspace{0.17em}}}=\frac{2\pi }{8\ln 2}{\theta }_{\text{beam\hspace{0.17em}}}^2,$$(2)

where θ_beam is the circularized FWHM of the beam, computed as the geometric mean of its major and minor axes: $${\theta }_{\text{beam\hspace{0.17em}}}=\sqrt{{\theta }_{\text{beam\hspace{0.17em}}}^{\text{maj}}\cdot {\theta }_{\text{beam\hspace{0.17em}}}^{\min }}.$$(3)

Lastly, a second and final filtering stage, which we specifically introduced in this improved version of CuTEx, is applied on the source catalog, requiring that the background-subtracted peak flux from photometry is higher than five times the map rms σ_rms (which is computed on the non-primary beam-corrected image), once σ_rms is also corrected for the local value of the primary beam. We imposed this rather conservative S/N threshold to obtain a more robust and reliable core catalog, which still provides very large statistics. Based on thorough visual inspection, lower S/N would lead to a significant number of spurious detections. The output of the procedure for the whole target sample, resulting in the catalog of ALMAGAL compact sources, is presented in detail in Sect. 4.

Figure 6 shows the outcome of the source extraction procedure for selected ALMAGAL fields, representing different fragmentation and emission patterns. We note that besides the sources that were detected, other relatively bright and centrally peaked structures might appear visually. While such objects seemingly did not satisfy the S/N and/or compactness requirements, they are in most cases recovered when allowing a lower S/N (e.g., 3 σ_rms, as shown in Fig. 6).

Fig. 5 Example of the source extraction procedure performed with CuTEx (see Fig. 4) for the ALMAGAL field AG349.6440-1.0948. Left panel: detection phase performed on curvature images along the reported four different directions; the blue crosses mark the candidate compact sources. Right panel: outcome of the photometry procedure performed on the corresponding intensity image; the blue ellipses represent the estimated shape of the compact sources (in terms of FWHMs and position angle, λPA) as a result of the 2D Gaussian fitting procedure (see Sect. 3.2 for details). Note: the color scale of the image was adjusted to highlight all the three sources at the same time, so their visual appearance in terms of emission extent might not correspond to their actual size (ellipses). The map beam is shown in red in the bottom-left corner of the plot. The angular and physical scales are also reported.

3.3 Flux completeness and photometric accuracy

Flux completeness levels for our compact source catalog were estimated through synthetic source extraction experiments. As a preliminary step, we visually divided the map sample in three groups according to their rms noise properties, as shown in Fig. 3: (i) low noise maps (σ_rms < 0.06 mJy/beam, 53 fields or 5% of the sample, hereafter Group 1), (ii) medium noise maps (0.06 ≤ σ_rms ≤ 0.1 mJy/beam, including the peak of the distribution, 769 fields or 76%, hereafter Group 2), and (iii) high noise maps (σ_rms > 0.1 mJy/beam, 191 fields or 19%, hereafter Group 3). This was done in view of performing a dedicated analysis for each of the groups to obtain more accurate flux completeness estimates. We then injected synthetic 2D elliptical Gaussian sources (20 per field) at multiple (from 10 to 15) integrated flux levels into 60 selected ALMAGAL fields representing a variety of emission morphologies, which we divided into three groups of 20 according to the rms noise classification presented above. In total, we generated ~4000 synthetic sources for each group. We ran our extracting algorithm on the whole set of generated images using the same parameter setup applied on native fields (see Sect. 3.2) and compared the photometry output obtained for the simulated sources with their “truth tables”, annotating the number of detections recorded for each generated flux level $\left(F_{\text{INT}}^{\text{inj}}\right)$. Detailed description and motivations for this methodology are provided in Appendix B. The results are shown in Fig. 7 for the three different groups of fields. The flux completeness level $\left(F_{\text{INT}}^{\text{compl}}\right)$ is defined as the integrated flux value at which a 90% rate of detection is reached. Completeness flux is 0.58 mJy for Group 1 (low noise fields, left panel), 0.94 mJy for Group 2 (medium noise fields, mid panel), and 2.05 mJy for Group 3 (high noise fields, right panel), respectively. A summary of the flux completeness analysis is given in Table 2. Corresponding mass completeness is addressed in Sect. 5.2.4.

The flux completeness estimation procedure also enables us to characterize the accuracy of our source extraction method in terms of flux and size recovery. Indeed, for injected synthetic sources the estimated integrated flux is on average within ~20% from the generated one, while the source size is within ~20%, both consistent with the true values given the adopted uncertainties. Dedicated plots are shown in Appendix B.

Fig. 6 Outcome of the source extraction procedure performed with CuTEx on selected ALMAGAL clumps (from left to right, fields AG316.4045- 0.4657, AG012.4977-0.2233, AG013.2428-0.0854). As in Fig. 5, the blue ellipses mark the compact sources detected at 5 σrms. The three clumps showcase the variety of fragmentation degree, morphologies of structures and patterns of emission revealed across the ALMAGAL sample. The CuTEx algorithm was optimized and tested in order to efficiently perform in all these different conditions. The green ellipses mark the sources that are recovered by allowing a signal-to-noise ratio (S/N) cut of 3 σrms. The map beam is reported in red in the bottom-left corner of each panel.

Fig. 7 Completeness levels (i.e., percentage of detection) as a function of the “true” integrated flux of the synthetic sources injected into the maps $\left(F_{\text{INT}}^{\text{inj}}\right)$ for the three different groups of targets selected based on their rms noise, as introduced in Sect. 2.3: Group 1 (left panel), 2 (middle panel), and 3 (right panel), respectively. The two dashed red lines in each panel identify the estimated completeness flux, $F_{\text{INT}}^{\text{compl}}$, which is reported in red.

4 Compact source catalog analysis

4.1 ALMAGAL catalog of compact sources and their photometric properties

Table 3 reports the catalog and the main photometric properties of the compact sources (cores) extracted from the ALMAGAL continuum images with the pipeline described in Sect. 3.2 (the full table is available at the CDS in machine-readable form). The cores are progressively ordered within the hosting clump, thus indicating the number of detections reported in each target, which is discussed in Sect. 4.2. All fluxes have been computed on the primary beam-corrected image. In particular, the background-subtracted source peak flux is reported in Col. 8. The background flux F_BKG can be intended as a tentative proxy of the local conditions where the compact source has been detected and measured. For example, a relatively high F_BKG may suggest the potential presence of a diffuse emission plateau or other nearby bright sources under or in close proximity to the detection location on the map, so that its flux estimate may be considered to be less robust. A few F_BKG values (<3%) turn out to be negative. These values could, in principle, affect the goodness of the source flux estimate, but this is not an issue in our catalog since we verified that they are not significant as they stay within the local fluctuations on the map σ_rms. The parameters FWHM_a and FWHM_b refer to the non-beam-deconvolved full sizes of the axes of the fitted Gaussian ellipse (rotated with respect to the x-y coordinate system of the map by the position angle λ_PA). Once converted to σ_a and σ_b, respectively, they are used in Eq. (1) to compute the source integrated flux (F_INT). The main properties of the corresponding map elliptical beam (major and minor axes extent, and orientation) are also reported to give an indication of the size of the source relative to the beam. The distribution of this ratio is shown and discussed in Sect. 4.3.2. As detailed in Sect. 3.2, detected sources were selected to have a minimum FWHM above 0.95 times the beam minor axis, so that it may happen that in some cases ${\theta }_{\text{beam\hspace{0.17em}}}^{\min }$ turns out to be greater than the source minimum extent between FWHM_a and FWHM_b. All included core peak fluxes are above the 5 σ_ms cut applied in the last photometric filtering stage. After a thorough visual inspection of the detection maps, however, five fields (namely AG025.8252-0.1777, AG029.9558-0.0161, AG274.0659-1.1488, AG305.3676+0.2128, AG331.5117-0.1024) were found to be likely slightly affected by some artifacts coming from the continuum map cleaning process (see Sánchez-Monge et al. 2025). These fields then required a dedicated tuning of the extraction parameters in order to optimize the robustness of the compact source extraction. As a consequence, the cores detected within these clumps present fluxes that were selected to be above 6 σ_rms (for the four targets AG025.8252-0.1777, AG029.9558- 0.0161, AG274.0659-1.1488, and AG331.5117-0.1024) or 8 σ_ms (for target AG305.3676+0.2128).

4.2 Detection and fragmentation statistics

The ALMAGAL compact source catalog features 6348 cores detected in 844 clumps from the total 1013 inspected (∼83%). Therefore, 169 fields (∼17%) reported no detections. Such evidence does not imply the lack of signal in the continuum maps, but rather that such emission did not satisfy the requirements of our source extraction procedure, for example, in terms of morphology, compactness, and contrast. This may be related to the physical properties of those targets (e.g., surface density, evolutionary stage, IR fluxes) as investigated in detail in other ALMAGAL works (Molinari et al. 2025, Elia et al., in prep.).

Figure 8 (left panel) shows the distribution of the number of cores (N_frag) revealed across the inspected clumps. The overall N_frag range is 1–49, while the average number of detections per clump is 8 and the median is 5, respectively. Moreover, 130 clumps (15%) showed only 1 detection, while 58 clumps (7%) appear to be rather crowded, reporting more than 20 detected cores. Table 4 summarizes the overall detection statistics. The right panel of Fig. 8 shows the distribution of the spatial densities of cores within the hosting clumps (computed within the entire FOV), which ranges between ~0.5 and ~500 sources pc⁻², with a median value of ∼10 pc⁻² .

The large statistics and broad range of physical parameters provided by ALMAGAL reveal the fragmentation properties of a wide variety of high-mass star-forming regions, which may be characterized by varied morphologies, physical properties and evolutionary stages. This most probably results in the large range of fragmentation degrees that we report in our catalog, much wider than that found by other recent ALMA studies, which feature much smaller samples. Sanhueza et al. (2019) (ASHES) detected ∼300 cores in 12 massive 70 µm-dark prestellar clump candidates (>500 M_⊙) observed with mosaics at ~4800 au average resolution, with a number of fragments per clump between 13–39. Later, Morii et al. (2023) extended the sample to 39 clumps, revealing ∼840 cores in total, with a range of 8–39 fragments per clump (median of 20). Svoboda et al. (2019), instead, found 67 cores in 12 massive (>400 M_⊙) starless clump candidates with a fragmentation degree ranging from 1 to 11, observed at ∼3000 au spatial resolution. In a sample of 146 late-stage clumps (HII regions), Liu et al. (2022a) (ATOMS) detected a total of 450 cores (candidate hot cores/HC-UCHIIs). Traficante et al. (2023) (SQUALO) inspected a representative sample of 13 massive clumps with masses above 170 M_⊙ at various evolutionary stages (0.1 ≤ L/M ≤ 100L_⊙/M_⊙, covering from prestellar to advanced phases) with a maximum resolution of ∼2000 au, within which they identified 55 dense fragments (ranging from 1 to 9 fragments per clump). Avison et al. (2023) (TEMPO) found 287 cores at ∼2000 au resolution in their 38 HMSFRs spanning ~0.5 ≤ L/M ≤ 700L_⊙/M_⊙. Most recently, Xu et al. (2024) (ASSEMBLE) detected 248 cores in their 11 evolved massive clumps. Furthermore, Beuther et al. (2018) (CORE) observed 20 clumps with distances below 5.5 kpc and masses above 40M_⊙ at a similar ∼1000 au maximum resolution with NOEMA, revealing diverse fragmentation properties, going from fields with a single fragment to crowded fields with up to 20 detections.

Therefore, overall, in our clump sample we find a wider range of fragmentation degrees with respect to other similar recent high-resolution studies. However, the statistical and physical significance (in terms of number and range of physical properties of the observed clumps) and spatial resolution achieved by these and other interferometric works are considerably lower than the ones provided by the ALMAGAL sample. The capability to observe highly fragmented clumps is indeed expected, thanks to the high resolution and high sensitivity of the ALMAGAL observations. Nevertheless, a significant difference with previous studies is the high fraction of clumps revealing a small number of fragments (28% with 1 or 2 fragments). We underline, in this respect, that we focused our source extraction strategy on detecting only compact substructures. However, this might be also related to the properties of the clump, for instance, in terms of mass, surface density, and evolutionary stage. The relationships between the observed fragmentation and the physical properties of the clumps and their embedded cores are investigated in Sect. 7, while a more extensive analysis is given in Elia et al. (in prep.).

Ultimately, in more general terms, it has to be taken into account that the capability of detecting dense fragments always depends to some extent on several aspects, such as modality (mosaic or single-pointing), sensitivity, and resolution of the observations, design of the source extraction algorithm, and uncertainties on source size and flux estimates (see, e.g., Sadaghiani et al. 2020; Louvet et al. 2021), in addition to the inherent physical properties of the targets (e.g., morphology, spatial emission distribution, evolutionary stage; see Molinari et al. 2025). In Sect. 4.2.1, we analyze the potential impact on the revealed fragmentation properties of the different resolutions and distances featured in our data.

Fig. 8 Fragmentation statistics deduced from the compact source catalog (Table 3). Left panel: overall distribution of the number of cores detected in each clump; detailed statistics are reported in Table 4. Right panel: overall distribution of the derived spatial density of the cores in each clump, computed within the entire FOV. The vertical dashed black line marks the median value of ∼10 pc−2 . In both panels, only clumps revealing at least 1 detection are included.

Fig. 9 Fragmentation statistics evaluated for near sample and far sample sources separately. Left panel: probability density distribution of the number of cores detected in each clump for near sample (red bars) and far sample (blue bars) sources, respectively. Right panel: same as left panel, but only including clumps containing core masses above 0.23 M⊙ (see text for interpretation). Detailed statistics for the individual distributions are reported in Table 5.

4.2.1 Sensitivity and physical effects on the observed fragmentation statistics

The angular resolution and target distance are critical parameters for the interpretation of our data, given the broad range of clump distances covered and the different specific observational setups used for near and far sample sources (see Sect. 2).

Figure 9 (left panel) shows the distributions of the probability density of the number of fragments within each clump, for the two near and far subsamples separately (i.e., observed with the two different observational setups). The plot displays clear discrepancies between the two distributions. Detailed statistics for each distribution are reported in Table 5. The total number of fragments detected in near sample targets is more than twice that found in far sample ones, the gap being not merely justified by the different number of clumps included in each subsample. As for the probability density distribution, although the overall range of N_frag is nearly the same, the near sample distribution is clearly more shifted towards higher N_frag values with respect to the far sample one. The latter is indeed strongly peaked towards N_frag below 5 and in particular, importantly, reports nearly twice the number of fields with only 1 detection. This trend is also confirmed by the different N_frag mean and median values of the two distributions (considerably higher for the near sample) and the percentage distribution of counts among different ranges of N_frag . Since the angular FOV of the maps is the same (∼35″) for both near and far sample targets, due to the different distances we are actually sampling a wider spatial scale in far sample fields (on average ~1 pc compared to ~0.6 pc of near sample fields). This could in principle result in a higher number of detections in far sample fields (assuming a uniform distribution of fragments in the sky). The fact that this expected trend is not seen by comparing the two distributions tells us that another effect is actually overpowering resolution in dominating the fragmentation properties that we observe as a function of distance, which we identify with the change in mass sensitivity across our sample. As it is directly shown in Sect. 5.2.4, the constant sensitivity in flux of the sample translates into an increasing mass sensitivity threshold with distance, which gradually reduces the possibility to reveal lower-mass cores going from the nearest to the farthest clumps. From the fragmentation viewpoint, the wider range of masses available for detection naturally leads to a higher mean number of detections within near sample targets with respect to the far sample. Indeed, we verified (right panel of Fig. 9 and Table 5) that near and far sample distributions of N_frag become almost identical when selecting only fields hosting derived core masses above the completeness limit (∼0.23 M_⊙, see Sect. 5.2), thus proving that mass sensitivity is indeed the predominant factor in the observed fragmentation statistics.

The left panel of Fig. 10 reports the two subsample distributions of spatial density of fragments. Again, a clear shift between the two profiles is found, largely inherited from the discrepancy in the number of fragments seen in Fig. 9 (left panel), with near sample counts oriented towards higher densities than far sample ones, as testified by the corresponding median values (21 pc⁻² for near sample sources, 5 pc⁻² for far sample sources). This behavior is for the most part ascribable to the same mass sensitivity effect discussed above, and the balance between the two subsamples is indeed nearly recovered when considering core masses above completeness (right panel of Fig. 10). Evidently, most of the counts of the near sample distribution that caused the considerable shift with the far sample in the overall distributions are now lost, proving that they mostly pertained to low-mass sources, which we are actually not able to reveal at larger distances. However, a slight shift between the two subsamples (quantified by the different median values) remains also when filtering the detections by mass. Therefore, the influence of an additional physical effect on this trend cannot be ruled out. If one assumes the fragments to be spatially distributed in a rather uniform fashion in the clump, one would expect to find constant spatial densities between near and far sample maps, even if they cover a different range of spatial scales. The fact that we observe a discrepancy, instead, could suggest that fragments are preferentially concentrated within the central, relatively spatially limited regions of the clumps, so that we find lower densities when expanding the sampled FOV.

The separation and spatial distribution of the fragments revealed within the ALMAGAL clumps will be investigated in detail in a following work (Schisano et al., in prep.), to make further considerations on the possible physical mechanisms governing the fragmentation process (see, e.g., Kainulainen et al. 2017; Sanhueza et al. 2019; Svoboda et al. 2019; Gieser et al. 2023a; Morii et al. 2023; Traficante et al. 2023; Ishihara et al. 2024).

Fig. 10 Analysis of the derived spatial density of the cores. Left panel: distribution of the spatial density of the cores within each clump for near sample (red bars) and far sample (blue bars) sources, respectively. Right panel: same as left panel, but only including clumps containing core masses above 0.23 M⊙. For comparison, the overall distributions of the left panel appear here in transparency in the background as pink (near sample) and cyan (far sample) bars, respectively. In both panels, overlapping areas between the two distributions are marked by the magenta bars. The median values of the near and far sample distributions are also reported in each panel.

4.3 Photometric properties

4.3.1 Integrated fluxes

The integrated fluxes of the extracted compact sources, as estimated through the CuTEx photometry procedure described in Sect. 3.2 by assuming a Gaussian source profile, are reported in Col. 12 of Table 3. The overall distribution of integrated fluxes is shown in Fig. 11 (left panel). The dashed cyan line indicates the flux completeness limit of 0.94 mJy, estimated for our catalog as described in Sect. 3.3. We decided to adopt the completeness limit computed for Group 2 targets as overall reference, since such group represents the vast majority of the ALMAGAL maps (76%, Fig. 3), characterized by a “normal” level of noise. The measured integrated fluxes range from ∼0.3 to ∼10³ mJy. The median of the integrated flux distribution is 2 mJy. The faintest source detected is ∼2 times above the requested flux sensitivity of ∼0.1 mJy, corresponding to a S/N of ∼5.5. The brightest source has a S/N of ∼1500 (see Sánchez-Monge et al. 2025). The uncertainty of estimated integrated fluxes is ∼20%, as deduced by comparing the fluxes of synthetic sources measured with CuTEx with the true, generated ones at different flux levels (see Sect. 3.3 and Appendix B for details).

The right panel of Fig. 11 reports the measured integrated fluxes of the sources as a function of the distance of their hosting clumps. The trend traced by the blue line shows that the median measured integrated flux is roughly constant (i.e., within a factor of ∼2) with the source distance. This was indeed designed for the ALMAGAL survey in order to ensure that consistent and comparable fluxes can be recovered across the whole clumps sample despite its wide spread in distance (see Molinari et al. 2025). A similar trend is found for 15th and 85th percentiles of the flux distribution within bins of distance, indicating that detected objects have roughly the same statistical distribution in flux independently from their distance. Moreover, the range of fluxes recovered at the different distances is roughly the same. The measured integrated fluxes have been used to estimate the mass of the extracted cores, as described in Sect. 5.2.

Fig. 11 Integrated fluxes of the extracted compact sources as estimated with CuTEx. Left panel: overall distribution of the integrated fluxes of the sources. The adopted flux completeness value of 0.94 mJy (see Sect. 3.3) is marked by the vertical dashed cyan line. The vertical dashed black line marks the median value of the distribution (2 mJy). Right panel: measured integrated fluxes of the sources as a function of the hosting clump distance. The solid blue line traces the median FINT values within given bins of distance. Light blue lines correspond to the 15th (bottom) and 85th (top) percentiles of the flux distribution in the same bins of distance. Note: for improved visualization, distance outliers (<3%) are not shown in this plot.

Fig. 12 Angular sizes of the extracted sources. Left panel: probability density distribution of the circularized FWHM of the sources. Red and blue bars include sources detected within the near sample and the far sample targets, respectively (see Sect. 2). Right panel: distribution of the ratio between the estimated angular size of the source and the map beam size (both in circularized form), providing a resolution-independent estimate of the source sizes. The vertical dashed black line marks the median value of 1.15 for the ratio.

4.3.2 Angular sizes

The circularized angular size of the extracted sources, FWHM_circ, was computed as the geometric mean of the FWHM of the fitted source ellipse reported in Cols. 9–10 of Table 3. The overall distribution of angular sizes ranges between ∼0.15–1.4″. The two different observing setups adopted in the ALMAGAL survey (described in Sects. 2 and 2.3) translate into two distinct distributions of angular sizes for detected sources belonging to each of the two samples. This can be noted from the histogram of Fig. 12 (left panel), where the two distributions overlap only marginally. Focusing on each of the individual subsamples, far sample sources show angular sizes within ∼0.15– 0.9″ (with approximate mean and median values of ∼0.25″), while near sample sources within ∼0.3–1.4″ (with mean and median values of ∼0.5″). Although the effective spatial resolution is rather uniform across our sample (right panel of Fig. 1), the dependence of the physical size estimate from the target distance is not completely removed, as it is shown in Sect. 5.1. A resolution-independent view of the source sizes can be obtained by normalizing the angular size by the beam size of the corresponding map, as reported in Fig. 12 (right panel). In these terms, detected sources have sizes from ∼0.8 to ∼2.4 times the beam, with a median value of 1.15 (i.e., 15% discrepancy). These values essentially match the expected outcome of the algorithm, as such compact source sizes are the ones CuTEx, by design, is more sensitive to (see Molinari et al. 2016b for details). The absence of a rigid lower limit at ratio ∼1 for the peak of the size distribution tells us that most of the detected compact sources are resolved. They turn out to mostly be slightly larger than the beam, implying that they are compact and centrally peaked objects. We note that for this ratio, even values below the imposed minimum source size threshold are allowed, since we are comparing with the circularized beam, while the lower boundary limit for the source size is set to 0.95 times the beam minor axis (see Sect. 3.2). We also caution that, since we are expressing quantities in circularized form for consistency with the subsequent physical analysis, ratios above 1 may in principle result, for some sources, from the combination of a non (or barely) resolved minimum FWHM, and a more extended maximum FWHM.

The overall distribution of source ellipticities (ϵ), defined as the relative discrepancies between the maximum and minimum axes of the fitted ellipse, i.e., (FWHM_max − FWHM_min)/FWHM_max, is shown in Fig. 13. A null ellipticity corresponds to perfectly circular sources, while values close to 1 indicate highly elongated sources. Observed ϵ range from just above 0 to ~0.7, meaning that the most elongated sources have a ratio of ∼3 between their major and minor extents. The median value of ϵ across the whole sample is 0.25 (i.e., a 1.25 ratio between the two axes). A physical interpretation can be made of these estimates, although some caveats related to the technical properties of the maps and the features of the source extraction method have to be considered. On the one hand, lower ellipticities could possibly be tracing mainly compact, centrally concentrated and isolated sources, while higher values could represent more diffuse ones, potentially residing in larger scale and/or filamentary-like structures. On the other hand, it must be considered that the map beam is always elliptical to some extent (mostly within ~1 and ~1.4 in terms of ratio between major and minor beam axes in our case, see Sánchez-Monge et al. 2025), and this has an effect on the photometry algorithm when searching for the best fitting of the source ellipse, so that we will never obtain perfectly circular sources.

Fig. 13 Histogram of the ellipticities of the sources (ϵ), defined as the relative discrepancy between the maximum and the minimum axes (i.e., the two FWHMs) of the fitted ellipse. Ellipticity ϵ = 0 means perfectly circular source, whereas large ϵ indicate highly elongated ones. The vertical dashed black line marks the median value ϵ = 1.25.

5 Physical properties of the cores

In this section, parameters obtained from the source detection and photometry procedure, reported in the compact source catalog of Sect. 4.1, are used, together with the properties of the hosting clumps, to estimate the main physical parameters of the extracted core population. In particular, the distributions of core physical sizes and masses are shown and discussed.

5.1 Core physical sizes

From the measured circularized angular size of the detected source (FWHM_circ, Sect. 4.3.2), its corresponding physical size (i.e., diameter) can be computed as $$D_{\text{core\hspace{0.17em}}}(\text{au})=FWH{M}_{\text{circ\hspace{0.17em}}}\left({}^{\prime \prime }\right)d(\text{pc}),$$(4)

where d is the clump heliocentric distance.

The estimated sizes of the cores are listed in Col. 3 of Table 7. Figure 14 (left panel) shows the overall distribution of sizes estimated for the cores. Core sizes range within ∼800–3000 au in most cases (∼90%), with a mean value of ∼1900 au and a median value of ∼1700 au, roughly where the distribution peaks. A few sources (∼2%) show sizes even below 800 au (down to ∼250 au), corresponding to the cases when the best possible resolution is achieved in the observations thanks to the combination of both good angular resolution and small target distance. Similarly, a small population of cores (∼3%) with estimated sizes above ∼4000 au is also present, corresponding to lower angular resolution and large distance cases. As explained in Sect. 2.3, however, some of these extreme cases could correspond to targets for which the distance estimates were recently corrected, thus switching between near and far samples. The uncertainty of estimated sizes coming from the photometry procedure is within ∼20%, as deduced by comparing the angular sizes of a large population of synthetic sources measured by CuTEx with their true simulated values (see Sect. 3.3 and Appendix B). The relative error on the clump heliocentric distance, originating from the V_LSR estimation at clump level (see Molinari et al. 2025), is on average 10%.

Recent subparsec resolution studies of the fragmentation of high-mass star-forming clumps, mostly performed with the ALMA, NOEMA or Very Large Array (VLA) interferometers both through dust continuum and line observations, have measured fragment sizes over a wide range of spatial scales, from ~0.1 pc (~20 000 au, e.g., Sánchez-Monge et al. 2013; Anderson et al. 2021; Avison et al. 2023; Morii et al. 2023; Traficante et al. 2023, all featuring maximum spatial resolutions within 2000–5000 au) down to ∼2000–3000 au (e.g., ALMA studies Sanhueza et al. 2019; Svoboda et al. 2019; Pouteau et al. 2022; Morii et al. 2023; Xu et al. 2024), and even ≲1000 au (Beuther et al. 2018, observed with NOEMA). In more detail, Sanhueza et al. (2019) and later Morii et al. (2023) measured fragment sizes between 3500–10 000 and 2000–20 000 au, respectively, in early stage 70 µm-dark clumps within 6 kpc distance at a >2000 au resolution. Pouteau et al. (2022) found ∼2500–10 000 au sizes in the single, large complex W43, hosting the young and more evolved protoclusters MM2 and MM3, observed at a 2500 au resolution. In their sample of evolved (L/M > 20 L_⊙/M_⊙) massive protostellar objects observed at high resolution (≲1000 au), Beuther et al. (2018) revealed a wide range of sizes, from ∼600 to ∼14 000 au. Traficante et al. (2023) reported fragments from ∼4000 to 20 000 au in size across their sample covering different evolutionary stages, from 70 µm-dark to HII regions, observed at ∼2000 au maximum resolution. Xu et al. (2024) (ASSEMBLE) measured core sizes between ∼2200–6200 au. It must be noted that, since our source extraction method is sensitive to compact sources up to nearly three times the beam wide (Sect. 3), we may be actually missing possible more extended structures which are instead revealed by other literature works.

The two distinct observational setups chosen for the near and far samples allowed us to obtain a rather uniform distribution of spatial resolutions across the whole clump sample, which resulted in uniform estimated core physical sizes, without clear separation between sources belonging to different subsamples. However, the spatial resolution still has an effect, as it can be seen from Fig. 14 (right panel), where cores of fields observed with a spatial resolution better than 1500 au are plotted as gold bars, while cores coming from less resolved fields (spatial resolution above 1500 au) appear as dark green bars. The two distributions partly overlap, buta slight shift is visible (a factor of 2 discrepancy between the two medians), proving that indeed the smallest cores revealed are potentially a bi-product of the better resolution available in those targets, whereas less resolved ones are responsible for the high-size end of the global distribution. Therefore, in less resolved clumps, we generally extract larger scale substructures than in better resolved fields. In other words, at higher spatial resolutions the centrally condensed/contrasted sources tend to match the beam size (see also right panel of Fig. 12), so that it is likely that most of the emission actually comes from a region smaller than the beam itself. In any case, from a statistical point of view, the vast majority of targets have effective spatial resolutions consistent within a factor of ∼2 (or 3 at most, see right panel of Fig. 1), so that within such limits we can be confident to be tracing similar kinds of physical objects. Nevertheless, it is still possible that at least in a few extreme cases (particularly at lower spatial resolutions) we could be actually probing different types of objects from one field to another. Such structures, in reality, may represent an upper level in a hierarchical framework of fragmentation, as well as unresolved smaller scale source(s).

Fig. 14 Distribution of the cores’ estimated physical sizes (i.e., diameters). Left panel: overall distribution of core sizes. The vertical dashed black line marks the median value of the distribution (1700 au). For better visualization purposes, the few outliers above 5000 au (∼2%) are not shown. Right panel: same distribution of left panel (black), but overplotted with the core size distributions resulting from better resolved fields (gold, below 1500 au) and less resolved ones (dark green, above 1500 au). The two histograms overlap in the lighter green area.

5.2 Core masses

From the measured source integrated flux (F_INT, Sect. 4.3.1) the core mass has been estimated assuming optically thin dust continuum emission through the formula (e.g., Kainulainen et al. 2017) $$M_{\text{core\hspace{0.17em}}}=\frac{R_{\text{gd}}{F}_{\text{INT}}{d}^2}{B_{\nu }(T){\kappa }_{1.3}},$$(5)

where R_gd = 100 is the widely assumed gas-to-dust mass ratio (see Sect. 5.2.3), and κ_1.3 = 0.9 cm² g⁻¹ is the dust absorption coefficient at 1.3 mm wavelength, describing the opacity of dust grains with thin ice mantles at the typical gas densities of these regions (≳10⁶ cm⁻³, Ossenkopf & Henning 1994; Sanhueza et al. 2019, see Sect. 5.2.3 for details). B_ν(T) is the Planck function at the core temperature T (see Sect. 5.2.1).

The uncertainties on the core mass estimate brought by the assumptions and parameters involved in the mass calculation are discussed in Sect. 5.2.3.

5.2.1 Core temperature model

A key parameter for the computation of the core mass is the temperature of the dust envelope enclosing the core, which determines its observed emission through the Planck function B_ν(T) (Eq. (5)). Ideally, direct measurements of the temperature for each source should be used, obtained from molecular line emission of high-density tracers such as NH₃ (ammonia), H₂CO (formaldehyde), CH₃OH (methanol), or CH₃CN (methyl cyanide) (see, e.g., Zhang et al. 2002; Sánchez-Monge et al. 2013; Lu et al. 2014; Beuther et al. 2018; Cesaroni et al. 2019; Olguin et al. 2021; Brouillet et al. 2022; Li et al. 2022; Gieser et al. 2023b; Taniguchi et al. 2023; Izumi et al. 2024, and Mininni et al. 2024). However, these estimates are not yet available for the whole target sample and all extracted cores at the time of this work, as they will be presented in a subsequent ALMAGAL work (Jones et al., in prep.).

Without such core-level estimates, educated assumptions on dust temperature for the cores must be made. Assigning a single, universal temperature to the whole core sample would represent a too coarse and inaccurate (although widely used) approximation, given the very wide range of different evolutionary stages and physical properties covered in our sample. Instead, we chose to base our assumptions on the properties of the respective hosting clumps, as it has been done in some other recent ALMA studies. In particular, several authors directly assigned to the cores the temperature estimated at larger scale for the hosting clumps. For example, Sanhueza et al. (2019) used 10– 15 K in their 70 µm-dark early stage high-mass clumps, whereas Svoboda et al. (2019) used 10–14 K and Zhang et al. (2021) 10–20 K in high-mass starless clump candidates. Csengeri et al. (2017) estimated dust temperature from the clump luminosity of their early stage massive clumps, then adopting a unified value of 25 K. Traficante et al. (2023), featuring a limited but varied sample of clumps in terms of evolutionary stages, split the sample in three groups according to the clump evolutionary indicator L/M, assigning increasing temperature levels (20, 30 and 40 K, respectively) to the cores belonging to more evolved clumps. For this work, we decided to adopt a similar but further developed approach, attributing to the cores of each clump a temperature according to its evolutionary stage, based on the calibration of the clump L/M with the rotational temperature of the high- density tracer CH₃C₂H performed in Molinari et al. (2016a) (see Appendix C). In our model we assign: (i) 20 K to cores belonging to less evolved targets (clump L/M ≤ 1 L_⊙/M_⊙); (ii) 35 K to core pertaining to intermediate sources (1 < L/M ≤ 10 L_⊙/M_⊙); iii) a temperature following the relation (with L/M expressed in L_⊙/M_⊙) $$T_{\text{core\hspace{0.17em}}}=21.1\cdot {(L/M)}^{0.22}\text{\hspace{0.17em}K},$$(6)

to cores in more evolved targets (L/M > 10 L_⊙/M_⊙), thus proportional to a power law of the clump L/M, scaled with a proportionality factor (21.1) to best fit to the observed points in the Molinari et al. (2016a) model. This corresponds to a temperature of 35 K for clumps having L/M = 10 L_⊙/M_⊙. The three groups of clump L/M and the corresponding reference temperatures, T_core , adopted for the cores are summarized in Cols. 1 and 2 of Table 6. The number of cores for each group is shown in Col. 5. The resulting distribution is reported in Appendix C.

With this approach, we are able to broadly account for the wide range of different evolutionary stages and physical conditions covered by our clump sample, which to some extent reflects on the properties of the internal population of cores. In the first group, i.e., with L/M ≤ 1 L_⊙/M_⊙, we do not expect the dust temperature to significantly vary from the ∼20 K that typically characterizes pristine clouds in prestellar phases (see, e.g., the above cited ALMA works). A slightly larger variance in temperature may be expected in the intermediate group, 1 < L/M ≤ 10 L_⊙/M_⊙, according to the number of protostellar sources already formed within the clump, but the fairly constant and scarcely scattered trend found in the Molinari et al. (2016a) model ensures a good robustness to our approximation. A wider range of temperatures is indeed possible in most evolved clumps, with L/M > 10 L_⊙/M_⊙, where formed protostars and young HII regions resulting from densest cores may be efficiently heating up the dust envelopes and clump molecular environment at large. By not setting a unique temperature for those sources, but rather a value proportional to the clump luminosity (and L/M), we intend to account for these effects.

Our model results in an overall range of core temperatures between 20 and 81 K, with a median value of the distribution of 35 K. Such values are consistent with those used in other recent ALMA surveys, such as Anderson et al. (2021) (18–76 K in IR- dark hub-filaments systems) and Pouteau et al. (2022) (19–65 K within W43-MM2&MM3). It must be noted nonetheless that all those works feature lower spatial resolutions (≥2500 au) with respect to our observations (1400 au on average).

Ultimately, a caveat has to be raised about these temperature assumptions. Cores belonging to the same clump will have a distribution of individual evolutionary stages (varying from still quiescent regions to active star-forming sites) and related physical properties (including temperature), which will not necessarily coincide with the overall behavior of the hosting clump, whose estimated physical parameters are mostly dominated by the brightest embedded object (see Elia et al. 2021, Elia et al., in prep.). It also remains possible that some hot cores in much evolved clumps could reach temperatures even higher than 80 K (∼100 K or even more, e.g., Beuther et al. 2005; Beltrán et al. 2011; Silva et al. 2017). Moreover, the presence of young HII regions may also have an influence, at least locally, on the properties of individual nearby cores, namely raising their effective temperature, potentially leading to misleading flux and mass estimates (see Sect. 6.1 and Appendix D).

The uncertainties associated with our temperature assumptions, and their impact on estimated core masses, are discussed in Sects. 5.2.3 and 6.

Fig. 15 Overall distribution of the estimated core masses. The vertical dashed magenta line marks the average mass completeness limit of 0.23 M⊙ estimated for our sample (see Sect. 5.2.4), while the vertical dashed black lines marks the median value of the distribution (0.4 M⊙).

5.2.2 Core mass distribution

The core masses, estimated with Eq. (5) adopting the temperature assumptions described in Sect. 5.2.1, are listed in Col. 6 of Table 7. Their overall distribution (CMF) is shown in Fig. 15, and ranges ~0.002–345 M_⊙, with a mean value of 1.5 M_o and a median of 0.4 M_⊙. The CMF peaks at ~0.3M_⊙. Thanks to the high spatial resolution and high sensitivity of the ALMAGAL observational setup, we are thus able to reveal dense, compact fragments across nearly 5 orders of magnitude in mass with large statistics. The sample is complete for core masses above 0.23 M_⊙ (see Sect. 3.3 and 5.2.4 for more details).

The observed mass range is wider than that found by Traficante et al. (2023) at slightly lower spatial resolutions (0.4–300 M_⊙ on >2000 au scales). As discussed in Sect. 4.2, clumps show a variety of degrees of fragmentation, and more crowded fields also report significantly wide internal distributions of core masses. We reveal massive cores (>30M_⊙) only in 28 cases (contained in 25 clumps), corresponding to 0.05% of the core catalog. The largest core masses we report are higher than those found by other recent ALMA and NOEMA studies with similar (or slightly lower) spatial resolutions, such as, Beuther et al. (2018) reporting up to 40 M_⊙, and Pouteau et al. (2022) with 70 M_⊙ at most. More generally, several interferometric studies have revealed core masses from ~0.1 to some tens (up to 70) of Solar masses (e.g., Sánchez-Monge et al. 2013; Beuther et al. 2018; Sanhueza et al. 2019; Svoboda et al. 2019; Pouteau et al. 2022; Morii et al. 2023; Xu et al. 2024). Other authors measured masses of 100 M_⊙ (e.g., Sadaghiani et al. 2020 in NGC6334 molecular cloud) or above, in particular up to ~150 M_⊙ (Gieser et al. 2023b, in a sample of 11 HMSFR from IRDC to HIIs), and ~300 M_⊙ (Traficante et al. 2023).

Detailed analysis of the uncertainties and observational biases influencing the measured core masses follows in Sect. 5.2.3. The variation of core masses with evolution is analyzed in Sects. 6 and Sect. 7.1.

5.2.3 Uncertainties and biases on core mass estimation

The mass estimate for the cores is potentially affected by different sources of uncertainties and observational biases (e.g., Battersby et al. 2010). The choice of the dust temperature which is assigned to the cores clearly plays a crucial role in this respect, representing in our case the largest source of uncertainty. To quantify its impact, similarly to what was done by Traficante et al. (2023), we set three ranges of temperatures around the reference values used for the three groups of sources described above (see Sect. 5.2.1). The three ranges are: (i) 10–40 K; (ii) 20–60 K; (iii) T_core ± 30 K (only allowing resulting values above 20 K for this group, for consistency). The overall range of temperatures explored was therefore 10–111 K. The three ranges were chosen in a rather conservative manner, being wide enough to allow for a proper exploration of the related mass variation while preserving reasonable fundamental assumptions on physical properties of objects given their evolutionary stage, that is, not allowing prestellar sources to go beyond T = 40 K or evolved ones to drop below T = 20 K (Traficante et al. 2023). Then, for each source, we computed the masses corresponding to the extreme temperature values of the range assigned to its group using the same Eq. (5) (M_lim), and evaluated the ratios between the reference core mass (M_core, computed using the reference T_core) and the two extreme mass values. We obtained mass discrepancies from a factor of ~1.8 to ~3, the latter cases occurring when the lowest temperature in the explored range is used, leading to higher masses than the reference ones (i.e., M_core/M_lim < 1). Columns 3 and 4 of Table 6 report, respectively, the range of temperatures explored to quantify the mass estimate variation and the corresponding range of obtained M_core/M_lim ratios, for each of the three source groups based on clump L/M.

Uncertainty associated with integrated flux measurement is ~20% (see Sect. 3.3 and Appendix B). It is worth noting that, at least in more evolved regions, continuum flux contamination from nearby HII regions (mostly from free-free radio emission) might occur. This issue is discussed in Sect. 6.1 and Appendix D.

The error in the clump distance (also included in Eq. (5)) is on average 10% (Sect. 5.1), which brings a 20% potential uncertainty in the computed core mass. Nevertheless, distance has the effect to practically delimit the range of masses that we are able to reveal across the whole sample of targets. This can be easily noted in Fig. 16, where the estimated core masses are plotted against the distance of the respective clumps, so that the various cores of a given clump are aligned along the vertical direction. While the highest core masses are almost constant across the whole sample, the lower masses are inferiorly bounded, and this limit, due to the dependency relation between flux and distance, increases with the square of distance. This behavior translates into a variation of the mass sensitivity of about 1 order of magnitude in total across the sample, implying that we are more sensitive to low masses (down to ~10⁻² M_⊙) in closer targets with respect to farther ones, where only M_core ≳ 0.1 M_⊙ are accessible (see Sect. 5.2.4 for details). Such trends are easily noticeable from the lines representing the median value, and the 15th and 85th percentiles of the core mass distribution within given bins of distance (see Fig. 16 for full description), which also give an idea of how objects are distributed over the covered mass range. In detail, the 15th percentile value (i.e, in the lower-mass range) increases by about 1 order of magnitude in total, whereas the 85th percentile (i.e., higher masses) by a factor of ~3. As a consequence, the core mass distribution within the far sample is more oriented towards higher masses, while the near sample will cover a broader mass range. Moreover, given the mass-distance trend of Fig. 16 and the linear dependence between the beam size and the distance, lower spatial resolutions lead on average to higher core masses, with a discrepancy of up to a factor ~100 between minimum masses obtained in most (≲1000 au) and least (>3500 au) resolved fields.

Another source of uncertainty in the mass estimate is the choice of dust opacity, which depends on the grain and gas density model adopted (see, e.g., Mathis et al. 1977; Ossenkopf & Henning 1994; Draine 2011). For our mass computation, we adopted the value κ_1.3 = 0.9 cm² g⁻¹ for the dust absorption coefficient at 1.3 mm wavelength of a distribution of coagulated MRN grains (Mathis, Rumpl & Nordsieck, Mathis et al. 1977, i.e., grains as homogeneous spheres of mixed composition) with thin ice mantles at gas densities ≳ 10⁶ cm⁻³ (Mathis et al. 1977; Ossenkopf & Henning 1994). This value allows for a consistent comparison with the results of several other interferometric studies of fragmentation at core scales, such as Palau et al. (2013), Beuther et al. (2018), Sanhueza et al. (2017), Sanhueza et al. (2019), Svoboda et al. (2019), and Morii et al. (2023). However, gas particle densities in our core sample are mostly found in the range ~10⁷–10⁸ cm⁻³ (see Sect. 5.3). Adopting a value of κ_1.3 ≃ 1.0–1.1 cm² g⁻¹, as obtained from the estimates reported in Ossenkopf & Henning (1994) for those densities, would cause a systematic decrease of ~10–20% of the estimated core masses. The uncertainty brought by opacity would be lower or at most similar to those due to temperature, flux, and distance. The correction factor would nevertheless change up to a factor of ~2 (used by, e.g., Palau et al. 2013 and Pouteau et al. 2022) if pre-coagulation or different grain composition models (such as naked grains or ice mantles with larger thickness, see Ossenkopf & Henning 1994) would be considered.

We assumed a widely used gas-to-dust mass ratio R_gd = 100 (e.g., Kainulainen et al. 2017; Sanhueza et al. 2017, 2019; Motte et al. 2022; Morii et al. 2023; Traficante et al. 2023), which at best constitutes only a rough estimate. While such value is quite well established for the Solar surroundings, its extension to the entire Galaxy is still under debate, as it has actually been found to vary as a function of the Galactocentric radius and local metal- licity (see Giannetti et al. 2017). Other values that can be found in literature range between ~110–150 (e.g., Beuther et al. 2018; Svoboda et al. 2019; Gieser et al. 2023b), and would imply an up to ~50% upward correction for calculated masses.

Ultimately, it is clear that, under our assumptions, mass uncertainty due to temperature is potentially significantly higher than the contributions of the other parameters involved. The true impact of temperature, but also of distance and spatial resolution, on the observed core mass function are further investigated in Sect. 6.2 and Appendix E.

Fig. 16 Estimated core masses as a function of the hosting clump heliocentric distance, highlighting the variation of mass sensitivity across distance. The red line draws the trend of the median value among core masses contained within 1 kpc bins of distance, and the two pink lines mark the position of the 15th (bottom) and 85th (top) percentiles of the mass distribution within the same distance bins. The horizontal dashed magenta line marks the mass completeness limit of 0.23 M⊙, while the solid black curve is drawn to approximately follow the trend of the minimum core mass revealed with target distance. The grey shaded area highlights the incomplete part of the sample, and, in particular, the region below the black curve represents the “forbidden” (i.e., not detectable) range of masses based on revealed objects at the different target distances. Note: for better visualization purposes, distance outliers (<3%) are not shown in this plot.

5.2.4 Mass sensitivity and completeness

The trend in Fig. 16 shows that our mass sensitivity (i.e., the minimum detectable fragment mass) varies by more than 1 order of magnitude due to heliocentric distance variations across our sample. While in the nearest targets (d ≃ 2 kpc) we are sensitive to masses down to ~10⁻² M_⊙, in the farthest ones (d ≳ 6 kpc) we are actually able to reveal only M_core ≳ 0.1 M_⊙. The dynamic range of recoverable masses is then higher for closer sources (about 4 orders of magnitude in mass) than for farther ones (~3 orders of magnitude in mass). An estimate of the mass corresponding to the 1σ_rms flux level can be derived by converting the average flux sensitivity of ALMAGAL maps of 0.05 mJy/beam (see Sect. 2.3) with Eq. (5) (using a median temperature of 35 K, assuming this flux as the emission of a point-like source). This leads to a mass range of ~2 × 10⁻³–4 × 10⁻²M_⊙ across the covered range of distances (~2–8 kpc). Dedicated estimates of the 1σ_rms mass for each field (M_1σ) can be achieved by employing the rms noise of the map σ_rms together with the proper distance and temperature of the clump. In these terms, we obtain a range of ratios between the minimum core mass $\left(M_{\text{core\hspace{0.17em}}}^{\min }\right)$ detected within each clump and its M_1σ of ~5–50, whereas the ratios between the maximum core mass for each clump and M_1σ are within ~8–1000. It must be noted that the ratios $M_{\text{core}}^{\min }/M_{1\sigma }$ do not exactly correspond to the S/N of the core detection measured in the photometry stage (see Sect. 3.2), due to the local primary beam correction and background subtraction operations computed by CuTEx.

We derived an average mass completeness limit for the entire ALMAGAL sample from the flux limits estimated from experiments with synthetic sources (see Sect. 3.3). Using the median temperature of 35 K and an average distance of 4.5 kpc, the corresponding mass completeness limits for the three groups of targets (classified based on their rms noise, see Sect. 2.3) are 0.14 (Group 1), 0.23 (Group 2), and 0.41 (Group 3) M_⊙, sorted by increasing map rms. Hereafter, being Group 2 by far the most statistically representative subsample, the value of 0.23 M_⊙ obtained for it will be used as the reference mass completeness limit, derived from the corresponding flux completeness value of 0.94 mJy. However, whereas integrated fluxes are essentially constant across the whole sample (as shown in Fig. 11, right panel) and thus flux completeness limits can be considered largely independent of clump distance, the trend between estimated masses and distance (Fig. 16) suggests to make specific mass completeness limit computations for near and far sample targets, separately, to account for the different sensitivity available at different distances within our sample. Using the most statistically relevant flux completeness value of 0.94 mJy (Group 2) and average distances of 3.5 and 6 kpc for near and far samples respectively, we obtain corresponding values of 0.13 and 0.37 M_⊙ in terms of mass completeness limit. These values are within a factor 2 from the average value of 0.23 M_⊙ adopted for the whole sample. Finally, we note that the effect of sensitivity in defining the range of detectable masses is substantially removed if we consider only masses above the average completeness limit.

5.3 Mass–size relation, mass surface density and molecular volume density

Figure 17 (left panel) shows the estimated core masses as a function of the core radius (see Sect. 5.1). Red points mark the median M_core values within given bins of R_core; they display a rather good correlation between the two quantities. Moreover, the best fit to the data points (solid blue line) corresponds to a power-law relation M ∝ R^α with α ≃ 2.6 (namely α = 2.612 ± 0.001), which is higher compared with the ones (~1.6–2.5) proposed in multiple literature studies for molecular clouds, clumps, and cores (e.g., Larson 1981; Kauffmann et al. 2010; Tan et al. 2014; Lombardi et al. 2010; Miville-Deschênes et al. 2017; Traficante et al. 2018; Urquhart et al. 2018; Ballesteros-Paredes et al. 2019). Starting from these two quantities, we computed the average mass surface density of the cores assuming a spherical shape as $${\text{Σ}}_{\text{core\hspace{0.17em}}}=\frac{M_{\text{core\hspace{0.17em}}}}{\pi R_{\text{core\hspace{0.17em}}}^2}.$$(7)

The resulting distribution of core surface densities is shown in the right panel of Fig. 17. Surface density Σ_core ranges from ~0.15 to ~500 g cm⁻², with a mean value of 4.5 and a median of 1.7 g cm⁻². This range is wider (especially in the upper end) than what was found in other similar studies, such as Palau et al. (2013) (~0.7–10 g cm⁻²), Sanhueza et al. (2019) (~0.8–5 g cm⁻²), Morii et al. (2023) (~0.05–10 g cm⁻²), and Traficante et al. (2023) (~0.1–30 g cm⁻²). It must be noted, however, that comparisons with those studies are not straightforward since, differently from the present work, they featured lower resolution observations of much more limited samples, and surface densities were in some cases computed using beam-deconvolved core sizes. Moreover, computed surface densities are affected by the uncertainties on temperatures used to derive the core masses. In more detail, ~32% of our cores have surface densities <1 g cm⁻², ~59% are within 1 and 10 g cm⁻², while only ~9% of them exceed 10 g cm⁻². Therefore, ~68% of the cores exceed 1 g cm⁻², which was suggested by Krumholz & McKee (2008) as the minimum threshold above which coherent collapse is allowed to form high-mass stars (although the influence of magnetic fields may overcome this, Tan et al. 2013). This percentage is considerably higher than the one found by Sanhueza et al. (2019) and Morii et al. (2023) (~10%). This is probably due to the fact that given our observational setup and target sample we are able to probe more dense and compact substructures on smaller scales down to ~800 au, in a broader variety of clumps in terms of evolutionary stage, including also more evolved and condensed protoclusters. Different levels of surface density are also displayed as a reference as solid gold lines in the left panel of Fig. 17, from bottom to top: 0.1, 1, 10, and 100 g cm⁻², respectively.

The gas volume density of the cores was also estimated assuming a spherical shape as (e.g., Zhang et al. 2018): $$n\left({\text{H}}_2\right)=\frac{M_{\text{core\hspace{0.17em}}}}{\overline{m}\left(4\pi R_{\text{core\hspace{0.17em}}}^3/3\right)},$$(8)

where $\overline{m}$ is the mean molecular mass, for which we assume 2.8 m_H (Kauffmann et al. 2008), with m_H the hydrogen mass. We found a range for n(H₂) of ~10⁶–10⁹ cm⁻³ (with a median of 3 × 10⁷ cm⁻³). Most of the cores (~64%) lie between ~10⁷ and 10⁸ cm⁻³. Our estimates are roughly consistent with what found by Beuther et al. (2018), while they are on average higher than in other recent works, which typically range from ~10⁵ to at most ~10⁸ cm⁻³ (e.g., Sanhueza et al. 2019; Zhang et al. 2021; Pouteau et al. 2022; Morii et al. 2023; Xu et al. 2024). This is likely due to the higher resolution achieved, and confirms that we are indeed revealing dense, compact substructures that represent the intermediate entities between clumps and (proto)stars in the star formation framework. The estimated mass surface densities and molecular volume densities of the cores are reported in Cols. 7–8 of Table 7, respectively.

Fig. 17 Mass–radius relation and mass surface density distribution for the core sample. Left panel: core masses as a function of estimated radius. Solid gold lines represent different values of core surface density: from bottom to top, 0.1, 1, 10, and 100 g cm−2, respectively. Percentages of data points within the corresponding bins of surface density are given in grey. Red points indicate the median Mcore values in given bins (red bars) of Rcore . The solid blue line is the best fit to the data points, corresponding to the power-law relation reported in the top-left corner. Right panel: distribution of the mass surface density of the cores. The vertical dashed black line marks the median value of 1.7 g cm−2.

5.4 Estimated physical parameters of the cores

Table 7 summarizes the main physical properties of the extracted population of cores (the full table can be accessed at the CDS in machine-readable form). The non-deconvolved core physical diameter has been derived from the circularized FWHM of the fitted source ellipse (see Sect. 5.1). Core temperatures have been assigned to the cores according to the reported evolutionary stage of the hosting clump as described in Sect. 5.2.1. The core masses and related surface and volume densities (Sects. 5.2 and 5.3) are also reported.

6 CMF variation with evolution

In this section, the distribution of estimated core masses (Sect. 5.2) is evaluated as a function of the clump evolutionary indicator L/M. The CMF is shown and discussed in terms of inverse cumulative number density of objects. Such representation has the advantage, with respect to the differential form, of being independent from the adopted binning, while preserving the possibility to reproduce and compare the slope of the high-mass tail (see, e.g., Testi et al. 1998; Olmi et al. 2013; Fiorellino et al. 2021; Pezzuto et al. 2023). To visualize and quantify the impact on the distribution of computed core masses of mass uncertainties due to the adopted model of core temperatures (Sect. 5.2.1), which represent the largest source of uncertainty (see Sect. 5.2.3), we performed a series of 1000 Monte Carlo (MC) simulations of masses of the whole population of extracted cores (see, e.g., Sadaghiani et al. 2020; Traficante et al. 2023). We assigned to each core a random temperature within the range ∆T reported in Table 6 and described in Sect. 5.2.3. The mass distributions shown in Sect. 6.2 will then have the corresponding uncertainties associated with temperatures reported as colored bands drawn by the 1000 simulated CMFs, lying in the background of the observed profile.

6.1 Core sample selection

For the present analysis, we filtered the core catalog in order to enhance the robustness of the results. First, we include here only cores with masses above the estimated average mass completeness limit for our catalog (0.23 M_⊙, Sect. 5.2.4). This allows us to avoid any sensitivity bias due to distance variations on the analysis of the shape of the CMF (see Fig. 16). We also underline that in our clump sample distance and L/M are not correlated (see Elia et al., in prep.). Moreover, we exclude cores that were positionally matched with strong radio continuum emission based on the analysis of the available ALMAGAL counterparts from the CORNISH (Purcell et al. 2008; Hoare et al. 2012) and CORNISH-South (Irabor et al. 2023) surveys at 5 GHz (or 6 cm). We identified 909 matched cores (~14% of the catalog, flagged in Col. 16 of Table 3), mostly associated, as expected, with the more evolved clumps (i.e., L/M > 10L_⊙/M_⊙ group). Such clumps are in fact compatible with the potential presence of HII regions, in which free-free continuum emission can bring a significant contribution to the measured core fluxes (e.g., Liu et al. 2022a; Motte et al. 2022; Pouteau et al. 2022; Louvet et al. 2024). For further safety, we also exclude cores belonging to clumps for which no CORNISH of CORNISH-South map was available (582 cores, ~9% of the catalog). This procedure was done to avoid possible mass overestimates for the cores resulting from flux contamination by free-free emission. The highest core masses are in fact the ones that influence the most the shape of the high-mass end of the CMF. Full details of this analysis can be found in Appendix D. Lastly, we exclude cores from clumps farther than 9 kpc, i.e., distance outliers (166 cores, <3%). After this selection process, we still end up with a large sample of 3324 cores (~53% of the catalog), almost equally distributed among the three evolutionary groups (from least to most evolved group, 991, 1181, and 1152 cores, respectively).

6.2 Evolutionary trend of the CMF

Figure 18 shows the distribution of estimated core masses within the three different subsamples of clumps selected based on their L/M ratio, according to the evolutionary classification presented in Sect. 5.2.1. We remark that the evolutionary stage of the clump will not necessarily reflect the one of the individual embedded cores. Assigning the same temperature to all the cores within the same clump according to its properties, although representing a reasonably good approximation, inevitably ignores the distribution of temperatures that embedded cores actually have due to their different individual evolutionary stage and local physical conditions (e.g., Hennebelle et al. 2020; Izumi et al. 2024).

The CMF profiles for the three different evolutionary groups appear as clearly shifted relative to each other, showing differences in shape and slope that go beyond the mass uncertainties due to temperature assumptions. Through the MC realizations, we overall cover a temperature range 10–111 K (see Table 6). We note that if we assumed the same temperature for all the cores, as done by other works, the shift between the three CMF profiles would be even larger (Eq. (5)). In more detail, the high-mass tail of the CMF appears to become flatter as evolution proceeds. In other words, more evolved clumps host more massive cores. We report core masses up to ~10M_⊙ in the L/M ≤ 1 group, up to ~30M_⊙ for 1 < L/M ≤ 10, while up to ~100 M_⊙ for L/M > 10 L_⊙/M_⊙.

Our result is, in principle, consistent with the findings of some recent ALMA surveys, which however feature much more limited statistics and different samples. For example, Pouteau et al. (2023) (ALMA-IMF) analyzed the variation in the observed CMF across different subregions within W43- MM2&MM3, which they broadly classified as corresponding to different evolutionary phases (from quiescent to post-burst). They found shallower high-mass end slopes with evolution, namely, more massive cores in more evolved regions. Nevertheless, it must be noted such study includes 205 cores extracted within a single large region, resolved at a ~2500 au resolution. Moreover, importantly, the inspected subregions are all evolved, so that only the more advanced part of the evolutionary sequence is traced. Analogous arguments can be made for similar results obtained in the W43-Main protocluster by Armante et al. (2024). Furthermore, Liu et al. (2023) selected 19 hub-filament systems from the samples of IRDCs of ASHES and ATOMS surveys and divided them into IR-dark and IR- bright sources based on the presence of IR emission in central hubs, to (non-uniformly) cover an evolutionary sequence spanning ~3 orders of magnitude in L/M (~0.1–100L_⊙/M_⊙). They found larger masses of the most massive cores (MMCs) in more evolved regions at ~3000 au average resolution. Morii et al. (2024))found higher mass dynamic range (i.e., higher MMC masses) in clumps with higher protostellar core fraction in ASHES. Xu et al. (2024) measured average core masses about 2– 3 times higher in ASSEMBLE evolved protoclusters with respect to ASHES earlier stage IRDCs (Sanhueza et al. 2019; Morii et al. 2023). Conversely, Anderson et al. (2021) (~2500 au resolution) found no clear indication of different behaviors at different evolutionary stages for the masses of the cores extracted from the 35 ATLASGAL IRDCs selected from Csengeri et al. (2017), which they subdivided into five evolutionary groups based on their mid-IR brightness.

Before focusing on the physical interpretation of this result (Sect. 6.3), we must take into account the fact that by constructing the overall CMF we are aggregating core populations extracted from a large and varied clump sample representing different physical conditions and environments. Such factors have an influence on core evolution, and thus on their mass distribution (e.g., Lee & Hennebelle 2018; Hennebelle & Grudić 2024). Because of this, we do not focus on detailed computations of the CMF slopes (as done, e.g., in Motte et al. 2018b; Sanhueza et al. 2019; Sadaghiani et al. 2020; Pouteau et al. 2022), which require an in-depth study that will be the subject of a dedicated subsequent ALMAGAL work. Besides, we find our CMFs to be reasonably consistent with a power law (see, e.g., Salpeter 1955; Kroupa 2001) only starting from relatively high masses of about 3 M_⊙ , while more resembling a log-normal in the lower-mass range.

Moreover, in addition to providing an excellent statistical relevance, the variety of our sample brings potential biases due to distance, resolution, and temperature effects. We already removed the influence of distance on mass sensitivity variation by cutting above the mass completeness limit. We also verified that the evolutionary trend of Fig. 18 stands when considering near and far sample sources separately, as expected from the fact that distance and L/M are not correlated in our targets (Sect. 6.1). We performed further analysis and tests to verify the reliability of the observed trend and reinforce our results, whose details are reported in Appendix E. First, we evaluated the trend of Fig. 18 for the subsamples of fields with spatial resolutions within a factor 2 (namely between 1000 and 2000 au, including most of the ALMAGAL targets). Second, we further analyzed the effect of temperature on the CMF (besides what already done through MC realizations) by adopting an alternative model based on Sadaghiani et al. (2020), in which the core temperature is proportional to its measured flux through a power law. Lastly, we considered a reduced version of the core mass formula Eq. (5) (essentially removing the temperature term), evaluating the evolutionary trend for the F_INTd² quantity. In all these cases (see Appendix E), the trend of Fig. 18 still holds, both in terms of shift and different slope among the three CMFs, confirming that none of these potential effects or biases are relevant in establishing the observed behavior.

Therefore, given the preliminary selection of the core catalog to avoid flux contamination (Sect. 6.1), and based on the above described tests we performed to exclude potential underlying biases, we consider the result of Fig. 18 to be definitely reliable and statistically robust.

Fig. 18 Observed CMFs (black lines, expressed as inverse cumulative number density distributions) for the three populations of ALMAGAL cores classified, as reported, according to the L/M evolutionary indicator of the hosting clump, as explained in Sect. 5.2.1. The three colored bands trace, for each population, the 1000 MC realizations of CMFs computed by assigning to each core a random temperature within the ranges defined in Table 6. Less evolved sources (i.e., with L/M ≤ 1 L⊙/M⊙) are marked by golden lines, the intermediate ones (1 < L/M ≤ 10 L⊙/M⊙) by green lines, and the more evolved sources (L/M > 10 L⊙/M⊙) by blue lines. Note: we include only cores above the mass completeness value of 0.23 M⊙, marked by the vertical dashed magenta line (see text for explanation).

6.3 Physical interpretation of the observed trend

The presented result potentially has several relevant physical implications. Importantly, the observed trend of the CMF, which thanks to the ALMAGAL data we are able to inspect with unprecedentedly large statistics, suggests that the shape of the CMF is not constant throughout the star formation process (e.g., Clark et al. 2007; Offner et al. 2014; Pouteau et al. 2022, 2023). The CMF rather builds with evolution, being limited to lower masses in the early stages, and then gradually accumulating mass (and flattening) with time, forming the most massive structures (>30 M_⊙, Fig. 18) only at advanced stages. If we would assume a core-to-star formation efficiency of ~30% (see, e.g., Matzner & McKee 2000; Alves et al. 2007; Elia et al. 2013, 2021), this would imply that only clumps with L/M > 1 L_⊙/M_⊙ could host cores already massive enough (~25 M_⊙) tobe capable of forming high-mass stars (M_* ≥ 8 M_⊙). According to this interpretation, in practice, we could be observing cores of similar initial (low) mass at a different time across their lifetime, so that most evolved ones would have simply had more time to grow in mass. It has to be noted nonetheless that the one-core to one-star correspondence is not straightforward, since massive cores have been observed to further fragment and form multiple stars (Cao et al. 2021; Li et al. 2024; Xu et al. 2024).

Interestingly, a compatible interpretation of the observed trend is a scenario in which more massive cores accrete mass at a higher rate than lower-mass ones, evolving over much shorter timescales when compared to the hosting clump, so that they soon dominate the global clump evolution (e.g., Schmeja & Klessen 2004; Contreras et al. 2018; Redaelli et al. 2022). Therefore, once an early stage clump (i.e., with low L/M) forms such massive fragment(s), it may accelerate its evolution, gradually (but rapidly) increasing its observed L/M (e.g., Sanhueza et al. 2019). As a consequence, in this scenario, high core masses (and late clump evolutionary stages) are reached sooner in time, so that, observationally, it becomes relatively unlikely to reveal a clump with a low L/M hosting massive fragments. Ultimately, the variation of the CMF with evolution implies that the stellar IMF is not a simple scaling of an observed CMF (see, e.g., Motte et al. 1998; André et al. 2014; Cao et al. 2021; Hennebelle & Grudić 2024).

In more general terms, such result is in favor of a clump- fed scenario for massive star formation, in which cores gain their mass in a competitive way, with the more massive ones gradually accreting mass more efficiently from the surrounding intraclump medium thanks to their high gravitational attraction.

Further results and discussions about these topics are presented in Sect. 7.

Fig. 19 Evolution of core masses. Left panel: median values of the maximum ($M_{\text{core}}^{\max }$, red) and the minimum ($M_{\text{core}}^{\text{min}}$, blue) core masses revealed within different bins of clump L/M (delimited by vertical dashed grey lines). The respective colored shaded areas mark the 25th and 75th percentiles of the two distributions within the same L/M bins. Right Panel: combined mass of the cores $\left(M_{\text{core\hspace{0.17em}}}^{\text{tot\hspace{0.17em}}}\right)$ minus the mass of the MMC within each clump, as a function of the latter. Data points are color-coded by the clump L/M. The solid black line draws the trend of the median value of the quantity on the y-axis in given bins of $M_{\text{core\hspace{0.17em}}}^{\max }$. The dashed grey line corresponds to the $M_{\text{core\hspace{0.17em}}}^{\max }=50\%M_{\text{core\hspace{0.17em}}}^{\text{tot}}$ relation, while the shaded grey area marks (from top to bottom) the 25–75% range.

7 Core mass growth and fragmentation

The result discussed in Sect. 6 suggests that compact cores grow in mass with time, so that massive cores (>30 M_⊙) are found only at more evolved stages. To investigate in more depth this scenario it is crucial to evaluate, in particular, whether it involves the entire population of cores embedded in a clump or just the more massive ones, and how this relates to (or affects) the fragmentation process in the region throughout the evolutionary sequence.

To maximize the outcome of such analysis, we employ a slightly different core sample with respect to Sect. 6. In detail, we now exclude all the cores embedded in clumps which reported at least one potentially free-free contaminated cores (see Appendix D), i.e., also the ones not directly associated with strong local radio emission. This is done in order to not artificially alter the internal mass distribution of the interested clumps, especially when it comes to evaluate the most and least massive cores. On the other hand, we waive the completeness cut in mass to preserve the information on the lower-mass core population. Ultimately, the considered sample amounts to 3973 cores (63% of the catalog).

7.1 Core mass growth with evolution

In Fig. 19 (left panel), we report the median value of the maximum and the minimum core masses ($M_{\text{core\hspace{0.17em}}}^{\max }$ and $M_{\text{core\hspace{0.17em}}}^{\text{min}}$, respectively) revealed across the clumps within different bins of L/M. The two quantities show a clearly different trend with the evolutionary tracer. The maximum core mass (i.e., the mass of the MMC) increases on average by a factor ~10 over about 3 orders of magnitude in L/M. The increasing trend essentially reflects the evidence of Fig. 18 and confirms recent observational findings (e.g., Liu et al. 2023; Pouteau et al. 2023; Morii et al. 2024; Xu et al. 2024), but also agrees with numerical simulations produced by Bonnell et al. (2004), assuming a subsequent competitive accretion model in the framework of clustered environment. The minimum core mass, instead, remains roughly constant on average around 0.2 M_⊙ (i.e., it slightly decreases within a factor of ~2, between ~0.15–0.3 M_⊙), roughly matching the mass completeness limit. As a consequence, the ratio between median $M_{\text{core\hspace{0.17em}}}^{\max }$ and $M_{\text{core}}^{\text{min}}$ increases from ~2 to ~40 with evolution (e.g., Morii et al. 2024). The roughly constant trend of $M_{\text{core\hspace{0.17em}}}^{\min }$ may be interpreted in different ways: (i) such objects might correspond to newly formed fragments that we reveal at all evolutionary stages before they start accreting mass in turn; or (ii) we are tracing the “same” core population that forms at early stages and then stays unaltered throughout the cluster evolution. In the first case, interestingly, these cores might then represent the low-mass seeds for massive star formation (e.g., Bonnell & Bate 2006; Traficante et al. 2023; Morii et al. 2024). According to this scenario, core mass growth and clump fragmentation would proceed essentially in parallel through the star formation process (see below for further analysis on this), so that a population of relatively low-mass cores is present throughout. We note that this scenario does not necessarily imply that early formed cores will survive throughout evolution (as they could for example further fragment). Nevertheless, the observed $M_{\text{core}}^{\min }$ trend might also be the result of a combination of scenarios (i) and (ii), as it is argued later.

Figure 19 (left panel) displays only the most and least massive cores within the clumps. In Fig. 19 (right panel), we report the total mass in cores within each clump minus the mass of the MMC, as a function of the latter. The color code represents the clump L/M. Although the quantities on the two axes inherently share a certain degree of correlation, the increasing trend on average between them, and the hint of a trend also with the evolutionary indicator, are particularly meaningful. They show that, on average, the whole core population within a clump (or at least the majority of it) undergoes mass growth during evolution (e.g., Morii et al. 2024; Xu et al. 2024). One would in fact expect a flatter (if not completely flat) trend in case where the MMC dominates over the companions to the point of monopolizing the mass reservoir available for accretion. In ~85% of the clumps, the MMC hosts between 25% and 75% of the combined core mass, while in some extreme cases (~2% of the targets) it reaches 90%. The overall increasing trend appears to be present regardless of the maximum core mass value, thus explaining the shift observed among the three CMF profiles of Fig. 18 throughout the mass range. Furthermore, we verified that the mass ratio between the MMC and its successive companions (in order of mass) increases with evolution, thus further suggesting a potential variation of the mass accretion rate of individual cores as a function of their mass, as discussed in Sect. 6. However, detailed estimations of mass accretion rates cannot be performed from these data, and require further dedicated study. The same Fig. 19 (right panel) trend between the mass of the MMC and the companions was predicted by the numerical realizations of Bonnell et al. (2004), while the growth of the total core mass in the clump with evolution agrees with MHD and HD models simulations by Hennebelle et al. (2022). We note that both works trace the formation and evolution of sinks, namely, (proto)stars.

7.2 Relation with clump fragmentation

To investigate whether such core mass growth is compatible with a continuous fragmentation occurring throughout evolution, as suggested by Fig. 19 (left panel), we relate the number of detected fragments (cores) with $M_{\text{core\hspace{0.17em}}}^{\max }$ (Fig. 20). An increasing trend of the median number of fragments (from ~2 to ~12) can be noticed between ~0.1 and ~10 M_⊙, so that, on average, more massive cores are found in more crowded clumps (e.g., Bonnell et al. 2004). This indicates that at least in such a mass range, core mass growth and further clump fragmentation can coexist. Above ~10M_⊙, the trend instead stabilizes at around 10–12 fragments. This might suggest that as soon as the most massive core(s) gets sufficiently massive, the fragmentation process within the clump is halted. Based on the previously discussed evidence from Figs. 18 and 19, this most likely occurs at later evolutionary stages. Indeed, we verified that most (~85%) of the MMCs above ~10 M_⊙ belong to clumps with L/M ~ 10 L_⊙/M_⊙ or higher (see also Elia et al., in prep.). According to this interpretation, the roughly constant trend of $M_{\text{core\hspace{0.17em}}}^{\min }$ in Fig. 19 (left panel) could indeed result from a combination the above-presented scenarios (i) and (ii), i.e., with new fragments forming until the MMCs reach the ~10 M_⊙ threshold, when the present low-mass population stabilizes in number and mass. We also verified that the overall trend of N_frag is confirmed when considering (instead of their masses) the integrated fluxes of the cores, so that temperature- and distance-related effects can be ruled out. The trend in Fig. 20 is compatible with the results of numerical simulations performed by Hennebelle et al. (2022).

However, the causality between the cores’ mass growth and the simultaneous interruption of the fragmentation is not straightforward, as various mechanisms might be in action. For instance, the MMC could be able to directly inhibit further fragmentation in the clump by monopolizing the mass reservoir and/or through feedback, although this is actually not expected based on the trend in Fig. 19 (right panel). On the other hand, such reservoir could be globally almost consumed at these late stages, so that there is not enough available material to form further fragments and/or to efficiently accrete them.

Furthermore, some caveats must be put on interpreting the trend in Fig. 20. First, the overall spread of the data points is rather wide, especially on the x-axis, so that in some cases we observe clumps with only a few fragments also at high $M_{\text{core\hspace{0.17em}}}^{\max }$ . Second, as it can be also noticed in Fig. 19 (right panel), the relatively poor statistics available for masses above ~10 M⊙ makes the trend less robust over that range, so we refrain from making conclusive assertions in this respect.

As a whole, the results discussed in Sects. 6 and 7 depict a scenario in which the formation of high-mass stars is linked to the cluster formation, with cores initially forming with relatively low masses and then growing in mass with evolution while new fragments keep forming, at least until $M_{\text{core\hspace{0.17em}}}^{\max }$ ≃ 10 M_⊙ , where further fragmentation might be inhibited. Cores might gain their mass from the material reservoir of the surrounding intraclump medium, thus favoring a clump-fed mechanism in a dynamic, multi-scale framework (see, e.g., Liu et al. 2023; Traficante et al. 2023; Morii et al. 2024). Gravity could play a major role in regulating the process, facilitating both the assembly of more fragments and the mass accretion onto cores (see Elia et al. in prep. and, e.g., Motte et al. 1998; Svoboda et al. 2019; Morii et al. 2023, 2024; Traficante et al. 2023).

Further exploration of the correlations among core- and clump-level properties (e.g., mass, surface density, evolutionary stage, age; see Elia et al., in prep.), as well as detailed study of the spatial distribution and mass segregation of the cores (see Schisano et al., in prep.), performed on the same ALMAGAL data, will be able to bring further, potentially crucial evidence in support of the scenario proposed in this work. Moreover, crucial confirmations will come from comparing the observational findings with the outcome of numerical simulations performed under conditions fully comparable to the ALMAGAL data (Nucara et al., in prep.).

Fig. 20 Number of detected compact fragments (cores) as a function of $M_{\text{core}}^{\max }$ for each clump. Points are drawn in transparency to highlight the local density of the data. The solid brown line draws the trend of the median Nfrag in given bins of $M_{\text{core}}^{\max }$ . The solid lighter brown lines correspond to the 25th (bottom) and 75th (top) percentiles of the Nfrag distribution within the same bins.

8 Summary and conclusions

In this paper, we present the results from an extensive physical analysis performed on the full sample of candidate high-mass star-forming clumps from the ALMAGAL survey. We investigated the fragmentation properties of a statistically relevant sample of 1013 clumps, representing wide ranges of Galactic locations, distances, masses, and the full evolutionary path. They were observed at 1.38 mm with ALMA to obtain a uniform high spatial resolution of ~1400 au (median value), enough to resolve the typical scales of dense cores. A compact source detection and photometry procedure was performed on combined continuum images with an optimized version of the CuTEx algorithm to obtain the ALMAGAL compact source catalog. We then characterized the fragmentation statistics and the photometric properties of the extracted core population, which were then used to estimate and analyze the main physical properties. The very large statistics obtained across the broad sample of clumps confers an unprecedented level of robustness and physical significance to our analysis.

The main results of this work can be summarized as follows:

1. The source extraction procedure revealed a population of 6348 dense fragments within 844 clumps (83% of the entire sample), with a mean number of detections per clump of 8 (and median of 5). A wide variety of clump fragmentation degrees is observed (from 1 to 49 fragments per clump), with 130 clumps (15%) reporting only 1 detection and 58 (7%) with more than 20 detected fragments. Compared with other recent surveys, we find a larger fraction of clumps with a low number of cores (28% with only 1 or 2 fragments);

2. The estimated core diameters are mostly (~90%) within ~800–3000 au, with a median value of ~1700 au (slightly above the median resolution of 1400 au). Based on the size distribution, most of the detected compact sources are resolved, and prove to be slightly larger than the beam (15% median discrepancy), implying that the revealed fragments are compact and centrally peaked;

3. The core masses were estimated from the measured integrated fluxes adopting for each core a dust temperature based on the evolutionary stage of the hosting clump as traced by the L/M indicator (leading to 20–81 K core temperatures). Thanks to the high spatial resolution and sensitivity of the ALMAGAL observations, we have been able to reveal compact fragments across about 5 orders of magnitude in mass with large statistics, a span that has not been achieved by any other recent high-resolution survey. The measured core masses range 0.002–345 M_⊙, with a median value of 0.4M_⊙. The average mass completeness limit for our sample is ~0.23 M_⊙. We revealed 28 cores with masses above 30M_⊙ (0.05% of the total). Moreover, the core masses show a rather good correlation with core radii, with a trend consistent with the widely proposed power-law relation M ∝ R^α with α ≃ 2.6. The estimated surface density of the cores ranges within ~0.15–500 g cm⁻² (median of 1.7 g cm⁻²). This range is wider (especially in the upper end) than that found in most of recent similar studies, likely due to our better spatial resolution allowing us to reveal denser, compact substructures. This is also confirmed by the rather high estimated core volume densities of ~10⁶–10⁹ cm⁻³ (median of 3 × 10⁷ cm⁻³);

4. We evaluated the variation in the observed CMF with evolution by comparing the distributions of three core subsamples selected according to the hosting clump L/M. Interestingly, the three CMF profiles are clearly shifted between each other, and the high-mass tail gets flatter with evolution, going from the least (L/M ≤ 1 L_⊙/M_⊙) to the most evolved (L/M > 10L_⊙/M_⊙) sources. In other terms, higher core masses are found at later evolutionary stages. In particular, cores above 30M_g are found only at L/M > 10 L_⊙/M_⊙. We verified that such results stand when taking into account uncertainties and biases brought by distance, temperature assumptions, mass sensitivity, and resolution. The observed trend suggests that the shape of the CMF is not constant throughout the star formation process, but rather builds (and flattens) with evolution, forming the most massive objects only at later stages. This result is also compatible with a scenario in which more massive cores accrete their mass at a higher rate than lower-mass ones, so that the higher masses are reached sooner in time. Such a trend of the CMF with evolution is consistent with what was found in some recent literature studies; however, we have broadly expanded the analysis in terms of available statistics and range of physical parameters explored;

5. We investigated the core mass growth as a function of evolution in detail. The median value of the mass of the most massive core in the clumps clearly increases with evolution, by nearly 1 order of magnitude over 3 orders of magnitude in L/M. The other cores within the clump also grow in mass together with the most massive core (MMC), and the mass ratio between successive companion cores is found to increase with evolution, further suggesting different accretion rates depending on mass. However, the mass of the least massive core stays roughly constant throughout evolution (~0.2 M_⊙). These trends suggest a scenario in which all cores grow in mass with evolution, while a population of lower-mass cores is present throughout, possibly corresponding to fragments newly formed at all stages. In this framework, such objects might represent low-mass seeds for high-mass star formation;

6. We related the maximum core mass to the clump fragmentation statistics. The number of detected compact fragments (cores) appears on average to grow together with the mass of the MMC, at least until (likely at later stages) the MMC reaches ~10 M_⊙, when clump fragmentation is seen to be apparently halted (although limited statistics prevents conclusive statements in this respect). Such a turning point might be directly due to the action of the MMC in the clump (e.g., feedback, monopolization of the mass reservoir) or related to more global effects, such as a gradual depletion of available reservoir in the clump environment that would subsequently reduce the efficiency for fragment formation and accretion.

Ultimately, our findings depict a likely clump-fed scenario for massive star formation, in which cores form as relatively low- mass seeds and then grow in mass with evolution, while further clump fragmentation proceeds in parallel (at least until late evolutionary stages). Cores might gain their mass from the intraclump medium in a dynamic, multi-scale framework, with more massive ones potentially accreting more efficiently than the others.

Our results are in agreement with the outcome of various recent surveys, which we have been able to extend thanks to a much larger and representative sample. Future ALMAGAL studies, together with systematic comparisons with a new set of numerical simulations, will allow for additional and possibly crucial evidence to support the proposed scenario.

Data availability

Full Tables 3 and 7 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/696/A151

Acknowledgements

We thank the anonymous referee for the warm reception of this paper and the insightful comments. AC and the Team at INAF-IAPS, RSK, P. Hennebelle, and LT, acknowledge financial support from the European Research Council via the ERC Synergy Grant ECOGAL (project ID 855130). AS-M acknowledges support from the RyC2021-032892-I grant funded by MCIN/AEI/10.13039/501100011033 and by the European Union ‘Next Gen- erationEU’/PRTR, as well as the program Unidad de Excelencia María de Maeztu CEX2020-001058-M, and support from the PID2020-117710GB-I00 (MCI-AEI-FEDER, UE). C. Battersby gratefully acknowledges funding from the National Science Foundation under Award Nos. 1816715, 2108938, 2206510, and CAREER 2145689, as well as from the National Aeronautics and Space Administration through the Astrophysics Data Analysis Program under Award No. 21-ADAP21-0179 and through the SOFIA archival research program under Award No. 09_0540. GAF gratefully acknowledges the Deutsche Forschungs- gemeinschaft (DFG) for funding through SFB 1601 “Habitats of massive stars across cosmic time” (sub-project B2), the University of Cologne and its Global Faculty Programme. P. Sanhueza was partially supported by a Grant-in-Aid for Scientific Research (KAKENHI Number JP22H01271 and JP23H01221) of JSPS. RSK acknowledges financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) “STRUCTURES”, and from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 50OO2206). RSK thanks for computing resources provided by the Ministry of Science, Research and the Arts (MWK) of the State of Baden-Württemberg through bwHPC and the German Science Foundation (DFG) through grants INST 35/1134-1 FUGG and 35/1597-1 FUGG, and also for data storage at SDS@hd funded through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). SDC is supported by the National Science and Technology Council (NSTC) through grants 112-2112-M-001-066 and 111-2112-M-001-064. TL acknowledges supports from the National Key R&D Program of China (No. 2022YFA1603100); the National Natural Science Foundation of China (NSFC), through grants No. 12073061 and No. 12122307; the international partnership program of the Chinese Academy of Sciences, through grant No. 114231KYSB20200009; and the Shanghai Pujiang Program 20PJ1415500. LB gratefully acknowledges support by the ANID BASAL project FB210003. The National Radio Astronomy Observatory and Green Bank Observatory are facilities of the U.S. National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.00195.L. 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.

Appendix A Source extraction tests using different algorithms

Source extraction runs were performed with the astrodendo (Rosolowsky et al. 2008) and CuTEx (Molinari et al. 2011) codes on ALMAGAL continuum maps, to characterize and compare their performance in terms of detecting compact sources and measuring their fluxes.

Because of this test and the following flux completeness analysis (see Appendix B), an IDL routine was developed to generate sets of "hybrid" continuum images of real ALMAGAL fields with synthetic compact sources injected. Synthetic sources are ideal tools to characterize the performance of a source extraction algorithm in a fully controlled and tunable environment (see, e.g., Molinari et al. 2011; Hopkins et al. 2015; Molinari et al. 2016b). Instead of producing fully simulated fields with artificial noise and emission patterns, we chose to use as a basis continuum images of ALMAGAL fields with different properties, on which we injected populations of synthetic sources at different flux levels, similarly to the approach recently used, for example, by Pouteau et al. (2022) and Nony et al. (2023). This was made to test the algorithms in conditions as similar as possible to the real maps.

To this end, a sample of 60 ALMAGAL fields was selected to properly represent the different features observed across the whole target sample: from relatively empty fields mainly characterized by diffuse emission to more complex ones showing a certain degree of structuring, to even crowded fields with multiple apparent compact bright objects. Synthetic sources are assumed to have a 2D elliptical Gaussian spatial brightness profile, and locations and position angles are randomly generated with uniform probability distribution within defined regions of the field. FHWMs are randomly generated allowing variation up to ±40% of the map beam size. A dedicated control is applied within the code on the potential blending among simulated sources, although we checked that, given the large statistics provided by the number of different map realizations and synthetic sources, and the broad region considered for source placement, source blending cases turn out to be extremely rare across our sample. In each simulation, a population of synthetic sources with a defined flux is generated. Different flux levels are also explored, so that multiple realizations for each ALMAGAL field are created.

For the source extraction test, 10 ALMAGAL fields were selected, in each of which 10 synthetic sources were injected at three different flux levels, thus generating a set of 30 maps with a total of 300 added simulated sources. Injected flux levels were 0.2, 0.5, and 0.8, in terms of fraction of the peak flux of the native map. To optimize their performance and provide a fairer comparison, a series of tests and tuning employing different combinations of source extraction parameters was preliminarily performed for both the astrodendro and CuTEx codes. They were then executed on all maps realizations, and their outputs were compared with the truth tables from the source simulation routine, adopting as metric the number of synthetic sources detected and the accuracy in measuring their integrated flux and size. The positional matching between extracted and simulated sources was performed using a 2 pixel radius as a threshold (nearly one- third of the beam size of 7M+TM2+TM1 maps). We point out that, differently from what done by CuTEx (which is based on the position of the local peak flux, see Sect. 3.2), astrodendro defines as source position the mean positions of the structure (leaf) along the x and y directions, since it does not assume any specific morphology for the sources to detect. This could make the positional matching with perfectly elliptical simulated sources less accurate in some cases, especially when the leaf is highly irregular in shape.

Figure A.1 (panel a) compares, for the three different generated flux levels, the detection maps obtained with astrodendro (left column panels) and CuTEx (right column panels) for a relatively complex field. Predictably, the number of synthetic sources recovered increases as they become brighter. However, CuTEx is able to detect all 10 sources already at the faintest source flux level, whereas astrodendro fails. The reason behind this is that CuTEx is able to deal with the presence of a strong background pattern by filtering it, while astrodendro, which works on plain thresholding of intensity image, can pick up a source only if it lies above the adopted threshold.

The two codes were then compared in terms of accuracy of the integrated flux measurement. Astrodendro does not include any background emission subtraction by default (see Sect. 3.1). This approach could still provide fairly accurate flux estimates in case of isolated compact sources in relatively simple fields. However, it is worth noting that in our continuum images of starforming regions we observe that compact substructures often appear to reside within or on top of diffuse emission zones, structured background, and/or filamentary structures (see, e.g., Fig. 6). As a consequence, based on CuTEx estimates, the contribution of the background emission to the total source flux is not negligible in our maps, being on average of about ∼25%, but varying from ∼10% to up to ∼90% in the most complex, extreme cases. To enhance its capability to recover correct fluxes and for a more significant comparison with CuTEx, an additional post-processing module was implemented in the astrodendro software, able to obtain a reliable estimate of the local background emission as the median intensity across a ring of pixels around each source, which was then subtracted from the total leaf flux. The compared outcome of the photometry procedure performed with CuTEx and this implemented version of astrodendro on the above described set of simulated maps is reported in Fig. A.1 (panel b). The implemented version of astroden- dro achieves a rather good accuracy in recovering the integrated fluxes of the detected synthetic sources overall (∼10% better than that achieved by CuTEx, in terms of median value of the relative discrepancy ∆), mostly due to the introduced background estimation. However, astrodendro records a relatively poor performance in terms of number of synthetic sources detected, going from the ∼50% for faintest sources to a ∼75% at most for brighter ones (for a combined ∼66% over the whole population). This suggests that the detection approach applied by astrodendro is likely not the most suited to efficiently reveal compact sources over a wide range of fluxes across our continuum maps sample. On the other hand, CuTEx performs perfectly in the source detection stage, where it records 100% of matched sources regardless of the flux level considered, while maintaining a rather high accuracy in measuring source fluxes (≲ 20% overall; see also Appendix B).

In addition, the accuracy in recovering correct source sizes was also compared between the two codes, with CuTEx reporting excellent results overall: measured FWHMs were consistent with simulated ones within 10 – 15% at most. While this is not surprising given the analogous assumption of elliptical Gaussian sources (see Sect. 3.2), it is worth noting that the sources are simulated over the local interferometric extended emission, which presents sufficient fluctuations to influence the observed source profile. Since astrodendro considers for the source size the ellipse best matching the leaf mask, in such cases even small distortions in the mask can substantially change the estimated size.

Ultimately, by compensating for a slightly worse flux recovery accuracy with a significantly higher source detection rate and size recovery accuracy, CuTEx reported the best overall performance in this series of compact source extraction tests, confirming to be more suitable for the aims of this work given the properties of the ALMAGAL continuum maps.

Appendix A.1 Aperture photometry tests with hyper

Additionally, we performed source extraction runs with hyper (Traficante et al. 2015) on the same set of simulated maps presented above, with the aim of testing the aperture photometry approach and its outcome compared to CuTEx.

The hyper code provides reliable compact source photometry through a combination of multi-Gaussian fitting and aperture photometry. Local background subtraction and source deblending are performed to account for variable backgrounds and crowded regions. Full details on its operation and application can be found in Traficante et al. (2015) and Traficante et al. (2023). For this work, we tuned the parameter setup to efficiently perform based on the properties of our images.

Figure A.2 shows the results obtained with hyper on the simulated maps in terms of compact source detection and flux estimation accuracy. The detection rate ranges from ∼50% (lower flux level) to ∼95% (highest flux level), for a combined ∼75% overall. The integrated flux recovery is rather accurate, being within 10% on average over the whole population. A direct comparison between hyper and CuTEx is shown in Fig. A.3, in terms of relative discrepancy of the estimated fluxes for the detected sources in common. The fluxes are generally lower for hyper, but overall compatible within at most 20% on average, a value which is consistent with the flux uncertainty we estimate for CuTEx (see Appendix B).

Fig. A.1 Comparison of the outcome of the source extraction tests performed with astrodendro and CuTEx. a) Comparison of the detections obtained with astrodendro (left column panels) and CuTEx (right column panels) on the continuum map of the ALMAGAL field AG304.5567+0.3273 added with 10 synthetic sources with three different flux levels (see text for details). The three flux levels reported to the side refer to fractions of the native map peak flux. Green contours/ellipses mark the detected sources that were matched with the truth tables from the simulation code (i.e., the recovered synthetic sources), while red/magenta contours/ellipses indicate detected but unmatched sources. b) Comparison of the outcome of the photometry procedure performed with the implemented version of astrodendro (left column panels) and CuTEx (right column panels) on the set of simulated maps (see text for details). In detail, histograms of the relative discrepancies (referred to as ∆) between measured integrated fluxes (F) and generated ones (F0) are shown for the matched synthetic sources detected at each flux level (100 injected sources in total, first three rows) and for all flux levels combined (300 injected sources in total, last row). Boxes within individual panels report the corresponding number of detections and the mean and median values of the relative discrepancy in the flux.

Fig. A.2 Same as panel b) of Fig. A.1, but for the source extraction tests performed with the hyper code.

Fig. A.3 Histogram of the relative discrepancies (referred to as ∆) between the integrated fluxes measured by hyper and CuTEx, for the matched synthetic sources combining all flux levels.

Appendix B Full description and outcome of the flux completeness and photometric accuracy analysis

With the aim of estimating the flux completeness levels for our extracted compact source catalog, and contextually determining the accuracy of the photometric algorithm applied, we tested the code on 60 fields from the ALMAGAL target sample, which we divided into three groups of 20 each representing different map noise properties (see Table 2). Based on the rms noise distribution shown in Fig. 3 (see Sect. 3.3), fields were selected within the ranges σ_rms < 0.06 mJy/beam (lower noise, Group 1), 0.06 ≤ σ_rms ≤ 0.1 mJy/beam (medium noise including the peak of the distribution, Group 2), and σ_rms > 0.1 mJy/beam (higher noise, Group 3). Moreover, through a visual inspection of the maps, we made sure i) to select fields showing only few possible evident compact sources at most, in order to avoid frequent cases of blending with the incoming simulated ones, and ii) that maps with various emission patterns (e.g., localized, filamentary, extended, etc.) were adequately represented, in order to increase the reliability of the completeness estimation on the selected sample with respect to the whole target sample.

For each field, using the same method presented in Appendix A, we generated a set of simulated images by injecting 20 synthetic elliptical 2D Gaussians with constant integrated flux, but with a distribution of angular sizes. The integrated flux levels to be injected were initially chosen in order to properly sample the measured flux distribution reported across our compact source catalog (Fig. 11), taken as reference to reproduce in the simulations the same properties of the observed core population for best results in code testing. Then, additional flux levels were added to the test runs for each group of fields to be able to more accurately sample the specific range where the 90% detection threshold defining the flux completeness level lied. Integrated flux levels ($F_{\text{INT}}^{\text{inj}}$ , Table 2) were generated in the range 0.2 – 2 mJy (11 levels) for Group 1 fields, 0.5 – 2.5 mJy (12 levels) for Group 2, and 0.7 – 2.4 mJy (15 levels) for Group 3. Considering all combinations, a total of about 4000 synthetic sources were generated for each group of fields, ensuring high statistical significance to the analysis. Source positions were randomly generated within the inner region of the map (up to a primary beam value of 0.6), and their size (i.e., FWHMs of the ellipse) was left free to randomly vary within 40% from the beam size. A limited control is applied to avoid full source overlapping. However, given the number of synthetic sources injected in each simulation, blending is not affecting the displayed results. In the case of real observations, there is a limited number of cases where sources are close enough that the partial blending affects their flux estimation.

We ran our CuTEx compact source extraction code on all the generated maps with the same parameters configuration used for the catalog extraction, and compared the photometry output with the "truth tables" of simulated sources adopting a matching radius of 2 pixels (about one-third of the map beam size). We then analyzed the detection statistics as a function of $F_{\text{INT}}^{\text{inj}}$ , in order to determine the flux completeness level for each of the three groups (see plots of Fig. 7). Determined completeness fluxes ($F_{\text{INT\hspace{0.17em}}}^{\text{compl\hspace{0.17em}}}$, Table 2) were 0.58 mJy for Group 1, 0.94 mJy for Group 2, and 2.05 mJy for Group 3, respectively.

To characterize the photometric accuracy of our algorithm, we then went to compare, for the matched detections, the estimated integrated fluxes (F_INT) with the generated ones, recording the relative discrepancy $\left(F_{\text{INT}}-F_{\text{INT}}^{\text{inj}}\right)/F_{\text{INT}}^{\text{inj}}$, for each of the simulated flux levels. Figure B.1 shows the peak value (mode) of the flux discrepancy distribution as a function of the injected flux level, for the three groups of fields. Figure B.2 instead reports the discrepancy distributions for the three groups aggregating all flux levels above the corresponding completeness limit. Recovered integrated fluxes are overall consistent with the simulations, presenting average values correct within ~20% for all groups across nearly all the flux levels.

Similarly, Figure B.3 shows the distributions of source angular size discrepancies between the generated and the recovered ones (F WHM^inj) expressed in terms of major and minor axes of the source ellipse, still considering all flux levels above completeness. Again, recovered sizes are correct within ~20% for all groups.

Fig. B.1 Photometric accuracy analysis: source integrated fluxes. Peak value (the mode) of the distribution of the relative discrepancy between measured (FINT) and simulated $\left(F_{\text{INT}}^{\text{inj}}\right)$ integrated fluxes, as a function of the injected flux level $F_{\text{INT}}^{\text{inj}}$ (Group 1 left panel, Group 2 middle panel, Group 3 right panel). In all panels, the vertical error bars mark the interquartile range of the distribution within each level, while the vertical dashed red line reports the estimated flux completeness limit for the current group. The horizontal dashed grey line marks the equality between FINT and $F_{\text{INT}}^{\text{inj}}$.

Fig. B.2 Photometric accuracy analysis: source integrated fluxes. Flux relative discrepancy distribution aggregating all flux levels above the completeness limit of the group (Group 1 left panel, Group 2 middle panel, Group 3 right panel). The vertical dashed grey line marks the equality between FINT and $F_{\text{INT}}^{\text{inj}}$ .

Fig. B.3 Photometric accuracy analysis: source angular sizes. Distribution of the relative discrepancy between measured (FWHM) and simulated (FWHMinj) sizes, in terms of major (top row) and minor (bottom row) axes of the detected source ellipse (Group 1 left panels, Group 2 middle panels, Group 3 right panels). As in Fig. B.2, only flux levels above the completeness limit are considered for each group. In all panels, the vertical dashed grey line marks the equality between FWHM and FWHMinj .

Appendix C Further details of the adopted core temperature model

The top panel of Fig. C.1 shows the calibration of the clump L/M with the rotational temperature of the high-density tracer CH₃C₂H performed by Molinari et al. (2016a), on which we based our core temperature model (see Sect. 5.2.1). Their study covers the range ~1 − 200 L_⊙/M_⊙ in L/M, so we took 35 K as the temperature value fitting the points within 1 − 10L_⊙/M_⊙ for the corresponding evolutionary group, and applied the powerlaw relation T_core (K) = 21.1 · (L/M)^0.22, representing the curve that fits the points at L/M > 10L_⊙/M_⊙, for the more evolved sources. For the L/M ≤ 1L_⊙/M_⊙ range, not covered by the Molinari et al. (2016a) model, we assumed a typical temperature of 20 K (see Sect. 5.2.1).

The distribution of core temperature resulting from our assumptions is shown in the bottom panel of Fig. C.1, and ranges between 20 and 81 K.

Fig. C.1 Details of the temperature model adopted for the cores. Top panel: Calibration of the clump L/M with the rotational temperature of CH3C2H, edited from Molinari et al. (2016a). The vertical dashed red lines mark different ranges of L/M as done in our work. The horizontal dashed green line mark the 35 K temperature value fitting the point within the 1 − 10 L⊙/M⊙ range, which has been adopted in this work for the correspondent evolutionary group. The full black line represents the power law fitting the points at L/M > 10L⊙/M⊙ (see Molinari et al. 2016a), also adopted in this work. Bottom panel: Overall distribution of the assumed core temperatures as a function of the clump L/M, resulting from the model used in this work. The vertical dashed red lines serve as in the top panel.

Appendix D Analysis of potential free-free emission contamination on ALMAGAL cores

Potential flux contamination from free-free continuum emission on detected ALMAGAL cores is to some extent expected in the ALMAGAL sample, given the very wide range of evolutionary stages that also includes evolved regions (L/M up to 500 L_⊙/M_⊙). Such clumps might in fact host HII regions, which are typically characterized by strong continuum emission at radio frequencies. This might bring significant contribution to the integrated fluxes of the sources that we measure on our maps (see, e.g., Liu et al. 2022a; Pouteau et al. 2022), ultimately leading to potentially overestimated core masses.

To evaluate this occurrence, we searched for counterparts of ALMAGAL fields in the CORNISH (COordinated Radio and Infrared Survey for High-mass star formation, Purcell et al. 2008; Hoare et al. 2012) and CORNISH-South (Irabor et al. 2023) surveys, providing radio continuum maps at 5 GHz (or 6 cm) with ~1″ angular resolution, obtained with the VLA. We found 894 correspondences (88% of the ALMAGAL target sample), for which we produced cutouts matching the ALMAGAL FOV (~35″). We computed the rms noise of the radio continuum maps $\left({\sigma }_{\text{rms}}^{\text{rad}}\right)$ through a sigma clipping procedure. In some regions, such estimation is not straightforward, since the noise can vary due to the presence of artifacts. In those cases, we visually inspected the maps to verify that the automatic estimation was reliable. We then generated flux contour levels at 3, 5, and 7 times the ${\sigma }_{\text{rms}}^{\text{rad}}$. Based on visual inspection, the latter ones proved to be the most realistic, i.e., to best represent the observed emission pattern on the maps, whereas lower levels included regions that did not seem to be really significant. Moreover, a stricter flux threshold leads to more reliable results in the source matching, reducing the impact of potential spurious correlations induced by the difficult estimation of the map noise. The clumps showing the presence of strong radio emission (i.e., above the $7{\sigma }_{\text{rms}}^{\text{rad}}$ threshold) were 192 (21% of the inspected maps) and are flagged in Col. 17 of Table 3. The coordinates of the ALMAGAL cores (Table 3) were matched with the $7{\sigma }_{\text{rms}}^{\text{rad}}$ radio contours allowing a tolerance of 1 radio continuum map pixel. Figure D.1 shows some examples of the analysis conducted. We noted the cores that fell within such contours, which we consider to be potentially significantly contaminated in their flux by radio continuum emission. The clumps containing at least 1 matched core were 110 (11% of the ALMAGAL sample), combining for a total of 909 matched cores (14% of the ALMA- GAL core catalog). Considering also the other unmatched cores pertaining to such clumps, we get a total of 1793 cores (28%). We underline that, as actually expected due to the nature of the regions, ~90% of matched cores belong to the more evolved clumps (i.e., with L/M > 10 L_⊙/M_⊙). To avoid contamination on the results coming from flux and mass overestimations, those cores were excluded from the physical analysis of Sects. 6 and 7. Furthermore, we also excluded cores detected within clumps for which no CORNISH or CORNISH-South counterpart was available for the above described analysis (582 cores or 9% of the catalog, coming from 119 clumps). We remark that this represents a rather conservative choice since, based on the above reported numbers, only a small fraction of such cores will actually be potentially contaminated. More generally, we emphasize the use of the term "potentially" throughout, since it is still possible that some of the positionally matched cores might actually not be significantly contaminated by the presence of local free- free emission, as well as some not matched nearby cores might be.

Overall, this analysis led to the identification of 1491 potentially contaminated cores (~23% of the catalog).

Fig. D.1 Cutouts of CORNISH/CORNISH-South radio continuum emission maps (5 GHz frequency, ~1′′ angular resolution) corresponding to the ALMAGAL fields (from top to bottom): AG305.5535-0.0114 (first panel), AG320.6753+0.2459 (second panel), AG326.4486-0.7484 (third panel), and AG333.2217-0.4020 (fourth panel). Red contour levels indicate $7{\sigma }_{\text{rms}}^{\text{rad}}$ emission regions. Red crosses mark the locations of the ALMAGAL cores that fell within such radio emission contours and are thus considered as potentially contaminated, whereas blue crosses mark the unmatched ALMAGAL cores.

Appendix E Further analysis on the effect of resolution and temperature on the observed trend of the CMF with evolution

In this appendix we report details of the further tests and analysis that were conducted to verify and reinforce the result of the evolution of the observed ALMAGAL CMF, discussed in Sect. 6. In particular, we inspect the potential influence of differences in resolution and of the temperature assumptions on determining the evolutionary trend shown in Fig. 18.

The spatial resolution of the observations can have a significant effect on the measured masses of the detected objects and thus on the shape of the CMF (see, e.g., Sadaghiani et al. 2020; Louvet et al. 2021). Although the distribution of the linear sizes of the beam is rather uniform across ALMAGAL maps (right panel of Fig. 1), the resolution can still vary in some cases by a factor up to ~3 with respect to the median value of ~1400 au. This might influence the shape of the CMF, especially in the higher-mass end where the tail can be driven by relatively few objects (e.g., unresolved resolution outliers). In principle, this should not be related to the evolutionary stage of the observed clump. We nevertheless evaluated the evolutionary trend of Fig. 18 for a more homogeneous subsample of cores observed with resolution between 1000 and 2000 au (thus within a factor 2 of discrepancy at most). Such range actually includes most of the ALMAGAL targets (see Fig. 2). The outcome of this is shown in Fig. E.1 (top panel). As for Fig. 18, we consider only cores above the 0.23 M_⊙ mass completeness limit. The evolutionary trend holds, with the more evolved profile actually appearing here to be more shifted with respect to the other two groups. We then consider the result to be not significantly affected by potential resolution effects across our sample.

To further evaluate the impact of temperature, we first applied an alternative model taken from Sadaghiani et al. (2020), for which the core temperature is assumed to be proportional to its measured flux through a power law, so that $${\tilde{T}}_{\text{core\hspace{0.17em}}}=T_{\text{core\hspace{0.17em}}}^{\min }\cdot {\left(\frac{F_{\text{PEAK}}}{F_{\text{PEAK}}^{\min }}\right)}^{\gamma },$$(E.1)

where $T_{\text{core\hspace{0.17em}}}^{\min }$ is the minimum assumed core temperature (20 K in our case), F_PEAK the measured core peak flux (Table 3), and $F_{\text{PEAK\hspace{0.17em}}}^{\min }$ its minimum value. A power law exponent γ = 0.238 is used from Sadaghiani et al. (2020) (found to reproduce the Salpeter IMF slope of −1.35, see Salpeter 1955). In this way, we are waiving our assumption of a single temperature for all cores within the same clump. We obtain a temperature range of 20− 137 K, and a range of corresponding masses (computed with the same Eq. 5) ${\tilde{M}}_{\text{core\hspace{0.17em}}}\simeq \text{\hspace{0.17em}}0.02-73{\text{M}}_{\odot }$ . We evaluated the evolutionary trend, which is shown in Fig. E.1 (middle panel). Again, the shift among the three CMF profiles is clearly present.

As a further, conclusive test, we analyzed the behavior of the F_INTd² quantity, thus basically considering a reduced version of the core mass expression (Eq. 5) in which the temperature term is removed. The result is shown in the bottom panel of Fig. E.1. The three profiles corresponding to the evolutionary groups are still clearly separated. The evolutionary trend is thus already well visible from the measured fluxes. Fluxes are of course physically related to the temperature of the sources, but the assumptions made to convert them into masses have no role in creating the evolutionary trend of the CMF.

Based on these tests, we can definitely rule out that the evolutionary trend that we observe for the overall ALMAGAL CMF (Fig. 18) is determined by the adopted temperature model. Ultimately, the result is then robust against both resolution and temperature effects.

Fig. E.1 Analysis conducted to test the result of Fig. 18. Top panel: Same as Fig. 18, but for a subsample of cores selected based on the spatial resolution of the corresponding continuum map (within 1000 – 2000 au). Middle panel: CMF for the three evolutionary groups with core masses obtained using the temperature model from Sadaghiani et al. (2020) (see text for explanation). Bottom panel: Distribution of the FINTd2 quantity for the three evolutionary groups, computed to obtained an analog of the CMF but without the influence of the core temperature assumption. The vertical dashed pink line marks the minimum threshold we imposed for the plotted quantity, calculated by converting our mass completeness limit.

References

All Tables

All Figures

Fig. 1 Beam properties of the ALMAGAL continuum maps. Left panel: probability density distribution of the circularized angular sizes of the map beam, separately for near (red) and far (blue) observational samples, respectively (see Sánchez-Monge et al. 2025 for more details). Right panel: distribution of the corresponding linear sizes (i.e., achieved spatial resolution). The grey histogram represents the overall distribution, with the vertical dashed black line marking its median value (1400 au). Near and far subsamples are overplotted with the same color coding of the left panel. It can be noted that, due to the distance variation across our sample, the bimodal distribution in the left panel converts into a fairly uniform beam linear size overall distribution (see text for details).

In the text

Fig. 2 Linear size of the map beam as a function of the heliocentric distance of the ALMAGAL clumps. Targets observed with near or far sample configurations appear as red or blue points, respectively. The vertical dashed grey line marks the nominal ~4.7 kpc distance limit separating the near and far subsamples (with the few exceptions described in Sect. 2.3). The horizontal dashed grey line separates higher spatial resolution fields (i.e., beam linear size <1500 au) from lower resolution ones (>1500 au). Note: for better visualization purposes, distance and/or resolution outliers are not shown in this plot (~30 targets, ≲3%).

In the text

	Fig. 3 Distribution of the map rms noise level σrms across the ALMA- GAL maps sample (see Sánchez-Monge et al. 2025 for more details). The vertical dashed lines mark the subdivision into the three reported groups of fields that are employed in Sect. 3.3.
In the text

	Fig. 4 Diagram of the main phases of the source detection and photometry pipeline developed and used in this work. Underlined black steps were already present in the original CuTEx version (but were still further implemented in this work), while those in red have been newly introduced in the context of the present work.
In the text

Fig. 5 Example of the source extraction procedure performed with CuTEx (see Fig. 4) for the ALMAGAL field AG349.6440-1.0948. Left panel: detection phase performed on curvature images along the reported four different directions; the blue crosses mark the candidate compact sources. Right panel: outcome of the photometry procedure performed on the corresponding intensity image; the blue ellipses represent the estimated shape of the compact sources (in terms of FWHMs and position angle, λPA) as a result of the 2D Gaussian fitting procedure (see Sect. 3.2 for details). Note: the color scale of the image was adjusted to highlight all the three sources at the same time, so their visual appearance in terms of emission extent might not correspond to their actual size (ellipses). The map beam is shown in red in the bottom-left corner of the plot. The angular and physical scales are also reported.

In the text

Fig. 6 Outcome of the source extraction procedure performed with CuTEx on selected ALMAGAL clumps (from left to right, fields AG316.4045- 0.4657, AG012.4977-0.2233, AG013.2428-0.0854). As in Fig. 5, the blue ellipses mark the compact sources detected at 5 σrms. The three clumps showcase the variety of fragmentation degree, morphologies of structures and patterns of emission revealed across the ALMAGAL sample. The CuTEx algorithm was optimized and tested in order to efficiently perform in all these different conditions. The green ellipses mark the sources that are recovered by allowing a signal-to-noise ratio (S/N) cut of 3 σrms. The map beam is reported in red in the bottom-left corner of each panel.

In the text

Fig. 7 Completeness levels (i.e., percentage of detection) as a function of the “true” integrated flux of the synthetic sources injected into the maps $\left(F_{\text{INT}}^{\text{inj}}\right)$ for the three different groups of targets selected based on their rms noise, as introduced in Sect. 2.3: Group 1 (left panel), 2 (middle panel), and 3 (right panel), respectively. The two dashed red lines in each panel identify the estimated completeness flux, $F_{\text{INT}}^{\text{compl}}$, which is reported in red.

In the text

Fig. 8 Fragmentation statistics deduced from the compact source catalog (Table 3). Left panel: overall distribution of the number of cores detected in each clump; detailed statistics are reported in Table 4. Right panel: overall distribution of the derived spatial density of the cores in each clump, computed within the entire FOV. The vertical dashed black line marks the median value of ∼10 pc−2 . In both panels, only clumps revealing at least 1 detection are included.

In the text

Fig. 9 Fragmentation statistics evaluated for near sample and far sample sources separately. Left panel: probability density distribution of the number of cores detected in each clump for near sample (red bars) and far sample (blue bars) sources, respectively. Right panel: same as left panel, but only including clumps containing core masses above 0.23 M⊙ (see text for interpretation). Detailed statistics for the individual distributions are reported in Table 5.

In the text

Fig. 10 Analysis of the derived spatial density of the cores. Left panel: distribution of the spatial density of the cores within each clump for near sample (red bars) and far sample (blue bars) sources, respectively. Right panel: same as left panel, but only including clumps containing core masses above 0.23 M⊙. For comparison, the overall distributions of the left panel appear here in transparency in the background as pink (near sample) and cyan (far sample) bars, respectively. In both panels, overlapping areas between the two distributions are marked by the magenta bars. The median values of the near and far sample distributions are also reported in each panel.

In the text

Fig. 11 Integrated fluxes of the extracted compact sources as estimated with CuTEx. Left panel: overall distribution of the integrated fluxes of the sources. The adopted flux completeness value of 0.94 mJy (see Sect. 3.3) is marked by the vertical dashed cyan line. The vertical dashed black line marks the median value of the distribution (2 mJy). Right panel: measured integrated fluxes of the sources as a function of the hosting clump distance. The solid blue line traces the median FINT values within given bins of distance. Light blue lines correspond to the 15th (bottom) and 85th (top) percentiles of the flux distribution in the same bins of distance. Note: for improved visualization, distance outliers (<3%) are not shown in this plot.

In the text

Fig. 12 Angular sizes of the extracted sources. Left panel: probability density distribution of the circularized FWHM of the sources. Red and blue bars include sources detected within the near sample and the far sample targets, respectively (see Sect. 2). Right panel: distribution of the ratio between the estimated angular size of the source and the map beam size (both in circularized form), providing a resolution-independent estimate of the source sizes. The vertical dashed black line marks the median value of 1.15 for the ratio.

In the text

	Fig. 13 Histogram of the ellipticities of the sources (ϵ), defined as the relative discrepancy between the maximum and the minimum axes (i.e., the two FWHMs) of the fitted ellipse. Ellipticity ϵ = 0 means perfectly circular source, whereas large ϵ indicate highly elongated ones. The vertical dashed black line marks the median value ϵ = 1.25.
In the text

Fig. 14 Distribution of the cores’ estimated physical sizes (i.e., diameters). Left panel: overall distribution of core sizes. The vertical dashed black line marks the median value of the distribution (1700 au). For better visualization purposes, the few outliers above 5000 au (∼2%) are not shown. Right panel: same distribution of left panel (black), but overplotted with the core size distributions resulting from better resolved fields (gold, below 1500 au) and less resolved ones (dark green, above 1500 au). The two histograms overlap in the lighter green area.

In the text

	Fig. 15 Overall distribution of the estimated core masses. The vertical dashed magenta line marks the average mass completeness limit of 0.23 M⊙ estimated for our sample (see Sect. 5.2.4), while the vertical dashed black lines marks the median value of the distribution (0.4 M⊙).
In the text

Fig. 16 Estimated core masses as a function of the hosting clump heliocentric distance, highlighting the variation of mass sensitivity across distance. The red line draws the trend of the median value among core masses contained within 1 kpc bins of distance, and the two pink lines mark the position of the 15th (bottom) and 85th (top) percentiles of the mass distribution within the same distance bins. The horizontal dashed magenta line marks the mass completeness limit of 0.23 M⊙, while the solid black curve is drawn to approximately follow the trend of the minimum core mass revealed with target distance. The grey shaded area highlights the incomplete part of the sample, and, in particular, the region below the black curve represents the “forbidden” (i.e., not detectable) range of masses based on revealed objects at the different target distances. Note: for better visualization purposes, distance outliers (<3%) are not shown in this plot.

In the text

Fig. 17 Mass–radius relation and mass surface density distribution for the core sample. Left panel: core masses as a function of estimated radius. Solid gold lines represent different values of core surface density: from bottom to top, 0.1, 1, 10, and 100 g cm−2, respectively. Percentages of data points within the corresponding bins of surface density are given in grey. Red points indicate the median Mcore values in given bins (red bars) of Rcore . The solid blue line is the best fit to the data points, corresponding to the power-law relation reported in the top-left corner. Right panel: distribution of the mass surface density of the cores. The vertical dashed black line marks the median value of 1.7 g cm−2.

In the text

Fig. 18 Observed CMFs (black lines, expressed as inverse cumulative number density distributions) for the three populations of ALMAGAL cores classified, as reported, according to the L/M evolutionary indicator of the hosting clump, as explained in Sect. 5.2.1. The three colored bands trace, for each population, the 1000 MC realizations of CMFs computed by assigning to each core a random temperature within the ranges defined in Table 6. Less evolved sources (i.e., with L/M ≤ 1 L⊙/M⊙) are marked by golden lines, the intermediate ones (1 < L/M ≤ 10 L⊙/M⊙) by green lines, and the more evolved sources (L/M > 10 L⊙/M⊙) by blue lines. Note: we include only cores above the mass completeness value of 0.23 M⊙, marked by the vertical dashed magenta line (see text for explanation).

In the text

Fig. 19 Evolution of core masses. Left panel: median values of the maximum ($M_{\text{core}}^{\max }$, red) and the minimum ($M_{\text{core}}^{\text{min}}$, blue) core masses revealed within different bins of clump L/M (delimited by vertical dashed grey lines). The respective colored shaded areas mark the 25th and 75th percentiles of the two distributions within the same L/M bins. Right Panel: combined mass of the cores $\left(M_{\text{core\hspace{0.17em}}}^{\text{tot\hspace{0.17em}}}\right)$ minus the mass of the MMC within each clump, as a function of the latter. Data points are color-coded by the clump L/M. The solid black line draws the trend of the median value of the quantity on the y-axis in given bins of $M_{\text{core\hspace{0.17em}}}^{\max }$. The dashed grey line corresponds to the $M_{\text{core\hspace{0.17em}}}^{\max }=50\%M_{\text{core\hspace{0.17em}}}^{\text{tot}}$ relation, while the shaded grey area marks (from top to bottom) the 25–75% range.

In the text

Fig. 20 Number of detected compact fragments (cores) as a function of $M_{\text{core}}^{\max }$ for each clump. Points are drawn in transparency to highlight the local density of the data. The solid brown line draws the trend of the median Nfrag in given bins of $M_{\text{core}}^{\max }$ . The solid lighter brown lines correspond to the 25th (bottom) and 75th (top) percentiles of the Nfrag distribution within the same bins.

In the text

Fig. A.1 Comparison of the outcome of the source extraction tests performed with astrodendro and CuTEx. a) Comparison of the detections obtained with astrodendro (left column panels) and CuTEx (right column panels) on the continuum map of the ALMAGAL field AG304.5567+0.3273 added with 10 synthetic sources with three different flux levels (see text for details). The three flux levels reported to the side refer to fractions of the native map peak flux. Green contours/ellipses mark the detected sources that were matched with the truth tables from the simulation code (i.e., the recovered synthetic sources), while red/magenta contours/ellipses indicate detected but unmatched sources. b) Comparison of the outcome of the photometry procedure performed with the implemented version of astrodendro (left column panels) and CuTEx (right column panels) on the set of simulated maps (see text for details). In detail, histograms of the relative discrepancies (referred to as ∆) between measured integrated fluxes (F) and generated ones (F0) are shown for the matched synthetic sources detected at each flux level (100 injected sources in total, first three rows) and for all flux levels combined (300 injected sources in total, last row). Boxes within individual panels report the corresponding number of detections and the mean and median values of the relative discrepancy in the flux.

In the text

	Fig. A.2 Same as panel b) of Fig. A.1, but for the source extraction tests performed with the hyper code.
In the text

	Fig. A.3 Histogram of the relative discrepancies (referred to as ∆) between the integrated fluxes measured by hyper and CuTEx, for the matched synthetic sources combining all flux levels.
In the text

Fig. B.1 Photometric accuracy analysis: source integrated fluxes. Peak value (the mode) of the distribution of the relative discrepancy between measured (FINT) and simulated $\left(F_{\text{INT}}^{\text{inj}}\right)$ integrated fluxes, as a function of the injected flux level $F_{\text{INT}}^{\text{inj}}$ (Group 1 left panel, Group 2 middle panel, Group 3 right panel). In all panels, the vertical error bars mark the interquartile range of the distribution within each level, while the vertical dashed red line reports the estimated flux completeness limit for the current group. The horizontal dashed grey line marks the equality between FINT and $F_{\text{INT}}^{\text{inj}}$.

In the text

	Fig. B.2 Photometric accuracy analysis: source integrated fluxes. Flux relative discrepancy distribution aggregating all flux levels above the completeness limit of the group (Group 1 left panel, Group 2 middle panel, Group 3 right panel). The vertical dashed grey line marks the equality between FINT and $F_{\text{INT}}^{\text{inj}}$ .
In the text

Fig. B.3 Photometric accuracy analysis: source angular sizes. Distribution of the relative discrepancy between measured (FWHM) and simulated (FWHMinj) sizes, in terms of major (top row) and minor (bottom row) axes of the detected source ellipse (Group 1 left panels, Group 2 middle panels, Group 3 right panels). As in Fig. B.2, only flux levels above the completeness limit are considered for each group. In all panels, the vertical dashed grey line marks the equality between FWHM and FWHMinj .

In the text

Fig. C.1 Details of the temperature model adopted for the cores. Top panel: Calibration of the clump L/M with the rotational temperature of CH3C2H, edited from Molinari et al. (2016a). The vertical dashed red lines mark different ranges of L/M as done in our work. The horizontal dashed green line mark the 35 K temperature value fitting the point within the 1 − 10 L⊙/M⊙ range, which has been adopted in this work for the correspondent evolutionary group. The full black line represents the power law fitting the points at L/M > 10L⊙/M⊙ (see Molinari et al. 2016a), also adopted in this work. Bottom panel: Overall distribution of the assumed core temperatures as a function of the clump L/M, resulting from the model used in this work. The vertical dashed red lines serve as in the top panel.

In the text

Fig. D.1 Cutouts of CORNISH/CORNISH-South radio continuum emission maps (5 GHz frequency, ~1′′ angular resolution) corresponding to the ALMAGAL fields (from top to bottom): AG305.5535-0.0114 (first panel), AG320.6753+0.2459 (second panel), AG326.4486-0.7484 (third panel), and AG333.2217-0.4020 (fourth panel). Red contour levels indicate $7{\sigma }_{\text{rms}}^{\text{rad}}$ emission regions. Red crosses mark the locations of the ALMAGAL cores that fell within such radio emission contours and are thus considered as potentially contaminated, whereas blue crosses mark the unmatched ALMAGAL cores.

In the text

Fig. E.1 Analysis conducted to test the result of Fig. 18. Top panel: Same as Fig. 18, but for a subsample of cores selected based on the spatial resolution of the corresponding continuum map (within 1000 – 2000 au). Middle panel: CMF for the three evolutionary groups with core masses obtained using the temperature model from Sadaghiani et al. (2020) (see text for explanation). Bottom panel: Distribution of the FINTd2 quantity for the three evolutionary groups, computed to obtained an analog of the CMF but without the influence of the core temperature assumption. The vertical dashed pink line marks the minimum threshold we imposed for the plotted quantity, calculated by converting our mass completeness limit.

In the text