Vibrational energy relaxation (VER) of diatomics following collisions with the surrounding medium is an important elementary process for modeling high-temperature gas flow. VER is characterized by two parameters: the vibrational relaxation time τvib and the state relaxation rates. Here the vibrational relaxation of CO(ν=0ν=1) in Ar is considered for validating a computational approach to determine the vibrational relaxation time parameter (pτvib) using an accurate, fully dimensional potential energy surface. For lower temperatures, comparison with experimental data shows very good agreement whereas at higher temperatures (up to 25 000 K), comparisons with an empirically modified model due to Park confirm its validity for CO in Ar. Additionally, the calculations provide insight into the importance of Δν>1 transitions that are ignored in typical applications of the Landau–Teller framework.

The relaxation of vibrational energy from an energized molecule to its environment is a key process for energy redistribution in gas and condensed phase systems.1–9 Specifically, the vibrational relaxation of diatomic molecules during collisions with the surrounding gas is an important elementary process in hypersonic flow typical of spacecraft entering an atmosphere.10–16 The key parameters for modeling these systems are the rate coefficients describing the collision of diatomic molecules with other atoms and molecules in the atmosphere. The vibrational relaxation process is essential for characterizing the energy flow from vibrational to rotational and translational degrees of freedom and chemical kinetics can depend on the distribution of energy between these modes.

CO in Ar is a suitable proxy for a generic investigation of vibrational energy relaxation (VER) at high temperatures because experimental data up to 3000 K are available for comparison17 from which an extrapolation model due to Millikan and White (MW) based on Landau-Teller (LT) theory was developed.10 It is also of interest to note that the atmosphere of Mars contains 3% argon and significant amounts of CO which make Ar–CO relevant for space exploration missions.12 For Ar–CO, a recent high-level spectroscopically accurate (CCSD(T)-F12b/aug-cc-pV5Z) potential energy surface (PES) is available.18 The ab initio grid covers CO diatomic distances r from 1.0 Å to 1.35 Å. This interval includes the classical turning points up to ν=4 for CO.19 In the present work, this PES was employed to generate a comprehensive grid of energies (3861 points) using a reproducing kernel Hilbert space (RKHS)20 representation, following the description of Ref. 21. One advantage of an RKHS is that the necessary derivatives can be computed efficiently and analytically.22,23 The general approach for quantifying the vibrational energy relaxation Evib is to start from the LT formalism,24 

dEvibdt=1τvib(EvibEvib),
(1)

where Evib is the equilibrium thermal energy of the oscillator in the bath and τvib is the vibrational relaxation time. This quantity is commonly reported as the vibrational relaxation time parameter pτvib, where p is the pressure and

pτvib=kBTk10(T)(1eω/kBT).
(2)

In Eq. (2), k10(T) corresponds to the rate coefficient for the vibrational transition 10 at temperature T, kB is the Boltzmann constant, is the reduced Planck constant, and ω is the harmonic frequency of the diatomic molecule. The function pτvib(T) initially decreases with temperature because the population of Ar atoms with sufficient energy to efficiently couple to the CO molecule increases with increasing T. At some temperature, the efficiency to move vibrational energy in a collision is maximal and pτvib increases as it is limited by the time between collisions which turns out to be at T15 000 K. However, for Ar–CO the experimental data only extend up to 3000 K. Computing pτvib requires the calculation of k10(T), and at high temperatures quasi-classical trajectory (QCT) calculations are a suitable and meaningful way to do so.

QCT simulations have provided valuable insight into processes such as rotational excitation in the charge exchange reaction N2+ + N2 → N2 + N2+25 or the state-to-state dynamics of the Br+H2 reaction.26 Here, QCT calculations are used to determine thermal rate coefficients. The dynamics of the system is followed by propagating the equations of motions for given initial conditions using the Velocity-Verlet algorithm27 with a time step of 0.02 fs to ensure the conservation of total energy. Initial conditions for CO were generated from a WKB-quantized periodic orbit28 of the corresponding rotating Morse oscillator.29 The symmetry axis of CO, the axis of its rotation, and the angular momentum were randomly oriented and the importance sampling Monte Carlo scheme was used.30 For the impact parameter b, stratified sampling was employed.31,32 At a given temperature, b was sampled on the interval [0;bmax] for five different values of bmax. The 2bdbbmax2 probability function was employed for sampling b in each interval. Finally, all trajectories were analysed jointly by grouping them in non-overlapping intervals of b with a weight

Vk=(bk+)2(bk)2bmax2,
(3)

where b and b+ are the minimum and maximum values of b, respectively. Then the reaction probability is

Preac=kVkNkNk,tot,
(4)

where k labels one particular stratum (interval), and Nk and Nk,tot are the number of reactive trajectories (or event of interest, e.g., Δν=1 for vibrational relaxation) and the total number of trajectories in stratum k. The complex is considered to be dissociated when the Ar and CO moieties are separated by more than the initial separation and the final state is determined from subtracting the translational and rotational energies from the total energy and using semiclassical quantization to determine an approximate vibrational quantum number ν. Starting from ν= 1 relaxation to ν= 0 was considered to occur if ν<0.5, and similarly when analyzing relaxation from higher excited states. Using such a criterion should be sufficient here as no final state analysis is sought but only the information whether in a particular trajectory vibrational relaxation has occurred is of interest. This conventional binning procedure has recently been compared33 for vibrationally induced photodissociation of HSO3Cl with more refined procedures based on Gaussian binning34 and found to yield satisfactory agreement. The rate coefficient from state i to state l is then determined from

kil(T)=8kBTπμπbmax2Preac.
(5)

For most of the present work, the initial vibrational quantum number was ν=1, and therefore the vibrational temperature is undefined. Several tests were performed by choosing bmax at a given T. The Nk,tot at each temperature varies from 9000 to 25 000. QCT calculations between 4000 K and 25 000 K were performed for nine temperatures. At T<4000 K, almost no trajectory showed vibrational relaxation.

Figure 1 reports the dependence of Δν on the range of impact parameter sampled. For a given Δν, the interval of b increases with temperature. Hence, the effective cross section for changing the vibrational state of CO through collision with Ar increases with T. Transitions with Δν=±1 can occur over a wide range of b, with 889 (Δν=1) and 897 (Δν=+1) events, respectively, at 25 000 K (see Figure 1). Furthermore, the probability of a vibrational transition has no angular preference due to the isotropic long range interaction potential.

FIG. 1.

Relationship between b and the probability to observe a particular change in ν of CO starting from ν=1. Results are reported for three different temperatures: T = 4000 K (red), T = 13 000 K (blue), and T = 25 000 K (green). The maximum value of b used in the sampling procedure at a given T is also shown, bmaxT.

FIG. 1.

Relationship between b and the probability to observe a particular change in ν of CO starting from ν=1. Results are reported for three different temperatures: T = 4000 K (red), T = 13 000 K (blue), and T = 25 000 K (green). The maximum value of b used in the sampling procedure at a given T is also shown, bmaxT.

Close modal

Figure 2 shows the rate coefficients for vibrational transitions of CO (ν=1) after collision with Ar between 4000 and 25 000 K. Below 10 000 K only processes with Δν=±1 are important and even at 25 000 K, Δν=±1 is the most probable process. Over the temperature range of interest, the rates increase by three orders of magnitude for Δν=±1. Another process that could become important at high T is full atomization (Ar+C+O). However, the number of trajectories with dissociated CO (DeCO=11.2 eV) even at 20 000 K is minute (<0.03%). At this temperature, very high rotational levels are populated which create a centrifugal barrier and largely prevent CO dissociation.

FIG. 2.

Rate coefficients for transitions 1ν with ν=0,2,3, and 4. Also the sum of the inelastic rate coefficients (black line) is included.

FIG. 2.

Rate coefficients for transitions 1ν with ν=0,2,3, and 4. Also the sum of the inelastic rate coefficients (black line) is included.

Close modal

At high temperatures, excited vibrational states of CO are populated (e.g., at 10 000 K, p(νCO=2)=13%). Running 10 000 trajectories from the initial state ν=2 at 10 000 K showed that the rate coefficients for Δν>1 are one order of magnitude smaller than Δν=±1 (e.g., k20=1.25× 1012 cm3 molecule−1 s−1 and k23=1.68×1011 cm3 molecule−1 s−1). The results for ν=2 show that the calculations satisfy detailed balance: from k21=1.36×1011 cm3molecule−1 s−1, k12=1.01×1011 cm3 molecule−1 s−1 is obtained which is in excellent agreement with the value for k12=1.09×1011 cm3 molecule−1 s−1 computed directly from the QCT calculations from initial ν=1. For testing the validity of the LT formalism, 10 000 trajectories for ν=6 were run at 10 000 K and 25 000 K. At the lower temperature, the sum of de-excitation coefficients for Δν2 was less than 20% of k65 but as high as 72% at the higher temperature. Hence, Δν±1 is an appropriate assumption at 10 000 K but not at 25 000 K.

Table I shows the rate coefficients k10 and the vibrational relaxation parameter pτvib. Furthermore, at each T the bmax used in the stratified sampling is reported. The computed pτvib from 4000 to 25 000 K are compared with the available data and fits to empirical expressions, see Figure 3. The “MW model” is the expression fitted to the experimental data (measured from 1000 to 3000 K) for Ar–CO. For O2–O2, the MW model faithfully describes available experiments up to 8000 K. With this in mind and accounting for different molecular parameters for the systems, we expect the MW model to be valid for Ar–CO well beyond the experimentally accessible temperature of 3000 K.

TABLE I.

Rate coefficients for the 10 vibrational transition k10 (in cm3 molecule−1 s−1) and vibrational relaxation time pτvib (in atm s). The standard deviation of the rate coefficients Δk10, the maximum value of the impact parameter bmax (in a0) used in the stratified sampling, the number of trajectories Ntot at each temperature, the sum of inelastic rate, and the removal efficiency pτvibrem from Eq. (7) are also included.

TbmaxNtotk10Δk10pτvibν1k1νpτvibrem
4 000 25 000 9.36×1014 4.01×1014 1.09×105 1.36×1013 4.07×106 
5 000 25 000 2.38×1013 7.49×1014 6.23×106 6.73×1013 1.03×106 
6 000 10 000 8.02×1013 1.85×1013 2.55×106 1.96×1012 4.22×107 
7 000 10 000 2.46×1012 4.28×1013 1.09×106 4.91×1012 1.97×107 
10 000 9 000 1.15×1011 1.19×1012 4.46×107 2.43×1011 5.69×108 
13 000 10 000 2.54×1011 1.72×1012 3.31×107 5.75×1011 3.12×108 
15 000 10 000 3.48×1011 2.10×1012 3.17×107 9.42×1011 2.20×108 
20 000 9 000 8.30×1011 4.37×1012 2.30×107 2.13×1010 1.30×108 
25 000 10 000 1.07×1010 4.61×1012 2.74×107 3.10×1010 1.11×108 
TbmaxNtotk10Δk10pτvibν1k1νpτvibrem
4 000 25 000 9.36×1014 4.01×1014 1.09×105 1.36×1013 4.07×106 
5 000 25 000 2.38×1013 7.49×1014 6.23×106 6.73×1013 1.03×106 
6 000 10 000 8.02×1013 1.85×1013 2.55×106 1.96×1012 4.22×107 
7 000 10 000 2.46×1012 4.28×1013 1.09×106 4.91×1012 1.97×107 
10 000 9 000 1.15×1011 1.19×1012 4.46×107 2.43×1011 5.69×108 
13 000 10 000 2.54×1011 1.72×1012 3.31×107 5.75×1011 3.12×108 
15 000 10 000 3.48×1011 2.10×1012 3.17×107 9.42×1011 2.20×108 
20 000 9 000 8.30×1011 4.37×1012 2.30×107 2.13×1010 1.30×108 
25 000 10 000 1.07×1010 4.61×1012 2.74×107 3.10×1010 1.11×108 
FIG. 3.

Comparison of computed pτvib for the collision of Ar–CO (open circles) from the simulations with the following: the MW model10 (blue, “experiment”), a fit12 to the MW model including Park’s empirical correction (black, Eq. (6)), and fitting of the QCT-data to the Park model (black, pτvib=Aeb/T1/3+cT5/2) (red line). At high T the computed data are found to follow the T5/2 behaviour. The green line is the collision time.

FIG. 3.

Comparison of computed pτvib for the collision of Ar–CO (open circles) from the simulations with the following: the MW model10 (blue, “experiment”), a fit12 to the MW model including Park’s empirical correction (black, Eq. (6)), and fitting of the QCT-data to the Park model (black, pτvib=Aeb/T1/3+cT5/2) (red line). At high T the computed data are found to follow the T5/2 behaviour. The green line is the collision time.

Close modal

Up to 10 000 K, the agreement between the QCT-computed and the experimental17 data is very good. For T> 10 000 K, the MW model (blue trace, “experiment”) and the QCT simulation results (open red symbols) differ. For high temperatures, the computed data decay more slowly while the MW model . The collision time, pτcol=(πμkBT)1/2/22πσp,q2 where σp,q=0.5(σp+σq) and σp and σq are the hard sphere diameters of the CO and Ar, respectively,13 was also included (green line Figure 3). For T> 15 000 K, the MW model predicts pτvib<pτcol.

Avoiding such unphysical behaviour (vibrational equilibrium is reached before collisions can occur) was the motivation to include a correction term that maintained pτvib>pτcol up to an arbitrarily high temperature of 50 000 K.12 For Ar–CO, pτvib assumes the form12 

pτvib=e215(T1/30.0302)+πμTkB8(T/50 000)23×1016.
(6)

The pτvib from the QCT simulations and Park’s parametrized correction (black line) agree very well at high T. In fact, Park’s parametrized model is close to a fit up to 15 000 K, see Figure 3. Finally, the QCT-computed pτvib data were also fitted (red solid line) to Park’s model12 pτvib=Aeb/T1/3+cT5/2 and confirm the T5/2-dependence at high temperature. The final fitting parameters are A=5.53×1011 atm s, b = −193.53 K1/3, and c=2.72×1018 atm s K−5/2.

Conventional application of LT theory assumes that collisions behave in a manner analogous to an optical vibrational transition. Within this approximation (truncation of the system-bath coupling at first order), it is considered that including transitions with Δν=1 is sufficient and those with Δν>1 can be ignored. Since simulations track all transitions, limitations of the Δν=1 assumption can be at least qualitatively tested. For CO the vibrational matrix elements (|Rνν|2)—i.e., the integral of the initial and final vibrational wavefunctions over the infrared dipole operator—are available.35 For ν=1ν=0, |R10|2=1.074×102 D2 and for ν=3ν=1, |R31|2=1.134×104 D2 which is almost 2 orders of magnitude lower.

In the QCT simulation however, the raw probabilities can not be directly compared across Δν. The probabilities include a bias for the availability of energy to drive the transition. In other words, the probability for Δν=2 may be low at temperatures for which the average collision energy is below two vibrational quanta. A simplistic way to remove the bias is to rescale the probabilities against an effective temperature, TTν, where Tν is the characteristic vibrational temperature, i.e., Tν=ν(CO)Δν/0.82/3. Plotting pτ1,ν against (TTν) and fitting it to pτ1ν=A(ν)[TTν]d, where d1 provides a means to extract A(ν) which is the likelihood of a collision leading to a change by Δν. Table II summarizes this data and reports the fitted A(ν) values. For the QCT simulations, the A(ν) indicate that Δν=2 is actually an order of magnitude more efficient than Δν=1 when the temperature bias is removed. Hence, for accurate pτ calculations, vibrational matrix elements are not necessarily useful predictors for the collisional vibrational relaxation efficiency.

TABLE II.

The vibrational matrix elements |Rνν|2 and parameters of the function A(ν)=pτ1ν(T)/[TTν]d, where pτ1ν(T) is from the QCT simulations at temperature T.

Δν|R1Δν+1|2 (D2)35 A(ν) (atm s K)d
2.13×102 1.76×102 −1.17 
1.21×104 2.75×101 −1.38 
6.65×107 5.78×102 −1.16 
 4.27×103 −0.88 
 9.80×103 −0.96 
Δν|R1Δν+1|2 (D2)35 A(ν) (atm s K)d
2.13×102 1.76×102 −1.17 
1.21×104 2.75×101 −1.38 
6.65×107 5.78×102 −1.16 
 4.27×103 −0.88 
 9.80×103 −0.96 

Recently, Andrienko and Boyd14 studied the O2–O vibrational relaxation and characterized the vibrational relaxation time by the total removal of population from ν=1 according to

pτvibrem=kBTν1k1ν,
(7)

.These values are also reported in Table I. Eqs. (2) and (7) are equivalent at low temperature (kBTω), where the exponential term in Eq. (2) can be neglected and the sum in Eq. (7) only includes the de-excitation rates. However, at high temperatures all terms are required. Therefore, both pτvib and pτvibrem are useful for modelling conditions at the hypersonic flight regime.

In summary, extrapolation of the QCT calculations using a high-quality PES correctly describe experimental VER data for Ar–CO for T 3000 K. Up to 10 000 K, the empirical MW analysis follows the simulation data which supports the validity of conventional LT theory (Δν=1). Above 10 000 K, Park’s correction is required and the QCT-computed data confirm the T5/2-dependence. The present simulations are consistent with all the available data, and extending the approach to investigate reactive systems such as N + NO N2 + O is the next step. In these systems,36 collisions and subsequent reactions can lead to non-Boltzmann distributions in the product state distributions which will profoundly affect the vibrational relaxation rates and other chemical processes relevant in high-temperature gas flows.

The authors acknowledge stimulating discussions with Professor Tom Schwartzentruber. Part of this work was supported by the United State Department of the Air Force which is gratefully acknowledged (to M.M.). Support by the Swiss National Science Foundation through Grant No. 200021-117810, the NCCR MUST (to M.M.), and the University of Basel is also acknowledged.

1.
M.
Gruebele
and
R.
Bigwood
,
Int. Rev. Phys. Chem.
17
,
91
(
1998
).
2.
M.
Gruebele
and
P.
Wolynes
,
Acc. Chem. Res.
37
,
261
(
2004
).
3.
V.
Wong
and
M.
Gruebele
,
J. Phys. Chem. A
103
,
10083
(
1999
).
4.
A.
Pakoulev
,
Z.
Wang
,
Y.
Pang
, and
D. D.
Dlott
,
Chem. Phys. Lett.
380
,
404
(
2003
).
5.
M.
Gruebele
,
J. Phys. Chem.
100
,
12183
(
1996
).
6.
C.
Lawrence
and
J.
Skinner
,
J. Chem. Phys.
119
,
1623
(
2003
).
7.
S.
Li
,
J.
Schmidt
, and
J.
Skinner
,
J. Chem. Phys.
125
,
244507
(
2006
).
8.
A.
Callegari
,
R.
Pearman
,
S.
Choi
,
P.
Engels
,
H.
Srivastava
,
M.
Gruebele
,
K.
Lehmann
, and
G.
Scoles
,
Mol. Phys.
101
,
551
(
2003
).
9.
T.
Uzer
and
W.
Miller
,
Phys. Rep.
199
,
73
(
1991
).
10.
R. C.
Millikan
and
D. R.
White
,
J. Chem. Phys.
39
,
3209
(
1963
).
11.
C.
Park
,
J. Thermophys. Heat Transfer
7
,
385
(
1993
).
12.
C.
Park
,
J. T.
Howe
,
R. L.
Jaffe
, and
G. V.
Candler
,
J. Thermophys. Heat Transfer
8
,
9
(
1994
).
13.
F.
Thivet
,
M.
Perrin
, and
S.
Candel
,
Phys. Fluids A
3
,
2799
(
1991
).
14.
D.
Andrienko
and
I. D.
Boyd
,
Chem. Phys.
459
,
1
(
2015
).
15.
D. A.
Andrienko
and
I. D.
Boyd
,
Phys. Fluids
27
,
116101
(
2015
).
16.
I. S.
Ulusoy
,
D. A.
Andrienko
,
I. D.
Boyd
, and
R.
Hernandez
,
J. Chem. Phys.
144
,
234311
(
2016
).
17.
W. J.
Hooker
and
R. C.
Millikan
,
J. Chem. Phys.
38
,
214
(
1963
).
18.
Y.
Sumiyoshi
and
Y.
Endo
,
J. Chem. Phys.
142
,
024314
(
2015
).
19.
S. M.
Kirschner
and
J. K.
Watson
,
J. Mol. Spectrosc.
47
,
234
(
1973
).
20.
T.-S.
Ho
and
H.
Rabitz
,
J. Chem. Phys.
104
,
2584
(
1996
).
21.
O. T.
Unke
,
J. C.
Castro-Palacio
,
R. J.
Bemish
, and
M.
Meuwly
,
J. Chem. Phys.
144
,
224307
(
2016
).
22.
J. C.
Castro-Palacio
,
T.
Nagy
,
R. J.
Bemish
, and
M.
Meuwly
,
J. Chem. Phys.
141
,
164319
(
2014
).
23.
J. C.
Castro-Palacio
,
R. J.
Bemish
, and
M.
Meuwly
,
J. Chem. Phys.
142
,
091104
(
2015
).
24.
E. E.
Nikitin
and
J.
Troe
,
Phys. Chem. Chem. Phys.
10
,
1483
(
2008
).
25.
X.
Tong
,
T.
Nagy
,
J. Y.
Reyes
,
M.
Germann
,
M.
Meuwly
, and
S.
Willitsch
,
Chem. Phys. Lett.
547
,
1
(
2012
).
26.
A. N.
Panda
,
D.
Herraez-Aguilar
,
P. G.
Jambrina
,
J.
Aldegunde
,
S. C.
Althorpe
, and
F.
Javier Aoiz
,
Phys. Chem. Chem. Phys.
14
,
13067
(
2012
).
28.
R. L.
Liboff
,
Introductory Quantum Mechanics
(
Addison-Wesley
,
2003
).
29.
R. N.
Porter
,
L. M.
Raff
, and
W. H.
Miller
,
J. Chem. Phys.
63
,
2214
(
1975
).
30.
F.
Daan
and
S.
Berend
,
Understanding Molecular Simulation: From Algorithms to Applications
(
Academic Press
,
San Diego, California
,
2002
).
31.
J. D.
Bender
,
P.
Valentini
,
I.
Nompelis
,
Y.
Paukku
,
Z.
Varga
,
D. G.
Truhlar
,
T.
Schwartzentruber
, and
G. V.
Candler
,
J. Chem. Phys.
143
,
054304
(
2015
).
32.
D. G.
Truhlar
and
J. T.
Muckerman
,
Atom-Molecule Collision Theory
(
Springer
,
1979
), pp.
505
566
.
33.
J.
Yosa Reyes
,
S.
Brickel
,
O. T.
Unke
,
T.
Nagy
, and
M.
Meuwly
,
Phys. Chem. Chem. Phys.
18
,
6780
(
2016
).
34.
L.
Bonnet
and
J.
Rayez
,
Chem. Phys. Lett.
277
,
183
(
1997
).
35.
R. C.
Herman
and
K. E.
Shuler
,
J. Chem. Phys.
22
,
481
(
1954
).
36.
O.
Denis-Alpizar
,
M.
Meuwly
, and
R.
Bemish
,
Phys. Chem. Chem. Phys.
19
,
2392
(
2017
).