Accurate transport properties for O(³P)–H and O(³P)–H₂

The computational modeling of combustion in practical devices such as internal combustion engines requires both a set of rate constants describing the chemistry of the fuel-oxidant mixture and also transport properties to represent the flow of heat and chemical species in the spatially inhomogeneous combustion environment. Progress on computational simulation of laminar flames has recently been reviewed by Smooke.¹ 

The conventional approach for the estimation of transport properties utilizes isotropic Lennard-Jones (LJ) (12-6) potentials, with the well depth and length parameters obtained through combination rules for unlike pairs.^2–4 Calculations of transport properties through classical trajectory calculations employing potential energy surfaces (PES’s) with full dimensionality of the nuclear coordinates have recently been carried out for several systems.^5–8 There have also been a number of time-independent quantum scattering (QS) calculations of transport properties that have employed ab initio PES’s. In some cases,^9–15 these calculations have employed angular averaging to generate an isotropic potential, to simplify the scattering calculations.

In recent years, our group has been engaged in the calculation of transport properties on collision pairs, mainly involving the light H atom.^16–21 In these studies, we have carried out time-independent quantum scattering calculations on state-of-the-art ab initio PES’s. In some cases, we have compared the transport properties with those obtained by employing the spherical average of the PES’s, to test the importance of the anisotropy of the PES’s for those systems.^16,18 For some of the collision pairs [e.g., H–H (Ref. 10), H–H₂ (Ref. 22), H–O₂ (Ref. 19), H–CO and H–CO₂ (Ref. 21)], the accurately calculated transport properties are significantly different in magnitude and have a different temperature dependence than those computed with assumed LJ (12-6) potentials.

In order to test whether the differences between the computed quantum scattering (QS) and conventional LJ transport properties are significant for combustion modeling, we have compared 1-dimensional combustion simulations on several flames using these two sets of transport properties.^20,21 In particular, we found for 1 atm CH₄/air flames that when the QS transport properties were substituted for the LJ values, the differences in the computed laminar flame speeds were significantly greater than when different chemistry models were employed with the same set of transport properties.²⁰ 

In the present work, we extend QS calculations of transport properties to two collision pairs involving the H atom and oxygen, namely, O(³P)–H(²S) and O(³P)–H₂. Yarkony and Parlant^23,24 have carried out multi-configuration self-consistent-field (MCSCF) calculations of the potential energy curves (PEC’s) emanating from the O(³P)–H(²S) asymptote and the spin-orbit matrix elements coupling these states. They used these data to compute predissociation rates for OH(A²Σ⁺) rovibrational levels. Krems and co-workers^25,26 used the results of the calculations by Yarkony and Parlant^23,24 to compute rate constants for fine-structure excitation of O(³P_j) and ³P → ¹D electronic excitation in collisions with H atoms. We have carried out new, internally contracted multi-reference configuration-interaction calculations with single- and double-excitation (MRCISD). In contrast to the previous MCSCF calculations,^23,24 our calculations show shallow van der Waals wells on the excited-state curves.

There have been a number of calculations of PES’s for the interaction of O(³P) with H₂. Rogers et al.²⁷ reported chemically accurate PES’s for the lowest ³A′ and ³A″ states of the O(³P) + H₂ → OH + H reaction. Attention has been paid to characterizing the van der Waals well of the entrance channel^28–32 since the analytic fit of the PES’s of Rogers et al. poorly describes this region of the PES’s.³¹ We report here explicitly correlated internally contracted multi-reference interaction calculations for the three states emanating from the O(³P)–H₂ asymptote.

This paper is organized as follows: we present a brief review of the calculation of transport properties in Sec. II. The calculation of the O–H and O–H₂ interaction energies is presented in Secs. III A and III B. The scattering calculations and computed transport properties are described in Secs. IV and V, respectively. The paper concludes with a discussion in Sec. VI.

We present in this section a brief review of the calculation of transport properties of molecule immersed dilutely in a bath of a collision partner. Transport properties can be computed from collision integrals Ω^(n,s)(T), which can be obtained by taking an appropriate integral over the collision energy and carrying out a Boltzmann state average, namely,^33–36 

In Eq. (1), ε_i is the energy of the ith O(³P) fine-structure level (and energy of the H₂ rotational level for O–H₂), q_R is the partition function, and k_B is the Boltzmann constant. For O(³P)–H₂, the degeneracy g_i includes the nuclear degeneracies associated with para- and ortho-H₂. The state-specific collision integrals are integrals of state-specific transport cross sections $Qi(n)(E)$ over the collision energy E,

Here, μ is the collision reduced mass.

The state-specific transport cross section in Eq. (2) is a sum over final levels f of state-to-state transport cross sections,

The state-to-state transport cross sections on the right-hand side of Eq. (3) are weighted angle averages of the i → f differential cross section,

where $Rˆ=(θ,ϕ)$ is the orientation of the Jacobi vector R. The transport cross sections $Qi→f(n)(E)$ can be expressed as a weighted sum over several low-order Legendre moments of the differential cross sections. The pertinent weighting factors in Eq. (4) for n = 1 and 2 have been presented previously.^17,34,35 The two systems under study here involve collision pairs for which both partners have internal structure. The  Appendix presents an appropriate expression for the transport cross section in this case.

Reduced collision integrals Ω^(n,s)∗ are usually employed in the calculation of transport properties. These are related in the following way to the collision integrals in Eq. (1):^4,33

where the reduced temperature T^∗ = k_BT/ϵ, ϵ and σ are the LJ well depth and length parameters, respectively, and the factor F(n, s) equals

The binary diffusion coefficient D is related to the (1,1) reduced collision integral,

where N is the total number density of the gas. The following ratios of collision integrals are required for the calculation of thermal diffusion coefficients:³⁷ 

The OH PEC’s were calculated using the multi-reference configuration interaction method with single and double excitation and the Davidson correction (MRCISD+Q), implemented in the MOLPRO suite of programs.³⁸ We used the augmented, correlation-consistent, quadruple zeta basis of Dunning, Kendall, and co-workers.^39,40 In the MRCISD+Q calculations, we took reference wave functions obtained from state-averaged complete active space self-consistent field (SA-CASSCF) calculations. We used MOLPRO’s built-in projector to obtain states with specified values of the quantum number Λ, which is the molecule-frame projection of the electronic orbital angular momentum. The calculations were carried out on a grid of R values over the range 1.1–23.5a₀.

For completeness, we computed PEC’s for the eight states emanating from the O(³P, ¹D, ¹S) + H(²S) asymptotes. A table of the interaction energies is presented in the supplementary material. The left-hand panel of Fig. 1 displays the PEC’s for these electronic states. The three excited states correlating with the ground asymptote [1⁴Σ⁻, 1²Σ⁻, 1⁴Π] possess shallow van der Waals wells, as shown in the expanded plot of interaction energies in the right-hand panel of Fig. 1. We took the R-dependent spin-orbit matrix elements from the work of Yarkony and Parlant,²⁴ checking that the phases of the matrix elements were consistent.²⁶ In the scattering calculations, the electronic energies and the spin-orbit matrix elements at a given value of R were determined by spline fits of the grid of values as a function of R.

Table I presents the well depths D_e and equilibrium internuclear separations R_e for all the calculated OH states. The computed dissociation energy D₀ is in reasonable agreement with the experimental value [35 580 cm⁻¹ (Ref. 41)]. The supplementary material presents calculated spectroscopic parameters for the two strongly bound OH electronic states [X²Π, A²Σ⁺]. The rotational constants can be seen there to be in good agreement with experimental values.

Warehime⁴² has computed 3-dimensional PES’s depending on the Jacobi coordinates (r, R, θ) to describe the interaction of O(³P) with H₂. Since we are interested in elastic and rotationally inelastic collisions only, we consider here the PES’s for a fixed H₂ internuclear separation r = 1.40a₀, near r_e(H₂) = 1.402a₀.³⁰ We define the O–H₂ complex to lie in the yz plane and the H–H bond to lie along the z axis.

The PES’s for the one ³A′ and two ³A″ states emanating from the ground O(³P) + H₂ asymptote were computed using the explicitly correlated variant^43–45 of the MRCISD+Q method^46–48 with single and double excitations and with the addition of the Davidson correction.⁴⁹ We used Dunning’s augmented correlation-consistent triple-zeta basis (aug-cc-pvtz)³⁹ coupled with the JKFIT/MP2FIT density fitting basis. The interaction energies were calculated by subtracting from the total energy at a given geometry the asymptotic total energy of O(³P)–H₂ in linear geometry (θ = 0°) at R = 40a₀. The calculations were carried out with the MOLPRO suite of programs.³⁸ The reference orbitals for this procedure were taken to be the orbitals in collinear geometry. The calculations were carried out on a grid of angles from θ = 0° to 90° in steps of 5° and over R = 2–12a₀.

We also require the mixing angle γ between the two states of A″ symmetry, in order to transform to a diabatic basis to facilitate the scattering calculations. In previous calculations on other systems,^50–54 the mixing angle was determined by a method, developed originally by Rebentrost and Lester,⁵⁵ based on the off-diagonal matrix element of the angular momentum operator between the two states of the same symmetry. Diabatization can be carried out with one or more molecular properties.⁵⁶ In the present work, we employed the spin-orbit operator to determine the mixing angle.

Table II presents calculated dissociation energies D_e and equilibrium separations R_e for the O(³P)–H₂ states in linear and perpendicular geometries, for which there is no 1³A″–2³A″ mixing. Atahan et al.³² reviewed calculations of O(³P)–H₂ van der Waals wells. We compare our calculations with MRCI calculations by Alexander³⁰ and the coupled cluster [RCCSD(T)] calculations of Atahan.³² We see in Table II that there is generally good agreement of the dissociation energies and equilibrium separations between the calculations. The most noticeable difference is that the dissociation energies for linear geometry are significantly larger than in the earlier calculations. No counterpoise correction was applied in the present calculation.

In principle, O(³P) + H₂ can react to form OH + H products in a slightly endothermic process. However, since we have fixed the H₂ bond length, we do not see any evidence toward reaction in our calculations of the interaction potential.

For the scattering calculations, we require diabatic PES’s. These were obtained from the adiabatic energies [E(1A′), E(1A″), E(2A″)] and the mixing angle γ through the procedure described previously.³⁰ We briefly review this procedure. The matrix of the interaction in the diabatic Cartesian basis can be written as

The angular dependence of these diabatic matrix elements can be expanded in terms of reduced rotation matrices,

Since the molecule is homonuclear, only even-λ_r terms are nonvanishing in the expansions of V_zz, V_s, V_d, and V_yz.

Alexander and Yang^30,57 relate these Cartesian expansion coefficients to the body-fixed (BF) V_λrλaμ(R) expansion coefficients defined by Dubernet and Hutson,⁵⁸ 

We have employed Eqs. (19)–(22) to compute the V_λrλaμ(R) expansion coefficients. With our angular fit, we required ten BF expansion coefficients: V₀₀₀, V₂₀₀, V₄₀₀, V₀₂₀, V₂₂₀, V₄₂₀, V₂₂₁, V₄₂₁, V₂₂₂, and V₄₂₂. Figure 2 presents these coefficients as a function of R. Also included in Fig. 2 are the coefficients from the potential computed by Alexander.³⁰ As can be seen, the counterpoise correction has little effect on the coefficients describing the potential. The anisotropy of the potential with respect to rotation of the H₂ molecule is described by the V_λr00 terms λ_r≠0; these terms are small for R ≳ 4.5a₀. The largest anisotropic terms are V₀₂₀, V₂₂₁, and V₂₂₀; the first term represents the variation of the interaction energy with respect to the orientation of the O(³P) orbital holes, while the latter describe the coupling between the orbital occupations and the H₂ rotation.³⁰ 

We employ a space-fixed (SF) scattering basis in our time-independent close-coupling scattering calculations. Dubernet and Hutson⁵⁸ showed how to expand the interaction potential in a SF basis. The relationship between the SF expansion coefficients, denoted as V^{λ_rλ_aλ₁₂}(R), and the BF expansion coefficients in Eqs. (19)–(22) is

where $...|...$ is a Clebsch-Gordan coefficient.⁵⁹ We employ the coefficients defined in Eq. (23) and the asymptotic atomic spin-orbit energies in our scattering calculations for O(³P)–H₂.

The quantum description of collisions between atoms in ³P and ²S states was first derived by Krems et al.²⁵ and is an extension of treatments of previous work of fine-structure-changing transitions and collision-induced electronic transitions in atoms.^60–64 Our scattering calculations employ a SF basis. For the collision of O(³P_jO) with H(²S_jH), the SF basis functions can be written as

where L is the orbital angular momentum, J, M, and Ω are the total angular momentum and its SF and BF projections, respectively, [x] = 2x + 1, and (:::) is a 3j symbol.⁵⁹ We also have j_H = 1/2 and j_O = 0, 1, 2; the angular momentum j is the vector sum of j_O and   j_H.

Matrix elements of the Hamiltonian between the BF scattering basis functions $jOjHjΩJM$ have a simple form,

To evaluate the matrix elements of H in the SF scattering basis, we use Eq. (24) to obtain

The electrostatic potentials and spin-orbit matrix elements, presented in Sec. III A, have been computed in the case (a) $ΛS ΣΩ$ basis, and we need to transform these to the $jOjHjΩ$ basis. We employ Eq. (II.10b) from Singer, Freed, and Band:⁶⁵ 

In Eq. (27), ⋮⋮⋮ is a 9j symbol.⁵⁹ In the $jOjHjΩ$ basis, H_SO is diagonal in the asymptotic, atomic limit. This was verified numerically with the large-R matrix elements of Parlant and Yarkony,²⁴ as corrected by Krems et al.²⁵ 

We see in Fig. 1 that the repulsive PEC’s correlating with the ground O(³P) + H(²S) asymptote cross the attractive OH(A²Σ⁺) PES. These crossings are responsible for predissociation of this electronic state.^23,24 Since the lowest of these crossings are ∼6000 cm⁻¹ above the energy of O(³P) + H(²S), we have ignored these crossings in our calculations.

The formalism for carrying out quantum scattering calculations for the collision of an atom in P electronic state with a diatomic molecule has been given previously.^30,57 We find it convenient to work in a SF scattering basis, and hence to use a SF expansion of the interaction potential [Eq. (23)]. Explicit expressions for the matrix elements of the interaction potential in a SF scattering basis have been given by Dubernet and Hutson [see Eq. (17) of Ref. 58].

We have carried out calculations employing the present potential and also the potential computed by Alexander.³⁰ This was done in order to see how the differences in the potentials, highlighted in Table II, affect the computed transport properties. The results reported in Sec. V were obtained with the present potential. Diffusion coefficients calculated with the two potentials are compared in Sec. VI.

Time-independent close-coupling calculations were carried out with the HIBRIDON suite of programs.⁶⁶ This program suite already contained a module to handle collisions of an atom in a ³P state with a homonuclear diatomic molecule. HIBRIDON was extended in this work to treat collisions between atoms in ³P and ²S states. The calculation of transport cross sections [Eq. (4)] was also extended to treat collisions of two species, both with internal structure, as described in the  Appendix. The cross sections were checked for convergence with respect to the inclusion of a sufficient number of partial waves and energetically closed channels. For O(³P)–H₂, the H₂ rotational basis included all levels with j ≤ 6.

The high-energy tails of the integrals over collision energy [Eq. (2)] were evaluated by exponential extrapolation of the integrand. For O(³P)–H₂ collisions, the sum over H₂ rotational levels in Eq. (1) included j = 0 and 2 for para-H₂ and j = 1 and 3 for ortho-H₂.

State-specific transport cross sections $Qi(1)(E)$ as a function of collision energy for O(³P)–H(²S) and O(³P)–H₂ are presented in Figs. 3(a) and 3(b). We see that, in general, the cross sections decrease monotonically with increasing collision energy. The major exception is for O(³P_j=2)–H(²S) collisions; in this case there is structure due to resonances from quasibound levels of OH(X²Π). Some structure is also observed in the energy-dependent O(³P_j=1)–H(²S) transport cross sections. For O(³P)–H(²S), the cross section for the j = 2 fine-structure level is somewhat larger than the j = 1 and 0 levels at intermediate collision energies [≲1000 cm⁻¹]. However, at higher collision energies, cross sections for the latter levels become slightly larger than for j = 2.

For O(³P)–H₂ collisions, there are also some differences in the magnitude of the cross sections for the different H₂ rotational levels associated with a given O-atom fine-structure level for O(³P)–H₂ collisions. Transport cross sections for collisions of O(³P_j=2) with the H₂ j = 0 and 1 rotational levels have similar magnitude. This contrasts with collisions involving O(³P_j=1,0), for which the cross sections for collision with H₂ j = 1 are significantly larger and smaller, respectively, than for collision with j = 0. Except at low collision energies the $Qi(2)(E)$ cross sections have the same collision energy dependence as that shown in Fig. 3 and are not plotted here.

We have employed the state-specific transport cross sections $Qi(n)(E)$ to compute state-specific collision integrals [Eq. (2)] for (n, s) = (1, 1), (1, 2), (1, 3), and (2,2). Figure 4 presents computed state-specific collision integrals $Ωi(1,1)$ for the O(³P_j) fine-structure levels in collisions with H atoms at several temperatures. We see that at low temperatures, the O(³P_j=2) level has the larger value of the state-specific collision integral, while at high temperatures this level has the smallest value, This is consistent with the collision energy dependence of the transport cross sections plotted in the left-hand panel of Fig. 3.

We plot in Fig. 5, the state-specific collision integrals $Ωi(1,1)$ for the O(³P) fine-structure levels in collisions with H₂ in specified rotational levels at two temperatures. We see that at both 300 and 1500 K the collision integrals for collision of O(³P) fine-structure levels with the lowest H₂ para and ortho rotational levels (j = 0 and 1, respectively) are very similar in magnitude, while those of the next higher rotational levels (j = 2 and 3) are significantly larger. We also see that for the higher H₂ rotational levels, the state-specific collision integrals at 300 K have a different dependence upon the O(³P) fine-structure level, with the collision integral for collision with O(³P_j=2) being the smallest and largest for H₂ j = 2 and 3, respectively. At the higher temperature (1500 K), the differences in the magnitudes of the collision integrals are much less.

We have used the state-specific collision integrals, such as those plotted in Figs. 4 and 5 to compute the collision integrals defined in Eq. (1). These are presented in tables in the supplementary material. For O(³P)–H(²S), the sum over internal levels [namely, the O atom fine-structure levels] in Eq. (1) is complete. As indicated in Sec. IV B, the sum over H₂ rotational levels covers j = 0-3. These levels include 99.5% and 66% of the total H₂ rotational population at 300 and 1500 K, respectively. The neglect of higher rotational levels in the sum in Eq. (1) for O(³P)–H₂ should not cause a significant error in the calculation of the collision integrals at high temperatures since the state-specific integrals do not show a strong dependence upon H₂ rotational level at high temperatures [see Fig. 5(b)].

The collision integrals reported in the supplementary material can be used to compute transport properties such as the diffusion coefficient D [Eq. (7)] and the ratios of collision integrals [ Eqs. (8)–(10)] required to describe thermal diffusion. We present in Fig. 6, diffusion coefficients for O(³P)–H and O(³P)–H₂ obtained from our QS calculations. We compare in Fig. 6(b), diffusion coefficients computed with the O–H₂ potential of the present study and that of Alexander.³⁰ Also included in Fig. 6 are conventional calculations based on LJ(12-6) isotropic potentials. For these calculations, we employed ϵ and σ parameters derived from the parameters in the Sandia compilation.² The well depths ϵ and length parameters σ were computed as the geometric and arithmetic means, respectively, of the values for the like pairs [O(³P): ϵ = 80 K, σ = 2.75 Å; H(²S): ϵ = 145 K, σ = 2.05 Å; H₂: ϵ = 38 K, σ = 2.92 Å ]. The following values were employed here: ϵ = 107.7 K and σ = 2.40 Å for O(³P)–H(²S) and ϵ = 55.1 K and σ = 2.84 Å for O(³P)–H₂.

We see in Fig. 6(a) that the LJ diffusion coefficient for O(³P)–H is slightly larger than that obtained in our quantum scattering calculations. This contrasts with previous comparisons of diffusion coefficients from the two types of calculations for other systems,²¹ for which the LJ diffusion coefficient is smaller, in some case dramatically smaller.^10,19,21,22 The one exception, aside from O(³P)–H, is H–H. For both these systems, one of the states emanating from the separated collision partners is very strongly bound [H$2(X1Σg+)$ and OH(X²Π)]. The fact that the quantum scattering calculations yield diffusion coefficients usually larger than the corresponding LJ estimates can be traced to a steeper repulsive wall of the LJ 12-6 potential than that of a typical ab initio potential.¹⁷ 

For the case of O(³P)–H₂, the quantum scattering calculations yield diffusion coefficients significantly larger than the LJ estimate. Moreover, we see that the calculations with the present potential and that of Alexander³⁰ yield diffusion coefficients in very good agreement. This indicates that the slight differences in these potentials have only a minor effect on the computed transport properties. In future work, we will extend our quantum scattering calculations to other collision pairs of importance in combustion.

See the supplementary material for a table of the O(³P)–H(²S) interaction energies as a function of R, tables of spectroscopic parameters for the OH X²Π and A²Σ⁺ states, a table of the SF expansion coefficients [Eq. (23)] as a function of R for the O(³P)–H₂ interaction, tables of collision integrals, and plots of the collision integral ratios A^∗, B^∗, and C^∗ for O(³P)–H(²S) and O(³P)–H₂.

This work was supported by the Chemical, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy, under Grant No. DESC0002323. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We gratefully acknowledge helpful correspondence with Roman Krems.

This appendix presents an expression for the transport cross section involving the collision of two collision partners both with internal structure. As discussed previously,¹⁶ we require for the calculation of transport cross sections the moments A_λ [Refs. 67 and 68] of the Legendre polynomials P_λ(cosθ) over the state-to-state differential cross section. In this case, the moment equals

where f is the scattering amplitude, which is a function of θ. In Eq. (A1), the initial and final rotational angular momenta of the collision partners are denoted j₁, j₂ and $j1′$⁠, $j2′$⁠, respectively, and the m′s are the corresponding SF projection quantum numbers; k is the initial wave vector, and [x] = (2x + 1). (It should be noted that the A_λ are defined slightly differently in Refs. 16, 67, and 68.)

The n = 1 and 2 transport cross sections can be expressed in terms of the A_λ as

where E and E′ are the initial and final collision energies, respectively.

After some angular momentum algebra, the moment in Eq. (A1) equals

In Eq. (A4), the L′s are the orbital angular momenta, the J′s are total angular momenta of the collision complex, j₁₂ = j₁ + j₂, $j12′=j1′+j2′$⁠, (:::) are 3j symbols, and $:::$ are 6j symbols. The indexed T-matrix elements are expressed in a space-frame basis; these are obtained from time-independent quantum scattering calculations.