## Abstract

We discuss recent progress in the calculation and identification of rotation–vibrational states of at intermediate energies up to 13 000 cm^{−1}. Our calculations are based on the potential energy surface of Cencek *et al*. which is of sub-microhartree accuracy. As this surface includes diagonal adiabatic and relativistic corrections to the fixed nuclei electronic energies, the remaining discrepancies between our calculated and experimental data should be due to the neglect of non-adiabatic coupling to excited electronic states in the calculations. To account for this, our calculated energy values were adjusted empirically by a simple correction formula. Based on our understanding of the adiabatic approximation, we suggest two new approaches to account for the off-diagonal adiabatic correction, which should work; however, they have not been tested yet for . Theoretical predictions made for the above-barrier energy region of recent experimental interest are accurate to 0.35 cm^{−1} or better.

## 1. Introduction

There appears to be no need to explicate here the importance of and its rotation–vibrational states in astrophysics. The role it plays in planetary atmospheres, dense and diffuse interstellar clouds or even in the early development of our Universe will surely be discussed and presented in great detail in this issue by those who are experts in these fields. The focus of our presentation will be the description and analysis of an accurate and precise calculation of the rotation–vibrational states of , the simplest triatomic molecule. This is done with the aim to extrapolate guidelines from our experience, which should be useful for similar *ab initio* calculations of the rotation–vibrational states of other triatomic or even more complex molecules in the future.

We apologize for not reviewing the many other earlier *ab initio* calculations of the rotation–vibrational states of , as this has been done in the past, and by other authors of this issue, e.g. Kutzelnigg & Jaquet (2006) and Tennyson *et al*. (2006). At the same time, we are grateful for the highly accurate, sub-microhartree potential energy surface of , which resulted from an extensive electronic structure calculation, with the imaginative explicit inclusion of electron correlation by Cencek *et al*. (1998). Without these data for the potential energy surface, including the relativistic and diagonal adiabatic corrections, the results obtained would have not been possible. Equally valuable for the accuracy we could achieve for many of the spectroscopically still unassigned states of was the interaction with the researchers of the group of T. Oka, where the spectrum of is and has been studied expertly and extensively with high precision, e.g. Lindsay & McCall (2001), Gottfried *et al*. (2003) and Gottfried (2006).

The rotation–vibrational term values were computed using the potential surface of Cencek *et al*. (1998). To this end, the potential surface is expressed using the ‘democratic’ hyperspherical coordinates *ρ* for the hyper-radius and the two angles *θ* and *ϕ*. The rovibrational wave functions are expanded in terms of the hyperspherical harmonics . This ensures that the rotation–vibrational coupling, affected via the kinetic energy operator, is taken care of both analytically and exactly. The resulting set of coupled differential equations in the remaining coordinate, the hyper-radius *ρ*, were integrated numerically.

Since the deviation of the computed term values should be of the order of fractions of wavenumbers, i.e. one part in 10^{6}, to be useful in the interpretation of the observed spectral lines, it is mandatory to control accuracy and precision of all theoretical and computational aspects carefully. An error analysis of the systematic and random errors is called for.

To be sure, the fundamental constants needed here, the value of the proton mass, *M*_{p}=1836.15267261(85)*m*_{e}, and the energy conversion factor, *E*_{h}=219474.6313705(15) cm^{−1} are known with sufficient accuracy, generally with nine significant figures (Mohr & Taylor 2005). The relative precision of the computed points of the potential surface and the diagonal adiabatic corrections is high enough to yield potentially a deviation of the computed rotation–vibrational term values of a few hundreds of wavenumbers. The global fit of the potential surface, however, introduces uncertainties of the order of 0.1 cm^{−1}.

Most important, however, is a careful scrutiny of the adiabatic separation between the electronic and nuclear coordinates. For one, the discussion whether atomic or nuclear masses should be used in the computation of the rovibrational states appears to linger on. This would introduce differences in the fourth significant figure of the computed term values. For two, the detailed computation of the off-diagonal adiabatic corrections for appears to be prohibitive. Thus, the significance of these corrections must be estimated and empirical adjustments should be made, provided they are judged reliable in the light of the analysis.

Another important aspect, if contact is to be made between theory and experiment, is the establishment of the correspondence of the computed term values with their exact quantum numbers to the experimentally determined term value, labelled with the approximate quantum numbers as used in spectroscopic work. This association is relatively straightforward for the lower levels; however, for term values above 10 000 cm^{−1}, where the molecule is quite floppy, the assignments can become ambiguous.

## 2. Transformation of the Hamiltonian

The non-relativistic Hamiltonian of a molecule with *N* nuclei and *n* electrons in field-free space is readily written down in terms of the Cartesian coordinates of the *N* nuclei and *n* electrons , in the laboratory frame. Here, upper case letters denote nuclear and lower case letters electronic coordinates.(2.1)where the potential *V*(** Q**,

**) depends only on the inter-particle distances. The kinetic energy operators with respect to Cartesian coordinates of the laboratory fixed system are defined as(2.2)**

*q*In order to separate the motion of the centre of mass(2.3)with the mass of the nuclei and the total mass given as *M*=*M*_{N}+*n*.*m*_{e}, from the internal molecular motions, it is appropriate to define molecular-centred coordinates. These types of transformations are discussed in detail in the text book of Zülicke (1984). Of the many possible reference points as origin for the molecular-centred particle vectors, we prefer here the centre of nuclear mass,(2.4)to yield the molecular-centred vectors(2.5)

This choice for the origin has the advantage to yield the simple Hamiltonian(2.6)where the translational motion of the centre of mass of the molecule is well separated, and where there are no coupling terms between the momenta of the nuclei and those of the electrons. Obviously, the inter-particle potential, depending only on the distances between the particles, is unchanged by this transformation, i.e. *V*(** Q**,

**)=**

*q**V*(

**,**

*R***). The kinetic energy operators for the motion of the centre of mass , the motion of nuclei , the motion of electrons and the operator for the mass polarization involving electronic coordinates are obtained as(2.7)**

*r*(2.8)

(2.9)

(2.10)

Actually the form of the kinetic energy operator given here for the nuclear motion, (equation (2.8)), is rather symbolic, just to indicate that it depends only on 3(*N*−1) independent degrees of freedom. However, it is clear from the form of and that nuclear masses must be used in the *ab initio* calculation of the nuclear motion. In general, one attempts to find coordinates for the nuclear motion, which are more suitable for the system than the 3(*N*−1) Cartesian coordinates. For diatomic as well as for triatomic molecules, this can be achieved readily. For larger polyatomic systems, this is still an art, documented in a large number of recent publications also (e.g. Sutcliffe 2000 and references therein). It is not even clear yet, whether it is wise for *ab initio* calculations to follow the route chosen by spectroscopists, to separate first the rotational motion of the nuclear framework, by transforming to a molecular fixed coordinate system. For this transformation, three principal axes of inertia need to be chosen, which can be achieved only approximately, e.g. by using the Eckart conditions. For floppy molecules, this approximation becomes unsatisfactory.

## 3. The hyperspherical harmonics approach

Here, for , we transform the nuclear Hamiltonian to ‘democratic’ hyperspherical coordinates, which are derived from the corresponding Jacobi coordinates. The definition of the hyperspherical coordinates as well as their interpretation has already been discussed in detail (Johnson 1980, 1983*a*,*b*; Hinze *et al*. 1995). In the hyperspherical approach, the internal Hamiltonian reads(3.1)where *Λ*^{2}(*Ω*) is the grand angular momentum operator, which depends on the three Euler and two hyperspherical angles *Ω*≡{*θ*,*ϕ*,*α*,*β*,*γ*} and the three-particle reduced mass is given by , where *M*_{N} is the total mass of the three nuclei. The potential *U*(*ρ*,*θ*,*ϕ*) is a function of the hyper-radius *ρ* and the two hyperspherical angles, which determine the size and the shape of the three-particle system. Using *u*=cos *θ*, one obtains for the grand angular momentum operator(3.2)where *J*_{x}, *J*_{y} and *J*_{z} are the components of the total angular momentum operator in the body fixed frame. Hyperspherical harmonics (Avery 1989), the eigenfunctions of the grand angular momentum operator(3.3)can be represented as(3.4)(3.5)where are normalized Wigner rotation matrix elements in the convention of Zare (1988). The explicit form of the hyperspherical harmonics depends on the particular parametrization chosen for the coordinates. Wolniewicz (1989) has derived an algorithm for their evaluation and has given the ranges of the hyperspherical quantum numbers *K*, *ν*, *s*. Of the hyperspherical harmonics defined in equation (3.4), linear combinations(3.6)may be formed that are simultaneous eigenfunctions of the grand angular momentum operator and of the operators of the three-particle permutation-inversion group *S*_{3}×*I*. In the above equation, *ϵ*=e^{−3πiν/4} and *Γ* are one of the irreducible representations of *S*_{3}×*I*. *α* is a composite index, *α*={*K*,*ν*,*s*,*p*}, which collects the indices of the hyperspherical functions compatible with *Γ*. Since *Γ* is a strictly conserved quantity, the symmetry adaptation leads to a significant reduction in the size of the Hamiltonian matrix because calculations can be performed independently for each irreducible representation. The rovibrational states of are the solutions of the six-dimensional Schrödinger equation (Wolniewicz *et al*. 1993; Wolniewicz & Hinze 1994; Alijah *et al*. 1995*a*)(3.7)

The product ansatz(3.8)transforms equation (3.7) into a set of coupled equations in the hyper-radius *ρ*,(3.9)Here, is the column vector of the *m* channel functions and **1** is the *m*×*m* unit matrix. The coupling is affected by the matrix elements of the potential only(3.10)

As in our previous work, we have used a contracted basis rather than the primitive basis in order to reduce the number of coupled equations. We then obtain a set of coupled equations analogous to equation (3.9) with channel functions and and *U*^{JΓ}(*ρ*) replaced by(3.11)

This set of equations was solved as discussed previously (Wolniewicz & Hinze 1994). The advantages of using the expansion of the eigenfunctions for the nuclear motion in terms of the hyperspherical harmonics (equation (3.6)) are twofold: (i) the rotation–vibrational coupling is accounted for both analytically and exactly and (ii) as the hyperspherical harmonics obey the correct boundary conditions as the hyperangle *θ* approaches 90°, the collinear structure of the triatomic, there is no problem in the transition of the triatomic to a linear structure.

## 4. Exact hyperspherical versus spectroscopic quantum numbers

The rotation–vibrational states of can be classified rigorously by the total angular momentum, *F*, its space-fixed projection, *F*_{z}, which is however of no importance in the absence of external fields, and by their symmetry index *Γ* with respect to the three-particle permutation-inversion group, *S*_{3}×*I*. This group is isomorphic with the point group *D*_{3h} and has six irreducible representations, *A*_{1}, *A*_{2} and *E*, either ‘prime’ or ‘double-prime’ according to their parity. For , hyperfine interactions can be neglected so that the nuclear spin, *I*, and the remaining angular momentum, *J*, are conserved separately. Hence, the quantum numbers (*I*,*J*,*Γ*,*n*), where *n* is a counting index, completely describe the rovibrational states. In our hyperspherical method, which makes full use of the permutation-inversion symmetry, we readily obtain such a classification. The symmetry index *Γ* is given in terms of the hyperspherical quantum numbers *ν* and *p* as shown in table 1.

To be sure, not all calculated states are permitted by the Pauli principle as this requires the total wave function to be antisymmetric with respect to a permutation of two protons. Thus, only states with rovibrational symmetry *A*_{2} and *E* exist, the former being combined with the totally symmetric quartet nuclear spin function and the latter with the degenerate doublet function.

As an alternative to the rigorous quantum numbers, approximate spectroscopic quantum numbers may be employed to characterize the rovibrational states. Such quantum numbers have been discussed in detail by Watson (1984). Following his work, a rovibrational state is characterized by (*I*,*v*_{1},*v*_{2},*J*,*G*,*U*,*s*), where *v*_{1} and *v*_{2} are the quantum numbers of the symmetric stretch and the degenerate vibration, respectively. *G* is defined as , with *k* denoting the internal projection of the angular momentum *J* and denoting the vibrational angular momentum associated with the degenerate vibration. As Hougen (1962) has shown, *k* and are not conserved separately. In spectroscopic representation, the permutation-inversion symmetry of a rovibrational state is completely determined by the quantum numbers *G*, *v*_{2} and *s*, see table 1. The latter quantum number *s* is needed in the case of *G*≡0(mod3) to distinguish states of symmetry *A*_{1} and *A*_{2}.

For the assignment of spectroscopic quantum numbers to the calculated ‘hyperspherical’ states, a semi-automatic procedure has been developed (Alijah *et al*. 1995*b*), which is based on the transformation properties shown in table 1 and on a comparison of rotational progressions in vibrational states that differ only in the symmetric stretch quantum number *v*_{1}. The algorithm has been applied to recently (Schiffels *et al*. 2003*b*).

## 5. The adiabatic coordinate separation

It is necessary for the practical solution of the Schrödinger equation using *H*_{Tot} (equation (2.6)) to reduce the number of coordinates such that the problems for the electronic and nuclear motion can be solved separately. To this end, the adiabatic coordinate separation(5.1)is introduced (Born 1952; Born & Huang 1954), where the electronic state functions *Φ*_{n}(** R**;

**) are solutions of the electronic Schrödinger equation(5.2)**

*r*It should be noted that an adiabatic coordinate separation is always better than an analogous simple product ansatz, irrespective of what the relative masses are. Here, where the nuclear masses are 2000 times larger than the electronic masses, the separation is especially good, yielding a rapidly converging series, such that even the single leading term is a good approximation.

The expansion of the total wave function (equation (5.1)) yields a set of coupled equations to be solved for the nuclear motion(5.3)

In general, the operators *C*_{nn′} and *D*_{nn′} will be off-diagonal in the basis of the electronic eigenfunctions(5.4)

(5.5)

The first derivative coupling matrix defined by(5.6)has the property because the electronic eigenfunctions are normalizable. The diagonal elements are thus purely imaginary quantities, and *W*_{nn,I}(** R**)=0 immediately follows for real electronic wave functions. In the one state (simple adiabatic) approximation and for real electronic wave functions, we thus obtain(5.7)

In order to describe the nuclear motion in the one-state approximation, the fixed nuclei potential *U*_{n}(** R**) and the diagonal adiabatic correction

*C*

_{nn}(

**) should be known. The coupling of the nuclear motion to different electronic states via equation (5.4) and (5.5) is called the off-diagonal adiabatic or non-adiabatic correction. It can be treated perturbational, especially if equation (5.3) is written in matrix form, i.e. using a basis function expansion for**

*R**Θ*

_{nα}(

**). To illustrate this, let us focus on the electronic groundstate**

*R**n*=1 and on the first-order off-diagonal coupling elements only. With this, we obtain equation (5.3) in matrix form(5.8)with and . Using partitioning yields(5.9)

With this perturbation approximation, there is no need to compute the rovibrational wave function components within the excited electronic states, only the matrix elements with the wave functions of the excited electronic states need to be computed.

Another more simple, though empirical approximation is based on our understanding of the fixed nuclei approximation (Hinze *et al*. 1998). With this the electrons are rigidly attached to the nuclei and their masses are hidden in the electronic potential *U*. Thus, they are over-correlated with the nuclear motion, while in fact they should lag behind if the nuclei are accelerated. Now consider Newton's equation(5.10)and transform it to obtain(5.11)This yields(5.12)

Here, we have multiplied by the electronic masses *nm*_{e}, where *n* is the number of electrons, and the force by a distance *R*−*R*_{e}, where *R*_{e} corresponds to the minimum of the potential, to get the dimension of an energy for the expectation value, to be taken over rovibrational wave functions. Note that Δ*E* goes to 0 if the prefactor *nm*_{e}/*M*_{N} goes to 0, because there should be no correction in the limits of zero electron or infinite nuclear masses. Equally, the integral goes to 0 at the energy of dissociation because the derivative of the potential becomes 0 for the large distances, where the vibrational wave function is large. In the above equation, *U*_{J}(*R*) denotes the effective potential, which is the sum of the electronic energy and the *J* dependent centrifugal term. Again, the nuclear masses *M*_{N} and the distance *R*−*R*_{e} are symbolic here: they will depend on the coordinate system chosen for the description of the nuclear motion.

This formula has been tested for H_{2} for which exact numerical data (Wolniewicz 1995) are available. Figure 1 shows the vibrational non-adiabatic correction as approximated by the negative of(5.13)together with Wolniewicz's results. Here, the reduced nuclear mass *μ*=*M*_{1}*M*_{2}/(*M*_{1}+*M*_{2})=918.076336*m*_{e} has been used. The numerical value of the dimensionless parameter is *c*=0.113430±3.11×10^{−3}. As can be inferred from figure 1, the empirical formula (5.13) describes the general form of the correction appropriately and yields results accurate to better than 0.5 cm^{−1}, i.e. better than 10% of the maximum adiabatic correction. We expect that this approach, with the present numerical value of the parameter, can be transferred to . There it should work particularly well as the adiabatic correction is much smaller than in H_{2}, of the order of only one wavenumber.

## 6. Results

Using the method of hyperspherical harmonics and nuclear masses, we have calculated and analysed the rovibrational states, for *J*≤10, up to the energy region of around 13 000 cm^{−1} (Schiffels *et al*. 2003*a*,*b*). The highest states are thus located well above the barrier to linearity at *ca* 9000 cm^{−1}. States above the barrier have been observed in recent experiments by the Oka group (Gottfried *et al*. 2003), and the ongoing experiments in this group are probing even higher excited states, see Gottfried (2006). Our theoretical work was performed to assist in the assignment of observed states. In turn, comparison of our calculated data with the very precise experimental data provides a rigorous assessment of the accuracy achieved in theoretical work.

We begin our analysis with the rovibrational states below the barrier to linearity. Comparison of our calculated energies, which are based on the Cencck–Rychlewski–Jaquet–Kutzelniyg (CRJK) potential energy surface (Cencek *et al*. 1998; Jaquet *et al*. 1998), to the comprehensive compilation of experimental data by Lindsay & McCall (2001) reveals deviations of up to 1.2 cm^{−1}. Since the *ab initio* electronic energies on which the potential energy surface is based are accurate to *ca* 0.02 cm^{−1}, as the diagonal adiabatic correction and relativistic effects are included, the discrepancies between theory and experiment must be due to the neglected off-diagonal adiabatic contributions, i.e. due to coupling effects to higher electronic states.

The analysis of the differences between the observed and calculated rovibrational energies shows strong vibrational, but relatively small rotational, shifts. Such a behaviour has been observed before, in our work on the isotopologues (Alijah *et al*. 1995*b*), H_{2}D^{+} (Alijah *et al*. 1995*a*) and D_{2}H^{+} (Alijah & Beuger 1996), where it led to empirical corrections of the band origins. In our 2003 work (Schiffels *et al*. 2003*a*,*b*), we considered both vibrational and rotational non-adiabatic corrections. Following this work, the correction is expressed analytically as(6.1)

The coefficients, which depend explicitly on the vibrational quantum numbers [** v**] were determined by least squares fitting of the differences between calculated and observed data. When this correction is subtracted from the calculated energy values, an estimate is obtained for the rovibrational energy including non-adiabatic effects(6.2)

A close agreement of within 0.1 cm^{−1} of the corrected energies with experimental data was found. This accuracy is consistent with the accuracy of the analytic formula for the representation of the *ab initio* electronic energies. A disadvantage of expression (6.2) is that it cannot be applied directly. Owing to the dependence of the coefficients on the vibrational quantum numbers, a prior assignment of the *ab initio* term values to spectroscopic term values would be required. However, it may be improved by the following measures. Firstly, since the coefficients of the rotational correction were found to be almost the same for all bands, they may be replaced by band-independent average values and . Secondly, the correction of the band origins may be expressed as a linear function of the energy, . With this, we obtain a more general correction formula, allowing extrapolation to higher energies(6.3)

The accuracy obtained by applying the above correction formula is expected to be 0.2 cm^{−1}. Even though is more general than and with little loss of accuracy, it can be used only if the vibrational (to identify the band origin) and rotational quantum numbers are known. Frequently, this is not the case, and moreover, the spectroscopic quantum numbers become less meaningful at high energies anyway. As an alternative to formula (6.3), we suggested the following correction formula,(6.4)in which the *b*_{1} scaling is applied directly to a calculated rovibrational state rather than just to the band origins. With this, we expect to overcorrect in particular rotationally highly excited states. At not too high values of the rotational quantum numbers, the predictions made by this rather simple formula should still be fairly accurate.

Since the applicability of our recommended correction formula (6.3) relies on the correct identification of calculated states in terms of spectroscopic quantum numbers, we attempted to assign as many such states as possible. To this end, the semi-automatic assignment algorithm was employed first in a test to recover the known assignments (Lindsay & McCall 2001) of the rovibrational states below the barrier. As the proposed assignments proved reliable, we extended our analysis to higher energy states up to 5*ν*_{2}. Some heavily perturbed states, however, could not be assigned as in these cases, the spectroscopic quantum numbers do no longer hold. The energy values of such states were adjusted with help of the more approximate correction formula (6.4).

Our calculated energy levels, corrected empirically mostly by formula (6.3) (Schiffels *et al*. 2003*b*), and their proposed assignments have proved helpful in the interpretation of recent experimental data (Gottfried 2006) as they led to the most accurate calculated transition frequencies. In the analysis, it became obvious, however, that the energy shifts produced by the correction formula for the states in the 13 000 cm^{−1} region are too big. The average error in our data, based on 34 experimentally observed states, was found to be −0.35 cm^{−1} compared to 0.80 cm^{−1} for the uncorrected case. Following the discussion of the adiabatic correction in the preceding section, this over-estimation of the non-adiabatic effects can be understood readily. Our extrapolation procedure assumes a linear dependence of the vibrational correction on the energy, but we know that this is not strictly true as at the dissociation limit such a correction must approach zero, see figure 1. At 13 000 cm^{−1}, i.e. after linear extrapolation over not less than 4000 cm^{−1}, nonlinear behaviour of the non-adiabatic effects as a function of the total energy becomes noticeable.

## 7. Outlook

Thus, our correction formulae (6.3) and (6.4) should be adapted to account for this. The new, very precise experimental data will form the basis for such an adaptation. With this, we expect to make accurate predictions of the rovibrational levels at even higher energies.

## Footnotes

One contribution of 26 to a Discussion Meeting Issue ‘Physics, chemistry and astronomy of ’.

- © 2006 The Royal Society