Analytic energy-level densities of separable harmonic oscillators including approximate hindered rotor corrections

Energy-level densities are key for obtaining various chemical properties. In chemical kinetics, energy-level densities are used to predict thermochemistry and microscopic reaction rates. Here, an analytic energy-level density formulation is derived using inverse Laplace transformation of harmonic oscillator partition functions. Anharmonic contributions to the energy-level density are considered approximately using a literature model for the transition from harmonic to free motions. The present analytic energy-level density formulation for rigid rotor-harmonic oscillator systems is validated against the well-studied CO + ˙OH system. The approximate hindered rotor energy-level density corrections are validated against the well-studied H 2 O 2 system. The presented analytic energy-level density formulation gives a basis for developing novel numerical simulation schemes for chemical processes. © 2016 Author(s). All article content, except where other-wise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). [http://dx.doi.org/10.1063/1.4963921]


I. INTRODUCTION
The concept of energy-level densities is used throughout physics to describe e.g. semiconductors in condensed-matter physics, 1 proteins in bio-physics, 2 and chemical reactions in gas-phase kinetics. 3 For chemical reactions in particular, the Rice-Ramsperger-Kassel-Marcus (RRKM) theory has emerged as state-of-the-art formulation for microscopic reaction probabilities. The RRKM theory formulates the microscopic rate constants of unimolecular reactions as the ratio of the number of states of the dividing surface and the energy-level density of the reactant at a constant energy-level. 3 Therein, N ‡ (E − E 0 ) is the sum of states of the dividing surface relative to the electronic barrier E 0 , ρ(E) is the energy-level density of the reactant, and h is Planck's constant. The sum of states and the energy-level density are computed from the energy levels accessible by the transition state and the reactant, respectively. These energy levels in turn are deduced from the transition state's and reactant's electronic states and molecular modes, namely translation, rotations, and vibrations. While for translational and rotational modes analytic energy-level densities can be easily obtained via inverse Laplace transformation, 3 the energy-level density of vibrational modes necessitates more sophisticated computations. First energy-level density computations for harmonic oscillators were based on combinatorial considerations, namely the straight-forward evaluation of all possible combinations of distributing energy among multiple harmonic oscillators. 4 Significant improvements in terms of computational efficiency were achieved by Beyer and Swinehart (BS) by making use of the regular nature of harmonic oscillator energy levels. 5 Stein  account for general anharmonic vibrational modes. 6 This modified version of the Beyer-Swinehart algorithm is state-of-the-art for computing energy-level densities of small molecules, within the limitations of discretization. With increasing molecule size, however, the number of internal degrees of freedom increases, leading to significant increases in the computational requirements of combinatorial methods. More recent approaches tried to avoid the combinatorial nature of the Beyer-Swinehart algorithm and its derivative by proposing approximate energy-level densities. [7][8][9][10][11][12] McClurg et al. 7 developed an approximate description for the hindered rotor energy-level density based on a partition function expansion. Wadi and Pollak approximated the energy-level density of large polyatomic molecules with a rather narrow energy-level spacing based on a short-time expansion of the thermal energy-level density. 8 For general systems, Wang and Landau proposed an approach for computing the energylevel density based on a random walk Monte-Carlo algorithm. 9 Recently, Kamarchik and Jasper formulated an approximation for the energy-level density of fully coupled anharmonic vibrational modes. 12 Bozzelli et al. 13 proposed a three-frequency model to extract energy-level density information from heat capacities, which has been used by Allen et al. 14 to compute energy-level densities in an automated manner. The most general approach to obtain energy-level densities is presumably the use of steepest descents of an arbitrary partition function, as described by Baer and Hase. 15 This approach gives the exact energy-level density, if the resolution of this numerical scheme is sufficiently high.
While energy-level densities are widely used in theoretical studies on elementary reactions within the RRKM/Master Equation (ME) framework, 16 the concept of energy-level densities is rarely used in kinetic modeling studies of complete reaction networks. Current kinetic models are typically documented based on the ChemKin format proposed by Kee et al. 17 In there, the temperatureand pressure-dependence of rate constants are represented using the modified Arrhenius equation and either the Lindemann 18 or Troe 19 falloff formulation, respectively. A magnificent exception is the reaction mechanism generator (RMG) by Green et al. 20 This software package estimates molecular properties for obtaining energy-level densities, which in turn are used to estimate the pressure-dependence of rate constants.
The thermochemistry representation in the ChemKin format is based on the NASA format proposed in 1971 21 and uses seven-parameter polynomials to model the enthalpy and entropy in the lowand high-temperature regimes. Energy-level densities, however, directly relate to molecular partition functions, 3 thus to the thermochemistry of the species. Therefore, it would be possible to describe zero-dimensional chemical processes using energy-level densities for RRKM and thermochemistry exclusively. Although this would increase the computational efforts required for simulations, using energy-level density-based simulation schemes would significantly improve the prediction of chemical processes. For this purpose, a closed analytic energy-level density expression is desirable which satisfies the continuous nature of energy-level densities 22 and allows to store a set of parameters rather than actual energy-level densities.
While for numerical simulations of the high-pressure limit the use of energy-level densities would not change the thermochemistry or kinetics, recent theoretical kinetic studies point towards the need for energy-resolved numerical simulations if not in the high-pressure limit. [23][24][25][26][27][28] The formation and consumption of non-thermal intermediates appears to be of high importance under high-temperature and low-pressure conditions. These effects, however, are mostly neglected presently and are hard to implement in state-of-the-art numerical simulation schemes.
Here, an analytic energy-level density formulation will be developed and validated against the well-studied CO +ȮH reaction network. This reaction network describes the conversion of carbon monoxide to carbon dioxide, which is one of the most important reactions in atmospheric and combustion sciences. Stationary points of the reaction network are well-described Previous studies on CO +ȮH mainly focused on evaluating the potential energy surface (PES) of the CO +ȮH reaction network. [29][30][31][32][33] Yu et al. 29 used full coupled-cluster (FCC) and complete basis set (CBS) extrapolations at the CCSD(T)/cc-pVTZ level of theory for updating the PES. Recently, Li et al. 30 proposed a PES based on approximately 35000 ab initio points at the CCSD(T)-F12/aug-cc-pVTZ level of theory. Based on this detailed PES, Xie et al. 31 and Corchado and Espinosa-Garcia 32 performed quasi-classical trajectory simulations to study the dynamics of reactions on the CO +ȮH PES. In addition to the PES studies, Senosiain et al. 34 used a RRKM/ME formulation to predict rate constants for reactions on the CO +ȮH PES at the B3LYP/cc-pVTZ level of theory, including variational effects, conservation of angular momentum, and tunneling. In the present study, the mathematical derivation of an analytic energy-level density formulation is provided for molecular systems with separable harmonic oscillators, including approximate hindered rotor corrections. The analytic formula is validated against reference energy-level densities computed via the Beyer-Swinehart algorithm. First, the well-studied CO +ȮH system is used for validation, which is accurately described by the rigid rotor harmonic oscillator (RRHO) approximation. 34 Second, the one-dimensional hindered rotor of the H 2 O 2 system is used for validation of the approximate hindered rotor corrections. A concluding discussion on the benefits from the presented data format provides insights to potential future numerical simulation frameworks.

II. THEORY
In general, the energy-level density is the inverse Laplace transform of the partition function. 3 Therein, L −1 {Q(s)}(E) denotes the inverse Laplace transformation operator, s = 1/k B T is the inverse Laplace transformation variable, k B is the Boltzmann constant, and T is the temperature. The rather simple expressions for the partition functions of translational and rotational modes allow for directly transforming the respective partition function to the energy-level density. In turn, the vibrational partition function Q ϑ,i is too complex for solving the inverse Laplace transformation directly, even for a single harmonic oscillator: 35 with the vibrational quantum number n, the vibrational energy ε i = hϑ i , Planck's constant h, and the harmonic frequency ϑ i of mode i. In the following, the function f (n) describes the exponential function exp (−nε i · s). Note that the inverse Laplace transformation is only applicable to continuous functions and exclusively produces continuous functions. An implementation of the following algorithms can be found in a PYTHON library attached as supplementary material.

A. Harmonic oscillators
The above expression for the harmonic oscillator partition function is replaced with an asymptotic expansion using the Euler-Maclaurin formula: 36 with the domain [α, β], the Bernoulli numbers B k , the above introduced function f (n) and its derivatives f (k -1) (n). Note that for k = 0 the function f (n) has to be integrated rather than differentiated. The partition function for a single harmonic oscillator is represented by the following asymptotic expansion: For a system of N separable harmonic oscillators, the energy-level density is the inverse Laplace transform of the partition function of the complete system. The total vibrational partition function is given as the product of the harmonic oscillator partition functions of the single vibrational modes: The partition function of the complete system is simplified by using the Cauchy product formula for converting the product of sums to a single sum (cf. Appendix). In there, the combinatorial nature of the complete system is condensed to the coefficients Λ N,j . Using the Cauchy product formulation of the separable harmonic oscillator partition function allows to solve the inverse Laplace transformation and to obtain the exact analytic energy-level density ρ N ϑ (E) according to: wherein the coefficients Λ N,j are defined in Appendix. Any s with an exponent bigger than or equal to zero will transform into the Dirac delta function δ(E) or a derivative of the delta function. Thus, the resulting expression for the energy-level density ρ N ϑ (E) is split into polynomial energy terms and terms which are related to the Dirac delta function: with the Gamma-function Γ(n). The second sum vanishes since the delta function and its derivatives are zero except for E = 0. Note that adding the translational and rotational partition functions increases the energy exponents in the energy-level density by 3/2 (NVT ) or 5/2 (NPT ) due to translation and by 2/2 (linear) or 3/2 (non-linear) due to rotation. The presented asymptotic expansion can be interpreted as the different 'levels' of energy exchange between the separate modes. The coefficients of this expansion are generated in a combinatorial manner. While the coefficient for the first term, the high-temperature limit, is simply the product of the coefficients of the separate modes, the coefficient of the second highest energypower is already the sum of different combinations. The number of combinations to be considered in a coefficient increases with decreasing energy-power. This agrees with the fact, that distributing many small energy amounts has more combinations compared to distributing few large energy amounts.

B. Internal rotations
When solving the nuclear Schrödinger equation for a general internal rotation potential energy profile, the resulting eigenvalues cannot be expressed analytically. Thus, the partition function of the respective mode cannot be represented by a closed analytic formula. Nevertheless, modeling the transition from harmonic to free motions is essential for obtaining accurate energy-level densities.
A simple model for converting the harmonic partition function Q ϑ,i to the free rotor partition function Q FR,i is given by Chuang and Truhlar via formulation of an approximate anharmonic partition function Q anh,i : 37 in which the ratio of the free-rotor (FR) and the high-temperature (HT) limit harmonic partition functions are used to judge whether the respective rocking mode is more like a free rotor or harmonic oscillator. The free rotor is described by the classical one-dimensional rotational partition function with the internal rotation symmetry number σ i and the internal rotational temperature Θ i . The harmonic oscillator high-temperature limit is given by the classical harmonic oscillator partition function Q HT,i . The partition function for a single harmonic oscillator with the above approximate hindered rotor corrections is represented by the following expansion of the hyperbolic tangent: wherein is the ratio of the temperature-independent parts of free rotor and classical harmonic oscillator partition functions.
For a system of N separable harmonic oscillators with M approximate hindered rotor corrections, the energy-level density is the inverse Laplace transform of the partition function of the complete system. The total vibrational partition function is given as the product of the anharmonic vibrational partition functions of the single vibrational modes: Again, the partition function of the complete system is simplified by using the Cauchy product formula for converting the product of sums to a single sum (analogous to Appendix). In there, the combinatorial nature of the separable hindered rotor corrections is condensed to the coefficients Φ N,j . Using the Cauchy product formulation of the separable anharmonic vibrational partition function allows to solve the inverse Laplace transformation and to obtain the energy-level density ρ N,M anh (E) according to: wherein the coefficients Φ M,j−k are defined analogous to the coefficients Λ N,k (cf. Appendix). Again, any s with an exponent bigger than or equal to zero will transform into the Dirac delta function δ(E) or a derivative of the delta function. The finite series containing the energy polynomials and the infinite series containing the Dirac delta functions are given by: Note that the Dirac delta function and its derivatives in the second sum vanish except for E = 0, and that the translational and rotational partition functions increase the energy exponents by 3/2 (NVT ) or 5/2 (NPT ) and 2/2 (linear) or 3/2 (non-linear), respectively.

III. VALIDATION FOR RRHO SYSTEMS
For the energy-level density of a single harmonic oscillator, an analytic formula is given by Müller-Kirsten. 38 In there, the energy-level density is defined as the reciprocal of the derivative of the harmonic oscillator energy levels E(n) = ε 1 (n + 1/2) with respect to n. This yields an energy-level density of 1/ε 1 , which equals Equation 8 for a single harmonic oscillator. Further, this result equals the so-called high-temperature approximation, which is exact for a single harmonic oscillator.
The analytic energy-level density formula for separable harmonic oscillators presented in Section A is validated against Beyer-Swinehart energy-level evaluations. Here, the well-studied CO+ȮH toḢ+CO 2 reaction network is used as a validation case. The present comparison of energylevel densities of separable harmonic oscillators is based on the harmonic frequencies reported in the study of Li et al. 30 A total of six stationary points are picked from the Li et al. 30  harmonic frequencies each. Each of these harmonic frequencies is assumed to be well-approximated by the harmonic oscillator model (ϑ min = 400.4 cm −1 ). The Beyer-Swinehart algorithm yields exact degeneracies which can be easily summed up to the exact number of states. In contrast, energy-level densities cannot be computed exactly via the Beyer-Swinehart algorithm, due to the continuous nature of energy-level densities. 22 In turn, integrating the present analytic energy-level density yields a continuous rather than discrete expression for the number of states. Thus, a fair comparison between energy-level densities or number of states is not possible. In Figure 1, the exact Beyer-Swinehart number of states N BS (E) is used as reference for validating the integrated and integer-rounded analytic energy-level density, i.e. the integer-rounded analytic number of states N round anly (E). The number of states ratio N round anly (E)/N BS (E) is not unity over the whole energy regime. This is due to comparing the exact Beyer-Swinehart number of states N BS (E), which is discrete, to the analytic integer-rounded number of states N round anly (E), as mentioned previously. For the present example systems, the deviation between N round anly (E) and N BS (E) amounts to a factor of two for low-energy, i.e. an absolute deviation of one. For high-energy, however, the absolute deviation increases to 16 at most.
The oscillations observed for the ratio N round anly (E)/N BS (E) correspond to the accessible energy levels, thus to the multiples of harmonic frequencies at low-energy. With increasing energy, the amplitude of the oscillating number of states ratio decreases, converging to unity. This is due to the increasing number of states-values, for which the relative deviations become less pronounced. In order to illustrate the qualitative agreement between exact and analytic number of states, the results for the stable stationary point trans-HOĊO are shown in Figure 2. Note that the exact Beyer-Swinehart number of states is represented by the thin solid lines and that the continuous analytic number of states (N cont anly (E)) are represented by the thick red dashed lines. Here, the analytic number of states are not integer-rounded, but directly compared to the Beyer-Swinehart results.
As can be seen in Figure 2(a), the curves representing the number of states computed via the Beyer-Swinehart and the present analytic formula perfectly agree. The high-temperature approximation in turn significantly overestimates the number of states with increasing energy. When zooming in, the differences between N BS (E) and N cont anly (E) becomes obvious (cf. Figure 2(b)). The continuous analytic number of states, however, is roughly the average of the Beyer-Swinehart number of states. Again, the high-temperature approximation shows significant deviation to the reference results, but this time underestimating them.

IV. VALIDATION FOR ONE-DIMENSIONAL HINDERED ROTOR SYSTEMS
The analytic energy-level density formula for separable harmonic oscillators with approximate hindered rotor corrections presented in Section B is validated against Beyer-Swinehart energy-level evaluations. Here, the well-studied H 2 O 2 system is used as a validation case. This species has a single anharmonic torsional motion around the O−O axis, which was studied in great detail. [39][40][41][42] Here, the energy levels reported by Ellingson et al. 41 are combined with the free rotor energy levels 43 above the torsional barrier. The resulting energy levels and the harmonic oscillator energy levels of H 2 O 2 are used as input for the modified Beyer-Swinehart algorithm. 6 In Figure 3, the exact Beyer-Swinehart number of states N BS (E) is used as reference for validating the integrated and integer-rounded analytic energy-level density, i.e. the integer-rounded analytic number of states N round anly (E). The first two figures compare the exact and analytic number of states of the H 2 O 2 torsional motion only, excluding contributions from other modes. Firstly, the H 2 O 2 torsional motion is treated as a free rotor (cf. Figure 3(a)) and secondly, the H 2 O 2 torsional motion is treated as a hindered rotor (cf. Figure 3(b)). Finally, the torsional and the harmonic modes are treated in a combined manner (cf. Figure 3(c)).
Although using the model of Chuang and Truhlar 37 for the transition from the harmonic oscillator to the free rotor regime, N BS (E) and N round anly (E) are equal over the whole energy-regime for the pure free rotor reference results (cf. Figure 3(a)). When considering the continuous analytic number of states N cont anly (E), however, the deviation to the exact Beyer-Swinehart number of states becomes obvious and is as expected (cf. dashed line in Figure 3(a)). Since the model of Chuang and Truhlar 37 is for all three vibrational regimes: Harmonic oscillator, hindered rotor, and free rotor, the exact Beyer-Swinehart free rotor number of states must exceed N round anly (E) for low-energies, and both converge to the free rotor number of states with increasing energy. For the one-dimensional hindered rotor in turn, the deviation between N BS (E) and N round anly (E) is obvious (cf. Figure 3(b)). While at energies up to roughly 2500 cm −1 the energy levels of the hindered rotor are dominating N BS (E), higher energies are dominated by the free rotor energy levels. At high energies, the deviation between N BS (E) and N round anly (E) is much larger than for the pure free rotor (cf. Figure 3(a)). This is presumably due to the energy-level density contributions of the hindered rotor regime, affecting the convergence even at high energies. With increasing energy, however, N BS (E) and N round anly (E) are converging to the free rotor number of states. When combining the contributions from the H 2 O 2 hindered rotor and the harmonic oscillators, the vibrational number of states of H 2 O 2 are obtained for the complete system (cf. Figure 3(c)). In there, the exponentially decaying deviations between N BS (E) and N round anly (E), and N cont anly (E), reported in Figures 3(a) and 3(b), merge with the oscillatory deviations of the pure harmonic oscillator system of H 2 O 2 (cf. Figures 1 as an example). For high energy, N round anly (E) does not match N BS (E) due to the same reason discussed for the hindered rotor in Figure 3(b): The hindered rotor energy levels affect the number of states in the free rotor energy regime. Again, N BS (E) and N round anly (E) are converging to the same number of states.

V. CONCLUSION
The present work provided an analytic formula for vibrational energy-level densities, satisfying their continuous nature. This formula was derived using the inverse Laplace transformation of the asymptotically expanded vibrational partition function of separable harmonic oscillators. Further, the formula was extended to approximate hindered rotor energy-level densities, based on a literature model for the transition from the harmonic oscillator to the free rotor partition functions. Finally, the analytic formulas were validated against the well-studied rigid rotor harmonic oscillator CO +ȮH system, and against the well-studied one-dimensional hindered rotor H 2 O 2 system. For any validation case, the reference results were obtained via the exact Beyer-Swinehart algorithm.
For the RRHO case, the analytic solution was derived without any approximation. The comparison between the Beyer-Swinehart and analytic solutions showed oscillatory deviations, resulting from the discretization error made when integrating the energy-level densities. Here, the formula for smooth ρ(E) is supplied, whose Laplace transform is the RRHO partition function Q(T ), which is expected to be convinient in applications where a smooth continuous ρ(E) is desired (such as for RRKM and Master Equation calculations). For the hindered rotor case, the analytic solution is based on an approximation, thus cannot be exact by definition. The deviations introduced by this approximation were found to mainly originate from an inaccurate approximation of the hindered rotor energy levels. For the H 2 O 2 test system, deviations above a factor of two were observed below 2000 cm −1 , thus in the hindered rotor regime of H 2 O 2 . With increasing energy, however, this deviation vanished and the solution slowly converged to the free rotor. As mentioned in the introductory discussion on the representation of kinetic models, both the thermochemistry and the kinetics of a reaction network can be represented by energy-level densities. The present closed analytic energy-level density formulation satisfies the continuous nature of energylevel densities and allows to store a set of parameters, namely Λ N,j and Ω N,j , rather than actual energy-level densities. These parameters are directly computed from molecular properties in the present study. In general, however, fitting these parameters to experimental data or estimating them with not yet existing rules is possible in principle. Using the present analytic continuous energy-level densities would satisfy the need for continuous functions in RRKM calculations and avoid the number of states-based finite difference approximation of energy-level densities.
Here, the closed analytic energy-level density formulation was derived using the inverse Laplace transform of the canonical partition function, i.e. the Boltzmann-weighted energy levels. The Beyer-Swinehart algorithm in turn is based on counting microscopic energy states, thus is independent from the energy-distribution. Since the present and the Beyer-Swinehart number of states were shown to agree, the closed analytic energy-level density formulation could be used to predict non-Boltzmann thermochemistry when being weighted with the non-Boltzmann energy-distribution of the respective species.
As touched in the introductory discussion, a closed analytic energy-level density expression, such as presented, can be used to study zero-dimensional chemistry.

APPENDIX: ASYMPTOTIC EXPANSION
The product of sums in Equation 6 can be expressed as a single sum using the Cauchy product formula: 44 For two general infinite series with elements a k and b k the product is formulated by an infinite series with elements c j : Using this formulation for the product of partition functions gives the following series: with parameters Λ N,j being completely independent from s. The parameters are defined recursively, starting with Λ 0,j = 1 and Λ 1,j = λ 1,j , with All subsequent Λ i,j 's are computed using the following expression: for i ranging up to N.

SUPPLEMENTARY MATERIAL
See supplementary material (rho E.py) for a python library of the present algorithms.

ACKNOWLEDGMENTS
Financial support from the Deutsche Forschungsgemeinschaft (German Research Association) through grant GSC 111 is gratefully acknowledged.