Modeling of the electronic excited states in high-temperature flows Open Access

In recent decades, the direct simulation Monte Carlo (DSMC) method has become one of the standard methods for simulating rarefied non-equilibrium flows.^1,2 One reason for this is the relatively simple implementation of the method itself, along with the various physical and chemical models that can be implemented. Thus, with DSMC, it is possible to model very complex non-equilibrium effects in a flow, such as non-equilibrium in the internal degrees of freedom of molecules, i.e., rotation and vibrational energies, or non-equilibrium chemical processes, depending on the chosen physical models. It has long been known that excitation of the internal states strongly influences chemical processes, which in the case of dissociation and vibrational energy has already been modeled using the vibrationally favored dissociation (VDF) model.³ The link between the vibrational energy and the dissociation probability became even more pronounced with the introduction of the quantum kinetic (QK) chemistry models⁴—and its numerous extensions.^5–7 In contrast, the electronic excitation of chemical species is sometimes disregarded, as some applications occur within temperature ranges in which the electronic excitation can still be neglected.^8–10 

The assumption of neglecting the electronic mode depends on the chemical species involved in the flow. For instance, quasi-classical trajectory (QCT) calculations for the dissociation of oxygen conducted by Jo et al.¹¹ reveal that dissociation reactions not only occur from the electronic ground state but also from electronic excited states. These QCT calculations demonstrate the importance of introducing a correction factor to account for the contribution of electronic excited states to reproduce the most representative experimental measurements. Similarly, the low electronic excitation energies observed in iron have substantial implications for flow physics, even at low temperatures, such as during meteorite reentry.¹² For applications involving typical hypersonic reentry, e.g., Hayabusa,¹³ Stardust,^14,15 and Fire II,^16,17 ionization reactions and radiation are predominantly governed by the excitation of the electronic states.^13–21 In these applications, the additional degrees of freedom of the electronic excitation modify the thermal conductivity, hence influencing the Prandtl number of the flow and the specific heats of the chemical species. Likewise, in applications involving plasma, e.g., electric propulsion systems or laser–plasma interactions, electronic excitation should not be disregarded, as it plays a significant role in the transport properties, chemical composition, and radiation processes.^22–26 

In the DSMC community, four methods have been developed to model electronic excitation in flows. The first method has been introduced by Liechty and Lewis^5,27 and treats the electronic excitation in a similar manner to the vibrational energies. In this method, each chemical species is initialized with an electronic quantum number sampled from the Boltzmann distribution. During collisions, energy exchange with the translational degrees of freedom occurs through the Larsen–Borgnakke (LB) model.²⁸ The primary advantage of this technique lies in the storage of only one additional energy value per particle, while many aspects of vibrational modeling can be incorporated. However, a drawback of this method is that it exhibits significant noise, especially at low temperatures, requiring a large number of particles to compensate for the noise.

The second method, initially introduced by Bird,²⁹ improved by Carlson et al.,^30,31 and later expanded by Burt and Josyula,^32,33 partially addresses this problem by equipping each particle with a complete distribution function of electronic excitation instead of an electronic quantum number. However, this leads to a significant increase in storage requirements.

The third method developed by Li et al.²¹ treats each transition between two electronic states as a chemical process during collisions, where each electronic excited state is regarded as a distinct species. However, to maintain a reasonable level of complexity, assumptions are made on the total number of electronic excited states involved. The merits of this method are that the transition rates can be tuned to reproduce any experimental measurements or high-fidelity calculations; albeit, noise remains a problem.

The fourth method has been designed by Gallis and Harvey.^34,35 While the preceding methods rely on the assumption of an equilibrium distribution of the electronic state after a collision, this method computes the probability of an electronic excitation resulting from a collision with the aid of electronic excitation cross sections. To achieve a realistic description of the electronic excitation, accurate cross sections are required. However, the experimental measurement of accurate cross sections is extremely challenging and covers only a fraction of the electronic transitions of interest. Consequently, electronic excitation cross sections are typically either calculated from theoretical considerations or approximated from other known cross sections for similar interactions.

A common aspect of all methods is that each internal mode, i.e., rotational, vibrational, and electronic, is regarded separately. However, real-gas effects involve coupling between all internal modes, such as rotational–vibrational coupling and rotational–vibrational–electronic coupling.³⁶ This becomes particularly significant in the analysis of molecular radiation, where achieving the correct distribution of vibrational and electronic excitation energy is crucial.

In recent years, coupling of DSMC methods with radiation transport solvers, e.g., PARADE,³⁷ NEQAIR,³⁸ or Specair,³⁹ has emerged for such investigations.^{19,22,24,25,40,41} In contrast to the DSMC solvers, these radiation solvers adopt a fully coupled approach in which all types of interactions between the internal modes are incorporated. To perform such coupling, previous coupled radiation-DSMC simulations introduced additional assumptions, i.e., Boltzmann distribution^19,41 or quasi-steady-state (QSS),^{14,21,40,42,43} to transfer uncoupled information from DSMC into the coupled models of radiation solvers. With an ever-increasing number of missions involving spectral studies of flows during in-flight or experimental measurements, a detailed description of the physics of the flow is fundamental to allow for spectral comparison. However, a limitation arises due to the assumption of decoupling each internal mode of the chemical species in DSMC simulations, preventing the achievement of this level of detail.

The novelty of the present article lies in the development of a model for coupling the electronic and vibrational modes of molecular species in DSMC. The model involves the development of new numerical techniques for the initialization of particles in the domain, the redistribution of the internal energy after a collision, and the measurement of the corresponding internal temperature. Considering the challenge of performing experiments to measure chemical processes involving electronic excited states, the new model is verified against an extensive compilation of theoretical studies for the reproduction of thermophysical properties of molecular oxygen. Finally, the model is applied for a canonical nonreactive oxygen hypersonic flow past an infinite cylindrical body at an altitude of 85 km in Earth's atmosphere.

The rest of the paper is organized as follows: Section II describes the new model for the coupling of the electronic and vibrational modes. In Sec. III, the numerical techniques for the implementation of the new model in a DSMC solver are developed. Section IV provides a comprehensive verification of the model and demonstrates its application under typical hypersonic flow conditions. The final section provides conclusions and perspectives for future work.

This approximation is traditionally applied in DSMC simulations, i.e., each mode of a molecular system is regarded as insensitive to any other excitation. In the current section, the traditional approach for DSMC simulations is reviewed, and a new model involving the coupling between the vibrational and electronic modes is presented.

The traditional approach imposes that each mode of the molecule is individually treated assuming no interactions with one another. This approach takes into account the translational, rotational, vibrational, and electronic modes.

For the most dominant chemical species in Earth's atmosphere, the characteristic rotational temperatures are typically very low, e.g., $θr,O2=2.064$ K, $θr,N2=2.869$ K, and $θr,NO=2.470$ K.⁴⁹ Therefore, the classical limit $T≫θr$ is applicable in most DSMC simulations, allowing for continuous treatment of the rotational mode.

The maximum vibrational quantum number, v_max, allowed by Eqs. (5) and (6), is calculated as the vibrational quantum number immediately before the vibrational energy reaches the dissociation limit. In this approach, it is important to emphasize that each electronic excited state can only excite the same set of vibrational levels as the ground electronic state.

The evaluation of the thermodynamic properties of such a system becomes a matter of addressing each mode individually; hence, it is denominated the uncoupled approach in the current work. Note that each internal mode, i.e., rotational, vibrational, and electronic, is treated separately resulting in three one-dimensional distribution functions. It is important to emphasize that the Boltzmann distribution is limited to discrete energy levels. In the present article, the translational and rotational modes are treated continuously, therefore, their distribution functions are described by the Hinshelwood function.⁵⁵ 

Considering an electronic state i, the maximum vibrational quantum level permissible is determined as the vibrational quantum immediately before the dissociation energy of electronic state i, i.e., $εve(i)<εd(i)$⁠. For some electronic configurations, all the vibrational quantum levels lie below the corresponding dissociation energy. In this specific scenario, the maximum vibrational quantum level is determined as the last vibrational quantum level before the gradient $∂εve(i)∂j$ changes sign.

As illustrated by Eqs. (8) and (13), the two approaches yield distinct distribution functions for the vibrational and electronic modes. The uncoupled approach uses a single one-dimensional distribution function for each mode, whereas the coupled approach introduces a two-dimensional distribution function that accounts for the coupling between the vibrational and electronic modes. This fundamental difference between the two approaches leads to significant consequences. Specifically, in the new coupled approach, a unique partition function, mean internal energy, and specific heat capacity can be defined. From a numerical perspective, it imposes the definition of a unique mean internal temperature, i.e., T_ve. This contrasts with the traditional uncoupled approach which requires separate internal temperatures for each mode, i.e., T_v and T_e.

As highlighted in Sec. II, the coupled approach has a two-dimensional distribution function that requires sampling of both the electronic and the vibrational quantum numbers simultaneously. The coupled approach is implemented in dsmcFoam+ based on the procedure reported by Liechty and Lewis.⁵ Three procedures are revisited to accommodate the vibrational excitation of the electronic excited states. These functions are the equilibrium sampling, the post-collision sampling, and the measurement of the corresponding mean temperature.

Figure 1 illustrates the general procedure to assign the initial electronic and vibrational quantum numbers. The first step consists of searching for the electronic and vibrational quantum numbers, i.e., (i_s, j_s), for which Eq. (13) has a maximum. Then, a pair of electronic and vibrational quantum numbers, i.e., $(i⋆,j⋆)$⁠, uniformly distributed between 0 and $imax−1$ and 0 and $jmax−1$⁠, respectively, are independently chosen randomly. Finally, an acceptance–rejection scheme is used to select a pair of electronic and vibrational quantum numbers, i.e., $(i⋆,j⋆)$⁠, from the distribution, Eq. (17) that satisfies $f′>R(0,1)$⁠.

In Reaction 1, the energy exchange occurs from a change in the vibrational quantum number. In Reaction 2, the energy exchange involves an electronic transition, see Table I, which, in turn, also necessitates a modification of the vibrational quantum number.

Typically, in DSMC simulations, a relaxation probability of $Pv=0.02$ is allowed to result in vibrational excitation/de-excitation while a relaxation probability of $Pe=0.002$ is allowed for electronic energy exchange.²⁹ For the coupled approach to reproduce the two separate redistributions of internal energy in the uncoupled approach, two relaxation mechanisms, Reaction 1 and Reaction 2, are implemented in dsmcFoam+. Each mechanism operates with a relaxation probability, namely, P₁ for Reaction 1 and P₂ for Reaction 2. Note that the relaxation probability P₂ is a conditional probability on P₁ to be true, see Fig. 2. Therefore, to reproduce the two separate redistributions of internal energy, P₁ must be set to 0.02 and P₂ to 0.1. Note that these relaxation probabilities are user-defined; hence, they can be modified to reproduce any baseline relaxation model.

The energy exchange involving vibrational energy, i.e., Reaction 1, follows the quantum LB approach. A detailed explanation for the vibrational excitation/de-excitation during an inelastic collision under the assumption of an anharmonic oscillator model can be found in Civrais et al.⁵¹ For brevity, the present section focuses on the redistribution of internal energy involving an electronic and vibrational excitation/de-excitation.

If an inelastic collision is accepted for Reaction 2 energy exchange, the first step is to search the maximum allowed vibrational quantum number for each electronic excited state, i.e., (i_s, j_s), that satisfies $εc<εis,js$ and to determine a pair of quantum numbers for which the distribution function, Eq. (22), is maximum. Then, an acceptance–rejection procedure is performed to sample $(i⋆,j⋆)$ from the distribution, Eq. (22). Finally, a pair $(i⋆,j⋆)$ is accepted if Eq. (22) satisfies $h′>R$⁠; otherwise, the procedure is repeated until values for $(i⋆,j⋆)$ are obtained.

Finally, a point should be made about the calculation of the vibrational-electronic temperature. In the uncoupled approach, the vibrational and electronic modes, i.e., T_v and T_e, are treated separately leading to the one mean temperature for each mode. To maintain mathematical consistency with the two-dimensional distribution function, Eq. (13), a unique temperature, T_ve, must be defined for the vibrational-electronic mode. The difference between the two approaches will be assessed in Secs. IV A–IV C.

The vibrational-electronic temperature ultimately represents a mean excitation of the vibrational-electronic mode. Since the electronic mode admits a set of vibrational levels, the vibrational-electronic temperature is unable to reflect the vibrational excitation of each electronic excited state but, rather, a mean excitation of the electronic excited states and their corresponding vibrational quantum levels. Despite its minimal physical insights, the vibrational-electronic temperature serves the purpose of retaining the mean energy of a molecular system, making it valuable for verification exercise.

Furthermore, in the derivation, Eqs. (13)–(16), no assumption has been made regarding thermal equilibrium between the vibrational and electronic modes. Thus, these modes can exhibit thermal non-equilibrium with each other. These considerations contrast with the definition of T_ve employed in the vast majority of multi-temperature extended Navier–Stokes equations solvers based on Gnoffo;⁵⁶ albeit, there exist high-speed and high-temperature computational fluid dynamics (CFD) methods that do not impose thermal equilibrium between modes.^57–59 

In the uncoupled approach, where vibrational excitation is modeled with an infinite harmonic oscillator, the mean vibrational degrees of freedom can be computed using the equipartition theorem, and the vibrational temperature derived analytically from the mean vibrational quantum number. However, with an anharmonic oscillator model, the vibrational temperature cannot be derived analytically.⁵¹ Specifically, it must be determined by solving an implicit equation. Similarly, to account for the vibrational excitation of the electronic excited states, the vibrational-electronic temperature cannot be directly calculated using the equipartition theorem. A similar approach to that described by Civrais et al.⁵¹ is herein adopted.

This implicit equation, Eq. (23), is numerically resolved with a Newton iterative approach. It involves initially evaluating Eq. (15) at a guessed temperature and comparing it to the mean vibrational-electronic energy in a cell. The procedure is repeated until a user-defined tolerance factor is reached.

A series of adiabatic reactor simulations is conducted to verify the derivation and implementation of the coupled approach. A single cubic cell with edge length $1.88×10−4$ m filled with 1 × 10⁶ DSMC simulator particles and periodic boundaries is used for this purpose. A fixed time step size of $1×10−9$ s is adopted. The working gas is molecular oxygen with the electronic excited states summarized in Table I. Collision partners are selected according to the No Time Counter (NTC) model.¹ The inter-molecular collisions are processed with the variable hard sphere (VHS) model¹ with the properties at a reference temperature T_ref of 273 K, i.e., $m=48.12×10−27$ kg, $d=4.07×10−10$ m, and $ω=0.77$⁠. The probabilities for a particle to experience Reaction 1 or Reaction 2 are those presented in Sec. III unless stated otherwise. The simulations are performed in serial on a machine equipped with 12 Intel^® Core™ i7-10850H central processing unit (CPU) cores (base: 2.70 GHz, max: 5.1 GHz).

In the present article, the seven lowest electronic excited states are considered.^60–62 Note that additional electronic excited states are reported by Liu et al.;⁶⁰ however, the decision has been made to be consistent with the number of electronic excited states engaged in previous DSMC studies.⁶³ The spectroscopy constants are extracted from the National Institute for Standards and Technology (NIST) database⁴⁹ with the exception of the ground state, i.e., $X3Σg−$⁠, which is extracted from Civrais et al.⁷ The spectroscopic constants of the electronic excited states of O₂ are summarized in Table I. In addition, it is assumed that the electronic excited states conserve the collision properties of the ground state; hence, the mass, m, diameter, d and viscosity exponent, ω, of the electronic excited states are those of the ground state.

The first verification involves the measurement of the vibrational-electronic temperature. The adiabatic reactor is initialized in thermal equilibrium conditions for temperatures ranging between 2000 and 20 000 K. The vibrational-electronic temperature is monitored and sampled for 10⁵ iterations.

Table II shows that the measurement technique implemented for the calculation of the vibrational-electronic temperature is accurate across the temperature range considered. The numerical scatter of the measured temperature is inversely proportional to the square root of the sampling size;² hence, for a sampling size of about 10⁶, the numerical scatter is about 0.1%. Table II shows that the difference between the initialized and measured temperatures reduces as expected. This demonstrates the correct implementation of both the initialization algorithm and the new measurement technique.

Figure 3 demonstrates that the DSMC results achieve excellent reproduction of the theoretical predictions for both the excitation and de-excitation scenarios. Quantitatively, for the excitation scenario, the theoretical equilibrium temperature is $Teq,extheo=13 740$ K and the equilibrium temperature measured in DSMC is $Teq,exDSMC=13 741$ K. For the de-excitation scenario, the theoretical equilibrium temperature is $Teqtheo=10 521$ K, and the equilibrium temperature measured in DSMC is $TeqDSMC=10 518$ K. Both scenarios show good reproduction of the theoretical prediction that demonstrates the successful implementation of the post-collision sampling technique.

Figure 4 shows further verification of the model for the population of the electronic excited states and their corresponding vibrational quantum numbers. The adiabatic reactor is initialized in thermal equilibrium conditions, i.e., $Tt=Tr=Tve=10 000$ K, and the population of each quantum number is sampled for 10⁵ time steps to reduce the scatter. The theoretical and DSMC two-dimensional distribution of energy across all electronic excited states and their corresponding vibrational quantum levels are depicted in Figs. 4(a) and 4(b), respectively. Additionally, the population of the electronic excited states and vibrational quantum levels are, respectively, depicted in Figs. 4(c) and 4(d).

Figures 4(a) and 4(b) show that the DSMC achieves a good reproduction of the two-dimensional theoretical distribution for all quantum levels. Figure 4(c) demonstrates an excellent agreement between the population of the electronic excited states and statistical mechanics theory. The consistency between DSMC results and the theoretical calculations, Eq. (13), serves as a further verification of the implementation of the model in the DSMC code. Similarly, Fig. 4(d) also confirms the accuracy of the model for the population of the vibrational quantum numbers in all seven electronic excited states. Some scatter is evident at the tails of the distributions. This is expected in a DSMC simulation, because the probability of finding a molecule in these higher vibrational levels is relatively small, resulting in a low signal-to-noise ratio.

The uncoupled and coupled approaches are examined for the calculation of the isochoric-specific heat capacity. The adiabatic reactor is initialized in thermal equilibrium for temperatures ranging from 2000 to 20 000 K, and the thermophysical properties are sampled for 10⁵ time steps to reduce the scatter. These results are compared against four numerical studies, i.e., Capitelli et al.,⁶⁴ Jaffe,⁶⁵ Qin et al.,⁶⁶ and McBride et al.⁶⁷ This represents a selection of many works available in the literature. Additionally, the specific heat capacity results presented in Civrais et al.⁵³ are included for reference.

Figure 5 compares the uncoupled and coupled approaches for the calculation of the isochoric specific heat capacity of molecular oxygen. The DSMC results are in excellent agreement with the theoretical calculations, Eq. (16), for the entire range of temperatures considered which serves as a verification of the implementation of the model in the DSMC code. Figure 5 illustrates the benefits of the coupled approach. This approach provides a more accurate description of the specific heat capacity compared to its counterparts. Figure 5 emphasizes the importance of including the electronic mode of molecular systems and modeling their corresponding vibrational excitation to yield accurate predictions in low to moderate temperatures.

For temperatures in excess of 7000 K, the coupled approach deviates from the prediction of past studies.^64–67 In the derivation of the coupled approach, the rotational–vibrational coupling effects of the electronic excited modes have been disregarded. While this assumption holds for low-to-moderate temperatures, it becomes evident that as the temperature increases the coupled approach deviates from the compilation of past studies due to the omission of rotational–vibrational coupling effects. To improve the accuracy of the coupled approach, further investigation should be directed toward incorporating the contribution of rotational-vibrational coupling effects.^68–71 

The present section aims to compare the uncoupled and coupled approaches for the simulation of a canonical nonreactive oxygen hypersonic flow past a cylindrical body entering Earth's atmosphere at an altitude of about 85 km. The free stream conditions are extracted from the U.S. Standard Atmosphere⁷² and summarized in Table III.

The geometry represents a two-dimensional slice of a cylindrical body with diameter d = 0.1 m, with a domain length equal to two and a half radii upstream of the stagnation point. The mesh is refined near the stagnation point to ensure that the cell size, $Δx$⁠, remains around one-quarter of the local mean free path, λ, throughout, resulting in a total of 32 400 cells. The time step is carefully chosen to be an order of magnitude smaller than both the mean collision time and the cell residence time, the latter of which corresponds to the time for a DSMC particle to cross a cell length under the freestream conditions. To maintain a minimal number of 100 particles per cell throughout, the numerical mesh is populated with a total number of $8.59×106$ DSMC particles at steady state. Similar to the adiabatic reactor simulations, the inter-molecular collisions are computed with the VHS model¹ for a reference temperature of T_ref = 273 K. Collision partners are selected according to the NTC model.⁷³ To reproduce the typical relaxation behavior of the uncoupled approach, a constant probability of $P1=0.02$ and $P2=0.1$ for a particle to undergo Reaction 1 and Reaction 2, respectively, is applied. The gas–surface interactions are fully diffusive with an isothermal wall at temperature T_w = 1000 K. The simulations are performed in parallel on 8 cores on the aforementioned machine and over 10⁵ samples are taken after steady-state to reduce the numerical scatter. To illustrate the structure of the flow, the internal temperatures and velocity magnitude along the stagnation streamline are presented in Fig. 6. The distance to the stagnation point is normalized by the diameter of the cylinder. A value of 0 refers to the inlet boundary conditions and 2.5 refers to the stagnation point.

As mentioned previously, the standard approach to couple DSMC and radiation transport solvers consists of transferring the uncoupled internal temperatures, i.e., rotational, vibrational, and electronic, calculated by the DSMC solver into a radiation solver. Note that these internal temperatures are calculated by assuming a Boltzmann distribution and resolving an implicit equation for both vibrational and electronic temperatures.⁵¹ However, radiation solvers, e.g., NEQAIR,³⁸ PARADE/PICLas,^19,37 and Specair,³⁹ integrate a fully coupled approach. Specifically, these solvers compute the radiative properties in a cell from the distribution function of the internal quantum levels. To compute these radiative properties, an assumption on the distribution function is made where it is assumed to follow Boltzmann statistics^19,41 or the non-equilibrium quasi-steady-state method.^21,40

The main objective of this section is twofold: to compare the two approaches for the distribution of the internal quantum levels measured by the DSMC method and to assess the assumption that the internal quantum levels are distributed according to the Boltzmann distribution. For both approaches, the vibrational and electronic density functions are calculated with the local internal temperatures, namely, $f=fi(Tv)fj(Te)$ for the uncoupled approach and $f=fi,j(Tve)$ for the coupled approach.

Figure 7 compares the DSMC results for the population of the internal quantum levels and the corresponding Boltzmann distribution along the stagnation streamline. The DSMC results are denoted by color markers, where solid symbols refer to the uncoupled approach and open symbols to the coupled approach, and the Boltzmann distributions are denoted by color lines where solid lines refer to the uncoupled approach and dashed lines refer to the coupled approach. Similar to Fig. 6, the distance to the stagnation point is normalized by the diameter of the cylinder. A value of 0 refers to the inlet boundary conditions and 2.5 refers to the stagnation point.

Downstream of the shock wave, Fig. 7 shows that the highest vibrational and electronic quantum levels exhibit large deviations from the assumption of a Boltzmann distribution. These high-lying quantum levels are of primary importance for the modeling of the radiative properties of molecular species. Furthermore, a noticeable discrepancy arises between assuming a Boltzmann distribution corresponding to the cell temperature and density, and the actual distribution observed in DSMC simulations, which represents a source of error in typical flow-radiation coupling.

Figure 7 also indicates a significant population of the high-lying quantum levels in the upstream region. This non-equilibrium effect manifests itself in Knudsen number regimes in which the degree of rarefaction is sufficient for a particle to collide with the surface of the body and travel backwards without undergoing sufficient collisions with any other particles to de-excite the quantum levels to return to an equilibrium distribution; hence, carrying post-shock information upstream of the bow-shock. This non-equilibrium effect is characterized by a strong deviation from a Maxwell velocity distribution to the extent of quasi-Maxwellian or bi-modal distributions.²¹ Note that this phenomenon has also been observed previously.^12,18,20,41

The deviation from the Boltzmann distribution is further investigated at four locations throughout the shock wave denoted by vertical red lines in Fig. 6 with thermal conditions summarized in Table IV. A series of adiabatic reactor simulations initialized with the thermal conditions summarized in Table IV is performed. Figure 8 shows the Boltzmann distribution and DSMC measurements of the coupled approach for the first three electronic states of O₂ obtained along the stagnation line. Contrastingly, Fig. 9 shows the Boltzmann distribution and DSMC measurements obtained from the series of adiabatic reactor simulations. The DSMC results are denoted by color markers and the Boltzmann distributions by color lines.

Figure 8 shows that the Boltzmann distributions at all four locations are significantly different to the DSMC measurements along the stagnation streamline, whereas Fig. 9 indicates that the DSMC simulations precisely reproduce the theoretical vibrational density function for the same thermal conditions in an adiabatic reactor. This demonstrates that the different distributions in the hypersonic case are due to non-equilibrium effects that a radiation solver would not be capturing using solely the temperatures.

So far, the comparison of the two approaches was limited to the population of the quantum levels along the stagnation streamline. In Fig. 10, the total number density of the first six electronic states of O₂, i.e., $ni=∑jni,j$⁠, measured in the DSMC simulations are regarded for the flow field.

As anticipated, Fig. 10 indicates that the ground electronic state is the most populated. Both approaches offer a similar description of the distribution of the ground state throughout the domain. However, the two approaches differ in the distribution for all the electronic excited states herein considered. Specifically, the coupled approach predicts lower populations of the electronic excited states in comparison to the uncoupled approach. This discrepancy arises from the difference in the modeling of the vibrational energies of the electronic excited states in each approach. The coupled approach implies that each electronic state can excite a unique set of vibrational levels, see Fig. 4. In contrast, the uncoupled approach assumes that each electronic state can only excite the same set of vibrational energies as the ground electronic state. In the context of a flow-radiation coupling, such a difference between the two approaches in terms of the distribution of the electronic excited states throughout the domain could lead to a significant difference in the radiation properties.

While the coupled approach has demonstrated good reproduction of theoretical calculations, see Figs. 3–5, for hypersonic flow past a cylindrical body, see Figs. 6–10, the approach has limitations that merit additional discussions.

The first limitation concerns the omission of the rotational–vibrational coupling. Figure 5 has evidenced the importance of the rotational–vibrational coupling effects for the reproduction of the most representative database.^64–71 There exist various approaches to leverage this limitation. One approach for modeling the rotational–vibrational coupling consists of a continuous treatment of the rotational mode and a quantized treatment of the vibrational and electronic modes.⁷⁴ This approach models the elongation of the bond length due to vibrational excitation in the rotational partition function. Specifically, it assumes that each vibrational quantum level has its own rotational and dissociation energies. The benefits of this approach lie in its capacity to provide a rotational–vibrational coupling without increasing the number of quantum levels considered. Another approach assumes a quantized treatment of the internal modes, i.e., rotational, vibrational, and electronic modes.⁷⁵ This approach models triplets of quantum levels that lie in the bounded and quasi-bounded areas.^64,65 While this approach offers a detailed representation of the dynamics of the internal quantum levels, it requires the addition of numerous quantum levels to the molecular systems, resulting in a significant increase in computational expense. For a typical air mixture simulation modeled with the five most dominant species, i.e., O₂, N₂, NO, O, and N, this approach requires ∼100 000 internal quantum numbers.

The second limitation concerns the chemical activity of the electronic excited states. In the current work, the chemical activity of the electronic excited states has consciously been disregarded. The inherent lack of experimental measurements of chemical processes involving electronic excited states of molecular species for temperatures relevant for hypersonic flow conditions does not allow the development of more sophisticated models. Such an issue would ideally be leveraged by QCT calculations on non-adiabatic potential energy surfaces, incorporating electronic excited state transitions. However, state-of-the-art QCT calculations of relevant species for Earth reentry are, as yet, limited to ground electronic configurations.¹¹ Therefore, non-adiabatic QCT calculations are required to allow for model development.

A novel model for modeling the vibrational excitation of the electronic excited states is developed and implemented in a DSMC solver. The model is verified against an extensive compilation of theoretical studies for both thermal equilibrium and non-equilibrium conditions. Additionally, it is evaluated against the traditional uncoupled approach utilized in DSMC for a canonical hypersonic non-reactive oxygen gas flow past a cylindrical body at an altitude of about 85 km.

The series of adiabatic reactor simulations have shown excellent agreement between theoretical predictions and DSMC results, thereby demonstrating the successful implementation in a DSMC solver. Furthermore, the coupled approach has been shown to improve the reproduction of the specific heat capacity of molecular oxygen compared to the uncoupled approach.

The hypersonic flow scenario has revealed that the assumption of the quantum levels being distributed according to the Boltzmann distribution is invalid. In high Knudsen number regimes, non-equilibrium effects lead to significant deviations between the local distribution of quantum levels and that predicted by the Boltzmann distribution calculated from internal temperatures. Consequently, one should benefit from the capability of DSMC solvers to keep track of the local distribution of energy across the internal quantum levels and solely rely on these distributions for the determination of the radiative properties. Furthermore, the coupled approach has demonstrated a lower population of electronic excited states compared to the uncoupled approach. This discrepancy is attributed to the possibility for each electronic excited state to allow vibrational excitation, hence providing additional channels for redistributing internal energy during collisions. Such a difference in the distribution of electronic excited states between the two approaches is expected to result in modifications in the radiation properties.

The limitations of the model have been discussed with the results that both rotational–vibrational–electronic coupling and chemical reactions involving electronic excited states have been disregarded. Further investigations will focus on incorporating rotational–vibrational–electronic coupling into the mathematical description of the model and ultimately in a DSMC solver. Additionally, computational chemistry studies¹¹ have evidenced that adding the contribution of electronic excited states to the global dissociation rates significantly improves the reproduction of experimental measurements. Therefore, future research will involve the modeling of the chemical activity of the electronic excited states in the context of hypersonic flows.

This work was initiated during a secondment at the Institute of Space Systems at the University of Stuttgart. The lead author gratefully acknowledges scholarship funding from the James Watt School of Engineering through the School of Engineering Scholarship 2020 and the Jim Gatheral Travel Grant 2023.

The authors have no conflicts to disclose.

C. H. B. Civrais: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Software (lead); Supervision (equal); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). M. Pfeiffer: Conceptualization (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). C. White: Supervision (supporting); Writing – review & editing (supporting). R. Steijl: Supervision (equal); Writing – review & editing (supporting).

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