Modeling deuterium chemistry in starless cores: full scrambling versus proton hop

We constructed two new models for deuterium and spin-state chemistry for the purpose of modeling the low-temperature environment prevailing in starless and pre-stellar cores. The fundamental difference between the two models is in the treatment of ion-molecule proton-donation reactions of the form XH + Y −→ X + YH, which are allowed to proceed either via full scrambling or via direct proton hop, that is, disregarding proton exchange. The choice of the reaction mechanism affects both deuterium and spin-state chemistry, and in this work our main interest is on the effect on deuterated ammonia. We applied the new models to the starless core H-MM1, where several deuterated forms of ammonia have been observed. Our investigation slightly favors the proton hop mechanism over full scrambling because the ammonia D/H ratios are better fit by the former model, although neither model can reproduce the observed NH2D ortho-to-para ratio of 3 (the models predict a value of ∼2). Extending the proton hop scenario to hydrogen atom abstraction reactions yields a good agreement for the spin-state abundance ratios, but greatly overestimates the deuterium fractions of ammonia. However, one can find a reasonably good agreement with the observations with this model by increasing the cosmic-ray ionization rate over the commonly adopted value of ∼10−17 s−1. We also find that the deuterium fractions of several other species, such as H2CO, H2O, and CH3, are sensitive to the adopted proton-donation reaction mechanism. Whether the full scrambling or proton hop mechanism dominates may be dependent on the reacting system, and new laboratory and theoretical studies for various reacting systems are needed to constrain chemical models.


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 followed by the analogous reactions The dissociative recombination of NH + 4 with free electrons then forms ammonia. The most important destruction pathway (H + 3 + NH 3 −→ NH + 4 + H 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 example for H 2 , the combination of two protons gives rise to two spin states where the nuclear spin wavefunction is either antisymmetric (para-H 2 ; rotational quantum number J = 0, 2, 4...) or symmetric (ortho-H 2 ; 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 H + 3 + HD ←→ H 2 D + + H 2 + 232 K to proceed also in the backward reaction near 10 K. It is known that H 2 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 gasphase 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 2 D and NHD 2 , 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 selfconsistent description of spin-state chemistry for multiplydeuterated 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 H + 3 + H 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 2 Cl + formation in diffuse 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.

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 NH 3 + D 2 H + −→ (NH 4 D + 2 ) * −→ NH 3 D + + HD, the intermediate complex contains four H atoms and two D atoms. The four H/D atoms required to form NH 3 D + can be picked in 4 1 = 4 equivalent combinations: HHHD, HHDH, HDHH, and DHHH. The individual probabilities for each combination are the same, for example for HHHD: 4 6 3 5 2 4 2 3 = 2 15 . So, in total, the branching ratio for NH 3 D + formation is 4 1 2 15 = 8 15 . There is only one combination that forms NH + 4 (probability 1 15 ), and so the probability to form NH 2 D + 2 is 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; atomexchanging reaction channels are closed off. Therefore, the probability of forming NH 3 D + in NH 3 + D 2 H + is simply 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 3 + D 2 H + reaction first forms NH + 4 (or one of its isotopologs), which subsequently dissociates mainly into NH 3 + 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 NH 2 D + 2 , which is not possible in the PH mechanism, and NH 2 D + 2 can then decay into NHD 2 , 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 note that 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 H + 3 + H 2 system by Hugo et al. (2009), whose data we use here as well.

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 the deuteron-donation reaction, We note that 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 = 1 2 ) and deuterons (I = 1) must be treated separately. The treatment of H 2 spin is trivial: the spin state is unaltered in this reaction. For deuterium we consider the subreaction 1 D + 2 + D 2 −→ D + 4 * −→ D + D + 3 . 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 (D + 2 ) breaks up: D + 2 −→ D + + D, and in the second step, the deuterium nucleus reacts with the neutral reactant: 3 . Therefore, reaction (6) can be formally broken down to NH 2 DD + + D 2 −→ NH 2 D + D + + D 2 −→ NH 2 D + D 2 D + , (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 2 spin state. We note that the reaction channel that forms pD + 3 is forbidden by selection rules when the neutral reactant is oD 2 . There are again fewer reaction pathways in the PH model than in the FS model. For the reaction involving oD 2 , the two schemes yield the same product channels with identical branching ratios, but when pD 2 is involved, the formation of mD + 3 is forbidden in the PH model. Ammonia formation is governed mainly by hydrogen atom abstraction reactions of the type NH + i + H 2 −→ NH + i+1 + H, such as the N + + H 2 −→ 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.

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 a single reactive bulk, that is, no separation between a chemically active layer and an inert mantle is made. A63, page 3 of 14 A&A 631, A63 (2019) 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 2 , H + 2 , H + 3 , 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  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.

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 zerodimensional (0D) physical models for this purpose, that is, we calculate the chemical evolution as a function of time in pointlike 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 −17 s −1 .
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 3 (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 twophase model to ammonia abundances that are much higher than observed ). The effect of modifying the initial H 2 o/p ratio is discussed in Sect. 3.2.  Table 1. Initial abundances (with respect to the total proton number density n H ≈ 2 × n(H 2 )) used in the chemical modeling.
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.

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 H + 3 and its isotopologs are practically unaffected by the branching mechanism. This is because the chemistry of these species is mainly determined by the H + 3 + H 2 reacting system, for which we adopt the rate coefficients calculated by Hugo et al. (2009)  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 2 D/NH 3 ratio, but decreases the NHD 2 /NH 2 D and ND 3 /NHD 2 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 2 D/NH 3 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 5 yr at n(H 2 ) = 10 5 cm −3 , for example. The difference between the FS and PH models increases at late times for the NH 2 D/NH 3 ratio, but decreases for NHD 2 /NH 2 D and ND 3 /NHD 2 . 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 3 spin ratios are affected the most when changing the model. This is expected, given that ND 3 has three distinct spin states, and the PH model includes significantly less reaction channels for ND 3 formation. Statistical ammonia D/H ratios obey (NH 2 D/NH 3 )/(NHD 2 /NH 2 D) = 3 and (NHD 2 /NH 2 D)/(ND 3 /NHD 2 ) = 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 2 O and D 2 O/HDO ratios between the two models is almost density-independent. , and water (right column), including their respective deuterated forms, at a medium density of n(H 2 ) = 10 4 cm −3 (top row), n(H 2 ) = 10 5 cm −3 (middle row), and n(H 2 ) = 10 6 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.
(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 2 O + and D 3 O + (in particular the most abundant meta (m) and ortho forms). For example, mD 3 O + + e − can only lead to oD 2 O formation while oD 3 O + + e − can form ortho and para D 2 O with equal probability. At high density in the FS model, the D 3 O + m/o ratio is above unity while in the PH model it is below unity, meaning that production of para-D 2 O is more efficient in the PH model, decreasing the D 2 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 2 D, pNHD 2 , and oND 3 , 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. 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 core model. The results correspond to t = 10 5 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 2 /NH 2 D and ND 3 /NHD 2 ratios are clearly lower in the core center in the PH model. The sharp decrease in the H 2 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 5 -10 6 yr; the current model predicts a factor of a few more ND 3 .

The H-MM1 core model
The D/H ratios again follow the same trends as in the discussion above -for example the NH 2 D/NH 3 ratio is enhanced in the PH model. Notably, the FS model can only fit the NH 2 D/NH 3 ratio for realistic timescales ( 10 5 yr) where the abundances are not below detectable levels. The PH model reproduces the observed NHD 2 /NH 2 D ratio within the observational uncertainties, and simultaneously fits the observed ND 3 /NHD 2 ratio A63, page 7 of 14 A&A 631, A63 (2019) within a factor of two. While the NH 2 D/NH 3 ratio is better reproduced by the FS model, we note that H17 only detected para-NH 3 and derived the NH 2 D/NH 3 ratio using an NH 3 o/p ratio of unity. If we instead assume NH 3 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 2 D/NH 3 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 2 o/p ratio is reproduced by both models, and both underestimate the NH 2 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 H + 3 isotopologs (see the discussion in Sipilä et al. 2015), and in this paper we use the same chemistry for the H + 3 + H 2 system (Hugo et al. 2009) in both the FS and PH cases. Hence a modification to the H + 3 + H 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 2 o/p ratio of 10 −3 , 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 2 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 2 o/p ratio of 0.1. This modification leads to very low abundances of deuterated species at early times (before ∼10 5 yr) because the presence of ortho-H 2 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 5 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.

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 2 column density as well as gas and dust temperature profiles (see H17 for details). The observational error in the H 2 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 2 D/NH 3 ratio, and lower NHD 2 /NH 2 D and ND 3 /NHD 2 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 2 D/NH 3 ratio (assuming NH 3 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 3 /NHD 2 ratio is also best fit by these parameter combinations; in those cases the NHD 2 /NH 2 D ratio is just at the observed lower limit at t = 10 6 yr. Regardless of the changes to the density and temperature, we never obtain an NH 2 D o/p ratio higher than approximately two in our models, marking a clear discrepancy with the observations.

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 2 D + and N 2 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 5 and 10 7 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 χ 2 -test on the simulated lines 2 where the summation is over spectral channels. The final χ 2 value at each time step is a simple mathematical average of seven species; we exclude para-ND 3 from the test because it was not detected by H17. The best-fit time that we obtain is t = 3.29 × 10 5 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 6 yr, but at these late times 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 2 D lines are overestimated by the model, which is unsuprising given that the model predicts a lower NH 2 D o/p ratio than what is deduced from the observations; (2) the NHD 2 and ND 3 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 observations could 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 5 yr the NH 2 D and NHD 2 lines are too weak. The results highlight the fact that line emission simulations are essential for comparing observations and chemical models.

Effect of species-to-species rate coefficients
In Sipilä et al. (2017) we presented a new set of rate coefficients for the H + 3 + H 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 6 cm −3 and temperatures above 10 K, but we also found that the ortho/para ratio of H 2 D + is affected already at a density of a few times 10 5 cm −3 because of the critical density of the ground-state 1 10 -1 11 rotational transition of oH 2 D + . Given the dependency of the ammonia o/p (and D/H) ratios on the spin states of the H + 3 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 2 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 H + 3 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 2 D + o/p ratio which indicates that including rotational excitation is important for the correct interpretation of H 2 D + line emission in starless and pre-stellar cores, as already clearly demonstrated in Harju et al. (2017b).

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  D 2 CO, CH 2 D, and CHD 2 . Of this list of species, N 2 D + and DCO + are not affected by the choice of the model. This is not surprising given that both are straightforward derivatives of H 2 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 in the PH model. Formaldehyde is produced (at high density) in the dissociative electron recombination of CH 2 OH + with free electrons, and in neutral-neutral reactions of atomic oxygen with CH 3 . The fraction of deuterated CH 3 is higher in the FS model than in the PH model -this derives from differences in the chemistry of deuterated CH + 4 -leading to more efficient production of HDCO and D 2 CO through the neutral-neutral pathway. We note that (deuterated) H 2 CO and CH 3 have spin states that could in principle be tracked using our spin-separation routine, but we leave such an analysis for future work.

Hydrogen atom abstraction reactions
The analysis presented above shows that neither the FS or the PH model can reproduce the NH 2 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 2 than with para-H 2 (Dislaire et al. 2012), while the substitution of deuterium in H 2 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 2 D + + H 2 can then only lead to NH 3 D + + H -the product pair NH + 4 + 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 2 D/NH 3 and NHD 2 /NH 2 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 2 D/NH 3 abundance ratio is now over a factor of two higher than the observed value if we account for the NH 3 o/p ratio in the PH alt model (∼1.2). The value of the NH 2 D o/p ratio (∼2.7) at t ∼ 10 6 yr is close to the observed value of 3.0 ± 0.2. Likewise, the o/p ratio of NHD 2 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 cosmicray ionization rate leads to less efficient production of the deuterated H + 3 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 −17 s −1 and 1.0 × 10 −16 s −1 . 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 5 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 −17 s −1 . Figure 12 shows as an example the simulated lines for ζ p (H) = 5.0 × 10 −17 s −1 ; the core model is otherwise identical to that used in Fig. 9. The bestfit time in this case is t = 2.50 × 10 5 yr. The overall fit is better than the one obtained with the PH model, and in particular the NH 2 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 2 D o/p ratio corresponds to the observed value.
As noted in Sect. 3.1, statistical ammonia deuteration ratios obey (NH 2 D/NH 3 )/(NHD 2 /NH 2 D) = 3 and (NHD 2 /NH 2 D)/ (ND 3 /NHD 2 ) = 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 2 /NH 2 D)/(ND 3 /NHD 2 ) ratio is higher than the (NH 2 D/ NH 3 )/(NHD 2 /NH 2 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 2 D/ NH 3 )/(NHD 2 /NH 2 D) ratio in particular. The observed values in H-MM1 are (NH 2 D/NH 3 )/(NHD 2 /NH 2 D) = 1.77 ± 0.19 and (NHD 2 /NH 2 D)/(ND 3 /NHD 2 ) = 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 NH + 3 + H 2 −→ NH + 4 + 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 2 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 2 D and NHD 2 spin ratios while simultaneously failing to match the D/H ratios (assuming ζ p (H) = 1.3 × 10 −17 s −1 ) is especially puzzling. The NH 3 o/p ratio has not been measured toward H-MM1 because the low-lying rotational lines (belonging to ortho-NH 3 ) 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 2 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 freezeout, 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.

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 2 D/NH 3 , NHD 2 /NH 2 D, and ND 3 /NHD 2 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 2 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 −17 s −1 , 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 2 D + and DCO + , are unaffected by the nature of the proton-donation reaction mechanism. However, species such as H 2 CO, H 2 O, and CH 3 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.