© ESO, 2015

1. Introduction

The chemical evolution in high-mass star formation is still poorly understood and a field of intense investigations. The question of which molecules are to be used to distinguish between different evolutionary stages is of great interest. These so-called chemical clocks could be used to derive lifetimes of the different stages and help to infer the absolute ages of these objects. In addition, studying deuterium chemistry is also very useful to constrain physical parameters, for example the ionization fraction, temperature, and density (e.g., Crapsi et al. 2005; Chen et al. 2011). In particular, deuterated molecules are very prominent candidates for probing this evolutionary sequence, since their chemistry highly depends on the temperature and the thermal history of an object (Caselli & Ceccarelli 2012; Albertsson et al. 2013). The deuteration fraction (the ratio between the column density of a deuterated molecule and its non-deuterated counterpart) is therefore an important parameter in order to study these evolutionary effects.

Theoretical and observational deuteration studies of low-mass star-forming regions revealed a large increase by several orders of magnitude in the deuteration fraction of starless cores with respect to the cosmic atomic D/H ratio of ~10^-5 (Linsky 2003; Oliveira et al. 2003) and discussed possible trends with the evolutionary state of the star-forming region (e.g., Caselli et al. 2002; Crapsi et al. 2005; Bourke et al. 2012; Friesen et al. 2013). Correlations of the deuteration fraction were seen, for instance with the dust temperature and the level of CO depletion (Emprechtinger et al. 2009) or with the density (Daniel et al. 2007). Whether the deuterium chemistry during high-mass star formation behaves similar to that of low-mass cores is so far poorly constrained by observations. The current studies mostly target single-deuterated species instead of a larger set of molecules, or are focused on a limited number of sources. Miettinen et al. (2011) found deuteration fractions in a sample of seven massive clumps associated with IRDCs that are lower than the values found in low-mass starless cores. Early studies of very young IRDCs by Pillai et al. (2007) and more evolved high-mass protostellar objects (HMPOs) by Fontani et al. (2006) indicated a trend of higher deuteration fractions for the younger, cooler sources. A recent attempt to systematically study a larger sample of different evolutionary stages by Fontani et al. (2011) revealed that the N₂D⁺/N₂H⁺ column density ratio is an indicator for the evolutionary stage in high-mass star formation. Chen et al. (2011) observed several dense cores covering different evolutionary stages in three massive star-forming clouds and studied the deuteration fraction of N₂H⁺ and the role of CO depletion in this context. They found a clear trend of decreasing deuteration fraction with increasing gas temperature tracing different evolutionary stages. They also found an increasing trend of the deuteration fraction with the CO depletion factor, which is similarly seen in low-mass protostellar cores.

To study the deuteration in high-mass star-forming regions in an evolutionary sense, we divided the high-mass star formation sequence into different stages (see also Gerner et al. 2014). Beuther et al. (2007a), Zinnecker & Yorke (2007), and Tan et al. (2014) divided the evolutionary sequence into different phases based on the physical conditions. We describe the evolutionary picture from an observational point of view and distinguish between four observationally motivated stages based on the underlying physical sequence. The first is an initially starless infrared dark cloud phase (IRDC). At this point, these objects are close to isothermal and consist of cold and dense gas and dust. In this approach we do not consider a long-living pre-IRDC phase, which is proposed in theoretical works (e.g., Narayanan et al. 2008; Heitsch et al. 2008) and also supported by observations (e.g., Barnes et al. 2011). This phase probably is less dense and we define the year zero of our evolutionary sequence in our model when the densities start to be higher than 10⁴ cm^-3. While starless IRDCs only emit in the (sub-)millimeter regime, places of beginning star-formation become detectable as point sources at μm-wavelengths. Eventually, the overdensities within the IRDC begin to collapse and form one or several accreting protostars with >8 M_⊙ in the next phase, i.e., a HMPO. The internal sources of HMPOs emit at mid-infrared wavelengths, and their radiation starts to heat up the environment, leading to non-isothermal temperature profiles. The higher temperatures boost the molecular complexity leading to the hot molecular core (HMC) phase. From a physical point of view, this phase is a sub-group of the HMPO phase but is clearly distinguished from a chemical point of view, because it is driven by the higher temperatures that liberate molecules from molecular-rich ices and increase the molecular complexity of the source. Finally, the UV-radiation of the central star(s) ionizes the ambient gas and creates an ultra-compact Hii region (UCHii). This is the last stage we considered in the evolutionary sequence. These objects presumably have stopped accreting, and complex molecules seen in the HMC phase are not longer detectable. It is possible that overlaps occur between these stages, leading to HMCs associated with UCHii regions and even still-accreting protostars within UCHii regions. High-mass star-forming sites are complex objects. To circumvent the problem of the coexistence of different stages in one object, we wish to statistically characterize the evolution during the different stages.

Here, we continue and extend an investigation of the chemical evolution in 59 high-mass star-forming regions in different evolutionary stages (Gerner et al. 2014) toward deuterated molecules. In the previous work we measured the beam-averaged column densities of 14 different molecular species and derived a chemical evolutionary picture across the evolutionary sequence in high-mass star formation starting with IRDCs via HMPOs to HMCs and finally UCHii regions. We found that overall, the chemical complexity, column densities, and abundances increase with evolutionary stage. We fitted the data with a 1D physico-chemical modeling approach and found good agreement with the observations. Here we want to measure the deuteration fractions of the four deuterated molecules DCN, DNC, DCO⁺, and N₂D⁺ and test their correlations with evolutionary stage and with physical parameters such as the luminosity of the objects. Furthermore, we model the derived column densities with an advanced 1D deuterium chemical model along the evolutionary sequence from IRDCs via HMPOs to HMCs and UCHii regions.

The structure of the paper is the following. In Sect. 2 we introduce the source sample, followed by a description of the observations in Sect. 3 and an introduction to deuterium fractionation in Sect. 4. In Sect. 5 we present the results of the analyzed observational data. In Sect. 6 we introduce the model used to fit the data and discuss the modeling results as well as their implications. We conclude with a summary in Sect. 7.

2. Source sample

The sources were taken from Gerner et al. (2014) and were initially selected from different source lists. The total sample contains 59 high-mass star-forming regions, consisting of 19 IRDCs and 20 HMPOs as well as 11 HMCs and 9 UCHiis. The sources were selected from well-known source catalogs of the literature without specific selection criteria such as spherical symmetry. The lists of the IRDCs were first presented in Carey et al. (2000) and Sridharan et al. (2005) and are part of the Herschel guaranteed time key project “The Early Phase of Star Formation” (EPOS, Ragan et al. 2012). This sample consists of 6 IRDCs showing no internal point sources below 70 μm and 13 IRDCs that have internal point sources at 24 μm and 70 μm. The HMPOs were taken from the well-studied sample by Sridharan et al. (2002) and Beuther et al. (2002b,c). HMC sources were selected from the line-rich sample of Hatchell et al. (1998), including a few additional well-known HMCs: W3IRS5, W3H₂O, and Orion-KL. For the UCHii regions, we selected line-poor high-mass star-forming regions from Hatchell et al. (1998), and additional sources from Wood & Churchwell (1989b). Recent studies toward single sources of the ones included in this work overall confirm the evolutionary classification given in the referenced papers. The full source list is given in Table A.1. The table also lists the distance of each source. When the distance ambiguity of near and far kinematic solutions could not be resolved, we used the near kinematic solution. In Fig. 1 we show the distance distributions for all four subsamples including their median and mean value. The median and mean distance for the IRDCs, HMPOs, and HMCs are very similar with values around 4 kpc. Thus, the wider range of distances of the total sample does probably not affect the comparison between the median or mean properties of the different subsamples, for example, their detection statistics or their derived median column densities (e.g., due to uncertain beam-filling fractions). Only the UCHii subsample shows significantly larger distances. However, since this subsample is not the main interest in this work, this does not dramatically influence the results of this study.

3. Observations

The 59 sources were observed with the Arizona Radio Observatory Submillimeter Telescope (SMT) in 2013 between February 12−15, March 10−13 and on March 31 and April 1 with ~100 h total observing time. For the observations we used the ALMA type 1.3 mm dual polarization sideband-separating heterodyne receiver and the filterbanks as backends with a resolution of 250 kHz, which corresponds to ~0.3 km s^-1 resolution in velocity. The beam size of the SMT at 1.3 mm is ~30″. One single integration took 5 min in position-switching mode, with 2.5 min on-source time. The emission in the deuterated lines was integrated two to three times, the emission in the non-deuterated lines one to two times, depending on the observing conditions. The data was calibrated using data of Jupiter from the same observation runs assuming a sideband rejection of 13 db. The mean system temperatures of the spectra T_sys is 380 K. The data reduction was conducted with the standard GILDAS¹ software package CLASS. All spectra from each source were baseline subtracted, calibrated to T_mb scale with typical beam efficiencies of 0.6, and averaged. In a few cases we used the 1 MHz filterbanks spectra because of the very broad lines (e.g., in Orion-KL) or when the line was not detected in the 250 kHz filterbanks (with 128 channels), but in the 1 MHz filterbanks (with 512 channels). The line integrals were measured by summing the line emission channel by channel for detected lines. The detection criteria was a signal-to-noise ratio S/N> 3.

The line parameters of the observed molecules are given in Table 1.

3.1. Rms

The median 1σ root mean square (rms) values for the deuterated molecules are 0.031 K for DCN (3−2), 0.024 K for DNC (3−2), 0.025 K for DCO⁺ (3−2), and 0.035 K for N₂D⁺ (3−2). DCN (3−2) was observed in the later shifts, and the higher rms value is due to poorer observing conditions during that period. The higher rms value for N₂D⁺(3−2) is due to problems at this specific transition frequency and is described in the following section. An overview of all rms values is given in Fig. 2. The strong outlier in DCN (3−2) is from Orion-KL, which is still a detection.

3.2. Problems with N₂D⁺ spectra

Reducing and analyzing the spectra of the N₂D⁺ (3−2) transition at 231.32 GHz was problematic because of blending by an overlaying pressure broadened ozone line in the atmosphere. This ozone line depends on the elevation of the source and time of the observation and affects the different sources with varying strength. As a result of this contamination, the mean 1σ rms value of the N₂D⁺ line is 40% higher than for the other deuterated spectral lines, and thus the threshold for a detection is higher. In the five sources IRDC18454.1, IRDC18454.3, Orion-KL, HMC010.47, and HMC031.41 the ozone line prevents distinguishing between a detection or a non-detection.

4. Deuterium fractionation

Under cold conditions, the deuteration fraction of molecules can be enhanced through reactions with, for instance deuterated H, H₂D⁺. While the chemical reaction that deuterates H proceeds without a thermal barrier, the backward reaction is endothermic and has a thermal energy barrier (Watson 1974; Caselli & Ceccarelli 2012; Albertsson et al. 2013). This energy barrier leads to an enhanced formation fraction at low temperatures, but becomes negligible for higher temperatures at which the backward reaction becomes as efficient as the forward reaction. The chemical reactions and their energy barriers that introduce deuterium into some of the key molecules are the following:

The temperature range for an effective D-enhancement for the pathway via H isotopologs is ~10−30 K, whereas it is ~10−80 K for pathways via light hydrocarbons (Millar et al. 1989; Albertsson et al. 2013). This enhances the deuteration fraction of the key molecules H and CH in different temperature regimes. Another effect that increases the H₂D⁺/H ratio is the depletion of neutral gas-phase species (e.g., CO, N₂; see Dalgarno & Lepp 1984; Roberts & Millar 2000). This enhancement is then imprinted on the molecules formed through these reaction partners.

In the literature different models exist that describe the formation routes of deuterated molecules and their relative importance. The dominant formation pathways of DCO⁺ and N₂D⁺ are via a low-temperature route through H isotopologs, whereas DCN can be formed via a fractionation route involving light hydrocarbons. The back-reaction of the route via H isotopologs sets in at temperatures around 30 K and their deuteration fractions decrease with rising temperature. According to Roueff et al. (2007), both DCN and DNC can be formed efficiently at low temperatures via deuteration of HNC or HCN. However, at higher temperatures DNC is destroyed by reactions with atomic oxygen, unlike DCN. The chemical models of Turner (2001) indicate that only DCN is also formed via reactions with light hydrocarbons involved, but not DNC. If the Turner (2001) scheme is correct, then molecules such as C₂D, HDCO or C₃HD should show a similar behavior as DCN with temperature, since they are also formed via CH₂D⁺ and C₂HD⁺.

Parise et al. (2009) found a low DCO⁺/HCO⁺ column density ratio but significant deuteration fractions for HCN and H₂CO under temperature conditions of ~70 K toward the Orion Bar. Model calculations by Roueff et al. (2007) predicted that in general the DNC/HNC, DCO⁺/HCO⁺ and N₂D⁺/N₂H⁺ column density ratios decrease with temperature, but are almost constant with density. The DCN/HCN column density ratio shows a more complex behavior with temperature, reaching the highest ratio for ~30 K and shows a stronger increase with density. They found that the reason for this behavior is twofold. First, it is due to the enhanced abundance of radicals (e.g., CHD and CD₂) that form DCN. Second, the main destruction pathways of DCN are reactions with the ions HCO⁺ and H₃O⁺, which leads to DCNH⁺ which subsequently returns to DCN via dissociative recombination. The deuteration fractions also strongly depend on the assumed elemental abundances.

5. Results

5.1. Detection fractions for DCN, DNC, DCO⁺, and N₂D⁺

The detection fractions of DCN, DNC, DCO⁺, and N₂D⁺ in the respective stages are shown in Fig. 3. In this figure, the IRDC stage is split into the two source categories mentioned in Sect. 2, without and with an embedded point source. Based on this figure, we identify three trends:

1) N₂D⁺(3−2) is only detected toward IRDCs and HMPOs. Since the N₂H⁺ abundance is almost constant over the full evolutionary sequence (Paper I), a strong temperature dependence of the N₂D⁺ production leading to its disappearance at T ≳ 30 K is most likely the answer to its apparent absence in the HMC and UCHii stages. Another complication are the problems with the particular bandpass of its (3−2) transition mentioned in Sect. 3.2, which lead to a higher mean 1σ rms value. This does not mean that N₂D⁺ is absent, since some of our sources with non-detections in N₂D⁺(3−2) were detected by Fontani et al. (2011) in N₂D⁺(2−1). This could indicate the importance of excitation conditions for the detection of the (3–2) transitions. Nonetheless, the detection rates of N₂D⁺(3−2) are clearly lower than the detection rates of the other observed molecular transitions of this work. 2) The only deuterated molecule detected in the presumably coldest IRDCs without an embedded point source at 24 μm (Spitzer/MIPSGAL) or 70 μm (Herschel/PACS) is DCO⁺, which also has the highest detection fraction in the more evolved IRDCs. 3) The detection fractions of DCN, DNC, and DCO⁺ in warm IRDCs up to UCHii regions are comparable and peak at the HMC stage. In general, the detection fractions toward the observed high-mass star-forming regions are high, ≳60%. It is important to mention here that the strength of a transition is the product of column density and temperature, and thus detections at the lower temperatures present in the earlier stages are more difficult. A lower detection fraction does not necessarily mean that the column densities in the earlier stages are lower than at the later stages. The differences in column densities during the evolutionary sequence are discussed in Sect. 6.1.

5.2. Excitation temperatures and final column densities

In this section we discuss the excitation temperatures used to derive the column densities and how we combined the data from Paper I and this work to obtain the final column densities. The column densities are derived in Appendix A.1.

The column densities for the non-deuterated species were partly taken from Gerner et al. (2014) and partly calculated from the new observations. A detailed overview of the exact combination of the previous work and this one is given below. In the previous study we observed the (1−0)-transitions of several molecules with the IRAM 30 m telescope. The size of the beam for these observed (1−0)-transitions is very close to the size of the beam of the (3−2)-transitions observed for this work with the SMT. In the previous work we used likely optically thin (1−0)-transitions of H¹³CO⁺ and HN¹³C and, because of their hyperfine structure, optical-depth-corrected N₂H⁺ and HCN to derive their column densities. The spectra of the first study were partly affected by source emission at the present off-positions and strong optical depth effects. In these cases the optical depth and the integrated intensity could not be reliably determined. Thus we refrained from using these lines in the analysis. In this work we have additional information from the (3−2) transitions to complement the missing column densities. To obtain consistent results from both observation runs we compared the column densities of H¹³CO⁺(1−0) with the column densities of H¹³CO⁺(3−2) derived with T_ex = 20.9 K (IRDC), T_ex = 29.5 K (HMPO), T_ex = 40.2 K (HMC) and T_ex = 36.0 K (UCH ii). These are the mean temperatures of the best-fit models from Gerner et al. (2014). The comparison is shown in Fig. 4.

The derived column densities for the IRAM data are all clearly higher. Excluding problems in the quality of the data and the calibration leads to the assumption that, while the (1−0) transitions are in local thermodynamical equilibrium (LTE), the upper levels are subthermally populated and thus the (3−2) transitions are not in local thermodynamical equilibrium (LTE). This implies that the average densities of the gas are below ~10⁶−10⁷ cm^-3. Since the two different lines trace gas with different excitation conditions, the spatial extent probed by the two transitions might be different. However, we do not know the beam-filling factors for the different molecules and their different transitions. In addition, we calculated the column densities from the H¹³CO⁺(1−0) and (3−2) transitions with the non-LTE radiative transfer code RADEX² (van der Tak et al. 2007). The derived column densities for H¹³CO⁺ agreed within a factor of 2 with the column densities derived under the assumption of LTE from the H¹³CO⁺(1−0) transition. This is well within the uncertainties and validates the assumptions made for the (1−0) transitions. To compensate for the differences in the excitation conditions between (1−0) and (3−2) transitions, we calculated the excitation temperature of the (3−2) transition for each source, which would be necessary to derive the total column densities derived from the (1−0) transition. From those excitation temperatures we computed the median value for each stage. This resulted in T_ex = 5.2 (IRDC), T_ex = 6.2 (HMPO), T_ex = 7.2 (HMC) and T_ex = 6.4 (UCH ii) for the (3−2) transitions. The comparison of the calibrated H¹³CO⁺(3−2) column densities derived with the lower excitation temperatures and the (1−0) data is shown in Fig. 4. This correction should also reduce the error in the beam-filling fraction and make the different transitions comparable.

For the four different molecules we derived the column densities as described below.

5.2.1. HCO⁺ and HNC

For HCO⁺ and HNC we used H¹³CO⁺(1−0) and HN¹³C transitions, assuming the relative isotopic ratio depending on the galactocentric distance described in Wilson & Rood (1994). In Gerner et al. (2014) we used the representative value for the Sun and the local ISM of ¹²C/¹³C = 89 (Lodders 2003). This decreases the derived column densities for this work.

5.2.2. N₂H⁺ and HCN

For N₂H⁺ we took the column densities from the (1−0) transitions from Gerner et al. (2014), for which we used, in case of optically thick lines, the optical depth and the excitation temperature from the hyperfine fits made with CLASS³. The hyperfine structure fitting routine (METHOD HFS) assumes that all components of the hyperfine structure have the same excitation temperature and the same width and that the components are separated in frequency by the laboratory values. From the comparison of the ratios of the line intensities of the hyperfine components to the theoretically expected ratios, the fit routine estimates the optical depth of the line. In the optically thick case, this allows us to determine the excitation temperature. In the optically thin case we assumed the excitation temperatures given in Sect. 5.2.

For HCN we used the optical-depth-corrected column densities from the (1−0) transitions when available. In all other cases we computed the mean difference between column densities derived from the (1−0) transitions and the (3−2) transitions and applied this factor of ~1.7 to the derived column densities from the (3−2) transitions.

All column densities are beam-averaged quantities. The resulting median abundances including all detections and upper limits for each subsequent evolutionary phase are given in Table 2. In deriving the column densities, we assumed the excitation temperatures of the deuterated and non-deuterated molecule of the same transition to be equal. This assumption is valid for high densities and high excitation temperatures, but might lead to an underestimation of the deuteration column density ratios D_frac by a factor of ~2−4 (Shirley, in prep.).

6. Discussion

6.1. Deuteration fractions

The spread of the deuteration fractions D_frac and the median value for each stage is shown in Fig. 5. These values include detections and upper limits. For N₂D⁺ the median 3σ-limits of non-detections for the HMC-stage and UCHii-stage are given because the detection fraction is zero. These values can be considered as a sensitivity limit.

The ratios of the column densities of DNC, DCO⁺, and N₂D⁺ with their non-deuterated counterparts HNC, HCO⁺, and N₂H⁺ all show decreasing trends with evolutionary stage, despite the wide spread of ratios within individual stages. Only HCN shows an increase in the ratio with the maximum reached in the HMC-stage. In some cases the consecutive stages appear not very different from each other within the error budget (e.g., HMPO-HMC of DCO⁺/HCO⁺, HMC-UCHii of DNC/HNC). The reason for this is probably that the evolution throughout the sequence is continuous, and there might be overlaps between consecutive stages, for example, some of the HMCs are associated with UCHii regions. One can distinguish between two different scenarios that can lead to an overlap between consecutive stages. For HMC029.96 (G29.96), for instance high spatial resolution observations show a separation of the actual UCHii region and the neighboring HMC of 2″ (Cesaroni et al. 1994; Beuther et al. 2007b). The HMC and the UCHii region are two (or more) distinct entities. The second scenario are cases like HMC009.62 (G9.62), where a separate UCHii region exists as well, but in addition, a small hypercompact Hii region is present at the location of the HMC itself (Testi et al. 2000). This is a general caveat in single-dish observations of HMCs that might be contaminated by UCHii regions in the beam. Here, the sources were classified based on their chemistry which is probably dominated by the HMC. However, this does not totally cancel this effect out, and differences in the measured ratios between different stages are weaker in some cases. Nevertheless, the global trends are clearly visible.

Furthermore, Fig. 6 shows the observed column density ratios of DCN/DNC and HCN/HNC for sources detected in both molecules. While the ratio of the non-deuterated molecules is almost constant, the ratio of the deuterated molecules clearly shows a peak at the HMC stage. reflected by the results of the one-way ANOVA tests that yield probability values p< 10^-4 for DCN/DNC and p = 0.17 for HCN/HNC. The KS2 test of the consecutive stages confirms these findings. We conclude that DCN can be formed more efficiently or is destroyed less efficiently than DNC in the more evolved and thus warmer sources. The observed behavior can be understood according to the dominant formation pathways scenario stated by Turner (2001). He claimed that all four deuterated molecules can be efficiently formed at low temperatures via H₂D⁺, but only DCN can be formed at higher temperature via CH₂D⁺. This would lead to the observed trends in the deuteration fractions, since DCN can still be formed in the more evolved stages, in contrast to the other observed deuterated molecules. However, newer chemical models (e.g., Roueff et al. 2007) or the updated model we use, do not support this scenario of a difference in dominant formation pathways between DCN and DNC. According to these new models, both isomers are formed from light hydrocarbons and thus exhibit a similar evolution in the abundance ratios with temperature. Hiraoka et al. (2006) studied the association reactions of CN with D in laboratory experiments and found an intensity ratio of DNC/DCN of ~3 at a temperature of 10 K. At temperatures >20 K, the formation of DNC and DCN became negligible. This is consistent with the derived value for the IRDC stage in this work on the order of unity, but contradicts the higher DCN/DNC ratios found in the HMC stage, suggesting another possible formation route. The HCN/HNC median ratios is between 0.4−1.4 in the different stages and slightly increases with stage. This is consistent within the observational uncertainties with theoretical expectations of the ratio of 1 in cold clouds and >1 in warmer clouds (Sarrasin et al. 2010).

6.2. Relation between deuteration and other parameters

We measured the full width at half maximum (FWHM) of HCN, HN¹³C, H¹³CO⁺ and N₂H⁺. The isotopologs were assumed to be optically thin. For N₂H⁺ and HCN we corrected the measured optically thick FWHM Δν_thick according to Phillips et al. (1979) using the relation

(4)where τ is the optical depth of the line center and FWHM_thin the optically thin FWHM. We tested for linear relationships using the Pearson correlation coefficient. By definition, the absolute value of ρ is ≤1 with a stronger correlation for larger ρ. ρ = 0 means no correlation. The sign indicates a positive or negative correlation between the two quantities. To estimate ρ we randomly drew sets of data points from the measured values with their corresponding errors and the upper limits, and found for each of these drawn data sets the underlying correlation coefficient. We then determined the most frequent correlation coefficient among all data sets and its 16% and 84% confidence values (Bailer-Jones et al., in prep.). The correlation plots are shown in Fig. 7, and the corresponding Pearson correlation coefficients ρ are given in Table 3. The plots show that earlier phases have a smaller FWHM.

The resulting coefficient is consistent with the data being uncorrelated for all four molecules. The FWHM of N₂H⁺ might be an exception and correlated, but is still <0.5 and no conclusive answer can be given. Although the correlation with the FWHM is weak, the decrease of the deuterium ratio with increasing FWHM is consistent with the picture of more quiescent, early stages with lower FWHM and higher deuteration fractions followed by more turbulent stages with higher FWHM and lower deuterium ratios.

Furthermore, we studied the correlation of the deuteration with the H₂ column density and luminosity L. The luminosities and H₂ column densities are shown in Tables A.2–A.3 and the corresponding correlation plots in Figs. 8 and 9. In general, the deuteration fraction shows similarly weak correlations with the luminosities as with the FWHM. Clearly, there is a lack of correlation with H₂ column densities. Of the four deuteration fractions, the strongest correlation is found for N₂D⁺, but the correlation is still very low. The luminosity of a source is not necessarily a tracer of its evolutionary stage, but Fig. 8 shows that the four different evolutionary stages are approximately separated with respect to the luminosity. Taking into account that the temperature of the objects increases with evolutionary stage, the luminosity can be used as a proxy for the temperature in the case of our source sample. Since higher luminosities indicate higher temperatures and smaller regions with T ≲ 20 K, where CO is frozen out, it reduces the overall abundance as well as the D/H ratio of N₂H⁺. Of all four molecular deuteration fractions, DCN and DCO⁺ show the weakest correlation with any of the shown parameters. However, none of the studied correlations have correlation coefficients ρ> 0.5.

6.3. Modeling the chemical evolution

With the observed column densities of 18 species (14 from Paper I) including 4 deuterated species in different stages of massive star formation at hand, we applied the iterative physico-chemical fitting model MUlti Stage CLoud codE with D-chemistry (MUSCLE-D) to these data. A list of all fitted species is given in Table 4. The model fits the evolution of the observed chemical data and thereby constrains basic physical properties in the assumed evolutionary path of high-mass star formation, such as mean temperatures and mean chemical ages. The chemical ages can be interpreted as typical lifetimes of the various stages. This model is the extended version of the iterative fitting model MUlti Stage CLoud codE (MUSCLE) already used and described in Gerner et al. (2014) and includes a physical and a chemical model. Here we briefly summarize the main characteristics of the two parts of MUSCLE-D, the physical model in Sect. 6.3.1 and the chemical model ALCHEMIC in Sect. 6.3.2.

6.3.1. Physical model

The physical model treats the star-forming region as spherically symmetric 1D-clouds with a fixed outer radius of r_out = 0.5 pc, based on the largest beam size of our observations of 30″ and the typical size of high-mass star-forming regions. High angular resolution studies revealed in some of the sources complex structures that deviate from the spherically symmetric case (e.g., Beuther et al. 2002a). However, in single-dish maps these sources appear to be symmetric (e.g., Beuther et al. 2002b). Since we analyzed single-dish data, a 1D physical model appears to be sufficient. More sophisticated models that take into account complicated substructures would increase the number of fitting parameters and might lead to an overinterpretation of the data. The radial density and temperature structure is modeled with

(5)and (6)respectively. The temperature and density profiles do not change with time. The fitted physical quantities are the inner radius r_in, the temperature and density at the inner radius, T_in and ρ_in and the power-law index of the density p. The inner radius is limited between 5 × 10^-5−5 × 10^-2 pc. While we modeled the IRDC as an isothermal sphere, the temperature structure of the more evolved stages is modeled by a inner flat plateau with T_in and a power law with slope q = 0.4 as a standard value for r>r_in (see van der Tak et al. 2000). The radial density profile within r_in is flat with ρ_in and decreases for r>r_in as a power law with slope p. The value of p was limited to values between 1.5 and 2.0 to save computing time. This range is supported by several observations, for example, Guertler et al. (1991), Beuther et al. (2002b), Mueller et al. (2002), and Hatchell & van der Tak (2003). The temperature and density profiles were fitted simultaneously. The model does not take into account radiative transport. In the model, the whole cloud is embedded in a larger diffuse low-density medium that shields the high-mass star-forming cloud from the interstellar FUV radiation.

6.3.2. Chemical model

The chemical model is an updated version of the time-dependent gas-grain chemical model ALCHEMIC described in Semenov et al. (2010) that we used and described in Gerner et al. (2014). In addition to this, the deuterium network from Albertsson et al. (2013) was added and extended with high-temperature reactions (Harada et al. 2010, 2012; Albertsson et al. 2014b) and ortho/para states of H₂, H and H and their isotopologs (Albertsson et al. 2014a). In total, the chemical network comprises 15 elements that can form 1260 different species from 38 500 reactions.

The initial abundances prior to the IRDC stage are taken from the low-metals set given in Lee et al. (1998) with a changed elemental abundance of Si (3 × 10^-9 with respect to H) and S (8 × 10^-7 with respect to H). The changes were needed to achieve proper fits to the IRDC phase. Initially, all metals (C, O, N, S, etc.) are in atomic form. Only H₂ is already in molecular form. The initial ortho-para ratio of H₂ is assumed to be the statistical value of 3:1. The model includes nuclear spin-state exchange reactions to account for the evolution of the ortho-para ratio and freeze out of CO. For the subsequent stages, the chemical outcome of the previous best-fit model is used as input initial abundances.

6.3.3. Fitting procedure

The different stages were fit iteratively using the physical and chemical model described above. We modeled the observed column densities for the IRDC, HMPO, and HMC stage by varying the parameters r_in, T_in, ρ_in and p. While keeping all other parameters fixed, we varied these four parameters and ran the model over 10⁵ yr for each of the realizations. Then we computed the χ²-value for each time step and model realization, given the observed mean column densities and computed model column densities. We assumed the standard deviation between modeled and observed values to be one order of magnitude as a typical value. Molecules that were detected in less than 50% of the sources within one stage were considered as upper limits. Finally, the model with the minimum χ²-value was found as the best-fit model that matches with the calculated mean column densities for each stage best.

As a result of the iterative fitting along the evolutionary sequence, the introduced uncertainties increase with evolutionary stage, leading to a limited confidence in the obtained results for the later stages.

6.4. Modeling results

We fitted the combined data from Gerner et al. (2014) and this work with this model assuming an ortho-para H₂ ratio o/p = 3:1. The resulting best-fit model parameters for the IRDC, HMPO, HMC, and UCHii stage are shown in Tables A.4−A.7. The evolution of the best-fits with time is shown in Fig. 10. These distributions give a good impression of the uncertainty on the chemical age. The best-fit lifetime is 16 500 yr for the IRDC stage. A good fit was also achieved with a chemical age of only 1000 yr. However, this seems to be a rather short time and might be interpreted as a lower limit. The HMPO best-fit age yielded 32 000 yr with probable values between 10 000−40 000 yr. For the HMC stage we found a best-fit age of 35 000 yr; values below 15 000 yr are improbable. The best-fit age for the UCHii stage was found to be rather young with 3000 yr, but quite unconstrained with probable values up to several 10 000 yr.

The observed and modeled column densities are shown in Tables A.8−A.11. The best-fit of the IRDC stage reproduces 18 of 18 molecules within the assumed combined (observational + chemical) uncertainty of one order of magnitude. In the HMPO stage the model is able to reproduce 16 of 18 molecules. SO is overproduced by a factor of ~20 and CH₃OH underproduced by a factor of ~40. The underproduction of methanol is possibly due to shock- or outflow-triggered enhanced desorption of methanol ice from grains, which is not taken into account in the model. The reason for the misfit of SO might be the poorly understood chemistry of sulfur-bearing species in modern astrochemical models in general. In the HMC stage the model could fit 14 of 18 species. In addition to C₂H, which is ~20 times underproduced, HNCO, SO, and CH₃OH are misfit by more than a factor of 100. The underproduction of CH₃OH and overproduction of SO in the HMPO stage continues in the HMC phase. The overproduction of HNCO might be connected to not well enough understood shock- and surface chemistry. This difference in C₂H between the model and observations might be influenced by UV-radiation of the central star(s) or a clumpy structure of the environment, which is not considered in the model. This is especially important for the UCHii regions, but it is also already present in some HMCs. For the last considered stage of an UCHii region, the model reproduces 13 of 18 species. As in the HMC stage, the molecules HNCO, SO, and C₂H are misfitted in the UCHii stage. In addition, SiO and DCO⁺ are slightly overproduced.

In general, the overall fit of all 18 molecules of the four fitted phases is good. The specific time dependent evolution of the best-fit abundances of molecules DCN, DNC, DCO⁺, N₂D⁺ and their non-deuterated counterparts are shown in Fig. 11. Between two consecutive stages, the physical parameters are instantly changed and the molecular species show a quick response to that change, followed by a slower evolution under the new constant conditions. Figures A.1–A.4 show the modeled column densities in each stage separately. In total, the best-fit ages add up to 85 000 yr, which is on the same order as typical models of high-mass star formation of 10⁵ yr (e.g., McKee & Tan 2003; Tan et al. 2014).

The evolution for the different molecular ratios are shown in Fig. 12. The modeled abundance ratios partly show larger deviations from the observed values because the uncertainties from both molecules add up. In Sect. 6.1 we discussed the ratio of DCN/DNC and a possible chemical explanation based on a chemical network by Turner (2001). However, the chemical network used in our work shows a different trend with a constant ratio rather than an increase of DCN/DNC with temperature. Thus the observed trend disagrees with our model predictions, and the possibility of different chemical formation pathways needs to be tested by future observations of more deuterated molecules sharing the same formation pathways as DCN (e.g., light hydrocarbons with two or three C-atoms). The model agrees with the constant ratio of HNC/HCN of unity with recent developments in the collisional rates by Sarrasin et al. (2010).

6.5. Comparison with best-fit models of Paper I

In Paper I we iteratively fitted our model to an observational dataset containing the median column densities of 14 different non-deuterated species. The column densities were derived with typical temperatures for each stage. The model followed the observed evolution and successfully fit most of the species with few exceptions in individual stages. The fits also constrained the physical structure of the model and yielded mean best-fit temperatures for each stage. The mean temperatures deviated from the temperatures used to derive the observed column densities for the later stages. In a second step we fitted the newly derived column densities. The temperature in the best-fit models decreased again and we could not find a converging solution.

In this work we used the best-fit temperatures of the first model from Paper I to recalculate the median column densities for the previously analyzed 14 species and an additional four deuterated species. The temperature in the best-fit model of the IRDC stage is lower by a factor of 2 compared to Paper I, but the best-fit temperatures of the IRDC stage for both Paper I and this work agree with typical IRDC temperatures. The mean density increased by about an order of magnitude. The best-fit lifetime in this work is a factor of 1.5 longer, but agrees within the assumed uncertainties. The HMPO best-fit mean temperature in this work is lower by 30%. The mean density increased by a factor 3−4. The best-fit lifetime is shorter by a factor of ~2. The HMC best-fit mean temperature in this work is higher by 25%. The mean densities are similar. In comparison to Paper I, the best-fit lifetime is about 20% lower. The UCHii best-fit mean temperatures are similar. The best-fit ages deviate by a factor of 4, but are consistent within their likely ranges. The mean density is lower by one order of magnitude.

The total lifetimes from the IRDC to the UCHii phase in Paper I add up to ~110 000 yr. This work yields a total lifetime that is slightly shorter with ~85 000 yr. The molecules that could not be reproduced in Paper I and this work are not the same. Possible reasons are the additional four deuterated molecules as free fit parameters that were treated with equal weight as well as the revised excitation temperatures used to derive column densities compared to Paper I (see Sect. 5.2). Although D-chemistry is seen as a tracer of the thermal history and conditions of an object and we see weak correlations with luminosity as an evolutionary tracer, the model fits did not improve substantially. Increasing the number of fitted molecules from 14 to 18 species did not improve the achieved results much, probably because the less statistical importance of the deuterated molecules within the total number of 18 species.

6.6. Comparison with literature

In low-mass star formation both observations and theory show that it is possible to use the deuteration fraction of a molecule as an evolutionary tracer (e.g., Crapsi et al. 2005; Aikawa et al. 2005; Emprechtinger et al. 2009). In a recent work, Fontani et al. (2011) observed 27 cores within high-mass star-forming regions and derived the deuteration fraction of N₂D⁺/N₂H⁺. They found a similar behavior as in the low-mass regime. The deuteration fraction is highest in the starless cores with D_frac = 0.26 and decreases during the formation of the protostellar objects to D_frac = 0.04. We find a similar trend in our work, but the median deuteration fraction we measure in our IRDC sample is lower by a factor of four than their high-mass starless core sample. This might be due to the fact that our IRDC sample contains starless as well as already more evolved objects inhabiting 24 μm sources. The same is true for the HMPO samples of Fontani et al. (2011) and this work. We find lower ratios by about a factor of four. The highest D/H ratio is that of N₂D⁺/N₂H⁺ in IRDC20081 of D_frac = 0.12. This source inhabits no point sources below 70 μm, but it is close to a nearby source with extended emission. It is presumably in a very early pre-stellar phase. Miettinen et al. (2011) studied seven massive clumps within IRDCs and derived deuteration fractions for N₂H⁺ between 0.002−0.028 and for HCO⁺ between 0.0002−0.014. While the range of values for HCO⁺ is comparable with the IRDCs in our sample, the deuteration fractions for N₂H⁺ are about one order of magnitude higher in this work. However, this agrees with the high number of upper limits due to non-detections of N₂D⁺ in this work. Chen et al. (2011) observed several cores in various evolutionary stages and found deuteration values of N₂H⁺ between 0.004−0.1. The source lists have three targets in common. For the HMPO18151 both works agree on a value of D_frac = 0.01. For the other two sources in common, IRDC18151 and IRDC18223, we find a slightly higher value of D_frac = 0.03 instead of D_frac = 0.02. While the highest value they found is similar to the highest value we observe, they detected N₂D⁺ in more sources and also found lower deuteration fractions, especially in the more evolved sources for which we have no detections of N₂D⁺. Crapsi et al. (2005) studyied the N₂D⁺/N₂H⁺ ratio in low-mass starless cores and revealed deuteration fractions on the percentage level up to almost 50% in the most extreme case of Oph D. The highest ratio observed in our survey is ~0.1. Sakai et al. (2012) measured deuteration fractions of DNC/HNC for IRDCs and HMPOs. They found D_frac = 0.003−0.03 for Midcourse Space Experiment (MSX) dark sources at 8 μm and D_frac = 0.005−0.01 for HMPOs. These ratios agree well with D_frac = 0.004−0.02 for our IRDC sample and D_frac = 0.0015−0.015 for our HMPO sample. In a recent study, Kong et al. (2013) explored the effect of different parameters on the deuteration fraction of N₂H⁺ in dense and cold environments. They found an increase of the deuteration fraction with decreasing temperature. However, they also found a positive correlation with the density, which is not clearly seen in our data. They employed a 0D model to study deuteration in star-forming clouds and found that the timescales to reach equilibrium in the abundances is on the order of several free-fall times (~10⁶ yr) for typical densities of high-mass cores. In contrast, the best-fit models in our work are not supposed to reach chemical equilibrium and thus our model predicts a total timescale of ~10⁵ yr for the best-fit models. This is probably because in a 1D model with a power-law density profile with higher densities in the inner region timescales in the chemical evolution are shorter than in a 0D model.

7. Conclusion

We extended the analysis from Gerner et al. (2014) of the chemical evolution in 59 high-mass star-forming regions for deuterated molecules. We measured beam averaged column densities and deuteration fractions of the four deuterated species DCO⁺, DCN, DNC, and N₂H⁺. We found an overall high detection fraction toward the high-mass star-forming regions, except for N₂D⁺, which is probably due to the limited sensitivity of our survey. The detection fraction of DCO⁺, DCN, and DNC increases from IRDCs to HMPOs and peaks at the HMC sample, where we detect all three molecules in all sources but one, and it drops again toward the UCHii stage.

The (3−2) lines are subthermally populated. The median deuteration fractions excluding upper limits are 0.02 for DCN, 0.005 for DNC, 0.0025 for DCO⁺, and 0.02 for N₂D⁺. The deuteration fractions of DNC, DCO⁺, and N₂D⁺ show a decreasing trend from IRDCs over HMPOs and HMCs to UCHii regions, which supports the hypothesis that deuterated molecules may be used as indicators of the evolutionary stage in high-mass star-forming regions.

In general, we find no correlation between the deuteration fraction of the various molecules and physical parameters for DCO⁺/HCO⁺ and DCN/DNC. Only N₂D⁺/N₂H⁺ shows a slight anticorrelation with the luminosity and the FWHM. The total measured column density of the gas does not correlate with the deuteration fraction. This result indicates that within the range of probed densities, deuteration depends more strongly on the temperature of the environment than on the column density and is enhanced in colder, less luminous regions. However, Albertsson et al. (2013), for example, found in their models differences in the deuteration with density at lower densities.

The slight anticorrelation with luminosity and FWHM, which is related to the evolutionary stage, indicates that evolutionary stage plays an important role for the deuteration fraction. But its huge scatter within single stages leads to the assumption that the evolution might not be the only factor.

Furthermore, we fit the observed data with a chemical model and find reasonable physical model fits. Combining observations of non-deuterated and deuterated species to obtain best-fits led to reasonable chemical and physical models. The best-fit models produce reasonably good results for all stages. As a result of the large uncertainties in the observations and the model and the wide spread of observed values within the four subsamples, we could not substantially improve the best-fit model results compared to Paper I. Another reason is that the combination of the four deuterated molecules with the 14 non-deuterated species probably reduced the statistical importance of each molecule and thus the effect on the best-fit of the model.

---

Acknowledgments

The authors thank the anonymous referee for helping to improve the paper. T.G. is supported by the Sonderforschungsbereich SFB 881 “The Milky Way System” (subproject B3) of the German Research Foundation (DFG). T.G. is member of the IMPRS for Astronomy & Cosmic Physics at the University of Heidelberg. This research made use of NASA’s Astrophysics Data System. D.S. acknowledges support by the Deutsche Forschungsgemeinschaft through SPP 1385: “The first ten million years of the solar system − a planetary materials approach” (SE 1962/1-2 and SE 1962/1-3). TA acknowledges funding from the European Community’s Seventh Framework Programme [FP7/2007-2013] under grant agreement No. 238258. In the reduction and analysis of the data we made use of the GILDAS software public available at http://www.iram.fr/IRAMFR/GILDAS. The SMT and Kitt Peak 12 m are operated by the Arizona Radio Observatory (ARO), Steward Observatory, University of Arizona, with support through the NSF University Radio Observatories program (URO: AST-1140030).

References

Online material

Appendix A

Appendix A.1: Molecular column densities

We calculated the molecular column densities of the upper level of a particular transition following the equation for optically thin emission (A.1)where the integrated intensity given in K km s^-1, and the Einstein coefficient A_ul is in s^-1, and with (A.2)where ν is the frequency of the observed transition (e.g., Mangum & Shirley 2015). Then the total column density (which we refer to from now on as column density) can be calculated: (A.3)When we were able derive the optical depth τ from hyperfine splitting, we corrected for the column density using (A.4)

To obtain abundances, we derived H₂ column densities either from dust maps obtained with Mambo with the IRAM 30 m telescope at 1.2 mm (Beuther et al. 2002b) which provides a resolution of 11″, or the galactic plane survey ATLASGAL (Schuller et al. 2009) at 870 μm with a resolution of 19.2″, or the SCUBA Legacy Catalog (Di Francesco et al. 2008) at 850 μm with a resolution of 22.9″ (see Table A.1). The continuum data were smoothed to 29″ resolution to be beam-matched with the IRAM 30 m observations at 3 mm and the SMT molecular line data. H₂ column densities were calculated from the observed peak intensities assuming optically thin emission and LTE following Eq. (A.5) (Schuller et al. 2009). The dust opacities used were κ_{850 μm} = 1.48, κ_{870 μm} = 1.42, κ_{1.2 mm} = 0.97, interpolated values from Ossenkopf & Henning (1994), assuming grains with thin ice mantles, gas densities of n = 10⁵ cm^-3, and a gas-to-dust mass ratio R = 100. With these assumptions the H₂ column density is calculated as (A.5)The uncertainties in the derived H₂ column densities are mainly based on the dust and temperature properties and are about a factor of 3. A more detailed description of the derivation is given in Gerner et al. (2014).

The molecular column densities are then divided by the H₂ column densities, and the averaged abundances are derived.

All Tables

All Figures