Astrophysical modeling of processes in environments that are not in local thermal equilibrium requires the knowledge of state-to-state rate coefficients of rovibrational transitions in molecular collisions. These rate coefficients can be obtained from coupled-channel (CC) quantum scattering calculations, which are very demanding, however. Here, we present various approximate but more efficient methods based on the coupled-states approximation (CSA), which neglects the off-diagonal Coriolis coupling in the scattering Hamiltonian in body-fixed coordinates. In particular, we investigated a method called NNCC (nearest-neighbor Coriolis coupling) [Yang et al., J. Chem. Phys. 148, 084101 (2018)] that includes Coriolis coupling to first order. The NNCC method is more demanding than the common CSA method but still much more efficient than full CC calculations, and it is substantially more accurate than CSA. All of this is illustrated by showing state-to-state cross sections and rate coefficients of rovibrational transitions induced in CO2 by collisions with He atoms. It is also shown that a further reduction of CPU time, practically without loss of accuracy, can be obtained by combining the NNCC method with the multi-channel distorted-wave Born approximation that we applied in full CC calculations in a previous paper.

In modeling protoplanetary disks and other interstellar media that are not in local thermal equilibrium (LTE) with the aid of spectroscopic data from ground and satellite based telescopes,1–5 the effects of molecular collisions are important. In such non-LTE environments, the populations of the rovibrational states of a molecule—and thereby the characteristics of its spectrum—are not only determined by absorption and emission of electromagnetic radiation but also by transitions induced by molecular collisions. By analyzing these spectra, one gets crucial information not only about the abundance of various molecules but also about the local conditions. An essential element in this analysis is the knowledge of the rate coefficients of rovibrational transitions induced by the collisions of the molecule with H2 molecules, He atoms, and electrons. Quantum scattering calculations, based on intermolecular potentials determined by ab initio electronic-structure calculations, can provide inelastic collision cross sections, from which the required transition rate coefficients and their temperature dependence can be derived.

An important molecule in these studies is carbon dioxide, CO2.4–6 Accurate cross sections and rate coefficients for rotationally inelastic CO2–He collisions were recently reported by Godard Palluet et al.7 CO2 has no permanent dipole moment, so its rotational transitions are forbidden and it cannot be observed in microwave or far-infrared spectra. However, it can be observed in mid- and near-infrared spectra through its vibrational transitions. It has three vibrational modes: a twofold degenerate bend mode with an experimental frequency of 667 cm−1, an asymmetric stretch mode at 2349 cm−1, and a symmetric stretch mode at 1333 cm−1. The latter mode is not infrared active by itself but becomes observable through a Fermi resonance with the bend overtone. In the pioneering theoretical studies of rate coefficients for vibrational transitions in CO2 induced by collisions with rare gas (Rg) atoms by Clary et al.,8–11 they used VCC-infinite-order sudden (IOS), a vibrational coupled-channel (CC) method for the vibrations, combined with the infinite-order sudden (IOS) approximation for the rotations. This method provided rate coefficients for vibrational transitions, without considering specific initial and final rotational states. The more advanced models currently being developed by astronomers4,5 and the availability of data from the James Webb space telescope (JWST) in the near future require rovibrational state-to-state collisional rate coefficients. These can nowadays be obtained from the numerically exact coupled-channel (CC) method with the use of accurate ab initio calculated intermolecular potentials. Full CC calculations are still time-consuming, however, especially at the higher collision energies needed to obtain rate coefficients for higher temperatures. In a previous paper,12 we have shown how one can reach the CC level of accuracy with a less time-consuming procedure that handles the coupling between rotational states by the CC method and the weaker coupling between different vibrational states with the multichannel distorted-wave Born approximation (MC-DWBA). Here, we investigate further possibilities to speed up the calculation of collisional rate coefficients for rovibrational transitions by using the coupled-states approximation (CSA) and an improvement of it that includes Coriolis coupling to first order. We apply various methods to CO2–He collisions with CO2 excited in the symmetric stretch mode. We consider both the efficiency of these methods and the accuracy of the results they provide.

The methods discussed in the present paper are based on the coupled-channels (CC)—also called close-coupling— method, formulated in body-fixed (BF) coordinates. These coordinates refer to a BF frame with its z-axis along the vector R that points from the center of mass of CO2 to the He nucleus and the CO2–He complex lying in the xz-plane. The coordinates are the length R of the vector R, the angle θ between the CO2 axis and the vector R, and the normal coordinate Q along which CO2 is deformed with respect to its linear equilibrium geometry with equal C–O bond lengths of 1.162 Å.

The 3D Hamiltonian of the vibrotor-atom system CO2–He over a monomer normal coordinate Q is given, in the BF frame, by

Ĥ=22μR2R2R+ĤCO2(Q)+Ĵ2+ĵ22ĵĴ2μR2+V(Q,R,θ),
(1)

where μ=mCO2mHe/(mCO2+mHe) is the reduced mass of the complex, ĵ is the CO2 monomer rotational angular momentum operator, Ĵ is the total angular momentum operator of the complex, and Ĵ2+ĵ22ĵĴ represents the end-over-end angular momentum operator L2 in the BF-frame.13 The monomer Hamiltonian ĤCO2(Q) is defined in Eq. (2) of Ref. 12 and also the computation of the normal modes Q—its eigenstates in the rigid-rotor harmonic-oscillator approximation— is described there. Here, we consider the symmetric stretch mode QQ1 = 0.17678Δz1 − 0.17678Δz3 (a0), where Δz1 and Δz3 are the displacements of the O atoms along the CO2 axis with respect to their equilibrium positions. The symmetric stretch mode does not displace the C atom.

The eigenfunctions of the CO2 monomer Hamiltonian are

|vjΩ=χvj(Q)YjΩ(θ,ϕ),
(2)

and the corresponding eigenvalues are ϵvj. The vibrational functions χvj(Q) are j-dependent, and the rotational functions YjΩ(θ, ϕ) are spherical harmonics. The latter are expressed with respect to the BF frame, the angle θ is the same as defined above, and the angle ϕ coincides with the third Euler angle for the overall rotation of the complex. The quantum number Ω is the projection of the CO2 angular momentum ĵ and of the total angular momentum Ĵ on the BF z-axis along the vector R.

The channel basis in coupled channels (CC) scattering calculations for rovibrationally inelastic CO2–He collisions is, in BF coordinates,

|vjΩ;JMJ=2J+14πχvj(Q)YjΩ(θ,0)DMJΩJ(α,β,ϕ)*.
(3)

The angles (β, α) are the polar angles of the vector R with respect to a space-fixed frame (SF), and the Euler angles (α, β, ϕ) in the Wigner D-functions describe the orientation of the BF frame relative to the SF frame. The angle ϕ has been moved from the spherical harmonics in Eq. (2) to the overall rotation functions, which is mathematically equivalent.

The Coriolis coupling operator 2ĵĴ in Eq. (1) can be written as

2ĵĴ=2ĵzĴz+ĵ+Ĵ++ĵĴ.
(4)

The so-called helicity quantum number Ω is an eigenvalue of both ĵz and Ĵz. It is an approximate quantum number; basis functions with different Ω are mixed by the ladder operators ĵ±Ĵ± in Eq. (4), which couple functions with Ω to those with Ω ± 1.

The ab initio calculation of the 3D CO2–He potential with CO2 deformed along the symmetric stretch coordinate Q1 is described in Ref. 12. The well depth of this potential for CO2 at its equilibrium geometry is 47.43 cm−1. When this potential is expanded in Legendre polynomials Pλ(cos θ) of order λ, as in Ref. 12,

V(Q,R,θ)=λCλ(Q,R)Pλ(cosθ),
(5)

its matrix elements over the BF basis are

VvjΩ;vjΩ(R)=vjΩ;JMJ|V(Q,R,θ)|vjΩ;JMJ=δΩΩλ(1)Ω[(2j+1)(2j+1)]1/2×jλj000jλjΩ0Ω×vj(Q)|Cλ(Q,R)|vj(Q).
(6)

Equation (6) shows the advantages of the BF basis: the potential V(Q, R, θ) does not couple functions with different Ω, and the expression for its remaining matrix elements is simpler than in the SF basis,14 which makes the calculations more efficient.

The overall angular momentum J and its projection MJ on the SF z-axis are exact quantum numbers and also the overall parity P under inversion of the system is a conserved quantity. The basis in Eq. (3) is not invariant under inversion; a parity adapted basis is

|vjΩ̃;PJMJ=|vjΩ̃;PJMJ+P(1)J|vjΩ̃;PJMJ/×2(1+δΩ̃0),
(7)

where Ω̃0 and P = ±1 is the overall parity.

Another valid symmetry operation is the interchange P13 of the O atoms in CO2. This operator affects only the monomer wave functions in the basis of Eq. (3). For the symmetric stretch mode that we consider here, we find

P̂13|vjΩ̃=(1)j|vjΩ̃.
(8)

Since 16O nuclei are bosons with spin zero, the wave functions must be symmetric under P̂13. This implies that only functions with even j are allowed.

The scattering wave functions in CC calculations are written in terms of the parity-adapted BF channel basis as

ΨPJMJ=1RvjΩ̃|vjΩ̃;PJMJψvjΩ̃PJMJ(R).
(9)

When these functions are substituted into the time-independent Schrödinger equation, it follows that the radial wave functions ψvjΩ̃PJ(R) must obey a set of coupled second order differential equations, the CC equations

2R2ψvjΩ̃PJ(R)=vjΩ̃WvjΩ;vjΩPJ(R)ψvjΩ̃PJ(R),
(10)

or in matrix form

ψ(R)=W(R)ψ(R).
(11)

The column vector ψ(R) contains the radial wave functions ψvjΩ̃PJ(R). The quantum number MJ has been omitted, since the solutions do not depend on it. The elements of the matrix W over the primitive basis in Eq. (3) are given by

WvjΩ;vjΩJ(R)=δvvδjjδΩΩkvj2+TvjΩ;vjΩJ(R)+2μVvjΩ;vjΩ(R),
(12)

with

kvj2=2μ(Eϵvj)
(13)

and E being the total energy. The matrix elements of the potential are defined in Eq. (6); this matrix is diagonal in Ω. Only the matrix T that originates from the angular kinetic energy operator is not diagonal in Ω. Its elements are

TvjΩ;vjΩJ(R)=R2δvvδjjδΩΩJ(J+1)+j(j+1)2Ω2δΩΩ±1J(J+1)Ω(Ω±1)1/2×j(j+1)Ω(Ω±1)1/2.
(14)

Hence, T, and therefore also W, contains a series of blocks diagonal in Ω and a series of neighboring blocks with Ω′ = Ω ± 1. All other elements of these matrices are zero, thanks to the use of a BF basis in which the potential matrix V is diagonal in Ω.

We solve these equations with the renormalized Numerov propagator method.15,16 This method implies that one defines an equidistant grid Ri, i = 1, …, n and propagates the matrix Q, which defines the ratio of the radial wave functions in subsequent grid points Ri−1 and Ri,

ψ(Ri1)=Qiψ(Ri).
(15)

The propagation starts at small R1, where the potential is sufficiently repulsive that the wave function—and therefore Q—is zero and continues to large Rn, where the potential has vanished. Then, we assume that the radial wave functions obey flux-normalized K-matrix boundary conditions at large R,

ψ(R)=F(R)G(R)K.
(16)

The symbol ψ(R) is here used for a matrix with column vectors that are the solutions of the CC equations in Eq. (11). The blocks of the matrices F(R) and G(R) for the open channels are diagonal,

FvjL;vjL(R)=δvvδjjδLLkvj1/2RjL(kvjR),GvjL;vjL(R)=δvvδjjδLLkvj1/2RyL(kvjR),
(17)

and contain asymptotic wave functions that are proportional to spherical Riccati–Bessel functions17 of the first and second kind jL(z)=12π/zJL+12(z) and yL(z)=12π/zYL+12(z). Similarly, closed channels are matched by modified spherical Bessel functions of the third kind 12π/zKL+12(z) in the matrix G(R).

The asymptotic wave functions are defined in the SF frame and depend on the partial wave index L, while the matrix Qn is obtained from the propagation in BF coordinates. Therefore, this matrix is first transformed to SF coordinates,

QnSF=UQnU.
(18)

The elements of U,

UΩLJj=jΩL0|JΩ2L+12J+1,
(19)

contain Clebsch–Gordan coefficients ⟨⋯|⋯⟩.12 The matrix K can then be obtained from QnSF by solving the linear equations,

G(Rn1)QnSFG(Rn)K=F(Rn1)QnSFF(Rn).
(20)

Finally, we use the open-channel block Koo of the matrix K to obtain the scattering matrix

S=(IiKoo)1(I+iKoo),
(21)

with I being the unit matrix, and compute state-to-state scattering cross sections,

σv,jv,j(E)=π(2j+1)kvj2PJ(2J+1)×L=|Jj|J+jL=|Jj|J+jδvvδjjδLLSvjL;vjLPJ(E)2.
(22)

Expressions for the temperature (T) dependent rate coefficients kv′,j′←v,j(T) of the transitions from rovibrational state v, j to state v′, j′ are given in Ref. 12. Since the total angular momentum J and the overall parity P are exact quantum numbers, the CC equations can be solved separately for parities P = ±1 and for all values of J required to obtain converged cross sections.

In the coupled-states approximation (CSA), introduced long ago,18 one neglects the Coriolis coupling terms off-diagonal in Ω that appear in the second line of Eq. (14). This makes Ω an exact quantum number so that the CC equations can be separated into subsets of equations for each value of Ω. The dimension of each subset is smaller by a factor of min(2J + 1, 2jmax + 1), where jmax is the maximum j-value in the basis, than the dimension of the full CC equations. Since the CPU time to solve the coupled equations is proportional to the third power of their dimension, this yields a large reduction in computer time. Moreover, since different Ω values are not coupled, the largest absolute value of Ω is limited to the smallest of the initial or final j value in the scattering process, which is smaller than jmax—and also much smaller than J in most cases—so the number of equation subsets to be solved is small also.

In the standard application of the CSA, the whole angular kinetic operator in Eq. (4) is replaced by an operator L̂2, with eigenvalues Leff(Leff + 1). The possible values of L range from |Jj| to J + j, and different choices of Leff have been investigated. One mostly uses Leff = J, which yields an angular kinetic energy proportional to J(J + 1). In our BF implementation with the renormalized Numerov propagator, we have three different options. First, we can follow the original CSA algorithm by using Leff = J also in the matching of the asymptotic scattering wave functions to obtain the scattering matrices SPJΩ for all values of P, J, and Ω. The cross sections can then be obtained from the equation

σv,jv,j(E)=π(2j+1)kvj2PJΩ(2J+1)δvvδjjSvj;vjPJΩ(E)2.
(23)

Instead of using each of the matrices QnΩ obtained after solving the coupled-states equations for different Ω values directly in an asymptotic matching procedure that yields matrices SPJΩ, one may follow an alternative procedure. This alternative implies that the individual matrices QnΩ are collected into a total matrix Qn, which consists of diagonal subblocks with the matrices QnΩ for all Ω values. All other elements of Qn are zero, since in CSA there is no coupling between functions with different Ω. Since this total Qn matrix involves the full BF channel basis containing all values of Ω, it can be transformed to its SF equivalent in the same way as in the full CC treatment [see Eq. (18)]. It can then be used in the same asymptotic matching procedure as described in Sec. II A to obtain the scattering matrices SPJ that occur in Eq. (22) for the cross sections. In the propagation of the individual matrices QiΩ with i = 1, …, n, one may either use the full diagonal angular kinetic energy J(J+1)+j(j+1)2Ω2/2μR2 from the first line of Eq. (14) or an effective angular kinetic energy Leff(Leff + 1) = J(J + 1) as in the original CSA method. Altogether, this yields three different variants of the CSA, of which we compare the results in Sec. III. In Figs. 1 and 2, we label these variants with Leff = J for the standard CSA application and with the angular kinetic energies J(J+1)+j(j+1)2Ω2/2μR2 and J(J + 1) that we take into account in the latter two CSA methods.

The error in the CSA relative to full CC calculations is sometimes not acceptable. Recently, Yang et al.19 introduced an improvement of the CSA method, called NNCC, which includes the Coriolis coupling between functions with Ω and Ω ± 1 to first order. It implies that the CC equations are still solved separately for each Ω value but with the off-diagonal Coriolis couplings to the neighboring blocks with Ω ± 1 included in the matrix W [cf. Eqs. (12) and (14)]. In the full CC equations, these neighboring blocks are again coupled to functions with Ω ± 2, but the latter are not directly coupled to the functions with the given Ω, so they only give rise to second and higher order contributions to the solutions. The NNCC method neglects these higher order Coriolis couplings.

The basis to solve the CC equations for each Ω in the NNCC method involves also the bases for Ω − 1 and Ω + 1. The propagation can then be done with the full angular kinetic energy from Eq. (14). Functions with different Ω are coupled by the off-diagonal Coriolis coupling terms in the second line of this equation, which has the effect that Ω is not a good quantum number anymore. We could still label the separate propagations with Ω, but in order to distinguish the individual propagations from the quantum numbers, we label them with the index Ω̄. The propagation does not have to be performed for all Ω values because the Ω = 0 basis is already included in the calculation for the parity-adapted basis with Ω̃=1 [see Eq. (7)], and the Ω̃=J basis is included in the calculation for Ω̃=J1. In the NNCC method, as in the CSA method, the calculations can be further restricted to Ω̃ values limited by the smallest of the initial and final j quantum numbers, augmented by one in this case. At the end of the propagations, we have two options. The first one is that we follow the NNCC method as proposed by Yang et al.19 At the end of each propagation Ω̄, they construct an equivalent SF basis by diagonalizing the operator L̂2 with the BF matrix elements from Eq. (14) in the space of functions Ω − 1, Ω, and Ω + 1. The eigenvalues L(L + 1) from this diagonalization are non-integer, and also L is therefore non-integer. The same procedure has been implemented previously in the reactive scattering program ABC,20 which has the possibility to truncate the Ω basis. The matrix BΩ̄ with the eigenvectors of L̂2 in the restricted Ω space as column vectors is then used to transform the matrix QnΩ̄ obtained at the end of propagation Ω̄ to its SF equivalent,

QnΩ̄,SF=BΩ̄QnΩ̄BΩ̄.
(24)

This numerical diagonalization of L̂2 leaves the overall sign of each eigenvector, i.e., of each column of the matrix BΩ̄, undetermined. This has no effect on the final results, however, because the inverse of BΩ̄ is used in Eq. (25) to transform the S-matrix back from the SF to the BF frame. Indeed, we found numerically that the final cross sections in Eq. (26) do not depend on these signs. Since the SF basis thus obtained corresponds to non-integer L values, the spherical Bessel functions used in the matching procedure described in Eqs. (16)(20) must be the corresponding functions with non-integer L values.17 

The second option is that we consider the matrices QnΩ̄ over the BF basis with Ω − 1, Ω, and Ω + 1 as part of the matrix Qn over the full BF basis and transform them to the SF basis with the aid of Eq. (18). The transformation matrix U contains Clebsch–Gordan coefficients, the L values are integers, and the asymptotic matching procedure to obtain SF matrices SPJΩ̄(SF) from each propagation Ω̄ can be done with the usual spherical Bessel functions.

Finally, in both options, one has to transform the matrices SPJΩ̄(SF) in the SF basis back to the BF basis with the aid of the inverse BF to SF transformation. In the first option, this transformation is done with the matrix BΩ̄ for the given Ω̄. Since BΩ̄ is a real orthogonal matrix, its inverse is simply its transpose, and one obtains

SvjΩ;vjΩPJΩ̄(BF)=LLiLLBΩLΩ̄SvjL;vjLPJΩ̄(SF)BLΩΩ̄.
(25)

In the second option, the BF to SF transformation is done with the matrix U from Eq. (18), which is also orthogonal, and one obtains the same formula, with BΩ̄ replaced by U. In order to avoid double counting, we select from each matrix SPJΩ̄(BF) only the blocks with Ω=Ω̄ and calculate the cross sections with the formula

σv,jv,j(E)=π(2j+1)kvj2ΩΩ=Ω1Ω+1PJΩ(2J+1)×δvvδjjSvjΩ;vjΩPJΩ̄(BF)(E)2.
(26)

In a previous paper,12 we have shown how the application of a multi-channel distorted wave Born approximation (MC-DWBA) in scattering calculations reduces the computer time by about a factor of three with respect to exact CC calculations but produces cross sections and rate coefficients for rovibrational transitions in CO2–He collisions that are about equally accurate. The MC-DWBA algorithm and the reason why it can be favorably applied to rovibrational transitions are explained in detail in Ref. 12. Here, we investigate the application of this algorithm in combination with the NNCC method.

We study state-to-state (v, jv′, j′) rovibrational transitions in CO2 by collisions with He in which CO2 is de-excited from the v = 1 symmetric stretch fundamental to the v′ = 0 ground state. The details of the calculation and the characteristics of the 3D intermolecular potential of CO2–He with CO2 deformed along the symmetric stretch coordinate are given in Ref. 12. The basis in the scattering calculations consisted of CO2 symmetric stretch functions with v = 0 and 1, and the rotational basis contained all functions with j ≤ 70 for both v = 0 and 1. Calculations were made with a larger vibrational basis including also v = 2 functions, but this made practically no difference for the results. It is interesting that the states with v, j = 1, 0 and v, j = 0, 58 have nearly the same energy, so that transitions between these states are nearly resonant. The radial grid for the renormalized Numerov propagator contained 502 equidistant points in the range 3 ≤ R ≤ 30 a0. Rovibrational state-to-state cross sections were calculated for collision energies up to 3000 cm−1, with steps of 0.1 cm−1 in the resonance regime from 1 ≤ E ≤ 20 cm−1, steps of 1 cm−1 for 20 ≤ E ≤ 50 cm−1, 2 cm−1 for 50 ≤ E ≤ 100 cm−1, 50 cm−1 for 100 ≤ E ≤ 1000 cm−1, and 200 cm−1 for 1000 ≤ E ≤ 3000 cm−1. The largest total J value included was 100 for both parities P = ±1. The corresponding rate coefficients were calculated for temperatures from 10 to 500 K by cubic spline interpolating the cross sections over the energy grid and calculating the integral in Eq. (14) of Ref. 12 with the trapezoidal rule. Since the rotational states of CO2 up to j = 50 with energy 991 cm−1 are populated at the highest temperature, we calculated the cross sections and rate coefficients for initial states up to this value of j.

Figures 1 and 2 show the integral cross sections (ICSs) from the three different CSA methods described in Sec. II B, in comparison with ICSs from full CC calculations, for quenching from initial states v, j = 1, 0 and v, j = 1, 6 to different final states v′ = 0, j′. The meaning of the legends referring to the different CSA methods is explained in the caption of this figure. The sharp peaks in the ICSs for collision energies below 20 cm−1 correspond to resonances, which are extremely sensitive to the details of the calculations. Hence, it is not surprising that they are different for the different scattering methods used. However, also for higher energies, one observes that the ICSs from the CSA methods differ substantially from the full CC results. Two of the CSA methods produce very similar results. In both of them, the propagation is done in the BF frame with the angular kinetic energy term J(J + 1), but the asymptotic matching procedures to obtain the S-matrix are different (see Sec. II B). Apparently, the way of matching is less important than the approximation in the angular kinetic energy term used in the propagation. The third CSA method, with the angular kinetic energy term J(J + 1) + j(j + 1) − 2Ω2, produces rather different ICSs. The deviations of the CSA ICSs from the CC results are rather erratic, and one cannot conclude that one of the CSA methods is definitely better than the others. The method with the term J(J + 1) + j(j + 1) − 2Ω2 contains the full angular kinetic energy when neglecting the Coriolis coupling between basis functions with different Ω, so in the following, we will show the results from this particular CSA approximation.

FIG. 1.

ICSs from the three different CSA methods described in Sec. II B compared with full CC results. The upper legend CC refers to ICSs from full CC calculations. The next legend CSA Leff = J refers to the conventional CSA method with angular kinetic energy Leff(Leff + 1) and Leff = J in which the ICS is calculated directly from the S-matrices for different Ω in the BF frame. The lower two legends refer to the CSA methods with angular kinetic energy terms J(J + 1) + j(j + 1) − 2Ω2 or J(J + 1) used in the BF propagation and the ICS obtained from the S-matrix in the SF frame after transformation of the BF matrix Qn to the SF frame with the aid of Eq. (18). The initial state is v = 1, j = 0, and the final states are v′ = 0 with j′ = 0 (a), j′ = 2 (b), and j′ = 4 (c).

FIG. 1.

ICSs from the three different CSA methods described in Sec. II B compared with full CC results. The upper legend CC refers to ICSs from full CC calculations. The next legend CSA Leff = J refers to the conventional CSA method with angular kinetic energy Leff(Leff + 1) and Leff = J in which the ICS is calculated directly from the S-matrices for different Ω in the BF frame. The lower two legends refer to the CSA methods with angular kinetic energy terms J(J + 1) + j(j + 1) − 2Ω2 or J(J + 1) used in the BF propagation and the ICS obtained from the S-matrix in the SF frame after transformation of the BF matrix Qn to the SF frame with the aid of Eq. (18). The initial state is v = 1, j = 0, and the final states are v′ = 0 with j′ = 0 (a), j′ = 2 (b), and j′ = 4 (c).

Close modal
FIG. 2.

Same as Fig. 1, with a different initial state: v = 1, j = 6 and final states v′ = 0 with j′ = 2 (a), j′ = 6 (b), and j′ = 10 (c).

FIG. 2.

Same as Fig. 1, with a different initial state: v = 1, j = 6 and final states v′ = 0 with j′ = 2 (a), j′ = 6 (b), and j′ = 10 (c).

Close modal

Figure 3 illustrates the effect of the different asymptotic matching procedures used in the NNCC methods described in Sec. II C. In the method proposed by Yang et al.,19 they numerically determine the eigenvectors of the operator L̂2 in the BF basis restricted to a given Ω and Ω ± 1, which yields non-integer eigenvalues and, hence, non-integer L values. These non-integer L values are then used in the asymptotic matching procedure that yields the S-matrix. In the second method, we use specific rows of the matrix U from Eq. (19) that transforms the full BF basis to the SF basis. This full BF basis corresponds to a SF basis with integer L values, which are used in the asymptotic matching procedure. Figure 3 shows that the ICSs from the two methods are very similar for collision energies higher than about 10 cm−1. For lower energies, between 0.1 and 10 cm−1, they become different, but this is the resonance regime where the ICSs are very sensitive to details of the potential and the computational method. The difference becomes very large for still lower energies approaching the Wigner regime. In this regime, the ICSs of inelastic collisions should depend on E−1/2 according to the Wigner threshold laws.21 The method with integer L values nicely obeys this relation, but the method with non-integer L values fails completely. This is not surprising because the Wigner threshold law21 assumes pure s-wave scattering at the limit of very low energies, which implies that L = 0. Only the method with integer L correctly reaches this limit.

FIG. 3.

Comparison of ICSs calculated with the two different NNCC methods described in Sec. II C. They differ in the asymptotic matching procedure using integer L or non-integer L values.

FIG. 3.

Comparison of ICSs calculated with the two different NNCC methods described in Sec. II C. They differ in the asymptotic matching procedure using integer L or non-integer L values.

Close modal

Figures 4 and 5 show product distributions from the NNCC methods with integer and non-integer L values, as described in Sec. II C, for different initial states and different collision energies, compared with results from CSA and full CC calculations. A similar comparison is made in Figs. 6 and 7, which show the ICSs for the same initial v = 1, j values and specific final v′ = 0, j′ values as functions of the collision energy. It is obvious from all these figures that the NNCC method, which includes the first-order Coriolis coupling between basis functions with Ω and Ω ± 1, produces results in much better agreement with full CC calculations than the CSA methods. The NNCC methods with integer and non-integer L values perform quite similarly, except, of course, for collision energies below 0.1 cm−1 displayed in Fig. 3, where the NNCC method with integer L values becomes better because it obeys the Wigner threshold law, as discussed above.

FIG. 4.

Product (v′ = 0, j′) distributions from the NNCC methods with integer and non-integer L values for initial state v = 1, j = 0, compared with CSA and full CC calculations at collision energies E = 100 cm−1 (a) and 500 cm−1 (b). Note that only even values of j and j′ are physically allowed. The CSA method contains the angular kinetic energy J(J + 1) + j(j + 1) − 2Ω2.

FIG. 4.

Product (v′ = 0, j′) distributions from the NNCC methods with integer and non-integer L values for initial state v = 1, j = 0, compared with CSA and full CC calculations at collision energies E = 100 cm−1 (a) and 500 cm−1 (b). Note that only even values of j and j′ are physically allowed. The CSA method contains the angular kinetic energy J(J + 1) + j(j + 1) − 2Ω2.

Close modal
FIG. 5.

Same as Fig. 4 for the initial state v = 1, j = 6 at collision energies E = 100 cm−1 (a) and 500 cm−1 (b).

FIG. 5.

Same as Fig. 4 for the initial state v = 1, j = 6 at collision energies E = 100 cm−1 (a) and 500 cm−1 (b).

Close modal
FIG. 6.

ICSs from the NNCC methods with integer and non-integer L values described in Sec. II C compared with CSA and full CC results. The ICSs are shown as functions of the collision energy E for the initial state v = 1, j = 0 and final states v′ = 0, j′ = 0 (a), v′ = 0, j′ = 2 (b), and v′ = 0, j′ = 4 (c). The CSA method contains the angular kinetic energy J(J + 1) + j(j + 1) − 2Ω2.

FIG. 6.

ICSs from the NNCC methods with integer and non-integer L values described in Sec. II C compared with CSA and full CC results. The ICSs are shown as functions of the collision energy E for the initial state v = 1, j = 0 and final states v′ = 0, j′ = 0 (a), v′ = 0, j′ = 2 (b), and v′ = 0, j′ = 4 (c). The CSA method contains the angular kinetic energy J(J + 1) + j(j + 1) − 2Ω2.

Close modal
FIG. 7.

Same as Fig. 6 for the initial state v = 1, j = 6 and final states v′ = 0, j′ = 2 (a), v′ = 0, j′ = 6 (b), and v′ = 0, j′ = 10 (c).

FIG. 7.

Same as Fig. 6 for the initial state v = 1, j = 6 and final states v′ = 0, j′ = 2 (a), v′ = 0, j′ = 6 (b), and v′ = 0, j′ = 10 (c).

Close modal

In Figs. 8 and 9, we display the total quenching cross sections for transitions from initial states v = 1, j = 0 and v = 1, j = 6 to all final v′ = 0, j′ states. The resonance structure in the ICSs is clearly more pronounced for the v = 1, j = 0 initial state than for the v = 1, j = 6 state. These resonances are so sensitive to the method of computation that none of the methods can reproduce the CC results. For energies higher than 20 cm−1 above the resonance regime both NNCC methods produce results in good agreement with full CC results for both initial states v = 1, j = 0 and v = 1, j = 6. It seems somewhat surprising that also the CSA method yields fairly good results for the v = 1, j = 0 initial state (Fig. 8), but as one can see in Fig. 4, CSA considerably overestimates the ICSs for some of the final j′ states and underestimates them for other j′ values, and these errors nearly compensate each other. For the initial v = 1, j = 6 state (Fig. 9), the CSA results are clearly inferior to the ICSs from the NNCC methods for energies below 200 cm−1.

FIG. 8.

ICSs for quenching from the initial state v = 1, j = 0 to the v′ = 0 state, summed over all final j′ values. The different curves are results from the NNCC methods with integer and non-integer L values, compared with CSA and full CC results. The CSA method contains the angular kinetic energy J(J + 1) + j(j + 1) − 2Ω2.

FIG. 8.

ICSs for quenching from the initial state v = 1, j = 0 to the v′ = 0 state, summed over all final j′ values. The different curves are results from the NNCC methods with integer and non-integer L values, compared with CSA and full CC results. The CSA method contains the angular kinetic energy J(J + 1) + j(j + 1) − 2Ω2.

Close modal
FIG. 9.

Same as Fig. 8 for the initial state v = 1, j = 6.

FIG. 9.

Same as Fig. 8 for the initial state v = 1, j = 6.

Close modal

A remarkable observation when comparing Figs. 8 and 9 is that the total quenching cross sections are almost the same for the v = 1, j = 0 and v = 1, j = 6 initial states. We also considered other initial j values, and we found, as in our studies in Ref. 12, that these total quenching cross sections hardly depend on the initial j value. This is important as it implies that one can apply our results for some specific initial j values more generally to all different initial j states and, thus, obtain total v = 1 → v′ = 0 quenching cross sections and rate coefficients as functions of the temperature.

So far, we discussed the ICSs from different methods as functions of the collision energy, but in astrophysical modeling, one needs collisional rate coefficients as functions of the temperature. Figures 10 and 11 show state-to-state rate coefficients calculated from the ICSs displayed in Figs. 6 and 7. It is clear that the conclusions discussed above for the ICSs apply also to the rate coefficients. That is, the NNCC method is clearly superior to the CSA method in reproducing the full CC results.

FIG. 10.

State-to-state rate coefficients as functions of the temperature calculated for the initial state v = 1, j = 0 and final states v′ = 0, j′ = 0 (a), v′ = 0, j′ = 2 (b), and v′ = 0, j′ = 4 (c), with the ICSs shown in Figs. 6(a)6(c).

FIG. 10.

State-to-state rate coefficients as functions of the temperature calculated for the initial state v = 1, j = 0 and final states v′ = 0, j′ = 0 (a), v′ = 0, j′ = 2 (b), and v′ = 0, j′ = 4 (c), with the ICSs shown in Figs. 6(a)6(c).

Close modal
FIG. 11.

State-to-state rate coefficients as functions of the temperature calculated for the initial state v = 1, j = 6 and final states v′ = 0, j′ = 2 (a), v′ = 0, j′ = 6 (b), and v′ = 0, j′ = 10 (c), with the ICSs shown in Figs. 7(a)7(c).

FIG. 11.

State-to-state rate coefficients as functions of the temperature calculated for the initial state v = 1, j = 6 and final states v′ = 0, j′ = 2 (a), v′ = 0, j′ = 6 (b), and v′ = 0, j′ = 10 (c), with the ICSs shown in Figs. 7(a)7(c).

Close modal

Finally, Fig. 12 demonstrates that state-to-state ICSs for rovibrationally inelastic collisions from the NNCC method are accurately reproduced by combining this method with MC-DWBA. This is very useful, as it shows that MC-DWBA cannot only be applied to the full CC method, as in Ref. 12, but also to more approximate and less expensive methods.

FIG. 12.

Product (v′ = 0, j′) distributions from the NNCC method with non-integer L values, compared with results from MC-DWBA applied to this method. Panels (a) and (b) are for the initial state v = 1, j = 0 at collision energies E = 100 and 500 cm−1, respectively, and panels (c) and (d) are for the initial state v = 1, j = 6 at the same collision energies.

FIG. 12.

Product (v′ = 0, j′) distributions from the NNCC method with non-integer L values, compared with results from MC-DWBA applied to this method. Panels (a) and (b) are for the initial state v = 1, j = 0 at collision energies E = 100 and 500 cm−1, respectively, and panels (c) and (d) are for the initial state v = 1, j = 6 at the same collision energies.

Close modal

Another issue to be discussed is that how much more efficient are the CSA and NNCC methods than the full CC method, both in terms of computer memory and CPU time. We checked this in calculations with a basis of CO2 rotational functions up to j = 100 for v = 0 and 1 and total angular momentum J = 100. As the initial state we took v = 1, j = 6, so in the CSA calculations, we needed to solve seven sets of coupled-channel equations for Ω̃ in the symmetry-adapted basis ranging from 0 to 6. The NNCC method includes blocks with Ω and Ω ± 1, and we had to make calculations for Ω̃ = [0, 1, 2], [1, 2, 3], [2, 3, 4], [3, 4, 5], [4, 5, 6], and [5, 6, 7]. The number of channels in the CC calculations was 5202; in CSA, it was 100; and in NNCC, it was 290. The size of the matrices involved in solving the coupled-channel problem depends quadratically on the number of channels, so it is clear that the CSA and NNCC methods allow one to handle much larger problems than the full CC method. Even more significant is the CPU time for the propagation, which depends roughly on the third power of the number of channels. Our renormalized Numerov propagation involved 501 steps and took 6633 CPU seconds for the CC method, 0.40 seconds for CSA, and 4.15 seconds for NNCC on a single Intel Xeon Platinum 8268 processor. So even with the more advanced NNCC approximation, the savings in CPU time relative to full CC calculations are substantial. We note here that the actual gain in CPU time and matrix size in CSA and NNCC with respect to full CC calculations depends on the characteristics of the channel basis. In our calculations of the cross sections in rovibrationally inelastic CO2–He collisions, we had to use a channel basis with large maximum values of the monomer rotational angular momentum j and of the total angular momentum J.

In Ref. 12, we explained that the MC-DWBA method applied in full CC calculations reduces the CPU time by about a factor of 3. Here, we find that also when combined with the NNCC method MC-DWBA leads to a further decrease in CPU time by a factor of 3.

The advanced models currently being developed by astronomers4,5 and the availability of data from the James Webb space telescope (JWST) in the near future require the knowledge of rovibrational state-to-state collisional rate coefficients. These rate coefficients can be obtained from coupled-channel (CC) scattering calculations, but these are very demanding. Here, we presented more efficient methods based on the coupled-states approximation (CSA) in which one neglects the off-diagonal Coriolis coupling in the scattering Hamiltonian in body-fixed coordinates. This makes Ω, the projection of the total angular momentum J on the intermolecular axis, a good quantum number so that scattering calculations can be performed independently for each Ω. In addition to CSA, we investigated a method called NNCC (nearest-neighbor Coriolis coupling)19 that includes Coriolis coupling to first order by simultaneously including basis functions with Ω and Ω ± 1. The NNCC method is more expensive than the CSA method, but still much more efficient than full CC calculations. We tested three versions of the CSA method and two versions of the NNCC method. The cross sections and rate coefficients from the two NNCC methods are similar and substantially better than all CSA results.

All of this is illustrated by showing state-to-state cross sections and rate coefficients of rovibrational transitions induced in CO2 by collisions with He atoms. Results from the CSA and NNCC methods are compared in detail with those from full CC calculations. In a recent paper,12 we have shown that the application of the multi-channel distorted-wave Born approximation (MC-DWBA) in CC calculations reduces the required CPU time by a factor of 3. Here, we show that a further increase of efficiency by about the same factor can be obtained by applying MC-DWBA to the NNCC method, with practically no loss of accuracy.

Finally, we note that rovibrational transitions in CO2 are probably more strongly induced by collisions with H2 than by collisions with He for several reasons. First, CO2–H2 interactions are stronger than CO2–He interactions because H2 is more polarizable than He and it has a quadrupole moment. Second, as was recently shown for rovibrationally inelastic H2O–H2 collisions,22 such transitions may be enhanced by simultaneous rotational excitation of H2. The channel basis needed in CO2–H2 scattering calculations will be even larger than the bases we used for CO2–He, and the efficient methods to compute rovibrationally inelastic collision cross sections and rate coefficients presented here will be very advantageous.

We thank Ewine van Dishoeck and Arthur Bosman for stimulating and useful discussions. This work was supported by The Netherlands Organisation for Scientific Research, NWO, through the Dutch Astrochemistry Network DAN-II. We also acknowledge useful discussions with David Manolopoulos.

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article.

1.
K. M.
Pontoppidan
,
G. A.
Blake
,
E. F.
van Dishoeck
,
A.
Smette
,
M. J.
Ireland
, and
J.
Brown
,
Astrophys. J.
684
,
1323
(
2008
).
2.
A. M.
Mandell
,
J.
Bast
,
E. F.
van Dishoeck
,
G. A.
Blake
,
C.
Salyk
,
M. J.
Mumma
, and
G.
Villanueva
,
Astrophys. J.
747
,
92
(
2012
).
3.
S.
Bruderer
,
D.
Harsono
, and
E. F.
van Dishoeck
,
Astron. Astrophys.
575
,
A94
(
2015
).
4.
A. D.
Bosman
,
S.
Bruderer
, and
E. F.
van Dishoeck
,
Astron. Astrophys.
601
,
A36
(
2017
).
5.
A. D.
Bosman
, “
Uncovering the ingredients for planet formation
,” Ph.D. thesis,
Leiden University
,
2019
.
6.
K. I.
Öberg
,
A. C. A.
Boogert
,
K. M.
Pontoppidan
,
S.
van den Broek
,
E. F.
van Dishoeck
,
S.
Bottinelli
,
G. A.
Blake
, and
N. J.
Evans
,
Astrophys. J.
740
,
109
(
2011
).
7.
A.
Godard Palluet
,
F.
Thibault
, and
F.
Lique
,
J. Chem. Phys.
156
,
104303
(
2022
).
8.
D. C.
Clary
,
J. Chem. Phys.
75
,
209
(
1981
).
10.
A. J.
Banks
and
D. C.
Clary
,
J. Chem. Phys.
86
,
802
(
1987
).
11.
C. T.
Wickham-Jones
,
C. J. S. M.
Simpson
, and
D. C.
Clary
,
Chem. Phys.
117
,
9
(
1987
).
12.
T.
Selim
,
A.
Christianen
,
A.
van der Avoird
, and
G. C.
Groenenboom
,
J. Chem. Phys.
155
,
034105
(
2021
).
13.
A.
van der Avoird
,
P. E. S.
Wormer
, and
R.
Moszynski
,
Chem. Rev.
94
,
1931
(
1994
).
14.
A. M.
Arthurs
and
A.
Dalgarno
,
Proc. R. Soc. London, Ser. A
256
,
540
(
1960
).
15.
B. R.
Johnson
,
J. Chem. Phys.
69
,
4678
(
1978
).
16.
B. R.
Johnson
, in
NRCC Proceedings
(
Lawrence Berkeley Laboratory, Berkeley, CA
,
1979
), Vol. 5, p.
86
. https://www.osti.gov/biblio/5603893.
17.
M.
Abramowitz
and
I. A.
Stegun
,
Handbook of Mathematical Functions
(
National Bureau of Standards
,
Washington, DC
,
1964
).
18.
P.
McGuire
and
D. J.
Kouri
,
J. Chem. Phys.
60
,
2488
(
1974
).
19.
D.
Yang
,
X.
Hu
,
D. H.
Zhang
, and
D.
Xie
,
J. Chem. Phys.
148
,
084101
(
2018
).
20.
D.
Skouteris
,
J. F.
Castillo
, and
D. E.
Manolopoulos
,
Comput. Phys. Commun.
133
,
128
(
2000
).
21.
E. P.
Wigner
,
Phys. Rev.
73
,
1002
(
1948
).
22.
L.
Wiesenfeld
,
J. Chem. Phys.
155
,
071104
(
2021
).