© The Authors 2024

1 Introduction

The Atacama large millimeter array – initial mass function (ALMA-IMF¹) large program surveyed massive protoclusters of the Milky Way (2.5−33 × 10³ M_⊙, see Paper I by Motte et al. 2022). With distances spanning from 2 to 5.5 kpc, performing interferometric observations with ALMA was paramount to trace the dust and gas emission at a scale that probes the formation and evolution of pre-stellar cores and protostars in the dense gas of dusty filaments (see Paper II by Ginsburg et al. 2022 and Paper VII by Cunningham et al. 2023). Measuring the mass and thus the temperature of these structures is essential to understand the conditions in which stars form, and to constrain the core mass function (hereafter CMF), since the estimated core masses may vary substantially depending on the adopted temperature. In addition to the spectrum of masses, measuring the bolometric luminosity is indispensable to build the luminosity-to-mass ratio, a fundamental quantity that can be related to the evolutionary stage of protostellar objects (Motte & André 2001; Elia et al. 2010; Csengeri et al. 2016; Mottram et al. 2017). One of the challenges faced by the ALMA-IMF program is that millimeter observations, on their own, are insufficient to measure the bolometric luminosity. This paper addresses all of these issues and provides high-resolution (2.5″) column density, temperature, and luminosity maps based on a multiwavelength approach.

Constraining the bolometric luminosity, dust column density, and temperature requires the spectral energy distribution (hereafter SED) of dust grains to be modeled. While it is possible to accurately describe the scattering, absorption, and reemission of starlight through the dusty interstellar medium, phenomeno-logical approximations are more practical in most cases (see Galliano et al. 2018, and references therein). The modified black-body (MBB) description, a widely used approximation, allows such measurements to be inferred from an analysis of the far-infrared (hereafter FIR, 70 ≤ λ ≤ 500 µm) and millimeter fluxes (see Fig. 1, bottom panel). Assuming that the majority of the dust mass is in large grains (r > 0.02 µm, Galliano et al. 2018), and that these grains are in thermal equilibrium (because of their large enthalpy), both the mass M_dust and temperature T_dust can be measured through SED fitting based on the following equation: $$L_{\lambda }=M_{\text{dust }}{\kappa }_0{\left(\frac{\lambda }{{\lambda }_0}\right)}^{-\beta }4\pi B_{\lambda }\left(T_{\text{dust }}\right),$$(1)

where L_λ (Wm⁻¹) is the monochromatic luminosity, κ₀ (kg⁻¹) the dust mass absorption coefficient, β the opacity index (both tied to the dust grains’ physical and chemical properties), λ the wavelength, and B_λ the Planck function. It should be noted, however, that the MBB description cannot reproduce the dust emission at shorter wavelengths (λ < 70 µm), since the stochastic heating of very small grains and the contribution of aromatic features (Duley & Williams 1981; Leger & Puget 1984; Allamandola et al. 1985) result in a departure from the Planck function (see Fig. 1, bottom panel). In this paper, we do not attempt to accurately model the mid- and near-infrared domain of the dust SED (1 ≤ λ ≤ 70 µm), and instead focus on the emission of dust grains in thermal equilibrium. To this end, a thorough sampling of the SED above 70 µm is required to constrain both M_dust and T_dust through SED fitting, in particular in the 70–250 µm range, since the peak of the SED traces the dust temperature and opacity index.

We acknowledge that this method is subject to biases, since line-of-sight variations of the temperature and measurement noise can induce a degeneracy between T_dust and β (Shetty et al. 2009; Kelly et al. 2012; Galliano et al. 2018). Including observations longward of 250 µm can help to alleviate this issue. Several observatories can provide such FIR and millimeter measurements, namely Herschel, SOFIA, APEX and ALMA (see Sect. 2 and Table 1). Relying on a diversity of instruments and bands immediately poses a problem, illustrated in the top panel of Fig. 1: while the angular resolution of ALMA observations such as performed for the ALMAIMF program lies between 0.29″ × 0.26″and 1.52″ × 1.30″, the angular resolution of Herschel/SPIRE ranges from 17.6″ (at 250 µm) to 35.2″ (at 500 µm). The standard procedure (e.g., Galametz et al. 2012; Aniano et al. 2012; Giannetti et al. 2013; Köhler et al. 2014; Guzmán et al. 2015) for SED fitting involves smoothing the observations to the same angular resolution, that is, to the coarsest resolution (35.2″, in our case). Applying this smoothing procedure would entirely defeat the purpose of high-angular ALMA observations, undermining the immense usefulness of high resolution long wavelength data. Alternatively, Fourier-space combination of Herschel images with ground-based single-dish bolometer data would allow to work at an intermediate resolution, but the improved angular resolution attained with this technique remains coarse with respect to ALMA observations (e.g., 10″, Lin et al. 2016, 2017; 18″, Palmeirim et al. 2013; Könyves et al. 2020; Ladjelate et al. 2020).

To address this issue and retain the high-angular resolution information from ALMA observations, we employ the point process mapping (PPMAP) algorithm developed by Marsh et al. (2015). PPMAP allows us to combine and reproduce multiwave-length observations using the MBB description while preserving the spatial information contained in the higher angular scale maps. Recently, PPMAP was applied to reveal and constrain dust structures in supernova remnants (Chawner et al. 2019, 2020), star-forming filaments (Howard et al. 2019, 2021) and in the Milky Way disk (Bates & Whitworth 2023). All these studies were based on Herschel observations, and some included SCUBA-2 data (Holland et al. 2013). For the first time, Motte et al. (2018b) advanced PPMAP so far as to simultaneously fit observations from a data set that spanned two orders of magnitudes in angular resolution (from 0.37″ × 0.53″ to 35.2″). The results they achieved with a 2.5″ resolution in the massive W43-MM1 protocluster compelled us to apply this novel procedure to the analysis of the 15 ALMA-IMF fields. Through PPMAP, we perform SED fitting for a large set of continuum observations, while preserving the high-angular resolution capabilities of ALMA, providing a first step toward constraining the luminosity, column density, and temperature of the population of candidate cores and protostars.

In Sect. 2, we present the ALMA and complementary continuum observations toward the ALMA-IMF protoclusters. We then proceed to the analysis in Sect. 3, where we describe the PPMAP algorithm and the methods used to apply it to our specific problem. The derived luminosity, dust temperature and column density maps are presented and compared with previous studies in Sect. 4. Lastly, the construction of a PPMAP luminosity peaks catalog, encompassing luminosity and mass measurements, is detailed in Sect. 5, in which we discuss our findings.

Fig. 1 Observational constraints. Bottom panel: illustrative spectral energy distribution (SED), produced using a THEMIS grain mixture (Jones et al. 2017, gray curve). The blue curve represents a single modified blackbody that best fits the far-infrared and millimeter range of the SED. Black markers are overlaid on the gray curve to indicate the SED coverage enabled by the observations listed in Table 1, with the horizontal bars representing the bandwidths. Top panel: beam size of the observations used in our analysis, with the wavelength of the observations indicated on top of each marker (in micron, except for ALMA markers). Vertical bars represent the distance-induced variation in physical beam size across the ALMA-IMF sample.

2 Observations

In addition to the new observations obtained with the ALMA, the modified blackbody SED fitting analysis requires far-infrared data. To cover the necessary wavelength range between 70 µm and 870 µm, we have gathered observations from five different instruments. Additional mid- and near-infrared observations are required to derive bolometric luminosities. The details of these observations, along with the corresponding instruments, are provided in Table 1 and shown in Fig. 2. To complete the sampling of the SED below 70 µm and estimate the bolometric luminosity (see Sect. 3.3), we have also collected archival observations between 3.6 µm and 24 µm. In the following subsections, we present these different observations and the corresponding instruments in detail.

2.1 ALMA-IMF images

We used the continuum images presented in Díaz-González et al. (2023). These images are the result of combining the 3 mm (Band 3) and 1 mm (Band 6) ALMA-IMF continuum images presented in Paper I (Motte et al. 2022) and Paper II (Ginsburg et al. 2022) with the pilot of the Mustang-2 Galactic plane survey (MGPS90, Ginsburg et al. 2020) at 3 mm and the Bolocam Galactic plane survey (BGPS, Aguirre et al. 2011, Ginsburg et al. 2013). The more evolved ALMA-IMF fields have a significant contribution of free-free emission in the continuum images. To aid our photometry measurements, we also used the estimates of pure dust emission presented in Galván-Madrid et al. (in prep.), which subtract the free-free contribution to the continuum at 1 mm using the H41α recombination line within the ALMA-IMF data set. Appendix A.1.3 provide further details about this procedure. The absolute flux calibration for the ALMA observations is estimated to have an uncertainty of 10% as reported by Ginsburg et al. (2022).

Fig. 2 Data coverage chart. The color scale represents the percentage of observed and unsaturated pixels in each pair of region and map. Spitzer/IRAC data, which has no saturated pixels in the regions studied, is not shown here. The numbers account for the amount of pixels replaced through astrofix (example: in G012.80, the Hi-GAL 160 µm map has 6 saturated pixels that were interpolated). Maps that are either missing or discarded from the analysis are hatched.

2.2 ATLASGAL & APEX/SABOCA observations

The ALMA-IMF regions were mapped by the APEX telescope large area survey of the Galaxy (ATLASGAL, Schuller et al. 2009). APEX/LABOCA observations provide a 870 µm data point at an angular resolution similar to Herschel/SPIRE observations at 250 µm (see Table 1). Furthermore, 12 out of 15 ALMA-IMF regions were mapped by APEX/SABOCA at 350 µm (Lin et al. 2019). The absolute flux calibration uncertainties on SABOCA and LABOCA observations are estimated to be 20% and 15%, respectively (Lin et al. 2019; Contreras et al. 2013).

2.3 Hi-GAL survey

The ALMA-IMF regions have also been extensively mapped by the Herschel infrared Galactic plane survey (Hi-GAL, Molinari et al. 2010). Moreover, W43 was also imaged in the high-gain mode of PACS and SPIRE by HOBYS, a key imaging survey with Herschel (Motte et al. 2010; Nguyen-Lu’o’ng et al. 2013), to correct the field for saturation. The combination of PACS and SPIRE observations provides data points from 70 µm to 500 µm. However, it is worth noting that some Hi-GAL maps contain “NaN” (Not-a-Number) values, that correspond to pixels that are saturated. For instance, in the PSW band (250 µm) observations of the G012.80 (W33) region, there are 154 saturated pixels according to Table B.1 in Molinari et al. (2016). To address this issue, we have applied interpolation to replace the values of these saturated pixels. This interpolation process was performed using Gaussian process regression, as implemented in the astrofix Python package (Zhang & Brandt 2021), described in Appendix A. The extent of pixel replacement through interpolation is visualized in Fig. 2. When observations reach saturation levels that preclude interpolation, an alternative approach is to substitute them with SABOCA and/or SOFIA observations (Sect. 2.4). This solution is generally applicable, except in the case of G333.60, where the Herschel/SPIRE image at 250 µm is saturated, but there are not SOFIA observations at 214 µm (as shown in Fig. 2, sixth row). This makes G333 the least constrained region within our study. For Herschel/PACS and Herschel/SPIRE observations, the absolute flux calibration uncertainties are estimated to be 10% and 7%, respectively, as reported by Galametz et al. (2014).

2.4 SOFIA/HAWC+ observations

G012.80, G351.77, W51-E, and W51-IRS underwent observations at 53, 89, and 214 µm, conducted by Vaillancourt (2016) and Pillai & Simplifi Team (2023), employing the 2.7 m stratospheric observatory for infrared astronomy (SOFIA) telescope (Temi et al. 2018). These observations used the High-resolution Airborne Wideband Camera-plus (HAWC+; Harper et al. 2018). In our analysis of the G012.80, G351.77, W51-E, and W51-IRS dust emission, the SOFIA/HAWC+ data were employed to better sample the mid-infrared portion of the SED and to replace saturated Herschel/SPIRE 250 µm maps. The absolute flux calibration uncertainties for SOFIA observations are estimated at 15% for 53 µm and 89 µm, and 20% for 214 µm (Chuss et al. 2019).

2.5 Spitzer surveys

The ALMA-IMF protoclusters were imaged with Spitzer instruments. The MIPSGAL survey (Benjamin et al. 2003) covered the Galactic plane with the Spitzer/MIPS camera, while the GLIMPSE survey (Carey et al. 2009) provides Spitzer/IRAC observations. The absolute flux calibration uncertainties for Spitzer observations are estimated at 4% for MIPS at 24 µm (Engelbracht et al. 2007), and 2% for IRAC at 3.6, 4.5, 5.8 and 8.0 µm (Reach et al. 2005).

3 PPMAP description and analysis

3.1 Point process mapping

The point process mapping procedure, denoted as PPMAP (Marsh et al. 2015, 2017), is an iterative Bayesian SED fitting algorithm that allows to account for the mixing of physical conditions along the line of sight (dust temperature gradients and variations of the opacity index β). PPMAP is grounded in the point process formalism (Richardson & Marsh 1987, 1992, Marsh et al. 2006). The point process formalism represents complex astrophysical systems as an arrangement of individual components, referred to as “points” (e.g., Marsh et al. 2015). Points do not correspond to individual astrophysical objects in the image, nor to pixels; rather, they serve as elements that facilitate image representation. These elements are defined by a set of physical parameters, corresponding to a specific position within a state space. For our present purposes, the state space includes the position on the celestial sphere (x, y) and the dust temperature T_dust. Hence, the state space has a dimensionality N_states = N_temp × Nx × N_y, where N_temp, N_x and N_y are the number of temperature and positional cells that the points may occupy. The system can then be characterized by a vector containing the occupation number of each cell within the state space, denoted Γ. In essence, the local density of points corresponds to the density of dust at a specific sky position and temperature. Therefore, the dust column density within the n^th cell is determined by the corresponding occupation number Γ_n, and the set of occupation numbers Γ can also be viewed as a probability density function. The underlying PPMAP measurement model is expressed by the equation: $$d=A\Gamma +\mu .$$(2)

Here, d is the vector of observational measurement. The m^th component of d pertains to the pixel value at the coordinates (x_m, y_m) within the observed image at the wavelength λ_m. The term µ represents the measurement noise, assumed to be a Gaussian random process. Lastly, A denotes the system response matrix. Within this matrix, the mn^th element corresponds to the response of the m^th measurement to a point situated in the n^th cell of the state space (characterized by spatial position x_n, y_n, and temperature T_n). The PPMAP algorithm aims to solve this equation for Γ, given a set of observations d and an a priori distribution of points. This is achieved by minimizing the mean square error, ensuring that the best estimate is the a posteriori expectation value of Γ. The initial distribution (or “prior”) across all positions and temperatures is a Gaussian random process, expressed as: $$P\left({\mathrm{\Gamma }}_n\right)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(\frac{-{\left({\mathrm{\Gamma }}_n-\eta \right)}^2}{2{\sigma }^2}\right),$$(3)

where Γ_n is the occupation number in the n^th element of the state space, and $\sigma =\sqrt{\eta \left(1-1/N_{\text{states }}\right)}$. The quantity η is referred to as the “dilution”, since it controls the number of points relatively to the number of cells in state space. The initial distribution expressed by Eq. (3) means that points are equally likely to occupy any x, y position, and any given value in the user-defined log(T) temperature distribution.

PPMAP operates under the assumption that the radiation emitted by dust across observed wavelengths is optically thin. Consequently, the system response matrix A takes the form of the MBB approximation, expressed as: $$A_{mn}=H_{{\lambda }_m}{K}_{{\lambda }_m}\left(T_n\right)B_{{\lambda }_m}\left(T_n\right){\kappa }_{{\lambda }_m}\Delta {\Omega }_m.$$(4)

Here, H_λ represents the convolution operator associated with the point spread function (PSF) at wavelength λ, K_λ(T) accounts for a possible color correction, pertaining to the finite bandwidth of observations (cf. Appendix A.2.2), B_λ(T) denotes the Planck function, κ_λ corresponds to the dust opacity law, and ∆Ω denotes the solid angle corresponding to a specific pixel in the output map. The H_λ operator enables PPMAP to function without downgrading the spatial resolution of input maps, provided that the model benefits from accurate beam profiles. The PPMAP dust opacity law, κ_λ, exhibits a wavelength dependence parametrized by the opacity index β: $$\beta =-\frac{\text{d}\ln \left({\kappa }_{\lambda }\right)}{\text{d}\ln (\lambda )}.$$(5)

Thus, for any given λ, the dust opacity law κ_λ can be represented as: $${\kappa }_{\lambda }={\kappa }_{300}{\left(\frac{\lambda }{300\mu \text{m}}\right)}^{-\beta }.$$(6)

Here, β denotes the opacity power-law index, and κ₃₀₀ = 0.1 cm²g⁻¹ represents the reference opacity, measured at λ₀ = 300 µm, encompassing both dust and gas mass contributions. The selection of the dust absorption coefficient employed by PPMAP (κ₃₀₀ = 0.1 cm²g⁻¹) is in line with a gas-to-dust mass ratio of 100 (Hildebrand 1983). The value of κ₃₀₀ is identical for all points, and remains unchanged during the iterative process. Employing Eq. (6) with β = 1.8 yields κ_1.3mm = 0.007 cm²g⁻¹, which is consistent with the value adopted by Armante et al. (2024) following Ossenkopf & Henning (1994) $\left({\kappa }_{1.3 \text{mm}}={0.01}_{-0.0033}^{+0.005}\hspace{0.17em}{\text{cm}}^2\hspace{0.17em}{\text{g}}^{-1}\right)$. The reference opacity is in fact not well known, and may vary across the protoclusters. Depending on the size distribution and the composition of the dust, a range κ_{1 3mm} = 0 002–0 03 cm² g⁻¹ is predicted by Ysard et al. (2019) for the diffuse interstellar medium (cf. “Mix 1” and “Mix 2”, power-law size distributions). Aggregated grain models better represent denser media, and in that case the reference opacity is predicted to increase by a factor 3 to 7, depending on the addition of ice mantles into the models (Köhler et al. 2015).

Unlike conventional modified blackbody fitting approaches, PPMAP circumvents the need to homogenize the input observational data to a common resolution. Instead, when provided with a collection of observational data pertaining to dust continuum emission at varying instrumental resolutions, PPMAP produces maps of column density and temperature that are simultaneously optimized with respect to each specific PSFs associated with each dataset. Figure 3 provides a schematic illustration of the stepwise approach employed by PPMAP. In practical terms, PPMAP initializes an array given an a priori distribution of points. Subsequently, it generates a corresponding synthetic map for comparison with actual maps, accounting for synthetic noise. In each pixel of the maps, PPMAP minimizes the reduced-χ² metric, calculated from the deviations between components of the measurement model d and the corresponding observations, that are initially resampled to a common pixel size. Employing a truncated hierarchy of integro-differential equations (described in details by Marsh et al. 2015, Sect. 2.3), the distribution of points in the state space is iteratively updated until the model converges to match the observations. Upon completion of the process, the dust column density and temperature are given by the expectation value E(Γ_n|d). The H₂ column density is derived assuming a reference opacity of κ₃₀₀ = 0 1 cm² g⁻¹, fractional abundance by mass of hydrogen and molecular hydrogen X_H = 0.7, $X_{{\text{H}}_2}=1$, and fractional abundance by mass of dust Z_D = 0.01 (Howard et al. 2019, 2021). The differential column density cube $N_{{\text{H}}_2}$(T_dust) represents the H₂ column density at different dust temperatures, such that: $$N_{{\text{H}}_2}(i,j)={\displaystyle \sum _{k=1}^{N_{\text{temp }}}{N}_{{\text{H}}_2}}\left(i,j,T_{\text{dust },k}\right).$$(7)

Similarly, the density-weighted dust temperature map is defined by the following average quantity, based on the differential column density cube: $$T_{\text{dust }}(i,j)=\frac{1}{N_{{\text{H}}_2}(i,j)}{\displaystyle \sum _{k=1}^{N_{\text{emp }}}\left\{T_{\text{dust },k}{N}_{{\text{H}}_2}\left(i,j,T_{\text{dust },k}\right)\right\}}.$$(8)

Lastly, in accordance with Eq. (4), the synthetic intensity produced by PPMAP at wavelength λ in any pixel is determined by the following sum: $$I_{\lambda }(i,j)={\displaystyle \sum _{k=1}^{N_{\text{temp }}}\left\{\left(\frac{N_{{\text{H}}_2}\left(i,j,T_{\text{dust },k}\right)}{2.1\times {10}^{24} {\text{cm}}^{-2}}\right){\left(\frac{\lambda }{300\mu \text{m}}\right)}^{-\beta }{B}_{\lambda }\left(T_{\text{dust },k}\right)\right\}}.$$(9)

Here, $N_{{\text{H}}_2}$(T_d,k)/2.1 × 10²⁴ cm⁻² = τ₃₀₀ represents the optical depth at 300 µm. The numerical constant depends on the adopted values of κ₃₀₀, X_H, $X_{{\text{H}}_2}$ and Z_D.

Fig. 3 Schematic representation of the PPMAP iterative process. G012.80 images (at λ1 = 350 µm and λ2 = 870 µm) are used for illustrative purpose. In the upper part, we represent how PPMAP distributes “points” in a continuous state space (X, Y, Tdust) that can be divided into finite cells (corresponding to PPMAP pixels, that is, with a size fixed by the user, independent of the pixel size of the observed images). This distribution is then translated into a synthetic continuum emission map through the MBB description, taking into account the PSF of the instruments. Synthetic observations are finally compared to true observations, allowing to update the distribution of points. These iterative steps are repeated until the model converges.

3.2 PPMAP analysis

We apply the PPMAP procedure to the continuum data set detailed in Sect. 2. A comprehensive account of the methodology employed for implementing PPMAP to the data is provided in Appendix A. Here, we present an overview of the principles adopted throughout our analysis. To prevent contamination by free-free emission, we exclude the ALMA Band 3 data from our analysis (that is, the ALMA 3 mm continuum map, see also Appendix A.1.3 on the impact of free-free emission on PPMAP products). Consequently, the input images for PPMAP encompass the wavelength range from 70 µm (Herschel/PACS) to 1.3 mm (ALMA Band 6).

The opacity indices predicted by dust models describing diffuse and dense interstellar media are β = 1.5 and β = 1.8 (Köhler et al. 2015), respectively. There is in fact a range of plausible values, and β may vary across the field of observations, but we adopt a fixed value to minimize effects from the degeneracy between the temperature and opacity index and thus better constrain the temperature (Shetty et al. 2009; Kelly et al. 2012; Galliano et al. 2018). We therefore fix the opacity index to β = 1 .8, and we employ 8 MBB components with temperature values ranging from 10 K to 50 K, consistent with prior PPMAP applications (e.g., Marsh et al. 2017; Howard et al. 2019, 2021; Whitworth et al. 2019; Chawner et al. 2020;). Adopting a number of temperature components higher than 8 would result in larger uncertainties without improving the temperature sampling significantly. We experimented with larger temperature ranges and found that 10–50 K is sufficient to reproduce the observations. The final temperature estimate can exceed this temperature range following on the a posteriori temperature correction (see Sect. 3.4). The pixel sizes of the Nyquist-sampled PPMAP arrays are 1.25″, corresponding to an angular resolution of 2.5″. We execute the PPMAP analysis twice. The initial run (“Run1”) uses the ALMA B6 data decontaminated from free-free emission (as released by Galván-Madrid et al., in prep.), primarily created for deriving column density and dust temperature. The subsequent run (“Run2”) employs the standard ALMA B6 data (as released by Díaz-González et al. 2023), that is more optimal for accurately determining the luminosity by taking into account the millimeter excess tied to free-free emission. Column density and temperature estimates are directly obtained through the application of PPMAP, but deriving the bolometric luminosity requires additional steps.

3.3 Bolometric luminosity measurements

Using the outcomes of PPMAP Run2 (in order to take into account the free-free contribution to the luminosity), we computed the luminosity for all observed ALMA-IMF fields. The “PPMAP Luminosity” (L_MBB) is defined as the integral 4πd² ∫ I_vdv, where Ι_ν is defined by Eq. (9), and d denotes the distance to the observed star-forming region. On the other hand, the bolometric luminosity (L_bol) accounts for the additional near-infrared flux estimated from background-subtracted Spitzer/IRAC, Spitzer/MIPS, and SOFIA/HAWC+ observations, extracted from the MIPSGAL, GLIMPSE, and SOFIA archives (Carey et al. 2009; Benjamin et al. 2003; Vaillancourt 2016; Pillai & Simplifi Team 2023; also see Table 1). We performed a pixel-per-pixel merge between MIPSGAL, GLIMPSE, and HAWC+ observations and the PPMAP output SED using a piecewise cubic Hermite interpolating polynomial from the scipy package in Python. This results in a composite SED model as follows: $$\begin{array}{l}{I}_{\lambda }^{\prime }(\lambda <70\mu \text{m})=\mathrm{PCHI}(\lambda ),\\ I_{\lambda }^{\prime }(\lambda \ge 70\mu \text{m})=\mathrm{MBB}(\lambda ).\end{array}$$(10)

Here, “PCHI” denotes the cubic spline function employed for interpolating near-infrared data points, and “MBB” represents the best-fit PPMAP model. The bolometric luminosity is then defined as $L_{\text{bol}}=4\pi d^2{\displaystyle {\int }_{3.6\mu \text{m}}^{1.3 \text{mm}}{I}_{\lambda }^{\prime }}\text{d}\lambda$. The presence of saturation in the Spitzer/MIPS band occasionally resulted in lower limits on the infrared flux at 24 µm (refer to the first column in Fig. 2). Additionally, our integration approach effectively merges 2.5″ products with observations acquired at coarser resolutions (Spitzer/MIPS at 5.6″ and SOFIA/HAWC+ at 4.85″), resulting in a composite angular resolution. These two limitations have a relatively low impact on the outcome, given that the dominant contribution to the luminosity arises from the MBB emission (L_MBB/L_bol ≥ 0.77 across the 15 star-forming regions studied, where L_MBB and L_bol are respectively the modified blackbody and bolometric luminosities integrated over the field of observations).

Figure 4 displays a representative sample of illustrative spectral energy distributions (SEDs) extracted from the ATLASGAL sources’ footprints encompassing the protoclusters (Contreras et al. 2013, Urquhart et al. 2014). All of the observations align within ±2 standard deviation (σ) of the PPMAP MBB model, and 86% of the observations maintain agreement within ±1σ. We note that the ALMA error bars appear larger as a result of the integration area being large with respect to the 1.3 mm sources sizes.

Fig. 4 SEDs extracted from the ATLASGAL sources’ footprints (Contreras et al. 2013, Urquhart et al. 2014, see Sect. 3.3) corresponding to the protoclusters mapped by ALMA-IMF. The actual observations are represented by black points, while the red curve depicts the PPMAP MBB that provides the best fit to the data, with the gray shaded area representing the ±2σ standard deviation of the best fit. The blue curve represents the total piecewise model described in Sect. 3.3. Downward and upward arrows respectively correspond to lower limits and saturated observations.

3.4 Temperature correction and final products

As reported in Sect. 3.1, PPMAP assumes that the observed astrophysical object is optically thin to the thermal radiation emitted by dust across all input wavelengths. This primarily leads to a bias in the estimate of the dust temperature of deeply embedded sources, since the infrared fluxes at λ ≤ 250 µm may significantly deviate from the MBB shape (Men’shchikov 2016). To mitigate this limitation, in our analysis we have systematically applied an a posteriori correction to the PPMAP-derived temperature maps. After running PPMAP on synthetic observations generated with a dust model that incorporates the effect of the optical depth, the PPMAP outcome is compared with the input model dust parameters. We derive a correction table from this comparison, that can then be applied to PPMAP products (the dust model used to build this correction table is described in Appendix A.3). Figure 5 illustrates the change in temperature following the correction of the W43-MM1 temperature image. The outcome of the temperature correction is evident in localized areas, where the temperature is increased. For instance, in G012.80 the 99^th percentile temperature raises to 38.0 K from 33.2 K after the opacity correction. The magnitude of this correction scales with the column density estimated by PPMAP, as shown in Fig. A.3, therefore the high-density pixels benefit the most from it. As a result, the correction allows a better representation of embedded protostars and hot cores.

The final products delivered with this study are the maps of the H₂ column density ($N_{{\text{H}}_2}$ , in cm⁻²), bolometric luminosity (L_bol, in L_sun/px) and dust temperature (T_dust, in Kelvin) (see Fig. A.1 for reference). The dust temperature maps are declined in two versions:

1. The direct output of PPMAP, denoted hereafter as $T_{\text{dust }}^{\prime }$, represents the best-fit MBB temperature derived under the assumption of optically thin emission.

2. The a posteriori correction of the temperature yields a new estimate, denoted hereafter as T_dust.

The opacity-corrected temperature T_dust generally provides a better representation of the dust temperature, except in instances where the foreground is heated. Therefore, the temperature of internally heated, optically thick dust cores are best estimated by the second version of the temperature map, that we hereafter consider the default. Impending studies will attempt to determine the best combination of both maps based on the identification of prestellar cores and candidate protostars in the ALMA-IMF fields (Motte et al., in prep.).

An other caveat tied to PPMAP-derived temperatures is the impact of the 2.5″ angular resolution. Emission of the ALMAIMF cores with a median size of ~2100 au (Motte et al. 2022) is diluted in a PPMAP beam of physical size 6000−14 000 au, depending on the distance of the region. This dilution of the signal may result in underestimating the temperature of proto-stellar cores and overestimating that of prestellar cores, since cold and warm dust are mixed in the 2.5″ beam, as well as along the stratified line of sight. Therefore, additional processing should be applied to PPMAP temperature maps prior to their use in constraining core temperatures. This endeavor is also being undertaken by Motte et al. (in prep.), to which we refer the reader for a detailed account of the methodology.

Fig. 5 Temperature correction of the PPMAP images illustrated. The left panels displays the W43-MM1 Main-West image before correction, while the right panels presents the post-corrected map. White ellipses outline continuum cores identified by Louvet et al. (2023) in the ALMA 1.3 mm images at 0.4–0.9″ angular resolution.

4 Results, validation, and caveats

Following the procedures described in Sect. 3 and expanded upon in Appendix A, we have applied PPMAP to multiwave-length observations and obtained luminosity, column density and temperature maps at a 2.5″ angular resolution. We here present these outputs, evaluate their reliability against an accepted reference and then discuss their uncertainties.

4.1 Results

Figure 6 illustrates the improvement in angular resolution achieved through PPMAP’s application. We used the temperature and column density images of G353.41 as examples to illustrate the importance of gaining angular resolution for the ALMA-IMF studies. In this figure, we contrast the outcomes of the PPMAP approach at a 2.5″ resolution with the more typical approach, that requires smoothing all input images to the same angular resolution. To make this illustration, we first substituted the Herschel/SPIRE image at 350 µm with the SABOCA image and excluded the Herschel/SPIRE image at 500 µm, thus preventing further smoothing to a 35.2″ resolution. We then smoothed all continuum images to the coarsest angular resolution, which is that of LABOCA observations, 19.2″. Finally we performed SED fitting using PPMAP.

Complete representations of PPMAP luminosity, column density and temperature maps obtained for all regions are shown in Fig. B.1. The spatial variations of the reduced χ² square metric are also shown in Fig. B.2 for all regions. We here present an overview of these PPMAP data products for a subset of regions, selected to provide one example for each evolutionary stage (young, intermediate, evolved, as outlined by Motte et al. 2022).

In Fig. 7, we present the column density, bolometric luminosity and dust temperature maps obtained for this specific subset. These temperature images correspond to those corrected for the opacity because we consider them to best represent the dust temperature in dense regions. In the following, we simply call them the temperature images (see Sect. 3.4 and Appendix A.3 for a description of the temperature correction procedure).

With a 2.5″ angular resolution, PPMAP captures the morphology of filamentary structures, pinpoint the location of cores, protostars, HII regions, and constrain their surrounding physical conditions. The angular resolution provided by PPMAP offers an insight into the relationships between the continuum cores identified by Louvet et al. (2023) in the ALMA-IMF 1.3 mm continuum images and the column density and dust temperature maps. These cores, with typical sizes of 0.4–0.9″ (equivalent to ~2000–4000 au), align with filaments and aggregate within central hubs, a trend depicted in Figs. 6–7. Additionally, within the temperature images we observe correlations between massive, hot protostars and warmer spots, while HII regions appear as extended areas of enhanced temperature (as evidenced in the evolved protocluster G012.80, shown in the bottom panel of Fig. 7, also see Armante et al. 2024). This underscores the pivotal role of PPMAP’s resolution-optimization capacity in achieving our goals.

Table 2 presents the luminosity measurements using two distinct approaches:

1. The bolometric luminosity (L_bol ) measured in the primary beam response of the ALMA 12 m array mosaics (hereafter referred to as the “ALMA-IMF mosaic footprint”, as outlined in Paper I, Motte et al. 2022, Fig. 1).

2. The bolometric luminosity $\left(L_{\text{bol}}^{\prime }\right)$ measured in the ATLASGAL sources’ footprints, as outlined in Contreras et al. (2013) and Urquhart et al. (2014). These sources were extracted from the ATLASGAL survey (Schuller et al. 2009) using the source extraction routine SEXtractor. Protoclusters imaged by ALMA-IMF are associated with either one or two ATLASGAL sources (two sources in the case of G008.67), that are always smaller than the ALMAIMF mosaic footprints. An example of single ATLASGAL source is shown in the bottom panel of Fig. 6.

Table 2 also provides a compilation of average dust temperatures, peak column densities and total masses. The regions are classified according to their evolutionary stages, as defined in Motte et al. (2022), categorized as “young”, “intermediate”, and “evolved”, in descending order in the table. The luminosity-to-mass ratio, L/M, exhibits a noticeable trend aligned with this classification: younger regions generally correspond to lower ratios (L/M = 20 on average for the young protoclusters listed in Table 2), and vice versa (L/M = 76 on average for the evolved protoclusters listed in Table 2). There is also a trend of increasingly high temperatures across the evolutionary stages, with a 99^th percentile temperature of 33.7 K, 35.4 K and 38.7 K respectively for the young, intermediate and evolved regions. The average temperatures presented in the table, ranging from 21 K to 29 K, are representative of broader regions rather than locally heated zones (see Sect. 4.1). The 99^th percentile temperature ranges from 28 K to 42 K, depending on the specific region. As such, the temperatures inferred by PPMAP within the vicinity of massive cores are higher than those measured by Wienen et al. (2012, 2018), König et al. (2017) and used in Motte et al. (2022) as the means to estimate core masses (T_dust = 20–30 K). The average column density measured in the 2.5 maps ranges from 2.5 × 10²² to 2.5 × 10²³ cm⁻², about one to two decades higher than those measured in the wide Herschel images of nearby star-forming regions (Arzoumanian et al. 2011, Hill et al. 2011).

Fig. 6 Resolution enhancement of the dust temperature (top panels) and column density (bottom panels) images, achieved through the application of PPMAP to the G353.41 dataset. The left panels display the maps derived by smoothing all input images to a uniform resolution of 19.2″, while the right panels represent the PPMAP images at an angular resolution of 2.5″. White ellipses outline continuum cores identified by Louvet et al. (2023) in the ALMA 1.3 mm images at 0.4−0.9″ angular resolution. The larger dashed circles represent the footprint of the ATLASGAL source AGAL353.409–00.361 (Contreras et al. 2013, Urquhart et al. 2014).

Fig. 7 PPMAP products illustrated for three example regions: the young W43-MM1 (top), intermediate G008.67 (center), and evolved G012.80 (bottom) protoclusters. From left to right: column density map (N(H2)), bolometric luminosity (Lbol), dust temperature (Tdust). White continuous contours outline the ALMA 1.3 mm mosaic areas. The luminosity peaks extracted from the PPMAP luminosity maps (see Sect. 5.2) are overlaid in gray. The continuum cores identified by Louvet et al. (2023) in the ALMA 1.3 mm images are overlaid in white. The size of the ellipses reflects the FWHM of the sources.

4.2 Comparison to results from ATLASGAL

In order to benchmark our results, Fig. 8 compares the bolo-metric luminosity and mass measured in the PPMAP images to those obtained by König et al. (2017). Based on a two-temperatures MBB description, they inferred bolometric luminosities and masses for a selected sample of 110 ATLASGAL sources through SED fitting of their mid-infrared to submillimeter flux densities, between 8 and 870 µm, that is, a slightly narrower range than ours (3.6 µm–1.3 mm). We performed aperture photometry on the 2.5″ PPMAP-derived bolometric luminosity and column density maps using the exact same apertures. The König et al. (2017) apertures were designed to ensure consistent flux extraction over the same area from the mid-infrared to the submillimeter range. Throughout the comparison we introduced a systematic correction to account for variations in the adopted distances for the ALMA-IMF regions². The comparison is made on a subset of 10 regions observed by both studies: G008.67, G010.62, G012.80, G327.29, G333.60, G337.92, G338.92, G351.77, G353.41, W43-MM1, and W51-E.

On the one hand, our comparison (cf. Fig. 8) indicates that our bolometric luminosity estimates are generally consistent with the measurements made by König et al. (2017). The only noticeable discrepancy, exceeding a standard deviation, is observed for G333.60. This discrepancy is likely due to the saturation of Spitzer/MIPS observations in our study, a limitation that was circumvented by König et al. (2017) using MSX observations (Egan et al. 2003). Incorporating MSX observations into our study was not feasible, as our focus is to achieve the most accurate representation of the finer scales in the images (below 5 ). This objective contrasts with the relatively coarse angular resolution of MSX (18″), and the fact that MSX observations can only be added to the SED model a posteriori, since PPMAP cannot reproduce the emission from out-of-equilibrium dust grains.

On the other hand, our mass estimates exhibit more significant discrepancies with König et al. (2017)’s results (see Fig. 8). Because the derived mass depends on the estimated temperature, using the MBB description can yield larger discrepancies in the mass estimates compared to the luminosity estimates. Measuring luminosities involves a robust and straightforward measurement of the area below the data points, whereas the derived mass depends on the estimated temperature. As a consequence, it is anticipated that a more substantial variability may arise in mass, in contrast to luminosity (see Fig. 8). Furthermore, Galliano et al. (2018) reported that mass estimates may depend on the spatial resolution, since the temperature structure can be hidden in poorly resolved images, while it is accounted for in higherresolution observations (e.g., Fig. 14 in Aniano et al. 2012). This interpretation is consistent with the fact the more distant proto-clusters display larger mass discrepancies, while the less distant protoclusters better align with the identity curve, with the only exception of G327.29.

Additionally, our use of 8 dust temperatures for reproducing the observations contrasts with König et al. (2017) use of at most two temperatures. When using more temperature components, a small amount of warm dust can largely contribute to the 70 µm emission, thereby recovering an amount of cold dust that would otherwise would be missed, and consequently increasing the total mass (e.g., Aniano et al. 2012). Moreover, an other effect is the fact that the longest wavelength used by König et al. (2017) is 870 µm, whereas we also used 1.3 mm observations. We checked for these two effects and found that reducing the number of temperature components to two does result in higher density-weighted temperatures, while removing the 1.3 mm data generally results in recovering less mass. Taking into account these two effects combined, 14 out of 15 protoclusters masses are consistent with König et al. (2017)’s measurement within 1σ (with the exception of G338.93, that is consistent within 2σ). Finally, the use of a slightly different opacity index (β = 1.75) and larger reference opacity (κ₃₀₀ ≃ 0.125 cm²g⁻¹) by König et al. (2017) marginally accounts for these mass discrepancies.

Fig. 8 Accuracy of PPMAP measurements. Comparison of estimates made by PPMAP (see Table 2) and König et al. (2017) for the bolometric luminosity $L_{\text{bol}}^{\prime }$ (top panel) and mass M′ (bottom panel), integrated in the aperture defined by König et al. (2017) for 11 of the ALMA-IMF regions covered by both studies. The red and gray shades respectively indicate a 50% and 100% uncertainty range around the identity curve, represented by a solid black line.

4.3 Uncertainties

4.3.1 Description

We estimate here the uncertainties inherent to the PPMAP-derived products and measurements. While the MBB description itself constitutes an approximation, our results are influenced by additional complexities. Table 3 provides an estimate of the mean values of these uncertainties across the 15 protoclusters studied. The primary sources of errors within the PPMAP process, that impact the determination of luminosity, mass, column density and dust temperatures, are as follows:

1. Saturated pixels: saturation in the continuum observations used as PPMAP inputs. Most importantly, the saturation of the Spitzer/MIPS map at 24 µm (see Fig. 2) may impact the bolometric luminosity estimate, like suggested for that of G333.60 in Fig. 8. In contrast, the saturation of far-infrared and submillimeter images and its effect on PPMAP products are mitigated by the Gaussian process regression we applied (cf. Appendix A.1.1).

2. Free-free emission: the MBB description does not account for the free-free emission that contributes to the ALMA 1.3 mm image of evolved and intermediate regions. We used images approximately corrected for contamination by free-free emission, as described in Appendix A.1.3.

3. Noise estimates in the FIR to millimeter maps: they determine the uncertainty in estimating the measurement error for the input maps, that in turn determines the relative weights of data points used throughout the PPMAP SED fitting process. All the PPMAP results are thus sensitive to the methods used to determine the noise level of input maps (Appendix A.2.1).

4. PPMAP SED fitting: the uncertainties associated with the PPMAP fitting process for determining dust parameters. Includes systematic errors arising from the adopted opacity index (β = 1.8 ± 0.2).

5. Correction of the optically thick emission: while this correction is crucial to estimate the dust temperature maps of ALMA-IMF protoclusters, it relies on a model of extinction that introduces uncertainties (cf. Sect. 3.4).

Finally, we identified four minor sources of error.

1. Pointing errors inherent to any map sets taken with different observatories, here leading to relative shifts between Herschel, APEX, SOFIA, and ALMA data, could bias the SED fitting (cf. Table 3).

2. Uncertainties on the distance to the Sun of ALMA-IMF protoclusters (see Table 2) lead to errors on their luminosity and mass images.

3. Errors introduced by the splines interpolation of the near-infrared observations, below 70 µm.

4. Finally, ring-like artifacts around sources are expected to have an impact on the PPMAP-derived measurements (cf. Appendices A.2.5 and A.4 for a description of these artifacts).

4.3.2 Quantification of errors

We here quantify the mean uncertainty of each sources of error mentioned above. Systematic and random errors are listed and quantified in Table 3. The methods we employed to estimate the errors are the following:

1. To estimate the uncertainties caused by the saturation, free-free emission, noise estimates, opacity index, artifacts, and pointing errors, we ran PPMAP with modified input parameters and maps. We then inferred the errors from the discrepancies between the outcomes. For instance, the uncertainty on the column density originating from the saturation of input maps is defined as $\sigma \left(N_{{\text{H}}_2}\right)=\Delta N_{{\text{H}}_2}\left({\text{RUN}}_{\text{a}},{\text{RUN}}_{\text{b}}\right)=\left|N_{{\text{H}}_2}\left({\text{RUN}}_{\text{a}}\right)-N_{{\text{H}}_2}\left({\text{RUN}}_{\text{b}}\right)\right|$, where $N_{{\text{H}}_2}\left({\text{RUN}}_{\text{a}}\right)$ and $N_{{\text{H}}_2}\left({\text{RUN}}_{\text{b}}\right)$ are respectively the column density maps obtained with and without including the saturated images (see Fig. 2). The characteristics of the different PPMAP runs performed to derive uncertainties are described below Table 3.

2. For the uncertainty inherent to the PPMAP SED fitting process, we used the uncertainty output, “sigtdens.fits” (in which the random errors obtained from the SED fitting are stored), to derive the column density and temperature uncertainties following Eqs. (7) and (8).

We found that the uncertainty caused by potential Herschel pointing errors are negligible (<1%). Meanwhile, the primary contributors to uncertainties in constraining the dust temperature include the errors associated with PPMAP’s SED fitting, noise estimates, the choice of the opacity index (β), and the influence of ring-like artifacts. Uncertainties are in fact variable across the field of observations, thus we provide uncertainty maps corresponding to the relevant data products. Values in Table 3 are an account of the spatially averaged uncertainties (measuring the mean value across the error maps). The determination of the total uncertainties for column density ($N_{{\text{H}}_2}$ ) and temperature (T_dust) employs the combined standard uncertainty $u_{\text{tot}}^2={\displaystyle {\sum }_i{u}_i^2}$, that is, a quadratic sum over the uncertainties listed in Table 3), where u pertains to the errors enumerated above. The resulting total uncertainties are $u_{N_{{\text{H}}_2}}=0.29N_{{\text{H}}_2}$ and $u_{T_{\text{dust }}}=0.16T_{\text{dust }}\pm 1 \text{K}$. It should be noted that these total uncertainties do not include i.) the uncertainty on κ₃₀₀, due to the unknown composition and size distribution of dust (Köhler et al. 2015, Ysard et al. 2019, Schirmer et al. 2020); ii.) the bias induced by the PPMAP assumption of optical thinness, that may significantly affect the temperature and mass estimates, in particular toward high column density pixels. These considerations, along with the beam dilution bias mentioned earlier, should be treated as additional uncertainties if their relevance arises in the context of using our PPMAP estimates. Finally, we propagated the quadratic sum of errors for parameters such as β N, and T_dust to infer errors for luminosity and mass estimates (as defined by Eq. (9)), employing the same combined standard uncertainty approach. Variable errors associated with distances, as outlined in the second column of Table 2, are also factored into the mass and luminosity uncertainties.

5 Discussion

The PPMAP products presented in this paper have the potential for further analyses. Firstly, high-resolution dust temperature maps are an essential prerequisite for deriving core masses and constructing core mass functions (CMFs). Although the PPMAP 2.5″ beam is roughly five-fold larger than that of ALMA observations, our temperature maps currently offer the most comprehensive coverage and the best resolution available for the ALMA-IMF survey. In fact, ongoing studies by Louvet et al. (2023) and Armante et al. (2024) are employing these PPMAP-derived temperature maps for CMF investigations. Moreover, our column density maps provide a means to characterize the structure of the protoclusters (as illustrated in Fig. 6), and could be used along with different molecular tracers (e.g., N₂H⁺) to derive abundances in different parts of the same protocluster, or between protoclusters, that has a potential use as an evolutionary indicator. Finally, the luminosity-to-mass ratio may also constitute a tracer of the evolutionary stage of these regions, even enabling the discernment of subregions within the ALMAIMF fields. These tools open up new possibilities for further exploration within the ALMA-IMF survey, a topic we discuss in the following sections, featuring a selection of examples. More comprehensive analyses will be the focus of future studies.

5.1 PPMAP-derived column density

5.1.1 Probability density functions

Figure 9 presents the mean probability density function (hereafter PDF) of the PPMAP-derived column density across the 15 regions studied. Column density PDFs of molecular cloud are well described by a lognormal distribution in addition to a power-law tail (Kainulainen et al. 2009; Schneider et al. 2015, 2022, refer to Pouteau et al. 2022 for an analysis of column density PDFs in W43-MM2&MM3). The low-density tail is always limited by the included sky area, therefore the lognormal shape and position may be biased by the relatively small size of ALMA fields (typically 1′ × 1′). Even though they span a wide range in column densities (from 10²¹ cm⁻² to 10²⁵ cm⁻²), distances to the Sun and evolutionary stages, the cumulative PDF of the 15 regions shown in the left panel of Fig. 9 can be roughly described by a lognormal distribution, although substructures are seen. As an example of an individual PDF measured over the extent of the ALMA footprint, we show the PDF of the evolved G012.80 protocluster in the right panel of Fig. 9. In this individual case the best-fit lognormal distribution better matches the measurements, and for higher column densities the departure from the lognormal distribution can be clearly defined. Deviations from the lognormal shape are predicted for gas structures governed by self-gravity (Schneider et al. 2015 and references therein). The flattening of the distribution at higher column densities is described by a power-law, $p\propto {\left(N_{{\text{H}}_2}\right)}^{-s}$. We measure s = 2.1 ± 0.2 in G012.80, a value that is consistent with gravitational collapse of an isothermal sphere (Schneider et al. 2015). A comprehensive analysis of the column density PDFs is not in the scope of this study, we refer to Díaz-González et al. (2023) for a systematic, high-resolution study of column density PDFs inferred with PPMAP temperature maps.

Fig. 9 Probability density functions (PDFs) of the PPMAP-derived column density, normalized with respect to the area. The cumulative PDF across the 15 ALMA-IMF regions is shown on the left, and the PDF measured in the evolved G012.80 protocluster is shown on the right. Solid black lines represent the lognormal distribution that best fits the PDFs, while the dashed black lines correspond to the power-law tail, with the power-law index s indicated in the upper-right corner.

5.1.2 Comparison of the PPMAP column density maps with a N₂H⁺ line

From the comparison with König et al. (2017)’s measurements described in section 4.2, we established that the PPMAP products, including the column density and temperature maps, are robust in terms of large of the large scale measurements. These parameters ($N_{{\text{H}}_2}$, T_dust) are vital to understanding the chemistry in massive star-forming protoclusters (column density to measure abundances, temperature in relation to potential energy barriers), and to constraining dust simulations. However, the comparison presented in Sect. 4.2 pertains to mean values measured within large apertures (on the order of 10″). To assess our results at the precise angular resolution enabled by PPMAP (2.5″), we made further comparisons with observations of comparable resolution, such as the ALMA-IMF spectral data (Cunningham et al. 2023). Through this comparison we aim to check that small features (~2.5″) found in the ALMA high-resolution data are reproduced to some extent in the PPMAP products.

Figure 10 compares our column density maps with the N₂H⁺ J = 1–0 integrated line emission (Stutz et al. in prep; Alvarez-Gutierrez et al. in prep.), a tracer of the dense and cold medium (Pety et al. 2017) that has the potential to correlate with dusty filaments. Our analysis across regions in different evolutionary stages (G327.29: young, G353.41: intermediate, G012.80: evolved) reveals a general consistency between the PPMAP-derived column density and the N₂H⁺ integrated intensity map. Toward G327.29 and G353.41, both the global filamentary morphology and a fraction of the local emission peaks peaks are coherent between the dust column density and N₂H⁺ maps. Particularly remarkable is the image in the right panel of Fig. 10, displaying W33 Main-West filament (Immer et al. 2014; Armante et al. 2024). Along the filament, faint column density peaks align with local maxima of N₂H⁺ emission, while N₂H⁺ fades toward the central, higher column density peak. Local protostellar heating could account for the absence of N₂H⁺ emission within the central source, where a heightened gas temperature must result in its chemical destruction following the desorption of CO from grain mantles (Lee et al. 2004; Busquet et al. 2011; Sanhueza et al. 2012). A hot core was detected at this location (Armante et al. 2024; Bonfand et al. 2024), and the PPMAP-derived temperature does register a local increase in this specific area, up to 37 K (see Fig. 10), a measurement that may be consistent with a localized temperature of a hundred Kelvin, if we account for beam dilution (Motte et al., in prep.). Consistent associations between the positions of hot cores and local temperature increases in the dust temperature maps strengthen the case that the PPMAP products are reliable at their native angular scale of 2.5″, and demonstrate that they offer opportunities for interpretations of the physical and chemical mechanisms at work in the ALMA-IMF fields of observations.

5.2 Luminosity and mass of PPMAP luminosity peaks

Here we discuss the (bolometric) luminosity-to-mass ratio measured over the full extent of the protoclusters, before we delve into the luminosity-to-mass ratio of smaller sources mapped at the 2.5″ angular resolution. The bolometric luminosity is directly measured from the bolometric luminosity maps described in Sect. 3.3 (hence, including free-free emission), and the total gas mass is derived from the H₂ column density maps. Figure 11 shows the pixel-per-pixel PDFs of the luminosity-to-mass ratio partitioned with respect to the evolutionary stage proposed by Motte et al. (2022).

As star-forming protoclusters evolve, the contribution of HII regions to the luminosity is expected to gradually increase, thus the luminosity-to-mass ratio should be enhanced in the more evolved regions. This trend is found, with the respective distributions shifting toward a higher L/M ratio across the “young”, “intermediate” and “evolved” protoclusters, although a significant dispersion is measured: for any pair of regions the mean values of L/M are consistent with each other at a ±1σ level, where σ is the standard deviation of the distribution. The large dispersion may be primarily attributed to intra-region variability: regions within the ALMA-IMF survey, despite their relatively small sizes of a few parsecs, can encompass subre-gions with differing characteristics. This may result in a blend of lower and higher L/M ratios across the field of observations, as we observe in Fig. 11. This interpretation is reinforced by the fact that evolved regions display the largest dispersions. Indeed, evolved protoclusters may harbor a combination of young and more evolved subregions, pertaining to inhomogeneous initial conditions. Consequently, the spatially averaged L/M ratio (over the extent of the ALMA footprint) may not systematically serve as a robust indicator of the global evolutionary stage of one protocluster.

In Fig. 12, we illustrate the pixel-per-pixel spatial variations of the luminosity-to-mass ratio, specifically in the evolved G012.80 protocluster. Through a luminosity map clustering approach based on an arbitrary threshold (log(L/M) = 1.8), we discern a broad correlation between areas exhibiting higher L/M ratios and HII regions, that are traced by H41α and NeII line emissions (Beilis et al. 2022). This illustrates the capability to map luminosity-to-mass ratio variations at a 2.5″ angular resolution and to constrain the nature and evolutionary stage of resolved structures. In the subsequent subsection, we outline a method for estimating the luminosity-to-mas s ratios of resolved luminosity peaks that may be extracted as individual sources.

Fig. 10 PPMAP column density maps compared with the N2H+ integrated intensity map (overlaid in white contours), for three specific regions (from left to right: G327.29, G353.41, and G012.80). Contour levels: logarithmically spaced between 50–225 Κ km s−1 (left panel), 50–150 Κ km s−1 (center panel), and 50–150 Κ km s−1 (right panel).

Fig. 11 Probability density functions (PDFs) of the luminosity-to-mass ratio, normalized with respect to the area. Cumulative PDFs are shown for the evolved (top panel: G010.62, W51-IRS2, G012.80, G333.60), intermediate (central panel: G351.77, G008.67, W43-MM3, W51-E, G353.41) and young (bottom: W43-MM1, W43-MM2, G338.93, G328.25, G337.92, G327.29) protoclusters. Black markers and horizontal bars represent the median, first and last decile L/M measurements within individual protoclusters.

5.2.1 Source extraction with getsf

Here, we aim to compile a systematic catalog encompassing the luminosity and mass estimates of luminosity peaks found in the PPMAP dataset. With a limited angular resolution of 2.5″, it is conceivable that certain sources may not correspond to single entities, but rather to clusters of luminous sources, that could include a mixture of protostars, pre-stellar cores, and ultra-or hyper-compact HII regions. Impending studies will attempt to perform eros s-identifications based on already established core and protostar catalogs (Motte et al., in prep.). To identify sources within the luminosity maps, we employed the getsf method described by Men’shchikov (2021). We direct interested readers to that publication for a detailed exposition of the procedure. Below, the underlying principles and application criteria pertinent to our dataset are briefly summarized.

The getsf method entails a spatial deconstruction of observed images, effectively separating structural constituents and their background. The technique aims to parse distinct spatial scales and segregating sources and filaments from both one another and the background. Characterized by a single parameter, namely an approximate maximum size of sources to be extracted, detection yields initial approximations of source footprints, dimensions, and fluxes. Subsequently, more precise measurements of the source sizes and fluxes are conducted on background-subtracted images and, if warranted, on auxiliary images.

We executed the getsf algorithm on the ΡΡΜΑΡ luminosity maps generated from Run2, that is, without taking into account the free-free subtraction performed by Galvan-Madrid et al. (in prep.). This choice is rooted in our goal to ensure the best representation of the bolometric luminosities of the ΡΡΜΑΡ luminosity peaks, by including the contribution of free-free emission. In addition, ΡΡΜΑΡ column density maps resulting from Run1 were simultaneously given to getsf as auxiliary data, in order to measure the mass of the sources. Our initial step involved a resampling of the luminosity maps targeted at achieving a three-pixel sampling of the PPMAP beam (as elaborated in Appendix A.5), because it improves the detection of luminosity peaks associated with protostars. Furthermore, we fixed the maximum source size to 5″, a value equivalent to twice the dimensions of the PPMAP beam.

Fig. 12 Luminosity-to-mass ratio unveiled by PPMAP. Top panel: we show the complete histogram of L/M across the 15 ALMA-IMF regions studied, separated into two samples by the equation log(L/M) = 1.8. Bottom panel: as an example, a decomposition of the evolved G012.80 protocluster’s luminosity map is performed. Pixels with higher luminosity-to-mass ratio (log(L/M) ≥ 1.8 are plotted with a red col-ormap, while their counterpart (log(L/M) ≤ 1.8) are shown with a blue colormap. Superimposed white contours illustrate the H41α line emission, that traces regions dominated by free-free emission (contour levels: logarithmically spaced between 0.025 and 0.075 Jy beam−1). Green contours illustrate the NeII line emission (contour levels: logarithmically spaced between 0.005 and 0.04 erg s−1 cm−2 sr−1).

5.2.2 Results of the getsf extraction

The outcomes of the getsf extraction process are presented in Table B.1. This table provides details including celestial coordinates, angular and spatial full width at half maximum (FWHM), background-subtracted luminosities and masses, as well as the luminosity-to-mass ratio of the PPMAP luminosity peaks. We estimated the completeness level of the catalog of 313 luminosity peaks to be ~60 L_⊙, with a tendency for a better completeness level in young regions (~30 L_⊙) than in evolved ones (~100 L_⊙). In some instances (18% of the peaks), the luminosity peaks have no counterparts in the column density maps, resulting in a mass below our detection threshold. These mass measurements are discarded and signalled by the symbol “-” (see Table B.1). We interpret the luminosity peaks with no massive counterpart as diffuse areas heated by evolved protostars or HII regions.

The spatial distribution of the PPMAP luminosity peaks is superimposed on the column density and luminosity maps, in Figs. 7 and B.1. Upon visual examination, we observe a general correspondence between the PPMAP luminosity peaks and features such as dusty filaments, HII regions, clusters of cores and protostars extracted from the ALMA images. Correlations are also observed between PPMAP luminosity peaks and individual cores and/or protostars extracted from the ALMA images. A detailed cross-analysis of the PPMAP source catalog and the getsf core catalog extracted from ALMA continuum images falls beyond the scope of this paper and will be covered in a subsequent work by Motte et al. (in prep.).

The relation between the bolometric luminosity and mass is a useful metric to constrain the evolutionary stage of star forming objects. The PPMAP source sample spans a considerable range of luminosity-to-mass ratios, across four orders of magnitude (from L/M ≃ 10⁻¹ to L/M ≃ 10³, in solar units). This range underscores the wide spectrum of physical conditions and object types encompassed within the sample, ranging from (bright) ultra- and hyper-compact HII regions to (faint) cold, massive star-forming cores. Due consideration must be given to the fact that a single PPMAP source may correspond to several blended objects, possibly of different nature, because of the limited angular resolution (2.5″).

In Fig. 13, we present the distribution of masses and luminosities for the PPMAP source catalog, and its comparison with evolutionary tracks from Motte et al. (2018a) and Duarte-Cabral et al. (2013). The path of individual objects within the M–L_bol diagram can be predicted by accretion models. From a given initial envelope mass, the mass is expected to decrease at a given rate through material accretion onto the central star, in addition to material ejection by molecular outflows. The luminosity, on the other hand, is a function of stellar mass, hence it grows over time. These mechanisms steer the progression within the M–L_bol diagram from the top-left extremity to the end of the evolutionary tracks, allowing to follow the time evolution of cores and protostars. The positions of luminosity peaks shown in Fig. 13 are generally consistent with the accretion models presented by Motte et al. (2018a) and Duarte-Cabral et al. (2013), although approximately 20% of the sources are above the 50 M_⊙ final stellar mass track. We interpret these massive sources as clusters of unresolved cores and/or protostars. Furthermore, we observe that although the distributions of luminosity peaks pertaining to evolved and young protoclusters are overlapping, their centroid diverge within the diagram, aligning with different segments of the evolutionary tracks. Indeed, the median values of the luminosity-to-mass ratio (L/M) of luminosity peaks demonstrate a trend across regions when classified based on their evolutionary stages, as outlined by Motte et al. (2022). Specifically, we measure mean L/M values of 253.5, 12.8, and 12.3 L_⊙/M_⊙ for the luminosity peaks of evolved, intermediate, and young regions, respectively (with a respective dispersion around the mean of 100.0, 4.2 and 8.0 L_⊙/M_⊙, estimated by the mean absolute deviation).

Fig. 13 Mass and luminosity distribution of PPMAP luminosity peaks extracted by getsf. Caveat: with a 2.5″ angular resolution and at a distance of 2–5.5 kpc, PPMAP sources may correspond to several cores and/or protostars. The size of the markers reflects the FWHM of the sources, while the color indicates the evolutionary stage of the proto-clusters they belong to. Typical error bars are shown in the bottom-right corner of the diagram. Solid black lines represent the evolutionary tracks from Motte et al. (2018a) and Duarte-Cabral et al. (2013) for final stellar masses of 2, 4, 8, 20 and 50 M⊙.

5.3 Summary & perspectives

Using the multiwavelength, multiresolution Bayesian algorithm PPMAP, we performed a SED analysis of the dust emission in the 15 ALMA-IMF protoclusters. Near-infrared to millimeter observations from 8 instruments were included in the MBB analysis, spanning angular resolutions from subarcsecond (ALMA) to 35.2″ (Herschel). Our results are the following:

1. We present new measurements of the bolometric luminosity, column density and dust temperature toward the 15 massive ALMA-IMF protoclusters (Table 2, Fig. 4). The PPMAP estimates are consistent with previous measurements at a coarser angular resolution (König et al. 2017, Fig. 8), and thus constitute a benchmarked mapping of the dust parameters at the best angular resolution currently attainable (2.5″).

2. We compared our column density and dust temperature maps with the continuum cores identified by Louvet et al. (2023), the hot cores identified by Bonfand et al. (2024) and the N₂H⁺ J=1−0 line, showing that the 2.5″ features found in the PPMAP products are consistent with sources and structures mapped at ALMA’s native resolution, 0.3″ –0.9″ (Figs. 7, B.1, 10).

3. The pixel-per-pixel analysis of the luminosity-to-mass ratio shows that more evolved regions have, on average, a larger luminosity-to-mass ratio, although intra-region variability is observed (Fig. 11). We show, with an example in the G012.80 protocluster, that subregions pertaining to different evolutionary stages can be separated by setting a luminosity-to-mass ratio threshold (Fig. 12).

4. Using getsf (Men’shchikov 2021) on the PPMAP bolometric luminosity maps, we have extracted a catalog of 313 PPMAP sources, with the associated bolometric luminosity and mass measurements (Table B.1). We find that the luminosity-to-mass ratio of PPMAP sources tends to be consistent with the evolutionary stage attributed to the different protoclus-ters (Motte et al. 2022, Galván-Madrid et al., in prep.). We compared the PPMAP sources with evolutionary tracks from Duarte-Cabral et al. (2013) and Motte et al. (2018a), although the comparison is biased by the fact that one PPMAP source may correspond to several unresolved cores and/or protostars (Fig. 13).

Numerous potential directions stem from our current findings. Our analysis resulted in the determination of luminosity-to-mass ratios across ALMA footprints, ATLASGAL sources’ footprints, and PPMAP luminosity peaks extracted via getsf. Our catalog of 313 luminosity peaks presents a resource for forthcoming investigations. The cross-identification of PPMAP luminosity peaks with HII regions offers a pathway to estimate their luminosities. The dataset of column density, temperature, and luminosity maps we provide opens doors to investigate relationships among these quantities. Furthermore, the PPMAP-derived column density, in conjunction with PPMAP source and core catalogs, lays the groundwork for examining core formation efficiency and its dependence with the local gas density at a 2.5″ angular resolution, following the original study by Louvet et al. (2014) in W43-MM1.

Finally, despite potential challenges related to associating single PPMAP source with multiple protostars, our catalog may help constraining the temperatures and luminosity-to-mass ratios of young stellar objects. Assuming that a robust cross-identification can be completed, the catalog enables a direct comparison with accretion models, based on the M_core versus L_bol evolutionary diagram (e.g., Duarte-Cabral et al. 2013). One such study, underway by Motte et al. (in prep.), focuses on cross-identification and further processing of PPMAP-derived temperatures to estimate core and protostar temperatures and luminosities.

Acknowledgements

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.01355.L, #2013.1.01365.S, and #2015.1.01273.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This project has received funding from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130), from the French Agence Nationale de la Recherche (ANR) through the project COSMHIC (ANR-20-CE31-0009), and the French Programme National de Physique Stellaire and Physique et Chimie du Milieu Interstellaire (PNPS and PCMI) of CNRS/INSU (with INC/INP/IN2P3). M.B. is a postdoctoral fellow in the University of Virginia’s VICO collaboration and is funded by grants from the NASA Astrophysics Theory Program (grant number 80NSSC18K0558) and the NSF Astronomy & Astrophysics program (grant number 2206516). A.S. gratefully acknowledges support by the Fondecyt Regular (project code 1220610), and ANID BASAL project FB210003. A.S. gratefully acknowledges support by the Fondecyt Regular (project code 1220610), and ANID BASAL project FB210003. R.G.M. and D.D.G. acknowledge support from UNAM-PAPIIT project IN108822 and from CONACyT Ciencia de Frontera project ID 86372.

Appendix A Technical insights into the PPMAP implementation

In this appendix, we provide a comprehensive overview of the systematic steps employed while using the PPMAP software for the data analysis. The following sections outline the process, starting from the preparation of input continuum maps and point spread functions (PSFs) to the post-correction procedures applied to the PPMAP data products. Fig A.1 represents a flow chart of these procedures.

A.1 Preprocessing

A.1.1 Correction of saturated pixels

Several Herschel/PACS and Herschel/SPIRE maps exhibit saturation in the regions around the brightest sources, particularly within the peak of the spectral energy distribution (around 160250 µm). As a result, some pixels in these regions contain “NaN” values, rendering the use of PPMAP impossible for these affected maps. To address this issue, we took the following steps:

1. In severe cases of saturation (11 out of 124 images in our study), where a significant portion of the map was affected (≥ 25%) and the data quality was compromised, we chose to entirely remove these maps from the pool of input data.

2. In cases where the saturation level was lower (less than 25% of the map is saturated), we applied interpolation to estimate the missing values. We used the astrofix Python package (Zhang & Brandt 2021), that employs Gaussian process regression for this purpose. In this process, the algorithm used the unsaturated areas of the available input images as a training set to learn and determine the missing values, using an optimized Gaussian interpolation kernel (as shown in Fig. 2, 25 out of 124 images underwent interpolation for the present study, with a median of 6 interpolated pixels per image).

Figure 2 illustrates the extent of the pixels replaced using interpolation in the affected Herschel maps. We acknowledge that interpolating missing pixels, especially in regions of high intensity peaks, can be an uncertain procedure. To quantify the uncertainties associated with this treatment, we performed two iterations of the PPMAP analysis: one with the interpolated maps and the other without any interpolation (excluding the saturated maps from the analysis). By comparing the results between these two iterations, we measured the discrepancies introduced by the interpolation process. These discrepancies were then propagated into the final uncertainty maps released with this study (see Sect. 4.3.2 and Table 3).

A.1.2 Preparation of point spread functions

To combine a set of input images with different angular resolutions, PPMAP requires accurate Point Spread Function (PSF) profiles. The current PPMAP version³ is shipped with Herschel PSFs projected on 256 × 256 frames. We replaced these files with new ones projected on a larger 2048 × 2048 frame. This adjustment was necessary to ensure correct sampling of both the higher-resolution (ALMA) and lower-resolution (Herschel) PSFs while using a single array size for all PSFs, as required by PPMAP. The PSFs we used for Herschel/SPIRE (at 250,350, and 500 µm) were downloaded from the Herschel science archive⁴, while those for Herschel/PACS (at 70 and 160 µm) were obtained from Bocchio et al. (2016). Traficante et al. (2011) reported that the beams in the Hi-GAL images are elongated because of the sampling procedure adopted in the Galactic plane survey, with an ellipticity lower than 15% across all bands. We produced our own synthetic beam profile for SABOCA, LABOCA, and ALMA observations using the Gaussian2DKernel function from the astropy.convolution package in Python.

A.1.3 Free-free emission in the ALMA 1.3 mm emission map

As introduced in Sect. 2, to mitigate the contamination by freefree emission, we have excluded the ALMA Band 3 data from the PPMAP inputs. In addition, we employed two versions of ALMA Band 6 observations : the “native” version, and a free-free subtracted version, following a process detailed by Galván-Madrid et al. (in prep.). The subtraction method relied on the H41α line emission to disentangle the thermal emission of dust grains from the free-free emission. This was achieved by assuming an electron temperature of 7000 K and variable helium abundances, depending on the specific region under consideration. Throughout the PPMAP analysis, we either used the free-free subtracted maps to obtain the dust temperature and column density maps (refered to as the PPMAP “Run1”), or the native ALMA Band 6 maps to estimate the luminosity (refered to as the PPMAP “Run2”). The motivations for these two runs are described in Sect. 3.2. In Fig. A.2, we emphasize the distinctions in the PPMAP-derived column density between the Run1 (without free-free emission) and Run2 (with free-free emission) outcomes. While the overall column density within the ALMA footprint registers a mere 12.5% increment between the two instances, the column density estimate is significantly enhanced toward localized areas (up to a factor 20). This increase is particularly pronounced around the G012.80 Main-Central and Main-North HII regions as outlined by Armante et al. (2024). Conversely, the impact of the free-free removal process on the Main-West filament is minimal (relative difference below 10%), reinforcing its classification as a primarily dusty filament. The disparity between Run1 and Run2 underscores that the free-free subtraction is necessary to precisely estimate the column density in evolved regions.

A.2 PPMAP SED fitting

A.2.1 Determination of measurement errors

As discussed in Sect. 3.1, the PPMAP Bayesian fitting procedure necessitates accurate error estimates, in order to attribute a correct weight to the input observations. We explored three distinct methods to determine the noise:

1. First method (applied to ALMA, LABOCA, SABOCA observations): Assuming a Gaussian distribution for the sky noise, we conducted a Gaussian fit of the data histograms to extract the standard deviation (referred to as σ₁).

2. Second method (applied to Herschel observations): due to the high signal-to-noise sky background associated with Herschel maps, the first method was not applicable. Instead, we employed the built-in PPMAP noise estimate function, getnoise, by setting the sigobs input parameter to 0 in the premap routine (Marsh et al. 2015). Using Gaussian kernels of variable scales, this function facilitates the measurement of sky noise through iterative subtraction of a smoothed background (referred to as σ₂).

3. Finally, in a third attempt to estimate the noise, we calculated the median absolute deviation (referred to as “MAD”, and denoted as σ₃) in empty areas of the observed fields (searching for such empty areas in a 30′ radius around the Herschel field center coordinates). MAD is a robust statistic that provides an estimate of the image noise, based on the median of the absolute deviations from the median pixel value in the image (Leys et al. 2013).

Fig. A.1 Flow chart of the data preprocessing, PPMAP calculations, and products’ post-processing steps applied in this study. “Postcorrection” and “Post-proc.” respectively refer to the treatments described in Appendix A.3 and Appendix A.4. “ASTROFIX” is the Gaussian process regression algorithm described in Appendix A.1.1. “VSG” stands for very small grains. The light blue boxes correspond to the final products released with this study.

Based on our investigation, when applied to ALMA, LABOCA, and SABOCA observations, all three noise measurement methods exhibited a general agreement, with relative differences within the range of 30-60% when compared to σ₂ (the noise estimate obtained from the PPMAP built-in subroutine). The specific values of these relative differences depended on the regions being analyzed. We always observed that σ₂ > σ₁ > σ₃, indicating that the third method consistently yields lower noise estimates compared to the other two methods. This outcome is particularly beneficial since it effectively increases the weight of the 1.3 mm image, thus allowing a better representation of the low signal-to-noise features in the observed protoclusters (such as faint cores).

Therefore, we incorporated the median absolute deviation measurement, along with the systematic error on the absolute flux calibration of the sky background $F_{\lambda }^{\text{sky }}$, to estimate the total errors for each input map at a given wavelength λ (e.g., σ_70µm, σ_350µm, σ_870µm, σ_1.3mm, etc). The error estimates σ_λ were computed as follows: $${\sigma }_{\lambda }=\mathrm{MAD}\left(F_{\lambda }\right)+u_{\lambda }{F}_{\lambda }^{\text{sky }},$$(A.1)

where, MAD(F_λ) represents the median absolute deviation measured in the map F_λ, and u_λ denotes the absolute flux calibration error (e.g., 15% for LABOCA, as discussed in Sect. 2). The resulting values of σ_λ were subsequently assigned to the input parameter sigobs for the corresponding wavelengths. The relative contribution of the two terms in Eq. (A.1) varies significantly depending on the map. In typical low signal-to-noise pixels (dominated by Gaussian noise in the ALMA observations), the error estimates for Herschel maps are primarily dominated by the sky background component, while the MAD component is the largest contributor to the ALMA Band 6 errors.

A.2.2 Herschel color correction

We applied color corrections to the Herschel/PACS and SPIRE observations (see Balog et al. 2014; Valtchanov et al. 2018 and references therein). These corrections are taken directly from the PACS and SPIRE user manuals, and consider the observed flux in each band as an integral over the product of the instrumental bandpass and the source’s spectrum, assuming a MBB profile with β = 2. The correction tables used in our analysis are consistent with the correction factors used in Motte et al. (2018b) for the PPMAP analysis of W43-MM1. These color corrections are provided to PPMAP as input tables and then automatically applied during the analysis. While it is ideally more accurate to include the color corrections, their effect on the PPMAP outcome is minimal (< 1% in relative variation of the column density and temperature estimates) and can be considered negligible in our analysis, with respect to other sources of uncertainties.

Fig. A.2 Demonstration of the impact of free-free subtraction, carried out by Galván-Madrid et al. (in prep.), on the PPMAP-derived column density (top panel) and luminosity (bottom panel) maps of G012.80. The left panel displays the map obtained from the subtracted data (Run1), while the right panel depicts the outcomes achieved with the standard ALMA data (Run2). Superimposed white contours illustrate the H41α line emission, that traces regions dominated by free-free emission (contour levels: logarithmically spaced between 0.045 to 0.5 Jy beam−1). Arrows indicate the positions of the main HII regions.

A.2.3 PPMAP parameters

The fixed PPMAP parameters adopted in our initial analysis are shown in Table A.1. Following previous applications of PPMAP (e.g., Marsh et al. 2017; Howard et al. 2019; Whit-worth et al. 2019; Chawner et al. 2020; Howard et al. 2021), we assumed that we can reproduce the observations in the FIR to submm (70 µm to 1.3 mm) with 8 MBB components between 10 K and 50 K (using a logarithmic spacing). We have systematically checked that allowing PPMAP to employ lower or higher temperature components does not provide a substantial change in column density $N_{{\text{H}}_2}$ and luminosity (within ± 5%). Ground-based input images inherently suppress low spatial frequencies. To address this, we enabled the input parameter highpass (see note on Table A.1 footnote (d)), allowing PPMAP to subtract a constant background from the internal model images. This subtraction is applied at wavelengths corresponding to the ground-based data (SABOCA, LABOCA, and ALMA). The goal is to enhance the alignment between the model and the spatial characteristics of these particular observations, ultimately improving the fitting process. We adopted a standard reference opacity κ₃₀₀ = 0.1 cm²g⁻¹, defined with respect to the total mass and corresponding to a gas-to-dust ratio of 100 (Marsh et al. 2015; Hildebrand 1983). Based on previous studies (Planck Collaboration XI 2014; Juvela et al. 2015; Köhler et al. 2015 and references therein for more diffuse gas), in a first approach we fixed the dust opacity to β = 1.8 for all regions. If we were to assume β = 2, the estimated dust temperatures would likely be lower (see Table 3). Finally, to balance between our highresolution requirements and potential artifact-inducing effects, we adopted a synthetic angular resolution of 2.5″ for the analysis. In fact, extensive testing of PPMAP indicated that further reducing the angular resolution could lead to stronger ring-like artifacts (see Appendix A.4).

A.2.4 Temperature and opacity index optimization

Based on our first series of PPMAP runs, performed with a fixed temperature range (Tdust = 10 – 50 K, using 8 MBB components) and fixed opacity index (β = 1.8), we refined the input temperature range. From the temperature distributions measured in the zero-order PPMAP output, we systematically lowered the amplitude of the temperature range (e.g., from T_dust = 10 – 50 K to T_dust = 21 – 39 K in G012.80). This adjustment allowed us to obtain more precise temperature estimates for the sources, by decreasing the ∆T interval between the temperature components. Additionally, we explored the possibility of making the opacity index, β, a free parameter in the PPMAP analysis. To achieve this, we conducted iterative runs, modifying the value of β in increments of ∆β = 0.1, while measuring the corresponding change in chi-square (χ²). Following this procedure, we found a best-fit value of the opacity index, that could be different from β = 1.8 (between β = 1.4 and β = 2.1, depending on the specific fields of observation). However, we acknowledge that the degeneracy between the temperature and opacity index makes our attempt at simultaneously constraining them uncertain. The variations in X² associated with the modification of β were negligible compared to the expected χ² variation at a 1σ level. More precisely, for any β_i and β_j, we found that ∆x²(β_i,β_j) ≪ [(1 + 1/(N − p)]min(x²), where N is the number of data points used to fit the modified blackbody and p is the number of parameters. As a result, we decided to fix β = 1.8 for all regions, and focused on the determination of the temperature. We acknowledge the potential uncertainty in the opacity index, that may differ significantly from the assumed value of β = 1.8 across the ALMA-IMF protoclusters. To explore the impact of different opacity indices, we assessed variations in the PPMAP-derived temperature. For instance, we observed average changes of ∆T_dust = +3 K and ∆T_dust = +12 K when adopting β = 1 .5 and β = 1 (representative of more evolved dust in denser media), respectively.

A.2.5 Criteria to stop PPMAP iterations

The quality of PPMAP outputs underwent a systematic validation process through an a posteriori chi-square test, where we manually compared the synthetic fluxes obtained from the PPMAP model to the actual observed fluxes. Based on this a posteriori validation of the product’s reliability, we monitored the convergence of the PPMAP iterative process to prevent over-fitting, that can lead to increasingly severe ring-like artifacts in the column density and temperature maps (Sect. A.4). We refer to this particular effect as “over-fitting”, since the associated change in synthetic flux is minimal with respect to measurement errors. We systematically controlled for and prevented over-fitting through minor tweaks of the PPMAP iteration number, while validating the outcome through the aforementioned chi-square test. These adjustments were necessary to strike a balance and ensure the best representation of the observed data without introducing excessive artifacts. Most importantly, we note that although it does mitigate the growth of the spurious features, this validation process does not prevent their development.

A.3 Dust temperature post-correction

PPMAP relies on the assumption that the dust emission is optically thin. However, this assumption is not met in regions with high-density spikes, such as cores in high-mass star-forming regions. To account for the effect of larger optical depths, we used a procedure established in previous work (first applied by Motte et al. 2018b in the analysis of W43-MM1) to correct the dust temperature maps generated by PPMAP. The temperature correction table used in this procedure was constructed based on a grid of model sources with varying temperatures $(T_{\text{dust }}^{\text{model }})$, optical depths, sizes, and measurement errors. The synthetic sources were represented by the simplest model possible: a uniform dust slab with specified temperature and optical depth. Using this simple dust model, we produced synthetic flux density maps between 70 µm and 3 mm. Then, we applied PPMAP to the synthetic flux density maps to provide an estimate for the column density N(H₂) and temperature $T_{\text{dust }}^{\text{PPMAP }}$ Using different dust model parameters, we sampled the relationship $\Delta T_{\text{dust }}=f\left(N\left({\text{H}}_2\right),T_{\text{dust }}^{\text{PPAP }}\right)$ between the temperature bias $\Delta T_{\text{dust }}=T_{\text{dust }}^{\text{PPMAP}}-T_{\text{dust }}^{\text{model }}$, the PPMAP-derived column density N(H₂) and the PPMAP-derived dust temperature $(T_{\text{dust }}^{\text{model }})$. This investigation enabled the implementation of an a posteriori correction. The established relationship is depicted in Fig. A.3, covering the initial range of explored temperature and column densities. In fact, the original correction table provided by Motte et al. 2018b is constrained within certain limits, insufficient to accommodate our PPMAP products. To extend the correction beyond these limits, we performed an extrapolation using cubic splines while assessing the quality and reliability of the extrapolated results. Subsequently, we applied the extrapolated correction table to the PPMAP outputs a posteriori to account for the effect of optical depth on the dust temperature estimates. In each pixel of the PPMAP outputs, temperatures are updated using the tabulated value ∆T_dust = f (N(H₂), T_dust). The impact of the correction on the column density is considered negligible, since the added warm dust represents a minor contribution to the total mass. It is essential to note that this correction is based on a simplified model and may not be fully suitable for the diverse environments observed in the ALMA-IMF fields, including HII regions, internally heated protostellar envelopes and prestellar cores (Motte et al., in prep.). Despite its limitations, this correction helps compensate for the assumption of optical thinness of PPMAP and should be considered when aiming to constrain the temperature in regions with potential optical depth effects such as ALMA-IMF protoclusters.

Fig. A.3 Temperature correction ∆Tdust as a function of the PPMAP column density N(H2). The correction is shown for different input temperatures, spanning 20.5 K to 39.3 K. Data points are interpolated using cubic splines. N.B.: 99% of the PPMAP-derived column densities are below 8 × 1023 cm−2.

A.4 PPMAP ring-like artifacts and post-processing

Temperature and column density maps generated with PPMAP often exhibit ring-like artifacts, as shown in Fig. A.4. These artifacts consistently appear around sources of significant brightness, and their intensity grows with the brightness of the source they are associated with. They generally occur either as complete or partial spurious rings around a source, and in the most extreme case are accompanied by secondary peaks (cf. Fig. A.4, first panel). The underlying cause of these ring-like artifacts may be attributed to several inter-twinned factors. Firstly, PSF errors and missing spatial frequencies between the ALMA images and the rest of the observations could produce these artifacts. Secondly, due to the strong dependence of emission on temperature, PPMAP may fail to reproduce steep temperature gradients around bright sources, thus causing ring-like artifacts. Lastly, the violation of PPMAP’s assumption of optically thin emission toward spikes in column density may also contribute to errors (Howard 2020). As reported in Appendix A.2.5, elevating the number of PPMAP iterations beyond a certain threshold amplifies these spurious patterns. This emphasizes the importance of accurately assessing measurement errors and controlling the number of iterations to prevent over-fitting. Most importantly, the presence of ring-like artifacts may introduce a bias in determining the temperatures of cores. To gauge the extent of this effect, we quantified the discrepancy between our actual maps and PPMAP outputs resulting from an increased number of iterations (N_iter = 25 000). This yielded an estimated upper limit of ± 1 Kelvin for potential temperature biases stemming from these artifacts (cf. Table 3). This estimate of ±1 Kelvin effectively corresponds to the maximum value of the temperature variation measured between the crest and anti-crest of the artifacts in the temperature maps.

We implemented a conservative post-processing treatment to improve the quality of the PPMAP column density and luminosity products, specifically to reduce spurious features around high-density peaks. The post-processing method is defined as follows:

1. We first convolved the original map with a Gaussian kernel of size 1.25″ (equivalent to 1 pixel).

2. Next, we flagged pixels (i, j) that met two specific conditions, identifying them as potentially spurious and in need of correction: $$\begin{array}{l}(a)\frac{\left|P(i,j)-P^{\prime }(i,j)\right|}{P(i,j)}\ge t,\\ (b)\text{B}6(i,j)\le {\sigma }_{\text{MAD}},\end{array}$$(A.2)

   where P and P′ are respectively the PPMAP output (column density, temperature and luminosity maps) and its convolved counterpart, B6 is the ALMA 1.3 mm map smoothed at the PPMAP’s angular resolution of 2.5″, σ_MAD is the median absolute deviation measured in B6, and t = 0.05 is a suitable threshold.

3. The flagged pixels were then corrected by replacing their values with the corresponding values from the convolved product (P′(i, j)).

The post-processing approach smoothed out artifacts selectively, ensuring that only spurious pixels were modified. Pixels associated with statistically significant signal, as measured in the ALMA B6 continuum map, remained unchanged. The result of post-processing is shown in Fig. A.4 for a small sample of hand-selected artifacts.

Fig. A.4 Zooming in on characteristic ring-like artifacts, we present four column density maps (G012.80, G327.29, W43-MM2, and W43-MM3) in a left-to-right sequence. The top panel illustrates the initial PPMAP outputs, while the bottom panel displays the maps after undergoing post-processing as outlined in Appendix A.4. The color scale is logarithmic, to enhance visibility.

A.5 Resampling

Finally, we resampled the PPMAP column density and luminosity products onto new grids using the DeForest 2004 adaptive, anti-aliased resampling algorithm implemented through the reproject_adaptive routine from the Python package reproject. In this step, the sampling of the beam is increased from 2 pixels per beam to 3 pixels per beam (that is, the pixel size is lowered from 1.25″ to 0.83″). The main motivation for this increased sampling was to achieve better source detection by the getsf algorithm (see Sect. 5.2). The adaptive resampling process leads to minimal map variations, typically around 1–2% on average, and potentially a few Kelvin toward warm sources in the temperature maps. To address this limitation, we employed the reproject_exact routine on the PPMAP temperature maps. This approach enabled us to conserve temperatures at the positions of point-like sources, such as cores and protostars, albeit at the expense of introducing (negligible) aliasing artifacts.

Appendix B Supplementary tables and figures

Fig. B.1 PPMAP products illustrated for the remaining protoclusters. From left to right: column density map (N(H2)), bolometric luminosity (Lbol), corrected dust temperature (Tdust). White continuous contours outline the ALMA 1.3 mm mosaic areas. The luminosity peaks extracted from the PPMAP luminosity maps (see Sect. 5.2) are overlaid in gray. The continuum cores identified by Louvet et al. (2023) in the ALMA 1.3 mm images are overlaid in white. The size of the ellipses reflects the FWHM of the sources.

Fig. B.1 (Continued)

Fig. B.1 (Continued)

Fig. B.1 (Continued)

Fig. B.2 Performance of PPMAP modeling. Reduced χ2 maps are shown for each region studied. Black contours represent the continuum emission at 1 mm, using logarithmically spaced levels between 1% and 10% of the maximum flux.

References

All Tables

All Figures

Fig. 1 Observational constraints. Bottom panel: illustrative spectral energy distribution (SED), produced using a THEMIS grain mixture (Jones et al. 2017, gray curve). The blue curve represents a single modified blackbody that best fits the far-infrared and millimeter range of the SED. Black markers are overlaid on the gray curve to indicate the SED coverage enabled by the observations listed in Table 1, with the horizontal bars representing the bandwidths. Top panel: beam size of the observations used in our analysis, with the wavelength of the observations indicated on top of each marker (in micron, except for ALMA markers). Vertical bars represent the distance-induced variation in physical beam size across the ALMA-IMF sample.

In the text

Fig. 2 Data coverage chart. The color scale represents the percentage of observed and unsaturated pixels in each pair of region and map. Spitzer/IRAC data, which has no saturated pixels in the regions studied, is not shown here. The numbers account for the amount of pixels replaced through astrofix (example: in G012.80, the Hi-GAL 160 µm map has 6 saturated pixels that were interpolated). Maps that are either missing or discarded from the analysis are hatched.

In the text

Fig. 3 Schematic representation of the PPMAP iterative process. G012.80 images (at λ1 = 350 µm and λ2 = 870 µm) are used for illustrative purpose. In the upper part, we represent how PPMAP distributes “points” in a continuous state space (X, Y, Tdust) that can be divided into finite cells (corresponding to PPMAP pixels, that is, with a size fixed by the user, independent of the pixel size of the observed images). This distribution is then translated into a synthetic continuum emission map through the MBB description, taking into account the PSF of the instruments. Synthetic observations are finally compared to true observations, allowing to update the distribution of points. These iterative steps are repeated until the model converges.

In the text

Fig. 4 SEDs extracted from the ATLASGAL sources’ footprints (Contreras et al. 2013, Urquhart et al. 2014, see Sect. 3.3) corresponding to the protoclusters mapped by ALMA-IMF. The actual observations are represented by black points, while the red curve depicts the PPMAP MBB that provides the best fit to the data, with the gray shaded area representing the ±2σ standard deviation of the best fit. The blue curve represents the total piecewise model described in Sect. 3.3. Downward and upward arrows respectively correspond to lower limits and saturated observations.

In the text

	Fig. 5 Temperature correction of the PPMAP images illustrated. The left panels displays the W43-MM1 Main-West image before correction, while the right panels presents the post-corrected map. White ellipses outline continuum cores identified by Louvet et al. (2023) in the ALMA 1.3 mm images at 0.4–0.9″ angular resolution.
In the text

Fig. 6 Resolution enhancement of the dust temperature (top panels) and column density (bottom panels) images, achieved through the application of PPMAP to the G353.41 dataset. The left panels display the maps derived by smoothing all input images to a uniform resolution of 19.2″, while the right panels represent the PPMAP images at an angular resolution of 2.5″. White ellipses outline continuum cores identified by Louvet et al. (2023) in the ALMA 1.3 mm images at 0.4−0.9″ angular resolution. The larger dashed circles represent the footprint of the ATLASGAL source AGAL353.409–00.361 (Contreras et al. 2013, Urquhart et al. 2014).

In the text

Fig. 7 PPMAP products illustrated for three example regions: the young W43-MM1 (top), intermediate G008.67 (center), and evolved G012.80 (bottom) protoclusters. From left to right: column density map (N(H2)), bolometric luminosity (Lbol), dust temperature (Tdust). White continuous contours outline the ALMA 1.3 mm mosaic areas. The luminosity peaks extracted from the PPMAP luminosity maps (see Sect. 5.2) are overlaid in gray. The continuum cores identified by Louvet et al. (2023) in the ALMA 1.3 mm images are overlaid in white. The size of the ellipses reflects the FWHM of the sources.

In the text

Fig. 8 Accuracy of PPMAP measurements. Comparison of estimates made by PPMAP (see Table 2) and König et al. (2017) for the bolometric luminosity $L_{\text{bol}}^{\prime }$ (top panel) and mass M′ (bottom panel), integrated in the aperture defined by König et al. (2017) for 11 of the ALMA-IMF regions covered by both studies. The red and gray shades respectively indicate a 50% and 100% uncertainty range around the identity curve, represented by a solid black line.

In the text

Fig. 9 Probability density functions (PDFs) of the PPMAP-derived column density, normalized with respect to the area. The cumulative PDF across the 15 ALMA-IMF regions is shown on the left, and the PDF measured in the evolved G012.80 protocluster is shown on the right. Solid black lines represent the lognormal distribution that best fits the PDFs, while the dashed black lines correspond to the power-law tail, with the power-law index s indicated in the upper-right corner.

In the text

	Fig. 10 PPMAP column density maps compared with the N2H+ integrated intensity map (overlaid in white contours), for three specific regions (from left to right: G327.29, G353.41, and G012.80). Contour levels: logarithmically spaced between 50–225 Κ km s−1 (left panel), 50–150 Κ km s−1 (center panel), and 50–150 Κ km s−1 (right panel).
In the text

Fig. 11 Probability density functions (PDFs) of the luminosity-to-mass ratio, normalized with respect to the area. Cumulative PDFs are shown for the evolved (top panel: G010.62, W51-IRS2, G012.80, G333.60), intermediate (central panel: G351.77, G008.67, W43-MM3, W51-E, G353.41) and young (bottom: W43-MM1, W43-MM2, G338.93, G328.25, G337.92, G327.29) protoclusters. Black markers and horizontal bars represent the median, first and last decile L/M measurements within individual protoclusters.

In the text

Fig. 12 Luminosity-to-mass ratio unveiled by PPMAP. Top panel: we show the complete histogram of L/M across the 15 ALMA-IMF regions studied, separated into two samples by the equation log(L/M) = 1.8. Bottom panel: as an example, a decomposition of the evolved G012.80 protocluster’s luminosity map is performed. Pixels with higher luminosity-to-mass ratio (log(L/M) ≥ 1.8 are plotted with a red col-ormap, while their counterpart (log(L/M) ≤ 1.8) are shown with a blue colormap. Superimposed white contours illustrate the H41α line emission, that traces regions dominated by free-free emission (contour levels: logarithmically spaced between 0.025 and 0.075 Jy beam−1). Green contours illustrate the NeII line emission (contour levels: logarithmically spaced between 0.005 and 0.04 erg s−1 cm−2 sr−1).

In the text

Fig. 13 Mass and luminosity distribution of PPMAP luminosity peaks extracted by getsf. Caveat: with a 2.5″ angular resolution and at a distance of 2–5.5 kpc, PPMAP sources may correspond to several cores and/or protostars. The size of the markers reflects the FWHM of the sources, while the color indicates the evolutionary stage of the proto-clusters they belong to. Typical error bars are shown in the bottom-right corner of the diagram. Solid black lines represent the evolutionary tracks from Motte et al. (2018a) and Duarte-Cabral et al. (2013) for final stellar masses of 2, 4, 8, 20 and 50 M⊙.

In the text

Fig. A.1 Flow chart of the data preprocessing, PPMAP calculations, and products’ post-processing steps applied in this study. “Postcorrection” and “Post-proc.” respectively refer to the treatments described in Appendix A.3 and Appendix A.4. “ASTROFIX” is the Gaussian process regression algorithm described in Appendix A.1.1. “VSG” stands for very small grains. The light blue boxes correspond to the final products released with this study.

In the text

Fig. A.2 Demonstration of the impact of free-free subtraction, carried out by Galván-Madrid et al. (in prep.), on the PPMAP-derived column density (top panel) and luminosity (bottom panel) maps of G012.80. The left panel displays the map obtained from the subtracted data (Run1), while the right panel depicts the outcomes achieved with the standard ALMA data (Run2). Superimposed white contours illustrate the H41α line emission, that traces regions dominated by free-free emission (contour levels: logarithmically spaced between 0.045 to 0.5 Jy beam−1). Arrows indicate the positions of the main HII regions.

In the text

	Fig. A.3 Temperature correction ∆Tdust as a function of the PPMAP column density N(H2). The correction is shown for different input temperatures, spanning 20.5 K to 39.3 K. Data points are interpolated using cubic splines. N.B.: 99% of the PPMAP-derived column densities are below 8 × 1023 cm−2.
In the text

	Fig. A.4 Zooming in on characteristic ring-like artifacts, we present four column density maps (G012.80, G327.29, W43-MM2, and W43-MM3) in a left-to-right sequence. The top panel illustrates the initial PPMAP outputs, while the bottom panel displays the maps after undergoing post-processing as outlined in Appendix A.4. The color scale is logarithmic, to enhance visibility.
In the text

Fig. B.1 PPMAP products illustrated for the remaining protoclusters. From left to right: column density map (N(H2)), bolometric luminosity (Lbol), corrected dust temperature (Tdust). White continuous contours outline the ALMA 1.3 mm mosaic areas. The luminosity peaks extracted from the PPMAP luminosity maps (see Sect. 5.2) are overlaid in gray. The continuum cores identified by Louvet et al. (2023) in the ALMA 1.3 mm images are overlaid in white. The size of the ellipses reflects the FWHM of the sources.

In the text

	Fig. B.1 (Continued)
In the text

	Fig. B.1 (Continued)
In the text

	Fig. B.1 (Continued)
In the text

	Fig. B.2 Performance of PPMAP modeling. Reduced χ2 maps are shown for each region studied. Black contours represent the continuum emission at 1 mm, using logarithmically spaced levels between 1% and 10% of the maximum flux.
In the text