Quantum and quasi-classical collisional dynamics of O₂–Ar at high temperatures

High-temperature gas dynamics trigger atom-molecule collisions and determine the efficiency of energy transfer between translational and internal modes and temperature profiles around aircraft in hypersonic aerospace applications. The high convective speed around a vehicle disrupts the distribution of states in ambient air—due to the shock wave occurring close to the vehicle’s body — thereby inducing a nonequilibrium distribution.^1,2 The reorganization of molecular states and the return to a new equilibrium distribution requires a characteristic amount of time. Energy is exchanged in molecular collisions between the collision partners, and between translational, rotational, and vibrational modes. The readjustment of translation and rotational modes typically happens much faster than the readjustment of the vibrational modes or chemical reactions. Thus, at moderate hypersonic speeds, i.e., less than 3 km/s, the main processes governing the energy dissipation are vibrational energy transfer and dissociation.^3–12 

An understanding of vibrational relaxation is therefore crucial for realistically modeling hypersonic air flows. Vibrational relaxation of molecular species is governed by two processes, vibration-vibration (VV) energy transfer and vibration-translation (VT) energy transfer. In the former process, vibrational energy from one diatomic is transferred to another diatomic, effectively exciting and de-exciting vibrational levels of the collision partners. The net energy exchange between vibrational degrees of freedom (DOF) strongly correlates with vibrational anharmonicity. Such energy exchange takes place rapidly and with a high probability at low translational temperatures.¹³ On the other hand, vibration-to-translation (VT) energy transfer is a dominant mechanism of relaxation at high temperature conditions. In the VT energy transfer, the energy corresponding to a vibrational (de-)excitation is exchanged with the translational degrees of freedom.

The processes investigated here take place in the upper atmosphere, where the gases N₂ and O₂ are the most predominant compounds. Here, we choose to investigate the O₂ + Ar scattering system for its simplicity and availability of experimental data on vibrational relaxation;^14,15 once we have established our methodology, the routine can be applied to more complex systems such as O₂ + O and O₂ + O₂. We use the quantum dynamics multi-configuration time-dependent Hartree method (MCTDH)^16–18 and a quasi-classical trajectory (QCT) method to derive the rate coefficients for vibrational transitions ν → ν′(ν≠ν′) in O₂–Ar scattering in a wide range of translational temperatures,

In this inelastic collision, there is an energy exchange between the vibrational degrees of freedom of the diatomic molecule and the translational DOFs of the atom. This process of energy transfer is characterized by its rate coefficient, and we can derive rate coefficients for the excitation as well as the de-excitation collisions. These rates are connected to the vibrational relaxation time, and to the energy-equivalent vibrational temperature of the gas. Under nonequilibrium conditions, the vibrational, translational, and rotational temperatures in gas flows around hypersonic vehicles are not equal, and the distribution of states can deviate from a Boltzmann distribution. Therefore, nonequilibrium effects in the dynamics need to be taken into account. The present work is a first step in this direction, since the transition probabilities calculated here can be used in any nonequilibrium distribution to derive the rate coefficients. Any desired distribution of states can be input into our code for the rate calculation, yielding a nonequilibrium rate in a manner of seconds. The generated rates can then be used in the computational fluid dynamics to aid the development of thermal protection systems and flow control strategies.

The quasi-classical methods are typically adopted for calculation of state-specific VT rates due to the significantly higher cost of more detailed quantum models.^3,6,9,19,20 Here, we present a fully quantum dynamical approach for obtaining VT rates using a corresponding summation over the quantum states. At low temperatures, the QCT method fails to predict the VT rates accurately for molecular systems with strong repulsive potentials. This is due to the lack of detailed tunneling events and interference effects that occur with sufficient probability near threshold that the classical trajectories necessarily miss. Using quantum simulations, we are able to bridge this gap and are able to estimate the relative accuracy of quantum and QCT rates. These rates can be converted into vibrational relaxation times, using derivations commonly applied in gas dynamics,¹ and thus can be compared to experimental data. In the temperature range discussed here (1000 K to 10 000 K), we obtain very good agreement between calculated and measured vibrational relaxation times. The important effects of vibrational energy relaxation (VER) between the colliding oxygens resolved earlier^21,22 at low temperatures do not arise here because our mixtures consist of dilute O₂ which necessarily suppresses VER. Meanwhile, Hase and co-workers²³ found agreement between theory and experiment for VER in N₂ + C₆F₆ using a similar approach as pursued here.

In the following, we discuss the theoretical basis for our quantum dynamics and QCT simulations in Sections II A and II B. The standard model relating vibrational relaxation times to rate coefficients is discussed in Section II C. We present transition probabilities, scattering cross sections, rate coefficients, and vibrational relaxation times in Section III. Our findings are discussed in Section IV.

The physical basis of any dynamics calculation is the representation of the forces and their gradients. We use a potential energy surface (PES) with a pairwise combination of the Hulburt-Hirschfelder analytic function²⁴ with fitted parameters for the O₂ potential energy curve^25,26 and the Buckingham potential²⁷ for the O–Ar interaction.²⁸ This methodology of constructing the many-body potential energy surface is widely adopted for extensive quasi-classical studies of diatom and inert atom collisions.^29,30 The validity of this approach is based on the electronic structure of the Ar ground electronic state, which has a closed-shell configuration.³¹ The parameter values are given in Table I.

The zero-point energy for the O₂–Ar system using this PES is 788 cm⁻¹ [exp E_ZPE = 787.4 cm⁻¹ (Ref. 32)], and the first vibrationally excited level lies at E_ν1 = 1556 cm⁻¹ [exp E_ν1 = 1580 cm⁻¹ (Ref. 33)], as calculated using the MCTDH method.³⁴ 

The accuracy of the dynamics calculations depends on the accuracy of the PES employed in the dynamical simulation. To validate the accuracy of the (parameterized) model PES in Eq. (2) with respect to the O₂–Ar interaction, we performed CCSD/6-311(d,p) calculations using the GAMESS-US program package.^35–37 Only CCSD results are shown in Figure 1 because the inclusion of triples led to negligible corrections. The agreement between the model PES and CCSD results is excellent at an impact angle of 0° (Fig. 1). The transition probability is the highest for this head-on collision, and vibrational (de-)excitation becomes less likely for increasing impact angle. The comparison between the PES and CCSD for other impact angles shows an increasing deviation in the onset of the repulsive wall with increasing impact angle. However, the shape of the wall is captured correctly and this is known to be important as it would otherwise lead to deviations in the dynamics. The error in the onset of the repulsive wall can be attributed to the neglect of three-body terms in the O₂–Ar interaction which is represented in Eq. (2) only through the sum of distinct O–Ar interactions. The ensuing shift to lower r₂ in the onset of the repulsive wall is exaggerated in our model PES, with the largest effect for the largest impact angle 90°. This in turn reduces the effective scattering cross section of the scattering process, and may therefore lead to slightly reduced VT cross sections and rates. However, as shown below, the agreement for the small angles is excellent indicating that the incorrect position of the onset is leading to a negligible error. Thus this study, featuring the use of the same PES in determining the rates by several methods provides a mechanism to clearly benchmark the quantum effects. It does not provide exact rates because it is limited by the possible errors in the approximate PES.

The implementation of the MCTDH method¹⁶ by the group of H.-D. Meyer³⁴ is used in the present work to compute the quantum-mechanical rates and cross sections. Here, the nuclear wave function is approximated as a product of single-particle functions (SPFs)

for f DOFs, and n SPFs per DOF. The SPFs are represented as a linear combinations of the primitive grid basis. The equations of motion (EOM) for the SPFs and the time-dependent coefficients A_j1…jf(t) involve the computation of the mean fields $〈 H R 〉 j l κ$ at each time step, which is the most computationally demanding step in the propagation. By rewriting the Hamiltonian such that h_κ operates only on the κ degree of freedom, and approximating the residual Hamiltonian H_R as a sum of products of single-particle operators,

the multi-dimensional integrals required in the evaluation of the mean fields are circumvented. The kinetic operator is usually separable between the DOF, but the potential term must be fitted to this form.³⁸ We achieve a product representation of the potential by fitting it using the potfit program, as part of the MCTDH package. To reduce the grid size required in the computation, and to compute the reactive and inelastic fluxes, complex absorbing potentials (CAPs) are added to the Hamiltonian for each DOF.^39,40 These CAPs have the form of

where in our calculations, the order of the CAPs is three (b = 3), and the grid point Q_c denotes the position at which the CAP is switched on, with θ(Q − Q_c) the Heaviside step function. The parameter η denotes the CAP strength.

The reactant Jacobi coordinates comprise our choice for the coordinate system and consist of the vibrational coordinate r₂ (O–O distance), translational coordinate r₁ (distance of Ar to the center of mass of the O₂ molecule), and the angle θ between r₁ and r₂. In this coordinate system, the center of mass is constant. We can define the quantum numbers as ν for vibrational quanta, j for rotation of the diatomic, J for the total angular momentum quantum number, and K for the projection onto the body-fixed axis. The Hamiltonian of the scattering system is defined in the coupled-states (CS) approximation, where J and K are “good” quantum numbers (neglecting J, K-coupling):

The reduced masses μ_r1 and μ_r2 are the reduced masses of the O₂–Ar and O₂-systems, respectively. The PES V(r₁, r₂, θ) is specified in Eq. (2). To obtain scattering cross sections and rate constants, inelastic transition probabilities need to be calculated. In principle, the inelastic transition probabilities can be obtained by two different methods. For state-specific transition probabilities P_if(E) the method developed by Tannor and Weeks⁴¹ is most suitable, where the S-matrix elements are obtained from the cross-correlation function C_if(t). Final-state averaged transition probabilities can then be obtained by summation over the final states f

A more direct way to obtain the final-state averaged transition probabilities is by evaluating the flux through the CAP W (see Eq. (6)). This directly yields P_i(E) as

The energy distribution of the initial and final wave packets, Δ_i(E) and Δ_f(E), are obtained by diabatic correction of the translational DOF, where the influence of the interaction potential along r₁ is corrected for. The energy distribution is evaluated as the overlap of the translational SPF with a distorted wave (rather than a plane wave for a Gaussian initial SPF),

where the distorted wave Ψ_dist is represented by a WKB-wave (Wentzel-Kramer-Brillouin) multiplied with the r₂ eigenstate SPF. The associated local momentum k(r₁) is defined as the momentum available for r₁,

with $V ̄ ν ( r 1 )$ the diabatic mean-field potential. In contrast to the initial wave packet, its energy distribution depends on the internal state of the scattering system (vibrational state ν, rotational state j, total rotational angular momentum J, and its projection K onto the body-fixed axis). Since for the scattering cross section, the transition probability needs to be expressed in terms of the collision energy E_col (purely translational energy), the energy shift due to the internal state of the molecule is corrected for.

The flux can be projected onto final vibrational states ν′ by applying a projection operator to the wave packets in Eq. (11). In this way, the required transition probabilities from initial to final vibrational state, $P ν j → ν ′ J K$⁠, can be obtained. Since we are interested in highly averaged properties such as cross sections and rate constants, we calculate only final-state averaged transition probabilities which already implicitly contain the summation over final rotational states j′, J′, and K′.

The integral cross sections σ_νj→ν′(E_col) are then obtained as⁴² 

The integral cross sections are required over an energy range of 5.0 eV. Within this energy range, there is a large amount of internal states that need to be included. In order to reduce the computational effort, we assume that the dependence of $P ν j → ν ′ J K$ on J can be obtained by an energy shift as for example in the J-shifting approximation,⁴³ or the quantum analog of the capture model,⁴⁴ 

Here, the shifting energy is the barrier associated with an effective potential. The distance r_1,barr is associated with the distance of the rotational barrier, and can be related to the QCT parameter b_max. The K-dependence of $P ν j → ν ′ J K$ is assumed to be negligible, so that

The integral cross sections are then rotationally averaged over the initial rotational states of the diatomic,

Here, β = (k_BT)⁻¹ with k_B the Boltzmann constant. Due to the O₂ nuclear spin state, only odd rotational states can be populated. Consequently the even states can be ignored in the sum. Translational averaging of the cross sections yields the vibration-translation transition rate coefficients

with g_s = 1/3 the factor of spin degeneracy, and k in units of cm³/s. The value of g_s comes from the ratio of the degeneracies of PES and reactants.²⁹ Below, we calculate vibrational transition rates from vibrational level ν = 1 to all vibrational levels ν = x, {x = 0, 2, 3, …, 36}, thus obtaining rates of vibrational deactivation, k_ν→ν′(ν > ν′), as well as rates of vibrational activation, k_ν′→ν. The rates of vibrational activation can be transformed into rates of vibrational deactivation using the principle of detailed balance, assuming the system is in an equilibrium state,

where E_vib is the difference in vibrational energy between the ν and ν′ states.

For another point of comparison with experimental and empirical data, we also compute temperature-dependent transition probabilities. First, the final-state averaged transition probabilities are averaged over J,

These probabilities are then averaged over rotational states and integrated over the energy range E_col while normalizing the integral,

The QCT method is used to reduce the computational cost of generating the O₂–Ar vibrationally resolved cross sections. Details of this approach are given elsewhere.⁴⁵ Here, we briefly discuss the basics of this method. The initial rovibrational state of the oxygen molecule, $ν , j$⁠, is fixed in each batch of trajectories. The internuclear separation in the isolated oxygen molecule is initialized by the random sampling of oscillation phase based on the period of vibrational motion. The initial separation of the target and projectile particles is given by 15 Å. Hamilton’s differential equations are solved by the Adams-Moulton method to 11^th order of accuracy as described in Ref. 45. The impact parameter, b, is sampled with a step size of 0.1 Å. Every batch contains 2 × 10³ trajectories. Analysis of the final state is performed according to Ref. 45. Each trajectory is integrated with an error in the total energy not exceeding 10⁻⁴%. A trajectory is terminated after the distance between products exceeds the initial separation. Each such trajectory is then classified into one of two possible channels: inelastic or elastic. The probability, P, of the state-specific transition $ν , j → ν ′ , j ′$ at a given collision energy E_col is defined as follows:

where b_max is the impact parameter at which only elastic collisions are observed, N_{νj→ν′j′} and N are the number of trajectories with the desired transition and the total number of trajectories in the current batch of impact parameter, respectively. The cross section of a bound-bound transition is calculated as follows:

where P is sampled over the five variables, not explicitly mentioned in its argument in Eq. (28): the impact parameter b, the initial azimuthal and polar orientations of the target molecule, the initial orientation of the target molecule angular momentum and the initial vibrational phase of the target molecule. The present QCT calculations study dynamics of O₂–Ar collisions in the range of collision energy from 0.1 to 5 eV. This range is covered by 35 intervals of uneven length. Stratified sampling of the impact parameter is adopted to increase accuracy of statistical modeling.

The present study features the analysis of rates assuming trans-rotational equilibrium. In order to do so, the cross section of transition from the initial rovibrational state to any final rotational state of the final vibrational level is computed using Eq. (20). The corresponding rate of transition is obtained using Eq. (21).

The rate coefficients for vibrational deactivation can be related to vibrational relaxation time as derived from experimental data by Millikan and White,⁴⁶ Camac¹⁴ and recent measurements by Hanson and co-workers.¹⁵ The vibrational relaxation parameter τ, is expressed in Landau-Teller theory as¹ 

with pressure p and energy difference E_vib of the vibrational states ν′ and ν. In this approach, the de-excitation from the first excited state (ν = 1) to the ground vibrational state (ν′ = 0) is considered. The exponential term in Eq. (29) corresponds to the contribution of the reverse excitation process. When the excitation of vibrational ground state is improbable (i.e., E_vib ≫ k_BT), Eq. (29) simplifies to

Equation (30) is based on the assumption that the population of vibrational level ν = 1 is very small compared to the vibrational ground state. However, at temperatures relevant to hypersonic flows, the resulting population of excited vibrational states and the contribution of multiquantum transitions lead to relaxation times that are significantly different from those obtained assuming a two-state vibrational model. Given the complete set of relevant transitions for each internal state, one can then define the average relaxation time for the entire rovibrational ensemble. This was previously done for N₂–N (Ref. 10) and O₂–O (Ref. 47) molecular systems using the database of rates generated by the QCT method. The details of such simulations, involving the master equation solution, are given elsewhere.⁴⁷ In the present work, the O₂–Ar vibrational relaxation time is obtained using the e-folding method⁴⁸ from the O₂–Ar set of vibrationally resolved transitions rates assuming rotational equilibrium. Unfortunately, due to the high cost of the fully quantum approach, such simulations using the MCTDH method are currently computationally intractable.

In experiments, the vibrational relaxation time is obtained as the slope of the straight line fit to the measured profile of $log 1 − E / E f$ over laboratory time t, with E the energy in the vibrational state, and E_f the final vibrational energy. This involves a number of approximations, such as vibrational harmonicity and single-state quantum transitions, as in Landau-Teller theory. An empirical model of the relation between vibrational relaxation time and temperature is detailed in Millikan and White⁴⁶ for various scattering systems. The same dependence is observed for a large number of species with a strong repulsion between projectile and target particles. Below, we compare calculated vibrational relaxation times (for the above equations) with experimental ones, and with vibrational relaxation times obtained using master equation calculations based on our computed rates.

For all quantum propagations, we use the constant-mean field (CMF) integration scheme, with a Bulirsch-Stoer (BS) extrapolation scheme for the propagation of SPFs and a short iterative Lanczos (SIL) algorithm for the A-vector propagation. This is implemented in the MCTDH package version release 84.8.⁴⁹ The DVR grids used to represent the primitive basis run from $r 1 min =3.89 a 0$ to $r 1 max =10.0 a 0$⁠, and $r 2 min =1.78 a 0$ to $r 2 max =6.00 a 0$⁠, where a Sine DVR is used for both DOFs, and the number of grid points is n_r1 = 512 and n_r2 = 112, respectively. The angular DOF is represented on a Legendre DVR with 235 basis functions, using symmetry restrictions and only odd rotational functions as the even ones are not allowed for O₂. The CAPs start at 7.50 a₀ (r₁), and 5.00 a₀ (r₂), with strength $0.004 a 0 − 1$⁠, and order three.

The initial wavepacket Ψ₀ is represented by a Gaussian SPF in the translational DOF r₁ (with center at Q_r1 = 6.50 a₀, width ΔQr₁ = 0.10 a₀, and initial momentum p₀ from −20 a.u. to −105 a.u.). The SPF for the bound DOF r₂ is initially setup as a vibrational eigenstate of the O₂ one-dimensional Hamiltonian, with ν = 1 as quantum number of the initial vibrational state. The angular DOF is represented by an associated Legendre function with angular momentum of the diatomic from j = 1 to j = 217. The Hilbert space in the vibrational DOFs, r₁ and r₂, is spanned by 14 SPFs, while the space in the angular DOF, θ, is spanned by 10 SPFs.

To cover the entire energy range of translational energies from 0.0 eV to 5.0 eV for each selected initial state ν, j, six different propagations are carried out. The resulting transition probabilities are then connected, interpolated, and smoothed to yield transition probabilities P_ν,j(E) over a fine grid supporting the desired translational energies. For very low impact energies, the transition probabilities are set to zero, since for these energies the deconvolution of translational and rovibrational contributions to the wave packet is difficult to achieve. To remedy this situation, one would need to implement an explicit adiabatic correction scheme for this system. However, for simplicity we use the related scheme already available in MCTDH⁴⁹ to address H+H₂ scattering so as to avoid significant recoding.

Convergence in the cross sections and rates required summation over rotational momenta J from J = 0 to J = 800. Using test calculations for J≠0, we also found that the distance associated with the barrier potential is approximately r_1,barr = 3.8 Å.

The calculated transition probabilities for vibrational deactivation ν = 1 → ν′ = 0, J = 0, and selected rotational momenta of the O₂ diatomic, j, are shown in Fig. 2. For low j (j < 41), the onset of the transition probabilities lies between 1-2 eV. The onset of the QCT probabilities is more abrupt than the quantum probabilities, and may be attributed to the difficulty in resolving the low-energy region accurately with QCT. Additionally, tunneling effects which may play a role in the low-energy region can only be resolved exactly by a quantum method. Both QCT and MCTDH transition probabilities level off at a value of P(E) ≈ 0.3. The transitions from these low-lying rotational states to the vibrational ground state are triggered by a fairly high amount of translational impact energy.

The rotational states of intermediate rotational energy, from j = 43 to j = 151, exhibit a shift of the onset of the transition probability to lower energies; also, their transition probability levels off at higher impact energies. These rotational states are still bound in the potential, and show a tendency to vibrational deactivation for low to intermediate translational energies. For high translational energies, the process of vibrational activation to the next-highest vibrational state, ν = 1 → ν′ = 2, becomes a competing process and therefore reduces the transition probability. In this region, the QCT and MCTDH transition probabilities compare very well.

For the rotational states of high rotational energy, 151 ≤ j ≤ 191, the onset shifts even more to lower impact energies, and the competition with vibrational activation becomes more significant. The transition probabilities level off quickly, and QCT and MCTDH results are in good agreement, while the low-energy region is resolved better in the quantum calculations.

The rotational states of highest rotational momentum, j ≥ 193, show a similar behavior as the ones of high rotational momentum. However, in the quantum calculations, it was difficult to resolve this region accurately. We attribute this to the nature of these rotational levels: only a few of these are still bound, and the rotational levels above j = 195 are unbound in our MCTDH calculations. Therefore, the molecule starts to dissociate even before the impact of the Ar atom, and we chose the propagation settings (i.e., numerical grids) to represent bound wave packets. To describe this region close to the dissociation energy, we would need to extend the vibrational grid and choose a different CAP. In the QCT calculations, this energy region could be resolved without any issues.

In summary, we obtain good agreement between the MCTDH and QCT transition probabilities, especially in the region of intermediate to high rotational momenta. The MCTDH calculations clearly show that it is necessary to carry out quantum calculations to resolve the low-energy impact energies accurately, while the high-lying rotational states can be accessed using the QCT method without requiring significant computing cost or special care.

Summing over the rotational momenta J and K (MCTDH), or integrating over the impact parameter b (QCT), we can obtain the scattering cross sections for selected initial rotational states j as shown in Fig. 3. The low-lying rotational states (j < 41) show a steady rise of the scattering cross sections and reach about 3-4 Ų, with one exception: in the quantum calculations, the lowest state j = 1 reaches a scattering cross section of almost 7 Ų at 5 eV. This can be attributed to the high number of J-states that are accommodated in the averaging. Again, in the quantum method, the low-energy region is resolved better than in the QCT method because of the issues stated earlier concerning low impact energies and rotational states.

The rotational states of intermediate rotational energy, from j = 43 to j = 151, again show a shift of the onset of the cross sections as well as a decreasing cross section for higher translational energies. This effect however is mitigated compared with the transition probabilities, due to the averaging. For the very high rotational states (j > 151), the decrease in the cross sections due to competing processes becomes more significant. Additionally, the QCT method predicts a higher cross section than the MCTDH method for higher j; almost double for the low impact energies and j = 181. We attribute this deviation to the differences in the low-impact energies between QCT and MCTDH. These differences were small as the level of transition probabilities become significantly more pronounced when the averaging over the body-fixed rotational states is carried out.

The initial-state selected cross sections can be rotationally averaged over initial rotational momentum j leading to temperature-dependent cross sections such as those shown in Fig. 4. With increasing temperature, the onset of the rotationally averaged scattering cross sections shifts to lower energies. Comparing the quantum to the quasi-classical calculations, the quantum scattering cross sections rise slower and level off at around 3.0-3.5 Ų. The QCT cross sections show a slightly more prominent rise and level off between 3.0-4.0 Ų, where with higher temperature the differences between the cross sections are more pronounced than for the MCTDH results.

The rotationally averaged scattering cross sections can be averaged over translational energies to yield rates of vibrational (de-)activation. The VT rates from the MCTDH and QCT calculations are shown in Fig. 5, for transitions ν = 1 → ν = x, {x = 0, 2, 3, 5, 10, 15}.

The earlier onset of transition probabilities and therefore scattering cross sections in the MCTDH calculations results in a higher VT rate at lower temperatures. Also, the low-lying rotational states were resolved more accurately in the quantum method; these states dominate the low-temperature region. However, QCT and MCTDH rates agree for temperatures of 2000-3000 K and higher. At very high temperatures, the QCT rates tend to be slightly larger than the MCTDH rates. The MCTDH and QCT deactivation rates ν = 1 → ν = 0, ν = 1 → ν = 2, and ν = 1 → ν = 3 show a very good agreement for higher T, while for the rates ν = 1 → ν = 5, ν = 1 → ν = 10, and ν = 1 → ν = 15, the deviation between quantum and classical rates becomes larger also at high temperatures. Generally, the differences between MCTDH and QCT scattering cross sections are less pronounced than the differences in the transition probabilities due to the strongly averaging procedure that converts molecular trajectory information into macroscopic quantities.

The VT rates cannot be measured directly in the experiment, the only way to compare the calculated data to the experiment is via the vibrational relaxation times. However, this involves a conversion of the experimentally measured data into relaxation times, as well as a conversion of the theoretical data; this introduces some uncertainties in the comparison, as it involves inherent assumptions such as the adsorption coefficient being independent of the rotational state and vibrational harmonicity. Nevertheless, the experimental measurements give insights into the accuracy of the theoretical calculation and the applicability of the quantum and quasi-classical methods. In the present studies, two experimental data sets are used to assess the adopted computational models. The vibrational relaxation times in O₂–Ar collisions were reported for the first time by Camac,¹⁴ where the ultraviolet absorption technique was applied to collect data in O₂–Ar mixtures with the post-shock translational temperatures between 1200 and 8000 K. The absorption in Schumann-Runge bands was measured at 147 nm wavelength. Recently, the O₂–Ar molecular system was revisited by Hanson and co-workers¹⁵ using similar methods of laser absorption spectroscopy in the range between 1000 and 4000 K. For this purpose, a picosecond pulsed Ti:Sapphire laser was used to generate tunable UV light with the wavelength between 210 and 230 nm. The new data have much less scatter compared to the earlier measurements by Camac.

The calculated and experimental vibrational relaxation times are shown in Fig. 6. The present theoretical data are constructed from two different equations, Eqs. (29) and (30), and for QCT also from the solution of master equations. The MCTDH data are shown with black curves, while the QCT results are shown in red. The curve attributed to Millikan and White is obtained from their empirical model using the specific parameters for the O₂–Ar system.⁴⁶ Symbols represent the experimental data from Ref. 14, for different admixtures of O₂ in Ar, and from Ref. 15.

The experimental and empirical relaxation times are in excellent agreement which is not surprising since the Millikan-White model is based on the experimental data from Ref. 14. The data by Camac are confirmed by the recent experimental measurements by the Hanson group.¹⁵ The relaxation time derived from the MCTDH calculations using Eq. (30) is in very good agreement with the experimental data. The QCT relaxation times strongly overestimate these results at temperatures below 4000 K. This is in agreement with the results shown in Fig. 5: the QCT VT rates are generally smaller than the MCTDH rates. This is due to the insufficient statistics accumulated by the QCT method. Because the true O₂–Ar vibrational relaxation time increases by orders of magnitude in the low temperature region, the QCT method requires a proportional increase of the number of trajectories in each batch. This makes the QCT calculations impractical for the system of interest at low translational temperatures. Recently, a similar failure of the quasi-classical methods regarding the O₂–N₂ system was described in Ref. 50. Notably, once the translational temperature is high enough to accumulate the statistics using a reasonable amount of trajectories, the QCT and MCTDH methods demonstrate a very good agreement.

In our dynamics treatment, we only involve the O₂ electronic ground state, 3Σ_g. The other low-lying electronic states of O₂ also become accessible at the high temperatures discussed here, and part of the deviation between theory and experiment may also result from the neglect of nonadiabatic couplings at high temperatures.

In addition to nonadiabatic couplings, failures in the assumption that the system is in an equilibrium state (in the sense that Boltzmann distributions are satisfied for translational, rotational and vibrational degrees of freedom at the same temperature), may also lead to a further deviation between theoretical and experimental results. To confirm this origin of the deviation, details of the state distributions of the system as well as the related temperatures would need to be experimentally available. Using such data, we would be able to generate nonequilibrium rates and lifetimes, and could investigate the impact of the nonequilibrium effects on the accuracy of the calculation. This will be a future point for studies, and can easily be implemented in our VT rate and lifetime calculations.

Finally, the approximation that J-K coupling can be ignored [centrifugal sudden (CS) approximation], could introduce uncertainties in the quantum dynamics simulations. On another system, the H₂ + D₂ scattering system, J-K coupling has been shown to enhance the reactivity.⁵¹ In line with this finding, a more accurate kinetic energy operator in which J-K coupling is treated accurately, could lead to slightly enhanced transition probabilities and therefore higher VT rates (lower lifetimes). At this point, we chose to neglect J-K coupling due to the additional computational demand it introduces in the calculation, and since we also do not treat the J-K dependence of the transition probabilities explicitly.

In summary, the use of Eq. (29) leads to underestimates in relaxation times, especially at high temperatures. This is due to the influence of the reverse excitation process. It affects both QCT and MCTDH methodologies. At low temperatures, the excitation events are rare; relaxation times via Eqs. (29) and (30) are nearly identical. The QCT master equation relaxation time lies in between the two approaches. This more accurate approach at calculating relaxation times takes into account multiquantum transitions, and was only carried out for the QCT rates, since the full set of ν = x → ν = y, {x, y = 0, 1, …, 36} transition rates are required. At this point, we only generated quantum-mechanical rates involving the first vibrational excited state. However, from the general comparison of master equation and approximated relaxation times using Eqs. (29) and (30), it can be expected that the quantum-mechanical master equation relaxation times would lie in between the relaxation times from Eqs. (29) and (30).

While both MCTDH and QCT methods can be used to predict highly averaged properties that can be measured in macroscopic experiments, they also provide detailed microscopic parameters and time-dependent characteristics that enter into classic gas-kinetic theory and other microscopic models. One such quantity is the single-state quantum probability of O₂–Ar deactivation from the first excited state, P_1→0. Using QCT methodology, P_1→0 can be obtained for each initial rotational quantum number and collisional energy directly using Eq. (27), where the right-hand side should be summed over all j′. The probability at a given trans-rotational temperature is the average of $P 1 → 0 j , E col$ over the Boltzmann ensemble of initial rotational states and Maxwell distribution of thermal velocities. The MCTDH probability is estimated in an equivalent manner using Eq. (23) instead of (27). The deactivation probability can be estimated from experimental measurements through

where Z is the collision frequency, and K_1→0 denotes the rate of deactivation similar to that in Eq. (21). The collision frequency is estimated as $Z=π r 2 v ̄$⁠, where r is the effective radius of collision and $v ̄$ is the expected value over the velocity distribution. In the original work by Camac, r was set to 3.43 Å.

The probabilities of deactivation at several trans-rotational temperatures are shown in Table II. As expected, the QCT method strongly underestimates P_1→0 at low translational temperatures. The low deactivation probability leads to overestimation of vibrational relaxation times as can be seen in Fig. 6. At T = 2400 K and higher, QCT and MCTDH both match the Camac estimates quite accurately. Below these temperatures, only the MCTDH results should be considered as accurate. The minimal impact parameter at which only elastic collisions are observed is also given in Table II. The impact parameter b_max increases toward higher temperatures, which can be understood from the larger O₂ interatomic distance in excited rotational states. The difference between the QCT value of b_max and the radius of collision, used by Camac, is not more than 10%, but can nevertheless be a source of discrepancy in the relaxation times at low temperatures. On the other hand, the relative accuracy in the value of b_max obtained in the QCT calculations in comparison with experiment, is indicative of the relative accuracy of the underlying PES used in the present studies.

We have found that the transition probabilities and scattering cross sections calculated with the MCTDH and QCT methods compare quite well, especially for moderate to high translational energies. The O₂–Ar scattering for high initial rotational states are better resolved with the classical method, while the low-lying rotational states and low impact energies are more accurately obtained using MCTDH. After averaging the cross sections over impact energies to generate rates of vibrational deactivation, we found that quantum-mechanical and classical rates agree quite well above 2000-3000 K. Below that temperature range, the QCT method fails to predict the rates accurately, as can be expected from a quasi-classical method that does not completely account for tunneling and other low-energy quantum effects. For the remaining temperature range, from 3000-10 000 K, both sets of rates show very good agreement.

We also determined a standard procedure for obtaining fully quantum-mechanical and quasi-classical VT rate calculations that are directly comparable at high temperatures. In this approach, we make use of the capture model to generate transition probabilities for rotationally excited states in the body-fixed frame. Both the MCTDH and QCT methods have significant computational requirements. For the quantum-mechanical calculations, six propagations per O₂ rovibrational state were carried out to cover the entire range of impact energies. Each propagation takes about 6 h on average on a single core on the Stampede cluster [within an Intel E5 8-core (Sandy Bridge) coprocessor] so that we used 36 h of computation for each rovibrational state. Thus, in total, for initial vibrational state ν = 1 and all 217 initial rotational states, the quantum-mechanical calculations took about 7800 h computing time on a single node. The QCT calculations for ν = 1 and all possible j require 790 h with the fixed size of a batch equal to 2000 trajectories. Hence, each rovibrational state takes approximately 4 h to obtain the data in the entire range of collision energies. The computational time is directly proportional to the size of a batch. Trajectory propagation for excited pre-collisional vibrational states requires higher computational time due to increased impact parameter b_max. Overall, the MCTDH approach appears to be an order of magnitude more expensive than the QCT calculations. However, in light of the much better resolution of the vibrational relaxation times at low temperatures, the MCTDH method is recommended over the QCT approach for conditions when vibrational deactivation events are rare.

The quantum-mechanical VT rates for transitions ν = 1 → ν = x, {x = 0, 2, 3, …, 36}, and the QCT VT rates for all the possible initial rovibrational states in the temperature range from 1000 K–10 000 K are available upon request. The inclusion of excited vibrational states and corresponding single- and multi-state quantum transitions is important for constructing reliable models of thermodynamics at temperatures relevant to hypersonic flows. Despite the fact that the quantum-mechanical simulations in the present work are limited to ν = 1, the MCTDH approach seems to be a valuable tool for resolving the collisional dynamics of the entire rovibrational ensemble.

Using our models, we obtain very good agreement of calculated, experimental, and empirical data over the entire temperature range. It is now possible to improve the data obtained by the QCT method so as to bridge the gap between low (1000 K) and high (10 000 K) temperatures, yielding very accurate data for future aerodynamics simulations.

I.U. and R.H. thank Hans-Dieter Meyer for helpful discussions and for providing the MCTDH package. This work has been partially supported by the Air Force Office of Scientific Research through Grant No. FA9550-12-1-0483. I.U. acknowledges the Alexander von Humboldt Foundation, Germany, for support through a Feodor Lynen Fellowship. This work used HPC resources (Stampede) at the Texas Advanced Computing Center (TACC) at the University of Texas at Austin through an XSEDE allocation under award No. TG-CTS090079.