A timeline for massive star-forming regions via combined observation of o-H$_2$D$^+$ and N$_2$D$^+$

TL;DR

This study investigates how the abundances of o-H$_2$D$^+$ and N$_2$D$^+$ vary with the evolutionary stages of high-mass star-forming regions, proposing their ratio as an effective tracer of star formation progression.

Contribution

It provides observational evidence linking the relative abundances of o-H$_2$D$^+$ and N$_2$D$^+$ to the evolutionary stages of massive clumps, enhancing understanding of chemical evolution in star formation.

Findings

o-H$_2$D$^+$ abundance decreases with evolution

N$_2$D$^+$ abundance increases with evolution

Estimated age of the youngest clump is about 10^5 years

Abstract

Context: In cold and dense gas prior to the formation of young stellar objects, heavy molecular species (including CO) are accreted onto dust grains. Under these conditions H$_3^+$ and its deuterated isotopologues become more abundant, enhancing the deuterium fraction of molecules such as N$_2$H$^+$ that are formed via ion-neutral reactions. Because this process is extremely temperature sensitive, the abundance of these species is likely linked to the evolutionary stage of the source. Aims: We investigate how the abundances of o-H$_2$D$^+$ and N$_2$D$^+$ vary with evolution in high-mass clumps. Methods: We observed with APEX the ground-state transitions of o-H$_2$D$^+$ near 372 GHz, and N$_2$D$^+$(3-2) near 231 GHz for three massive clumps in different evolutionary stages. The sources were selected within the G351.77-0.51 complex to minimise the variation of initial chemical conditions, and to remove distance effects. We modelled their dust continuum emission to estimate their physical properties, and also modelled their spectra under the assumption of local thermodynamic equilibrium to calculate beam-averaged abundances. Results: We find an anticorrelation between the abundance of o-H$_2$D$^+$ and that of N$_2$D$^+$, with the former decreasing and the latter increasing with evolution. With the new observations we are also able to provide a qualitative upper limit to the age of the youngest clump of about 10$^5$ yr, comparable to its current free-fall time. Conclusions: We can explain the evolution of the two tracers with simple considerations on the chemical formation paths, depletion of heavy elements, and evaporation from the grains. We therefore propose that the joint observation and the relative abundance of o-H$_2$D$^+$ and N$_2$D$^+$ can act as an efficient tracer of the evolutionary stages of the star-formation process.

Tables

Table 1. Table 1: Properties of the clumps derived from dust continuum emission. The clumps go from the least- to the most evolved.

Source	$T_d$	$N\left(H_2\right)$	a𝑎aa𝑎a$a$Within the central $28″$, see text.	a𝑎aa𝑎a$a$Within the central $28″$, see text.	$M$b𝑏bb𝑏b$b$Rescaled for the new dust temperatures from the values in Leurini et al. (2011).	$Diameter$c𝑐cc𝑐c$c$From Leurini et al. (2011).	$n\left(H_2\right)$
	K	${10}^{22}\text{usk}{\mathrm{cm}}^{-2}$	$L_{\odot }$	$M_{\odot }$	$M_{\odot }$	pc	${10}^5\text{usk}{\mathrm{cm}}^{-3}$
Clump 7	13.0	7.3	23	31	120	0.18	1.3
Clump 5	15.5	8.7	111	35	100	0.10	2.9
Clump 2	20.0	10.5	331	44	200	0.19	1.8

Table 2. Table 2: Results of the line-fitting procedure with MCWeeds.

Source	$N\left(\text{o-}{H}_2{D}^{+}\right)$		$N\left(N_2{D}^{+}\right)$		$V_{\mathrm{LSR}}\left(\text{o-}{H}_2{D}^{+}\right)$		$V_{\mathrm{LSR}}\left(N_2{D}^{+}\right)$		$\Delta V\left(\text{o-}{H}_2{D}^{+}\right)$		$\Delta V\left(N_2{D}^{+}\right)$
	${10}^{11}\text{usk}{\mathrm{cm}}^{-2}$		${10}^{10}\text{usk}{\mathrm{cm}}^{-2}$		$\mathrm{km}\text{usk}{s}^{-1}$		$\mathrm{km}\text{usk}{s}^{-1}$		$\mathrm{km}\text{usk}{s}^{-1}$		$\mathrm{km}\text{usk}{s}^{-1}$
	med.	$95\%$ CI	med.	$95\%$ CI	med.	$95\%$ CI	med.	$95\%$ CI	med.	$95\%$ CI	med.	$95\%$ CI
Clump 2	$\dots$		20.4	$17.0;24.6$	–3.0	$\dots$	–2.4	$-2.5;-2.3$	1.0	$\dots$	1.6	$1.3;1.9$
Clump 5	21.4	$13.9;29.6$	12.0	$8.6;16.0$	–2.6	$-2.9;-2.3$	–2.7	$-3.0;-2.3$	1.8	$1.1;2.6$	2.3	$1.7;3.2$
Clump 7	33.3	$22.4;47.2$	$\dots$		–3.8	$-4.0;-3.5$	–3.0	$\dots$	1.5	$1.0;2.6$	1.0	$\dots$

Equations

$${\mathrm{N}}_2+{\mathrm{H}}_2{\mathrm{D}}^{+}\to {\mathrm{N}}_2{\mathrm{D}}^{+}+{\mathrm{H}}_2,$$

$${\mathrm{H}}_3^{+}+\mathrm{HD}\leftrightarrows {\mathrm{H}}_2{\mathrm{D}}^{+}+{\mathrm{H}}_2+\Delta E$$

Full text

*Aims. *We investigate how the abundances of o-$\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$ vary with evolution in high-mass clumps.

*Methods. *We observed with APEX the ground-state transitions of o-$\mathrm{H}_{2}\mathrm{D}^{+}$ near $372\usk\mathrm{G}\mathrm{Hz}$, and $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2) near $231\usk\mathrm{G}\mathrm{Hz}$ for three massive clumps in different evolutionary stages. The sources were selected within the G351.77–0.51 complex to minimise the variation of initial chemical conditions, and to remove distance effects. We modelled their dust continuum emission to estimate their physical properties, and also modelled their spectra under the assumption of local thermodynamic equilibrium to calculate beam-averaged abundances.

*Results. *We find an anticorrelation between the abundance of o-$\mathrm{H}_{2}\mathrm{D}^{+}$ and that of $\mathrm{N}_{2}\mathrm{D}^{+}$, with the former decreasing and the latter increasing with evolution. With the new observations we are also able to provide a qualitative upper limit to the age of the youngest clump of about 105 yr, comparable to its current free-fall time.

*Conclusions. *We can explain the evolution of the two tracers with simple considerations on the chemical formation paths, depletion of heavy elements, and evaporation from the grains. We therefore propose that the joint observation and the relative abundance of o-$\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$ can act as an efficient tracer of the evolutionary stages of the star-formation process.

††offprints: A. Giannetti (agianne@ira.inaf.it)

1 Introduction

Depletion of heavy elements proceeds progressively in dense and cold material ($n_{\mathrm{H_{2}}}>\mathrm{few}\times 10^{4}\usk\mathrm{c}\mathrm{m}^{-3}$, $T\lesssim 25\usk\mathrm{K}$), rendering heavy molecular species unobservable, and thus not viable to trace gas in the earliest phases of the star-formation process. It has been demonstrated that deuterated species become more abundant under these conditions, and that there is a correlation between the degree of deuteration and the depletion of CO (Caselli et al., 1999, 2002; Bacmann et al., 2003). Large values for the deuterium fraction are also found towards protostellar sources (e.g. Ceccarelli et al., 1998; Parise et al., 2002).

Pivotal to the deuteration process is the exothermic reaction chain $\mathrm{H_{3}^{+}}\rightarrow\mathrm{H_{2}D^{+}}\rightarrow\mathrm{D_{2}H^{+}}\rightarrow\mathrm{D_{3}^{+}}$ (e.g. Walmsley et al., 2004; Flower et al., 2004). $\mathrm{H}_{2}\mathrm{D}^{+}$, $\mathrm{D_{2}H^{+}}$, and $\mathrm{D_{3}^{+}}$ react with other species, producing the vastly enhanced deuterium fractions in molecules at the core of the ion-neutral chemistry. $\mathrm{N_{2}}$ in particular can react with $\mathrm{H}_{2}\mathrm{D}^{+}$ to form $\mathrm{N}_{2}\mathrm{D}^{+}$via

[TABLE]

and more efficiently with $\mathrm{D_{2}H^{+}}$ and $\mathrm{D_{3}^{+}}$. The molecule $\mathrm{N}_{2}\mathrm{D}^{+}$ is frequently used as a reliable tracer for gas affected by heavy CO depletion (e.g. Caselli et al., 2002; Barnes et al., 2016). Because $\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$ are chemically linked, and because they are abundant only under specific circumstances, they could be used to determine the evolutionary stage of a source in early phases (e.g. Caselli et al., 2008; Emprechtinger et al., 2009; Fontani et al., 2015), something of particular interest in the high-mass regime.

Pillai et al. (2012) have obtained maps for part of the DR21 complex of the o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}) and $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2) transitions, observed with the James-Clerk-Maxwell Telescope (JCMT) and the Submillimeter Array (SMA), respectively. They find very extended o-$\mathrm{H}_{2}\mathrm{D}^{+}$ emission as do Vastel et al. (2006) in the low-mass source L1544. They also find that this species mainly traces gas that is not seen in dust continuum emission or in the interferometric $\mathrm{N}_{2}\mathrm{D}^{+}$ data. $\mathrm{H}_{2}\mathrm{D}^{+}$ may therefore be sensitive to gas that would elude detection in the most commonly used tracers, and can represent an even earlier stage in the process of star formation.

These results raise two important questions that we try to address in this work: Firstly, is the lack of correlation between o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}) and $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2) found by Pillai et al. (2012) real, or is it an effect of interferometric filtering? Secondly, how do the abundances of these molecules evolve with time, and why? In this Letter we first describe our observations, present the results, and finally discuss the importance of our findings.

2 Source selection and observations

G351.77–0.51 (hereafter G351) is the closest and most massive filament ($D=1\usk\mathrm{kpc}$, $M\sim 2000\usk\mathrm{M_{\odot}}$, Leurini et al., subm.) identified in the $870\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ APEX Telescope Large Survey of the Galaxy (ATLASGAL) (Schuller et al., 2009). A three colour image of the region is shown in Fig. 1, emphasising the large amount of dense and cold gas in its massive ridge.

G351 offers several advantages for our purposes: 1) it is close-by, 2) it hosts massive clumps that exceed the Kauffmann & Pillai (2010) threshold for high-mass star formation, and 3) the clumps are in different evolutionary stages. Therefore, sources along the spine of the G351 filament are extremely well suited to studying the evolution of chemistry in high-mass clumps, not only because they are nearby and all at the same distance, but also because they have similar properties and share the same initial chemical conditions.

We select three clumps in this complex (Clumps 2, 5, and 7, following Leurini et al., 2011) in different evolutionary stages, based on their IR properties: Clump 7 is the least evolved, being still quiescent at $70\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$; Clumps 5 and 2 are both bright at this wavelength, with the former being weaker than the latter (Fig. 1). These sources are the ones with the largest column densities and masses found in G351 in each evolutionary stage. Their mass and peak column densities are also within a factor of two, using a temperature of $25\usk\mathrm{K}$ for Clump 2, and $10\usk\mathrm{K}$ for Clumps 5 and 7 (cf. Table 4 in Leurini et al. 2011).

The dust emission peak for each of the three clumps has been targeted by APEX 12m submillimeter telescope observations (project M-099.F-9508A), using the PI230 and FLASH+ receivers, to cover the $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2) and o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}) lines at $231.3216$ and $372.42134\usk\mathrm{G}\mathrm{Hz}$, respectively. Observations were performed between 2017 May 31 and 2017 September 19, and completed on 2018 July 1-2. The rms noise on the main-beam brightness temperature scale, $T_{\mathrm{{MB}}}$, is $\sim 20\usk\mathrm{m}\mathrm{K}$ at $231\usk\mathrm{G}\mathrm{Hz}$ and $30\usk\mathrm{m}\mathrm{K}$ at $372\usk\mathrm{G}\mathrm{Hz}$, for a spectral resolution of $0.45\usk\mathrm{k}\mathrm{m}\usk\mathrm{s}^{-1}$. We converted the antenna temperature $T_{A}^{*}$ to $T_{\mathrm{{MB}}}$, using $\eta_{\mathrm{MB}}=0.60$ for o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}) and $\eta_{\mathrm{MB}}=0.69$ for $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2).

3 Results

In our APEX observations Clumps 5 and 7 were detected in o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}), while Clumps 2 and 5 were detected in $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2). Figure 2 shows the spectra for all sources observed, highlighting their different features in line emission.

To give a physical explanation of this behaviour, we first characterise the properties of the clumps.

3.1 Dust spectrum

We perform aperture photometry on the clumps to extract dust continuum fluxes. We use images from ATLASGAL ($870\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ Schuller et al., 2009), HiGAL ($350\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, $250\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, $160\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and $70\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ Molinari et al., 2010), MIPSGAL ($24\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ Carey et al., 2009) and MSX (Egan et al., 2003). All images were smoothed to $28\arcsec$, the resolution of the APEX 230 GHz observations, and the fluxes were extracted in the central beam, after removing the median background estimated in an annulus with $2\times\theta_{beam}<R_{annulus}<5\times\theta_{beam}$, for each clump. The apertures and annuli are centred on the positions of the clumps, as determined in Leurini et al. (2011).

The dust temperature and column density are estimated via a greybody fit at $\lambda\geq 70\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$. The bolometric luminosity is computed as the sum of the integral of the best-fit greybody curve and the integral at $\lambda\leq 70\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, computed with the trapezoidal rule in log-log space. For the greybody we use $\kappa_{870\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}}=1.85\usk\mathrm{c}\mathrm{m}^{2}\usk\mathrm{g}^{-1}$ and $\beta=1.75$, and consider the $70\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ flux as an upper limit.

To convert the dust column density to that of molecular hydrogen, we use a gas-to-dust ratio $\gamma=120$, as calculated from Eq. 2 in Giannetti et al. (2017a), and $R_{GC}=7.4\usk\mathrm{kpc}$ (Leurini et al., subm.). The impression that Clump 7 is the least evolved of the three sources considered, and that Clump 2 is the most evolved one is confirmed by $T_{\mathrm{{d}}}$ and the $L/M$ ratio (see Table 1), two efficient indicators of evolution (e.g. Saraceno et al., 1996; Molinari et al., 2008; König et al., 2017; Urquhart et al., 2018).

The average volume density along the line of sight listed in Table 1 is estimated as $n_{\mathrm{H_{2}}}=N_{\mathrm{H_{2}}}/\ell$, where $\ell$ is the size of the clump obtained by Leurini et al. (2011). We note however that the filament width is nearly constant ($\sim 0.2\usk\mathrm{pc}$; Leurini et al., subm.), even at the location of the clumps, so the difference in mean volume density may be overestimated.

3.2 Molecular column densities and abundances

To calculate the column densities, a spectral-line fit is performed, under the assumption of local thermodynamic equilibrium (LTE) (cf. Vastel et al., 2012), with MCWeeds (Giannetti et al., 2017b). This allows the uncertainty on this parameter to be obtained as well (Table 2). We use the partition function for o-$\mathrm{H}_{2}\mathrm{D}^{+}$ (in the relevant temperature range, $9.375\usk\mathrm{K}$: $10.3375$, $18.750\usk\mathrm{K}$: $12.5068$, $37.500\usk\mathrm{K}$: $15.5054$) from CDMS (Müller et al., 2001). Because the observed clumps are sufficiently dense to attain thermal coupling between gas and dust (cf. Table 1), and because we lack an estimate of the excitation temperature for o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}) and $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2), we use $T_{\mathrm{{d}}}$, assuming that $T_{\mathrm{{d}}}=T_{\mathrm{{ex}}}$. The fit is performed using Monte Carlo Markov Chains. We used an adaptive Metropolis-Hastings sampler (Haario et al., 2001), with 100,000 total samples, a burn-in period and a delay for the adaptive sampling of 10,000 and 5,000 iterations, respectively, and a thinning factor of 20 (see Giannetti et al., 2017b, for more details on these parameters). Convergence and independence of the samples are ensured with the Raftery-Lewis (Raftery & Lewis, 1995), Geweke (Geweke, 1992a) and Gelman-Rubin (Geweke, 1992b) tests.

We include in the budget a Gaussian calibration uncertainty with $\sigma=5\%$. Multiple tests have been performed on the priors to make sure that the choice of the latter is not crucial for the results.

Beam-averaged abundances can be estimated from the peak column density of H2 derived from the dust continuum emission. In the following, the same value of $N_{\mathrm{H_{2}}}$ is used for both o-$\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$, despite the difference in the angular resolution of the spectral-line observations ($18-28\arcsec$, respectively). Smoothing the continuum maps to the resolution of the ATLASGAL data (excluding the HiGAL $350\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ image) leads to an increase in the peak H2 column density of $\lesssim 15\%$, well within the uncertainties; this has the effect of lowering $X(\mathrm{H}_{2}\mathrm{D}^{+})$ in Clump 5, thus increasing the difference with Clump 7 (see Fig. 3, left).

If $X(\mathrm{o-}\mathrm{H}_{2}\mathrm{D}^{+})$ does not vary significantly over scales of 18000-28000 AU ($18-28\arcsec$ at $1\usk\mathrm{kpc}$), one can also estimate the relative abundance of the two species. In Fig. 3 one can see that o-$\mathrm{H}_{2}\mathrm{D}^{+}$ becomes rarer with evolution, while $\mathrm{N}_{2}\mathrm{D}^{+}$ shows the opposite behaviour. Considering the abundance variation and its uncertainty, o-$\mathrm{H}_{2}\mathrm{D}^{+}$ is much more sensitive to the clump evolution than $\mathrm{N}_{2}\mathrm{D}^{+}$. The relative abundance (panel (b)) of the two species offers another way of looking at the same finding, but with the advantage of being independent of the H2 column density, then removing an additional source of uncertainty. From the figure, we see that the o-$\mathrm{H}_{2}\mathrm{D}^{+}$/$\mathrm{N}_{2}\mathrm{D}^{+}$ ratio progressively decreases in time and o-$\mathrm{H}_{2}\mathrm{D}^{+}$ is more abundant than $\mathrm{N}_{2}\mathrm{D}^{+}$ by a factor $\gtrsim 10$ in the youngest Clump 5 and 7.

4 Discussion

In the previous section we show that o-$\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$ have opposite trends in abundance as a function of evolution. In the following we discuss a few possible explanations for the abundance differences between the three clumps.

4.1 Deuteration of H${}_{3}^{+}$ and availability of gas-phase N2

$\mathrm{H}_{2}\mathrm{D}^{+}$ rapidly forms in cold and dense gas, where the reaction

[TABLE]

only proceeds in the forward direction, unless there is a substantial fraction of ortho-$\mathrm{H_{2}}$ (e.g. Gerlich et al., 2002), the only case in which $\Delta E$ can be negative. $\mathrm{H}_{2}\mathrm{D}^{+}$ then drives the deuteration of the neutral molecules that have a deuteron affinity larger than molecular hydrogen (Millar et al., 1989). The o-$\mathrm{H}_{2}\mathrm{D}^{+}$ abundance derived in Clumps 5 and 7 is of the order of $\sim 3\times 10^{-11}$, similar to the values reported by Pillai et al. (2012) for their o-$\mathrm{H}_{2}\mathrm{D}^{+}$ peaks in the DR 21 region, and by Vastel et al. (2006) for the outskirts of L1544. For the temperatures of our sources, $13\usk\mathrm{K}$ and $15.5\usk\mathrm{K}$ for Clumps 7 and 5, respectively, the ortho-to-para ratio (OPR) predicted by Flower et al. (2004) for $\mathrm{H}_{2}\mathrm{D}^{+}$ is in the range $0.05-0.1$. Recently, we also performed 3D magnetohydrodynamic simulations of collapsing clumps and filaments, which included chemistry under the assumption of total depletion (Körtgen et al., 2017, 2018). In these simulations we also follow the OPR for $\mathrm{H}_{2}\mathrm{D}^{+}$, finding that it rapidly decreases to values below unity, in agreement with Flower et al. (2004), on time scales of the order of $10^{5}\usk\mathrm{yr}$. This OPR can be used to estimate the total $\mathrm{H}_{2}\mathrm{D}^{+}$ abundance, which would be in the range $10^{-10}-10^{-9}$. Caselli et al. (2003) show how in L1544, $\mathrm{H}_{2}\mathrm{D}^{+}$ is nearly as abundant as electrons ($10^{-9}\mathrm{\,vs.\,}2\times 10^{-9}$), indicating that this species is a major ion, which only happens when the depletion is extremely high. Our results, when compared to these findings, also suggest that regions exist in the clump where heavy elements are vastly depleted, which is also supported by our ongoing study of CO depletion for the entire filament (Sabatini et al., in prep.).

Clump 7, however, not only has a high abundance of o-$\mathrm{H}_{2}\mathrm{D}^{+}$, but also shows low values of $X(\mathrm{N}_{2}\mathrm{D}^{+})$. One possible explanation for this is the time lag that is needed to form $\mathrm{N}_{2}\mathrm{D}^{+}$ from $\mathrm{H}_{2}\mathrm{D}^{+}$. As an example, Figure 1 in Sipilä et al. (2015) shows that in the first $10^{5}\usk\mathrm{yr}$, according to their model, o-$\mathrm{H}_{2}\mathrm{D}^{+}$ is much more abundant than $\mathrm{N}_{2}\mathrm{D}^{+}$, and the difference progressively decreases. This provides, qualitatively, an upper limit to the age of the clump of 105 years, comparable to the current free-fall time ($t_{ff}\sim 1.2\times 10^{5}$ yr). The clump is indeed sub-virial ($\alpha\sim 0.4$333Assuming a homogeneous clump (see MacLaren et al., 1988), the radius from Table 1, and the line width from o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}).), and shows a very weak stellar activity in its centre. Clump 7 is also the only one to show a blue-skewed profile in $\mathrm{HCO^{+}}$(3–2), which indicates that the clump is collapsing (see Fig. 4). We note that if the $\mathrm{H_{2}}$ ortho-to-para conversion occurs on the surface of dust grains, the chemical timescale could be even shorter (Bovino et al., 2017).

In Clump 5, on the other hand, $X(\mathrm{N}_{2}\mathrm{D}^{+})$ is two to three times larger than in Clump 7. In this clump, being more evolved, $\mathrm{N}_{2}\mathrm{D}^{+}$ has had time to form, while the stellar activity is not yet sufficient to warm up the gas to temperatures high enough to significantly alter the chemistry on large scales.

Interestingly, Clump 2 still has a large fraction of $\mathrm{N}_{2}\mathrm{D}^{+}$ with respect to $\mathrm{H_{2}}$, but at the same time $X(\mathrm{H}_{2}\mathrm{D}^{+})$ is at least an order of magnitude lower than $X(\mathrm{N}_{2}\mathrm{D}^{+})$. Because $\mathrm{N}_{2}\mathrm{D}^{+}$ forms via deuteron exchange from $\mathrm{H}_{2}\mathrm{D}^{+}$, one would expect to see a correlation in the abundance of the two molecules, after they have had time to form. Furthermore, Clump 2 hosts luminous young stellar objects (YSOs), possibly close to the zero-age main sequence (ZAMS), as indicated by its $L/M$ ratio (Giannetti et al., 2017b), and is therefore a relatively evolved source.

What is keeping the level of $\mathrm{N}_{2}\mathrm{D}^{+}$ high, even higher than in Clump 5, and at the same time destroying o-$\mathrm{H}_{2}\mathrm{D}^{+}$? The beam- and line-of-sight-averaged depletion of CO does not vary significantly from Clump 7 to Clump 2, and a relatively high depletion is one of the key ingredients for the efficient formation of $\mathrm{N}_{2}\mathrm{D}^{+}$ because carbon monoxide acts as a destroyer of both $\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$ (e.g. Walmsley et al., 2004; Emprechtinger et al., 2009). The other fundamental ingredient for the formation of $\mathrm{N}_{2}\mathrm{D}^{+}$ is the presence of $\mathrm{N_{2}}$ in the gas phase, which can react with $\mathrm{H}_{2}\mathrm{D}^{+}$, $\mathrm{D_{2}H^{+}}$, and $\mathrm{D_{3}^{+}}$. Close to the YSOs, where the temperature is high and CO is in the gas phase, $\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$ are both actively destroyed and Eq. 2 is also fast in the backward direction. However, at larger distances from luminous YSOs, the material is still dense and cold, and CO remains frozen out onto dust grains, contrary to $\mathrm{N_{2}}$, either because the latter has not been significantly depleted (due to its longer formation timescale compared to CO and the fact that atomic nitrogen could have a lower sticking coefficient than $\mathrm{N_{2}}$ Flower et al., 2006), or because it evaporates faster. In fact, according to the most recent estimates (Wakelam et al., 2017) $\mathrm{N_{2}}$ and CO binding energies differ by a few hundred Kelvin (1100 K vs 1300 K), enough to cause a significant difference in the evaporation timescales. The latter, evaluated by employing the standard $\tau_{\mathrm{evap}}=\nu_{0}^{-1}\mathrm{e}^{E_{D}/k_{B}T_{\mathrm{d}}}$, are $\sim 5\times 10^{8}\usk\mathrm{yr}$, and $\sim 2\times 10^{4}\usk\mathrm{yr}$, for CO and N2, respectively. In the equation above, $\nu_{0}=10^{12}\mathrm{s}^{-1}$ is the typical harmonic frequency (Hasegawa et al., 1992), $E_{D}$ is the binding energy, $k_{B}$ is the Boltzmann’s constant, and $T_{\mathrm{d}}$ is the dust temperature (see Table 1). More importantly, with time, a large amount of $\mathrm{H}_{2}\mathrm{D}^{+}$ is converted to $\mathrm{D_{2}H^{+}}$ and $\mathrm{D_{3}^{+}}$, decreasing its abundance. $\mathrm{N}_{2}\mathrm{D}^{+}$, on the other hand, is more efficiently formed by the reaction of $\mathrm{D_{2}H^{+}}$ and $\mathrm{D_{3}^{+}}$ with $\mathrm{N_{2}}$, compared to $\mathrm{H}_{2}\mathrm{D}^{+}$. The combination of these processes could succeed in reproducing the observed reduction in abundance of o-$\mathrm{H}_{2}\mathrm{D}^{+}$, and in maintaining a high column density of $\mathrm{N}_{2}\mathrm{D}^{+}$. In even more evolved sources, $X(\mathrm{N}_{2}\mathrm{D}^{+})$ decreases, as expected (Fontani et al., 2015).

If this is the case, we propose that the combined abundances of o-$\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$ represent a more efficient evolutionary indicator for the first stages of the high-mass star-formation process, compared to their individual values (e.g. Caselli et al., 2008; Fontani et al., 2015), in the same way that they are both needed to trace the full reservoir of gas (Pillai et al., 2012).

4.2 Alternative explanations

One alternative possibility to explain the emission in Clump 7 could be that the region in the clump where heavy elements (including $\mathrm{N_{2}}$) are completely depleted is extended, strongly affecting the abundance measurement. Contamination in the beam and along the line-of-sight would explain the weak emission of high-density tracers towards this source (see also Leurini et al., 2011, for $\mathrm{N_{2}H^{+}}$ $J=1\rightarrow 0$).

The critical density of the o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}) is $\sim 10^{5}\usk\mathrm{c}\mathrm{m}^{-3}$, roughly one order of magnitude lower than that of $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2), and similar to the mean values of $n(\mathrm{H_{2}})$ reported in Table 1. Although this is not always a direct indication that emission is coming from denser gas (see e.g. Kauffmann et al., 2017), it is possible that in Clump 7 the region dense enough to excite $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2) is smaller than in the other two clumps.

In both cases, the abundance of o-$\mathrm{H}_{2}\mathrm{D}^{+}$ relative to $\mathrm{N}_{2}\mathrm{D}^{+}$ keeps its potential as an evolutionary indicator, on the one hand because if a region of complete depletion exists at the centre of the clump, the warm-up by the YSOs will progressively erode it, and on the other hand because density increases as a result of the clump collapse, before feedback becomes important at the clump scale and starts to dissipate it.

For Clumps 2 and 5, the $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2) emission is likely to arise in cold gas surrounding the star-forming cores, as observed by Fontani et al. (2009) for example. If o-$\mathrm{H}_{2}\mathrm{D}^{+}$ is abundant in these compact structures, its emission may be completely saturated and diluted in the APEX beam for Clump 2; o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}) must also be weak in the lower-density gas enveloping these cores. This is in contrast with the observations of extended emission in DR 21, and those of L1544 (Vastel et al., 2006), although the envelope of high-mass star forming regions could be warmer due to the stronger external illumination from close-by high-mass stars; this increases the ortho-to-para $\mathrm{H_{2}}$ ratio, which suppresses D-fractionation (e.g. Kong et al., 2015). Additional high-resolution observations are needed to investigate these possibilities.

5 Summary and conclusions

We observed the o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}) and $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2) lines in three clumps along the spine of G351, the most massive filament within $1\usk\mathrm{kpc}$ (Leurini et al., 2011) in ATLASGAL. These observations provide for the first time the possibility to investigate the variation with evolution of the abundances of o-$\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$ in high-mass sources, from a clump that is still quiescent at $70\usk\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, to one which hosts luminous YSOs, likely close to the ZAMS. The selected clumps not only belong to the same complex, ensuring the best conditions for a comparison of the chemistry, but also have a comparable peak column density of $\mathrm{H_{2}}$.

The abundance of $\mathrm{N}_{2}\mathrm{D}^{+}$ progressively increases with evolution, while an opposite trend is found for o-$\mathrm{H}_{2}\mathrm{D}^{+}$. We propose that the chemical evolution of the clumps causes this behaviour, and that the relative abundance of these species is a good indicator of time evolution.

First, after the clump reaches a high-enough density to allow for an efficient CO depletion, $\mathrm{H}_{2}\mathrm{D}^{+}$ starts to form, but to attain a considerable amount of $\mathrm{N}_{2}\mathrm{D}^{+}$ would require a longer time (e.g. Sipilä et al., 2015). Clump 7 still has $X(\mathrm{N}_{2}\mathrm{D}^{+})\lesssim 9\times 10^{-13}$, which also suggests an age of no more than 105 yr, and a fast collapse time. As time passes, more $\mathrm{N}_{2}\mathrm{D}^{+}$ forms via reaction 1 in cold and dense gas. High degrees of deuteration are commonly observed in relatively young and massive clumps (e.g. Caselli et al., 2002; Fontani et al., 2011). This is the situation of Clump 5, where both $\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$ are abundant, despite the presence of YSOs (Fig. 1), which are not yet luminous enough to warm up a significant portion of the gas in the clump. In Clump 2, $X(\mathrm{N}_{2}\mathrm{D}^{+})$ increases by a factor of $\sim 2$ and $X(\mathrm{H}_{2}\mathrm{D}^{+})$ decreases by a factor $\gtrsim 10$ compared to Clump 5. Because this source is the most evolved, $\mathrm{N_{2}}$ has had more time to form. Also, the line-of-sight and beam-averaged dust temperature in Clump 2 is $\sim 20\usk\mathrm{K}$, not high enough to release CO into the gas phase on large scales (see Sabatini et al., in prep.), but sufficient to reduce the evaporation timescale of $\mathrm{N_{2}}$ to $\sim 2\times 10^{4}\usk\mathrm{yr}$. Therefore, $\mathrm{N_{2}}$ is abundant in the gas phase throughout the clump, whether or not it was depleted from the gas phase earlier. $\mathrm{H}_{2}\mathrm{D}^{+}$ may be progressively transformed into $\mathrm{D_{2}H^{+}}$ and $\mathrm{D_{3}^{+}}$, reducing the abundance of $\mathrm{H}_{2}\mathrm{D}^{+}$, and, at the same time, boosting the production of $\mathrm{N}_{2}\mathrm{D}^{+}$. This behaviour could be confirmed with observations of p-$\mathrm{D_{2}H^{+}}$.

A weak $\mathrm{N}_{2}\mathrm{D}^{+}$(3–2) line in Clump 7 may alternatively be caused by excitation effects, because it has a critical density an order of magnitude larger than the o-$\mathrm{H}_{2}\mathrm{D}^{+}$ line (see also Pillai et al., 2012) and larger than the mean $n(\mathrm{H_{2}})$ of the source. A similar effect could result from a large part of the gas being affected by complete depletion of heavy elements. In both cases the potential of the relative abundance of o-$\mathrm{H}_{2}\mathrm{D}^{+}$ and $\mathrm{N}_{2}\mathrm{D}^{+}$ in tracing evolution remains unaltered, because density increases with time, and the region of complete depletion is progressively eroded by the warm-up caused by the YSOs.

Without high-resolution observations, however, we cannot exclude that $\mathrm{H}_{2}\mathrm{D}^{+}$ is also abundant in compact $\mathrm{N}_{2}\mathrm{D}^{+}$-emitting regions for sources like Clump 2. Here, o-\mathrm{H}_{2}\mathrm{D}^{+}$$(1_{1,0}-1_{1,1}) could be completely saturated, and dilution in the APEX beam could cause the non-detection of the line.

Acknowledgements.

SB is financially supported by CONICYT Fondecyt Iniciación (project code 11170268), CONICYT programa de Astronomia Fondo Quimal 2017 QUIMAL170001, and BASAL Centro de Astrofisica y Tecnologias Afines (CATA) AFB-17002. DRGS is financially supported by CONICYT Fondecyt regular (project code 1161247) and the international collaboration project PII20150171. AG thanks the Astronomy Department, University of Concepción for having supported his visit during July 2018 (via Redes Internacionales project number REDI170093). B.K. acknowledges funding from the German Science Foundation (DFG) within the Priority Programm ”The Physics of the ISM” (SPP 1573) via the grant BA 3706/3-2. This work was partly supported by the Collaborative Research Council 956, sub-project A6, funded by the Deutsche Forschungsgemeinschaft (DFG). This paper is based on data acquired with the Atacama Pathfinder EXperiment (APEX). APEX is a collaboration between the Max Planck Institute for Radioastronomy, the European Southern Observatory, and the Onsala Space Observatory. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013, http://www.astropy.org), of NASA’s Astrophysics Data System, and of Matplotlib (Hunter, 2007). MCWeeds makes use of the PyMC package (Patil et al., 2010).

Appendix A Additional spectra