Formation of the NH molecule and its isotopologues through radiative association

ABSTRACT

The rate coefficients and the cross-sections for the formation of imidogen (NH) molecule (and its isotopologues: ¹⁵NH and ND) through radiative association are determined by employing quantum mechanical perturbation theory, classical Larmor formula, and Breit–Wigner theory. We suggest the radiative association process as possible route for NH production in diffuse interstellar clouds.

molecular processes, scattering, ISM: atoms, ISM: clouds, ISM: molecules

1 INTRODUCTION

The imidogen (NH) molecule is one of the most simple and basic building block of more complex nitrogen bearing molecules. In interstellar diffuse clouds it has been detected towards ζ Persei and HD 27778 (Meyer & Roth 1991), Sagittarius B2 (Cernicharo, Goicoechea & Caux 2000), ζ Ophiuchi (Crawford & Williams 1997), in the ultraluminous infrared galaxy Arp 220 (González-Alfonso et al. 2004), and towards the star-forming regions W49N and G10.6−0.4 (Persson et al. 2012). The spectral lines of NH frequently can be observed in the spectra of comets (Swings, Elvey & Babcock 1941; Litvak & Kuiper 1982) and cool stars (Lambert & Beer 1972; Lambert et al. 1984; Smith & Lambert 1986; Aoki & Tsuji 1997). Its emission lines are used to determine the abundance of nitrogen in the interstellar space (Grevesse et al. 1990).

NH might form in a multistep mechanism in a sequence of photodissociation and reactive collisions of ions and atoms (Prasad & Huntress 1980). This established gas phase mechanism contains H₂, NH₂, NH₃, and CN molecules and the corresponding ions and atoms (see fig 1 from Wagenblast et al. 1993). The main formation paths of NH may be the recombinations of NH|$_{2}^{+}$| and NH|$_{3}^{+}$| ions with electrons. Furthermore, the associative detachment channel N + H⁻ → NH + e⁻ or in the presence of strong radiation fields the photodissociation of NH₂ might have significant contribution to the gross NH production in the diffuse clouds (Prasad & Huntress 1980). According to Wagenblast et al. (1993) this gas phase chemistry model is not able to reproduce the abundance of the observed species, like NH and CN in clouds toward ζ Persei and HD 27778. They established a new non-equilibrium model that contains grain surface production channels for NH and OH. This model could reproduce the NH and also the CN abundances in the mentioned clouds. In a recent study by Weselak et al. (2009), the spatial correlations between column densities of diatomic species have been examined. According to their results the abundance of NH molecules more significantly correlates with the neutral species than the ionic.

The radiative association mechanism might have a contribution to NH formation in the interstellar medium, but its possible role has not been discussed in the literature yet. The rate constant for the radiative association of NH molecule is still required, even though the rate coefficient for several diatomic molecules with this mechanism was calculated in the last decades with high-level quantum dynamical methods (Barinovs & van Hemert 2005; Antipov et al. 2009; Franz, Gustafsson & Nyman 2014; Gustafsson & Nyman 2014; Gustafsson, Monge-Palacios & Nyman 2014; Nyman, Gustafsson & Antipov 2015; Babb & McLaughlin 2017; Forrey et al. 2018).

There is only one potential energy curve, which of the ground state (³Σ⁻, see Fig. 1), that can provide an accessible path for dipole allowed radiative stabilization in the collision of ground state atoms, when spin–orbit and non-adiabatic couplings are neglected. Under those circumstances NH may stabilize in the following direct reaction:

\begin{eqnarray*} \mathrm{ N}\left(^{4}\mathrm{ S}\right) + \mathrm{ H}\left(^{2}\mathrm{ S}\right) \rightarrow \left[\mathrm{ NH}\left(^{3}\Sigma ^{-}\right) \right]^{*} \rightarrow \mathrm{ NH}\left(^{3}\Sigma ^{-}\right) + \hbar \omega. \end{eqnarray*}

(1)

Figure 1.

Potential energy curves of NH molecule that correlate with the ground and the two lowest excited states of N atom. The term symbols of the six lowest electronic states are assigned. Singlet manifold is represented with solid lines, triplet manifold with full circles, and the quintet with full triangles. (See Computational Details section for the further details of the employed quantum chemical level.)

Figure 2.

Potential energy and electric dipole curve for reaction (1) used in the dynamical calculations.

Our purpose in this work is to calculate the cross-sections and rate constants of formation for NH molecule (and its isotopologues: ¹⁵NH and ND) in reaction (1) in absence of spin–orbit and non-adiabatic couplings.

2 THEORY

The classical cross-section of radiative association can be obtained by the collision theory:

\begin{eqnarray*} \sigma \left(E_{\mathrm{ coll}}\right)=2 {\pi} \int_{0}^{\infty }b P_{\mathrm{ rad}}\left(E_{\mathrm{ coll}},b\right)\mathrm{ d}b , \end{eqnarray*}

(2)

where the opacity function (the probability of the radiation as a function of impact parameter, b, and collision energy, E_coll) can be computed through a combination of the classical trajectory method and the Larmor formula, where the emission is a result of a time-dependent dipole. According to this recently developed method by Gustafsson (2013) the probability of radiation for the collision of two structureless neutral particles is

\begin{eqnarray*} P_{\mathrm{ rad}} = \frac{2 f_{\mathrm{ stat}}}{3c^{3} {\pi} \hbar \left(4 {\pi} \epsilon _{0}\right)} \int_{\omega _{\rm coll}}^{\omega _{\rm max}}\omega ^{3}\left| \int_{-\infty }^{\infty }\mathrm{ e}^{\mathrm{ i}\omega t} \boldsymbol {D}\left(b,t,E_{\rm coll}\right)\mathrm{ d}t\right|^{2}\mathrm{ d}\omega, \end{eqnarray*}

(3)

where |$\boldsymbol{D}(t)$| is the time-dependent dipole vector along the trajectory, ω is the frequency of radiation, ω_coll = E_coll/ℏ and ω_max = E_max/ℏ, where E_max is maximally possible energy loss of the formed molecule corresponding to the total energy of the system and the maximum depth of the centrifugal energy corrected potential well. f_stat is the statistical weight factor that can be calculated from the spin- and orbital angular momentum multiplicity of the electronic quantum states of the reactants and of the system formed by interacting fragments (Nyman et al. 2015).

The quantum mechanical perturbation theory treatment of the radiative stabilization yields a golden rule-like formula for the cross-section (Babb & Dalgarno 1995):

\begin{eqnarray*} \sigma {\left(E_{\rm coll}\right)} = \frac{64 {\pi} ^{5} f_{\rm stat}}{3k_{\rm ini}^{2} \left(4 {\pi} \epsilon_{0}\right) } \sum _{J,j^{\prime },v^{\prime }} \frac{ S_{J,j^{\prime }}}{\lambda ^{3}_{J,j^{\prime },v^{\prime },E_{\rm coll}}} {\left| M_{J,j^{\prime },v^{\prime },E_{\rm coll}} \right|^{2}}, \end{eqnarray*}

(4)

where the sum runs over the initial angular momentum (J) and final rotational (⁠|$j^{\prime}$|⁠) and vibrational (⁠|$v^{\prime}$|⁠) quantum numbers. |$k_{\rm ini}=\sqrt{\left(2 \mu E_{\rm coll} \right)}/\hbar$| is the wavenumber, where μ is the reduced mass of the colliding fragments, |$S_{J,j^{\prime }}$| is the Hönl–London factor, and |$\lambda _{J,j^{\prime },v^{\prime },E_{\rm coll}}$| is the wavelength of the emitted photon. The transition dipole matrix elements are defined as

\begin{eqnarray*} M_{J,j^{\prime },v^{\prime },E_{\rm coll}} = \left\langle \varphi _{j^{\prime },v^{\prime }}^{\rm fin}\left(r\right)\left|\hat{\boldsymbol{D}}\left(r\right)\right|\chi _{J,E_{\rm coll}}^{\rm ini}\left(r\right)\right\rangle, \end{eqnarray*}

(5)

where |$\hat{\boldsymbol {D}}(r)$| is the operator of the dipole moment vector, |$\chi _{J,E_{\rm coll}}^{\rm ini}(r)$| is the radial part of the energy normalized scattering wavefunction of the initial state, and |$\varphi _{j^{\prime },v^{\prime }}^{\rm fin}(r)$| is the radial part of the final rovibrational wavefunction, normalized to unity.

The thermal rate constant can be calculated from the cross-section according to collision theory as

\begin{eqnarray*} k\left(T\right)=\left(\frac{8}{ {\pi} \mu k_{\rm B}^{3}T^{3} } \right)^{1/2} \int_{0}^{\infty } \sigma \left(E_{\rm coll} \right) E_{\rm coll}\, \mathrm{ e}^{-\frac{E_{\rm coll}}{k_{\rm B}T}} \mathrm{ d}E_{\rm coll}, \end{eqnarray*}

(6)

where T is the temperature and k_B is the Boltzmann constant.

The complex resonance structures can make impractical the evaluation of the rate coefficient by numerical integration in equation 2. Furthermore, the established perturbation theory approach overestimates the heights of the resonances, for which the tunnelling width is smaller than or comparable to the radiative width (Bennett el al. 2003). A possible approach to the problem is the application of the Breit–Wigner theory (Breit & Wigner 1936; Bain & Bardsley 1972). The total rate coefficient can be decomposed as

\begin{eqnarray*} k\left(T\right)=k_{\rm dir}\left(T\right) + k_{\rm res}\left(T\right), \end{eqnarray*}

(7)

where k_res(T) is the resonance contribution and k_dir(T) is the direct contribution, which can be identified as the baseline of the cross-section. This partitioning of the rate constant means that the possible interference effects between the direct and resonance mechanisms are neglected. According to the Breit–Wigner theory the resonance contribution to the rate coefficient can be obtained as (Nyman et al. 2015)

\begin{eqnarray*} k_{\rm res}\left(T\right)=\hbar ^{2}f_{\rm stat} \left(\frac{2 {\pi} }{ \mu k_{\rm B} T } \right)^{3/2} \sum _{j,v} \frac{\left(2J+1 \right) \mathrm{ e}^{-\frac{E_{j,v}}{k_{\rm B}T}}}{1/\Gamma ^{\rm tun}_{v,j}+ 1/\Gamma ^{\rm rad}_{v,j}}, \end{eqnarray*}

(8)

where |$v$| and j are vibrational and rotational quantum numbers of the quasi-bound states. |$E_{j, v}$| is the energy of the quasi-bound state that corresponds to the peak position of resonances. |$\Gamma ^{\rm tun}_{v,j}$| is the width associated with tunnelling and |$\Gamma ^{\rm rad}_{v,j}$| is the width due to the radiative decay to all lower lying levels (bound and quasi-bound) of the final electronic state.

3 COMPUTATIONAL DETAILS

The molpro package (Werner et al. 2015) was used for the quantum chemical calculations. The potential energy curve used in the dynamical calculations and the permanent dipole moment is calculated with the explicit correlated internally contracted multireference configuration interaction method (icMRCI-F12) with Davidson correction using the aug-cc-pVQZ-F12 basis set as implemented in molpro. All calculations were carried out in the C_2|$v$| symmetry group. The molecular orbitals were constructed using the state averaged complete active space self-consistent field (CASSCF) method with an active space consisting of eight electrons on 10 orbitals (5a₁, 2b₁, 2b₂, 1a₂). Three states with A₂ irreducible representation were included in the state average. Excited state potential energy surfaces are also calculated to present the manifold of the singlet, triplet, and quintet states that correlate with the ground state and the two lowest excited states of N atom (see Fig. 1). For these calculations we used the icMRCI-F12 method with aug-cc-pVTZ-F12 basis set. An active space consisting of six electrons on eight orbitals (4a₁, 2b₁, 2b₂, 0a₂) is used for the reference CASSCF wavefunction.

We have also taken the dipole function from a previous electronic structure study made by Brook et al. (2014). A discrepancy can be observed between our and the previous dipole below 0.9 Å . We have tested both data sets in perturbation theory calculations of the cross-sections in the entire energy range. The cross-sections obtained with our dipole have overestimated the other one with 10–15 per–cent. The dipole obtained by Brook et al. (2014) seems more accurate below 0.9 Å. Thus, we have used their data in the final calculations (see Fig. 2).

We have carried out the dynamical calculations of radiative association cross-sections according to the quantum mechanical perturbation theory (PT), classical (CL), and Breit–Wigner (BW) theories.

The CL mechanical simulations were done as described in Gustafsson (2013). In the PT treatment the vibrational Hamiltonian is represented on the Sinc-DVR basis in the calculations of bound states. Maximum value of the total angular momentum was set to J = 50 and the investigated energy range was E_coll = [10⁻⁵, 7.0] eV. The Numerov method is used for the calculation of the scattering wavefunctions. The energy resolution of ΔE = 10⁻⁶ eV is employed in the entire energy range for the excitation function of radiative association to properly represent the narrow resonances. The parameters needed in the BW formulas are computed with the program level (Le Roy 2007). The value of the statistical weight factor is f_stat = 3/8 for reaction (1).

4 RESULTS AND DISCUSSION

The cross-section for formation of NH and its isotopologues, obtained with the CL and quantum mechanical PT, is shown in Fig. 3. The CL and the PT curves show good agreement in the whole energy range: The CL one can be considered as the baseline of the PT curve. Similar behaviour has been observed in the previous studies of diatomic molecules (Nyman et al. 2015). Because of this agreement the pure quantum effects can be isolated. The quantum mechanical curves have a rich resonance structure and in the cases of NH and ¹⁵NH a tunnelling effect can be observed below 10⁻³ eV. The isotope substitution also influences the dynamics. As expected from the reduced masses, the deuterium exchange has bigger effect than the ¹⁵N substitution. The CL curve of ¹⁵NH is almost the same as the original NH curve, but the cross-section of ND is smaller. The effect of isotope substitution is more significant in the PT results: the position and magnitude of resonances strongly depend on the masses of the colliding partners.

Figure 3.

CL and quantum mechanical PT cross-sections for formation of NH and its isotopologues.

Because of the narrow resonances at high energies we should use a small energy step and account for partial waves with high J values. Under those circumstances the PT method can be cumbersome even for a diatomic system. As mentioned in Section 2, the PT method overestimates the heights of the resonances. Moreover, the transitions from the quasi-bound to lower-lying quasi-bound states are not included in PT cross-sections. That is why the cross-sections and rate constants are also calculated with the BW theory to provide a complete description of the resonances in the whole energy range. The CL cross-sections are used as the direct contribution to supplement the BW theory. Fig. 4 shows the comparison of the BW and the PT cross-sections for the NH formation. The BW result almost perfectly matches the position of resonances that are obtained with the PT method. At high energies (above 0.5 eV) the difference between the PT and BW+CL cross-sections is more significant: there are a lot of very narrow BW resonances that are not predicted by PT method. These resonances correspond to the quasi-bound to quasi-bound transitions. Although these transitions result in metastable states, they may have a significant contribution to the stable molecule formation (Bennett et al. 2003). That is because of the decay into the continuum can be slower than the radiative relaxation to the stable rovibrational states.

Figure 4.

Radiative association cross-sections for NH molecule obtained by Breit–Wigner theory and quantum mechanical perturbation theory.

We tested the validity of recently established isotope-scaling relation for the CL cross-sections (Gustafsson 2013). The isotope substitution significantly can change the dynamics of the collisions. Although the dipole and the potential energy surface are independent of isotopes in the Born–Oppenheimer approximation, but the time-scale of a collision process is strongly depends on masses. Based on the CL analysis of mass-dependent time-scales in collisions of two atoms, the following isotope-scaling relation can be obtained for the cross-sections (Gustafsson 2013):

\begin{eqnarray*} \sigma _{1}\left(E_{\rm coll}\right) \approx \frac{\mu _{2}}{\mu _{1}}\sigma _{2}\left(\sqrt{\frac{\mu _{1}}{\mu _{2}}}E_{\rm coll}\right), \end{eqnarray*}

(9)

where σ₁ and σ₂ are the cross-sections and μ₁ and μ₂ are the reduced masses of the isotopologues. Fig. 5 shows the comparison of the cross-sections for the isotopologues that are calculated with the CL method and obtained from the scaling of the original NH curve according to formula (equation 9). In both cases the isotope-scaled curves excellently match the dynamical results. This means that the isotope-scaling relation is fulfilled for the direct contribution of this reaction.

Figure 5.

Comparison of the CL cross-sections and those obtained by the isotope-scaling relation.

The rate constants for the formation of NH molecule and its isotopologues are presented in Fig. 6. As expected from the cross-sections, the CL rate constants of NH and ¹⁵NH are almost the same in the whole energy range and the rate constant of the ND formation is smaller. Because of the resonance contribution the PT rates coefficients are much larger than the CL ones, mainly at low temperatures. For the ND and ¹⁵NH molecules we simply used the CL and BW rate constants for the direct and resonance contribution, respectively, but for the NH molecules we used a different approach. Only in the case of NH molecule there is no resonance peak below 10⁻³ eV and the baseline of the PT is significantly larger than the CL curve. It means that only the tunnelling into the repulsive wall of the potential is responsible for the increased reactivity in the quantum mechanical calculations. To consider this tunnelling effect in the case of NH molecule we combined the CL and PT curve to calculate the direct contribution to the BW theory. The PT cross-section was used at low energies until it first crossed the CL one (at E_coll = 7 × 10⁻⁴ eV). Above that energy the CL curve was used in the rate coefficient calculation for the BW+CL results. Although in the case of ¹⁵NH there is a big tunnelling contribution at low energies, we could not use a similar combined cross-section for the direct contribution in BW theory, because there is a resonance below 10⁻³ eV.

Figure 6.

Rate constants as a function of temperature for the formation of NH and its isotopologues.

Following the recommendation of Antipov et al. (2009), the BW+CL cross-section is used to compute the rate constant in this work. The extended Arrhenius equation

\begin{eqnarray*} k\left(T\right) = A\left(\frac{T}{300} \right) ^{\alpha } \mathrm{ e}^{-\beta / T} \end{eqnarray*}

(10)

is used to fit the BW+CL rate constants. The temperature dependence of the rate coefficients shows a complex structure with multiple inflexion points. This complex structure is caused by the quantum mechanical effects like the resonances and the tunnelling at low temperatures. Because of this structure we split the curves into seven or eight temperature ranges and fit them individually to get applicable fitting parameters in a wide temperature range (2–10 000 K). The fitting parameters are summarized for the formation of NH and its isotopologues in Table 1. These rate constants will be made available in the data base KIDA (Wakelam et al. 2012).

Table 1.

Fitting parameters of the extended Arrhenius curves for the formation of NH and its isotopologues.

	T (K)	A (cm3 s−1)/10−21	α	β (K)
NH	2–3	352 413.8	2.7452	−4.3864
NH	3–10	0.6077	−0.9286	6.5983
NH	10–25	1.9275	−0.3847	−0.3896
NH	25–70	13.7884	0.8420	−27.5819
NH	70–500	8.5618	0.1772	6.5126
NH	500–1000	13.9902	−0.1602	167.9334
NH	1000–2500	44.6853	−0.6699	721.3719
NH	2500–10 000	292.6699	−1.2685	2273.757
ND	2–10	1.3541	−0.1361	−2.0892
ND	10–25	0.9766	−0.2607	−1.3074
ND	25–100	2.7236	0.4116	−17.6110
ND	100–600	2.8390	0.0960	22.9811
ND	600–1000	5.1260	−0.2797	225.4916
ND	1000–2500	13.8915	−0.7144	703.7595
ND	2500–10 000	75.8599	−1.2525	2119.963
15NH	2–10	0.1725	−1.3963	3.9834
15NH	10–25	1.9086	−0.3950	−6.1760
15NH	25–60	9.4229	0.5593	−25.8539
15NH	60–300	7.6309	0.3402	−17.35167
15NH	300–1000	11.6185	−0.0696	111.6673
15NH	1000–2500	43.5454	−0.6648	725.3351
15NH	2500–10 000	285.5763	−1.2630	2281.373

The magnitude of the rate coefficients is similar to formation of other diatomic molecules (Nyman et al. 2015), but it is much smaller than rates of exchange reactions. Despite this small rate constant this process might be important, because it is a direct route to the NH formation, rather than through a complex chemical network (see fig. 1 from Wagenblast et al. 1993) where multiple consecutive exchange reactions with polyatomic reactants are needed for the NH production.

5 CONCLUSIONS

In this work the formation of NH molecules and its isotopologues by radiative association mechanism is studied. The cross-sections and thermal rate coefficients for the radiative stabilization are obtained by quantum mechanical PT, BW theory, and by a CL method. Some distinct isotope effects are observed. For example, the rate constant for the formation of ¹⁵NH is at least five time larger than that of NH at a few Kelvins. We suggest the radiative association process as possible route for NH formation in the dust-poor interstellar medium. The numerical data set of the presented results is attached as supplementary material.

SUPPORTING INFORMATION

NH_paper_support_info.tar.gz

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

Support from the Kempe Foundation is gratefully acknowledged. MG acknowledges support from the Knut and Alice Wallenberg Foundation. Computational resources provided by the Swedish National Infrastructure for Computing (SNIC) at HPC2N are acknowledged.