© ESO, 2015

1. Introduction

Theory and observations suggest that the abundance of deuterated molecules in dense star-forming cores is related to the core evolution. The formation of deuterated molecules is favoured by the combination of low temperatures (T ≤ 20 K) and high densities (n ≥ 10⁴ cm^-3), which on one hand boosts the depletion of CO and other neutrals and, on the other, makes the relative abundance between a species containing D and its hydrogenated counterpart (the so-called deuterated fraction, D_frac) higher by three to four orders of magnitude with respect to the [D/H] interstellar abundance (~10^-5, e.g. Oliveira et al. 2003), due to the endothermicity of their backward reactions (see e.g. Millar et al. 1989; Gerlich et al. 2002). After protostellar birth, the young stellar object formed at the core centre heats up its surrounding material, and the temperature enhancement favours the progressive destruction of deuterated species and, consequently, makes D_frac decrease (see e.g. Caselli et al. 2002). Observations of low-mass star-forming cores have confirmed this theoretical scenario: both the column density ratio D_frac(N₂H⁺) and the column density of ortho-H₂D⁺, the parent species of most of the deuterated molecules formed in the gas (e.g. DCO⁺, N₂D⁺, DNC), increase in starless cores close to the onset of gravitational collapse, and then, after the formation of the protostar, they decrease as the core evolves (Crapsi et al. 2005; Emprechtinger et al. 2009; Caselli et al. 2008). Growing observational evidence suggests that high values of D_frac are also typical in high-mass star-forming cores (e.g. Fontani et al. 2006; Pillai et al. 2007, 2011; Miettinen et al. 2011), and that D_frac of some species could be also an evolutionary indicator in the intermediate- and high-mass regime (e.g. Busquet et al. 2010; Fontani et al. 2011; Sakai et al. 2012).

To investigate the relation between D_frac and core evolution in the high-mass regime in a systematic way, our team started a survey of deuterated molecules in about 30 dense cores that were carefully selected and almost equally divided among the three evolutionary phases in which we can roughly divide the high-mass star formation process observationally (see e.g. Beuther et al. 2007 and Tan et al. 2014): high-mass starless cores (HMSCs), high-mass protostellar objects (HMPOs) and ultracompact Hii regions (UC Hiis). In brief, the targets were selected as follows: The HMSCs had to be dense molecular cores not associated with indicators of star formation; the HMPOs had to show outflows, infrared sources, and/or faint (S_{3.6 cm}< 1 mJy) radio continuum emission; the UCHIIs had to be associated with stronger (S_{3.6 cm} ≥ 1 mJy) radio continuum. In selecting the sources, we rejected cores that are blended with nearby cores to avoid confusion and make the emission of the targeted core dominant.

In the first study, we (Fontani et al. 2011, hereafter Paper I) presented the results obtained from spectroscopic observations of millimetre rotational transitions of N₂H⁺ and N₂D⁺ obtained with the IRAM-30 m telescope, where we showed that D_frac(N₂H⁺) is ~0.26 in HMSCs, and drops by about an order of magnitude in the HMPO and UC Hii stages. These results are consistent with the fact that deuteration of N₂H⁺ starts from the reaction H₂D⁺ + N₂→ N₂D⁺ + H₂, which is efficient only at temperatures ≤20 K (Gerlich et al. 2002). In a subsequent study focused on DNC/HNC, Fontani et al. (2014, Paper II) showed that D_frac(HNC) also decreases from the pre- to the protostellar phase, but much more moderately, indicating that the ratio N₂D⁺-to-N₂H⁺ is more appropriate to identifying massive starless cores. This is consistent with the prediction that DNC can also easily form when the gas gets warmer, because the route reaction for the deuteration of HNC is linked to CH₂D⁺, which can stay abundant up to temperatures of 70 K (e.g. Leurini et al. 2006). However, N₂H⁺, HNC and their deuterated isotopologues can form mainly (HNC, DNC) or solely (N₂H⁺, N₂D⁺) in the gas phase. Other important molecules, such as NH₃, H₂CO, CH₃OH and their deuterated forms, can be produced on dust grain surfaces (e.g. Aikawa et al. 2005), and theoretical models show that this can make relevant differences in their D_frac, especially during the protostellar phase in which grain mantles evaporate (Aikawa et al. 2012).

In this work we investigate the role of surface chemistry by means of measurements of D_frac(NH₃) and D_frac(CH₃OH). Because methanol and its deuterated forms can be produced only on grain surfaces (see e.g. Parise et al. 2002; Garrod et al. 2007), D_frac(N₂H⁺) and D_frac(CH₃OH) represent the two “extreme conditions” under which deuteration can occur: only on grain surfaces (CH₃OH) and only in gas (N₂H⁺). Therefore, the results obtained in this work and in Paper I can be used as reference for the deuteration process of any other species formed potentially both in the gas and on dust grains. In Sect. 2 we present the source sample and give an overview of the technical details of the observations. The main results are presented and discussed in Sects. 3 and 4, respectively. A summary with the main conclusions of the work is given in Sect. 5.

2. Source list, observations and data reduction

2.1. Source list

We targeted the same sources as studied in Paper I, to avoid possible biases due to the source selection when comparing the deuterated fractions. Table 1 contains the list of the observed sources selected as explained briefly in Sect. 1. In particular, three HMSCs have been classified as “warm” cores because they show evidence of heating from external sources (see Paper I for details). More information extracted from the literature about the star-forming regions in which the sources lie are given in Table A.1 of Paper I. To the list of HMSCs reported in Paper I, we have added the source G028–C3, selected by applying the same selection criteria as for the other HMSCs.

2.2. IRAM-30 m observations

Run-1: towards all sources in Table 1, observations of the ortho- and para-NH₂D(1_1,1−1_0,1) line were obtained simultaneously to the N₂D⁺ and N₂H⁺ observations described in Paper I. Table 2 lists the main observational parameters. We refer to Sect. 2 of Paper I for any other technical detail related to these observations.

Run-2: we performed CH₃OH and CH₂DOH observations towards all sources observed in Paper I from 6 to 9 February 2013. We simultaneously observed two bands at 3 and 1.3 mm, covering some important rotational transitions of CH₃OH, ¹³CH₃OH and CH₂DOH. Table 2 presents the observed spectral windows and some main technical observational parameters. The atmospheric conditions were very stable during the whole observing period, with precipitable water vapour usually below ~2 mm. The observations were made in wobbler-switching mode. Pointing was checked almost every hour on nearby quasars or bright Hii regions. The data were calibrated with the chopper wheel technique (see Kutner & Ulich 1981), with a calibration uncertainty of ~20%. The spectra were obtained in antenna temperature units, , and then converted to main beam brightness temperature, T_MB, via the relation ). The spectra were obtained with the fast Fourier transform spectrometers (FTS), thereby providing a broad band of ~8 GHz simultaneously at 3 and 1.3 mm (see Table 2 for details). All calibrated spectra were analysed using the GILDAS¹ software developed at the IRAM and the Observatoire de Grenoble. The rest frequencies used for the line identification have been taken from the Cologne Molecular Database for Spectroscopy² (CDMS; Müller et al. 2001, 2005)

2.3. GBT observations

The ammonia (1,1) and (2,2) inversion transitions (rest frequencies 23.6944955 and 23.7226336 GHz, respectively) were observed with the 100 m Robert C. Byrd Green Bank Telescope³ (GBT) during 13 and 21 March and 4 and 21 April 2013. The GBT spectrometer backend was configured to simultaneously observe the two transitions in separate spectral windows, using bands of 50 MHz and spectral resolution of 12.2070 kHz, which corresponds to 0.154 km s^-1 for both lines. The main observational parameters are listed in Table 2. The data were taken using in-band frequency switching with a throw of 7.5 MHz. The beam FWHM was approximately 32′′. The pointing was checked at hourly intervals on a nearby quasar, with corrections approximately 2′′–3′′. Flux calibration was performed on 3C 123 and NGC 7027. The absolute flux accuracy is 10% to 20%. Data reduction and calibrations were done using the GBTIDL⁴ package and subsequently converted to CLASS format.

3. Results and derivation of physical parameters

3.1. NH₃ and NH₂D

3.1.1. Detection summary and parameters derived directly from the fits

NH₃: the NH₃(1,1) and (2,2) inversion lines were detected with excellent signal-to-noise ratio in all sources observed. The spectra of all HMSCs, HMPOs, and UC Hiis are shown in Figs. A.1–A.3, respectively. Both transitions consist of 18 hyperfine components, grouped in five lines: the main one at the centre of the spectrum, and four satellites symmetrically placed in frequency with respect to the main one (see Ho & Townes 1983 for details). The spectra have been fit considering this hyperfine structure when the satellites are detected. When they are not, we adopted a simplified approach in which we fitted the main line with a Gaussian curve. This simplified method was used for eight of the (2,2) spectra observed, in which the satellites have not been detected. The fit procedure has given good results (with very low residuals, see Figs. A.1–A.3) using both methods. The simplified approach tends to overestimate the intrinsic line width, because the main line is in reality a blending of several hyperfine components. To quantify this, we have taken a (2,2) spectrum with hyperfine structure fit nicely (one spectrum per evolutionary group), applied the simplified method, and compared the derived line width with the one obtained from the accurate method (the hyperfine fit method). From this comparison, we quantify an overestimate of at most the 10% of the true intrinsic line width. Nevertheless, because the column density in this approach is computed from the integral of the line (see Sects. 3.1.2 and 3.1.3), this overestimation does not influence the calculation of either the column density or the deuterated fraction.

The line parameters derived from these fit procedures are listed in Tables 3 and 4. The accurate method has given a well-constrained value of the optical depth of the main component of the (1,1) line (τ_m(1,1)/Δτ_m(1,1) ≥ 3) for all objects except for I04579–VLA1, in which the line is optically thin. The average τ_m(1,1) is ~1 with no significant differences between the three evolutionary groups, while τ_m(2,2) is usually smaller than 1. The average line widths of the (1,1) lines are ~1.7, ~2.3, and ~2.6 km s^-1 for the HMSC, HMPO, and UC Hii groups (standard deviations 0.5, 0.6, and 0.7 km s^-1, respectively), and tend to increase with evolution, as expected (Sánchez-Monge et al. 2013).

NH₂D: the ortho-NH₂D(1_1,1−1_0,1) line has been detected in all the observed sources except in two HMPOs (I04579–VLA1 and I21307). The detection rate is thus ~92%. The para-NH₂D(1_1,1−1_0,1) line has been detected in 13 out of 26 sources observed (detection rate of 50%). The spectra of both lines are shown in Figs. A.4 and A.5. In this work we use the ortho-NH₂D(1_1,1−1_0,1) line to derive the physical parameters of our interest because of its higher detection rate and signal-to-noise ratio. The line of the para- species will be used to test whether the ortho-/para- ratio assumed to derive the total column density in Sect. 3.1.3 is correct. Like ammonia, the NH₂D lines have been fit by taking their hyperfine structure driven by the quadrupole moment of the deuterium and nitrogen nuclei into account (see Olberg et al. 1985, for details).

In general, the procedure has provided good fits to the spectra, except in a few cases in which deviations from the LTE (symmetric) pattern are seen (e.g. G028–C1, I20293–WC, I20293–MM1, 23033+5951, see Fig. A.4). To check if (and how) our simplified local thermodynamic equilibrium (LTE) approach gives results that are different from those of a non-LTE analysis, we ran the non-LTE radiative transfer code RADEX⁵ (Van der Tak et al. 2007) in order to reproduce the measured line ratios of the two lines (ortho-NH₂D and para-NH₂D). The molecular data were taken from the LAMDA database (Schöier et al. 2005) using the collisional rate coefficients with H₂ of Daniel et al. (2014). We built grids of models with kinetic temperatures in the range 8–25 K, H₂ volume densities in the range 10³−10⁸ cm^-3, and total column densities of 10¹²−10¹⁵ cm^-2. We assumed line widths of 1.5 km s^-1 and an ortho-to-para- ratio of 3. The “best estimate” of the column densities that we find are roughly consistent with the values measured from the LTE approach, but since we only have one line ratio, we cannot properly distinguish between the different non-LTE models. Therefore, with only one line ratio, all we can say is that the column densities of NH₂D derived assuming LTE conditions are consistent with the values expected from a non-LTE approach.

The average optical depth of the main hyperfine component derived from this fitting procedure is about one in all three evolutionary groups. For the sources for which the mentioned fitting procedure did not give good results (because of poor signal-to-noise ratio), we fit the lines with Gaussians. As for the NH₃ lines, this simplified method could overestimate the line widths by at most ~10%, and we find yet an increasing trend of the line widths going from the HMSC phase to the HMPO and UC Hii phases, for which mean values (and standard deviations) are: 1.4(0.6), 2.5(1.3) and 2.4(1) km s^-1, respectively. All line parameters are listed in Table 5.

3.1.2. NH₃ rotation temperature and total column density

From the NH₃(1,1) and (2,2) line parameters, we obtained rotation temperatures, T_rot, adopting two methods: for the nine sources having τ_(2,2)/Δτ_(2,2) ≥ 3 and τ_(2,2)> 0.1, we have first derived the excitation temperature of the (1,1) and (2,2) lines (T_ex1,1 and T_ex2,2, respectively) independently using Eq. (A.2) of Busquet et al. (2009), and, from these, the column densities of the two levels, N_(2,2) and N_(1,1), from the relations given in Anglada et al. (1995). Although Eq. (A.2) in Busquet et al. (2009) is derived for the (1,1) line, it is also valid for the (2,2) line given the small difference in frequency between the two transitions. Then, the rotation temperature has been derived from the relation: (1)For sources with an optically thin (2,2) line, or with τ_(2,2) not determined because the satellites are undetected, T_ex2,2 is assumed to be equal to T_ex1,1. This hypothesis is justified by the good agreement between the two excitation temperatures in the sources in which they can both be measured (see Sect. 4.2). Under this assumption, we have applied Eq. (A.4) in Busquet et al. (2009), which utilizes the peak intensity of the main hyperfine component of the (2,2) line.

In both methods, the total NH₃ column density, N(NH₃), has been calculated from Eq. (A.6) in Busquet et al. (2009). Both T_rot and N(NH₃) are listed in Table 7. Rotation temperatures range from 11.7 to 29 K, and on average they are ~17, ~22, and ~22 K for HMSCs, HMPOs, and UC Hiis, respectively (standard deviations are 2.6, 3.5, and 4 K, respectively). Separately, quiescent and “warm” HMSCs have mean temperatures of 16 and 20 K (standard deviations of 2.6 and 1.2 K, respectively), which confirms the higher gas temperature in the “warm” cores. Total NH₃ column densities range from 5.6 × 10¹³ to 3.6 × 10¹⁵ cm^-2, and the average values are 9.4 × 10¹⁴, 9.3 × 10¹⁴ and 1.6 × 10¹⁵ cm^-2 in the HMSC, HMPO and UC Hii groups, respectively. We assumed a filling factor of one because available VLA interferometer ammonia maps of some of the targets show that the ammonia emission is extended and fills most of the GBT beam. Nevertheless, we stress that the emission from the target cores is clearly dominant with respect to the emission of nearby condensations (see Sanchez-Monge et al. 2013).

3.1.3. NH₂D total column density

The NH₂D column densities were computed from the line parameters of the ortho-NH₂D line following Eq. (1) in Busquet et al. (2010), which assumes the same T_ex for all the hyperfine components. The value of T_ex was computed as described in Sect. 3.1.2 for sources where the opacity of the main component is well-constrained. For the others, we have assumed T_ex = 7.5 K, which is the average value derived from the sources with well-constrained opacity and obtained the column density from Eq. (A4) of Caselli et al. (2002), which is valid for optically thin lines.

Again, we assumed a unity filling factor because there are few high angular resolution observations of this line towards the targets from which the emitting region of NH₂D can be determined. This assumption is critical, since the ortho-NH₂D line has a critical density of ~10⁶ cm^-3, so higher than that of the inversion transitions of NH₃ (~10³⁻⁴ cm^-3). However, while by neglecting the beam dilution the absolute values of the column densities can certainly be affected, the evolutionary trend of the column density ratio should not be affected by this assumption because the beam dilution is expected to be almost constant, and it should thus introduce only a systematic correction (see also Paper I). Also, observations at high angular resolution towards massive star forming regions (Busquet et al. 2010; Pillai et al. 2011) indicate that the emission of NH₂D can be as extended as that of NH₃, despite the different critical density. For example, the emitting regions of NH₂D(1_1,1−1_0,1) and NH₃ (1,1) in I20293–WC and I20293–MM1, both included in our survey, are approximately the same (Busquet et al. 2010). N(NH₂D) is listed in Table 7.

3.2. Methanol and deuterated methanol lines

We focus the attention on the deuterated fraction of CH₃OH and on the physical quantities relevant to deriving it (i.e. temperature and total column density). Therefore, in what follows we present the approach we adopted to identify the lines from which D_frac(CH₃OH) will be derived (Sect. 3.2.1), the method for computing rotation temperature and total column density from the line parameters (Sect. 3.2.2), and the deuterated fraction in the sources detected in CH₂DOH (Sect. 4.4).

3.2.1. Lines detected and fit procedure

Multiple CH₃OH lines are detected in the observed spectral windows (Col. 1 of Table 2) towards all sources, while ¹³CH₃OH lines are clearly detected in four HMSCs, seven HMPOs, and six UC Hii regions. CH₂DOH lines are detected towards 6 objects only, three HMSCs and three HMPOs, and two out of the three HMSCs are “warm” cores (see Sect. 2). Moreover, in two HMPOs (AFGL5142–MM and 18089–1732), the CH₃OD (5_1,5−4_1,4A++) line at 1.3 mm has been detected, although in 18089–1732 this could be blended with emission of (CH₂OH)₂. Tables B.1 and B.2 give the line parameters obtained from Gaussian fits to the lines detected at the 3σ level and not probably blended with other transitions.

The detection of the deuterated lines was double-checked by comparing observed and synthetic spectra. For this purpose, the observed spectra were smoothed to 1.0 km s^-1 at 1 mm and to 2.5 km s^-1 at 3 mm, to improve the signal-to-noise ratio. The CH₂DOH synthetic spectra were computed by assuming LTE and optically thin emission as in Palau et al. (2011) and using the molecular data from the Jet Propulsion Laboratory (Pickett et al. 1998). To build the synthetic spectra, we adopted a line width of 1.5 km s^-1 at 1 mm and 2.5 km s^-1 at 3 mm and used the rotational temperature listed in Table 8 derived from CH₃OH. Examples of the synthetic spectra of CH2DOH overplotted on the observed spectra are shown in Fig. 1. The figure shows that several transitions are (marginally) detected at 1 mm in each of the four cases shown.

3.2.2. Derivation of molecular column densities and rotation temperatures

From the line parameters in Tables B.1 and B.2, we derived rotation temperature (T_rot) and total column densities, N, of CH₃OH, ¹³CH₃OH and CH₂DOH from the rotation diagram method. As an example, in Fig. 2 we show the rotation diagrams obtained for 18089–1732. We included all rotation diagrams in Appendix C. The method was applied when the number of transitions detected was sufficient to build a “reliable” rotation diagram: for example, we rejected the results obtained from this method for sources in which rotation diagrams provide meaningless negative temperatures or for objects in which a few lines associated with large uncertainties and/or similar energy of the upper levels have been detected. Specifically, for CH₂DOH the rotation diagram method has given acceptable results for only two sources, AFGL5142–MM and 18 089–1732. However, because in AFGL5142–MM we only have two lines, and in 18089–1732 the fit results are not very accurate (bottom panel in Fig. 2), the column densities have also been derived from Eq. (A4) of Caselli et al. (2002), taking the strongest line detected and assuming the gas temperature equal to T_rot derived from CH₃OH.

For ¹³CH₃OH, we assumed that all transitions are optically thin; for CH₃OH, we first derived a rough estimate of the opacity by comparing two identical lines (specifically, we compared the 2_{(− 1,2)}−1_{(− 1,1)} and the 2_(0,2)−1_(0,1) transitions) of CH₃OH and ¹³CH₃OH and assumed an LTE abundance ratio of [¹²C]/[¹³C] = 77 (Wilson & Rood 1994). From this check, we derived low opacities (values below 1) in all sources, so that we decided to compute N and T_rot assuming optically thin conditions too. As for NH₃ and NH₂D, the source sizes of CH₃OH and CH₂DOH are unknown, but they are expected to be smaller than the beam size and to have a comparable extent, based on observations at high angular resolution in Orion (Peng et al. 2012). Therefore, to take the beam dilution into account, the column densities in the rotation diagrams were corrected by assuming the same source sizes as in Paper I, namely 6.5, 4.1, and 5.5′′ for HMSCs, HMPOs, and UC Hiis, assuming that methanol and its deuterated forms trace approximately the same material. This is a reasonable general assumption from a theoretical point of view if deuterated methanol is formed from methanol through H–D substitution reactions on dust grains. In principle, CH₂DOH could be formed following other pathways, but the H–D substitution reaction on solid ices remains the most efficient one (Nagaoka et al. 2005). Moreover, because of the lack of direct measurements in the cores, assuming a different source size for methanol and their deuterated forms would be an arbitrary choice not supported by observations.

For the deuterated species for which only one line is detected, and for sources in which N(¹³CH₃OH) cannot be derived from rotation diagrams, we derived N using Eq. (A4) of Caselli et al. (2002) from one transition only assuming optically thin conditions and adopting the rotation temperature derived from CH₃OH, available in all sources as excitation temperature. The partition functions of all species have been calculated from the approximated expressions valid for asymmetric rotors provided, e.g., in Ratajczak et al. (2011; see also Parise 2004). The results of this analysis are presented in Table 8.

4. Discussion

4.1. The ortho-/para- ratio of NH₂D

The total column density of NH₂D has been derived from lines of ortho-NH₂D taking the statistical ortho-/para- ratio (3:1) into account. In the sources also detected in the para-NH₂D line at ~110 GHz (see Fig. A.5), we have verified whether the assumption is correct: First, we fitted the hyperfine structure of the para-NH₂D line, and found that all detected lines are optically thin. Because most of the detected ortho-NH₂D lines are either optically thin or have τ ≤ 1, we decided to compare the integrated areas of the two transitions under the channels with signal. These are reported in Table 6. As we see, the ratio is consistent with three within the errors in most of the sources: the mean value is 2.6, with standard deviation 0.6, hence consistent with three.

Shah & Wootten (2001) have found similar results in a sample of protostellar cores, in which they compare the integrated intensity of the same two transitions and derived a mean value of the ortho-/para- ratio of 3.2 (with a larger standard deviation of ~1.3). Comparable values have also been found by Pillai et al. (2007) in infrared dark clouds and by Tiné et al. (2000) in the two dark molecular clouds L134N and TMC1.

4.2. NH₃ and NH₂D excitation temperatures

The excitation temperatures of the three lines examined in the previous sections have very similar mean values: 7.8, 8.2 and 7.5 K for NH₃(1,1), NH₃(2,2) and ortho-NH₂D, respectively. The NH₃(1,1) and (2,2) lines are also well correlated (see upper panel of Fig. 3), while the excitation temperatures of the ortho-NH₂D line and of NH₃(1,1) are not correlated owing to the different dispersion around the average values: in fact, T_ex of NH₃(1,1) spans a range from ~3.5 K to 15 K, while T_ex of the ortho-NH₂D lines is distributed tightly around the average value. This may indicate either that ortho-NH₂D is in sub-thermal conditions, as suggested by the asymmetric pattern of the hyperfine structure observed in some spectra (see Sect. 3.1.1), or to our neglect of the correction for beam dilution. The former hypothesis seems plausible for the ortho-NH₂D line, which has a high critical density (~10⁶ cm^-3). About the possible different beam dilution: as stated in Sect. 3.1.3, in the few cores in which both NH₃(1,1) and ortho-NH₂D(1_{− 1,1}−1_0,1) have been mapped at high angular resolution, the emissions have comparable extensions, despite the different critical densities. Therefore, sub-thermal conditions of the ortho-NH₂D lines seem the most likely explanation of the different excitation temperatures.

4.3. Deuterated fraction of NH₃

By dividing N(NH₂D) for N(NH₃), we computed D_frac(NH₃). The three parameters are given in Table 7. The average values of D_frac(NH₃) for HMSCs, HMPOs, and UC Hiis are 0.26 (0.23 if one excludes the “warm” HMSCs, see Sect. 2), 0.34, and 0.21, respectively. These values are consistent with those obtained by Pillai et al. (2007) in a sample of infrared-dark clouds, for which, however, the evolutionary stage of the embedded sources has not been determined. The mean D_frac(NH₃) is thus at the maximum in the HMPO stage, although the large disperion of the data does not allow finding a statistical difference between the three groups. This is apparent in Fig. 4, where we compare the total column densities of NH₂D and NH₃: the plot shows that the three groups are not clearly separated. Kolmogorov-Smirnov tests applied to the data confirm that the difference is not statistically significant. If one compares Fig. 4 with the same plot in Paper I for N₂H⁺, we note clearly that, unlike D_frac(N₂H⁺), D_frac(NH₃) does not decrease with core evolution. Thus, it is not a tracer of pre-protostellar or young protostellar objects, because it keeps above 0.1 even in the evolved stage of UC Hii region. Moreover, because the deuteration in the gas-phase for both N₂H⁺ and NH₃ is linked to H₂D⁺, our results would confirm that the formation of NH₂D is largely influenced by surface chemistry.

Furthermore, D_frac(NH₃) does not show any clear anti-correlation with the typical indicators of evolution. This is suggested by Figs. 5 and 6, where we plot D_frac(NH₃) against the gas temperature and the line widths of the (1,1) transition, both known to increase with time (e.g., Sánchez-Monge et al. 2013): by applying statistical tests, we even find that D_frac(NH₃) could be slightly correlated with both the ammonia rotation temperature (Pearson’s linear correlation coefficient ρ ~ 0.2) and line widths (Pearson’s linear correlation coefficient ρ ~ 0.4). In contrast, D_frac(N₂H⁺) is anti-correlated with both parameters, as shown in Fig. 2 of Paper I. We stress, however, that the p-value (measure of the probability of chance correlation) is 0.12 for D_frac(NH₃) vs. T_rot and 0.22 for D_frac(NH₃) vs. Δv(1,1). Therefore, because the significance level under which one can reject the null hypothesis is typically p ~ 0.1, both correlations are very weak from a statistical point of view. Nevertheless, the relevant result provided by Figs. 5 and 6 is the absence of anti-correlation, contrary to what found for D_frac(N₂H⁺).

4.4. Deuterated fraction of methanol

In the six objects detected in CH₂DOH, we computed D_frac(CH₃OH) = N(CH₂DOH)/N(CH₃OH). The results are listed in the last column of Table 8. The average D_frac(CH₃OH) in the three HMPOs detected is ~0.04 if one uses N(CH₂DOH) derived from rotation diagrams, and ~0.01 if we use the simplified approach from only one line (see Sect. 3.2.2). In the two “warm” HMSCs is ~0.0025. G034–G2 is the unique quiescent HMSC detected in CH₂DOH, and in this core D_frac(CH₃OH) is ~0.015. For the cores undetected in CH₂DOH, the large majority of the targets observed, we have estimated upper limits of D_frac(CH₃OH) in this way: we computed the 3σ level in the spectrum of the (2_0,2−1_0,1)e0 line, which is the transition with the lowest energy of the upper level (E_u ~ 6.5 K) at 3 mm, and estimated the upper limit to the integrated area from the relation , valid for a Gaussian line having peak temperature =3σ. We assumed ΔV = 1 km s^-1, which is the average value of the detected CH₂DOH lines both in the HMSCs and the HMPOs (see Table B.1); then, the upper limit on the CH₂DOH column density was computed using the same equations as for the detected sources. We followed the same method so as to compute the upper limits on the ¹³CH₃OH lines, using this time the (2_0,2−1_0,1 + +) line.

The case of the HMSC G034–G2 is quite peculiar, because it is the only quiescent starless core detected in CH₂DOH (in one line only), and its D_frac(CH₃OH) exceeds 0.01, while the upper limits found in the other quiescent HMSCs are lower. Its detection is thus quite doubtful. We checked for possible contamination from other species by running synthetic spectra (see Sect. 3.2.1) of molecules that possess transitions at a similar frequency, and concluded that indeed contributions from lines of CH₃OOH, and HCCCH₂OH are possible. Therefore, this detection remains doubtful.

Although the number of detections is low, and the results are affected by the faintness of the CH₂DOH lines, these findings suggest that high values of D_frac(CH₃OH) tend to be associated with “warm” HMSCs and HMPOs rather than with cores that are very young (quiescent HMSCs) or evolved (UC Hiis), although the remarkable value derived in G034–G2 (if confirmed) suggests that the story could not be so simple. Parise et al. (2006) measured values of D_frac(CH₃OH) higher than ours by at least an order of magnitude in a sample of low-mass protostellar cores. Nevertheless, due to the higher linear resolution of their observations (most of their cores are in Perseus, at a distance of ~200 pc), their measurements should be less affected than ours by non-deuterated gas along the line of sight. Moreover, our D_frac(CH₃OH) are consistent with the upper limits found by Loinard et al. (2003) in high-mass protostellar objects (where, however, they observed D₂CO and derived [D₂CO] / [H₂CO] < 0.5%), as well as with observations of deuterated methanol in the intermediate-mass protostar NGC 7129-FIRS2 (Fuente et al. 2014) and in Orion BN/KL (Peng et al. 2012).

4.5. Deuteration and core evolution: the role of surface chemistry

In Fig. 7 we report the mean values of D_frac obtained in HMSCs, HMPOs and UC Hiis for the four molecular species investigated so far towards our source sample: N₂H⁺ (Paper I), HNC (Paper II), NH₃ and CH₃OH (this work). We show the values derived for quiescent HMSCs and “warm” HMSCs separately to underline the effect of nearby star formation. We also include the mean values (with standard deviation) of the ammonia rotation temperatures derived in this work to highlight a possible (anti-)correlation between D_frac and gas temperature. Inspection of Fig. 7 leads to these immediate results: (i) only D_frac(N₂H⁺) shows a net decrease from the HMSC stage to the HMPO stage, associated to a temperature enhancement; (ii) D_frac(NH₃) remains nearly constant in all stages; (iii) D_frac(CH₃OH) is at its maximum in the HMPO stage, although this result must be interpreted carefully due to the low number of detections and the caveats on the methods to derive D_frac (see Sects. 3.2.2); (iv) the behaviour of D_frac(HNC) is something in between D_frac(N₂H⁺) and D_frac(NH₃), because its maximum value is found in the HMSC phase, like D_frac(N₂H⁺), but the statistically significant decrease when going to the HMPO stage is not seen. In Paper II we have already discussed this difference, and attributed it to a slower process of destroying DNC into the warm gas with respect to N₂D⁺.

As stated in Sect. 1, D_frac(N₂H⁺) and D_frac(CH₃OH) should represent the two “extreme” situations where deuteration can occur: only in the gas and only on grain mantles, respectively. In the classical framework, both ammonia and methanol (and their deuterated isotopologues) are produced efficiently on grain mantles during the pre-stellar phase through hydrogenation of N and CO, respectively. Specifically, hydrogenation of CO sequentially forms formaldehyde first and then methanol: thus, as time proceeds, the formation of methanol and their deuterated isotopologues is boosted, until the energy released by the nascent protostellar object in the form of radiation increases the temperature of its environment, causing the evaporation of the grain mantles and the release of these molecules into the gas. As the temperature increases and the protostar evolves towards the UC Hii region phase, the deuterated species are expected to be gradually destroyed because of the higher efficiency of the backward endothermic reactions (see Caselli & Ceccarelli 2012, for a review). The trends shown in Fig. 7 are consistent with this classic framework, and they clearly show that high deuterated fractions of ammonia cannot be used as an evolutionary indicator of a high-mass star-forming core. On the other hand, D_frac(CH₃OH) may potentially be a tracer of the very early stages of the protostellar evolution, where the evaporation or sputtering of the grain mantles is most efficient. Our results, however, suffer from statistics that are too low and need to be reinforced by other observations of deuterated methanol at higher sensitivity.

Chemical models of low-mass star-forming cores predict how the abundance of several deuterated species varies during the evolution, including the amount formed on ices during the early cold phase (e.g. Taquet et al. 2012; Aikawa et al. 2012). Aikawa et al. (2012) show that the relative abundance ratios [NH₂D]/[NH₃] and [CH₂DOH]/[CH₃OH] in the ices during the pre–stellar phase are both in between 0.01 and 0.1. These values are consistent with the D_frac(CH₃OH) measured in this work, and they confirm that methanol and its deuterated forms are products of the evaporation of grain mantles. On the other hand, the D_frac(NH₃) measured in our work (≥0.1) is larger than the values predicted on ices by Aikawa et al. (2012), suggesting that the emission we see must include a contribution from material formed through gas-phase reactions. Awad et al. (2014) modelled the deuterium chemistry of star-forming cores using both gas-phase and grain-surface reactions, but focus on the protostellar phase, when the evaporation of the icy mantles of dust grains is at its maximum. The model that best reproduces a HMPO predicts D_frac(NH₃) ~10^-3−10^-2 and D_frac(CH₃OH) ≤4 × 10^-3, both of which are smaller than our observed values. However, the abundance of deuterated species strongly depends on the density of the gas: lower density cores have lower abundances of deuterated species, owing to a lower degree of CO depletion. Therefore, higher core densities could be able to reproduce the larger deuterated fractions that we measure.

In any case, the huge dispersion of the data does not allow us to derive firmer quantitative conclusions and push us to interpret any comparison with chemical models with caution. Moreover, the chemical models of Taquet et al. (2012) and Aikawa et al. (2012) neglect the spin states of the deuterated species, which can significantly influence the deuterium fractionation depending on the ortho-to-para H₂ ratio (Flower et al. 2006). Nevertheless, the clearly different trends between D_frac(NH₃) and D_frac(N₂H⁺) undoubtedly indicate that gas-phase chemistry cannot play a dominant role in the production of NH₂D.

5. Summary and conclusions

The deuterated fraction of species that can be formed on dust grains (in part, like NH₃, or uniquely, like CH₃OH) has been investigated towards a sample of dense cores harbouring different evolutionary stages of the high-mass star formation process. As expected, the deuterated fraction of these species and those of molecules totally or predominantly formed in the gas, such as N₂H⁺ and HNC, evolve differently with time and with temperature: D_frac(NH₃) does not show statistically significant changes with evolution, unlike D_frac(N₂H⁺) and D_frac(HNC), which decrease (especially D_frac(N₂H⁺)) when temperature increases. Very few lines of CH₂DOH and CH₃OD are clearly detected, and only towards protostellar cores or externally heated starless cores. Only one line of CH₂DOH could have been detected in a quiescent starless core, but the detection is doubtful. No lines of deuterated methanol species are detected in UC Hii regions. This work clearly supports the scenario in which the contribution of surface chemistry to the formation of deuterated forms of ammonia is relevant, hence D_frac(N₂H⁺) remains the best indicator of massive starless cores. High values of D_frac(CH₃OH) seem suitable to tracing the earliest protostellar phases, at which the evaporation or sputtering of the grain mantles is most efficient, but this result needs to be supported by further, higher sensitivity observations. The data presented in this work represent an excellent starting point for higher angular resolution studies to address further questions. In particular: if the various deuterated molecules are formed with different mechanisms, do we expect a different distribution of the emission, too?

Appendix A: NH₃ and NH₂D spectra

Appendix B: Tables

Appendix C: Boltzmann plots for CH₃OH and ¹³CH₃OH

---

Acknowledgments

F.F. and A.P. are grateful to the IRAM-30 m staff for their help with the observations at the IRAM-30 m telescope. G.B. is grateful to Amanda Kepley for her help during the GBT observations and to Jeff Magnum for providing the procedures for converting the GBT spectra from GBTIDL to CLASS format. A.P. acknowledges the financial support from UNAM, and CONACyT, México. G.B. is supported by the Spanish MICINN grant AYA2008–06189–C03–01 (co-funded with FEDER funds) and by the Italian Space Agency (ASI) fellowship under contract number I/005/07/01. P.C. acknowledges the financial support of the European Research Council (ERC; project PALs 320620). A.S.-M. is supported by the Deutsche Forschungsgemeinschaft (DFG) through the collaborative research grant SFB 956 “Conditions and Impact of Star Formation”, project area A6. The research leading to these results has received funding from the European Commission Seventh Framework Programme (FP/2007-2013) under grant agreement Na 283393 (RadioNet3).

References

Online material

All Tables

All Figures