Full-dimensional quantum dynamics of rovibrationally inelastic scattering between CN and H₂

The hydrogen molecule is the dominant species in the vast majority of interstellar environments, while the cyano radical (CN) is also relatively abundant. Both may play important roles in the physics and chemistry of diffuse clouds, for example. CN is mainly observed via visible absorption in diffuse and translucent clouds,^1,2 and it is also detected in the millimeter where it can be important in the chemistry of dense molecular clouds. Ever since it was identified by McKellar in 1940,³ CN has been detected in a variety of molecular clouds.^4–8 For example, Wootten et al.⁴ reported the first observations of the $N=2−1$ line at 226.9 GHz in the circumstellar envelope of the carbon-rich star IRC + 10216 and in the Orion Molecular Cloud (OMC-1). CN is also believed to be an excellent tracer of photodissociation regions (PDRs) which are environments exposed to intense ultraviolet (UV) radiation.^9,10 To model the infrared radiation from CN in such environments requires a rich variety of molecular data, including inelastic rotational and vibrational collision rate coefficients, primarily due to H₂. However, because of the difficulty of measurements of these quantities, numerical models often rely on theoretical rates. While state-to-state inelastic rotational (de)excitation rate coefficients due to H₂ collisions were measured for the vibrationally excited state v₁ = 2 by Brunet et al.,¹¹ no data are available for vibrationally inelastic transitions in CN induced by H₂.

Due to its fundamental and astrophysical importance, the CN–H₂ collisional system has been the subject of many theoretical^12–18 and experimental^11,19 studies. An early four-dimensional (4D) potential energy surface (PES) within the rigid-rotor model was obtained by Kaledin et al.¹² and applied to a study of the spectroscopy of the CN–H₂ van der Waals complex. More recently, Kalugina et al.^14,15 calculated a 4D PES using the partially spin-restricted coupled-cluster with single, double, and perturbative triple excitations (RCCSD(T)) method and an aug-cc-pVTZ basis set augmented by 3s, 2p, 1d bond functions. Fine- and hyperfine-resolved collisional excitation calculations of CN (v₁ = 0) with para-H₂ (v₂ = 0,j₂ = 0) were performed, but using a two-dimensional (2D) PES obtained by averaging the 4D PES over the H₂ angular degrees of freedom. Later, the 4D PES of Kalugina et al.^14,15 was used in close-coupling (CC) calculations of rotational de-excitation cross section of CN (v₁ = 0) by para- and ortho-H₂.¹⁵ In both PES and scattering calculations,^14,15 fixed bond distances (i.e., equilibrium bond distance of CN at r₁ = 2.2144 a₀ and vibrationally averaged bond distance of H₂ at r₂ = 1.4487 a₀) were adopted. Recently, the 4D PES of Kalugina and Lique was also applied to calculate fine- and hyperfine-resolved rotational excitation rate coefficients of CN with para- and ortho-H₂¹⁷ and of CN isotopologues, ¹³CN and C¹⁵N, by para-H₂.¹⁸ 

To go beyond the rigid-rotor approximation, in this work we have constructed a full dimensional (6D) PES for the CN–H₂ system and carried out comprehensive CC calculations of state-to-state cross sections in an effort to initiate the development of a complete database of rovibrational quenching rate coefficients. In Secs. II–V, the 6D PES calculations, the adopted scattering approach, and the results are discussed. Conclusions are presented in Sec. VI. Transitions at the fine- or hyperfine-resolved level are not considered here.

The electronic ground state interaction potential of CN–H₂ was computed using the explicitly correlated coupled-cluster (RCCSD(T)-F12B) method,^20,21 with an aug-cc-pVTZ basis,²² as implemented in MOLPRO 2010.1.²³ The interaction energies were calculated using the supermolecular approach in which the counter-poise (CP)²⁴ correction was employed to reduce the basis set superposition error (BSSE).²⁵ To construct the PES, the computations were performed on a 6D grid using Jacobi coordinates as shown in Fig. 1. The radial coordinate R denotes the distance between the center-of-masses of CN and H₂, while r₁ and r₂ represent the bond lengths of CN and H₂, respectively. The angular coordinates $θ1$ and $θ2$ de angle between $r→1$ and $R→$ and $r→2$ and $R→$⁠, while ϕ denotes the out-of-plane dihedral angle. For the potential energy calculations, the bond lengths of H₂ and CN spanned the range 1.1011–1.8011 a₀ and 2.0184–2.5184 a₀, respectively. R is scanned from 5 to 21 a₀ with a step-size of 0.5 a₀ for $R<8$ a₀ and 1.0 a₀ for $R>8$ a₀. The angle $θ1$ ranges from 0° to 180° with a step-size of 22.5°, and $θ2$ and ϕ have range 0° to 180° and 0° to 90°, respectively, with a step-size of 30°. The range of r₁ and r₂ allows for the treatment of vibrational excitation up to v₁ = 2 for CN, but only for v₂ = 0 for H₂, and as the bond lengths do not approach dissociation, rearrangement channels (e.g., to form HCN) are excluded from consideration on this PES. While the barrier to reaction is around 1260 cm⁻¹ (3.6 kcal/mol), which is comparable to the CN (v₁ = 1) excitation energy, CN plays the role of a spectator in the reaction dynamics because the effective barrier at the saddle point increases by the same energy as the CN (v₁ = 1) excitation energy. As a consequence, reaction channels are not incorporated into the PES or scattering.

An invariant polynomial method^26,27 has been used to fit the 6D CN–H₂ interaction PES, which is expanded in the form

where y_i = e^−0.5d_i is a Morse-type variable, d_i are the internuclear distances (in atomic units) between two atoms, $d1=dHH′$⁠, $d2=dNH′$⁠, $d3=dCH′$⁠, $d4=dCH$⁠, $d5=dNH$⁠, and $d6=dCN$⁠. The total power of the polynomial, $N=∑i=17ni$⁠, was restricted to 7. The fit for the expansion coefficients $cn1…n6$ was performed using a weighted least-squares method for interaction energies up to 3000 cm⁻¹. The root-mean-square (RMS) error in the PES fit was 3.67 cm⁻¹, which included 84 734 ab initio points. The 6D interaction PES, $V(R,r1,r2,θ1,θ2,ϕ)$⁠, hereafter referred to as V6ES, can be averaged over the ground-state or excited vibrational wave functions of CN and H₂ to obtain 4D vibrationally averaged PESs, $V(R,θ1,θ2,ϕ)$⁠. However, V6ES was used in the majority of the scattering calculations. The comparison of the global and local minima on the vibrationally averaged V6ES-4D PES and the 4D PESs of Kalugina et al.¹⁵ and Kaledin et al.¹² is given in Table I. Note that this is not an exact comparison as the 4D PESs were constructed differently.

Some features of the fitted V6ES PES are illustrated in Figs. 2–4. Fig. 2 displays the R dependence of V6ES for various angles and bond distances fixed at r₁ = r_e(CN) = 2.2144 a₀ and r₂ = r_e(H₂) = 1.4011 a₀. The computed ab initio potential energy points are also displayed as symbols. Fig. 3 shows two-dimensional contour plots in $θ1$ and $θ2$ for the V6ES PES. The CN and H₂ bond lengths are fixed at their equilibrium distances for R = 3.0, 3.5, 4.0, and 6.0 Å, ϕ= 0°, 45°, and 90°. This figure can be compared to Fig. 3 of Kalugina et al.¹⁵ and Fig. 2 of Kaledin et al.¹⁴ The current Fig. 4 depicts the anisotropy of the V6ES potential in $θ1$ for variousϕ. Here $θ1$ describes the orientation of CN; the anisotropy with respect to $θ1$ is large.

In the quantum-mechanical description of inelastic molecular scattering, the theory for the collision of two ¹Σ diatomic molecules within the rigid rotor and vibrating rotor models was developed^28–31 based on the close-coupling (CC) formulation of Arthurs and Dalgarno.³² Recently the theory has been implemented for the full-dimensional study of the vibrationally inelastic scattering between two diatomic molecules.^33–35 Here, we apply this formalism for the CN + H₂ system by approximating CN as a ¹Σ molecule, hence the rotational quantum number N₁ = j₁ since we neglect spin. In the scattering calculations, the interaction PES $V(R,r1,r2,θ1,θ2,ϕ)$ was expanded in the form

where $Yλ1,λ2,λ12$ is a bi-spherical harmonic function given by

where $0≤λ1≤8$⁠, $0≤λ2≤4$⁠, and only even values of $λ2$ are retained due to the symmetry of H₂.

We use the term combined molecular state (CMS) to describe a combination of rovibrational states of the CN–H₂ scattering system. A CMS, denoted as (v₁j₁v₂j₂), represents a unique quantum state before or after a collision: v = v₁v₂ and j = j₁j₂ denote the vibrational and rotational quantum numbers of the two molecules.

The rovibrational state-to-state cross section as a function of collision energy E is given by

where (v₁j₁v₂j₂) and (⁠$v1′j1′v2′j2′$⁠) denote the initial and final CMSs, respectively, and l is the orbital angular momentum quantum number. The total angular momentum of the collision system is given by $J→=l→+j→12$ with $j→12=j→1+j→2$ and $k=2μE/ℏ2$ is the wave vector with respect to the initial channel.

For a given H₂ transition, $j2→j2′$⁠, the total vibrational quenching cross section of CN from $v1→v1′$ can be obtained by summing over the final rotational levels $j1′$ of CN in vibrational state $v1′$⁠,

In this work, H₂ is taken to remain in its ground vibrational state, $v2=v2′=0$⁠. However, the CN total vibrational quenching cross sections were calculated for both elastic (⁠$j2′=j2$⁠) and inelastic (⁠$j2′≠j2$⁠) H₂ transitions.

All CC scattering calculations were carried out using the TwoBC code³⁶ with the coupled-channel equations propagated using the log-derivative method of Johnson³⁷ and Manolopoulos³⁸ with a radial step-size of $ΔR=0.05$ a₀ and the range of R propagation from 5 to 21 a₀. The number of discrete variable representation points N_r1 and N_r2; the number of points in $θ1$ and $θ2$ for Gauss-Legendre quadrature, $Nθ1$ and $Nθ2$⁠; and the number of points in ϕ for Chebyshev quadrature, $Nϕ$⁠, used to project out the potential expansion coefficients are presented in Table II. The basis sets used in the pure rotational and rovibrational scattering calculations are also presented in Table II.

The state-to-state rate coefficients at a temperature T can be obtained by thermally averaging the corresponding integral cross sections over a Maxwellian kinetic energy distribution,

where $μ$ is the reduced mass of the CN–H₂ complex, $β=(kBT)−1$⁠, and k_B is Boltzmann’s constant.

We first consider the rotational transitions of CN in collisions with H₂ and D₂ when the molecules remain in their ground vibrational states $v1=v2=v1′=v2′=0$⁠,

and

Full-dimensional calculations using TwoBC with V6ES and rigid rotor approximation calculations using MOLSCAT³⁹ with vibrationally averaged 4D PESs have been performed to obtain the rotational excitation and de-excitation cross sections of CN in collisions with H₂. For CN with D₂, only rigid-rotor approximation calculations using vibrationally averaged 4D PESs are carried out. Three vibrationally averaged 4D PESs in the CN ground and excited vibrational states v₁ = 0, 1, and 2, referred to as V4H2-0, V4H2-1, and V4H2-2, are computed for CN–H₂. For CN–D₂, the vibrationally averaged 4D PESs are referred to as V4D2-0, V4D2-1, and V4D2-2. In constructing these 4D PESs, H₂ and D₂ are in their ground vibrational state v₂ = 0.

The initial rotational state of CN was taken to be j₁ = 4 and 11 to make comparisons with available measurements as described below. For para-H₂ and ortho-D₂, rotational levels j₂ = 0 and 2 were included in the basis set while for ortho-H₂ and para-D₂, j₂ = 1 and 3 were included.

Rotational excitation and de-excitation cross sections of CN (v₁ = 0,j₁ = 4) in collisions with para-H₂ (v₂ = 0,j₂ = 0, 2) and ortho-H₂ (v₂ = 0,j₂ = 1) have been computed for collision energies ranging from 1.0 to 1500 cm⁻¹. In Figs. 5–7, the calculated de-excitation cross sections using V6ES and V4H2-0 are compared with the 4D results of Kalugina et al.¹⁵ It can be seen that the 6D V6ES and 4D V4H2-0 PESs yield nearly identical results. The cross sections for de-excitation from j₁ = 4 are seen to be in good agreement with the results of Kalugina et al.

Brunet et al.¹¹ measured the rotational inelastic, but vibrationally elastic, state-to-state rate coefficients of CN (⁠$X2Σ+,v1=2$⁠) due to collisions with normal H₂ and D₂. For comparison, the state-to-state rotational excitation and de-excitation rate coefficients of CN (v₁ = 0,j₁ = 4 and 11) in collisions with H₂ and D₂ with v₂ = 0, j₂ = 0, 1, 2 are computed using Eq. (6). The Boltzmann average, over the H₂ and D₂ rotational levels, of the state-to-state rate coefficients k(j₁j₂) was then obtained with the partition function Q of the first three rotational states of H₂ and D₂ (j₂ = 0, 1, 2) and given by

where

k(j₁j₂) is the rate coefficients from initial CMS (0j₁0j₂), $ϵj2=B0j2(j2+1)$⁠, the rotational constant B₀ is 59.3219 cm^$−$1 for H₂ and 29.9043 cm⁻¹ for D₂, and $w$_j2 is the spin degeneracy, $w$_j2 = 3 for ortho-H₂ and para-D₂ and 1 for para-H₂ and ortho-D₂.

In Figs. 8 and 9 we compare the present rate coefficients with measurements¹¹ and the theoretical results of Kalugina et al.¹⁵ For collisions with normal H₂ and j₁ = 4, the results of full-dimensional calculations using V6ES and rigid-rotor approximation calculations using vibrationally averaged PESs V4H2-0, V4H2-1, and V4H2-2 are displayed in Fig. 8(a). The 4D PESs yield very similar scattering results, which are also comparable to the 6D results. Due to the degeneracy factors and assumed relative rotational populations given by Equation (7), the rate coefficients for normal-H₂ are dominated by ortho-H₂ transitions. Therefore, while Fig. 5 displays an even $Δj1$ propensity rule for para-H₂, applicable for a near-homonuclear molecule like CN, an exponential energy-gap-law behavior is evident for ortho-H₂ cross sections in Fig. 6. As a consequence, any evidence of an even $Δj$ propensity is washed-out in Fig. 8 with the result that $Δj1=±1$ dominates. For j₁=11, presented in Fig. 8(b), the even $Δj1$ propensity is partially recovered, though $Δj1=±1$ transitions give the largest rate coefficients.

For both j₁ = 4 and 11, the comparison of the present rate coefficients with experiment and the theoretical results of Kalugina et al.¹⁵ shows good agreement particularly in light of the fact that this is an absolute value comparison. The exception occurs for $Δj1$ = $±$1 transitions where noticeable differences are found between the experiment and all calculations. The calculations of Kalugina et al.¹⁵ are somewhat closer to the measurements which are mostly due to their slightly larger (0401) to (0301) cross section (see bottom right panel of Fig. 6).

Fig. 9 shows similar comparisons between our theoretical results of CN in collision with D₂ and measurements.¹¹ The results from three vibrationally averaged 4D PESs, V4D2-0, V4D2-1, and V4D2-2 are very similar for both j₁=4 and 11, except that V4D2-2 gives smaller rate coefficients for $Δj1=±1$ transitions.

While the measurement was performed for vibrationally excited CN (v₁ = 2), calculations which used v₁ = 0, though for vibrationally averaged PESs, are not direct comparisons. We explored whether this was an issue by repeating our 4D calculations on v₁ = 1 and 2 vibrationally averaged PESs. Figs. 8 and 9 demonstrate that rotational transitions are nearly independent of v₁, except as noted above, the $Δj1=±1$ transitions are somewhat smaller for D₂ collisions on the V4D2-2, i.e., v₁ = 2 surface. This might be related to the partial recovery of an even $Δj1$ propensity for ortho-D₂ collisions. Nevertheless, it appears that CN is sufficiently harmonic whose dependence on v₁ for pure rotational transitions is small.

Returning to the experiment, they were state-to-state measurements performed by laser probing the CN level populations and then using a master equation analysis to extract rate coefficients. For $v1>0$⁠, the level populations may also be quenched by reactive scattering to H + HCN and vibrational de-excitation, which are not considered in the experimental modeling. Vibrational quenching is typically many orders of magnitude smaller than rotational quenching but may become competitive for a quasi-resonant transition as found for CO–H₂³⁵ (e.g., the 1002→0006 + 83 cm⁻¹ transition, in CMS notation). For CN–H₂, the quasi-resonant transitions 1002 →0006 −18 cm⁻¹ and 2002→1006 − 45 cm⁻¹ are potentially relevant to the Brunet et al.¹¹ measurements. Both processes are slightly endoergic and will be studied in future work.

Full-dimensional quantum CC calculations of rovibrational cross sections were carried out for initial CMSs (100j₂), for para-H₂ (j₂ = 0) and ortho-H₂ (j₂ = 1) with collision energies ranging from 1 to 1000 cm⁻¹. In Figs. 10 and 11 the state-to-state quenching cross section to final states (0$j1′$0j₂), $j1′$ = 0, 2, 4, …, 24 are shown for para- and ortho-H₂, respectively, where $j2=j2′$⁠. The cross sections have similar behavior and display resonances at low energies due to quasibound states of the collision complex. Small $|Δj1|=|j1′−j1|$ transitions dominate the quenching with the cross sections generally decreasing with increasing $j1′$ with that for $j1′=0$ being the largest. This indicates that the cross section is dominated by the strength of the interaction terms $Aλ1,λ2,λ12$ and not by the asymptotic channel energy gaps.

The total cross sections $σv1j1v2j2→v1′v2′j2′T(E)$ for CN vibrational quenching from $v1=1→v1′=0$ are given in Fig. 12. Fig. 12(a) shows the total cross section from (1000) to CN (⁠$v1′=0$⁠) + para-H₂ (⁠$j2′$ = 0, 2), while Fig. 12(b) displays the total cross section from (1001) to CN (⁠$v1′=0$⁠) + ortho-H₂ (⁠$j2′$=1, 3). The total vibrational quenching cross section display orbiting resonances between collision energies of 1.0 and 100 cm^$−$1. However, the magnitude of the resonances are suppressed for H₂ elastic processes, $j2=j2′$⁠. For both the para- and ortho-H₂ cases, the total vibrational quenching cross sections for elastic H₂ processes are nearly an order of magnitude larger than that for collisions involving inelastic H₂ transitions.

Total rate coefficients for the vibrational quenching of CN from (1000) to CN (⁠$v1′=0$⁠) + para-H₂ (⁠$j2′$ = 0 and 2) and from (1001) to CN (⁠$v1′=0$⁠) + ortho-H₂ (⁠$j2′$ = 1 and 3) are displayed in Fig. 13 for temperatures ranging from 5 to 300 K. The trends are similar to that given for the cross sections in Fig. 12. The total vibrational quenching rate coefficients of CN for $Δj2=0$ are nearly an order of magnitude larger than the results for $Δj2=2$⁠. For the quenching due to para-H₂ collisions, as shown in Fig. 13(a), between 5 and ∼100 K, the rate coefficients generally decrease weakly with increasing temperature, while weak oscillatory temperature dependences are exhibited due to the presence of resonances. For temperatures above ∼100 K, the rate coefficients increase with increasing temperature. The trends for ortho-H₂ are very similar to those noted for para-H₂ collisions. CO–H₂ rate coefficients³⁵ for the same transitions are typically factors of ∼5-10 smaller, as shown in Fig. 13. The smaller CO–H₂ rate coefficients are likely related to the fact that the CN–H₂ potential is more anisotropic, e.g., with global and local minima deeper by 38.69 and 11.12 cm⁻¹, respectively.

As highlighted in the Introduction, in astrophysical environments with intense UV fields, the radiation can drive the chemistry and internal level populations out of equilibrium. In such situations, a photodissociation region (PDR) resides at the interface of the hot H ii region and the cold molecular region. Sternberg and Dalgarno⁴⁰ investigated the production of gas-phase atomic and molecular species in dense molecular clouds and identified molecular diagnostics of photon-dominated chemistry. It was illustrated that CN was formed by photodissociation of HCN, HCN + $ν→$CN + H, and by the reactions NO + C → CN + O and C + N → CN + ν. Minor production of CN can also be obtained from dissociative recombination H₂CN⁺ + e → CN + H₂. The CN abundance strongly depends on the UV radiation and increases in PDRs as the HCN photodissociation rate increases. Hence, CN emission lines may serve as useful diagnostics of PDRs.

Rodriguez-Franco and Martin-Pintado⁹ reported the observations of the CN 1-0 and 2-1 rotational emission towards the Orion A molecular complex and showed that the ionized gas of the H II region dominated the morphology of the CN emission. Aalto et al.¹⁰ detected CN 1-0 and 2-1 line emission in IR luminous galaxies and found that CN emission required high density, $n>104$ cm⁻³. Therefore CN is expected to be a good tracer of high density regions and cold dense molecular clouds. This was also confirmed by Hakobian and Crutcher⁴¹ in their multi-species study of structure and composition of molecular clouds through the CN Zeeman transitions. Finally, Podio et al.⁴² argue that high j₁ lines, particularly the 5-4, are unique tracers of embedded disks of young stellar objects, such as T Tau N.

We are, however, unaware of any CN vibrational transition observations though it is predicted to be abundant in models of protoplanetary disks⁴³ and in the outflows of giant C- and O-rich evolved stars.^44,45 It is anticipated that vibrational emission lines can probe the hot inner wind regions of the outflows and therefore provide information on the stellar mass-loss rates.⁴⁴ With the upcoming launch of NASA’s James Webb Space Telescope, the CN vibrational bands will become observable allowing astronomers to see deep into warm regions of these objects, particularly where the vibrational lines can be pumped by intense ultraviolet radiation. The collisional data computed here will be critical for models of such near infrared observations.

Accurate collisional rate coefficients for CN due to H₂, needed for astrophysical modeling, are currently lacking for vibrationally excited CN. Quantum scattering calculations are the primary source of such rate coefficients, but accurate values that explicitly consider all six internal degrees of freedom of the CN–H₂ complex have not been reported due to the difficulty of the computations. We performed the first full-dimensional quantum dynamics calculations for inelastic rovibrational transitions of CN in collisions with H₂ using a high-level 6D PES. The PES was fitted with an invariant polynomial approach with an RMS error of less than 3.8 cm⁻¹. Cross sections and rate coefficients for rotational (de)excitation of CN due to H₂ and D₂ collisions were studied using the new 6D PES for excited initial rotational levels of CN (j₁ = 4 and 11). Good agreement with the results of previous theoretical and experimental studies confirms the accuracy of the current PES fit and scattering calculations. We have carried out the first calculations of rovibrational quenching of CN (v₁=1,j₁ = 0) in collisions with H₂ (v₂ = 0,j₂ = 0, 1, 2) on the new 6D PES. To aid in astrophysical modeling, computations are in progress for a more highly excited vibrational level of v₁ = 2 of CN.

See supplementary material⁴⁶ for a Fortran subroutine to generate V6ES PES.

Work at UGA and Emory was supported by NASA Grant No. NNX12AF42G from the Astronomy and Physics Research and Analysis Program, at UNLV by NSF Grant No. PHY-1505557, and at Penn State by NSF Grant No. PHY-1503615. This study was supported in part by resources and technical expertise from the UGA Georgia Advanced Computing Resource Center (GACRC), a partnership between the UGA Office of the Vice President for Research and Office of the Vice President for Information Technology. We thank Shan-Ho Tsai (GACRC) and Jeff Deroshia (UGA Department of Physics and Astronomy) for computational assistance. We also thank Dr. Yulia Kalugina for providing their potential subroutine and cross-sectional data.