A vibrational state-specific model for dissociation and recombination reactions within the direct simulation Monte Carlo method is introduced to study the energy level dynamics of the O_{2} + O system. The state-resolved cross sections for vibrational relaxation and dissociation reactions are obtained from a rotationally averaged quasi-classical trajectory database based on the Varandas and Pais $O2(\u2009\u20093\Sigma g\u2212)+O(\u2009\u20093P)$ potential energy surface. A two-step binary collision framework is outlined to characterize the vibrational state-resolved recombination probabilities, which are constrained by detailed balance for orbiting pair formation, and microscopic reversibility applied to the dissociation cross sections for orbiting pair stabilization. The vibrational state-to-state (STS) model is compared to the phenomenological total collision energy (TCE) and quantum kinetic (QK) models through a series of 0-d non-equilibrium relaxation calculations. A quasi-steady state (QSS) region is established in the vibrational temperature profiles of the TCE, QK, and STS models under non-equilibrium heating. This QSS region is a result of the competition between vibrational relaxation by vibrational-translational (VT) transitions and O_{2} dissociation. The duration of QSS predicted by the STS model is approximately ten and four times that of the TCE and QK model predictions, respectively, and the total time to reach equilibrium is approximately 3.5 times that of the TCE model and 1.5 times that of the QK model. A distinct QSS region is not observed in the non-equilibrium cooling case. This is attributed to the relatively rapid VT transitions that work to equilibrate the vibrational energy distribution upon recombination, which is comparatively slow. The total time to reach equilibrium by the STS model in the non-equilibrium cooling case is five times and three times greater than those of the QK and TCE models, respectively.

## I. INTRODUCTION

Modeling the flowfield surrounding hypersonic vehicles poses a significant challenge owing to the complex thermochemical processes in the shock-heated region and in the wake. These regions are characterized by both thermal and chemical non-equilibrium owing to the disparity in characteristic relaxation times. Models based on state-to-state (STS) kinetics have received considerable attention and have been employed to study relaxation behind strong shock waves^{1,2} in nozzle flows^{3,4} and in boundary layer flows.^{5} In particular, the $N2(\u2009\u20091\Sigma g+)+N2(\u2009\u20091\Sigma g+)$ system,^{6–10} the $N2(\u2009\u20091\Sigma g+)+N(\u2009\u20094Su)$ system,^{11–14} and the $O2(X3\Sigma g\u2212)+O(\u2009\u20093P)$ system^{15–17} have been studied extensively owing to their importance in Earth atmospheric reentry. These efforts have demonstrated the inadequacy of legacy multi-temperature models, but most importantly, these studies have provided a wealth of information regarding the complex evolution of the internal energy distributions under thermal and chemical non-equilibrium. A recent study by Macdonald *et al.*^{8} examined the non-equilibrium processes during dissociation of the $N2(\u2009\u20091\Sigma g+)+N2(\u2009\u20091\Sigma g+)$ system. It was established that dissociation processes are driven by the high-energy tail of the rovibrational distribution. Furthermore, nearly all dissociation occurs under quasi-steady state (QSS) conditions at high temperatures, while approximately half of the dissociation occurs after QSS at low (<13 000 K) temperatures. This is in contrast to dissociation processes for the $N2(\u2009\u20091\Sigma g+)+N(\u2009\u20094Su)$ system, where more than half of dissociation occurs before QSS conditions are established.^{13} Panesi *et al.*^{13} also found that rotational-translational (RT) and vibrational-translational (VT) relaxation times converge at high temperatures in nitrogen, again contrary to the conventional idea that vibrational relaxation occurs over comparatively longer characteristic time scales.

Recent work by Andrienko and Boyd^{15} outlines a comprehensive master equation study of the O_{2} + O system based on the $O2(\u2009\u20093\Sigma g\u2212)+O(\u2009\u20093P)$ potential energy surface by Varandas and Pais.^{18} The STS rovibrationally resolved cross sections for inelastic (non-reactive), exchange, and dissociation collisional processes employed in the master equation calculation were obtained from quasi-classical trajectory (QCT) calculations.^{15} The cross sections were computed on a single surface corresponding to the ground electronic state (accounting for 1/27 degeneracy). VT processes resulting from this surface were thus assumed to represent processes occurring on the other eight surfaces, while a constant scaling factor of 16/3 proposed by Nikitin was used to account for dissociation from electronic excited states.^{19} A number of important results emerged from the analysis of this system. Exchange reactions were found to play a significant role in the thermalization of the internal energies during rovibrational heating, particularly for highly excited quasi-bound states. Furthermore, equilibration of rotational and vibrational degrees of freedom was found to occur over similar time scales for a broad range of temperatures, which is in contrast to recent studies of N_{2}–N and O_{2}–Ar molecular systems.^{13,20} Thermochemical relaxation under dissociation, recombination, and exchange reactions was also studied in a 0-d heat bath calculation for a range of translational temperatures. A QSS region is established owing to the competition between internal energy relaxation (via RVT non-reactive transitions and exchange reactions) and dissociation. Remarkably, however, the characteristic QSS rotational and vibrational temperatures exhibit a maximum at 10 000 K and decrease for higher heat bath temperatures. This is in stark contrast to both the N_{2}–N and N_{2}–N_{2} systems, which show a monotonic increase in the internal temperatures under QSS.

State-specific models have also been developed for use in the direct simulation Monte Carlo (DSMC) method for $N2(\u2009\u20091\Sigma g+)+N(\u2009\u20094Su)$, $N2(\u2009\u20091\Sigma g+)+N2(\u2009\u20091\Sigma g+)$, and $O2(\u2009\u20093\Sigma g\u2212)+O(\u2009\u20093P)$^{21–25} in place of the widely used phenomenological models, such as the Larsen-Borgnakke (LB) model for internal energy relaxation^{26} and the total-collision-energy (TCE) and quantum-kinetic (QK) models for chemical reactions.^{27–29} The state-specific cross sections for internal energy exchange favor monoquantum transitions, in contrast to the Larsen-Borgnakke model in which transitions to the ground state are most favored.^{24,25} A recent study by Borges *et al.*^{23} compared the state-specific dissociation (SSD) rate coefficients obtained from $O2(\u2009\u20093\Sigma g\u2212)+O(\u2009\u20093P)$ QCT calculations with rate coefficients from QCT and the TCE model in which the TCE model parameters were calibrated with the QCT data.^{23} It was found that in non-equilibrium conditions, the TCE model overestimates the QCT dissociation rate coefficient by more than an order of magnitude in vibrationally cold conditions and underestimates the dissociation rate coefficient in vibrationally hot conditions.

Although it is necessary to include recombination reactions to capture the correct mixture composition and temperature in the shock layer and in the wake of hypersonic flows,^{27,30} a formal treatment of the recombination reaction in state-resolved hypersonic DSMC calculations has been lacking. Recombination modeling in DSMC is included in both TCE and QK models,^{31} in which the probability of recombination is computed from the dissociation rate coefficient and an equilibrium constant. The internal state of the recombining molecule is then selected based on an assumed Boltzmann distribution at a given temperature, which cannot capture the up-pumping effect observed in detailed state-resolved calculations.^{4,32,33} Recently, however, Gimelshein and Wysong proposed a new chemistry model in DSMC, inclusive of both dissociation and recombination reactions, that utilizes the vibrational bias model of dissociation and an empirical microscopic recombination model.^{34–36} The dissociation probability, *P*_{d}, is computed for each dissociation collision based on the translational and internal energies of the colliding pair and includes two parameters, *λ* and *A*_{d}. *A*_{d} is calibrated to match a known rate constant of the specific reaction, and *λ* is used to specify the degree of vibrational favoring. The recombination probability, *P*_{r}, is computed from the equilibrium recombination rate coefficient *k*_{r} of Arrhenius form for a rovibrational state (*v*, *j*), which is obtained as the ratio of the corresponding dissociation rate coefficient *k*_{d} to the equilibrium constant (further details may be found in Ref. 34). An important feature of this model is that by its construction, both the principle of detailed balance and microscopic reversibility are satisfied inherently. Detailed balance is achieved through the determination of *k*_{r} as mentioned above, while microscopic reversibility is satisfied indirectly through the selection of the internal energy states of recombining pairs; the internal energy state of dissociated molecules is stored in a list in each computational cell, and the internal state of each recombining molecule is selected from this list. *Thus, the inelastic recombination (backward) cross section for each rovibrational state (v, j) is inherently reconstructed through the sampling of this list, which is populated according to the dissociation (forward) processes*. A similar approach was also employed for the TCE model.^{1} After assigning the internal state to the recombining molecule, the relative translational energy between the recombined product and the third body is determined based on energy conservation. While the chemistry model by Gimelshein and Wysong addresses both dissociation and recombination reactions in an elegant way that satisfies detailed balance and microscopic reversibility, the computed dissociation (and, therefore, recombination) probability relies on the calibration of *A*_{d} and *λ* to match the known equilibrium or non-equilibrium reaction rate constants.

This work introduces a new vibrational state-to-state DSMC chemistry model that has been developed for dissociation/recombination and internal energy relaxation of the O_{2} + O system. This model is based on vibrational state-specific VT and dissociation cross sections obtained from QCT calculations^{16} involving the Varandas and Pais potential energy surface (PES),^{18} which has been previously used for studying the O_{2} + O system.^{15,23,37,38} The database utilized in this work was developed by Esposito and Capitelli which includes rotationally averaged cross sections.^{38} As will be detailed later in Sec. III, this database is a suitable first step for the extension of the two-step binary collision (TSBC) model to state-resolved cross sections,^{35} which is the aim of the present work. It is noted that a rovibrationally resolved cross-sectional database has recently been developed by Andrienko and Boyd^{15} based on the Varandas and Pais PES and new accurate PES by Varga *et al.*^{39} and Paukku *et al.*^{40,41} for the O_{3} and O_{4} systems is also recently available. This paper aims to first assess the vibrational state-resolved DSMC chemistry model based on the O_{3} database from Esposito and Capitelli. This will enable important comparisons to the fully resolved rovibrational state-to-state model based on the O_{3} database of Andrienko and Boyd and to the new O_{3} PES by Varga *et al.* Recent comparisons of the Varandas and Pais PES to the 1^{1}*A*′ surface by Varga show significant differences in the structure of the potential as well as the vibrational relaxation rates. A full discussion of these comparisons may be found in Grover *et al.*^{42}

The present model employs the two-step binary collision (TSBC) approach for DSMC recombination, which was first introduced by Koura^{35} for a modified forced harmonic oscillator. The TSBC approach treats recombination through two binary interaction processes, including orbiting pair formation (constrained by detailed balance) followed by subsequent stabilization involving a third body collision. As will be presented later, the internal energy of the recombined molecule formed through the stabilization collision is determined through the recombination cross section via microscopic reversibility, but in such a way that this principle is satisfied directly from the known vibrational state-resolved QCT dissociation cross sections. The internal energy level dynamics of the O_{2} + O system are investigated for 0-d vibrational heating and cooling cases in the presence of both dissociation and recombination reactions, and these results are compared against existing phenomenological models. To the authors’ best knowledge, this is the first time that both dissociation and recombination reactions have been included in the study of energy level dynamics of the O_{2} + O system in DSMC. The results from the relaxation simulations capture a distinct QSS region in the vibrational heating cases examined, but a QSS region is not established during vibrational cooling.

The rest of the paper is organized as follows: a brief overview of the Larsen-Borgnakke, TCE and QK phenomenological reaction models is presented in Sec. II, followed by a detailed outline of the TSBC approach employed in the current chemistry model in Sec. III. In Sec. IV A, a comparison of the dissociation and recombination rate coefficients computed from QCT with those of the TCE, QK, and STS models in equilibrium (*T* = *T*_{v}) and non-equilibrium (*T* ≠ *T*_{v}) conditions is presented. Results of the 0-d vibrational heating and cooling cases are presented in Sec. IV B, and the observed differences in the QSS regions predicted by each of the models are discussed, followed by conclusions.

## II. PHENOMENOLOGICAL CHEMISTRY MODELS IN DSMC

In this section, a brief overview of the phenomenological DSMC chemistry models for internal energy exchange, dissociation, and recombination reactions is presented. Although many chemistry models have been developed for DSMC, we limit our discussion to the Larsen-Borgnakke (LB), TCE, and QK models for later comparison to the STS chemistry model.

### A. Larsen-Borgnakke model

The Larsen-Borgnakke model has been widely used for describing energy transfer in non-reactive inelastic collisions^{26} owing to its simplicity and its ability to reproduce relaxation with the Landau-Teller form. A comparison of the LB cross sections with the QCT cross sections for VT transitions, however, shows significant differences, as seen in Fig. 1. The LB cross sections were evaluated using the variable hard sphere (VHS) model and a vibrational collision number based on the Millikan-White form with Park’s high-temperature correction^{43} and the Gimelshein DSMC collision number correction factor.^{44} The graphs in the figure show a comparison of the LB and QCT cross sections for VT transitions (including inelastic and exchange collisions; excluding elastic collision) of the O_{2} + O system from initial vibrational states *v* = 0, *v* = 20, and *v* = 40 to all possible final states. A rotational temperature of *T*_{rot} = 300 K and a relative collision energy of *E*_{tr} = 2 eV are specified. It can be seen that the LB cross section (and hence the probability for transition into the final vibrational state, *v*′, indicated by the x-axis) is the largest for transitions into the lowest state, regardless of the initial vibrational state. In contrast, QCT cross sections favor mono-quantum transitions over multi-quantum transitions. Thus, the most probable post-collision vibrational state is in general not the ground state and depends upon the initial vibrational state. This QCT cross-sectional behavior is well-established in the literature and has been observed by others for the O_{2} + O PES employed in this work.^{17}

### B. Quantum kinetic model

One of the main advantages of the QK model is that it does not require information about the rates of reaction.^{45,46} Starting from the LB model, the post-collision vibrational states are chosen from a Boltzmann distribution corresponding to a collision energy equivalent temperature. The collision energy is the sum of the relative translational energy and the vibrational energy. For instance, in a collision of O_{2} and O, the collision energy used to redistribute the vibrational energy of O_{2} is given by *E*_{coll} = *E*_{tr} + *E*_{v}, where *E*_{v} is the pre-collision vibrational energy of O_{2}. The post-collision vibrational state *v*′ is chosen using a uniform random number from 0 to $vmax\u2032$ = ⌊*E*_{coll}/(*k*_{B}*θ*)⌋ and accepted with a probability of

In the phenomenological models considered here, a simple harmonic oscillator (SHO) assumption is used, meaning the vibrational levels are spaced evenly by a characteristic vibrational temperature *θ*_{vib}. A dissociation event occurs if the maximum possible vibrational state $vmax\u2032$ satisfies $vmax\u2032$ > *θ*_{d}/*θ*_{vib}, where *θ*_{d} is the characteristic dissociation temperature. The QK dissociation rate can be expressed explicitly for VHS gas at equilibrium temperature *T*,

where *σ*_{ref} is the cross section at temperature *T*_{ref} and *ϵ* is the symmetry factor. Recombination reactions occur when there is a third body particle within a ternary collision volume *V*_{coll}. The corresponding recombination probability can be expressed as

where *n*_{T} is the number density of the third body particles and *V*_{coll} is suggested by Bird^{28} in the form $Vcoll=a(T/\theta vib)bVref$, where the reference volume *V*_{ref} is conveniently set to the volume of the sphere with the radius equal to the sum of the reference radii of the three atoms.^{28} For an equilibrium VHS gas, the recombination rate coefficient can be expressed as

### C. Total collision energy model

The TCE model^{27} expresses the probability of dissociation in terms of an Arrhenius dissociation rate constant of the form

with Arrhenius parameters Λ and *η*. *E*_{d} denotes the dissociation energy. Unlike the QK model, the collision energy of the TCE model takes the rotational energy into account, which is given as a sum of the internal and relative translational energies; i.e., for the O_{2} + O collision, *E*_{coll} = *E*_{tr} + *E*_{v} + *E*_{j}, where *E*_{v} and *E*_{j} are the vibrational and the rotational energy of pre-collision O_{2}, respectively. The probability of dissociation is given by

where $\zeta \xaf$ is the average number of internal degrees of freedom contributing to the collision energy. Again, an equilibrium constant *K*_{eq} can be determined to cast the recombination rate in terms of the dissociation rate. Since it is assumed that the recombination reactions have zero activation energy, the recombination rate can be expressed as

The recombination probability which yields the rate in Eq. (7) is

Similar to the QK model, the TCE model was recently improved to satisfy detailed balance.^{1}

## III. VIBRATIONAL STATE-TO-STATE (STS) DSMC CHEMISTRY MODEL

To provide a more detailed description of dissociation and recombination in DSMC, the phenomenological QK and TCE models are replaced by a high-fidelity vibrationally resolved state-to-state (STS) description using QCT collision cross sections. The cross sections are computed by performing scattering simulations of individual atoms or molecules for a prescribed pre-collision state to determine the resulting post-collision state. This provides a set of state-resolved dissociation cross sections for direct use in DSMC, while the recombination cross sections are inferred according to the TSBC model. In this section, we first introduce the STS cross-sectional database employed in this work. We then provide an overview of the TSBC model for the O_{3} system and assess the appropriateness of an indirect (two step) recombination mechanism for O + O interactions, including O + O orbiting pairs and their associated lifetime. We then develop in detail the mathematical descriptions of dissociation and recombination kinetics for the TSBC model, followed by their implementation in DSMC.

This work employs the available cross-sectional database from Esposito and Capitelli^{37,38} for the $O2(\u2009\u20093\Sigma g\u2212)+O(\u2009\u20093P)$ system.^{18} The potential energy surface was constructed using the double many-body expansion (DMBE) method for the ground electronic state of ozone. Rovibrational levels were computed using the WKB approximation for the O_{2} potential.^{18} The cross sections are available at pre-collision vibrational states: *v*_{init} = 0, 1, 3, 5, 7, 10, 15, 20, 25, 30, 35, 39, 40–46; rotational temperatures: *T*_{rot} = 300, 3000, 20 000 K; and relative translational energies *E*_{tr} = 0.006–9 eV in increments of 0.01 eV. Details of this database may be found in Refs. 16, 18, 37, and 38.

It should be noted that the reactions in these QCT calculations proceed on a single surface (accounting for 1/27 of the degenerate states) corresponding to the ground electronic state. VT processes involving the other eight surfaces, representing the remaining degeneracies, were not considered in the generation of this database. Dissociation processes, however, are unlikely to proceed adiabatically because the required dissociation energy (>5 eV) is much higher than the ^{1}Δ first excited molecular state, which is about 1 eV above the molecular ground state. Accordingly, Esposito and Capitelli introduced a variable multisurface factor^{38} in the current database to account for dissociation from electronically excited states. This variable multisurface factor is an extension of the constant 16/3 factor introduced by Nikitin,^{19} in which the dissociation cross section at an electronically excited state is estimated as its degeneracy divided by the degeneracy of the electronic ground state. The variable multisurface factor employed in this work is computed as the summation of the degeneracy of each excited electronic state which lies at or below the energy of the rovibrational level under consideration, divided by the degeneracy of the ground state. Further details of this variable multisurface factor and of the QCT calculations are provided in Ref. 38. In summary, there is only one PES considered in the current work, which is the representative of 1/27 of O_{2} + O collisions, and the variable multisurface factor proposed by Esposito and Capitelli is employed to account for dissociation from electronically excited states.

As mentioned earlier, the database employed in this work provides rotationally averaged cross sections. The QCT trajectories are carried out for a specified relative translational energy, *E*_{tr}, and an initial vibrational state, *v*. The initial rotational state, *j*, is drawn from a Boltzmann distribution at temperature *T*_{rot}, based on the WKB approximation of rovibrational levels for the O_{2} potential. The resulting cross sections are obtained by summing over all final rotational states, *j*, corresponding to a given final vibrational state, *v* (full details may be found in Ref. 38),

In this model, the influence of rotational energy on the reaction probability is limited to cell-averaged *T*_{rot}. Because of this, only *E*_{tr} and *v* have a collision-specific influence on the collision probability, and the reaction of interest is expressed in terms of *v* only,

Since the cross sections contained in this QCT database have significant scatter, a polynomial with the form below is used to fit the cross sections as a function of relative translational energy,

A bilinear interpolation technique is applied to determine cross sections for vibrational states (or rotational temperatures) not explicitly available in the database. Although the cross sections employed in this work are rotationally averaged, the study by Andrienko and Boyd suggests that O_{2} rotational nonequilibrium is an important factor at low temperatures, where recombination is most significant. Further studies should investigate the role of rotational nonequilibrium in recombination using a set of rotationally resolved cross sections.

Although the current database does not provide rotationally resolved cross sections, it is well suited for Koura’s two-step binary collision (TSBC) model,^{35} which considers only vibrational non-equilibrium. Our objective is to thus extend Koura’s approach for the current QCT database. The TSBC model describes dissociation and recombination reactions in the following two steps:

The first step of the dissociation reaction [the forward direction of Eq. (12)] is the collision of O_{2} with a third-body oxygen atom, O^{(3)}, which forms an orbiting pair, OP, in the vibrational free state *f*. The OP is formed from the atoms of the O_{2} molecule and is assumed as a unique state in the chemistry model. The second step of dissociation [the forward direction of Eq. (13)] is the splitting of the OP into its two oxygen atoms, referred to as O^{(1)} and O^{(2)}.

Recombination reactions are modeled as the reverse of the dissociation reaction detailed above. The first step of recombination [the backward direction of Eq. (13)] is the formation of an OP from oxygen atoms O^{(1)} and O^{(2)}, and the third body O^{(3)} remains unchanged. The second step [the backward direction of Eq. (12)] is the stabilization of the OP through collision with the third-body oxygen atom O^{(3)}. This produces the recombined O_{2}, while the third-body O^{(3)} remains unchanged. A diagram outlining the TSBC process is shown in Fig. 2.

Before discussing the details of the TSBC model as applied to the QCT database, we first assess the appropriateness of an indirect (two step) mechanism for the recombination of O_{2}. We consider specifically the evidence for O + O orbiting and the associated lifetime of these orbiting pairs. It is well-known that orbiting occurs if there exists an attractive portion of the interaction potential,^{48–50} and such orbiting has been observed by Weaver and Alexeenko^{51} and Subramaniam and Stephani^{52} for the O_{3} system based on the Varandas and Pais PES considered in this work. A recent study by Elchinger and co-workers^{53} also demonstrates long-lived O + O orbiting interactions based on the Hulbert-Hirschfelder pairwise potential. Although this potential is different from the pairwise potential employed in the Varandas and Pais PES, both feature an attractive portion of the interaction potential, and thus, long-lived orbiting interactions are observed.

The time spent in orbit prolongs the interaction between the O atoms, making stabilization of the orbiting pair (i.e., O_{2} recombination) by a third body more likely as the orbiting time increases. Obviously, as the orbiting time increases, the range of impact parameters (i.e., the cross section) that ultimately results in appreciable orbiting decreases. If appreciable orbiting is defined as greater than three revolutions, then the corresponding orbiting pair cross section is approximately *σ*_{OP} = *O*(10^{−3})Å^{2}. In what follows, by employing Eq. (29) with *E*_{tr} = 0.1 eV, the corresponding orbiting lifetime is approximately *τ*_{OP} = *O*(10^{−12})s. The orbiting time is thus very short, which is important for the development of the TSBC model below.

In order to compute the recombination probability in DSMC, several assumptions are introduced in the development of the TSBC model: (1) the orbiting pair is assumed as a unique state, and the mechanisms in Eqs. (12) and (13) are assumed as the only possible recombination channels. Therefore, only the two-step (indirect) recombination pathway is considered; other possible channels, such as the direct three-body recombination mechanism, are not considered. (2) It is assumed that the OP splitting and its formation [forward and backward reactions of Eq. (13)] are sufficiently rapid such that the formation and splitting processes may be assumed to be in equilibrium. This assumption implies that the excitation cross section from the vibrational state *v* to the OP free state may be approximated by the dissociation cross section. (3) The velocity distribution function of oxygen atoms follows the Maxwellian distribution, which allows us to derive the analytical form of expressions used in the TSBC model. (4) The orbiting lifetime is assumed to be constant in this work as suggested by Koura.^{35} As will be discussed below, the value of orbiting lifetime itself is never evaluated in the current TSBC model; rather, only the product of the orbiting pair cross section and its lifetime is evaluated.

### A. Dissociation reactions

In Koura’s TSBC model, the dissociation reaction can be thought of as the following reaction:

Note that the reaction in Eq. (14) is simply the result of the two consecutive reactions in Eqs. (12) and (13). Energy conservation is applied to dissociation by

where *E*_{tr} is the relative translational energy of O_{2} with O^{(3)} [Eq. (12)]. Although the rotational energy does not directly influence the dissociation probability in this model, *E*_{rot} goes into the post-reaction scattering of O^{(1)} and O^{(2)} (i.e., the splitting of the OP) through energy and momentum conservation.

The description of dissociation through Eqs. (12) and (13) is necessary in order to establish microscopic reversibility for the reverse recombination process described in Sec. III B. The approach outlined above, in which the rotational energy has no influence on the dissociation probability but is used to split the dissociating molecule, is identical to the treatment of the QK model and similar to the TCE model. Elastic collisions are not included in the current QCT database and are accounted for using the phenomenological VHS model. Finally, the probability of a DSMC dissociation for each collision is simply the ratio of the dissociation cross section to the total cross section,

The total cross section *σ*_{tot} is the sum of all possible interaction cross sections for O + O_{2} collisions, including elastic, inelastic (non-reactive), exchange, and dissociation reactions.

Finally, the state-specific dissociation rate coefficient can be computed by integrating the state-specific dissociation cross section over the velocity space,

where *f*_{M} and *f*_{A} are the velocity distribution functions of the colliding molecule and atom, respectively, and are assumed to be Maxwellian for a given temperature *T*. It should be noted that Eq. (17) is only computed in post-processing for the purpose of analysis; the dissociation rate coefficient is not required as input in DSMC simulations.

### B. Recombination reactions

Because the cross sections for the reverse recombination reaction 3O → O_{2} + O are not included in the QCT database, these cross sections must be backed out from the forward dissociation reaction through detailed balance and microscopic reversibility. Again following the TSBC method by Koura,^{35} recombination reactions occur through two steps: (i) OP formation [the backward direction of Eq. (13)] and (ii) OP stabilization [the backward direction of Eq. (12)], as shown in Fig. 2. The stabilization cross sections are determined from the state-specific dissociation cross sections through microscopic reversibility, while the orbiting pair cross sections are derived from detailed balance using equilibrium rates. The derivations of OP formation and stabilization cross sections are introduced next, followed by the state-specific recombination probability for DSMC.

The forward process of OP splitting and the backward process of OP formation are assumed to be rapid and in equilibrium,^{35} where

Therefore, the sum of the stabilization rate over all vibrational states can be referred to as the overall recombination rate according to the above assumption, which gives a recombination rate

Substituting for *n*_{OP} from Eq. (18) into Eq. (19) gives the recombination rate coefficient at a given equilibrium temperature *T*,

The lifetime of the orbiting pair (*τ*_{OP}) is introduced in Eq. (20) where $\tau OP=ksplit\u22121$, which is assumed to be constant.

In order to determine the recombination probability, the recombination rate coefficient must be expressed in terms of the dissociation rate coefficient. The equilibrium dissociation rate coefficient is expressed by the sum of the vibrational state-specific dissociation rate coefficient weighted by the Boltzmann vibrational distribution at temperature *T*,

Recalling the second assumption outlined in the beginning of this section, the rate coefficient from state *v* → *f* is assumed to be approximately equal to the dissociation rate coefficient ($kdv\u2248kvf$). The law of mass action can then be used to relate the dissociation and recombination rate constants as^{54}

where *Q*_{M} and *Q*_{A} are partition functions of molecular and atomic oxygen, respectively. The superscripts *tr*, *r*, *v*, and *el* denote the translational, rotational, vibrational, and electronic mode, respectively. Combining Eqs. (22), (20), and (21) and employing detailed balance between *k*_{vf}(*T*) and *k*_{fv}(*T*),

we obtain

It is important to note that since the forward and backward reactions of Eq. (13) are assumed to be in equilibrium, the TSBC model does not capture the pressure dependence for the recombination rate coefficient and displays the standard fall-off region with respect to pressure. Such a pressure dependence can however be captured if the steady-state approximation is used, which will be addressed in future work. It is worth noting that the recombination rate coefficient in the current work is consistent with the low-pressure limit, for which the equilibrium and steady-state approximations yield the same rate coefficient. To compute the probability of recombination in DSMC, consider the expression for *k*_{OP} from kinetic theory,

where the OP formation rate coefficient is obtained from the OP formation cross section when integrated over the velocity distribution functions for the reactant atoms. Assuming that the oxygen atoms follow a Maxwellian distribution at the temperature T, we have the following:

where *μ*^{O–O} is the reduced mass of two oxygen atoms. Here, we introduce the average value of the quantity *E*_{tr}*σ*_{OP}, which can be expressed as

Combining Eqs. (24) and (26) with the ratio of the electronic partition function $\eta el=(QAel)2/QMel$, the translational partition function for species s, $Qstr=(2\pi kTms/h2)3/2$, and the rotational partition function $QMr=T/2\theta r$ and applying Eq. (27), we obtain

where *m*_{O} is the molecular weight of an oxygen atom and *h* is Planck’s constant. To implement in DSMC, the product *σ*_{OP}*τ*_{OP} for each O + O collision pair is computed based on its relative collision energy *E*_{tr},

The total cross section of OP stabilization [the backward direction of Eq. (12)], *σ*_{stab}, is the sum of the cross section of all possible *f* → *v* transitions,

The cross section of stabilization into each vibrational state *v* (*σ*_{fv}) can be determined in terms of the corresponding excitation cross section *σ*_{vf} (from QCT database) via microscopic reversibility and energy conservation,^{55}

where $EtrOP\u2212O$ and $EtrO2(v)\u2212O$ are the relative translational energy between OP and O and between O_{2}(*v*) and O, respectively. After computing $EtrO2(v)\u2212O$ by Eq. (32) for the specific vibrational state *v* and substituting this value into Eq. (31), *σ*_{fv} can then be solved. Note that the rotational energy of the stabilized O_{2} is assumed equal to the orbiting energy (the relative translational energy between two atoms with respect to the center-of-mass frame), consistent with the dissociation process,^{35} while the bond energy goes into vibration.

### C. DSMC algorithm

Previous studies have proposed various approaches to three-body collisions in DSMC,^{56} including NTC-based methods.^{57} The collision frequencies outlined in Koura’s TSBC method use a time counter approach, which may include null collisions (candidates which are not selected for collision).^{35} The no-time-counter (NTC) scheme for the TSBC model is introduced instead for use in this model. For a binary collision between species *p* and *q*, the number of collision candidates selected is

where *g* denotes the relative velocity of the collision pair, *F* is the number of real molecule represented by a simulation particle, and *ϵ* is the symmetric factor [(1) if *p* and *q* are the same species; (2) otherwise]. The candidates are selected for collision with probability

To extend the NTC scheme for a three-body collision based on the TSBC approach, the number of collision candidates of OP + O is computed first, which is expressed as

with N_{OP} defined^{35} as

and *n*_{O} = *N*_{O}*F*, we substitute Eq. (36) into Eq. (35) to obtain the number of recombination collision candidates,

This form is similar to the NTC scheme for binary collisions but involves a third *N* term to account for the three oxygen atoms involved in the recombination reaction. These recombination collision candidates are then selected for recombination with probability

The vibrational and rotational energies of the recombined molecule are determined according to

The relative translational energy of the recombined species and the third body is then determined through energy conservation.

The no-time-counter scheme of the TSBC algorithm is shown schematically in Fig. 3, and the steps for implementation are as follows:

Three particles (atoms) are randomly selected from the DSMC particle list (first two particles are OP candidates; third particle is a stabilization candidate).

Compute

*σ*_{OP}*τ*_{OP}[Eq. (29)] from the first two particles.Compute

*σ*_{fv}for each possible recombined vibrational state*v*[Eqs. (31) and (32)]. For example, for*v*∈ [0–46], take*v*= 0, solve Eq. (32) for $EtrO2(v)\u2212O$ (all other quantities known), and then solve Eq. (31) for*σ*_{fv}(all other quantities known); repeat for other*v*.Compute

*σ*_{stab}[Eq. (30)].Compute

*P*_{rec}[Eq. (38)]; use acceptance-rejection to determine if recombination occurs.Select vibrational state of the recombined molecule [Eq. (39) via acceptance-rejection].

Assign rotational energy of the recombined molecule equal to orbiting energy [Eq. (40)].

Assign relative translational energy of the recombined molecule and third body through energy conservation.

## IV. RESULTS

The results presented in this section are focused on the relaxation of the state-specific vibrational mode through both thermal and chemical relaxation processes. Although the rotationally averaged cross sections are obtained by assuming that the pre-collision rotational energy is Boltzmann distributed at a temperature *T*_{rot}, the rotational mode is not constrained to be in equilibrium with either translational or vibrational modes, both in the cross-sectional calculation and in the DSMC simulations.

### A. Comparison of equilibrium and non-equilibrium rate coefficients

The TCE, QK, and STS models are first tested in equilibrium conditions, where *T* = *T*_{v} = 4000–20 000 K. In Fig. 4, the QCT rate coefficients are computed by integrating the QCT cross-sectional data directly. All of the TCE parameters for dissociation and recombination reactions are calibrated from the QCT data for a consistent comparison. The experimental data from Ibraguimova^{58} are also plotted in Fig. 4, in which the low-temperature (2000–6000 K) and high-temperature (6000–11 000 K) ranges of the dissociation rate coefficient are provided for comparison in Fig. 4(a). We note that the rates reported by Ibraguimova *et al.* include contributions from O_{2}–O_{2} collisions, and contributions from O_{2}–O interactions are considered negligible. The corresponding recombination coefficients [Fig. 4(b)] are computed based on the experimental dissociation rate coefficients and the equilibrium constant via detailed balance. The dashed line of the experimental data denotes the extrapolation from the valid temperature region.

The agreement between the STS implementation in DSMC and the QCT data is excellent. The phenomenological TCE model also shows good agreement with the QCT results in the equilibrium conditions since the TCE parameters are calibrated from the QCT data. The QK model shows reasonable agreement with the dissociation coefficients from the QCT calculations, but comparison with the recombination coefficients is relatively poor, particularly at high temperatures [Fig. 4(b)]. Recall that the QK parameters are not explicitly calibrated to match the reaction rate coefficients; rather, the parameters *a* and *b* are calibrated to match the equilibrium constant (as discussed earlier in Sec. II). Therefore, it is not surprising to have some discrepancy in the equilibrium dissociation/recombination rate coefficients between QK and STS models.

In comparison with the experimental data, the dissociation rate coefficients of the STS model show good agreement for *T* < 7000 K. The STS recombination rate coefficients show reasonable agreement with the low-temperatureexperimental data. At high temperatures (*T* > 7000 K), the STS model shows a higher dissociation rate coefficient than the experiment. This is also observed by Kulakhmetov *et al.*^{17} Similar discrepancies are observed for the recombination rate coefficient at high temperatures. This overprediction may be a result of the use of rotationally averaged cross sections from the Varandas and Pais PES as well as the nonadiabatic correction factor that is used to account for dissociation through electronically excited reaction pathways. It is also noted that the experiments do not directly measure dissociation rates. Rather, the rates are reconstructed from the absorption data for O_{2}–O_{2} collisions, using gas dynamic relations and determination of a vibrational temperature from the low-lying vibrational states using the Landau-Teller (harmonic oscillator) model.^{58} To obtain the final dissociation rate coefficient of O_{2}–O, a factor is employed based on the O_{2}–O_{2} rate coefficient, $kdO2\u2212O=3.5kdO2\u2212O2$.^{58}

Figure 5 shows the vibrational state-specific dissociation and recombination coefficients at *T* = 7000 K for the STS model compared with the TCE and QK models. This comparison provides a more detailed picture of how the total rate coefficients (shown previously in Fig. 4) compare by examining the contribution to the equilibrium total rate coefficient by each vibrational state. At low vibrational states, the STS model is characterized by a relatively low dissociation and recombination rate coefficient compared to the TCE and QK models, while at higher vibrational states, the STS rate coefficients exceed those of the TCE/QK models. This indicates a preference for dissociation/recombination from high-lying internal (vibrational) states.

The rate coefficients of the TCE, QK, and STS models are compared under non-equilibrium conditions in Fig. 6. The coefficients are computed for *T* = 10 000 K, with the vibrational temperature ranging from *T*_{v} = 4000–20 000 K. Results from the QCT data are also shown for comparison. It is observed in general that the rate coefficients of the TCE and QK phenomenological models have weak vibrational temperature dependence for the dissociation rate coefficient, and these models perform poorly in strong non-equilibrium conditions. The recombination rate coefficients show no vibrational dependence and are approximately constant across the vibrational temperature range based on the current model.

### B. 0-d non-equilibrium relaxation calculations

In this section, the dynamics of the vibrational distribution function and the state-resolved reaction rates arestudied through a series of 0-d non-equilibrium heat bath calculations. The simulations track the chemical relaxation of the gas in an isochoric reactor, in which the translational temperature is held constant. The vibrational temperature is initialized at the beginning of the simulation below or above the translational temperature such that either vibrational heating or cooling occurs as a result of the relaxation process. Dissociation, recombination, and VT transitions are modeled in the STS DSMC simulations using the QCT cross-sectional data and the TSBC recombination model, while RT transitions and elastic collisions are modeled using the LB and VHS models, respectively. It should be noted that the vibrational collision number, *Z*_{v}, which controls the rate of vibrational relaxation in the TCE and QK model implementations was fit to a constant value to match the relaxation rate of the state-based model for a consistent comparison. This step was taken in efforts to minimize the influence of differences in non-reactive inelastic collisions when interpreting the current results. The chemical reaction parameters for the TCE and QK models are shown in Table I. The parameters for the TCE model are computed by fitting the reaction rate coefficients from the QCT database with the Arrhenius form, and the correction suggested by Gimelshein and co-workers^{59} is applied for employing the discrete vibrational energy in DSMC. The parameters *a* and *b* in the QK model^{28} are used to compute the collision volume for the third-body particle using the same equilibrium constant as the STS model.

#### 1. Vibrational heating case (dissociation regime)

The first case presented is the 0-d vibrational heating case, in which the system is held at a constant translational temperature, *T* = 7000 K, and the initial vibrational and rotational temperatures are set to *T*_{v} = *T*_{rot} = 2000 K. The molecular/atomic oxygen gas mixture is initialized with number density 6.79 × 10^{24} and 1.48 × 10^{24}, respectively, such that the system undergoes net dissociation during relaxation toward chemical equilibrium. Approximately 40 000 simulation particles were employed to obtain good statistics.

The results of the simulation are shown in a series of figures to highlight the important stages of relaxation and the energy level dynamics corresponding to the observed macroscopic behavior. Figures 7(a) and 7(b) show the evolution of vibrational temperature and mole fraction of atomic and molecular oxygen as a function of time. The system eventually reaches thermal and chemical equilibrium at *T* = 7000 K, where the equilibrium mole fraction of the molecular oxygen is about 0.6%. Figures 8(a) and 8(b) indicate the simulation times [(a)–(e)] at which vibrational distribution functions (Fig. 9) and state-specific dissociation/recombination rates (Fig. 10) are sampled. From Fig. 7(a), the evolution of the system can be divided into 4 distinct stages: (i) vibrational heating, (ii) quasi-steady state (QSS), (iii) final relaxation, and (iv) equilibrium. In stage (i), the energy added to the system by the isothermal reservoir goes into heating the vibrational mode. This is seen as a rapid initial increase in T_{v} at the beginning of the simulation. This increase in vibrational temperature is also seen in Figs. 9(a)–9(c) as distribution broadening from t = 0 s to t = 1.2 × 10^{−7} s.

Figure 7(a) shows that all models spend approximately the same duration of time in the initial vibrational heating stage, which is reasonable since Z_{v} in both the TCE and QK models is determined from the STS vibrational relaxation rate. The system then enters into the quasi-steady state (QSS) region in stage (ii), where the dissociation rate becomes comparable to the vibrational relaxation rate. This results in a competition between dissociation and vibrational heating, which work to deplete and replenish the vibrational states, respectively, and a nearly stationary vibrational distribution from t = 1.2 × 10^{−7} s to t = 2.2 × 10^{−7} s is observed in Fig. 9(e). It is noted that the TCE and QK phenomenological models exhibit strongernet dissociation, resulting in comparatively short QSS regions compared to the STS model, whose QSS region is longer than TCE and QK by a factor of ten and four, respectively. These results are consistent with the fact that the phenomenological models overestimate the dissociation coefficient for *T*_{v} < *T*, as shown in Fig. 6(a). The total time to reach equilibrium for the STS model is approximately 3.5 times that of the TCE models and 1.5 times that of the QK model.

The evolution of the vibrational distribution function and the state-resolved dissociation and recombination rates of STS model are plotted in Figs. 9 and 10. The local vibrational temperature is computed by the energy-equivalent approach.^{13} Again, the instances in time indicated in the figure captions correspond to times indicated by labels (a)–(e) in Fig. 8. At *t* = 4.0 × 10^{−9} s, the VT transition pumps the molecules from low to high vibrational states and leads to theoverpopulation at the high-lying vibrational states [Fig. 9(a)]. At *t* = 1.2 × 10^{−8} s, the vibrational distribution continues to broaden and becomes closer to the Boltzmann distribution. The start and end of the QSS region correspond with *t* = 1.2 × 10^{−7} s and *t* = 2.2 × 10^{−7} s, respectively. At this stage, dissociation becomes more important and starts to deplete the tail population of the distribution. At the same time, VT transitions are working to populate the tail, resulting in a competition between the two processes. As seen in Figs. 9(c) and 9(d), the populations of the vibrational states at these times are nearly stationary, and the vibrational temperature [indicated by the slope of the green and blue lines fit to the data in (c) and (d)] is nearly constant [Fig. 9(f)]. It is also noted that the vibrational distribution function can be represented by two distinct temperatures, where the vibrational temperature at high-lying states is about 2000 K lower than the low-lying states during QSS. The corresponding state-specific dissociation and recombination rates during QSS from the STS model are shown in Figs. 10(c) and 10(d). During QSS, the state-specific dissociation and recombination rates are nearly stationary. It is also noted that at these times, the dissociation/recombination rates of the high-lying vibrational states are already balanced during QSS (indicated by the overlapping of the red and blue symbols). As expected, these rates are balanced across all vibrational states at equilibrium as required by detailed balance [Fig. 10(e)].

#### 2. Vibrational cooling case (recombination regime)

The final case presented is the 0-d vibrational cooling case, in which the system is held at a constant translational temperature, *T* = 2000 K, and the initial vibrational and rotational temperatures are set to *T*_{v} = *T*_{rot} = 7000 K. The molecular/atomic oxygen gas mixture is initialized with mole fractions of 0.2 and 0.8, respectively, such that the system undergoes net recombination during relaxation toward chemical equilibrium. The results of the simulation are again shown in a series of figures to highlight the energy level dynamics corresponding to the observed relaxation behavior. Figures 11(a) and 11(b) show the temperature and mole fraction as a function of time during the vibrational cooling process. The evolution of the vibrational distribution function and the state-resolved dissociation/recombination rates are plotted in Figs. 12 and 13.

Unlike the heating case, there is no apparent plateau in the vibrational temperature profile during cooling [Fig. 11(a)], and a QSS region is not observed for any of the models. The competing processes which would contribute to QSS in the cooling case are VT transitions, which work to deplete the vibrational distribution tail, and recombination, which works to replenish the tail population. The recombination rate, however, is not large enough to balance the rate of depletion by VT transitions, which is comparatively faster. Recall that for the vibrational heating case, QSS was observed in which the dissociation rates shown in Fig. 10 were on the order of 10^{30} s^{−1}, while the recombination rates shown in Fig. 13 during vibrational cooling are more than an order of magnitude lower. This suggests that the VT transitions dominate the thermal relaxation process, and the temperature profile observed in Fig. 11(a) resembles a relatively rapid thermal relaxation, with a comparatively slow chemical relaxation. This is also reflected in Fig. 11(b), which shows the mole fractions of molecular and atomic oxygen as a function of time. Relaxation toward chemical equilibrium requires more time in the cooling case, owing to slow recombination rates.

It is observed that the STS model reaches chemical equilibrium more slowly than the phenomenological TCE and QK models. This is due to the vibrational favoring of dissociation/recombination in the STS model, as shown in Fig. 5. Themolecules at high-lying vibrational states tend to be dissociated, while the recombining molecules have a preference to enter into the high-lying vibrational states for the STS model. Therefore, the molecules produced by recombination are very readily re-dissociated. For the phenomenological models, the trend of the state-resolved dissociation coefficients is similar to the STS model; however, the recombining molecules tend to enter into low-lying vibrational states. Therefore, the net recombination rate of STS model is lower than those of the phenomenological models and requires more time to reach chemical equilibrium.

The evolution of the vibrational distribution function and the state-resolved dissociation and recombination rates are plotted in Figs. 12 and 13, respectively. The labels (a)–(d) correspond to the results at time 2.5 × 10^{−8} s, *t* = 3.5 × 10^{−8} s, *t* = 9.5 × 10^{−7} s, and (d) *t* = 8.5 × 10^{−4} s, respectively. During the vibrational cooling stage, the vibrational temperature drops rapidly due to VT transitions which drive the tail population to lower-lying vibrational energy states. In the first stages of rapid vibrational cooling [Figs. 12(a) and 12(b)], the low-lying vibrational states relax to lower temperatures first. At these early times, both dissociation/recombination reactions are important due to the relatively high vibrational temperature and the high mole fraction of atomic oxygen, respectively, as shown in Fig. 13(a).

As the vibrational temperature is further decreased and approaches the translational temperature, the vibrational distribution becomes closer to Boltzmann; however, the very high-lying vibrational states are in strong non-equilibrium, as shown in Figs. 12(b) and 12(c). For instance, in Fig. 12(c), the vibrational distribution function can be represented by three distinct temperatures. The vibrational temperature of the intermediate energy states is much hotter (≈14 000 K) due to recombination, while the vibrational temperature of the tail (≈2000 K) is close to the translational temperature due to the relatively strong dissociation in the tail. The overpopulated high-lying states assist the dissociation; therefore, even though the mole fraction of molecular oxygen is low (*χ*_{O} ≈ 0.2) and the vibrational temperature is low (≈2000 K), there are still appreciable amounts of dissociation. There is an obvious overpopulation at high-lying vibrational states in Figs. 13(b) and 13(c), which is mainly due to the up-pumping by the recombination. In Fig. 13(d), the gas has nearly reached chemical equilibrium.

## V. CONCLUSION

In this work, a vibrational state-specific model for dissociation and recombination reactions in DSMC is presented based on QCT cross-sectional calculations. The two-step binary collision framework is applied for recombination reactions in which both microscopic reversibility and state-specific detailed balance are used to constrain the orbiting pair and stabilization cross sections. The dissociation and recombination rate coefficients are discussed in both equilibrium and non-equilibrium conditions. In comparing the STS model with the phenomenological models, in the equilibrium condition, TCE and STS show good agreement in the overall reaction rate coefficients.

The TCE, QK, and STS models were compared in 0-d vibrational heating and cooling cases. The vibrational heating case shows a distinct QSS region during relaxation, while QSS is not established during vibrational cooling. In all cases, the vibrational distribution function exhibited strong non-equilibrium behavior during the vibrational heating and cooling processes. In the heating case, the phenomenological models start to dissociate ahead of the STS model because of the overestimation of the dissociation coefficient at low-lying vibrational states, which also leads to a shorter time to reach thermal and chemical equilibrium. In the cooling case, the phenomenological models again reach chemical equilibrium faster than STS model. A stronger dissociation rate is observed in the STS model which is due to the assistance of the recombination, resulting in a longer time to reach equilibrium. Finally, the balance between dissociation and recombination in the STS model is first achieved in the high vibrational states, which indicates that the net dissociation/recombination during the heating/cooling processes is driven by the low-lying vibrational states. Because of its importance at low temperatures for O_{2} relaxation,^{15} rotational nonequilibrium should be addressed by extending the current work to a rotationally resolved cross-sectional database.

## ACKNOWLEDGMENTS

This work was supported by an Early Career Faculty grant from NASA’s Space Technology Research Grants Program and the NASA Space Technology Research Fellowship, and by the Air Force Office of Scientific Research under award number FA9550-17-1-0127. The authors gratefully acknowledge Dr. F. Esposito for providing the QCT database and for many useful discussions.

## REFERENCES

_{2}/O flows

_{2}–N

_{2}system

_{2}+ N

_{2}dissociation reactions

_{2}(

^{1}Σ

_{g}

^{+})–N(

^{4}S

_{u}) system in hypersonic flows

_{2}-N collisions using direct molecular simulation

_{2}–O collisions

_{2}state-to-state vibrational relaxation and dissociation rates based on quasiclassical calculations

_{2}+ O system

_{3}from a multiproperty fit to ab initio calculations, and to experimental spectroscopic, inelastic scattering, and kinetic isotope thermal rate data

_{2}

_{2}

_{2}–O mixtures

_{2}molecules by atomic oxygen

_{2}(v) → 3O

_{2}collisions

_{4}

_{4}

_{2}–Ar, Kr, and Xe

_{2}collisions

_{P})–O(3

_{P}) interaction