Rotational Excitations in CO–CO Collisions at Low Temperature: Time

Article pubs.acs.org/JPCA

Rotational Excitations in CO−CO Collisions at Low Temperature: Time-Independent and Multiconfigurational Time-Dependent Hartree Calculations Steve A. Ndengué,† Richard Dawes,*,† and Fabien Gatti‡ †

Department of Chemistry, Missouri University of Science and Technology, 142 Schrenk Hall, 400 West 11th Street, Rolla, Missouri 65409, United States ‡ CTMM, Institut Charles Gerhardt, UMR 5253, Univeristé de Montpellier II, Place Eugène Bataillon, 34095 Montpellier, France S Supporting Information *

ABSTRACT: The rotational excitation in collisions between two carbon monoxide molecules was studied while combining the use of both time-independent and time-dependent formalisms. All of the calculations made use of a recently published four dimensional PES for CO dimer. Timeindependent scattering calculations were performed in the lower part of the collision energy range, thus limiting the number of open channels and computational cost. The PES features a low-energy path for geared motion that appears to affect the excitation propensities in low-energy collisions. For reactants colliding without initial rotational excitation, symmetric excitations (both monomers excited equally) are strongly favored. This behavior deviates significantly from an exponential gap model based on endo- or exoergicity. Comparable timedependent calculations were performed in an extended energy range made feasible by the lower cost of those calculations. The wave packet propagation in the time-dependent approach was performed with the multiconfiguration time-dependent hartree (MCTDH) method and analyses via the flux method, and the Tannor and Weeks approach was used to calculate the transition probabilities in the energy range up to 1000 cm−1. We deduce from the cross sections the corresponding reaction rates for temperatures between 10 and 250 K. MCTDH was found to yield well-converged results, where the methods overlap, validating the use of MCTDH as an efficient tool to study collision processes.

■

INTRODUCTION Carbon monoxide is the second most abundant molecule in the universe after molecular hydrogen. Its importance has been demonstrated in processes occurring in the atmosphere, hydrocarbon combustion, and interstellar space. In astrophysics, the CO emission spectrum is frequently used to detect H2 and other molecules that do not have a spectrum visible from Earth. Usually, molecular level populations are not in equilibrium, and the prediction of spectral line intensities depends on the magnitude of collisional excitation rate coefficients with dominant species (H2, He, H and also CO). The scattering of CO by various colliders (H,1−5 H2,6−10 He,11−14 Ne,15−18 Ar19−23) has been intensively studied experimentally and theoretically over the past few decades. Despite the wide range of studies that have been performed on these other CO−X systems, to our knowledge, no rigorous quantum-mechanical calculations have been reported to date on the nonreactive scattering of CO−CO even with simplifications such as the coupled states (CS) or infinite order sudden (IOS) approximations despite the availability of accurate potential energy surfaces (PESs)24−26 for more than a decade. The lack of nonreactive scattering predictions for CO dimer can be explained by the difficulty of the calculations caused by the “heavy” masses © XXXX American Chemical Society

of the collision partners (small rotational constant), therefore requiring the use of very large basis sets to converge the quantum-mechanical calculations. This requirement is further stressed by the heteronuclear nature of the diatom, preventing the use of higher symmetry, which allows one to reduce the size of the basis, as has been done for other low-temperature inelastic scattering calculations on “heavy” homonuclear dimers O227,28 and N2.29,30 The most common approach to accurately describe rotational (de)excitations in nonreactive collisions of dimers is the timeindependent approach within the close coupling (CC) formalism. Although for collisions involving light−light or light−heavy partner these calculations can be done almost routinely with publicly available codes, heavy−heavy collisions pose a more significant problem because they typically become much more computationally expensive. Time-dependent approaches have Special Issue: 100 Years of Combustion Kinetics at Argonne: A Festschrift for Lawrence B. Harding, Joe V. Michael, and Albert F. Wagner Received: January 31, 2015 Revised: April 7, 2015

A

DOI: 10.1021/acs.jpca.5b01022 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 1. Plot of the 4D intermolecular PES for planar geometries using the extended angles described in ref 26. V(R = ∞) = 0. At each point in the plot, the energy is minimized with respect to the intermonomer distance R. A low-energy channel along a disrotatory path (corresponding to geared motion) connects the global minimum nonpolar isomer (NC, E = −135.1 cm−1) through a polar minimum (P1, E = −121.1 cm−1) to the other nonpolar isomer (NO, E = −119.6 cm−1). The path continues through P2 arriving back at NC. Images of the NC and NO isomers are centered about the corresponding minima on the plot. The extended angles are defined: θ̃1 = θ1 and θ̃2 = θ2 (quadrant I); θ̃1 = −θ1 and θ̃2 = θ2 (quadrant II); θ̃1 = −θ1 and θ̃2 = −θ2 (quadrant III); θ̃1 = θ1 and θ̃2 = −θ2 (quadrant IV).

dimensionality41,42 on the molecular hydrogen dimer benchmark system. More recently, Malenda et al.43 showed on a “heavy− heavy” atom−diatom collision system that the MCTDH approach is able, even at low energies, to produce cross sections of similar accuracy to numerically exact benchmarks obtained using time-independent CC calculations. In the same work, it was shown that the MCTDH calculations provide increased computational efficiency for larger systems over other timedependent approaches. Those results support the growing acceptance of MCTDH as the method of choice to treat quantum dynamics of high-dimensional systems and led us to consider its application for the titled diatom−diatom nonreactive collision process. The aim of this paper is thus to present quantum-mechanical results of the CO−CO nonreactive collision relevant for astrophysics applications that complement previous results on CO−X complexes and to discuss related questions such as the propensities of specific rovibrational transitions observed in nearly homonuclear diatoms. The work will explore the effectiveness of using the MCTDH method for the study of

been less commonly applied to study nonreactive scattering. Some examples are reported in the literature but mainly for light−light or heavy−light collisions.31−37 The two approaches differ in how the results are obtained. Using a time-dependent approach, for each set of collider initial states and each total angular momentum, the partial reaction probability is obtained in a broad range of collision energies. Time-dependent approaches are also known to be more easily applicable to higher ranges of collision energies. In contrast, time-independent approaches produce accurate cross sections (in the limit of the completeness of the basis) for particular well-defined collision energies and have been successfully applied for very low collision energies. The choice of method is therefore determined by the quantity of interest, and, in general, time-independent methods are most often used to describe very low collision energies and timedependent approaches for higher collision energies and reactive processes. The MCTDH method38−40 is an approach to solve the timedependent Schrödinger equation and has been used in the past decade to study inelastic scattering in full and reduced B

DOI: 10.1021/acs.jpca.5b01022 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

above the ZPE level. The rovibrational states are naturally sorted into stacks45,46 (sets of levels with similar rotational constants), and close correspondence between calculated and observed levels has permitted unambiguous assignment of numerous levels. The CBS_PL PES used in this study produces calculated levels in very close agreement with experiments. For 68 levels with J ≤ 6, the rms error is only 0.29 cm−1. The symmetry labels and many details of the states in various stacks have been discussed previously. Of interest here, anticipated to be reflected in the scattering cross sections and possibly even the state-tostate propensities, is the very low frequency of the intermolecular modes, particularly the geared bend mode (3.73 cm−1). Because of the low-frequency modes and multiple isomers the density of states is very high.

heavy diatom−diatom collisions and how these results could be confirmed and complemented by time-independent calculations at low energy. We will discuss aspects of the recently developed 4D PES used in our calculations and the nature of the bound states of the CO dimer. Next, the MCTDH and timeindependent calculations presented in this work are described. The last sections of this work present Results and Discussion, including comparisons between time-independent (with reduced basis) and time-dependent MCTDH calculations.

■

INTERMOLECULAR POTENTIAL ENERGY SURFACE The 4D (rigid monomer) PES used here was developed recently for a computational study of the intermolecular rovibrational bound states.26 In that study, several ab initio PESs were made using an automated algorithm and the method of interpolating moving least-squares (IMLS). The final PESs are composed of 2226 automatically determined data points, sufficient to reach an rms fitting error below 0.1 cm−1. This level of fidelity to the ab initio methods made it possible to explore the accuracy of various coupled-cluster-based schemes. As discussed in that paper (and references within), it is very challenging to obtain an ab initio PES of spectroscopic accuracy for the (CO)2 system due to the unusual importance of high-order correlation. Remarkably, the perturbative (T) triples correction increases the binding energy by >30%. Here we only use the best PES resulting from the previous study denoted CBS_PL, which uses data from explicitly correlated F12 coupled-cluster theory (with all-electrons correlated) at the (AE)-CCSD(T)-F12b/CVnZ-F12 (n = 3,4) levels extrapolated to the complete basis limit (CBS) using an optimized power law. The CBS_PL PES is available for download as supporting information to that paper. The global minimum on the PES is at −135.1 cm−1 relative to separate monomers and corresponds to a nonpolar slippedparallel isomer with the C atoms facing in labeled NC on the PES plot in Figure 1. The center-of-mass separation for the NC isomer is Re = 4.33 Å. (More details and geometric parameters of the isomers are found in ref 26.) One of the most notable features of the PES for planar geometries is a low-energy path for geared motion. A disrotatory path with very low barriers passes through a polar isomer (labeled P1) connecting to the other nonpolar isomer (O atoms facing in, labeled NO, E = −119.6 cm−1) and back again. Tunneling is important between NC and NO, and even the lowest energy bound states have wave functions that are quite delocalized along the geared-motion disrotatory path. The NO isomer has a significantly shorter center-of-mass distance of Re = 3.63 Å. The two lowest rovibrational bound states can be assigned to the NC and NO isomers (and have significantly different rotational constants reflecting the different values of Re). Although the two minima differ energetically by 15.6 cm−1 on the PES, due to the shape of the PES and tunneling, the calculated levels are only 1.2 cm −1 apart (0.9 cm − 1 experimentally).44 The lowest bound states are delocalized along the disrotatory path. (See figure 6 of ref 26 for a plot of the probability densities.) Interestingly, the ZPE level (assigned as the NC isomer) is at −80.5 cm−1, well above many features of the PES including the entire disrotatory path, yet retains the signature character of the well. Similarly, the next level (only ∼1 cm−1 above the ZPE) assigned to the NO isomer is also delocalized and affected by tunneling and also exhibits character imparted by the NO local minimum feature of the PES. Figure 7 of ref 26 shows that a state corresponding to the polar isomer was also identified localized about the P1 minimum but significantly higher (17.6 cm−1)

■

TIME-DEPENDENT CALCULATIONS WITH MCTDH The multiconfiguration time-dependent hartree 3 8−4 0 (MCTDH) is an algorithm to solve the time-dependent Schrödinger equation, which can be considered as a timedependent version of the complete active-space self-consistent field (CASSCF). Within this method, the wave function Ψ(Q,t) of the system is written as a sum of products of single-particle functions (SPFs), forming a time-dependent orthonormal basis set. SPFs are low-dimensional functions: when they contain more than one degree of freedom (DOF), the combined coordinates Qκ = q1,κ, ..., qdκ,κ that comprise dκ physical DOF are introduced. The ansatz of the MCTDH wave function reads Ψ(q1 ,..., qf , t ) ≡ Ψ(Q 1 ,..., Q p , t ) n1

=

j1

=

np

∑ ··· ∑ A j ,..., j (t )Πκp= 1φ j(κ)(Q κ , t ) p

1

jp

κ

∑ AJ ΦJ (1)

J

where f denotes the number of degrees of freedom and p is the number of MCTDH particles (combined modes). The AJ ≡ Aj1,···,jp denotes the MCTDH expansion coefficients, and the configuration or Hartree products ΦJ are products of SPFs defined in relation 1. The SPFs are finally represented by linear combinations of time-independent primitive basis functions N1, κ

φ j(κ)(Q κ , t ) = κ

N1, κ

∑ ··· ∑ c (j κl)··· l (t )χl(κ) (q1, κ )··· χl(κ) (qd ,κ ) κ 1

l1= 1

d

1

d

l1= 1

(2) 47,48

usually within a discrete variable representation (DVR), here χli(κ)(qi,κ) with the time-dependent coefficients cjκl1··· ld. Propagation of the Wave Packet with the MCTDH Method. Applying the Dirac−Frenkel variational principle on the ansatz (eq 1) leads to equations of motion for the expansion coefficients as well as the SPF. Evolving the wavepacket at each time step requires the calculation of multidimensional integrals over all DOF. It is therefore helpful to express the highdimensional terms within the Hamiltonian as a sum of products of low-dimensional ones. When the Hamiltonian is of product form, matrix elements can be expressed by a sum of products of monomode integrals. While the kinetic energy operator (KEO) often satisfies this requirement (especially when using the socalled polyspherical coordinates as in the present work), PESs are C

DOI: 10.1021/acs.jpca.5b01022 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A often given as complicated, nonanalytic, nonseparable expressions. Then, the PES cannot be used directly but instead needs to be converted to a product form. In the next section, we will provide more details of the sum of products representation of the CO dimer PES. The MCTDH method differs from conventional propagation methods in the fact that, in MCTDH, the wave function is expressed in an intermediate time-dependent basis set that is optimized by a variational principle: generally n1 × ... × np is significantly smaller than N1 × ... × Nf. The MCTDH equations of motion conserve the norm and the total energy (for timeindependent Hamiltonians), which follows directly from the variational principle. However, when a complex absorbing potential is added to the Hamiltonian, the norm of the wave packet is reduced as the outgoing wavepacket is absorbed by the wall. Kinetic Energy Operator and the Potential Energy Surface Representation. The treatment of the CO+CO collision requires six coordinates for a full-dimensional treatment. For technical reasons, it is advantageous to incorporate one external degree of freedom, the overall rotation about the intermolecular axis, into the Hamiltonian for the internal motion. We therefore begin with a 7D Hamiltonian, before reducing it to 5D with the rigid rotor approximation. We use in this work, for practical reasons related to the MCTDH implementation, a Jacobi coordinates parametrization in the E2 frame49,50 (see Figure 2), as used in previous studies on the H2+H2 collision41,42

where ji ̂ = −

ki2 ∂ 1 ∂ + sin θi ∂θi sin θi ∂θi sin 2 θi

(4)

ji ±̂ = ±

∂ − ki cot θi ∂θi

(5)

2

(with additional lowering and raising shifts of ki, ki → ki ± 1) and C±(J , K ) =

J(J + 1) − K (K ± 1)

(6)

Here J is the quantum number associated with the total angular momentum quantum number and K is the quantum number associated with the projection of the total angular momentum on the internuclear axis. jî , i = 1,2 are the angular momenta of the monomers and k1, k2 are the projections of the monomers angular momenta on the BF frame. μR is the reduced mass of the system. The rigid rotor approximation has been applied to the system and consists of removal from the 7D KEO terms involving the ∂2/∂r2i and replacing 1/2 μir2i by the CO groundstate rotational constant Brot = 1.931281 cm−1. As previously mentioned, to be optimal for the MCTDH calculations, the potential energy operator should be expressed in the sum of products form. The PES used in this work,26 however, is not represented in that particular form, at least in the E2 coordinates frame where the wave packet propagation will be performed. A convenient way to express the PES in a form readily usable by the MCTDH code is through the use of the Potfit algorithm proposed by Jäckle and Meyer.51,52 The PES in the rigid rotor approximation is given in four coordinates, the vibrational coordinates of the diatoms being frozen at their equilibrium value V = V (R , θ1 , θ2 , φ1BF)

(7)

The procedure used for the representation of the potential and repeated below follows essentially the same steps as performed in previous non reactive collisions studies41,42 with MCTDH in the E2 frame. The PES expressed in the four coordinates is fitted with the Potfit procedure after replacing the φBF 1 by a momentum variable Ω

Figure 2. Coordinates description of the four atom system in the E2 frame. θi and φEi 2, i = 1,2, are the polar and azimuthal of the monomer which orient those with respect to the E2 plane. The Z axis of the E2 frame is the same as the Z axis of the BF frame.

VΩ̃ (R , θ1, θ2) = (2π )−1

∫0

2π

BF

e−iΩφ1 V (R , θ1, θ2 , φ1BF) dφ1BF (8)

or rather than the familiar body-fixed frame. Note that the coordinates shown in Figure 2 (and used to perform the dynamics) are different than those from Dawes et al. used to plot the potential in Figure 1. The KEO in these coordinates has been described previously41 and takes the form

V (R , θ1 , θ2 , φ1BF) =

Ω

=

− −

− 2k1k 2) −

C −(J , k1 + k 2) μR R2

μR R2

(j1̂ − + j2̂ − )

(j ̂

1+

E2

e−iΩφ2

(9)

⎛ ⎛ 1 ∂ 1 ⎞ 2̂ 1 ⎞ ⎟j + ⎜⎜2Brot + ⎟⎟ + ⎜⎜2Brot + 2T̂ = − 2 2⎟1 μR ∂R μR R ⎠ μR R2 ⎠ ⎝ ⎝ 1 2 j2̂ + (j ̂ j ̂ + j1̂ − j2̂ + + J(J + 1) − 2k12 μR R2 1 + 2 − C+(J , k1 + k 2)

E2 1

∑ VΩ̃ (R , θ1, θ2) × eiΩφ Ω

E2 E2 because φBF 1 =φ1 − φ2 . The action of the potential on the wave function then writes

2

2k 22

BF 1

∑ VΩ̃ (R , θ1, θ2) × eiΩφ

(V̂ Ψ)(R , θ1, k1, θ2 , k 2) =

∑ VΩ̃ (R , θ1, k1, θ2 , k 2) × Ψ(R , θ1, k1 − Ω, θ2 , k 2 + Ω) Ω

(10)

+ j2̂ + )

In the present work the potential is real and symmetric and it suffices to only perform calculations with Ω ≥ 0; here Ω = 0,···,4. A very generous set of parameters for the primitive basis set as well as the SPF basis (reported in Table 1) was used for the

(3) D

DOI: 10.1021/acs.jpca.5b01022 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A with

Table 1. Parameters of the Primitive Basis Used for the Nonreactive Scattering Calculationsa

ζj , m(θ , k) =

R

(θ1,φ1), (θ2,φ2)

primitive basis number of points/functions

FFT 768

range size of spf basis

5.0−60.0 a0 30−75

KLeg-DVR j = 0, 1, ···, 25 for θ and K = −25, −24, ···, 25 for φ 0 to π for θ, 0 to 2π for φ 30−75

FFT representation and KLeg-DVR stands, respectively, for fast Fourier transform and extended Legendre DVR.53

Pfi(E) =

present calculation, although the calculations could be performed confidently with a smaller set. This choice was made to ensure convergence of our calculations but also because of the high density of states of the system under study, which leads to higher correlations compared with the H2+H2 case. The Dawes CO−CO PES is provided as a Fortran routine in four coordinates R, θ1, θ2, and ϕ (φBF 1 ). We first perform a Fourier transform of the potential, projecting it according to relation 9. Then, the Ṽ Ω(R,θ1,θ2), Ω = 0, ..., 4 are fitted as sum-of-products using the Potfit program, which comes with the Heidelberg MCTDH package. The quality of the projection as well as the fit rely on the size of the primitive basis used. Table 1 presents the primitive basis used for the dynamics calculations. A fast Fourier transform representation is used for the R coordinate: the extended-Legendre DVR (K-DVR),53 which is a combined θ,φDVR, for the polar and azimuthal angle of the diatoms. For the MCTDH calculations, the fitted potential has to be expressed in the same primitive basis as for the dynamics calculations, for each relevant coordinate. The primitive basis presented in Table 1 gives the mean error in the projection of the potential of 0.00187 meV (0.015 cm−1). The errors in the sum-of-products form fit of the Ṽ Ω(R,θ1,θ2) are then, respectively, for Ω = 0, ..., 4, 0.132, 0.087, 0.075, 0.043, and 0.032 meV.

1 | 2 (2π ) |Δf (E)|2 |Δi (E)|2

∫0

∞

dt eiEt ⟨Ψf |Ψ(t )⟩|2 (16)

where |Δ(E)| are the energy distributions of the initial wave function and the reference wave function, Ψ(t) is the propagated wave packet, and Ψf is the reference wave function: a direct product of an internal eigenstate of the diatom after the scattering and an outgoing Gaussian. The final internal eigenstate and Gaussian have the same form as equations 14 and 13, respectively, except that the momentum of the Gaussian has a positive value for the outgoing Gaussian instead of the negative value for the incoming Gaussian. The total angular-momentumdependent energy distribution is given by |Δ(E)|2 = ⟨χi |δ(Ĥ i − E)|χi ⟩

(17)

where Ĥ i = ⟨ψi|Ĥ |ψi⟩ is a mean field Hamiltonian and the total angular momentum dependence of the Hamiltonian is hidden to simplify the notation. In the flux formalism, the transition probability is given by P fiJ = |SfiJ(E)|2 =

■

Ffi(E) |Δi (E)|2

(18)

where the projected quantum flux is obtained from the flux operator57 F̂̂= i[Ĥ ,Θ]as

ANALYSIS OF THE WAVE PACKET PROPAGATION Transition Probabilities. The first step to obtain the cross section is to determine the transition probabilities, Pf i, which gives the probability that a molecule starting in the state i has after collision the quantum numbers f. It is related to the S-matrix by

̂ f̂ δ(E − Ĥ )|Ψ⟩ Ffi(E) = 2π ⟨Ψ|i δ(E − Ĥ )Pf̂ FP i

(19)

The term Θ in the flux operator expression is a characteristic function in the scattering coordinate, R, having value 1 if the scattering coordinate was greater than some radius Rc and zero otherwise. The working equation in the flux approach thus reads

(11)

To determine the S-matrix elements in the present work, we used both the Tannor and Weeks42,54 approach based on timecorrelation functions and the flux method39,55,56 based on analysis of the flux operator57 for comparison purposes. The initial wavepacket Ψi for either method is a product of an internal eigenstate ψi of the dimer and an incoming Gaussian χi in R Ψi(R , θ1 , k1 , θ2 , k 2) = χi (R )ψi(θ1 , k1 , θ2 , k 2)

(15)

In the above relations, the vibrational part of the dimer has been removed as we work in the rigid rotor approximation, and consider our diatoms to be in their ground vibrational state. The justification of eq 14 can be found in Appendix B of Gatti et al. In the Tannor and Weeks formulation, the energy-dependent transition probabibilities are given by

a

Pfi(E) = |Sfi(E) − δfi|2

2j + 1 (j − m)! m P j (cos θ )δkm 2 (j + m )!

P fiJ(E) =

2 π |Δ(E)|2

∫0

T

[g W (τ ) + g fiΘ(τ )]eiEτ dτ fi

(20)

where gW (τ ) = fi

(12)

∫0

T−τ

⟨Ψi(t )|Pf̂ WPf̂ |Ψi(t + τ )⟩ dt

(21)

and

with χi (R ) =

⎡ ⎛ R − R ⎞2 ⎤ 1 0⎟ + ip0 (R − R 0)⎥ exp⎢ −⎜ 2πd ⎣ ⎝ 2d ⎠ ⎦

g fiΘ(τ ) =

ψi(θ1 , k1 , θ2 , k 2) = (2 + 2δ j , j δm1, −m2) 1 2

× [ζ j , m1(θ1 , k1)ζ j , m2(θ2 , k 2) + ( −1) J ζ j , −m2(θ1 , k1) ζ j , −m1(θ2 , k 2)] 1

2

(22)

where T is the final propagation time and P̂f is a projection onto the final rotational states. The characteristic function in the flux formalism serves as a dividing surface that separates the internal region from the asymptotic region. A complex absorbing potential (CAP), written as −iW(R), is used to annihilate the wavepacket that crossed the dividing surface and thus helps to reduce the grid size and prevent unwanted reflection. The W function is a non-

(13)

and

1

1 ⟨Ψi(T − τ )|Pf̂ ΘPf̂ |Ψi(T )⟩ 2

2

(14) E

DOI: 10.1021/acs.jpca.5b01022 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A negative real function that vanishes in the internal region and increases smoothly in the asymptotic region. Cross Sections and Rate Coefficients. After computing the transition probabilities for all relevant values of the total angular momentum in the energy range of interest, the energyresolved cross sections can be expressed in terms of the probabilities as σfi(E) =

π 2μR Ecoll

Table 2. CO Monomer Rotational Energies

∞

∑ [1 + (−1)J δ j j δm1, −m2] 12

J=0

× [1 + (− 1) J δ j′j′δm1′, −m2′](2J + 1)P fiJ(E) 12

(23)

that we can rewrite j1

j2

∑ ∑

σ(j1′j2′ ← j1 j2 ; E) =

m1 =−j1 m2 = 0

(2 − δ0, m2)

with π 2μR Etrans

j1′

j2′

Jmax

∑ ∑

∑

×

m1′=−j1′ m2′=−j2′ J ′=|m1′+ m2′|

(2J + 1)[1 + ( −1) J δ j , j δm1, −m2] × P j′m1′, j′ m2′← j m1, j m2(E) 1 2

1

2

1

2

(25)

The cross sections for rotational collisional excitation, starting from the ground state, are then simply written as π (1 + δ j′j′)

σ(j1′j2′ ← 00; E) = j1′

j2′

∑

∑

m1′=−j1′ m2′=−j2 v

μR Ecoll

∞

∑ (2J + 1) J=0

P j′m1′, j′ m2′← 00,00(E) 1

2

(26)

The thermal rate coefficients for the state-to-state transitions are obtained by averaging the cross sections over the Boltzmann distribution of the collision energy k fi(T ) =

8β 4 πμR β

∫0

∞

j

ε (cm−1)

0 1 2 3 4 5 6 7 8 9 10 11 12

0.0 3.8 11.5 22.9 38.2 57.3 80.2 106.9 137.4 171.8 209.9 251.9 297.7

13 14 15 16 17 18 19 20 21 22 23 24 25

347.3 400.7 457.9 518.9 583.7 652.2 724.6 800.8 880.8 964.5 1052.0 1143.4 1238.4

Alexander and Manolopoulos.61 The propagation is carried out with the diabatic modified log-derivative method from a minimum distance of 4.5 bohrs to an intermediate one of 33 bohrs and with the Airy method up to a maximum intermolecular distance R of 60 bohrs. The 4D intermolecular (CO)2 potential was expanded in θ and ϕ in terms of Legendre polynomials up to the fifth order. The range of the PES was extended from its original 42 bohrs to 60 bohrs to converge the calculations. The rotational parameters62,63 and the reduced mass of CO used in this study are Be = 1.93128087 cm−1, De = 6.12 × 10−6 cm−1, αe = 1.7504 × 10−2 cm−1, and μ = 13.9974575 AMU. The first few (rigid) rotational levels64 of carbon monoxide monomer obtained from these parameters are shown in Table 2. They can be used to determine the endo- or exoergicity of particular state-to-state processes and the number of energetically accessible channels for each collision energy. The rotational basis functions for the time-independent calculation were selected using j1max = j2max = 7. For these basis parameters, at ET = 100 cm−1 there are 36 levels with 23 open, 13 closed, and none missing. At ET = 200 cm−1, 35 channels are open and 1 is closed, with 10 missing open channels. The maximum total angular momentum was set to Jtot = 100 for energies below 100 cm−1 and increased up to Jtot = 150 at the end of our energy range, in this case 150 cm−1. The main issue with these calculations is that the computational resources, that is, the computational time as well as the memory usage, quickly becomes prohibitive. A calculation using a smaller basis with j1max = j2max = 6 at ET = 100 cm−1 typically takes 2 days to run for all of the J values and parities. The same type of calculation with j1max = j2max = 7 takes typically 3 weeks, and for j1max = j2max = 8 it takes a few months. The number of channels basis for j1max = j2max = 6, 7, and 8, is, respectively 588, 1008, and 1620. It is therefore not currently feasible for us to extend the study to higher energies using this method (or even perform as many convergence tests as might be desired). To test the convergence of our time-independent calculations below 100 cm−1, we performed calculation of inelastic scattering cross sections for j1max = j2max = 6 and j1max = j2max = 7 for the same maximal Jtot = 100. The results for some transitions are plotted in Figure 1 of the Supporting Information. We also tested calculations of the cross sections with j1max = j2max = 7 comparing maximal Jtot = 100 and Jtot = 120: results are shown in Figure 2 of Supporting Information. The results for other transitions show similar behavior. This supports the idea that the time-independent calculations at low energy correctly capture the dynamics of the system. The change in the cross sections

(24)

12

ε (cm−1)

(2j1 + 1)(2j2 + 1)

× σm1m2(j1′j2′ ← j1 j2 ; E)

σm1m2(j1′j2′ ← j1 j2 ; E) =

j

dEcoll e−βEcollEcollσfi(E) (27)

where β = 1/kT with k being the Boltzmann constant.

■

TIME-INDEPENDENT CALCULATIONS: CROSS SECTIONS AND COLLISIONAL RATES The quantum-mechanical close-coupling approach of Arthurs and Dalgarno58 was used to calculate the cross sections in carbon monoxide collisions. The theory describing the collision of two 1 Σ diatoms59 is well known and is therefore not detailed explicitly here. The reaction rates of the scattering are then obtained from the cross sections from relation eq 27. The details of the numerical calculations and convergence tests will be given in the next section.

■

RESULTS AND DISCUSSION Numerical Details. The time-independent CC calculations reported in this work were done with the general purpose scattering program MOLSCAT.60 The propagation was carried out with the hybrid modified log-derivative Airy propagator of F

DOI: 10.1021/acs.jpca.5b01022 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 3. Time-independent state resolved CO+CO inelastic cross sections for various j′1,j′2 ← 0,0 transitions.

cm−1. For this energy range, the relevant maximum total angular momentum necessary to achieve converged results was found to be Jmax = 200. We performed propagations for J = 0, 1, 2, 3, 4, 6, 8, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, and 200. The partial probabilities for the other values were produced with the J-interpolation algorithm, an improved version of the J-shifting approximation,65−67 presented in appendix A of Gatti et al. The parameters of the primitive basis set and grid used for the present calculations are listed in Table 1. A dense primitive basis in a large grid was used for the present calculations to minimize numerical convergence errors. Several tests were performed to select the various parameters of the calculations. The primitive grid, that is, the integration range and number of functions, was analyzed for increasing number of functions and spatial range. The initial position of the wavepacket was chosen so that it does not “feel” the effect of the attractive potential until long after its start. The convergence of the propagation was confirmed by inspecting selected reaction probabilities with increasing basis size (number of single particle functions). Partial Cross Sections and Reaction Rates. The energy dependence of the state-resolved cross sections for j1′ ,j2′ ← 0,0 obtained from the time independent calculations are depicted in Figure 3. The state-resolved cross sections for other j′1,j′2 ← j1,j2 transitions are given in the Supporting Information. The calculations presented in these plots focus on collision energies below 150 cm−1. Because the rotational basis was limited to j1max

Figure 4. CO+CO inelastic cross sections for the 1,2 ← 0,0 transition using the Tannor and Weeks, flux, and time-independent methods.

going from j1max = j2max = 6 to j1max = j2max = 7 is on average