Nuclear spin ratios of deuterated ammonia in prestellar cores LAsMA observations

Context. Molecules containing two or more hydrogen or deuterium atoms have different nuclear spin states which behave as separate chemical species. The relative abundances of these species can give clues to their origin. Formation on grains is believed to yield statistical spin ratios whereas gas-phase reactions are predicted to result in clear deviations from them. This is also true for ammonia and its deuterated forms NH 2 D, NHD 2 , and ND 3 . Aims. Here we aim to determine the ortho/para ratios of NH 2 D and NHD 2 in dense, starless cores, where their formation is supposed to be dominated by gas-phase reactions. Methods. The Large APEX sub-Millimeter Array (LAsMA) multibeam receiver of the Atacama Pathfinder EXperiment (APEX) tele-scope was used to observe the prestellar cores H-MM1 and OphD in Ophiuchus in the ground-state lines of ortho and para NH 2 D and NHD 2 . The fractional abundances of these molecules were derived employing three-dimensional radiative transfer modelling, using different assumptions about the abundance profiles as functions of density. We also ran gas-grain chemistry models with different scenarios concerning proton or deuteron exchanges and chemical desorption from grains to find out if one of these models can reproduce the observed spin ratios. Results. The observationally deduced ortho/para ratios of NH 2 D and NHD 2 are in both cores within 10% of their statistical values 3 and 2, respectively, and taking 3 σ limits, deviations from these of about 20% are allowed. Of the chemistry models tested here, the model that assumes proton hop (as opposed to full scrambling) in reactions contributing to ammonia formation, and a constant efficiency of chemical desorption, comes nearest to the observed abundances and spin ratios. Conclusions. The nuclear spin ratios derived here are in contrast with spin-state chemistry models that assume full scrambling in proton donation and hydrogen abstraction reactions leading to deuterated ammonia. The efficiency of chemical desorption strongly influences the predicted abundances of NH 3 , NH 2 D, and NHD 2 , but has a lesser effect on their ortho/para ratios. For these the proton exchange scenario in the gas is decisive. We suggest that this is because of rapid re-processing of ammonia and related cations by gas-phase ion-molecule reactions.


Introduction
Astrochemical models that distinguish between different nuclear spin modifications assume that ion-molecule reactions in interstellar gas proceed via the formation of a relatively long-lived intermediate complex, where the light H and D nuclei can be completely mixed.In this scenario, one of the main deuteration reactions of ammonia in dense gas can be written as NH 3 + H 2 D + → (NH 5 D + ) † → NH 3 D + + H 2 , where the reaction complex is indicated with a dagger.This reaction is followed by electron recombination of NH 3 D + , which can yield NH 3 or NH 2 D. The increasing abundances of H + 3 , H 2 D + , HD + 2 , and D + 3 in the interiors of dense, starless cores, and the possibility of mixing between A&A proofs: manuscript no.lasma_arxiv The few observational determinations of the o/p ratios of NH 2 D or NHD 2 in dense cores are consistent with statistical ratios (Gerin et al. 2006;Lis et al. 2006;Daniel et al. 2016;Harju et al. 2017;Wienen et al. 2021).The uncertainties of these estimates are, however, large -typically greater than 30%, and one cannot confidently exclude the deviations from statistical values predicted by chemistry models.
Here we attempt to determine the o/p ratios of NH 2 D and NHD 2 in two dense and cold interstellar gas clouds with greater accuracy than achieved previously, using the Large APEX sub-Millimeter Array (LAsMA) of the Atacama Pathfinder EXperiment (APEX).The ground state transitions of the ortho and para modifications of NH 2 D and NHD 2 (hereafter oNH 2 D, pNH 2 D, oNHD 2 , and pNHD 2 ) at λ ∼ 0.9 mm were observed with this instrument towards the starless cores H-MM1 and Oph D in the region of the dark cloud L1688 in Ophiuchus.The large-scale structure and general properties of the cores in the L1688 region have been discussed in, for example, Pattle et al. (2015), Friesen et al. (2017), Choudhury et al. (2020), and Ladjelate et al. (2020).Detailed studies of the structures and the ammonia distributions in the two target cores have been presented in Harju et al. (2020) and Pineda et al. (2022) (for H-MM1), and in Ruoskanen et al. (2011) (for Oph D).The two cores can be classified as prestellar, meaning that they are on the verge of collapse, as they have central densities significantly above 10 5 cm −3 , and high levels of deuterium fractionation (e.g., Crapsi et al. 2005;Harju et al. 2017).
In Sect. 2 of this paper we describe the LAsMA observations, and in Sect. 3 we present the data and estimate the abundances and the o/p ratios of NH 2 D and NHD 2 in the two cores.Spin ratios expected from the proton hop and full scrambling scenarios, and the results of previous spin-state chemistry models are discussed in Sect. 4. In Sect.5, we show results from an updated version of our gas-grain chemistry model.Implications of the observed spin ratios and NH 2 D/NH 3 and NHD 2 /NH 2 D abundance ratios are discussed in Sect.6.Finally, in Sect.7 we draw our conclusion from this study.

Observations
The ground-state rotational lines of oNH 2 D, pNH 2 D, oNHD 2 , and pNHD 2 near 333 and 335 GHz were observed with the LAsMA array (Güsten et al. 2008) on the APEX telescope (Güsten et al. 2006) towards the prestellar cores H-MM1 and Oph D in Ophiuchus.The observations were done between April 10 and September 28, 2022, under the project number M-0109.F-9503A-2022.The APEX telescope is located at an altitude of 5100 m at Llano de Chajnantor, in Chile.
The LAsMA array has 7 pixels separated by approximately 40 ′′ .Figure 1 shows the footprint of the instrument superposed on the 8 µm maps of the cores measured by the InfraRed Array Camera (IRAC) of the Spitzer Space Telescope.The beamsize of APEX is 18 ′′ (FWHM) at 335 GHz.The pattern was unintentionally tilted by 19 • .1 during observations in the Spring and Summer 2022.This was corrected during additional observations in September towards Oph D (solid circles on the right).With two tilt angles, the total number of positions observed towards Oph D is 13.
The backend was composed of Fast Fourier Transform Spectrometers (FFTS4G) that covered two 4 GHz intermediate frequency (IF) bands resolved into 65 536 spectral channels of 61.03 kHz in width (corresponding to ∼ 55 m s −1 ).The lines of ortho and para NH 2 D and NHD 2 were observed in the lower sideband (LSB) of the receiver which was centered at  1. Also listed there are the energies of the upper transition level and the Einstein A coefficients of the transitions.The spectroscopic data are adopted from Melosso et al. (2021).At 10 K, the optically thin critical densities of the four transitions (estimated using Eq. ( 4) of Shirley 2015) lie in a narrow range between 6.9 × 10 4 cm −3 and 7.9 × 10 4 cm −3 .The NH 2 D and NHD 2 spectra towards the density peaks of the two cores are shown in Fig. 2. The positions correspond to pixel 1 for H-MM1 and pixel 7 for Oph D in Fig. 1.The spectra are shown on the main-beam brightness temperature, T MB , scale.In the conversion from T * A scale, we assumed the following main-beam and forward-scattering and spill-over efficiencies: η MB = 0.73, η fss = 0.97.In the course of the measurements, we noticed line intensity changes up to 15%, which probably are caused by occasional rapid variations in the water vapour content of the atmosphere.We estimate that the absolute calibration is accurate to 20%.Because the four lines were observed in the same spectral band, this uncertainty does not concern their relative intensities, and does not propagate to the derived o/p ratios.

LTE analysis of the spectra
The line parameters and column densities derived from Gaussian fits to the hyperfine structure of the lines observed towards the density peaks of H-MM1 and Oph D are presented in Table 2.The fits are shown with thin red lines in Fig. 2. In the column density estimates we have assumed line-of-sight homogeneity and that the excitation temperature, T ex , is constant for all rotational transitions of the molecule.Furthermore, the calculation involves the assumption that populations of hyperfine states in a certain rotational level are proportional to their statistical weights.The transitions have three groups which consist of several unresolved components (Coudert & Roueff 2006;Melosso et al. 2021).The linewidth listed in Table 2 refers to the intrinsic width of an individual hyperfine component.One can see in this Table that the NH 2 D lines are broader than those of NHD 2 , and the difference in H-MM1 clearly exceeds that expected from thermal broadening, ∼ 5 m s −1 at 10 K.The NHD 2 abundance is likely to peak at slightly higher densities and deeper in the core as compared to NH 2 D, and the linewidth difference could possibly be understood in terms of lower temperatures and smaller non-thermal velocity dispersion in the interior parts of the cores.The 3D models described below include realistic temperature gradients but assume a constant non-thermal velocity dispersion.
Despite rather high signal-to-noise ratios especially in H-MM1, the optical thicknesses and the excitation temperatures of the pNH 2 D and pNHD 2 lines have large relative uncertainties which propagate to the column densities.The total optical thicknesses (sums of three separate groups) of the pNH 2 D and pNHD 2 lines derived towards H-MM1 are below 1, suggesting that one could use for the para lines the optically thin approximation, where the column density is proportional to the integrated intensity.In this method, the T ex must be assumed.We examined the excitation of the observed transitions by applying radiative transfer calculations to the 3D model of H-MM1 used previously by Pineda et al. (2022), and described in more detail in Sect.3.2.Here, the molecules are assumed to have constant fractional abundances, corresponding to the column densities derived from the hyperfine fits described above.The T ex profiles along a line crossing the density peak (with n(H 2 ) ∼ 1.2 × 10 6 cm −3 ) are shown in Fig. 3.According to this calculation, the T ex of pNH 2 D line is consistently ∼ 10% lower than that of the oNH 2 D line, whereas for the para-and ortho-NHD 2 lines the T ex ratio is 1.The higher T ex of the oNH 2 D line compared to that of the pNH 2 D line is likely to be caused by an increased effect of radiative trapping with larger optical thickness (Shirley 2015;Roueff et al. 2021).In the last column of Table 2 we give the pNH 2 D and pNHD 2 column densities derived from the integrated intensities assuming that T ex,para = 0.9 × T ex,ortho for NH 2 D, and T ex,para = T ex,ortho for NHD 2 .
Using para column densities from optically thin approximation, the values listed in Table 2 give the following o/p-NH 2 D ratios: 3.0 ± 0.2 (H-MM1) and 3.5 ± 0.6 (Oph D).The corresponding o/p-NHD 2 ratios are 2.2 ± 0.5 (H-MM1) and 2.4 ± 1.6 (Oph D).If we use the assumption T ex,para = T ex,ortho also for NH 2 D, the o/p-NH 2 D ratios increase by ∼ 20% from the values given above.To distinguish between different formation scenarios for these molecules, the o/p ratios should be determined more accurately, which can be achieved via non-LTE modelling as we will discuss below.
Besides the uncertainty involved in the determination of the optical thickness, τ, from only three resolved components, the assumption of a constant T ex along the line of sight in the presence of large density gradients is an additional source of uncertainty, as demonstrated in Fig. 3.In what follows, we construct three-dimensional physical models for the H-MM1 and Oph D cores, and estimate the fractional abundances of the observed species using radiative transfer calculations.

3D radiative transfer modelling
We derived the fractional abundances of the observed molecules by comparing simulated spectra with the observations.Physical models for H-MM1 and Oph D were constructed using the 8 µm surface brightness maps shown in Fig. 1.Both cores appear as absorption features in the maps, and this was used to derive highresolution H 2 column density maps.The method is described in Appendix A of Harju et al. (2020) and in Pineda et al. (2022).The density structures of the cores were estimated by fitting a Plummer-type function to the cross-sectional H 2 column density profiles in different positions along the ridges of the cores.The method is adopted from Arzoumanian et al. (2011).We assumed that the density distribution has circular symmetry in the plane perpendicular to the spine of the core.The latter is defined by a spline fit through local N(H 2 ) peaks.A more detailed description of the method is presented in Appendix D of Pineda et al. (2022).
The H 2 column densities of the models of H-MM1 and Oph D are comparable to those that can be obtained from the high-resolution N(H 2 ) map of L1688 derived from multiwavelength Herschel data (Ladjelate et al. 2020), which has an effective resolution of 18 ′′ .2. When the model column density distributions are smoothed to the same resolution, the density peaks have the following H 2 column densities (values from the Herschel map of Ladjelate et al. (2020) are given in brackets): H-MM1 (0, 0) N(H 2 ) = 3.7 × 10 22 cm −2 (4.4 × 10 22 cm −2 ); Oph D (−19 ′′ , −37 ′′ ) N(H 2 ) = 3.2 × 10 22 cm −2 (2.2 × 10 22 cm −2 ); Oph D (+20 ′′ , +35 ′′ ) N(H 2 ) = 2.2 × 10 22 cm −2 (2.2 × 10 22 cm −2 ).The fact that the total N(H 2 ) column density derived from Herschel is lower than that of the core model at Oph D (−19 ′′ , −37 ′′ ) can be understood by the relatively high dust colour temperature (T C = 12.9 K) derived from the far-infrared data.The mass-averaged physical temperature is probably lower towards this position (see below).The high-resolution N(H 2 ) and T C Fig. 2. NH 2 D and NHD 2 spectra from the peak positions of H-MM1 and Oph D (positions 1 and 7 in Fig. 1, respectively).Gaussian fits to the hyperfine structure are shown with thin red lines.The gas temperature was assumed to be equal to the dust temperature, T dust , in the inner parts of the cores.Based on the T kin distribution in H-MM1 derived from NH 3 by Pineda et al. (2022) (see their Figs.5 and 11) we assumed, however, that the gas temperature does not rise above 11 K in the outer parts of the cores owing to molecular line cooling.The dust temperature was computed using the continuum radiative transfer program CRT (Juvela 2005).In the case of H-MM1, the external radi-ation field was assumed to consist of an isotropic component and a point source located on the western side of the core.The point source represents the subgiant binary HD 147889 located about 1.2 pc west of H-MM1.The density model of H-MM1 and the assumptions about the external radiation field are the same as those used in Pineda et al. (2022).For Oph D, which is known to be cooler than H-MM1 (∼ 6 − 7 K in the densest parts; Harju et al. 2008;Ruoskanen et al. 2011), the external radiation field was assumed isotropic and weaker by a factor of two than that around H-MM1.The spectrum of the unattenuated interstel- lar radiation field (ISRF) was taken from Black (1994).We used the dust opacity data from Ossenkopf & Henning (1994) for unprocessed dust grains with thin ice coatings 1 .The density and dust temperature models of H-MM1 and Oph D are shown in Fig. 4. The assumed levelling off of the gas temperature at 11 K only affects the model of H-MM1, where T dust rises to ∼ 14 K in the outer parts.
The fractional abundances were determined by applying minimum chi-square estimation.We ran three sets of radiative transfer calculations with different assumptions about the abundance distributions: (a) one assuming that the fractional NH 2 D and NHD 2 abundances are constant throughout the core, (b) another assuming that the abundances are constant in the outer parts but start to decrease when the H 2 density exceeds 2 × 10 5 cm −3 , and (c) a third one assuming that the abundances are constant in the density range 5×10 4 cm −3 < n(H 2 ) < 2×10 5 cm −3 but decrease both towards the outer edge and the core centre.In models (b) and (c), the fractional abundances in the central parts were assumed to be inversely proportional to the density, that is, X ∝ n −1 , owing to freezing onto the grains.This tendency was recently found in H-MM1 for NH 3 (Pineda et al. 2022), and it is possible that also the deuterated forms of ammonia behave that way (see, e.g., Caselli et al. 2022).In model (c), the abundance in the outer part was assumed to be, somewhat arbitrarily, directly proportional to the density, X ∝ n.At the core edge, the abundances of NH 2 D and NHD 2 are expected to increase towards higher densities as the production of these molecules is favoured by CO depletion.This behaviour is also predicted by our chemistry models (see Fig. 11 below).
1 hera.ph1.uni-koeln.de/ossk/Jena/tables.html For each test value of the fractional abundance and in each position included, we evaluated the statistic where T obs,i is the observed intensity and T mod,i is the computed intensity in channel i, and σ is the RMS noise of the observed spectrum.The summation is over approximately 100 channels including the three hyperfine groups.The sum of the obtained χ 2 values over two or three positions shows a parabolic distribution as a function of the fractional abundance.The x-coordinate of the minimum corresponds to the best-fit fractional abundance.
We fitted a parabola to this distribution, and denote the fitted function with χ2 .The best-fit values, X, and their 1σ errors, σ X , were calculated from where P χ ∝ e χ2 /2 is the normalised probability distribution.The method is based on the assumption that the error distribution of the spectra is Gaussian (Andrae 2010).We tested the error estimates of the X values using Monte Carlo.To this end, we generated a few thousand realisations of the spectra by adding normally distributed noise (σ = RMS) to each channel, and calculated the 'best-fit' X values by minimum chi-squared estimation.The distribution of these X values is Gaussian and the range containing 68 % of the values agrees with X ± σ X derived above.The resulting maximum and average fractional abundances and their 1 σ uncertainties are presented in Table 3.By 'maximum' abundance we mean the undepleted, baseline abundance, and by 'average' abundance we mean the column density ratio N(mol)/N(H 2 ) in the direction of the density peak of the core model.The average abundances are used for comparison with the predictions of the chemistry model in Sect. 5.
The assumption of molecular depletion at high (or low) densities results naturally in higher baseline (undepleted) abundances and lower average abundances towards the density peak compared to the assumption of a constant X.However, the relative abundances of the four species are approximately the same for all three X profiles.As the main interest of this study is in the o/p ratios, we do not attempt to optimise the depletion power laws or the density thresholds.The simulated spectra using the X=constant model are shown in Appendix A together with the observed spectra.
One can see from Table 3 that the best-fit o/p-NH 2 D and o/p-NHD 2 ratios derived towards both cores are close to their statistical values 3 and 2, respectively, and that these ratios do not depend much on the abundance profile assumed.For H-MM1, the 3 σ limits are within 12 % from the statistical ratios.For Oph D, both o/p-NH 2 D and o/p-NHD 2 are higher than those derived in H-MM1, and the 3 σ limits allow deviations up to 15 % and 24 % in excess of the statistical values 3 and 2, respectively.
noise level of the present 0.8 mm spectra is a factor of two lower than those used in the earlier study, but perhaps a more important difference is that the analysis is now based on several positions and a realistic 3D model is used for the core.

Spin ratios from the proton hop and full scrambling scenarios and the predictions of previous chemistry models
Deuterium enrichment of NH 3 in dense interstellar gas is supposed to occur primarily through deuteron donation from H 2 D + , D 2 H + , or D + 3 to ammonia.Deuterium can also enter the ammonia formation sequence in reactions of other nitrogen hydrides, such as NH, with the cations H 2 D + , D 2 H + , and D + 3 .These exother-mic reactions can be thought to proceed via a direct proton or deuteron hop (hereafter PH) or via the formation of an intermediate complex where H and D nuclei can be exchanged.In the latter mechanism, called full scrambling (hereafter FS), the nuclear spin state of the resulting cation, for example NH 3 D + , depends on the initial spin states of both reactants according to the selection rules that take into account the conservation of the spin and symmetry (e.g., Rist et al. 2013;Faure et al. 2013;Sipilä et al. 2015;Hily-Blant et al. 2018).The spin states of the H 2 D + , D 2 H + , and D + 3 ions, in turn, are determined by spin-state conversions in the H + 3 + H 2 reaction system and its deuterated variants.The state-to-state and ground-state-to-species rate coefficients of this system were derived by Hugo et al. (2009), and these have been used in spin-state chemistry models ever since The letters indicate three different assumptions about the abundance distributions: (a) X is constant throughout the core; (b) X decreases in the densest parts with n(H 2 ) > 2 × 10 5 cm −3 , following the power law n −1 ; (c) X increases proportional to n(H 2 ) at densities below 5 × 10 4 cm −3 , and decreases again as in model (b) above 2 × 10 5 cm −3 .For models including depletion, 'maximum' means the undepleted fractional abundance and 'average' means the column density ratio N(mol)/N(H 2 ) through the density peak of the core.‡The X values are given in units of 10 −9 .
To illustrate how the reaction details affect the nuclear spin ratios, we discuss the outcome of a basic sequence leading to NH 2 D, assuming that this follows either the PH or the FS scenario.In the PH mechanism, H 2 D + simply donates the deuteron to ammonia, and the spin symmetry of the hydrogen complex does not change.The H 3 complex of E ("para") symmetry dissociates to form oH 2 + H and pH 2 + H at equal probabilities, whereas the A 1 ("ortho") symmetry dissociates to form oH 2 + H (if not H + H + H).Assuming equal amounts of oNH 3 and pNH 3 , this sequence produces three times more oNH 2 D than pNH 2 D. In the FS scenario, the first step forms the reaction complex (NH 5 D + ) † .With reactants oNH 3 and pH 2 D + , the resulting H 5 complex is of G 1 symmetry, which gives H 3 (A 1 ) ("ortho") and H 3 (E) ("para") at equal probabilities when dissociating (Hugo et al. 2009, Table III;Hily-Blant et al. 2018, Table C2).The H 5 complex formed in the reaction pNH 3 +pH 2 D + has H 1 symmetry with strong preference to H 3 (E) in dissociation.Assuming again equal amounts of oNH 3 and pNH 3 , the FS sequence results in an o/p-NH 2 D ratio of 27/13 ≈ 2.1.
The difference between the PH and FS scenarios is similar for other routes to NH 2 D. For example, according to the PH model, the sequence NHD + H 2 → NH 2 D + H 2 → NH 3 D + e − → NH 2 D produces again oNH 2 D and pNH 2 D in the statistical ratio 3:1, whereas the same sequence in the FS model results in an o/p-NH 2 D ratio of 7/5 = 1.4,assuming that the reacting H 2 is always in the more common para form.Finally, the formation of NHD 2 involves the combination of two deuterium nuclei.If this happens one by one, as in the PH model, the resulting o/p-NHD 2 ratio is always 2. In the FS model, the outcome depends on the deuteron donor.For example, in a reaction between a singly deuterated species and oD 2 H + , the o/p ratio of the D 2 complex formed is 13/5=2.6.The branching ratios of the reactions discussed here can be verified by inspecting the nuclear spin symmetry induction and subduction matrices for systems with multiple H and D atoms presented in Hugo et al. (2009), Sipilä et al. (2015), andHily-Blant et al. (2018).
The fact that the FS assumption in proton-donation reactions leads to deviations of the o/p ratios of NH 2 D and NHD 2 from their high-temperature statistical values was noted by Sipilä et al. (2015), Daniel et al. (2016), and Hily-Blant et al. (2018) (see Table 5 of Hily-Blant et al. 2018).This disagreed with observations which suggested statistical ratios, albeit with A&A proofs: manuscript no.lasma_arxiv large uncertainties.Daniel et al. (2016) pointed out that the observed o/p ratios of NH 2 D and NHD 2 are consistent with formation on grains, while Harju et al. (2017) argued that statistical spin ratios could be explained by gas-phase reactions if they can be characterised as single-particle hops, in analogy with H and D atom additions on grains.The latter suggestion was examined by Sipilä et al. (2019a) who considered the effects of PH vs. FS in proton-donation reactions, such as NH 3 + D 2 H + , and in the hydrogen abstraction chain leading to the ammo- , and its deuterated variants.It turned out that replacing the FS with PH in proton-donation reactions brought the predicted o/p ratios a little closer to the statistical values but the difference was still greater than 30%.On the other hand, by extending the PH scenario also to H abstraction reactions, the o/p ratios agreed with those observed previously towards H-MM1, but the modelled fractionation ratios NH 2 D/NH 3 and NHD 2 /NH 2 D exceeded those derived from the observations.

Gas-grain chemistry model with dynamic chemical desorption
One way to explain the discrepancy between the observed and modelled spin ratios is that the desorption of ammonia from grain surfaces is in reality more efficient than what has been considered in the previous simulations.Ammonia does not desorb thermally in low-temperature conditions owing to its high binding energy on ice; usually values around 4000 − 5000 K are assumed for the NH 3 binding energy on water ice (Minissale et al. 2022;Ferrero et al. 2023).Although recent studies indicate that we should rather speak of a broad distribution of binding energies on amorphous water ice (Grassi et al. 2020;Tinacci et al. 2022), the desorption of ammonia most likely occurs through non-thermal mechanisms.
There are several non-thermal desorption mechanisms that are typically considered in astrochemical models, namely photodesorption, cosmic-ray induced desorption, and chemical desorption (also known as reactive desorption).The first two mechanisms do not release ammonia into the gas phase efficiently enough to cause significant "contamination" of the gas-phase ammonia reservoir.Chemical desorption (CD) on the other hand has been shown to have a very large impact on ammonia abundance in starless core conditions (Sipilä et al. 2019b), assuming that the desorption efficiency is constant at 1 % (Garrod et al. 2007).Hence, out of the three options mentioned here, CD is the most promising candidate for pushing the gas-phase ammonia spin ratios toward statistical values.However, the constant efficiency assumption for CD was found to lead to gas-phase ammonia abundances much higher than what is indicated by observations for volume densities around 10 4 cm −3 (Sipilä et al. 2019b).
It is reasonable to posit that in reality the efficiency of CD is in fact not universally constant, but would instead depend not only on the desorbing molecule but also on the type of surface on which the desorption is taking place.An alternative description of CD exploring this scenario has been presented by Minissale et al. (2016), and applied to a chemical model by Vasyunin et al. (2017).An important conclusion in these works is that the CD efficiency tends to increase when switching from a water ice surface to a CO ice surface, which in a starless core would correspond to the transition from outer regions to the dense inner core where the freeze-out timescale is short and the grains are covered by a multitude of species, and not only water.Motivated by these findings, we investigate here whether the dynamic nature of the CD process influences the spin ratios of deuterated ammonia -we recall that the observed ammonia lines originate in volume densities approaching 10 5 cm −3 where the ices on the grains are made up of a mixture of molecules.
We use for the simulations the chemical model described in Sipilä et al. (2019a), and adopt the full scrambling and proton hop chemical networks presented therein.In the proton hop model, proton donation and abstraction reactions proceed via the exchange of one proton (or deuteron) between the reactants, while in the full scrambling model multiple atom exchanges are possible.Also, for the sake of simplicity, we employ a one-dimensional physical model for H-MM1 recently presented in Pineda et al. (2022;their Fig. 11).The adoption of a one-dimensional spherically symmetric model means that we can with reasonable accuracy compare the simulation results to the observations obtained toward the dust peak, but we cannot attempt to reproduce the morphology of the deuterated ammonia distributions beyond a distance of a few tens of arcseconds from the core center.
We consider two different approaches to simulating CD: one where the CD efficiency is set to a constant value of 1 % (Garrod et al. 2007), and another where the CD efficiency changes dynamically.The latter is an updated version of the model of Vasyunin et al. (2017) -in the present work, the energy released in a chemical desorption event is partitioned among the reaction products according to their total degrees of freedom.In the dynamic model, the CD efficiency of a given species depends on the ice composition, and hence on the location in the core and on the simulation time, as well as the properties (mass, binding energy) of the species itself.At late simulation times, however, we obtain nearly constant CD efficiencies as a function of radius.For example for NH 2 D, the dynamic model predicts a CD efficiency of 4.1% at the center of the core, and 4.9% at the edge, at a simulation time of 5 × 10 5 yr.A detailed account of the new dynamic CD model will be given in Riedel et al. (2023).
We employ the same initial abundances as in Sipilä et al. (2019a), not recounted here for brevity, and most of the simulation parameters are also the same as in that paper.There are two notable differences: here we assume an external visual extinction A V = 3 mag following Pineda et al. (2022) (2 mag in Sipilä et al. 2019a), and we use the new description for cosmic-ray induced desorption from Sipilä et al. (2021).
In Figs. 5 to 10, we compare the average fractional abundances, that is, the column density ratios N(mol)/N(H 2 ) towards the density peaks of the cores, and similarly derived average o/p and deuterium fractionation ratios predicted by the chemical models to those derived by fitting spectra observed towards H-MM1 and Oph D. The fractional NH 2 D and NHD 2 abundances relative to H 2 and the abundance ratios derived from observations are column density ratios through the density peak of the 3D model with abundance profiles (b) and (c) provided in Table 3.The ranges indicated with hatched areas correspond to the inverse-variance weighted averages and the standard deviations of the values derived towards H-MM1 and Oph D. For the observed para-NH 3 ( hereafter pNH 3 ) abundance, we adopt the result from Pineda et al. (2022) that the fractional (ortho+para) NH 3 abundance is 2 × 10 −8 (±10%) in the outer core, but decreases with the density when n(H 2 ) > 2×10 5 cm −3 .Pineda et al. (2022) could only derive the pNH 3 abundance, and assumed o/p-NH 3 = 1.
One can see from Figs. 5, 6, and 7 that the two CD scenarios result in large differences in the average NH 3 , NH 2 D, and NHD 2 abundances.The fractional pNH 3 , NH 2 D, and NHD 2 abundances predicted by the constant CD models are close to Fig. 5. Evolution of the average fractional ortho-and para-NH 3 abundances according to different varieties of our chemistry model: "const CD" means chemical desorption (CD) with a constant (1 %) efficiency, "dyn CD" means changing efficiency of CD depending on the composition of the grain surface, "FS" means full scrambling of protons in proton donation and hydrogen abstraction reactions, and "PH" means proton hop in these reactions (see text).The averages are column density ratios N(mol)/N(H 2 ) taken through the density maximum of the core model used in chemical simulations.The shaded area in the lower panel shows the range of similarly taken average of the fractional pNH 3 abundance through the 3D model of H-MM1 based on the results of Pineda et al. (2022).3).
those derived from the observations, whereas the dynamic CD models seem to produce too much of these molecules.On the other hand, the assumed gas-phase reaction mechanism -either PH or FS -has a minor effect on the abundances, and for most species the predicted abundances from the two models are contained within or lie close to the range of observationally derived values.The only exception is oNH 2 D, for which the model assuming constant CD with PH seems to give a better agreement than the FS model.
It is evident from Fig. 8 that both CD models assuming PH agree with the o/p ratios of NH 2 D and NHD 2 derived from the observations, whereas the FS models under-predict o/p-NH 2 D Fig. 7. Same as Fig. 6 but for ortho-and para-NHD 2 .Fig. 8. Evolution of the average o/p-NH 2 D and o/p-NHD 2 ratios according to our chemistry model.The labelling is the same as in Fig. 5.The ranges derived from observations using the 3D models for H-MM1 and Oph D with abundance profiles (b) and (c) (see Table 3) are indicated with shaded areas.5.In the top panel, the shaded area indicates the ratio derived from observations towards H-MM1, where the fractional pNH 3 abundance has been determined by Pineda et al. (2022).In the bottom panel, the "observed" ratios include both H-MM1 and Oph D. ratio by approximately 30% and over-predict the o/p-NHD 2 ratio by a similar amount.Including results for both H-MM1 and Oph D, the observationally derived ratios are approximately o/p-NH 2 D ∼ 2.9 ± 0.2 and o/p-NHD 2 ∼ 2.1 ± 0.1.As expected, the PH models predict statistical spin ratios.In the FS models, the o/p-NH 2 D ratios settle around 1.9-2.0,and the o/p-NHD 2 ratios are found in the range 2.6-2.8 at late times.According to Fig. 10, the NHD 2 /NH 2 D fractionation ratio is reproduced by the constant CD PH model at simulation times around 5×10 5 yr, and also the constant CD FS model approaches the observed values at late times of the simulation.The constant CD PH model over-predicts the NH 2 D/pNH 3 ratio except at late times (> 5.5 × 10 5 yr), whereas the prediction of the constant CD FS model agrees with the observed ratio only at early times (< 4 × 10 5 yr).The dynamic FS models seem to work better for this ratio.
To summarize, the models assuming constant CD (with PH or FS) come closer to the observed abundances than the dynamic CD models.The models assuming PH in the gas-phase reactions (combined with constant or dynamic CD) can reproduce the observed spin ratios, whereas the FS models predict too low o/p-NH 2 D ratios and too high o/p-NHD 2 ratios.Furthermore, the NHD 2 /NH 2 D ratio can be reproduced by the constant CD PH model at a certain time interval (∼ 3 − 6 × 10 5 yr), but this model over-predicts the NH 2 D/pNH 3 ratio before the simulation time ∼ 5.5 × 10 5 yr, reaching the observed ratio only later.Both deuterium fractionation ratios predicted by the constant CD PH model are within the observed range only during a short period before 6 × 10 5 yr.
In an attempt to quantify the overall agreement or disagreement of the four models with the observations, we present in Table 4 the relative deviations of the predicted NH 2 D/pNH 3 and NHD 2 /NH 2 D ratios, o/p ratios, and fractional NH 2 D and NHD 2 abundances from the values derived from observations.The relative deviation of the quantity q is defined by ∆q ≡ |q mod − q obs |/q obs .The modelled values are taken at the time 5 × 10 5 yr.The observed values are variance-weighted averages of those derived towards the centres of H-MM1 and Oph D. One can see that both PH models give small relative deviations for the o/p ratios, but the constant CD PH model provides the smallest mean deviation when all quantities are included.As discussed above, the dynamic CD models over-predict the fractional abundances, and this causes the large relative deviations in the last two rows of Table 4 for these models.As can be expected from Figs. 6, 7, 8, and 10, the relative deviations change with time, especially at early times (< 3 × 10 5 yr), except for the o/p ratios predicted by the PH models which are nearly constant.The constant CD PH model gives the lowest mean deviation at all times.We note that in general, the abundances of NH 2 D and NHD 2 , and thus the deuterium fractionation ratios, depend on the physical conditions (temperature and density), whereas the spin ratios are directly related to the chemical reaction mechanisms.
Spectra calculated from the 3D cloud model for H-MM1 with abundances interpolated from the 1D chemical model for three simulation times are shown in  and 7, abundances from the dynamic CD models overpredict the line intensities of all the observed lines, and the difference is particularly stark for the NHD 2 lines.The simulated spectra, which depend on the beam-averaged column densities towards the centre of the core model, support the idea conveyed by Table 4 that the constant CD model with PH agrees better with the observations than the other three chemical models tested here.

Discussion
The observations and their analysis presented here indicate that the o/p ratios of NH 2 D and NHD 2 in two cold, dense cores deviate at most 20% from their high-temperature statistical values o/p-NH 2 D = 3 and o/p-NHD 2 = 2.The result conforms with previous observational determinations, with a higher accuracy than previously attained.Statistical spin ratios do not agree with previously published astrochemical models that assume complete scrambling of light nuclei in gas-phase reactions.These typically predict an o/p-NH 2 D ratio of ∼ 2 or slightly below (Sipilä et al. 2015;Hily-Blant et al. 2018).For o/p-NHD 2 , the predictions diverge depending on some differences in complexforming reactions.In the model of Hily-Blant et al. (2018), o/p-NH 2 D remains below 3, whereas it approaches 4 in the original spin-state chemistry model of Sipilä et al. (2015); see Sect. 3.3.4 and Fig. 6 of Hily-Blant et al. (2018).
Using chemical models that include both gas-phase and grain-surface reactions, and two different chemical desorption scenarios, we tested two suggestions that have been put forward to explain the discrepancy between observations and models: 1) gas-phase ammonia has a substantial contribution from grains where H and D addition reactions result in statistical spin ratios (Daniel et al. 2016), and 2) gas-phase reactions that form ammonia and its deuterated isotopologues occur through proton or atom hops instead of particle exchange in a reaction complex (Harju et al. 2017;Sipilä et al. 2019a).According to nuclear spin selection rules, reactions where particles are added one by one should produce species in proportion to their nuclear spin statistical weights.
It turns out that while the adopted chemical desorption mechanism (constant or dynamic) affects strongly the average abundances of ammonia and its deuterated isotopologues (Figs. 5,6,and 7), the effect of this mechanism on their spin ratios is less marked (Figs. 8 and 9).Switching from constant to dynamic CD increases the fractional NH 3 , NHD 2 , and NHD 2 abundances by a factor of 2-3.In the FS scenario, dynamic CD decreases the o/p-NH 2 D and o/p-NHD 2 ratios by ∼ 10 %, and increases the o/p-NH 3 ratio by a similar amount, but has hardly any effect on the spin ratios predicted by the PH model.The poor agreement of the dynamic CD model with the observed fractional abundances suggests that desorption is not as efficient as assumed there, rendering the constant CD model the more realistic of the two.
An input of ammonia from grains seems plausible in view of the fact that the fractional abundance of NH 3 residing on grains is probably several orders of magnitude higher than that in the gas in starless cores (e.g., Tinacci et al. 2022).The recent quantum mechanical simulations of Ferrero et al. (2023) indicate that direct chemical desorption of ammonia from water ice is unlikely, because the reaction energy of hydrogenation is efficiently absorbed by the ice surface.However, the situation may be different in CO-rich ices, which are supposed to be characteristic of dense, starless cores, where also deuterated ammonia is found.The dynamic CD model discussed above is based on these ideas, but it seems to over-predict the desorption efficiency from COrich ices.
Our modelling results suggest that the observed statistical o/p ratios of NH 2 D and NHD 2 are primarily determined by gas phase reactions and are not an effect of vigorous desorption from grains.Firstly, the spin ratios depend strongly on the assumed gas-phase proton exchange scenario (FS or PH) also in the dynamic CD model with very efficient desorption.Secondly, the predicted deuterium fractionation ratios in the grain mantles (NH 2 D/NH 3 < 0.1 and NHD 2 /NH 2 D ∼ 0.1) are lower than those observed in starless dense cores, including the targets of the present study.The observed fractionation ratios agree with the predictions for gas-phase formation.Inspection of the rates of the principal reactions producing and destroying ammonia in a given simulation cell at a certain time shows that desorption and accretion typically account for a few percent of the production and destruction of ammonia at densities (∼ 10 5 cm −3 ) where we believe the observed lines to originate, while the formation through NH + 4 + e − and destruction through NH 3 + H + 3 (and isotopologues), and to a lesser extent through charge exchange between NH 3 and H + , dominate (see also Fig. 2 in Sipilä et al. 2015).This means that, according to our simulations, ammonia and its deuterated forms are rapidly re-processed by ionmolecule reactions in the gas phase, even though desorption from grain surfaces may have a significant contribution to their total abundances.By re-processing we mean constant conversion of nitrogen hydrides to cations through charge or proton exchange reactions, and their return to neutral species in recombination reactions with electrons, for example, NH 3 Besides the fact that many uncertainties are involved in chemical modelling, also observational determination of molecular abundances is non-trivial.There, we are dealing with changing abundances and excitation conditions along the line of sight.To demonstrate this, we show in Fig. 11 the fractional abundances of NH 3 , NH 2 D, and NHD 2 as functions of the radial distance from the core centre and the density, predicted by the constant CD model with PH at a certain time of the simulation.When comparing the best-fit abundances derived using the 3D model of H-MM1 to the model predictions, we used 1) column density ratios N(mol)/N(H 2 ) along a line going through the core centre (Figs. 5 to 10 and relative deviations given in Table 4), and 2) simulated and observed spectra towards the core centre, sampling a cylindrical region including the densest parts of the core (Figs.C.1 to C.4).In both comparisons, the constant CD model with PH gave the best overall agreement with the observations.The interpretation of the four spectral lines used here is facilitated by the fact that they lie close in frequency, meaning that they can be observed in the same spectral band and have similar critical densities.
The PH scenario was previously tested by Sipilä et al. (2019a) with a core model resembling H-MM1.They found that PH reactions favour the singly deuterated NH 2 D at the cost of NHD 2 and ND 3 , but also that all the predicted deuterium fractionation ratios, NH 2 D/NH 3 , NHD 2 /NH 2 D, and ND 3 /NHD 2 , were higher than observed in H-MM1 when the ammonia production only involved proton hop reactions.Here, we do not find large discrepancies between the modelled and observed NH 2 D/pNH 3 and NHD 2 /NH 2 D ratios when using the constant CD PH model.The main difference between the present and the former modelling efforts is that Sipilä et al. (2019a) did not include CD.The replenishment of gas by ammonia formed on grains can diminish the deuterium fractionation as the NH 2 D/NH 3 and NHD 2 /NH 2 D ratios on grains are predicted to be substantially lower than those in the gas.For example, in all four models examined here, the NH 2 D/NH 3 ratio on grains always stays below 0.1.

Conclusions
We have determined the o/p ratios of NH 2 D and NHD 2 in two prestellar cores by observing the ground-state rotational lines of these molecules with the LAsMA multi-beam receiver on APEX.By simulating line emission from realistic 3D models of the cores, and using different assumptions regarding the abundance distributions, we could establish that the o/p ratios deviate at most 20% from their statistical values, o/p-NH 2 D = 3, o/p-NHD 2 = 2.
We also ran chemical models to test if the observed statistical nuclear spin ratios can be explained by enhanced chemical desorption (CD) in the core, employing two different scenarios for proton/deuteron exchanges in gas-phase chemical reactions.An increase in the efficiency of CD in the core interior can be envisioned because the grain surface composition is expected to change from H 2 O-dominated to CO-rich ice.However, models assuming constant CD, with 1% of exothermic reactions leading to desorption, can reproduce the observed abundances better than the so called dynamic CD models with changing efficiency.Despite diverging predictions for the fractional abundances, both constant and dynamic CD models assuming single-particle hop in proton donation and hydrogen abstraction reactions gave the better agreement with the observed o/p ratios than the corresponding models with full scrambling.
The efficiency of CD naturally affects the abundances of NH 3 and its deuterated isotopologues in the gas, but according to our modelling results the effect of the proton exchange scenario in gas-phase reactions is more important to the spin ratios.This is because the accretion and desorption processes are outcompeted by rapid conversion of gas-phase nitrogen hydrides to cations through charge or proton exchange reactions, and the return of cations to neutral species in recombination with electrons.

Fig. 1 .
Fig. 1.Footprint of the APEX/LAsMA array on the 8 µm surface brightness images of H-MM1 (left) and Oph D (right) observed by the Spitzer satellite.The dashed and solid circles on the right represent two orientations of the array in observations towards Oph D. The positions of the central pixels are RA 16 h 27 m 58 s .5,Dec. −24 • 33 ′ 40 ′′ (J2000.0)for H-MM1, and RA 16 h 28 m 30 s .2,Dec. −24 • 18 ′ 42 ′′ for Oph D.

Fig. 3 .
Fig. 3. Excitation temperature profiles of the ground-state lines of ortho and para NH 2 D and NHD 2 in a core model resembling H-MM1.
The sequence is initiated by NH 3 + H 2 D + and leads to NH 2 D through dissociative electron recombination of NH 3 D + .Each reaction step has several possible products, but here we concentrate on the branch ending with NH 2 D. The spinseparated sequences according to the PH and FS scenarios are shown in Fig. B.1.

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

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

Fig. 10 .
Fig.10.Evolution of the average deuterium fractionation ratios NH 2 D/pNH 3 and NHD 2 /NH 2 D according to the chemistry model.The labelling is the same as in Fig.5.In the top panel, the shaded area indicates the ratio derived from observations towards H-MM1, where the fractional pNH 3 abundance has been determined byPineda et al. (2022).In the bottom panel, the "observed" ratios include both H-MM1 and Oph D.
Figs. C.1 to C.4 in Appendix C. The interpolation uses fractional abundances as functions of density in the 1D model.The spectra observed towards the centre of H-MM1 are superposed on the simulated spectra.The abundances from the constant CD model with FS can approximately reproduce the pNH 2 D and pNHD 2 lines, but predict too weak oNH 2 D lines and too strong oNHD 2 lines (Fig. C.1).The constant CD model with PH (Fig. C.2) reaches rather a good agreement in all four spectra at late times.As expected from Figs. 6

Fig. 11 .
Fig. 11.Fractional abundances of the ortho and para modifications of NH 3 , NH 2 D, and NHD 2 as functions of radial distance from the centre (left) and density (right) in the 1D cloud model used for chemical calculations.The chemistry model assumes constant CD with PH.The abundances are taken at the simulation time 5 × 10 5 yr.

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

Fig
Fig. B.1.Spin-state branching ratio diagrams for the reaction sequence NH 3 H 2 D + → NH 3 D + e − → NH 2 D in the PH (left) and FS (right) models.Assuming that o/p-NH 3 = 1, the PH scenario leads to o/p-NH 2 D = 3, whereas the FS scenario results in o/p-NH 2 D = 27/13.

Fig. C. 1 .
Fig. C.1.Simulated NH 2 D and NHD 2 spectra towards the centre of the 3D model for H-MM1 with abundance distributions interpolated from the chemistry model assuming constant CD and full scrambling (FS, red lines).The spectra are convolved with the APEX beam.Spectra observed with APEX towards the centre of H-MM1 are shown with blue lines.

Table 2 .
Line parameters and column densities from Gaussian fits to the hyperfine structure.
* Column density using optically thin approximation with T ex estimated from the ortho line (see text).maps of L1688 from Herschel are accessible via http://gouldbeltherschel.cea.fr.

Table 4 .
Relative deviations, |q mod − q obs |/q obs , of the predicted abundance ratios from the observed values.
†The deviations are given for the simulation time 5 × 10 5 yr.