A hypersonic vehicle traveling at a high speed disrupts the distribution of internal states in the ambient flow and introduces a nonequilibrium distribution in the post-shock conditions. We investigate the vibrational relaxation in diatom-atom collisions in the range of temperatures between 1000 and 10 000 K by comparing results of extensive fully quantum-mechanical and quasi-classical simulations with available experimental data. The present paper simulates the interaction of molecular oxygen with argon as the first step in developing the aerothermodynamics models based on first principles. We devise a routine to standardize such calculations also for other scattering systems. Our results demonstrate very good agreement of vibrational relaxation time, derived from quantum-mechanical calculations with the experimental measurements conducted in shock tube facilities. At the same time, the quasi-classical simulations fail to accurately predict rates of vibrationally inelastic transitions at temperatures lower than 3000 K. This observation and the computational cost of adopted methods suggest that the next generation of high fidelity thermochemical models should be a combination of quantum and quasi-classical approaches.

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.13 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 N2 and O2 are the most predominant compounds. Here, we choose to investigate the O2 + 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 O2 + O and O2 + O2. 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 O2–Ar scattering in a wide range of translational temperatures,

O 2 ( ν ) + Ar O 2 ( ν ) + Ar .
(1)

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,1 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 earlier21,22 at low temperatures do not arise here because our mixtures consist of dilute O2 which necessarily suppresses VER. Meanwhile, Hase and co-workers23 found agreement between theory and experiment for VER in N2 + C6F6 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 function24 with fitted parameters for the O2 potential energy curve25,26 and the Buckingham potential27 for the O–Ar interaction.28 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.31 The parameter values are given in Table I.

V r O a , O b , r O a , Ar , r O b , Ar = V O 2 r O a , O b + i = a , b V O i , Ar r O i , Ar
(2a)
V O 2 r O a , O b = D 1 exp a ( r O a , O b r e ) 2 + c a 3 ( r O a , O b r e ) 3 exp 2 a ( r O a , O b r e ) 1 + a b ( r O a , O b r e )
(2b)
V O i , Ar r O i , Ar = A exp B r O i , Ar C exp 1 . 28 r m ( r O i , Ar 1 . 0 ) 2 / r O i , Ar 6 .
(2c)
TABLE I.

Potential parameters used in the dynamics study. Parameters for the O2-potential are taken from Refs. 25 and 26. Parameters for the Oi–Ar potential are taken from Ref. 28.

VO2(rOa,Ob) VOi,Ar(rOi,Ar)
Parameter Value Unit Parameter Value Unit
a  2.653 645  Å−1  A  209.935 7  Eh 
b  1.670 039    B  2.016 90  a 0 1  
c  6.141 914 × 10−2    C  43.851 40  E h a 0 6  
D  4.206 738 × 104  cm−1  rm  6.427 00  a0 
re  1.207 520  Å       
VO2(rOa,Ob) VOi,Ar(rOi,Ar)
Parameter Value Unit Parameter Value Unit
a  2.653 645  Å−1  A  209.935 7  Eh 
b  1.670 039    B  2.016 90  a 0 1  
c  6.141 914 × 10−2    C  43.851 40  E h a 0 6  
D  4.206 738 × 104  cm−1  rm  6.427 00  a0 
re  1.207 520  Å       

The zero-point energy for the O2–Ar system using this PES is 788 cm−1 [exp EZPE = 787.4 cm−1 (Ref. 32)], and the first vibrationally excited level lies at Eν1 = 1556 cm−1 [exp Eν1 = 1580 cm−1 (Ref. 33)], as calculated using the MCTDH method.34 

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 O2–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 O2–Ar interaction which is represented in Eq. (2) only through the sum of distinct O–Ar interactions. The ensuing shift to lower r2 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.

FIG. 1.

Comparison between the model PES of Eq. (2) with CCSD calculations, for six different impact angles, as indicated in the legend as a function of the distance r2 between Ar and the O2 center of mass.

FIG. 1.

Comparison between the model PES of Eq. (2) with CCSD calculations, for six different impact angles, as indicated in the legend as a function of the distance r2 between Ar and the O2 center of mass.

Close modal

The implementation of the MCTDH method16 by the group of H.-D. Meyer34 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)

Ψ ( Q 1 , , Q f , t ) = j 1 = 1 n 1 j f = 1 n f A j 1 j f ( t ) κ = 1 f ϕ j κ ( κ ) ( Q κ , t )
(3)

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 Aj1jf(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 HR as a sum of products of single-particle operators,

H = κ = 1 f h κ + H R
(4)
H R = r = 1 s c r κ = 1 f h r ( κ )
(5)

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.38 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

i W ( Q ) = i η ( Q Q c ) b θ ( Q Q c ) ,
(6)

where in our calculations, the order of the CAPs is three (b = 3), and the grid point Qc denotes the position at which the CAP is switched on, with θ(QQc) 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 r2 (O–O distance), translational coordinate r1 (distance of Ar to the center of mass of the O2 molecule), and the angle θ between r1 and r2. 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):

H ˆ = 1 2 μ r 1 2 r 1 2 1 2 μ r 2 2 r 2 2 + j ˆ 2 2 μ r 2 r 2 2 + j ( j + 1 ) 2 μ r 2 r 2 2 + 1 2 μ r 1 r 1 2 J ( J + 1 ) 2 K 2 + j ˆ 2 + V ( r 1 , r 2 , θ ) ,
(7)
j ˆ = 1 sin θ θ sin θ θ K 2 sin 2 θ .
(8)

The reduced masses μr1 and μr2 are the reduced masses of the O2–Ar and O2-systems, respectively. The PES V(r1, r2, θ) 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 Pif(E) the method developed by Tannor and Weeks41 is most suitable, where the S-matrix elements are obtained from the cross-correlation function Cif(t). Final-state averaged transition probabilities can then be obtained by summation over the final states f

P i ( E ) = f P i f ( E ) .
(9)

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 Pi(E) as

P i ( E ) = 2 π | Δ ( E ) | 2 Re 0 g ( τ ) e i E τ d τ ,
(10)
g ( τ ) = 0 d t Ψ ( t + τ ) W Ψ ( t ) .
(11)

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 r1 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),

Δ ( E ) = Ψ dist | Ψ 0 = μ r 1 2 π 0 exp i r 1 k ( x ) d x k ( r 1 ) χ 0 ( r 1 ) d r 1 ,
(12)

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

k ( r 1 ) = 2 μ r 1 ( E V ̄ ν ( r 1 ) )
(13)

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 Ecol (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ν(Ecol) are then obtained as42 

σ ν j ν ( E col ) = π ħ 2 2 μ r 1 E col J = 0 J max K = min ( J , j ) K = min ( J , j ) ( 2 J + 1 ) ( 2 j + 1 ) P ν j ν J K ( E col ) .
(14)

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,43 or the quantum analog of the capture model,44 

P ν j ν J K ( E ) = P ν j ν 0 , K ( E E shift ) ,
(15)
E shift ( r 1 ) = ν j | V | ν j + ħ 2 2 μ r 1 r 1 , barr 2 J ( J + 1 ) .
(16)

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

P ν j ν J K ( E col ) = P ν j ν J , 0 ( E col ) .
(17)

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

σ ν ν , T ( E col ) = j = 1 j max p ( j , T ) σ ν j ν ( E col ) ,
(18)
p ( j , T ) = ( 2 j + 1 ) exp β E rot Q rot ( T ) ,
(19)
Q rot ( T ) = j = 1 j max ( 2 j + 1 ) exp β E rot .
(20)

Here, β = (kBT)−1 with kB the Boltzmann constant. Due to the O2 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

k ν ν ( T ) = g s 0 E max 8 π μ r 1 β 3 e β E col σ ν ν , T ( E col ) E col d E col
(21)

with gs = 1/3 the factor of spin degeneracy, and k in units of cm3/s. The value of gs comes from the ratio of the degeneracies of PES and reactants.29 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,

k ν ν = k ν ν exp β E v i b ,
(22)

where Evib 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,

P ν j ν ( E col ) = J J max ( 2 J + 1 ) ( 2 min ( j , J ) + 1 ) P ν j ν J ( E col ) J J max ( 2 J + 1 ) ( 2 min ( j , J ) + 1 ) .
(23)

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

P ν ν , T ( E col ) = j = 1 j max p ( j , T ) P ν j ν ( E col ) ,
(24)
P ν ν ( T ) = g s 0 E max 8 π μ r 1 β 3 e β E col P ν ν , T ( E col ) E col d E col / N ,
(25)
N = 0 E max 8 π μ r 1 β 3 e β E col E col d E col .
(26)

The QCT method is used to reduce the computational cost of generating the O2–Ar vibrationally resolved cross sections. Details of this approach are given elsewhere.45 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 11th 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 × 103 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−4%. 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 Ecol is defined as follows:

P ν j ν j E col = 2 b max 2 0 b max N ν j ν j N b d b ,
(27)

where bmax 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:

σ E col , ν , j = π b max 2 P E col , ν , j ,
(28)

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 O2–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,46 Camac14 and recent measurements by Hanson and co-workers.15 The vibrational relaxation parameter τ, is expressed in Landau-Teller theory as1 

p τ ν ν ( T ) = 1 β k ν ν ( T ) 1 e β E vib
(29)

with pressure p and energy difference Evib 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., Evib ≫ kBT), Eq. (29) simplifies to

p τ ν ν ( T ) = 1 β k ν ν ( T ) .
(30)

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 N2–N (Ref. 10) and O2–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.47 In the present work, the O2–Ar vibrational relaxation time is obtained using the e-folding method48 from the O2–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 Ef 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 White46 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.49 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 nr1 = 512 and nr2 = 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 O2. The CAPs start at 7.50 a0 (r1), and 5.00 a0 (r2), with strength 0 . 004 a 0 1 , and order three.

The initial wavepacket Ψ0 is represented by a Gaussian SPF in the translational DOF r1 (with center at Qr1 = 6.50 a0, width ΔQr1 = 0.10 a0, and initial momentum p0 from −20 a.u. to −105 a.u.). The SPF for the bound DOF r2 is initially setup as a vibrational eigenstate of the O2 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, r1 and r2, 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 MCTDH49 to address H+H2 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 r1,barr = 3.8 Å.

The calculated transition probabilities for vibrational deactivation ν = 1 → ν′ = 0, J = 0, and selected rotational momenta of the O2 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.

FIG. 2.

(a) Final-state averaged transition probabilities obtained from MCTDH, for v = 1 → ν = 0, with the body-fixed and projected rotational quantum numbers J = K = 0 and initial j as specified. (b) Final-state averaged transition probabilities, QCT, for v = 1 → ν = 0 and initial j as specified.

FIG. 2.

(a) Final-state averaged transition probabilities obtained from MCTDH, for v = 1 → ν = 0, with the body-fixed and projected rotational quantum numbers J = K = 0 and initial j as specified. (b) Final-state averaged transition probabilities, QCT, for v = 1 → ν = 0 and initial j as specified.

Close modal

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 Å2, with one exception: in the quantum calculations, the lowest state j = 1 reaches a scattering cross section of almost 7 Å2 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.

FIG. 3.

(a) Scattering cross sections as calculated with MCTDH, for v = 1 → ν = 0 and initial j as specified. (b) Scattering cross sections calculated with the QCT method, for v = 1 → ν = 0 and initial j as specified.

FIG. 3.

(a) Scattering cross sections as calculated with MCTDH, for v = 1 → ν = 0 and initial j as specified. (b) Scattering cross sections calculated with the QCT method, for v = 1 → ν = 0 and initial j as specified.

Close modal

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 Å2. The QCT cross sections show a slightly more prominent rise and level off between 3.0-4.0 Å2, where with higher temperature the differences between the cross sections are more pronounced than for the MCTDH results.

FIG. 4.

(a) Rotationally averaged scattering cross sections as calculated with MCTDH, for v = 1 → ν = 0 and rotational temperature T as specified. (b) Rotationally averaged scattering cross sections as calculated wit the QCT method, for v = 1 → ν = 0 and rotational temperature T as specified.

FIG. 4.

(a) Rotationally averaged scattering cross sections as calculated with MCTDH, for v = 1 → ν = 0 and rotational temperature T as specified. (b) Rotationally averaged scattering cross sections as calculated wit the QCT method, for v = 1 → ν = 0 and rotational temperature T as specified.

Close modal

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}.

FIG. 5.

Vibration-translation rates from MCTDH and QCT calculations, for vibrational deactivation from ν = 1 → ν = x, {x = 0, 2, 3, 5, 10, 15}.

FIG. 5.

Vibration-translation rates from MCTDH and QCT calculations, for vibrational deactivation from ν = 1 → ν = x, {x = 0, 2, 3, 5, 10, 15}.

Close modal

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 O2–Ar collisions were reported for the first time by Camac,14 where the ultraviolet absorption technique was applied to collect data in O2–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 O2–Ar molecular system was revisited by Hanson and co-workers15 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 O2–Ar system.46 Symbols represent the experimental data from Ref. 14, for different admixtures of O2 in Ar, and from Ref. 15.

FIG. 6.

Vibrational relaxation times, experiment data, and MCTDH and QCT methods. The theoretical relaxation times are obtained using different approximations (MCTDH (a) and QCT (a) using Eq. (29), MCTDH (b) and QCT (b) using Eq. (30), and QCT (c) using master equation calculations). The Millikan-White data are generated according to the empirical model in Ref. 46. The experimental data included were obtained for different gas mixtures (Exp(a): 0.25% O2 and 99.75% Ar, and 1.0%O2 and 99.0% Ar (Ref. 14); Exp(b): 21.5% O2 and 78.5% Ar (Ref. 14); Exp(c): 1.0% O2 and 99.0% Ar (Ref. 15)).

FIG. 6.

Vibrational relaxation times, experiment data, and MCTDH and QCT methods. The theoretical relaxation times are obtained using different approximations (MCTDH (a) and QCT (a) using Eq. (29), MCTDH (b) and QCT (b) using Eq. (30), and QCT (c) using master equation calculations). The Millikan-White data are generated according to the empirical model in Ref. 46. The experimental data included were obtained for different gas mixtures (Exp(a): 0.25% O2 and 99.75% Ar, and 1.0%O2 and 99.0% Ar (Ref. 14); Exp(b): 21.5% O2 and 78.5% Ar (Ref. 14); Exp(c): 1.0% O2 and 99.0% Ar (Ref. 15)).

Close modal

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.15 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 O2–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 O2–N2 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 O2 electronic ground state, 3Σg. The other low-lying electronic states of O2 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 H2 + D2 scattering system, J-K coupling has been shown to enhance the reactivity.51 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 O2–Ar deactivation from the first excited state, P1→0. Using QCT methodology, P1→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

P 1 0 Z = K 1 0
(31)

where Z is the collision frequency, and K1→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 P1→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 bmax increases toward higher temperatures, which can be understood from the larger O2 interatomic distance in excited rotational states. The difference between the QCT value of bmax 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 bmax obtained in the QCT calculations in comparison with experiment, is indicative of the relative accuracy of the underlying PES used in the present studies.

TABLE II.

Probability of O2–Ar deactivation from O 2 ν = 1 .

T, K 1000 1200 1800 2400 3000 4000 6000 8000 10 000
Camac14   …  1.6 × 10−6  1.2 × 10−5  5 × 10−5  1.6 × 10−4  6.3 × 10−4  3.5 × 10−3  …  … 
MCTDH  9.28 × 10−8  4.40 × 10−7  7.10 × 10−6  3.50 × 10−5  1.05 × 10−4  3.72 × 10−4  1.69 × 10−3  3.99 × 10−3  6.77 × 10−3 
QCT  3.24 × 10−9  4.54 × 10−8  4.22 × 10−6  4.39 × 10−5  1.84 × 10−4  7.95 × 10−4  3.53 × 10−3  7.41 × 10−3  1.14 × 10−2 
bmax, Å (QCT)  3.68  3.70  3.72  3.73  3.74  3.76  3.79  3.81  3.82 
T, K 1000 1200 1800 2400 3000 4000 6000 8000 10 000
Camac14   …  1.6 × 10−6  1.2 × 10−5  5 × 10−5  1.6 × 10−4  6.3 × 10−4  3.5 × 10−3  …  … 
MCTDH  9.28 × 10−8  4.40 × 10−7  7.10 × 10−6  3.50 × 10−5  1.05 × 10−4  3.72 × 10−4  1.69 × 10−3  3.99 × 10−3  6.77 × 10−3 
QCT  3.24 × 10−9  4.54 × 10−8  4.22 × 10−6  4.39 × 10−5  1.84 × 10−4  7.95 × 10−4  3.53 × 10−3  7.41 × 10−3  1.14 × 10−2 
bmax, Å (QCT)  3.68  3.70  3.72  3.73  3.74  3.76  3.79  3.81  3.82 

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 O2–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 O2 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 bmax. 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.

1.
W. G.
Vincenti
,
J.
Charles
, and
H.
Kruger
,
Introduction to Physical Gas Dynamics
, Reprint 1977 ed. (
Robert E. Krieger Publishing Company
,
Huntington, New York
,
1967
).
2.
C.
Park
,
Nonequilibrium Hypersonic Aerothermodynamics
(
Wiley
,
New York
,
1989
).
3.
J. G.
Kim
and
I. D.
Boyd
, “
State-resolved master equation analysis of thermochemical nonequilibrium of nitrogen
,”
Chem. Phys.
415
,
237
246
(
2013
).
4.
I.
Armenise
and
E. V.
Kustova
, “
On different contributions to the heat flux and diffusion in non-equilibrium flows
,”
Chem. Phys.
428
,
90
104
(
2014
).
5.
F.
Esposito
and
M.
Capitelli
, “
Quasiclassical trajectory calculations of vibrationally specific dissociation cross sections and rate constants for the reaction O + O2 = 3O
,”
Chem. Phys. Lett.
364
,
180
187
(
2002
).
6.
F.
Esposito
,
I.
Armenise
,
G.
Capitta
, and
M.
Capitelli
, “
O + O2 state-to-state vibrational relaxation and dissociation rates based on quasiclassical calculations
,”
Chem. Phys.
351
,
91
98
(
2008
).
7.
R.
Jaffe
,
D.
Schwenke
, and
G.
Chaban
, “Theoretical analysis of N2 collisional dissociation and rotation-vibration energy transfer,” AIAA Paper 2009-1569, 2009.
8.
J.
Bender
,
I.
Nompelis
, and
P.
Valentini
, “Quasiclassical trajectory analysis of the N2 + N2 reaction using a new ab initio potential energy surface,” AIAA Paper 2014-2964, 2014.
9.
D. A.
Andrienko
and
I. D.
Boyd
, “
High fidelity modeling of thermal relaxation and dissociation of oxygen
,”
Phys. Fluids
27
,
116101
(
2015
).
10.
M.
Panesi
,
R. L.
Jaffe
,
D. W.
Schwenke
, and
T. E.
Magin
, “
Rovibrational internal energy transfer and dissociation of N 2 ( 1 Σ g + ) N ( 4 S u ) system in hypersonic flows
,”
J. Chem. Phys.
138
,
044312
(
2013
).
11.
M.
Kulakhmetov
,
M.
Gallis
, and
A.
Alexeenko
, “
Effect of O2 + O ab initio and morse additive pairwise potentials on dissociation and relaxation rates for nonequilibrium flow calculations
,”
Phys. Fluids
27
,
087104
(
2015
).
12.
J. D.
Bender
,
P.
Valentini
,
I.
Nompelis
,
Y.
Paukku
,
Z.
Varga
,
D. G.
Truhlar
,
T.
Schwartzentruber
, and
G. V.
Candler
, “
An improved potential energy surface and multi-temperature quasiclassical trajectory calculations ofN2 + N2 dissociation reactions
,”
J. Chem. Phys.
143
,
054304
(
2015
).
13.
M.
Capitelli
,
C. M.
Ferreira
,
B. F.
Gordiets
, and
A. I.
Osipov
,
Plasma Kinetics in Atmospheric Gases
(
Springer Science & Business Media
,
2013
), Vol.
31
.
14.
M.
Camac
, “
O2 vibration relaxation in oxygen-argon mixtures
,”
J. Chem. Phys.
34
,
448
459
(
1961
).
15.
K. G.
Owen
,
D. F.
Davidson
, and
R. K.
Hanson
, “
Oxygen vibrational relaxation times: Shock tube/laser absorption measurements
,”
J. Thermophys. Heat Transfer
(published online).
16.
H.-D.
Meyer
,
U.
Manthe
, and
L. S.
Cederbaum
, “
The multi-configurational time-dependent Hartree approach
,”
Chem. Phys. Lett.
165
,
73
78
(
1990
).
17.
H.-D.
Meyer
, “
Studying molecular quantum dynamics with the multiconfiguration time-dependent Hartree method
,”
WIREs: Comput. Mol. Sci.
2
,
351
374
(
2012
).
18.
Multidimensional Quantum Dynamics: MCTDH Theory and Applications
, edited by
H.-D.
Meyer
,
F.
Gatti
, and
G. A.
Worth
(
Wiley-VCH
,
Weinheim
,
2009
).
19.
F.
Esposito
,
I.
Armenise
, and
M.
Capitelli
, “
N–N2 state to state vibrational-relaxation and dissociation rates based on quasiclassical calculations
,”
Chem. Phys.
331
,
1
8
(
2006
).
20.
D.
Andrienko
and
I. D.
Boyd
, “
Investigation of oxygen vibrational relaxation by quasi-classical trajectory method
,”
Chem. Phys.
459
,
1
13
(
2015
).
21.
K. F.
Everitt
and
J. L.
Skinner
, “
Vibrational energy relaxation of oxygen in liquid mixtures with argon
,”
J. Chem. Phys.
110
,
4467
4470
(
1999
).
22.
X.
Tong
,
T.
Nagy
,
J. Y.
Reyes
,
M.
Germann
,
M.
Meuwly
, and
S.
Willitsch
, “
State-selected ionmolecule reactions with coulomb-crystallized molecular ions in traps
,”
Chem. Phys. Lett.
547
,
1
8
(
2012
).
23.
A. K.
Paul
,
S. C.
Kohale
,
S.
Pratihar
,
R.
Sun
,
S. W.
North
, and
W. L.
Hase
, “
A unified model for simulating liquid and gas phase, intermolecular energy transfer: N2 + C6F6 collisions
,”
J. Chem. Phys.
140
,
194103
(
2014
).
24.
D.
Steele
and
E. R.
Lippincott
, “
Comparative study of empirical internuclear potential functions
,”
Rev. Mod. Phys.
34
,
239
251
(
1962
).
25.
J. T.
Vanderslice
,
E. A.
Mason
,
W. G.
Maisch
, and
E. R.
Lippincott
, “
Ground state of hydrogen by the Rydberg-Klein-Rees method
,”
J. Mol. Spectrosc.
3
,
17
29
(
1959
).
26.
J. G.
Kim
and
I. D.
Boyd
, “Thermochemical nonequilibrium modeling of electronically excited molecular oxygen,” 11th AIAA Aviation/ASME Joint Thermophysics and Heat Transfer Conference, 2014.
27.
A.
Gross
and
G. D.
Billing
, “
Rate constants for ozone formation and for isotopic exchange reactions
,”
Chem. Phys.
173
,
393
406
(
1993
).
28.
G. J.
Kroes
and
R. P. H.
Rettschnick
, “
Vibrational relaxation of glyoxal in collisions with He and Ar
,”
Chem. Phys.
156
,
293
307
(
1991
).
29.
A.
Gross
and
G. D.
Billing
, “
Rate constants for ozone formation
,”
Chem. Phys.
187
,
329
335
(
1994
).
30.
K.
Mizobata
, “An analysis of quasiclassical molecular collisions and rate processes for coupled vibration-dissociation and recombination,” AIAA Paper No. 97–0132, 1997.
31.
A. J. C.
Varandas
, “
Four-atom bimolecular reactions with relevance in environmental chemistry: Theoretical work
,”
Int. Rev. Phys. Chem.
19
,
199
245
(
2000
).
32.
K. K.
Irikura
, “
Experimental vibrational zero-point energies: Diatomic molecules
,”
J. Phys. Chem. Ref. Data
36
,
389
397
(
2007
).
33.
K. P.
Huber
and
G.
Herzberg
,
Molecular Spectra and Molecular Structure. IV. Constants of Diatomic Molecules
(
Van Nostrand Reinhold Co.
,
1979
).
34.
M. H.
Beck
,
A.
Jäckle
,
G. A.
Worth
, and
H.-D.
Meyer
, “
The multiconfiguration time-dependent Hartree (MCTDH) method: A highly efficient algorithm for propagating wavepackets
,”
Phys. Rep.
324
,
1
105
(
2000
).
35.
M. W.
Schmidt
,
K. K.
Baldridge
,
J. A.
Boatz
,
S. T.
Elbert
,
M. S.
Gordon
,
J. H.
Jensen
,
S.
Koseki
,
N.
Matsunaga
,
K. A.
Nguyen
,
S.
Su
,
T. L.
Windus
,
M.
Dupuis
, and
J. A.
Montgomery
, “
General atomic and molecular electronic structure system
,”
J. Comput. Chem.
14
,
1347
1363
(
1993
).
36.
J. L.
Bentz
,
R. M.
Olson
,
M. S.
Gordon
,
M. W.
Schmidt
, and
R. A.
Kendall
,
Comput. Phys. Commun.
176
,
589
600
(
2007
).
37.
M. M.
Francl
,
W. J.
Pietro
,
W. J.
Hehre
,
J. S.
Binkley
,
M. S.
Gordon
,
D. J.
DeFrees
, and
J. A.
Pople
,
J. Chem. Phys.
77
,
3654
3665
(
1983
).
38.
A.
Jäckle
and
H.-D.
Meyer
, “
Product representation of potential energy surfaces
,”
J. Chem. Phys.
104
,
7974
(
1996
).
39.
U. V.
Riss
and
H.-D.
Meyer
, “
Calculation of resonance energies and widths using the complex absorbing potential method
,”
J. Phys. B
26
,
4503
(
1993
).
40.
U. V.
Riss
and
H.-D.
Meyer
, “
Investigation on the reflection and transmission properties of complex absorbing potentials
,”
J. Chem. Phys.
105
,
1409
(
1996
).
41.
D. J.
Tannor
and
D. E.
Weeks
, “
Wave packet correlation function formulation of scattering theory: The quantum analog of classical S-matrix theory
,”
J. Chem. Phys.
98
,
3884
3893
(
1993
).
42.
N.
Bulut
,
J. F.
Castillo
,
F. J.
Aoiz
, and
L.
Bañares
, “
Real wave packet and quasiclassical trajectory studies of the H+ + LiH reaction
,”
Phys. Chem. Chem. Phys.
10
,
821
827
(
2008
).
43.
Q.
Sun
,
J. M.
Bowman
,
G. C.
Schatz
,
J. R.
Sharp
, and
J. N. L.
Connor
, “
Reduced-dimensionality quantum calculations of the thermal rate coefficient for the Cl + HCl → ClH + Cl reaction: Comparison with centrifugal-sudden distorted wave, coupled channel hyperspherical, and experimental results
,”
J. Chem. Phys.
92
,
1677
1686
(
1990
).
44.
S. K.
Gray
,
E. M.
Goldfield
,
G. C.
Schatz
, and
G. G.
Balint-Kurti
, “
Helicity decoupled quantum dynamics and capture model cross sections and rate constants for O(1D) + H2 → OH + H
,”
Phys. Chem. Chem. Phys.
1
,
1141
1148
(
1999
).
45.
D. G.
Truhlar
and
J. T.
Muckerman
, “
Reactive scattering cross sections III: Quasiclassical and semiclassical methods
,” in
Atom-Molecule Collision Theory
(
Springer
,
1979
), pp.
505
566
.
46.
R. C.
Millikan
and
D. R.
White
, “
Systematics of vibrational relaxation
,”
J. Chem. Phys.
39
,
3209
3213
(
1963
).
47.
D. A.
Andrienko
and
I. D.
Boyd
, “High fidelity modeling of thermal relaxation and dissociation of oxygen,” AIAA Paper 2016-0736, 2016.
48.
C.
Park
, “
Rotational relaxation of N2 behind a strong shock wave
,”
J. Thermophys. Heat Transfer
18
,
527
533
(
2004
).
49.
G. A.
Worth
,
M. H.
Beck
,
A.
Jäckle
, and
H.-D.
Meyer
, The MCTDH Package, version 8.2, 2000; H.-D. Meyer, version 8.3, 2002; version 8.4, 2007; current version: 8.4.11, 2015, see http://mctdh.uni-hd.de/.
50.
E.
Garcia
,
A.
Kurnosov
,
A.
Lagana
,
F.
Pirani
,
M.
Bartolomei
, and
M.
Cacciatore
, “
Efficiency of collisional O2 + N2 vibrational energy exchange
,”
J. Phys. Chem. B
120
,
1476
1485
(
2016
).
51.
S.
Sukiasyan
and
H.-D.
Meyer
, “
On the effect of initial rotation on reactivity. A multi-configuration time-dependent Hartree (MCTDH) wave packet propagation study on the H + D2 and D + H2 reactive scattering systems
,”
J. Phys. Chem. A
105
,
351
374
(
2001
).