© ESO, 2017

1. Introduction

Ammonia belongs to the most useful probes of the dense cores of molecular clouds owing to its energy spectrum and chemical properties (Ho & Townes 1983; Benson & Myers 1989; Tafalla et al. 2002; Friesen et al. 2009). The molecule can survive in the gas phase also in the cold, dense interior parts of starless and pre-stellar cores, and in these regions, reactions with deuterated ions convert part of NH₃ to NH₂D, NHD₂ and ND₃ (Rodgers & Charnley 2001). Ammonia and its deuterated isotopologues are also formed on grain surfaces through H/D-atom addition reactions to N atoms (Brown & Millar 1989; Fedoseev et al. 2015). Substantial deuteration in both phases occurs after the disappearance of CO from the gas phase, and needs some time to take effect, so the relative abundances of the mentioned molecules can give an idea of the evolutionary stage of a dense core (Roueff et al. 2005; Flower et al. 2006a; Roueff et al. 2015). The N₂D⁺/ N₂H⁺ abundance ratio has already been used extensively for this purpose (Crapsi et al. 2005; Pagani et al. 2009b).

Besides enabling nuanced investigation into deuteration, spectral line observations of NH₃, NH₂D, NHD₂ and ND₃ test our understanding of spin-state chemistry, that is, selection rules in chemical reactions that determine the relative abundances of different nuclear spin symmetry species or “modifications”. Each spectral line observed from these molecules belongs to one of the two spin modifications of NH₃, NH₂D or NHD₂, or one of the three spin modifications of ND₃. The abundance ratios of these modifications, which must be treated as separate chemical species, can deviate from their nuclear spin statistical ratios and are predicted to change with time (Faure et al. 2013; Sipilä et al. 2015a,b).

The spin symmetry species of a molecule is defined based on how the nuclear spin wave function transforms under symmetry operations such as the interchange or the permutation of identical nuclei. For molecules containing only H nuclei, there is a one to one correspondence between the symmetry and the nuclear spin angular momentum, but this is not true for molecules with multiple D nuclei (Hugo et al. 2009). Schmiedt et al. (2016) have recently shown that the spin angular momentum and the permutation symmetry are, after all, inherently coupled, and that it is possible to construct an unequivocal and practical representation for each spin angular momentum – symmetry species using the Young diagrams. In the present study, chemical species are distinguished solely by their nuclear spin symmetries.

The statistical weight of a symmetry species is the number of possible nuclear spin functions that have that symmetry. When there are no more than three spin modifications, it is customary to call the species with the largest nuclear spin statistical weight “ortho”, and the one with the lowest weight “para”. If there is a third one, this is called “meta”. For example, the three elementary spin functions of the deuteron, | 1,1 ⟩, | 1,0 ⟩ and | 1,−1 ⟩, can combine in ND₃ in 27 different ways, and these combinations can be arranged into a set of 27 orthogonal functions, which form an irreducible representation of the appropriate permutation group S₃. Of these functions, 10 have symmetry A₁ (“meta”), 1 has symmetry A₂ (“para”), and 16 have symmetry E (“ortho”) in S₃ (Bunker & Jensen 2006; Sipilä et al. 2015b; Daniel et al. 2016b, see their Appendix A).

The abundances of all three deuterated forms of ammonia have been previously estimated in the protostellar core “B1-b” of Barnard 1 in Perseus (Lis et al. 2002b, 2006; Roueff et al. 2005), and in the starless core (“I16293E”) of the dark cloud L1689N in Ophiuchus (Roueff et al. 2005; Gerin et al. 2006). The fractionation ratios are similar in these two objects. The recent analysis of these observations by Daniel et al. (2016b) gives [ NH₂D ] / [ NH₃ ] ≳ [ NHD₂ ] / [ NH₂D ] ≈ 0.2, [ ND₃ ] / [ NHD₂ ] ≈ 0.05−0.10. In both cores, the ortho/para ratios of NH₂D and NHD₂ were found to be close to their statistical values, 3:1 and 2:1, respectively. The interpretation of chemical data from these two, bright sources of deuterated ammonia is complicated by the possibility of enhanced evaporation from grains due to shocks. The Barnard 1 core contains two, albeit very young and low-mass protostars driving outflows (B1-bS and B1-bN; Hatchell et al. 2007; Hirano & Liu 2014; Gerin et al. 2015). The core in L1689N may be compressed by an outflow from the adjacent protostellar source IRAS 16293-2422 (Lis et al. 2002a; Stark et al. 2004; Gerin et al. 2006).

In the present work, we use the 100-m Robert C. Byrd Green Bank telescope, the 12-m APEX and the IRAM 30-m telescope to determine the abundances of para-NH₃, meta-ND₃ and both ortho and para modifications of NH₂D and NHD₂ in the starless core H-MM1 in Ophiuchus (Johnstone et al. 2004; Parise et al. 2011). Through these observations, we obtain a homogenous data set pertaining to a quiescent region, presumably characterised by a simple physical structure, where comparison with chemistry models is more straightforward than in a star-forming core. The observations are interpreted by modelling the chemical evolution of a hydrostatic core resembling H-MM1, and simulating observations towards this model core. Besides the previously published collisional coefficients of NH₂D with para-H₂ (Daniel et al. 2014), we use newly calculated coefficients for NHD₂ and ND₃ (Daniel et al. 2016b), in conjunction with the radiative transfer program of Juvela (1997).

The paper is organised as follows. In Sect. 2, we describe the target core and give details of the observations. The direct observational results are presented in Sect. 3. In Sect. 4, we construct a physical model of the core and describe the modelling tools, the collisional coefficients and the chemistry model used in this work. In Sect. 5, we make predictions for the NH₃, NH₂D, NHD₂ and ND₃ abundances and the observable line emission from the model core. In Sect. 6, we compare the outcome of the modelling with the observations and discuss the implications of this comparison. Finally, in Sect. 7, we draw our conclusions.

2. Observations

2.1. The target: H-MM1

The dense core H-MM1 in Ophiuchus was discovered by Johnstone et al. (2004) using the Submillimeter Common-User Bolometric Array (SCUBA) on the James Clerk Maxwell Telescope (JCMT). The core was also covered by the JCMT Gould Belt Survey with SCUBA-2 at 450 and 850 μm (Pattle et al. 2015). H-MM1 lies in relative isolation in the eastern part of Lynds 1688, far from sites of active star formation. Parise et al. (2011) detected extended para-D₂H⁺ emission towards this core. Based on the analysis of the para-D₂H⁺ and ortho-H₂D⁺ lines towards its centre, Parise et al. suggested that the average density of the core is high compared to typical starless cores, a few times 10⁵ cm^-3. In accordance with the high abundances of the deuterated ions, the deuterium fraction in N₂H⁺ is also extremely high: N₂D⁺/ N₂H⁺ = 0.43 ± 0.11 (Punanova et al. 2016).

In Fig. 1, we show the dust colour temperature (T_C) and the H₂ column density (N(H₂)) maps of the core derived from far-infrared images observed by Herschel (Pilbratt et al. 2010). Contours of the 850 μm emission map from the SCUBA-2 survey of Pattle et al. (2015) are superposed on the N(H₂) map.

The Herschel/SPIRE images were extracted from the pipeline-reduced images of the Ophiuchus complex made in the course of the Herschel Gould Belt Survey (André et al. 2010). The data are downloaded from the Herschel Science Archive (HSA)¹. We calculated the T_C and N(H₂) distributions using only the three SPIRE (Griffin et al. 2010) bands at 250 μm, 350 μm and 500 μm, for which the pipeline reduction includes zero-level corrections based on comparison with the Planck satellite data. A modified blackbody function with the dust emissivity spectra index β = 2 was fitted to each pixel after smoothing the 250 μm and 350 μm images to the resolution of the 500 μm image (~ 40″), and resampling all images to the same grid. For the dust emissivity coefficient per unit mass of gas, we adopted the value from Hildebrand (1983), κ_{250 μm} = 0.1 cm² g^-1 (1 /C₂₅₀ in Table 1 in their paper). Suutarinen et al. (2013) derived a similar value for κ_{250 μm} in the starless core CrA C. According to the derived maps, the dust colour temperature minimum and the column density maximum of the core is found at RA ${\mathrm{16}}^{\mathrm{h}}{\mathrm{27}}^{\mathrm{m}}\mathrm{59}\stackrel{\mathrm{s}}{.}\mathrm{0}$, Dec. −24°33′33″ (J2000). The line observations presented here were done towards this position. The obtained minimum colour temperature is 12.4 K and the maximum column density is N(H₂) = 5.7 × 10²² cm^-2. These values are averages over the 40″ beam. The fitted T_C overestimates the mass-averaged dust temperature because of line-of-sight temperature variations. This effect is more marked towards the centre of a starless core than on the outskirts of the core (Nielbock et al. 2012; Suutarinen et al. 2013).

Our position lies approximately 13″ northeast from the centre position used by Parise et al. (2011), and approximately 7″ east of the 450 and 850 μm peaks observed with SCUBA-2. In a later section we use the SCUBA-2 maps to derive a simple spherically symmetric model of the core for the purpose of radiative transfer modelling.

2.2. GBT observations

The observations were carried out using the 7-beam K-Band Focal Plane Array (KFPA) at the GBT, with the Versatile GBT Astronomical Spectrometer (VEGAS) backend, as part of the Gould Belt Ammonia Survey (GBT15A-430, PIs: Friesen & Pineda). VEGAS was configured in Mode 20, which uses 8 separate spectral windows per KFPA beam, each with a bandwidth of 23.44 MHz and 4096 spectral channels. The spectral resolution is 5.7 kHz, corresponding to ~ 0.07 km s^-1. Observations were performed using in-band frequency switching with a frequency throw of 4.11 MHz. Here we use the NH₃(1,1) and (2,2) line maps of a 6′ × 6′ region centered on the column density peak of H-MM1. The integrated NH₃(1,1) intensity map of this region is shown in Fig. 2.

These observations were part of a much larger area map of the entire L1688 region during the 15A semester, which will be presented by Friesen & Pineda et al. (in prep.), and carried out on 10′ × 10′ boxes scanned in right ascension with rows separated by 13″ in declination. The scanning rate was 6.2″ s^-1, with a data dump every 1.044 s. A fast frequency switching rate of 0.348 s was used, which results in a root mean square noise (RMS) of 0.1 K (on the T_MB scale).

The data are calibrated using the GBT KFPA data reduction pipeline (Masters et al. 2011). The data were calibrated to the T_MB scale using the gain factors for each beam calibration derived from the Moon observations. The final cubes are created by a custom-made gridder using a tapered Bessel function for the convolution following Mangum et al. (2007)².

The observed ammonia transitions are indicated in the energy level diagram shown in Fig. 3.

2.3. APEX observations

The centre position of H-MM1 was observed using the upgraded version of the First Light APEX Submillimeter Heterodyne instrument (FLASH; Heyminck et al. 2006) on APEX (Güsten et al. 2006). This instrument, FLASH⁺ (Klein et al. 2014), operates in the 345 GHz and the 460 GHz atmospheric windows simultaneously, and can record two 4-GHz-wide sidebands separated by 12 GHz in both windows, that is, altogether 4 × 4 GHz frequency bands. The receivers were connected to MPIfR Fast Fourier Transform Spectrometers (XFTTS, Klein et al. 2012) with spectral resolutions of ~ 0.03 and ~ 0.05 km s^-1 at 345 and 460 GHz, respectively.

The sky subtraction was done by position switching, using an absolute reference position (RA 16^h28^m32^s, Dec. −24°31′00″, J2000) which, judging from Herschel far-infrared maps is void of dense gas. In the lower frequency window we used two frequency settings, which covered 1) the ground-state lines of para- and meta-ND₃ (hereafter pND₃ and mND₃) at 306.7 and 309.9 GHz, and 2) the ground-state lines of ortho- and para-NH₂D (oNH₂D and pNH₂D) at 332.8 GHz, and the ground-state lines of ortho- and para-NHD₂ (oNHD₂ and pNHD₂) at 335.5 GHz. The first tuning also covered the N₂D⁺(4−3) line at 308.4 GHz. The observations in the 460 GHz window were aimed at the N₂D⁺(6−5) and N₂H⁺(5−4) lines. Because these two lines could be measured simultaneously, the tuning of the 460 GHz receiver was kept constant during the whole observing run.

A list of transitions covered is given in Table 1, mentioning only the most significant to the present study. Here we give the centre frequencies of transitions, upper state energies and Einstein coefficients for spontaneous emission. These parameters are obtained from the Cologne Database for Molecular Spectroscopy, CDMS³. The observed transitions of NH₃, NH₂D, NHD₂ and ND₃ are indicated in the energy level diagrams in Figs. 3.

The APEX beamsize (Full width at half maximum, FWHM) is ~ 20″ at 310−330 GHz, and ~ 14″ at 465 GHz. The main beam efficiency, η_MB, is 0.73 at 310−330 GHz, and 0.6 at 465 GHz (Güsten et al. 2006). The observations were carried out between the 29th and 31st May, 2015. The total observing time was 11.4 h. The weather conditions were stable and fairly good (PWV 0.7−1.2 mm). The absolute calibration, pointing and focus were checked by observing Saturn. The system temperatures were in the following ranges: 180−190 K (310 GHz), 230−250 K (335 GHz) and 570−600 K (465 GHz). The resulting RMS noise levels at the mentioned frequencies at a velocity resolution of 0.1 km s^-1 were 18, 27 and 33 mK, respectively, on the T_MB scale.

2.4. IRAM observations

The column density peak of H-MM1 was observed with the IRAM 30-m telescope on the 5th July, 2015, in acceptable weather conditions (PWV 8−10 mm). Pointing and focus were checked towards QSO 1253-055. The following transitions were observed: oNH₂D(1₁₁−1₀₁) at 85.9 GHz, pNH₂D(1₁₁−1₀₁) at 110.2 GHz and N₂D⁺(2−1) at 154.2 GHz. The measurements were obtained with the EMIR 090 and 150 receivers⁴ and the VESPA auto-correlator with a spectral resolution of 20 kHz; the corresponding velocity resolutions were 0.04−0.07 km s^-1. The beam sizes were 29″, 23″ and 16″ for oNH₂D(1₁₁−1₀₁), pNH₂D(1₁₁−1₀₁) and N₂H⁺(2−1), respectively. The system temperatures ranged from 166 to 310 K depending on the frequency (see Table 2 for the details). The spectra were taken in the position switching mode, using an off position 400′′ East of the target. The integration time for each line was between 15 and 22 min, which resulted in RMS noise levels of 0.08−0.14 K on the T_MB scale. The intensity scale was converted to the main-beam temperature scale using the beam efficiency values given in the IRAM 30-m report on a beam pattern (Kramer et al. 2013); see Table 2 for details.

3. LTE Analysis of the observed spectra

In this section, we present the observed spectra and the results of the standard hyperfine structure fitting to the detected lines. The method assumes line-of-sight homogeneity and that the populations of hyperfine states in a certain rotational manifold are proportional to their statistical weights, according to the assumption of local thermodynamic equilibrium (LTE). The total column density estimates rely, furthermore, on the assumption that the excitation temperature, T_ex, is constant for all rotational transitions of a molecule.

The observed spectra are shown in Figs. 4 (NH₃), 5 (NH₂D, NHD₂ and ND₃) and 6 (N₂D⁺ and N₂H⁺). The spectra were reduced using the GILDAS software package⁵. All the observed transitions have hyperfine structure. The detected lines were analysed using the HFS method implemented in the CLASS software (part of GILDAS) or fitting routines written in the IDL language.

The oNH₃(3,3) and N₂D⁺(6−5) lines were not detected, while N₂H⁺(5−4) shows perhaps a weak line with an integrated intensity of ~ 30 mK km s^-1. The upper limits for the intensities of these lines are 0.1 K, 0.05 K and 0.1 K, respectively, on the T_MB scale. For the N₂D⁺(6−5) and N₂H⁺(5−4) spectra we have only used scans (20 s integrations) with smooth baselines, showing no visible disturbances, which means approximately half of the measurements.

The hyperfine patterns of the detected lines are dominated by the splitting caused by the electric quadrupole moment of the ¹⁴N nucleus. In the NH₃ inversion lines, the structure owing to the magnetic moments of the N and H nuclei can be partially resolved (Kukolich 1967; Ho & Townes 1983). Furthermore, for the NH₂D lines we use the hyperfine patterns calculated by Daniel et al. (2016a), including the effects of both N and D nuclei. As shown recently by Daniel et al. (2016a), quadrupole coupling of the D nucleus substantially broadens the observed 1₁₁−1₀₁ lines of NH₂D at 86 and 110 GHz, and disregarding this effect results in overestimates of the line widths by ~ 50% for a cold cloud. For the ground-state lines at 333 GHz, and for any lines at higher frequencies, the effect of the hyperfine structure owing to D is, however, negligible compared with the Doppler broadening.

For other molecules observed here, we use line lists that only take the splitting due to N into account. The data are from Coudert & Roueff (2006) (NHD₂ and ND₃), and from Caselli et al. (1995), Dore et al. (2004) and Pagani et al. (2009a) (N₂H⁺ and N₂D⁺). The effect of other interactions is likely to be small compared with the thermal broadening (Kukolich 1969; Dore et al. 2004). In the hyperfine fitting, we have assumed that the frequencies listed in Table 1, which are adopted from the CDMS, represent weighted averages of the hyperfine components.

The results of Gaussian fits to the hyperfine structure are presented in Table 3. This table contains the peak main-beam brightness temperatures, T_MB, the radial velocities, V_LSR, the line widths, Δυ, and the sums of the peak optical thicknesses, τ_sum, of the hyperfine components. The hyperfine fit also gives an estimate for the T_ex of the transition, provided that the spectrum is on the brightness temperature scale. We have approximated this by the T_MB scale, which is equivalent to assuming that the source fills the telescope beam uniformly. Finally, the two last columns of Table 3 give the total column densities, N_total, of the molecules, and the fractional abundances, X, using the H₂ column density derived from Herschel observations. The derived total optical thicknesses of the detected lines range from ~ 1 to ~ 5. This means that the satellites are mainly optically thin, and we are in the regime where the integrated intensity is proportional to the column density of the molecule.

While the accuracy of the LSR velocities and line widths from the hyperfine fits is high, the optical thicknesses, and consequently the column densities derived using this method have relatively large uncertainties, except for pNH₃ and oNH₂D, which are the brightest lines. The oNH₂D and pNH₂D column densities derived from the lines observed with APEX and IRAM are consistent when the different beam sizes (see Table 2) and the uncertainties owing to noise are taken into account. In contrast, the N₂D⁺ column densities derived from the N₂D⁺(2−1) (IRAM) and N₂D⁺(4−3) (APEX) spectra differ by a factor of five. We consider the value obtained from the N₂D⁺(2−1) spectrum more reliable because this line has several resolved hyperfine components, and because this transition connects rotational levels that are much more densely populated than J = 3 and J = 4. Furthermore, the N₂D⁺ column density obtained from N₂D⁺(2−1) is consistent with that derived from N₂D⁺(1−0) by Punanova et al. (2016) (3.8 ± 0.9 × 10¹² cm^-2, beamsize 32″) towards a position lying 13″ away from our centre position.

Assuming that ortho- and para-NH₃ have equal abundances (which would correspond to their nuclear spin statistical weights), we obtain a total ammonia column density of N(NH₃) = (3.3 ± 0.4) × 10¹⁴ cm^-2. The total column densities for NH₂D and NHD₂ are N(NH₂D) = (1.5 ± 0.1) × 10¹⁴ cm^-2, N(NHD₂) ~ 6.7 × 10¹³ cm^-2. As in the case of NH₃, we only detected one spin modification of ND₃. Assuming that the relative abundances of ortho-, meta- and para-ND₃ correspond to their nuclear spin statistical weights, 16:10:1, we get an estimate for total ND₃ column density with large error margins: N(ND₃) = (2.7 ± 1.4) × 10¹² cm^-2. These column density estimates imply the following fractionation ratios: NH₂D/NH₃ = 0.45 ± 0.09, NHD₂/ NH₂D ~ 0.45 and ND₃/ NHD₂ ~ 0.04. The spin ratios can only be estimated for NH₂D and NHD₂ for which we get o/pNH₂D = 3.0 ± 1.1 and o/pNHD₂ ~ 1.5. We see that both ortho/para ratios are close to their statistical values, 3 and 2, respectively.

We used the NH₃(1,1) and (2,2) maps from GBT to derive the kinetic temperature, T_kin, and the pNH₃ column density, N(pNH₃), distributions in the vicinity of the H-MM1. The standard analysis described in Ho & Townes (1983), Walmsley & Ungerechts (1983) and Ungerechts et al. (1986) was used. According to the modelling results of Juvela et al. (2012), the ammonia spectra faithfully trace the real mass-averaged gas temperature. In Fig. 7, we show the average N(pNH₃) and T_kin as functions of distance from the core centre. The diagrams are derived by averaging spectra over concentric, 10″ wide annuli, and calculating the parameters and their errors from these averaged spectra. Also shown in this figure are the corresponding distributions of N(H₂) and T_C derived from the Herschel/SPIRE maps. Two things are perhaps worthy of notice in the diagrams: the NH₃ abundance seems to decrease and the T_kin seems to increase towards the edge of the core, although both quantities have large errors away from the centre of the core.

The line widths of the spectra observed towards the core centre range from 170 to 260 m s^-1. Assuming that the average kinetic temperature in the core is ~ 11 K, the non-thermal velocity dispersions obtained from various transitions are between ~ 50 and ~ 100 m s^-1.

4. Abundances from radiative transfer modelling

4.1. Physical model of H-MM1

In order to account for inhomogeneities in the density and temperature distributions along the line of sight, we constructed a spherically symmetric physical model of the core. Besides providing a realistic description of molecular line excitation conditions for the radiative transfer modelling, the core model is also needed for making a connection between the observations and the theory of interstellar chemistry. In this model, the core structure is described by a modified Bonnor-Ebert sphere (MBES; Evans et al. 2001; Zucconi et al. 2001; Sipilä et al. 2011, 2015c), which is a pressure-bound, hydrostatic sphere of gas and dust with temperature decreasing towards the centre.

We fixed the outer radius of the core to be R_out = 9600 AU (80″). At this distance, the core was assumed to be merged with the ambient cloud. We assumed a constant non-thermal velocity dispersion of σ_NT = 100 m s^-1 inside the core, which increases the internal pressure slightly. The dust temperature profile was calculated using a Monte Carlo program for continuum radiative transfer, CRT⁶, developed by M. Juvela (Juvela 2005). The spectrum of the unattenuated interstellar radiation field (ISRF) was taken from Black (1994). We used the dust opacity data from Ossenkopf & Henning (1994) for unprocessed dust grains with thin ice coatings⁷, which agree with the opacities at 250, 350 and 500 μm used in the derivation of the T_C and N(H₂) maps in Sect. 2.

The MBES model was constructed using the following constraints: 1) the dust temperature at the boundary should agree with the T_C outside the core derived from Herschel; 2) the model should approximately reproduce the 450 and 850 μm emission profiles derived from SCUBA-2 maps of Pattle et al. (2015) and 3) the mass-averaged gas kinetic temperature profile, smoothed to the angular resolution of the GBT, should agree with the observed T_kin profile shown in Fig. 7. In order to achieve the large temperature difference between the edge (~ 15 K) and the centre (< 11 K) with the model where the external heating is dominated by the ISRF, we set the visual extinction to A_V = 2^mag at the outer boundary of the core. With this choice, we assume that the core lies near the edge or in a protrusion of the ambient cloud. The total hydrogen column density in the neighbourhood of the core is of the order of 10²² cm^-2 (Fig. 1), which corresponds to A_V ~ 10^mag. The strongly peaked sub-millimetre emission observed with SCUBA-2 could only be reproduced with central H₂ densities of the order of 10⁶ cm^-3.

The iteration was started from the density distribution of an isothermal Bonnor-Ebert sphere at 11 K with a central density of n(H₂) = 10⁶ cm^-3. The intensity of the ISRF and the central density were adjusted to reach agreement with the constraints 1) and 2) above. The sub-millimetre intensity profiles of the model core were calculated by evaluating the integrals ${{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{L}}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathit{,z}\mathrm{\right)}\mathit{\rho }\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{\kappa }}_{\mathit{\nu }}\mathrm{d}\mathit{z}$, where the integration is along the line of sight, at different distances from the core centre. The resulting intensity profiles at 450 μm and 850 μm, smoothed to appropriate resolutions, were compared with corresponding circularly averaged profiles from the SCUBA-2 maps.

The gas and dust temperatures are assumed to be equal at densities above ~ 10⁵ cm^-3 (Goldsmith 2001). The observed T_kin profile shown in Fig. 7 suggests, however, that the (3-dimensional) gas temperature distribution inside the core does not show the steep gradient characteristic of T_dust caused by the attenuation of the ISRF. In fact, the ammonia observations can be explained with a model where the core is mostly isothermal, but the temperature rises steeply near the edge. We note that this conclusion is probably influenced by the limited angular resolution (cf. temperature determination in L1544 by Crapsi et al. 2007). To conform with the prediction that T_gas ~ T_dust at the highest densities, we assumed that this is true above a certain, adjustable density threshold, but that below this threshold, T_gas is constant up to the transition layer, where it rises abruptly. The adopted gas temperature distribution required a slight adjustment of the density profile to keep the core in hydrostatic equilibrium. The iteration converged after three to four rounds. A reasonable agreement with the dust continuum observations was found with a model where the central density of the core is n(H₂) = 1.2 × 10⁶ cm^-3 and the standard IRSF is scaled up by the factor 1.7. The observed T_kin distribution could be reproduced assuming that T_gas separates from T_dust at the density n(H₂) ~ 4 × 10⁵ cm^-3. At this point, ~ 15″ from the centre, the temperature is ~ 11 K.

The distributions of T_gas and T_dust as functions of the radial distance from the core centre for the best fit model are shown in Fig. 8, together with the observable mass-averaged temperature profiles smoothed to the angular resolutions of GBT and Herschel. Figure 9 shows the circularly averaged 450 μm and 850 μm intensity profiles of H-MM1 derived from the SCUBA-2 maps of Pattle et al. (2015), together with predictions from the MBES model.

4.2. Average abundances from line modelling

We derived the fractional abundances of the observed molecules in H-MM1 by applying the Monte Carlo radiative transfer program of Juvela (1997) to the physical model described in Sect. 4.1. The collisional rate coefficients for NH₂D were adopted from Daniel et al. (2014), and for NHD₂ and ND₃ we used the newly calculated coefficients from Daniel et al. (2016b). For N₂D⁺, we used the collisional rate coefficients for N₂H⁺ from the recent work of Lique et al. (2015).

In this calculation, we assumed constant fractional abundances throughout the core, that is, that there is no dependence on the distance from the core centre. The obtained values can thus be taken as averages over the core. Starting from the estimates presented in Table 3 of Sect. 3, we varied fractional abundances until the modelled spectra produced the same integrated intensities as the observed ones. The best-fit fractional abundances are presented in Table 4. The uncertainties correspond to the 1σ errors of the integrated intensities. The para-NH₃ lines are broader than those of its deuterated isotopologues, and for this molecule we had to increase the assumed non-thermal velocity dispersion to σ_NT = 150 m s^-1 in order to reproduce the integrated intensities. This suggests that NH₃ emission has contribution from the ambient cloud having a larger velocity dispersion than the core, and that the derived fractional para-NH₃ abundance is an upper limit for the core.

For the different isotopologues of ammonia, the agreement reached between the predicted spectra and observations is equally good as for the hyperfine fits shown in Figs. 4 and 5. The two N₂D⁺ lines detected cannot be reproduced by a single abundance. The value listed in Table 4 is a compromise that over-predicts the N₂D⁺(4−3) intensity but gives an overly weak N₂D⁺(2−1) line. The situation is thus opposite to what is expected from the results of the hyperfine fits. This indicates that either the assumption of a constant abundance is unrealistic for N₂D⁺ or that the physical model is inaccurate.

The abundances listed in Table 4 imply the following total fractional abundances: X(NH₃) = (7.6 ± 0.1) × 10^-9 (assuming o:p=1:1), X(NH₂D) = (2.9 ± 0.1) × 10^-9, X(NHD₂) = (6.4 ± 0.5) × 10^-10 and X(ND₃) = (3.7 ± 0.4) × 10^-11 (assuming o:m:p = 18:10:1). Here we have assumed statistical spin ratios for species without a line detection. The corresponding deuterium fractionation ratios for ammonia, and the ortho/para ratios for NH₂D and NHD₂ are listed in the bottom part of Table 4.

For the molecules with the brightest lines, the abundances from the radiative transfer modelling agree relatively well with those from the LTE analysis (Table 3). The most glaring discrepancy (by a factor of 2.5) is found for para-NHD₂ with the weakest detection and a very large uncertainty of the optical thickness from the LTE method. The statistical errors of the fractional abundances listed in Table 4 are small in most cases, but the values are subject to systematic errors depending on accuracy of the physical model, and on the validity of the assumption of constant abundances. The present physical model is, however, consistent with dust continuum observations, and provides a more realistic description of the excitation conditions in the core than the assumption of line-of-sight homogeneity. Therefore, we consider the abundance ratios listed in Table 4 to be more accurate than those implied by the values presented in Table 3, and use the former ratios to assess the validity of the chemistry model described below.

5. Chemical modelling

5.1. Model description

We model the chemistry of H-MM1 using the pseudo-time-dependent gas-grain chemical code presented in earlier papers (Sipilä 2012; Sipilä et al. 2013, 2015a,b) where the details of the code (e.g., the expressions of the various reaction rate coefficients) can be found. The model includes gas-phase chemistry, adsorption onto and (non-thermal) desorption from grain surfaces, and grain-surface chemistry. Tunnelling diffusion of H and D atoms on grains is not considered in the present calculations, whereas tunnelling through activation energy barriers in surface reactions has been included.

The program can be instructed to include various desorption mechanisms: thermal desorption (negligible in the physical conditions explored here), cosmic ray desorption, reactive desorption and photodesorption. Cosmic ray desorption is treated following Hasegawa & Herbst (1993). Exothermic association reactions on the surface can result, when this option is turned on, in desorption of the reaction product with an efficiency of 1% (Garrod et al. 2007). Finally, it is possible to include photodesorption of water, CO and ammonia caused by secondary UV photons created by H₂ excitation (Prasad & Tarafdar 1983). In the simulations presented here, however, the photodesorption and reactive desorption options have been turned off. According to extensive testing, the inclusion of these processes increases the ammonia production, which in turn makes it necessary to lower the elemental N abundance to ensure compliance with the observed line intensities, but fractionation and spin ratios remain largely unchanged.

The spin-state chemical model presented in Sipilä et al. (2015a) describes the spin states of species involving multiple protons. Recently, we upgraded this model to include a self-consistent description of the spin states of multiply-deuterated species (Sipilä et al. 2015b). This is achieved by considering nuclear spin selection rules arising from molecular symmetries, assuming full scrambling of nuclei in reactive collisions. The model gives the necessary information for the present application, that is, the calculation of simulated line emission from deuterated ammonia. In the beginning of the simulation, all elements are in the atomic form, with the exceptions of hydrogen and deuterium, which are initially locked in H₂ and HD, respectively.

The descriptions of deuterium and spin-state chemistry adopted here apply to both gas-phase chemistry and grain-surface chemistry. However, the formation mechanism of ammonia (as a particular example) is different in the gas and on the grains. In the gas phase, (deuterated) ammonia forms through a network of ion-molecule reactions, while on the surface, the formation mechanism is hydrogen/deuterium addition. Therefore we expect non-statistical deuterium and spin-state ratios in the gas phase, and statistical ratios on the grain surface, although complete scrambling is assumed in both cases. The assumed binding energies on grain surfaces are the same as listed in Table 2 of Sipilä et al. (2015a).

For the elemental abundances, we have adopted the set of low-metal elemental abundances labelled EA1 in Wakelam & Herbst (2008), except for the nitrogen abundance for which we needed to increase the EA1 abundance by a factor of 2.5 to N/H = 5.3 × 10^-5 to reproduce the observed line intensities of NH₂D, NHD₂, ND₃ and N₂D⁺. According to the recent compilation of Jenkins (2009), the appropriate value for the diffuse ISM is N/H = 6.2 × 10^-5. The adopted carbon and oxygen abundances are O/H = 1.8 × 10^-4 and C/H = 7.3 × 10^-5.

We derive abundance profiles for the various molecules by separating the core model (see Sect. 4.1) into a series of concentric spherical shells, each associated with unique values of density, gas/dust temperatures and visual extinction A_V. Chemical evolution is then calculated separately in each shell, which leads to simulated abundances for each chemical species as functions of time and radial distance from the core centre. The integration of the modelled abundance gradients into the radiative transfer model is discussed in Sect. 5.3.

5.2. Chemical evolution of the core

We calculate the evolution of chemical abundances in the core model with the aim being to examine if the model can reproduce, at a certain stage of the simulation, the intensities of the lines of deuterated ammonia and N₂D⁺ observed towards H-MM1 in the present study, as well as the intensities of the oH₂D⁺ and pD₂H⁺ lines observed previously by Parise et al. (2011).

We assume that the spin temperature of H₂ has been thermalised during the initial contraction phase of the cloud down to ~ 20 K (Flower et al. 2006b; Sipilä et al. 2013), and accordingly, set the initial o/pH₂ to 1 × 10^-3. For the cosmic ray ionisation rate of H₂, we assume ζ_H2 = 1.3 × 10^-17 s^-1, and the average grain radius is set to a = 0.1 μm. These are kept unchanged in the present simulations. We note, however, that deuteration can be delayed by increasing the initial o/pH₂ ratio, and that the fractionation ratios can be lowered by decreasing the average grain size or by increasing the cosmic ionisation rate (Sipilä et al. 2010). On the other hand, an increase of the cosmic ray ionisation rate would generally increase the abundances and the line intensities of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, NH₃ and their deuterated isotopologues, and decrease the abundances and line intensities of N₂H⁺ and N₂D⁺.

We first discuss the predictions for some of the most common species. The gas-phase abundances of the H, D and N atoms, and the HD, CO and N₂ molecules, relative to the total hydrogen abundance, are plotted in Fig. 10, as functions of time for our fiducial core model discussed in Sect. 4.1. The abundances are averages over the line of sight through the centre of the core, weighted by the density. The freeze-out of CO is followed by an increase of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, which in turn results in an enhanced ortho-para conversion of H₂. At the same time, deuterium is efficiently transferred from HD to deuterated ions in the gas phase (to H₂D⁺, D₂H⁺ and ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ in the first place). The H and D atoms released in the dissociative recombination of deuterated ions mainly accrete onto grains, where they can combine to give back H₂ or HD, but also react with heavier atoms or radicals. At late stages of chemical evolution, deuterium becomes increasingly incorporated into icy compounds. In the gas phase, this is reflected by the reduction of the HD abundance.

The rapid decrease of atomic nitrogen at the beginning of the simulation is mainly caused by accretion onto grains. A fraction of the nitrogen atoms in the gas phase is converted to N₂ through N + OH → NO + H, N + NO → N₂ + O (Flower et al. 2006a; Le Gal et al. 2014). This reaction is important at early times, when the N and OH abundances are high. Also, the N₂ molecules accrete onto grains, but their desorption from grains is more significant than for N atoms, which quickly react with other atoms or radicals, for example, N^∗ + H^∗ → NH^∗, N^∗ + O^∗ → NO^∗ or ${\mathrm{N}}^{\mathrm{\ast }}\mathrm{+}{\mathrm{N}}^{\mathrm{\ast }}\mathrm{\to }{\mathrm{N}}_{\mathrm{2}}^{\mathrm{\ast }}$. Species attached to grains are indicated here with asterisks. The surface species ${\mathrm{N}}_{\mathrm{2}}^{\mathrm{\ast }}$ is destroyed by only two processes; desorption or photodissociation, ${\mathrm{N}}_{\mathrm{2}}^{\mathrm{\ast }}\mathrm{+}\hspace{0.17em}\mathrm{photon}\hspace{0.17em}\mathrm{\to }{\mathrm{N}}^{\mathrm{\ast }}\mathrm{+}{\mathrm{N}}^{\mathrm{\ast }}$, by cosmic-ray induced UV photons. This is the main difference from CO^∗ for which hydrogenation, CO^∗ + H^∗ → HCO^∗, competes strongly against desorption (see also Sect. 6.1).

Owing to desorption, the N₂ abundance remains high in the gas phase until late stages of the simulation. Molecular nitrogen is prerequisite to N₂H⁺, which forms through ${\mathrm{N}}_{\mathrm{2}}\mathrm{+}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\to }{\mathrm{N}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}$. In case reactive desorption is included, the most important source of ammonia at very early stages of simulation is formation on grains, $\mathrm{N}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{\ast }}\mathrm{+}{\mathrm{H}}^{\mathrm{\ast }}\mathrm{\to }\mathrm{N}{\mathrm{H}}_{\mathrm{3}}$, where ammonia is supposed to be released into the gas phase at the probability of 1%. After a few thousand years, gas-phase formation takes over. In the simulations presented here, the ammonia production is always dominated by the well-known chain of gas-phase reactions, terminating in $\mathrm{N}{\mathrm{H}}_{\mathrm{4}}^{\mathrm{+}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}}\mathrm{\to }\mathrm{N}{\mathrm{H}}_{\mathrm{3}}\mathrm{+}\mathrm{H}$ (e.g. Le Gal et al. 2014; Roueff et al. 2015). As discussed by Sipilä et al. (2015b), the dissociative ionisation of HNC by He⁺ helps the initiation of this chain by producing NH⁺, which would otherwise be solely dependent on N⁺ + oH₂ → NH⁺ + H (Dislaire et al. 2012).

The evolution of the abundances of NH₃ and N₂H⁺ as well as their deuterated isotopologues in the gas phase are shown in Fig. 11 (top panel). The other two panels show the fractionation ratios and the spin ratios for these species. In Fig. 12, we show the abundances, fractionation ratios and the spin ratios of the grain-surface species $\mathrm{N}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{\ast }}$, NH₂D^∗, $\mathrm{NH}{\mathrm{D}}_{\mathrm{2}}^{\mathrm{\ast }}$ and $\mathrm{N}{\mathrm{D}}_{\mathrm{3}}^{\mathrm{\ast }}$. The middle panel of this latter figure also shows the atomic D^∗/H^∗ on grains.

In the present model, the gas-phase NH₃ abundance is built up early and does not change significantly at later times. The largest variations are seen in the abundances of NHD₂, ND₃ and N₂D⁺, which grow rapidly in the beginning, and decay slowly after the deuteration peak. This behaviour seems to reflect the variations in the D₂H⁺ and ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ abundances that are shown in Fig. 13. The ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ ion reaches its maximum before D₂H⁺, which in turn peaks before H₂D⁺. Likewise, the maximum fractionation ratio ND₃/ NHD₂ occurs earlier than the maximum in the NHD₂/ NH₂D ratio, which again takes place long before the NH₂D/NH₃ peak.

In the bottom panel of Fig. 11, one can see that while the o/p-NH₂D ratio decreases slightly with time, o/p-NHD₂, m/p-ND₃ and m/o-ND₃ have increasing tendencies. The predicted ratios are close to their statistical values in the beginning of the simulation. The o/p-NHD₂ and m/p-ND₃ ratios mimic the corresponding ratios of D₂H⁺ and ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ shown in Fig. 13, bottom panel. The relationship between the spin modifications of NH₃ and ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ is discussed in Sect. 6.

The ammonia production on grain surfaces is very efficient in the present model. The deuteration of ammonia occurs more slowly than in the gas phase, and never reaches fractionation ratios as large as those seen there. The spin ratios on grains stay close to their statistical values at all times, except that m$\mathrm{N}{\mathrm{D}}_{\mathrm{3}}^{\mathrm{\ast }}$ is enhanced at the cost of o$\mathrm{N}{\mathrm{D}}_{\mathrm{3}}^{\mathrm{\ast }}$ and p$\mathrm{N}{\mathrm{D}}_{\mathrm{3}}^{\mathrm{\ast }}$.

5.3. Predicted spectra

The radial distributions of the density, temperature and the (time-varying) chemical abundances in the gas phase are used as input for a Monte Carlo radiative transfer program (Juvela 1997) to predict observable rotational line profiles. As in simulations described in Sect. 4.2, we use a larger non-thermal velocity dispersion (σ_N.T. = 150 ms^-1) for NH₃ than for the deuterated species (for which σ_N.T. = 100 ms^-1) in order to reach agreement with the observed integrated intensities.

The line intensities depend, besides the (mass averaged) abundances of the species within the telecope beam, on their radial distributions, which change with time. At early stages, deuterated species are concentrated on the core centre with high densities, whereas later on, when freezing onto grains reduces the abundances in the centre, line emission is dominated by lower-density outer parts of the core. The changes of the abundance profiles mean that at early times, lines with large transition dipole moments, such as N₂D⁺(4−3) and mND₃(1₀−0₀), are much stronger than the o/pNH₂D(1₁₁−1₀₁) lines, for example; the reverse is true at later times, however. The fractional abundances of selected species as functions of the radius are shown in Fig. 14. The distributions are taken at the time when most of the predicted spectra agree reasonably well with the observations.

The intensities of the simulated NH₂D, NHD₂, ND₃ and N₂D⁺ spectra are comparable with those of the observed ones during rather a short period around the NHD₂/ NH₂D peak, occurring at ~ 3 × 10⁵ yr from beginning of the simulation. At this time, however, the predicted beam-averaged pNH₃ abundance, X(pNH₃) ~ 5 × 10^-9, is approximately 30% higher than needed to reproduce the NH₃(1,1) and (2,2) line intensites observed at the GBT (see Table 4). The simulated NH₃ lines agree with the observed ones at an early stage, around 10⁵ yr, and much later, around 2 × 10⁶ yr, when ammonia depletion has finally started to take effect. The predicted spectra at t = 3 × 10⁵ yr are shown in Figs. 15 (NH₃), 16 (NH₂D, NHD₂ and ND₃), 17 (N₂D⁺) and 18 (oH₂D⁺ and pD₂H⁺).

The model predicts that the abundances of the different isotopologues of NH₃ and N₂H⁺ grow at different rates. This implies that the fractionation ratios, NH₂D/NH₃, NHD₂/ NH₂D, ND₃/ NHD₂, do not necessarily reach their maxima at the same time. This is illustrated in the middle panel of Fig. 11. In this model, the peak fractionation ratios range from 0.15 to 0.30, depending on the pair of species considered. The ND₃/ NHD₂ ratio mimics the N₂D⁺/ N₂H⁺ ratio divided by 10, and these two fractionation ratios are the first to peak in all the models we have run. This tendency is related to the rapid growth of the ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ abundance occurring at early stages of the simulation (Fig. 13). The fact that the average temperature exceeds 10 K has a favourable effect on deuteration, but its fast advancement is made possible by the low initial o/pH₂ ratio assumed in the simulation. A higher abundance of oH₂ would delay the deuterium peak by obstructing the primary deuteration through reaction ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{HD}\mathrm{\leftrightarrow }{\mathrm{H}}_{\mathrm{2}}{\mathrm{D}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}$ (Flower et al. 2006b; Pagani et al. 2011, 2013; Kong et al. 2015).

In Figs. 18 and 19, we compare our model predictions with the previous observations of Parise et al. (2011) and Punanova et al. (2016). The modelled oH₂D⁺(1₁₀−1₁₁) and pD₂H⁺(1₁₀−1₀₁) spectra shown in Fig. 18 are on the ${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}$ scale as observed with APEX to allow comparison with the spectra shown in Fig. 3 of Parise et al. (2011). These spectra are reproduced in Fig. 18. The 13″ offset from the supposed core centre position is taken into account. The hyperfine patterns of the lines have been adopted from Jensen et al. (1997). The line profiles are dominated, however, by the thermal broadening. While the simulated pD₂H⁺ line approximately agrees with the observations, the simulated oH₂D⁺ line is brighter than the observed one by a factor of three. A good agreement with both observations would be found at very early times by setting the initial o/pH₂ ratio to 10^-4. On the other hand, that model cannot reproduce the observed line ratios for the other molecules. After t ~ 4 × 10⁵ yr, the pD₂H⁺ line intensity decreases rapidly below the observed ${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}\mathrm{~}\mathrm{0.1}$ K while oH₂D⁺ remains relatively strong. This behaviour is determined by the close correlation between o/pH₂D⁺ and o/pH₂, and the increase of the o/pD₂H⁺ ratio with time (Fig. 13).

The N₂H⁺(1−0), N₂D⁺(1−0), and C¹⁷O(1−0) spectra observed at the IRAM 30-m telescope by Punanova et al. (2016) are shown in Fig. 19 together with the predicted spectra at the time 3 × 10⁵ yr. The H-MM1 spectra of Punanova et al. (2016) were obtained towards the (0, 0) of Parise et al. (2011), and also here the 13″ offset from the core centre is taken into account. The model reproduces the observed N₂D⁺(1−0) spectrum reasonably well, but gives a vastly undervalued C¹⁷O(1−0) intensity. According to the model, CO is heavily depleted in the core at the time 3 × 10⁵ yr (Fig. 10). The peak velocity and velocity dispersion of the observed C¹⁷O(1−0) line are, however, different from those of the N₂H⁺ and N₂D⁺ lines. The C¹⁷O emission is probably dominated by the ambient cloud, which is not included in our model. Also the shape of the N₂H⁺(1−0) line suggests that part of the emission originates in the ambient cloud. According to the radial abundance distributions shown in Fig. 14, N₂H⁺ belongs to species that are not confined to the core. Nevertheless, the model seems to under-predict the N₂H⁺(1−0) emission from the core.

To summarise comparisons with observations, the core model predicts the observed intensities of the N₂D⁺, NH₂D, NHD₂, ND₃ and para-D₂H⁺ lines rather well, overproduces para-NH₃ and ortho-H₂D⁺ and gives overly low N₂H⁺(1−0) and C¹⁷O(1−0) line intensities. Besides these discrepancies, the present model has problems in reproducing the observed spin ratios. Inspection of Fig. 16 reveals that the model underestimates the o/pNH₂D ratio and overestimates the o/pNHD₂ ratio.

6. Discussion

6.1. Ubiquity of ammonia

Our chemical network predicts that the ammonia abundance builds up quickly in the gas phase. The predicted fractional abundances are similar to those found previously in molecular clouds, in particular, in the Ophiuchus complex (Friesen et al. 2009), even though the model over-predicts the pNH₃ abundance in H-MM1 during the deuteration peak. During the very beginning of the simulation, the production is dominated by surface reactions followed by desorption, but after a few thousand years, ion-molecule reactions in the gas phase take over as the main source of gaseous ammonia. At early stages, the rapidly evolving carbon chemistry comes to aid as the dissociative ionisation reaction HNC + He⁺ provides a short cut to NH⁺, past the slow formation of N₂ and the famous bottle-neck reaction N⁺ + H₂.

The ammonia abundance remains high until very late stages of the simulation. One of the reasons is the slow decrease of the molecular nitrogen abundance, which is sustained by desorption. In this respect, N₂ acts differently from CO, which freezes out quickly. In the presence of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and H⁺, N₂ replenishes the gas with the NH⁺ ion through the sequence ${\mathrm{N}}_{\mathrm{2}}\stackrel{{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}}{\mathrm{\to }}{\mathrm{N}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\stackrel{{\mathrm{e}}^{\mathrm{-}}}{\mathrm{\to }}\mathrm{NH}\stackrel{{\mathrm{H}}^{\mathrm{+}}}{\mathrm{\to }}\mathrm{N}{\mathrm{H}}^{\mathrm{+}}$. NH⁺, in turn, fuels the ammonia production through successive reactions with H₂ (see, e.g. Flower et al. 2006a; Le Gal et al. 2014).

In the present model, N₂ attached to a grain can either photodissociate or be desorbed. N₂ and CO, even if isoelectronic, have very different surface chemistries. The two molecules have approximately the same adsorption energies on amorphous ice: approximately 1100 K for CO and approximately 1000 K for N₂ (Hama & Watanabe 2013; Fayolle et al. 2016). Hydrogenation of CO (CO^∗ → HCO^∗ → H₂CO^∗ → H₃CO^∗ → H₃COH^∗) has been experimentally observed and characterised (Hama & Watanabe 2013, and references therein) at temperatures down to 3 K, where the intermediate HCO could be detected (Pirim & Krim 2011). The reaction CO + H → HCO proceeds by tunnelling through an activation barrier computed to be approximately 2000 K in the gas phase, and probably lower on the surface of amorphous ice (Peters et al. 2013; Rimola et al. 2014). The exothermicity of the reaction is approximately 6700 K, while the endothermicity of the reaction CO + H → COH is approximately 10 000 K (Zanchet et al. 2007; note that this theoretical paper overestimates the exothermicity of the reaction leading to HCO). The case of N₂ is different: the reaction N₂ + H → N₂H is endothermic by approximately 4400 K (Bozkaya et al. 2010) and does not occur on ice surfaces. The neutral chemistry of CO and N₂ proceeds thus in very different ways, with no surface hydrogenation of N$\begin{array}{c}\mathrm{\ast }\\ \mathrm{2}\end{array}$ towards ammonia NH₃ or hydrazine N₂H₄. However, the ionic chemistry on low-temperature ice surfaces is not fully characterised.

6.2. Fractionation ratios

The abundance of NH₂D starts to increase gradually, first through the deuteron transfer to ammonia, primarily by HCND⁺ or DCNH⁺, e.g. HCND⁺ + NH₃ → NH₃D⁺ + HCN, followed by dissociative recombination of NH₃D⁺. The depletion of CO boosts the abundance of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, which in turn is efficiently deuterated to H₂D⁺, D₂H⁺ and ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ in successive reactions with HD. This stage is characterised by a rapid increase of NHD₂, ND₃ and N₂D⁺. The most important reactions contributing to the formation of NH₂D during the deuteration peak are shown in Fig. 20. These comprise the deuteron transfer from H₂D⁺, or some other deuterated ion, to NH₃ giving NH₃D⁺, dissociative recombination of NH₃D⁺, charge transfer between NH₂D and H⁺ and hydrogen abstraction from H₂ to the NH₂D⁺ ion.

Characteristic of the reaction scheme is the circulation between neutral and ionic species generated by the charge transfer reaction with H⁺, one of the ions which increases in number after the disappearance of CO. Precursors of doubly and triply deuterated ammonia, $\mathrm{N}{\mathrm{H}}_{\mathrm{2}}{\mathrm{D}}_{\mathrm{2}}^{\mathrm{+}}$ and $\mathrm{NH}{\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$, are mainly formed from reactions between NH₃ and D₂H⁺ or ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ in the present model. This reaction network is discussed in more detail in Sipilä et al. (2015b).

The N₂D⁺ ion, which is mainly produced in reactions between N₂ and H₂D⁺, D₂H⁺ or ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$, is strongly favoured by the successive deuteration of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$. At the time of the most vigorous deuteration, the NH₃ abundance decreases slightly owing to accretion onto grains and enhanced charge exchange reactions caused by the increase of H⁺. After that, until very late times, probably exceeding the lifetime of the core, the NH₃ abundance remains almost constant, and so does the abundance of NH₂D. In contrast, the abundances of NHD₂, ND₃, N₂H⁺ and N₂D⁺ rise and fall in the time range shown in Fig. 11. At late times, these species are most strongly influenced by the depletion of nitrogen and deuterium in the gas phase.

The best overall agreement between the modelled and observed NH₂D, NHD₂ and ND₃ spectra is achieved at t ~ 3 × 10⁵ yr. At this stage of the model, the fractionation ratios are NH₂D/NH₃ ~ 0.25, NHD₂/ NH₂D ~ 0.25 and ND₃/ NHD₂ ~ 0.1 (while the observed fractionation ratios are NH₂D/NH₃ ~ 0.4, NHD₂/ NH₂D ~ 0.2 and ND₃/ NHD₂ ~ 0.06). The time is coincident with the NHD₂ maximum, whereas ND₃ is already decreasing at that time. The obtained fractionation ratios are reasonably close to what was found in previous observational studies (Roueff et al. 2005; Daniel et al. 2016b), but it should be noted that our results suggest that large temporal variations are possible, also when the physical conditions remain constant. The spectral line simulations show that all three deuterated forms of ammonia should be easily detectable from a core such as H-MM1 even if the fractionation ratios were reduced to half of those derived here. In particular, the early formation of ND₃ and the large transition dipole moment of the rotation-inversion transition of mND₃ at 309.9 GHz makes this line a useful signpost of the deuterium peak.

According to the present chemistry model, the fractionation ratios on grain surfaces are lower than those in the gas phase by a factor of two. The atomic D^∗/H^∗ ratio on the grain surfaces, which determines the overall degree of deuteration there, reaches a high value of ~ 0.4 a little before 10⁵ yr in the present simulation (Fig. 12). Because of competition between various viable addition reactions for H^∗ and D^∗ (for example with NO^∗, HCO^∗ and HS^∗), the abundances of deuterated forms of ammonia build up slowly. Towards the end of the simulation, the abundances settle, however, almost exactly to the values expected from the statistical rule $\mathrm{N}{\mathrm{H}}_{\mathrm{2}}{\mathrm{D}}^{\mathrm{\ast }}\mathit{/}\mathrm{N}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{\ast }}\mathrm{=}\frac{\mathrm{3}}{\sqrt{\mathrm{2}}}{\mathrm{D}}^{\mathrm{\ast }}\mathit{/}{\mathrm{H}}^{\mathrm{\ast }}$, $\mathrm{NH}{\mathrm{D}}_{\mathrm{2}}^{\mathrm{\ast }}\mathit{/}\mathrm{N}{\mathrm{H}}_{\mathrm{2}}{\mathrm{D}}^{\mathrm{\ast }}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}{\mathrm{D}}^{\mathrm{\ast }}\mathit{/}{\mathrm{H}}^{\mathrm{\ast }}$, $\mathrm{N}{\mathrm{D}}_{\mathrm{3}}^{\mathrm{\ast }}\mathit{/}\mathrm{NH}{\mathrm{D}}_{\mathrm{2}}^{\mathrm{\ast }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}\sqrt{\mathrm{2}}}{\mathrm{D}}^{\mathrm{\ast }}\mathit{/}{\mathrm{H}}^{\mathrm{\ast }}$, as predicted by Brown & Millar (1989).

6.3. Spin ratios

In the gas phase, the ortho/para ratio of NH₂D is largely determined by the cycle consisting of reactions with H⁺, H₂ and e⁻ shown in Fig. 20. In the present chemistry model, full scrambling of H nuclei is assumed to take place in these reactions, and, owing to nuclear spin selection rules, o/pNH₂D should settle to approximately 2.3 (Sipilä et al. 2015b). The full reaction set predicts that the ratio decreases to approximately 2.0 at late times (Fig. 11).

Similar cycles involving NHD₂ and ND₃ preserve the spin states of D₂ and D₃. Consequently, the spin ratios of doubly and triply deuterated ammonia are determined by the primary deuteration reactions D₂H⁺ + NH₃ and ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{N}{\mathrm{H}}_{\mathrm{3}}$. The reaction scheme is discussed in detail in Sipilä et al. (2015b). The fact that o/pNHD₂ follows o/pD₂H⁺ closely can be seen in Figs. 11 and 13. A tight correlation between m/p-ND₃ and m/p-${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ is also evident from these figures.

By comparing the modelled and observed spectra shown in Fig. 16, one finds that while the modelled o/pNH₂D ratio is too low, the corresponding ratio for NHD₂ is too high. The discrepancy is more pronounced in the case of pNHD₂ for which the predicted spectrum underestimates the observed intensity by approximately 40%. As mentioned in Sect. 5.2, the timing and the strength of the deuterium peak can be affected by the selection of the initial o/p-H₂ ratio, the cosmic ray ionisation rate and the average grain size, but these modifications do not change the fact that the o/pNHD₂ ratio given by our model during the deuterium peak is larger than the observed ratio. If the adopted deuteration scheme is correct, the implication is that also the o/pD₂H⁺ ratio is overestimated in the model, because o/pNHD₂ is directly related to o/pD₂H⁺.

However, the spin ratio of D₂H⁺ is determined by the well-studied ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}$ isotopic system (Hugo et al. 2009), and the prediction of a high o/pD₂H⁺ ratio seems to be well-founded. The lower energy, ortho form of D₂H⁺ is favoured over pD₂H⁺ both in the primary production through H₂D⁺ + HD → D₂H⁺ + H₂ and in the “backward” reactions $\mathrm{m}\mathit{/}\mathrm{o}{\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{o}{\mathrm{H}}_{\mathrm{2}}\mathrm{\to }\mathrm{o}\mathit{/}\mathrm{p}{\mathrm{D}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}\mathrm{HD}$ and o/pD₂H⁺ + oH₂ → o/pH₂D⁺ + HD. In addition, para - ortho conversion of D₂H⁺ is viable in cold clouds through reaction with HD (see discussion in Flower et al. 2006b; and Sipilä et al. 2010).

6.4. Statistical abundance ratios

In their analysis of the deuterated ammonia observations towards Barnard 1 and L1689N, Daniel et al. (2016b) concluded that the observational data (considering the error margins) are consistent with the assumption that both the fractionation ratios and the spin ratios are equal to the corresponding statistical ratios. Daniel et al. (2016b) point out that statistical fractionation and spin ratios would be expected if the production of ammonia molecules were dominated by surface reactions. On grain surfaces, where ammonia formation takes place through H/D atom additions to N, the fractionation ratios, NH₂D/NH₃, NHD₂/ NH₂D and ND₃/ NHD₂, should successively diminish by a factor of three (Brown & Millar 1989; Rodgers & Charnley 2001), and the spin ratios should follow the ratios of the corresponding nuclear spin statistical weights.

Also in the present study, the nuclear spin ratios derived directly from the observed lines agree with their statistical values, and NHD₂/ NH₂D ~ 3ND₃/ NHD₂ as expected from combinatorial principles (see Table 4 in Sect. 4.2). The NH₂D/NH₃ ratio falls, however, below the value expected from the other two fractionation ratios. On the other hand, the derived NH₃ abundance for the core may be an overestimate as discussed in Sect. 4.2.

The observed spectra originate, however, in the gas phase, and also sublimated species are likely to be exposed to rapid processing by ion-molecule reactions there. The fractionation ratios on grain surfaces depend on the atomic D^∗/H^∗ ratio on grains, which, according to the present simulation, is, at any time, much higher than the atomic D/H ratio in the gas phase (Figs. 10 and 12). As mentioned in Sect. 6.2, the statistical fractionation ratios for ammonia on grains are only reached at very late times of the simulation.

The observations suggest, therefore, that the adopted gas-phase deuteration scheme for ammonia (illustrated in Fig. 2 of Sipilä et al. 2015b) is not correct. The model assumes complete scrambling of H or D nuclei in the intermediate reaction complexes. This assumption is highly uncertain for several reactions involved in the production of ammonia (Rist et al. 2013).

An argument against a full scrambling of the reaction forming the ammonium ion, $\mathrm{N}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\to }\mathrm{(}\mathrm{N}{\mathrm{H}}_{\mathrm{5}}^{\mathrm{+}}{\mathrm{)}}^{‡}\mathrm{\to }\mathrm{N}{\mathrm{H}}_{\mathrm{4}}^{\mathrm{+}}\mathrm{+}\mathrm{H}$, comes from the energetics of this reaction (an intermediate reaction complex is indicated here with ()^‡). The reaction is likely to occur through two minima on the potential energy surfaces, $\mathrm{(}\mathrm{N}{\mathrm{H}}_{\mathrm{3}}\mathrm{·}\mathrm{·}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}{\mathrm{)}}^{‡}$ and (NH₄··H⁺)^‡. According to the calculations of Ischtwan et al. (1992), some of the conceivable interchanges of two H nuclei between different parts of these complexes are facile, but most of them involve high energy $\mathrm{(}\mathrm{N}{\mathrm{H}}_{\mathrm{5}}^{\mathrm{+}}{\mathrm{)}}^{‡}$ transition states with different geometries. The D-substituted cases should be very similar. However, since it is known experimentally that, for example, D/H exchange occurs between $\mathrm{N}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and D₂, Ischtwan et al. (1992) proposed that it is possible to circumvent a high-energy transition state by a process where hydrogen transfer is followed by internal rotation and reverse transfer. In none of the configurations of $\mathrm{(}\mathrm{N}{\mathrm{H}}_{\mathrm{5}}^{\mathrm{+}}{\mathrm{)}}^{‡}$ do all H (or D) nuclei occupy equivalent positions. Consequently, the spin symmetry rules for this reaction are far from obvious, even less so given that internal rotations are likely to occur.

For the reaction NH₃ + H₂D⁺ and its doubly deuterated analogue (Fig. 21 below), the situation is even less clear. While the reaction proceeds at 30 K (Lindinger et al. 1975; Marquette et al. 1989), there is no theoretical computation on the (NH₅D⁺)^‡ complex; see for example Rist et al. (2013). There is no current experimental or theoretical discrimination between scrambling or proton hop mechanisms for this reaction.

In the present model, we assume that the reaction NH₃ + D₂H⁺ forms the complex $\mathrm{(}\mathrm{N}{\mathrm{H}}_{\mathrm{4}}{\mathrm{D}}_{\mathrm{2}}^{\mathrm{+}}{\mathrm{)}}^{‡}$, which can dissociate to $\mathrm{N}{\mathrm{H}}_{\mathrm{4}}^{\mathrm{+}}$, NH₃D⁺ or $\mathrm{N}{\mathrm{H}}_{\mathrm{2}}{\mathrm{D}}_{\mathrm{2}}^{\mathrm{+}}$. If this is true, the reaction is one of the main sources of NHD₂. If the reaction complex is not formed, but D₂H⁺ simply donates the proton or one of the deuterons to ammonia, the outcome after dissociative recombination could be either NH₃ or NH₂D, as illustrated in Fig. 21. The statistical branching ratios are indicated in this figure.

It seems that several important reactions, which in the model have been assumed to proceed through long-lived intermediate complexes where nuclei can be scrambled, should rather be described as proton/deuteron hops or hydrogen/deuterium abstractions. In these reactions, H and D nuclei are effectively added one by one, as in surface reactions, and they produce each nuclear spin modification according to its statistical weight.

The assumption that ammonia is primarily processed in reactions with the isotopologues of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, such as the one described in Fig. 21, can eventually lead to a ratio of 3 between successive levels of deuteration, although now the fractionation ratios would not depend on D/H nor D^∗/H^∗, but on the relative abundances of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, H₂D⁺, D₂H⁺ and ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$. Using combinatorics, one can show that in steady state NH₂D/NH₃ = 3γ, NHD₂/ NH₂D = γ, and $\mathrm{N}{\mathrm{D}}_{\mathrm{3}}\mathit{/}\mathrm{NH}{\mathrm{D}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathit{\gamma }$, where $\begin{array}{ccc}\mathit{\gamma }\mathrm{=}\frac{\mathrm{\left[}{\mathrm{H}}_{\mathrm{2}}{\mathrm{D}}^{\mathrm{+}}\mathrm{\right]}\mathrm{+}\mathrm{2}\mathrm{\left[}{\mathrm{D}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right]}\mathrm{+}\mathrm{3}\mathrm{\left[}{\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\right]}}{\mathrm{3}\mathrm{\left[}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\right]}\mathrm{+}\mathrm{2}\mathrm{\left[}{\mathrm{H}}_{\mathrm{2}}{\mathrm{D}}^{\mathrm{+}}\mathrm{\right]}\mathrm{+}\mathrm{\left[}{\mathrm{D}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right]}}\mathrm{·}& & \end{array}$Nyman (2016) has recently discussed a simple statistical model for the deuteration of interstellar ammonia that neglects the reaction kinetics. In its simplest form, this model leads to the same rule for the fractionation ratios, NH₂D/NH₃ = 3 NHD₂/ NH₂D = 9 ND₃/ NHD₂, as the models discussed above, but now the parameter γ = NHD₂/ NH₂D depends on the elemental N/D ratio through γ = 1/(3 N/D−1). The fractionation ratio NHD₂/ NH₂D ~ 0.2 observed in H-MM1 and in L1689N (Roueff et al. 2005) would imply N/D ~ 2, which is approximately 30% lower than the ratio assumed in the chemistry model used here. The inclusion of energetics into this model through the molecular partition functions makes the fractionation ratios correspond to what would be obtained at LTE, giving NH₂D/NH₃ ≤ NHD₂/ NH₂D. While this agrees with the earlier results presented in Roueff et al. (2005), the fractionation ratios found in the present study, and in the re-analysis of the B1b results by Daniel et al. (2016b), suggest NH₂D/NH₃ ~ 2−3 NHD₂/ NH₂D. Therefore, the assumption of thermal equilibrium between different deuterated isotopologues of ammonia does not seem to be universally valid. This situation is also unlikely because large deviations from thermal equilibrium are usually found in cold interstellar gas. On the other hand, as recommended by Nyman (2016), the energetics, and especially the differences in the vibrational zero-point energies between different ammonia isotopologues, should be taken into account in kinetic models. This would particularly affect reactions working against deuteration.

7. Conclusions

Ammonia and its three deuterated isotopologues were detected towards the starless core H-MM1 in Ophiuchus. By modelling the observed spectra, we derived the following fractionation ratios: NH₂D/NH₃ ~ 0.4, NHD₂/ NH₂D ~ 0.2 and ND₃/ NHD₂ ~ 0.06. The relative line intensities of the ortho and para modifications of NH₂D and NHD₂ are consistent with the statistical spin ratios o/pNH₂D = 3 and o/pNHD₂ = 2. The fractionation and spin ratios are similar to those obtained towards two young cores (L1689N and B1), which are the only objects previously observed in all four molecules (Roueff et al. 2005; Roueff et al. 2015; Daniel et al. 2016b).

The observations towards H-MM1 were simulated using a gas-grain chemistry model in conjunction with a Monte Carlo radiative transfer program. In the chemistry model, which includes cosmic-ray-induced desorption and tunnelling through activation energy barriers in surface reactions, ammonia forms early and stays long in the gas phase. The model over-predicts pNH₃ and oH₂D⁺, but can approximately reproduce the observed NH₂D, NHD₂, ND₃, N₂D⁺ and pD₂H⁺ lines. The oH₂D⁺ and pD₂H⁺ observations used here are from Parise et al. (2011). According to the simulation, the longevity of NH₃ and N₂H⁺ in the gas phase can be traced back to the chemical inertness of N₂ on the grain surfaces. Unlike CO, the nitrogen molecule is not supposed to be hydrogenated on grain surfaces, which makes it susceptible to desorption.

The present chemistry model cannot account for the observed o/pNH₂D and o/pNHD₂ ratios satisfactorily. In conditions prevailing in H-MM1 and in other starless, dense cores, characterised by low temperature, high obscuration and inefficient desorption, the spin ratios should be determined by gas-phase ion-molecule reactions. When complete scrambling of H and D nuclei in these reactions is assumed, the spin ratios are predicted to settle to o/pNH₂D ~ 2 and o/pNHD₂ ~ 3−4 by the time all three deuterated forms of ammonia become detectable, instead of their statistical ratios 3 and 2, respectively. The fact that the observed spin ratios nevertheless correspond to the nuclear spin statistical weights suggests that full scrambling in reactions forming deuterated ammonia is not a valid assumption. At the moment, there are very little experimental data and few theoretical calculations concerning the probability of proton/deuteron scrambling in the principal deuteration reactions of ammonia.

In contrast to what our chemistry model predicts, the currently available observational data suggest that the nuclear spin ratios of deuterated ammonia isotopologues do not depend strongly on physical conditions or time. On the other hand, the degree of deuterium fractionation does show clear variations from source to source, and is also predicted to have a temporal pattern, depending on CO depletion, the ortho/para-H₂ ratio and eventually on HD depletion. The deuterium fractionation of ammonia has therefore the potential to be used as a probe of dense core evolution. While the present core model is static, chemical reaction rates depend strongly on the density, that is, the dynamical evolution. Therefore, definitive conclusions on the effects of cloud evolutionary stages on deuterium fractionation await models where deuterium chemistry is coupled to dynamical evolution.

---

Acknowledgments

We thank the anonymous referee for insightful comments which helped to improve the manuscript. We thank Malcolm Walmsley and Juris Kalvāns for helpful discussions. J.H. and L.W. thank the Max-Planck-Institute for Extraterrestrial Physics for generous support. J.H. acknowledges financial support from the Academy of Finland grant 258769. P.C., J.P., and A.P. acknowledge the financial support of the European Research Council (ERC; project PALs 320620). F.D., L.W., A.F., and C.R. thank the Agence Nationale de la Recherche (ANR-HYDRIDES), contract No. ANR-12-BS05-0011-01, and the CNRS national program “Physique et Chimie du Milieu Interstellaire”.

References

All Tables

All Figures