Michael Rey
email: michael.rey@univ-reims.fr

Novel methodology for systematically constructing global effective models from ab initio-based surfaces: a new insight into high-resolution molecular spectra analysis

Keywords: Block-diagonalization procedure, ab initio global effective model, Symmetry, Irreducible tensor operators, Infrared spectra, Analysis

In this paper, a novel methodology is presented for the construction of ab initio effective rotation-vibration spectroscopic models from potential energy and dipole moment surfaces. Non-empirical effective Hamiltonians are obtained via the blockdiagonalization of selected variationally-computed eigenvector matrices. For the first time, the derivation of an effective dipole moment is carried out in a systematic way.

This general approach can be implemented quite easily in most of the variational computer codes and turns out to be a clear alternative to the rather involved Van Vleck perturbation method. Symmetry is exploited at all stages to translate firstprinciples calculations into a set of spectroscopic parameters to be further refined on experiment. We demonstrate on H 2 CO, PH 3 , CH 4 , C 2 H 4 and SF 6 that the proposed effective model can provide crucial information to spectroscopists within a very short time compared to empirical spectroscopic models. This approach brings a new insight into high-resolution spectra analysis of polyatomic molecules and will be also of great help in the modelling of hot atmospheres where completeness is important.

I. INTRODUCTION

Precise knowledge of high-energy molecular states and absorption spectra is of primary importance because it gives access to the determination of the physical properties of various planetary objects 1,2 and clearly demonstrates the necessity of having consistent line-byline molecular databases (e.g. like HITRAN 3 or GEISA 4 for the modelling of the Earth's atmosphere). The interpretation of strong spectral features requires the use of sophisticated and robust theoretical models (i) for an accurate quantum-mechanical description of highlyexcited molecular states and (ii) for the prediction of line intensities for reliable opacity calculations. Two approaches are commonly used in spectroscopy for the modelling of rotation-vibration spectra:

Effective models. For the analysis of absorption spectra of molecules whose energy levels are organized as small groups of strong interacting vibrational levels, called polyads, the idea of introducing effective, or phenomenological, models (EM) is now well established [START_REF] Aliev | Higher-order effects in the vibration-rotation spectra of semirigid molecules[END_REF][START_REF] Papousek | Molecular vibrational-rotational spectra[END_REF][START_REF] Champion | Spherical Top Spectra[END_REF] . Such models are defined by a set of empirical parameters, either directly fitted to experiment or fixed manually by the user. By definition, each polyad is formed by groups of nearly degenerate vibrational states characterized by the same number, called polyad number

P = c 1 v 1 + c 2 v 2 + • • • + c Nm v Nm (c i ∈ R >0 ), (1)
where v i are the vibrational quantum numbers and N m the number of vibrational modes.

In a basis of the type | γ; J, C P , every block of the Hamiltonian matrix will be properly labelled by the total angular momentum J and by the symmetry C. As stated in Ref. [START_REF] O'brien | [END_REF] , the polyad numbers P are approximate good quantum numbers which will be also used to label every block of the effective Hamiltonian. γ denotes all other labels or quantum numbers (v i , etc). The first polyad P = 0, or simply denoted as P 0 , corresponds to the ground vibrational state. Though the second polyad does not necessarily equal P = 1, it will be denoted as P 1 , and so on. The full problem is thus divided into a series of much smaller problems making EMs. The success of the effective approach in spectra analyses will be partly governed by a convenient choice of the polyad vector c = (c 1 , c 2 , • • • , c Nm ) t . Krasnoshchekov and Stepanov 9 introduced a number Θ which is interpreted as a measure of proportionality between the polyad coefficients c i and the vibrational harmonic frequency ω i . Molecules like methane or phosphine exhibit a clear polyad structure (see Section IV), at least for the first vibrational states, while many others do not have such a regular structure.

initio effective Hamiltonian and dipole moment operators. This method combines the best of both worlds, namely (i) manipulation and diagonalization of block-diagonal matrices of small dimension for fast calculations, even for very high J values (> 100), with the possibility of refining some molecular parameters on experiment without much computational effort and (ii) completeness of the variational method. As a prerequisite, we start from a complete nuclear motion ab initio Hamiltonian, assumed known here, and compute the variational solutions of the stationary Schrödinger equation for low-J values only.

In this work, we focus on the study of semirigid molecules in their ground electronic state belonging to arbitrary Abelian or non-Abelian point groups. The motivations of this work as well as the methodology for the derivation of ab initio EMs are presented in Section II.

Computational details are given in Section III and the validation of the method is presented in Section IV on the calculation of rovibrational energy levels and infrared spectra of H 2 CO, PH 3 , CH 4 , C 2 H 4 and SF 6 . It is highly likely that this work will be of great help in current and futures analyses by accompanying spectroscopists in the assignment of rotationally resolved infrared spectra, even for complex polyads. As to the construction of an effective model for more flexible molecules, it will be briefly discussed in Conclusion.

II. CONSTRUCTION OF AB INITIO EFFECTIVE MODELS: MOTIVATIONS AND METHODOLOGY

A. Motivations

Undoubtedly, the effective Hamiltonians have greatly contributed to the "golden age" of the high-resolution molecular spectroscopy, but the current researches in the study of various planetary atmospheres require knowledge of increasingly complex molecular systems over wide wavenumber and temperature ranges 2 . Unfortunately, empirical EMs are beginning to reach their limits for studying molecules with complex rovibrational energy-level structures and for which the successive polyads contain many vibrational bands and numerous degeneracies and quasi-degeneracies. Modelling of the "dark states" which are not directly observable is one of the major obstacles in the empirical effective approach. Another challenging problem concerns the modelling of line intensities, both for cold and hot band transitions. Indeed, the density of states may rapidly increase for many molecules resulting in strongly congested spectra where almost no individual lines can be extracted easily due to the overlap of thousand experimental transitions, even at room temperature. This lack of information may lead to a poor determination of both resonance coupling and effective dipole moment parameters with possibly wrong intensity transfers between weak and strong lines.

It is evident that the complete, or even partial line-by-line analysis of very crowded spectral regions can take years or even decades using "traditional" EMs (see Section IV).

Extrapolation to high temperatures (say >1000 K) is another limiting factor because too many hot bands are still missing in the available EMs. Consequently, it is highly desirable to propose an alternative model which is capable of both dealing with complex polyads within a very short time and accompanying spectroscopists in the assignment of dense and rich rotationally resolved infrared spectra.

The aim of this paper is to build an ab initio EM whose the key features can be summarized as follows:

-This model must contain most of the resonance coupling terms up to a given polyad where the polyad scheme defined in Eq. ( 1) or (2) has to be conveniently chosen by the user. Contrary to empirical EMs where some vibrational bands in polyads are sometimes voluntary omitted to simplify calculations, our ab initio EM includes the major contributions for computing line positions and line intensities over a large spectral range.

-This model takes full advantage of the symmetry and provides to spectroscopists an initial set of physically meaningful "ab initio" effective parameters. Except for the zero order vibrational and the rotational constants, it is common to start an analysis from parameters initialized to zero whereas those we are able to provide are initialized to "good" values. An illustrative example on the determination of the ground state (P = 0) and dyad (P = 1) parameters of PH 3 will be given in Section IV.

-A part of the initial set of Hamiltonian parameters can be further optimized by fine tuning to experimental spectra in order to compute energy levels at the "spectroscopic accuracy", say of 10 -3 cm -1 which is for example the typical resolution of Fourier transform rotationally resolved spectra in the infrared. This procedure takes few seconds and is thus much less demanding than refining some PES parameters.

-This model allows computing high-J rovibrational energy levels in only few minutes from the diagonalization of matrices whose the size does not usually exceed 1000.

-If the DMS is available, an ab initio effective dipole moment is derived almost automatically for line intensity calculations.

-Finally, this model attempts changing the time-scale from months or years to some days or weeks for understanding and modelling the main spectral features of polyatomic molecules, even for transitions involving high energy rovibrational states.

Usually, the formal derivation of ab initio EMs in quantum mechanics is based on perturbation theory which has a long history (see for example the review by Watson 73 ). The oldest approach is known as the Rayleigh-Schrödinger perturbation method 74 which has been generalized to tackle the case of quasi-degenerate states 75,76 . Another approach was suggested by Van Vleck 77 and is known as canonical perturbation theory or contact transformation (CT) method [78][79][80][81][82][83][84][85][86][87] falling in the domain of the Lie algebra 88 . The basic idea is to a apply a series of unitary transformations T CT to the nuclear motion Hamiltonian in order to transform it to a new, block-diagonal representation. Very briefly, in perturbation theory the untransformed Hamiltonian is commonly expanded as a power series

H ≡ H (Γ 0 ) = H 0 + λH 1 + λ 2 H 2 + • • •
where H 0 is the zero order model, λ a small formal parameter and Γ 0 the totally symmetric irreducible representation of G. Following Van Vleck, the initial Hamiltonian H (Γ 0 ) is transformed to an effective one H(Γ 0 ) from successive unitary transformations as

H(Γ 0 ) = T -1 CT H (Γ 0 ) T CT , T CT = e -iλS (Γ 0 ) 1 e -iλ 2 S (Γ 0 ) 2 • • • , (3)
where the S (Γ 0 ) n are the generators of the CTs which insure a block-diagonal form. In order to accelerate calculation of high-order contributions, the modern algorithms employ the super-operator technique proposed by Primas 88 (see e.g. the MOL -CT computer code 89 by Tyuterev & Tashkun designed for rigid asymmetric top molecules of symmetry C s and C 2v ).

In this work, we present a numerical method which is not based on perturbation theory. Some tedious algebraic calculations involved in the derivation of effective Hamiltonian and transition dipole moment operators are clearly obviated.

B. Methodology for the derivation of "polyad" Hamiltonian and dipole moment operators

Before giving more details in Section III, our procedure for building an ab initio EM up to a polyad P N P can be summarized in three steps as follows.

(i) First, we start from the complete nuclear motion Hamiltonian H (Γ 0 ) and compute the full matrix

H (J,C) = γ ; J, C | H (Γ 0 ) | γ; J, C , (4)
in a basis set {| γ; J, C }. This matrix is diagonalized for each block (J, C) for some low-J values (see Section IV) and the corresponding variational eigenpairs are stored. As already mentioned in Introduction, a large number of theoretical developments and associated computer codes based on sophisticated algorithms [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27] are now able to perform accurate variational calculations, even for molecules with more than 10 atoms.

(ii) Then, we search for a unitary transformation T P ≡ T (J,C) P that brings H (J,C) into block diagonal form up to a maximum polyad P N P following a polyad scheme defined by (1) or (2). In other words, the new matrix H (J,C,P )

Polyad = T -1 P H (J,C) T P = [H (J,C) P 0 ⊕ • • • ⊕ H (J,C) P N P ] ⊕ H (J,C) R , (5)
has the same eigenvalues as H (J,C) but is now composed of N P + 2 blocks. The last block H R includes all the "remaining" polyads that are not of interest for the present study (see Section III B).

(iii) In a last step, we assume that there exists a phenomenological Hamiltonian H(s) with a set of unknown effective parameters {s} = {s P 0 , sP 1 , • • • , sP Np } associated with each polyad P k . We propose to determine these parameters such that the matrix representation of H(s) in a basis set {| γ; J, C P k ≡| m P k } restricted to a polyad P k matches the corresponding block in Eq. ( 5). So, a set of parameters sP k is determined by solving

P k m | H(s) | m P k = (H (J,C) P k ) m m =⇒ {s P k }. (6)
The present procedure amounts to (a) searching for the matrix representation of T CT in Eq. (3) without knowing explicitly the generators S n and (b) computing the block-diagonal representation of an effective operator H(Γ 0 ) before knowing its parameters. A non-empirical EM is thus deduced without performing a perturbative expansion of H and T CT . This clearly obviates the need to compute thousands or even millions of multiple commutators of the type [iS

(Γ 0 ) n , • • • [iS (Γ 0 ) n , H (Γ 0 ) j ] • • • ],
as both the order of expansion and the number of atoms N increase. Moreover, our procedure allows a full account of the symmetry properties through the use of irreducible tensor operators adapted to Abelian and non-Abelian point groups. This generalizes the work by Sadovskii & Zhilinskii 90 on the dyad of methane to arbitrary point groups and polyads. Most importantly, we will also see that the construction of effective operators other than Hamiltonians (e.g. dipole moment, polarizability, etc.) is carried out for the very first time in a systematic and direct manner from the transformation T P defined in Eq. ( 5). To our knowledge, no effective dipole moment operator obtained by CT in a systematic manner has been published so far, except for some simple asymmetric rotors [START_REF] Camy-Peyret | Vibration-rotation dipole moment operator for asymmetric rotors[END_REF][START_REF] Lamouroux | [END_REF] .

III. COMPUTATIONAL DETAILS

Let us now focus with more details on the three steps presented in Section II. We just assume that the Hamiltonian H (Γ 0 ) is written as a sum-of-product of irreducible tensor operators to deal with arbitrary D n , C nv , D nh , D nd , T d or O h point groups.

A. Variational calculation (step 1)

For calculating rovibrational energy levels and eigenstates, it is common to first solve the J = 0 stationary Schrödinger equation H

(Γ 0 ) v Ψ (Cv) v = E v Ψ (Cv) v
for each symmetry block

C v where H (Γ 0 ) v = T v + V refers to the vibrational Hamiltonian. To this end, we consider a basis {Φ (Cv) v,j } of a subspace F r spanned by M Cv primitive functions selected through the pruning condition F λ (p) = Nm i=1 λ i v i ≤ p with v i = 0, • • • , p
and where λ i are weight coefficients. p and λ i are chosen to properly converge the vibrational levels. So far the largest molecule our computer code was able to treat using this direct-product basis set is SF 6 48 . Undoubtedly, our procedure would strongly benefit from the use of reduced-dimension

Hamiltonians 28 to construct contracted bending, stretching or torsional functions, as done in many other theoretical developments and computer codes. An updated version of our TENSOR computer code is currently in progress to tackle big molecules (N > 7) by using contracted functions in conjunction with irreducible tensor operators for a full account of symmetry.

Once the J = 0 problem solved, the standard procedure consists in retaining

N Cv vi- brational eigenvectors Ψ (Cv) v,i = M Cv j=1 U (Cv) j,i φ (Cv) v,j , (i = 1, • • • , N Cv ), for making J > 0 calculations 35 . H v is diagonal in the product contracted basis (Ψ (Cv) v,i ⊗ Ψ (Cr) r
) ≡| ir where Ψ (Cr) r are symmetry-adapted rotational functions. Computation of the rovibrational matrix elements i r | H vr | ir is much more costly because it requires making sums over j and j with the evaluation of the matrix elements j r | H vr | jr which cannot be stored in memory when M v is large. Several dedicated algorithms and optimized computer codes (e.g. GENIUSH 15 or TROVE 12 ) have been developed to compute variational solutions for

J > 0.
In this work, we follow the strategy of Refs. 26,93 which consists in considering a smaller

vibrational basis set F λ (p ) ∈ F p of dimension M Cv << M Cv such that the full vibrational eigenvectors (i = 1, • • • , M Cv ) can be decomposed as Ψ (Cv) v,i = j∈F p U (Cv) j,i φ (Cv) v,j + k / ∈F p U (Cv) k,i φ (Cv) v,k . (7)
A set of M v eigenvectors is then selected from analysis of the Gram matrix [U (Cv) ] t U (Cv) before introducing approximate or reduced (p → p ) eigenvectors

Ψ (Cv) app,v,i = j∈F p U (Cv) app,j,i φ (Cv) v,j (8)
where U (C) app ∈ R M v ×M v contains a set of vectors orthonormalized using the Gram-Schmidt algorithm. In order to further reduce the dimension of the problem, N Cv eigenvectors can be retained using the so-called "VSS parameter" (see e.g. Fig. 5 of Ref. 94 as an illustration).

Recently, this procedure allowed computing rovibrational energy levels up to J = 120 for a seven atomic molecule 48 . However, Ψ app,v are not eigenvectors of H

(Γ 0 ) v
, strictly speaking but we can show that U

(Cv) app,j,i ≈ U (Cv) j,i (j = 1, • • • , M v )
in the energy range of many molecules where observation are available, making the second term in the right-hand side of Eq. ( 7) very small. Obviously, if p = p then Ψ app,v = Ψ v and the "usual" procedure is applied.

Finally, it can be shown that the rovibrational Hamiltonian matrix (4) reads

H (J,C) ≈ g ⊕ i=1 q i =0 E (C i ) v ⊗ I q i + T -1 v h rv T v , (9)
where h rv is the matrix representation of

H (Γ 0 ) -H (Γ 0 ) v
built in the rovibrational primitive basis of dimension M , g the number of classes of G and T v a unitary transformation given by

T v = g ⊕ i=1 q i =0 U (C i ) app ⊗ I q i . (10)
In Eqs. ( 9) and ( 10), q i is a multiplicity index, that is the number of allowed vibrational basis functions Φ

(C i ) v
which appear in the rovibrational functions of symmetry C. E

(C i ) v
is a diagonal matrix which contains the variationally-computed energy levels E v . To improve the "quality" of the eigensolutions of H (J,C) , namely

{E vr , Ψ (J,C) vr = [U (J,C) ] -1 Φ vr }, (11)
the E v values in Eq. ( 9) can be directly replaced by the observed levels, when available.

As an illustration, let us focus on the calculation of the energy levels of PH 3 for the block

(J = 3, E) using either (a) Ψ (Cv) v or (b) Ψ (Cv)
app,v . More details about the construction of the model can be found in Section IV. The vibrational problem is solved using the basis F (14) which leads to symmetry blocks of dimensions M A 1 = 6945, M A 2 = 5985 and M E = 12915.

(a) Without "reduction", the number of rovibrational primitive functions using the pruned basis F (p = 14) is of 90435 for the block (J = 3, E). In the variational calculation, we have retained only 5% of eigenvectors, that is N v = 0.05M v , leading to a final symmetry block of dimension 4522. The corresponding eigenvalues, which are all converged within 10 -4 cm -1 up to 4800 cm -1 with respect to the calculation with N v = M v , are taken as benchmark. The variational calculation took 80 min on a computer with 28 processors.

(b) The vibrational eigenvectors are now "reduced" using a pruned basis F (p ), with p =8, 10 and 11, of dimensions (M A 1 , M A 2 , M E )=(590,415,999), (1506,1170,2666) and (2290,1842,4122) and

N v = 0.65M v , N v = 0.25M v and N v = 0.15M v "approximate" eigen-
vectors are retained, respectively. The final dimensions of the block (J = 3, E) are of 4553, 4675 and 4330, respectively, and are similar to that in (a). Fig. 1 shows the rotational errors between with respect to calculation (a). The variational calculation using the reduction 14→8 with the basis F (p = 8) took 20 seconds and leads to errors below 0.003 cm -1 up to 2300 cm -1 and of ∼0.06 cm -1 up to 4800 cm -1 using very few basis functions. The calculation using F (p = 10) took 2 min and introduced errors of 0.003 cm -1 up to 3500 cm -1 and of 0.01 cm -1 up to 4800 cm -1 . Finally, the calculation using F (p = 11) took 5 min with errors less than 0.0008 cm -1 up to 3500 cm -1 and less than 0.003 cm -1 up to 4800 cm -1 . We can conclude that the reduction procedure proposed in Eqs. (7) and (8) gives quite good results and requires computation of much less matrix elements than in the "full" problem. This procedure is particularly suited when performing variational calculations with very high J values.

B. Block-diagonalization of the Hamiltonian matrix (step 2)

Now, we search for a way to transform the H (J,C) matrix ( 9) into a new matrix H (J,C,P ) block given by Eq. ( 5) following a specific polyad scheme defined by (1) or (2). There exist many published papers dedicated to this task (see e.g. Refs. [95][96][97][98][99][100] ) and infinitely many transformations bringing a given matrix into block diagonal form. The approach proposed by Cederbaum et al. 101,102 to block-diagonalize Hermitian matrices turns out very relevant for this work. They required that the transformation T P "changes the initial matrix as little as possible" which amounts to imposing the condition || T P -I ||=min where || . || is the Euclidean norm. This latter condition is sufficient to uniquely determine T P which is finally given by 101

T P = U (J,C) (U (J,C) BD ) t U (J,C) BD (U (J,C) BD ) t -1/2 , ( 12
)
where U (J,C) is the eigenvector matrix of Eq. ( 11) and U (J,C) BD its non-singular, block-diagonal part related to the choice (1) or (2). The success of the method is mainly governed by the good extraction of U

(J,C) BD from U (J,C
) which has to be made with care for polyads that are not well isolated. In order to build (5) up to the polyad P N P , we proceed in two steps.

1. Extraction of a single block containing all vibrational states up to P N P In a first step, we want to block-diagonalize the Hamiltonian matrix in two blocks: a block including all vibrational states and coupling terms from the polyad P 0 to the polyad P N P and another one including all the "remaining" states beyond P N P . In order to properly define the block-diagonal part U (J,C) BD leading to this structure, we have to organize U (J,C) in a convenient manner. To this end, the M primitive rovibrational basis functions | γ; J, C P which were used to build (9) are sorted in increasing values of P . The U (J,C) matrix is thus rearranged via a permutation P row of its rows as

P row U (J,C) = U (J,C) r
. If M denotes the number of basis functions up to the polyad P N P , then we extract from U r the first M rows and form the new matrix U r ∈ R M ×M . The M relevant vectors of U r associated with

P 0 • • • P N P are chosen such that the trace M i=1 [U t r U r ]
ii of the Gram matrix is maximal. This property can be obtained by applying a permutation P col of the columns. A new sorted matrix is written as

U (J,C) r P col = U (J,C) s . (13)
Finally, we can decompose the sorted matrix (13) as the sum of a block-diagonal and antidiagonal block part as

U (J,C) s =   U P 0-Np s 0 0 U R s   +   0 U 12 s U 21 s 0   = U (J,C) BD + U (J,C) ABD , (14)

Complete nuclear-motion Hamiltonian (input)

Eigenpair, Eq. (11) Effective Hamiltonian, Eq. ( 18) parameters Hamiltonian matrix, Eq. ( 9) Eq. ( 13 BD in the transformation (12). One can show that the transformed eigenvector matrix reads Ũ(J,C)

s = T -1 P 0-Np U (J,C) s = U l V t r =   ŨP 0-Np s 0 0 ŨR s   , (15)
where U l and V r are the left and right orthogonal matrices involved in the singular de-

composition U BD = U l DV t r .
Here, D is a diagonal matrix with positive singular values as entries. ŨR s ∈ R M -M ×M -M is the block which describes all the "remaining" polyads that are not of great interest for our purposes. The corresponding block-diagonal Hamiltonian is obtained from H (J,C,P )

Polyad = Ũ(J,C) s E c vr [ Ũ(J,C) s ] -1 = H (J,C,P ) P 0-Np ⊕ H (J,C) R , (16)
where E c vr is a diagonal matrix composed by the energy levels E vr (11) after applying the permutation vector (13).

Extraction of each individual block P k

At this stage, H (J,C,P ) P 0-Np can be further decomposed as H (J,C,P 0 ) ⊕ H (J,C,P 1 ) ⊕ • • • . To this end, we first extract H (J,C,P 0 ) from a transformation T P 0 that block-diagonalizes ŨP 0-Np

s as ŨP 0 s ⊕ ŨP 1-Np s
, then ŨP 1-Np s will be block-diagonalized from T P 1 , and so on. The process ( 13)-( 16) is repeated until getting the desired equation ( 5) which is finally obtained by applying an unitary transformation formed by the N p successive matrix products

T (J,C) P = T P 0-Np T P 0 T P 1 • • • T P N P -1 . (17)
T (J,C) P will completely block-diagonalize the Hamiltonian matrix (9) up to the polyad

P N P . The last block, namely H (J,C) R
, is not considered in this work. By analogy with perturbative CTs, Eq. ( 17) could be thus seen as the matrix representation of T CT = exp(-iλS 1 )exp(-iλ 2 S 2 ) • • • , without needing to perform a perturbative expansion in λ.

The other advantage of the proposed approach is the possibility of substituting E vr by E obs vr , when available, directly in Eq. ( 16). Consequently, the eigenvalues of the final block-diagonal matrix (5) can match observation without any fit.

At this stage, we have obtained the block-diagonal representation (5) of the full nuclear Hamiltonian matrix (4) using the transformation (17). Usually, we start from an Hamiltonian operator with known parameters and we compute the matrix elements to deduce its energy spectrum after diagonalization. Conversely, we start in this work from a known matrix H (J,C) P k for a polyad P k and we search for the corresponding effective Hamiltonian.

Mathematically speaking, we have to deal with an inverse problem for determining a set of unknown parameters {s} = {s P 0 , sP 1 , • • • , sP Np } of an effective operator H(s) by solving Eq. ( 6) for each P k .

In the effective Hamiltonian theory, it is common to use creation-annihilation operators because the terms of the type a +m i i a n i i must satisfy the resonance condition 84 c i (m i -n i ) = 0 where the c i 's are defined in (1). Within this condition, it is possible to write a formal Hamiltonian for each polyad by forming the Γ 0 -invariant polynomials as

H(Γ 0 ) = N P i=0 h i j=1 sj,P i V Ωv(Γ) {α,P i } ⊗ R Ωr(Kr,αrΓ) (Γ 0 ) j = N P i=1 h i j=1 sj,P i Õij , (18)
where V and R are vibrational and rotational operators of degree Ω v and Ω r in a i , a + i and in the total angular momentum components J α , respectively. K r is the rank of the tensor in SO(3) and the parity in the conjugate momenta with = (-1) Ωr because of the time reversal invariance. α is a set of vibrational labels associated with each vibrational mode.

More details about the construction of these tensor operators can be found elsewhere [START_REF] Champion | Spherical Top Spectra[END_REF]103 .

Using the so-called vibrational extrapolation scheme 7 , the effective Hamiltonian (18) for a given polyad contains the contributions of all lower polyads. In that case, the parameters specific to the next polyad should be, in principle, only small corrections. The major limitation of the vibrational extrapolation is that spectra analysis has to be carried out polyad by polyad. For example, if we want to study the polyad P 3 of a molecule, the polyads P 0 , P 1 and P 2 need generally to be studied before. In our approach, all the s parameters are determined simultaneously.

Finally, from Eqs. ( 6) and ( 18) the h k parameters sj,P k specific to a polyad P k can be determined iteratively by solving the overdetermined system of equations using the dgelsy routine of LAPACK

h k j=1 sj,P k Õik m m = H (J,C,P k ) m m - k-1 i=1 h i j=1 sj,P i Õij m m , (19)
where Õij m m are the matrix elements of Õij computed using the Wigner-Eckart theorem adapted to the point group. Only some low-J values will be sufficient to determine the full set of parameters. As a rule, we have to compute the variational eigenpairs at least up to J = Ω r /2 (resp. (Ω r + 1)/2) if Ω r even (resp. odd) such that number of equations≥ number of unknowns. In this work, the system ( 19) is first solved for the purely vibrational part to deduce the corresponding parameters. The eigenvalues of the Hamiltonian (18) will reproduce exactly the variationally-computed or observed levels if Ω v has been conveniently chosen. For example, if a polyad contains the band 4ν i , then we should consider Ω max v = 8.

As to the rovibrational part, the J > 0 levels of the Hamiltonian (18) whose parameters have been obtained by solving (19) will not match exactly those in Eq. ( 11) because convergence of the variational calculation is rarely fully achieved. This is in principle not an issue because the major part of the resonance parameters will be fixed to their consistent ab initio values while only few ones will be slightly refined to observation. This strategy was already applied with success in the framework of CT (see e.g. Ref. 104 ).

Effective dipole moment operator

Last but not least, our approach is not restricted to the construction of effective Hamiltonians but can be also applied to the derivation of the transformed dipole moment for line intensity calculations. Let M Θ (Θ = X, Y , Z) be the matrices of the laboratory-fixed frame dipole moment components (C (Γ ) ⊗ M (Γ) ) ( Γ) computed in the same primitive basis as the Hamiltonian. Here, M

Such a construction turns out very tedious in the framework of CTs because of the rotational algebra which is larger. Indeed, the direction cosines λ Θα do not commute with J β , making calculation significantly more complicated. This mainly explains why only few studies were published on triatomic molecules these past few decades [START_REF] Camy-Peyret | Vibration-rotation dipole moment operator for asymmetric rotors[END_REF][START_REF] Lamouroux | [END_REF] , and not in a systematic manner.

In this work, we adopt the same strategy as for the Hamiltonian for building an effective dipole moment, namely we start by forming all the Γ-covariant tensors allowed by the symmetry

M (Γ) α = N P i,k=0 h ik j=1 μj,P i P k V Ωv(Γv) {α,P i P k } ⊗ R Ωr(Kr,αrΓr) (Γ) j , (21)
up to a given order Ω = Ω v + Ω r where Γ is not necessary the totally symmetric irrep. Then, we determine the effective parameters μ such that the matrix elements of (

C (Γ ) ⊗ M (Γ) ) ( Γ)
match MΘ in Eq. ( 20). Again, only few transitions between P i and P k will be required to determine the whole set of parameters. Typically, all the dipole moment parameters can be obtained from transitions with J ≤ 2. Finally, the methodology presented in Section III is schematically depicted in Fig. 2.

IV. CONVERGENCE STUDY AND SPECTRA CALCULATION: VALIDATION AND DISCUSSION

In order to validate theory, we focus on five candidates whose analysis of high-resolution spectra in the infrared is not complete: H 2 CO, PH schemes, convergence studies, spectra comparisons and vibrational extrapolation are now given.

• H 2 CO: Starting from the PES of Ref. [START_REF] Nikitin | [END_REF] , whose geometry has been slightly refined for this work using the strategy of Ref. 26 , we have considered the reduced 18→9 normal-mode model and the reduced 18→10 basis for a variational calculation in the C 2v point group (see e.g. Ref. 26 for the various definitions). A brief inspection of the harmonic frequencies would suggest a polyad scheme (1) with c = {2.7, 1.8, 1.5, 1.2, 2.7, 1.2}. Up to 2000 cm -1 , such a scheme works quite well but many "effective" polyads overlap beyond with for instance the low-lying energies of the polyad P 9 far below those of the polyad P 7 . In this work, we choose c = {2, 1, 1, 1, 2, 1} which obviously leads to bigger polyads but in turn to a more consistent treatment, at least up to 4000 cm -1 . We have built two effective Hamiltonians (18) up to the polyad P 3 using the variational eigenpairs (11) up to J = 5 to solve (19). The first one was expanded up to Ω v = 6 and Ω r = 6 and the second one up to Ω v = 8 and Ω r = 6. They are composed respectively of 722 and 2740 irreducible tensor operators. Fig. 3a displays the convergence error of the energy levels up to J = 8 for these two models with respect to the variational calculation, taken as the benchmark. We can see that the errors between effective and variational calculation are below 10 -3 cm -1 up to 3000 cm -1 and ∼10 -3 up to 5000 cm -1 when Ω max v = 8. The results are thus quite similar, except that energy level calculation took less than 1 min up to J = 40 using the effective model against few hours from the variational approach. A comparison between the reduced energy levels obtained from our global model and from those used to build HITRAN 3 is given in Fig. 3 (panels b & c). Unlike traditional empirical EMs, we can see that our model can predict all bright and dark states, here up to the polyad P 3 , and could predict other ones far beyond if necessary.

• PH 3 : Previous calculations were published for this molecule 107 . Here, we have considered the PES and DMS of Refs. 108,109 where the quadratic force constants have been slightly refined for this work and the geometry fixed to r e = 1.4119002 Å, α e = 93.39082 degrees using the optimization procedure described in Ref. 26 . A variational calculation was carried out in the C 3v point group using the normal mode model 18→10 and the basis 17→10 up to J = 5. Unlike H 2 CO, the PH 3 molecule exhibits a clear polyad structure with a resonance scheme 2:1:2:1. Using this polyad model, the effective Hamiltonian ( 18) was expanded at order 10 up to P 5 with Ω v ≤ 10 and Ω r ≤ 8 and the resulting 10415 rovibrational parameters have been determined by solving (19). Tab. I gives a comparison between some effective TABLE I. Comparison between the effective parameters a s (in cm -1 ) obtained from Eq. ( 19) and from a fit for the ground state and dyad of PH 3 . All the coupling parameters between ν 2 and ν 4

were fixed to their "ab initio" values.

Operator s, Eq. ( 19) s, fit a

Ground state (P = 0)

R 2(2,0A1) -1.08859121×10 -1 -1.08856(2)×10 -1 R 2(0,0A1) -1.85097033 -1.850969(1) R 4(4,0A1) -4.21792347×10 -6 -4.224(3)×10 -6 R 4(4,3A1) 6.35064638×10 -6 6.393(9)×10 -6 R 4(2,0A1) -4.54280135×10 -6 -4.530(5)×10 -6 R 4(0,0A1) -1.91476610×10 -5 -1.9143(3)×10 -5 RMS=0.0018 RMS=0.00013 b Dyad(P = 1) R 0(0,0A1) (a + 2 ⊗ a 2 ) (A 1 ) 992.134856 992.1354(1) R 2(2,0A1) (a + 2 ⊗ a 2 ) (A 1 ) 6.11040167×10 -3 5.877(5)×10 -3 R 2(0,0A1) (a + 2 ⊗ a 2 ) (A 1 ) -2.85659105×10 -3 -2.863(6)×10 -3 R 0(0,0A1) (a + 4 ⊗ a 4 ) (A 1 )
1581.52444 1581.5238(1)

R 2(2,0A1) (a + 4 ⊗ a 4 ) (A 1 ) -1.05612963×10 -2 -1.0614(5)×10 -2 R 2(0,0A1) (a + 4 ⊗ a 4 ) (A 1 ) -1.12134942×10 -3 -1.119(6)×10 -3 R 2(2,1E) (a + 4 ⊗ a 4 ) (E) 2.18084810×10 -2 2.176(3)×10 -2 R 2(2,2E) (a + 4 ⊗ a 4 ) (E) -5.13139342×10 -5 -3.52(7)×10 -5
RMS=0.019 RMS=0.00057 b a 6 parameters fitted among 24 at order 10 for the ground state. 8 parameters fitted among 150 at order 10 for the dyad. The parameters that are not shown were fixed to their "ab initio" values.

b Fit of 171 levels up to J = 15 for the ground state and of 288 levels up to J = 11 for the dyad.

The observed energy levels were taken from HITRAN.

ground state and dyad parameters s derived from Eq. ( 19) and obtained from a fit. We can see that the initial parameters s only need to be slightly refined to reproduce high-resolution data.

Concerning the line intensity calculation, the effective dipole moment ( 21) was built at order 6 with Ω r ≤ 2 to account for all transitions P i ← P j (i, j=0,• • • ,5). The corresponding 8519 parameters were determined from the variationally-computed transitions up to J = 2. We thus obtain the first global effective model for PH 3 up to P 5 whereas the current spectroscopic models derived from high-resolution analysis these past three decades are available only up to P 3 [110][111][112] . Due to its completeness, our model will be used for the analysis of P 4 and P 5 as well as for the construction of hot line lists for astrophysical applications.

In order to validate this model, Fig. 4(a) shows the convergence of the rovibrational levels up to J = 8 using the Hamiltonian (18) at order 6, 8 and 10, with respect to the variational energies in the P 4 region. The two other panels in Fig. 4 show a comparison of both line positions and line intensities up to J = 15 in the P 4 region between our effective model and the variational calculation. We can see a very good agreement between these two approaches, even for small transitions (say < 10 -23 cm/molecule) with errors on the line intensities below 0.1%. Fig. 5 gives an overview of the spectrum up to P 5 computed from our global model ( 18) and (21). In this figure, the consistency of our model is clearly seen from a direct comparison of the theoretical absorption spectrum for the polyad P 4 with experimental PNNL records 105 .

• 12 CH 4 : We started from the ab initio model described in Ref. 114 using the PES and DMS of Refs. 115,116 . As for PH 3 , the methane polyad scheme is quite clear and is also given by the resonance scheme 2:1:2:1. In this work, a very first global effective model beyond tetradecad (P 4 ) 117 was built using an effective Hamiltonian expanded at order 8, up to Ω r = 4 (7912 parameters) and an effective dipole moment at order 5 with Ω r ≤ 2 (522 parameters), only for cold band transitions. A higher-order expansion in Ω r will be easily carried out for practical studies. The first global predictions of 12 CH 4 (T d ) spectra in the icosad (P 5 ) and triacontad (P 6 ) regions using an effective model are given in Fig. 6.

Some observed energy levels for low-J values in the icosad (P 5 ) determined from the partial assignment 29 of the Grenoble WKLMC experimental line list 113 have been directly replaced in Eq. ( 16). In the next polyad (P 6 ), the big deviations on line positions are partly explained by the "poor" accuracy (of ∼0.1-1 cm -1 ) of the PES 115 in this region. Though most of the line will deserve to be refined in position, the agreement with high-resolution experimental spectra remains quite good, in particular if we remember that the complete modelling of room-temperature methane spectra up to the tetradecad (P 4 ) took more than forty years.

Undoubtedly, our model will allow to complete the current assignments in the regions above 12 CH 4 up to the polyad P 6 using our model and that used to build HITRAN these past forty decades (up to the lower edge of P 5 ). Calculated spectra compared to the Grenoble WKLMC line list 113 at 80 K and to HITRAN at 296 K in the P 5 and P 6 region are also given.

P 5 118,119 in a shorter time.

• 12 C 2 H : Several studies were focused on ethylene 36,120-122 and an accurate line list was built 44 . Here, we start from the variational calculation performed in the D 2h point group which is based on the PES and DMS of Refs. 123,124 and described in detail in Ref. 125 . This molecule is the typical case where the choice (1) is unclear beyond 2000 cm -1 , so the energy criterion (2) was used to define the successive polyads up to 5000 cm -1 . For example, the polyad P 1 falls in the region [800, 1050] cm or remove vibrational states quite easily. The panel (a) in Fig. 7 gives an overview of the ethylene spectrum up to 5000 cm -1 between our calculated EM and those used to generate HITRAN2020. The panels (b-e) show the very good agreement between our calculated spectra in the P 2 , P 5 , P 8 and P 13 regions and the PNNL databases 105 and Fig. 8 illustrates the so-called vibrational extrapolation. In this figure, the "global EM P 13 " corresponds to the full effective model up to P 13 whereas the "global EM P 8 " contains spectroscopic parameters only up to P 8 but was used to generate the P 13 spectrum by extrapolation. As expected, the agreement of the P 13 model with experiment is better though the P 8 model gives reasonable results despite fewer parameters (16081 for P 13 vs. 4157 for P 8 ).

Once again, a first "global" effective model is presented for this molecule, both for line positions and line intensities. For example, our model takes into account 52 interacting vibrational bands in the 3 µm region, which is to be compared to the 4 bands included in the current HITRAN database. In the present EM, the resonances between all bright and dark states induce intensity transfers that cannot be properly described using the available empirical effective Hamiltonians. We can note that the variational calculation presented in Ref. 125 up to 6400 cm -1 and J = 71 took five days using 96 processors while the computed spectra up to J = 40 in Fig. 7 were built on a laptop in 1 hour only.

• 32 SF 6 : We conclude this work by an illustration on the 32 SF 6 (O h ) molecule for which completeness is crucial. Here, the PES and DMS of Ref. 126 were employed and the Hamiltonian model as well as the choice of the basis functions were described in Ref. 48 . According to the six harmonic frequencies, we have fixed the polyad vector (1) as c = {2.3, 1.8, 2.8, 1.5, 1}.

In order to take into account all the cold and hot band transitions falling in the strongest absorption region around 10 µm, a polyad model up to P 24 has been constructed. In order to consider 5-quanta vibrational bands, we have developed the effective Hamiltonian (18) at order 10 with Ω r ≤ 8. The resulting 9773 parameters have been determined from the variational eigenpairs up to J = 10. The dipole moment has been expanded at order 5 and contains all operators corresponding to the cold and hot bands having the strongest integrated intensities, as tabulated in Supplementary Material of Ref. 48 . A set of 4311 effective dipole moment parameters has been determined from variationally-computed transitions up to J = 2. Fig. 9a shows a comparison of both line positions and line intensities up to J = 80 in the P 7 region between our effective model and the variational calculation of Ref. 48 .

The most extensive empirical "global" model to date contains only 13 vibrational bands 127,128 and has been used to generate the SF 6 HITRAN database (see Fig. 9b,c).

In this work, hundreds of cold and hot bands have been generated in 5 hours up to J = 120.

The very good agreement with PNNL in Fig. 9c validates the construction of our effective Hamiltonian (18) and dipole moment (21). We can conclude that all the missing hot bands in HITRAN-like databases cannot provide good opacity calculations required for accurate atmospheric applications.

V. CONCLUSION AND PROSPECTS

To summarize, our TENSOR computer code initially designed for variational calcula- we have shown that the present numerical approach obviates the need to know explicitly the rotation-vibration algebra, with sometimes very involved calculations.

We have demonstrated that the proposed effective model can provide crucial information to spectroscopists (e.g. a good set of initial parameters) within a very short time and is a clear alternative to more traditional spectroscopic models if, for a given molecule, the ab initio surfaces are available. So far, only variational calculations were able to provide complete molecular line lists for different temperature and spectral ranges contrary to empirical EMs which are able to make accurate but quite "localized" predictions. Finally, the contribution of the proposed effective approach is threefold: (1) it is much less demanding than performing a full variational calculation, (2) it allows computing very rapidly and simultaneously all the cold (P i ← P 0 ) and hot (P i ← P k ) band transitions up to a given polyad and (3) the parameters can be refined much more easily. Undoubtedly, this model brings a new insight into high-resolution spectra analysis and will be of great help, not only in current or future infrared spectra analyses of polyatomic molecules but also in the modelling of hot atmospheres for which completeness is crucial.

For future studies, the case of nonrigid molecules exhibiting one or several inversion or internal rotation motions could be considered. Numerous papers dealing with the construction of empirical effective models for nonrigid molecules have been published these past four decades (see e.g. Refs. [129][130][131][132][133] and references therein). We could imagine that the nonrigid counterpart of ( 18) can be also derived in the same fashion using a convenient polyad scheme and parametrization 134 . We will show in a forthcoming study that the tensorial Hamiltonian ( 18) can be formulated in the Hougen-Bunker-Johns (HBJ) formalism 135 by making the substitution s → s(ρ) where ρ is a large amplitude curvilinear coordinate. The effective parameters for the rigid motions will be determined on a numerical grid ρ i . Coordinate transformations in the framework of HBJ as well as choice of the axis system where the angular momentum due to the large amplitude motion vanishes or is of constant magnitude have been recently discussed 136,137 . The first candidate will be the NH 3 molecule for which accurate PESs have been recently constructed 138 .

Another point concerns the construction of a "full" curvilinear EM. This could be done for example by using the curvilinear creation-annihilation operators A and A + introduced in Ref. 139 . Other formalisms like those implemented in GENIUSH 15 or TROVE 12 could possibly benefit from the present study while anharmonic basis functions (Morse, etc.) associated with generalized ladder operators 140,141 could replace the harmonic oscillator basis functions. so that any Hamiltonian H(X i , Y i , J α , J β , J γ ; t) can be transformed like (18) and hundreds of rovibrational parameters t and s can be linked together in few seconds by solving linear systems of equations. We write ti =

ZOOM PH 3 FIG. 1 .

 31 FIG. 1. Rotational errors for the block (J = 3, E) of PH 3 introduced when using the reduced vibrational eigenfunction (14 → p ) of Eq. (8) with p =8, 10 and 11. The variational calculation using the basis F (14) is taken as reference (see text).

FIG. 2 .

 2 FIG. 2. Schematic representation of the methodology presented in Section III for a systematic construction of global, ab initio effective Hamiltonian and dipole moment operators.

  x, y, z) are the molecular-fixed frame dipole moment components where the µ parameters are determined from the ab initio dipole moment surfaces and C (Γ ) is the tensor counterpart 7 of the direction cosines λ Θα . Using the unitary transformation (17) that block-diagonalizes the Hamiltonian, the transformed dipole moment matrices involved in line strength calculations for transitions of the type J , C , P ← J, C, P are given by MΘ = [T (J ,C ) P ] -1 M Θ T (J,C) P .

FIG. 3 .

 3 FIG. 3. (Panel a) Convergence of the energy levels of H 2 CO using the effective Hamiltonian (18) expanded at order 6 (Ho6) and order 8 (Ho8) with respect to the variational energies (V) up to J = 8. (Panels b & c) Comparison of reduced energy levels E -B 0 J(J + 1) for H 2 CO up to the polyad P 3 using empirical effective Hamiltonians fitted to observation and extracted from the HITRAN database 3 and using our global effective Hamiltonian (18). The different colours in (c) represent all the eigenvector mixings.

FIG. 4 .

 4 FIG. 4. (Panel a) Convergence of the P 4 energy levels of PH 3 using the effective Hamiltonian (18) expanded at order 6 (Ho6), order 8 (Ho8) and order 10 (Ho10) with respect to the variational energies (V) up to J = 8. (Panels b & c) Comparison of line positions and line intensities in the polyad P 4 and P 5 between the effective model and variational calculations.

FIG. 5 .

 5 FIG. 5. Overview spectrum of PH 3 up to the polyad P 5 and comparison between first-principles calculation (black circles) and global effective model (red dots). A comparison of the theoretical absorption spectrum for the polyad P = 4 with experimental medium resolution PNNL records 105 at T = 298 K is also given..

FIG. 7 .

 7 FIG. 7. (a) Overview spectrum of 12 C 2 H 4 up to the 5000 cm -1 and comparison between the empirical effective (black dots) and global effective (red dots) models. A comparison of the theoretical absorption spectrum for the P 8 (panel b) and P 13 (panel c) polyad with experimental medium resolution PNNL records 105 at T = 298 K is given. In panels (d) and (e), detailed portions of the ethylene spectra in the P 2 and P 5 regions, respectively, and comparison with the PNNL experimental database 105 .

FIG. 8 .

 8 FIG. 8. Comparison between the PNNL experimental database 105 and calculated spectra of ethylene using the full P 13 model and the P 8 model extrapolated up to P 13 (see text).

FIG. 9 .

 9 FIG. 9. (a) Comparison of line positions and line intensities in the ν 3 region (P 7 ) of 32 SF 6 between the global effective (this work) and variational calculations 48 . (b) Overview spectrum of 32 SF 6 in the ν 3 region using our model and that used to build HITRAN. (c) Comparison with the PNNL database 105 (black curve) is given at 298K. Note that only 5 hours were necessary in our approach to get a global EM and compute the red spectrum up to J = 120 against two decades to obtain the blue spectrum.

  tions is now able to manage EMs for spectra analyses. More generally, the steps 2 & 3 in Section III, namely the block-diagonalization procedure and the construction of effective Hamiltonians, could be also implemented in most of the variational computer codes. For users who do not feel comfortable with the theory of irreducible tensor operators, the Hamiltonian(18) and dipole moment (21) could be easily replaced either by the standard Watson formalism[START_REF] Aliev | Higher-order effects in the vibration-rotation spectra of semirigid molecules[END_REF] or by another one. Alternatively, the parameters involved in different formalisms could be linked together, as shown in Appendix A. Unlike Van Vleck perturbation theory,

ic

  ij si , si = i d ij ti . (A4)Note that a recent paper142 provides relations between some rotational parameters involved in the tensorial and Watson formalisms for a limited number of point groups whereas illustrative examples were already given elsewhere for the vibration-rotation Hamiltonian (see Eqs. (31)-(33) of Ref. 143 and Tab. 1 and Eqs. (28)-(29) of Ref. 144 ). The same holds for the derivation of dipole moment parameters (see e.g. Eqs. (54)-(55) of Ref. 145 ). Obviously, the transformation (A2) or (A3) can be provided upon request. Sometimes, calculations are not performed in the highest symmetry point group but in one of its subgroups. If for example methane is treated in C 2v instead of T d our procedure allows converting any C 2v operator to a better suited symmetry-adapted form. The components of the total angular momentum in the C 2v (x , y , z ) frame can be thus expressed as a linear combination of rotational T d tensor operators as . The vibrational part could be also transformed in this way and other point groups could be considered (see e.g. Ref. 145 where the link between C 3v and T d was clearly established).

  3 , 12 CH 4 , 12 C 2 H 4 and 32 SF 6 . We recall that the aim of this work is to present a new methodology and give a general overview of what our global ab initio EM can bring to improve the current models. Future spectra analyses will be carried out with the help of specialists in the domain (i) who will use our parameters as input data in dedicated programs for spectra analysis and (ii) who will publish

refined spectroscopic parameters in dedicated papers in due time. Note that global effective model is often used in this paragraph instead of ab initio effective model in order to highlight the difference with the empirical models which usually focus on small spectral regions. All details concerning the construction of the ab initio and effective models, the various polyad

  -1 and contains 4 vibrational states, P 8 belongs to [2785, 3390] cm -1 and has 52 vibrational states while P 13 ∈ [4100, 4850] cm -1 and has 251 vibrational states. Such a polyad scheme can be customized by the user who can add

ACKNOWLEDGMENTS

The author acknowledges support from the Romeo computer center of Reims-Champagne-Ardenne, from the French-Russian collaboration program IRP "SAMIA2", and from the French ANR TEMMEX project (Grant 21-CE30-0053-01). The author warmly acknowledges Dr. Andrei Nikitin for providing the ab initio surfaces.

Appendix A: Link between the tensorial formalism and other ones Irreducible tensor operators are great tools to deal with symmetry for both Abelian and non-Abelian point groups but in turn may appear as abstract mathematical objects for many users. We can show that any rovibrational tensor operator T

iσ can be expressed in terms of "standard" vibrational (X, Y ) = (q, p), (a + , a) or (S, -id/dS) and rotational (J x , J y , J z ) or (J + , J + , J z ) operators of degree Ω v and Ω r , respectively. It is common to write these operators in a normally-ordered form. Typically, if we know the structure constants c k ij of a Lie algebra (e.g. Heinsenberg-Weyl, so(n), su(n), su(m, n), etc.), defined by a set of generators G 1 , G 2 , etc. satisfying the commutation rules [G i , G j ] = c k ij G k , any linear combination of arbitrary powers of G i can be normally ordered as

Such ordering is greatly simplified with the aid of symbolic calculation where non commutative objects can be now easily manipulated. Eq. (A1) helped us to express the rotational tensors R Ωr(Kr,nrΓr) in terms of angular momentum components. For example, a rovibrational coupling term for the ν 4 band of a D 3h molecule reads

Eq. (A1) could be also used for the derivation of the structure constants of the rotational algebra supplemented by the direction cosines λ θ,α , involved for example in dipole moment contact transformations.

Finally, if M i denotes an operator of the type

, with (α, β, γ) = (x, y, z) or (+, -, z), the TENSOR computer code is able to expand each rovibrational tensor in elementary operators (see Eq. iσ and M i , respectively. In order to properly invert Eq. (A2), we must consider all the Γ-covariant tensors. For example, in the case of C 3v species the vector T is composed of T