Inelastic cross sections and rate coefficients for collisions between CO and H 2 Christina Castro,1 Kyle Doan,1 Michael Klemka,1 Robert C. Forrey,1,a) B. H. Yang,2 Phillip C. Stancil,2 and N. Balakrishnan3 1)Department of Physics, Penn State University, Berks Campus, Reading, PA 19610-6009 7 2)Department of Physics and Astronomy and the Center for Simulational Physics, 1 0 University of Georgia, Athens, Georgia 30602 2 n 3)Department of Chemistry, University of Nevada, Las Vegas, a J NV 89154 8 1 (Dated: 20 January 2017) ] M A five-dimensional coupled states (5D-CS) approximation is used to compute cross I . h sections and rate coefficients for CO+H collisions. The 5D-CS calculations are 2 p - benchmarked against accurate six-dimensional close-coupling (6D-CC) calculations o r t for transitions between low-lying rovibrational states. Good agreement between the s a two formulations is found for collision energies greater than 10 cm−1. The 5D-CS [ 1 approximation is then used to compute two separate databases which include highly v 3 excited states of CO that are beyond the practical limitations of the 6D-CC method. 1 2 The first database assumes an internally frozen H molecule and allows rovibrational 2 5 0 transitions for v ≤ 5 and j ≤ 30, where v and j are the vibrational and rotational . 1 0 quantum numbers of the initial state of the CO molecule. The second database 7 1 allows H rotational transitions for initial CO states with v ≤ 5 and j ≤ 10. The two 2 : v databases are in good agreement with each other for transitions that are common i X to both basis sets. Together they provide data for astrophysical models which were r a previously unavailable. PACS numbers: 34.10.+x, 34.50.Ez a)Electronic mail: [email protected] 1 I. INTRODUCTION H and CO are the most abundant molecules in most astrophysical environments and 2 have been the focus of numerous astrophysical studies and observations1. CO is easily excited by collisions with H and other species in interstellar gas, and the resulting emission 2 lines are commonly used to provide important diagnostics of gas density and temperature. When these environments are irradiated with an intense UV field, the radiation drives the chemistry and internal level populations out of equilibrium. Models which aim to interpret the emission lines arising from these environments must account for all mechanisms which can excite or de-excite the molecules. Simulation packages have been developed2–5 which account for inelastic transitions of both H and CO in various astrophysical environments 2 such as photodissociation regions (PDRs) which occur between hot HII (ionized hydrogen) and cold molecular regions. Observations of star-forming regions and protoplanetary disks (PPDs) of young stellar objects have shown evidence of rovibrational transitions involving states of CO which are highly excited. In particular, vibrational transitions have been detected in the near infrared by the Infrared Space Observatory, the Gemini Observatory, and the Very Large Telescope (VLT). A recent VLT survey6 of 69 PPDs found that 77% of the sources showed CO vibrational bands, including vibrational level v as high as 4 and rotational level j as high as 32. CO vibrational bands in the 1-5 micron regime will soon be accessible by NASA’s James Webb Space Telescope which will allow the warm inner regions of young stellar objects to be probed. To reliably model these environments requires an extensive set of state-to-state rate coefficients for rovibrational transitions induced by H 2 collisions. Many of these rate coefficients are either unavailable or else are estimated with potentially unreliable methods7. Toprovideaccuratedataforastrophysical modelsgenerallyrequires alargeandconcerted effort between several communities within atomic and molecular physics. Experimental data are sparse but are critical for benchmarking theoretical methods8–10. Numerical solution of the Schro¨dinger equation offers the best hope for generating the bulk of the needed data. However, to make these computations possible, it is often necessary to invoke decoupling approximations which reduce the dimensionality of the full scattering problem. Even when full-dimensional quantum dynamics calculations are feasible, their accuracy depends on the underlying electronic structure calculations and potential energy surface (PES) which is 2 used as input. Recently, a high-level full-dimensional PES was developed for CO+H and 2 used in a six-dimensional close-coupling (6D-CC) formulation10 which is essentially an exact treatment of the dynamics. The calculations were in excellent agreement with experiment for vibrational de-excitation of CO, whereas prior calculations which consisted of various levels of approximation varied by more than two orders of magnitude11–14. Subsequently, the 6D-CC approach was used to provide additional high-quality cross sections for CO+H 2 collisions15. These calculations are computationally demanding and have been completed for only a small subset of rovibrational levels needed for the astrophysical models. Nevertheless, these calculations provide the best benchmarks for transitions where experimental data are lacking, and may serve as a foundationfor testing approximate schemes which aim to further extend the available data. One such scheme is the coupled states or centrifugal sudden (CS) approximation which requires significantly less computational effort due to the decoupling of orbital and internal angular momenta. A recent study16 compared 6D-CS and 6D-CC cross sections for CO+H and found encouraging results for energies above 10 cm−1. Additional 2 decoupling of the “twist” angle between the two molecules led to a 5D-CS approximation which also yielded encouraging results. In the present work, we compute a large database of cross sections and rate coefficients using the 5D-CS approximation. The 5D-CS calculations are compared with existing 6D-CC results in order to assess the expected level of accuracy of the approximation. The database provides a balance between accuracy and computational efficiency and is made available in a convenient format for use in astrophysical models. II. THEORY Thequantum mechanical CCand CSformulations fordiatom-diatomcollisions have been given previously17–19. In order toclarifythe differences between the 5Dand6Dformulations, we provide a brief summary of the theory. The Hamiltonian of the four-atom system may be written H(~r ,~r ,R~) = T(~r )+T(~r )+T(R~)+V(~r ,~r ,R~), (1) 1 2 1 2 1 2 where T(R~) is a radial kinetic energy term describing the center-of-mass motion, T(~r ) and 1 T(~r ) are kinetic energy terms for the diatomic molecules, and V(~r ,~r ,R~) is the potential 2 1 2 3 energy of the system. It is convenient to define V(~r ,~r ,R~) = U(~r ,~r ,R~)+V(~r )+V(~r ), (2) 1 2 1 2 1 2 where V(~r ) and V(~r ) are the two-body potential energies of the isolated CO and H 1 2 2 molecules, and U(~r ,~r ,R~) is the four-body interaction potential which vanishes at large 1 2 separations. The 6D Jacobi coordinate system in Figure 1 is used where R is the distance between the centers-of-mass of the diatomic molecules, θ is the angle between ~r and R~, 1 1 θ is the angle between ~r and R~, and φ is the out-of-plane dihedral angle or “twist” angle. 2 2 The interaction potential may be expanded as U(~r ,~r ,R~) = A (r ,r ,R)Y (rˆ ,rˆ ,Rˆ) (3) 1 2 λ1,λ2,λ12 1 2 λ1,λ2,λ12 1 2 all λ X with Y (rˆ ,rˆ ,Rˆ) = hλ m λ m |λ m iY (rˆ )Y (rˆ )Y∗ (Rˆ), (4) λ1,λ2,λ12 1 2 1 λ1 2 λ2 12 λ12 λ1mλ1 1 λ2mλ2 2 λ12mλ12 all m X where h...|...i represents a Clebsch-Gordan coefficient and Y (rˆ) is a spherical harmonic. λm The total wave function for the four-atomsystem is expanded in terms of a diabatic basis set which contains products ofmolecular wave functions χ (r ) withvibrational androtational viji i quantum numbers v and j , respectively. We describe the combined molecular state (CMS) i i comprised of CO(v ,j ) and H (v ,j ) using the notation n = (v ,j ,v ,j ) so that the basis 1 1 2 2 2 1 1 2 2 functions for vibrational motion may be written χ (r ,r ) = χ (r )χ (r ) . (5) n 1 2 v1j1 1 v2j2 2 Therotationalwave functions aregiven interms ofproductsofspherical harmonicsinatotal angular momentum representation. The basis sets are defined by the number of vibrational levels and the maximum rotational level jmax for a given v. The radial interaction potential v matrix elements are obtained by integrating over the internal coordinates ∞ ∞ Bnλ;1n,λ′2,λ12(R) = χn(r1,r2)Aλ1,λ2,λ12(r1,r2,R)χn′(r1,r2)r12r22dr1dr2 . (6) Z0 Z0 The full potential matrix depends on the scattering formulation. For the CC method, the channels are defined by the set {n,j ,l}, where l is the orbital angular momentum quantum 12 number and ~j =~j +~j . The total angular momentum quantum number J is defined by 12 1 2 J~ =~l +~j , and the interaction potential matrix for 6D-dynamics is given by 12 4 UnJj12l;n′j1′2l′(R) = (4π)−3/2 (−1)j1′+j2′+j12+J [j1][j2][j12][l][j1′][j2′][j1′2][l′][λ1][λ2][λ12]2 1/2 aXll λ (cid:16) (cid:17) j j j l′ λ l j′ λ j j′ λ j l l′ λ 12 2 1 × 12 1 1 1 2 2 2 12  j′ j′ j′ Bλ1,λ2,λ12(R) (7)  0 0 0 0 0 0  0 0 0  j1′2 j12 J  12 2 1  n;n′     λ λ λ     12 2 1 which is diagonal with respect to J and independent of the projection ofJ~ in the space-fixed frame. The notations (===), {===}, and {≡≡≡} are the usual 3j, 6j, and 9j symbols, and [j] = (2j +1). For the 5D-CS formulation, the potential matrix is given by Um1,m2(R) = (4π)−3/2 (−1)λ1+λ2+m1+m2 [j ][j ][j′][j′][λ ][λ ][λ ]2 1/2 n;n′ 1 2 1 2 1 2 12 aXll λ (cid:16) (cid:17) j′ λ j j′ λ j j′ λ j j′ λ j λ λ λ × 1 1 1 2 2 2 1 1 1 2 2 2 1 2 12 Bλ1,λ2,λ12(R)       n;n′ 0 0 0 0 0 0 −m 0 m −m 0 m 0 0 0 1 1 2 2             (8) which is diagonal with respect to m and m , the projection quantum numbers of~j and~j , 1 2 1 2 respectively. The 5D-CS formulation averages the twist-angle dependence of the PES and yields a potential matrix which is independent of j . Both dynamical formulations require 12 the solution of a set of coupled radial equations derived from the Schro¨dinger equation. Cross sections may be expressed in terms of the appropriate T-matrix by π 2 σ6D-CC(E ) = (2J +1) TJ (E ) , (9) n→n′ c (2j +1)(2j +1)2µE nj12l;n′j1′2l′ c 1 2 c j12Xj1′2ll′J (cid:12) (cid:12) (cid:12) (cid:12) and (cid:12) (cid:12) σ5D-CS(E ) = π (2¯l+1) T¯lm1m2(E ) 2 , (10) n→n′ c (2j +1)(2j +1)2µE n;n′ c 1 2 c ¯lmX1m2 (cid:12) (cid:12) (cid:12) (cid:12) where E is the collision energy, and µ is the reduced mass of (cid:12)the CO-H (cid:12)system. The CS c 2 approximationassumes theoff-diagonalCoriolismatrixelements ofˆl2 withrespect tom and 1 m may beneglected, andthediagonal elements maybeapproximated byaneffective orbital 2 angular momentum quantum number ¯l which replaces l. For all calculations in the present work, l ≡ J which is its average value between |J−j | and J+j . Rate coefficients may be 12 12 obtained by thermally averaging the cross sections over a Maxwellian velocity distribution 8k T ∞ kn→n′(T) = B (kBT)−2 σn→n′(Ec)e−Ec/kBTEcdEc , (11) s πµ Z0 where T is the temperature and k is Boltzmann’s constant. B 5 III. RESULTS Previouswork16 showedthat5D-CSand6D-CSresultsforCO+H werevirtuallyidentical 2 for collision energies above 10 cm−1. The agreement with 6D-CC calculations was also found tobegoodattheseenergies forthelimitedamountofdatathatwasavailableforcomparison. Recently, new 6D-CC results have been reported15 which allow for more extensive testing of the 5D-CS approximation. All calculations in the present work were performed using a modified version of the TwoBC code20 which replaces equations (7) and (9) with equations (8) and (10). The radial coordinates of both molecules were represented as discrete variables with 20 points each. Gauss-Legendre quadratures were used for θ and θ with 14 points 1 2 each, and a Chebyshev quadrature was used for φ with 8 points. The maximum values for λ and λ were 10 and 6, respectively. The log-derivative matrix propagation method21,22 1 2 was used to integrate the set of coupled equations derived from the Schro¨dinger equation from R = 4 − 18 a.u. in steps of 0.05 a.u. The calculations were performed in unit steps for four decades of collision energy on a logarithmic energy grid. The maximum effective orbital angular momentum for each set of calculations is given by 20, E = 1−10 cm−1 , c  lmax = 4800,, EEcc == 11000−−1100,0c0m0−c1m,−1 , (12) Excellent agreement was seen at t1h6e0,EEbco=un1d,a0r0ie0s−, a1n0d,0co0n0vcemrg−e1nc.e tests at low E verified c c that the results were converged to within 5% or better with respect to size of the basis set. The 5D-CS calculations were performed for collision energies between 10-10,000 cm−1 and compared against 6D-CC results obtained perviously15. Figure 2 shows elastic and rotationally inelastic cross sections for the (1,1,0,0) initial state. The curves for rotational excitation (1,2,0,0) and de-excitation (1,0,0,0) are in near perfect agreement with the corresponding 6D-CC points for energies above 100 cm−1. The agreement remains good for energies down to 10 cm−1, however, the 5D-CS results do not show the resonant structure that was found in the 6D-CC calculations. This is partly due to the energy step size, which is too large to resolve the resonances, however, it is not expected that the effective l used in the CS calculations would accurately describe the resonances even if a finer energy grid were used. The elastic (1,1,0,0) cross section shows a similar pattern of agreement with the 6 exception of small discrepancies at 500 and 1000 cm−1. This is probably due to incomplete convergence of the 6D-CC elastic cross section with respect to higher partial waves. Figure 3 shows a set of similar comparisons for the (1,2,0,0), (1,3,0,0), (1,4,0,0), and (1,5,0,0) initial states. Only the rotational de-excitation cross sections are shown. The agreement between the 5D-CS and 6D-CC results is again very good for energies between 10-1000 cm−1. Apart from a few small discrepancies, the approximate 5D-CS formulation is able to reproduce the results of the 6D-CC calculation and extend them to higher energies. The agreement is not quite as good for the vibrationally inelastic cross sections shown in Figure 4. Cross sections for transitions from the same initial states as Figures 2 and 3 are shown for selected final states. The total cross section represents the sum over all final CO rotational level contributions in the v = 0 manifold. The state-to-state cross 1 sections show a larger number of discrepancies between the 5D-CS and 6D-CC calculations. However, these discrepancies aregenerally small, and the 5D-CSresults areable to provide a good approximation to the overall shape and relative magnitudes of the state-to-state cross sections for vibrationally inelastic collisions. The agreement between the two computational methods forthetotalcross sectionappearsto worsen athighenergies astheinitial rotational level of the CO molecule increases. This is due to the contributions from transitions to high rotational levels. For example, the (1,5,0,0) to (0,10,0,0) cross section in Figure 4(f) shows a relatively large discrepancy compared to the (0,0,0,0) and (0,2,0,0) cross sections which are in very good agreement at high energies. This pattern was found for all of the vibrationally inelastic cross sections and may be attributed to basis set truncation error in the 6D-CC calculations which used jmax = 22 and jmax = 20. These limits provide good v1=0 v1=1 convergence for vibrationally inelastic cross sections involving low-lying rotational levels at all energies. However, for transitions involving high rotational levels, the truncation limits are insufficient, particularly at high collision energies. It should benoted that the truncation in the 6D-CC calculations was a necessity due to the computational resource limitations. For comparison, the 5D-CS calculations in Figures 2-4 used jmax = jmax = 40. v1=0 v1=1 Thesize oftheCO vibrational basis set isconsidered inFigure5. Forboththe5D-CSand 6D-CC calculations, the vibrational basis set included v = 0−2. The vibrational quenching 1 cross section for the (2,0,0,0) initial state is dominated by the ∆v = −1 contribution, and 1 the comparison between the 5D-CS and 6D-CC results is similar to the (1,0,0,0) initial state shown in Figure 4. The ∆v = −2 contribution is at least two orders of magnitude 1 7 smaller than the ∆v = −1 contribution for both sets of calculations at the energies shown. 1 The agreement between the 5D-CS and 6D-CC results is not as good for the ∆v = −2 1 contribution above 300 cm−1, but given the relatively small magnitude of this cross section, the discrepancy is not very significant. Expanding the basis set to include closed vibrational levels yielded no change to these results. Figure 5b shows similar results for the (5,0,0,0) initial state. These cross sections are again dominated by the ∆v = −1 contribution, 1 however, the ∆v = −2 contribution increases relative to the ∆v = −1 contribution due 1 1 to the decreasing difference between vibrational energy levels. The small contribution from ∆v = −2 transitions allows a compact representation of the vibrational motion to be 1 obtained using basis sets with v equal to v and v−1, where v is the vibrational level of the 1 initial state. This enables large rotational basis sets to be employed (see Table I). The size of the para-H rotational basis set is considered in Figure 6. The solid curves 2 correspond to a rotational basis set which includes j = 0 − 4. In order to improve the 2 computational efficiency of these calculations, the truncation limit jmax of the initial CO v1=1 rotational basis was reduced to 20. The figure shows that the cross sections for vibrational relaxation of CO accompanied by rotational excitation of H are comparable in magnitude 2 to those with no change in H . Expanding the H basis further to allow more rotational 2 2 excitation yielded small contributions and provided little change to these results. The points shown in the figure correspond to the jmax = 40 basis set with H restricted to its ground v1 2 rovibrational state. The agreement in the H (0,0) cross sections is nearly perfect for the two 2 sets of calculations for energies below 2000 cm−1. At higher energies, it is not clear which basis set yields the more reliable result due to the different truncation limits on H and CO. 2 In order to explore rovibrational transitions for a large sample of internal states and collision energies, we constructed two separate databases for state-to-state cross sections and rate coefficients. The “frozen H ” database uses jmax = 40 and includes transitions 2 v for initial CO(j = 0−30) due to collisions with rovibrationally rigid para-H (j = 0) and 1 2 2 ortho-H (j = 1). The “flexible H ” database uses jmax = 20 and includes transitions for 2 2 2 v initial CO(j = 0−10) due to collisions with rotationally flexible para-H (j = 0,2,4) and 1 2 2 ortho-H (j = 1,3,5). Additional details of the basis sets are given in Table I. 2 2 Figures7-9showratecoefficientsforasampleofinitialstatesfromthefrozen-H database. 2 The pure rotational relaxation rate coefficients in Figure 7 are very similar to the rigid rotor results reported previously23. The present results are for v = 1, however, similar plots for 8 higher v are nearly identical, which confirms the validity of the rigid rotor approximation for pure rotational transitions of CO. For the rovibrationally inelastic transitions shown in Figures 8 and 9, the curves are well-ordered with respect to the initial vibrational level and generally exhibit a step-like structure with increasing temperature. The height of the steps increases with final state rotational level j′, and the rate coefficients are found to be 1 relatively large for j′ near the basis set truncation limit at the highest temperatures shown. 1 This indicates that the rovibrationally inelastic calculations are still not fully converged at the highest temperatures (T > 1000 K) even with the large truncation limit jmax = 40. v−1 Figure10showsresultsfromtheflexible-H databaseforthe(v,10,0,0)initialstate. Rate 2 coefficients for transitions involving no change in H are nearly identical to those presented 2 in Figure 8 and are not shown. Vibrational relaxation of CO(v = 1−5) for transitions which leave H in a rotationally excited state yield rate coefficients which are well-ordered with 2 respect to v as shown. The solid black and dashed red curves represent ∆j = 2 and ∆j = 4 2 2 transitions for H , respectively, with both sets of curves increasing with v. The step-like 2 structure is again seen in the rate coefficients and is a general feature for both databases. ThedatabasesincludeCOcollisionswithbothpara-H andortho-H . Figure11showsthe 2 2 ratio of total CO vibrational quenching rate coefficient for ortho-H and para-H collisions 2 2 for several initial states of CO. The ratio is generally found to approach unity for T > 100 K. As the initial rotational level of CO is increased, the ratio is close to unity for all vibrational levels and temperatures shown. These results suggest that para-H rate coefficients may be 2 directly used for ortho-H rate coefficients for T > 100 K with an error of no more than 10% 2 and a considerable savings in computational expense. IV. CONCLUSIONS Full-dimensional dynamics of CO interacting with H is a challenging computational 2 problem for quantum mechanical scattering formulations. The closely spaced rotational levelsforCOleadtosubstantialinternalangularmomentumcouplingbetweenthemolecules. Whenorbitalangularmomentumisproperlytakenintoaccountwithinthenumerically exact CC formulation, the additional angular momentum coupling creates a bottleneck for many rovibrationally excited initial states, which makes it impossible to bring the calculations to completion on a practical timescale. In such cases, it is necessary to use approximate 9 methods to obtain the desired data. The CS approximation offers a good compromise between accuracy andcomputational effort. Inthis approximation, the exact orbital angular momentum quantum number is replaced by an effective value, and the off-diagonal Coriolis coupling between states with different projection quantum numbers is neglected. Averaging the PES over the twist angle to obtain the 5D-CS approximation further improves the computational efficiency without introducing a significant additional lossof accuracy beyond the CS approximation itself. Using the5D-CSapproximation, we constructed two separate databases for state-to-state cross sections and rate coefficients for CO+H collisions. The “frozen-H ” database assumes 2 2 the collision takes place with the H molecule frozen in the lowest rotational level of a given 2 symmetry. This restriction allows the bulk of the computational effort to be used for the CO basis sets which include jmax = jmax = 40 for each initial v. The database then includes v v−1 transitions for each basis set with initial CO(v = 0 − 5,j = 0 − 30). The “flexible-H ” 2 database allows rotational excitation of H during the collision. To compensate for the 2 enlarged H basis sets, the CO basis sets are reduced to jmax = 20 with jmax = 40, and a 2 v v−1 smaller set of transitions from initial CO(v = 0−5,j = 0−10) is included in the database. The two databases are in good agreement with each other for transitions that are common to both basis sets. Both databases generally contain cross sections and rate coefficients for E = 1 − 10,000 cm−1 and T = 10 − 3,000 K. It would be desirable to improve the c databases by replacing the 5D-CS cross sections with 6D-CC values for the lower energies, and by expanding the 5D-CS basis sets to include additional CO rotational levels at the higher energies. These extensions will be implemented in future versions of the project. The present version provides a large amount of data which may be used as a starting point for further approximations (e.g. scaling) and provides important data for astrophysical models. The databases may be downloaded from the supplementary material associated with this article and from the UGA Excitation Database24 or PSU website25. ACKNOWLEDGMENTS This work was supported at Penn State by NSF Grant No. PHY-1503615, at UGA by NASA grant NNX12AF42G, and at UNLV by NSF Grant No. PHY-1505557. 10