© The Authors 2024

1 Introduction

Astrochemical models that distinguish between different nuclear spin modifications assume that ion-molecule reactions in interstellar gas proceed via the formation of a relatively long-lived intermediate complex, where the light H and D nuclei can be completely mixed. In this scenario, one of the main deuteration reactions of ammonia in dense gas can be written as NH₃ + H₂D⁺ → (NH₅D⁺)^† NH₃D⁺ + H₂, where the reaction complex is indicated with a dagger. This reaction is followed by electron recombination of NH₃D⁺, which can yield NH₃ or NH₂D. The increasing abundances of ${\text{H}}_3^{+}$, H₂D⁺, ${\text{HD}}_2^{+}$, and ${\text{D}}_3^{+}$ in the interiors of dense, starless cores, and the possibility of mixing between nuclei in gas-phase reactions, leads to a situation where the nuclear spin ratios of ammonia and its deuterated forms deviate from their statistical values. This is caused by the conservation of spin angular momentum and permutation symmetry. In the models of Sipilä et al. (2015) and Hily-Blant et al. (2018), the ortho-to-para ratio (hereafter o/p) of NH₂D approaches 2 or gets slightly below it (instead of the statistical value 3), whereas o/p-NHD₂ reaches up to ~3 (Hily-Blant et al.) or even 4 (Sipilä et al.) at high densities (instead of the statistical 2).

While interesting from a theoretical point of view, this phenomenon also provides a means of determining the origin of species that possess different nuclear spin modifications: hydrogenated species formed on grains should always show statistical spin ratios because they are formed through H or D atom addition reactions, whereas species formed in the gas phase can have non-statistical spin ratios (e.g., Brown & Millar 1989; Sipilä et al. 2015, 2019a; Harju et al. 2017; Hily-Blant et al. 2018).

The few observational determinations of the o/p ratios of NH₂D or NHD₂ in dense cores are consistent with statistical ratios (Gerin et al. 2006; Lis et al. 2006; Daniel et al. 2016; Harju et al. 2017; Wienen et al. 2021). The uncertainties of these estimates are, however, large – typically greater than 30%, and one cannot confidently exclude the deviations from statistical values predicted by chemistry models.

Here we attempt to determine the o/p ratios of NH₂D and NHD₂ in two dense and cold interstellar gas clouds with greater accuracy than achieved previously, using the Large APEX subMillimeter Array (LAsMA) of the Atacama Pathfinder Experiment (APEX). The ground state transitions of the ortho and para modifications of NH₂D and NHD₂ (hereafter oNH₂D, pNH₂D, oNHD₂, and pNHD₂) at λ ~ 0.9 mm were observed with this instrument towards the starless cores H-MM1 and Oph D in the region of the dark cloud L1688 in Ophiuchus. The large-scale structure and general properties of the cores in the L1688 region have been discussed in, for example, Pattle et al. (2015), Friesen et al. (2017), Choudhury et al. (2020), and Ladjelate et al. (2020). Detailed studies of the structures and the ammonia distributions in the two target cores have been presented in Harju et al. (2020) and Pineda et al. (2022; for H-MM1), and in Ruoskanen et al. (2011; for Oph D). The two cores can be classified as prestellar, meaning that they are on the verge of collapse, as they have central densities significantly above 10⁵ cm⁻³, and high levels of deuterium fractionation (e.g., Crapsi et al. 2005; Harju et al. 2017).

In Sect. 2 of this paper, we describe the LAsMA observations, and in Sect. 3 we present the data and estimate the abundances and the o/p ratios of NH₂D and NHD₂ in the two cores. Spin ratios expected from the proton hop and full scrambling scenarios, and the results of previous spin-state chemistry models are discussed in Sect. 4. In Sect. 5, we show results from an updated version of our gas-grain chemistry model. Implications of the observed spin ratios and NH₂D/NH₃ and NHD₂/NH₂D abundance ratios are discussed in Sect. 6. Finally, in Sect. 7 we draw our conclusion from this study.

Fig. 1 Footprint of the APEX/LAsMA array on the 8 μm surface brightness images of H-MM1 (left) and Oph D (right) observed by the Spitzer satellite. The dashed and solid circles on the right represent two orientations of the array in observations towards Oph D. The positions of the central pixels are RA 16h27m58s.5, Dec – 24°33′40″ (J2000.0) for H-MM1, and RA 16h28m30s.2, Dec – 24°18′42″ for Oph D.

2 Observations

The ground-state rotational lines of oNH₂D, pNH₂D, oNHD₂, and pNHD₂ near 333 and 335 GHz were observed with the LAsMA array (Güsten et al. 2008) on the APEX telescope (Güsten et al. 2006) towards the prestellar cores H-MM1 and Oph D in Ophiuchus. The observations were done between April 10 and September 28, 2022, under the project number M-0109.F-9503A-2022. The APEX telescope is located at an altitude of 5100 m at Llano de Chajnantor, in Chile.

The LAsMA array has 7 pixels separated by approximately 40″. Figure 1 shows the footprint of the instrument superposed on the 8 μm maps of the cores measured by the InfraRed Array Camera (IRAC) of the Spitzer Space Telescope. The beamsize of APEX is 18″ (FWHM) at 335 GHz. The pattern was unintentionally tilted by 19°.1 during observations in the Spring and Summer 2022. This was corrected during additional observations in September towards Oph D (solid circles on the right). With two tilt angles, the total number of positions observed towards Oph D is 13.

The backend was composed of Fast Fourier Transform Spectrometers (FFTS4G) that covered two 4 GHz intermediate frequency (IF) bands resolved into 65 536 spectral channels of 61.03 kHz in width (corresponding to ~55m s⁻¹). The lines of ortho and para NH₂D and NHD₂ were observed in the lower sideband (LSB) of the receiver which was centered at 334.145 GHz. The frequencies of the observed transitions are listed in Table 1. Also listed there are the energies of the upper transition level and the Einstein A coefficients of the transitions. The spectroscopic data are adopted from Melosso et al. (2021). At 10 K, the optically thin critical densities of the four transitions (estimated using Eq. (4) of Shirley 2015) lie in a narrow range between 6.9 × 10⁴ cm⁻³ and 7.9 × 10⁴ cm⁻³.

The NH₂D and NHD₂ spectra towards the density peaks of the two cores are shown in Fig. 2. The positions correspond to pixel ① for H-MM1 and pixel ⑦ for Oph D in Fig. 1. The spectra are shown on the main-beam brightness temperature, T_MB, scale. In the conversion from $T_{\text{A}}^{\ast }$ scale, we assumed the following main-beam and forward-scattering and spill-over efficiencies: η_MB = 0.73, η_fss = 0.97. In the course of the measurements, we noticed line intensity changes up to 15%, which probably are caused by occasional rapid variations in the water vapour content of the atmosphere. We estimate that the absolute calibration is accurate to 20%. Because the four lines were observed in the same spectral band, this uncertainty does not concern their relative intensities, and does not propagate to the derived o/p ratios.

Fig. 2 NH2D and NHD2 spectra from the peak positions of H-MM1 and Oph D (positions ① and ⑦ in Fig. 1, respectively). Gaussian fits to the hyperfine structure are shown with thin red lines.

3 Results

3.1 LTE analysis of the spectra

The line parameters and column densities derived from Gaussian fits to the hyperfine structure of the lines observed towards the density peaks of H-MM1 and Oph D are presented in Table 2. The fits are shown with thin red lines in Fig. 2. In the column density estimates we have assumed line-of-sight homogeneity and that the excitation temperature, T_ex, is constant for all rotational transitions of the molecule. Furthermore, the calculation involves the assumption that populations of hyperfine states in a certain rotational level are proportional to their statistical weights. The transitions have three groups which consist of several unresolved components (Coudert & Roueff 2006; Melosso et al. 2021). The linewidth listed in Table 2 refers to the intrinsic width of an individual hyperfine component. One can see in this Table that the NH₂D lines are broader than those of NHD₂, and the difference in H-MM1 clearly exceeds that expected from thermal broadening, ~5 m s⁻¹ at 10 K. The NHD₂ abundance is likely to peak at slightly higher densities and deeper in the core as compared to NH₂D, and the linewidth difference could possibly be understood in terms of lower temperatures and smaller non-thermal velocity dispersion in the interior parts of the cores. The three-dimensional models described below include realistic temperature gradients but assume a constant non-thermal velocity dispersion.

Despite rather high signal-to-noise ratios especially in H-MM1, the optical thicknesses and the excitation temperatures of the pNH₂D and pNHD₂ lines have large relative uncertainties which propagate to the column densities. The total optical thicknesses (sums of three separate groups) of the pNH₂D and pNHD₂ lines derived towards H-MM1 are below 1, suggesting that one could use for the para lines the optically thin approximation, where the column density is proportional to the integrated intensity. In this method, the T_ex must be assumed. We examined the excitation of the observed transitions by applying radiative transfer calculations to the three-dimensional model of H-MM1 used previously by Pineda et al. (2022), and described in more detail in Sect. 3.2. Here, the molecules are assumed to have constant fractional abundances, corresponding to the column densities derived from the hyperfine fits described above. The T_ex profiles along a line crossing the density peak (with n(H₂) ~ 1.2 × 10⁶ cm⁻³) are shown in Fig. 3. According to this calculation, the T_ex of pNH₂D line is consistently ~10% lower than that of the oNH₂D line, whereas for the para- and ortho-NHD₂ lines the T_ex ratio is 1. The higher T_ex of the oNH₂D line compared to that of the pNH₂D line is likely to be caused by an increased effect of radiative trapping with larger optical thickness (Shirley 2015; Roueff et al. 2021). In the last column of Table 2 we give the pNH₂D and pNHD₂ column densities derived from the integrated intensities assuming that T_ex,para = 0.9 × T_ex,ortho for NH₂D, and T_ex,para = T_ex,ortho for NHD₂.

Using para column densities from optically thin approximation, the values listed in Table 2 give the following o/p-NH₂D ratios: 3.0 ± 0.2 (H-MM1) and 3.5 ± 0.6 (Oph D). The corresponding o/p-NHD₂ ratios are 2.2 ± 0.5 (H-MM1) and 2.4 ± 1.6 (Oph D). If we use the assumption T_ex,para = T_ex,ortho also for NH₂D, the o/p-NH₂D ratios increase by ~20% from the values given above. To distinguish between different formation scenarios for these molecules, the o/p ratios should be determined more accurately, which can be achieved via non-LTE modelling as we will discuss below.

Besides the uncertainty involved in the determination of the optical thickness, τ, from only three resolved components, the assumption of a constant T_ex along the line of sight in the presence of large density gradients is an additional source of uncertainty, as demonstrated in Fig. 3. In what follows, we construct three-dimensional physical models for the H-MM1 and Oph D cores, and estimate the fractional abundances of the observed species using radiative transfer calculations.

Fig. 3 Excitation temperature profiles of the ground-state lines of ortho and para NH2D and NHD2 in a core model resembling H-MM1.

3.2 Three-dimensional radiative transfer modelling

We derived the fractional abundances of the observed molecules by comparing simulated spectra with the observations. Physical models for H-MM1 and Oph D were constructed using the 8 μm surface brightness maps shown in Fig. 1. Both cores appear as absorption features in the maps, and this was used to derive highresolution H₂ column density maps. The method is described in Appendix A of Harju et al. (2020) and in Pineda et al. (2022). The density structures of the cores were estimated by fitting a Plummer-type function to the cross-sectional H₂ column density profiles in different positions along the ridges of the cores. The method is adopted from Arzoumanian et al. (2011). We assumed that the density distribution has circular symmetry in the plane perpendicular to the spine of the core. The latter is defined by a spline fit through local N(H₂) peaks. A more detailed description of the method is presented in Appendix D of Pineda et al. (2022).

The H₂ column densities of the models of H-MM1 and Oph D are comparable to those that can be obtained from the high-resolution N(H₂) map of L1688 derived from multi-wavelength Herschel data (Ladjelate et al. 2020), which has an effective resolution of 18″.2. When the model column density distributions are smoothed to the same resolution, the density peaks have the following H₂ column densities (values from the Herschel map of Ladjelate et al. 2020 are given in brackets): H-MM1 (0, 0) N(H₂) = 3.7 × 10²² cm⁻² (4.4 × 10²² cm⁻²); Oph D (−19″, −37″) N(H₂) = 3.2 × 10²²cm⁻² (2.2 × 10²² cm⁻²); Oph D (+20″, +35″) N(H₂) = 2.2 × 10²² cm⁻² (2.2 × 10²² cm⁻²). The fact that the total N(H₂) column density derived from Herschel is lower than that of the core model at Oph D (−19″, −37″) can be understood by the relatively high dust colour temperature (T_C = 12.9 K) derived from the far-infrared data. The mass-averaged physical temperature is probably lower towards this position (see below). The high-resolution N(H₂) and T_C maps of L1688 from Herschel are accessible via a website set up by the CEA-Saclay¹.

The gas temperature was assumed to be equal to the dust temperature, T_dust, in the inner parts of the cores. Based on the T_kin distribution in H-MM1 derived from NH₃ by Pineda et al. (2022, see their Figs. 5 and 11) we assumed, however, that the gas temperature does not rise above 11 K in the outer parts of the cores owing to molecular line cooling. The dust temperature was computed using the continuum radiative transfer program CRT (Juvela 2005). In the case of H-MM1, the external radiation field was assumed to consist of an isotropic component and a point source located on the western side of the core. The point source represents the subgiant binary HD 147889 located about 1.2 pc west of H-MM1. The density model of H-MM1 and the assumptions about the external radiation field are the same as those used in Pineda et al. (2022). For Oph D, which is known to be cooler than H-MM1 (~6 – 7 K in the densest parts; Harju et al. 2008; Ruoskanen et al. 2011), the external radiation field was assumed isotropic and weaker by a factor of two than that around H-MM1. 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². The density and dust temperature models of H-MM1 and Oph D are shown in Fig. 4. The assumed levelling off of the gas temperature at 11 K only affects the model of H-MM1, where T_dust rises to ~14 K in the outer parts.

The fractional abundances were determined by applying minimum chi-square estimation. We ran three sets of radiative transfer calculations with different assumptions about the abundance distributions: (a) one assuming that the fractional NH₂D and NHD₂ abundances are constant throughout the core, (b) another assuming that the abundances are constant in the outer parts but start to decrease when the H₂ density exceeds 2 × 10⁵ cm⁻³, and (c) a third one assuming that the abundances are constant in the density range 5 × 10⁴ cm⁻³ < n(H₂) < 2 × 10⁵ cm⁻³ but decrease both towards the outer edge and the core centre. In models (b) and (c), the fractional abundances in the central parts were assumed to be inversely proportional to the density, that is, X ∝ n⁻¹, owing to freezing onto the grains. This tendency was recently found in H-MM1 for NH₃ (Pineda et al. 2022), and it is possible that also the deuterated forms of ammonia behave that way (see, e.g., Caselli et al. 2022). In model (c), the abundance in the outer part was assumed to be, somewhat arbitrarily, directly proportional to the density, X ∝ n. At the core edge, the abundances of NH₂D and NHD₂ are expected to increase towards higher densities as the production of these molecules is favoured by CO depletion. This behaviour is also predicted by our chemistry models.

We varied the undepleted baseline abundances, and calculated from the three-dimensional model described above spectral cubes that were smoothed to the 20″ resolution of APEX. Thereafter spectra corresponding to the pixels of LAsMA were extracted. The simulations were done using the LOC radiative transfer program (Juvela 2020). We compared the simulated and observed spectra in the brightest positions where the signal-to-noise ratios of the pNH₂D spectra exceed 10. For H-MM1, we included positions (0, 0) and (+8″, +39″), and for Oph D, we used positions (0, 0), (+20″, +35″), and (−19″, −37″) for fitting the fractional abundances.

For each test value of the fractional abundance and in each position included, we evaluated the statistic

$${\chi }^2=\underset{i}{{{\displaystyle \sum }}^{\text{}}}\frac{{\left(T_{\text{obs},i}-T_{\text{mod},i}\right)}^2}{{\sigma }^2},$$

where T_obs,i is the observed intensity and T_mod,i is the computed intensity in channel i, and σ is the RMS noise of the observed spectrum. The summation is over approximately 100 channels including the three hyperfine groups. The sum of the obtained χ² values over two or three positions shows a parabolic distribution as a function of the fractional abundance. The x-coordinate of the minimum corresponds to the best-fit fractional abundance. We fitted a parabola to this distribution, and denote the fitted function with ${\tilde{\chi }}^2$. The best-fit values, $\overline{X}$, and their 1σ errors, σ_X, were calculated from

$$\overline{X}={{\displaystyle \int }}^{\text{}}X{P}_{\tilde{\chi }}(X)dX,{\sigma }_X^2={{\displaystyle \int }}^{\text{}}{(X-\overline{X})}^2{P}_{\tilde{X}}(X)\text{d}X,$$

where $P_{\tilde{\chi }}\propto {\text{e}}^{{\tilde{\chi }}^2/2}$ is the normalised probability distribution. The method is based on the assumption that the error distribution of the spectra is Gaussian (Andrae 2010). We tested the error estimates of the X values using Monte Carlo. To this end, we generated a few thousand realisations of the spectra by adding normally distributed noise (σ = RMS) to each channel, and calculated the ‘best-fit’ X values by minimum chi-squared estimation. The distribution of these X values is Gaussian and the range containing 68% of the values agrees with $\overline{X}\pm {\sigma }_X$ derived above. The resulting maximum and average fractional abundances and their 1 σ uncertainties are presented in Table 3. By ‘maximum’ abundance we mean the undepleted, baseline abundance, and by ‘average’ abundance we mean the column density ratio N(mol)/N(H₂) in the direction of the density peak of the core model. The average abundances are used for comparison with the predictions of the chemistry model in Sect. 5.

The assumption of molecular depletion at high (or low) densities results naturally in higher baseline (undepleted) abundances and lower average abundances towards the density peak compared to the assumption of a constant X. However, the relative abundances of the four species are approximately the same for all three X profiles. As the main interest of this study is in the o/p ratios, we do not attempt to optimise the depletion power laws or the density thresholds. The simulated spectra using the X=constant model are shown in Appendix A together with the observed spectra.

One can see from Table 3 that the best-fit o/p-NH₂D and o/p-NHD₂ ratios derived towards both cores are close to their statistical values 3 and 2, respectively, and that these ratios do not depend much on the abundance profile assumed. For H-MM1, the 3 σ limits are within 12% from the statistical ratios. For Oph D, both o/p-NH₂D and o/p-NHD₂ are higher than those derived in H-MM1, and the 3 σ limits allow deviations up to 15% and 24% in excess of the statistical values 3 and 2, respectively.

The fractionation ratio NHD₂/NH₂D is lower in H-MM1 (0.23–0.25) than in Oph D (0.27–0.29), where an exceptionally high N₂D⁺/N₂H⁺ column density ratio, 0.44 ± 0.08, has been derived previously by Crapsi et al. (2005). The baseline fractional NH₂D abundance of ~9 × 10⁻⁹ (ortho+para) obtained here for H-MM1 assuming molecular depletion, together with the (undepleted) fractional NH₃ abundance 2 × 10⁻⁸ (ortho+para) derived by Pineda et al. (2022), imply an NH₂D/NH₃ abundance ratio of 0.45. The o/p ratios of NH₂D and NHD₂, and the abundance ratios NH₂D/NH₃ and NHD₂/NH₂D derived in H-MM1 are consistent with the previous estimates by Harju et al. (2017). The noise level of the present 0.8 mm spectra is a factor of two lower than those used in the earlier study, but perhaps a more important difference is that the analysis is now based on several positions and a realistic three-dimensional model is used for the core.

Fig. 4 Density and dust temperature models for H-MM1 and Oph D. Top: H-MM1 The isodensity surfaces (left) are 104 cm−3 (green), 105 cm−3 (purple), 5 × 105 cm−3 (red), and 106 cm−3 (blue). The isotemperature Tdust surfaces (right) are 14 K (green), 12 K, (purple), 11 K (red), and 10 K (blue). Bottom: OphD The isodensity surfaces (left) are 104 cm−3 (green), 105 cm−3 (purple), 2 × 105 cm−3 (red), and 4 × 105 cm−3 (blue). The isotemperature Tdust surfaces (right) are 11 K (green), 10 K (purple), 9 K (red), and 8 K (blue). The positive X-axis points to the observer, and the Y-axis points to the equatorial west. The scale is in arcseconds. In both cores, the gas kinetic temperature, Tkin, is assumed to be equal to Tdust in the inner parts, and settle to 11 K in the outer parts.

4 Spin ratios from the proton hop and full scrambling scenarios and the predictions of previous chemistry models

Deuterium enrichment of NH₃ in dense interstellar gas is supposed to occur primarily through deuteron donation from H₂D⁺, D₂H⁺, or D⁺ to ammonia. Deuterium can also enter the ammonia formation sequence in reactions of other nitrogen hydrides, such as NH, with the cations H₂D⁺, D₂H⁺, and D⁺. These exothermic reactions can be thought to proceed via a direct proton or deuteron hop (hereafter PH) or via the formation of an intermediate complex where H and D nuclei can be exchanged. In the latter mechanism, called full scrambling (hereafter FS), the nuclear spin state of the resulting cation, for example NH₃D⁺, depends on the initial spin states of both reactants according to the selection rules that take into account the conservation of the spin and symmetry (e.g., Rist et al. 2013; Faure et al. 2013; Sipilä et al. 2015; Hily-Blant et al. 2018). The spin states of the H₂D⁺, D₂H⁺, and ${\text{D}}_3^{+}$ ions, in turn, are determined by spin-state conversions in the ${\text{H}}_3^{+}$ + H₂ reaction system and its deuterated variants. The state-to-state and ground-state-to-species rate coefficients of this system were derived by Hugo et al. (2009), and these have been used in spin-state chemistry models ever since (e.g., Pagani et al. 2009; Sipilä et al. 2010; Le Gal et al. 2014; Hily-Blant et al. 2018).

To illustrate how the reaction details affect the nuclear spin ratios, we discuss the outcome of a basic sequence leading to NH₂D, assuming that this follows either the PH or the FS scenario. The sequence is initiated by NH₃ + H₂D⁺ and leads to NH₂D through dissociative electron recombination of NH₃D⁺. Each reaction step has several possible products, but here we concentrate on the branch ending with NH₂D. The spin-separated sequences according to the PH and FS scenarios are shown in Fig. B.1.

In the PH mechanism, H₂D⁺ simply donates the deuteron to ammonia, and the spin symmetry of the hydrogen complex does not change. The H₃ complex of E (“para”) symmetry dissociates to form oH₂ + H and pH₂ + H at equal probabilities, whereas the A₁ (“ortho”) symmetry dissociates to form oH₂ + H (if not H + H + H). Assuming equal amounts of oNH₃ and pNH₃, this sequence produces three times more oNH₂D than pNH₂D. In the FS scenario, the first step forms the reaction complex (NH₅D⁺)^†. With reactants oNH₃ and pH₂D⁺, the resulting H₅ complex is of G₁ symmetry, which gives H₃(A₁) (“ortho”) and H₃(E) (“para”) at equal probabilities when dissociating (Hugo et al. 2009, Table III; Hily-Blant et al. 2018, Table C2). The H₅ complex formed in the reaction pNH₃ + pH₂D⁺ has H₁ symmetry with strong preference to H₃(E) in dissociation. Assuming again equal amounts of oNH₃ and pNH₃, the FS sequence results in an o/p-NH₂D ratio of 27/13 ≈ 2.1.

The difference between the PH and FS scenarios is similar for other routes to NH₂D. For example, according to the PH model, the sequence ${\text{NHD}}^{+}\stackrel{{\text{H}}_2}{\to }{\text{NH}}_2{\text{D}}^{+}\stackrel{{\text{H}}_2}{\to }{\text{NH}}_3{\text{D}}^{+}\stackrel{{\text{e}}^{-}}{\to }{\text{NH}}_2\text{D}$ produces again oNH₂D and pNH₂D in the statistical ratio 3:1, whereas the same sequence in the FS model results in an o/p-NH₂D ratio of 7/5 = 1.4, assuming that the reacting H₂ is always in the more common para form. Finally, the formation of NHD₂ involves the combination of two deuterium nuclei. If this happens one by one, as in the PH model, the resulting o/p-NHD₂ ratio is always 2. In the FS model, the outcome depends on the deuteron donor. For example, in a reaction between a singly deuterated species and oD₂H⁺, the o/p ratio of the D₂ complex formed is 13/5 = 2.6. The branching ratios of the reactions discussed here can be verified by inspecting the nuclear spin symmetry induction and subduction matrices for systems with multiple H and D atoms presented in Hugo et al. (2009), Sipilä et al. (2015), and Hily-Blant et al. (2018).

The fact that the FS assumption in proton-donation reactions leads to deviations of the o/p ratios of NH₂D and NHD₂ from their high-temperature statistical values was noted by Sipilä et al. (2015), Daniel et al. (2016), and Hily-Blant et al. (2018) (see Table 5 of Hily-Blant et al. 2018). This disagreed with observations which suggested statistical ratios, albeit with large uncertainties. Daniel et al. (2016) pointed out that the observed o/p ratios of NH₂D and NHD₂ are consistent with formation on grains, while Harju et al. (2017) argued that statistical spin ratios could be explained by gas-phase reactions if they can be characterised as single-particle hops, in analogy with H and D atom additions on grains. The latter suggestion was examined by Sipilä et al. (2019a) who considered the effects of PH vs. FS in proton-donation reactions, such as NH₃ + D₂H⁺, and in the hydrogen abstraction chain leading to the ammonium ion, ${\text{N}}^{+}\stackrel{{\text{H}}_2}{\to }{\text{NH}}^{+}\stackrel{{\text{H}}_2}{\to }{\text{NH}}_2^{+}\stackrel{{\text{H}}_2}{\to }{\text{NH}}_3^{+}\stackrel{{\text{H}}_2}{\to }{\text{NH}}_4^{+}$, and its deuterated variants. It turned out that replacing the FS with PH in proton-donation reactions brought the predicted o/p ratios a little closer to the statistical values but the difference was still greater than 30%. On the other hand, by extending the PH scenario also to H abstraction reactions, the o/p ratios agreed with those observed previously towards H-MM1, but the modelled fractionation ratios NH₂D/NH₃ and NHD₂/NH₂D exceeded those derived from the observations.

5 Gas-grain chemistry model with dynamic chemical desorption

One way to explain the discrepancy between the observed and modelled spin ratios is that the desorption of ammonia from grain surfaces is in reality more efficient than what has been considered in the previous simulations. Ammonia does not desorb thermally in low-temperature conditions owing to its high binding energy on ice; usually values around 4000 – 5000 K are assumed for the NH₃ binding energy on water ice (Minissale et al. 2022; Ferrero et al. 2023). Although recent studies indicate that we should rather speak of a broad distribution of binding energies on amorphous water ice (Grassi et al. 2020; Tinacci et al. 2022), the desorption of ammonia most likely occurs through non-thermal mechanisms.

There are several non-thermal desorption mechanisms that are typically considered in astrochemical models, namely photodesorption, cosmic-ray induced desorption, and chemical desorption (also known as reactive desorption). The first two mechanisms do not release ammonia into the gas phase efficiently enough to cause significant “contamination” of the gas-phase ammonia reservoir. Chemical desorption (CD) on the other hand has been shown to have a very large impact on ammonia abundance in starless core conditions (Sipilä et al. 2019b), assuming that the desorption efficiency is constant at 1% (Garrod et al. 2007). Hence, out of the three options mentioned here, CD is the most promising candidate for pushing the gas-phase ammonia spin ratios toward statistical values. However, the constant efficiency assumption for CD was found to lead to gas-phase ammonia abundances much higher than what is indicated by observations for volume densities around 10⁴ cm⁻³ (Sipilä et al. 2019b).

It is reasonable to posit that in reality the efficiency of CD is in fact not universally constant, but would instead depend not only on the desorbing molecule but also on the type of surface on which the desorption is taking place. An alternative description of CD exploring this scenario has been presented by Minissale et al. (2016), and applied to a chemical model by Vasyunin et al. (2017). An important conclusion in these works is that the CD efficiency tends to increase when switching from a water ice surface to a CO ice surface, which in a starless core would correspond to the transition from outer regions to the dense inner core where the freeze-out timescale is short and the grains are covered by a multitude of species, and not only water. Motivated by these findings, we investigate here whether the dynamic nature of the CD process influences the spin ratios of deuterated ammonia–we recall that the observed ammonia lines originate in volume densities approaching 10⁵ cm⁻³ where the ices on the grains are made up of a mixture of molecules.

We use for the simulations the chemical model described in Sipilä et al. (2019a), and adopt the full scrambling and proton hop chemical networks presented therein. In the proton hop model, proton donation and abstraction reactions proceed via the exchange of one proton (or deuteron) between the reactants, while in the full scrambling model multiple atom exchanges are possible. Also, for the sake of simplicity, we employ a one-dimensional physical model for H-MM1 recently presented in Pineda et al. (2022; their Fig. 11). The adoption of a one-dimensional spherically symmetric model means that we can with reasonable accuracy compare the simulation results to the observations obtained toward the dust peak, but we cannot attempt to reproduce the morphology of the deuterated ammonia distributions beyond a distance of a few tens of arcseconds from the core center.

We consider two different approaches to simulating CD: one where the CD efficiency is set to a constant value of 1% (Garrod et al. 2007), and another where the CD efficiency changes dynamically. The latter is an updated version of the model of Vasyunin et al. (2017) – in the present work, the energy released in a chemical desorption event is partitioned among the reaction products according to their total degrees of freedom. In the dynamic model, the CD efficiency of a given species depends on the ice composition, and hence on the location in the core and on the simulation time, as well as the properties (mass, binding energy) of the species itself. At late simulation times, however, we obtain nearly constant CD efficiencies as a function of radius. For example for NH₂D, the dynamic model predicts a CD efficiency of 4.1% at the center of the core, and 4.9% at the edge, at a simulation time of 5 × 10⁵ yr. A detailed account of the new dynamic CD model will be given in Riedel et al. (2023).

We employ the same initial abundances as in Sipilä et al. (2019a), not recounted here for brevity, and most of the simulation parameters are also the same as in that paper. There are two notable differences: here we assume an external visual extinction A_V = 3 mag following Pineda et al. (2022, 2 mag in Sipilä et al. 2019a), and we use the new description for cosmic-ray induced desorption from Sipilä et al. (2021).

In Figs. 5 to 10, we compare the average fractional abundances, that is, the column density ratios N(mol)/N(H₂) towards the density peaks of the cores, and similarly derived average o/p and deuterium fractionation ratios predicted by the chemical models to those derived by fitting spectra observed towards H-MM1 and Oph D. The fractional NH₂D and NHD₂ abundances relative to H₂ and the abundance ratios derived from observations are column density ratios through the density peak of the three-dimensional model with abundance profiles (b) and (c) provided in Table 3. The ranges indicated with hatched areas correspond to the inverse-variance weighted averages and the standard deviations of the values derived towards H-MM1 and Oph D. For the observed para-NH₃ (hereafter pNH₃) abundance, we adopt the result from Pineda et al. (2022) that the fractional (ortho+para) NH₃ abundance is 2 × 10⁻⁸ (± 10%) in the outer core, but decreases with the density when n(H₂) > 2 × 10⁵ cm⁻³. Pineda et al. (2022) could only derive the pNH₃ abundance, and assumed o/p-NH₃ = 1.

One can see from Figs. 5–7 that the two CD scenarios result in large differences in the average NH₃, NH₂D, and NHD₂ abundances. The fractional pNH₃, NH₂D, and NHD₂ abundances predicted by the constant CD models are close to those derived from the observations, whereas the dynamic CD models seem to produce too much of these molecules. On the other hand, the assumed gas-phase reaction mechanism – either PH or FS – has a minor effect on the abundances, and for most species the predicted abundances from the two models are contained within or lie close to the range of observationally derived values. The only exception is oNH₂D, for which the model assuming constant CD with PH seems to give a better agreement than the FS model.

It is evident from Fig. 8 that both CD models assuming PH agree with the o/p ratios of NH₂D and NHD₂ derived from the observations, whereas the FS models under-predict o/p-NH₂D ratio by approximately 30% and over-predict the o/p-NHD₂ ratio by a similar amount. Including results for both H-MM1 and Oph D, the observationally derived ratios are approximately o/p-NH₂D ~ 2.9 ± 0.2 and o/p-NHD2 ~ 2.1 ± 0.1. As expected, the PH models predict statistical spin ratios. In the FS models, the o/p-NH₂D ratios settle around 1.9–2.0, and the o/p-NHD₂ ratios are found in the range 2.6–2.8 at late times.

According to Fig. 10, the NHD₂/NH₂D fractionation ratio is reproduced by the constant CD PH model at simulation times around 5 × 10⁵ yr, and also the constant CD FS model approaches the observed values at late times of the simulation. The constant CD PH model over-predicts the NH₂D/pNH₃ ratio except at late times (>5.5 × 10⁵ yr), whereas the prediction of the constant CD FS model agrees with the observed ratio only at early times (<4 × 10⁵ yr). The dynamic FS models seem to work better for this ratio.

To summarise, the models assuming constant CD (with PH or FS) come closer to the observed abundances than the dynamic CD models. The models assuming PH in the gas-phase reactions (combined with constant or dynamic CD) can reproduce the observed spin ratios, whereas the FS models predict too low o/p-NH₂D ratios and too high o/p-NHD₂ ratios. Furthermore, the NHD₂/NH₂D ratio can be reproduced by the constant CD PH model at a certain time interval (~3 – 6 × 10⁵ yr), but this model over-predicts the NH₂D/PNH₃ ratio before the simulation time ~5.5 × 10⁵ yr, reaching the observed ratio only later. Both deuterium fractionation ratios predicted by the constant CD PH model are within the observed range only during a short period before 6 × 10⁵ yr.

In an attempt to quantify the overall agreement or disagreement of the four models with the observations, we present in Table 4 the relative deviations of the predicted NH₂D/PNH₃ and NHD₂/NH₂D ratios, o/p ratios, and fractional NH₂D and NHD₂ abundances from the values derived from observations. The relative deviation of the quantity q is defined by Δq = |q_mod – q_obs|/q_obs. The modelled values are taken at the time 5 × 10⁵ yr. The observed values are variance-weighted averages of those derived towards the centres of H-MM1 and Oph D. One can see that both PH models give small relative deviations for the o/p ratios, but the constant CD PH model provides the smallest mean deviation when all quantities are included. As discussed above, the dynamic CD models over-predict the fractional abundances, and this causes the large relative deviations in the last two rows of Table 4 for these models. As can be expected from Figs. 6, 7, 8, and 10, the relative deviations change with time, especially at early times (< 3 × 10⁵ yr), except for the o/p ratios predicted by the PH models which are nearly constant. The constant CD PH model gives the lowest mean deviation at all times. We note that in general, the abundances of NH₂D and NHD₂, and thus the deuterium fractionation ratios, depend on the physical conditions (temperature and density), whereas the spin ratios are directly related to the chemical reaction mechanisms.

Spectra calculated from the three-dimensional cloud model for H-MM1 with abundances interpolated from the one-dimensional chemical model for three simulation times are shown in Figs. C.1 – C.4. The interpolation uses fractional abundances as functions of density in the one-dimensional model. The spectra observed towards the centre of H-MM1 are superposed on the simulated spectra. The abundances from the constant CD model with FS can approximately reproduce the pNH₂D and pNHD₂ lines, but predict too weak oNH₂D lines and too strong oNHD₂ lines (Fig. C.1). The constant CD model with PH (Fig. C.2) reaches rather a good agreement in all four spectra at late times. As expected from Figs. 6 and 7, abundances from the dynamic CD models over-predict the line intensities of all the observed lines, and the difference is particularly stark for the NHD₂ lines. The simulated spectra, which depend on the beam-averaged column densities towards the centre of the core model, support the idea conveyed by Table 4 that the constant CD model with PH agrees better with the observations than the other three chemical models tested here.

Fig. 5 Evolution of the average fractional ortho- and para-NH3 abundances according to different varieties of our chemistry model: “const CD” means chemical desorption (CD) with a constant (1%) efficiency, “dyn CD” means changing efficiency of CD depending on the composition of the grain surface, “FS” means full scrambling of protons in proton donation and hydrogen abstraction reactions, and “PH” means proton hop in these reactions (see text). The averages are column density ratios N(mol)/N(H2) taken through the density maximum of the core model used in chemical simulations. The shaded area in the lower panel shows the range of similarly taken average of the fractional pNH3 abundance through the three-dimensional model of H-MM1 based on the results of Pineda et al. (2022).

Fig. 6 Same as Fig. 5 but for ortho- and para-NH2D. The shaded areas show the ranges of averages through the three-dimensional model with the abundance profiles (b) and (c) used in the fitting of the observed spectra towards H-MM1 and Oph D (see Table 3).

Fig. 7 Same as Fig. 6 but for ortho- and para-NHD2.

Fig. 8 Evolution of the average o/p-NH2D and o/p-NHD2 ratios according to our chemistry model. The labelling is the same as in Fig. 5. The ranges derived from observations using the three-dimensional models for H-MM1 and Oph D with abundance profiles (b) and (c) (see Table 3) are indicated with shaded areas.

Fig. 9 Same as Fig. 8 but for NH3. There is no observational estimate of the o/p-NH3 ratio towards H-MM1 or OphD.

Fig. 10 Evolution of the average deuterium fractionation ratios NH2D/pNH3 and NHD2/NH2D according to the chemistry model. The labelling is the same as in Fig. 5. In the top panel, the shaded area indicates the ratio derived from observations towards H-MM1, where the fractional pNH3 abundance has been determined by Pineda et al. (2022). In the bottom panel, the “observed” ratios include both H-MM1 and Oph D.

6 Discussion

The observations and their analysis presented here indicate that the o/p ratios of NH₂D and NHD₂ in two cold, dense cores deviate at most 20% from their high-temperature statistical values o/p-NH₂D = 3 and o/p-NHD₂ = 2. The result conforms with previous observational determinations, with a higher accuracy than previously attained. Statistical spin ratios do not agree with previously published astrochemical models that assume complete scrambling of light nuclei in gas-phase reactions. These typically predict an o/p-NH₂D ratio of ~2 or slightly below (Sipilä et al. 2015; Hily-Blant et al. 2018). For o/p-NHD2, the predictions diverge depending on some differences in complex-forming reactions. In the model of Hily-Blant et al. (2018), o/p-NH₂D remains below 3, whereas it approaches 4 in the original spin-state chemistry model of Sipilä et al. (2015); see Sect. 3.3.4 and Fig. 6 of Hily-Blant et al. (2018).

Using chemical models that include both gas-phase and grain-surface reactions, and two different chemical desorption scenarios, we tested two suggestions that have been put forward to explain the discrepancy between observations and models: 1) gas-phase ammonia has a substantial contribution from grains where H and D addition reactions result in statistical spin ratios (Daniel et al. 2016), and 2) gas-phase reactions that form ammonia and its deuterated isotopologues occur through proton or atom hops instead of particle exchange in a reaction complex (Harju et al. 2017; Sipilä et al. 2019a). According to nuclear spin selection rules, reactions where particles are added one by one should produce species in proportion to their nuclear spin statistical weights.

It turns out that while the adopted chemical desorption mechanism (constant or dynamic) affects strongly the average abundances of ammonia and its deuterated isotopologues (Figs. 5–7), the effect of this mechanism on their spin ratios is less marked (Figs. 8 and 9). Switching from constant to dynamic CD increases the fractional NH₃, NHD₂, and NHD₂ abundances by a factor of 2–3. In the FS scenario, dynamic CD decreases the o/p-NH₂D and o/p-NHD₂ ratios by ~10%, and increases the o/p-NH₃ ratio by a similar amount, but has hardly any effect on the spin ratios predicted by the PH model. The poor agreement of the dynamic CD model with the observed fractional abundances suggests that desorption is not as efficient as assumed there, rendering the constant CD model the more realistic of the two.

An input of ammonia from grains seems plausible in view of the fact that the fractional abundance of NH₃ residing on grains is probably several orders of magnitude higher than that in the gas in starless cores (e.g., Tinacci et al. 2022). The recent quantum mechanical simulations of Ferrero et al. (2023) indicate that direct chemical desorption of ammonia from water ice is unlikely, because the reaction energy of hydrogenation is efficiently absorbed by the ice surface. However, the situation may be different in CO-rich ices, which are supposed to be characteristic of dense, starless cores, where also deuterated ammonia is found. The dynamic CD model discussed above is based on these ideas, but it seems to over-predict the desorption efficiency from CO-rich ices.

Our modelling results suggest that the observed statistical o/p ratios of NH₂D and NHD₂ are primarily determined by gas phase reactions and are not an effect of vigorous desorption from grains. Firstly, the spin ratios depend strongly on the assumed gas-phase proton exchange scenario (FS or PH) also in the dynamic CD model with very efficient desorption. Secondly, the predicted deuterium fractionation ratios in the grain mantles (NH₂D/NH₃ < 0.1 and NHD₂/NH₂D ~ 0.1) are lower than those observed in starless dense cores, including the targets of the present study. The observed fractionation ratios agree with the predictions for gas-phase formation. Inspection of the rates of the principal reactions producing and destroying ammonia in a given simulation cell at a certain time shows that desorption and accretion typically account for a few percent of the production and destruction of ammonia at densities (~ 10⁵ cm⁻³) where we believe the observed lines to originate, while the formation through ${\text{NH}}_4^{+}+{\text{e}}^{-}$ and destruction through ${\text{NH}}_3+{\text{H}}_3^{+}$ (and isotopologues), and to a lesser extent through charge exchange between NH₃ and H⁺, dominate (see also Fig. 2 in Sipilä et al. 2015). This means that, according to our simulations, ammonia and its deuterated forms are rapidly re-processed by ion-molecule reactions in the gas phase, even though desorption from grain surfaces may have a significant contribution to their total abundances. By re-processing we mean constant conversion of nitrogen hydrides to cations through charge or proton exchange reactions, and their return to neutral species in recombination reactions with electrons, for example, ${\text{NH}}_3\stackrel{{\text{H}}^{+}}{\to }{\text{NH}}_3^{+}\stackrel{{\text{H}}_2}{\to }{\text{NH}}_4^{+}\underset{{\text{H}}_3^{+}}{\stackrel{{\text{e}}^{-}}{\rightleftarrows }}{\text{NH}}_3$.

Besides the fact that many uncertainties are involved in chemical modelling, also observational determination of molecular abundances is non-trivial. There, we are dealing with changing abundances and excitation conditions along the line of sight. To demonstrate this, we show in Fig. 11 the fractional abundances of NH₃, NH₂D, and NHD₂ as functions of the radial distance from the core centre and the density, predicted by the constant CD model with PH at a certain time of the simulation. When comparing the best-fit abundances derived using the three-dimensional model of H-MM1 to the model predictions, we used 1) column density ratios N(mol)/N(H₂) along a line going through the core centre (Figs. 5 – 10 and relative deviations given in Table 4), and 2) simulated and observed spectra towards the core centre, sampling a cylindrical region including the densest parts of the core (Figs. C.1–C.4). In both comparisons, the constant CD model with PH gave the best overall agreement with the observations. The interpretation of the four spectral lines used here is facilitated by the fact that they lie close in frequency, meaning that they can be observed in the same spectral band and have similar critical densities.

The PH scenario was previously tested by Sipilä et al. (2019a) with a core model resembling H-MM1. They found that PH reactions favour the singly deuterated NH₂D at the cost of NHD₂ and ND₃, but also that all the predicted deuterium fractionation ratios, NH₂D/NH₃, NHD₂/NH₂D, and ND₃/NHD₂, were higher than observed in H-MM1 when the ammonia production only involved proton hop reactions. Here, we do not find large discrepancies between the modelled and observed NH₂D/pNH₃ and NHD₂/NH₂D ratios when using the constant CD PH model. The main difference between the present and the former modelling efforts is that Sipilä et al. (2019a) did not include CD. The replenishment of gas by ammonia formed on grains can diminish the deuterium fractionation as the NH₂D/NH₃ and NHD₂/NH₂D ratios on grains are predicted to be substantially lower than those in the gas. For example, in all four models examined here, the NH₂D/NH₃ ratio on grains always stays below 0.1.

Fig. 11 Fractional abundances of the ortho and para modifications of NH3, NH2D, and NHD2 as functions of radial distance from the centre (left) and density (right) in the one-dimensional cloud model used for chemical calculations. The chemistry model assumes constant CD with PH. The abundances are taken at the simulation time 5 × 105 yr.

7 Conclusions

We have determined the o/p ratios of NH₂D and NHD₂ in two prestellar cores by observing the ground-state rotational lines of these molecules with the LAsMA multi-beam receiver on APEX. By simulating line emission from realistic three-dimensional models of the cores, and using different assumptions regarding the abundance distributions, we could establish that the o/p ratios deviate at most 20% from their statistical values, o/p-NH₂D = 3, o/p-NHD₂ = 2.

We also ran chemical models to test if the observed statistical nuclear spin ratios can be explained by enhanced chemical desorption (CD) in the core, employing two different scenarios for proton or deuteron exchanges in gas-phase chemical reactions. An increase in the efficiency of CD in the core interior can be envisioned because the grain surface composition is expected to change from H₂O-dominated to CO-rich ice. However, models assuming constant CD, with 1% of exothermic reactions leading to desorption, can reproduce the observed abundances better than the so called dynamic CD models with changing efficiency. Despite diverging predictions for the fractional abundances, both constant and dynamic CD models assuming single-particle hop in proton donation and hydrogen abstraction reactions gave the better agreement with the observed o/p ratios than the corresponding models with full scrambling.

The efficiency of CD naturally affects the abundances of NH₃ and its deuterated isotopologues in the gas, but according to our modelling results the effect of the proton exchange scenario in gas-phase reactions is more important to the spin ratios. This is because the accretion and desorption processes are out-competed by rapid conversion of gas-phase nitrogen hydrides to cations through charge or proton exchange reactions, and the return of cations to neutral species in recombination with electrons.

Acknowledgements

We thank the anonymous referee for helpful comments, and Mika Juvela for his advice on error analysis and on the use of the LOC radiative transfer program. We also thank the Max Planck Society for financial support.

Appendix A Observed and simulated spectra

The NH₂D and NHD₂ spectra observed towards H-MM1 and Oph D with LAsMA are shown in Figs. A.1 to A.6. The simulated spectra (shown with red) are calculated from the three-dimensional models assuming that the fractional abundance is constant.

Fig. A.1 NH2D spectra obtained with the LAsMA array towards H-MM1. Simulated spectra assuming constant fractional abundances in the three-dimensional core model are shown with red. The assumed ortho- and para-NH2D abundances are those given in the left column of Table 3, that is, X(oNH2D) = 4.520 × 10−9, X(pNH2D) = 1.469 × 10−9.

Fig. A.2 Observed and simulated NHD2 spectra towards H-MM1. The model abundances X(oNHD2) = 9.25 × 10−10, X(pNHD2) = 4.55 × 10−10.

Fig. A.3 NH2D spectra observed towards Oph D in the Spring and Summer 2022. The model abundances X(oNH2D) = 4.640 × 10−9, X(pNH2D) = 1.469 × 10−9.

Fig. A.4 NH2D spectra towards Oph D observed in September 2022, when the orientation of the LAsMA array was different from that used in the Spring and Summer. The model abundances are the same as in Fig. A.3.

Fig. A.5 NHD2 spectra observed towards Oph D in the Spring and Summer 2022. The model abundances X(oNHD2) = 1. 122 × 10−9, X(pNH2D) = 5.02 × 10−10.

Fig. A.6 NHD2 spectra observed towards Oph D in September 2022. The model abundances are the same as in Fig. A.5.

Appendix B Spin-state branching ratio diagrams for the reaction sequence NH₃→NH₃D⁺ → NH₂D

Fig. B.1 illustrates the spin selection rules of the reaction sequence ${\text{NH}}_3\stackrel{{\text{H}}_2{\text{D}}^{+}}{\to }{\text{NH}}_3{\text{D}}^{+}\stackrel{{\text{e}}^{-}}{\to }{\text{NH}}_2\text{D}$, assuming that this follows either the proton hop (PH) or the full scrambling (FS) scenario. This is one of the principal reaction sequences leading to NH₂D. In the FS diagram, we have assumed that H₂D⁺ is the more common para modification; in the PH mechanism, the spin state of the proton or deuteron donor does not matter.

Fig. B.1 Spin-state branching ratio diagrams for the reaction sequence ${\text{NH}}_3\stackrel{{\text{H}}_2{\text{D}}^{+}}{\to }{\text{NH}}_3{\text{D}}^{+}\stackrel{{\text{e}}^{-}}{\to }{\text{NH}}_2\text{D}$ in the PH (left) and FS (right) models. Assuming that o/p-NH3 = 1, the PH scenario leads to o/p-NH2D = 3, whereas the FS scenario results in o/p-NH2D = 27/13.

Appendix C Simulated spectra using abundances from the chemistry model

Figures C.1 to C.4 show NH₂D and NHD₂ spectra towards the centre of the three-dimensional model of H-MM1 with abundances interpolated from the one-dimensional core model used in the chemical model calculations. The spectra represent three simulation times with different assumptions about the desorption efficiency and the viability of proton or deuteron exchanges in chemical reactions. Also shown are the observed spectra towards the centre of H-MM1.

Fig. C.1 Simulated NH2D and NHD2 spectra towards the centre of the three-dimensional model for H-MM1 with abundance distributions interpolated from the chemistry model assuming constant CD and full scrambling (FS, red lines). The spectra are convolved with the APEx beam. Spectra observed with APEx towards the centre of H-MM1 are shown with blue lines.

Fig. C.2 Same as Fig. C.1 but for a chemistry model assuming constant CD and proton hop (PH).

Fig. C.3 Same as Fig. C.1 but for a chemistry model assuming dynamic CD and full scrambling (FS).

Fig. C.4 Same as Fig. C.1 but for a chemistry model assuming dynamic CD and proton hop (PH).

References

All Tables

All Figures

Fig. 1 Footprint of the APEX/LAsMA array on the 8 μm surface brightness images of H-MM1 (left) and Oph D (right) observed by the Spitzer satellite. The dashed and solid circles on the right represent two orientations of the array in observations towards Oph D. The positions of the central pixels are RA 16h27m58s.5, Dec – 24°33′40″ (J2000.0) for H-MM1, and RA 16h28m30s.2, Dec – 24°18′42″ for Oph D.

In the text

	Fig. 2 NH2D and NHD2 spectra from the peak positions of H-MM1 and Oph D (positions ① and ⑦ in Fig. 1, respectively). Gaussian fits to the hyperfine structure are shown with thin red lines.
In the text

	Fig. 3 Excitation temperature profiles of the ground-state lines of ortho and para NH2D and NHD2 in a core model resembling H-MM1.
In the text

Fig. 4 Density and dust temperature models for H-MM1 and Oph D. Top: H-MM1 The isodensity surfaces (left) are 104 cm−3 (green), 105 cm−3 (purple), 5 × 105 cm−3 (red), and 106 cm−3 (blue). The isotemperature Tdust surfaces (right) are 14 K (green), 12 K, (purple), 11 K (red), and 10 K (blue). Bottom: OphD The isodensity surfaces (left) are 104 cm−3 (green), 105 cm−3 (purple), 2 × 105 cm−3 (red), and 4 × 105 cm−3 (blue). The isotemperature Tdust surfaces (right) are 11 K (green), 10 K (purple), 9 K (red), and 8 K (blue). The positive X-axis points to the observer, and the Y-axis points to the equatorial west. The scale is in arcseconds. In both cores, the gas kinetic temperature, Tkin, is assumed to be equal to Tdust in the inner parts, and settle to 11 K in the outer parts.

In the text

Fig. 5 Evolution of the average fractional ortho- and para-NH3 abundances according to different varieties of our chemistry model: “const CD” means chemical desorption (CD) with a constant (1%) efficiency, “dyn CD” means changing efficiency of CD depending on the composition of the grain surface, “FS” means full scrambling of protons in proton donation and hydrogen abstraction reactions, and “PH” means proton hop in these reactions (see text). The averages are column density ratios N(mol)/N(H2) taken through the density maximum of the core model used in chemical simulations. The shaded area in the lower panel shows the range of similarly taken average of the fractional pNH3 abundance through the three-dimensional model of H-MM1 based on the results of Pineda et al. (2022).

In the text

	Fig. 6 Same as Fig. 5 but for ortho- and para-NH2D. The shaded areas show the ranges of averages through the three-dimensional model with the abundance profiles (b) and (c) used in the fitting of the observed spectra towards H-MM1 and Oph D (see Table 3).
In the text

	Fig. 7 Same as Fig. 6 but for ortho- and para-NHD2.
In the text

	Fig. 8 Evolution of the average o/p-NH2D and o/p-NHD2 ratios according to our chemistry model. The labelling is the same as in Fig. 5. The ranges derived from observations using the three-dimensional models for H-MM1 and Oph D with abundance profiles (b) and (c) (see Table 3) are indicated with shaded areas.
In the text

	Fig. 9 Same as Fig. 8 but for NH3. There is no observational estimate of the o/p-NH3 ratio towards H-MM1 or OphD.
In the text

Fig. 10 Evolution of the average deuterium fractionation ratios NH2D/pNH3 and NHD2/NH2D according to the chemistry model. The labelling is the same as in Fig. 5. In the top panel, the shaded area indicates the ratio derived from observations towards H-MM1, where the fractional pNH3 abundance has been determined by Pineda et al. (2022). In the bottom panel, the “observed” ratios include both H-MM1 and Oph D.

In the text

	Fig. 11 Fractional abundances of the ortho and para modifications of NH3, NH2D, and NHD2 as functions of radial distance from the centre (left) and density (right) in the one-dimensional cloud model used for chemical calculations. The chemistry model assumes constant CD with PH. The abundances are taken at the simulation time 5 × 105 yr.
In the text

	Fig. A.1 NH2D spectra obtained with the LAsMA array towards H-MM1. Simulated spectra assuming constant fractional abundances in the three-dimensional core model are shown with red. The assumed ortho- and para-NH2D abundances are those given in the left column of Table 3, that is, X(oNH2D) = 4.520 × 10−9, X(pNH2D) = 1.469 × 10−9.
In the text

	Fig. A.2 Observed and simulated NHD2 spectra towards H-MM1. The model abundances X(oNHD2) = 9.25 × 10−10, X(pNHD2) = 4.55 × 10−10.
In the text

	Fig. A.3 NH2D spectra observed towards Oph D in the Spring and Summer 2022. The model abundances X(oNH2D) = 4.640 × 10−9, X(pNH2D) = 1.469 × 10−9.
In the text

	Fig. A.4 NH2D spectra towards Oph D observed in September 2022, when the orientation of the LAsMA array was different from that used in the Spring and Summer. The model abundances are the same as in Fig. A.3.
In the text

	Fig. A.5 NHD2 spectra observed towards Oph D in the Spring and Summer 2022. The model abundances X(oNHD2) = 1. 122 × 10−9, X(pNH2D) = 5.02 × 10−10.
In the text

	Fig. A.6 NHD2 spectra observed towards Oph D in September 2022. The model abundances are the same as in Fig. A.5.
In the text

	Fig. B.1 Spin-state branching ratio diagrams for the reaction sequence ${\text{NH}}_3\stackrel{{\text{H}}_2{\text{D}}^{+}}{\to }{\text{NH}}_3{\text{D}}^{+}\stackrel{{\text{e}}^{-}}{\to }{\text{NH}}_2\text{D}$ in the PH (left) and FS (right) models. Assuming that o/p-NH3 = 1, the PH scenario leads to o/p-NH2D = 3, whereas the FS scenario results in o/p-NH2D = 27/13.
In the text

	Fig. C.1 Simulated NH2D and NHD2 spectra towards the centre of the three-dimensional model for H-MM1 with abundance distributions interpolated from the chemistry model assuming constant CD and full scrambling (FS, red lines). The spectra are convolved with the APEx beam. Spectra observed with APEx towards the centre of H-MM1 are shown with blue lines.
In the text

	Fig. C.2 Same as Fig. C.1 but for a chemistry model assuming constant CD and proton hop (PH).
In the text

	Fig. C.3 Same as Fig. C.1 but for a chemistry model assuming dynamic CD and full scrambling (FS).
In the text

	Fig. C.4 Same as Fig. C.1 but for a chemistry model assuming dynamic CD and proton hop (PH).
In the text