© O. Sipilä et al. 2019

1 Introduction

Ion-molecule proton-donation reactions are at the heart of low-temperature chemistry in the interstellar medium (ISM; Herbst & Klemperer 1973). One particular molecule that is formed through such reactions is ammonia, which is an important diagnostic tool for probing the gas temperature in the ISM, and has been extensively used for this purpose (e.g., Friesen et al. 2017, and references therein). Ammonia is a relatively simple molecule and its main gas-phase formation and destruction pathways have long since been identified: its formation begins with the slightly endothermic hydrogen atom abstraction reaction $${\text{N}}^{\text{+}}{\text{+H}}_{\text{2}}\to {\text{NH}}^{\text{+}}\text{+H},$$(1)

followed by the analogous reactions $$\begin{array}{ll}{\text{NH}}^{\text{+}}{\text{+H}}_{\text{2}}\hfill & \to {\text{NH}}_{\text{2}}^{\text{+}}\text{+H}\hfill \\ {\text{NH}}_{\text{2}}^{\text{+}}{\text{+H}}_{\text{2}}\hfill & \to {\text{NH}}_{\text{3}}^{\text{+}}\text{+H}\hfill \\ {\text{NH}}_{\text{3}}^{\text{+}}{\text{+H}}_{\text{2}}\hfill & \to {\text{NH}}_{\text{4}}^{\text{+}}\text{+H}.\hfill \end{array}$$

The dissociative recombination of ${\text{NH}}_{\text{4}}^{\text{+}}$ with free electrons then forms ammonia. The most important destruction pathway (${\text{H}}_{\text{3}}^{\text{+}}+{\text{NH}}_{\text{3}}\to {\text{NH}}_{\text{4}}^{\text{+}}{\text{+H}}_{\text{2}}$) actually leads back to the ammonia formation ladder.

Molecules with multiple hydrogen (or deuterium) atoms can exist in two or several forms owing to the nuclear spin of the proton (deuteron). For examplefor H₂, the combination of two protons gives rise to two spin states where the nuclear spin wavefunction is either antisymmetric (para-H₂; rotational quantum number J = 0, 2, 4...) or symmetric (ortho-H₂; J = 1, 3, 5...); these states behave as distinct chemical species that in practice cannot be converted from one form to the other by radiation or by inelastic collisions. The ~170 K energy difference between the J = 0 and 1 rotational levels is an important source of chemical energy in dense clouds, and allows in particular the exothermic reaction $${\text{H}}_{\text{3}}^{\text{+}}\text{+HD}\to {\text{H}}_{\text{2}}{\text{D}}^{\text{+}}{\text{+H}}_{\text{2}}\text{+232}\text{K}$$(5)

to proceed also in the backward reaction near 10 K. It is known that H₂ D⁺ is the main driver of deuteration in cold clouds (e.g., Rodgers & Charnley 2001). Spin-state and deuterium chemistry are thus very closely tied, and to understand observations of deuterated species through chemical modeling, one needs to take spin-state chemistry into consideration. Several such chemical models have been introduced in the past decades, describing gas-phase or gas-grain chemistry (e.g., Pagani et al. 1992; Walmsley et al. 2004; Sipilä et al. 2010; Furuya et al. 2015; Roueff et al. 2015; Hily-Blant et al. 2018).

Ammonia can exist in a total of nine variants when one takes deuterium and spin states into account. Almost all of these exhibit emission lines that are observable from the ground. Deuterated ammonia traces the innermost areas of starless cores, where many otherwise common molecules such as CO are frozen onto dust grains, as demonstrated in the case of H-MM1 by Harju et al. (2017a; hereafter H17) who carried out an extensive modeling exercise in order to reproduce the observed D/H and spin-state abundance ratios. Their model however failed to reproduce in particular the ortho-to-para (hereafter o/p) ratios of NH₂ D and NHD₂, and it remains unclear why the model prediction deviates from the observed values.

Harju et al. (2017a) used in their study the gas-grain chemical model of Sipilä et al. (2015) that includes a self-consistent description of spin-state chemistry for multiply-deuterated species derived using group theory. The underlying assumption in that model is that the relevant chemical reactions proceed through full scrambling (FS), where it is believed that the reactants form a relatively long-lived intermediate complex in which several atom interchanges can take place (Gerlich & Horning 1992). Full scrambling actually consists of three possible outcomes (Oka 2004): identity, proton hop (PH), and proton exchange. It has been shown that the reactions belonging to the ${\text{H}}_{\text{3}}^{\text{+}}{\text{+H}}_{\text{2}}$ system proceed preferentially through FS even at low temperature (Hugo et al. 2009; Suleimanov et al. 2018). It is however not at all clear that other reacting systems involving more massive molecules also prefer FS over PH at low temperature; Le Gal et al. (2017) have shown that the reactions relevant to H₂ Cl⁺ formation indiffuse cloud conditions do in fact prefer PH, but overall there is unfortunately very little data on this topic in the literature.

In this paper we investigate the effect of proton exchange on ion-molecule chemistry in the ISM by modifying the chemical model of Sipilä et al. (2015), allowing proton-donation reactions, in which proton exchanges were previously allowed, to proceed instead through PH only. From here on, we simply refer to these two types of model either as the FS or PH model, respectively. We concentrate our analysis on deuterated ammonia and reproduce some of the modeling results of H17, using our gas-grain chemical model in combination with radiative transfer simulations of ammonia emission lines. We also investigate the effect of the reaction mechanism on other species besides ammonia, with the aim of quantifying whether FS or PH is preferred. A small part of the results presented here have already been introduced in Caselli et al. (2019).

The paper is organized as follows. In Sect. 2, we introduce the new FS and PH chemical networks, and discuss the details of our modeling. Section 3 presents the results of our analysis, which are further discussed in Sect. 4. We give our conclusions in Sect. 5.

2 Model

2.1 Full scrambling versus proton hop: deuteration

We consider in the FS model that reactions proceed through an intermediate complex that constitutes a pool of atoms from which the atoms are drawn one by one. The branching ratios are constructed by multiplying simple probabilities. For example, in the reaction ${\text{NH}}_{\text{3}}{\text{+D}}_{\text{2}}{\text{H}}^{\text{+}}\to {({\text{NH}}_{\text{4}}{\text{D}}_{\text{2}}^{\text{+}})}^{\text{*}}\to {\text{NH}}_{\text{3}}{\text{D}}^{\text{+}}\text{+HD}$, the intermediate complex contains four H atoms and two D atoms. The four H/D atoms required to form NH₃ D⁺ can be picked in $\left(\begin{array}{c}4\\ 1\end{array}\right)=4$ equivalent combinations: HHHD, HHDH, HDHH, and DHHH. The individual probabilities for each combination are the same, for example for HHHD: $\frac{4}{6}\frac{3}{5}\frac{2}{4}\frac{2}{3}=\frac{2}{15}$. So, in total, the branching ratio for NH₃D⁺ formation is $\left(\begin{array}{c}4\\ 1\end{array}\right)\frac{2}{15}=\frac{8}{15}$. There is only one combination that forms ${\text{NH}}_{\text{4}}^{\text{+}}$ (probability $\frac{1}{15}$), and so the probability to form ${\text{NH}}_{\text{2}}{\text{D}}_{\text{2}}^{\text{+}}$ is $\frac{6}{15}$.

The PH mechanism instead assumes that atom exchanges do not occur during the reaction, and the donating ion simply transfers one proton or deuteron to the other reactant; atom-exchanging reaction channels are closed off. Therefore, the probability of forming NH₃ D⁺ in NH₃ + D₂H⁺ is simply $\frac{2}{3}$, for example. We obtain in general fewer product pathways when assuming PH instead of FS.

Figure 1 illustrates the difference between the models. Here we show the branching ratios arising from the PH and FS assumptions, when the NH₃ + D₂H⁺ reaction first forms ${\text{NH}}_{\text{4}}^{\text{+}}$ (or one of its isotopologs), which subsequently dissociates mainly into NH₃ + H (Öjekull et al. 2004) in a dissociative recombination reaction with an electron. The PH reaction scheme is straightforward. The FS mechanism however, can result in the formation of ${\text{NH}}_{\text{2}}{\text{D}}_{\text{2}}^{\text{+}}$, which is not possible in the PH mechanism, and ${\text{NH}}_{\text{2}}{\text{D}}_{\text{2}}^{\text{+}}$ can then decay into NHD₂, which is not produced in the PH scheme either in these particular reactions. The FS model strongly favors the production of multiply deuterated species.

The PH mechanism as described above can be applied to reactions where one proton is transferred between the reactants, such as proton-donation reactions of the form XH⁺ + Y→X + YH⁺, or hydrogen abstraction reactions such as the main ammonia formation reactions (1)–(4). Fiducially, we applied the PH mechanism to ion-molecule proton-donation reactions only, but we also considered a case where abstraction reactions are modified in similar fashion (see Sect. 4.5). For nearly all other reactions we assume that FS is valid, where applicable. A notable exception is charge-transfer reactions, where the H and D atoms cannot be mixed during the reaction.

We notethat the approach described above to calculate the branching ratios in reactions involving deuterium is extremely simplistic, and is based purely on statistical considerations. It neglects, for example, any zero-point energy effects that come into play when introducing different counts of hydrogen or deuterium in the reactions. Detailed experimental and/or theoretical data are available only for a limited set of reacting systems including deuterated species. An example of this is the study of the ${\text{H}}_{\text{3}}^{\text{+}}{\text{+H}}_{\text{2}}$ system by Hugo et al. (2009), whose data we use here as well.

2.2 Full scrambling versus proton hop: spin states

Assuming PH instead of FS also changes the way that spin-state chemistry is treated in proton-donation reactions. This point has been discussed for example by Oka (2004) in the context of hydrogen. Here we extend the discussion to deuterium. Consider as an example thedeuteron-donation reaction, $${\text{NH}}_{\text{2}}{\text{D}}_{\text{2}}{}^{+}+{\text{D}}_{\text{2}}\to {\text{NH}}_{\text{2}}\text{D}+{\text{D}}_{\text{3}}^{\text{+}}.$$(6)

We notethat this reaction is endothermic and hence inefficient at low temperature. It was selected to illustrate the effect of the adopted reaction mechanism on product spin states when multiple deuterons are present in both reactants. In the FS scenario, the reaction proceeds through a complex that contains in this case two protons and four deuterons. The nuclear spin symmetries of protons ($I=\frac{1}{2}$) and deuterons (I = 1) must be treated separately. The treatment of H₂ spin is trivial: the spin state is unaltered in this reaction. For deuterium we consider the subreaction¹ ${\text{D}}_{\text{2}}^{\text{+}}{\text{+D}}_{\text{2}}\to {\left({\text{D}}_{\text{4}}^{\text{+}}\right)}^{\text{*}}\to {\text{D+D}}_{\text{3}}^{\text{+}}$. The reactions between different spin states in this subreaction can be calculated using the group-theoretical method presented in detail in Sipilä et al. (2015). In brief, the governing principle is the conservation of nuclear spin, and the task is to determine the symmetry representations of the reactants, products, and the intermediate complex, and to derive correlation tables between the group of the complex and its subgroups which represent the reactants and products. The correlation tables yield all of the necessary information, that is, which product spin species are possible given a particular pair of reactants. Example tables are given in Sipilä et al. (2015) and we omit them here for brevity; see also Hily-Blant et al. (2018) for systems not covered in Sipilä et al. (2015).

In the PH scenario, reaction (6) is divided into two subreactions that take place in succession. First, the deuteron-donating ion (${\text{D}}_{\text{2}}^{\text{+}}$) breaks up: $${\text{D}}_{\text{2}}^{\text{+}}\to {\text{D}}^{\text{+}}\text{+D},$$

and in the second step, the deuterium nucleus reacts with the neutral reactant: $${\text{D}}^{\text{+}}{\text{+D}}_{\text{2}}\to {\text{D}}_{\text{3}}^{\text{+}}.$$

Therefore, reaction (6) can be formally broken down to $${\text{NH}}_{\text{2}}{\text{DD}}^{\text{+}}{\text{+D}}_{\text{2}}\to {\text{NH}}_{\text{2}}{\text{D+D}}^{\text{+}}{\text{+D}}_{\text{2}}\to {\text{NH}}_{\text{2}}{\text{D+D}}_{\text{2}}{\text{D}}^{\text{+}}\text{},$$(7)

where the colors are added to help track the placements of the individual D atoms in the reaction. The branching ratios for each of the two subreactions are calculated using the method outlined above, that is, full scrambling applies for the subreactions separately, but not for the reaction as a whole since it proceeds in two distinct steps. The resulting branches for reaction (6) are shown in Fig. 2, omitting the H₂ spin state. We note that the reaction channel that forms p${\text{D}}_{\text{3}}^{\text{+}}$ is forbidden by selection rules when the neutral reactant is oD₂. There are again fewer reaction pathways in the PH model than in the FS model. For the reaction involving oD₂, the two schemes yield the same product channels with identical branching ratios, but when pD₂ is involved, the formation of m${\text{D}}_{\text{3}}^{\text{+}}$ is forbidden in the PH model.

Ammonia formation is governed mainly by hydrogen atom abstraction reactions of the type ${\text{NH}}_i^{+}+{\text{H}}_{\text{2}}\to {\text{NH}}_{i+1}^{+}+\text{H}$, such as the N⁺ + H₂→NH⁺ + H reaction, and it is clear that considering these reactions as abstraction reactions – as opposed to assuming them to proceed through scrambling as is our standard approach – would have an effect on our results. We discuss this topic in detail in Sect. 4.5.

Fig. 1 Branching ratio diagrams for the NH3 + D2H+ reaction in the FS (left) and PH (right) models, showing the re-distribution of hydrogen and deuterium in the reaction followed by dissociative recombination with electrons.

Fig. 2 Spin-state branching ratio diagrams for the ${\text{NH}}_{\text{2}}{\text{oD}}_{\text{2}}^{\text{+}}\text{+}(\text{o/p}){\text{D}}_{\text{2}}$ reaction in the FS (left) and PH (right) models, omitting the H2 spin state.

2.3 Chemical model

We use the gas-grain chemical code described in Sipilä et al. (2019). In short, the code uses the rate-equation approach to calculate chemical time-evolution in the gas phase and on grain surfaces. Gas-phase and grain-surface chemistry are connected via adsorption and (non-thermal) desorption processes. We consider a two-phase model where the ice on dust grains consists of asingle reactive bulk, that is, no separation between a chemically active layer and an inert mantle is made.

As in our recent modeling papers (e.g., Sipilä et al. 2019), we took the latest public release of the KIDA gas-phase network (kida.uva.2014; Wakelam et al. 2015) as our basis network, which was deuterated as described in Sect. 2.1. Spin-state chemistry was then introduced in the deuterated network following either the FS or PH methods discussed above. As in Sipilä et al. (2019), we explicitly track the spin states of H₂, ${\text{H}}_{\text{2}}^{\text{+}}$, ${\text{H}}_{\text{3}}^{\text{+}}$, and their deuterated isotopologs, as well as the spin states of any multiply protonated or deuterated species involved in the water and ammonia formation networks. We introduced deuteration and spin states in our grain-surface network (Sipilä et al. 2019) in an analogous way, although we assume that no ions are present on or in the ice, and consider only the FS mechanism for those surface reactions where atom interchanges occur. The rate coefficients of electron recombination reactions involving the protonated ions in the main ammonia formation pathway are taken from the KIDA database. The products of the dissociative pathways follow statistical ratios.

The final number of gas-phase reactions depends on whether we adopt the FS or the PH model. In the former case, the network contains ~75 000 reactions, while in the latter case the number is reduced to ~64 000. The number of grain-surface reactions is similar in both cases: ~2400 versus ~2200 in the FS and PH models, respectively.

2.4 Physical models

In this paper we demonstrate the differences arising from using either the FS or PH mechanism on the D/H and spin-state abundance ratios of various chemical species. We use zero-dimensional (0D) physical models for this purpose, that is, we calculate the chemical evolution as a function of time in point-like fixed physical conditions. We adopt in the 0D models a constant temperature of T_dust = T_gas = 10 K, a visual extinction of A_V = 10 mag, and a primary cosmic ray ionization rate per hydrogen atom of ζ_p(H) = 1.3 × 10⁻¹⁷ s⁻¹.

We also reproduce some results presented in H17, where we attempted to fit observations of deuterated ammonia in the starless core H-MM1 using a previous version of our chemical model. The model results presented therein were derived using the FS assumption. Here we compare the FS and PH mechanisms to see which one better matches the observations. We therefore used the core model derived in H17 and divided it into concentric spherical shells in which we calculated the chemical evolution. We then combined the results to produce radius-dependent abundance profiles, and compared them to the observed tendencies. The density and temperature profiles of the model core are displayed in Fig. 3. We refer the reader to H17 for details on how these profiles were derived, but highlight here two prominent features of the gas temperature profile: (1) The profile is based on NH₃ (1,1) and (2,2) observations. The sharp increase in the temperature at around 7500 au is accompanied with a transition from subsonic to supersonic turbulence regime; (2) The gas temperature lies below the dust temperature outside the inner region where gas-dust coupling is important, which is probably due to efficient line cooling (see e.g., Sipilä & Caselli 2018). The initial abundances, reproduced in Table 1, correspond to the EA1 set of Wakelam & Herbst (2008), except that the elemental nitrogen abundance was multiplied by a factor of 2.5 (H17). We do not consider chemical desorption because it leads in our two-phase model to ammonia abundances that are much higher than observed (Sipilä et al. 2019). The effect of modifying the initial H₂ o/p ratio is discussed in Sect. 3.2.

We also carried out radiative transfer simulations of the ammonia lines in order to compare our chemical modeling results with the observations of H17 in a meaningful way. These simulations are described in Sect. 4.2.

Fig. 3 Density and temperature profiles of the H-MM1 core model.

3 Results

3.1 Zero-dimensional models

Figures 4 and 5 show the time-evolution of various D/H and spin-state ratios at different medium densities in the FS and PH models at T = 10 K. First, it is evident that ${\text{H}}_{\text{3}}^{\text{+}}$ and its isotopologs are practically unaffected by the branching mechanism. This is because the chemistry of these species is mainly determined by the ${\text{H}}_{\text{3}}^{\text{+}}{\text{+H}}_{\text{2}}$ reacting system, for which we adopt the rate coefficients calculated by Hugo et al. (2009) regardless of what we assume for any other proton-donation reactions. The feedback of reactions such as those illustrated in Figs. 1 and 2 on (deuterated)${\text{H}}_{\text{3}}^{\text{+}}$ is evidently very small.

Both the D/H and spin ratios of ammonia are in turn clearly affected by the branching mechanism, although the spin ratios experience a lesser effect. Adopting the PH mechanism increases the NH₂D∕NH₃ ratio, but decreases the NHD₂∕NH₂D and ND₃∕NHD₂ ratios by about the same factor. While we recover similar time-dependence in the D/H ratios as in our previous work (Sipilä et al. 2015), we obtain generally higher values for the ratios, and notably the NH₂D∕NH₃ ratio can exceed unity in the current model. There is an evident peak for the ratios which depends on time; this occurs at ~ 2 × 10⁵ yr at n(H₂) = 10⁵ cm⁻³, for example. The difference between the FS and PH models increases at late times for the NH₂D∕NH₃ ratio, but decreases for NHD₂∕NH₂D and ND₃∕NHD₂. The ammonia spin ratios show some density dependence in the difference between the FS and PH models, although this is always a factor of two at most. The ND₃ spin ratios are affected the most when changing the model. This is expected, given that ND₃ has three distinct spin states, and the PH model includes significantly less reaction channels for ND₃ formation. Statistical ammonia D/H ratios obey (NH₂D∕NH₃)∕(NHD₂∕NH₂D) = 3 and (NHD₂∕NH₂D)∕(ND₃∕NHD₂) = 3 (Brown & Millar 1989; H17). The PH model predicts for these ratios (at late times at high density) values of ~2.9 and ~2.2, respectively. The corresponding values in the FS model are ~1.45 and ~2.3. We discuss these ratios in more detail in Sect. 4.5.

The water D/H ratios present similar behavior to that of the ammonia ones: the PH model predicts more deuteration in singly substituted molecules, but less deuteration in multiply substituted molecules. The difference in the HDO/H₂O and D₂ O/HDO ratios between the two models is almost density-independent. The H₂O o/p ratio is largely unaffected by the choice of branching mechanism. The D₂ O o/p ratio presents small deviations (of the order of a few tens of percent; note the y-axis scale in Fig. 5) depending on the model, which are attributed to the spin ratios of its precursors HD₂O⁺ and D₃O⁺ (in particular the most abundant meta (m) and ortho forms). For example, mD₃O⁺ + e⁻ can only lead to oD₂ O formation while oD₃ O⁺ + e⁻ can form ortho and para D₂ O with equal probability. At high density in the FS model, the D₃O⁺ m/o ratio is above unity while in the PH model it is below unity, meaning that production of para-D₂O is more efficient in the PH model, decreasing the D₂O o/p ratio with respect to the FS model.

The general conclusion from these results is that the PH model favors singly deuterated species over multiply deuterated ones and produces deuterated species slightly more efficiently in the higher-energy spin states, that is, pNH₂D, pNHD₂, and oND₃, when compared to the FS model. (The energy level diagrams of the various ammonia isotopologs are shown in Fig. 3 in H17.) The former conclusion was already expected based on our discussion in Sect. 2, where the example reaction showed the tendency of the PH model to omit multiply deuterated pathways.

Fig. 4 D/H abundance ratios of ammonia (left column), ${\text{H}}_{\text{3}}^{\text{+}}$ (middle column), and water (right column), at a medium density of n(H2) = 104 cm−3 (top row), n(H2) = 105 cm−3 (middle row), and n(H2) = 106 cm−3 (bottom row). The abundances represent sums over the spin states where applicable. Solid lines represent the FS model, while dashed lines represent the PH model.

Fig. 5 Spin-state abundance ratios of ammonia (left column), ${\text{H}}_{\text{3}}^{\text{+}}$ (middle column), and water (right column), including their respective deuterated forms, at a medium density of n(H2) = 104 cm−3 (top row), n(H2) = 105 cm−3 (middle row), and n(H2) = 106 cm−3 (bottom row). Solid lines represent the FS model, while dashed lines represent the PH model. The statistical ammonia spin-state abundance ratios are indicated with colored arrows.

3.2 The H-MM1 core model

Figure 6 shows the D/H and spin ratios of the same species as in Figs. 4 and 5, but as radial profiles calculated using the H-MM1 coremodel. The results correspond to t = 10⁵ yr since the start of the simulation. (The time was arbitrarily chosen for the sake of illustration.) The same trends already discussed above are present here as well: for example, the NHD₂∕NH₂D and ND₃∕NHD₂ ratios are clearly lower in the core center in the PH model. The sharp decrease in the H₂ D⁺ o/p ratio near the edge of the core is due to the sudden increase in the gas temperature (Fig. 3). A corresponding feature is present in some of the other spin ratios as well, although it is far less prominent for the other ratios.

The D/H and spin ratios are heavily time-dependent, and it is difficult to draw any firm conclusions on the difference between the model and observations based on the radial profiles presented in Fig. 6. To complement that figure, we show in Fig. 7 the density-weighted abundances and abundance ratios of the ammonia isotopologs as functions of time. The observed D/H and spin ratios derived by H17 are overlaid. The early-time behavior of the various ammonia abundances is in the present model very different from that of H17. This is because of several updates to the chemical model since Sipilä et al. (2015) that affect the early-time (due to ion-molecule reactions) and late-time (due to cosmic rays) chemistry in particular (see Sipilä et al. 2019). However, we recover similar abundances compared to H17 in the time interval 10⁵ –10⁶ yr; the current model predicts a factor of a few more ND₃.

The D/H ratios again follow the same trends as in the discussion above – for example the NH₂D∕NH₃ ratio is enhanced in the PH model. Notably, the FS model can only fit the NH₂D∕NH₃ ratio for realistic timescales (≳10⁵ yr) where the abundances are not below detectable levels. The PH model reproduces the observed NHD₂∕NH₂D ratio withinthe observational uncertainties, and simultaneously fits the observed ND₃∕NHD₂ ratio within a factor of two. While the NH₂D∕NH₃ ratio is better reproduced by the FS model, we note that H17 only detected para-NH₃ and derived the NH₂D∕NH₃ ratio using an NH₃ o/p ratio of unity. If we instead assume NH₃ o/p = 0.5 as suggested by our chemical model, the PH model can match the observed ratio (orange hatched box) relatively well.

To summarize, the PH model can reproduce all three D/H ratios on a similar timescale, whereas the FS model can only reproduce the NH₂D∕NH₃ ratio. We note that H17 derived the observed abundance ratios assuming constant abundances throughout the core, and this is in stark contrast with the chemical model which predicts strong abundance gradients. Indeed, neither the FS nor the PH model can match the observed lines well, as we show in Sect. 4.2.

The observed NHD₂ o/p ratio is reproduced by both models, and both underestimate the NH₂D o/p ratio by about 50%. The similarities in the spin ratios of the ammonia isotopologs in the two models is unsurprising given that the ratios are dominated by reactions with ${\text{H}}_{\text{3}}^{\text{+}}$ isotopologs (see the discussion in Sipilä et al. 2015), and in this paper we use the same chemistry for the ${\text{H}}_{\text{3}}^{\text{+}}{\text{+H}}_{\text{2}}$ system (Hugo et al. 2009) in both the FS and PH cases. Hence a modification to the ${\text{H}}_{\text{3}}^{\text{+}}{\text{+H}}_{\text{2}}$ system would be required to obtain a better agreement with the observed ammonia spin state ratios. This is beyond the scope of the present work; see Sects. 4.3 and 4.5 for some additional discussion.

We adopted in our models a low initial H₂ o/p ratio of 10⁻³, assuming that the ratio has been thermalized to 20 K before the dense core stage (H17). Low-temperature deuterium chemistry is very sensitive to the H₂ o/p ratio, and hence it is clear that changing the initial value would influence our results, at least at early times in the simulation. To test this, we ran the FS and PH models assuming instead an initial H₂ o/p ratio of 0.1. This modification leads to very low abundances of deuterated species at early times (before ~ 10⁵ yr) because the presence of ortho-H₂ suppresses deuteration, and moves the abundance peaks displayed in Fig. 7 forward by a factor of approximately two. The results of these test models and those of our fiducial FS and PH models are very close to each other for t ≳ a few × 10⁵ yr, and in fact at late times the test models present very slightly higher D/H ratios (of the order of a few percent) for the various species than the corresponding fiducial models do.

Fig. 6 D/H ratios (top row) of ammonia (left column), ${\text{H}}_{\text{3}}^{\text{+}}$ (middle column), and water (right column) and the corresponding spin-state ratios (bottom row), as functions of radius in the H-MM1 core model. The curves correspond to t = 105 yr. Solid lines represent the FS model, while dashed lines represent the PH model.

Fig. 7 Abundances of the ammonia isotopologs summed over the spin states (left), corresponding D/H ratios (middle), and spin-state ratios (right), as functions of time in the H-MM1 core model. Solid lines represent the FS model, while dashed lines represent the PH model. Observational ratios derived by H17 are overlaid (hatched boxes). The orange hatched box in the middle panel represents the NH2 D∕NH3 ratio calculated using an NH3 o/p ratio of 0.5 (see text for details).

4 Discussion

4.1 Variation in the physical core model

The results discussed above suggest that the PH model is better than the FS model in reproducing the ammonia D/H ratios observed toward H-MM1. This result holds within the assumption that our physical model for H-MM1 accurately represents the actual core. The core model is based on observationally derived H₂ column density as well as gas and dust temperature profiles (see H17 for details). The observational error in the H₂ column density is less than a factor of two, but for the gas temperature the error ranges from less than one kelvin to several kelvin depending on the distance from the center of the core, which can affect the results of our modeling. An important question in the present context is: can the uncertainty in the core model affect the goodness of the fit in the PH model versus FS model; in other words, is it possible that for another physical model the PH model would no longer be the better match to the observations? Our 0D calculations show that the difference between the two models does not present strong density variations, but to further quantify this issue we ran additional core models, varying the density and gas temperature profiles. In these tests our fiducial density profile was scaled by a factor of 0.5, 1, or 2, and the fiducial gas temperature profile was shifted by −1, 0, or 1 K for a total of nine models including our fiducial model. We note that although the observational error in the gas temperature may be several kelvin, in the inner core where most of the line emission originates (depending on the species and the line) the error is of the order of one kelvin only (H17), and hence shifting the gas temperature by this conservative value is justified. We do not alter the dust temperature profile for the sake of simplicity.

Figure 8 shows the results of these calculations, with all nine cases overlaid. There is some spread in the ratios owing to the changes in the physical conditions, but an analysis of the individual curves shows that the PH model always yields a higher NH₂D∕NH₃ ratio, and lower NHD₂∕NH₂D and ND₃∕NHD₂ ratios, than the FS model does, no matter which parameter combination we use. This strengthens the view that the PH model is indeed the better match to the observed abundance ratios.

The observed NH₂D∕NH₃ ratio (assuming NH₃ o∕p = 0.5) can be reached by the PH model if the density of the model core is multiplied by a factor of 0.5, and the best match is obtained if the fiducial temperature profile is either unaltered or shifted down by one kelvin. The ND₃∕NHD₂ ratio is also best fit by these parameter combinations; in those cases the NHD₂∕NH₂D ratio is just at the observed lower limit at t = 10⁶ yr. Regardless of the changes to the density and temperature, we never obtain an NH₂D o/p ratio higher than approximately two in our models, marking a clear discrepancy with the observations.

4.2 Line emission profiles

The match between the models and observations is best quantified by simulating line emission and comparing it with the observed lines, so that possible optical depth and excitation effects can be taken into account. Harju et al. (2017a) presented a detailed analysis of ammonia in H-MM1 including the lines of several other species, such as H₂ D⁺ and N₂ D⁺, and searched for the best match between models and observations by tuning the chemical and physical model parameters (e.g., interstellar radiation field strength, diffusion-to-binding energy ratio on grain surfaces, etc.). In this paper we perform only a reduced version of their analysis, concentrating on the ammonia lines. We adopt the radiative transfer pipeline used by H17 and apply it to the results from the up-to-date chemical models presented above.

We extracted the abundance profiles of the eight distinct (deuterated) ammonia isotopologs (taking into account spin states) observed by H17 at 51 time steps spaced between 10⁵ and 10⁷ yr, and carried out line simulations for all of the eight species at each time step using the non-local-thermal-equilibrium Monte Carlo radiative transfer code of Juvela (1997). We then searched for the best match between the observations and models by performing Pearson’s χ² -test on the simulated lines² $${\chi }^2={\displaystyle \sum _i}\frac{{\left(T_{\text{MB}}^{\text{obs}}(i)-T_{\text{MB}}^{\text{mod}}(i)\right)}^2}{T_{\text{MB}}^{\text{mod}}(i)},$$(8)

where the summation is over spectral channels. The final χ² value at each time step is a simple mathematical average of seven species; we exclude para-ND₃ from the test because it was not detected by H17. The best-fit time that we obtain is t = 3.29 × 10⁵ yr for both the FS and the PH model. The various modeled D/H ratios are closest to the observed ratios at around t = 1 × 10⁶ yr, but at these latetimes the abundances of the deuterated ammonia isotopologs are too low (Fig. 7) to match the intensities of the observed lines.

Figure 9 shows the best-fit simulated lines along with the lines observed by H17. We find that: (1) the para-NH₂D lines are overestimated by the model, which is unsuprising given that the model predicts a lower NH₂D o/p ratio than what is deduced from the observations; (2) the NHD₂ and ND₃ lines are clearly overestimated by the current chemical model, but the PH model reproduces the relative strength of the ortho and para lines slightly better; and (3) the PH model is overall the better match to the observations, although the difference is not great. The line simulation results confirm the overall tendencies noted for the abundance distributions discussed above.

The agreement with the observations obtained by H17 using the chemical model of Sipilä et al. (2015) is better than what we obtain with the present model. In this paper we did not attempt to modify the chemical model parameters for the purposes of obtaining the best possible match to the observations, and the present results do not imply that the new model performs worse. Indeed, if we explored the parameter space in the same way as H17 did, the match between the model and the observationscould be improved. We did carry out another set of calculations using one of our alternative core models where the fiducial density profile is multiplied by 0.5 and the gas temperature profile is shifted down by one kelvin. We found that in this case the overall fit is worse than in the regular FS or PH models, and already at t = 3.29 × 10⁵ yr the NH₂D and NHD₂ lines are too weak. The results highlight the fact that line emission simulations are essential for comparing observations and chemical models.

Fig. 8 Left: D/H and spin-state abundance ratios in the PH model, calculated using nine different physical core models as explainedin the text. Right: corresponding data from the FS model. Hatched boxes correspond to observed ranges as in Fig. 7.

4.3 Effect of species-to-species rate coefficients

In Sipilä et al. (2017) we presented a new set of rate coefficients for the ${\text{H}}_{\text{3}}^{\text{+}}{\text{+H}}_{\text{2}}$ reacting system, taking rotational excitation into account in the derivation of the rate coefficients (see Hily-Blant et al. 2018 for a recent similar work). We found that the inclusion of excited states leads in general to reduced deuteration levels. This effect is prominent for densities above 10⁶ cm⁻³ and temperatures above 10 K, but we also found that the ortho/para ratio of H₂ D⁺ is affected already at a density of a few times 10⁵ cm⁻³ because of the critical density of the ground-state 1₁₀ –1₁₁ rotational transition of oH₂D⁺. Given the dependency of the ammonia o/p (and D/H) ratios on the spin states of the ${\text{H}}_{\text{3}}^{\text{+}}$ isotopologs, it is conceivable that rotational excitation could also modify the ammonia ratios indirectly. We tested this issue and found no significant effect on our results apart from a decrease in the H₂ D⁺ o/p ratio of a few tens of percent across the core. The density and temperature of H-MM1 are too low for rotational excitation of the ${\text{H}}_{\text{3}}^{\text{+}}$ isotopologs to affect ammonia chemistry. This is in line with the discussion in Sipilä et al. (2017) where we pointed out that the new rates are more important in protostellar disks than in starless or pre-stellar cores. However, we do obtain a noticeable effect on the H₂ D⁺ o/p ratio which indicates that including rotational excitation is important for the correct interpretation of H₂ D⁺ line emission in starless and pre-stellar cores, as already clearly demonstrated in Harju et al. (2017b).

4.4 Abundances of additional species

We concentrated mainly on ammonia in the above analysis because of the direct comparison with our previous observations toward H-MM1. However, the nature of the proton-transfer mechanism has the potential to affect a wide variety of species other than ammonia, either directly or indirectly. In Fig. 10, we show the abundances and D/H ratios N₂ D⁺, DCO⁺, HDCO, D₂ CO, CH₂D, and CHD₂. Of this list of species, N₂D⁺ and DCO⁺ are not affected by the choice of the model. This is not surprising given that both are straightforward derivatives of H₂ D⁺ which is itself virtually unaffected by the transfer mechanism as explained above. Deuterated formaldehyde shows behavior that is different from water for example: both the singly and doubly deuterated variants are produced less efficiently inthe PH model. Formaldehyde is produced (at high density) in the dissociative electron recombination of CH₂OH⁺ with free electrons, and in neutral-neutral reactions of atomic oxygen with CH₃. The fraction of deuterated CH₃ is higher in the FS model than in the PH model – this derives from differences in the chemistry of deuterated ${\text{CH}}_{\text{4}}^{\text{+}}$ – leading to more efficient production of HDCO and D₂CO through the neutral-neutral pathway. We note that (deuterated) H₂CO and CH₃ have spin states that could in principle be tracked using our spin-separation routine, but we leave such an analysis for future work.

Fig. 9 Simulated emission lines for deuterated ammonia (specific species and line indicated in the panels) in the FS model (red) and PH model (blue). The black histograms show the observed lines (H17). The results correspond to the best-fit time t = 3.29 × 105 yr (see text).

Fig. 10 As the left and middle panels of Fig. 7, but showing the abundances and D/H ratios of some additional species.

Fig. 11 As Fig. 7, but showing the averaged abundances in H-MM1 in the FS model (solid lines) and in the PHalt model discussed in the text (dotted lines). The hatched boxes showing the observed ranges in the right panel are omitted here for better visibility of the model results.

4.5 Hydrogen atom abstraction reactions

The analysis presented above shows that neither the FS or the PH model can reproduce the NH₂D o/p abundance ratio observed toward H-MM1 regardless of variations in the physical core model, and that the various modeled D/H ratios also differ from the observed values. Although we did not perform an exhaustive parameter-space exploration, it seems clear that the observed ratios cannot be reproduced by modifying the reaction mechanism in proton-donation reactions only. The main ammonia formation reactions are abstraction reactions, and hence we constructed another chemical model where we assume that all ion-molecule hydrogen abstraction reactions analogous to reactions (1)–(4) proceed through PH instead of FS. The efficiencies of some of these reactions depend on the inclusion of deuterium and/or spin states. For example, reaction (1) proceeds much faster with ortho-H₂ than with para-H₂ (Dislaire et al. 2012), while the substitution of deuterium in H₂ slows down the reaction (Marquette et al. 1988). Reaction (4) is slow at low temperature, although it does increase somewhat in efficiency toward 10 K (Barlow & Dunn 1987). The branching ratios for the modified abstraction reactions are calculated in a way analogous to what is described in Sect. 2 for the proton-donation reactions. For example, the reaction NH₂D⁺ + H₂ can then only lead to NH₃D⁺ + H – the product pair ${\text{NH}}_{\text{4}}^{\text{+}}\text{+D}$ is not accessible like it is in the FS model, meaning that this alternative model (hereafter PH_alt) is expected to produce larger deuterium fractions as many formation pathways leading to hydrogenated species are no longer available, compared to the FS model. Changing the reaction mechanism naturally also leads to repercussions in spin-state ratios. Proton-donation reactions and hydrogen abstraction reactions are both treated under the PH assumption in the PH_alt model.

Figure 11 shows the averaged abundances in H-MM1 calculated with the FS model and the PH_alt model. Switching to a direct hydrogen atom hop mechanism for the abstraction reactions greatly increases the D/H ratios, as expected, and indeed the NH₂D∕NH₃ and NHD₂∕NH₂D ratios are now both higher than the observed ones. PH_alt produces ortho and para ammonia in roughly equal proportions, and in fact the modeled NH₂D∕NH₃ abundance ratio is now over a factor of two higher than the observed value if we account for the NH₃ o/p ratio in the PH_alt model (~1.2). The value of the NH₂D o/p ratio (~2.7) at t ~ 10⁶ yr is close to the observed value of 3.0 ± 0.2. Likewise, the o/p ratio of NHD₂ approaches two at late times. This indicates that the modified abstraction reaction mechanism may give a better representation of ammonia spin-state ratios if the problem of excessively high D/H ratios can be countered.

In Sipilä et al. (2010), we found that increasing the cosmic-ray ionization rate leads to less efficient production of the deuterated ${\text{H}}_{\text{3}}^{\text{+}}$ isotopologs, hence decreasing the efficiency of deuteration in general. Motivated by this result, we tested the influence of the cosmic-ray ionization rate on the H-MM1 model results by varying the ionization rate between ζ_p (H) = 1.0 × 10⁻¹⁷ s⁻¹ and 1.0 × 10⁻¹⁶ s⁻¹. We found that increasing ζ_p(H) does indeed decrease the ammonia D/H ratios, but it also causes the abundances of the ammonia isotopologs to drop below detectable levels in an increasingly shorter timescale, of the order of a few times 10⁵ yr. However, the spin-state abundance ratios of the ammonia isotopologs are insensitive to changes in ζ_p (H). This suggests that it may be possible for the PH_alt model to reproduce the observed emission lines if the ionization rate is increased from our fiducial value of 1.3 × 10⁻¹⁷ s⁻¹. Figure 12 shows as an example the simulated lines for ζ_p(H) = 5.0 × 10⁻¹⁷ s⁻¹; the core model is otherwise identical to that used in Fig. 9. The best-fit time in this case is t = 2.50 × 10⁵ yr. The overall fit is better than the one obtained with the PH model, and in particular the NH₂D o/p ratio is correctly reproduced, but the better agreement comes with a cost. The steep decline of the ammonia abundances at relatively short timescales is incompatible with observations toward the starless core L1544 which show no ammonia depletion even though the core is in an evolved dynamical state (Crapsi et al. 2007; Caselli et al. 2017; Sipilä et al. 2019); it would be rather serendipitous if the observations caught H-MM1 during the short-lived ammonia abundance peak, assuming that the two cores have evolved similarly. This problem may be solved once it is understood how ammonia can maintain a high gas-phase abundance even at high gas densities (see Sipilä et al. 2019 for more discussion on this issue). Figure 12 also shows the simulated lines in the PH_alt model for our fiducial cosmic-ray ionization rate, and it is clear that the agreement with the observations is not better when compared to the PH model, except for the fact that the NH₂D o/p ratio corresponds to the observed value.

As noted in Sect. 3.1, statistical ammonia deuteration ratios obey (NH₂D∕NH₃)∕(NHD₂∕NH₂D) = 3 and (NHD₂∕NH₂D)∕ (ND₃∕NHD₂) = 3. Model predictions of these ratios depend on the deuterated ammonia formation mechanism (Rodgers & Charnley 2002), and thus they may be useful in discriminating between the FS and PH scenarios. Figure 13 shows the time-evolution of these ratios in the H-MM1 core model in the three cases studied above (FS, PH, and PH_alt). We recover the same tendency as for point models (Sect. 3.1), that is, that the (NHD₂∕NH₂D)∕(ND₃∕NHD₂) ratio is higher than the (NH₂D∕NH₃)∕(NHD₂∕NH₂D) ratio in the FS model, and vice versa for the PH model. Switching from the PH to the PH_alt model increases both ratios by a similar factor. These results highlight the strong model-dependent behavior of the (NH₂D∕NH₃)∕(NHD₂∕NH₂D) ratio in particular. The observed values in H-MM1 are (NH₂D∕NH₃)∕(NHD₂∕NH₂D) =1.77 ± 0.19 and (NHD₂∕NH₂D)∕(ND₃∕NHD₂) = 3.67 ± 0.70, as calculated from Table 4 in H17. Neither one of these ratios is reproduced by the models. However, we note that the comparison depends on how the column densities are derived; H17 simply assumed a constant abundance throughout the core when determining the best-fit column density. Comparison of line profiles remains the best method of constraining models with the aid of observations, because this allows to take into account the effects of kinematics and optical depth, for example.

We did not consider any corrections to rate coefficients due to differences in mass or zero-point energy between hydrogenated and deuterated species, and so deuterated reactions proceed with the same rate coefficients as their hydrogenated counterparts. We did test the relative rates for the deuterated ${\text{NH}}_{\text{3}}{}^{+}+{\text{H}}_{\text{2}}\to {\text{NH}}_{\text{4}}^{\text{+}}\text{+H}$ reactions tabulated by Roueff et al. (2015; their Table 3), but this change had a negligible influence on our results, which is readily explained by the fact that collisions of ammonia with HD or D₂ are relatively rare. Unfortunately, in general little theoretical or laboratory data are available for the deuterated reactions involved in the main ammonia formation pathway; our results show that more work on these reactions is urgently needed to reach a satisfactory agreement between chemical models and observations. The fact that the PH_alt model shows such good agreement with the observations for the NH₂D and NHD₂ spin ratios while simultaneously failing to match the D/H ratios (assuming ζ_p (H) = 1.3 × 10⁻¹⁷ s⁻¹) is especially puzzling. The NH₃ o/p ratio has not been measured toward H-MM1 because the low-lying rotational lines (belonging to ortho-NH₃) cannot be observed from the ground, but the ratio that the PH_alt model predicts (~1.2) is clearly higher than the value (~0.7) inferred toward L1544 by Caselli et al. (2017).

It is possible that nonthermal desorption of ammonia from grain surfaces plays a part in setting the abundance ratios in the gas phase. Ammonia is formed on the grain surfaces by sequential association reactions, and in these reactions the products are formed in statistical ratios; for example, NH + H leads to ortho and para NH₂ in a 3:1 ratio. Thus, if ammonia was released in significant quantities from the grain surfaces and was mixed with the nonstatistical gas-phase distribution, the resultant abundance ratios could be close to statistical as observed by H17. Our chemical model however predicts strong ammonia freeze-out, and standard thermal or nonthermal mechanisms do not lead to appreciable ammonia desorption in starless core conditions as we have recently shown in Sipilä et al. (2019). We do not delve into this matter further in this paper, and refer the reader to Sipilä et al. (2019) for a complete discussion including suggested solutions to the problem.

The results presented in this paper highlight the fact that multiply deuterated species are very clearly affected by the nature of the proton-transfer process, and that our modeling of the starless core H-MM1 seems to support the PH mechanism as the favored one, but it is difficult to draw definitive conclusions based on the limited data available. It is possible that the prevalence of PH or FS is dependent on the reacting system. We propose that pending further theoretical and laboratory work, a systematic observational programme targeting pairs of (multiply) deuterated species toward different objects could help to constrain the reaction mechanism, when paired with detailed modeling.

Fig. 12 As Fig. 9, but for the PHalt model assuming either ζp(H) = 1.3 × 10−17 s−1 (blue) or ζp(H) = 5.0 × 10−17 s−1 (red). The corresponding best-fit chemical times are t = 4.33 × 105 yr and t = 2.50 × 105 yr.

Fig. 13 Abundance ratios (NH2D∕NH3)∕(NHD2∕NH2D) and (NHD2∕NH2D)∕(ND3∕NHD2) as functions of time in the H-MM1 core model. The blue, red, and green curves represent the FS model, the PH model, and the PHalt model, respectively.

5 Conclusions

We investigated the effect of modifying the reaction mechanism in gas-phase proton-donation reactions on the abundances of (deuterated) ammonia and other key species in the starless core H-MM1, limiting the reactions to proceed through a direct proton hop and disallowing atom-exchange processes. We quantified the difference between the full scrambling and hop models by investigating density-averaged abundance profiles and by carrying out line emission simulations using a radiative transfer code. We also constructed an alternative hop model that includes the treatment of hydrogen atom abstraction reactions as a direct hop process. Our results apply to hydrogenated as well as deuterated species, and the effect of the reaction mechanism on spin-state chemistry is consistently taken into account.

We found that the proton hop mechanism favors the formation of singly deuterated species over multiply deuterated ones, because in this case many reaction channels that would lead to multiply deuterated species in the FS model are closed off. Switching the reaction mechanism also has implications for spin-state abundance ratios. The PH model reproduces somewhat better the NH₂D∕NH₃, NHD₂∕NH₂D, and ND₃∕NHD₂ ratios observed toward H-MM1 by Harju et al. (2017a) than the FS model does, but in particular the observed o/p ratio of NH₂D is not reproduced by either model – both models predict a ratio of approximately two while the observed value is close to the statistical value of three. We carried out radiative transfer simulations of the observed ammonia emission lines to confirm that these tendencies found in the modeled abundance profiles are present in the lines as well. In the alternative model where also hydrogen abstraction reactions proceed through proton hop, the observed ammonia spin-state abundance ratios are reproduced but the D/H ratios are very strongly overestimated. This problem can be mitigated by considering cosmic-ray ionization ratios higher than our fiducial value of ζ_p (H) = 1.3 × 10⁻¹⁷ s⁻¹, but achieving a good match with the observations would also require fine-tuning of the physical and chemical model parameters, which is beyond the scope of the present work.

Chemical species produced directly through proton-donation reactions, such as N₂ D⁺ and DCO⁺, are unaffected by the nature of the proton-donation reaction mechanism. However, species such as H₂ CO, H₂ O, and CH₃ whose gas-phase formation pathways depend on proton-donation reactions, are strongly affected by the reaction mechanism.

Our models cannot constrain definitively whether the proton hop process prevails over full scrambling (i.e., allowing for proton exchange) at low temperature in the ISM. The results presented here do however present a preference for the former, at least for ammonia. Further theoretical and laboratory work on multiple reacting systems is required to shed more light on this problem. In the meantime, we propose that observations of multiple pairs of (multiply) deuterated species toward different starless and pre-stellar cores could provide constraints for models of low-temperature deuterium and spin-state chemistry.

Acknowledgements

We thank Stephan Schlemmer for insightful comments on an early version of the manuscript, and the anonymous referee for a constructive report that improved the paper.

References

All Tables

All Figures

	Fig. 1 Branching ratio diagrams for the NH3 + D2H+ reaction in the FS (left) and PH (right) models, showing the re-distribution of hydrogen and deuterium in the reaction followed by dissociative recombination with electrons.
In the text

	Fig. 2 Spin-state branching ratio diagrams for the ${\text{NH}}_{\text{2}}{\text{oD}}_{\text{2}}^{\text{+}}\text{+}(\text{o/p}){\text{D}}_{\text{2}}$ reaction in the FS (left) and PH (right) models, omitting the H2 spin state.
In the text

	Fig. 3 Density and temperature profiles of the H-MM1 core model.
In the text

	Fig. 4 D/H abundance ratios of ammonia (left column), ${\text{H}}_{\text{3}}^{\text{+}}$ (middle column), and water (right column), at a medium density of n(H2) = 104 cm−3 (top row), n(H2) = 105 cm−3 (middle row), and n(H2) = 106 cm−3 (bottom row). The abundances represent sums over the spin states where applicable. Solid lines represent the FS model, while dashed lines represent the PH model.
In the text

Fig. 5 Spin-state abundance ratios of ammonia (left column), ${\text{H}}_{\text{3}}^{\text{+}}$ (middle column), and water (right column), including their respective deuterated forms, at a medium density of n(H2) = 104 cm−3 (top row), n(H2) = 105 cm−3 (middle row), and n(H2) = 106 cm−3 (bottom row). Solid lines represent the FS model, while dashed lines represent the PH model. The statistical ammonia spin-state abundance ratios are indicated with colored arrows.

In the text

	Fig. 6 D/H ratios (top row) of ammonia (left column), ${\text{H}}_{\text{3}}^{\text{+}}$ (middle column), and water (right column) and the corresponding spin-state ratios (bottom row), as functions of radius in the H-MM1 core model. The curves correspond to t = 105 yr. Solid lines represent the FS model, while dashed lines represent the PH model.
In the text

Fig. 7 Abundances of the ammonia isotopologs summed over the spin states (left), corresponding D/H ratios (middle), and spin-state ratios (right), as functions of time in the H-MM1 core model. Solid lines represent the FS model, while dashed lines represent the PH model. Observational ratios derived by H17 are overlaid (hatched boxes). The orange hatched box in the middle panel represents the NH2 D∕NH3 ratio calculated using an NH3 o/p ratio of 0.5 (see text for details).

In the text

	Fig. 8 Left: D/H and spin-state abundance ratios in the PH model, calculated using nine different physical core models as explainedin the text. Right: corresponding data from the FS model. Hatched boxes correspond to observed ranges as in Fig. 7.
In the text

	Fig. 9 Simulated emission lines for deuterated ammonia (specific species and line indicated in the panels) in the FS model (red) and PH model (blue). The black histograms show the observed lines (H17). The results correspond to the best-fit time t = 3.29 × 105 yr (see text).
In the text

	Fig. 10 As the left and middle panels of Fig. 7, but showing the abundances and D/H ratios of some additional species.
In the text

	Fig. 11 As Fig. 7, but showing the averaged abundances in H-MM1 in the FS model (solid lines) and in the PHalt model discussed in the text (dotted lines). The hatched boxes showing the observed ranges in the right panel are omitted here for better visibility of the model results.
In the text

	Fig. 12 As Fig. 9, but for the PHalt model assuming either ζp(H) = 1.3 × 10−17 s−1 (blue) or ζp(H) = 5.0 × 10−17 s−1 (red). The corresponding best-fit chemical times are t = 4.33 × 105 yr and t = 2.50 × 105 yr.
In the text

	Fig. 13 Abundance ratios (NH2D∕NH3)∕(NHD2∕NH2D) and (NHD2∕NH2D)∕(ND3∕NHD2) as functions of time in the H-MM1 core model. The blue, red, and green curves represent the FS model, the PH model, and the PHalt model, respectively.
In the text