QUANTUM CALCULATION OF INELASTIC CO COLLISIONS WITH H. I. ROTATIONAL QUENCHING OF LOW-LYING ROTATIONAL LEVELS

The current understanding of the chemical and physical conditions in diffuse clouds (Herbst 2001), and for the majority of other low-temperature astrophysical environments, is deduced primarily from observations of molecular hydrogen and CO, which are typically the most abundant molecules. CO can be detected in absorption in the ultraviolet (UV) and near infrared (IR) and in emission at wavelengths from the far IR (FIR) to the submillimeter (submm). As examples, CO IR emission lines, due to rovibrational transitions, toward Orion Peak 1 and 2 were observed by the Infrared Space Observatory (ISO) with the Short Wavelength Spectrometer (González-Alfonzo et al. 2002). The Long Wavelength Spectrometer on ISO detected CO pure rotational lines in the C-rich planetary nebulae NGC 7027 (Liu et al. 1996; Cernicharo et al. 1997), in the C-rich objects AFGL 2688 and AFGL 618 (Justtanont et al. 2000), and in the reflection nebula NGC 1333 (Molinari et al. 2000). Using Herschel/HIFI, Bujarrabal et al. (2010) observed the protoplanetary nebula CRL 618, detecting high-j transitions of CO in the FIR/submm. Herschel/HIFI detection of high-j lines was also reported by Yildiz et al. (2010) for CO rotational levels up to j =10 toward the stars NGC 1333 IRAS 2A, IRAS 4A, and IRAS 4B. Recently, Castro-Carrizo et al. (2012) carried out subarcsecond-resolution observations of CO line emission from the early planetary nebula M 2-9, while Lupu et al. (2012) performed CO redshift measurements in the H-ATLAS survey. A CO excitation analysis, considering departures from local thermodynamic equilibrium (LTE) using the RADEX modeling package (van der Tak et al. 2007), was also carried out by Lupu et al. (2012) for the 4 → 3 and 10 → 9 pure rotational lines. Carbon monoxide was also detected in the brown dwarf Gliese 229B, but the predicted CO abundance obtained from LTE spectral modeling was ∼10³ times larger than the value computed using local chemical equilibrium models (Noll et al. 1997).

In the majority of low-density environments, the level populations of molecules typically depart from LTE. As a consequence, to accurately predict spectral line intensities, rate coefficients for collisional excitation due to the impact of the dominant species, H₂, H, and He are necessary. For CO rotational excitation in collision with H atoms, one of the earliest quantum calculations was carried out by Chu & Dalgarno (1975) using semiempirical potentials. They presented rate coefficients for 5–150 K. Green & Thaddeus (1976) performed another early quantal scattering calculation for the CO–H system over the same temperature range, but using different semiempirical potentials. From their calculations, it was found that the collisional rate coefficients for CO with H were significantly smaller than those of CO by H₂, so that the contribution of CO–H collisions was usually neglected in later CO emission modeling.

Later theoretical investigations on the CO–H collision system have been motivated by the availability of ab initio potential energy surfaces (PESs) with increasing sophistication. An early ab initio surface was obtained by Bowman et al. (1986), the so-called BBH PES. BBH has a barrier, with a calculated height of 0.087 eV (Wang et al. 1973), to the formation of the radical HCO. It was demonstrated from a later scattering calculation by Lee & Bowman (1987), using the BBH surface, that for collision energies below the barrier, a homonuclear-like behavior was observed for the state-to-state rotational excitation cross sections, i.e., the cross sections displayed an even-Δj propensity.

Keller et al. (1996) constructed another ab initio PES for the CO–H system using the multireference configuration-interaction (MRCI) method. They reported that this PES, called the WKS surface, gave better agreement with experimental spectroscopy than the BBH surface. Using both the BBH and WKS potentials, Green et al. (1996) performed scattering calculations for CO rotational excitation. Cross sections on the two potentials revealed significant discrepancies, particularly for the j = 0 → 1 transition. These calculations were revisited by Balakrishnan et al. (2002, hereafter BYD), who verified the differences in the cross sections obtained in the BBH and WKS surfaces observed by Green et al. (1996). BYD, however, used the WKS PES to compute a large set of CO–H rotational excitation data. Subsequently, the BYD rate coefficients have been used in a variety of astrophysical models. For example, Liszt (2006), adopted the CO–H rotational excitation rate coefficients of BYD to investigate CO rotational populations in the diffuse interstellar medium (ISM). Based on comparisons to observations, Liszt (2006) argued that H impact excitation makes an important contribution to the CO rotational populations and should be included in models. As argued by Lee & Bowman (1987), an even-Δj propensity is likely to be valid for CO. However, it was later pointed out by Shepler et al. (2007) that the excitation rate coefficients of BYD from initial j = 0 failed to obey the homonuclear-like propensity for even-Δj transitions, because the computed j = 0 → 1 rate coefficients are larger than the j = 0 → 2 rate coefficients.

In order to investigate this apparent departure from the expected behavior, two completely independent CO–H rigid-rotor PESs were constructed in Shepler et al. (2007) using state-of-the-art methodologies in quantum chemistry. One of the PESs was computed using the coupled-cluster singles-and-doubles-with-perturbative-triples method (CCSD(T)) together with the frozen-core approximation (Purvis & Bartlett 1982), hereafter referred as the CCSD potential. In the CCSD potential calculation, the doubly augmented correlation consistent basis set of Woon & Dunning (1994), d-aug-cc-pVnZ, n = T, Q, and 5, was adopted. The second surface was constructed using the complete active space self-consistent field method (Werner & Knowles 1985) with the internally contracted MRCI approach (Werner & Knowles 1989) based on aug-cc-pVQZ basis sets, hereafter referred to as the MRCI potential. Quantum-mechanical close-coupling (CC) computations of the cross sections for initial state j = 0 were carried out by Shepler et al. (2007) on the two PESs at kinetic energies of 400 and 800 cm⁻¹. The state-to-state cross sections on both PESs showed even-Δj propensity with odd-Δj transitions being significantly suppressed. In a subsequent investigation of the formation and excitation of CO in diffuse clouds (Liszt 2007), the rotational excitation temperature for the j = 1 → 0 transition was recomputed using the CO–H rotational excitation rate coefficients of Green & Thaddeus (1976). It was also pointed out by Liszt (2007) that the observations are better reproduced with these smaller CO–H rotational excitation rate coefficients.

As discussed in Shepler et al. (2007), a deficiency in the long-range part of the WKS PES results in the inaccurate rate coefficients. For odd Δj, the rate coefficients are typically too large, and if implemented in astrophysical models, can lead to questionable astrophysical deductions. However, due to the lack of more accurate rate coefficients for the CO–H system, the BYD values are still adopted. Examples include a study of molecular excitation in interstellar shock waves by Flower & Pineau des Forêts (2010), an investigation of the molecular regions of the supernova remnants IC443C, W28, W44, and 3C391 and the Herbig-Haro objects HH7 and HH54 (Yuan & Neufeld 2011), and in radiation thermo-chemical models of protoplanetary disks (Thi et al. 2013). Therefore, to provide accurate rate coefficients of CO rotational excitation due to impact by H atoms, we extend here the calculations of Shepler et al. (2007), carrying out computations of state-to-state quenching cross sections and rate coefficients for CO from initial states j = 1–5 in the ground vibrational state to all lower j' levels using the potentials presented in Shepler et al. (2007).

The quantum scattering theory of a linear rigid-rotor with a structureless atom was developed by Arthurs & Dalgarno (1963). In the current calculations, we treat CO as a rigid-rotor with bond-length fixed at its equilibrium distance. Calculations were carried out using the quantum CC method and the coupled-states (CS) approximation (see, for example, Flower 2007). The two-dimensional CO–H interaction potential is given by V(R, θ), where R is the separation between the CO center of mass and the H atom, and θ is the angle between R and the CO molecular axis. The PES was expressed as

$$\begin{equation} V(R,\theta)=\sum _{\lambda }^{\lambda _{\rm max}}v_{\lambda }(R)P_{\lambda }(\hbox{cos} \theta), \end{equation} \tag{ 1 }$$

where P_λ are Legendre polynomials of the order of λ.

For a rotational transition from an initial state j to a final state j', the degeneracy-averaged-and-summed integral cross section can be given, within the CC method, as

$$\begin{eqnarray} \sigma _{j\rightarrow j^{\prime }}(E_{j}) &=& \frac{\pi }{(2j+1)k_{j}^2}\sum _{J=0}(2J+1)\sum _{l=|J-j|}^{J+j} \sum _{l^{\prime }=|J-j^{\prime }|}^{J+j^{\prime }} \nonumber\\ && \times |\delta _{jj^{\prime }}\delta _{ll^{\prime }} -S_{jj^{\prime }ll^{\prime }}^J(E_j)|^2, \end{eqnarray} \tag{ 2 }$$

where $S_{jj^{\prime }ll^{\prime }}^J$ is an element of the scattering matrix and j and l denote the rotational angular momentum of the CO molecule and the orbital angular momentum of the collision complex, respectively. J is the total angular momentum given by J = l + j. $k_{j}=\sqrt{2\mu (E-\epsilon _{j})}/\hbar$ is the wave vector for the initial channel, E the total energy, _j the rotational energy of CO, E_j the relative translational energy for the initial channel, and μ the reduced mass of the CO–H system. In the CS formulation, the integral cross section is

$$\begin{eqnarray} \sigma _{j\rightarrow j^{\prime }}(E_{j}) &=& \frac{\pi }{(2j+1)k_{j}^2}\sum _{J=0}(2J+1)\sum _{\Omega =0}^{\Omega _{\bf max}} \nonumber\\ && \times (2-\delta _{\Omega 0})|\delta _{jj^{\prime }} -S_{jj^{\prime }}^{J\Omega }(E_j)|^2, \end{eqnarray} \tag{ 3 }$$

where Ω is the projection of the angular momentum quantum number of the diatom along the body-fixed axis.

The nonreactive scattering code MOLSCAT (Hutson & Green 1994) was used in all our calculations. The modified log-derivative Airy propagator of Alexander & Manolopoulos (1987), with a variable step-size, was adopted to solve the coupled-channel equations. The propagation was carried out from R = 1 a₀ to a maximum distance of R = 100 a₀, where a₀ is the atomic unit of length (the Bohr radius). Cross sections were calculated for collision energies between 10⁻⁵ and 15,000 cm⁻¹. The CC calculations have been carried out for initial kinetic energies up to 11,000 cm⁻¹ on the MRCI potential. Computations with the CS approximation were performed for the MRCI, CCSD, and WKS potentials above 2000 cm⁻¹, but CC calculations were performed for the latter two surfaces for collision energies less than 2000 cm⁻¹. The angular dependence of the potential expressed in Equation (1) was expanded to λ_max  =  20 with a 22-point Gauss–Legendre quadrature. In order to ensure convergence of the state-to-state cross sections, at least 5–10 closed channels in the basis and a sufficient number of angular momentum partial waves were included in our calculations. The CO rotational constant of 1.9225 cm⁻¹ (Lovas et al. 1979) was used to calculate the rotational energy levels.

3.1. State-to-state Deexcitation Cross Sections

Calculations of state-to-state rotational quenching cross sections have been carried out for initial levels j = 1, 2, ..., 5 for the deexcitation of CO(v = 0, j) by H.⁵ To illustrate the different behavior of the cross sections on the WKS, CCSD, and MRCI potentials, as an example we present in Figure 1, a comparison of the j = 1 → 0 deexcitation cross section. Due to the limited R range (3–20 a₀) of the CCSD potential, the quenching cross sections are given for collision energies greater than 0.1 cm⁻¹. From Figure 1, one can see that the cross sections on the WKS PES are significantly larger than that obtained using the CCSD and MRCI potentials, especially at ultracold collision energies. In the energy region between 0.1 and ∼20 cm⁻¹, the cross sections computed on all three potentials exhibit resonances caused by the van der Waals interaction, but due to differences in the potential well depths (see Shepler et al. 2007) very different structures are displayed, in particular the resonances on the WKS potential are largely suppressed. There exists generally good agreement between the cross sections on the CCSD and MRCI PESs above ∼100 cm⁻¹, while the cross sections on all three surfaces appear to converge above 10⁴ cm⁻¹. From Figure 1, it can also be seen that the CS approximation is in excellent agreement with the CC method.

Zoom In Zoom Out Reset image size

For initial rotational states j = 2, 3, 4, and 5, the state-to-state quenching cross sections were calculated using the MRCI PES of Shepler et al. (2007) only. Figure 2 presents the cross sections from initial j = 2, 3, and 4 into the individual final states j'. Over the entire collision energy range considered, similar structure is observed for the different initial rotational j states. At collision energies below ∼10⁻³ cm⁻¹, only s-wave (l = 0) scattering in the incident channel contributes, with the cross section for inelastic collisions varying inversely with the relative velocity, indicative of Wigner's threshold law (Wigner 1948). Similar to that shown in Figure 1, a number of resonances at collision energies between ∼1 and ∼30 cm⁻¹ are evident, but the resonances become suppressed compared with the j = 1 → 0 resonances (see also Yang et al. 2006). In general, the cross sections increase with increasing j', with that for j' = 0 being the smallest for initial levels j = 3 and 4. For each initial j state, the |Δj| = |j' − j| = 2 transition is seen to dominate the quenching. This finding is consistent with the fact that odd-Δj transitions in homonuclear molecules are forbidden and since CO is nearly homonuclear, it should follow an even-Δj propensity (e.g., Lee & Bowman 1987; Shepler et al. 2007).

Zoom In Zoom Out Reset image size

3.2. State-to-state Deexcitation Rate Coefficients

By thermally averaging state-to-state quenching cross sections over a Maxwellian distribution of collision energy, the state-to-state quenching rate coefficients for deexcitation from specific initial rotational states can be obtained at a temperature T. For the j = 1 → j' = 0 rate coefficients in Figure 3, we compare the present results on the CCSD and MRCI potentials with previous calculations of Chu & Dalgarno (1975), Green & Thaddeus (1976), and BYD. Significant differences between the present results and those of Chu & Dalgarno (1975) and BYD are evident. In particular, the BYD rate coefficients obtained using the WKS potential are over an order of magnitude larger than the present results. However, the results of Green & Thaddeus (1976) are in agreement with the present rate coefficients on the MRCI PES, though their adopted semiempirical potential is expected to be less accurate. As one can see from Figure 3, the current rate coefficients on the CCSD and MRCI potentials agree very well for temperatures above 150 K, though differences occur for lower temperatures. The cross section resonances, as shown in Figure 1, are responsible for the undulations in the rate coefficients.

Zoom In Zoom Out Reset image size

For j = 2, 3, 4, and 5, the state-to-state deexcitation rate coefficients were computed using the MRCI potential only and are presented in Figures 4 and 5 for temperatures ranging from 1 to 3000 K. For temperatures below ∼20 K, the rate coefficients display a less prominent temperature dependence due to suppression of the resonances in the cross sections. For temperatures above ∼20 K, the rate coefficients show a general increasing trend for all transitions. Further, the rate coefficients increase with decreasing Δj, except for Δj = −2 transitions which dominate. For those dominant transitions out of each initial level, 2 → 0, 3 → 1, 4 → 2, and 5 → 3, the rate coefficients of Chu & Dalgarno (1975), Green & Thaddeus (1976), and BYD are also displayed in Figures 4 and 5. Similar to Figure 3 for j = 1 → 0, Figures 4 and 5 show large differences between the present rate coefficients and those of BYD as well as those of Green & Thaddeus (1976). However, the rate coefficients of Chu & Dalgarno (1975) approach the current results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Most of the observed FIR/submm absorption and emission lines originate from rotational transitions due to simple molecular species such as CO. To accurately model molecular spectral features, in addition to spectroscopic data, a knowledge of collisional excitation rate coefficients due to the dominant impactors is required. As discussed in the Introduction, carbon monoxide is nearly ubiquitous, being detected in a vast number of objects with varying physical, chemical, and radiative properties. As a consequence, CO plays an important role in the thermal balance and chemistry and is a key target of many modeling efforts. For example, the radiative and collisional rates of CO were used by Martin et al. (2004) in estimating the gas density and kinetic temperature in a model of the Galactic center molecular region. In a study of the observed properties of Galactic diffuse gas, Pety et al. (2008) used the j = 1 → 0 and 2 → 1 rotational transitions of CO to perform non-LTE models of CO brightness, while Srianand et al. (2008) detected CO absorption in a damped Lyα system at high redshift. From their modeling with the radiative transfer code RADEX (van der Tak et al. 2007), Srianand et al. (2008) found that the cosmic microwave background dominates the CO excitation and that the rotational excitation temperatures of CO are higher than measured in the Galactic ISM. More highly excited CO transitions, from j = 14 to 30, have been observed with Herschel/PACS in the extragalactic source NGC 1068 (Hailey-Dunsheath et al. 2012).

To investigate the excitation mechanism of CO rotational and rovibrational emission seen in Herbig Ae disks (for example, the Herschel/PACS observations of Meeus et al. 2012), Thi et al. (2013) have constructed a complete CO rovibrational model which adopted existing rate coefficients of CO in collision with H, He, H₂, and electrons into the radiative-photochemical code ProDiMo. Due to the incompleteness of the available rate coefficients, the missing values were estimated using various scaling laws. For the case of CO–H, the rate coefficients of BYD were used, though they are not reliable as pointed by Shepler et al. (2007) and as indicated in Figures 3–5. From modeling of protoplanetary disks, Thi et al. (2013) found that the abundance of atomic H in the CO line-emitting region is significant so that CO–H collisions are important.

As CO collisional rate coefficients are needed for CO line intensity and abundance profile modeling, the current CO–H rate coefficients will be very useful for interpreting observed CO emission lines. In particular, recent observations with Herschel/PACS and Herschel/HIFI provide a large amount of CO FIR/submm data for warm gas around evolved stars. Modeling the CO emission with accurate collisional rate coefficients will enhance our ability to understand the evolution of these objects. The availability of the new rotational CO–H rate coefficients provided here, along with our earlier CO–H₂ (Yang et al. 2010) results, should improve CO modeling for the environments discussed above and many others. Further, the excitation rates due to H and H₂ become comparable for abundance ratios $n_{\rm H} /n_{\rm H_2}$ of 2–3 at 10 K and 0.5–1 at 3000 K, highlighting the importance of H collisions in regions with a significant atomic hydrogen fraction.

We have performed explicit quantum-mechanical calculations of rotational quenching of CO due to H collisions using the CC approach and the CS approximation on two different interaction potential surfaces. Deexcitation cross sections reveal even-Δj propensity for the CO–H system over the collision energy range considered as expected on physical grounds. State-to-state quenching rate coefficients for initial CO rotational levels j =1, 2, ..., 5 to all lower j' levels were computed for temperatures ranging from 1 to 3000 K and are available in tables formatted for astrophysical applications. Discrepancies in the cross sections and rate coefficients on the two new potentials are due to the differences in the well depth and anisotropy of the surfaces. Comparison of the present rate coefficients with previous theoretical results show significant differences primarily due to the quality of the adopted potentials. The present CO–H rate coefficients can be used to aid astrophysical modeling, while quenching from higher rotational excited states as well as vibrational transitions will be studied in the future.

B.H.Y., P.C.S., and J.M.B. acknowledge support from NASA under grant No. NNX12AF42G, while the work of N.B. and R.C.F. was partially supported by the NSF through grant PHY-1205838 and grant PHY-1203228, respectively.