Construction of a coarse-grain quasi-classical trajectory method. I. Theory and application to N₂–N₂ system

Modeling the environment surrounding vehicles flying at hypersonic speeds into Earth’s atmosphere is a very challenging task due to the complex interaction between multiple physical processes. In hypersonic flows, the shock heated gases can reach temperatures on the order of tens of thousands of degrees Kelvin. At these extreme temperatures, high energy collisions between the atoms and molecules can lead to excitation of the internal energy modes, dissociation and eventually ionization of the atoms and molecules in the gas.^1,2 Thus, at high temperatures, air comprises a large number of chemical components including both neutral and ionized species. This work focuses on the modeling of energy transfer and dissociation of molecular nitrogen in its ground electronic state. Molecular nitrogen is the most abundant species in Earth’s atmosphere and must be accurately modeled for the design of high-speed vehicles. In particular, understanding the chemical processes due to N₂–N₂ interactions is critical to predicting the initial excitation and dissociation in air flows.

Because the time scale of the flow and the characteristic time of the chemical reactions are comparable, flow calculations require the use of a coupled approach, where chemical kinetics and flow governing equations are solved simultaneously.^2,3 Early work on non-equilibrium modeling focused on the empirical multi-temperature models, which were constructed on the assumption that the population of each internal (rotational, vibrational, or electronic) energy mode can be described by means of a Maxwell-Boltzmann distribution at a specific temperature.^4–7 The vast majority of multi-temperature models rely on the assumption of equilibrium among the rotational and translational modes at T and equilibrium among the vibrational, electronic, and electron temperature at T_V. This assumption stems from the observation that, at low temperatures, rotational excitation occurs very quickly, given the reduced energy spacing among the rotational states. Therefore, the rotational relaxation has reached equilibrium with translation long before the vibrational, electronic, and electron energy modes are significantly excited. Finally, it is worth mentioning that these multi-temperature models require the knowledge of numerous parameters which are calibrated by matching experimental data, collected at ground-based facilities.^8–11 

State-to-state (StS) models provide a physics-based alternative to multi-temperature models. In this approach, internal states of atoms and molecules are treated as individual pseudo-species, and their densities are evaluated by considering the detailed kinetic processes.¹² This can be done to varying degrees of accuracy, e.g., electronic,^13–18 vibrational, and rovibrational state-to-state.

Vibrational StS models generally assume either rotationless molecules or equilibrium of the rotational mode at the local translational temperature and explicitly solve for the individual vibrational state populations.^19–30 These models have been used extensively for studying non-equilibrium excitation and dissociation processes. However, the accuracy of these models is limited by the validity of the assumptions made about the rotational mode and the availability of accurate kinetic data for the individual processes. Previous work on computing the kinetic data focused on semi-empirical models such as the ladder-climbing model,^31,32 or semi-classical models such as the Schwartz-Slawsky-Herzfeld (SSH) model,^33–35 and the Forced Harmonic Oscillator (FHO) model.^36–43 Early work on characterizing StS kinetics for three-body systems made use of surprisal analysis based on the work of Procaccia and Levine.^44–46 

Finally, rovibrational StS models treat all internal rovibrational states individually as separate species. These models are prohibitively expensive for two or three-dimensional simulations but can be used to study excitation and dissociation in zero-dimensional reactors and 1D Euler flows for validation purposes.^47–49 Moreover, a full rovibrational simulation for even simple systems (e.g., atom-diatom interaction) requires a huge amount of kinetic data, which are largely unavailable.

Recent advances in computational quantum chemistry have enabled work detailing the microscopic kinetics for several systems.^{23–25,27,28,47,48,50–67} This is accomplished in two steps: first determining the interaction potential between two particles and then running scattering calculations to compute transition probabilities. The interaction potential, also known as the potential energy surface (PES), is computed by solving the electronic Schrödinger equation at many geometric arrangements of the atoms.^68–71 Next, the transition probability can be computed by simulating many collisions between two particles, using the PES computed from quantum mechanics. In this work, the collisions are simulated using the quasi-classical trajectory (QCT) method, which assumes that the collision occurs classically, solving Hamilton’s equations of motion for the nuclei, but making use of the PES from quantum mechanics to obtain the force between the particles throughout the course of the collision.^71–74 The outcome of many collisions is analyzed to obtain the probability that a given reaction occurred. This approach has been used to determine the rovibrational StS kinetics for atom-diatom systems, including the $N2(X1Σg+)−N(Su4)$ system. Recent work by Esposito et al.⁷⁵ coupled QCT with scaling laws to determine the kinetics for the $N2(X1Σg+)−N2(X1Σg+)$ system. However, a full QCT characterization of the kinetic data for diatom-diatom systems is impractical because the number of possible transition is on the order of 10¹⁵.

As an alternative to the computationally expensive rovibrational StS model, the coarse-grain model (CGM) was developed with the objective of capturing the same physics of StS models at a fraction of the computational cost.^76–86 The coarse-grain model, first developed for studying the kinetics of the electronic states of atomic species,^14,87 has been recently extended to molecular species and applied to study the non-equilibrium energy transfer and dissociation processes for the $N2(X1Σg+)−N(Su4)$ system.^80–82,86 In the CGM, energy states are lumped together, and the actual population distribution is approximated with a local reconstruction function derived by maximizing the entropy within the group. This approach is similar to that proposed by Haug et al.,^44,45 which lumps energy states, assuming states in the same group to be in equilibrium. Previous work on the CGM has focused on systems with known microscopic StS kinetics. However, for more complex systems, such as interactions between two diatoms, detailing all the microscopic reactions becomes impossible due to the large number of kinetic channels.

The objective of this work is twofold: first to present the Coarse-Grain Quasi-Classical Trajectory (CG-QCT) method, a novel approach to constructing the coarse-grain model for diatom-diatom systems, and second to use the reduced order model constructed to study non-equilibrium energy transfer and dissociation of N₂ molecules. The CG-QCT method couples the coarse-grain model with the quasi-classical trajectory calculations to directly obtain the grouped kinetic properties for the CGM, providing a framework for constructing reduced order models for more complex chemical systems. In the proposed method, the quasi-classical trajectory calculations require only information about the energy level grouping and reconstruction within each group to determine the group average reaction rates. The CG-QCT method is used to study the energy transfer and dissociation processes of molecular nitrogen, considering only the interaction between two nitrogen molecules in their ground electronic states, the $N2(X1Σg+)−N2(X1Σg+)$ system. Two grouping strategies are considered in this work: energy-based bins and vibrational-based bins. The energy-based grouping strategy lumps energy states with similar internal energy together, regardless of their rotational and vibrational quantum numbers. Vibrational-based binning lumps states together which share the same vibrational quantum number. This is identical to a vibrational StS model assuming the rotational states are in equilibrium at a rotational temperature. The two grouping strategies are compared in a zero-dimensional isothermal and isochoric reactor simulation under severe non-equilibrium conditions.

This work is divided in two parts: The present paper focuses on the model development and application of grouping strategies to the N₂–N₂ system. In particular, the differences in the results due to the grouping scheme are highlighted. Insights on the strengths and shortcomings of each model are obtained by direct comparison against experimental data. Paper II⁹⁵ presents a comparison of the results obtained with the CG-QCT method with the ones obtained with the direct molecular simulation (DMS) method.^88–95 

This paper is organized as follows: first the physical models, the CGM and CG-QCT methods, will be presented in Sec. II, followed by the results and comparison against the experimental data in Sec. III, a discussion of the physics captured in the results in Sec. IV, and finally the conclusions will be discussed in Sec. V.

In this section, the model reduction framework is described in detail. In Subsection II A, the coarse-grained modeling approach is recalled. This is followed by the Coarse-Grain Quasi-Classical Trajectory (CG-QCT) method in Subsection II B. Finally, the details on the potential energy surface for the $N2(X1Σg+)−N2(X1Σg+)$ system are presented in Subsection II C.

State-to-state modeling requires the solution of a large system of equations since different internal energy levels of reacting particles are treated as distinct chemical species. This description although exceedingly accurate has limited applicability to multi-dimensional computational fluid dynamics (CFD) problems due to its prohibitive cost. The present work overcomes this bottleneck through the use of reduced-order models based on the coarse-grain method. The coarse-grain method has been previously applied to different flow problems and chemical systems. The model reduction procedure can be broken down into two key steps:^80,82

1. Local representation and reconstruction: the rovibrational energy states of $N2(X1Σg+)$ are grouped together according to their energy, and the actual distribution within each group is prescribed assuming an equilibrium distribution within the bin at the local translational temperature.

2. Macroscopic governing equations: the macroscopic governing equations are obtained by taking moments of the microscopic master equations and applying the reconstruction function based on the assumption of equilibrium within each bin.

In the present work, the distribution of states within each bin is assumed to be an equilibrium Boltzmann distribution at the local translational temperature,

where quantities n_i, E_i, and g_i are, respectively, the population, energy, and degeneracy of state i, n_p represents the population of group p, $Fpi(T)$ denotes the distribution of states i within group p at translational temperature T, k_B represents the Boltzmann constant, and Q_p(T) denotes the partition function of group p, given by

where the set $Ip$ denotes the set of states contained in group p. The general coarse-grain method framework introduces a bin internal temperature, T_p by imposing constraints to maximize the entropy within the group.^76,78

One of the most important considerations in the construction of the CGM is how to group the energy states. Previous work has focused on an energy-based grouping method, by which states were grouped according to increasing energy.^76,82 In this work, two grouping strategies were considered: uniform width energy-based grouping and vibrational-based grouping.

Uniform width energy-based grouping considers the generic diatomic molecule, M, with internal structure, made up of atoms A and B. The energy states of molecule M are first split into bound and quasi-bound (or pre-dissociated) states. Then, given a number of bound and quasi-bound groups, denoted by N_B and N_QB, respectively, the energy width of the bins can be determined by

where ΔE_B and ΔE_QB denote the energy width of the bound and quasi-bound bins, respectively, E_A and E_B represent the formation energy of atoms A and B, and E_⋆ represents the energy of the rovibrational state with the largest energy. A schematic of the energy-based grouping strategy is shown in Fig. 1(a), considering only three energy bins. In this schematic, the different colors denote different vibrational states, which in the energy-based grouping strategy are lumped together. In this work, the grouping method will be applied to the $N2(X1Σg+)−N2(X1Σg+)$ system using 60 bins (40 of bound states, and 20 of quasi-bound or pre-dissociated states). The resulting bin widths are ΔE_B = 0.24 eV and ΔE_QB = 0.26 eV. However, it is worth recalling that the proposed framework is general and, as such, it can be applied to any diatom-diatom system.

Vibrational-based grouping instead considers the vibrational quantum numbers when constructing the groups. In this approach, all the rotational states which share a vibrational quantum number are grouped together. A schematic of this grouping strategy is shown in Fig. 1(b) for three example vibrational states, where the red states belong to $v=0$⁠, the blue states belong to $v=1$⁠, and the green states belong to $v=2$⁠. The resulting bins contain states which span the entire energy spectrum (e.g., the $v=0$ bin will contain states with all possible rotational quantum numbers or energies).

The second step of the CGM consists of constructing the macroscopic governing equations to determine group properties n_p. This is accomplished by summing the microscopic moment equations of order m for states (i, j, k, l) contained in groups (p, q, r, s),

where the k symbols stand for the rate coefficients for excitation, $kij−klE$⁠, de-excitation, $kkl−ijE$⁠, combined excitation-dissociation, $kij−kED$⁠, combined excitation-recombination, $kk−ijER$⁠, dissociation, $kijD$⁠, and recombination, $kijR$⁠. Quantities n_A and n_B denote the number densities of the atomic constituents of molecule M, sets (⁠$Ip$⁠, $Iq$⁠, $Ir$⁠, $Is$⁠) contain the states within groups (p, q, r, s), and the set $I$ stores all the groups of molecule M. Taking the zeroth-order energy moment of Eq. (4) yields the conservation equations of group mass,

Equation (5) constitutes the general conservation equations for the linear CGM. The capitalized symbols K denote the group reaction rate coefficients. The solution of Eq. (5) allows for the determination of the group number densities, which can be then used to reconstruct the population of elementary states via Eq. (1).

The QCT method was used in this work for the scattering calculations. This method makes use Hamilton’s equations of motion to describe the dynamics of the collision but uses the ab initio PES generated by solving the Schrödinger equation at various geometries (i.e., configurations of the nuclei). The QCT code used in this work is a modified version of the Vectorized Variable Timestep Trajectory (VVTC) code written by Schwenke at NASA Ames Research Center.⁷¹ The initial states can either be explicitly defined, to determine rate coefficients for specific rovibrational states, or sampled from a thermal distribution at a specific internal temperature. The translational energy of the collisions was sampled from a thermal distribution at prescribed translational temperature values. All other collision parameters (e.g., impact parameter, orientation, vibrational phase, etc.) were sampled from an appropriate distribution using Monte Carlo averaging. The maximum impact parameter was taken to be 7.41 Å, conservatively chosen based on previous work by Valentini et al.⁹² Further details on the QCT procedure can be found in Karplus et al.⁷² and Jaffe et al.⁷⁴ 

After simulating a large number of collisions, the reaction cross section can be determined from the probability that a given event occurs. For a general excitation reaction from states (i, j) to states (k, l), the cross section, $σij−klE$⁠, can be written as

where $Pij−klE$ is the probability of the transition (i, j) → (k, l). This quantity is computed from the ratio between the number of trajectories which resulted in the transition, $Nij−klE$⁠, and the total number of trajectories simulated, N_tot. The error due to sampling is reduced by the square root of the number of samples (i.e., $dσij−kl∝1/Nij−kl$⁠). Therefore, increasing the number of samples resulting in a given outcome improves the aleatoric uncertainty of a given rate. Quantity b_max denotes the maximum allowed value for the impact parameter. Once the cross section is known, the reaction rate coefficient can be determined by taking a thermal average over the relative velocities of the reaction cross section

where quantities g and χ denote, respectively, the relative velocity of the colliding particles and the symmetry factor (which equals 1 when the colliding partners are different and 2 when they are identical). The symbol f^M(T) represents a Maxwellian velocity distribution at temperature T.

The CG-QCT method provides the link between the CGM and the QCT simulations and enables the construction of the CGM for diatom-diatom systems. Previous work on the CGM obtained the group-average kinetic parameters by summing the microscopic data and applying the appropriate weighting to obtain the correct distribution within the group. For example, the grouped cross section for an excitation reaction from groups (p, q) to groups (r, s), $σpq−rsE$⁠, can be written as the sum of the microscopic probabilities, $Pij−klE$⁠, weighted by the distribution of states within groups p and q, $Fpi$ and $Fqj$⁠, multiplied by the area of the impact ring, $πbmax2$⁠,

In the CG-QCT method, the summation over the initial states (i, j) is evaluated by sampling these states from groups (p, q) populated according to the distribution functions $(Fpi,Fqj)$⁠. This results in the rate coefficients from pair of groups (p, q) populated according to the prescribed distribution functions $(Fpi,Fqj)$ going to states (k, l). Then the summation over the final states is performed after the QCT calculations are complete. This eliminates the need to directly compute all the elementary cross sections, $σij−klE$⁠, which are impossible for a rovibrational StS model for the $N2(X1Σg+)−N2(X1Σg+)$ system. From the grouped cross sections, the grouped reaction rate coefficients can be obtained in the same way as the corresponding microscopic rate coefficients

For the remaining reactions considered in the diatom-diatom CGM, the cross sections for combined dissociation-excitation and dissociation reactions can be written as

The reverse reaction rate coefficients can be determined from the forward grouped rates by applying micro-reversibility. Details on computing the reverse reaction rate coefficients can be found in the  Appendix.

In this work, the CG-QCT method has been applied to the $N2(X1Σg+)−N2(X1Σg+)$ system to study the dynamics of excitation and dissociation in molecular nitrogen. The N₂ molecule in its ground electronic state, $N2(X1Σg+)$⁠, has 9390 rovibrational states, with vibrational quantum numbers $v$ and rotational quantum numbers J. The rovibrational states were determined using quantum mechanics calculations using the WKB approximation⁷¹ with the $N2(X1Σg+)$ potential developed by Le Roy et al.⁹⁶ After sorting the rovibrational states according to increasing energy, these can be assigned a global index i, where the energy and degeneracy of the state can be written E_i and g_i. Of these states, most of them have energy below dissociation energy of 9.75 eV, referred to as the bound states. The remaining levels have energy above the dissociation energy but below the J-dependent centrifugal barrier.

The potential energy surface for the $N2(X1Σg+)−N2(X1Σg+)$ system used in this work was developed by NASA Ames Research Center.^59–61 The nuclear positions were divided into regions where both $N2(X1Σg+)$ molecules had bond lengths near equilibrium, and where one or both $N2(X1Σg+)$ molecules had a bond length far from equilibrium. In the first group, where both molecules have bond lengths near equilibrium, electronic structure calculations were performed using the closed-shell coupled-cluster-single-double method to parametrize the wavefunction.^68,69 For this region, calculations were performed for 3821 nuclear geometries. In the second region, where one or both bond lengths are far from equilibrium, calculations were performed for 325 nuclear geometries using the multi-reference averaged-coupled-pair-function method.⁷⁰ Further details on the calculation of the PES can be found in Refs. 59–61.

In this section, an isochoric reactor model is used to study the non-equilibrium dissociation process of nitrogen modeled using the CG-QCT method. In all the simulations, the gas is initially composed of cold nitrogen molecules, populated according to a Boltzmann distribution at T_I = 300 K with a density of 1.2 kg/m³ corresponding to a pressure of 1 atm. Under these assumptions, the initial population of the groups is given by

where $nN2$ denotes the total number density of nitrogen molecules given by $nN2=∑p∈Inp$⁠, E_p is the formation energy of group p, $ΔEip$ is the energy of state i relative to the formation energy of group p such that $Ei=Ep+ΔEip$⁠, T_I is the initial internal temperature, and $Q̃p(T)$ is the partition function of a group relative to the group formation energy: $Q̃p(T)=∑i∈Ipgiexp−ΔEip/(kBT)$⁠. At the beginning of the simulation, the translational temperature of the reactor is instantaneously raised and held constant, driving the gas out of equilibrium. Four different translational temperatures are considered: 10 000 K, 13 000 K, 20 000 K, and 25 000 K.

All the simulations are carried out using the energy-based and vibrational-based grouping strategies discussed in Sec. II. The former comprises 60 groups, 40 bound and 20 quasi-bound groups of equal energy width, while the latter considers the 61 vibrational energy levels of N₂ as separate groups. In both cases, the collisional processes between two nitrogen molecules considered are excitation, combined excitation-dissociation, and double dissociation. Collisions between atomic and molecular species are ignored.

The section is organized as follows: First, the time evolution of the distribution function predicted by the CG-QCT model will be presented, followed by an analysis of the dissociation and energy transfer processes. Finally, the results obtained with the two different grouping strategies are compared against the available experimental data.

At the beginning of the simulation, only the low-lying energy groups are significantly populated, given the initial values of T_I. With time, the random motion of molecules brings about collisions, thus enabling the transfer of kinetic energy into rotational and vibrational energy. Figure 2(a) shows the distribution of energy groups at various times in the 10 000 K reactor simulation. Early in the relaxation process, t = 10⁻¹⁰ s, the distribution is still significantly colder than the final equilibrium distribution, indicating that the gas is still in the midst of the relaxation process. This phase is completed at about t = 10⁻⁹ s. Between t = 10⁻⁸ s and t = 10⁻⁷ s, the distribution is frozen, indicating that the gas has reached the quasi-steady-state (QSS) distribution. After t = 10⁻⁷ s, the high energy groups are replenished through recombination, and the gas approaches the final equilibrium distribution after t ≈ 10⁻⁵ s. Figure 2(b) shows the time evolution of the population for a subset of the groups. In this figure, the different phases of the thermochemical relaxation can be clearly observed along with the plateau in the population densities which represent the QSS of the gas.

Figure 3(a) shows the distribution of vibrational groups at various times in the 10 000 K reactor simulation. Because of the difference in grouping strategies, the distribution looks very different from the energy binned results. Early in the relaxation process, t < 10⁻⁹ s, the low energy vibrational states (⁠$Ev<1.5$ eV) appear nearly frozen at a colder temperature than the higher energy states. This bi-modal distribution persists until the QSS is reached at t = 10⁻⁷ s. The QSS is clearly observable in Fig. 3(b), which shows a narrow plateau in the population of the high vibrational energy level populations around this time. As observed for the other grouping strategy, the distribution shows significant deviation from equilibrium. It is important to note that the initial energy transfer process predicted by the vibrational-based grouping model is significantly slower than the one predicted by the energy-based grouping method.

To better illustrate the differences between the two grouping strategies on the microscopic level, the reconstructed distributions are shown in Fig. 4 at two different times in the 10 000 K reactor. Figure 4(a) shows the distribution during the relaxation process, at t = 1.173 × 10⁻¹⁰ s. At this instant, the total internal energy contained in the molecules is the same; however, due to different assumptions in the grouping strategy, the rovibrational distributions are very different. Since the energy-based grouping assumes equilibrium over a small range of energies (i.e., each energy bin is only 0.25 eV in width), the distribution is mostly continuous across the energy spectrum, with each bin approximating only a reduced number of levels in a narrow energy range. In contrast, the vibrational binning strategy assumes equilibrium at 10 000 K among all rotational states within a given vibrational state. This results in the strand like structure observed in Fig. 4(a).

Figure 4(b) shows the reconstructed distribution of states during the dissociation process for both grouping strategies. In both groupings, the distribution of low energy states is very similar and approaches equilibrium. However, nearing the dissociation energy (9.75 eV), the distribution predicted by the energy binning strategy deviates from the Boltzmann distribution. The population predicted by the vibrational grouping model is more complicated. The high-lying vibrational levels contain a reduced number of rotational levels and therefore will capture the depopulation of the high-lying rovibrational states. On the contrary, the low-lying vibrational levels contain a large number of rotational levels characterized by large rotational quantum number. These levels are forced to be in equilibrium with the low-lying rovibrational states by the averaging procedure, thus creating an artificial overpopulation in the tail region of the distribution.

The two models provide a very different representation of rovibrational relaxation. The vibrational grouping appears adequate for the description of the low-lying vibrational level, characterized by a small rotational quantum number, where the mode separation is clearly significant. On the contrary, the energy binning strategy seems more adequate for the description of the relaxation of the high-lying states characterized by low vibrational quantum numbers and high rotational energy content.

The mole fraction of atomic nitrogen predicted by both grouping strategies at various temperatures is shown in Fig. 5. At 25 000 K, the two grouping strategies predict similar dissociation rates, even if the onset of dissociation occurs significantly earlier with the vibrational specific model. As the temperature decreases, the two grouping strategies diverge. At 10 000 K, the energy-based groups predict faster dissociation, with a shorter incubation period. The discrepancy between the two grouping models is due to the treatment of the rotational energy mode. In the vibrational grouping model, the rotational levels are assumed to be populated according to a Boltzmann distribution at the translational temperature. This assumption hinders dissociation at low temperature because it artificially imposes equilibrium between low and high energy states [e.g., (⁠$v$⁠, J) = (0, 0) and (⁠$v$⁠, J) = (0, 270)], which is clearly incorrect. On the contrary, the energy-based grouping is able to capture the non-equilibrium between high and low-lying rotational levels since they belong to different groups.

The importance of the high energy states for predicting dissociation is highlighted in Fig. 6. This shows the fraction of molecules which dissociate from a given group when the molecules are in the QSS condition, given by

where p and q are the groups of the dissociating and exciting molecules, respectively, and P_diss(E_p) is the probability that a molecule in group p dissociates. Figure 6(a) presents the results obtained with the energy-based grouping model. The distribution clearly shows that molecules climb nearly to the dissociation energy before dissociating. This effect is more pronounced at low temperatures, where there is less energy available in the translational mode to facilitate dissociation from the lower energy groups. Thus, the importance of the quasi-bound states increases with decreasing temperature. Table I gives the percentage of dissociation occurring from quasi-bound states as predicted by the model. At low temperatures, nearly 50% of the dissociation occurs from quasi-bound states, highlighting their importance for the prediction of the dissociation process. This conclusion is similar to that reached by Bender et al.⁹⁷ for the $N2(X1Σg+)−N2(X1Σg+)$ system who observed that up to 58% of dissociation events come from trajectories with at least one quasi-bound molecule, as well as Panesi et al. for the $N2(X1Σg+)−N(Su4)$ system.⁴⁸ 

The vibrational grouping method, shown in Fig. 6(b), predicts very different behavior. At the lowest temperature, 10 000 K, the molecules generally climb to higher vibrational states (≈7–8 eV) before dissociating, as predicted by the ladder climbing model. However, at the highest temperature, 25 000 K, the molecules are much more likely to dissociate from low energy vibrational states (≈1–4 eV). This behavior can be explained as follows: at high temperatures, the low vibrational states are rotationally excited (given the assumption of rotational equilibrium), and the high lying rotational states are therefore significantly populated, thus significantly contributing to dissociation. In other words, the dissociation energy of a molecule decreases as J increases. However at low temperatures, the Boltzmann weighting favors the lower lying rotational energy levels, characterized by significantly lower dissociation probability, thus hindering dissociation from the high-lying rotational states of a given group.

To gain more insight on the relative contribution of rotational and vibrational energy to dissociation, Fig. 7 shows the fraction of energy lost by the two modes during dissociation at 10 000 K, normalized by the total internal energy lost, given by

where $ΔErot,pdiss$ and $ΔEvib,pdiss$ denote, respectively, the energy lost in dissociation from the rotational and vibrational modes of group p, $ΔEtotdiss$ denotes the total energy lost during dissociation, and $Eirot$ and $Eivib$ denote, respectively, the rotational and vibrational energies of state i. It is interesting to observe that the energy-based bins predict a similar contribution between rotation and vibration, with the latter contributing slightly more. The percentage of energy lost from each mode in the energy-based grouping model is also summarized in Table I, highlighting that at high temperature, the energy loss from each mode becomes comparable. In contrast, the vibrational grouping provides a completely different picture: the contribution of rotation does not exceed 10% of the total internal energy lost. In other words, as a result of the assumptions made, the importance of rotation on the kinetics is downplayed at low temperatures in the vibrational-based averaging. The work of Bender et al.⁹⁷ obtained similar results to what is observed in the vibrational grouping method, which finds that the rotational energy contribution to dissociation is significantly smaller than the vibrational contribution. However, this may be an artifact of the assumptions made in the model (e.g., rotation and vibrational modes are decoupled).

The group distribution during QSS for the energy and vibrational-based bins is shown in Fig. 8 for various temperatures. In both cases, significant departures from the QSS temperature Boltzmann distributions are observed. The low energy groups appear to be close to equilibrium, whereas the distribution of the high energy tail appears significantly depopulated. The behavior of the distribution is consistent with the dissociation probability functions discussed above. The high energy molecules, more likely to dissociate, are responsible for the departures of the distribution from equilibrium. Between the two models, the vibrational energy model seems to predict lower non-equilibrium of the high energy tail.

Figure 9 shows the local dissociation rate profiles during the relaxation. These have been calculated from the production rate of atomic nitrogen as follows:

The rate profiles are monotonically increasing in the early stage of the relaxation for both models, until the onset of the QSS distribution, responsible of the formation of a plateau in the profiles. A comparison of the nitrogen concentration profiles and the dissociation rates shows that, across a range of conditions analyzed, the entire dissociation process proceeds under QSS conditions for the N₂–N₂ processes. This was not the case for the N₂–N system, where significant part of the dissociation occurred in non-QSS conditions.⁴⁸ The comparison of the rate profiles obtained with the two different grouping strategies shows that the QSS dissociation rate predicted at low temperatures by the vibrational grouping is approximately half of those predicted by the energy-based grouping. Moreover, the onset of the QSS region is significantly delayed for the vibrational grouping model. This indicates that not only is the dissociation process different between the two groupings but also the energy transfer process has significant differences. As the temperature increases, the vibrational grouping QSS dissociation rate overshoots the rate predicted by the energy groups. This is due to the over-population of high energy rovibrational states in the vibrational grouping model discussed earlier.

The internal energy of the molecules, shown in Fig. 10, can be computed as

The energy-based grouping starts with very low internal energy and excites quickly to the QSS, by 10⁻⁸ s at all temperatures. At the low temperatures, the total internal energy of the molecules is the same as the thermal internal energy, indicating that the QSS distribution is not much different from the thermal Boltzmann distribution. This is not surprising since most of the energy is contained in the low energy states, whose population is close to equilibrium at low temperatures, and the high energy states do not significantly contribute to the internal energy content. At higher temperatures, the QSS has significantly lower internal energy from the final equilibrium value, indicating the presence of stronger non-equilibrium effects.

The vibrational specific bins exhibit very different behavior throughout the relaxation process. Despite the initial excess of internal energy, due to the assumption of T = T_R, the vibrational grouping model predicts significantly slower relaxation across the range of temperatures, if compared with the energy-based model.

Further insights on the energy transfer processes can be obtained by computing the second order moment of the transition rates for both models. To this aim, the master equation is expressed as diffusion or Fokker-Planck type equation,^2,98 for which the dynamics of relaxation is controlled by the diffusion coefficients, expressed in the function of the second order moment of the transition rates. The expression of the second order moment for the group p is

where K(p, r) is the effective rate of energy relaxation from group p to group r,

where $fq=nqnN2$ denotes the fraction of particles in group q.

Figure 11(a) shows the second order moment for the energy-based bins formulation at 10 000K. Initially, the moment exhibits a monotonically increasing behavior with the bin energy. At later times, however, the efficiency of the first few groups increases creating a shallow bottleneck between 1 and 2 eV. At high energies, the dependence of the moment on the group energy is nearly exponential and appears to be unchanged during the relaxation.

Figure 11(b) shows the second order moment for the vibrational bins at 10 000 K. In general, the moments appear several orders of magnitude smaller across the entire energy spectrum. In particular, below 4 eV the coefficients are extremely small. This explains the formation of the bimodal distribution shown in Fig. 3(a) and in general slow vibrational relaxation observed in Fig. 10. It is important to stress that, contrary to what observed by other researchers, the vibrational relaxation does not exhibit a bottleneck.⁹⁹ This is due to the effect of rotation, that provides additional channels to vibrational processes, hence enhancing the relaxation.

Figure 12 shows the excitation rate contour plots, as a function of the energies of the products, for given pairs of initial groups. The analysis was repeated for the energy and vibrational grouping models. To facilitate the comparison between the models, we selected reactants with similar energy groups. In both cases, the energy of the second reactant was increased to analyze the behavior of the rates for the lower and the upper part of the distribution.

The contours for the energy-based grouping, shown in Figs. 12(a) and 12(c), exhibit a maximum corresponding to the “elastic” processes in which neither of the groups change. The magnitude of the rates decays exponentially for increasing energies of the products. The behavior is nearly isotropic, indicating that both reactants have equal probability of being excited. It is interesting to note that all the possible energy transfer reactions follow this behavior, and it is possible to fit the rates with a unique exponential function with three parameters (A, σ, γ), based on the work of Barker et al.,¹⁰⁰ 

In Figs. 12(b) and 12(d), the reaction rates for the vibrational energy groups exhibit a different behavior. At low energies, the vibrational-translational (VT) energy transfer processes are very inefficient compared with the high energy states. This is consistent with the results shown in Fig. 11(b). In contrast, due to the anharmonicity of the vibrational states, multi-quantum jumps are very probable for the high lying vibrational energy levels [Fig. 12(d)]. This justifies the establishment of a multimodal distribution in the early part of the relaxation. Moreover, at low energies, vibrational-vibrational (VV) energy transfer appears very efficient, thus favoring the thermalization of the distribution.

Sections III A– III C focused on the analysis of the dynamics of dissociation and energy transfer obtained with the two different grouping strategies. In order to assess the validity of each model, the results obtained are now compared against the available experimental data. Two different observables are used: phenomenological dissociation rate coefficient and an energy transfer relaxation time.

The thermal and QSS dissociation rate coefficients predicted by both the grouping models are compared against experimental data in Fig. 13. The numerical values of the rates in function of temperature can be found in Table II.

The thermal dissociation rate, shown in Fig. 13(a), is computed by assuming an equilibrium distribution between all the groups

Although in excellent agreement with the work of Jaffe et al.¹⁰¹ and Bender et al.,⁹⁷ both predictions grossly overestimate the experimental rates from Appleton et al.,⁹ and Kewley and Hornung.¹⁰² This is not surprising, since the population of the high-lying energy levels for both models was found to strongly deviate from equilibrium.

The QSS rates shown in Fig. 13(b) for both models are more consistent with the experimental data. In particular, the energy-based groups show excellent agreement with the Appleton⁹ data across the range of experimental conditions. The vibrational specific grouping predicts a lower rate coefficient at low temperatures (e.g., 10 000 K), and the slope of the rate does not agree with the experimental fit. At higher temperatures, both grouping strategies predict similar rate coefficients. The Kewley and Hornung data¹⁰² predict a significantly different temperature dependence, inconsistent with what was predicted by the CG-QCT method, demonstrated by the different slope in the data in Fig. 13. It is the authors’ opinion, based on preliminary analysis, that the discrepancy in the Kewley and Hornung data compared with the CG-QCT method is due to the non-equilibrium model adopted by the authors to interpret the experimental data (not shown). Additional experimental data, not included in the figure, have been used in the comparison. Hanson and Baganoff¹⁰³ predict the same temperature dependence as Kewley and Hornung, but their rates are significantly higher than all other experimental results. Park¹⁰⁴ reinterpreted the experimental data of Appleton, by including non-equilibrium effects in the post-processing. As a result, the Park dissociation rate coefficient appears larger than the original value given by Appleton; however, the temperature dependence is similar.

The vibrational relaxation time was computed using the e-folding method (the time required for the mode to reach 63.2% of its steady state energy) from the vibrational bins. The vibrational relaxation time is shown in Fig. 14, compared with the Millikan-White correlation,¹⁰⁵ and high-temperature correlations developed by Park⁶ and Boyd.¹⁰⁶ At low temperatures, the vibrational relaxation matches well with the Millikan-White correlation. The high-temperature corrections to Millikan-White match well at high temperatures with the vibrational binning data.

The results have clearly shown how the choice of grouping strategy has a profound impact on the characteristics of the thermochemical relaxation. The two grouping strategies adopted in this work are based on two fundamentally different assumptions: the vibrational grouping is constructed on the assumption of rigid separation between the rotational and vibrational energy modes, while the energy-based grouping, by lumping the levels independently of their vibrational and rotational characteristics, assumes exactly the opposite. Each of the governing assumptions has its merits: the separation of modes was clearly demonstrated experimentally and theoretically for the low lying energy levels. This justifies the use of a vibrational-based grouping strategy. On the contrary, the high energy states exhibit much stronger rovibrational coupling, which implies that energy-based binning provides a more accurate description of their behavior.

The energy-based grouping model predicts that molecules must climb to high energy rovibrational levels before dissociating and that a significant amount of energy lost in dissociation comes from rotationally excited molecules. By dissecting the energy ladder into narrow rovibrational energy groups, the contribution of vibration and rotation to dissociation can be separated without introducing significant artificial bias. The energy grouping model predicts nearly equal contribution of rotation and vibration to dissociation. This is not the case for the vibrational grouping method, which strongly underpredicts the role of rotation. Moreover, the narrow width of the energy bins (0.25 eV) limits the negative influence of the assumption of thermal equilibrium between translation and the group internal temperature.

In contrast, the vibrational binning strategy assumes equilibrium among rotational states across a wide range of energy, especially for the low vibrational energy groups (e.g., $v=0$ spans nearly 15 eV). As a result of the averaging process, the importance of dissociation from the high energy states is overwhelmed by the improbability of the low energy states dissociating. The weakness of the vibrational grouping strategy originates from the fact that low and high energy states are governed by different kinetics. To be successful, rovibrational states with similar rates should be grouped together (or averaged together).⁸⁶ 

The energy transfer process predicted by the energy level grouping was very fast across the range of temperatures. Because states with very different quantum configurations are lumped together in this approach, the mode separation known to exist for low energy states during the relaxation is not captured. This problem is not present in the vibrational-based grouping because the difference in the dynamics of slow and fast processes (i.e., vibrational and rotational relaxation, respectively) is correctly captured. Since internal energy relaxation and dissociation do not overlap in the thermochemical relaxation process, errors in the modeling of the relaxation are unlikely to affect the dynamics of dissociation.

The objective of this work is twofold: first to present a novel approach to constructing reduced order models (CG-QCT) for diatom-diatom interactions and second to study the non-equilibrium relaxation of nitrogen molecules due to N₂–N₂ kinetic processes. The CGM is constructed by lumping rovibrational states together into groups and by determining the group averaged kinetic properties by means of the quasi-classical trajectory calculations. Moments of the master equations are used to compute the population of the groups during the relaxation. Two very different grouping strategies, energy and vibrational-based grouping, were used to investigate the non-equilibrium relaxation of molecular nitrogen in zero-dimensional isothermal and isochoric reactors. It was found, based on comparison with experimental data, that both approaches captured the dissociation process,⁹ while the vibrational bins better captured the energy transfer process.^6,105,106

The comparison of the energy transfer relaxation rates, predicted by the vibrational grouping method, against available experimental data,^6,105,106 validates the classical picture of relaxation based on the rigid separation among the rotational and vibrational energy modes. It was found that the relaxation data provided by the Millikan and White correlation are in good agreement with the theoretical predictions provided by the vibrational grouping method for temperatures below 10 000 K. On the contrary, any attempt to ignore the mode separation by grouping rovibrational levels characterized by different vibrational quantum number has lead to an overestimation of the relaxation rate.

The dissociation process is controlled by the tail of the rovibrational distribution. For the high-lying states, the strong rovibrational coupling is well captured by the energy-based grouping. So molecules have to climb the rovibrational energy ladder before dissociating, thus acquiring a significant vibrational and rotational energy. In all the cases analyzed, rotational energy accounts for more than 40% of the total energy required to dissociate the N₂ molecules. While vibrational bias is introduced in the conventional non-equilibrium models, no models include bias to account for effect of rotational excitation of the molecules.⁹⁷ Finally, it is important to note that during the dissociation process, the distribution of the N₂ molecules is in QSS. This is very different from what is observed in the N₂–N case, where significant part of the dissociation takes place in non-QSS conditions.

To complete the analysis of these two grouping strategies, the second part of this paper presents a detailed comparison of the CG-QCT method with the Direct Molecular Simulation (DMS)^92–94 method.⁹⁵ 

The authors would like to thank AFOSR and in particular Dr. Ivett Leyva for supporting this research. Ms. R. Macdonald was supported by the Department of Defense (DoD) through the National Defense Science and Engineering Graduate (NDSEG) Fellowship Program. Dr. M. Panesi was supported by the Air Force Office of Scientific Research Young Investigators Program No. FA9550-15-1-0132 with Program Office Dr. Ivett Leyva. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the AFOSR or the U.S. government.

The reverse reaction rates are obtained by applying microreversibility. Microreversibility holds for the StS reactions

where the ⋆ denotes the equilibrium composition at the local translational temperature of each species. The internal temperature of the group is assumed to be in equilibrium with the translational mode, giving an equilibrium distribution at the local translational temperature. Rearranging this expression and using the definition of the grouped rate, $Kr−pqER$⁠,

Multiplying and dividing the entire expression by the equilibrium population of the groups p and q allows the expression to be re-written in terms of the forward group rate, $Kpq−rED$⁠,

The equilibrium composition is obtained using the partition functions

For the molecules in groups p, q, and r, the partition function can be written as the product of the internal and translational partition functions

where m_p denotes the mass of the molecule containing group p, and h_P denotes Planck’s constant. For the atoms A and B, assuming that they both are in the ground electronic state, which have degeneracy given by g_A and g_B, respectively, the partition function is written as

where m_A denotes the mass of atom A. The partition function of atom B is of the same form. This framework can also be applied to determine the recombination rates of the $KpqR$ reaction

Similar to the combined excitation-dissociation reaction, the partition functions are used to compute the equilibrium composition of groups p and q and atoms A and B.