Rovibrational energy transfer and dissociation in O₂–O collisions

Simulation of the environment around a spacecraft or hypersonic vehicle during flight in the upper atmosphere presents a challenging problem due to the abundance of physical and chemical effects that take place in the shock-heated air. Among the most important processes are nonequilibrium excitation of internal degrees of freedom, chemical transformations, ionization, and radiation.^1–4 The flow is often found in a state of chemical and thermal nonequilibrium because of the rapid conversion of energy from translational to internal degrees of freedom. A number of empirical models have been proposed to deal with such complex interaction between physical and chemical processes.^1,5–8 However, recent works in aerothermodynamics rely on the gradual rejection of empirical models in favor of accurate state-resolved chemical kinetics based on ab initio methods and quantum calculations.^9–13 

The chemistry of oxygen is one of the most important subjects in the study of cruising flight at hypersonic speed. Among other collisional pairs, the interaction of molecular oxygen with the parent atom is interesting from the physical point of view due to the fast vibrational quenching¹⁴ and a strong chemical effect¹⁵ in the O₃ molecular system. The dynamics of O₂–O collisions have been studied previously theoretically^11,16–20 and experimentally.^21–26 Owing to growing computational resources, the statistical analysis of molecule-atom interaction became computationally tractable. Particularly, the quasi-classical trajectory (QCT)²⁷ method was applied^11,16,17 to obtain the rates of vibrational transitions in O₂–O collisions using an semi-empirical potential energy surface (PES)²⁸ at temperatures observed in hypersonic flows. Recently, the kinetics of the O₂–O molecular system was investigated by means of a system of master equations, using a database of energy transfer rates, obtained for each vibrational state.²⁹ While that work relied on the assumption of trans-rotational equilibrium, the importance of multiquantum transition rates and fast vibrational relaxation at low temperatures in O₂–O collisions was demonstrated in contrast with the O₂–Ar system that has a large repulsive potential. To the authors’ best knowledge, no other master equation studies of the O₃ system using the non-empirical energy transfer rates have been performed.

The assumption of trans-rotational equilibrium is questionable under the strong non-equilibrium conditions in re-entry flows.^30–32 With ever-increasing computational resources, the rovibrational models of energy transfer became available.^32,33 In these simulations, the collisional cross sections are generated for each rovibrational state of a target molecule. Such collisional models of rotational and vibrational relaxation are free of empiricism. Due to the large number of rovibrational states and the cost of trajectory simulations, these models became computationally manageable only recently.

The present paper extends the study of the $O 2 ( X 3 Σ g − ) – O 3 P$ system by introducing a new set of transition rates, generated for each rovibrational state of molecular oxygen using the double many-body potential surface by Varandas and Pais.²⁸ The cross sections of nonreactive, exchange, and dissociation collisions are generated in a wide range of kinetic energies, relevant to hypersonic flows. The corresponding energy transfer rates are employed in a system of master equations that simulates the thermal relaxation of molecular oxygen in a heat bath of oxygen atoms. Each rovibrational state is treated as a separate species. Following this approach, vibrational and rotational relaxation times are obtained by the folding of average internal energy of molecular ensemble. A quasi-steady phase of relaxation is assessed by simulating the dissociation and recombination in the absence of thermal equilibrium. Some aspects of nonequilibrium relaxation, such as the dissociation-rovibrational coupling and the importance of pre-dissociated states, receive special attention. The present paper investigates the O₂–O collision dynamics with the emphasis on the macroscopic and highly averaged parameters that can be utilized in the future for verification of reduced order models.

The paper is organized as follows. The details about the O₂–O molecular structure and trajectory simulations are given in Section II. The system of master equations is described in Section III. Discussion and results, given in Section IV, are divided into three subsections: Section IV A overviews the O₂–O transition rates and internal energy transfer in the absence of dissociation, Section IV B presents the master equation study of rotational and vibrational relaxation, and Section IV C investigates the O₂ thermalization in the presence of depletion mechanisms. Section IV D assesses the validity of the multi-temperature (MT) model using present dissociation rate coefficients. Finally, Section V provides a summary of findings and concludes the present paper.

The ozone molecule plays a crucial role in atmospheric chemistry, demonstrating strong absorption properties of harmful UV radiation. For this reason, the O₃ molecule has been extensively studied experimentally^34–38 and theoretically.^28,39–41 Some interesting phenomena associated with ozone include the high enrichment of the atmosphere by heavy ozone isotopomers,⁴² negative temperature dependence of isotope exchange reaction,³⁸ the dependence of ozone rate formation on pressure,⁴³ and the nonstatistical behavior of reactive scattering in the exchange reaction.⁴⁴ These effects are relevant, first of all, to the chemistry of the atmosphere and, thus, the detailed discussion of this matter is outside the scope of the present paper. Refer to the review of current problems by Schinke et al.⁴⁵ 

The computational study of molecular systems with multiple open-shell atoms presents a great challenge due to the multiple PESs that correlate with the reactants. The interaction of triplet molecular oxygen with $O 3 P$ results in 27 different electronic states if spin-orbit coupling is neglected.⁴⁶ If spin-orbit coupling is considered, this manifold leads to 15-, 9-, and 3-fold degenerate states; however, only one state adiabatically correlates with the ground electronic state of ozone. Tashiro and Schinke⁴⁷ estimated the effect of spin-orbit coupling for these 27 electronic states and compared it to that of the wave packet calculations on a single PES of O₃ ground electronic state. This research was motivated by a disagreement in the rate coefficient at room temperature between global ab initio PES⁴⁸ and experimental measurements. Since full-scale quantum mechanical calculations including all relevant PESs are impractical, the influence of fine structure levels was estimated by calculating the probability of O₃ complex formation after oxygen atom crosses the transition-state region. These results indicated that the major contribution to the complex formation comes from the lowest state $O 3 P 2$⁠, and the spin-orbit coupling might not be the source of disagreement with experiment. It is worth to note that the electronic Coriolis and Renner-Teller coupling effects were omitted from this study.

Another problem associated with the O₃ molecule is the conical intersections of low-lying electronic states.⁴⁹ This matter has received special attention due to the impact of electronic structure on photodissociation.^50–52 In fact, the two lowest electronic states of the same symmetry, namely, 1¹ A′ and 2¹ A′, are connected by a conical intersection seam passing near the location of the transition state.^53,54 The presence of the conical intersection results in the complication of the Born-Oppenheimer approximation that is fundamental for spectroscopy and collision theory. As shown by Herzberg,⁵⁵ the electronic wave function changes sign when one traverses along the closed path around a conical intersection. In order to keep the total wave function continuous and single valued, Mead and Truhlar introduced a vector potential in the Schrödinger equation for nuclear motion.⁵⁶ Meanwhile, whenever the nuclear coordinate space encloses the conical intersection, the electronic wave function gains an additional phase. This phenomenon is known as the geometric or Berry phase effect.⁵⁷ The Hamiltonian of nuclear motion must be modified accordingly. For a review of computational methods regarding this problem, one may refer to the work of Yarkony.⁵⁸ 

One of the aspects closely associated with the study of the ozone complex is the vibrational spectrum and the dissociation limit of $O 2 3 Σ g −$⁠. The recent very accurate ab initio calculations of all valence shell electron correlations⁵⁹ indicate the existence of 42 vibrational states, while the experimental data on Schumann-Runge band predict only 36 vibrational states.⁶⁰ The most recent experimental measurement of O₂ dissociation energy in the ground state was reported by Ruscic⁶¹ with a scatter of 10 cm⁻¹ between previously reported data.⁵⁹ The correlation of this data with the energies levels used in the present work is discussed in Sec. II B.

A number of $O 3 X 1 A ′$ potential energy surfaces have been previously proposed in the literature.^28,62,63 A high quality PES, based on ab initio calculations, is also available.⁴⁸ Owing to the large cost of trajectory calculations, the present paper adopts the PES by Varandas and Pais²⁸ that presents a good compromise between accuracy and required computational time. This PES is obtained via the double many-body expansion method by the multi-property fitting of ab initio energies, existing experimental data on total scattering cross sections and the measurements of kinetic thermal rates. The Varandas PES reproduces the ab initio value of energy difference between open chain stable and cyclic metastable conformers, has no dissociation barrier, and reproduces the R⁻ⁿ function of long-range forces in the dissociation channel. The energy transfer between the kinetic motion of a projectile and the internal state of a target particle is efficient in the wide range of collision energy due to the absence of a potential barrier in this PES. This leads to a relatively weak dependence of rovibrational transition rates on temperature and a large probability of multiquantum transitions due to energy scrambling in the exchange reaction.^15,16

The potential energy surface by Varandas generates 48 vibrational levels and a maximum of 236 rotational levels for molecular oxygen in the ground electronic state. The vibrational energies and turning points of each rovibrational state are calculated by the Wentzel-Kramers-Brillouin (WKB) method⁶⁴ and are given elsewhere.⁶⁵ The total number of rovibrational levels in the ground electronic state of oxygen is 6245; however, even numbered rotational states are forbidden due to the selection rules of nuclear spin statistics. Since, during a bound-bound rovibrational transition, the symmetry of the initial and final states cannot change, transitions of the O₂ molecule in the ground electronic state are allowed only between odd-numbered rotational levels.⁶⁶ The dissociation energy with respect to the ground vibrational state, given by the Varandas PES, is 41 259.167 cm⁻¹, which is within the range between 41 256.6 and 41 269.6 cm⁻¹ of previously reported data⁵⁹ and only 9 cm⁻¹ less than the most recent estimation by Ruscic.⁶¹ The comparison of the present vibrational energy levels²⁹ with those obtained on the O₂ ab initio potential curve⁵⁹ reveals a difference within approximately 10 cm⁻¹ and 200 cm⁻¹ for low and high lying vibrational states, respectively. Taking into account that the characteristic temperatures, simulated in the present work, are between 5000 and 10 000 K, the Varandas PES and the WKB method seem to satisfy the accuracy of the present simulation. One should note that the Varandas PES slightly underestimates the vibrational energy levels, leading to an excess of 6 vibrational states, compared to the vibrational ladder by Bytautas.⁵⁹ 

The trajectory of a target molecule and a projectile atom is described by a system of six first order partial differential equations, formulated for the relative and internal systems of coordinates.⁶⁷ Solution of the governing equations is obtained by the 11th-order accurate Adams-Bashforth-Moulton method⁶⁸ with a variable time step between 10⁻¹³ and 10⁻⁷ s. The initial rotational and vibrational quantum numbers are assigned directly for each batch of trajectories, rather than sampling according to the Boltzmann distribution assuming some internal temperature. As discussed previously, this step makes calculations significantly more expensive, compared to the simulation with the Boltzmann sampling of rotational states; however, the present collisional model, based on the direct sampling, allows the investigation of thermal relaxation avoiding the assumption of trans-rotational equilibrium.

Due to the centrifugal forces arising from the rotation of a molecule, the vibrational and rotational levels of energy are coupled. In order to perform an appropriate sampling of trajectories in the QCT method, the energy of the internal degrees of freedom should be decoupled. In the present work, the “vibration-prioritized” framework^10,69 is adopted. The energy of a vibrational state is calculated as the energy of the rovibrational state with j = 0, while the energy of a rotational state is obtained as the difference between the rovibrational and vibrational energies.

In order to describe chemical kinetics at hypersonic temperatures, the range of collision energy is described by 43 unevenly spaced points and spans between from 5 × 10⁻³ to 10 eV. Trajectory calculations are carried out for each rovibrational state at a varying impact parameter of collision at each energy point. The increment of impact parameter is set to 0.1 Å, and stratified sampling is terminated when only elastic collisions are observed. Each batch in the given interval of impact parameter contains 2000 trajectories per single rovibrational level per energy point. The intermediate data are stored in order to improve the statistical accuracy of calculations by computing additional trajectories in the future. Each trajectory is integrated with the relative error in the total energy not exceeding 10⁻⁴%. The database of cross sections contains nearly 9.6 × 10⁶ bound-bound transitions, more than 6000 bound-free transitions and supports 690 quasi-bound (QB) states. Nearly 12 × 10⁹ trajectories have been processed on the supercomputer of the Advanced Research Computing Center in the University of Michigan.

The O₃ system is known for the close spacing of electronic levels. The energy difference between the ground $1 A 1$ and first excited $3 B 2$ states is nearly 3 eV.³⁹ Strictly speaking, the electronic excitation is probable at high collisions energies. One of the most rigorous techniques that deals with the nonadiabaticity of collisions is developed by Sholl and Tully⁷⁰ and involves the calculations of PES hopping probability and trajectory branching. In the present case, the simulation of electronic excitation would require to take into account hopping between multiple potential energy surfaces and an exhaustive trajectory branching, since some collisions occur with a large number of interatomic interactions.⁷¹ Such an approach would substantially increase the cost of the QCT method without significantly influencing the results in the range of temperature between 1000 and 5000 K, where nonadiabatic collisions play the dominant role. Moreover, the global PES of electronically excited ozone is not currently available. Instead, in the present work, all trajectories are forced to stay on the ground PES of ozone. At the post-processing step, the rates of bound-bound transitions are scaled to account for the spin and orbital degeneracy of the reactants. The dissociation rates are modified to take into account electronic nonadiabaticity, in a manner proposed by Nikitin.⁷² The correction procedure is discussed below.

Details of calculations of transition probabilities and cross sections were given previously.¹⁶ In the present work, the rates of bound-bound and bound-free rovibrational transitions in the range of temperatures between 500 and 20 000 K are generated. The former rates have to be modified for the spin degeneracy of reactants and PES. The O(³P) and $O 2 ( 3 X g − )$ reactant asymptote is 27-fold degenerate if spin-orbit coupling is neglected. Tashiro and Schinke demonstrated that the multistate calculations including the spin-orbit coupling result in similar O₃ formation probabilities compared to that of the single-surface approximation. In fact, the adiabatic calculations tend to underestimate the true rate coefficients.⁷³ The underestimation of O₃ formation probability using only $O 3 X 1 A ′$ was estimated as 10%-30%,⁴⁷ which is a reasonable level of uncertainty for the present model.

Since the complete coupling of all relevant states is computationally impractical, the present work neglects the spin-orbit coupling and adopts the following form of the degeneracy factor:³⁹ 

where the factor of 3 corresponds to the $O 2 ( X 3 Σ g − )$ triplet state. The expression in parentheses describes the degeneracy of $O 3 P J$⁠, with the angular momentum J = 0, 1, and 2. Now, the degeneracy factor of the bound-bound transition is found as a ratio of the PES degeneracy to the degeneracy of reactants. Taking into account that the $O 3 X ̃ 1 A ′$ PES is nondegenerate, the total degeneracy factor is given by Eq. (1). At high temperatures, the degeneracy factor asymptotically approaches a value of 1/27. This simple form of degeneracy factor should be reconsidered in more accurate calculations of reaction rates due to possible intersections of PESs and arising complications in the Born-Oppenheimer approximation.⁷⁴ In the present work, as the first step towards generation of a complete set of reaction rates for each O₂ rovibrational state, all trajectories are forced to stay on the ground PES of ozone.

The dissociation of molecular oxygen in a shock wave with a temperature above a few thousands of K is unlikely to proceed in the adiabatic manner. One way to account for alternative channels of dissociation from excited electronic states is to assume equilibrium between lower vibrational levels of an excited electronic state and high vibrational levels of the ground state.⁷² The contribution of dissociation from the excited electronic level should be added to the dissociation rate from the ground electronic level. It is possible to calculate the degeneracy factor of the dissociation reaction for each rovibrational level of the ground electronic state individually.⁷⁵ In this case, the degeneracy is calculated by including only electronic levels with a minimum of energy below the energy of the current rovibrational level. The degeneracy of the bound-free transition for the rovibrational level $v , j$ is then computed as follows:

where g^X = 3 is the degeneracy of the ground electronic state, and g^E is the degeneracy of the electronically excited state. The assumption of equilibrium between the vibrational levels of the ground and excited states is valid if the following condition holds:

Assuming that Eq. (3) holds for every rovibrational level of the ground electronic state, one can obtain the degeneracy factor of the bound-free transition, originally proposed by Nikitin,⁷² 

The system of master equations includes the processes of energy transfer in the inelastic (nonreactive) and exchange collisions, impact dissociation of bound and quasi-bound states, as well as the tunneling dissociation of the latter. The general appearance of the kinetic equation that describes these mechanisms has the following form:

where n_i is the population of the i rovibrational state, K_i′→i is the rate of transition from i′ to i rovibrational states, n_O is the number density of atom oxygen, D_i and R_i are the state-specific dissociation and recombination rate coefficients of the i state, $T i f$ and $T i b$ are the forward and backward tunneling rate coefficients, and N denotes the total number of rovibrational states. The summation in Eq. (5) is performed over the entire rovibrational ladder. The principle of detailed balance is invoked to generate rates of endothermic transitions, in order to reduce the statistical error of the QCT method. The rate of change of number density of atomic oxygen is given by the following equation:

In the present paper, rotational and vibrational relaxation is studied using two different thermodynamic models. The first model assumes rotational equilibrium in a chemical reactor. It is assumed that the rotational and translational temperatures are equal through out the calculations. This model is referred to as the VT thermodynamic model. However, rotational nonequilibrium may take place at some conditions, observed in hypersonic flows. The second thermodynamic model, adopted in the present work, treats the rotational and vibrational degrees of freedom in a similar manner. This approach is referred to as the RVT thermodynamic model.

Initial conditions must be specified for Eqs. (5) and (6). The number density of rovibrational level i at t = 0 is evaluated as follows:

where $Q i = 2 j + 1 exp − e i , vib / k T vib exp − e i , rot / k T rot$ is the rovibrational factor of state i, e_i,vib and e_i,rot are the vibrational and rotational energies of state i, and T_vib and T_rot are the energy-equivalent vibrational and rotational temperatures of O₂ gas, defined below. The VT and RVT models differ by initial conditions, specified for Eq. (5). For the former, it is assumed that T_rot = T and T_vib = T₀, where T is the translational temperature of a heat bath, set to a constant value for all calculations, T₀ is the initial temperature of O₂ gas, set to 100 K. For the RVT model, it is assumed that T_rot = T_vib = T₀.

The solution of Eq. (5) describes the time-dependent population of each rovibrational level during the relaxation to thermal equilibrium. Because the large number of states complicates further analysis, the energy-equivalent vibrational and rotational temperatures are used to describe the population of the entire vibrational ladder. The average vibrational and rotational energies are evaluated in the following manner:

where x stands for either vibrational or rotational mode. Due to the “vibrationally prioritized” paradigm of decoupling of rovibrational energy, the vibrational temperature is defined as a function of e_vib only, while the rotational temperature is defined as a function of both e_vib and e_rot. The following implicit equations define T_vib and T_rot,

where $e vib =∑exp − e i , vib / k T vib e i , x /∑ − e i , vib / k T vib$ and e_rot = ∑Q_ie_i,rot/∑Q_i. The solution of Eq. (9) is obtained by the bisection method. It is important to note that the internal temperatures, defined by Eqs. (8) and (9), are evaluated only at the post-processing step, after the solution of master equations has been obtained.

It is possible to define the rotational and vibrational relaxation times based on the temporal evolution of internal energy under heat bath conditions. There are several ways to define a characteristic relaxation time.⁷⁶ One approach is to calculate the relaxation time as if the equilibration process follows the Landau-Teller (LT) equation, the so-called e-folding method:⁷⁶ 

where x corresponds to either rotational or vibrational mode. The internal energy, given by Eq. (11), evolves exponentially in time. The relaxation time, τ_x, is then defined by the corresponding energy e_x,efold,

where e_x,efold is the average energy at the time τ_x and $e x T 0$ and $e x T$ are the average energies evaluated at the initial and equilibrium temperatures, respectively.

The attractive component in the O₃ potential energy surface has a significant influence on the rates of mono- and multiquantum state-to-state transitions.^25,77,78 Further understanding of the influence of the $O 3 1 A 1$ PES on the dynamics of O₂ thermal relaxation comes from the comparison of the temperature dependence of transition rates in the O₂–O system with that of a collisional pair with a large repulsive potential. Recently, the QCT method was applied to generate the complete set of O₂–Ar bound-bound transition rates.³² In the present paper, the comparison is performed for the rates of vibrational deactivation, averaged at rotational temperature T_rot = T in the range between 1000 and 10 000 K. The rates of v = 1 → v′ = 0, v = 10 → v′ = 9, and v = 10 → v′ = 5 transitions are shown in Fig. 1. For reference, the equilibrium dissociation rate coefficients, obtained from averaging of state-specific rate coefficients assuming T_vib = T_rot = T, are shown by the long dashed symboled lines.

At low temperatures, the O₂–O vibrational transition rates are substantially higher than those for the O₂–Ar system. The deactivation in the former molecular system demonstrates a very weak temperature dependence, while in O₂–Ar collisions the probability of vibrational quantum jump quickly increases at highly energetic conditions. Moreover, the rates of O₂–O multiquantum transitions are comparable to those of a single vibrational jump, while in O₂–Ar collisions the multiquantum jumps are much less probable. At high temperatures, the O₂–O vibrational deactivation occurs slower, compared to that in O₂–Ar collisions. This can be explained by a diminishing role of the O₂–O exchange mechanism¹⁶ and a less pronounced influence of the O₃ potential well⁷¹ at high energies. At the same time, the dissociation rate coefficients in O₂–O and O₂–Ar systems increase in a similar manner. For the O₃ system, the dissociation at high temperatures occurs faster than the relaxation. This phenomenon is related to the efficient scrambling of internal pre-collisional states via the exchange reaction and subsequent internal excitation that precedes the dissociation.⁷⁹ This observation suggests that the coupling of relaxation and dissociation processes, known as the quasi-steady state (QSS), proceeds in the O₂–O system in a different manner, compared to that in a collisional pair with a strong internuclear repulsion.

An important insight on the process of relaxation can be gained by analyzing the state-specific populations of the rovibrational ladder during thermalization to the equilibrium state. To simplify the analysis, the dissociation is artificially excluded from this simulation. The solution of master equations for the constant heat bath conditions of 10 000 K is shown in Fig. 2. All populations n_i are normalized to the degeneracy of rotational angular momentum, 2j_i + 1, and to the total number density of oxygen molecules, which is set to a constant of 10¹⁵ cm⁻³. At t = 0, most molecules are in the ground vibrational state and occupy the first few low-lying rotational levels. As the relaxation continues, excited states become populated. The bump in the population of molecules with the energy close to the dissociation threshold corresponds to the large probability of dissociation and a relatively low probability to remain in a bound state. Overall, the scattering of populations is narrow due to the effectiveness of relaxation via multiquantum jumps. The majority of excited states are at a similar internal temperature, which is, as follows from the slope of the normalized population distribution function, significantly higher than the temperature of the ground vibrational state. The scattering of populations of quasi-bound states is only slightly broader than that of the bound states. The mechanism of population of QB states is discussed below.

The conventional behavior of the population distribution function at the late phase of relaxation is a dissection into separate strands for each low-lying vibrational state.^10,80 The close-up view of these strands is shown in Fig. 3 at t = 10⁻⁷ s. Different colors represent low-lying vibrational states v = 0…8. The normalized population distribution function of the ground vibrational state, v = 0, significantly differs from the rest of the vibrational ladder. Each subsequent vibrational state is populated less densely than the previous one; however, they all demonstrate a very narrow scattering of normalized populations. The tails of the vibrational ladders nearly overlap, indicating that the relaxation proceeds at nearly identical rotational temperatures in these vibrational states.

It is possible to calculate the state-specific rotational temperature, $T rot v$ for each vibrational state using the energy-equivalent definition of temperature, similar to the one defined in Eq. (9). In the present work, the following equation for $T rot v$ is solved by the bisection method:

where the summation is performed over the range of rotational quantum numbers j, compatible with the current vibrational state v. The evolution of state-specific rotational temperature with time is shown in Fig. 4 for T = 10 000 K. In the early phase of relaxation, T_rot of the ground vibrational state is frozen at the initial conditions. The majority of vibrational states are, indeed, populated at nearly the same rotational temperature, which explains the narrow distribution of populations in Fig. 2. Finally, the top-lying vibrational states are populated at much lower rotational temperature, due to the relatively slow excitation of these vibrational states. The flat distribution of $T rot v$ can be attributed to the efficient energy randomization via multiquantum state-to-state transitions in O₂–O collisions.¹⁵ As relaxation proceeds, $T rot v = 0$ rapidly increases, while the rotational temperatures of the other vibrational states are nearly constant and slowly increasing toward the thermal equilibrium. In the late phase of relaxation, the highly excited states eventually become thermalized.

The exchange channel plays an important role in the energy randomization in molecule-atom collisions. The role of exchange reactions is especially significant in the O₃ molecular system, since the projectile can be easily trapped by the potential well, effectively converting the kinetic energy of motion into the internal energy of products in the excited rovibrational state.⁸¹ During the insertion of a projectile atom, the system has a short memory about the initial rovibrational state, and the final state can take a wide range of possible quantum numbers, defined by the kinetic energy of the collision. The database of cross sections, generated in the present work, separately contains the data about the contribution of inelastic (nonreactive) and exchange channels of collisions, which can be converted into a set of transition rates accounting for either inelastic channel or both channels of collisions. Following this approach, it is possible to estimate the contribution of the exchange channel on the relaxation at different temperatures of a heat bath.

The evolution of population of the rovibrational ladder is shown in Figs. 5 and 6 for the heat bath conditions of 3000 and T = 10 000 K, respectively. Dissociation is not considered. The populations are shown at time between 10⁻¹¹ and 10⁻⁵ s for 3000 K and 10⁻¹¹ and 10⁻⁷ s for 10 000 K due to the different pressures in the chemical reactor. In the early stage of relaxation, the majority of molecules are in the low-lying rovibrational states, while the tail of the population is sparsely populated. The spectrum of populations is significantly broadened when only nonreactive collisions are considered. More importantly, the population of quasi-bound states is governed mostly by the exchange channel of interaction, which is explained by a large probability of energy randomization in exchange collisions. Failure to account for the exchange reaction can result in strong underestimation of the highly excited states.

Further study of O₂–O relaxation dynamics analyzes the average internal energies and temperatures of the rovibrational ensemble. The variation of average vibrational and rotational energies is shown in the range of translational temperatures between 500 and 20 000 K in Fig. 7. No dissociation is considered in these calculations. The initial number density of oxygen atoms is set to 9.99 × 10¹⁷ cm⁻³ and the number density of molecular oxygen is equal to 10¹⁵ cm⁻³. The initial population of the rovibrational ladder is computed according to the Boltzmann distribution at a temperature of 100 K.

At the initial phase of relaxation, most molecules occupy the ground vibrational state with an energy of 0.097 62 eV and low-lying rotational states, which explains lower e_rot compared to e_vib. Thermalization of the rotational and vibrational degrees of freedom occurs within the same time scale at these temperatures, which corresponds to the fact that rotational nonequilibrium in a heat bath of oxygen atoms takes place at all temperatures observed in hypersonic flows. This result is profoundly different from the behavior of the rotational mode in other molecular systems at low temperatures. The master equation simulations of other molecular systems revealed a much faster rotational relaxation in N₂–N and O₂–Ar collisional pairs.^10,32

Due to the enormous number of rovibrational states, it is convenient to describe the relaxation process of different modes by their internal temperature. An average energy of the ensemble can be used to define such energy-equivalent O₂ vibrational and rotational temperatures. In the present work, these temperatures are calculated only at the post-processing step. The average vibrational and rotational temperatures, defined by Eq. (9), are shown in Fig. 8 for T between 500 and 20 000 K. In this range of translational temperatures, the vibrational and rotational temperatures are closely coupled due to the very efficient energy randomization in the O₃ complex.⁸² 

The average internal energies, shown in Fig. 7, can be used to define the mean relaxation times of the vibrational and rotational modes. In the present work, τ_vib and τ_rot are calculated by the e-folding method.⁸³ Vibrational and rotational relaxation times, obtained from the solution of master equations using the RVT thermodynamic model, are shown in Fig. 9 with circular and diamond symbols, respectively. The vibrational relaxation time, assuming the VT model, is shown by the dashed line. Available experimental data by Breen et al.²⁵ are shown by empty symbols and the vibrational relaxation time, derived by Park⁸³ from the Millikan-White relation,⁸⁴ is shown by the dashed-dotted line. The ratio of τ_vib and τ_rot, obtained from the RVT model, is shown by deltas.

The RVT rotational and vibrational relaxation times are comparable to each other in the considered temperature range. This result suggests that the modeling of O₂–O collisions should include individual simulation of both rotational and vibrational modes. The ratio of τ_vib and τ_rot increases at high temperatures, indicating that at these conditions the vibrational and rotational relaxations are strongly intertwined and proceed at almost the same rate. Since the vibrational and rotational temperatures closely follow each other, as can be seen from Fig. 8, the population of the rovibrational ladder can be approximately described by a unique internal temperature. Previously, the existence of such an internal temperature for N₂–N relaxation at very high temperatures was pointed out by Panesi et al.¹⁰ The vibrational relaxation time, derived from the VT thermodynamic model, is smaller than that from the RVT model. Since the vibrational and rotational relaxations proceed within the same time scale, the rotational nonequilibrium delays the overall thermalization, extending the O₂–O vibrational relaxation. The relaxation parameter, obtained via the VT model, describes the existing experimental data by Breen et al.²⁵ best, which implies the presence of rotational equilibrium in the oxygen gas during experiments. Nevertheless, the present data on $p τ vib R V T$ and $p τ rot R V T$ should be useful in the construction of reduced order models.^12,30

The present relaxation times can be fit by a polynomial function of translational temperature in the range from 500 to 20 000 K. Due to the slow variation of relaxation times in this temperature range, the polynomial function is used,

where P is the pressure of the gas of oxygen atoms, Pτ_x is in atm s, T is in K.

Vibrational and rotational relaxation times in the presence of only inelastic collisions are shown in Fig. 10 with the dashed lines and circular and diamond symbols, respectively. This simulation excludes the relaxation via the exchange channel. The relaxation times, obtained when both nonreactive and exchange channels are present, are shown by solid lines. Thermalization in the absence of exchange channels proceeds more slowly, which confirms the importance of including this reaction in the temperature range of interest. To compare the contribution of the exchange mechanism at different temperatures, the ratio of relaxation times, $p τ nonreact / p τ total$⁠, is shown by the dashed-dotted lines in Fig. 10. Red and black colors correspond to vibrational and rotational modes, respectively. For the vibrational degrees of freedom, this ratio slowly varies in the range between 1.7 and 1.9. However, the influence of the exchange channel on the rotational mode is more pronounced: at higher temperatures, the rotational relaxation in the presence of the exchange channel proceeds substantially faster than that at low temperatures.

Characteristic relaxation times, obtained via the RVT thermodynamic model, can be used to verify the possibility of describing the relaxation process in O₂–O collisions by the simple LT model. For the latter, the evolution of the average vibrational and rotational energies can be described by Eq. (11). Using the e-folding RVT relaxation times from Table I as the relaxation parameter τ_x, the evolution of the average energy of ensemble, e_x, is obtained. The average rotational and vibrational energies, computed from the system of master equations under the heat bath conditions of 1000, 10 000, and 16 000 K, are compared to those from Eq. (11) in Fig. 11. The solution of master equations is shown with solid, dashed, and dashed-dotted lines and symbols correspond to Eq. (11). At first glance, the agreement is satisfactory for both rotational and vibrational modes of relaxation. However, the relative difference between average energies, shown in Fig. 12, indicates that these discrepancies have a persistent pattern. The deviation of average rotational energy predicted by the Landau-Teller model from the solution of master equations is typically larger than that for the vibrational mode and increases with temperature. The LT model has better accuracy at low temperatures, which is an expected conclusion.⁸³ The maximum discrepancy is as high as 14% in the considered temperature range, which is, generally, an acceptable inaccuracy of reduced order models.

Thermal relaxation of molecular oxygen in the presence of dissociation and recombination processes is studied in this section. The tunneling of quasi-bound states is also included in the simulations of heat bath at constant translational temperatures between 1000 and 20 000 K.

The energy-equivalent vibrational and rotational temperatures with dissociation and recombination processes enabled are shown in Fig. 13 for translational temperatures between 3000 and 16 000 K. Several important observations follow from these simulations. First, because at low temperatures the vibrational and rotational relaxation occurs much faster than dissociation, the latter takes place from the rovibrational manifold, populated nearly at the equilibrium temperature. As follows from Fig. 13, this regime takes place at translational temperatures below 5000 K. At higher temperatures, the vibrational and rotational degrees of freedom do not completely equilibrate before dissociation occurs. During this phase, known as the QSS regime, T_vib and T_rot are substantially lower than the temperature of the heat bath. It is worth to note that T_rot is higher than T_vib, which possibly means that the bath of molecules preferably dissociates from higher vibrational and lower rotational states. The energy rate coefficients, that define the average loss of internal energy due to dissociation, are calculated below.

At a translational temperature of 10 000 K, the quasi-steady values of T_vib and T_rot experience a maximum. The explanation of this fact lies in the behavior of the relaxation and dissociation rate coefficients in O₂–O collisions at high temperatures. While the latter has a conventional, Arrhenius-type, temperature dependence, the former only weakly increases with temperature. Such temperature dependence is different from what is typically observed in other molecular systems.³³ Since, at some translational temperature, the dissociation becomes substantially faster than the relaxation, the depletion of the rovibrational manifold takes place at lower internal temperature, compared to that in a heat bath with lower translational temperature. Moreover, the QSS phase becomes shorter and the quasi-steady plateau is smeared. At temperatures higher than 10 000 K, T_rot and T_vib change significantly during the QSS phase. This means that the relaxation and dissociation proceed within the same time scale and should be modeled only in a concurrent manner. The non-monotonicity of O₂ vibrational and rotational temperatures during the QSS phase is a remarkable feature of the O₂–O system, different from the steady increase of internal temperatures in N₂–N collisions.¹⁰ These results should be carefully reviewed when collisional models of electronic excitation become available.

The QSS global dissociation rate coefficient, D^QSS, is calculated from the solution of master equation in a manner adopted in Ref. 32. The present D^QSS, obtained from the RVT and VT thermodynamic models, is shown in Fig. 14. Both RVT and VT dissociation rate coefficients are substantially lower than the equilibrium rate coefficient, D^eq. The ratio of D^QSS,RVT and D^eq increases with temperature from a factor of 3–39. The lower range of D^QSS,RVT/D^eq is in agreement with the statement made by Park,⁸⁵ while at higher temperatures the QSS rate coefficient is much lower than expected. This is explained by the declining internal temperature during the QSS phase, as can be seen in Fig. 13. Meanwhile, the dissociation rate coefficient, obtained from the VT model, is higher than D^QSS,RVT, which is due to the rotational nonequilibrium over the entire temperature range of interest.

The present equilibrium and quasi-steady dissociation rates are compared with the existing experimental data in Fig. 15. The red curve corresponds to D^eq estimated from the QCT data assuming T_rot = T_vib = T. Black and blue curves describe the QSS dissociation rate coefficients, obtained via the RVT and VT models, respectively. In these simulations, the variable g^BF factor is utilized. The experimental data are taken from the work of Shatalov⁸⁶ and Byron.⁸⁷ In the former work, the equilibrium dissociation rate was originally reported; hence circular symbols should be compared to D^eq, given by the red curve. Byron reported the equilibrium rate, given by the dashed line, as well. Since the experimental measurements describe the dissociation rate during the quasi-steady phase, it is of interest to compare the actual experimental measurements with the rates computed for the QSS regime. In the present work, Shatalov’s data points for equilibrium dissociation rate are converted to D^QSS using the methodology described therein. Triangular symbols in Fig. 15, that describe measurements of D^QSS, should be compared with the present D^QSS,RVT and D^QSS,VT dissociation rate coefficients.

The agreement of equilibrium dissociation rate coefficient with experimental data is very good. The measurements of D^QSS, which lie below D^eq due to the incomplete thermal relaxation prior to the quasi-steady regime, demonstrate excellent agreement with the present D^VT,QSS dissociation rate. The rate coefficient, computed under the assumption of rotational nonequilibrium, D^RVT,QSS, falls in the lower range of the experimental data, possibly indicating that the oxygen gas was in rotational equilibrium in the original experimental setup. Present equilibrium and quasi-steady dissociation rate coefficients are given in Table II.

Normalized populations of the rovibrational ensemble at different stages of thermal relaxation are shown in Figs. 16-18 for the heat bath conditions of 5000, 10 000, and 20 000 K. In these calculations, the quasi-bound states are assumed to have a finite lifetime, defined by the probability of tunneling. The initial total number density of particles in the chemical reactor is 10¹⁸ cm⁻³; the molar fraction of atomic oxygen is set to 5%. The distribution of internal states is shown at times between 10⁻¹² and 6 × 10⁻⁵ s. The population distribution function in the midst of the QSS phase and at equilibrium is shown by orange and maroon symbols, respectively.

The scattering of the population distribution function is very narrow throughout the relaxation process in this temperature range, similarly to what is observed in Fig. 2. A fraction of molecules in quasi-bound states is depleted before the QSS state begins, due to their short lifetimes. At low temperatures, the depletion mechanism does not play a significant role until the majority of states are thermalized. The equilibration for low-lying rovibrational states at T = 5000 K is nearly complete before the gas becomes significantly dissociated. During the QSS phase, a significant conversion of chemical energy takes place. This energy is removed from the excited states, as follows from Fig. 16. The tail of the rovibrational ladder, that corresponds to highly excited bound and quasi-bound states, is strongly underpopulated, indicating the presence of a preferential mechanism of dissociation at these conditions.⁶ 

Due to the strong dependence of the dissociation rate coefficient on temperature, the incubation period prior to the QSS regime is shorter at higher temperatures. In the case of 10 000 K, thermalization of low-lying states at the beginning of the QSS phase is not complete, and the population distribution function of bound states is broadened, compared to that at T = 5000 K. As follows from Fig. 17, the ground rovibrational states are overpopulated due to incomplete relaxation to the heat bath conditions. Rovibrational states with energy more than 3 eV are significantly underpopulated. As follows from the slope of the population distribution function, the vibrational and rotational temperatures during the QSS phase are significantly lower than their equilibrium values. At translational temperature of 20 000 K, the nonequilibrium effects are even more pronounced than at T = 10 000 K. The low-lying states are strongly overpopulated which corresponds to a lower internal temperature during the QSS phase, compared to that at T = 10 000 K. The rovibrational states with energies more than 2 eV are underpopulated due to the incomplete relaxation. Overall, the underpopulated branch of the ensemble shifts toward the lower internal energies at higher temperatures.

More insight can be obtained by analyzing $T rot v$ in the late phase of relaxation and during the QSS phase. The evolution with time of rotational temperatures at T = 5000, 10 000, and 20 000 K is shown in Figs. 19–21. Thermalization of the rotational degrees of freedom proceeds differently at low and high translational temperatures of the heat bath. At T = 5000 K, the majority of states, excluding those with v = 0, have nearly the same rotational temperatures in the late phase of the incubation period, shown by the blue and green lines. As the QSS phase begins, the equilibration occurs initially for low-lying states. The dip in rotational temperature in the middle of the vibrational manifold is caused by the maximum dissociation from these states, as will be shown later. At translational temperature of 10 000 K, the equilibration of low- and high-lying states occurs nearly at the same time. The minimum in $T rot v$ is now shifted to the lower energies, compared to what is observed at T = 5000 K. Finally, rotational thermalization at T = 20 000 K occurs first for highly excited states, while the low-energetic states are delayed.

In order to describe the effect of rovibrational relaxation and dissociation that proceed concurrently in the O₂–O system at high temperatures, it is convenient to define the average loss of energy in collisions that lead to dissociation. In the present paper, such coupling coefficients between the relaxation and dissociation processes are calculated for the RVT thermodynamics models. The average loss of energy in a single dissociation event is calculated as follows:⁷ 

where x stands for the vibrational, rotational, or rovibrational mode. By substituting the corresponding internal energy, the dissociation-vibration, dissociation-rotation, and dissociation-rovibration coupling coefficients are obtained. The summation in Eq. (14) is performed over the entire rovibrational ladder. In the present work, energy rates coefficients are normalized to the O₂ classical energy of dissociation, D_e = 5.212 93 eV.

The state-specific dissociation rate coefficients are strongly dependent on the factor of nonadiabaticity, given by either Eq. (2) or (4). Since the energy rate coefficients directly depend on the dissociation rate, the average loss of internal energy during depletion is sensitive to the choice of g^BF. The influence of the nonadiabaticity factor on the energy rate coefficients is shown in Fig. 22. The solid lines correspond to the constant g^BF and the dashed lines describe the results obtained using the variable g^BF. The data, reported in Fig. 22, correspond to the QSS phase.

Several important observations follow from Fig. 22. First, one should notice a relatively large drop of energy rate coefficients with temperature, compared, for example, to similar calculations for the O₂–Ar³² and N₂–N¹⁰ systems. At T = 10 000 K, which corresponds to mild hypersonic conditions, the loss of rovibrational energy is only 0.6 and the loss of vibrational energy is less than 0.4, compared to C^DV = 0.7 in N₂–N and 0.6 in O₂–Ar collisions. The O₂–O energy rate coefficients quickly diminish at high temperatures. One of the possible explanations lies in the fact of relatively inefficient rovibrational relaxation in O₂–O collisions in highly energetic collisions, while the dissociation rate coefficient quickly increases with temperature.⁶⁵ During the QSS phase, dissociation takes place from the rovibrational ladder populated at a temperature that is substantially lower than the equilibrium one. In fact, as can be seen from Fig. 13, this effect is stronger at higher temperatures. As a result, the average energy, removed due to dissociation, is essentially lower than the classical dissociation threshold.

Second, the energy rate coefficients are very sensitive to the nonadiabaticity factor. At lower temperatures, the influence of g^BF is less pronounced due to the smaller probability of electronic excitation for low-lying states. In fact, variable g_BF affects only low-lying states, reducing their contribution to the global dissociation rate coefficient. Another reason, which is discussed below, is the preferential dissociation of highly excited states at low temperatures. However, at high temperatures, the difference between energy rate coefficients, computed with various g^BF, is more than 40%.

Finally, the rotational mode demonstrates a non-monotonic behavior of energy rate coefficient, which is different from similar calculations for other molecular systems.^10,32 At low temperatures, the vibrational coupling coefficient, C^DV, is substantially higher than C^DR, which partially explains the higher rotational temperature during the QSS phase, as follows from Fig. 13. Since the internal energy is removed mostly from the vibrational mode, it is expected for T_vib to be lower than T_rot. At high temperatures, the rovibrational relaxation lags behind the dissociation process, and the depletion of the ladder occurs at low vibrational and rotational temperatures during the quasi-steady regime. At these conditions, the vibrational and rotational energy exchanges are strongly intertwined, and C^DV and C^DR both have a negative temperature dependence. Note that the vibrational and rotational energy rate coefficients do not converge to a common asymptote, as stated in Ref. 10. In fact, as shown in Ref. 32, the C^DV and C^DR curves intersect at high temperatures. In the present work, C^DR is only slightly larger than C^DV at temperatures of 14 000 K and higher. In the light of the higher rotational temperature during the QSS phase at these conditions, one may conclude that besides the energy removal from internal modes due to chemistry, other processes, such as inelastic and exchange collisions, play an important role here. Again, these results should be carefully reviewed in future, when accurate models of collisional electronic excitation become available.

A comparison of coupling coefficients, obtained for the O₂–O, N₂–N,¹⁰ and O₂–Ar³² systems, is shown in Fig. 23. Similar to the O₂–O calculations, results for other molecular system are obtained using the RVT thermodynamic model. Since in the original paper by Kim and Boyd³² the values of C^DRV were not presented, the missing information is computed in the present work in a similar manner as for the O₂–O calculations. The vibrational energy rate coefficient demonstrates similar behavior for all three collisional pairs. However, for O₂–O collisions, C^DV is substantially lower over the entire temperature range. This can be attributed to a faster dissociation in O₂–O collisions, as can be seen from Fig. 1. The second reason, which plays an important role at high temperatures, is a relatively slow relaxation in the O₂–O collisional pair. Combination of these two factors leads to a smaller loss of vibrational energy during the quasi-steady phase of relaxation, compared to the N₂–N and O₂–Ar systems. The incomplete thermalization of the vibrational manifold certainly reduces the average loss of internal energy, as follows from comparison of C^DRV in Fig. 23. It is interesting to note that C^DR for O₂–O is higher than that in other collisional pairs at low temperatures, which can be partially explained by the efficient excitation of highly excited bound and quasi-bound states due to the exchange mechanism. It was recently shown¹⁶ that the exchange reaction is important for multiquantum excitation at low temperatures.

The energy rate coefficients are important parameters in reduced order models. Results, reported in Fig. 22, can be included in multi-temperature thermodynamic models^3,12,30,83 in order to describe the rotational and vibrational relaxation behind the shock front. In the present work, C^DV, C^DR, and C^DRV, obtained with the variable g^BF, are presented as a function of translational temperature with coefficients given in Table III. The bi-exponential form is adopted to describe the temperature dependence,

where T is in K and all other parameters are dimensionless. This curve fit data reproduce the temperature dependence of energy rate coefficients in the interval between 500 and 20 000 K. The data, reported in Table III, should be multiplied by the oxygen dissociation threshold in order to obtain energy units.

The state-specific energy rate coefficients are calculated in a manner, similar to Eq. (14), taking into account only contributions from particular rovibrational states. The sum of all state-specific contributions, $C i D x$⁠, exactly equals to C^Dx. The state-specific $C i DV$ and $C i DR$ are shown in Figs. 24-26 for the heat bath conditions of 5000, 10 000, and 20 000 K, respectively. Red and black solid lines represent the cumulative C^DR and C^DV, respectively. At each energy, the cumulative curve gives the sum of $C i D x$ with rovibrational energies less than or equal to the current value of interest.

The bulk contribution to C^DR and C^DV comes from the highly excited bound states in the temperature range between 5000 and 20 000 K. The analysis of cumulative curves indicates the minor contribution of quasi-bound states: nearly 8% in the rotational mode and less than 5% in the vibrational mode. The maximum of $C i D x$ shifts from 4.3 eV at 5000 K to 1.67 eV at 20 000 K. In other words, at high temperatures, the low-lying states lose most of the internal energy during dissociation. The maximum energy rate coefficient has a direct influence on the rotational temperature of individual vibrational states, as can be seen from Figs. 19-21. Those rovibrational states with the largest loss of internal energy due to chemical reactions have the lowest rotational temperature.

Cumulative functions of energy rate coefficients, shown in Figs. 24-26, imply a small influence of quasi-bound states on the chemical and thermal relaxation in O₂–O collisions. This observation differs from the result of N₂–N analysis,¹⁰ where the contribution of quasi-bound states on C^DR and C^DV is 80 and 40%, respectively. In the present work, the distribution of internal states, shown in Figs. 16-18, demonstrates a strong underpopulation of the quasi-bound manifold during the QSS state of relaxation. This can serve as a partial explanation of the relative unimportance of QB states. Another factor, that plays a crucial role in thermal relaxation, is the dissociation probability density function (PDF), which can be estimated individually for each rovibrational state. In the present work, the following equation is used to calculate the dissociation PDF:

where p_i is the dissociation PDF of state i and n_O2 is the instantaneous total number density of molecular oxygen. The evolution of dissociation PDFs with time is shown in Figs. 27 and 28 for heat bath temperatures of 5000 and 10 000 K, respectively. The dissected strands describe the strong dependence of dissociation probability on the rotational quantum number. Each strand corresponds to a single vibrational state. Low-lying rotational states have PDFs severals orders of magnitude lower compared to that of states with high j. This is valid for the majority of bound states, excluding the ground vibrational state. Since, at the early stage of relaxation, nearly all molecules are in the ground vibrational state and occupy the first few low-lying rotational levels, the dissociation PDF of these states is substantially higher compared to that of other states. When depletion is starting to affect the internal distribution of the ensemble, the population of low-lying rovibrational states changes drastically, which is consistent with the high dissociation PDF of these states.

During the QSS stage, the maximum of the dissociation PDF shifts toward the highly energetic bound rovibrational states. Note that quasi-bound states have a relatively low PDF during the QSS phase, which is explained by their relatively low population, as follows from Fig. 17. In fact, the pre-dissociated states contribute only 5% to the cumulative PDF at translational temperatures of 5000 and 10 000 K. The cumulative PDFs shift toward the low rovibrational energies at high temperatures, which is consistent with the shift of state-specific energy rate coefficients, shown in Figs. 24-26.

This section investigates the accuracy of the MT model for the description of O₂ dissociation in the presence of vibrational nonequilibrium. For this purpose, Park’s model of vibration-dissociation coupling is used. The governing equations of the zero-dimensional MT model are as follows:

In Eq. (17), e_v and $e v ⋆$ are the O₂ vibrational energy evaluated at T_vib and T, respectively, ρ and ρ_O2 are the density of O₂–O mixture and partial density of O₂, D_e is the classical dissociation energy, and τ_vib is the relaxation time, taken from Fig. 9, dashed line. The only type of collision considered is between O₂ and O. Initially, the total number density is set to 10¹⁸ cm⁻³ with 5% of atomic oxygen molar fraction. The governing temperature T_a is calculated as $T T vib$⁠. The global recombination rate coefficient, R, is estimated from D via the principle of detailed balance. Vibrational energy coupling coefficient corresponds to the loss of internal energy in O₂–O collisions in the presence of rotational equilibrium.⁸⁸ 

The global dissociation rate, D, is calculated according to the recommendation by Park,⁸³ as well as using the present, thermal equilibrium rate, shown in Fig. 15 by the red line. For the latter, the data are curve fitted to the Arrhenius form in order to enable straightforward coupling with Park’s model. For the purpose of comparison, instead of Park’s model, the actual QSS dissociation rate coefficient, given in Fig. 15 by the black line, is used. When the quasi-steady rate coefficient is utilized, Park’s model should not be introduced in Eq. (17).

Comparison of T_vib and n_O via the MT and STS approaches for the constant translational temperature of 10 000 K is shown in Figs. 29-31. Park’s model coupled to Park’s dissociation rate⁸³ is shown in Fig. 29. Park’s model coupled to the present thermally equilibrium dissociation rate is shown in Fig. 30. The solution of Eq. (17), using the present QSS rate without Park’s model, is shown in Fig. 31. Both simulations, using Park’s model, predict much later dissociation of oxygen, compared to STS approach. The STS approach indicates the depletion at the early stage of vibrational relaxation. In this situation, the QSS assumption is not valid: the time scales of thermal relaxation and dissociation are similar, and these processes should be considered simultaneously. In fact, the instantaneous and QSS dissociation rate coefficients have similar values at T = 10 000 K.⁸⁸ This explains the satisfactory accuracy of the MT model that utilizes the QSS rate with no vibration-dissociation coupling, as can be seen from Fig. 31.

In the situation, when thermalization is nearly complete prior to the onset of dissociation, the relaxation and dissociation mechanisms can be virtually decoupled from each other. The evolution of T_vib and n_O at constant heat bath conditions of 5000 K is shown in Figs. 32-34. The STS approach demonstrates later dissociation compared to the MT model. Again, the amount of atomic oxygen is most accurately described in the case when the MT model is used with the constant QSS dissociation rate coefficient. One of the reasons, resulting in a difference, is the overestimation of dissociation rate when the thermal equilibrium rate coefficient is applied. For example, the present QSS and equilibrium rate coefficients, estimated at T_a = 5000 K, are different by a factor of 3.5. This difference is explained by the depopulation of excited vibrational state during the QSS phase that is responsible for O₂ depletion.

A comprehensive study of the $O 2 ( X 3 Σ g − ) – O 3 P$ molecular system is presented. Thermal relaxation and chemical transformation in the oxygen gas heated up to 20 000 K are simulated by means of master equations. The extensive QCT simulation of O₂–O collisions using the O₃ DMBE potential energy surface by Varandas and Pais precedes this master equation study. The investigation is performed assuming the existence of only relaxation processes in a chemical reactor as well as studying the relaxation and dissociation in a concurrent manner. The e-folding O₂–O vibrational and rotational relaxation times are obtained over a wide range of temperature taking into account excitation of the entire rovibrational ladder. The present results indicate a temperature dependence of vibrational relaxation time that is different from the widely used Millikan-White formula. Namely, the relaxation proceeds with extreme efficiency at low temperatures, while at high temperature the vibrational relaxation time slowly varies with temperature. Rotational relaxation is tightly intertwined with the vibrational energy transfer in the studied temperature range. These features are attributed to the efficient energy randomization in the O₂–O system, particularly due to the absence of a potential barrier. Further simulations of O₂–O collisions should be performed in a manner that accounts for rotational nonequilibrium, if no alternative channels of relaxation are present.

The weak temperature dependence of the rovibrational transition rates strongly affects the properties of the O₂–O quasi-steady phase. The QSS vibrational and rotational temperatures do not increase monotonically with the temperature of the heat bath. Instead, there is a distinct maximum of internal temperature, which is defined by the simultaneous relaxation and dissociation processes in the O₂–O mixture. At lower translational temperatures, the rovibrational relaxation is significantly faster than dissociation, while at higher temperatures the relaxation lags behind.

The inefficient rovibrational relaxation in highly energetic collisions leads to a rapid drop of energy rate coefficients at high temperatures. Due to the strong vibrational-rotational coupling at these conditions, the rotational energy rate coefficient demonstrates a non-monotonic behavior with temperature. The analysis of cumulative curves of energy rate coefficient indicates a minor influence of quasi-bound states during the process of equilibration to thermal and chemical equilibrium. Instead, at high temperatures, the maximum of the dissociation probability density function shifts toward the low-lying bound rovibrational states.

The authors gratefully acknowledge funding for this work through Air Force Office of Scientific Research Grant No. FA9550-12-1-0483.