Identification of prestellar cores in high-mass star forming clumps via $\rm H_2D^+$ observations with ALMA

Context. The different theoretical models concerning the formation of high-mass stars make distinct predictions regarding their progenitors, i.e. the high-mass prestellar cores. However, so far no conclusive observation of such objects has been made. Aims. We aim to study the very early stages of high-mass star formation in two infrared-dark, massive clumps, to identify the core population that they harbour. Methods. We obtained ALMA observations of continuum emission at 0.8mm and of the ortho-$\rm H_2D^+$ transition at 372GHz towards the two clumps. We use the SCIMES algorithm to identify cores in the position-position-velocity space, finding 16 cores. We model their observed spectra in the LTE approximation, deriving the centroid velocity, linewidth, and column density maps. We also study the correlation between the continuum and molecular data, which in general do not present the same structure. Results. We report for the first time the detection of ortho-$\rm H_2D^+$ in high-mass star-forming regions performed with an interferometer. The molecular emission shows narrow and subsonic lines, suggesting that locally the temperature of the gas is less than 10K. From the continuum emission we estimate the cores' total masses, and compare them with the respective virial masses. We also compute the volume density values, which are found to be higher than $10^{6}\, \rm cm^{-3}$. Conclusions. Our data confirm that ortho-$\rm H_2D^+$ is an ideal tracer of cold and dense gas. Interestingly, almost all the $\rm H_2D^+$-identified cores are less massive than 13M_sun , with the exception of one core in AG354. Furthermore, most of them are subvirial and larger than their Jeans masses. These results are difficult to explain in the context of the turbulent accretion models, which predict massive and virialised prestellar cores.


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, these play a key role in the energetic budget of the ISM. Furthermore, there is evidence that our Sun was born in a cluster containing also 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 socalled high-mass prestellar 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 highmass star forming clumps, 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, due to several reasons. High-mass stars are intrinsically rarer and more short-lived with respect to low-mass counterparts, as predicted by stellar evolution theories. As a consequence, they are on average more distant, which in turn affect the achievable spatial resolution of the observations. Moreover, they form in crowded and dense environments, heavily affected by extinction (e.g. Zinnecker & Yorke 2007). This means that the lack of infrared emission detected with single-dish facilities is not a conclu-Article number, page 1 of 25 arXiv:2104.06431v1 [astro-ph.GA] 13 Apr 2021 sive evidence of prestellar stage (e.g. Motte et al. 2018). In this context, interferometers such as the Atacama Large Millimeter/submillimeter Array (ALMA) represent a powerful tool, being able to provide the necessary angular resolution and sensitivity to resolve substructures in distant high-mass clumps.
In the context of the search for HMPCs, targeting 70µmdark clumps with interferometric 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 investigate with ALMA band 6 observations two HMPCs candidates in the rich high-mass star forming regions W43-MM1, previously investigated by Nony et al. (2018). One of the two cores does not show clear evidence of protostellar 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 9 cm −3 ) and mass (M 30 M ), it lacks line emission in several, abundant species such as N 2 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 prestellar, 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 ASHES survey (Alma Survey of 70 µm dark High-mass clumps in Early Stages; Sanhueza et al. 2019;Li et al. 2020) targeted twelve IRDCs dark 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 1 ≈ 70% of them are classified as prestellar. None of them appears 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-/proto-stellar. Continuum observations at mm/sub-mm wavelengths (especially when only one frequency is available) cannot provide information on the mass and on the temperature independently, and they cannot unveil the sources 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 4 cm −3 ), 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 reaction: which is particularly efficient at low temperatures (assuming a low ortho-to-para H 2 ratio; see e.g. Pagani et al. 1992), due to the lower zero-point energy of H 2 D + . From this, and from its doubly and triply deuterated forms D 2 H + and D + 3 , all other deuterated molecules in the gas phase (such as N 2 D + and DCO + ) are formed (Dalgarno & Lepp 1984).
In a hunt for HMPCs, Tan and collaborators have used an extensive dataset of N 2 D + observations with ALMA. Tan et al. (2013) reported initially the detection of six cores, the most massive of which was revealed to be in fact composed by several cores, some of which already in a protostellar stage Kong et al. 2018). The sample was increased to 141 N 2 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 protostellar 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 2 D + /N 2 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. (2019Feng et al. ( , 2020 reported similar findings using N 2 H + and HCO + isotopologues. Despite the fact that N 2 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 molecule containing carbon and/or oxygen, but at high densities (n 10 6 cm −3 ) N 2 H + and N 2 D + show sign 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 prestellar gas, since they can persist in the gas phase at high densities. H 2 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 2 D + ) has a critical density of n cr ≈ 10 5 cm −3 (Hugo et al. 2009). Furthermore, H 2 D + is very sensitive to the temperature, due to the chemical reactions with its main destroyer, carbon monoxide. When -due for instance to protostellar activity-the gas temperature raises beyond 20 K, the CO molecules frozen onto the dust grains evaporate back into the gas phase, thus lowering the H 2 D + abundance. Moreover, above 30 K reaction (1) start to proceed backwards, thus reducing the deuteration level. Hence, detecting H 2 D + is an unambiguous sign of prestellar conditions in the sense that core structures which show significant emission in o-H 2 D + are trustfully in a prestellar phase.
Studies of H 2 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. (2003Caselli et al. ( , 2008) produced a small sample of low-mass cores observed in o-H 2 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 2 D + emission in the star forming region DR21, in Cignus X, using JCMT. Brünken et al. (2014) demonstrated the use of H 2 D + as a chemical clock for low-mass star forming regions. More recently, Giannetti et al. (2019) detected o-H 2 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 anticorrelates with that of N 2 D + . In a following study, Sabatini et al. (2020) published a detailed census of this tracer in 16 highmass 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 2 D + with interferometers is present in literature.
In this work, we report the first interferometric observation of o-H 2 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. Us-E. Redaelli et al.: Prestellar cores in HM clumps identified with H 2 D + ing a core-identifying algorithm, we find a total of 16 cores in H 2 D + . We study their properties using a spectral fitting analysis. Interestingly, a significant fraction of the gas emitting in o-H 2 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 part of the clumps is still in a very cold, prestellar stage. We study 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 2 D + emission and continuum emission. They represent ideal candidates to be prestellar 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 2 D + data-cubes and we perform the line spectral analysis. We discuss the results in Sect. 5: we first focus on the properties of the gas traced by the o-H 2 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.

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 (∼ 10000) of massive clumps in dif-ferent evolutionary stages (Schuller et al. 2009). In particular, the two chosen clumps belong also to the ATLASGAL TOP100 subsample (Giannetti et al. 2014;König et al. 2017), and their ID 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 (∼ 500M ) obtained for the 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 several infrared 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 µmdark, since they are not associated to detected point-like sources at 24 and 70 µm. We have furthermore 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 protostellar 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 prestellar phases, even though we highlight that several works showed that 70µm-dark clumps were later discovered to contain protostars with intereferometric observations Li et al. 2019Li et al. , 2020.
Both clumps are included in the sample investigated by Sabatini et al. (2020), who detected the o-H 2 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 −1 . In AG351 and AG354 the line peak temperature is T peak ≈ 200 mK and the linewidth is ≈ 1.0 − 1.2 km s −1 . Using the H 2 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 2 D + ) ≈ 10 −10 in both clumps.
The spectral setup consists of four spectral windows (SPWs). One is centred on the o-H 2 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 −1 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 2 and methanol lines , with a total bandwidth of 938 MHz and spectral resolution of 564 kHz (≈ 0.47 km s −1 ). These lines however were not detected. At these frequencies, and with the used configuration, the maximum recoverable scale is θ MRS ≈ 20 , the primary beam size is θ FoV ≈ 26 2 , and the resolution is ≈ 1.0 (corresponding to ≈ 1600 AU at the sources average distance of 1.6 kpc). The total observing time were 180 min (ACA) and 28 min (12m-array) for each source. During the observations, the precipitable water vapour was typically 0.3 mm < PWV < 0.7 mm, and in general lower than 1 mm. The average system temperature values are found in the range 300−400 K for the SPW containing the o-H 2 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 datacubes of the o-H 2 D + line, we focused on the central 25 MHz 2 At 372 GHz, the primary beams of the main array and of ACA are 17 and 30 , respectively (≈ 20 km s −1 ) around the source centroid velocity. The datacubes hence present 100 channels in the frequency axis, and the final spectral resolution is 0.22 km s −1 . Finally, in order to avoid oversampling, both the continuum and the line images (already primary-beam corrected) have been re-gridded in order to ensure 3 pixel per beam minor axis, in agreement with the Nyquist theorem. The molecular line data have been converted into the brightness temperature T b scale, computing the corresponding gain, G = 11 mK/(mJy beam −1 ) 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. 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 2 D + (1 1,0 − 1 1,1 ) transition, obtained integrating the cubes over the range [−4.3; −2.0] km s −1 (AG351) and [−7.5; −4.0] km s −1 (AG354). The continuum maps present a sensitivity of 0.8 mJy/beam for both sources, whilst the median root-mean-squared (rms) of the line data, computed over emission-free channels, is 300 mK per 0.22 km s −1 . These values refer to already primary-beam corrected data. The observation resolutions and sensitivities for each source are summarised in Table 1.

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

Core identification with scimes
In order to identify substructures in the molecular emission we use the Spectral Clustering for Interstellar Molecular Emission Segmentation package (scimes, Colombo et al. 2015). scimes is based on dendrograms algorithms (Rosolowsky et al. 2008), and it is 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 2 D + emission. We therefore focus on the leaves, which in the context of this paper correspond to prestellar cores.
In order to first build the dendrogram, a number of parameters is needed. We performed multiple tests, followed by visual inspection of the results, to find the best choice for our o-H 2 D + data. First of all, following Rosolowsky & Leroy (2006), we find the regions in the emission characterised by signal-tonoise ratio (S/N) higher than a given threshold (S/N lim ), which contain peaks of emission brighter than a second higher threshold (S/N peak ). This implementation maximises the recovered information in case of moderate-to-low signal-to-noise-ratio data (S/N = 3−10), such as our o-H 2 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 E. Redaelli et al.: Prestellar cores in HM clumps identified with H 2 D +  Notes. (a) The beam size is expressed as: major axis × minor axis, position angle (PA).
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/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 apply 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 signal-to-noise ratio map. We set the minimum number of channels that a leaf must span to N min chan = 2, due to the low linewidth of the o-H 2 D + line. Finally, we exclude 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 identify seven cores in AG351 and nine cores in AG354, which we have labelled in order of decreasing peak value of o-H 2 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 cores overlap (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 which overlap spatially but present different centroid velocities (i.e. distinct velocity components along the line of sight). The first four columns of Table regular; ii) the limited S/N of the o-H 2 D + data, which translates in 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 in simulations cores do not always appear regular in shapes (see e.g. Smullen et al. 2020). Figure 2 shows that the o-H 2 D + integrated intensity distribution does not always follow the structure seen in continuum emission. The cores just identified confirm this scenario: in some cases, dust thermal emission peaks are found within molecularidentified cores (e.g. core 3 in AG351, or core 2 in AG354). In other cases, bright continuum spots do not overlap with H 2 D + emission. In Sec. 5.2 and Appendix B we further investigate this point.

Spectral fitting with MCWeeds
In order to derive the physical parameters of each core from the o-H 2 D + data, we perform a spectral fitting of the lines pixel-perpixel 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 use 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 3 . This improves the code speed when dealing with large datasets, allowing to fit simultaneously spectra from multiple positions. For each spectrum, the code performs 100000 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 was in any case checked.
Since our ALMA data comprise only one o-H 2 D + transition, they do not allow to constrain independently the molecular column density (N col ) and the excitation temperature (T ex ). We therefore assume 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 prestellar and protostellar cores; Friesen et al. (2014) adopted T ex = 12 K when studying prestellar cores in Ophiuchus A. The remaining free parameters in the fit are hence the molecular column density N col (o-H 2 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 are chosen to improve the code convergence. In order to estimate correctly the parameters distributions, the noise level of the spectral line data has to be indicated. We adopt rms = 180 mK (AG351) and rms = 170 mK (AG354) 4 .
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 component on the line of sight, to self-absorption, or to large scale emission filtering are present in the observed o-H 2 D + spectra. Even though they represent a small minority of the observations, these spectra cannot be successfully fitted with a single velocitycomponent, LTE model. Since these positions are characterised by an overestimation of the linewidth, we find that a successful masking strategy is to discard fits with FWHM higher than a given threshold. This threshold, which is chosen for each core individually, is in the range 1 − 1.2 km s −1 . The fraction of rejected pixels is 3% for both clumps, and we checked that no singlecomponent 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. We however convert 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: 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, whilst 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 Fig. 4 and 5 show that the cores are coherent structures in velocity, with dispersions around mean values of usually ≈ 0.2 km s −1 . At the clump scale, stronger velocity gradients of ≈ 1 km s −1 and ≈ 2 km s −1 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 to infer the large-scale gas kinematics (see Appendix A for more details).
Typical values for the molecular column density are in the range N col = (0.6 − 4) × 10 13 cm −2 . These are higher than those found by Sabatini et al. (2020) using APEX observations, most likely due the dilution of the o-H 2 D + emission in the large single-dish beam (16.8 ). The observed values are on average higher also than the maximum column density reported by Friesen et al. (2014) with ALMA at a higher spatial resolution (N col = 1.2 × 10 13 cm −3 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 have estimated the optical depth τ maps, using the spectral parameters listed in the Cologne Database for Molecular Spectroscopy (CDMS 5 ), through the equation: where A ul = 1.20 × 10 −8 s −1 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 2 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 is not affecting the line shapes strongly, even though we cannot completely exclude this hypothesis given the scarce spectral resolution of ALMA observations.
The velocity dispersion maps of the clumps show that ALMA detects very narrow o-H 2 D + lines. In several cores, especially in AG351, large portions of the gas traced by H 2 D + present total linewidths narrower than the isothermal sonic speed, which at T gas = 10 K is 6 : 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 2 and 10% of helium). In particular, in AG351 36% of positions show subsonic linewidths 7 , whilst 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 linewidths. This has profound implication on the physical conditions traced by o-H 2 D + (1 1,0 − 1 1,1 ), as discussed in detail in Sect. 5.1. Locally, the line velocity dispersion increases up to 0.4−0.5 km s −1 , without however showing a clear correlations with other quantities, such as the emission in continuum.
We have compared the distribution of the velocity dispersion and column density in the two clumps in Fig. 6. These plots show that AG351 presents in general narrower lines and lower column density values than AG354, which could hint to differences in their evolutionary stage and/or physical conditions.

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 have computed the average o-H 2 D + (1 1,0 − 1 1,1 ) spectrum in each core. As highlighted in the previous subsection, the clumps present velocity gradients of the order of 1 − 2 km s −1 , while the individual cores show smaller variations of the centroid velocity. Nevertheless, given the narrow linewidths of the o-H 2 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 have then been re-analysed using MCWeeds, as previously described, setting T ex = 10 K and fitting the linewidth, centroid velocity, and column density of o-H 2 D + . We have 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 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 selfabsorption) 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 turns lead 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 subor super-sonic. The total velocity dispersion values derived from the average spectra (Table 2, 7th column) show that about half of cores present subsonic or trans-sonic gas motions, within uncertainties. We can further disentangle the thermal and non-thermal components of the lines. In fact, the velocity dispersion of any line is due to a combination of its thermal broadening and nonthermal one, 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. In the assumption that the two components are independent, they sum in quadrature to compose the total, observed velocity dispersion: σ V = σ 2 V,NT + σ 2 V,th (Myers 1983). The thermal component of the o-H 2 D + line at the gas temperature T gas is: where m H 2 D + is the H 2 D + molecular mass in g. We have computed the ratio of non-themal 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 −1 . 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).

The gas traced by o-H 2 D +
As explained in Sect. 1, the o-H 2 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 to determine the gas temperature and density independently at the same time. In order to perform the spectral analysis, we have 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 2 D + is abundant when the temperature is 20 K. However, dense prestellar gas can reach temperatures lower than 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 linewidth of the o-H 2 D + (1 1,0 − 1 1,1 ) transition can give us hints on the gas temperature. As mentioned in the previous Sect., the velocity dispersion of any line can be separated in thermal and non-thermal broadening, which are summed in quadra-11. ture. It is therefore straightforward that the total linewidth of a line cannot be smaller than its thermal component. The thermal broadening of H 2 D + at 10 K is σ V,th = 0.14 km s −1 . From the maps of this quantity obtained in Sec. 4.2 it results that 17% of pixels in AG351 and 7% in AG354 show linewidths smaller than this level (see Fig. 6). This is a clear indication that at least some parts of the gas traced by o-H 2 D + have a temperature lower that 10 K. We can reverse this argument, and derive the gas temperature T H 2 D + at which the observed velocity dispersion values are purely due to thermal broadening: We have hence derived maps of T H 2 D + pixel-per-pixel in each source. This quantity represents an upper limit on the real gas temperature, since 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 T H 2 D + = 6 K from APEX observations of o-H 2 D + (1 1,0 − 1 1,1 ) in the starless clump Ophiuchus D. In our clumps, we constrain T H 2 D + to be in the range 5 − 20 K. The lower limit is due to the fact that i) temperature 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 linewidth (in FWHM) equals the spectral resolution of our observations. The upper limit is instead due to the fact that the H 2 D + abundance is expected to drop at temperatures higher than 20 K (which corresponds to σ V,th = 0.21 km s −1 ). Linewidths broader than this limit most likely involve non-thermal contributions, such as turbulence. A significant fraction of positions presents T H 2 D + < 20 K (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): where G is the gravitational constant. Equation (6) holds in the assumption of a spherical core of radius R core and uniform density, and σ dyn represents the total linewidth of the gas, and it can be derived from the observed σ V of o-H 2 D + assuming that the non-thermal component of the latter is the same for all the gas components: In a first step, we estimate M vir from the o-H 2 D + linewidths and the effective radii listed in Table 2 (R core = R eff ), using a constant gas temperature of 10 K. For each core we compute the distribution of the virial mass from the corresponding distribution of σ V , and we then determine the median and the 95% HPD intervals. The resulting virial masses are summarised in Table 3.
As previously discussed, we have reasons to believe that at least part of the gas traced by o-H 2 D + is colder than 10 K. As a test, we have therefore computed the virial masses using the minimum T min H 2 D + value in each core as derived from Eq. (5). With some little algebra it can be shown that σ dyn (and hence M vir ) is positively correlated with temperature. Since T H 2 D + 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 T H 2 D + 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 % from the previous ones. It is interesting to notice that the estimated virial masses are in the range 0.3 − 2.2 M , and hence these cores are essentially lowmass, in the hypothesis that they are virialised by the kinetic energy. We further discuss this point in Sec. 5.2 and 5.4.

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 2 D + transition. As previously mentioned, the two sets of data do not present the same morphology. We investigate this point in de-  0.6 core 7 -3.3 km s 1 Fig. 7. The 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 V lsr value to which the spectra have been aligned, before averaging them. tail in Appendix B, where we identify core-like structures in the continuum maps using dendrogram analysis. Out of the 16 H 2 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 2 D + -identified cores the continuum flux peak is found outside the core boundaries. It is hence clear that the two datasets are not tracing in general the same material in the clumps. In order to further investigate this point, we have used the continuum maps to estimate the total gas mass and density. In conditions of high densities (n 10 4−5 cm −3 ), dust and gas are well coupled, which means that they are thermalised at the same temperature (Goldsmith 2001). We can hence assume T dust = T gas = 10 K in order to compute the total gas column density N(H 2 ) and hence the mass of the identified cores. We compute pixel-by-pixel the quantity: 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 −1 ) and A pix is the pixel area; µ H 2 = 2.8 is the mean molecular weight per hydrogen molecule , and κ ν is the dust opacity at the frequency of the observations. For the latter, we use the power-law expression: in which we employ β = 1.5 for the dust emissivity index (Mezger et al. 1990;Walker et al. 1990) and the κ 0 = 10 cm 2 g −1 for the dust opacity at the reference frequency corresponding to the wavelength λ 0 = 250 µm (Hildebrand 1983;Beckwith et al. 1990). Errors on N(H 2 ) are computed with standard propagation from the uncertainties on continuum flux, and the average uncertainty is rms = 2.9 × 10 22 cm −2 . We have computed the average N(H 2 ) values for each core, which span the range (1 − 4) × 10 23 cm −2 . From the gas column density we derive the molecular abundance X mol The scatterplots of the o-H 2 D + linewidth and molecular abundance with respect to H 2 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 2 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 2 ) 5 × 10 23 cm −2 ) are characterised by narrow linewidths (σ V 0.35 km s −1 ). The molecular abundance seems anti-correlated with respect to N(H 2 ) for AG354, whilst 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 suggested also by the fact that AG354 shows on average larger velocity dispersion and o-H 2 D + column density values. However, the long tail towards high gas densities in AG354 arises from the data points of one core mainly (core 2), and hence we cannot derive strong conclusions regarding the whole clump based on this.
The anticorrelation between X mol and N(H 2 ) can be explained by the fact that at later evolutionary stages several factors can contribute to lowering the H 2 D + abundance: depletion of HD onto the dust grains (Sipilä et al. 2013), conversion of H 2 D + into its doubly and triply deuterated isotopologues, or a rise in temperature due to the protostellar feedback. Regarding the lat-ter point, we lack the interferometric data to verify whether any protostellar 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 2 D + , using: 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 have 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, as for instance the dust opacity. A variation of 30% of κ ν translates in a variation of 33% in the mass. We have therefore increased the error on M core to this level, in order to take into account the parameters uncertainties. The obtained masses range between 0.4-6 M , with larger values observed in AG354. However, these values might not represent the entire mass budget of the H 2 D + -identified cores, particularly for the ones where the continuum emission is weak and Table 3. Properties of the identified cores: virial mass, total mass from dust continuum emission, and average volume density. The first set of values is computed in the assumption that T gas = T dust = 10 K; in the second half of the table, the quantities are computed using the minimum temperature derived from the H 2 D + thermal broadening. The used temperature values are listed in the 5th column.
The uncertainties on the masses estimated from the continuum emission are 33% (relative error), to take into account uncertainties on parameters such as the dust opacity, on top of the flux error from the ALMA observations. (b) T H 2 D + is the minimum T H 2 D + in each core computed through Eq. 5. It still represents, however, an upper limit on the real gas/dust temperature, since other components (non-thermal, instrumental,...) can contribute to the line broadening.
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 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 2 D + emission traces part of the gas colder than 10 K, locally also the dust temperature can be lower than this value. We have therefore re-computed the N(H 2 ) (and hence the o-H 2 D + abundance) pixel-per-pixel using the gas temperature maps T H 2 D + obtained from the o-H 2 D + (1 1,0 − 1 1,1 ) linewidths, as discussed in Sec. 5.1. The average uncertainties on the new N(H 2 ) values are rms = 2.1 × 10 22 cm −2 for AG351 and rms = 1.4×10 22 cm −2 for AG354, and they are computed as standard error propagation from the flux uncertainty. The core average column density values are N(H 2 ) = 0.3 − 4 × 10 23 cm −2 . The new scatterplots showing the correlation of σ V and X mol (H 2 D + ) with respect to N(H 2 ) 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 however expected, since the narrower the line, the lower the derived T H 2 D + , and hence the higher N(H 2 ) (since, for a given flux, dust temperature and mass/density are anti-correlated). On the other hand, the anti-correlation between the molecular abundance and the H 2 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 2 D and dust thermal emission in several high-mass star forming regions. Concerning o-H 2 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 2 D + (1 1,0 − 1 1,1 ) in a sample of ATLASGAL sources to fit the relation between X mol (o-H 2 D + ) and N(H 2 ), which in the log-log scale reads: log 10 [X mol (o-H 2 D + )] = −1.06 × log 10 [N(H 2 )] + 13.93. We plot the predictions from this relation in Fig. 10. We have fit a linear relation in the log-log space between the molecular abundance and the H 2 column density for our data. We find slopes of m ≈ −0.8, −0.9 with high correlation coefficient (|R| > 0.85), with the exception of the X mol (o-H 2 D + )-N(H 2 ) 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: 16.8 2 /(2.5 2 × 10) ≈ 5.
It is important to notice that X mol (o-H 2 D + ) and N(H 2 ) 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 2 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 have checked the correlation between the molecular column density N col (o-H 2 D + ) and N(H 2 ) (not shown here). These two quantities are in fact independent 8 . We do not find any significant correlation between these variables, as previously noted by Sabatini et al. (2020). The o-H 2 D + column density is flat with respect to N(H 2 ). This confirms that higher N(H 2 ) values present lower molecular abundances.
We have computed the core masses using the T min H 2 D + of each core, as we did for the virial masses. We stress again that this temperature value is the minimum T H 2 D + value found within each core, but it represents an upper limit, due both to possible nonthermal components to the linewidth, 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 in the assumption of very low dust temperatures (T dust ≈ 5 K), 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 in the assumption of T dust = 5 K, one of which however do not present significant o-H 2 D + emission, hence preventing us from classifying it as truly prestellar. We further discuss these findings in Sec. 5.4. It is worth noting though the lack of total-power observations means that the ALMA data might be filtering out the large scale emission associated with the gas surrounding the cores here identified, hence leading to 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, in the assumption of constant distribution and spherical geometry, using the equation: We have evaluated n(H 2 ) for both cases: assuming T dust = 10 K and using T dust = T min H 2 D + derived from the o-H 2 D + (1 1,0 − 1 1,1 ) linewidth. 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 6 cm −3 . This confirms that the ALMA data are tracing high-density material, and justifies the hypothesis of dust-and-gas coupling. Furthermore, at such high densities the o-H 2 D + (1 1,0 − 1 1,1 ) is expected to be thermalised, and hence LTE is a good approximation (Harju et al. 2008). The found volume densities correspond to surface densities in the range Σ = (0.1 − 10.0) g cm −2 , 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 is not enough alone to balance the gravitational pull, and if no other source of pressure is present, they will collapse. On the other hand, supervirial sources (α vir > 1) are unlikely to undergo gravitational contraction. We have computed the virial ratio in both cases, with a constant temperature of 10 K and with T = T min H 2 D + , and we report the two sets of values in Table 3. Most cores are subvirial, or consistent with being subvirial within the 95% HPD intervals in AG351, regardless of the assumed temperature. In AG354, the supervirial cases are represented by core 5, 6, and 7 (i.e. 33% of the sample). This is consistent with what found by Sabatini et al. (2020) at the clump scales, who reported that AG351 is more subvirial than AG354 (α vir = 0.4 and 0.8, respectively). We cannot exclude however than 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.

A population of cores at different evolutionary stages
We can speculate that the different correlation between the dust continuum emission and the H 2 D + distribution in the cores is linked to their evolutionary stage. In fact, H 2 D + forms when the gas become 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 2 D + abundance. H 2 D + can be transformed in D 2 H + and D + 3 or deplete as a consequence of HD depletion (Sipilä et al. 2013). Secondly, one has to take into account opacity effects: if the volume density grows much larger than the line critical density, the transition becomes optically thick (τ > 1), which means that it does not trace anymore the whole bulk of gas, but only the layer at τ ≈ 1. In Sect. 4.2 we have derived the optical depths from the derived N col maps, and despite being a minority, there are positions where lines are moderately optically thick. This could contribute to explain 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 prestellar core collapses and forms a protostellar object, the abundance of H 2 D + is finally reduced by the rising temperature.
We can therefore speculate that cores with strong H 2 D + emission, either corresponding to a continuum-identified structure or not, are early in their evolution and truly prestellar. The lack of detected continuum peaks and associated continuumidentified cores could be due to observational limits. In fact, if the dust temperature in the correspondence of the o-H 2 D + emission is lower than in the surroundings -as expected in case of prestellar 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 which do not have a correspondence in continuum-identified structures could still be centrally peaked in reality.
As the cores evolve, the H 2 D + abundance at the dust continuum peak is lowered for the aforementioned reasons (e.g. the conversion to D 2 H + and D + 3 ). H 2 D + -identified cores which partially overlap with continuum structure, but which present a con- Fig. 9. Scatterplots of o-H 2 D + (1 1,0 − 1 1,1 ) velocity dispersion σ V with respect to H 2 column density N(H 2 ) in AG351 (blue, top panels) and AG354 (red, bottom panels). In the left panels, N(H 2 ) is computed assuming a constant temperature of 10 K. In the right panels, we use the temperature derived pixel-per-pixel from the o-H 2 D + (1 1,0 − 1 1,1 ) linewidth, 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 H 2 D + -identified cores, for which the continuum flux is significantly detected. The σ V values are taken from Table 2, whilst the N(H 2 ) values are computed as averages in each core. tinuum flux peak outside their boundaries, may belong to this evolutionary stage. Cores 4 and 6 in AG354 are representative of this case (as shown in Appendix B). The fact that they do not present significant increase of the velocity dispersion in the proximity of the continuum peaks suggest that the decrease of X mol (o-H 2 D + ) is most likely due either to HD depletion onto dust grains or to conversion into other isotopologues. Finally, cores which are only visible in continuum are the most evolved. The lack of detection of H 2 D + in their surroundings suggests that either their temperature is overall higher than 20 K, hinting to an already protostellar stage, or that H 2 D + is depleted onto dust grains and/or converted into its doubly and triply deuterated forms.
In order to test our hypothesis, we have searched for correlations of the suggested evolutionary phases and the cores dynamical properties reported in Table 2 and 3, such as linewidths, virial masses, and average densities. We did not find any significant trend. Nevertheless, Fig. 10 shows a decrease of X mol (o-H 2 D + ) with an increase of N(H 2 ), that could resemble a situation where the collapse is more advanced at least in some of the cores. Further observations of tracers of protostellar activity (such as SiO jets or CO outflows) or kinematics tracers could confirm this hypothesis.
According to our analysis, core 3 and 7 (AG351) and core 2 (AG354) represent good candidates to be truly prestellar cores in high-mass clumps, since they overlap significantly with continuum-identified structures with masses M ≥ 14 M . Fig.  11 shows the distribution of the o-H 2 D + abundance in core 2 in AG354, the most massive structure we identify, with contours from the continuum (thus representative of the H 2 column density, in the assumption of constant dust temperature). It can be seen that the molecular distribution seems anticorrelated with the continuum one: at the dust peak X mol (o-H 2 D + ) present the lowest values, whilst the molecular abundance in- Fig. 10. Scatterplots of o-H 2 D + molecular abundance with respect to H 2 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-H 2 D + (1 1,0 − 1 1,1 ) linewidth, 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 H 2 D + -identified cores, for which the continuum flux is significantly detected. The N col values to compute X mol (o-H 2 D + ) are taken from Table 2, whilst the N(H 2 ) values are computed as averages in each core. creases moving away from it, in a similar fashion to the general trend shown by the scatterplots in Fig. 10. This core shows particularly low linewidths ( σ V = 0.20 km s −1 ), which become as low as 0.10 km s −1 towards the core centre, where hence 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 hence represents an ideal candidate of HMPCs.
The evolution here discussed can either be chemical or physical (i.e. a density evolution). The first case holds if the cores are virialised by other 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 2 D + is efficiently formed, and hence the cores show bright molecular emission. As time passes by, the o-H 2 D + abundance decreases, because the molecule is depleted onto the dust grains and/or converted into other isotopologues, until the o-H 2 D + (1 1,0 − 1 1,1 ) at the continuum peaks is not detectable anymore.
In the case the cores are instead subvirial, and hence they are experiencing gravitational contraction, the density is evolving and increasing with time. The evolution of X mol (o-H 2 D + ) is however 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 2 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. E. Redaelli et al.: Prestellar cores in HM clumps identified with H 2 D + Fig. 11. The X mol (o-H 2 D + ) distribution in core 2 in AG354, the most massive structure identified in this work. X mol (o-H 2 D + ) has been computed adopting a constant T dust = 10 K. The contours show the ALMA continuum flux in band 7, at levels [5, 10, 15] mJy beam −1 . In the assumption of constant dust temperature, the contours represent also the distribution of N(H 2 ).

Our results in the context of star formation theories
The analysis of our data shows that the H 2 D + -identified cores have masses in the range 0.1 − 13 M depending on the assumed temperature. These values, which are probably accurate within a factor of ≈ 5 (see Appendix A), would hence suggest that several of the H 2 D + -identified cores are essentially low-mass. The only clear exception is core AG354-2 for which, in 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 2 D + -identified cores, 80% of them are consistent with α vir 1, regardless of the assumed temperature. The virial parameter computed in the previous section however takes into account only the kinetic energy content (T ) and the gravitational energy (U), whilst in a more complete fashion one has to take into account also the magnetic energy term, which in case of a homogeneous and spherical core threaded by a uniform magnetic field B reads: where the equation is expressed in cgs units. We can determine the strength of the magnetic field needed to halt the gravitational collapse by imposing the 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 5 − 10 6 cm −3 the authors report values of B ≈ a few hundreds µG, but never higher than 1 mG. However, analysing the polarised dust thermal emission, several authors found magnetic fields of the order of the 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 provides 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 of further observations to investigate the level of magnetisation of the clumps, and their large-scale kinematics. We have also estimated the Jeans masses for the H 2 D +identified cores. For volume densities in the range n = 10 6 − 10 7 cm −3 and temperatures of 5−10 K, the Jeans mass is M Jeans = 0.03 − 0.3 M . Hence, almost the totality of the H 2 D + -identified cores contains 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 analysis we have performed on the linewidths show that the motions of the gas traced by o-H 2 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 highmass stars as in the competitive accretion model (Bonnell & Bate 2006), growing during the low-mass protostellar stage similarly to what also proposed by Tigé et al. (2017); Motte et al. (2018). We have highlighted how AG351 and AG354 sits in the lowend of the clump masses in the ATLASGAL sample. Kauffmann & Pillai (2010) investigated empirically the threshold for highmass star formation in galactic IRDCs, 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 M lim clump = 260 M . 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 clumps population (M ≈ 500 M ), can host several high-mass prestellar cores. Further studies focused not only on the continuum emission but also on H 2 D + on more massive clumps are then needed to assess the existence of massive prestellar cores and disentangle between the different theories.
Up to now, in literature the best candidates of HMPCs have been the C1S core dentified 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 2 D + and in a follow-up paper a non-detection of o-H 2 D + . 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, protostellar cores seen by Tan et al. (2016). The source investigated by Cyganowski et al. (2014), on the other hand, has been identified only in continuum, which prevents from 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 subvirial dynamical states, in agreement with what we report in this study.
The main advantage of this work is the detection of o-H 2 D + , which traces the very early stages, and we can then state these are possibly the first prestellar 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 the their dynamical states. In addition, numerical studies including chemistry are needed to shed light on the different proposed theories and disentangle among the different physical processes that can affect the star-formation process. Examples of these works, including detailed deuteration chemistry, have recently started to be developed (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.

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 2 D + (1 1,0 − 1 1,1 ) in this kind of sources with interferometric observations. Our molecular line data show that o-H 2 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 have identified 16 cores in the o-H 2 D + datacubes. We have 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 2 D + (1 1,0 − 1 1,1 ) lines, with linewidths lower than their thermal broadening at 10 K. This indicates that this species traces very cold and quiescent region in the analysed sources.
We have 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 2 D + is quite abundant, and tracing the high-density gas. The presence of cores 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 Sec. 5.3 for more details). Later on, as the gas becomes denser, or as the chemistry evolves, the X mol (o-H 2 D + ) starts to decrease, but the molecule is still detectable. Eventually, the H 2 D + abundance drops, due either to depletion at very high densities, or to protostellar feedback effects.
Our results highlight how the continuum emission alone is generally not a good probe of prestellar gas. Bright cores in continuum flux which do not show significant emission in cold-gas tracers (such as H 2 D + ) are likely in a more evolved, possibly protostellar stage. This suggests that particular care must be taken when doing surveys of "prestellar" cores seen only in continuum, especially when only one frequency band is available, since these data do not allow 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 2 D + represents a good tracer of cold gas at densities of n ≈ 10 6 cm −3 . This gas is in a prestellar phase, in the sense that it has not be influenced by protostellar activity, and -being cold and dense-it has the potential to form new protostars; its evolution, however, is determined by its dynamic state.
At densities higher than 10 6 cm −3 even o-H 2 D + could be not ideal anymore, as shown by its abundance drop as a function of N(H 2 ), due to opacity effects, depletion, or a combination of both. In those physical conditions, D 2 H + and D + 3 represent probably the only good tracer available. Whilst the latter is not observable, the former could represent, in combination with o-H 2 D + , the best choice to investigate the more evolved and denser regions.
Most of the H 2 D + -identified cores are less massive that 10 M , even in 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, hence leading to 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, protostellar phase. Two of the cores contains peaks of both molecular and continuum emission, which allows us to reliably estimate their total masses. Core AG354-2 shows particularly narrow linewidths, 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 hence represents an ideal candidate of HMPCs, but further investigation is needed to better constrain its temperature and thus mass.
We have investigated the dynamic state of the H 2 D +identified cores by means of the virial analysis. Most of the cores appear subvirial if we take into account the gravitational energy and the kinetic energy only. We however do not have information on the magnetic properties of the sources. We estimate that magnetic field strengths of the order of several hundreds of mi-croGauss or a few milliGauss are needed to virialise the cores. Further observations aimed to recover the magnetic field properties are needed to make a definite conclusion on the core dynamic states.
As future perspective, we plan to recover the missing flux in the large scale emission of o-H 2 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 in turn will provide key information about the cores' dynamics. Furthermore, ALMA observations at a similar spatial resolution of molecules which are good tracers of the protostellar activity (SiO, CO,...) will help us disentangle unambiguously prestellar cores from protostellar ones. also filtered out, hence leading to mass underestimation, as observed in the nearby low-mass prestellar 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. Hence also the clump total masses could be underestimate, 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 2 − 5.

Appendix B: Continuum-identified cores
In the main text we have focused on the properties of the H 2 D + -identified cores. However, it is worth looking also to the continuum-identified structures, which is the purpose of this Appendix. scimes is developed to work in position-position-velocity space, and hence it is suitable to analyse molecular line data only. We have hence used the python package astrodendro (http: //www.dendrograms.org/), 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 have 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 −1 as the rms level of the continuum maps before the correction for the primary beam response.
We find in total 15 cores, which are shown in Fig. B.1, together with the H 2 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 2 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, or just outside, of the H 2 D +identified cores. In Table B.1 we report the label of the corresponding H 2 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 range 0 − 70%.
Similarly to what done for the H 2 D + -identified cores, we have 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 have hence computed M core via Eq. (10) at three different temperatures, i.e. 5, 10, and 20 K, which cover the range of T H 2 D + 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 2 D +identified ones, and present an average effective radius of 2.2 × 10 3 AU. In general, the cores present low masses, unless low dust temperature values are assumed, which however are supported by our o-H 2 D + analysis. In AG351, three cores are expected to be more massive than 10 M in the assumption of T dust = 5 K. Core c-4 and core c-5 were identified also in the o-H 2 D + data. Core c-6, instead, do not present bright emission in the molecular line, which may be due either to high depletion of H 2 D + , or to the presence of protostellar feedback from an unidentified protostellar object, as discussed in Sect 5.3.
Continuum cores in AG354 are on average more massive than in AG351. Six out of nine cores present in fact M core ≥ 10 M at 5 K. Core c-9, which correspond to the H 2 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 they however 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 2 D + -identified structures is poor, since the majority of cores either do not overlap completely in the two datasets, or they do 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 2 D + emission, represents an ideal HMPC candidate.