© E. Redaelli et al. 2021

1 Introduction

It is generally known that stars form from the fragmentation and subsequent collapse of the cold molecular phase of the interstellar medium (ISM). However, while the star formation process in the low-mass regime is fairly well understood, we still lack a comprehensive view of how high-mass stars (M > 8–10 M_⊙) are born. On the other hand, they play a key role in the energetic budget of the ISM. Furthermore, there is evidence that our Sun was born in a cluster also containing massive stars (Adams 2010; Pfalzner & Vincke 2020). For all these reasons, the study of high-mass star formation is one of the crucial and still unanswered questions of modern astrophysics.

In this context, several competing theories have been developed (McKee & Tan 2003; Bonnell & Bate 2006; Smith et al. 2009; Motte et al. 2018; Padoan et al. 2020). Important differences among these models concern the very early stages of the process, characterised by the formation and evolution of so-called high-mass pre-stellar cores (HMPCs). In particular, the different theoretical models make distinct predictions on the masses, accretion modes, and in general on the dynamical and physical initial stages of HMPCs. Observations targeting high-mass star formingclumps, aimed to investigate the mass distribution and kinematic structure of such objects, could therefore provide crucial constraints to the theory.

Such observations are, however, difficult to perform for several reasons. High-mass stars are intrinsically rarer and more short-lived with respect to their low-mass counterparts, as predicted by stellar evolution theories. As a consequence, they are on average more distant, which in turn affects the achievable spatial resolution of the observations. Moreover, they form in crowded and dense environments, which are heavily affected by extinction (e.g. Zinnecker & Yorke 2007). This means that the lack of infrared emission detected with single-dish facilities is not conclusive evidence of a pre-stellar stage (e.g. Motte et al. 2018). In this context, interferometers such as the Atacama Large Millimeter/submillimeter Array (ALMA) represent a powerful tool, capable of providing the necessary angular resolution and sensitivity to resolve substructures in distant high-mass clumps.

In the context of the search for HMPCs, targeting 70 μm-dark clumps withinterferometric studies has been a powerful tool (Sanhueza et al. 2017; Contreras et al. 2018; Li et al. 2019; Pillai et al. 2019). In Molet et al. (2019), the authors investigated two HMPC candidates in the rich, high-mass star forming regions W43-MM1, previously investigated by Nony et al. (2018), with ALMA Band 6 observations. One of the two cores does not show clear evidence of proto-stellar activity, but even the ALMA resolution is not sufficient to draw clear conclusions on its nature. Using Submillimeter Array (SMA) and Very Large Array (VLA) data, Cyganowski et al. (2014) studied the massive starless core G11.92-0.61-MM2 and found that, despite its high density (n ≳ 10⁹ cm⁻³) and mass (M ≳ 30 M_⊙), it lacks line emission in several abundant species such as N₂H⁺, HCO⁺, HCN (down to the sensitivity of their observations). Pillai et al. (2019) investigated two infrared dark clouds (IRDCs), selected to be 70 μm-dark and therefore believed to be pre-stellar, with the SMA. The authors found several substructures in the dust thermal emission data. However, using CO (2–1) observations, they found that both IRDCs host a large population of molecular outflows, indicating that several protostars are already active. The Alma Survey of 70 μm dark High-mass clumps in Early Stages survey (ASHES; Sanhueza et al. 2019; Li et al. 2020) targeted twelve IRDCs in the wavelengths 3.6–70 μm with ALMA at 1.3 mm continuum and with several molecular tracers. They identified ≈300 cores, and based on the detection of outflows and/or of warm transitions¹, ≈ 70% of them are classified as pre-stellar. None of them appear more massive than 10–30 M_⊙.

These examples highlight the difficulties not only in observing cores in the high-mass regime, but also in correctly classifying them as pre- or proto-stellar. Continuum observations at millimetre and submillimetre wavelengths (especially when only one frequency is available) cannot provide information on the mass and on the temperature independently, and they cannot unveil the source’s kinematics either. Other wavelengths could provide useful information in this sense, as done, for instance, with X-ray observations in Yu et al. (2020). On the other hand, chemistry is a powerful tool to probe the evolution of star forming regions. In the cold molecular phase of the ISM, at low temperatures (T ≲ 20 K) and high densities (n ≳ 10⁴ cm⁻³), most C- and O-bearing species are frozen out onto dust grains (Caselli et al. 1999; Bacmann et al. 2002; Giannetti et al. 2014; Sabatini et al. 2019), and they are therefore depleted from the gas phase. In these physical conditions, processes such as deuteration are greatly enhanced (Caselli & Ceccarelli 2012; Ceccarelli et al. 2014). In particular, the first deuterated species in the gas phase is produced by the following reaction: (1)

which is particularly efficient at low temperatures (assuming a low ortho-to-para H₂ ratio; see e.g. Pagani et al. 1992), due to the lower zero-point energy of H₂D⁺. From this, and from its doubly and triply deuterated forms D₂H⁺ and , all other deuterated molecules in the gas phase (such as N₂D⁺ and DCO⁺) are formed (Dalgarno & Lepp 1984).

As partof the search for HMPCs, Tan and collaborators used an extensive dataset of N₂D⁺ observations with ALMA. Tan et al. (2013) initially reported the detection of six cores, the most massive of which was revealed to be in fact composed of several cores, some of which were already in a proto-stellar stage (Tan et al. 2016; Kong et al. 2018). The sample was increased to 141 N₂D⁺-identified cores in Kong et al. (2017). Interestingly, only one core appears to be more massive than 10 M_⊙. However, the lack of tracers of proto-stellar activity leaves the question about the evolutionary stage of this core still open.

Deuteration of simple molecules has indeed raised huge interest recently, since it has been shown to be a useful ‘chemical clock’, correlated with the evolutionary stage of star forming regions, as shown by Emprechtinger et al. (2009) in the low-mass regime. Fontani et al. (2011) showed that the N₂D⁺∕N₂H⁺ isotopic ratio decreases from ≈0.2–0.4 to less than 0.1 when the temperature rises above 20 K in a sample of high-mass star forming regions in different evolutionary stages. Feng et al. (2019, 2020) reported similar findings using N₂H⁺ and HCO⁺ isotopologues. Despite the fact that N₂D⁺ is an abundant ion in cold and dense gas, it may not be the ideal tracer of HMPCs. N-bearing species are less affected by freeze-out onto the dust grain surfaces than molecules containing carbon or oxygen, but at high densities (n ≳ 10⁶ cm⁻³) N₂H⁺ and N₂D⁺ show signs of depletion as well (Pagani et al. 2007; Redaelli et al. 2019).

Transitions of lighter, deuterated molecules with high critical densities then represent the only reliable tracers of pre-stellar gas, since they can persist in the gas phase at high densities. H₂D⁺ is an ideal candidate, as shown by the complete depletion model of Walmsley et al. (2004). The (1_1,0 −1_1,1) transition of its ortho form (hereafter o-H₂D⁺) has a critical density of n_cr ≈ 10⁵ cm⁻³ (Hugo et al. 2009). Furthermore, H₂D⁺ is very sensitive to the temperature, due to the chemical reactions with its main destroyer, carbon monoxide. When (e.g. as a result of proto-stellar activity) the gas temperature rises beyond 20 K, the CO molecules frozen onto the dust grains evaporate back into the gas phase, thus lowering the H₂D⁺ abundance. Moreover, above 30 K, reactions (1) start to move backwards, thus reducing the deuteration level. Hence, detecting H₂D⁺ is an unambiguous sign of pre-stellar conditions in the sense that core structures that show significant emission in o- H₂D⁺ are undoubtedly in a pre-stellar phase.

Studies of H₂D⁺ in star forming regions are, however, rare and mainly performed with single-dish facilities, due to the weakness of its transitions. Caselli et al. (2003, 2008) produced a small sample of low-mass cores observed in o- H₂D⁺(1_1,0−1_1,1), and Miettinen (2020) targeted dense cores in Orion B9. In the high-mass regime, Pillai et al. (2012) spatially resolved the H₂D⁺ emission in the star forming region DR21 in Cygnus X using JCMT. Brünken et al. (2014) demonstrated the use of H₂D⁺ as a chemical clock for low-mass star forming regions. More recently, Giannetti et al. (2019) detected o- H₂D⁺(1_1,0−1_1,1) in three clumps belonging to the same filament with the Atacama Pathfinder EXperiment telescope (APEX; Güsten et al. 2006), showing that the abundance of this molecule anti-correlates with that of N₂D⁺. In a following study, Sabatini et al. (2020) published a detailed census of this tracer in 16 high-mass clumps with the APEX telescope. To our knowledge, only one interferometric study of this molecule has been performed: Friesen et al. (2014) used ALMA Band 7 observations to target the low-mass star forming region Ophiuchus. No high-mass dedicated study of H₂D⁺ with interferometers is present in literature.

In this work, we report the first interferometric observation of o- H₂D⁺(1_1,0−1_1,1) in two high-mass clumps performed with ALMA at ≈1″ resolution. The line was successfully detected in both sources, and its emission appears extended. Using a core-identifying algorithm, we find a total of 16 cores in H₂D⁺. We studied their properties using a spectral fitting analysis. Interestingly, a significant fraction of the gas emitting in o- H₂D⁺(1_1,0−1_1,1) presents low, subsonic velocity dispersion, often lower than the thermal broadening of the line at 10 K. This indicates that some of the clumps are still in a very cold, pre-stellar stage. We studied the correlation of the molecular emission and the dust thermal emission, also detected with ALMA at 0.8 mm. The line-integrated intensity usually does not correlate with the continuum data, and we speculate that this is due to different evolutionary stages of the identified cores. Among them, two present both bright o- H₂D⁺ emission and continuum emission. They represent ideal candidates to be pre-stellar cores embedded in high-mass clumps.

The observations and the resulting data are described in Sects. 2 and 3. We present the data analysis in Sect. 4, where we identify cores in position-position-velocity space in the o- H₂D⁺ data cubes, and we describe the line spectral analysis we performed. We discuss the results in Sect. 5: we first focus on the properties of the gas traced by the o-H₂D⁺ line (Sect. 5.1), and then on the comparison between the molecular and the continuum data (Sect. 5.2), using the dust thermal emission to estimate the gas column density and the sources’ total masses. In Sect. 5.3, we speculate on the possible evolutionary sequence of the identified cores. We discuss our findings in the context of star formation theories in Sect. 5.4. Section 6 contains the summary and concluding remarks of this work.

Fig. 1 Maps of the continuum emission at 870 μm seen with APEX by the ATLASGAL survey in AG351 (left) and AG354 (right). The FoV of the ALMA observations is indicated with a red circle. The beam size (18.2″) is shown in the bottom left corner, while the scale bar is in the bottom right.

2 Source selection and observations

2.1 The targeted sources

The sources selected for this work belong to the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL), which comprises a large number (~10 000) of massive clumps in different evolutionary stages (Schuller et al. 2009). In particular, the two chosen clumps also belong to the ATLASGAL TOP100 sub-sample (Giannetti et al. 2014; König et al. 2017), and their IDs are AGAL351.571+00.762 (hereafter AG351) and AGAL354.944-00.537 (hereafter AG354). Figure 1 shows the dust thermal emission at 870 μm as seen with APEX at a resolution of ≈18″. The clumps have similar masses and distances: M_clump = 170 M_⊙ and D = 1.3 kpc for AG351, M_clump = 150 M_⊙ and D = 1.9 kpc for AG354 (see König et al. 2017 and references therein). It is worth noting that these clumps have masses below the median value of the distribution (~ 500 M_⊙) obtained forthe quiescent stage classified in the ATLASGAL sample (see Urquhart et al. 2018). König et al. (2017) investigated the dust emission properties in the TOP100 sample using continuum data in the wavelength range 3–870 μm from severalinfrared surveys. From the fit of the spectral energy distribution, they estimated dust temperatures of 17–19 K, which must be considered average values over the source sizes (80–100″).

Concerning their evolutionary stage, both AG351 and AG354 are classified by the ATLASGAL catalogue as 24 μm-dark, since they are not associated with detected point-like sources at 24 and 70 μm. Furthermore, we inspected the available archive data at mid-infrared wavelengths: the Herschel-PACS maps at 70 μm and the Spitzer-IRAC maps at 3.6, 4.5, 5.8, and 8.0 μm (published in the GLIMPSE survey; Benjamin et al. 2003; Churchwell et al. 2009). No point source is detected within a radius of ≈ 30″. In addition, Kuhn et al. (2020) produced a Spitzer-IRAC catalogue of proto-stellar candidates in the Galactic midplane, and they did not find young stellar objects associated with the positions of the investigated clumps. This supports the scenario according to which AG351 and AG354 are in early pre-stellar phases, even though we highlight that several works showed that 70 μm-dark clumps werelater discovered containing protostars with intereferometric observations (Pillai et al. 2019; Li et al. 2019, 2020).

Both clumps are included in the sample investigated by Sabatini et al. (2020), who detected the o- H₂D⁺(1_1,0−1_1,1) transition with APEX in 16 ATLASGAL clumps. Their data have a beam size of 16.8″ and spectral resolution of 0.55 km s⁻¹. In AG351 and AG354, the line peak temperature is T_peak ≈ 200 mK and the linewidth is ≈ 1.0–1.2 km s⁻¹. Using the H₂ column densities from literature data (König et al. 2017; Urquhart et al. 2018), Sabatini et al. (2020) estimated an abundance with respect to molecular hydrogen of X_mol(o-H₂D⁺) ≈ 10⁻¹⁰ in both clumps.

2.2 Observations and data reduction

AG351 and AG354 were observed during Cycle 6 as part of the ALMA project #2018.1.00331.S (PI: Bovino) in two runs, from November 2018 to April 2019. The observations made use of both the Main Array (12 m-array 45 antennas) and the Atacama Compact Array (ACA, 12 antennas), with baselines ranging from 9 to 313 m. They were acquired as single pointings, centred at the sources’ coordinates: RA = 17^h20^m51.0^s, Dec = −35°35′23.29″ for AG351 and RA = 17^h35^m12.0^s, Dec = −33°30′28.97″ for AG354 (J2000). The quasars J1700-2610, J1717-3342, J1924-2914, and J1517-2422 were used as calibrators.

The spectral setup consists of four spectral windows (SPWs). One is centred on the o- H₂D⁺ (1_1,0 − 1_1,1) line at the frequency ν_rest = 372421.3558 MHz (Jusko et al. 2017). This SPW has a resolution of 244 kHz (corresponding to 0.20 km s⁻¹ at 372 GHz) and a total bandwidth of 500 MHz. A second spectral window is dedicated to the continuum, with a total bandwidth of 2.0 GHz, centred at the frequency of 371 GHz. The remaining two SPWs were centred on SO₂ and methanol lines, with a total bandwidth of 938 MHz and spectral resolution of 564 kHz (≈ 0.47 km s⁻¹). These lines, however, were not detected. At these frequencies, and with the configuration used, the maximum recoverable scale is θ_MRS ≈ 20″, the primary beam size is θ_FoV ≈ 26″², and the resolution is ≈1.0″ (corresponding to ≈1600 AU at the sources average distance of 1.6 kpc). The total observing times were 180 min (ACA) and 28 min (12 m-array) for each source. During the observations, the precipitable water vapour was typically 0.3 mm < PWV < 0.7 mm, and in general it was lower than 1 mm. The average system temperature values are found in the 300–400 K range for the SPWcontaining the o-H₂D⁺(1_1,0−1_1,1) line.

The data were calibrated by the standard pipeline (CASA, version 5.4; McMullin et al. 2007). The quality assessment outputs were checked, and the automatic calibration results are satisfactory. From a first inspection of the dirty maps, the emission both in continuum and in line appear very extended in the whole field of view (FoV). We therefore applied a modified weight of 2.4 to the ACA observations, chosen to maximise the recovery of large-scale flux without downgrading the resolution. The ACA and 12m-array observations were then concatenated to proceed with the imaging.

Both continuum and line data were imaged with the CASA tclean task interactively (CASA version 5.6), using a natural weighting and the multiscale deconvolver algorithm (Cornwell 2008). In order to maximise the bandwidth used to produce the continuum maps, we combined the continuum-dedicated SPW with the line-free channels in the other windows, obtaining a total bandwidth of ≈ 4.2 GHz. In order to obtain the data cubes of the o-H₂D⁺ line, we focused on the central 25 MHz (≈ 20 km s⁻¹) around the source centroid velocity. The data cubes hence present 100 channels in the frequency axis, and the final spectral resolution is 0.22 km s⁻¹. Finally, in order to avoid oversampling, both the continuum and the line images (already primary-beam-corrected) were re-gridded in order to ensure three pixels per beam minor axis, in agreement with the Nyquist theorem. The molecular line data were converted into the brightness temperature T_b scale, computing the corresponding gain, G = 11 mK mJy⁻¹ beam⁻¹ at a beam size of 1″ × 0.8″. In Appendix A, we discuss the possible missing flux problem affecting our data by comparing single-dish and interferometric observations.

3 Results

Figure 2 shows the continuum maps obtained as described in Sect. 2.2 for AG351 and AG354, respectively. The contours show the integrated-intensity maps of the o- H₂D⁺(1_1,0−1_1,1) transition, obtained integrating the cubes over the range [−4.3;−2.0] km s⁻¹ (AG351) and [−7.5;−4.0] km s⁻¹ (AG354). The continuum maps present a sensitivity of 0.8 mJy beam⁻¹ for both sources, while the median rms of the line data, computed over emission-free channels, is 300 mK per 0.22 km s⁻¹. These values refer to already primary-beam-corrected data. The observation resolutions and sensitivities for each source are summarised in Table 1.

4 Analysis

Our goal is to identify substructures seen in the o-H₂D⁺ emission in each clump and to analyse their physical and chemical properties. To reach our aim, we made use of an automated algorithm that makes it possible to identify such structures, fit the molecular line spectra, and reconstruct properties such as their column density, velocity dispersion, and centroid velocity. In the following subsections, we describe the steps of our analysis.

4.1 Core identification with SCIMES

In order to identify substructures in the molecular emission, we used the Spectral Clustering for Interstellar Molecular Emission Segmentation package (SCIMES, Colombo et al. 2015). SCIMES is based on dendrogram algorithms (Rosolowsky et al. 2008), and it was developed specifically to analyse molecular emission data in the form of position-position-velocity cubes. It identifies the substructures at the basis of the dendrogram (leaves), and it then finds the parental structures to which they belong (branches and trunks). To our scopes, we are interested in the smallest, significant substructures visible in o- H₂D⁺ emission. We therefore focussed on the leaves, which in the context of this paper correspond to pre-stellar cores.

In order to first build the dendrogram, a number of parameters are needed. We performed multiple tests, followed by visual inspection of the results, to find the best choice for our o- H₂D⁺ data. First of all, following Rosolowsky & Leroy (2006), we find the regions in the emission characterised by a signal-to-noise ratio (S∕N) higher than a given threshold (), which contain emission peaks that are brighter than a second higher threshold (). This implementation maximises the recovered information in case of moderate-to-low S/N data (S∕N = 3–10), such as our o- H₂D⁺ ALMA observations. In fact, the lines are locally bright (T_b > 0.5 K), but the peak temperature rapidly drops. We therefore set (S/N_)peak = 2, and (S/N)_lim = 1.5 for AG351, and (S/N)_peak = 2, and (S/N)_lim = 1.3 for AG354. By doing so, we set the minimum value (min_val) of the dendrogram algorithm to zero.

Another key parameter to build the dendrogram is the minimum height (in flux or brightness) that a structure must have to be catalogued as an independent leaf (Δ_min). In our case, we set Δ_min = 2.0 × rms. SCIMES performs better with data cubes with constant noise, since Δ_min has to be set as a single value. We therefore applied the algorithm to the cubes before applying the primary-beam correction. Their average noise is hence rms = 150 mK. We highlight that this correction applies a corrective factor dependent on the position within the primary beam, but uniform in the frequency axis. Therefore, it alters the absolute value of the rms, which becomes larger at the edges of the primary beam, but it does not affect the S/N map. We set the minimum number of channels that a leaf must span to , due to the low line width of the o-H₂D⁺ line. Finally, we excluded leaves smaller than two times the beam area. All the identified cores are found within regions where the integrated intensity is detected above a 3σ level ((S/N)_II > 3.0).

At the end of the process, we identified seven cores in AG351 and nine cores in AG354, which we labelled in order of decreasing peak value of o-H₂D⁺-integrated intensity within each clump. Figure 3 shows the identified structures on top of the continuum maps. It is worth noting that there are rare cases of core overlaps (e.g. cores 3 and 5 in AG354). This is due to the fact that SCIMES works in position-position-velocity space, and it is therefore able to disentangle structures that overlap spatially but present different centroid velocities (i.e. distinct velocity components along the line of sight). The first four columns of Table 2 report the cores’ identification labels, centre positions, and effective radii (R_eff, defined as the radius of a circle with the same area as the core). Several cores appear irregular in shape. This is due to i) the limited spatial resolution of our observations, combined with a Nyquist sampling of the beam size, which hence tends to make the borders of the identified structures more irregular: ii) the limited S/N of the o- H₂D⁺ data, which translates as the fact that the identified cores, despite being significant at a 3σ level in integrated intensity, are often at the limit of detection in peak brightness temperature; and iii) the lack of total-power observations, which should recover the most extended emission surrounding the cores. Moreover, we highlight that even insimulations, cores do not always appear regular in shape (see e.g. Körtgen et al. 2018; Smullen et al. 2020).

Figure 2 shows that the o- H₂D⁺-integrated-intensity distribution does not always follow the structure seen in continuum emission. The identified cores confirm this scenario: in some cases, dust thermal emission peaks are found within molecular-identified cores (e.g. core 3 in AG351, or core 2 in AG354). In other cases, bright continuum spots do not overlap with H₂D⁺ emission. We further investigate this point in Sect. 5.2 and Appendix B.

Fig. 2 Map of the continuum emission as seen by ALMA in Band 7 in AG351 (left panel) and AG354 (right panel). The contours show the o-H2D+-integrated intensity emission from the original not re-gridded and not primary-beam-corrected data, at levels = [5, 7, 9, 11]σ, where 1σ = 60 mK km s−1 (AG351) and 70 mK km s−1 (AG354). For the continuum data, we show the primary-beam-corrected maps. The beam size and the scale bar are shown in the bottom left and bottom right corners, respectively.

Fig. 3 Contours show the leaves identified by SCIMES, overlaid on the primary-beam-corrected continuum map in AG351 (left panel) and AG354 (right panel).

4.2 Spectral fitting with MCWEEDS

In order to derive the physical parameters of each core from the o- H₂D⁺ data, we performed a pixel-per-pixel spectral fitting of the lines using MCWEEDS (Giannetti et al. 2017). Briefly, the code combines the WEEDS package from GILDAS (Maret et al. 2011), which is optimised to produce synthetic spectra assuming local thermodynamic equilibrium conditions (LTE), with Bayesian statistical models implemented using PYMC (Patil et al. 2010) to perform the actual fit. We used a Markov chain Monte Carlo (MCMC) algorithm to sample the parameter space, using uninformative flat priors over the models’ free parameters. We recently upgraded MCWEEDS, parallelising the code using the methods of the python package MULTIPROCESSING³. This improves the code speed when dealing with large datasets, allowing us to fit spectra from multiple positions simultaneously. For each spectrum, the code performs 100 000 iterations and excludes the first 1000 steps (burn-in), which represents a reliable choice, as previously investigated (see Appendix B in Sabatini et al. 2020). The stability and convergence of the chains were checked in any case.

Since our ALMA data comprise only one o-H₂D⁺ transition, they do not allow us to constrain the molecular column density (N_col) and the excitation temperature(T_ex) independently. We therefore assumed T_ex = 10 K, a choice consistent with previous works; Caselli et al. (2008) found T_ex= 7–14 K in a sample of low-mass pre-stellar and proto-stellar cores, and Friesen et al. (2014) adopted T_ex = 12 K when studying pre-stellar cores in Ophiuchus A. The remaining free parameters in the fit are hence the molecular column density N_col(o-H₂D⁺), the line local-standard-of-rest velocity V_lsr, and the line full width at half-maximum (FWHM). For each core, individual initial guesses for the free parameters were chosen to improve the code convergence. In order to estimate the parameter distributions correctly, the noise level of the spectral line data has to be indicated. We adopted rms = 180 mK (AG351) and rms = 170 mK (AG354)⁴.

The parallelised version of MCWEEDS returns the maps of the best fit parameters, together with those of the lower and upper 95% high probability density (HPD) intervals. The best-fit models have been visually inspected, in order to assess their quality. In both clumps double-peak features, either due to multiple components on the line of sight, self-absorption, or large-scale emission filtering are present in the observed o- H₂D⁺ spectra. Even though they represent a small minority of the observations, these spectra cannot be successfully fitted with a single-velocity-component LTE model. Since these positions are characterised by an overestimation of the line width, we find that a successful masking strategy is to discard fits with a FWHM higher than a given threshold. This threshold, which was chosen for each core individually, is in the 1–1.2 km s⁻¹ range. The fraction of rejected pixels is ≲3% for both clumps, and we checked that no single-component pixels were masked. In Figs. 4 and 5, we show the resulting best-fit parameter maps in AG351 and AG354, respectively. Individual maps for each core are presented in Appendix C. As mentioned, MCWEEDS use the line FWHM as free parameter. However, we converted the obtained values of this quantity into velocity dispersion σ_V, in order to allow a more straightforward comparison with, for instance, the isothermal sound speed.

The median values for the relative errors derived from the 95% HPD intervals obtained with MCWEEDS are as follows: 57% in σ_V, 43% in N_col, and 3.5% in V_lsr for AG351; 38% in σ_V, 33% in N_col, and 1.9% in V_lsr for AG354. Uncertainties for AG351 are on average larger, since the spectra in this source are less bright. The centroid velocity is usually well constrained, while the velocity dispersion is characterised by the highest uncertainties. These are particularly large at low σ_V values, due to the limited spectral resolution of our observations, which is discussed further later on.

The V_lsr maps in Figs. 4 and 5 show that the cores are coherent structures in velocity, with dispersions around mean values of usually ≈ 0.2 km s⁻¹. At the clump scale, stronger velocity gradients of ≈ 1 km s⁻¹ and ≈ 2 km s⁻¹ are seen in AG351 and AG354, respectively. We avoid detailed discussion on the centroid velocity gradients seen in the clumps, since the lack of zero-spacing observations prevents us from inferring the large-scale gas kinematics (see Appendix Afor more details).

Typicalvalues for the molecular column density are in the N_col = (0.6–4) × 10¹³ cm⁻² range. These are higher than those found by Sabatini et al. (2020) using APEX observations, most likely due the dilution of the o- H₂D⁺ emission in the large single-dish beam (16.8″). The observed values are also, on average, higher than the maximum column density reported by Friesen et al. (2014) with ALMA at a higher spatial resolution (N_col = 1.2 × 10¹³ cm⁻³ at ~150 AU). However, those authors targeted a low-mass star forming region, focusing on a core significantly smaller and less massive than ours. From the column density maps, we estimated the optical depth τ maps using the spectral parameters listed in the Cologne Database for Molecular Spectroscopy (CDMS⁵) through the following equation: (2)

where A_ul = 1.20 × 10⁻⁸ s⁻¹ is the Einstein coefficient for spontaneous emission, g_u = 9 is the statistical weight, E_u = 17.9 K is the upper level energy, Q(T_ex) is the partition function (which at 10 K is Q = 10.48), k_B is the Boltzmann constant, and c is the light speed. We report the average value of τ and its dispersion around the mean value in each core in Table 2. Most of the cores present τ ≲ 1.0, indicating that the o-H₂D⁺(1_1,0−1_1,1) line is optically thin. Less than 7% of pixels in both clumps show optical depths higher than unity. This suggests that self-absorption does not strongly affect the line shapes, even though we cannot completely exclude this hypothesis given the low spectral resolution of ALMA observations.

The velocity dispersion maps of the clumps show that ALMA detects very narrow o- H₂D⁺ lines. In several cores, especially in AG351, large portions of the gas traced by H₂D⁺ present total line widths narrower than the isothermal sonic speed, which at T_gas = 10 K is⁶ : (3)

where m_H is the mass of the hydrogen atom, and μ is the mean molecular weight of the gas (μ = 2.33 for a gas composition of H₂ and 10% of helium). In particular, in AG351 36% of positions show subsonic line widths⁷, while the percentage for AG354 is 23%. The fraction of subsonic gas could be even higher, if we consider that the limited spectral resolution of our observations may cause an overestimation of the line widths. This has profound implication on the physical conditions traced by o- H₂D⁺(1_1,0−1_1,1), as discussedin detail in Sect. 5.1. Locally, the line velocity dispersion increases up to 0.4−0.5 km s⁻¹, without, however, showing a clear correlation with other quantities, such as the emission in continuum.

We compared the distribution of the velocity dispersion and column density in the two clumps in Fig. 6. These plots show that AG351 presents generally narrower lines and lower column density values than AG354, which could hint at differences in their evolutionary stage and/or physical conditions.

Fig. 4 Maps of o-H2D+ column density (top left panel), σV (top right), and centroid velocity (bottom) obtained from the MCWEEDS analysis in the individual cores in AG351. The beam size and scale bar are shown in the bottom left and bottom right corners, respectively. The black contours show the continuum emission in Band 7 at levels [1, 5, 10, 15, 20] mJy beam−1.

Fig. 5 Same as Fig. 4, but for AG354.

Fig. 6 Normalised kernel density distribution of σV and Ncol in AG351 (blue) and AG354 (red) obtained from the MCWEEDS pixel-per-pixel fit. The contour levels are [0.1, 0.3, 0.5, 0.7, 0.9]. The horizontal solid line shows the sound speed at 10 K, while the dot-dashed line represents the thermal broadening of H2D+ at the same temperature.

4.3 Average cores’ parameters

Important dynamical quantities of cores, such as the virial mass (see Sect. 5.1 for further discussion), depend on their average properties. In order to derive these parameters, we computed the average o- H₂D⁺(1_1,0−1_1,1) spectrum in each core. As highlighted in the previous subsection, the clumps present velocity gradients on the order of 1–2 km s⁻¹, while the individual cores show smaller variations of the centroid velocity. Nevertheless, given the narrow line widths of the o- H₂D⁺ spectra, even these small gradients could affect the line shapes of the averaged spectra if not taken into account. Before computing the mean spectrum from a core, therefore, we first aligned the spectra pixel-per-pixel to the mean centroid velocity obtained from the spectral fitting in Sect. 4.2.

The averaged spectra were then re-analysed using MCWEEDS, as previously described, setting T_ex= 10 K and fitting the line width, centroid velocity, andcolumn density of o-H₂D⁺. We computed the rms values of the average spectra considering line-free channels. Table 2 summarises the best-fit values, together with their HPD intervals at the 95% level, obtained with the method just described. We show the average spectra, overlaid with the best-fit model obtained with MCWEEDS, in Figs. 7 and 8. It is worth noting that all the spectra present Gaussian line shapes. This confirms first of all that secondary velocity components are, in general, negligible. Furthermore, it also suggests that opacity effects (such as self-absorption) are likely negligible as confirmed by the optical depth values computed in Sect. 4.2 (see Table 2), even though we cannot exclude them completely. In fact, due to the limited spectral resolution of the ALMA observations, the observed velocity dispersion may be locally overestimated, which in turn leads to the underestimation of the optical depth.

Dynamical parameters estimated from the observed velocity dispersions can provide useful information on the level of turbulence and its impact on the core dynamics. In particular, it is interesting to investigate whether the gas motions are subsonic or supersonic. The total velocity dispersion values derived from the average spectra (Table 2, seventh column) show that about half of the cores present subsonic or trans-sonic gas motions, within uncertainties. We can further disentangle the thermal andnon-thermal components of the lines. In fact, the velocity dispersion of any line is due to a combination of its thermal and non-thermal broadening, the latter comprising all terms which are not due to the gas temperature, such as turbulence, multiple components on the line of sight if spectrally unresolved, and instrumental broadening. Under the assumption that the two components are independent, they sum in quadrature to compose the total observed velocity dispersion: (Myers 1983). The thermal component of the o-H₂D⁺ line at the gas temperature T_gas is (4)

where is the H₂D⁺ molecular mass in g. We computed the ratio of non-thermal velocity dispersion and sound speed for each core from the observed velocity dispersion, as obtained with the spectral fit, assuming T_gas = 10 K, at which σ_V,th = 0.14 km s⁻¹. The one-dimension turbulent Mach number (σ_V,NT∕c_s) values are summarised in Table 2. In AG351, all cores present σ_V,NT∕c_s values lower than unity or consistent with unity within 95% HPD. Three cores in AG354 instead present slightly supersonic turbulent motions within the confidence intervals. In a similar way, Sabatini et al. (2020) found that at the clump scales the one-dimension turbulent Mach number in AG354 is higher than in AG351, even though according to their results both sources are mildly supersonic (σ_V,NT∕c_s = 1.4 and 1.8 in AG351 and AG354, respectively).

Fig. 7 Mean signal obtained averaging in each core the observed spectra is shown in black. The red curve shows the best-fit model obtained with MCWEEDS. In the top-right corner of each panel we report the core id number and the Vlsr value to which the spectra have been aligned, before averaging them.

5 Discussion

5.1 Gas traced by o-H₂D₊

As explained in Sect. 1, the o-H₂D⁺(1_1,0−1_1,1) transition probes a cold and dense component of the interstellar medium. However, having access to only one molecular transition prevents us from determining the gas temperature and density independently at the same time. In order to perform the spectral analysis, we assumed an excitation temperature of T_ex = 10 K, which is consistent with previous works (Harju et al. 2006; Pillai et al. 2012; Friesen et al. 2014; Sabatini et al. 2020) and with the fact that H₂D⁺ is abundant when the temperature is ≲ 20 K. However, dense pre-stellar gas can reach temperatures below 10 K, as shown both by theoretical works (Zucconi et al. 2001; Evans et al. 2001; Ivlev et al. 2019) and by observational evidence (Crapsi et al. 2007; Pagani et al. 2007).

The line width of the o-H₂D⁺(1_1,0−1_1,1) transition can give us hints concerning the gas temperature. As mentioned in the previous section, the velocity dispersion of any line can be separated in thermal and non-thermal broadening, which are summed in quadrature. It is therefore straightforward that the total line width of a line cannot be smaller than its thermal component. The thermal broadening of H₂D⁺ at 10 K is σ_V,th = 0.14 km s⁻¹. From the maps of this quantity obtained in Sect. 4.2, it results that 17% of pixels in AG351 and 7% in AG354 show line widths smaller than this level (see Fig. 6). This is a clear indication that at least some parts of the gas traced by o- H₂D⁺ have a temperature lower than 10 K. We can reverse this argument and derive the gas temperature at which the observed velocity dispersion values are purely due to thermal broadening: (5)

We hence derived maps of pixel-per-pixel in each source. This quantity represents an upper limit on the real gas temperature; from which point, when any contribution from the non-thermal component is present it holds that σ_V > σ_V,th. This approach has already been used, for instance, by Harju et al. (2008), who estimated a gas temperature of from APEX observations of o- H₂D⁺(1_1,0−1_1,1) in the starless clump Ophiuchus D. In our clumps, we constrain to be in the 5–20 K range. The lower limit is due to the fact that i) temperatures lower than 5 K are not predicted, even at very high volume densities (see, for instance, Hocuk et al. 2017), and ii) this is the temperature at which the thermal line width (in FWHM) equals the spectral resolution of our observations. The upper limit is instead due to the fact that the H₂D⁺ abundance is expected to drop at temperatures above 20 K (which corresponds to σ_V,th = 0.21 km s⁻¹). Line widths broader than this limit most likely involve non-thermal contributions, such as turbulence. A significant fraction of positions present (50% in AG351 and 31% in AG354).

The quantities derived from the MCWEEDS fit of the average spectrum in each core allow us to derive another important property, the virial mass (M_vir), which can be used to assess the dynamical state of a source. This is the mass at which the kinetic energy content of a system equals its gravitational potential energy, and it can be expressed in observational units as (Bertoldi & McKee 1992): (6)

where G is the gravitational constant. Equation (6) holds under the assumption of a spherical core of radius R_core and uniform density, and σ_dyn represents the total line width of the gas, and it can be derived from the observed σ_V of o- H₂D⁺ assuming that the non-thermal component of the latter is the same for all the gas components: (7)

In a first step, we estimated M_vir from the o-H₂D⁺ line widths and the effective radii listed in Table 2 (R_core = R_eff), using a constant gas temperature of 10 K. For each core, we computed the distribution of the virial mass from the corresponding distribution of σ_V, and we then determined the median and the 95% HPD intervals. The resulting virial masses are summarised in Table 3.

As previously discussed, we have reason to believe that at least part of the gas traced by o- H₂D⁺ is colder than 10 K. As a test, we therefore computed the virial masses using the minimum value in each core as derived from Eq. (5). Using algebra, we can show that σ_dyn (and hence M_vir) is positively correlated with temperature. Since T_H2D⁺ represents an upper limit for the gas temperature, the virial masses derived in this way also represent upper limits. In Table 3, we also report the M_vir values obtained with this approach, together with the corresponding values used for the computation. The virial mass does not strongly depend on the gas temperature, and the new M_vir values differ by only a few percent from the previous ones. It is interesting to notice that the estimated virial masses are in the 0.3–2.2 M_⊙ range, and hence these cores are essentially low mass, in the hypothesis that they are virialised by the kinetic energy. We further discuss this point in Sects. 5.2 and 5.4.

Fig. 8 Same as Fig. 7, but for AG354.

5.2 The correlation between molecular and continuum emission

Our ALMA data allow us to compare the gas properties as traced by the dust thermal emission and by the o- H₂D⁺ transition. As previously mentioned, the two sets of data do not present the same morphology. We investigate this point in detail in Appendix B, where we identify core-like structures in the continuum maps using a dendrogram analysis. Out of the 16 H₂D⁺-cores, five do not overlap at all with continuum-identified structures, and, conversely, three structures identified in continuum emission do not find correspondence in the SCIMES analysis. Furthermore, for two H₂D⁺-identified cores, the continuum flux peak is found outside the core boundaries. It is hence clear that, in general, the two datasets do not trace the same material in the clumps. In order to further investigate this point, we used thecontinuum maps to estimate the total gas mass and density. In high-density conditions (n ≳ 10⁴ – 10⁵ cm⁻³), dust and gas are well coupled, which means that they are thermalised at the same temperature (Goldsmith 2001). We can thus assume T_dust = T_gas = 10 K in order to compute the total gas column density N(H₂) and hence the mass of the identified cores. We computed pixel-per-pixel the quantity as follows: (8)

where f is the gas-to-dust ratio (assumed to be 100, Hildebrand 1983), D is the source’s distance; B_ν(T_dust) is the Planck function at the frequency ν and temperature T_dust, S_pix is the pixel flux (in Jy pix⁻¹) and A_pix is the pixel area, is the mean molecular weight per hydrogen molecule (Kauffmann et al. 2008), and κ_ν is the dust opacity at the frequency of the observations. For the latter, we used the following power-law expression: (9)

in which we employed β = 1.5 for the dust emissivity index (Mezger et al. 1990; Walker et al. 1990) and the κ₀ = 10 cm² g⁻¹ for the dust opacity at the reference frequency corresponding to the wavelength λ₀ = 250 μm (Hildebrand 1983; Beckwith et al. 1990). Errors on N(H₂) were computed with standard propagation from the uncertainties on continuum flux, and the average uncertainty is rms = 2.9 × 10²² cm⁻². We computed the average N(H₂) values for each core, which span the (1–4) × 10²³ cm⁻² range. From the gas column density, we derive the molecular abundance X_mol(o-H₂D⁺) = N_mol(o-H₂D⁺)∕N(H₂). The derived values are found in the (0.2–3) × 10⁻¹⁰ range for AG351 and (0.1–10) × 10⁻¹⁰ for AG354. These values are consistent with previous measurements. For instance, in low-mass pre-stellar and proto-stellar cores, Caselli et al. (2008) found X_mol(o-H₂D⁺) = (0.1–5) × 10⁻¹⁰.

The scatterplots of the o-H₂D⁺ line width and molecular abundance with respect to H₂ are shown in Figs. 9 and 10 (left panels), where we also show the average values for each core with star symbols. Concerning the relation between σ_V and the H₂ column density, no clear correlation is visible for AG351, while a tentative trend can be seen in AG354, where it seems that positions at high column densities (N(H₂) ≳ 5 × 10²³ cm⁻²) are characterised by narrow line widths (σ_V≲ 0.35 km s⁻¹). The molecular abundance seems anti-correlated with respect to N(H₂) for AG354, while the data do not show a clear trend in AG351. These differences in the scatterplots for the two clumps suggest that they are in two distinct evolutionary stages, as possibly also suggested by the fact that AG354 shows, on average, larger velocity dispersion and o-H₂D⁺ column density values. However, the long tail towards high gas densities in AG354 mainly arises from the data points of one core (core 2), and hence we cannot derive strong conclusions regarding the whole clump based on this.

The anti-correlation between X_mol and N(H₂) can be explained by the fact that at later evolutionary stages several factors can contribute to lowering the H₂D⁺ abundance: depletion of HD onto the dust grains (Sipilä et al. 2013), conversion of H₂D⁺ into its doubly and triply deuterated isotopologues, or a rise in temperature due to the proto-stellar feedback. Regarding the latter point, we lack the interferometric data to verify whether any proto-stellar source is embedded in the clumps. However, the fact that the observed velocity dispersion values in AG354 are higher than those of AG351 could be consistently explained by higher gas temperatures.

From the continuum emission, one can derive the total gas mass of the cores identified in H₂D⁺ using the following: (10)

where the notation is the same of Eq. (8), with the exception that S_tot is now the total flux (in Jy) integrated over the area of the core. The obtained values are summarised in Table 3. We computed M_core for cores where the peak flux is detected above the 5σ level. The ALMA absolute flux uncertainty in Band 7 is ≈10%. However, several parameters involved in the computation of M_core are not well constrained (i.e. the dust opacity). A variation of 30% of κ_ν translates to a variation of 33% in the mass. We therefore increased the error on M_core to this level, inorder to take into account the parameters’ uncertainties.

The obtained masses range between 0.4 and 6 M_⊙, with larger values observed in AG354. However, these values might not represent the entire mass budget of the H₂D⁺-identified cores, particularly for the ones where the continuum emission is weak and anti-correlated with the molecular one. Furthermore, the interferometer may filter out the larger scale envelope emission, hence underestimating the cores’ total masses.

Similarly to what was discussed for the gas temperature, the assumption of a constant dust temperature T_dust = 10 K is questionable. We do not have access to other data to better constrain this quantity. However, based on the assumption of dust-gas coupling, it holds that if the o- H₂D⁺ emission traces part of the gas colder than 10 K, the dust temperature can also be lower than this value locally. We therefore re-computed the N(H₂) (and hence the o- H₂D⁺ abundance) pixel-per-pixel using the gas temperature maps obtained from the o-H₂D⁺(1_1,0−1_1,1) line widths, as discussed in Sect. 5.1. The average uncertainties on the new N(H₂) values are rms = 2.1 × 10²² cm⁻² for AG351 and rms = 1.4 × 10²² cm⁻² for AG354, and they are computed as standard error propagation from the flux uncertainty. The core’s average column density values are N(H₂) = 0.3–4 × 10²³ cm⁻². The new scatterplots showing the correlation of σ_V and X_mol(H₂D⁺) with respect to N(H₂) are shown in the right panels of Figs. 9 and 10. With the updated values, a trend between the observed velocity dispersion values and the gas total column density is visible for both sources. This is expected, however, since the narrower the line, the lower the derived , and the higher N(H₂) (for a given flux, dust temperature and mass or density are anti-correlated). On the other hand, the anti-correlation between the molecular abundance and the H₂ column density is now clear in both clumps.

The anti-correlation between continuum and deuterated species has already been observed in the literature. For instance, Zhang et al. (2020) found a shift between NH₂D and dust thermal emission in several high-mass star forming regions. Concerning o- H₂D⁺, Friesen et al. (2014) similarly reported a shift between the molecular emission peak and the continuum peak. Sabatini et al. (2020) used the detections of o- H₂D⁺(1_1,0−1_1,1) in a sample of ATLASGAL sources to fitthe relation between X_mol(o-H₂D⁺) and N(H₂), which in the log–log scale reads: . We plot the predictions from this relation in Fig. 10. We fitted a linear relation in the log–log space between the molecular abundance and the H₂ column density for our data. We find slopes of m ≈−0.8, −0.9 with high correlation coefficient (), with the exception of the X_mol(o-H₂D⁺)-N(H₂) relation in AG351 with T_dust = 10 K, where the data do not present a significant trend. The relations for our datasets show similar slopes to that found by Sabatini et al. (2020), but they are shifted upwards by a factor of 5–10. Interestingly, the cores identified by SCIMES are on average 2.5″ wide in angular size. Hence, the smaller abundances of Sabatini et al. (2020) can be explained by beam dilution, if one takes into consideration the ratio of the source sizes with respect to the APEX beam, assuming ≈ 10 cores in each clump: .

It is important to notice that X_mol(o-H₂D⁺) and N(H₂) are not independent variables. This point has been investigated deeply by Sabatini et al. (2020), who performed a robust statistical analysis of the correlation between these two quantities. In this work, we are more interested in the general trend that X_mol(o-H₂D⁺) shows in relation to the gas total density and in comparing our results with the literature ones. However, to further test whether the correlations that we see in Fig. 10 are real, we checked the correlation between the molecular column density N_col(o-H₂D⁺) and N(H₂) (not shown here). These two quantities are in fact independent⁸. We do not find any significant correlation between these variables, as previously noted by Sabatini et al. (2020). The o- H₂D⁺ column density is flat with respect to N(H₂). This confirms that higher N(H₂) values present lower molecular abundances.

We computed the core masses using the of each core, as we did for the virial masses. We stress again that this temperature value is the minimum value found within each core, but it represents an upper limit, due both to possible non-thermal components to the line width, and to the instrumental broadening caused by the limited spectral resolution of the observations. Hence, in reality the cores’ masses could be higher than these new M_core values, which are listed in Table 3. It is interesting to notice that according to the dust masses computed at 10 K, the cores are essentially low mass. Only under the assumption of very low dust temperatures (T_dust ≈ 5 K) do the estimated masses become larger, but in any case M_core ≲ 13 M_⊙, with the only exception of core 2 in AG354. This holds also for the cores identified in the continuum emission maps in Appendix B. Out of the 15 continuum-identified structures, only two are more massive than 20 M_⊙ under the assumption of T_dust = 5 K, one of which does not present significant o-H₂D⁺ emission, hence preventing us from classifying it as truly pre-stellar. We further discuss these findings in Sect. 5.4. It is worth noting that the lack of total-power observations means the ALMA data might be filtering out the large-scale emission associated with the gas surrounding the cores identified here, hence leading to an underestimation of their masses based on the dust thermal emission (see Appendix A for more details).

From the dust masses and effective radii, we can estimate the average volume density of the cores, under the assumption of constant distribution and spherical geometry, using the following equation: (11)

We evaluated n(H₂) for both cases: assuming T_dust = 10 K and using derived from the o-H₂D⁺(1_1,0−1_1,1) line width. The resulting values are also summarised in Table 3. Given the uncertainties of this estimation (e.g. due to the assumption of uniform, spherical distribution), we expect these values to be accurate within a factor of 2. In all cases, even when the dust mass is computed at a high temperature of 20 K, the average density is higher than 10⁶ cm⁻³. This confirms that the ALMA data trace high-density material, and it justifies the hypothesis of dust-and-gas coupling. Furthermore, at such high densities the o-H₂D⁺(1_1,0−1_1,1) is expected to be thermalised, and hence LTE is a good approximation (Harju et al. 2008). The volume densities found correspond to surface densities in the Σ = (0.1−10.0) g cm⁻² range, depending on the core and on the assumed temperature.

The ratio between the virial mass and the total mass, known as virial parameter (α_vir = M_vir∕M_core), provides indication on the dynamical state. For sub-virial sources (α_vir < 1), the kinetic energy content alone is not enough to balance the gravitational pull, and if no other source of pressure is present, they will collapse. On the other hand, super-virial sources (α_vir > 1) are unlikely to undergo gravitational contraction. We computed the virial ratio in both cases, with a constant temperature of 10 K and with , and we report the two sets of values in Table 3. Most cores are sub-virial, or consistent with being sub-virial within the 95% HPD intervals in AG351, regardless of the assumed temperature. In AG354, the super-virial cases are represented by cores 5, 6, and 7 (i.e. 33% of the sample). This is consistent with what was found at the clump scales by Sabatini et al. (2020), who reported that AG351 is more sub-virial than AG354 (α_vir = 0.4 and 0.8, respectively). We cannot exclude, however, that the cores are virialised by sources of pressure other than the kinetic one, and in particular by the presence of magnetic fields. This point is further discussed in Sect. 5.4.

Fig. 9 Scatterplots of o-H2D+(11,0−11,1) velocity dispersion σV with respect to H2 column density N(H2) in AG351 (blue, top panels) and AG354 (red, bottom panels). In the left panels, N(H2) is computed assuming a constant temperature of 10 K. In the right panels, we use the temperature derived pixel-per-pixel from the o- H2D+(11,0−11,1) line width, in the hypothesis of only thermal contributions. In order to allow an easier comparison between the two clumps, the plot ranges are kept equal in the corresponding panels. For the sake of readability,we only show data-points for which the relative 95% HPD interval derived with MCWEEDS is < 50%. The black stars show the average values referring to the H2D+-identified cores, for which the continuum flux is significantly detected. The σV values are taken from Table 2, while the N(H2) values are computed as averages in each core.

Fig. 10 Scatterplots of o-H2D+ molecular abundance with respect to H2 column density in AG351 (blue, top panels) and AG354 (red, bottom panels). In the left panels, the quantities are computed assuming a constant temperature of 10 K. In the right panels, we use pixel-per-pixel the temperature derived from the o- H2D+(11,0−11,1) line width, in the hypothesis of only thermal contributions. Scatterplot ranges are fixed to allow an easier comparison.For the sake of readability, we only show data points for which the relative 95% confidence interval is < 50%. The dashed black curve shows the correlation found by Sabatini et al. (2020). The black stars show the average values referring to the H2D+-identified cores, for which significant continuum flux is detected. The Ncol values to compute Xmol(o-H2D+) are taken from Table 2, and the N(H2) values are computed as averages in each core.

5.3 A population of cores at different evolutionary stages

We can speculate that the difference in correlation between the dust continuum emission and the H₂D⁺ distribution in the cores is linked to their evolutionary stages. In fact, H₂D⁺ forms when the gas becomes dense and cold enough for reaction (1) to be efficient. As the gas becomes denser due to contraction motions, several factors could contribute to lowering the H₂D⁺ abundance. H₂D⁺ can be transformed into D₂H⁺ and or deplete as a consequence of HD depletion (Sipilä et al. 2013). Secondly, one has to take opacity effects into account: if the volume density grows much larger than the line critical density, the transition becomes optically thick (τ > 1), which means that it no longer traces the whole bulk of gas, only the layer at τ ≈ 1. In Sect. 4.2, we described the optical depths from the derived N_col maps, and despite their being a minority, there are positions where lines are moderately optically thick. This could contribute to explaining why the observed velocity dispersion values in AG354, which we believe to be more evolved, are broader than in AG351. In fact, optically-thick lines deviate from Gaussian shapes, which leads to the overestimation of their widths. When the pre-stellar core collapses and forms a proto-stellar object, the abundance of H₂D⁺ is finally reduced by the rising temperature.

We can therefore speculate that cores with strong H₂D⁺ emission, either corresponding to a continuum-identified structure or not, are early in their evolution and truly pre-stellar. The lack of detected continuum peaks and associated continuum-identified cores could be due to observational limits. In fact, if the dust temperature in the correspondence of the o- H₂D⁺ emission is lower than in the surroundings – as expected in the case of pre-stellar cores – we may not be able to detect a centrally concentrated structure at 0.8 mm continuum emission, as shown, for instance, by Di Francesco et al. (2004). Hence, cores that do not have a correspondence in continuum-identified structures could still be centrally peaked in reality.

As the cores evolve, the H₂D⁺ abundance at the dust continuum peak is lowered for the aforementioned reasons (e.g. the conversion to D₂H⁺ and ). H₂D⁺-identified cores that partially overlap with continuum structure but present a continuum flux peak outside their boundaries, may belong to this evolutionarystage. Cores 4 and 6 in AG354 are representative of this case (as shown in Appendix B). The fact that they do not present significant increases of the velocity dispersion in the proximity of the continuum peaks suggest that the decrease of X_mol(o-H₂D⁺) is most likely due either to HD depletion onto dust grains or to conversion into other isotopologues. Finally, cores that are only visible in continuum are the most evolved. The lack of detection of H₂D⁺ in their surroundings suggests that either their temperature is overall higher than 20 K, suggesting an already proto-stellarstage, or that H₂D⁺ is depleted onto dust grains and/or converted into its doubly and triply deuterated forms.

In order to test our hypothesis, we searched for correlations of the suggested evolutionary phases and the cores’ dynamical properties reported in Tables 2 and 3 such as line widths, virial masses, and average densities. We did not find any significant trend. Nevertheless, Fig. 10 shows a decrease of X_mol(o-H₂D⁺) with an increase of N(H₂), which could resemble a situation where the collapse is more advanced, at least in some of the cores. Further observations of tracers of proto-stellar activity (such as SiO jets or CO outflows) or kinematics tracers could confirm this hypothesis.

According to our analysis, cores 3 and 7 (AG351) and core 2 (AG354) represent good candidates to be truly pre-stellar cores in high-mass clumps, since they overlap significantly with continuum-identified structures with masses M ≥ 14 M_⊙. Figure 11 shows the distribution of the o-H₂D⁺ abundance in core 2 in AG354, the most massive structure we identify, with contours from the continuum (thus representative of the H₂ column density, with the assumption of constant dust temperature). It can be seen that the molecular distribution appears anti-correlated with the continuum one: at the dust peak, X_mol(o-H₂D⁺) present the lowest values, while the molecular abundance increases moving away from it, in a similar fashion to the general trend shown by the scatterplots in Fig. 10. This core shows particularly low line widths (⟨σ_V⟩ = 0.20 km s⁻¹), which become as low as 0.10 km s⁻¹ towards the core centre, where the temperature sinks to ≈ 5 K. If we assume this dust temperature, the core mass is M_core = 39 M_⊙, significantly more massive than all the other cores studied in this work, and the mass of the associated continuum-identified structure (c-9 in Appendix B) could be as high as 60 M_⊙. Core AG354-2 thus represents an ideal HMPC candidate.

The evolution discussed here can either be chemical or physical (i.e. a density evolution). The first case holds if the cores are virialised byother sources of pressure other than thermal and turbulent motions, and hence their density is kept approximately constant while the chemistry evolves. In the initial stages, due to the high volume densities and low temperatures, o- H₂D⁺ is efficientlyformed, and hence the cores show bright molecular emission. With time, the o- H₂D⁺ abundance decreases, because the molecule is depleted onto the dust grains and/or converted into other isotopologues until the o- H₂D⁺(1_1,0−1_1,1) at the continuum peaks is not detectable anymore.

If the cores are instead sub-virial, and hence they experience gravitational contraction, the density is evolving and increasing with time. The evolution of X_mol(o-H₂D⁺) is similar, because the high initial abundance is then lowered, for instance, when the increasing density causes more depletion onto dust grains. In this case, at a certain point the collapse will form a central protostar, whose feedback will reduce the o- H₂D⁺ abundance below the detection limit. Lastly, we cannot exclude a combination of the two cases, where cores experience periods of time in equilibrium, chemically evolving at a constant density, followed by contraction phases.

Fig. 11 Xmol(o-H2D+) distribution in core 2 in AG354, the most massive structure identified in this work. Xmol (o-H2D+) was computed adopting a constant Tdust = 10 K. The contours show the ALMA continuum flux in Band 7, at levels [5, 10, 15] mJy beam−1. Under the assumption of constant dust temperature, the contours also represent the distribution of N(H2).

5.4 Our results in the context of star formation theories

The analysis of our data shows that the H₂D⁺-identified cores have masses in the 0.1–13 M_⊙ range, depending on the assumed temperature. These values, which are probably accurate within a factor of ≈ 5 (see Appendix A), would suggest that several of the H₂D⁺-identified cores are essentially low mass. The only clear exception is core AG354-2 for which, under the assumption of a dust temperature value T_dust = 5 K, we estimate M_core = 39 ± 13 M_⊙. We highlight that in order to form a 8 M_⊙ star, a core of at least 30 M_⊙ is needed, assuming a star formation efficiency of 30%.

Furthermore, both at core and clump scales the analysed sources show sub-virial conditions, in contrast with the predictions of turbulent-core accretion models of HMPCs. In our sample of H₂D⁺-identified cores, 80% are consistent with α_vir ≲ 1, regardless of the assumed temperature. The virial parameter computed in the previous section, however, only takes into account the kinetic energy content (T) and the gravitational energy (U), whereas to be more comprehensive the magnetic energy term should be taken into account, which in the case of a homogeneous and spherical core threaded by a uniform magnetic field B reads: (12)

where the equation is expressed in cgs units. We can determine the strength of the magnetic field neededto halt the gravitational collapse by imposing the following virialisation: U +2T + Ω_B = 0. In our sample of cores, the estimated B values range from hundreds of μG to several mG. Direct magnetic field measurements via the Zeeman effect are rare, but we note that in Crutcher & Kemball (2019), at densities of n ≈10⁵−10⁶ cm⁻³ the authors report values of B ≈a few hundreds μG, but never higherthan 1 mG. However, analysing the polarised dust thermal emission, several authors found magnetic fields of the order of a milliGauss in massive star-forming regions (Girart et al. 2013; Zhang et al. 2014; Liu et al. 2020). In conclusion, it is plausible that magnetic fields provide the extra support needed to virialise the cores, even though for those with very low virial ratios (α_vir < 0.3), magnetic strengths of several mG would be needed. We highlight the need for further observations to investigate the level of magnetisation of the clumps and their large-scale kinematics.

We also estimated the Jeans masses for the H₂D⁺-identified cores. For volume densities in the n =10⁶–10⁷ cm⁻³ range and temperatures of 5–10 K, the Jeans mass is M_Jeans = 0.03–0.3 M_⊙. Hence, almost the totality of the H₂D⁺-identified cores contain several tens, up to hundreds, of Jeans masses. This will cause a large degree of further fragmentation in the more massive cores, unless magnetic pressure or a high-level of turbulence are assumed. The analyses we have performed on the line widths show that the motions of the gas traced by o- H₂D⁺(1_1,0−1_1,1) is at most slightly supersonic (σ_V,NT∕c_s < 2 in all cases), similarly to what is found at the clump level (see Sabatini et al. 2020), pointing towards the necessity of magnetic fields to prevent fragmentation.

If we exclude that the mass values derived from the ALMA observations are significantly underestimated, it looks reasonable to evaluate two alternative possibilities: i) these cores are going to form low-to-intermediate-mass stars if no further accretion is considered, with perhaps the exception of core AG354-2, or ii) these cores can be clump-fed and eventually form high-mass stars as in the competitive accretion model (Bonnell & Bate 2006), growing during the low-mass proto-stellar stage similarly to what was also proposed by Tigé et al. (2017) and Motte et al. (2018). We highlight how AG351 and AG354 sits in the low end of the clump masses in the ATLASGAL sample. Kauffmann & Pillai (2010) investigated the threshold for high-mass star formation in Galactic IRDCs empirically, and according to their Eq. 1 our clumps are above the mass limit, and in principle they can form high-mass stars. On the other hand, Sanhueza et al. (2017) computed the minimum clump mass needed to form at least one massive star (M > 8 M_⊙), assuming a Kroupa initial mass function (IMF) and a star formation efficiency of 30%, which is . This threshold is however strongly dependent on the assumed IMF shape and star-formation efficiency, and it is hence very uncertain. We therefore cannot exclude that more massive clumps, representative of the high-mass clump population (M ≈ 500 M_⊙), can host several high-mass pre-stellar cores. Further studies focused not only on the continuum emission but also on H₂D⁺ on more massive clumps are then needed to assess the existence of massive pre-stellar cores and disentangle different theories.

Up to now, in literature the best candidates of HMPCs have been the C1S core identified by Tan et al. (2013; M = 10–50 M_⊙ depending on the assumed temperature), the one by Cyganowski et al. (2014; 30 M_⊙), and the core W43-MM1#6 studied by Nony et al. (2018). In the first study, the authors reported high abundances of N₂D⁺, and in a follow-up paper a non-detection of o- H₂D⁺ (Kong et al. 2016). Based on the results from Giannetti et al. (2019) (see Sect. 1 for more details), the C1S core could hence already host an embedded protostar in its very initial stages, or it could be perturbed by the activity of the two nearby, proto-stellar cores seen by Tan et al. (2016). The source investigated by Cyganowski et al. (2014), on the other hand, has only been identified in continuum, which prevents a conclusive assessment of its evolutionary stage. Core 6 in W43-MM1, with M = 60 M_⊙ and no outflow detected, remain a good candidate to be an isolated and massive HMPCs. Future observations targeting molecular tracers of dense and cold gas are needed to validate this hypothesis. Recent studies (see e.g. Sanhueza et al. 2019; Contreras et al. 2018) pursued with high-resolution ALMA observations revealed a large population of low-mass cores in high-mass clumps, showing sub-virial dynamical states, in agreement with what we report in this study.

The main advantage of this work is the detection of o-H₂D⁺, which traces the very early stages, and we can then state that these are possibly the first pre-stellar cores identified unambiguously in high-mass star forming clumps. Further observations focusing on the magnetic properties of the observed cores will provide us with the detailed analysis of the magnetic fields needed to assess conclusively their dynamical states. In addition, numerical studies including chemistry are needed to shed light on the different proposed theories and disentangle the different physical processes that can affect the star formation process. Examples of these works, including detailed deuteration chemistry, are beginning to emerge (Goodson et al. 2016; Körtgen et al. 2018; Bovino et al. 2019; Hsu et al. 2021), and a proper comparison of these simulations with observations is a viable way to find an answer to this longstanding problem.

6 Summary and conclusions

In this work, we report ALMA observations in Band 7 at a resolution of ≈ 1″ in two quiescent, intermediate-to-high-mass clumps. For the first time, we report the detection of o- H₂D⁺(1_1,0−1_1,1) in this kind of sources with interferometric observations. Our molecular line data show that o- H₂D⁺ is very extended, and its distribution does not correlate with the one of the dust thermal emission in the same ALMA band (λ = 0.8 mm).

Using the algorithm SCIMES, we identified 16 cores in the o- H₂D⁺ data cubes. We fitted their spectra pixel-per-pixel using a Bayesian approach implemented in the code MCWEEDS in order to derive the line velocity dispersion, the molecular column density, and the centroid velocity map. The first important conclusion is the detection of narrow o- H₂D⁺(1_1,0−1_1,1) lines, with line widths lower than their thermal broadening at 10 K. This indicates that this species traces a very cold and quiescent region in the analysed sources.

We investigated the general lack of correlation between the dust continuum peaks and the molecular line peaks, which has profound implications. We suggest that this is due to a possible physical and/or chemical evolution. In the initial stages o- H₂D⁺ is quite abundant, and tracing the high-density gas. The presence of cores that are bright in the molecular emission, but lacking a continuum peak nearby, can be due to the fact that the ALMA Band 7 observations are not able to trace the peaked dust distribution due to temperature effects (see Sect. 5.3 for more details). Later on, as the gas becomes denser, or as the chemistry evolves, the X_mol(o-H₂D⁺) starts to decrease, but the molecule is still detectable. Eventually, the H₂D⁺ abundance drops, due either to depletion at very high densities, or to proto-stellar feedback effects.

Our results highlight how the continuum emission alone is generally not a good probe of pre-stellar gas. Bright cores in continuum flux that do not show significant emission in cold-gas tracers (such as H₂D⁺) are likely in a more evolved, possibly proto-stellar stage. This suggests that particular care must be taken when doing surveys of ‘pre-stellar’ cores seen only in continuum, especially when only one frequency band is available since these data do not allow us to determine the temperature, the evolutionary stage, and the kinematic properties of the sources. Complementary molecular data should be used to distinguish pre- and proto-stellar objects. In this context, we find that o- H₂D⁺ represents a good tracer of cold gas at densities of n ≈ 10⁶ cm⁻³. This gas is in a pre-stellar phase, in the sense that it has not been influenced by proto-stellar activity, and – being cold and dense – it has the potential to form new proto-stars; its evolution, however, is determined by its dynamic state.

At densities higher than 10⁶ cm⁻³, even o- H₂D⁺ is no longer ideal, as shown by its abundance drop as a function of N(H₂) due to opacity effects, depletion, or a combination of both. In those physical conditions, D₂H⁺ and represent probably the only good tracer available. While the latter is not observable, the former could represent, in combination with o- H₂D⁺, the best choice to investigate the more evolved and denser regions.

Most of the H₂D⁺-identified cores are less massive than 10 M_⊙, even under the assumption of low dust temperatures (5 K). It is important to highlight, however, that the ALMA observations could be filtering out the most extended core envelopes, and hence leading to the underestimation of their total masses. Furthermore, we cannot exclude that they are still in the process of accreting material from the parental clump. On the other hand, our data could support a scenario in which high-mass clumps fragment in a population of low-to-intermediate cores, which can continue the accretion to larger masses during the later, proto-stellarphase. Two of the cores contain peaks of both molecular and continuum emission, which allows us to reliably estimate their total masses. Core AG354-2 shows particularly narrow line widths, consistent with temperature values as low as 5 K. With this dust temperature, we estimate a core mass of 39 M_⊙, which could represent a lower limit for the aforementioned reasons. Furthermore, it is associated with a continuum-identified core, which contains several tens of solar masses. To our knowledge, it represents an ideal candidate of HMPCs, but further investigation is needed to better constrain its temperature and thus mass.

We investigated the dynamic state of the H₂D⁺-identified cores by means of the virial analysis. Most of the cores appear sub-virial if we take into account the gravitational energy and the kinetic energy only. We do not have information on the magnetic properties of the sources. We estimate that magnetic field strengths of the order of several hundreds of microGauss or a few milliGauss are needed to virialise the cores. Further observations aimed at recovering the magnetic field properties are needed to make a definite conclusion on the core dynamic states.

As a future perspective, we plan to recover the missing flux in the large-scale emission of o- H₂D⁺(1_1,0−1_1,1) (see Appendix A). This will allow us to use this transition to investigate the kinematics of the gas at the clump scales, which will in turn provide key information about the cores’ dynamics. Furthermore, ALMA observations at a similar spatial resolution of molecules that are good tracers of the proto-stellar activity (SiO, CO, etc.) will help us unambiguously disentangle pre-stellar cores from proto-stellar ones.

Acknowledgements

We thank the anonymous referee, for her/his suggestions to improve the manuscript. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.00331.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. In addition, publications from NA authors must include the standard NRAO acknowledgement: The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research made use of SCIMES, a Python package to find relevant structures into dendrograms of molecular gas emission using the spectral clustering approach. ER is thankful to Dominique M. Segura-Cox, for the help provided in the data reduction. SB acknowledges the BASAL Centro de Astrofisica y Tecnologias Afines (CATA) AFB-17002. The data analysis has been performed with resources provided by the KULTRUN Astronomy Hybrid Cluster at Universidad de Concepción. DC acknowledges support by the Deutsche Forschungsgemeinschaft, DFG project number SFB956A.

Appendix A ALMA and APEX comparison

In Fig. A.1, we show the comparison of the o- H₂D⁺(1_1,0−1_1,1) spectra obtained with APEX (from Sabatini et al. 2020) and with ALMA (this work) in AG351 and AG354. In order to compare the two datasets, we converted the temperature scale of the APEX observations in flux, using the gain G = 40 Jy K⁻¹ (from the APEX telescope efficiency webpage⁹). The ALMA data were integrated over an area equal to the APEX beam size, and then convolved to the same velocity resolution (≈ 0.55 km s⁻¹).

Fig. A.1 Comparison of the o-H2D+(11,0−11,1) spectra observed by APEX (red) and ALMA (black) in AG351 and AG354, from left to right. For a fair comparison, the APEX spectra from Sabatini et al. (2020) were converted in flux units using the APEX gain (40 Jy K−1). The ALMA spectra were computed integrating the signal over an area corresponding to the APEX beam size (16.8″), and they were smoothed to the APEX spectral resolution.

The ALMA spectra in both sources present a peak flux of ≈ 1∕3−1∕5 of the APEX peak. This is due to the so-called missing-flux problem, which is caused by the fact that the molecular emission is diffuse over scales larger than the maximum recoverable scale θ_MRS ≈ 20″. As a consequence, a significant amount of the emission is filtered out. However, it is important to notice that the ALMA integrated spectra do not present anomalous line shapes, and their centroid velocity is consistent with the one from the single-dish observations. Furthermore the core apparent sizes are smaller than the maximum recoverable scale θ_MRS ≈ 20″. All this considered, we do not expect that the missing-flux problem severely affects the line profiles in the cores, and hence we believe that the kinematics parameters (σ_V, V_lsr) obtained from the spectral fitting are reliable. However, we do highlight the fact that future observations recovering the zero-spacingflux are needed (see e.g. Henshaw et al. 2014; Sokolov et al. 2018) if we want to use the ALMA o- H₂D⁺(1_1,0−1_1,1) data to discuss the large-scale kinematics of the clumps.

Concerning the continuum data, we can estimate the fraction of flux filtered out by the interferometer comparing the ALMA and the APEX observationsat 870 μm of the clumps from the ATLASGAL project. Integrating over an area equal to the ALMA FoV, we obtain S_{870 μm} = 2.6 Jy (AG351) and S_{870 μm} = 1.8 Jy (AG354), while the total flux seen in the ALMA data is S_{0.8 mm} = 1.0 Jy (AG351) and S_{0.8 mm} = 0.93 Jy (AG354). This is just a rough comparison, since we are not taking into account the ≈ 10% difference in the observed wavelength, but it shows that the interferometer is filtering out ≈ 50−60% of the emission.

Despite most of the missing flux coming from the large clump scales, we cannot exclude that part of the cores’ envelopes are also filtered out, hence leading to mass underestimation, as observed in the nearby low-mass pre-stellar core L1544 by Caselli et al. (2019). It is, however, not straightforward to determine by which factor the masses are possibly underestimated, since the APEX single-dish observations themselves could underestimate the flux due to filtering of the large-scale emission. Therefore, the clump total masses could also be underestimated, and the total mass budget available for the cores could be higher than ≈ 150 M_⊙. In a conservative approach, we estimate that the mass values could be underestimated of a factor of 2− 5.

Appendix B Continuum-identified cores

In the main text, we focus on the properties of the H₂D⁺-identified cores. However, it is also worth looking to the continuum-identified structures, which is the purpose of this appendix. SCIMES was developed to work in position-position-velocity space, and hence it is suitable to analyse molecular line data only. We hence used the python package ASTRODENDRO¹⁰, on which SCIMES is also built, to identify structures in the 0.8 mm continuum maps. After a few tests followed by visual inspection of the results, we selected the dendrogram parameters as it follows: min_val = 3σ, Δ_min = 1σ, and min_area = 2 beam (i.e. we exclude structures smaller than two beams), a choice consistent with similar works (see e.g. Barnes et al. 2021). As previously stated, dendrogram analysis algorithms perform better on constant-noise maps, and hence we input 1σ = 0.5 mJy beam⁻¹ as the rms level of the continuum maps before the correction for the primary beam response.

We find 15 cores in total, which are shown in Fig. B.1, together with the H₂D⁺ cores identified in Sect. 4.1. As previously discussed, there is no one-to-one correspondence in the identified structures. Some of the brightest figures seen in dust thermal emission do not correspond to cores visible in o- H₂D⁺, and vice versa. Several structures overlap, but it is interesting to notice that for some of them (c-2 and c-4 in AG354) the continuum peak is found at the border of (or just outside) the H₂D⁺-identified cores. In Table B.1, we report the label of the corresponding H₂D⁺-identified core that overlaps the most with each of the continuum-identified ones, together with the fraction of overlap. One can see that several structures overlap only partially. For nine of the continuum-identified cores, the overlapping fraction is in the 0–70% range.

Similarly to what was done for the H₂D⁺-identified cores, we computed the sizes in term of effective radii and the masses of the cores identified in continuum emission. As extensively discussed in Sects. 5.1 and 5.2, the assumption of constant dust temperature T_dust = 10 K might not be appropriate. We hence computed M_core via Eq. (10) at three different temperatures, that is, 5, 10, and 20 K, which cover the range of values. Similarly to Table 3, we estimate a 33% relative error on the masses. We highlight again that the lack of total-power data (which is, however, not offered in continuum observations with ALMA) may cause a partial underestimation of the total core masses, in the hypothesis that part of the extended emission is filtered out. The resulting parameters are summarised in Table B.1.

The continuum-core sizes are similar to those of the H₂D⁺-identified ones and present an average effective radius of 2.2 × 10³AU. In general, the cores present low masses, unless low dust temperature values are assumed, which, however, are supported by our o- H₂D⁺ analysis. In AG351, three cores are expected to be more massive than 10 M_⊙ under the assumption of T_dust = 5 K. Core c-4 and core c-5 were also identified in the o-H₂D⁺ data. Core c-6, instead, does not present bright emission in the molecular line, which may be due either to high depletion of H₂D⁺, or to the presence of proto-stellar feedback from an unidentified proto-stellar object, as discussed in Sect. 5.3.

Fig. B.1 Background image shows the dust thermal emission in Band 7, with the continuum-identified structures overlaid as grey contours. The black crosses show the position of the flux peak within each core. The white contours show the H2D+ -identified cores resulting from the SCIMES analysis.

Continuum cores in AG354 are on average more massive than in AG351. In fact, six out of nine cores present M_core ≥ 10 M_⊙ at 5 K. Core c-9, which corresponds to the H₂D⁺-identified core 2, appears as massive as 60 M_⊙. It is worth noticing that core c-9 contains two emission peaks, separated by 2.5″, the fainter of which has a peak flux approximately one-third of the brighter one. In the dendrogram analysis, however, they are never separated, despite the choice of input parameters, most likely due to the fact that they lay closer than three ALMA beams within each other. Core c-8 is in a similar situation as core c-6 in AG351, since it does not correspond to bright molecular emission, and hence could be in a later evolutionary stage.

Overall, the dendrogram analysis confirms our findings: most cores are essentially low mass, even though we stress again that we may be underestimating masses from the continuum emission due to the lack of zero-spacing observations. Furthermore, the correlation between continuum- and H₂D⁺-identified structures ispoor, since the majority of cores either do not overlap completely in the two datasets, or they do so for less than 70% of their projected area. There are, however, a few exceptions, the most important of which is represented by core c-9, which, with an estimated mass budget of several tens of solar masses and detectable o- H₂D⁺ emission, represents an ideal HMPC candidate.

Appendix C Parameter maps of individual cores

Figure C.1 shows the maps of the best-fit parameters obtained with MCWEEDS in each core.

Fig. C.1 Left panels: maps of σV, in km s−1 units. Central panels: Maps of the molecular column density of o-H2D+, in unit of log10(cm−2). Right panels: maps of Vlsr, in km s−1 units. Each row shows a different core, labelled at the top of the central panel. The beam size of the H2D+ observations is shown in the bottom right corners. The coordinates are shown in the ICRS system. The dashed contours show the continuum emission at levels [1, 5, 10, 15, 20] mJy beam−1.

References

All Tables

All Figures

	Fig. 1 Maps of the continuum emission at 870 μm seen with APEX by the ATLASGAL survey in AG351 (left) and AG354 (right). The FoV of the ALMA observations is indicated with a red circle. The beam size (18.2″) is shown in the bottom left corner, while the scale bar is in the bottom right.
In the text

Fig. 2 Map of the continuum emission as seen by ALMA in Band 7 in AG351 (left panel) and AG354 (right panel). The contours show the o-H2D+-integrated intensity emission from the original not re-gridded and not primary-beam-corrected data, at levels = [5, 7, 9, 11]σ, where 1σ = 60 mK km s−1 (AG351) and 70 mK km s−1 (AG354). For the continuum data, we show the primary-beam-corrected maps. The beam size and the scale bar are shown in the bottom left and bottom right corners, respectively.

In the text

	Fig. 3 Contours show the leaves identified by SCIMES, overlaid on the primary-beam-corrected continuum map in AG351 (left panel) and AG354 (right panel).
In the text

	Fig. 4 Maps of o-H2D+ column density (top left panel), σV (top right), and centroid velocity (bottom) obtained from the MCWEEDS analysis in the individual cores in AG351. The beam size and scale bar are shown in the bottom left and bottom right corners, respectively. The black contours show the continuum emission in Band 7 at levels [1, 5, 10, 15, 20] mJy beam−1.
In the text

	Fig. 5 Same as Fig. 4, but for AG354.
In the text

	Fig. 6 Normalised kernel density distribution of σV and Ncol in AG351 (blue) and AG354 (red) obtained from the MCWEEDS pixel-per-pixel fit. The contour levels are [0.1, 0.3, 0.5, 0.7, 0.9]. The horizontal solid line shows the sound speed at 10 K, while the dot-dashed line represents the thermal broadening of H2D+ at the same temperature.
In the text

	Fig. 7 Mean signal obtained averaging in each core the observed spectra is shown in black. The red curve shows the best-fit model obtained with MCWEEDS. In the top-right corner of each panel we report the core id number and the Vlsr value to which the spectra have been aligned, before averaging them.
In the text

	Fig. 8 Same as Fig. 7, but for AG354.
In the text

Fig. 9 Scatterplots of o-H2D+(11,0−11,1) velocity dispersion σV with respect to H2 column density N(H2) in AG351 (blue, top panels) and AG354 (red, bottom panels). In the left panels, N(H2) is computed assuming a constant temperature of 10 K. In the right panels, we use the temperature derived pixel-per-pixel from the o- H2D+(11,0−11,1) line width, in the hypothesis of only thermal contributions. In order to allow an easier comparison between the two clumps, the plot ranges are kept equal in the corresponding panels. For the sake of readability,we only show data-points for which the relative 95% HPD interval derived with MCWEEDS is < 50%. The black stars show the average values referring to the H2D+-identified cores, for which the continuum flux is significantly detected. The σV values are taken from Table 2, while the N(H2) values are computed as averages in each core.

In the text

Fig. 10 Scatterplots of o-H2D+ molecular abundance with respect to H2 column density in AG351 (blue, top panels) and AG354 (red, bottom panels). In the left panels, the quantities are computed assuming a constant temperature of 10 K. In the right panels, we use pixel-per-pixel the temperature derived from the o- H2D+(11,0−11,1) line width, in the hypothesis of only thermal contributions. Scatterplot ranges are fixed to allow an easier comparison.For the sake of readability, we only show data points for which the relative 95% confidence interval is < 50%. The dashed black curve shows the correlation found by Sabatini et al. (2020). The black stars show the average values referring to the H2D+-identified cores, for which significant continuum flux is detected. The Ncol values to compute Xmol(o-H2D+) are taken from Table 2, and the N(H2) values are computed as averages in each core.

In the text

	Fig. 11 Xmol(o-H2D+) distribution in core 2 in AG354, the most massive structure identified in this work. Xmol (o-H2D+) was computed adopting a constant Tdust = 10 K. The contours show the ALMA continuum flux in Band 7, at levels [5, 10, 15] mJy beam−1. Under the assumption of constant dust temperature, the contours also represent the distribution of N(H2).
In the text

Fig. A.1 Comparison of the o-H2D+(11,0−11,1) spectra observed by APEX (red) and ALMA (black) in AG351 and AG354, from left to right. For a fair comparison, the APEX spectra from Sabatini et al. (2020) were converted in flux units using the APEX gain (40 Jy K−1). The ALMA spectra were computed integrating the signal over an area corresponding to the APEX beam size (16.8″), and they were smoothed to the APEX spectral resolution.

In the text

	Fig. B.1 Background image shows the dust thermal emission in Band 7, with the continuum-identified structures overlaid as grey contours. The black crosses show the position of the flux peak within each core. The white contours show the H2D+ -identified cores resulting from the SCIMES analysis.
In the text

Fig. C.1 Left panels: maps of σV, in km s−1 units. Central panels: Maps of the molecular column density of o-H2D+, in unit of log10(cm−2). Right panels: maps of Vlsr, in km s−1 units. Each row shows a different core, labelled at the top of the central panel. The beam size of the H2D+ observations is shown in the bottom right corners. The coordinates are shown in the ICRS system. The dashed contours show the continuum emission at levels [1, 5, 10, 15, 20] mJy beam−1.

In the text