Open Access
Issue
A&A
Volume 631, November 2019
Article Number A63
Number of page(s) 14
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201936416
Published online 22 October 2019

© O. Sipilä et al. 2019

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Open Access funding provided by Max Planck Society.

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 (1)

followed by the analogous reactions

The dissociative recombination of with free electrons then forms ammonia. The most important destruction pathway () 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 H2, the combination of two protons gives rise to two spin states where the nuclear spin wavefunction is either antisymmetric (para-H2; rotational quantum number J = 0, 2, 4...) or symmetric (ortho-H2; 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 (5)

to proceed also in the backward reaction near 10 K. It is known that H2 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 NH2 D and NHD2, 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 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 H2 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 , the intermediate complex contains four H atoms and two D atoms. The four H/D atoms required to form NH3 D+ can be picked in equivalent combinations: HHHD, HHDH, HDHH, and DHHH. The individual probabilities for each combination are the same, for example for HHHD: . So, in total, the branching ratio for NH3D+ formation is . There is only one combination that forms (probability ), and so the probability to form is .

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 NH3 D+ in NH3 + D2H+ is simply , 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 NH3 + D2H+ reaction first forms (or one of its isotopologs), which subsequently dissociates mainly into NH3 + 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 , which is not possible in the PH mechanism, and can then decay into NHD2, 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 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, (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 () and deuterons (I = 1) must be treated separately. The treatment of H2 spin is trivial: the spin state is unaltered in this reaction. For deuterium we consider the subreaction1 . 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 () breaks up:

and in the second step, the deuterium nucleus reacts with the neutral reactant:

Therefore, reaction (6) can be formally broken down to (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 H2 spin state. We note that the reaction channel that forms p is forbidden by selection rules when the neutral reactant is oD2. There are again fewer reaction pathways in the PH model than in the FS model. For the reaction involving oD2, the two schemes yield the same product channels with identical branching ratios, but when pD2 is involved, the formation of m is forbidden in the PH model.

Ammonia formation is governed mainly by hydrogen atom abstraction reactions of the type , such as the N+ + H2→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.

thumbnail 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.

Open with DEXTER
thumbnail Fig. 2

Spin-state branching ratio diagrams for the reaction in the FS (left) and PH (right) models, omitting the H2 spin state.

Open with DEXTER

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 H2, , , 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 Tdust = Tgas = 10 K, a visual extinction of AV = 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 NH3 (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 H2 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.

thumbnail Fig. 3

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

Open with DEXTER
Table 1

Initial abundances (with respect to the total proton number density nH ≈ 2 × n(H2)) used in the chemical modeling.

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 and its isotopologs are practically unaffected by the branching mechanism. This is because the chemistry of these species is mainly determined by the 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) 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 NH2D∕NH3 ratio, but decreases the NHD2∕NH2D and ND3∕NHD2 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 NH2D∕NH3 ratio can exceed unity in the current model. There is an evident peak for the ratios which depends on time; this occurs at ~ 2 × 105 yr at n(H2) = 105 cm−3, for example. The difference between the FS and PH models increases at late times for the NH2D∕NH3 ratio, but decreases for NHD2∕NH2D and ND3∕NHD2. 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 ND3 spin ratios are affected the most when changing the model. This is expected, given that ND3 has three distinct spin states, and the PH model includes significantly less reaction channels for ND3 formation. Statistical ammonia D/H ratios obey (NH2D∕NH3)∕(NHD2∕NH2D) = 3 and (NHD2∕NH2D)∕(ND3∕NHD2) = 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/H2O and D2 O/HDO ratios between the two models is almost density-independent. The H2O o/p ratio is largely unaffected by the choice of branching mechanism. The D2 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 HD2O+ and D3O+ (in particular the most abundant meta (m) and ortho forms). For example, mD3O+ + e can only lead to oD2 O formation while oD3 O+ + e can form ortho and para D2 O with equal probability. At high density in the FS model, the D3O+ m/o ratio is above unity while in the PH model it is below unity, meaning that production of para-D2O is more efficient in the PH model, decreasing the D2O 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, pNH2D, pNHD2, and oND3, 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.

thumbnail Fig. 4

D/H abundance ratios of ammonia (left column), (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.

Open with DEXTER
thumbnail Fig. 5

Spin-state abundance ratios of ammonia (left column), (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.

Open with DEXTER

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 = 105 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 NHD2∕NH2D and ND3∕NHD2 ratios are clearly lower in the core center in the PH model. The sharp decrease in the H2 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 105 –106 yr; the current model predicts a factor of a few more ND3.

The D/H ratios again follow the same trends as in the discussion above – for example the NH2D∕NH3 ratio is enhanced in the PH model. Notably, the FS model can only fit the NH2D∕NH3 ratio for realistic timescales (≳105 yr) where the abundances are not below detectable levels. The PH model reproduces the observed NHD2∕NH2D ratio withinthe observational uncertainties, and simultaneously fits the observed ND3∕NHD2 ratio within a factor of two. While the NH2D∕NH3 ratio is better reproduced by the FS model, we note that H17 only detected para-NH3 and derived the NH2D∕NH3 ratio using an NH3 o/p ratio of unity. If we instead assume NH3 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 NH2D∕NH3 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 NHD2 o/p ratio is reproduced by both models, and both underestimate the NH2D 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 isotopologs (see the discussion in Sipilä et al. 2015), and in this paper we use the same chemistry for the system (Hugo et al. 2009) in both the FS and PH cases. Hence a modification to the 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 H2 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 H2 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 H2 o/p ratio of 0.1. This modification leads to very low abundances of deuterated species at early times (before ~ 105 yr) because the presence of ortho-H2 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 ta few × 105 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.

thumbnail Fig. 6

D/H ratios (top row) of ammonia (left column), (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.

Open with DEXTER
thumbnail 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).

Open with DEXTER

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 H2 column density as well as gas and dust temperature profiles (see H17 for details). The observational error in the H2 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 NH2D∕NH3 ratio, and lower NHD2∕NH2D and ND3∕NHD2 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 NH2D∕NH3 ratio (assuming NH3 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 ND3∕NHD2 ratio is also best fit by these parameter combinations; in those cases the NHD2∕NH2D ratio is just at the observed lower limit at t = 106 yr. Regardless of the changes to the density and temperature, we never obtain an NH2D 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 H2 D+ and N2 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 105 and 107 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 lines2 (8)

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-ND3 from the test because it was not detected by H17. The best-fit time that we obtain is t = 3.29 × 105 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 × 106 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-NH2D lines are overestimated by the model, which is unsuprising given that the model predicts a lower NH2D o/p ratio than what is deduced from the observations; (2) the NHD2 and ND3 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 × 105 yr the NH2D and NHD2 lines are too weak. The results highlight the fact that line emission simulations are essential for comparing observations and chemical models.

thumbnail 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.

Open with DEXTER

4.3 Effect of species-to-species rate coefficients

In Sipilä et al. (2017) we presented a new set of rate coefficients for the 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 106 cm−3 and temperatures above 10 K, but we also found that the ortho/para ratio of H2 D+ is affected already at a density of a few times 105 cm−3 because of the critical density of the ground-state 110 –111 rotational transition of oH2D+. Given the dependency of the ammonia o/p (and D/H) ratios on the spin states of the 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 H2 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 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 H2 D+ o/p ratio which indicates that including rotational excitation is important for the correct interpretation of H2 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 N2 D+, DCO+, HDCO, D2 CO, CH2D, and CHD2. Of this list of species, N2D+ and DCO+ are not affected by the choice of the model. This is not surprising given that both are straightforward derivatives of H2 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 CH2OH+ with free electrons, and in neutral-neutral reactions of atomic oxygen with CH3. The fraction of deuterated CH3 is higher in the FS model than in the PH model – this derives from differences in the chemistry of deuterated – leading to more efficient production of HDCO and D2CO through the neutral-neutral pathway. We note that (deuterated) H2CO and CH3 have spin states that could in principle be tracked using our spin-separation routine, but we leave such an analysis for future work.

thumbnail 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).

Open with DEXTER
thumbnail Fig. 10

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

Open with DEXTER
thumbnail 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.

Open with DEXTER

4.5 Hydrogen atom abstraction reactions

The analysis presented above shows that neither the FS or the PH model can reproduce the NH2D 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-H2 than with para-H2 (Dislaire et al. 2012), while the substitution of deuterium in H2 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 NH2D+ + H2 can then only lead to NH3D+ + H – the product pair is not accessible like it is in the FS model, meaning that this alternative model (hereafter PHalt) 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 PHalt model.

Figure 11 shows the averaged abundances in H-MM1 calculated with the FS model and the PHalt model. Switching to a direct hydrogen atom hop mechanism for the abstraction reactions greatly increases the D/H ratios, as expected, and indeed the NH2D∕NH3 and NHD2∕NH2D ratios are now both higher than the observed ones. PHalt produces ortho and para ammonia in roughly equal proportions, and in fact the modeled NH2D∕NH3 abundance ratio is now over a factor of two higher than the observed value if we account for the NH3 o/p ratio in the PHalt model (~1.2). The value of the NH2D o/p ratio (~2.7) at t ~ 106 yr is close to the observed value of 3.0 ± 0.2. Likewise, the o/p ratio of NHD2 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 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 105 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 PHalt 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 best-fit time in this case is t = 2.50 × 105 yr. The overall fit is better than the one obtained with the PH model, and in particular the NH2D 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 PHalt 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 NH2D o/p ratio corresponds to the observed value.

As noted in Sect. 3.1, statistical ammonia deuteration ratios obey (NH2D∕NH3)∕(NHD2∕NH2D) = 3 and (NHD2∕NH2D)∕ (ND3∕NHD2) = 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 PHalt). We recover the same tendency as for point models (Sect. 3.1), that is, that the (NHD2∕NH2D)∕(ND3∕NHD2) ratio is higher than the (NH2D∕NH3)∕(NHD2∕NH2D) ratio in the FS model, and vice versa for the PH model. Switching from the PH to the PHalt model increases both ratios by a similar factor. These results highlight the strong model-dependent behavior of the (NH2D∕NH3)∕(NHD2∕NH2D) ratio in particular. The observed values in H-MM1 are (NH2D∕NH3)∕(NHD2∕NH2D) =1.77 ± 0.19 and (NHD2∕NH2D)∕(ND3∕NHD2) = 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 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 D2 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 PHalt model shows such good agreement with the observations for the NH2D and NHD2 spin ratios while simultaneously failing to match the D/H ratios (assuming ζp (H) = 1.3 × 10−17 s−1) is especially puzzling. The NH3 o/p ratio has not been measured toward H-MM1 because the low-lying rotational lines (belonging to ortho-NH3) cannot be observed from the ground, but the ratio that the PHalt 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 NH2 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.

thumbnail 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.

Open with DEXTER
thumbnail 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.

Open with DEXTER

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 NH2D∕NH3, NHD2∕NH2D, and ND3∕NHD2 ratios observed toward H-MM1 by Harju et al. (2017a) than the FS model does, but in particular the observed o/p ratio of NH2D 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 N2 D+ and DCO+, are unaffected by the nature of the proton-donation reaction mechanism. However, species such as H2 CO, H2 O, and CH3 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


1

In this paper we refer to any systems that represent a part of a complete reaction, such as reaction (6), as “subreactions”.

2

For the summation we only included channels for which the main beam temperature is greater than 3σ, as determined from a line-free part of the observed spectrum.

All Tables

Table 1

Initial abundances (with respect to the total proton number density nH ≈ 2 × n(H2)) used in the chemical modeling.

All Figures

thumbnail 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.

Open with DEXTER
In the text
thumbnail Fig. 2

Spin-state branching ratio diagrams for the reaction in the FS (left) and PH (right) models, omitting the H2 spin state.

Open with DEXTER
In the text
thumbnail Fig. 3

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

Open with DEXTER
In the text
thumbnail Fig. 4

D/H abundance ratios of ammonia (left column), (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.

Open with DEXTER
In the text
thumbnail Fig. 5

Spin-state abundance ratios of ammonia (left column), (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.

Open with DEXTER
In the text
thumbnail Fig. 6

D/H ratios (top row) of ammonia (left column), (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.

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail Fig. 10

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

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.