ROTATIONAL QUENCHING OF CO DUE TO H₂ COLLISIONS

Molecular hydrogen (H₂) and carbon monoxide (CO) are the most abundant molecular species in the majority of interstellar environments, and they play important roles in determining the physics and chemistry of diffuse clouds (Herbst 2001). CO is observed in diffuse clouds in absorption in the ultraviolet (UV) and visible and in emission at infrared (IR) wavelengths (Snow & McCall 2006), revealing the physical and chemical complexity along these lines of sight. As an example, Liszt (2006, 2007) recently investigated the formation, fractionation, and rotational excitation of CO in H₂ bearing diffuse H i clouds. Due to its small rotational constant, CO can be easily rotationally excited by collisions with other species in interstellar gas, mostly H₂, with CO rotational lines providing important diagnostics of gas density and temperature. The abundance ratio, CO/H₂, is assumed to be roughly constant in dense molecular gas, and as such, CO is often used as a tracer of H₂, as the latter is difficult to detect in emission in cold environments. The observed CO abundance is therefore used to estimate the total H₂ content (Sonnentrucker et al. 2007).

In environments with an intense UV field, 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. In recent years, various codes that model the physics and chemistry of PDRs have been developed with an emphasis on the excitation mechanisms of both H₂ and CO (van Dishoeck & Black 1988; Warin et al. 1996; Le Petit et al. 2006; Shaw et al. 2005; Röllig et al. 2007).

In more dense environments, such as low-mass dwarf stars (i.e., M, L, or T dwarfs), CO is an important opacity source and its infrared line list can be obtained quite accurately (Pavlenko & Jones 2002; Jones et al. 2005). It has also been observed in the day side spectrum of the transiting hot-Jupiter extrasolar giant planet (EGP) HD 189733b (Swain et al. 2009). While molecular level populations are typically treated in local thermodynamic equilibrium (LTE) in stellar and planetary atmospheres, Schweitzer et al. (2000) have pointed out that non-LTE conditions may exist in very cool stellar atmospheres or in cases where the temperature of the incident radiation, i.e., close-in EGPs, is much different from the kinetic temperature of the gas.

Due to the astrophysical importance of H₂ and CO, the CO–H₂ collisional system has been the subject of numerous experimental (Kudian et al. 1967; Butz et al. 1971; Kupperman et al. 1973; Bouanich & Brodbeck 1973; Nerf & Sonnenberg 1975; Brechignac et al. 1980; Andresen et al. 1982; Picard-Bersellini et al. 1983; Schramm et al. 1991; Drascher et al. 1998; McKellar 1990, 1991, 1998; Antonova et al. 2000) and theoretical (Green & Thaddeus 1976; Prisette et al. 1978; Flower et al. 1979; Poulsen 1982; Schinke et al. 1984; Bacic et al. 1985a, 1985b; Jankowski & Szalewicz 1998, 2005; van Hemert 1983; Parish et al. 1992; Danby et al. 1993; Salazar et al. 1995; Reid et al. 1997; Mengel et al. 2001; Baker & Flower 1984; Flower & Launay 1985; Flower 2001; Antonova et al. 2000; Wernli et al. 2006; Yang et al. 2006a, 2006b) studies. Quantitative determinations of state-to-state cross sections and rate coefficients for CO–H₂ collisions are crucial to numerical models of various astrophysical environments, such as those highlighted above. However, as measurements of these quantities are difficult, numerical models often rely on cross sections and rate coefficients derived from theoretical methods. In a highly cited work, Green & Thaddeus (1976) performed close-coupling (CC) calculations of rate coefficients based on an approximate CO–H₂ potential surface. Since then, a number of quantum scattering calculations were carried out on various potential energy surfaces (PESs). In this work, we have performed comprehensive calculations of state-to-state cross sections in an effort to develop a complete database of rotational quenching rate coefficients, but utilizing a more recent, and presumably more accurate PES. In the following sections, the choice of the potential surface, the adopted scattering approach, and the results are discussed. We briefly highlight important astrophysical applications of the current CO–H₂ rate coefficients for modeling far IR (FIR) observations.

A number of PESs (Prisette et al. 1978; Flower et al. 1979; Poulsen 1982; Schinke et al. 1984; Bacic et al. 1985a; Jankowski & Szalewicz 1998, 2005) have been developed for the CO–H₂ complex. One of the most accurate surfaces was given by Jankowski & Szalewicz (1998) who calculated a four-dimensional, rigid-rotor PES referred to as V₉₈. This surface was constructed using the symmetry-adapted perturbation theory (SAPT) method with high-level electron correlation effects. An analytical fit of the ab initio PES has a global minimum of 109.3 cm⁻¹ at the intermolecular separation of 7.76 a₀ for the linear geometry with the C atom pointing toward H₂. Antonova et al. (2000) adopted V₉₈ in their calculations of state-to-state cross sections for CO rotational excitation at collision energies between 795 and 991 cm⁻¹ which were found to be in good agreement with their measurements. The V₉₈ PES was also utilized in second virial coefficient calculations for mixtures of hydrogen and CO (Gottfried & McBane 2000). A comparison with experimental data suggested that the van der Waals well of the V₉₈ PES is too deep by 4%–9%, though it represents an improvement over earlier surfaces. The V₉₈ surface has also been used in full coupled-channel cross section and rate coefficient calculations for rotationally inelastic scattering of CO by para- and ortho-H₂ (Flower 2001; Mengel et al. 2001). However, because the attractive well of V₉₈ was presumed to be too deep, Mengel et al. (2001) modified the PES by merely multiplying the surface by a constant factor of 0.93 which was subsequently used in their scattering calculations.

A newer CO–H₂ PES was reported by Jankowski & Szalewicz (2005), referred to as V₀₄. To improve the surface accuracy, they adopted the coupled-cluster method with single, double, and non-iterative triple excitations [CCSD(T)] and the supermolecular approach. The V₀₄ surface was obtained on a five-dimensional grid including the dependence on the H–H separation, but with the CO molecule taken to be rigid with the C–O separation set to the value of the C–O distance averaged over the CO ground state vibrational wave function. The PES was then obtained by averaging over the intramolecular vibration of H₂ to yield a four-dimensional rigid-rotor surface. Similar to V₉₈, an analytical fit of V₀₄ was presented which has a global minimum of 93.05 cm⁻¹, or 16.22 cm⁻¹ shallower than V₉₈. Jankowski & Szalewicz (2005) used V₀₄ to calculate the rovibrational energy levels for the para- and ortho-H₂–CO complex as well as the second virial coefficient. The rovibrational energies were found to agree with the experimental values of McKellar (1998) for para-H₂ to within 0.1 cm⁻¹ while a scaling of V₀₄, resulting in a deepening of the well depth by 4 cm⁻¹, was needed to match the measured second virial coefficients. V₀₄ may be too shallow though it is possible that the experimental second virial coefficients are systematically too low (Jankowski & Szalewicz 2005). Comparative studies of collisional cross sections using the two surfaces have been presented in our earlier work (Yang et al. 2006b) and by Wernli et al. (2006). In Yang et al. (2006b), calculations with the V₀₄ surface gave cross sections in excellent agreement with the state-to-state measurements of Antonova et al. (2000).

The theory for the scattering of two linear rigid rotors has been developed in Green (1975) and Heil et al. (1978). The calculations presented here were performed for collision energies between 10⁻⁵ and 10,000 cm⁻¹ by applying both the CC method and the coupled-states (CS) approximation. All the CC and CS calculations were performed using the non-reactive scattering code MOLSCAT (Hutson & Green 1994). In the quantum scattering calculations, the coupled-channel equations were integrated using the modified log-derivative Airy propagator of Alexander & Manolopoulos (1987) with a variable step size. The largest Legendre terms adopted in the potential expansion for H₂ and CO were, respectively, 8 and 10. The numbers of Gauss integration points used in projecting angular components of the potential were 10, 11, and 12 for integration in θ₁, θ₂, and ϕ, respectively. The propagation was carried out to a maximum intermolecular separation of R = 100 Å. At each energy, a sufficient number of total angular momentum partial waves was included to ensure convergence of the state-to-state cross sections to within 1%. The maximum value of the total angular momentum quantum number J employed in the calculations was 300. Since the computation time for CC calculations scales as j⁶_max, the CS method is adopted for energies greater than 800 cm⁻¹ which results in considerable time savings with a scaling proportional to j⁴_max, where j_max is the maximum size of the rotational basis. Further details about the scattering calculations can be found in Yang et al. (2006a, 2006b). The adopted rotational constants for H₂ and CO are 60.853 cm⁻¹ (Huber & Herzberg 1979) and 1.9225 cm⁻¹ (Lovas et al. 1979), respectively. Hereafter, j₁ denotes the rotational quantum number for H₂, and likewise j₂ for CO. In the present study, we adopted the four-dimensional PES for the CO–H₂ complex, V₀₄, of Jankowski & Szalewicz (2005).

Calculations of the collision energy dependence of state-to-state quenching cross sections for initial rotational states of CO, j₂ = 1, 2, ..., 40, by collisions with both para- and ortho-H₂ (i.e., j₁ = 0, 1, respectively) were performed. The cross sections were thermally averaged over the kinetic energy distribution to yield state-to-state rate coefficients from specific initial rotational states j₂ as functions of the temperature T,

$$\begin{equation} k_{j_2\rightarrow j_2^{\prime }}(T)= \left(\frac{8}{\pi \mu \beta } \right)^{1/2}\!\!\beta^{2}\!\!\int ^{\infty }_0 E_k \sigma _{j_2\rightarrow j_2^{\prime }}(E_k)\exp (-\beta E_k)dE_k{,} \end{equation} \tag{ 1 }$$

where $\sigma _{j_2\rightarrow j_2^{\prime }}$ is the rotational transition cross section with j₂ and j'₂ being respectively the initial and final rotational quantum number of CO, μ is the reduced mass of the CO–H₂ complex, E_k the kinetic energy, and β = (k_BT)⁻¹, where k_B is the Boltzmann constant.

All state-to-state cross sections and rate coefficients for quenching are available on the UGA Molecular Opacity Project Web site (http://www.physast.uga.edu/ugamop/). The rate coefficients are also available in the format of the Leiden Atomic and Molecular Database (LAMDA, Schöier et al. 2005) on our Web site and the LAMDA site (http://www.strw.leidenuniv.nl/~moldata/). Figure 1 shows typical examples of state-to-state cross sections for j₂ = 5 and 20 for collisions with para-H₂. A number of resonances are evident between 10⁻² and 10² cm⁻¹ which result in significant modulation in the rate coefficients. However, the magnitude of the resonances are seen to decrease with j₂. Cross sections for some other j₂ can be found in Yang et al. (2006a, 2006b). The agreement between the CC and CS calculations is shown to be excellent with a difference typically better than ∼20%, justifying the adoption of the CS approximation at the higher energies.

Zoom In Zoom Out Reset image size

Quenching rate coefficients from initial rotational states, j₂ = 5, 10, 20, and 40, at temperatures ranging from 1 to 3000 K are shown in Figures 2–9 for CO scattering with para- and ortho-H₂. Unfortunately, we are unaware of any experimental rate coefficient data for these initial states. Therefore, the current results are compared with the theoretical rate coefficients of Flower (2001), which were obtained over the limited temperature range, 5–400 K, with the V₉₈ PES. Similar results with V₀₄ were obtained by Wernli et al. (2006) for rotational levels of CO up to 5 and temperatures in the range 5–70 K. Wernli et al. also presented an analytic fit valid in the same temperature range. The rate coefficients calculated using their analytic relation are almost identical to the current results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the quenching of j₂ = 5 due to para-H₂ collisions (Figure 2 shows that for temperatures between ∼1 and 100 K, which is the van der Waals interaction-dominated regime) the rate coefficients exhibit an undulatory temperature dependence for Δj₂ = −1, −2, and −3 transitions due to the presence of resonances,⁴ though the magnitude of the undulations decrease with increasing j₂. At temperatures above ∼100 K, the rate coefficients generally increase with increasing temperature. Comparison with the rate coefficients of Flower shows that at temperatures higher than ∼50 K, there is generally good agreement, except for j'₂ = 0 and 2, where Flower's results get larger at temperatures above about 200 K. The results of Flower are smaller than the present rate coefficients for lower temperatures with the discrepancy increasing with decreasing temperature. Exceptions to this behavior are the transitions j₂ = 5 → j'₂ = 2 and 4, in which Flower's rate coefficients are generally larger than the present results. For quenching from initial levels j₂ = 10 and 20 given in Figures 3 and 4, respectively, the state-to-state rate coefficients have a structure similar to that shown for j₂ = 5. The rate coefficients are dominated by small Δj₂. For Δj₂ = −1 transitions, a slight bump can be seen at temperatures between ∼1 and 100 K for j₂ = 10 and 20 as shown in Figures 3 and 4. Figure 5 displays quenching rate coefficients for j₂ = 40. No previous results are available as the calculations of Flower (2001) stopped at j₂ = 29 for para-H₂ collisions.

For scattering by ortho-H₂, the trends noted for para-H₂ are also evident as displayed in Figures 6–9. The quenching rate coefficients are of similar magnitude as those obtained for para-H₂ and show very similar structures. The highest level considered by Flower (2001) for ortho-H₂ collisions was j₂ = 20.

Comparing Figures 2–9, some general trends can be noted. With few exceptions, the state-to-state rate coefficients appear to attain constant values as they decrease to 1 K. However, this is deceptive as an orbiting resonance is typically present in the cross sections for collision energies between ∼0.01 and 1 cm⁻¹ (see Figure 1) computed on the V₀₄ PES (Yang et al. 2006b). For temperatures less than shown, the rate coefficients typically decrease (increase) prior to the Wigner regime (Wigner 1948) for intermediate (high) j₂ before eventually attaining finite limiting values, typically below ∼1 mK.

Disregarding the effect of the resonances, one would expect two possible trends in the relative ordering of the rate coefficient values as a consequence of energy and angular momentum gaps. The exponential energy gap law (Steinfeld et al. 1991) predicts that the rate coefficients should decrease with increasing internal energy difference between the initial and final levels (i.e., with increasing |Δj₂|). However, CO is nearly homonuclear and rotational transitions in homonuclear molecules follow an even Δj₂ propensity rule (e.g., Shepler et al. 2007), though the effect may diminish for large values of |Δj₂|. Both effects appear to be important for CO–H₂. For example, Δj₂ = −2 is the dominant transition for quenching for temperatures greater than ∼1000–3000 K. It is also the dominant transition for temperatures less than ∼10 K, but only for j₂ = 5 and 10. For other temperatures, Δj₂ = −1 transitions are predominant, while for some cases Δj₂ = −3 is secondary, but with a rate coefficient larger than the Δj₂ = −2 transitions. However, for Δj₂ < −3, the exponential energy gap law appears to be fulfilled in nearly all considered cases.

As discussed in the Introduction, CO plays an important role in a variety of astrophysical and atmospheric environments and has been detected in countless sources under a variety of excitation conditions. We discuss a few of these below focusing primarily on observed and predicted transitions from intermediate- and high-j₂ in the FIR.

Rotational emission lines of CO have been observed in high-redshift (z) objects including the quasar SDSS J1148+5251 at z = 6.42 using the Very Large Array. Narayanan et al. (2008) modeled the emission from CO(v = 0, j₂) with j₂ = 1–10 with a non-LTE radiative transfer code and considered excitation due to H₂ with rate coefficients taken from the LAMDA database. They found that the CO spectral energy distribution peaks at j₂ = 8 at the beginning of the quasar phase, z ∼ 7–8. Observations and modeling of the CO emission can serve to probe the hierarchical buildup of the host galaxy and the quasar phase and can be used to estimate the halo mass and quasar host morphology.

Spaans & Meijerink (2008) proposed that CO emission from j₂ as high as 30 could be observed from high-z massive accreting black holes using the Atacama Large Millimeter/submillimeter Array (ALMA). Such observations could provide a diagnostic of the radiation source: PDR or X-ray dominated region (XDR), for example. However, as the densities of such objects are typically 10⁵ cm⁻³, collisions with the dominant species, H₂, will determine the CO rotational populations. In Figure 10, the critical densities for CO due to para-H₂ collisions are plotted (neglecting optical depth effects) where the critical density of rotational level j₂ is defined as

$$\begin{equation} n_{\rm cr}(j_2) = {{\sum _{j_2^{\prime }<j_2} A_{j_2\rightarrow j_2^{\prime }}}\over {\sum _{j_2^{\prime }\ne j_2} k_{j_2\rightarrow j_2^{\prime }}(T)}} \end{equation} \tag{ 2 }$$

(e.g., Osterbrock & Ferland 2006) and $A_{j_2\rightarrow j_2^{\prime }}$ are the spontaneous transition probabilities for dipole transitions. Rotational levels j₂ ≳ 10 are seen to be clearly out of thermal equilibrium for a gas density of 10⁵ cm⁻³ requiring a non-LTE analysis. As the LAMDA database incorporates the CO–H₂ rate coefficients of Wernli et al. (2006) and Flower (2001) with extrapolations for para- and ortho-H₂ collisions above j₂ = 29 and 20, respectively, the predicted CO line intensities of Spaans & Meijerink (2008) for emission from levels larger than j₂ = 20 would likely be improved if the current rate coefficients are adopted. Interestingly, their predicted CO SEDs peak at j₂ = 10 and j₂ = 25 for PDR and XDR environments, respectively.

Zoom In Zoom Out Reset image size

The j₂ = 6 → 5 and j₂ = 7 → 6 lines have been observed in the starburst galaxy NGC 253 (Hailey-Dunsheath et al. 2008; Bradford et al. 2003). Using a non-LTE analysis, Bradford et al. (2003) concluded that the cooling in the CO lines is so large that it must be balanced by a cosmic ray heating rate ∼800 times greater than that in the Milky Way. Previously, Krolik & Lepp (1989) suggested that highly excited j₂ lines (e.g., j₂ = 16–58) could be detectable from Seyfert galaxies giving diagnostics of the internal pressure of the dust-obscured torus.

The Infrared Space Observatory studied a large variety of pre-main sequence objects (i.e., protostars) and detected numerous FIR CO rotational lines. For example, emission lines from j₂ = 14–26 were observed in the Class 0 object L1448-mm (Nisini et al. 1999), j₂ = 14–31 in the Class I object IC 1396N, j₂ = 14–16 in the Class I object W28 A2, j₂ = 14–20 in the Class II object R CrA (Saraceno et al. 1999), and j₂ = 14–25 in T Tauri (Spinoglio et al. 1999). Non-LTE analyses of these pre-main sequence objects suggested multiple components with various excitation conditions. Future studies of these protostars, the infrared sources discussed above, and future FIR and submillimeter observations with Herschel and ALMA can benefit from the CO–H₂ quenching rate coefficients presented here.

Rotational quenching of CO due to para- and ortho-H₂ collisions has been studied using an explicit quantum-mechanical CC approach and the CS approximation on the potential surface, V₀₄, of Jankowski & Szalewicz (2005). State-to-state quenching rate coefficients for initial rotational levels j₂ = 1, 2, ..., 40 of CO are obtained over a wide temperature range and available in tables formatted for astrophysical applications. Resonances result in undulations in the temperature dependence of the rate coefficients with the amplitudes of the undulations typically decreasing with j₂. For temperatures less than ∼50 K, the current state-to-state rotational quenching rate coefficients obtained with V₀₄ are found to depart from the results of Flower (2001) obtained with the earlier V₉₈ PES. The discrepancies are likely related to the differences in the well depth and anisotropy of the two potentials.

B.H.Y. and P.C.S. acknowledge support from NASA grant NNX07AP12G. N.B. acknowledges support from the National Science Foundation (NSF) through grant PHY-0855470. R.C.F. acknowledges support from NSF through grant PHY-0854838. We thank Dr. Piotr Jankowski for providing his CO–H₂ potential subroutine and Prof. Gary Ferland and Dr. Donghui Quan for discussions on astrophysical applications. Part of calculations were performed on EITS/RCC pcluster at the University of Georgia.