The High Temperature Vibrational Partition Function of Stretching of Triatomic Systems

The statistical thermodynamic model for the vibrational partition function with separated stretching and bending is developed. The model is studied on the example of CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CO}_{2}$$\end{document} molecule for temperature up to 20,000 K with the aim to describe efficient dissociation by deposition of energy mainly to the stretching modes of vibration. The observed separation of bending mode at lower temperatures suggest that it is possible to construct such kinetic model of plasma in which the high vibrational temperature of stretching and the low vibrational temperature of bending are obtained resulting in an efficient dissociation. In particular, the proposed model is extended to ideal-gas version where all the interactions between atoms are taken into account. The idea behind such approach is to eliminate contributions to partition functions stemming from non-interacting dissociated fragments of the molecule. The application areas of such partition functions are discussed and the full vibrational partition functions based on that model are compared with the known data.


Introduction
The vibrational partition function is a basic quantity discussed in statistical mechanics books [1], but its thermodynamical applications are limited. If the temperature is high enough separation of rotations and vibrations fails [2] because of ro-vibrational coupling, theoretical mechanics do not allow to separate rotations from vibrations exactly.
The multi-temperature approaches to thermodynamics are sometimes discussed [3] (in the present manuscript the situation of high vibrational and low rotational temperatures is of interest) but such approach has never been successfully applied. This topic is discussed in the review paper on non-equilibrium phenomena in thermal plasma flows [4].
In general idea of vibrational partition function in the context multi-temperature thermodynamics should be considered incorrect [5] but still it finds its applications in kinetics. The vibrational partition function is used for rate constants in the theory of gas and 1 3 plasma flows. The specific physical situations are hypersonic flows in aerodynamics or shock waves [6][7][8][9]. In those applications the recurring model for vibration-dissociation coupling is the Treanor-Marrone (or Marrone-Treanor) model where the rate constant expression contains vibrational partition function [10]. Note that the Treanor-Marrone model can be extended for state specific rates [11]. The Treanor-Marrone model is discussed in the standard plasma chemistry texts [12].
In particular hypersonic flows of CO 2 , with the methods that needs vibrational partition functions, are studied in [13] and other studies by this group. The recent study for the Mars entry problem [14] describe ionized shock wave (the vibrational excitation of CO 2 by electrons is considered) in which the vibrational temperatures equilibrate to 5000 K and exceed this value close to the shock front. The rate constants in Annaloro [14] study are given up to 10,000 K. The earlier work of those authors [15] describes models for CO 2 plasma jets generated in high enthalpy wind tunnels, in this models all vibrational modes have the same temperature equal to electron temperature which is 11,000 K. This theory needs the separate vibrational partition functions of electronic states (state-to-state model). One more example [16] assumes temperature of asymmetric vibrational mode to be 8000 K. The example of experimental study of vibrations in microwave plasma at very high (translational and vibrational) temperatures is Ref. [17].
Finally, another one application of the vibrational partition function is for the spectral lines intensities in non-equilibrium conditions [18].
All those applications may require very high temperatures (for example in hypersonic flows temperatures exceeding 10,000 K are considered) which is the perfect situation for a fully classical approach. At high temperatures the ro-vibrational partition functions are known to disagree between various studies [19], this happens even for the simplest H + 2 molecule [20]. The same can be expected for vibrational partition functions too. At high temperature the molecules are highly dissociated so that the most important are the simplest ones-namely diatomic and triatomic. The author already dealt with the high temperature partition functions for diatomic molecules [2, [20][21][22][23] and this work is an extension of that theory towards triatomic systems.
For diatomic molecules, the present approach amounts to one dimensional integration with no numerical problems. It will be shown that the triatomic model of stretches, with two or four dimensional integration only, proves to be numerically challenging.
The main aim of the present study is to calculate the vibrational partition function of stretching on example of the CO 2 molecule. The separated stretching model of triatomic molecules presented here gives thermodynamical description of the situation when the energy in the system is distributed mainly into the stretching modes (what results in an effective dissociation) and very little energy is distributed into the bending mode which do not lead to the molecule dissociation. Note that dissociation to C + O 2 (not covered by the model presented here) is a result of symmetric stretching at certain bending angles (only stretches cause C-O bond elongation and, as a result, dissociation of CO 2 ).
A lot of research has been done recently into the CO 2 dissociation [24][25][26], but in those studies there is distinction between symmetric and asymmetric modes which is not the most effective with respect to the thermal decomposition. Despite the presence of the Fermi resonance, the separation of bending observed in the present study suggests that it is possible to construct a kinetic model to obtain high temperature of stretching modes and low temperature of bending mode for CO 2 or other molecules. In principle is seems theoretically to be the most effective model for thermal decomposition. Note also usefulness of reduced-dimensionality models in the research of larger molecules [27].
The recent study [14], which was already mentioned, supports the separation of vibrational modes in plasma, the non-conventional treatment of CO 2 vibrations based on distinguishing the sets of vibrational levels with lower and higher energy (the latter set has one mode excited) is used. This model uses separate vibrational temperature for each mode.
Apart from possible applications to dissociation of CO 2 in plasmas, the linear CO 2 model for vibrational partition function was used recently at very high temperatures for equation of state. [28]. That model assumes the ideal-gas approach and uses harmonic approximation which cannot be valid at extreme temperatures of the study.
It is worth to mention the other methods of phase space integration for partition functions. The classical approach to vibrational partition function due to Varandas group [29,30] is the Monte Carlo based but deals with the bound states only, in such cases kinetic energy should be lower than potential energy and that condition restricts the integration volume to the regions near the equilibrium position of nuclei.
The other approaches, the quantum ones based on path integrals, are applied to ro-vibrational partition functions and allows convenient multi-dimensional integration. The novel one by Szekeres et al. [31] (this publication reviews other numerical methods) does not apply boundedness condition but because of that is applicable to relatively low temperatures (but still high from the point of view of conventional chemistry).
In the present study there is no restriction on nuclei positions and integration spans over large volume of phase space. The application of Hill decomposition [22,32,33] of the partition function to the partition function of interacting atoms (i.e. partition function of actual system) and the partition functions of non-interacting products of thermal decomposition allows to calculate the partition function at arbitrary high temperature; the Monte Carlo simulations (when Hills decomposition is not taken into account) will then cause problems because of sampling non-interacting species and the breakup of molecular system would be observed. To sum up, the present method is complementary to the mentioned path integrals ones-at high temperatures the present classical approach is correct (also no problem with system breakup) but at lower temperatures the Monte Carlo based on path integrals approaches are applicable (or even simple harmonic oscillator approximations).
Finally, it is worth mentioning that in addition to the fact that that the high temperature thermodynamic data are not accurate, the 2017 Plasma Roadmap [34] identifies some thermodynamic data as missing and points out importance of the non-equilibrium conditions as well as some situations where the temperatures are very high.
The "Classical Vibrational Partition Function of Stretching" section of this article describes the classical method of phase space integration for the vibrational partition function. The "Results" section presents the results for the CO 2 molecule. Finally "Conclusions" section concludes the present research.

Classical Vibrational Partition Function of Stretching
The classical partition function was calculated in valence coordinates (two bond lengths and angle between them), the kinetic energy of triatomic system is then given by [35] for the molecule X1 − X3 − X2 , where p 1 , p 2 , p 3 are momenta of respective atoms, m 1 , m 2 , m 3 are masses of respective atoms, 's are the reduced masses, and the valence coordinates are the bond lengths ( r 1 , r 2 ) and the angle between them ( ). Note that other coordinate systems can be used in principle but they are less convenient regarding physical picture; their main advantage seems to be mathematical convenience because of kinetic energy simplification.
This function can be integrated over momenta (with the inverse temperature = 1∕(k B T )): and then the vibrational partition function is then given by where V(r 1 , r 2 , ) is the potential energy surface in valence coordinates.
The vibrational partition function of stretching of the triatomic molecule (i.e. separated bending because of the lower temperature of bending mode) is obtained by setting p = 0 (no bending) and = 0 (linear molecule; for non-linear molecule it would be generalized as = 0 ). In this case the kinetic energy of CO 2 stretching is ( 13 = 23 = CO , m 3 = m C ; subscript S indicates stretches only expression): After integration, as in Eq. 2 but without p momentum, the expression independent of bond lengths is obtained Finally, the vibrational partition function of stretching can be calculated note (2 ) 2 in this expression instead of (2 ) 3 because of one integration (degree of freedom) less. For the carbon dioxide molecule the potential energy surface in valence coordinates according to Cerezo et al. [36] was used. The influence of higher electronic states on partition function will increase with temperature, it can be estimated by calculation of e − E factor to account for the energy shift ( E ) of the lowest excited state that for the 10,000 K the partition function of the lowest excited state calculated without energy shift e − T(r 1 ,r 2 , ,p 1 ,p 2 ,p ) dp 1 dp 2 dp , would have to be multiplied by a 0.0085 factor. There are two reasons to do the calculations with only one potential energy surface-the theoretical one to study the effect on a given potential energy surface (not a mixed effect from more potential energy surfaces) and the practical one based on the idea to treat every state (electronic, and in most detailed approaches also vibrational) as a separate "molecule" (state-to-state models) [13].
In order to get the full vibrational partition function, the above expression has to be complemented by partition function of bending Q b in classical or quantum harmonic oscillator approximation, where the vibrational frequencies ( ss = 0.00606 a.u. = 1330 cm −1 , as = 0.0107 a.u. = 2349 cm −1 , b = 0.003039 a.u. = 667 cm −1 ; all calculations were done in atomic units) were taken form Ref. [37] and the second powers are due to bending mode degeneracy.
Because it is customary to give partition functions without zero point energy in order to transform all the above expressions which contain zero point contribution, the multiplication by exp( ∕2) factors appropriate for each vibrational mode is done and all the results are presented in such form.
If the temperature is high enough, meaning that the high energy states (over dissociation threshold) contribute significantly to the partition function, there is a possibility to include all the interactions (scattering and resonance) in the partition function what allows to perform the calculations at temperature as high as needed with no difficulties. According to the above mentioned theory of Hill (the resulting expression will be called Hill decomposition and the partition function will be indicated by HD superscript) the following expression is the ideal-gas counterpart of Eq. 6 where infinity sign is understood as the appropriate limit. In the above expression the second term was subtracted twice in order to remove the partition function of non-interacting CO + O species. Stretching one of the bonds ( r 1 = ∞ ) causes dissociation and the stretching of the other one ( r 2 = ∞ ) also, because of the symmetry the resulting physical situations are the same and can be treated by subtracting twice the same expression. In the dissociation product CO the bond can have arbitrary length, this means that both expressions describing CO + O system describe also dissociated CO (i.e. describe O + C + O noninteracting atoms), because CO 2 molecule can dissociate on one set of O + C + O atoms it have to be compensated adding the third terms in Eq. 9 which describes O + C + O ( r 1 = ∞ = r 2 , that is both bonds are dissociated).
The above analysis explaining how to write Eq. 9 is based on valence coordinates (i.e. bond lengths) and would be much more difficult in other coordinate systems if possible at all.

3
The Hill decomposition expression can be verified by expanding integration limits to infinity at some very high temperature-in case of Eq. 6 the result will tend to infinity (because it describes non-interacting products of CO 2 thermal decomposition) and the value given by Eq. 9 will remain constant with expanding integration range, so that only this expression leads to well defined ideal-gas high temperature result. The well defined result will give also the classical approach using only bound states [38] (the contributions from resonance and scattering states are not taken into account).
There is also the possibility to take into account resonance and scattering states with the scattering theory [39], the results for water molecule were obtained by considering dissociation to H + OH but not the total dissociation to atoms. The temperatures at which those effects are important allow for classical treatment. In general ideal-gas treatments of triatomic molecules are very scarce because from the viewpoint of standard chemistry it is not so important but in high temperature situations may become crucial.

Results
In the calculation of partition functions the numerical problems are known so that a lot of research is done into devising new numerical methods [31] which are often based on Monte Carlo methods. In light of that fact, the numerical errors have to be considered, it turns out that even only two dimensional integration can cause problems so that the various methods of integration were considered.
In the present study the numerical integration methods implemented in the Mathematica computation system [40] and library of numerical multidimensional integration procedures Cuba [41] (a very reliable free of charge alternative) were used. The results in temperature range of 1000-20,000 K were given in Table 1 where NIntegrate is the default function of Mathematica, NIntegrate+mod is the improved modification of the default function (Min-Recursion → 2, Method → { "GlobalAdaptive", "MaxErrorIncreases" → 100,000, "Symbol-icProcessing" → 0, "SingularityHandler" → None} , PrecisionGoal → 2), AMC is the Adaptive Monte Carlo procedure also implemented in Mathematica, Sauve is one the Monte Carlo procedures implemented in the Cuba library, finally for the check of both theory and numerical integration, the integration with modified NIntegrate of Mathematica of both momenta and positions (without analytical M function) was used (NIntegrate + mod + p). All the methods usually yield very similar results but sometimes some can differ (there are others numerical integration methods implemented both in Mathematica and Cuba). Moreover, not all the methods behave reliably with increasing the upper integration limit, for the subsequent calculations the modified Mathematica integration will be used. In general there are some differences between methods so that it is good to check more than one (in Table 3 it can be seen that AMC method is not good-values at 3000 K, 6000 K and 15,000 K are significantly different than all the other).
Various theoretical methods for the vibrational partition function of stretching are compared in Table 2. To confirm at what temperatures the classical approximation is valid, the quantum (QHO) and classical (CHO) harmonic oscillator approximations based on the frequencies given in the previous section are calculated (for each vibrational mode the expressions as in Eqs. 7 and 8 were used). It can be concluded that above 5000 K the classical approximation gives reasonable estimate and over 10,000 K its accuracy is very good.
The following two columns are the Q vib,s according to Eq. 6 and the ideal-gas version Q HD vib,s of Eq. 9. The ideal-gas approach starts to differ from the conventional expression over 10,000 K, so that below that temperature simpler Eq. 6 is sufficient. In general, the present methodology is applicable at temperatures over 5000 K.
As mentioned before, convergence with respect to the integration limit should be verified, indeed there is a constant value of Q HD vib,s at high distances r 1 = r 2 which confirms that contributions from non-interacting dissociated fragments were correctly removed. For all calculations the integrations are performed up to 10a.u. which was found to be sufficient.
Even though the method is proposed for the vibrational partition function of high temperature stretching complemented with the vibrational partition function of low temperature bending, it is interesting to compare the results in case of equal stretching and bending temperature to the other known results. Table 3 shows the following: quantum harmonic oscillator approximation (separated normal modes), Q vib,s Q b with Q vib,s based on Eq. 6 and quantum harmonic oscillator bending, Q HD vib,s Q b with Q HD vib,s based on ideal-gas expression Eq. 9 and quantum harmonic oscillator bending, Capitelli [3] data and HITRAN [42] data calculated as the respective internal partition function divided by the rotational partition function Q rot = I∕ = 280181∕ (I-moment of inertia).
Interestingly, it turns out that for the CO 2 molecule the approximation with the quantum harmonic oscillator below 5000K gives reasonable agreement with Capitelli or HITRAN data. This confirms that the both separation of vibrational modes, ro-vibrational coupling, and anharmonicity have very minor effect on vibrational partition function up to relatively high temperatures. This should be considered before undertaking summation of energy levels. Good performance of harmonic oscillator approximation which is based on separated modes suggests that this separation at lower temperatures can be the basis for separation of bending also at higher temperatures.

Conclusions
The possibility to separate stretching and bending in linear CO 2 molecule is found to be possible from the viewpoint of statistical thermodynamics and the appropriate classical mechanics methodology applicable for arbitrary high temperatures is developed to calculate the vibrational partition function.
The presented classical approach allows calculation of vibrational partition function which is not becoming more difficult with increasing temperature. For some applications the present method is the most convenient because it can be used at arbitrary high temperature.
The present considerations are applicable to systems with a distinct bending mode, it would not be possible to use it for a molecule like H + 3 where bending is not separated.