Vibrational relaxation and triggering of the non-equilibrium vibrational decomposition of CO₂ in gas discharges

Plasma chemical conversion of CO₂ into CO is intensively studied now as a potential upstream process for production of synthetic fuels and chemicals (power-2-X) [1–3]. Initially those studies were inspired by, presumably, very high energy efficiency of the CO₂ conversion in microwave plasmas reported in the past [4, 5]. The chemical energy efficiency η_chem up to 80% achieved at relatively low translational–rotational temperatures T ≲ 1000 K was ascribed to dissociation from high vibrational states in conditions of strong vibrational non-equilibrium [4–6]. The reported high efficiencies at low T were not reproduced later. In the modern day experiments the highest η_chem of 40%...50% are observed at T > 2000 K and are largely explained by thermal quenching [7–10]. Nevertheless, activating the non-equilibrium conversion mechanism remains an attractive option for increasing the efficiency of plasma processes to make them in this respect competitive with electrolysis.

The non-equilibrium vibrational conversion can be triggered when the kinetic energy of electrons is deposited predominantly into excitation of vibrational states. The electron kinetics calculations had undefined shown [11–13] that this shall be indeed the case for reduced electric fields typical for the microwave plasmas. This condition is necessary, but not sufficient. In order for the dissociation from high vibrational states to be eventually activated the speed of relaxation of vibrational energy with the translational–rotational degrees of freedom must be slower than the speed of the chemical process. Vibrational relaxation in CO₂ and its connection with activation of the non-equilibrium CO₂ dissociation is the subject of the present paper.

Since the non-equilibrium mechanism of dissociation has not been undoubtedly detected in plasma experiments the theoretical and computational analysis is indispensable for understanding how its activation could be achieved. A simple consideration based on the semi-empiric vibrational relaxation equation known from shock tube and sound absorption experiments [14, 15] suggests that the activation of the non-equilibrium dissociation is controlled by the governing parameter $Q/{n}_{0}^{2}$, where Q is the specific (volumetric) power input into vibrational states, and n₀ is the initial number density of molecules. The non-equilibrium process can be triggered when the parameter $Q/{n}_{0}^{2}$ exceeds some critical value. To obtain quantitative estimates and verify the result of the semi-empiric theory this latter is supplemented by a computational study performed with a 0D state-to-state vibrational kinetics model.

There is no generally accepted 'standard' way of modelling the vibrational kinetics of CO₂. Unlike for diatomic molecules one cannot describe each single vibrational state with an individual particle balance (master) equation because the number of states is prohibitively large. Therefore, all existing state-to-state CO₂ models use one or the other kind of 'coarse graining' [16–19]. An extensive account on the subject can be found in [20]. For the present study the 'two-modes model' [21] is applied where the treatment of vibrational levels is simplified by combining symmetric stretching and bending vibrations into one effective mode. This is an approximate model with restricted quantitative accuracy which main purpose is verification and extension of the semi-empiric approach. The modelling study is divided into three steps, or into consideration of three model systems:

• (a)  
  A test which mimics the conditions behind a shock front. Its purpose is to check that the vibrational kinetics model is properly calibrated to reproduce the empiric relaxation times.
• (b)  
  Thermal dissociation of CO₂. The results are compared with shock tube measurements of the rate of the process CO₂ + M → CO + O + M for M = Ar, and with literature data for M = CO₂. This benchmark indicates that at least for T > 2000 K the model is able to reproduce the correct orders of magnitude of the dissociation rate.
• (c)  
  The model calibrated, verified and partly validated in the two previous steps is applied to study the onset of the low T non-equilibrium vibrational dissociation in conditions relevant for microwave plasmas.

On the step (c) in a series of numerical experiments it is shown that the rate of the primary dissociation process CO₂ → CO + O rapidly increases starting from zero when $Q/{n}_{0}^{2}$ is increased above some critical value. The critical value of $Q/{n}_{0}^{2}$ as well as that of the critical vibrational temperature are found to be not sensitive with respect to model uncertainties. The ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ calculated by the semi-empiric formula agrees well with the results of numerical simulations.

The rest of the paper is organized as follows. The next section summarizes the analytic and semi-empiric treatment of vibrational relaxation in CO₂, introduces the governing parameter $Q/{n}_{0}^{2}$ and the criterion of activating the non-equilibrium dissociation. A brief overview of the two-modes model, underlying assumptions and constrains is given in section 3. Sections 4–6 describe the numerical experiments with the model systems listed above, including the comparison of the two-modes simulations with experimental data and analytic formulas. Section 7 puts the theoretical results presented in the paper into the context of the current CO₂ plasma conversion research. Last section gives a summary of the main findings and a brief outlook on the prospective of future work.

2.1. Summary on vibrational relaxation equations

The basic information about vibrational relaxation in molecular gases is obtained in shock tube and acoustic experiments [14, 15, 22]. Here for compactness the consideration is structured around relaxation in shock waves. Same relaxation processes take place in sound waves and the relaxation rates deduced from both types of experiments agree [15, 22]. Behind shock wave fronts a very rapid increase of translational temperature is followed by its fast equilibration with rotational temperature. The characteristic time of translational–rotational relaxation is $\sim 1{0}^{-9}$ sec at atmospheric pressure (see [14], sections 48 and 49). The thermal excitation of vibrational states is much slower: the typical characteristic times at 1 bar are 10⁻⁶...10⁻⁵ s. At the same time, unless the temperatures are not very high, the change of the gas composition due to chemical reactions takes place on a significantly longer time scale.

It is well known that for a gas of diatomic molecules with unchanging number density the relaxation of vibrational energy after instant increase of translational–rotational temperature T can be described by a simple differential equation, see e.g. [14], chapter 19:

$$\begin{equation}\frac{\mathrm{d}{E}_{\text{vibr}}}{\mathrm{d}t}=\frac{1}{\tau }\left[{E}_{\text{vibr}}^{\mathrm{e}\mathrm{q}}\left(T\right)-{E}_{\text{vibr}}\right]\end{equation} \tag{ 1 }$$

Here E_vibr is the average specific vibrational energy per molecule, ${E}_{\text{vibr}}^{\mathrm{e}\mathrm{q}}\left(T\right)$ is the average vibrational energy per molecule at Boltzmann equilibrium with temperature T, τ is the relaxation time which depends only on T and on the number density of collision partners M (see below) n_M :

$$\begin{equation}{\tau }^{-1}={R}_{10}\left(T\right){n}_{\mathrm{M}}\left(1-{\mathrm{e}}^{-\frac{\hslash \omega }{T}}\right)\end{equation} \tag{ 2 }$$

where R₁₀ is the rate coefficient of transition from the first vibrationally excited state to vibrational ground state, ℏω is the energy of one vibrational quantum. In this simplest case of one-component gas n_M equals to the number density n of the molecules themselves. Here and below, unless specified explicitly, the temperatures are always expressed in energy units.

Equations (1) and (2) are derived on the following assumptions: (i) the molecules are linear (harmonic) oscillators; (ii) the probabilities of vibrational–translational (VT) transitions can be calculated by a first order perturbation theory: e.g. SSH (Schwartz–Slawski–Herzfeld) or Landau–Teller theory; (iii) the number of vibrational states is infinite. For CO₂ which has three modes of oscillations (1) and (2) have to be modified. Two more assumptions are added: (iv) exchange of vibrational energy with translational–rotational modes proceeds only via transitions between the bending mode oscillations:

$$\begin{equation}{\mathrm{C}\mathrm{O}}_{2}\left({v}_{1}{v}_{2}^{l}{v}_{3}\right)+\mathrm{M}\to {\mathrm{C}\mathrm{O}}_{2}\left({v}_{1},{v}_{2}^{l\pm 1}\pm 1,{v}_{3}\right)+\mathrm{M}\end{equation} \tag{ 3 }$$

(v) exchange of energy between bending modes ${v}_{2}^{l}$ and stretching modes v₁, v₃ is very fast. Here Herzberg notation of vibrational states is used, see [23]. Rigorous derivation under assumptions above can be found in appendix A. The resulting equation is same as (1), but with E_vibr on the right-hand side replaced by E_vibr-2—the specific vibrational energy per molecule accumulated in the bending modes only:

$$\begin{equation}\frac{\mathrm{d}{E}_{\text{vibr}}}{\mathrm{d}t}={R}_{10}\left(T\right){n}_{\mathrm{M}}\left(1-{\mathrm{e}}^{-\frac{\hslash {\omega }_{2}}{T}}\right)\left[{E}_{\text{vibr}-2}^{\mathrm{e}\mathrm{q}}\left(T\right)-{E}_{\text{vibr}-2}\right]\end{equation} \tag{ 4 }$$

$$\begin{equation}{E}_{\text{vibr}-2}^{\mathrm{e}\mathrm{q}}\left(T\right)=\frac{2\hslash {\omega }_{2}}{{\mathrm{e}}^{\frac{\hslash {\omega }_{2}}{T}}-1}\end{equation} \tag{ 5 }$$

here ℏω₂ is the fundamental energy of bending oscillations v₂. ${E}_{\text{vibr}-2}^{\mathrm{e}\mathrm{q}}\left(T\right)$ is calculated by applying the known formula for harmonic oscillators, see e.g. [24], equation (49.3) there. Factor 2 takes into account that the mode v₂ of CO₂ is double degenerate. R₁₀ is the rate coefficient of the process:

$$\begin{equation}{\mathrm{C}\mathrm{O}}_{2}\left(0{1}^{1}0\right)+\mathrm{M}\to {\mathrm{C}\mathrm{O}}_{2}\left(0{0}^{0}0\right)+\mathrm{M}\end{equation} \tag{ 6 }$$

In [25] it had been shown that when a system of linear oscillators initially has Boltzmann distribution with temperature equal to the translational–rotational temperature T, and T is instantly increased, then vibrational relaxation towards final equilibrium state proceeds via a sequence of Boltzmann distributions with vibrational temperature T_vibr ⩽ T. The result of [25] is not directly applicable to CO₂, but numerical experiments presented in section 4 below confirm its validity in that case as well. For the Boltzmann vibrational distribution equation (4) is transformed into equation for T_vibr which can be easily integrated:

$$\begin{align}\hfill \frac{\mathrm{d}{T}_{\text{vibr}}}{\mathrm{d}t}& =\frac{{R}_{10}\left(T\right){n}_{\mathrm{M}}\left(1-{\mathrm{e}}^{-\frac{\hslash {\omega }_{2}}{T}}\right)}{{c}_{\text{vibr}}\left({T}_{\text{vibr}}\right)}\hfill \\ \hfill & \quad \times \left[{E}_{\text{vibr}-2}^{\mathrm{e}\mathrm{q}}\left(T\right)-{E}_{\text{vibr}-2}^{\mathrm{e}\mathrm{q}}\left({T}_{\text{vibr}}\right)\right]\hfill \end{align} \tag{ 7 }$$

Here ${c}_{\text{vibr}}\left({T}_{\text{vibr}}\right)=\frac{\mathrm{d}{E}_{\text{vibr}}}{\mathrm{d}{T}_{\text{vibr}}}$ is the (dimensionless) heat capacity due to vibrational modes which is calculated by applying the formulas for linear oscillators, see e.g. [24], equation (49.4):

$$\begin{equation}{c}_{\text{vibr}}\left({T}_{\text{vibr}}\right)={c}_{1}+2{c}_{2}+{c}_{3},\hspace{25.0pt}{c}_{i}=\frac{{\left(\hslash {\omega }_{i}\right)}^{2}{\mathrm{e}}^{\frac{\hslash {\omega }_{i}}{{T}_{\text{vibr}}}}}{{T}_{\text{vibr}}^{2}{\left({\mathrm{e}}^{\frac{\hslash {\omega }_{i}}{{T}_{\text{vibr}}}}-1\right)}^{2}}\end{equation} \tag{ 8 }$$

ℏω_i are the fundamental energies of the three modes of CO₂ oscillations.

The assumption of Boltzmann distribution of vibrational levels allows to derive one more equation where the rate of vibrational relaxation appears explicitly—by transforming (4) into equation for E_vibr-2:

$$\begin{align*}\hfill \frac{\mathrm{d}{E}_{\text{vibr}}}{\mathrm{d}t}& =\frac{\mathrm{d}{E}_{\text{vibr}}}{\mathrm{d}{E}_{\text{vibr}-2}}\frac{\mathrm{d}{E}_{\text{vibr}-2}}{\mathrm{d}t}=\frac{\frac{\mathrm{d}{E}_{\text{vibr}}}{\mathrm{d}{T}_{\text{vibr}}}}{\frac{\mathrm{d}{E}_{\text{vibr}-2}}{\mathrm{d}{T}_{\text{vibr}}}}\frac{\mathrm{d}{E}_{\text{vibr}-2}}{\mathrm{d}t}\hfill \\ \hfill & =\frac{{c}_{\text{vibr}}\left({T}_{\text{vibr}}\right)}{2{c}_{2}\left({T}_{\text{vibr}}\right)}\frac{\mathrm{d}{E}_{\text{vibr}-2}}{\mathrm{d}t}\hfill \end{align*}$$

$$\begin{align}\hfill \frac{\mathrm{d}{E}_{\text{vibr}-2}}{\mathrm{d}t}& =\frac{2{c}_{2}\left({T}_{\text{vibr}}\right)}{{c}_{\text{vibr}}\left({T}_{\text{vibr}}\right)}{R}_{10}\left(T\right){n}_{\mathrm{M}}\left(1-{\mathrm{e}}^{-\frac{\hslash {\omega }_{2}}{T}}\right)\hfill \\ \hfill & \quad \times \left[{E}_{\text{vibr}-2}^{\mathrm{e}\mathrm{q}}\left(T\right)-{E}_{\text{vibr}-2}\right]\hfill \end{align} \tag{ 9 }$$

One can see that the rate of relaxation (subsequently, the instant relaxation time) is not independent of T_vibr, but this dependency is weak. The factor $\frac{2{c}_{2}}{{c}_{\text{vibr}}}$ only changes from 1 at low T_vibr where c₂ ≫ c₁, c₃ (because ℏω₂ < ℏω₁, ℏω₃) to $\frac{2{c}_{2}}{{c}_{\text{vibr}}}=\frac{1}{2}$ at high T_vibr where c₁, c₂, c₃ → 1.

Equation (9) suggests the following approximation. One can assume that the total vibrational energy E_vibr decays with exactly same rate as E_vibr-2, and then apply equation (1) with τ defined as:

$$\begin{equation}{\tau }^{-1}=\frac{2{c}_{2}\left(\bar{T}\right)}{{c}_{\text{vibr}}\left(\bar{T}\right)}{R}_{10}\left(T\right){n}_{\mathrm{M}}\left(1-{\mathrm{e}}^{-\frac{\hslash {\omega }_{2}}{T}}\right)\end{equation} \tag{ 10 }$$

with $\bar{T}$ being some average temperature, e.g. $\bar{T}=\frac{1}{2}\left({T}_{0}+T\right)$ where T₀ is the initial Boltzmann temperature before the vibrational relaxation starts. Of note is that a relation similar to (10) appears also in [22], equation (24) there.

To estimate the accuracy of this approximation a comparison is made with numerical solution of equation (7) with constant T. Both solutions are expressed in terms of the dimensionless time:

$$\begin{equation}{t}^{\prime }=t\cdot {R}_{10}\left(T\right){n}_{\mathrm{M}}\left(1-{\mathrm{e}}^{-\frac{\hslash {\omega }_{2}}{T}}\right)\end{equation} \tag{ 11 }$$

That is, the explicit dependence on n_M and R₁₀ is eliminated. The temperature T_vibr obtained by solving (7) is translated into E_vibr by the linear oscillator formula:

$$\begin{equation}{E}_{\text{vibr}}\left({T}_{\text{vibr}}\right)=\frac{\hslash {\omega }_{1}}{{\mathrm{e}}^{\frac{\hslash {\omega }_{1}}{{T}_{\text{vibr}}}}-1}+\frac{2\hslash {\omega }_{2}}{{\mathrm{e}}^{\frac{\hslash {\omega }_{2}}{{T}_{\text{vibr}}}}-1}+\frac{\hslash {\omega }_{3}}{{\mathrm{e}}^{\frac{\hslash {\omega }_{3}}{{T}_{\text{vibr}}}}-1}\end{equation} \tag{ 12 }$$

The fundamental energies of the CO₂ oscillations [26]: ℏω₁ = 0.167 83 eV, ℏω₂ = 0.083 427 eV, ℏω₃ = 0.297 10 eV. The time evolution ${E}_{\text{vibr}}\left({t}^{\prime }\right)$ is compared with the exponential function:

$$\begin{equation}{E}_{\text{vibr}}^{\tau }\left(t\right)={E}_{\text{vibr}}\left(T\right)+\left[{E}_{\text{vibr}}\left({T}_{0}\right)-{E}_{\text{vibr}}\left(T\right)\right]{\mathrm{e}}^{-t/\tau }\end{equation} \tag{ 13 }$$

where t is replaced by t' and τ is replaced by ${\tau }^{\prime }=\frac{{c}_{\text{vibr}}\left(\bar{T}\right)}{2{c}_{2}\left(\bar{T}\right)}$, see (10) and (11). The relative difference between two solutions is defined as:

$$\begin{equation}{\epsilon}=\underset{{t}^{\prime }}{\mathrm{max}}\left[\frac{\left\vert {E}_{\text{vibr}}\left({t}^{\prime }\right)-{E}_{\text{vibr}}^{\tau }\left({t}^{\prime }\right)\right\vert }{{E}_{\text{vibr}}\left({t}^{\prime }\right)}\right]\end{equation} \tag{ 14 }$$

where the maximum is taken over the time interval from t' = 0 to t' = 10 ⋅ τ'.

For the test with T₀ = 300 K and $\bar{T}=\frac{1}{2}\left({T}_{0}+T\right)$ the quantity is found to be always ⩽5% for T ⩽ 2500 K. This outcome means that, provided that all the assumptions underlying (7) are fulfilled, the time evolution of E_vibr can be approximated with good accuracy by equation (1) with τ defined by (10). This characteristic time, in turn, depends on the rate coefficient of only one process (6). The coefficient R₁₀ then serves as a single parameter which fits the vibrational kinetics model into experimental data.

2.2. Threshold of the non-equilibrium vibrational dissociation

There is a direct link between the rate of vibrational relaxation and the activation of the non-equilibrium vibrational dissociation in a gas discharge. This link can be found by the following simple consideration. In the most general terms the balance of vibrational energy of the CO₂ molecules is written as:

$$\begin{equation}\frac{d\left(n{E}_{\text{vibr}}\right)}{dt}=Q-{Q}_{\mathrm{V}\mathrm{T}}-{Q}_{\text{diss}}\end{equation} \tag{ 15 }$$

where n is the number density of CO₂ molecules, Q is the power deposited into their vibrational modes, Q_VT is the rate of the energy losses into translational–rotational modes, and Q_diss is the rate of vibrational energy losses into chemical transformations. For the dissociation from the upper vibrational states to be significant their population has to be high enough, that is, T_vibr—or E_vibr—has to increase up to a certain high value. For this to happen at low translational–rotational temperature T the right-hand side of (15) has to be larger than zero (n ≈ const. before the dissociation starts).

In the near threshold regimes where Q_diss < Q_VT the last term can be neglected. The term Q_VT can be written in the form suggested by equation (1):

$$\begin{equation}{Q}_{\mathrm{V}\mathrm{T}}\left(T,{E}_{\text{vibr}}\right)={R}_{\mathrm{V}\mathrm{T}}\left(T\right){n}_{\mathrm{M}}n\left[{E}_{\text{vibr}}-{E}_{\text{vibr}}^{\mathrm{e}\mathrm{q}}\left(T\right)\right]\end{equation} \tag{ 16 }$$

where ${R}_{\mathrm{V}\mathrm{T}}\left(T\right)={\left(\tau {n}_{\mathrm{M}}\right)}^{-1}$ is the empiric energy relaxation rate coefficient.

At constant T the term (16) monotonically increases with increased E_vibr so that with Q_diss = 0 the solution of (15) would finally arrive at some nearly stationary value ${E}_{\text{vibr}}^{\ast }$ required for efficient dissociation to start. For this critical ${E}_{\text{vibr}}^{\ast }$ the equality $Q={Q}_{\mathrm{V}\mathrm{T}}\left(T,{E}_{\text{vibr}}^{\ast }\right)$ is approximately valid. Since only the very beginning of the conversion process is considered it can be taken that n = n_M = n₀, where n₀ is the initial CO₂ density before the process starts. Then the threshold (critical) value of $Q/{n}_{0}^{2}$ which corresponds to the critical ${E}_{\text{vibr}}^{\ast }$ is readily found as:

$$\begin{equation}{\left(\frac{Q}{{n}_{0}^{2}}\right)}_{\text{crit}}={R}_{\mathrm{V}\mathrm{T}}\left(T\right)\left[{E}_{\text{vibr}}^{\ast }-{E}_{\text{vibr}}^{\mathrm{e}\mathrm{q}}\left(T\right)\right]\end{equation} \tag{ 17 }$$

The criterion (17) has a simple physics meaning. For the vibrationally non-equilibrium process to start the rate of producing the vibrationally excited states has to be at least as fast as the rate with which vibrational energy is transferred into translational–rotational modes.

To apply the equation (17) in practice ${R}_{\mathrm{V}\mathrm{T}}\left(T\right)$ and ${E}_{\text{vibr}}^{\ast }$ must be known. R_VT can be extracted from experimental data, but for ${E}_{\text{vibr}}^{\ast }$ there is no simple estimate. The main goal of the model system studies which will be presented below in sections 4–6 is to find plausible estimate of ${E}_{\text{vibr}}^{\ast }$ and to verify the simple formula (17) by comparing it with the results of detailed vibrational kinetics simulations. It will be shown that empiric R_VT obtained in shock tube experiments are valid approximately for the conditions of microwave plasmas as well. One will also see that the term in square brackets $\left[\dots \right]$ in (17) which comes out of the simulation results does not vary a lot from case to case. The magnitude of the corresponding vibrational temperature ${T}_{\text{vibr}}^{\ast }$—around 3000...4000 K—coincides with the expectations based on the shock tube dissociation data (see section 5).

Prior to that main study the assumptions behind equation (7) will be verified (section 4) by the same state-to-state vibrational kinetics model which is later used for the microwave discharge study in section 6. It will be shown that they hold, and in a certain temperature range the state-to-state model calibrated as prescribed by equation (10) does reproduce the experimental relaxation times τ with good accuracy. The experimental data which are used for both calibrating the state-to-state model and defining ${R}_{\mathrm{V}\mathrm{T}}\left(T\right)$ in equation (17) are outlined in the next subsection.

2.3. Experimental data on the rate of vibrational relaxation in CO₂

The experimental data on VT-transfer in CO₂ had been commonly published in form of relaxation times τ defined in terms of equation (1). The experimental τs could be reduced to functions of translational–rotational temperature T only [15, 22]. The most recent review found on the subject is the publication [27]. When comparing that paper with the older report [28] one finds no references on process (6) with M = CO₂ which had not been taken into account in [28]. Therefore, here the rate coefficient of process (6) as a function of T is taken from [28] (process $\left(0{1}^{1}0\right)\to \left(0{0}^{0}0\right)$, M = CO₂, in table IIIa there) as the best fit of the all available experimental data. Apparently, this rate coefficient was obtained by applying equation (2) with ℏω = ℏω₂ rather than equation (10). The correctness of this assumption is confirmed by translating the rate coefficient back into $\tau \left(T\right)$, and comparing the result with selected primary publications on the relaxation time measurements. This comparison is shown in figure 1 where $\tau \left(T\right)$ 'reverse engineered' as described above ('Blauer–Nickerson 1973') is plotted together with results of the shock tube [15, 29] and sound absorption [30, 31] experiments. The agreement is very good.

Zoom In Zoom Out Reset image size

The rate coefficient of process (6) applied in the present work is calculated by using equation (10) with τ derived from the rate coefficient ${R}_{10}^{\text{BN}73}$ of that process in [28]. That is:

$$\begin{equation}{R}_{10}\left(T\right)=\frac{{c}_{\text{vibr}}\left(\bar{T}\right)}{2{c}_{2}\left(\bar{T}\right)}{R}_{10}^{\text{BN}73}\left(T\right),\hspace{25.0pt}\bar{T}=\frac{1}{2}\left(300\,\mathrm{K}+T\right)\end{equation} \tag{ 18 }$$

The coefficient ${R}_{10}^{\text{BN}73}$ corrected in that way is fitted by the equation:

$$\begin{equation}{R}_{10}\left(T\right)=4.48\times 1{0}^{-13}\mathrm{exp}\left(-198.8{T}^{-1/3}+517.6{T}^{-2/3}\right)\end{equation} \tag{ 19 }$$

Here T is in Kelvin and the resulting R₁₀ is in m³ s⁻¹. For T from 250 K to 2000 K (19) approximates (18) with maximum relative deviation 0.5%. The magnitude of the correction factor was already discussed in section 2: the original and corrected rate coefficients are nearly equal at low T, and R₁₀ is twice as large as ${R}_{10}^{\text{BN}73}$ at high T (⩾2000 K).

Finally, the coefficient ${R}_{\mathrm{V}\mathrm{T}}\left(T\right)$ for equation (17) is calculated from the coefficient ${R}_{10}^{\text{BN}73}\left(T\right)$ by applying (2) backwards:

$$\begin{equation}{R}_{\mathrm{V}\mathrm{T}}\left(T\right)={R}_{10}^{\text{BN}73}\left(T\right)\left(1-{\mathrm{e}}^{-\frac{\hslash {\omega }_{2}}{T}}\right)\end{equation} \tag{ 20 }$$

There are at least two readily available vibrational kinetics models of CO₂ developed primarily for the spacecraft atmospheric entry problems which include full manifold of vibrational levels. The model of Kustova et al [17, 18, 32–34] and the model of Vargas et al [19]. In the model [17, 18, 32–34] the states ${\mathrm{C}\mathrm{O}}_{2}\left({v}_{1}{v}_{2}^{l}{v}_{3}\right)$ with same rotational number l are combined into effective states ${\mathrm{C}\mathrm{O}}_{2}\left({v}_{1}{v}_{2}{v}_{3}\right)$. When all vibrational states up to dissociation limit are included in the model there are in total approximately 8000 master equations to integrate. A 1st order perturbation theory—Schwartz–Slawski–Herzfeld (SSH) model [14]—is applied for calculation of the probabilities of transitions between effective states. The model [19] is a modern development which uses completely different approach. In [19] the master equations are written for the individual modes of oscillations rather than individual or effective vibrational states. Their number is reduced to 250. To calculate the transition probabilities the non-perturbative forced harmonic oscillator (FHO) model [35] is applied. Moreover, in the model [19] the process of CO₂ dissociation itself is modelled in detail. Whereas in [17, 18] Treanor–Marrone model is applied to calculate the dissociation rate in [19] the transition from singlet ground state X¹Σ to electronically excited triplet state ³ B₂ via potential surfaces crossing is modelled explicitly as well as subsequent vibrational kinetics of the state ³ B₂. According to shock tube studies the 'forbidden' transition X¹Σ → ³ B₂ is the dominant channel of CO₂ dissociation [36].

The model of Vargas et al is physically much more accurate than the models based on the 1st order perturbation theories and is likely to substitute them in future. Nonetheless, those latter have the advantage that they are based on straightforward textbook physics theories and their numerical implementation is relatively simple. Thus, they still suit well for preliminary investigations preceding more advanced computational studies. The model applied in the present work is of such 'preliminary' kind and is based on the 1st order perturbation theories similar to the model of Kustova et al. The main issue of that latter is too large number of equations which have to be integrated making its implementation computationally challenging [34]. The solution of that problem is further coarse graining by grouping the vibrational states. One of such methods is the 'adaptive energy bin' technique proposed in [37, 38] which is applied to couple the full vibrational kinetics including all levels up to dissociation limit with hydrodynamics.

In the present work the study is restricted to 0D simulations and conceptually more simple approach is applied—the so called 'two-modes approximation' [21]. This approximation originates from the well known fact that the ${\mathrm{C}\mathrm{O}}_{2}\left({v}_{1}{v}_{2}^{l}{v}_{3}\right)$ states with the same good quantum number v_s = 2v₁ + v₂ have very close vibrational energies. Thus, it is plausible to assume that the relaxation between those states is fast. In addition, some of them are in Fermi resonance, see [23, 28], which adds to increasing the speed of mutual relaxation. Experimental information about the real speed of energy exchange between symmetric stretching and bending modes v₁ and ${v}_{2}^{l}$ is very scarce [27], and in two-modes approximation it is taken as postulate that this relaxation is always faster than all other relevant time-scales of the system. Therefore, all vibrational states ${\mathrm{C}\mathrm{O}}_{2}\left({v}_{1}{v}_{2}^{l}{v}_{3}\right)$ with same v_s can be gathered into combined states ${\mathrm{C}\mathrm{O}}_{2}\left[{v}_{\mathrm{s}},{v}_{\mathrm{a}}\right]\ ({v}_{\mathrm{a}}={v}_{3})$. Moreother, to derive the formulas for effective probabilities of transitions between combined states it is assumed that the energies and populations of the elementary states within one combined state ${\mathrm{C}\mathrm{O}}_{2}\left[{v}_{\mathrm{s}},{v}_{\mathrm{a}}\right]$ are exactly equal. This is an extreme assumption, but it greatly simplifies the calculations when 1st order perturbation theories are applied for the probabilities of transitions between elementary states and leads for combined states to simple analytic formulas in this case [39].

The total number of combined states included in the model is determined by the prescribed and fixed dissociation energy E_diss = 5.5 eV which represents the enthalpy change of the transformation CO₂ → CO + O calculated according to [40]. All states ${\mathrm{C}\mathrm{O}}_{2}\left[{v}_{\mathrm{s}},{v}_{\mathrm{a}}\right]$ with total vibrational energies E_v smaller than E_diss are included, plus for each ${\mathrm{C}\mathrm{O}}_{2}\left[{v}_{\mathrm{s}},{v}_{\mathrm{a}}\right]$ closest to the boundary ${E}_{\mathit{v}}\left({v}_{\mathrm{s}},{v}_{\mathrm{a}}\right)={E}_{\text{diss}}$ the states ${\mathrm{C}\mathrm{O}}_{2}\left[{v}_{\mathrm{s}}+1,{v}_{\mathrm{a}}\right]$ and ${\mathrm{C}\mathrm{O}}_{2}\left[{v}_{\mathrm{s}},{v}_{\mathrm{a}}+1\right]$. Those latter are treated as 'unstable' states and the dissociation process CO₂ → CO + O is modelled by assuming their infinitely fast dissociation (technically that means setting for them a very large dissociation rate coefficient). In total 816 master equations for the number densities of regular and unstable combined states are solved.

Especially the extension of the two-modes approximation to high vibrational states is based to large extent on physics intuition rather than on solid theoretical justification. Indeed, for large v_s and v_a the basic assumption that the energies of the elementary states within one ${\mathrm{C}\mathrm{O}}_{2}\left[{v}_{\mathrm{s}},{v}_{\mathrm{a}}\right]$ state are equal is expected to get less and less realistic due to anharmonicity. The way how the dissociation itself is modelled is also apparently only a crude representation of the real physical process. Thorough analysis presented in [19, 41] shows that the potential energy crossings which drive the singlet–triplet transition have highest probability for specific combinations of quantum numbers, and not at vibrational energy of exactly 5.5 eV. The extension of the two-modes approach to high energy states, as well as the preposition that the dissociation is determined only by the total vibrational energy indirectly make use of the concept of 'vibrational quasi-continuum' [42]: that for large vibrational energies precise distinguishing between individual states has no meaning anymore because of very large number of states and very fast transitions between them. Besides the anharmonicity the detailed kinetics of the electronic excitation and of the 'unstable' triplet state, see [19], could also be of importance. In [41] a hypothesis was proposed that the singlet–triplet transition could be one of the rate determining steps which substantially affects the speed of the dissociation process. If this hypothesis is correct the assumption of infinitely fast dissociation of unstable states would be wrong. The model of the present work still can be amended to reflect more accurately the kinetics of the dissociation process by replacing this assumption by finite dissociation rate coefficients dependent on the specific combinations of the numbers v_s, v_a. Despite the above uncertainties and objections comparing with shock tube data confirms that the two-modes model is able to yield realistic rates of CO₂ dissociation for CO₂ diluted in Ar. For dissociation in pure CO₂, however, the reference model has shown to overestimate the literature data by an order of magnitude, although this discrepancy can be explained by the shortcoming of the 1st order perturbation method (see below). These benchmarks will be discussed in section 5.

The full list of transitions between ${\mathrm{C}\mathrm{O}}_{2}\left[{v}_{\mathrm{s}},{v}_{\mathrm{a}}\right]$ states included in the model as well as the details of calculating the corresponding rate coefficients can be found in the open access publication [21]. The rates of vibrational transitions between combined states are calculated by applying the scalings based mainly on the Herzfeld (SSH) theory [43]. A combination of SSH and Sharma–Brau [44] scalings is used for vibrational–vibrational (VV) transitions between asymmetric modes. The absolute value of the rate coefficient R₁₀ of the vibrational–translational (VT) process (6) is adjusted to experimental data by applying equation (10). The absolute values of the rate coefficients of inter-mode transitions between symmetric and asymmetric modes are adjusted to the laser florescence data. The same applies to VV transitions between asymmetric modes. Only the electronic ground state of the molecule is considered.

Calculations of the transition rate coefficients in the model contain two external parameters: f_VV and ${p}_{\text{SSH}}^{\,\mathrm{max}}$. Experimental data on VV-transitions between symmetric modes were found to be very scattered, and therefore for those transitions purely theoretical (SSH) rate coefficients are used. In order to evaluate the consequences of modification of the absolute values of those coefficients they are multiplied by the factor ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$. Its reference (nominal) value is ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}=1$. The range within this parameter may vary is from 10⁻² to 10 estimated approximately on the basis of the experimental data reviewed in [27]. The parameter ${p}_{\text{SSH}}^{\,\mathrm{max}}$ had to be introduced because in SSH theory—as in any 1st order perturbation theory—it is implicitly assumed that the transition probabilities are small, but fulfillment of this assumption is not guaranteed by the theory itself. As a result the SSH formulas can grossly overestimate transition probabilities, and can even formally yield probabilities larger than 1. Therefore, in the model all transition probabilities for which SSH scalings are applied are not allowed to be larger than ${p}_{\text{SSH}}^{\,\mathrm{max}}$. The reference value of this parameter is ${p}_{\text{SSH}}^{\,\mathrm{max}}$ = 1.

Normally for each series of numerical experiments the sensitivity with respect to the parameter ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ and the consequence of the reduction of ${p}_{\text{SSH}}^{\,\mathrm{max}}$ are always checked. In particular, as it will be shown in section 5 the agreement with literature data for the dissociation rate of pure CO₂ could only be achieved with ${p}_{\text{SSH}}^{\,\mathrm{max}}$ = 0.1. This fact is a manifestation of the general limitation of 1st order perturbation theories which can be fully resolved only with non-perturbative methods such as FHO [19, 35]. The two-modes model of the present work is not expected to give exact quantitative predictions, and the result of simulations is rather an expectation range of solutions. The main purpose of the model calculations is to estimate the intervals within which the quantities of interest—such as ${T}_{\text{vibr}}^{\ast }$—may lay. As well as to verify some assumptions behind the semi-empiric models of section 2.

The goal of the present section is to demonstrate that for conditions which mimic those behind a shock front the applied vibrational kinetics model does reproduce the experimental relaxation times. The model is benchmarked using the same test problem as in section 2.1. The initial gas has Boltzmann vibrational distribution with T₀ = 300 K, the translational–rotational temperature T > T₀ is prescribed and fixed. The system relaxes by thermal excitation of vibrational states until the equilibrium vibrational distribution is reached which corresponds to the temperature T. For T ⩽ 2500 K considered here the change of the CO₂ number density due to dissociation on the time-scale of the test is in all cases $<$0.01%.

In a very general sense a state-to-state model can be written as initial value problem for the set of ordinary differential equations:

$$\begin{equation}\frac{\mathrm{d}{n}_{k}}{\mathrm{d}t}=\sum\limits _{i,j}\,{a}_{ij}^{k}{n}_{i}{n}_{j}\end{equation} \tag{ 21 }$$

where n_k,i,j are the number densities of individual species—excited states, and coefficients ${a}_{ij}^{k}$ are functions of the temperature T only. It can be easily shown that changing of the gas pressure p at fixed T does not change the shape of the solution of (21). Indeed, replacing densities with concentrations (molar fractions) c: n_k = c_k n₀, where n₀ is the total initial density of molecules, yields ¹ :

$$\begin{equation*}\frac{1}{{n}_{0}}\frac{\mathrm{d}{c}_{k}}{\mathrm{d}t}=\sum\limits _{i,j}\,{a}_{ij}^{k}{c}_{i}{c}_{j}\end{equation*}$$

The solution of this equation written as a function of n₀ t—the integration time scaled proportional to n₀—depends only on T. Therefore, all the reference model runs here are made for one nominal pressure p = 1 bar. Two tests with p = 0.1 bar and p = 10 bar yield, as expected, exactly same results within relative difference 10⁻⁷ when the time-traces ${E}_{\text{vibr}}\left({n}_{0}t\right)$ are compared.

The time evolution of specific vibrational energy E_vibr calculated by the two-modes model is compared with exponential function (13) with τ defined by (10). To remind, since R₁₀ is defined by (18) this τ is the fit of the primary experimental data. Selected time-traces of E_vibr are shown in figure 2. One can see a very good agreement up to T = 1500 K, the agreement deteriorates at higher temperatures. Quantitative comparison of the numerical solution and the analytic formula is given in table 1 in terms of the relative deviation defined by (14). One can see that for T ⩽ 1500 K the exponential function fits the numerical solution within 4%. That is, the conclusion on the applicability of the simple formula (10) made in section 2.1 is confirmed by the two-modes model. Or, the other way around, the test demonstrates that in conditions which mimic those behind a shock front the state-to-state model calibrated as prescribed by equation (10) yields at T ⩽ 1500 K (to some extent up to 2000 K) the same relaxation behavior as the empiric equation (1).

Zoom In Zoom Out Reset image size

Table 1. Results of the tests of section 4. Comparison of the numerical solution and exponential function in terms of the relative error defined by (14).

	${p}_{\text{SSH}}^{\mathrm{max}}$ = 1		${p}_{\text{SSH}}^{\mathrm{max}}$ = 0.1
T, (K)	, %	ηV1%	, %	ηV1%
350	0.2 (0.2)	100.0	0.2	100.0
400	0.6 (0.5)	99.9	0.6	99.9
500	1.6 (1.5)	99.8	1.6	99.8
600	2.6 (2.6)	99.6	2.6	99.6
800	3.9 (3.9)	98.9	3.9	98.9
1000	4.0 (4.0)	97.7	4.0	97.7
1200	3.1 (3.2)	95.5	3.1	95.5
1500	1.0 (0.6)	89.6	1.0	89.6
2000	12.2 (10.4)	72.0	12.2	72.1
2500	31.9 (28.4)	48.8	31.3	50.4

The numerical simulations allow to verify some assumptions behind the equations (7) and (10). In particular, the linear oscillators assumption which is applied in the two-modes model only for calculation of the matrix elements of transitions. The vibrational energies are calculated with the full anharmonic expression [26]. To estimate the importance of anharmonicity the simulations were repeated with energies calculated in linear (harmonic) approximation. The resulting values of are shown in the second column of table 1 in parentheses, the difference in the outcome of the benchmark is negligible.

The most non-obvious assumption underlying the derivation of (7) is that of the dominance of transition (3) in the thermal excitation of vibrational energy. For its verification in two-modes calculations the relative contribution η_V1 of process (3) in the total energy exchange is calculated:

$$\begin{equation*}{\eta }_{V1}=\frac{{\int }_{t}{Q}_{V1}\left(t\right)\mathrm{d}t}{{\int }_{t}{Q}_{\text{tot}}\left(t\right)\mathrm{d}t}\end{equation*}$$

here ${Q}_{\text{tot}}\left(t\right)$ is the total power transferred from translational–rotational to vibrational modes, ${Q}_{V1}\left(t\right)$ is the power transferred by the process (3) only, the integral is calculated over the whole time interval where the equations are solved. One can see from table 1 that η_V1 ⩾ 90% for T ⩽ 1500 K, and this is exactly the temperature range where the agreement between the two-modes simulations and the semi-empiric fit is good. Above 2000 K the assumption of the dominance of (3) does not hold any more, and the relaxation equations of section 2 are not applicable. At high T the transfer of translational–rotational energy into vibrational modes is taken over by the inter-mode transitions ${\mathrm{C}\mathrm{O}}_{2}\left({v}_{1},{v}_{2}^{l},{v}_{3}\right)+\mathrm{M}\to {\mathrm{C}\mathrm{O}}_{2}\left({u}_{1},{u}_{2}^{m},{v}_{3}\pm 1\right)+\mathrm{M}$ (processes V2 in [21]). At T ⩾ 3000 K this is the dominant mechanism of thermal excitation of vibrational energy. This conclusion drown for high temperatures has to be taken with caution because the absolute values of the rate coefficients of inter-mode transitions in the model are based on the laser fluorescence data available only for T ⩽ 1000 K.

As it was already mentioned in section 3 it is a well known issue of the SSH theory that it can significantly overestimate transition probabilities, especially for high excited states and high T. In order to check which influence this issue could have on the tests of the present section they were repeated with the upper boundary of the SSH transition probabilities (see section 3) reduced to ${p}_{\text{SSH}}^{\mathrm{max}}$ = 0.1. The results are shown in the two last columns of table 1. The difference is only visible starting from T = 2000 K, and does not affect the conclusions drawn above. The issue of the SSH theory is solved in the non-perturbative FHOs model [35]. Recently, FHO calculations of the probabilities of vibrational transitions in CO₂ became available [19, 45, 46]. A comparison between the rate coefficients R of VT-transitions (3) calculated by the SSH and FHO methods has shown that the FHO correction has only a small impact on the calculated rate of this process. This comparison was performed for the set of the FHO rate coefficients of the transitions:

$$\begin{equation}{\mathrm{C}\mathrm{O}}_{2}\left(0{v}_{2}^{{v}_{2}}0\right)+\mathrm{M}\to {\mathrm{C}\mathrm{O}}_{2}\left(0{\left({v}_{2}-1\right)}^{{v}_{2}-1}0\right)+\mathrm{M},\quad \mathrm{M}={\mathrm{C}\mathrm{O}}_{2}\end{equation} \tag{ 22 }$$

taken from [46], file K_CO2v2_CO2_VT.mat, following the description in [45]. The SSH calculations are described in [21, 39], the vibrational energies of the molecules are calculated according to [26]. The rate coefficients of selected transitions calculated by the two methods are plotted in figure 3. Both SSH and FHO calculations are normalized such that for v₂ = 1 the rate coefficient is equal to that of the process $\left(0{1}^{1}0\right)\to \left(0{0}^{0}0\right)$ in [28]. One can see that in the temperature range from 300 to 1500 K the deviation between SSH and FHO scalings is small. For v₂ ⩽ 10 which corresponds to the total vibrational energy of the molecule ⩽0.84 eV the relative difference ${{\epsilon}}_{\text{SSH}}=\left\vert {R}_{\text{FHO}}-{R}_{\text{SSH}}\right\vert /{R}_{\text{FHO}}$ is always ⩽12%. The upper boundary of _SSH is increased for higher v₂, e.g. for v₂ = 20 (energy 1.70 eV) _SSH ⩽ 20%, but even for the highest level v₂ = 63 (energy 5.54 eV) _SSH stays smaller than 32%.

Zoom In Zoom Out Reset image size

The FHO model also allows to calculate the probabilities of multi-quantum transitions which are not accounted for in first order perturbation theories. Comparison between the rate coefficients of the single-quantum transitions (22) and the corresponding two-quantum rates in the data set [46] shows that at room temperature T = 300 K the two-quantum transitions are always by more than a factor 100 less probable. The importance of multi-quantum transitions is increased at higher T. However, the ratio ${R}_{1-q}/\left(2{R}_{2-q}\right)$—where R_1−q is the rate of a single-quantum transition, and R_2−q is the rate of the two-quantum transition from the same v₂—stays $>$10 at T ⩽ 1000 K even for the largest v₂ = 63 ($>$5 at T ⩽ 1500 K). With decreasing v₂ the ratio ${R}_{1-q}/\left(2{R}_{2-q}\right)$ is increased.

For transitions other than (3), in particular for VV transfer, the difference between FHO and SSH scalings can be much large. As has been shown above, at T < 2000 K those other transitions are relatively unimportant for thermal excitation of vibrational energy. At the same time, the inter-mode and VV transitions are the processes which govern the fast energy exchange between different kinds of vibrational modes and bring them to the Boltzmann distribution with single temperature T_vibr. Although the absolute values of their coefficients are adjusted to experimental data the SSH scaling can lead to very inaccurate values for high excited states. Nevertheless, as long as those processes in the model are fast enough this inaccuracy must not have large impact on the final result. Strong coupling between different types of vibrations in the two-modes simulations agrees with the fact that in shock tube experiments only one or two very close relaxation times were detected [15, 22].

The 'thermal dissociation' test serves as experimental check of the calculation of dissociation rate in the vibrational kinetics model. The simulations are performed for conditions of temperature equilibrium—translational–rotational temperature T is equal to vibrational temperature T_vibr. Their results are compared with the dissociation rates measured is shock tube experiments. Behind a shock front vibrational relaxation is much faster than the dissociation, and the assumption that this latter takes place under T = T_vibr conditions is valid. The most reliable experimental data for the speed of reaction CO₂ + M → CO + O + M are available for CO₂ diluted in Ar. Therefore, the first benchmark here is made for M = Ar rather than for M = CO₂ which will be discussed later. The experimental data of [36, 47] are used which were obtained for molar fraction of CO₂ in Ar in the range from 10⁻⁵ to 10⁻².

A comparison of the model with [36, 47] was already shown in the previous paper [21], see section 3.1 and figure 4 there. However, in [21] M = Ar was modelled by M = CO₂ with disabled VV transitions. Here the benchmark of [21] is improved by removing that inconsistency. The vibrational transitions of CO₂ in collisions with Ar which are taken into account in the present model are listed in table 2. They are added to the transitions between CO₂ molecules, table 1 in [21]. The process tags in the first column of table 2 correspond to the notation of [21].

Zoom In Zoom Out Reset image size

The rate coefficients of process V1 are based on the shock tube data of [15]. The relaxation times measured in [15] are in a good agreement with independent measurements of [48]. In [15] vibrational relaxation in pure CO₂ and in Ar/CO₂ mixture with 10% molar fraction of CO₂ was investigated. To reconstruct the rate R₁₀ of (6) with M = Ar $({R}_{10}^{\mathrm{A}\mathrm{r}})$ it is assumed that the measured relaxation time τ of the mixture can be expressed as:

$$\begin{equation*}{\tau }^{-1}=cn{R}_{\mathrm{V}\mathrm{T}}^{{\mathrm{C}\mathrm{O}}_{2}}+\left(1-c\right)n{R}_{\mathrm{V}\mathrm{T}}^{\mathrm{A}\mathrm{r}},\end{equation*}$$

where ${R}_{\mathrm{V}\mathrm{T}}^{\mathrm{M}}$ is the energy relaxation rate coefficient in terms of equation (16) for the corresponding type of 2nd particles M, n is the total number density of Ar and CO₂ (the density which determines the total pressure), c is the concentration—molar fraction—of CO₂ in the mixture. ${R}_{\mathrm{V}\mathrm{T}}^{{\mathrm{C}\mathrm{O}}_{2}}$ is reconstructed using the data of the pure CO₂ experiment in the same paper [15]. After that the ${R}_{10}^{\mathrm{A}\mathrm{r}}$ is found by applying equation (10) with ${\tau }^{-1}={R}_{\mathrm{V}\mathrm{T}}^{\mathrm{A}\mathrm{r}}n$. The resulting least square fit is:

$$\begin{equation*}{R}_{10}^{\mathrm{A}\mathrm{r}}\left(T\right)=6.318\times 1{0}^{-14}\mathrm{exp}\left(-162.3{T}^{-1/3}+271.6{T}^{-2/3}\right)\end{equation*}$$

In this section in equations for rate coefficients the temperature T is always in Kelvin, and the resulting rate coefficients are in m³ s⁻¹.

The rate coefficients of the processes V2 (see table 2) are calculated in exactly same way as it was done for M = CO₂ in [21], see section 2.1 there. The absolute values of the basis rate coefficients for the lowest vibrationally excited states are obtained from the experimental decay rate of the state ${\mathrm{C}\mathrm{O}}_{2}\left(0{0}^{0}1\right)$. The SSH theory calculations are only used to find the branching ratios. The rate coefficients R_V2 of the ${\mathrm{C}\mathrm{O}}_{2}\left(0{0}^{0}1\right)$ decay determined by laser diagnostics are taken from [49, 50]. The experimental data in [49, 50] are only available for T up to 1000 K, and for the benchmark in question they have to be extrapolated up to 4500 K. The rate coefficient ${R}_{10}^{\mathrm{A}\mathrm{r}}\left(T\right)$ had to be extrapolated to high T as well, but for the process V1 this is a smaller issue because: (i) the data are available for T up to 2600 K; (ii) this process was found to be less important for thermal dissociation than the group of processes V2. The small systematic difference between the two references [49, 50] gets significant when they are extrapolated to high T. Therefore, two cases are considered. The decay rate coefficient ${R}_{V2}^{\mathrm{I}}$ which is the least square fit of the data of [50], and ${R}_{V2}^{\mathrm{I}\mathrm{I}}$ which is the fit of [49]:

$$\begin{equation*}{R}_{V2}^{\mathrm{I}}\left(T\right)=1.046\times 1{0}^{-8}\mathrm{exp}\left(-371.0{T}^{-1/3}+1169.2{T}^{-2/3}\right)\end{equation*}$$

$$\begin{equation*}{R}_{V2}^{\mathrm{I}\mathrm{I}}\left(T\right)=6.490\times 1{0}^{-12}\mathrm{exp}\left(-235.4{T}^{-1/3}+581.9{T}^{-2/3}\right)\end{equation*}$$

Calculation of the transition rate coefficients of individual combined states applying SSH scaling is described in detail in [21], section 2.1 there, and in [39]. The only difference is that the reduced mass μ and the characteristic interaction length a⁻¹ must be taken for the pair CO₂ + Ar.

The consistency of the final CO₂/Ar vibrational kinetics model and the relaxation data of [15] is checked by the same 'behind a shock front' test as in section 4. The test with 10% CO₂ in CO₂/Ar mixture was performed for the temperatures which correspond exactly to the experimental points of [15]. The time trace ${E}_{\text{vibr}}\left(t\right)$ calculated for each experimental case is compared with exponential function (13) where τ is the experimental relaxation time for that temperature point. The maximum relative difference defined by (14) is 17%. As it was seen in section 4 one would expect that at high T the excitation of vibrational energy is not dominated by the processes V1, and R₁₀ calculated by equation (10) is not valid anymore. Hence, the relatively good agreement at the highest T is to some extent a coincidence. Although for M = Ar even in the case with the highest T = 2600 K the channel V1 is still responsible for more than 60% of vibrational energy gain, and (10) is still valid to a certain degree of accuracy.

The results of the 'thermal dissociation' simulations with the updated CO₂/Ar model are compared with experimental data [36, 47] in figure 4. The rate coefficients R_diss are calculated by a least square fit of the time trace ${n}_{{\mathrm{C}\mathrm{O}}_{2}}\left(t\right)$ with the function ${n}_{{\mathrm{C}\mathrm{O}}_{2}}\left(t=0\right)\mathrm{exp}(-{R}_{\text{diss}}{n}_{0}t)$, where ${n}_{{\mathrm{C}\mathrm{O}}_{2}}\left(t\right)$ is the number density of CO₂ molecules as a function of time t and n₀ is the initial total number density of Ar and CO₂. It was already shown in the previous work [21] that the calculated R_diss does not explicitly depend on initial pressure (or on n₀). Those tests were not repeated here, all calculations are performed for the nominal pressure 1 bar. The results in figure 4 are all obtained with c = 10⁻⁴. Only a little variation in the results was seen when c is increased. The maximum relative difference between the solutions obtained with c = 10⁻⁴ and c = 10⁻² is 36%, and the relative difference of the solutions obtained with c = 10⁻⁴ and c = 10⁻³ is always smaller than 4%.

As one can see in figure 4 the difference between the simulations made with ${R}_{V2}^{\mathrm{I}}$ and ${R}_{V2}^{\mathrm{I}\mathrm{I}}$ is not large. In both cases at temperatures below 3000...3500 K the model tends to underestimate the rate of dissociation, and at higher temperatures tends to overestimate it. As it was already discussed above in section 3 in the simulations the effective transition probabilities calculated by the SSH theory scalings are always constrained by the prescribed upper boundary ${p}_{\text{SSH}}^{\,\mathrm{max}}$. As it is shown in figure 4 reducing ${p}_{\text{SSH}}^{\,\mathrm{max}}$ from its nominal value 1 to 0.1 slightly improves the agreement at high temperatures. The dashed lines in figure 4 give the estimated uncertainty of the experimental data [36, 47] provided by the least square fit:

$$\begin{equation}{R}_{\text{diss}}\left(T\right)=2.388\times 1{0}^{5}{T}^{-5.2576}\mathrm{exp}\left(-\frac{68\,715}{T}\pm 2\times 0.54\right)\end{equation} \tag{ 23 }$$

where 0.54 is the standard deviation of residuals. Formally the numerical results obtained with ${p}_{\text{SSH}}^{\,\mathrm{max}}$ = 0.1 are almost within this uncertainty. Nevertheless, the calculated dissociation rate apparently grows steeper with temperature than in the experiments. This discrepancy is likely a consequence of the oversimplified way in which the method of 'unstable states' applied in the model to simulate the dissociation reflects the real two-step process as discussed in section 3. The fact that the steepness, hence, effective activation energy, is overestimated by the model speaks in favor of the hypothesis proposed in [41] that the speed of dissociation is also constrained by the speed of the singlet-triplet transition and kinetics of subsequent vibrational excitation of the electronic state ³ B₂.

In CO₂ heavily diluted in Ar VV transitions are effectively disabled and are not involved in the benchmark. In order to test the whole model the rates of dissociation in pure CO₂ have to be compared. Opposite to M = Ar there are no direct experimental data for M = CO₂. There are two main reasons for that. First is that when processing the results of shock tube experiments there is no straightforward way to separate the primary process CO₂ + M → CO + O + M from the oxygen exchange reaction CO₂ + O → CO + O₂. Second is that the reconstruction of the flow parameters in molecular gases is less certain and more error prone than in atomic gases. Because of that the rate coefficients of the process CO₂ + M → CO + O + M with M = CO₂ which can be found in the literature are those for M = Ar multiplied by an enhancement factor. This latter is estimated roughly on the basis of experimental data and its value varies from 2.5 to 11 depending on the (literature) source. In figure 5 the model calculations are compared with the fits of dissociation rate coefficients obtained by this method taken from two data reviews [51, 52].

Zoom In Zoom Out Reset image size

One can see that the reference model overestimates the literature data by an order of magnitude. Good agreement can only be achieved with reduced ${p}_{\text{SSH}}^{\,\mathrm{max}}$. (Agreement with the data from [52] in figure 5 appears to be exact in this case, but this is coincidence). The results are found to be much less sensitive with respect to the parameter ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ (see section 3). Its increase up to ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ = 10 leads to only a very moderate increase of the rate coefficient, and its decrease down to ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ = 10⁻² does not lead to its reduction. Unless ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ is much larger than 10 the VV-transitions between asymmetric modes dominate the transitions to unstable states and thus the rate of dissociation. The fact that the result is sensitive with respect to reduction of ${p}_{\text{SSH}}^{\,\mathrm{max}}$ is manifestation of the general shortcoming of the 1st order perturbation theories which was already discussed in section 3. In the microwave discharge studies in section 6 below the ${p}_{\text{SSH}}^{\,\mathrm{max}}$ = 0.1 case is always part of sensitivity tests. Proper solution of this issue would be application of the FHO model [19, 35] instead of SSH at least for the VV-transitions between asymmetric states.

The benchmark discussed in this section still can be seen at best only as a partial validation of the model because the temperatures T > 2000 K are not in the range relevant for investigating the non-equilibrium vibrational dissociation which is the main goal of the present work. Nevertheless, this test shows that despite its very approximate nature especially for high excited states the applied vibrational kinetics model is capable of producing the correct orders of magnitude of the CO₂ dissociation rate.

The vibrational kinetics model [21] calibrated, verified and partially validated as described in the two previous sections has been applied to study the onset of the non-equilibrium vibrational dissociation of CO₂ in microwave gas discharges. To simulate the electron impact excitation of vibrational states a model system similar to that of [21] has been adopted with one important modification. Instead of degree of ionization the total specific (volumetric) input power Q into electrons is prescribed and fixed in a model run. In practice this means that the electron density n_e is calculated from the power balance:

$$\begin{equation}{n}_{\mathrm{e}}=\frac{Q}{{\sum }_{k,i}\,{R}_{ei}^{k}{E}_{i}^{k}{n}_{i}}\end{equation} \tag{ 24 }$$

where n_i is the number density of the excited state i, ${R}_{ei}^{k}$ is the total rate coefficient of all the electron impact processes which transfer the species i into species k, and ${E}_{i}^{k}$ is the energy spend on average by electrons in each transformation from i to k.

The term Q only includes vibrational excitations and does not include other types of the electron processes such as direct electron impact dissociation or ionization. This approximation is applicable when the electron kinetic energy goes predominantly into production of vibrationally excited states. In weakly ionized plasmas the electron energy distribution function (EEDF) as well as the breakdown of electron energy over different channels is determined by reduced field E/n (see e.g. [11] or [53]) where E is the strength of electric field, and n in the case in question is the number density of CO₂ molecules. Numerical calculations [11–13] had shown that for CO₂ there is a certain range of E/n around 10 Td—approximately between 5 Td and 50 Td—in which most of the electron energy is indeed spent on the excitation of vibrational states. The order of magnitude of E/n indicated above is typical for the microwave sustained plasmas [54, 55]. Hence, the model system considered here is relevant for the conditions of microwave discharges.

The rate coefficients ${R}_{ei}^{k}$ are calculated assuming Maxwellian EEDF. This assumption is shown to be approximately valid in microwave plasmas for the processes for which the high energy tail of EEDF is of minor importance, see [56–58]. A more accurate calculation of EEDF was considered to be superfluous since no excited state resolved data are readily available for the electron impact transitions anyway, thus, only a very basic model of those processes is applied. Only one-quantum transitions are taken into account, and same rate coefficients are applied for all excited states, see [21]. The range of reduced field 5 Td ≲ E/n ≲ 50 Td corresponds to the average electron energies of the order of 1 eV [11, 12]. Therefore, electron temperature T_e = 1 eV is assumed in the reference simulations here. The effect of variation of T_e near this value was found to be very small, see below.

Applicability of the model to gas discharge conditions—although not exactly to microwave plasmas—is checked by simulations performed for the pulsed glow discharge experiment reported in [59]. Only pure CO₂ is included in the model. Therefore, as a test case the 'single pulse' experiment presented in figure 7(b) of [59] is taken in which the off-time is long enough to purge the gas from CO and O₂ produced due to electron impact processes during the discharge. The duration of the pulse is 5 ms. The pulse voltage is 1.5 kV, discharge current is 50 mA—total discharge power 75 W. The length of the plasma is approximately estimated to be equal to the distance between electrodes—17 cm, and the plasma diameter is the diameter of the discharge tube—2 cm. This yields specific input power 1.4 × 10⁶ W m⁻³. The gas pressure is 6.7 mbar. At the beginning of discharge when the gas temperature T_g is around 300 K the reduced field is E/n = 55 Td. According to [11, 12] for this value of E/n 70%...80% of input power goes into vibrational excitations. Therefore, in the calculations the specific power which goes into vibrational states is set to Q = 10⁶ W m⁻³. The electron temperature for the value of E/n above is set to T_e = 1.5 eV according to [11].

To reflect the gas heating in the model a time-profile of translational–rotational temperature T is prescribed according to the time dependency of T_g shown in [59] (figure 7(b)). The master equations are solved in constant pressure form. The dissociation of CO₂ from high vibrational states in the conditions in question is found to be zero. The results of simulations are shown in figure 6 as time evolution of the temperatures of symmetric T_s and asymmetric T_a modes. The profile of T_g and the definitions of T_s and T_a can be found in appendix B, as well as the constant pressure form of the equations (B.3). The time traces of T_s, T_a in figure 6 can be compared with T_1,2 and T₃ respectively in figure 7(b) of [59]. Qualitatively the behavior of the both temperatures is reproduced: T_a is rapidly increased and saturates, and the increase of T_s nearly follows that of T_g. Quantitatively the model overestimates the peak T_a and both T_s, T_a in the end of discharge phase by more than 100 K.

Zoom In Zoom Out Reset image size

At least part of this discrepancy can be ascribed to the lack of self consistency in the present model. As T_g is increased at constant pressure the reduced field E/n is increased too. Thus, the fraction of input power which goes into vibrationally excited states has to get smaller, see [11, 12], and Q has to be reduced as well. Indeed, calculations with Q = 5 × 10⁵ W m⁻³ shown in figure 6 by dashed lines yield temperatures in the end of discharge phase which are closer to the experiment. The parameters ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ and ${p}_{\text{SSH}}^{\mathrm{max}}$ of the two-modes model, see section 3, were found to have very weak influence on the outcome of the test. Variation of the parameter ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ from 10⁻² to 10, and reduction of ${p}_{\text{SSH}}^{\,\mathrm{max}}$ to 0.1 have found not to affect the results shown in figure 6. Two more factors not included in the simulations may be responsible for the remaining disagreement with experiment besides the vibrational kinetics itself. First are the losses of vibrational energy due to spontaneous emission which can become significant at low pressures. Second is the hydrodynamics of the gas expending due to heating approximated here by the constant pressure assumption.

The set of master equations solved by the model in the discharge zone Q > 0 can be written in the same form as (21) with the electron excitation term added:

$$\begin{equation*}\frac{\mathrm{d}{n}_{k}}{\mathrm{d}t}=\sum\limits _{i,j}\,{a}_{ij}^{k}\left(T\right){n}_{i}{n}_{j}+\sum\limits _{i}\,{R}_{ei}^{k}\left({T}_{\mathrm{e}}\right){n}_{\mathrm{e}}{n}_{i}\end{equation*}$$

The following transformation applies when Q = const. Equation (24) is substituted for n_e, the number densities are replaced by concentrations c_k : n_k = c_k n₀. Finally, specific energy input (SEI) per one initial molecule = Qt/n₀ is introduced, bringing the equations above into the form:

$$\begin{equation}\frac{\mathrm{d}{c}_{k}}{\mathrm{d}{\epsilon}}=\frac{{n}_{0}^{2}}{Q}\sum\limits _{i,j}\,{a}_{ij}^{k}\left(T\right){c}_{i}{c}_{j}+\frac{{\sum }_{i}{R}_{ei}^{k}\left({T}_{\mathrm{e}}\right){c}_{i}}{{\sum }_{k,i}{R}_{ei}^{k}\left({T}_{\mathrm{e}}\right){E}_{i}^{k}\left({T}_{\mathrm{e}}\right){c}_{i}}\end{equation} \tag{ 25 }$$

As one can see from (25) for the constant T and T_e the solution expressed in terms of c_i and depends only on the parameter $Q/{n}_{0}^{2}$. This is exactly same parameter as in the approximate threshold criterion (17) in section 2.2.

The outcomes of the model calculations for selected fixed temperatures T are shown in figure 7. The results are presented in terms of $Q/{n}_{0}^{2}$ and the 'process rate' χ defined as the fraction of initial CO₂ molecules which dissociate in the process CO₂ → CO + O. For all model runs this quantity is taken at the time instants which correspond to = 3 eV:

$$\begin{equation}\chi =1-\frac{{n}_{{\mathrm{C}\mathrm{O}}_{2}}\left({\epsilon}=3\enspace \mathrm{e}\mathrm{V}\right)}{{n}_{{\mathrm{C}\mathrm{O}}_{2}}\left({\epsilon}=0\right)}\end{equation} \tag{ 26 }$$

Conversion in the afterglow (Q = 0) is not taken into account since in the previous work [21] the conversion on that phase was found to be insignificant. The chosen final value of is close to the net enthalpy change of the total process CO₂ → CO + $\frac{1}{2}$O₂, ΔH_f = 2.9 eV [40] which is the ideal cost of producing one CO molecule. That is, which is larger than that value will lead to inevitable losses of energy not into the target process even at 100% conversion of CO₂ into CO and is therefore not considered.

Zoom In Zoom Out Reset image size

The CO₂ kinetics model applied for the simulations here is incomplete in terms of chemistry—it does not include the secondary processes with O-atoms and reverse reactions. Also, interaction of CO₂ with reaction products CO, O₂, O is not taken into account. Therefore, the calculated χ is not expected to be a quantitatively correct prediction, but rather an indicator that the process has been started. At the same time, the model can be considered as valid in the vicinity of the threshold values of $Q/{n}_{0}^{2}$ there the total CO₂ conversion is yet low and the influence of products on the primary dissociation process is yet small. That is, the determination of the $Q/{n}_{0}^{2}$ threshold above which the efficient dissociation from high vibrational states can start is the main result of the simulations, not the absolute values of χ at higher $Q/{n}_{0}^{2}$.

All the simulations presented in figure 7 were performed technically for the nominal pressure p = 100 mbar. To confirm that, as suggested by (25), there is no explicit pressure dependence the calculations were repeated with different values of p. The relative differences of χ calculated with p = 50 mbar and p = 200 mbar and those with p = 100 mbar for the same values of $Q/{n}_{0}^{2}$ (only χ > 0.1% are taken into account) are found to be always $<$0.7%.

The solid lines in figure 7 present the solutions obtained with the reference model, and the dashed lines are the results of the simulations where the model uncertainties are evaluated. In particular, the impact of the parameters ${p}_{\text{SSH}}^{\,\mathrm{max}}$ and ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ described in section 3. The effect of reducing the parameter ${p}_{\text{SSH}}^{\,\mathrm{max}}$ from 1 to 0.1 is visible although relatively small for T = 300 K, but tend to increase at higher T. For the parameter ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ the values 10⁻², 10⁻¹ and 10 were tested. One can see a relatively large modification of the solution, especially when ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ is increased. At the same time, one can also say that both the variations of ${p}_{\text{SSH}}^{\,\mathrm{max}}$ and ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$—although they have substantial impact on the steepness of the χ increase with increased $Q/{n}_{0}^{2}$—have a very limited effect on the threshold value of this parameter. Finally, the effect of changing the electron temperature T_e was also investigated. In the reference model (T_e = 1 eV) 69%...81% of Q is deposited into asymmetric modes v_a. Two model runs have been added with T_e = 0.5 eV and 2 eV. In the first case the fraction of Q which goes into asymmetric modes is slightly increased up to 80%...92%, in the second case it is reduced down to 57%...70%. As one can see in figure 7 those variations have negligible impact on the solution.

6.1. Process threshold ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$

In this subsection the results of the vibrational kinetics simulations above are compared with the simple formula (17). First, it has to verified that equation (16) with R_VT determined from gas dynamic experiments is approximately valid for the microwave discharge conditions as well. For that the magnitude of the term Q_VT as it appears in the simulations is compared directly with Q_VT calculated by (16). Examples of this comparison are shown in figure 8. Instead of E_vibr the more illustrative quantity T_vibr is used. To translate the simulated E_vibr into T_vibr equations (A.2) and (A.3) from [21] are used, E_vibr in (16) is calculated as a function of T_vibr by applying (12). The ${R}_{\mathrm{V}\mathrm{T}}\left(T\right)$ in (16) is calculated by (20), with R₁₀ taken from [28]. To remind, the two-modes model [21] which has been applied for the simulations here is adjusted to match the experimental rates of vibrational relaxation which correspond to exactly this choice of ${R}_{\mathrm{V}\mathrm{T}}\left(T\right)$.

Zoom In Zoom Out Reset image size

One can see in figure 8 that the agreement is relatively good. This statement is quantified by examining the ratio of ${Q}_{\mathrm{V}\mathrm{T}}\left({T}_{\text{vibr}}\right)$ taken from the simulations and ${Q}_{\mathrm{V}\mathrm{T}}\left({T}_{\text{vibr}}\right)$ calculated with (16). For T = 300...600 K this ratio always lies between 0.4 and 1.5. These maximum and minimum are taken over the all modelling runs including those with varied T_e, ${p}_{\text{SSH}}^{\mathrm{max}}$ and ${f}_{\mathrm{V}\mathrm{V}}^{\mathrm{s}}$ (only points with T_vibr > T + 100 K are taken into account). For higher T the upper boundary (maximum) of the ratio above is increased: up to 2.4 for T = 800 K and up to 3.4 for T = 1000 K. Nonetheless, even in those cases (16) can always be taken as a good order of magnitude estimate. This is an expected result since vibrational distribution function of the lower vibrational levels in the gas discharge simulations was found to be close to Boltzmann, see an example in [21]. The populations of the upper levels near the dissociation limit deviate strongly from the Boltzmann distribution. However, since their absolute values are relatively small the high energy tail they form has no overwhelming impact on the total energy exchange rate.

The critical ${E}_{\text{vibr}}^{\ast }$ (or ${T}_{\text{vibr}}^{\ast }$) which enters (17) can be formally determined from the condition that at this ${T}_{\text{vibr}}^{\ast }$ the dissociation term Q_diss starts to be comparable with Q_VT, e.g. by taking ${Q}_{\text{diss}}\left(T,{T}_{\text{vibr}}^{\ast }\right)=0.1{Q}_{\mathrm{V}\mathrm{T}}\left(T,{T}_{\text{vibr}}^{\ast }\right)$. Unlike Q_VT there is no simple way of calculating ${Q}_{\text{diss}}\left(T,{T}_{\text{vibr}}^{\ast }\right)$, in particular because the dissociation proceeds from the non-Boltzmann high energy tail of vibrational distribution. Therefore, the determination of ${T}_{\text{vibr}}^{\ast }$ relies on the results of numerical simulations. The model uncertainties lead to a relatively large scattering of that quantity, but, as it will be seen below, this scattering has only minor impact on the final result.

The ratio ${Q}_{\text{diss}}\left({T}_{\text{vibr}}\right)/{Q}_{\mathrm{V}\mathrm{T}}\left({T}_{\text{vibr}}\right)$ obtained in the simulations is plotted for selected T in figure 9. The maximum and minimum of T_vibr which correspond to ${Q}_{\text{diss}}\left({T}_{\text{vibr}}\right)/{Q}_{\mathrm{V}\mathrm{T}}\left({T}_{\text{vibr}}\right)\approx 0.1$ for each T are quoted in table 3. The numbers are rounded off because no high accuracy is required here. The corresponding values of ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ are obtained by applying equation (17) with ${E}_{\text{vibr}}^{\ast }\left({T}_{\text{vibr}}^{\ast }\right)$ calculated using (12). One can see that the variation of ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ caused by the uncertainty of ${T}_{\text{vibr}}^{\ast }$ is not large. The increase of ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ with increased T is mainly due to the increase of ${R}_{\mathrm{V}\mathrm{T}}\left(T\right)$ rather than due to the term in [...] in (17).

Zoom In Zoom Out Reset image size

Table 3. Semi-empiric estimate of the critical value of $Q/{n}_{0}^{2}$, see section 6.1.

T, (K)	${T}_{\text{vibr}}^{\ast }$, (K)	${\left(\frac{Q}{{n}_{0}^{2}}\right)}_{\text{crit}}$, (W m3)	${\left(\frac{Q}{{p}^{2}}\right)}_{\text{crit}}$, $(\mathrm{W}\enspace {\mathrm{m}}^{-3}\enspace {\mathrm{P}\mathrm{a}}^{-2})$
300	2500..3200	(5..7) × 10−40	30..40
500	2900..4300	(2..3) × 10−39	40..60
800	3300..5000	(6..11) × 10−39	50..90
1200	3700..5500	(2..3) × 10−38	60..100

The range of ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ from table 3 is shown as shaded rectangles in figure 7. The agreement with the process threshold as it appears in the results of the simulation runs is very good. That is, for practical purposes the accuracy of the formula (17) together with the approximate values of ${T}_{\text{vibr}}^{\ast }$ provided in table 3 in determining the process threshold is notworse than the accuracy of a complicated vibrational kinetics model. The estimated magnitude of ${T}_{\text{vibr}}^{\ast }\approx$ 3000...4000 K agrees with the temperatures at which efficient dissociation is observed experimentally under temperature equilibrium T = T_vibr conditions, see figure 4.

Besides the restrictions discussed above straightforward energetic constrains on SEI have to be fulfilled for efficient dissociation of CO₂. SEI is the specific energy input per initial CO₂ molecule. SEI has to be at least larger than the energy per molecule which corresponds to the target vibrational temperature T_vibr: e.g. for T_vibr = 3000 K this is 0.8 eV. At the same time, as it was already mentioned in section 6 SEI cannot be allowed to be much larger than the net enthalpy change of the net process CO₂ → CO + $\frac{1}{2}$O₂, ΔH_f = 2.9 eV [40].

To summarize, for the non-equilibrium vibrational dissociation at low T to be possible at least three conditions have to be fulfilled simultaneously. SEI and the reduced field E/n must lie in a certain—relatively narrow—range, and the parameter $Q/{n}_{0}^{2}$ must be larger than some critical value. As it was shown in section 6.1 the semi-empiric estimate yields ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ = 6 × 10⁻⁴⁰ W m³ (averaged) at T = 300 K. The numerical experiments suggest that within variation of the modelling result no non-equilibrium vibrational process is possible when $Q/{n}_{0}^{2}< 0.5{\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$. At the same time, $Q/{n}_{0}^{2}={\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ merely means that the non-equilibrium process starts to be measurable. In this regime still most of the power Q goes directly into translational–rotational modes. Efficient conversion of CO₂ into CO by the non-equilibrium process can be expected only when Q is further increased, possibly by a factor 5–10, that is, when $Q/{n}_{0}^{2}$ is larger than ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ within an order of magnitude.

The critical $Q/{n}_{0}^{2}$ calculated theoretically can be compared with the peak values observed in real microwave plasmas. In order to do that it is more convenient to replace n₀ with the pressure p = n₀ T (T in energy units):

$$\begin{equation*}\frac{Q}{{p}^{2}}=\frac{1}{{T}^{2}}\frac{Q}{{n}_{0}^{2}}\end{equation*}$$

The corresponding values of ${\left(Q/{p}^{2}\right)}_{\text{crit}}$ can be found in the last column of table 3. Despite T² in the denominator the ${\left(Q/{p}^{2}\right)}_{\text{crit}}$ increases with increased T. This is mainly because ${R}_{\mathrm{V}\mathrm{T}}\left(T\right)$ in (17) grows faster than ∼T².

Very accurate reconstructions of the profiles of the volumetric specific power Q have been performed for the surfaguide CO₂ plasma conversion experiments at the DIFFER Institute. In figure 10 peak values of Q/p² taken from the publications [55, 60, 61] are shown. The largest set of data is taken from Wolf 2020 [60], from figure 4 there. Those data represent the experiments with power absorbed in plasma varying from 435 to 1300 W and flow rate 12 slm (standard liter per minute). The parameters correspond to global SEI = 0.5...1.5 eV. The data from Wolf 2020 are supplemented by the pressure scan from Viegas 2020 [55] (figure 3 there): input power 1400 W, flow rate 18 slm, global SEI = 1.1 eV. Two more points have been added from the recent publication Viegas 2021 [61] where the accuracy of the reconstruction of Q was improved. The peak values of Q are taken from figures 9 and 11 in [61] which present respectively the experiments with pressure 110 mbar, power 950 W, flow rate 9 slm (SEI = 1.5 eV), and 150 mbar, 860 W, 12 slm (SEI = 1 eV).

Zoom In Zoom Out Reset image size

In all those experiments T > 2000 K are observed, and the measured CO₂ conversion is explainable by thermal quenching mechanism. At the same time, one can see that almost all Q/p² in figure 10 lie below even the room temperature ${\left(Q/{p}^{2}\right)}_{\text{crit}}$. That is, from the point of view of the theory of the present work no non-equilibrium vibrational dissociation process was possible in those experiments, and most of the absorbed microwave power should indeed eventually go into the gas heating. Two exceptions are the p = 150 mbar point from Viegas 2021 and the lowest pressure p = 10 mbar point from Wolf 2020. In the first case Q/p² is only marginally larger than ${\left(Q/{p}^{2}\right)}_{\text{crit}}$ for T = 300 K. In the second case Q/p² is by more than a factor 2 larger than ${\left(Q/{p}^{2}\right)}_{\text{crit}}$. The theory does not exclude that in this case the non-equilibrium process could start at T = 300 K, although with yet low process rate, see figure 7(a). In the experiment strong gas heating was observed in this case same as in other cases, with peak temperature 2300 K. So a weak non-equilibrium effect even if it occurred in the beginning could be masked by subsequent heating of the gas and transition to thermal process.

For the IPP Garching plasma torch experiment [10] no detailed reconstructions of the Q profiles have been published. In contracting mode at high pressures the specific input power can be approximately calculated from the estimated plasma diameter d_pl and length L_pl presented in figures 9(b) and 11(b) of [10]. Measurements for the experiments with the following parameters are shown there: pressure 200, 500, 900 mbar, flow rate 20 slm, absorbed power from 900 to 2700 W, global SEI from 0.6 to 1.9 eV. The boundary of the plasma region is approximately defined as the region where the intensity of radiation is above 15% of its maximum. Dividing the absorbed power by the plasma volume calculated as $\frac{\pi }{4}{d}_{\mathrm{p}\mathrm{l}}^{2}{L}_{\mathrm{p}\mathrm{l}}$ gives Q/p² which are always smaller than 6 W (m^{− 3} Pa^{− 2}). That is, always well below ${\left(Q/{p}^{2}\right)}_{\text{crit}}$ even at T = 300 K. In homogeneous mode at low pressures before the plasma contracts the extreme assumption that the plasma occupies the whole discharge tube (inner diameter 26 mm) leads to the following numbers. For L_pl = 40 mm (approximately the high of the waveguide) and the largest absorbed power 2700 W Q = 1.3 × 10⁸ W m⁻³. For the smallest reported pressure 20 mbar that yields Q/p² = 32 W (m^{− 3} Pa^{− 2}) which does not exceed the critical value. One has to keep in mind that the local peak values can be higher than the averaged ones.

To conclude, for the selection of experiments reviewed here the parameter Q/p² in most cases formally smaller than the theoretical ${\left(Q/{p}^{2}\right)}_{\text{crit}}$. This agrees with the fact that no significant dissociation due to the non-equilibrium vibrational process could be identified in those experiments. Exceptions are the experiments with very low pressure p ≲ 20 bar where Q/p² can be marginally higher than the room temperature critical value. In particular under low pressure conditions the reduced field E/n constrain may be more critical than the one due to the parameter Q/p². At low pressures the reduced field E/n may become too large and out of the optimal range for preferential excitation of vibrational states. Besides that, at very low pressures one more issue occurs which was left entirely out of the scope of the present work. Here total opacity for spontaneous emission from vibrational levels was always implicitly assumed. At sufficiently low gas density the radiative losses of vibrational energy may become significant. Proper account of them would require solving a radiative transport problem. Analysis of glow discharge experiments in [62] indicated that at p = 5 Torr spontaneous emission from asymmetric vibrational states could already have non-negligible impact on their population.

In this paper activation of the non-equilibrium vibrational dissociation of CO₂ at low translational–rotational temperatures T in microwave gas discharges is investigated by means of semi-empiric and computational models. A governing parameter $Q/{n}_{0}^{2}$ has been introduced, where Q is the volumetric specific power coupled into vibrational modes, and n₀ is the initial number density of CO₂. It has been shown that the non-equilibrium process can only be triggered when $Q/{n}_{0}^{2}$ exceeds some critical (threshold) value. To back and extend the semi-analytic results the approximate two-modes state-to-state model of CO₂ vibrational kinetics [21] is applied.

The model is calibrated and benchmarked to match at T ⩽ 1500 K the characteristic vibrational relaxation times measured in shock tube and acoustic experiments [15, 22, 28]. The dissociation process CO₂ + M → CO + O + M is included in the model via introduction of unstable states whose vibrational energy is larger than prescribed dissociation limit. This approach is partially validated by comparing with shock tube measurements [36, 47] of dissociation rates of CO₂ diluted in M = Ar at T > 2000 K, as well as with literature data [51, 52] on dissociation in M = CO₂. In this latter case good agreement could only be achieved with artificial upper limit applied to transition probabilities calculated by the SSH theory. This result emphasizes drawbacks of the first order perturbation theories and strongly supports the application of non-perturbative models such as FHO [19, 35] despite their higher complexity. Applicability of the model for the non-equilibrium gas discharge conditions with T < 1000 K has been checked by comparing with the pulsed glow discharge experiment [59]. The model is shown to reproduce qualitatively the time evolution of the temperatures of symmetric and asymmetric vibrational modes, although there is no accurate quantitative much with the experimental results.

A series of 0D two-modes simulations has been performed for a model system which mimics vibrational excitation of CO₂ in microwave discharges. In agreement with the analytic consideration the simulations indicate a rapid increase of the rate of the process CO₂ → CO + O with increased $Q/{n}_{0}^{2}$ when this parameter exceeds some critical value. The solutions are found to be sensitive with respect to variation of the uncertain model parameters. However, the position of the process threshold ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ as it appears in the simulation results is relatively stable. A simple estimate of ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ based on the experimental rates of vibrational relaxation agrees well with the modelling results. At T = 300 K this estimate provides ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$ = 6 × 10⁻⁴⁰ W m³ or ${\left(Q/{p}^{2}\right)}_{\text{crit}}$ = 35 W (m⁻³ Pa⁻²) (averaged), where p is the gas pressure. In the temperature range from 300 to 1200 K the ${\left(Q/{p}^{2}\right)}_{\text{crit}}$ is shown to always increase with increased T. Analysis of selected experimental publications [10, 55, 60, 61] has shown that in the real microwave discharges the parameter Q/p² is typically below the estimated critical value. At best the Q/p² has found to be only marginally (within a factor 2–3 at maximum) above the room temperature ${\left(Q/{p}^{2}\right)}_{\text{crit}}$ at very low pressures p ≲ 20 mbar. This observation agrees with the fact that no significant CO₂ conversion due to the non-equilibrium vibrational dissociation in the experiments in questions is reported.

The state-to-state model used in the present paper is applicable only near the $Q/{n}_{0}^{2}$ threshold where the rate of conversion via the primary process CO₂ → CO + O is not yet too high. In order to further extend its applicability the secondary reactions with oxygen, reverse processes, reaction products and their influence on vibrational kinetics must be added. It can be expected that such extensions will not significantly affect the position of ${\left(Q/{n}_{0}^{2}\right)}_{\text{crit}}$, but they might have a large impact on the steepness of the conversion rate growth with increased $Q/{n}_{0}^{2}$.

Adding CO to the model can be seen as a purely technical matter since all the data and models for calculation of the transition rate coefficients are readily available. Adding the oxygen chemistry, on the opposite, can face large difficulties because of the exchange reaction CO₂ + O → CO + O₂. Too little is known about the influence of vibrational non-equilibrium on that process. To which exactly extent vibrational excitation of CO₂ speeds up this reaction? The process has a relatively low enthalpy change, where the vibrational energy of CO₂ will then go? Will it go into vibrational energy of the products or into their translational–rotational energy? Those questions have no definitive answers. The energy efficiency of the overall conversion process CO₂ → CO + $\frac{1}{2}$O₂ triggered by the non-equilibrium dissociation can be larger than 53% only if the reaction CO₂ + O → CO + O₂ is faster than the three-body recombination O + O + M → O₂ + M, and if vibrational energy of CO₂ which takes part in the former goes predominantly into the products vibrational energy.

It is unlikely that all the unknowns will be resolved soon. Thus, it is unlikely that models which correctly predict the efficiency of the overall non-equilibrium CO₂ conversion process can be developed in the near future. At the same time, accurate model predictions of the conditions at which the non-equilibrium process could be triggered are realistic. In particular, the models for CO₂/CO/O₂/O mixtures without exchange reaction CO₂ + O → CO + O₂ are feasible. Computational models of that kind may help to find the parameters under which the low T non-equilibrium process can be obtained and studied experimentally. The $Q/{n}_{0}^{2}$ criterion alone introduced in the present work only makes sense when most of the power absorbed in plasma goes into excitation of vibrational states, that is, only for a certain range of reduced field E/n. To account for this constrain automatically it is mandatory that a model includes detailed electron kinetics. Models which couple self-consistently the electron kinetics and the CO₂ vibrational kinetics and chemistry do exist [13, 57, 58]. Their main issue is that the cross sections of the electron processes with vibrational states of CO₂ are not known to a high degree of accuracy.

The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.26165/JUELICH-DATA/IHO9NP. The up-to-date source code of the numerical model (Fortran) can be found at the following link: https://jugit.fz-juelich.de/reacflow/co2-20.

Derivation of the equation which describes thermal excitation of vibrational energy in diatomic molecules can be found in [14], chapter 19. Here the derivation of [14] is modified for CO₂ which has three modes of oscillation v₁, v₂, v₃. The mode v₂ is double degenerate and has rotational momentum l: ${v}_{2}^{l}$. It is assumed that the molecules are linear (harmonic) oscillators and their vibrational energy E_v can be expressed as:

$$\begin{equation*}{E}_{\mathrm{v}}=\hslash {\omega }_{1}{v}_{1}+\hslash {\omega }_{2}{v}_{2}+\hslash {\omega }_{3}{v}_{3}\end{equation*}$$

Further assumptions are that vibrational–translational (VT) transfer takes place only from the mode ${v}_{2}^{l}$ via process (3), and the number of vibrational states is infinite.

Since in the linear oscillator approximation E_v does not depend on l all states with same v₂ can be gathered into a combined state. There are v₂ + 1 different values of l possible for a fixed v₂. The normalized matrix elements X of all possible VT-transitions from the mode ${v}_{2}^{l}$ are summarized in [39], see equation (27) there:

$$\begin{equation*}{X}_{{v}_{2}^{l}\to {\left({v}_{2}-1\right)}^{l+1}}={v}_{2}-l;\hspace{25.0pt}{X}_{{v}_{2}^{l}\to {\left({v}_{2}-1\right)}^{l-1}}={v}_{2}+l\end{equation*}$$

$$\begin{equation*}{X}_{{v}_{2}^{l}\to {\left({v}_{2}+1\right)}^{l+1}}={v}_{2}+l+2;\hspace{25.0pt}{X}_{{v}_{2}^{l}\to {\left({v}_{2}+1\right)}^{l-1}}={v}_{2}-l+2\end{equation*}$$

The averaged probability p of transition between combined states v₂ and v₂ − 1 calculated by the first order perturbation theory (SSH) is proportional to:

$$\begin{equation}{p}_{{v}_{2}\to {v}_{2}-1}\sim \frac{{\sum }_{l}\left({X}_{{v}_{2}^{l}\to {\left({v}_{2}-1\right)}^{l+1}}+{X}_{{v}_{2}^{l}\to {\left({v}_{2}-1\right)}^{l-1}}\right)}{{\sum }_{l}1}\sim {v}_{2}\end{equation} \tag{ A.1 }$$

The sum ∑_l is taken over all values of l allowed for the corresponding v₂. Same for the transition between combined states v₂ and v₂ + 1:

$$\begin{equation}{p}_{{v}_{2}\to {v}_{2}+1}\sim \frac{{\sum }_{l}\left({X}_{{v}_{2}^{l}\to {\left({v}_{2}+1\right)}^{l+1}}+{X}_{{v}_{2}^{l}\to {\left({v}_{2}+1\right)}^{l-1}}\right)}{{\sum }_{l}1}\sim {v}_{2}+2\end{equation} \tag{ A.2 }$$

Some transitions which appear in the sums are not possible. For example the transition from ${v}_{2}^{{v}_{2}}$ to ${\left({v}_{2}-1\right)}^{{v}_{2}+1}$ is not possible because this latter state does not exist. Such transitions must not be excluded explicitly from the sums because the formulas for X automatically give zero in those cases.

Taking into account the assumption that only VT-transitions between modes ${v}_{2}^{l}$, process (3), are possible the particle balance equations for the number densities n of the combined vibrational states $\left({v}_{1},{v}_{2},{v}_{3}\right)$ read as follows:

$$\begin{align}\hfill \frac{\mathrm{d}{n}_{{v}_{1}{v}_{2}{v}_{3}}}{\mathrm{d}t}& ={k}_{{v}_{2}+1\to {v}_{2}}{n}_{{v}_{1}{v}_{2}+1{v}_{3}}+{k}_{{v}_{2}-1\to {v}_{2}}{n}_{{v}_{1}{v}_{2}-1{v}_{3}}\hfill \\ \hfill & \quad -\left({k}_{{v}_{2}\to {v}_{2}-1}+{k}_{{v}_{2}\to {v}_{2}+1}\right){n}_{{v}_{1}{v}_{2}{v}_{3}},\quad {v}_{2} > 0\hfill \end{align} \tag{ A.3 }$$

$$\begin{equation}\frac{\mathrm{d}{n}_{{v}_{1}0{v}_{3}}}{\mathrm{d}t}={k}_{10}{n}_{{v}_{1}1{v}_{3}}-{k}_{01}{n}_{{v}_{1}0{v}_{3}}\end{equation} \tag{ A.4 }$$

Here k are the transition rates: transition rate coefficients multiplied by the total density of particles M in the process (3). The coefficients k₁₀ and k₀₁ are connected by the detailed balance relation:

$$\begin{equation}\frac{{k}_{01}}{{k}_{10}}=\frac{{w}_{1}}{{w}_{0}}\mathrm{exp}\left(-\frac{{E}_{1}-{E}_{0}}{T}\right)=2\,\mathrm{exp}\left(-\frac{\hslash {\omega }_{2}}{T}\right)=2K\end{equation} \tag{ A.5 }$$

Here w_j is the degeneracy of the vibrational level $\left(0,j,0\right)$ which is equal to j + 1, and E_j is its vibrational energy.

All the other rate coefficients in (A.3) and (A.4), can be connected with k₁₀ using the expressions (A.1), (A.2) and (A.5):

$$\begin{equation}{k}_{{v}_{2}\to {v}_{2}-1}={v}_{2}{k}_{10},\hspace{25.0pt}{k}_{{v}_{2}\to {v}_{2}+1}=\left({v}_{2}+2\right)K{k}_{10}\end{equation} \tag{ A.6 }$$

$$\begin{equation}{k}_{{v}_{2}+1\to {v}_{2}}=\left({v}_{2}+1\right){k}_{10},\hspace{25.0pt}{k}_{{v}_{2}-1\to {v}_{2}}=\left({v}_{2}+1\right)K{k}_{10}\end{equation} \tag{ A.7 }$$

Equations (A.6) and (A.7) are substituted into equations (A.3) and (A.4). For convenience of notation indexes v₁, v₂, v₃ are replaced by i, j, k respectively. The resulting master equations read:

$$\begin{align}\hfill \frac{\mathrm{d}{n}_{\mathit{ijk}}}{\mathrm{d}t}& ={k}_{10}\left(j+1\right){n}_{ij+1k}+{k}_{10}K\left(j+1\right){n}_{ij-1k}\hfill \\ \hfill & \quad -{k}_{10}j{n}_{\mathit{ijk}}-{k}_{10}K\left(j+2\right){n}_{\mathit{ijk}},\quad j > 0\hfill \end{align} \tag{ A.8 }$$

$$\begin{equation}\frac{\mathrm{d}{n}_{\mathit{i0k}}}{\mathrm{d}t}={k}_{10}{n}_{\mathit{i1k}}-2K{k}_{10}{n}_{\mathit{i0k}}\end{equation} \tag{ A.9 }$$

To obtain the energy relaxation equation the (A.8) and (A.9) are multiplied by:

$$\begin{equation*}{E}_{\mathrm{v}}=\hslash {\omega }_{1}i+\hslash {\omega }_{2}j+\hslash {\omega }_{3}k={E}_{ik}+\hslash {\omega }_{2}j\end{equation*}$$

and then their sum over indexes i, j, k is calculated. This sum reads:

$$\begin{align}\hfill & \frac{\mathrm{d}\left[{\sum }_{i,j,k=0}^{\infty }\left({E}_{ik}+\hslash {\omega }_{2}j\right){n}_{\mathit{ijk}}\right]}{\mathrm{d}t}\hfill \\ \hfill & \hspace{25.0pt}={k}_{10}\hslash {\omega }_{2}\sum\limits _{i,k=0}^{\infty }\,\sum\limits _{j=1}^{\infty }\left\{j\left(j+1\right){n}_{ij+1k}-{j}^{2}{n}_{\mathit{ijk}}\right.\hfill \\ \hfill & \hspace{25.0pt}\quad +\left.Kj\left(j+1\right){n}_{ij-1k}-Kj\left(j+2\right){n}_{\mathit{ijk}}\right\}\hfill \\ \hfill & \hspace{25.0pt}\quad +{k}_{10}{E}_{ik}\sum\limits _{i,k=0}^{\infty }\left[\sum\limits _{j=1}^{\infty }\left\{\left(j+1\right){n}_{ij+1k}-j{n}_{\mathit{ijk}}\right.\right.\hfill \\ \hfill & \hspace{25.0pt}\quad +\left.K\left(j+1\right){n}_{ij-1k}-K\left(j+2\right){n}_{\mathit{ijk}}\right\}\hfill \\ \hfill & \hspace{25.0pt}\quad +\left.{n}_{\mathit{i1k}}-2K{n}_{\mathit{i0k}}\right]\hfill \end{align} \tag{ A.10 }$$

The second term in (A.10) equals zero. Indeed, combining the terms without K in the innermost sum yields:

$$\begin{align*}\hfill & \sum\limits _{j=1}^{\infty }\left[\left(j+1\right){n}_{ij+1k}-j{n}_{\mathit{ijk}}\right]+{n}_{\mathit{i1k}}=\left\vert {j}^{\prime }=j+1\right\vert \hfill \\ \hfill & \hspace{25.0pt}=\sum\limits _{{j}^{\prime }=2}^{\infty }\,{j}^{\prime }{n}_{i{j}^{\prime }k}+{n}_{\mathit{i1k}}-\sum\limits _{j=1}^{\infty }\,j{n}_{\mathit{ijk}}\hfill \\ \hfill & \hspace{25.0pt}=\sum\limits _{j=1}^{\infty }\,j{n}_{\mathit{ijk}}-\sum\limits _{j=1}^{\infty }\,j{n}_{\mathit{ijk}}=0\hfill \end{align*}$$

Same for the terms with factor K:

$$\begin{align*}\hfill & \sum\limits _{j=1}^{\infty }\left[\left(j+1\right){n}_{ij-1k}-\left(j+2\right){n}_{\mathit{ijk}}\right]-2{n}_{\mathit{i0k}}\hfill \\ \hfill & \hspace{25.0pt}=\left\vert {j}^{\prime }=j-1,\,j={j}^{\prime }+1\right\vert =\sum\limits _{{j}^{\prime }=0}^{\infty }\left({j}^{\prime }+2\right){n}_{i{j}^{\prime }k}\hfill \\ \hfill & \hspace{25.0pt}\quad -\sum\limits _{j=1}^{\infty }\left(j+2\right){n}_{\mathit{ijk}}-2{n}_{\mathit{i0k}}=\sum\limits _{j=0}^{\infty }\left(j+2\right){n}_{ijk}\hfill \\ \hfill & \hspace{25.0pt}\quad -\sum\limits _{j=0}^{\infty }\left(j+2\right){n}_{\mathit{ijk}}=0\hfill \end{align*}$$

The sum:

$$\begin{equation*}\sum\limits _{i,j,k=0}^{\infty }\left({E}_{ik}+\hslash {\omega }_{2}j\right){n}_{\mathit{ijk}}={E}_{\text{vibr}}n,\quad n=\sum\limits _{i,k,j=0}^{\infty }\,{n}_{ijk}\end{equation*}$$

is the total specific (per unit volume) vibrational energy; E_vibr is the specific vibrational energy per molecule, n is the total number density of molecules. The whole equation (A.10) is reduced to the following form:

$$\begin{align}\hfill \frac{\mathrm{d}\left(n{E}_{\text{vibr}}\right)}{\mathrm{d}t}& ={k}_{10}\hslash {\omega }_{2}\sum\limits _{i,k=0}^{\infty }\sum\limits _{j=1}^{\infty }\left[j\left(j+1\right){n}_{ij+1k}-{j}^{2}{n}_{\mathit{ijk}}\right]\hfill \\ \hfill & \quad +{k}_{10}\hslash {\omega }_{2}K\,\sum\limits _{i,k=0}^{\infty }\,\sum\limits _{j=1}^{\infty }\left[j\left(j+1\right){n}_{ij-1k}\right.\hfill \\ \hfill & \quad -\left.j\left(j+2\right){n}_{\mathit{ijk}}\right]\hfill \end{align} \tag{ A.11 }$$

Calculation of the first sum of (A.11):

$$\begin{align*}\hfill & \sum\limits _{j=1}^{\infty }\left[j\left(j+1\right){n}_{ij+1k}-{j}^{2}{n}_{\mathit{ijk}}\right]=\left\vert {j}^{\prime }=j+1,\,j={j}^{\prime }-1\right\vert \hfill \\ \hfill & \hspace{25.0pt}=\sum\limits _{{j}^{\prime }=1}^{\infty }\left({j}^{\prime }-1\right){j}^{\prime }{n}_{i{j}^{\prime }k}-\sum\limits _{j=1}^{\infty }{j}^{2}{n}_{\mathit{ijk}}\hfill \\ \hfill & \hspace{25.0pt}=\sum\limits _{j=1}^{\infty }\left[\left(j-1\right)j-{j}^{2}\right]{n}_{ijk}=-\sum\limits _{j=0}^{\infty }\,j{n}_{ijk}\hfill \end{align*}$$

Calculation of the second sum of (A.11):

$$\begin{align*}\hfill & \sum\limits _{j=1}^{\infty }\left[j\left(j+1\right){n}_{ij-1k}-j\left(j+2\right){n}_{\mathit{ijk}}\right]=\left\vert {j}^{\prime }=j-1,\,j={j}^{\prime }+1\right\vert \hfill \\ \hfill & \hspace{25.0pt}=\sum\limits _{{j}^{\prime }=0}^{\infty }\left({j}^{\prime }+1\right)\left({j}^{\prime }+2\right){n}_{i{j}^{\prime }k}-\sum\limits _{j=0}^{\infty }\,j\left(j+2\right){n}_{\mathit{ijk}}\hfill \\ \hfill & \hspace{25.0pt}=\sum\limits _{j=0}^{\infty }\left[\left(j+1\right)\left(j+2\right)-j\left(j+2\right)\right]{n}_{ijk}=\sum\limits _{j=0}^{\infty }\left(j+2\right){n}_{ijk}\hfill \end{align*}$$

Substituting the calculated sums back into (A.11) yields:

$$\begin{align}\hfill \frac{\mathrm{d}\left(n{E}_{\text{vibr}}\right)}{\mathrm{d}t}& ={k}_{10}\hslash {\omega }_{2}\left[K\sum\limits _{i,k,j=0}^{\infty }\left(j+2\right){n}_{ijk}-\sum\limits _{i,k,j=0}^{\infty }\,j{n}_{ijk}\right]\hfill \\ \hfill & ={k}_{10}n\left[K{E}_{\text{vibr}-2}+2\hslash {\omega }_{2}K-{E}_{\text{vibr}-2}\right]\hfill \end{align} \tag{ A.12 }$$

where:

$$\begin{equation*}\sum\limits _{i,k,j=0}^{\infty }\,\hslash {\omega }_{2}j{n}_{ijk}=n{E}_{\text{vibr}-2}\end{equation*}$$

is the total vibrational energy stored in the mode ${v}_{2}^{l}$, E_vibr-2 is its specific value per molecule. Transforming the right-hand side of (A.12) and substituting K, equation (A.5), yields finally:

$$\begin{align*}\hfill \frac{\mathrm{d}\left(n{E}_{\text{vibr}}\right)}{\mathrm{d}t}& ={k}_{10}n\left[2\hslash {\omega }_{2}K-\left(1-K\right){E}_{\text{vibr}-2}\right]\hfill \\ \hfill & ={k}_{10}n\left(1-K\right)\left[\frac{2\hslash {\omega }_{2}K}{1-K}-{E}_{\text{vibr}-2}\right]\hfill \\ \hfill & ={k}_{10}n\left(1-{\mathrm{e}}^{-\frac{\hslash {\omega }_{2}}{T}}\right)\left[\frac{2\hslash {\omega }_{2}}{{\mathrm{e}}^{\frac{\hslash {\omega }_{2}}{T}}-1}-{E}_{\text{vibr}-2}\right]\hfill \end{align*}$$

For n = const. and with k₁₀ replaced by ${R}_{10}\left(T\right)n$ this equation gives (4) and (5).

For a constant pressure system it is more convenient to write the master equations in terms of total number of particles N_i of each component i rather than their number densities n_i :

$$\begin{equation}\frac{\mathrm{d}{N}_{i}}{\mathrm{d}t}={S}_{i}V\end{equation} \tag{ B.1 }$$

Here S_i is the specific (per unit volume) particle source of the component i and V is the variable volume of the system determined from the equation:

$$\begin{equation}\sum \,{N}_{i}=V\sum \,{n}_{i}=V\frac{p}{T}\quad {\Rightarrow}\quad V=\frac{T}{p}\sum \,{N}_{i}\end{equation} \tag{ B.2 }$$

Here the sum ∑ is taken over all components i, T is the translational temperature (in energy units), p is the total pressure. Equation (B.1) is readily translated into equation for n_i when one can assume that the total number of particles is conserved, ∑N_i = const. This assumption is valid for the pulsed glow discharge test of section 6 where vibrational temperatures are very low and vibrational dissociation does not take place. In this case substituting N_i = n_i V and (B.2) into (B.1), and recalling that p = const. yields:

$$\begin{equation}\frac{\mathrm{d}{n}_{i}}{\mathrm{d}t}={S}_{i}-{n}_{i}\frac{\mathrm{d}\,\mathrm{ln}\,T}{\mathrm{d}t}\end{equation} \tag{ B.3 }$$

The dependence of the gas temperature T_g shown in figure 7(b) of [59] from time t is fitted by the function:

$$\begin{equation*}{T}_{\mathrm{g}}=\begin{cases}{a}_{0}+{a}_{1}t+{a}_{2}{t}^{2}+{a}_{3}{t}^{3}+{a}_{4}{t}^{4},\quad \hfill & \quad t\leqslant {t}_{\mathrm{p}\mathrm{l}}\hfill \\ {b}_{0}+{b}_{1}t+{b}_{2}{t}^{2},\quad \hfill & \quad t > {t}_{\mathrm{p}\mathrm{l}}\hfill \end{cases}\end{equation*}$$

where a₀ = 300, a₁ = 3.479 6276 × 10⁴, a₂ = 6.289 9233 × 10⁷, a₃ = −1.809 2153 × 10¹⁰, a₄ = 1.425 3195 × 10¹²; b₀ = 675.767 76, b₁ = −1.005 9133 × 10⁵, b₂ = 9.708 9841 × 10⁶; t_pl = 5 × 10⁻³ sec, t is in seconds, T_g in Kelvin.

The temperatures (in energy units) T_s and T_a are defined as:

$$\begin{align*}\hfill {T}_{\mathrm{s}}& ={E}_{[1,0]}{\left[\mathrm{ln}\left(\frac{2{n}_{[0,0]}}{{n}_{[1,0]}}\right)\right]}^{-1},\hfill \\ \hfill {T}_{\mathrm{a}}& ={E}_{[0,1]}{\left[\mathrm{ln}\left(\frac{{n}_{[0,0]}}{{n}_{[0,1]}}\right)\right]}^{-1},\hfill \end{align*}$$

where ${n}_{[{v}_{\mathrm{s}},{v}_{\mathrm{a}}]}$ are the number densities of the corresponding two-modes states, and ${E}_{[{v}_{\mathrm{s}},{v}_{\mathrm{a}}]}$ are their vibrational energies (E_[0,0] = 0). For the case shown in section 6 the temperatures T_a defined for the first 4 levels [0, v_a] are equal within relative difference of 3%, the same applies to T_s defined for the first 10 levels [v_s, 0].