Rovibrational internal energy transfer and dissociation of high-temperature oxygen mixture

Behind strong shock waves, diatom–diatom collisions are essential as they initiate the dissociation process and determine the chemical production rates of atomic species. Once a significant number of atoms are created, atom–diatom processes govern the thermochemical nonequilibrium state of the gas.^1–4 In hypersonic flows within Earth’s atmosphere, O₂ + O₂,^5–7 O₂ + N₂,^8,9 and N₂ + N₂¹⁰ are the relevant diatom–diatom systems. Among them, O₂ + O₂ interaction is the focus of the present work owing to its dominant contribution to the early stages of thermochemical relaxation, justified by the weak bond energy of molecular oxygen.

The O₂ + O₂ system has been investigated extensively both experimentally^7,11 and numerically.^5,6,12–16 Recently, shock tube experiments were carried out by Ibraguimova et al.¹¹ and Streicher et al.⁷ to measure vibrational relaxation times and nonequilibrium dissociation rate coefficients. Those measurements are mainly used as a reference for validating models. Theoretical calculations^5,6,12–16 of the O₂ + O₂ system include approximate and high-fidelity estimates. Lino da Silva et al.¹² proposed a vibrational-specific (VS) state-to-state (StS) model constructed using the Forced Harmonic Oscillator (FHO) theory and one-dimensional Morse intermolecular potential. Andrienko and Boyd¹³ devised another VS model using the quasi-classical trajectory (QCT) method. In their work,¹³ a double many-body expansion (DMBE) potential energy surface (PES) developed by Varandas and Pais¹⁴ was used to describe molecular interactions. The major shortcoming of these studies^12,13 is the crude treatment of rotational degrees of freedom, assumed to be in equilibrium at the translational energy. Previous studies^1,2,17 demonstrated that rotational nonequilibrium is essential to model hypersonic shock layers accurately. Grover et al.⁶ performed direct molecular simulation (DMS) using the ab initio PESs proposed by Paukku et al.^15,16 The results of Grover et al. accounted for the aforementioned rotational nonequilibrium and showed good agreement with the available experimental data. Although the DMS captured general aspects of nonequilibrium kinetics in the O₂ + O₂ system, the method could not capture the details of the rovibrational-specific population dynamics, which is still unknown to the author’s best knowledge. Understanding the rovibrational distributions and their dynamics in time is crucial to developing an advanced reduced-order model for nonequilibrium thermochemistry.^2,18

Numerical simulation of rovibrational-specific molecular interactions in the O₂ + O₂ system is computationally prohibitive. Several efforts have been made to model the rovibrational energy transfer and nonequilibrium dissociation occurring in four-atom systems. Kim and Boyd¹⁹ employed an ordinary Kriging model to estimate fully rovibrational-specific StS rate coefficients in the H₂ + H₂ system. In another paper, Kim et al.²⁰ proposed an approximated treatment of the internal states, assumed to be populated according to a Boltzmann distribution of one of the molecules in the H₂ + H₂ system to reduce the computational complexity. Macdonald et al.¹⁰ developed a coarse-grain QCT (CG-QCT) method to investigate nonequilibrium rovibrational energy transfer and dissociation in N₂ + N₂ systems by relying on a binning strategy among internal states. Following Kim et al.,²⁰ the present work proposes to model rovibrational energy transfer and dissociation in the O₂ + O₂ system, assuming a thermal distribution for a collision partner.

Toward this end, the present study proposes a hybrid strategy to model the rovibrational-specific interactions in the O₂ + O₂ system by extending the modeling concepts by Kim et al.²⁰ This procedure imposes a single energy bin on the colliding molecule that is Boltzmann distributed at a prescribed internal temperature. At the same time, the target O₂ contains complete information on the rovibrational states. The present model is built on the framework consisting of the QCT method and master equation analyses. The QCT simulations employ the ab initio PESs used in the work of Paukku et al.^15,16 to construct a rovibrational-specific StS kinetic database. Then, the master equation analysis was performed in the kinetic temperature range of 7500–20 000 K. The rovibrational-specific population dynamics, energy transfer, and dissociation are thoroughly investigated to understand the fundamental physics behind the O₂ + O₂ collisional system. A direct comparison with a VS model, which this work has also constructed, is made. In addition, the underlying physics in the rovibrational distributions of the O₂ + O₂ system is directly compared with that of the O₂ + O system to make clear the differences between these systems. Finally, the master equation study is extended to combined O₂ + O and O₂ + O₂ systems (i.e., O₂ + O and O₂ + O₂) to identify their relative contributions.

This paper is organized as follows: A cost-effective modeling strategy to evaluate rovibrational distributions in the O₂ + O₂ system is proposed in Sec. II. In Sec. III A, the constructed StS kinetic database is presented. Then, the rovibrational energy transfer and nonequilibrium dissociation are discussed in Secs. III B and III C, respectively. Section III D provides the master equation results for the combined O₂ + O and O₂ + O₂ systems. In Sec. III E, the present macroscopic results are compared with existing data from the literature. Finally, Sec. IV provides the conclusions from the present study.

This work assumes that all species relevant to the O₂ + O₂ system remain in their electronic ground state during collisional interaction. Then, the rovibrational energy transfer and dissociation processes of the O₂ + O₂ system can be expressed as

where the indices (i, m, j, l) denote the rovibrational states of $O2X3Σg−$ that are stored by increasing energy in the set $IO2$⁠. Each rovibrational level corresponds to a unique vibrational (v) and rotational (J) quantum number [e.g., i = i(v, J)]. The reaction in Eq. (1) includes both inelastic and exchange kinetics. The processes in Eqs. (2)–(4) contain dissociation from both direct and exchanged pairs. In the present master equation study, the double-dissociation reaction in Eq. (4) is neglected since it was found to have little influence for the conditions of present interest.²¹ Then, the system of master equations describing the evolution of the rovibrational states of O₂ and the concentration of O atom reads

Here, n_X identifies the number density of X species/level and t denotes time. The rate coefficient k describes the corresponding bound–bound (e.g., k_im→jl) and bound-free (e.g., k_im→cl) transitions where the subscript c denotes the dissociated state. It is important to note that the rate coefficient k is a dependent variable of the kinetic temperature T, although its dependency is omitted for brevity.

Description of the rovibrational energy transfer and dissociation of O₂ using Eqs. (5) and (6) requires a prohibitive amount of computational cost since $O2X3Σg−$ has 6115 different rovibrational levels (see Sec. II B), resulting in about 6115⁴ collisional processes to be considered. This constraint in the computational cost motivated the development of coarse-graining strategies.^{2,10,18,22,23} However, those order reduction strategies provide an approximated description of the StS population dynamics rather than the complete rovibrational-specific distributions, which are crucial ingredients in developing advanced reduced-order models.^2,18,24 This work, thus, presents a novel and simple reduced-order model, which allows investigating the rovibrational distribution in the O₂ + O₂ system by combining the concept of coarse-grain QCT method¹⁰ with the model proposed by Kim et al.²⁰ assuming a Boltzmann distribution among the collision partner’s internal states.

Derivation of the proposed model can start from Eqs. (1)–(4). We enforced the rovibrational states m of the collision partner being populated as a Boltzmann distribution at a prescribed internal temperature T_int by introducing a state-specific distribution function $Fm$⁠,

The partition functions Q_m and $QO2$ are defined as

where e_m and k_B denote the mth state energy and the Boltzmann constant, respectively. It is worth mentioning that the physical meaning of T_int corresponds to the rovibrational temperature, T_RV, since the excitation of electronic state is excluded. The rovibrational degeneracy $gm=gnuc2J+1$ includes the nuclear spin contribution as g_nuc = 0.5 for all the levels based on a semiclassical approximation.² Then, Eqs. (1)–(4) can be recast as

where O₂(T_int) denotes the collision partner Boltzmann-populated at T_int. It is important to note that Fig. S1 in the supplementary material verifies the insignificant impact of the double-dissociation in Eq. (12) by comparing the thermal dissociation rate coefficients from the mechanisms in Eqs. (10)–(12). By neglecting the double-dissociation, the set of master equations corresponding to Eqs. (9)–(11) leads to

The form of rate coefficients k_i→j, k_i→c, and $ki→j⋆$ can be defined as

where k_i→j, k_i→c, and $ki→j⋆$ correspond to the reactions in Eqs. (9)–(11), respectively. It is important to note that dependence of the rate coefficients on collision partner’s internal states (i.e., m and l) is eliminated by assuming that the collision partner is thermalized at the given T_int, represented by $FmTint$ for the mth level for example. This implies that the computational cost for the $O2X3Σg−$ + $O2X3Σg−$ system for the given T_int and T can be reduced to the same order of that for the $O2X3Σg−$ + $OP23$system. Consequently, solutions from Eqs. (13) and (14) describe the rovibrational energy transfer and dissociation in $O2X3Σg−$ + $O2X3Σg−$ system with much less computational expense compared with the exact formulation in Eqs. (5) and (6).

The rovibrational-specific state-to-state rate coefficients in Eqs. (15)–(17) are computed by means of the QCT method discussed in Sec. II B. At a given T, the QCT calculations are carried out for several T_int that range from 300 K to T (see Table I for more details). To integrate the set of master equations in Eqs. (13) and (14), the rate coefficients k and k^⋆ depending on the T_int are updated at each time step by interpolating them in T_int-space. For this purpose, T_int is evaluated by equating the internal energy extracted from the computed population to an equivalent average energy as defined in  Appendix C. Figure S4 shows an example of the evolution of the so-calculated T_int in time at T = 10 000 K. Furthermore, a linear interpolation is subsequently carried out for log scale of the rate coefficients over the T_int-space. This strategy is adopted to perform the master equation study at T = 10 000 K, whereas for other T the rate coefficients are averaged over the T_int-space and kept as constant values for further reduction of computational costs. It is important to note that the latter averaging approach does not affect to the key findings of the present work. The aforementioned reduced-order modeling strategy is denoted as thermal collider (TC) approach in the following sections.

Integration of the master equations is carried out using plato (PLAsmas in Thermodynamic nOn-equilibrium),²⁵ a library for nonequilibrium plasmas developed within the Center for Hypersonics and Entry System Studies (CHESS) at the University of Illinois at Urbana-Champaign. In the master equation analysis to be presented, the chemical composition is initialized to have 100% of O₂. The initial gas pressure and internal temperature are set to 1000 Pa and 300 K, respectively. Those conditions are applied for all simulations unless otherwise stated. As per the heat-bath temperatures, the values of 7500, 10 000, 15 000, and 20 000 K are considered.

To compute the rovibrational-specific StS rate coefficients for Eqs. (9)–(12), QCT calculations are performed using CoarseAIR.^2,26 It is an in-house QCT code that is a modernized version of the original VVTC code developed at NASA Ames Research Center by Schwenke.²⁷ The three different ab initio PESs proposed by Paukku et al.^15,16 are employed to consider the adiabatic interactions occurring in the electronic ground state of the O₂ + O₂ system through the singlet, triplet, and quintet spin couplings. The statistical weights for each spin interaction are approximated as 1/9, 3/9, and 5/9, respectively, based on the high-temperature limit. The 6115 rovibrational energy levels are obtained by solving Schr$ö$dinger’s equation based on the WKB semiclassical approximation.^27,28 The resultant dissociation energy of the rovibrational ground state is found as 5.113 eV. The initial state of thermal collision partner $O2Tint$ is sampled from the prescribed distribution function $FmTint$⁠. Consequently, the rate coefficients in Eqs. (15)–(17) are computed by summing up the statistics of individual trajectories with respect to the final state of the collision partner.

For a given bound–bound transition, both endothermic and exothermic trajectories can exist as the result of QCT calculations. In this case, one of them must be selected as a reference direction and the other one needs to be reconstructed based on micro-reversibility to guarantee appropriate equilibration of the initial heat-bath simulations. In the present QCT calculations, the sampled trajectories are investigated to decide which direction (i.e., endothermic or exothermic) has better statistics as presented in Sec. II A, where a detailed discussion is provided. As a result, we utilize the endothermic sampled trajectories as the reference forward direction with higher priority than the exothermic ones to construct the StS kinetic database described in Sec. III A. Exceptionally, the exothermic trajectory is used in a case that the given reaction channel does not have the endothermic trajectory. The reverse rate coefficients are then reconstructed based on micro-reversibility. Influence of the choice of reference direction can be found in Fig. S16 in terms of the vibrational and rotational relaxation times. In the considered temperature range, the maximum variations due to the choice of direction occur at the lowest temperature. They are within factors of 2.2 and 1.4 for the vibrational and rotational relaxation times, respectively.

The constructed StS kinetic database in this work covers the kinetic temperature range 7500–20 000 K. The combination of T and T_int for the QCT calculations are summarized in Table I. The present analysis at T = 10 000 K ensures the highest fidelity in terms of accuracy since it covers the largest number of QCT databases along the T_int-space. The results in this work at other T, such as 7500, 15 000, and 20 000 K, demonstrate the extension of the present modeling strategy to various kinetic temperature ranges. This section investigates the influence of the sampling temperature T_int on the constructed StS kinetic database, especially for the inelastic and target particle’s dissociation of Eqs. (15) and (16), respectively.

Figure 1 shows the inelastic rate coefficient distributions at T = 10 000 K. The magnitude of rates is overlayed on the diatomic potential of $O2X3Σg−$⁠. The overall trends observed from both T_int are qualitatively similar to each other: In the low-lying states, transition pathways are mostly aligned to the same vibrational strands, whereas they are spread across different vibrational levels in the high-lying energy regime. This may imply that the assumption for separation between rotational and vibrational states is less feasible in the high-energy area. A distinct difference according to the variation of T_int is the rate coefficient for vibrational ladders. At higher T_int, the vibrational transition rates are larger particularly in the low-v and low-J region that governs majority of the rovibrational energy transfer. This difference is quantitatively highlighted in Fig. 2. The rotational transition shown in Fig. 2(a) is not sensitive to variations in T_int. The magnitudes of the rotational transition rates are almost identical to each other. On the other hand, the vibrational excitation rates change by more than an order of magnitude with variations in T_int as observed in Fig. 2(b). As T_int increases, the vibrational energy transfer is significantly raised. This increase starts to be asymptotic around T_int = 5000 K. A similar trend is consistently observed for other initial target states as shown in Figs. S6 and S7. From the investigations in Figs. 1 and 2, it is found that the sampling temperature T_int has a significant influence on the vibrational transition rates, whereas their impact on the rotational transitions is insignificant.

The target particle’s dissociation rate coefficient distributions are presented in Fig. 3. The figure compares the rates calculated at the three different T_int. The dependence of the dissociation rate coefficient on T_int is a function of the distance from the centrifugal barrier. This behavior was observed by Venturi et al.² In the case of T_int around 10 000 K, the rate coefficients for O₂ + O₂ are indistinguishable from the ones computed for the O₂ + O system (See Fig. S8).

The influence of variations in T_int is confined to the low-lying energy levels, which do not have a significant impact on the dissociation kinetics. By relating the rates with the concept of energy deficit,² defined as a distance from the centrifugal barrier in energy space, it is clear that the T_int of the colliding molecule works as a stretching factor for the magnitude of rates, especially for states far away from the centrifugal barrier. This is because as the collision partner has higher internal energy (i.e., higher T_int), the low-lying states of target O₂ can more easily gain the required barrier energy to dissociate during collisions.

It is important to note that in the high-lying energy states, the rates are very close to each other and, consequently, the difference in thermal rate coefficients between T_int = 300 and 10 000 K is less than a factor of 1.5 at T = 10 000 K. Figure S9 supports this by showing a more quantitative distribution of the rates for selected quantum pairs.

In this section, internal energy transfer by inelastic and exchange kinetics is studied in detail by neglecting the dissociation mechanisms while solving the set of master equations. The influence of sampling temperature T_int on the rovibrational population dynamics by the inelastic process is first investigated at T = 10 000 K. Then, the master equation analysis is extended to other kinetic temperatures along with examination of the exchange process.

Figure 4 presents snapshots of the rovibrational distributions by inelastic process at t = 10⁻⁷ s and T = 10 000 K. The time instant 10⁻⁷ s is selected to compare the distributions in the early stage of energy transfer. Overall, structures at the different T_int values are similar to each other. At the higher T_int, the excited vibrational levels are more populated than can be verified from the slope among the heads of vibrational strands which implies the vibrational temperature T_V. At the same time instant, the rate of average vibrational energy transfer at T_int = 10 000 K occurs about two times faster compared to that at T_int = 300 K. At both T_int, the thermalization process induced by the inelastic kinetics slowly occurs for the high-lying bound levels around the dissociation threshold (i.e., high-v and low-J states). This is a part of direct consequences of the dissociative properties of the molecule and has been consistently observed in previous studies.^1,2,10 The rovibrational levels near the threshold are characterized by their significantly high dissociation probability as shown in Figs. S5 and S8. This results in that their populations are much lower than those of the other rovibrational levels during relaxation via inelastic transition.

The master equation study is extended to other kinetic temperatures to investigate global trend of rovibrational energy transfers in the O₂ + O₂ system. Figure 5 shows the temporal evolution of vibrational (T_V) and rotational (T_R) temperatures with the variation of T. From the master equation results, the characteristic rotational–translational (τ_RT) and vibrational–translational (τ_VT) relaxation times are determined using the e-folding method.²⁹ Then, they are fed into the linear Landau–Teller (L–T) relaxation equation³⁰ to compare with the internal energy transfers from the master equation study. The L–T relaxation equation can be expressed as

where E_V and E_R are the average energy of vibration and rotation, respectively. The subscript eq indicates the quantity evaluated at the equilibrium temperature T. It is important to note that the L–T model reasonably reproduces the master equation results in the considered T range as shown in the figure. Although the difference between the vibrational and rotational characteristic relaxation times is reduced as T increases (See Fig. 16), the rotational energy transfer is distinctly faster than the vibration in the considered temperature range.

The internal energy transfer processes in O₂ + O₂ system are analyzed in detail by investigating the time evolution of rovibrational distributions as shown in Fig. 6. The figure presents the distributions at T = 10 000 K. In earlier times, the low-lying levels form several distinct vibrational strands whereas the high-lying states near the dissociation limit are preferentially excited than the others around (see t = 1 × 10⁻⁸ and t = 1 × 10⁻⁷ s) due to the exchange kinetics that can be inferred by comparing Figs. 4 and 6 for the time instant t = 1 × 10⁻⁷ s. The quasi-bound levels having higher energy than 6.5 eV are quickly equilibrated to each other that forms a Boltzmann distribution in the tail region. From t = 1 × 10⁻⁸ to t = 1 × 10⁻⁷ s, the excitation of high-lying states is saturated and thermalization among the low-lying levels is slowly followed. One of important findings of the present study is that in the O₂ + O₂ system, the thermalization among different vibrational levels is relatively slow compared to that of the O₂ + O system. This finding is attributed to two facts: In O₂ + O₂, (1) the inelastic process occurs in more rotation-favorable ways rather than vibration and (2) the exchange kinetics is much less effective in comparison with those of the O₂ + O system. The aforementioned points are discussed in detail through Figs. 7–9.

Figure 7 highlights the influence of exchange kinetics on the internal energy transfer of the O₂ + O₂ and O₂ + O systems. The lower and upper bounds of T are considered. Results obtained from a present VS model for the O₂ + O₂ system are also presented for comparison. The rovibrational-specific kinetic database for the O₂ + O system used in the present study is taken from the work by Venturi et al.² Obviously, the exchange kinetics has little impact on the internal energy transfer of the O₂ + O₂ system, whereas it significantly hastens the relaxation of the O₂ + O system. This fact results in different rates of thermalization among the internal energy states. In the O₂ + O₂ system, the equilibration among different vibrational strands slowly occurs as observed from Fig. 6. On the contrary, the exchange process advances the mixing of the internal levels in the O₂ + O system. Figure S10 presents more details of this in terms of the rovibrational distributions of the O₂ + O system. The exchange reaction more favorably occurs in the O₂ + O system compared to the O₂ + O₂ system since the three-atom system has no barrier for atom exchange.^31–33 

Figure 8 provides an alternative interpretation regarding the results shown in Fig. 7 on a more microscopic scale. The figure compares the rovibrational distributions in the O₂ + O₂ and O₂ + O systems at T = 10 000 K. Both distributions are extracted from 20% of vibrational relaxation for direct comparison. The two systems commonly have strand-like structures in their low-lying energy regime. In the O₂ + O₂ system, impact of the exchange process is only confined to levels near the dissociation limit that can be inferred by comparing Figs. 4 and 8 (top). On the other hand, in the O₂ + O system most of the vibrational strands are already thermalized, except for some low-lying levels having the state energy less than 1.5 eV. The result shown in Fig. 8 implies that the fundamental structure of the rovibrational distributions between the two systems is different. It is important to note that the difference observed in Fig. 8 is not only due to exchange kinetics but also attributed to the inelastic process. Even though with the inelastic kinetics only, the rovibrational distributions of O₂ + O₂ are more broadly spread than those of O₂ + O, especially in the low-lying energy region in which the energy transfer is dominated by the inelastic kinetics.

The different natures in the inelastic kinetics can be more effectively distinguished by investigating the vibrational transition probability $Pvi→vj$ due to the inelastic process, where v_i and v_j are the initial and final vibrational quantum numbers of the target O₂, respectively. Figure 9 shows the distributions of $Pvi→vj$ in the O₂ + O₂ and O₂ + O systems at T = 10 000 K. The probability is computed by summing up the individual probabilities, which are obtained from the QCT calculations, along with the initial and final rotational states. In O₂ + O₂, the $Pvi→vj$ distribution is more clustered around the diagonal line (i.e., v_i = v_j) compared to O₂ + O, implying that the vibrational transition is less favorable. In the same vein, distributions of the rotational transition probability $PJi→Jj$ are also provided to Fig. S11, confirming that the O₂ + O system has more diagonally clustered probability, resulting in a faster rate of vibrational transitions through the inelastic process.

In this section, dissociation and recombination processes in the O₂ + O₂ system are investigated. Figure 10 shows temporal evolution of the internal temperatures and O₂ mole fraction with T. In the considered T range, the plateau region of the temperatures that implies the quasi-steady-state (QSS) period is clearly observed as shown in Fig. 10(a). In the QSS region, the rotational and vibrational temperatures are distinctly separated from each other, especially at T = 15 000 and 20 000 K. The slow thermalization among the vibrational strands discussed in Fig. 8 contributes to that separation by allowing the two energy modes to have different gradients of population distribution (see Fig. S12 for more details). As shown in Fig. 10(b), it is confirmed that O₂ is completely dissociated in the given T range. It is important to note that in Eq. (11), the influence of combined inelastic and dissociation process is negligible, whereas the contribution due to combined exchange and dissociation kinetics is substantial. This is because when exchange occurs, in O₂ + O₂ collision, bond breaking happens, resulting in the higher probability of the dissociation. In Fig. S13, the O₂ composition profiles obtained from the TC model are compared with those obtained from the VS approach. This demonstrates that the two methods show reasonable agreement in the considered T range.

Accuracy of the present modeling strategy for O₂ + O₂ dissociation can be verified by comparing QSS vibrational populations with the existing result obtained from the DMS by Grover et al.⁶ The normalized cumulative vibrational population f_v is defined as

where J_max implies the maximum rotational quantum number at the given v. It is important to note that comparison of f_v can function as a necessary, but not a sufficient, condition.^2,34

Figure 11 shows the comparison of f_v at QSS for T = 10 000 K. The distributions are computed from the three different models. In the QSS region, the low-lying vibrational levels form a Boltzmann distribution whereas the high-lying ones are in strong non-Boltzmann state. The dissociation process is mostly governed by the high energy levels, confirming the O₂ + O₂ dissociation occurs in a nonequilibrium manner during the QSS period. As shown in the figure, the results from the different models are in reasonable agreement, which demonstrates the capability of the present TC approach to predict the QSS vibrational distribution of the molecule.

The rovibrational distribution in the molecular QSS period is one of the most important ingredients for studying physics behind dissociation kinetics. In this work, the QSS rovibrational distribution of O₂ + O₂ is compared with that of O₂ + O as shown in Fig. 12 for T = 10 000 K. By comparing Figs. 12(a) and 12(b), it is found that structures of the distributions are different from each other; overall, the O₂ + O₂ system has a more scattered rovibrational population. This difference is mainly attributed to dissimilarity in the behavior of energy transfers (i.e., inelastic and exchange) rather than the dissociation dynamics. As discussed in Sec. III B, the O₂ + O₂ system represents a significantly slower thermalization rate among the internal states, especially for the vibrational strands. This ends up as the more scattered population distribution at the time the system enters the QSS period in which majority of O₂ molecule dissociates.

In Figs. 12(c) and 12(d), the QSS distributions are compared with those from the VS model. The rotational populations predicted by the VS model are reconstructed by assuming a Boltzmann distribution at T. The comparisons are made in the middle of the QSS regions, determined by the TC model in the O₂ + O₂ system and by the fully StS approach in the O₂ + O system, respectively. The inset figures present normalized O₂ mole fractions, where $χO2,0$ implies the initial value of $χO2$⁠. In O₂ + O₂, the QSS distributions are well matched for the bound levels, whereas substantial discrepancy exists for the quasi-bound states. The VS model overestimates the quasi-bound populations. Despite that difference, both models predict almost identical mole fraction profiles. In O₂ + O, the mole fraction evolution calculated by the VS model is distinctly different from that calculated by the fully StS approach: The VS model overestimates the amount of dissociation. In a previous study by Venturi et al.,² a similar issue was addressed and an improved coarse-grain model for the O₂ + O dissociation kinetics was designed by leveraging the physics behind the diatomic potential and the centrifugal barrier. Since the faster thermalization quickly squeezes the rovibrational distribution, the high-lying levels around the dissociation limit have relatively larger contributions to the amount of dissociation in O₂ + O, compared with O₂ + O₂ in which the equilibration among the vibrational states is much slower in nature.

In Fig. 13, the normalized rovibrational energy loss ratio ϵ_i/ϵ_I is compared between O₂ + O₂ and O₂ + O systems. The state-specific internal energy loss ratio ϵ_i is calculated using the definition proposed by Panesi et al.¹, where ϵ_I is the accumulated quantity of ϵ_i over the internal states. By using the present TC model for the O₂ + O₂ system, the ϵ_i can only be defined for the dissociation process in Eq. (10) but not for Eq. (11). This is because the present TC model has lost access to the state-specific energy loss ratio regarding the dissociation of the collision partner while applying the thermalized assumption to Eq. (3) to obtain Eq. (11). This is one of the drawbacks of the present TC approach. In O₂ + O₂, the distribution has a more gradual shape, whereas it shows a spike-like structure in O₂ + O centered around the dissociation limit. From the comparison, it is found that in O₂ + O₂, the quasi-bound states have much less contribution to the dissociation-energy coupling than that of O₂ + O. This is attributed to lower amount of rovibrational populations of the quasi-bound levels in O₂ + O₂ as discussed in Fig. 12.

In this section, the master equation analysis is extended to complete the oxygen mixture that includes both O₂ + O₂ and O₂ + O interactions. The purpose of this extension is to investigate relative contributions from those two collisional systems to the internal energy transfer and dissociation of O₂ in the full oxygen mixture. The set of master equations corresponding to the combined O₂ + O and O₂ + O₂ system is presented in  Appendix B. The present results are compared with those of a study by Grover et al.⁶ obtained from the DMS method. For this comparison, the chemical components are initialized to have 100% of O₂ and initial number density of 2.414 × 10¹⁹ cm⁻³. The initial rotational and vibrational temperatures are set to 1394.3 and 985.1 K, respectively. Then, the kinetic temperature T is suddenly raised to 10 000 K and kept constant to study the nonequilibrium energy transfer and dissociation processes.

Figure 14 presents the normalized average internal energy and O₂ mole fraction profiles. It is important to note that for the present results shown in the figure, the recombination process is turned off for consistency with the DMS result obtained by Grover et al.⁶ Among the results, Present-Rovib implies a rovibrational-specific approach that is a hybrid of the TC model for O₂ + O₂ and the fully StS method for O₂ + O. For both average vibrational energy and O₂ mole fraction profiles, the present result obtained by the rovibrational-specific approach is in better agreement than the VS model compared to the reference DMS data.⁶ This improvement in the accuracy is mostly attributed to the rovibrational-specific treatment of the O₂ + O rather than O₂ + O₂ system since we observed that the VS and the TC approaches provide very similar macroscopic mole fraction profiles as shown in Fig. 12(c). The present average rotational energy profile shows a larger discrepancy from the DMS result than the vibrational energy profile during nonequilibrium energy transfer. According to the previous studies conducted by Venturi et al.² and Macdonald et al.,³⁵ it was demonstrated that the StS master equation approach yields very similar results to those obtained from the DMS provided that an identical PES is adopted. Therefore, the discrepancy observed in the average rotational energy profile in Fig. 14(a) requires further investigation.

To the authors’ best knowledge, so far no detailed rovibrational distribution has been presented for the complete oxygen system. In Fig. 15, the evolution of rovibrational distributions and the time-cumulative dissociation production rate are presented. It should be noted that the recombination is now turned on to simulate the results shown in Fig. 15. Figure 15(a) confirms that the O₂ + O₂ collision governs the early stage of energy transfer since the rovibrational distribution at 10% of vibrational relaxation represents a structure very similar to that of isolated O₂ + O₂ shown in Fig. 8. In the QSS, however, the O₂ + O collision takes control of the oxygen mixture. This can be inferred from the fact that the QSS distribution shown in Fig. 15(a) is almost identical to that of isolated the O₂ + O system shown in Fig. 12(b). One more important takeaway from Fig. 15(a) is that the QSS distribution of the complete oxygen system (i.e., bottom) is clearly correlated with each rovibrational level’s energy deficit $eiD$⁠. This implies that the coarse-grain strategy,² which was developed for diatomic dissociation by correlating the energy deficit with the multi-group maximum entropy approach,³⁶ can also be applied to the case of a complete oxygen mixture.

The quantity presented in Fig. 15(b) is defined by summing up and normalizing the mass production rate for the dissociation and recombination in Eq. (13) at given time step with the total amount of chemical production rate relevant to dissociation and recombination. This definition was taken from the previous study by Jo et al.³ From the figure, it is found that before the onset of QSS most of the dissociation is dominated by O₂ + O₂ collision. During the QSS period, about 35% of dissociation occurs in which O₂ + O system has larger contribution than O₂ + O₂. This observation is in line with the rovibrational distribution in Fig. 15(a). It is important to note that in total, O₂ + O₂ collision has a larger contribution to the cumulative dissociation production rate. This is because the four-atom system has more varied dissociation pathways compared to the three-atom system.

In this section, the macroscopic quantities, such as rovibrational relaxation time, dissociation rate coefficients, and average energy loss ratio, obtained from the present analyses are compared with existing data.

Figure 16 compares the rotational and vibrational relaxation times of the O₂ + O₂ system with the literature data.^{6,7,11,37–41} The symbol p denotes the pressure of the collision partner. The present vibrational relaxation time obtained by the TC model is in good agreement with the DMS data.⁶ The result obtained from the VS approach, however, shows the distinctly shorter relaxation time at T = 15 000 and 20 000 K and represents closer agreement with the Park model, including high-temperature correction.⁴⁰ In comparison with the measured data,^7,11,37 the present TC result shows slightly improved agreement than the VS model. It is worth mentioning that the pτ_VT data measured from the recent experiments^7,11 represent longer relaxation time than the present results. The rotational relaxation time pτ_RT calculated by the present TC model shows significant discrepancy compared to that obtained from the Parker model,³⁸ which was derived by classical mechanics using the rigid-rotor approximation. From the present result obtained from the TC model, it is found that in O₂ + O₂ collision the convergence of the rotational and vibrational relaxation time to a common asymptote in high-temperature range⁴² is relatively weak in comparison to other three-atom systems, such as N₂ + N,^1,43 O₂ + O,⁴⁴ N₂ + O, and NO + N.³ 

Figure 17 compares the present QSS dissociation rate coefficients in the O₂ + O₂ system with data from literature.^{6,7,11,40,45–48} Since some of excited electronic states of O₂ are closely spaced near the ground level, the multi-surface correction factor⁴⁹ 16/3 is multiplied to the present QSS rate coefficients to consider the dissociation contribution from the excited electronic states. Those results are labeled as w Corr. Fac. in the figure. The present results obtained from using the TC and VS models are in very close agreement to each other, and they are well matched with the DMS data,⁶ although slightly larger discrepancy between the present results and the DMS data is observed in the lower temperature region. It the multi-surface correction factor is considered, the predicted data from this work show better agreement with the measurements^7,45–48 except for the nonequilibrium experiment conducted by Ibraguimova et al.¹¹ This measured data is closer to the present result that does not include the correction factor. The trend of inconsistency implies uncertainty in both measurement and modeling to consider dissociation from the excited electronic states. It is worth mentioning that the Park model⁴⁰ represents somewhat different temperature dependency compared to the present data in the high-temperature region. In Fig. S15, the overall QSS dissociation rates predicted by the TC model shown in Fig. 17 are decomposed into individual components by the two reaction channels in Eqs. (10) and (11).

In Fig. 18, comparison of the average energy loss ratio for the combined oxygen mixture (i.e., O₂ + O and O₂ + O₂) is presented. It is important to note that unlike in Fig. 17, the multi-surface correction factor here is multiplied to the rovibrational-specific dissociation rate coefficients to consider the nonlinear nature in master equations [see Eq. (13)]. In Fig. 18(a), the predicted results are compared with the measured vibrational energy loss ratio reported by Streicher et al.⁷ The present rovibrational-specific result shows closer agreement with the measured data than the VS approach, which overpredicts the measured quantity. Including the multi-surface correction factor improves the agreement between the experimental data and the present rovibrational-specific calculations by lowering the predicted value. This implies that the contribution of the excited electronic states to dissociation is commonly captured by both the experiment and numerical modeling. It is found that the average rotational energy loss ratio is relatively less sensitive to the multi-surface correction, especially at T = 7500 and 10 000 K. This might be attributed to the fact that the dissociation is mostly controlled by the high-v and low-J levels at those temperatures. Then as T increases, the high-J levels at low-v region have larger contribution that ends up the greater sensitivity to the multi-surface correction. It should be noted that the QSS value of the average energy loss ratio of the combined O₂ + O and O₂ + O₂ system is almost identical with that of O₂ + O (i.e., O₂ + O) only. It is because the dissociation of O₂ + O and O₂ + O₂ during the QSS period is mostly governed by the O₂ + O system as observed from Fig. 15. This fact is supported by Fig. 18(b), showing the decomposition of the average vibrational energy loss ratio. During the QSS period, the overall curve in red is very close to that of O₂ + O. From the present investigation, we conclude that the quantity measured by Streicher et al.⁷ most likely corresponds to the O₂ + O collision rather than O₂ + O₂.

This work proposes a reduced-order rovibrational kinetic model of the O₂ + O₂ system for applications relevant to high-enthalpy flows. To overcome the exponential number of possible reactive channels that characterizes diatom–diatom systems, the population of internal levels of one of the reactants is assumed to follow a Boltzmann distribution at a prescribed internal temperature. Detailed QCT calculations are performed to determine the relevant inelastic and reactive reaction rate coefficients, carefully selecting between the endothermic and exothermic directions, the one characterized by the lower statistical error. The rovibrational-specific kinetic database constructed covers a wide temperature range (7500–20 000 K) and is used in the master equation solver to study the kinetics of thermochemical relaxation in isothermal isochoric conditions. To the authors’ best knowledge, this work constitutes the first application of the rovibrational master equation approach to a diatom–diatom system (except for the H₂ + H₂ system, which is far more straightforward). The results obtained have been compared with the DMS calculations, showing excellent agreement in regard to population distribution and average energy and composition.

The main findings are summarized below:

• The rovibrational population distribution of the O₂ molecules is found to depart from equilibrium throughout the dissociation process. The extent of nonequilibrium is more pronounced if compared to the O₂ + O case. Furthermore, the vibrational–rotational mixing is less efficient than in the O₂ + O case, leading to a more distinct separation of the rotational strands characterized by a constant vibrational quantum number. The strands exhibit significant curvature due to less efficient inelastic and exchange energy processes.

• The internal temperature used to prescribe the population of the collision partner has an impact on the results. It leads to a parametrization of the reaction rate parameters as a function of T and T_int.

• Analysis of the combined O₂ + O and O₂ + O₂ system shows how O₂ + O₂ collision governs the energy transfer in the early stage, whereas O₂ + O takes control of the dissociation physics in the QSS period. This is an important finding as it facilitates the construction of reduced-order models for oxygen chemistry.

• Finally, the predicted macroscopic quantities, such as the rovibrational relaxation times, QSS dissociation rate coefficients, and dissociation–vibration coupling constants, agree with the available literature data. This provides validity to the physical model constructed in the present work.

See the supplementary material for additional QCT and master equation results.

This work was supported by AFOSR Grant No. FA9550-22-1-0039 with Dr. Sarah Popkin as Program Officer. The authors would like to thank Dr. R. L. Jaffe (NASA AMES Research Center) for useful discussions on QCT calculations.

The authors have no conflicts to disclose.

Sung Min Jo: Conceptualization (equal); Methodology (equal); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Simone Venturi: Conceptualization (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Jae Gang Kim: Validation (supporting). Marco Panesi: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In Fig. 19(a), the number of sampled bound–bound trajectories are presented at T = 10 000 K and T_int = 300 K. A bound–bound trajectory denotes that a given initial collision pair ends up with another one in which both target and collider are in the bound or quasi-bound states. For each energy level, 30 000 trajectories are simulated. Majority of the trajectories lead to bound final states for the low-lying initial target level (i.e., low level index i). As the initial target state is aligned to higher v and J, the final energy level is more probable to end up in the dissociated state.

Figure 19(b) shows the sampled bound–bound trajectories decomposed into the endothermic and exothermic forward directions. Determination of the energy direction is done using the individual trajectories that are not summed up for the collision partner’s final state. As a result, all four energy levels in Eq. (1) are considered. In the sampled trajectories, the endothermic direction overwhelms the exothermic one, especially in the low-lying energy levels that govern the rovibrational energy transfer. This also happens for the trajectory calculations at T_int = 10 000 K as shown in Fig. S2 in the supplementary material. Fig. 19 shows that the endothermic transitions of the low-lying states are more probable to occur in description of the rovibrational–translational (RV-T) and RV–RV energy transfers when the initial states of collision partner are sampled from the Boltzmann distribution at a given T_int.

In a similar vein, Fig. S3 shows the ratio of the number of sampled trajectories by Eq. (9) to the exact one from Eq. (1) for both exothermic and endothermic directions. It is clear that the endothermic trajectories have better sampling quality than those of the exothermic direction, especially for the heads of each v [i.e., i = i(v, J = 0)] that governs the vibrational energy transfer.

The set of master equations corresponding to the combined O₂ + O and O₂ + O₂ system leads to

where $k̃i→j$⁠, $k̃j→i$⁠, $k̃i→c$⁠, and $k̃c→i$ are the excitation/de-excitation and dissociation/recombination rate coefficients of O₂ + O system.

The average rotational, vibrational, and internal (i.e., rovibrational) energy and temperature are, respectively, denoted as E_R, E_V, E_int and T_R, T_V, T_int. By solving the master equations in Eqs. (13) and (14), the state population n_i of O₂ is obtained at each time t. The corresponding average energy can then be calculated as

where $eri$⁠, $evi$⁠, and $einti$ denote the rotational, vibrational, and rovibrational energy of the ith level, respectively. The temperatures then can be derived as

where E^B for each energy mode denotes the average energy based on a Maxwell–Boltzmann distribution of the internal states defined by the partition function. The temperatures are evaluated by solving Eqs. (C4)–(C6) using a Newton–Raphson method.