Ab initio-informed maximum entropy modeling of rovibrational relaxation and state-specific dissociation with application to the O₂ + O system

The kinetics of highly nonequilibrium gas flows, found in hypersonic flight or plasmas, cannot be characterized by single-temperature Arrhenius rates. Due to the difficulty of measuring state-specific high energy collisions experimentally, simplified phenomenological reaction and relaxation models have been developed, based on computational convenience, for nonequilibrium flow modeling. These models are calibrated to sparse equilibrium measurements and extrapolated to nonequilibrium conditions.^1–4 Phenomenological models, however, are not guaranteed to be accurate in nonequilibrium regimes.⁵ At the same time, experimental dissociation and relaxation measurements are not available at high temperatures, relevant to hypersonic flight, so the low temperature experiments are often extrapolated to flight-relevant conditions producing large uncertainties.^6,7 Accurate and efficient nonequilibrium energy exchange and dissociation models are needed in order to properly model nonequilibrium flows.

The quasi-classical trajectory (QCT) method^8,9 has been used by the chemistry community since the 1960s and by the aerospace community for over two decades to study nonequilibrium gas dynamics at high enthalpy conditions.^10–12 The QCT calculations yield state-specific reaction cross sections, state-to-state energy exchange cross sections, and scattering differential cross sections. However, there are two main challenges in adapting QCT results to flowfield simulations. First, the calculations require potential energy surfaces (PESs) for all collision pairs and rarely all potentials needed for mixture systems, such as air, are available. In an earlier work, it was shown that for high temperatures (T > 5000 K), when ab initio surfaces are not available, it is possible to construct Morse additive pairwise (MAP) potentials based on spectroscopic measurements.^13,14 The simplified potentials can produce reaction rates within experimental uncertainty bounds. The second challenge is in coupling QCT results to flow simulations using direct simulation Monte Carlo (DSMC) and computational fluid dynamics (CFD) methods. Air molecules have a high number of internal states (>10 000). QCT calculations for over 10⁴ ⋅ N_Et initial states, where N_Et is number of translational energy bins, would be required to construct a complete set of diatom-atom state-to-state cross sections, with over 10⁸ ⋅ N_Et elements. Diatom-diatom collision cross sections sets, which are more commonly needed for air flows, further grow to >10¹⁶∗N_Et elements. Since memory access is becoming more expensive than flop calculations in high-performance scientific computations, building and accessing large datasets need to be avoided.¹⁵ This paper addresses the second challenge by proposing new efficient models for internal energy exchange and dissociation based on a compact set of QCT calculations. These models can be directly implemented in DSMC codes or integrated to produce rates for CFD simulations.

There are currently four approaches available in the literature for coupling QCT results to flow calculations. In the first approach, QCT cross sections are integrated to obtain equilibrium dissociation rates that are then used with phenomenological models in DSMC simulations.¹⁶ This approach does not require any modification to the DSMC method but it does not guarantee accurate nonequilibrium results. Some researchers tabulated the complete set of state-to-state cross sections and rates, which are sampled directly within DSMC or master equation calculation.^11,17–19 Although this approach produces the highest fidelity, it cannot be easily scaled up to full flowfield simulations due to the large cross section and rate database sizes. This approach is best suited for testing fidelity of other models. QCT-DSMC hybrid codes that perform a QCT calculation, for all DSMC-selected collisions, on the fly have also been developed.^20–23 These implementations avoid large state-to-state tables at the expense of increased computational cost. QCT-DSMC hybrid codes are ideally suited for testing sensitivity of PES to reaction and relaxation phenomena, but they have not been adapted to flowfield simulations yet. Furthermore, on the fly trajectory, calculations cannot be used within CFD simulations. The fourth approach combines cross sections and rates into larger energy bins.^24–26 These bins can then be used in master and CFD calculations or sampled within DSMC. The size and internal population of bins need to be calibrated to properly reproduce global phenomena.

The present work proposes new internal energy exchange and dissociation models based on QCT cross sections and maximum entropy (ME) considerations. Maximum entropy considerations were originally proposed by Levine and Bernstein²⁷ and they have been applied to hypersonic flows before; however, ME models require an accurate surprisal functional form that was not available until now.^28–30 Previous investigators used a simpler surprisal form postulated by Levine and Bernstein and fit their models to macroscopic data. This work develops new surprisal functional form based on QCT observations. The developed ME-QCT models are fit to a reduced number of QCT calculations, and they do not require creation of massive state-to-state tables. These two advantages make it possible to upscale ab initio calculations to DSMC and CFD simulation of full flowfields. Implementation of the ME-QCT and reaction models requires only simple changes to DSMC codes and they are extendable to diatom-diatom collisions systems.

The state-specific dissociation models and the ME-QCT energy exchange models proposed in this work are fit to our previous O₂ + O collision calculations.¹³ Oxygen is found in most hypersonic flows and oxygen dissociation can dominate in the scramjet flight regimes (2-10 km/s). A survey of O₃ chemical dynamics studies performed by the chemistry community is presented by Schinke et al.;³¹ however, most of the earlier work considered gas temperatures below 1000 K. Esposito et al.,^32,33 Kulakhmetov et al.,^13,14 and Andrienko and Boyd^34,35 used the O₃ double many-body expansion (DMBE) PES by Varandas and Pais³⁶ to calculate O₂ + O vibrational relaxation and dissociation rates at temperatures up to 20 000 K. These works showed good agreement with experimentally measured equilibrium dissociation rates. Therefore, we use the same PES in our calculation. Although the new models are applied to the O₂ + O system in this work, because many high enthalpy diatom-atom collisions dynamic features are similar, the proposed models can also be applied to other chemical systems and extended to diatom-diatom systems.

The remainder of this paper is organized as follows. QCT calculations of the O₂ + O system and the ME formulation are presented in Sec. II. The ME-QCT-VT energy exchange model is then developed in Sec. III and then extended to the ME-QCT-VRT energy exchange model in Sec. IV. An O₂ + O dissociation model is then presented in Sec. V.

In this work, state-to-state energy exchange and dissociation cross sections are calculated using the quasi-classical trajectory (QCT) method. This method is described in detail by Karplus et al.⁸ and Truhlar and Muckerman,⁹ while a detailed description of our implementation can be found in our previous work.¹³ In the QCT method, initial quantum ro-vibrational energies are first related to classical position and momentum coordinates. The molecular trajectories are then propagated using the variable time step, 5th-6th order Runge-Kutta algorithm. Once the colliding particles move sufficiently far apart, the integration stops and we determine if and what type of reaction happened. For stable product molecules, we quantize final rotational and vibrational energies and place the excess energy into the post-collision translational energy, $E t ′$⁠. To construct cross sections, QCT calculations are sampled using the Monte Carlo method with varied initial impact parameters, molecular orientation angles, and phase angles. Integration time step tolerance and the maximum impact parameter are determined from a cross section convergence study. A maximum impact parameter, b_max, is maintained at a constant value of 6.5 A for all ro-vibrational and translational energies. Although a smaller b_max may be used at higher translational energies, using the larger value does not affect calculated inelastic transitions. State-to-state transition rates and state-specific reaction rates are calculated directly by Monte Carlo averaging QCT over pre-collision translational energies, E_t at a specified temperature. Therefore, cross sections and rates described later are produced by independent QCT calculations.

The double many-body expansion (DMBE) O₃ potential energy surface (PES) obtained by Varandas and Pais³⁶ is used to study $O 2 ( 3 Σ g − ) +O ( 3 P )$ collisions. The PES for a collinear collision is shown in Fig. 1. This surface is fit to ab initio calculations by Shih et al.³⁷ using self-consistent field (SCF) and configuration interaction (CI) methods, and an extended set of experimental data for ground state ozone. Most of the potential terms and parameters are found in Ref. 36. An O₂ potential surface is created by moving one of the oxygen atoms, in the O₃ potential, 100 A away. The vibrational levels, E_v(v), of O₂ have been calculated and tabulated by Esposito et al.,³³ using the Varandas and Pais potential in the limit of large O₂–O distances, using the WKB method. The ro-vibrational levels in our work are modeled as

where

The spectroscopic constants a_e, B_e, and D_e are found in Table I. The potential has 47 vibrational levels and, due to nuclear spin statistics, 2826 ro-vibrational states.

Collisions between oxygen molecules at temperatures approaching 20 000 K may experience nonadiabatic transitions and electronic excitations. Although nonadiabatic transitions can be included in QCT calculations with the trajectory surface hopping (TSH) techniques,³⁹ these techniques require a complete set of excited electronic surfaces and coupling potentials. These are currently not available for the O₂ + O chemical system. For cases when only the ground surface is available, Truhlar⁴⁰ suggested running trajectory calculations on just the ground surface but then correcting the cross sections and rates by the ratio of PES degeneracy to reactant degeneracies. Transitions via other surfaces are then assumed to not take place. For our system, this correction factor is 1/27. This correction was previously used by the authors of the PES³⁶ and Esposito and Capitelli³² in their O₂ + O dynamics studies. Tashiro and Schinke tested this approach in their formation of ozone calculations.⁴¹ The tests revealed that at room temperature, the thermal rate constant based on single surface calculations, corrected by the common electronic degeneracy factor, is comparable to the constant based on multisurface calculations within 30%. At higher temperatures, transitions to other surfaces may have a more significant effect on collision dynamics and internal excitation; however, this effect cannot be easily quantified due to the unavailability of electronically excited and coupling surfaces. In addition, more accurate oxygen dynamic calculations that better handle nonadiabatic transitions can still be fit to the new models proposed in this work, once they become available.

The dissociation cross sections and rates are corrected for possible nonadiabatic transitions by following Nikitin’s suggestions.⁴² In this approach, trajectories are again calculated using a single surface but bound-free transitions are scaled using the factor

The summation of electronic degeneracies, g_i in Eq. (3), is over all stable electronic levels with energies, E_el(i), below the ro-vibrational energy, E_rv(v, J). As is shown in Fig. 10(b), calculated dissociation rates match equilibrium experiments when this correction is applied. Additional discussion on our trajectory calculations and scaling can be found in our previous work.¹³ The state-to-state and dissociation cross sections used in the present work match previous calculations by Esposito et al..³³ Proper scaling is applied to all fitted ME-QCT and dissociation model parameters presented in this work. The reader may use the parameters directly without needing to apply further scaling.

There are approximately 8 × 10⁶ possible ro-vibrational transitions for the considered O₂ + O system. With 10 translational energy bins, the database of state-to-state cross sections grows to 80 × 10⁶, which is too large for practical implementation in CFD or DSMC. Instead, these cross sections are reduced to compact models based on the maximum entropy (ME) formalism that was originally proposed by Levine and Bernstein.²⁷ The ME theory proposes that molecular collisions try to produce post-collision distributions, P, with N nondegenerate post-collision states, j, that maximize entropy,

The distribution of post-collision states can be obtained by maximizing Eq. (4) subject to the normalizing constraint,

The Lagrange function for this constrained problem is

Solving Eq. (6) for the final states,

yields a result that does not depend on j. Therefore all nondegenerate states are equally likely. Normalizing Eq. (7) leads to P_j = 1/N.

The distribution functions for degenerate post-collision states can be obtained by considering degeneracies of these states, g(v′, E_c). In their original work, Levine and Bernstein used the translational and vibrational degeneracies of the bulk gas in equilibrium and calculated the distribution of post-collision vibrational states as²⁷ 

where E_v(v′) is the energy of the post-collision vibrational level v′ and E_c is the total collision energy that includes both translational and internal energies. This distribution can be identified as a distribution of Maxwell particles and is also referred to as the prior. When Marriott and Harvey first applied the ME model to DSMC, they used this distribution as well.²⁸ The prior described in Eq. (8) is not fundamental to the ME theory. For example, Gallis and Harvey proposed using an equilibrium distribution of the collision-pairs, which is more physically appropriate, as the prior.²⁹ This latter distribution is also known as the Larsen-Borgnakke (LB)¹ distribution for vibrational-translational (VT) energy exchange in diatom-atom collisions,

Equation (9) has an additional parameter, ν, that is related to the number of translational degrees of freedom involved in the collision. For variable hard sphere (VHS) or variable soft sphere (VSS) collision models in which collision cross sections decrease with relative translational energy, this parameter is calculated as ν = 3/2 − w, where w is the viscosity exponent that is fit to available experimental data.¹ In general, this parameter depends on the potential surfaces and it is treated as an additional fitting parameter in the present work.

Absent of additional influences molecular collisions would produce an equilibrium post-collision distribution, Eq. (9), that would maximize the entropy generation rate. Ab initio PES and details of the collisions, described in Sec. II B, introduce additional constraints on the post-collision distributions. To capture this effect, Levine and Bernstein proposed a second constraint that fixes average post-collision vibrational energy as 〈E_v〉,

where E_v(j) is the vibrational energy of a nondegenerate state j. With the additional constraint, the Lagrange function becomes

Solving Eq. (11) for nondegenerate states, P_j, yields a distribution presented in Eq. (7) multiplied by exp(λ₁E_v(j)). Therefore, the second constraint produces a degenerate post-collision distribution,

where P₀ is the prior distribution presented in Eq. (9). The first constant, λ₀, in Eq. (12) is obtained by normalizing the distribution and the 2nd constant is obtained by fixing the average energy for the final distribution. Positive λ₁ values favor higher states while negative values favor lower states. This multiplier describes how an actual post-collision distribution function, P, deviates from an equilibrium distribution, P₀. Marriott and Harvey suggested that distributions of products of endothermic and exothermic reactions can be modeled by adjusting this parameter.²⁸ 

Additional constraints show up in Eq. (12) inside the exponential term that is multiplied by the equilibrium distribution, P₀. A more general deviation from equilibrium can be described by a surprisal function,

For the model originally proposed by Levine and Bernstein, the surprisal function is S = λ₀ + λ₁E_v(v′). After analyzing available experimental and trajectory calculations, Procaccia and Levine developed a new surprisal function, $S= λ v | E v − E v ′ |$⁠, that matched the data better and satisfied micro-reversibility.⁴³ This surprisal form has been used in CFD^44,45 and master equation^46,47 calculations. However, at the time data were only available for de-excitation cross sections at bulk gas temperatures below 1000 K. As is shown later, at low collision energies, the surprisal function is linear but as collision energies are increased, this is no longer true. We extend the work of Levine and Bernstein, propose new ME-QCT models by comparing QCT-calculated cross sections, σ(v, E_c), to the equilibrium cross sections, σ₀(v, E_c), and develop new surprisal functions. This is discussed in more detail in Secs. III A and III B. The same approach can be used to extend ME-QCT models further to diatom-diatom systems and to exothermic or endothermic exchange reactions.

State-to-state energy exchange models are needed in order to properly describe flowfield relaxation. The first energy exchange model proposed in this work considers exchange between vibrational and translational (VT) modes and is referred to as ME-QCT-VT model. This model is based on the ME considerations from Sec. II, and it takes the form of a perturbed equilibrium distribution. This perturbation is described by the surprisal function, which is also introduced in Sec. II as Eq. (13). The surprisal function for the O₂ + O system is obtained from QCT calculations and is plotted in Fig. 2. The VT transition cross sections (QCT-VT), in Eq. (13), are obtained by summing QCT-calculated ro-vibrational (QCT-VRT) cross sections over all final rotational levels,

The LB model is used to produce equilibrium cross sections, σ₀(v → v′, E_c), in Eq. (13). The surprisal function, plotted in Fig. 2 for collisions with (v = 5, J = 51) initial states, shows strong favoring toward small ΔE_v transitions. This favoring drops off exponentially for ΔE_v ≤ 1 eV and linearly for ΔE_v > 1 eV. The favoring is also sensitive to translational energy. The surprisal function has higher peaks and widens as the translational energy increases. Collisions with higher translational energy typically have shorter collision duration and therefore less energy exchange.^26,48 When the total collision energy exceeds dissociation energy (orange and red lines), the surprisal function plots collapse into one.

In addition, the ME-QCT-VT model has to satisfy micro-reversibility,

In this equation, g is the state degeneracy, E_t is the translational energy involved in the collision, and σ(E_c, v, J → v′, J′) is the cross section for transition from the initial level (v, J) to the final level (v′, J′). Primed variables in this and all subsequent equations refer to post-collision states. By satisfying Eq. (15), the relaxation model will automatically achieve detail balance and relax toward the correct equilibrium distribution.⁴⁹ In the ME-QCT-VT model, the degeneracy stays constant since only v to v′ transitions are considered. Proper degeneracies are used in the ME-QCT-VRT model.

The proposed form for the ME-QCT-VT model is

where E_t is the translational energy, E_v is the vibrational energy, E_c is the collision energy, and S is the surprisal function. In the ME-QCT-VT model, the collision energy only includes vibrational and translational energies, E_c = E_t + E_v = E_t′ + E_v′. Reference cross section, σ_ref, and the exponent, ν, parameters are fit to QCT-calculated cross sections. Without the surprisal function (S = 0), Eq. (16) can be recognized as the LB model and it satisfies micro-reversibility under all conditions. To maintain micro-reversibility, the additional surprisal function has to be invariant to exchange of final and initial vibrational states (S(E_c, v, v′) = S(E_c, v′, v)). A surprisal function,

is invariant and reproduces QCT observations. This function has linear and exponential favoring terms. Following the discussion in Sec. II C, the first surprisal term constrains the VT distribution to an average vibrational energy change of 〈ΔE_v〉, while the second surprisal constrains a higher moment of the vibrational distribution function. In contrast, previous ME models only constrained average final vibrational energies, 〈E_v′〉. To capture the relationship between vibrational favoring and collision energy, the coefficients (A through C) that control the Δv favoring in the surprisal function are allowed to depend on the total collision energy,

When collision energy exceeds dissociation energy, E_c > E_d, the surprisal function, seen in Fig. 2, and state-to-state cross sections, seen in Fig. 3, should collapse into the same shape. Although it is possible to reproduce this observation by setting E_c = E_d in Eq. (16) and placing the excess energy in post-collision translational energy, doing so would break micro-reversibility in Eq. (15). Instead, we force the surprisal function to decay as E_v(v′) approaches E_d using the third surprisal function term that is still invariant to exchange of v with v′.

The ME-QCT-VT model has eleven parameters that are fit to QCT calculations. The fitting procedure is discussed in Sec. III A and their values for the O₂ + O system are presented in Table II. The new model can be used directly within DSMC calculations or integrated to produce rates for CFD calculations. For DSMC implementation, all state-to-state transitions in Eq. (16) need to be calculated for all collisions and then summed to produce the total VT transition cross section, σ_{V T,T}. The σ_{V T,T} could then be used to determine which collisions undergo VT transitions. If the surprisal function in Eq. (16) could be normalized then σ_{V T,T} could be calculated without the additional summation, improving DSMC efficiency; however, doing so breaks micro-reversibility.

The ME-QCT-VT energy exchange model, presented in Eq. (16), has eleven parameters that need to be adjusted to match QCT calculations. These parameters are optimized using the simulated annealing algorithm implemented in MATLAB with the following global objective function:

The objective function considers both the model error in cross sections, E_XSection, and error in model-predicted state-to-state rates, E_Rate. The weight assigned to the rate error is determined by the W_Rate parameter. A value of 10 is used for W_Rate in this work. The error in cross sections is determined by first calculating the log ratio between model-predicted cross sections, σ_{ME−QCT−V T}, and QCT-calculated cross sections, σ_{QCT−V T}. This difference is then summed over all final vibrational states and available initial vibrational, v_i, rotational, J_j, and translational energy, E_t,k states,

Elastic cross sections, v = v′, are not included in the summation. In Eq. (20), N_vib is the number of non-zero post-collision states. This factor ensures that the cross section fitting is not biased toward high energy collisions that have more possible post-collision states. The initial states considered in Eq. (20) are

The QCT calculation of rovibrational cross sections has been described in our previous work¹³ and it is presented in Fig. 3.

The error in the model rates is calculated as

Again, we do not consider elastic (v_i = v_j) rates in Eq. (22). Since calculating ∼2800² state-to-state rates within the parameter optimization routine would be too expensive, we only consider the rates for the following transitions:

We calculate k_{QCT−V T} rates directly within our QCT code instead of first calculating the entire set of state-to-state cross sections and then integrating the cross sections to produce rates. Although this approach requires us to run independent cross section and rate QCT calculations, because Monte-Carlo integration converges faster for high dimensional integrals than deterministic integration, we end up needing fewer total trajectory calculations to obtain k_{QCT−V T}.

Model-predicted state-specific rates are calculated by integrating model-predicted cross sections,

where X is the final state of interest and σ(v, J, E_t → X) is the cross section for the process of interest. In addition, μ_a,bc, is the reduced mass for an atom a colliding with a molecule bc, and k_b is the Boltzmann constant. For the ME-QCT-VT model, state-to-state rates are calculated by using $σ ME − QCT − V T ( v i , E t → v f ′ )$⁠, with $X= v f ′$ in Eq. (24). The optimized ME-QCT-VT model parameters are presented in Table II. The reference cross section, σ_ref in this table is divided by 27 to compensate for nonadiabatic transitions, which were discussed in Section II B. Cross sections have the units of Angstrom squared and all energies are presented in electron volts. The dissociation energy for this system is 5.212 75 eV.

The ME-QCT-VT model cross sections are compared to QCT-calculated cross sections in Fig. 3, where the former are shown with lines and the latter with symbols. QCT calculations show that mono-quantum transitions are the most likely transitions at translational energies above 1.4 eV. At lower translational energies, shown by black and blue symbols, multi-quantum de-excitations cross sections become as large as mono-quantum de-excitation cross sections. For excitation, mono-quantum transitions remain the most likely transitions; however, the value for the average quantum transition increases with translational energy. The average quantum transition from the ground state in collisions with 5 eV translational energy is 5. The model is able to reproduce major cross section features including Δv = 0 peak, decay in Δv > 0 cross sections, and relatively flat Δv < 0 region. The VT model does not consider initial rotational energy even though the true post-collision vibrational cross sections are sensitive to initial rotational levels. As the result, the model over predicts excitation cross sections, Δv ≫ 0, in Fig. 3, when it is compared against QCT-calculations with J = 11 initial level. When the total collision energy exceeds dissociation energy, as shown by green and orange squares in Fig. 3(d), the QCT-calculated post-collision Δv > 0 cross sections collapse on top of each other. The model tries to reproduce this trend but it again overpredicts the QCT calculations. The model is currently fit to collision energies that span over 2 orders of magnitude. If the ME-QCT-VT model is fit to a smaller range of translational, rotational, or vibrational energies, then the comparison between the model and QCT calculations improves significantly.

VT transition rates within 2500–20 000 K temperature range are presented in Fig. 4. De-excitation (v → v′ = 0) transitions are only weakly sensitive to temperature. The rate for the Δv = − 1 transition changes by less than 30% within this temperature range. At the same time, multi-quantum transitions are important in the de-excitation process. For instance, Δv = − 1 transition rates are only 3 times higher than Δv = − 5 transitions and only 7 times higher than Δv = − 20 transitions. The relaxation rates are high because vibrational transitions are possible via both inelastic collisions and exchange reactions. In addition, O₃ PES has a weak reaction barrier.

The quality of the ME-QCT-VT model is quantified by comparing relaxation rates it predicts to the relaxation rates calculated directly by QCT. This comparison is presented in Fig. 4. Within 5000–20 000 K temperature range, the relaxation, Δv < 0, and excitation, Δv > 0, rates predicted by the model are within 30% of rates calculated directly by QCT. The highest deviation in rates appears for large Δv excitation at low temperatures. For instance, the model-predicted and QCT-calculated rates differ by 60% for the v = 0 → v′ = 20 transition at 2500 K. However at 2500 K this transition is 10 orders of magnitude less likely than the v = 0 → v′ = 1 transition and the error is not expected to affect global flowfield calculations. The current model is sufficiently accurate for modeling flowfields but, if needed, a better comparison may be obtained by optimizing the surprisal function used in Eq. (16), the objective function used for optimization in Eq. (19), or by fitting to a smaller energy range.

Although the ME-QCT-VT model developed in Section III reproduces energy exchange between vibrational and translational modes, a comprehensive model that includes rotational modes is also necessary. In this section we propose a new ME-QCT-VRT model based on observations made from QCT calculations and ME considerations. State-to-state cross sections from the v = 5 initial vibrational level, E_t = 0.5 eV translational level and J = 11, 51 rotational levels are shown in Fig. 5. When the vibrational level does not change during the collision, Δv = 0, the most likely post-collision rotational level is the same as the initial, ΔJ = 0, and there is a strong favoring for small ΔJ transitions. However, when collisions experience vibrational energy exchange, Δv≠0, colliding molecules forget their initial rotational levels and the post-collision rotational distribution can be modeled with an equilibrium distribution,

In Eq. (25), C is a constant that scales the rotational transition to QCT results,

This supports the observation that rotational modes often equilibrate faster than vibrational modes. However, when the translational energy in collision increases post-collision rotational distributions begin to favor ΔJ = 0 despite vibrational transitions, Δv≠0. This favoring appears for translational energies ≈2.7 eV, as shown in Fig. 5(c), and it increases with increasing translational energies, as shown in Fig. 5(d). The rotational favoring leads to a slower rotational relaxation rate at higher temperatures, which was observed in other recent ab initio high temperature gas relaxation studies.^19,34 Similar favoring appears when the initial rotational level is increased as shown in Fig. 6(a), but it does not appear when the initial vibrational level is increased as shown in Fig. 6(b). In all cases, rotational favoring disappears when Δv > 5 transitions take place.

A new ME-QCT-RVT model is proposed based on the discussed observations and ME considerations,

This model considers ro-vibrational energies, E_rv, and degeneracies of post-collision states, (2J′ + 1). Although the ME-QCT-VRT model is similar to the ME-QCT-VT model, the rotational mode could not be added to the ME-QCT-VT model simply by multiplying it by the rotational distribution term. Doing so would violate micro-reversibility, Eq. (15), when rotational levels are coupled to vibrational levels. The vibrational surprisal function from the ME-QCT-VT model is not changed and coefficient A through C still depend on collision energy, E_c, as shown in Eq. (18). A new rotational surprisal function,

is introduced in the ME-QCT-VRT model. This function favors small ΔJ transitions for vibrationally elastic collisions and provides an equilibrium rotational distribution for vibrationally inelastic collisions. In the future, this model should include a smooth transitions between these two conditions as a function of translational and rotational energies. The third term in Eq. (28) forces the post-collision distribution function to decay when product internal energy approaches the dissociation energy. Parameters H and I in the rotational surprisal function also depend on collision energy as

The ME-QCT-VRT model parameters fit to the O₂ + O system are presented in Table III. This model describes approximately 80 × 10⁶ state-to-state cross sections using just 19 parameters.

The ME-QCT-VRT model uses the global objective function presented in Eq. (19) but the error in ro-vibrational cross sections is now calculated as

The new form considers rotational levels for Δv = 0 transitions. The additional parameter, W_rot, assigns weight to the rotational function and its value is set to 0.1 in this work. N_vib and N_rot in Eq. (30) are the numbers of non-zero cross sections to ensure that the fit is not biased toward large collision energies that have more non-zero cross sections. The ME-QCT-VRT model is also fit using the simulated annealing method and the optimized parameters are presented in Table III.

The vibrational and rotational cross sections fit to the ME-QCT-VRT model are shown in Fig. 7. The first column of figures has v = 0 initial vibrational level, while the second column has v = 10 initial vibrational level. The top figure of each column presents the vibrational post-collision cross sections, while the next two figures present ro-vibrational post-collision cross sections. The vibrational cross sections are obtained by summing ro-vibrational cross sections over all final rotational levels, as was shown in Eq. (14). QCT-calculated cross sections for the J = 11 initial level are shown with squares while cross sections for the J = 51 level are shown with circles. Similarly, model predictions for the J = 11 initial level are shown with dashed lines while those for the J = 51 level are shown with solid lines. Figs. 7(a) and 7(b) show that the model is able to accurately reproduce QCT-calculated vibrational cross sections and capture the sensitivity of cross sections to initial translational energy and vibrational and rotational levels. The model over predicts elastic cross sections (v′ = v); however, these cross sections do not affect relaxation rates.

The ME-QCT-VRT model also reproduces the sensitivity of rovibrational cross sections to initial translational energy and vibrational and rotational levels. Largest deviations between model-predicted and QCT-calculated cross sections are found when translational energy exceeds 3 eV. As was discussed in Section IV, rotational cross sections become increasingly sensitive to translational energy above 3 eV. The difficulty in matching the rotational cross sections is caused by the over-simplified rotational surprisal function in Eq. (28). This function can be further improved.

The vibrational relaxation rates predicted by the ME-QCT-VRT model are shown in Fig. 8. Within the 5000–15 000 K temperature range, the model-predicted and QCT-calculated rates differ by less than 60%. At 10 000 K rates differ less than 20% but at 20 000 K the rates differ by up to 120%. Currently the ME-QCT-VT model is not recommended for modeling vibrational relaxation outside the 5000–10 000 K range. The ME-QCT-VRT model, however, can be improved by using a more accurate rotational surprisal function.

A consistent dissociation model based on the same ab initio surface and QCT calculations as the relaxation model is also required. A new dissociation model is proposed in this work by assuming reaction cross sections of the form

where A, α₁, and α₂ are coefficients that are fit at fixed vibrational, v, and rotational, J, levels and linearly interpolated for intermediate levels. In the future, the linear interpolation should be replaced by a parametric function as presented for ME-QCT-VT and ME-QCT-VRT models. Doing so will reduce the number of parameters that need to be tabulated. The form for Eq. (31) is obtained by modifying the total collision energy (TCE) model¹ to comply with QCT observations. Recombination cross sections are not calculated but they can be constructed by using the micro-reversibility considerations, which were presented in Eq. (15).

The reaction model is fit using the simulated annealing algorithm with the following objective function:

where the summation is over all available translational energies that produce non-zero cross sections. The fitting is performed for initial levels,

and the fitted parameters are found in Table IV.

The dissociation cross sections for J = 11 and J = 101 levels and v = 0, 10, 20, 30 levels are shown in Fig. 9. The cross sections predicted by the model match QCT-calculated cross sections within 10% in the entire energy range and they capture both vibrational and rotational favoring. Dissociation from v > 10 levels are more than 10 times more likely than dissociation from v = 0 level. Therefore accurate relaxation models, proposed in Secs. III and IV, are needed in order to properly resolve dissociation rates in a flowfield.

The quality of the reaction model is tested further by calculating state-specific and equilibrium reaction rates. The state specific reaction rates are shown in Fig. 10(a). In these plots, the rotational and translational temperatures are equilibrated although the model does not require this. The vibrational levels, for which the model was fit, are shown with red lines in the same figure. Even though the reaction model is fit to only 5 vibrational levels, the new model is able to reproduce state-specific rates for vibrational levels up to v = 40 by interpolating between available points. The model under predicts rates above the 40th vibrational level because higher levels were not fit. Although state-specific rates for levels above 40 are high, these levels have low equilibrium populations and do not contribute significantly to equilibrium dissociation rates. In compression nonequilibrium flows, the highest levels are even further underpopulated. Of course, the implementation can be improved by including higher levels in the fit.

Equilibrium reaction rates calculated by the new model are shown in blue in Fig. 10(b). The new model is able to reproduce rates calculated directly by QCT, shown by orange squares, within 25% at 2500 K and within 10% above 10 000 K. This is a good agreement considering that the reaction model is fit to only 160 cross sections (compared to ∼2800² initial states that were considered in the rates calculated directly by QCT). Both the model and direct-QCT rates fall within the range of available experimental measurements made by Ibraguimova⁵⁰ and Shatalov.³³ In contrast to the 10% disagreement between the model and QCT calculations, the experimental data disagree from each other by 300% at 10 000 K. Nonequilibrium dissociation measurements for the O₂ + O → 3O reaction do not exist; however, the popular two-temperature phenomenological model⁶ can disagree from QCT-calculated rates by a few orders of magnitude at the vibrationally cold conditions.¹³ 

Recent quasi-classical trajectory (QCT) calculations began to probe non-equilibrium phenomena found in hypersonic and other energetic flows. In this work, we adapt QCT calculations to study energy exchange and dissociation in the O₂ + O collision system. We observed that collisions favor small Δv transitions with translational energies below 1.4 eV. The favoring drops off at first exponentially and then linearly as Δv increases. Based on these observations and maximum entropy (ME) considerations, we develop the ME-QCT-VT vibrational-translational energy exchange model. The model is then fit to the O₂ + O system using only 11 parameters. The state-to-state vibrational rates calculated by this model reproduce rates calculated directly by QCT within 30% in the 5000 K to 20 000 K temperature range.

Quasi-classical trajectory calculations show that in absence of vibrational energy transfer, post-collision rotational energies favor small ΔJ transitions. When large vibrational energy exchange takes place, collisions forget their original J state and post-collision rotational levels follow the equilibrium distribution. The rate at which favoring decreases depends on pre-collision translational and rotational energies. Based on these observations and ME considerations, we propose a new ME-QCT-VRT vibrational-rotational-translational energy exchange model. The model reproduced main features of post-collision ro-vibrational distribution functions and is able to predict vibrational excitation and relaxation rates within 60%. Methods to improve this model further are also discussed.

In addition to the energy exchange models, we develop a new dissociation model based on observations made in QCT calculations. Dissociation cross sections and rates show strong vibrational and rotational favoring. For the same total collision energy, dissociation cross sections increase by over an order of magnitude when vibrational levels increase. The new model is fit to QCT-calculated O₂ + O dissociation cross sections at discrete ro-vibrational levels and interpolated between these levels. The model reproduces main state-specific dissociation cross section and rate features. The equilibrium rates predicted by the model and direct-QCT calculation match within 25% in the 2500 K to 20 000 K temperature range.

Although this work applies the energy exchange and dissociation models to the O₂ + O system, the proposed models can be applied to other diatom-atom systems. There is also a clear path for extending ME-QCT models to diatom-diatom collision systems. The models reduce the number of required QCT calculations and parameters that need to be stored. This approach makes it possible to upscale ab initio calculations to full flowfield DSMC and CFD simulations. All parameters needed to model internal energy exchange and dissociation for the O₂ + O system are presented in this work.

This work has been supported by Sandia National Laboratory through the Sandia Excellence in Science and Engineering Research Fellowship. In particular, we would like to acknowledge Dr. Dan Rader at Sandia. We would also like to thank Professor Antonio Varandas at University of Coimbra for providing an O₃ potential energy surface.