Collisional excitation of doubly and triply deuterated ammonia NDH and ND by H.
Abstract
The availability of collisional rate coefficients is a prerequisite for an accurate interpretation of astrophysical observations, since the observed media often harbour densities where molecules are populated under non–LTE conditions. In the current study, we present calculations of rate coefficients suitable to describe the various spin isomers of multiply deuterated ammonia, namely the NDH and ND isotopologues. These calculations are based on the most accurate NH–H potential energy surface available, which has been modified to describe the geometrical changes induced by the nuclear substitutions. The dynamical calculations are performed within the close–coupling formalism and are carried out in order to provide rate coefficients up to a temperature of = 50K. For the various isotopologues/symmetries, we provide rate coefficients for the energy levels below 100 cm. Subsequently, these new rate coefficients are used in astrophysical models aimed at reproducing the NHD, NDH and ND observations previously reported towards the prestellar cores B1b and 16293E. We thus update the estimates of the corresponding column densities and find a reasonable agreement with the previous models. In particular, the ortho–to–para ratios of NHD and NHD are found to be consistent with the statistical ratios.
keywords:
molecular data – molecular processes – scattering1 Introduction
Ammonia is an ubiquitous molecule in space, which has been observed in a large variety of media. In some of these places, the temperature can fall to a few kelvins, as in low mass star forming regions. Others, as for example shocked media, have temperatures as high as a few thousand kelvins. Depending on the temperature, the chemistry that leads to ammonia formation can follow different paths. Given the importance of this molecule, in particular in the context of nitrogen chemistry, much effort has been made in the past in order to characterize the corresponding formation or destruction mechanisms. This proved to be successful and nowadays, most of the observations can be, at least qualitatively, understood. In some cases, however, astrochemical models and observations still show discrepancies, as in the case of circumstellar envelopes where the models under–predict the observed molecular abundances by a few orders of magnitude (menten2010). Additionally, in recent years, it was found that the nuclear–spin states of a given molecule are not necessarily populated according to simple statistical considerations. This is the case for ammonia, for which some observations report ortho–to–para ratios (OPR) in the range 0.5–0.7 (hermsen1988; persson2012), while a statistical ratio of 1 was expected in the past, for a formation in the gas phase. If ammonia is formed in the solid phase, and subsequently released from the grains, this ratio would be 2. Recently, astrochemical gasphase models have been updated by including rigorous nuclearspin selection rules relevant to describe the ammonia formation. It was found that ratios lower than 1 can be well accounted for by such models in the low temperature regime, where the H gas is paraenriched (rist2013; faure2013; sipila2015).
While NH is observed in media that span a wide range of physical conditions, the detection of its multi–deuterated counterparts, i.e. NDH and ND, is restricted to cold regions. The reasons behind this are now well established and lie in the fact that the deuterium enrichment comes from the difference of zero point energies associated with the hydrogenated or deuterated isotopologues. These differences create a dichotomy in the chemical network. Indeed, the reactions that lead to deuterated species have, in general, exothermicities which are enhanced by comparison to the equivalent reactions that involve their hydrogenated counterparts. This effect creates a significant deuterium enrichment at low temperatures if H is mainly in its para form. Indeed, in the high temperature regime, the dichotomy between deuterated and hydrogenated molecules vanishes, since the difference of exothermicities become negligible with respect to the kinetic temperature. In the particular case of ammonia, the main formation route consists of a sequence of reactions involving H, which is initiated by the N + H NH + H reaction and ends with the dissociative recombination reaction NH + NH + H (see e.g. dislaire2012). The initiating reaction has a small endothermicity (in the range 100–300K, see e.g. gerlich1993). On the other hand, the endothermicity of the equivalent reaction involving HD and leading to ND is reduced (20–70K, see marquette1988; sunderlin1994), which hence favors the deuteration processes at low temperatures. Finally, the abundances of the NH isotopologues decrease with the number of deuterium atoms. Given the physical conditions characteristics of low–mass star forming regions, which harbour low gas temperatures (T 20K) and relatively low gas densities (n(H) cm), the low abundance implies that the signal associated with these isotopologues is weak, which makes the NDH and ND molecules hard to detect with the current facilities. Hence, so far, the detection of all four deuterated istopologues of ammonia is restricted to two astrophysical objects, the prestellar cores 16293E^{1}^{1}1This prestellar core is sometimes referred as L1689N in the literature and Barnard 1 (roueff2005). This census has still to be enlarged to other objects to ensure statistically meaningful conclusions, starting with objects where ND has already been observed, as in NGC 1333 (vandertak2002).
The current work is aimed at providing collisional rate coefficients for the NDH and ND molecules, with p–H as a collisional partner. For these molecules, earlier calculations were performed with He or H with a reduced basis set (machin2006; machin2007; yang2008; wiesenfeld2011). Such rate coefficients are essential in order to interpret the emission observed in astrophysical objects and imply solving a problem of quantum dynamics. Recently, similar calculations (ma2015; tkac2015) were performed independently for NH and ND and using the Potential Energy Surface (PES) reported by maret2009, which is also used in the current work. However, rather than providing collisional rate coefficients, the focus of ma2015 or tkac2015 was to assess the possibility of testing the accuracy of the PES by comparing theoretical predictions based on the PES with experimental data. The comparisons performed so far show a good agreement between experiments and theory at collision energies above 430 cm (tkac2015).
The present paper is organized as follows. The potential energy surface is described in Section 2 and the collisional dynamics based on this surface in Section 3. We then describe the rate coefficients in Section 4. In Section 5, we discuss the current results with respect to other related collisional systems. In Section 6, we use our new rate coefficients in order to re–interpret the observations previously reported toward 16293E and Barnard 1. Finally, we present our conclusions in Section 7.
2 Potential Energy Surface (PES)
The collisional excitation of an isotopic homologue takes place on the same BornOppenheimer PES as the main isotopologue. All differences are therefore contained in the dynamical treatment. Within the rigidrotor approximation, these differences involve changes in i) the internal stateaveraged geometry, ii) the center of mass position, iii) the reduced mass of the total system, iv) the energy level spacing and v) the rotation of the principal axes of inertia. All these effects were considered here for NDH and ND, except the change of internal geometry: the averaged geometry of NH was employed for both deuterated isotopologues, thereby neglecting the deuterium substitution effect on the NH bond length and angle. This assumption was previously adopted for the NHDH and NDHH systems by daniel2014 and wiesenfeld2011, respectively. The (rigidrotor) PES for NDH was thus again derived from the CCSD(T) NHH PES computed by maret2009, but in the principal inertia axes of the isotopologue. We note that internal geometry effects in the estimate of the rate coefficients are expected to be only moderate ( 2030%) at the temperatures investigated here, i.e. K, as demonstrated by scribano2010 for the DOH system.
The transformation between the reference frame of an isotopologue and the original frame of the NHH system can be found in wiesenfeld2011 (Eqs. 23) where the angle is zero for ND and the coordinates of the ND center of mass in the NHH reference frame are =0.0 bohr and =0.088548436 bohr (where and were calculated for the average geometry of NH used by maret2009). As for NHDH and NDHH, the NDH PES was generated on a grid of 87 000 points consisting of 3000 random angular configurations combined with 29 intermolecular distances in the range 315 bohr. This PES was expanded in products of spherical harmonics and rotation matrices (see Eq. 4 in Wiesenfeld et al. 2011), using a linear leastsquares fit procedure detailed in rist2011. The fit actually employs the 120–terms expansion optimized for NHH by maret2009, where details can be found. The final expansion thus includes anisotropies up to =11 for ND and =4 for H. The root mean square residual was found to be lower than 1 cm for intermonomer separations larger than 4.5 bohr, with a corresponding mean error on the expansion coefficients smaller than 1 cm.
3 Collisional dynamics
For the various symmetries of NDH and ND^{2}^{2}2The NDH molecule is an asymmetric top which can be described by the quantum numbers , and . In what follows, we use the pseudo quantum number =  and adopt the notation to describe the rotational structure. The ND molecule is a symmetric top which is described by the and quantum numbers and we hence label the levels as . For both molecules, the umbrella vibrational motion further splits every rotational level in two energy levels. We then introduce the quantum number in order to distinguish the levels in every doublet. With , the NDH and ND energy levels are then respectively labeled as and (see details in appendix A). the calculations aim at providing rate coefficients up to = 50 K and for rotational energy levels below 100 cm. More precisely, in the case of NDH, we thus provide rate coefficients for the levels up to . For the para symmetry of ND, we consider the levels up to = and for the ortho symmetry, up to = (see Table 1).
Depending on the deuterated isotopologue, we either solved the collisional dynamics using the MOLSCAT^{3}^{3}3J. M. Hutson and S. Green, MOLSCAT computer code, version 14 (1994), distributed by Collaborative Computational Project No. 6 of the Engineering and Physical Sciences Research Council (UK) computer code or HIBRIDON^{4}^{4}4HIBRIDON is a package of programs for the timeindependent quantum treat ment of inelastic collisions and photodissociation written by M. H. Alexander, D. E. Manolopoulos, H.J. Werner, B. Follmeg, Q. Ma and P. J. Dagdigian, with contributions by P. F. Vohralik, D. Lemoine, G. Corey, R. Gordon, B. Johnson, T. Orlikowski, A. Berning, A. DegliEsposti, C. Rist, B. Pouilly, G. van der Sanden, M. Yang, F. de Weerd, S. Gregurick, J. Klos and F. Lique. More information and/or a copy of the code can be obtained from the website http://www2.chem.umd.edu/groups/alexander/hibridon/hib43. package which includes the treatment of both inversion motion and non–spherical linear collider . The former was used to treat the NDHH system, since we neglected the inversion motion for this particular isotopologue. The inversion motion splits every rotational level in two rotation–inversion energy levels, the splitting being of the order 0.4 cm for NDH. However, the two levels of every doublet are associated to different NDH nuclear–spin symmetries, either the ortho or para species, which are not connected through inelastic collisions. Hence, by neglecting the inversion motion we consider a system which is suitable to describe both the ortho and para species, if we omit the small error (0.4 cm) made in the estimate of the energy of the rotation–inversion state. The methodology of the corresponding calculations is similar to the recent study of the NHDH collisional system (daniel2014). In the case of NDH, we followed the description adopted in wiesenfeld2011, which is based on the spectroscopic parameters given by coudert1986; delucia1975. Hence, we described the energy levels adopting the rotational constants A = 5.3412, B = 7.4447 and C= 3.7533 cm and the centrifugal terms = , = and = , these values being given in the II representation (gordy1984) suitable to MOLSCAT.
In the case of ND, two levels of a doublet can pertain to the same symmetry and thus can be connected through collisions. The MOLSCAT code could be used to treat such a problem, but only in the case of a spherical collider, like He or p–H restricted to its fundamental rotational level . Since our calculations aim at including the first excited state of p–H (i.e. =2), we used the HIBRIDON package which includes the possibility of treating the inversion motion in the case of a non–spherical collider. We checked the consistency of the results obtained with the two codes, for the collisions of o–ND with H () and we could verify that the cross sections agree with an accuracy typically of the order of the per cent.
o–ND  p–ND  m–ND  

8.27  0.05  0.00  
8.32  10.29  10.34  
22.78  30.91  30.86  
22.83  43.55  43.55  
28.84  43.60  43.60  
28.89  61.71  61.76  
53.64  84.69  84.69  
53.69  84.74  84.74  
59.69  102.90  102.85  
59.74  136.12  136.12  
70.56  136.17  136.17  
70.61  143.34  143.34  
94.78  143.39  143.39  
94.83  154.28  154.33  
100.84  197.83  197.83  
100.89  197.88  197.88  
103.83  215.34  215.34  
103.88  215.39  215.39  
121.99  216.04  215.99  
122.04  269.83  269.83 
To define the ND energy structure, we adopted the rotational constants and cm and set the inversion splitting at a constant value of 0.05 cm. The corresponding energy levels are given in Table 1. As can be seen in this table, the energy levels associated with the para and meta symmetries of ND are mostly identical. The only differences between the two sets of energies are found for the =0 levels. This similarity entails that the cross sections associated to the two symmetries will be roughly identical for a given transition. Indeed, the coupling terms that enter the dynamic equations are identical for the two symmetries and the only differences in the two sets of equations will thus come from the differences in the energy levels. These similarities are verified in our state–to–state cross sections performed at a few energies. Cross sections involving levels with both initial and final differ at most by 1%. For transitions that involve levels with , the differences may be higher but never exceed 5%. Hence, we chose to perform calculations for only one symmetry, i.e. for p–ND. Full details on the spectroscopy of ND are presented in Appendix A.
The accuracy of the cross sections is assessed by a few test calculations where we checked the convergence against the step of integration used to propagate the wavefunctions, and the number of energy levels considered during the propagation. These parameters are summarized in Table 2 for NDH and Table 3 for ND. In the case of ND, the integration step is not reported since we used a single value of 0.1 at every energies. The parameters quoted in these tables ensure a convergence better than 5% at all energies. In these tables, we also report the step between two consecutive energies. These values ensure a good characterization of the resonances, for the transitions between levels under consideration. For both NDH and ND, all the calculations were performed within the close–coupling formalism. Recently, tkac2015 reported quantum dynamical calculations for the NDH collisional system, with the aim of testing the PES against experimental differential cross sections. These comparisons were performed at total energies above 430 cm. They assessed the convergence of their calculations with respect to the size of the rotational basis of the two rotators and, in particular, they found that omitting the energy level of p–H can lead to misestimate the cross sections by 25%, for transitions associated with the H state–to–state transitions. Hence, we checked the accuracy of our calculations by testing the influence of the size of the H rotational basis, at various total energies between 200 and 400 cm. We find that below 400 cm, the cross sections were affected by less than 5% due to the truncation of the p–H rotational basis to . At 400 cm, we indeed found that some transitions were more largely affected, with inaccuracies of up–to 35%. However, less than 5% of the total number of transitions are affected by errors greater than 5%. Moreover, the corresponding cross sections are the lowest in magnitude, i.e. lower than 0.001 , the highest cross sections being of the order of 10 at this total energy. Hence, given the low influence of the energy level of p–H for the range of energy considered in the current work, and since the inclusion of this level would considerably raise the calculation time, we chose to restrict the rotational levels of H to = 0, 2 .
Energy range (cm)  step in energy (cm)  STEPS  JMAX 

25  0.1  40  7 
25 60  0.1  20  7 
60 80  0.1  10  7 
80 110  0.1  10  8 
110 155  0.2  10  8 
155 250  0.5  10  8 
250 400  1  10  8 
Energy range (cm)  step in energy (cm)  JMAX  E 

30  0.1  6  + 
25 65  0.1  7  + 
65 105  0.1  8  650 
105 150  0.1  9  700 
150 200  0.1  9  750 
200 300  0.5  10  800 
300 450  2.0  10  800 
4 Rate coefficients
The deexcitation rate coefficients for the rotationinversion transitions were calculated by averaging the above cross sections with a Maxwell–Boltzmann distribution that describes the distribution of velocity of the molecules in the gas (see e.g. eq. (5) in wiesenfeld2011). The rate coefficients were thus obtained for temperatures in the range 5–50 K and for all transitions involving NHD levels below 82 cm (i.e. up to ) and ND levels below 103 cm (i.e. up to and for the para and ortho species, respectively).
Since the hyperfine structure due to the nitrogen nuclei (N) is resolved observationally (see Section 6), rate coefficients for the hyperfine transitions are also required. The recoupling theory now routinely employed for linear species (see e.g. daniel2004; daniel2005; kalugina2015) has not yet been extended to nonlinear molecules for which only the statistical approximation is possible. In this approach, it is assumed that the hyperfine deexcitation rate coefficients are proportional to the degeneracy (2 + 1) of the final hyperfine level and completely independent of the initial hyperfine level. This corresponds to a statistical reorientation of the quantum number after collision and it is also called the randomizing limit or proportional approach. It has been shown to be inaccurate by large factors when predicting hyperfine individual rate coefficients (see e.g. daniel2005; faure2012). On the other hand, when the opacity of the individual hyperfine lines is low (1), the hyperfine levels are populated according to their statistical weights and, in this regime, the statistical approximation is well adapted (daniel2005; faure2012). As the hyperfine transitions of the deuterated isotopologues are usually optically thin, the statistical approximation should be accurate enough to model the spectra of NHD, NHD and ND. Note that in radiative transfer models, the quasi–elastic rate coefficients are also needed. Unfortunately, these rate coefficients cannot be guessed from elastic rate coefficients, which are typically of the order cm s. Indeed, faure2012 showed that applying the proportional approach to elastic rate coefficients lead to overestimate the quasi–elastic rates by large factors, with respect to rate coefficients calculated with a recoupling method. Hence, in the radiative transfer calculations described in Sect. 6, we applied the proportional approach to a canonical rate of cm s in order to obtain an estimate of the quasi–elastic rate coefficients.
The whole set of hyperfine rate coefficients will be available through the LAMDA (schoier2005) and BASECOL (dubernet2013) databases.
5 Discussion
To date, two sets of rate coefficients were computed for NDH. The first set considered He as a collider (machin2007), while the second set was calculated for p–H using a basis set (wiesenfeld2011). Quite often, the collisions with He are used to emulate the rate coefficients with p–H. However, as for other molecular systems, the comparison between these two sets showed that the differences can be quite large, the rates with p–H(=0) being higher by a factor 3–30 depending on the transition.
In order to perform some comparison with these previous works, we considered the NDHH rate coefficients determined by wiesenfeld2011, where data were computed for the nine first rotational energy levels and for temperatures in the range 5–35K. These calculations only differ from the current ones in two points. First, in wiesenfeld2011, the cross sections were computed in part with the coupled–states approximation. In the current study, all the calculations were performed within the close–coupling formalism. Secondly, the present H rotational basis includes the level, while it was reduced to in wiesenfeld2011. As already commented in daniel2014, increasing the H basis can have a non negligible effect on the estimate of the magnitude of the rate coefficents. In the case of NHD, it was found that using a basis rather than a basis can lead to underestimating by a factor 3 the cross sections of highest magnitude (see Fig. 2 in daniel2014). Moreover, in that study, it was already mentioned that the rates given by wiesenfeld2011 should be accurate to within a factor 3. Such a factor was guessed from a few preliminary calculations and considering the behaviour of a reduced set of cross sections (see Fig. 9 in daniel2014). With the current data, the comparison can be extended and done on more quantitative grounds. In Fig. 1, we show the ratio between the two sets of rate coefficients as a function of the magnitude of the rate coefficients. From this figure, it can be seen that the largest differences are found for the rates of highest magnitude. This is similar to what was already observed in the NHDH collisional system. Additionally, the alteration of the rates of highest magnitude is at most of a factor 2.5, which is of the same order of magnitude as for the NHDH collisional system.
The calculation of NDHe rate coefficients was a goal of the PhD thesis of L. Machin (machin2006)^{5}^{5}5https://tel.archivesouvertes.fr/tel00128133 (these unpublished rate coefficients are given in Table 5.3 of the thesis manuscript). As already evidenced for the NHD and NHD isotopologues, the differences between the He and H rate coefficients are large (see Fig. 8 in daniel2014). In Fig. 2, we plot the ratio of the H and He rate coefficients, at a temperature =10K and for the levels up to . As for NHD and NHD, the ratios are in the range 1–100 and the largest differences are obtained for the highest rate coefficients. Interestingly, the calculations of Machin with He were performed for both the para and meta symmetries (however, the meta symmetry was misleadingly referenced as ortho in that work). Machin could thus check that the rate coefficients for the para and meta species agree with a good accuracy, as expected. In fact, he found that at temperatures higher than 50K, the agreement is better than a few per cent. The comparison at =10K shows larger differences, but most of the rate coefficients agree to within 20% (in Table 5.3 of the manuscript, the transition is the only transition that shows a higher ratio, i.e. 2). These differences in the para and meta rate coefficients were attributed to differences in the resonance structure. In Sect. 3, we mentioned that we performed some test calculations at some specific energies and we could verify that at these energies, the meta and para cross sections agree within a few per cent. However, as in the case of the collisions with He, we cannot discard that the resonance structure for the two symmetries could induce larger differences, but not above 20%, which is the typical accuracy of the present rate coefficients.
Recently, ma2015 reported dynamical calculations for the NDH system that were performed on the same PES as the current calculations. These two sets of calculations were performed independently and with a different focus. While in the current case, we are interested in providing collisional rate coefficients, ma2015 put emphasis on the description of the resonance features in order to probe the possibility of characterizing these resonances in future experiments. Since the energy, strength and shape of the resonances depend strongly on the accuracy of the PES, this would provide a test of the accuracy of the PES of maret2009 and hence, would give an estimate of the accuracy we could expect for the rate coefficients. In Fig. 3, we give the cross sections of o–ND that originate from the level. This figure is similar to Fig. 13 of ma2015 and in the current figure, the cross sections obtained by ma2015 are reported as dotted lines. The overall agreement between the two sets of calculations is very good, even if there remains some small differences in the description of the resonances. These differences do not exceed a few per cent and presumably originate from the respective accuracy of the two collision dynamics calculations .
6 Astrophysical modeling
The rate coefficients described in the previous sections are a key ingredient for astrophysical models interpreting the deuterated ammonia inversion/rotational lines observed from far–infrared to centimeter wavelengths. Such observations of the NHD, NDH and ND isotopologues were reported and analyzed by tine2000; roueff2000; saito2000; lis2002; vandertak2002; roueff2005; lis2006; gerin2006, but the lack of rate coefficients hampered the analysis, which thus relied on simplifying assumptions. In some cases, the recourse to such simplification leads to discrepancies, as outlined by lis2006. Indeed, it was then shown that an LTE analysis applied to two rotational lines of o–NDH, at 336 and 389 GHz, gives inconsistent column density estimates. In the case of the observations towards the prestellar cores 16293E or B1b, it was found that the two estimates disagree by factors of 3 and 6, respectively.
In what follows, we reconsider the observations available towards B1b and 16293E, in light of the newly calculated NDH and ND rate coefficients presented here, as well as the NHD / H rate coefficients described in daniel2014. The molecular excitation and radiative transfer is solved with the 1Dart code described in daniel2008. As in daniel2013, we report column densities, which are calculated by counting the number of molecules across the diameter of the spheres that describe the sources. These column densities are not tied to any particular molecular transition and hence or not averaged over a particular telescope beam size. Such beam averaged column densities are often introduced when estimating molecular column densities and in what follows, we will sometimes introduce it for the sake of comparison with previous results.
6.1 Barnard 1b
B1b, located at 235 pc in the Perseus molecular cloud, is a region of active low–mass star formation. The region around B1b has been the subject of many observational studies (see references in daniel2013) and consists of two objects, B1b–N and B1b–S (hirano1999), whose evolutionary states are still under debate. Currently, B1b–N is thought to be in the first hydrostatic core stage while B1b–S seems to be a slightly more evolved object (hirano2014; gerin2015).
6.1.1 NhD
The physical structure adopted for this source (i.e. dust and gas temperatures, H density) is taken from daniel2013, where it was inferred from the analysis of 850 m and 1.3 mm continuum observations. In the latter study, the map of the o–NHD – emission was already analyzed by solving the molecular excitation. However, the rate coefficients for the NHD–H system were not available at that time and the analysis thus relied on the NHD–He rate coefficients (machin2006), that were subsequently scaled. We refer the reader to the discussion in daniel2013 for the details that concern the scaling of the rates. More recently, daniel2014 reported rate coefficients for the NHD–H system. Using these new rate coefficients with the o–NHD abundance profile given in daniel2013, we find that it leads to overestimate the line intensity by a factor 2. This illustrates the inaccuracy introduced by the scaling of the rate coefficients. We thus re–analyze the o–NHD and o–NHD data available towards B1b with the same methodology as in daniel2013. In particular, we note that the B1b is a complex region, with two cores separated by 20”, which have V that differ by 1 km s. Hence, to model the spectra, we follow the methodology described in daniel2013: we use a 1D spherical model, adopt a V6.7 km s and constrain the models by fitting the blue part of the spectra. The estimate of the NHD abundance profile across the core is constrained thanks to the o–NHD – emission map and the comparison between the model and observations is given in Fig. 4. The corresponding o–NHD abundance profile is reported in red in Fig. 5. The confidence zone, as defined in daniel2013, is indicated as grey area. In this figure, the abundance profile derived by daniel2013 is shown in blue. Finally, the o–NHD data observed towards the central position was also reanalyzed. The comparison of the models and observations for the two isotopologues at this position are given in Fig. 6.
The previous estimate for the o–NHD column density made by daniel2013 was log() = . With the new rate coefficients, this value is revised to log() = . For o–NHD, the column density was previously estimated to log() = while the current estimate is . The corresponding isotopologue column density ratio, previously estimated to is now .
In the case of the p–NHD spin isomer, the abundance was constrained using a mini–map of the – line at 110 GHz, observed at the IRAM 30m telescope, and the – line at 494 GHz, observed with Herschel. These latter data were obtained using the HIFI instrument (degraauw2010) on Herschel (pilbratt2010), in the context of the open time program OT2_dlis_3 (Herschel observation ID 1342248917). The data were taken in the frequencyswitching mode, were processed using the standard HIFI data reduction pipeline and the resulting spectra were subsequently reduced using the GILDAS CLASS software package. The FWHM HIFI beam size is 44 at 492 GHz and the main beam efficiency is 62%^{6}^{6}6http://herschel.esac.esa.int/twiki/pub/Public/HifiCalibrationWeb/HifiBeamReleaseNote_Sep2014.pdf.
By using an overall scaling of the oNHD abundance profile, we managed to obtain a reasonable fit of the pNHD observations, although the line shape of the – transition is not correctly reproduced. Indeed, the best fit model predicts self–absorption features which are not observed (see right panel in Fig. 6). Some alternative models were tried in order to improve the fit, but we could not find a model that would considerably improve the comparison with the observations. We thus kept the simplest model, i.e. the model where the p–NHD abundance is scaled from the o–NHD abundance profile. The OPR derived from the modeling is . This ratio is higher than the ratio predicted through statistical considerations, i.e. 3, but, compatible given the large error bars. Additionally, we note that examining the possible scenarios relative to the NHD formation, a ratio of 3 would be expected in the case of formation on grains (with T 10K) while a ratio 3 is expected if the formation occurs in the gas phase (sipila2015). Note, however, that it was shown by persson2015 that, in the case of NH, the interconversion between ortho and para states due to reactive collisions with H atoms is an efficient process that may substantially alter the NH OPR. Such process is, however, generally not included in the gasphase formation models due to the lack of either theoretical or experimental value for the reactive rate coefficients. Finally, to date, OPR higher than 3 are not expected for NHD. Although shah2001 reported similar values, we believe that the accuracy of the analysis is influenced by the high opacity of the line, which introduces uncertainties in the estimate of the abundance profile. This is illustrated by the models shown in Fig. 7. In this figure, we show a model for the 110 GHz line where the p–NHD abundance profile is set identical to the o–NHD profile (blue curve). The corresponding intensities are compared to the model adopted for B1b (red dashed curve). The fact that the intensity does not scale linearly with the column density is emphasized through the comparison with the model where the intensities are scaled according to the OPR (red curve). In particular, it can be seen that the hyperfine satellites are in a regime of opacity where they still scale proportionally to the column density. The central component is however opaque and the two models differ by a factor 2. From this comparison, it can be inferred that if we consider the ratio of integrated intensities of the 86 and 110 GHz lines, this ratio will only give a lower limit to the OPR. For the (0”,0”) position, this ratio is 2.9. Finally, if we only consider the most optically thin hyperfine satellites, their integrated intensity ratio should be closer to the actual OPR. Considering these components, we obtain an integrated intensity ratio 3.50.2. To conclude, the various ways of deriving the OPR of NHD rule out a ratio lower than 3 and therefore favor the grain surface formation processes.
Considering these various estimates and since, to date, the highest value expected for the OPR of NHD is 3, we use this latter value, in what follows, to derive the NHD total column density. The total NHD column densities is thus log() = .
6.1.2 Nhd
Using the same methodology as for NHD, we re–analyzed the o–NHD (single position) observations published in lis2006. The spectroscopy of o–NHD was obtained from the splatalogue database where the hyperfine structure spectroscopy corresponds to the data published in coudert2006; coudert2009. As a first guess, we used the NHD abundance profile discussed in the previous section, and tried to scale it by a constant value in order to reproduce the two observed lines at 336 and 389 GHz. By doing so, we could not find a reasonable model. Indeed, a model that fits the 336 GHz line could be obtained by scaling down the o–NHD abundance by a factor 4.5. However, this would lead to a model where the integrated area of the 389 GHz line would be overestimated by a factor 3. As stated by lis2006, the 389 GHz line has a higher critical density than the 336 GHz line. Hence, in order to simultaneously fit the two o–NHD lines, the NHD abundance profile has to be modified with a global shift of the abundance towards lower densities, i.e. larger radii. Such a model corresponds to the abundance profile given in Fig. 5 (right panel) and the corresponding fit of the observations is shown in Fig. 8. As an alternative to this modification, we checked other possibilities that would enable to define NHD and NHD abundance profiles that would only differ by a scaling factor. We found that assuming a constant gas temperature around 10K would allow such a solution. However, given the constraints obtained on the dust temperature across the B1b region (see daniel2013) and given the high densities in the innermost part of the core, this solution was discarded.
The NHD observations in B1b were already discussed in lis2006. As mentioned earlier, in the latter study, the column density estimated separately from the 336 or 389 GHz lines, using an LTE approximation, were found to disagree by a significant factor. Indeed, these beam averaged column densities were respectively found to be and cm, i.e.log() = and . Here, the o–NHD column density is estimated to log() = . Finally, a previous analysis by lis2006 have shown that in B1b, the OPR is consistent with the ratio of the nuclear spin statistical weights, i.e. 2:1, within error bars. With such a ratio, the total NHD column densities is thus log() = .
6.1.3 Nd
To model the ND transition, we scaled the o–NHD abundance profile by an overall factor. In fact, the exact shape of the abundance profile is not constrained since the observed line is optically thin and because observations are available only towards a single position. As a consequence, alternate abundance profiles would provide similarly good fits to the observations. Hence, the only parameter which is well constrained by the model is the column density. The m–ND column density is estimated to log() = . In a previous study, lis2002 reported a ND column density of cm, ie. log() = . This value corresponds to the total ND column density, i.e. which consider all the spin isomers, and is averaged over the beam size. With the current model, the column density averaged over a beam of 25” lead to a column density for m–ND of cm. Assuming, that the ortho, meta and para spin isomers are populated according to their spin statistical weights, i.e. 16:10:1, we obtain a total column density of cm which is a factor 3.5 higher than the estimate of lis2002. Note however that in lis2002, the column density estimate assumes that the various spin isomers are thermalized at the gas temperature. At 5K and 10K, this implies that 60% and 45% of the molecules are, respectively, in the meta modification. In the current estimate, we assume that the proportion of molecules in the meta modification is independent of the temperature and is set at 37%. Hence, a factor 1.5 in the two estimates comes from the hypothesis respectively made on the partition functions.
Finally, assuming that the ortho, para and meta ND states are populated according to the ratios 16:1:10, we derive log() = .
6.1.4 Discussion
Deviations of the OPR from the statistical weights ratio are expected if the molecules form in the gasphase owing to nuclear–spin selection rules. Nuclear–spin effects in the gasphase have, however, a complex dependence on density, temperature and time (rist2013; faure2013; sipila2015). Here, we have thus assumed that the various nuclear–spin species of all four ammonia isotopologues are populated according to their spin statistical weights, in order to translate the observed column densities to total column densities. In the case of the main isotopologue of ammonia, the p–NH column density was estimated to be log() = (daniel2013), which translates to a total column density of log() = assuming a statistical OPR of 1:1.
[NHD]/[NH]  [NDH]/[NHD]  [ND]/[NHD]  [ND]/[NH]  

old  new  old  new  old  new  old  new  
B1b  0.230.05  0.69  0.150.03  0.16  0.0330.01  0.11  1.10.5  1.30.6 
16293E  0.190.05  …  0.220.05  0.21  0.0240.01  0.05  1.00.5  … 
The fractionation ratios derived from the current analysis are reported in the first entry of Table 4. Note that with respect to the previous estimate reported by roueff2005 and displayed in this table, the current analysis leads to larger ratio, the updated values being higher by factors in the range 1–10. We thus obtain for Barnard 1b NHD/NH0.69, NHD/NHD0.16 and ND/NHD0.11. These ratios provide very useful constraints for astrochemical models and they should help, in particular, to disentangle between a gasphase and a grain surface formation of the ammonia isotopologues. The gasphase model of roueff2015 was shown to reproduce the observed column density ratios of all four ammonia isotopologues in B1b and 16293E within a factor of 3. The updated observed value for NHD/NH would be however underpredicted by a factor of . On the other hand, it is interesting to note that a purely statistical model^{7}^{7}7A purely statistical approach predicts the following ratios: NHD/NH=3D/H; NHD/NHD=D/H; ND/NHD=D/H of grain surface formation (i.e. through statistical hydrogenation/deuteration of solid atomic nitrogen) reproduces the observed ratios within the error bars, provided that the accreting atomic D/H ratio is 0.2. Such a high value for the D/H ratio is actually predicted in regions with high density and heavy depletion, such as prestellar cores (roberts2003). Clearly, more observations of the four ammonia isotopologues, in several other sources, as well as other related species such as ND and NHD, are necessary to make progress.
6.2 16293e
The IRAS 16293E source is a prestellar core located close to the IRAS 162932422 protostar. Hence, its physical and chemical properties are influenced by this source, mainly because of the interaction with the ouflows that emanate from it (castets2001; lis2002b).
To model the emission of the ammonia isotopologues observed towards the 16293E core, we adopt the source structure described in Bacmann et al. (2015). Based on ground–based and space–based observations of the dust continuum emission at various wavelengths, the center of the density profile is assumed to be at the position , (bacmann2015). The reference position of the observations is , (roueff2005; lis2006).
In order to constrain the behaviour of the abundance, we start by considering the o–NHD isotopologue, for which we have maps of the 336 GHz line emission obtained with APEX (gerin2006), as well as CSO observations of the 336 and 389 GHz lines obtained close to the dust emission peak (lis2006). As for the B1b modeling, it is confirmed that these two lines are good probes of the radial variations of the abundance profile, because of their respective critical densities and because of the difference of opacity of the lines.
6.2.1 Nhd
The comparison between the model and observations for the 336 and 389 GHz lines is given in Fig. 10. The map of the 336 GHz line obtained with APEX is reported in Fig. 11. The abundance profile, which corresponds to these models, is given in Fig. 12.
As outlined in gerin2006, the o–NHD and the dust emission peak are offset by 10”. With the current spherical model, such an asymmetry cannot be taken into account and the model shown in Fig. 11 thus only reproduces qualitatively the overall behaviour of the o–NHD emission. Finally, assuming an OPR of 2:1, i.e. given by the nuclear spin statistical weights, we obtain a good fit to the APEX p–NHD map. The comparison between model and observations is reported in Fig. 13. The OPR is similar to the one obtained by gerin2006 where the abundance ratio was derived from the line intensity ratio, which is accurate in the case of optically thin lines.
The o–NHD column density is log() = and with a 2:1 ratio between the ortho and para species, this gives a total column density of log() = .
6.2.2 NhD
The o–NHD (–) line was observed at the single position , , which is offset by 15” from the center of our model. This o–NHD observation can be satisfactorily reproduced by applying an overall scaling factor to the o–NHD abundance profile discussed in the previous section. The comparison between the model and observations is shown in Fig. 14.
The o–NHD column density is log() = and with a 3:1 ratio between the ortho and para species, this gives a total column density of log() = .
6.2.3 Nd
As for NHD, we model the ND transition by applying an overall scaling factor to the o–NHD abundance profile. The m–ND column density is log() = and with a 16:10:1 ratio between the ortho, meta and para species, this gives a total column density of log() = . The comparison between the model and observations is shown in Fig. 15.
6.2.4 Discussion
The column densities derived for the four ammonia isotopologues in 16293E are reported in the second entry of Table 4. Again, our new analysis does not change much the fractionation ratios with respect to the study by roueff2005. We thus obtain for 16293E the ratios NHD/NHD0.21 and ND/NHD0.05, which are in fact very similar to those derived in B1b. Such small variations in the fractionation ratios suggest similar physical conditions in both sources and, therefore, similar chemical pathways. We note that these fractionation ratios can be reproduced within a factor of 3 by the gasphase model of roueff2005, but also by a purely statistical grain surface model. Once again, more observations are required to make further progress in our understanding of the deuterium fractionation of ammonia and to rule out one of the two possible scenario.
7 Conclusions
We used the state–of–the–art potential energy surface that describes the NH–H interaction (maret2009) in order to determine collisional rate coefficients for the NDH and ND isotopologues. We subsequently used the Close–Coupling method to calculate a set of collisional rate coefficients for NDH, which applies to both the para and ortho symmetries of this molecule. For ND, we made specific calculations for the ortho and para species. A few test calculations showed that this latter set was also suitable to treat the meta symmetry of ND.
The NDH–H rate coefficients were compared to earlier calculations based on the same PES but using a reduced basis for H (wiesenfeld2011). We found that the new set shows differences of up to a factor 2.5 with respect to the old data. Moreover, the rates which are the most affected are those of highest magnitude. In the case of ND, we compared our calculations with ND–He calculations (machin2006). As previously found for NDH (wiesenfeld2011) or NHD (daniel2014), we find large differences between the H and He rate coefficients with variations of up to two orders of magnitude.
Finally, we used these new rate coefficients, as well as the NHD–H calculations published earlier (daniel2014), in order to re–interpret the observations available in the literature of the various deuterated isotopologues of NH (roueff2000; lis2002; roueff2005; gerin2006; lis2006). We focused on the B1b and 16293E prestellar cores and we used a comprehensive radiative transfer analysis based on the latest estimates of the physical structure of these objects (daniel2013; bacmann2015). By comparison with the earlier estimates of the column densities, based on approximate methods which did not solve the coupled set of radiative transfer and statistical equilibrium equations, we obtained a modest revision of the column densities. Additionally, the column density ratios of the various deuterated isotopologues agree, within a factor of 2–3, with the previous values published for these two objects. This shows that the LTE method, which is often used to carry out a first order analysis independent of collisional rates, is a reliable alternative to analyze the emission of the NH deuterated isotopologues, at least if one is concerned with obtaining the order of magnitude of the deuterium enrichment. However, in the case a factor of 2 accuracy is required, as when constraints are to be set on the abundance ratios of molecular spin isomers, the new rate coefficients have to be used in a model where the molecular excitation is solved.
Acknowledgments
All (or most of) the computations presented in this paper were performed using the CIMENT infrastructure (https://ciment.ujfgrenoble.fr), which is supported by the RhôneAlpes region (GRANT CPER07_13 CIRA: http://www.cira.org). D.L. support for this work was provided by NASA (Herschel OT funding) through an award issued by JPL/Caltech. This work has been supported by the Agence Nationale de la Recherche (ANRHYDRIDES), contract ANR12BS05001101 and by the CNRS national program “PhysicoChimie du Milieu Interstellaire”. The authors thanks Q. Ma, P. Dadgigian and A. Van der Avoird for providing some ND cross sections. We also thank J. Harju for constructive discussions.
References
Appendix A Nd rotational levels
ND rotational levels can be inferred from ammonia spectroscopy (Townes). In their ground electronic state, both NH and ND are symmetric top with a C pyramidal symmetry. The symmetric top rotational levels are labeled where stands for the rotational momentum quantum number and is its absolute projection along the symmetry axis. The rotation states are eigenstates of the symmetric top Hamiltonian, namely where is the projection quantum number of the angular momentum along any spacefixed axis (Green76). All rotational levels are degenerate in .
As for ammonia, the nitrogen atom can tunnel back and forth through the deuterium plane. In this umbrella vibrational mode () the potential presents a double minimum along the inversion coordinate. This results in the splitting of all rotational levels into two inversion levels (states) respectively symmetric and antisymmetric with respect to the inversion coordinate. Although we neglect any vibrational coupling during the collision, the weak energy splitting between inversion levels is set to a constant cm (Tkac2014; Fusina1986397; Fusina1994464) independantly of the rotational level. The lower inversion level is symmetric whereas the upper level is antisymmetric with respect to the inversion coordinate.
a.1 Deuterium permutation symmetry
Because D nuclei have a nuclear spin I = 1, they obey Bose–Einstein statistics and the total internal wavefunction of ND is unchanged by permutation of two identical D. Thus the determination of all allowed internal states requires the description of the deuterium nuclearspin states together with the symmetry under permutation of the three identical D atoms for both motional and nuclearspin states. From BoseEinstein statistics one derives new statistical weights of the rotationinversion levels compared to NH (Bunker; Hugo2009).
In the ground electronic state, the molecular symmetry group of ND is C and the complete nuclear spin permutation group is S (Bunker). The associated irreducible representations correspond to A, A and E symmetry species. The one dimension symmetry adapted basis of A and A symmetry are wavefunctions respectively symmetric and antisymmetric under permutation of two nuclei whereas in any two dimension basis that spans the E representation, the permutation matrix is non diagonal. In the following we recall the symmetry species of each rotational, inversion and nuclearspin state.
a.2 Permutation symmetry adapted motional states
As for ammonia, the symmetric top rotational states transform as follows under permutation of two deuterons (Townes; Bunker):
Thus symmetry adapted rotational states are combinations of and states (Townes). Symmetric and antisymmetric combinations of and , for which is a multiple of , span the A and A representations :
For each rotational level the permutation symmetry species A or A are given by the sign of . Rotational states of both symmetries exist with the exception of for which and even states span the A symmetry whereas odd states span the A symmetry. Any basis ( ; ) of rotational states for which is not a multiple of spans a representation of E symmetry.
For the inversion motion the symmetric and antisymmetric vibrational states and span respectively the A and symmetry (Townes):
a.3 Symmetry adapted nuclearspin states
The total nuclearspin wavefunction is obtained from the coupling of three deuterium nuclearspin momentum with nuclearspin quantum number . The direct product of the three deuteron nuclearspin states (with ) where is the nuclearspin projection quantum number, gives deuterium spin states for ND. In terms of nuclearspin angular momentum, from the single deuteron rotation group irreducible representation, D one builds the complete nuclearspin representation as the direct product: which decomposes as follows on the total nuclearspin irreducible representations D of total quantum number I (Atkins):
Each of the angular momentum basis characterized by spans the permutation symmetry species as follows (Hugo2009; machin2006): , , and . Thus the angular momentum representation characterized by the total nuclearspin generate the permutation symmetry representations (Hugo2009):
(1)  
(2)  
(3)  
(4) 
In terms of permutation group, the total nuclearspin representation decomposes on the irreducible representations as follows:
For the total nuclear spin basis (of 27 nuclearspin states), the symmetry degeneracies are respectively . We note that the E representation is of dimension 2 and therefore the number of nuclearspin basis states of E symmetry is 16. From this decomposition we identify the three ND nuclearspin modifications as defined by Maue, in decreasing order of the symmetry degeneracy : ortho (E, ), meta (A, ) and para (A, ).
a.4 Nuclearspin statistics
The total nuclearspin rotation inversion wavefunctions are obtained from the direct product of nuclearspin states with the rotation states and vibration state. The direct product of two S irreducible representations generates a new representation shown in table 5 (Bunker; Atkins).
A  A  E  

A  A  A  E 
A  A  A  E 
E  E  E  A+A+E 
Because of BoseEinstein principle the complete wavefunction is unchanged by the interchange of two deuterium. In terms of permutation symmetry it belongs to the A symmetry of S. Therefore given the symmetry species A and A of ND inversion states and , each inversion state can only combine to nuclearspinrotation states of the same symmetry. The symmetry species of the nuclearspinrotation states can also be deduced from the direct product of Table 5. The para and meta nuclearspin state of symmetry A or A can combine to any A symmetry rotational state with , whereas the ortho nuclearspin states of symmetry E can combine to any E symmetry rotational state with . The symmetry of the resultant direct product nuclearspinrotation state is either A or A.
In reference to previous description of ammonia (Green76; Offer89; Rist93; ma2015), for each modification, the symmetry adapted ND nuclearspinrotationinversion states are commonly labeled by where specifies the symmetry of the nuclearspinrotation state: for the meta and ortho modifications and for the para modification. Therefore each nuclearspin rotation inversion state is uniquely defined by the label . Conversely, for a given spin isomer, it is only necessary to specify either or the inversion state since these quantities are related through the relations : for p–ND and for o–ND or m–ND. We note that this definition of is consistent with the previous description of NH spinrotation states but leads to an opposite inversion symmetry for the meta and ortho modifications (ma2015).
For rotational states with , as the spinrotation states are restricted to only one of the inversion state fulfills BoseEinstein principle; for the A meta spin states and for the A para spin states. For all other rotational states with , both inversion states can combine to the rotationnuclearspin states. A summary of the permutation symmetry of all nuclearspinrotationinversion states with their nuclearspin weights is given in table 6. The spin statistics of all rotationinversion states is deduced in table 7. The consequent rotational energy level diagram is given by Tkac2014.
nuclear spin ()  g  

Para (A)  A  A  A  odd  1  
A  A  A  even  1  
Meta (A) 
A  A  A  even  10  
A  A  A  odd  10  
Ortho (E) 
E  A  A    8  
E  A  A    8 
spin modification  g  

even  Meta  
even  Para  
odd  Para  
odd  Meta  
all  Meta  
all  Para  
all  Meta  
all  Para  
all  Ortho  
all  Ortho  
