Bruno Lepetit
email: bruno.lepetit@irsamc.ups-tlse.fr

Computation and analysis of bound vibrational spectra of the Neon tetramer using row orthonormal hyperspherical coordinates

Keywords:

Computation and analysis of bound vibrational spectra

of the Neon tetramer using row orthonormal hyperspherical coordinates.

I. INTRODUCTION

The hyperspherical coordinate method has been used extensively over the past decades to study few-body problems. It relies on the parametrization of a N-body system with a hyperradius ρ which gives to the global size of the system and a set of 3N-4 angles which describe the shape of the system. Since the method has the capability to describe the fragmentation channels of the system, it has been used extensively for collisional problems which involve continuum states, like reactive scattering (for example, ref. 1,2), ultracold physics [3], high energy nuclear collisions [4], as well as for bound state computations, for instance in atomic clusters [5] or nuclei [4]. It is surprising however that an overwhelming majority of these works concerns 3-body problems. This is an indication that it remains even today technically challenging to deal with 4 bodies without drastic approximations.

Yet, on rare occasions like the study of 4-body recombinations, the use of such methods was an essential ingredient to achieve significant progress in our understanding of ultra-cold fermionic [6] and bosonic [3,7] physics.

In the field of rare gas clusters Rg N also, almost all studies of (ro-)vibrational spectra concern trimers (for a review, ref. 5). As an exception, a pioneering study [8,9] of the vibrational spectra for Rg=He, Ne, Ar was performed using hyperspherical coordinates for N up to N=6, but within the frame of the drastic adiabatic approximation which restricts the validity of the results at most to a very limited number of low lying states. However, this study provided useful qualitative insight in the cluster dynamics. In fact, the only two accurate computations of (ro-)vibrational spectra for rare gas clusters beyond trimers which we are aware of were obtained using Jacobi coordinates and potential-optimized discrete variable representation [START_REF] Light | Advances in chemical physics[END_REF] functions. The first one is the vibrational spectrum computation of Ar 4 described in ref. [START_REF] Wang | [END_REF], the second one the full rovibrational spectrum computation of Ne 4 described in ref. 12. The recent availability of these latter Ne 4 spectra provides a good opportunity to test for the first time the capability of the hyperspherical coordinate method to obtain the vibrational spectrum of a tetratomic system beyond the drastic adiabatic approximation.

The hyperspherical harmonics [START_REF] Avery | Hyperspherical Harmonics: Applications in Quantum Theory[END_REF] are eigenfunctions of the fixed-hyperradius kinetic energy operator and form an orthogonal basis in which the complicated angular part of the kinetic energy operator has a very simple diagonal matrix representation. They can be used in 2 ways in the expansion of the wavefunction of the system [START_REF] Suzuki | [END_REF]. One way popular in nuclear sciences [4] is to combine them with hyperradial functions in product functions on which the total wavefunction of the system is expanded directly. Another way popular in chemical physics is to use them to expand the eigenfunctions of the fixed ρ-Hamiltonian which are called channel functions [6] or surface functions or, following ref. 15, local-hypersphericalsurface-functions (LHSF). These LHSF form a compact basis set on which to expand the total wavefunction of the system. This set of surface functions provides a hyperspherical adiabatic representation of the system containing useful physical information on its dynamical properties. The hyperspherical adiabatic representation becomes a hyperspherical adiabatic approximation [8,9] when a single surface function is kept in the expansion of the full wavefunction. In the present paper, we will show how analysis of the surface functions and of their ρ-dependent energies provides insight into the vibrational dynamics of the Ne 4 cluster.

In particular, we will identify the pyramid, rhombus and trigonal geometries which play an important role. For completeness, we note that correlated Gaussians have been used as a successful alternative to hyperspherical harmonics in the expansion of LHSF [6,[START_REF] Suzuki | [END_REF].

There are 2 types of hyperspherical coordinates and harmonics. The most standard one is the Delves type [6,16,17] obtained from an extension of the vector spherical representation to a set of several Jacobi vectors associated to a given clustering scheme of the particles.

The corresponding hyperspherical harmonics are known analytically [6,[START_REF] Avery | Hyperspherical Harmonics: Applications in Quantum Theory[END_REF]. They may not be always the optimum choice of coordinates in practical applications. Their angles are defined with respect to a space fixed frame and do not take advantage of the isotropy of space, as Euler angles would. Furthermore, an identical particle symmetrization operation is difficult to implement in this coordinate system because it potentially changes all angular coordinates [18]. Consequently, giving up implementing permutations on the wavefunction may be preferred in some applications even if the price to pay is larger computations [19]. In order to overcome these difficulties, symmetrized coordinates, which we call, following ref. 20, row-orthonormal-hyperspherical coordinates (ROHC), were developped in ref. 20,21 following the pioneering work of ref. 22,23. In this coordinate system, the 9 dimensional parameter space can be considered as the superposition of 3 subspaces. One is generated by the 3 Euler angles, which define the orientation in space of a body frame tied to the principal axes of inertia of the molecular system. The second subspace is the 3-dimensional kinematic invariant subspace [24] S KI . It contains the hyperradius ρ and 2 angles which provide a parametrization of the 3 moments of inertia of the system. Finally, the kinematic rotation subspace S KR is defined by 3 angles which parametrize the shape of the system. This coordinate system has distinctive advantages. The Euler angles allow to take full advantage of the isotropy of space and when the total angular momentum of the system is 0, the wavefunction becomes independent of them. The effects on coordinates of identical particle permutations is confined to S KR , which makes their implementation in the wavefunction easier. On the other hand, the hyperspherical harmonics are more difficult to express in ROHC than in Delves' coordinates. In early works, only a limited number of harmonics could be obtained analytically [22][23][24][25][26]. A general recursion relation was implemented in a mathematical formal calculation program to enlarge the set of harmonics known analytically [15], at the price of practical difficulties. For instance, the recursion provides over-complete sets of linearly-dependent harmonics which have to be culled to reduce them to linearlyindependent sets using a formal but computationally intensive algorithm. An alternative generation approach entirely numerical was therefore proposed in [27], it relies on a selection of simple trigonometric basis functions in which to expand the hyperspherical harmonics according to their symmetry properties. However, the resulting harmonics have never been used in physical applications. Therefore, the present paper represents to our knowledge the first successful use of 4-body ROHC harmonics to a molecular problem. It presents the results of the computation of the vibrational spectrum of the Neon tetramer using these symmetrized hyperspherical coordinates and harmonics. Section II presents the formalism, first the coordinates (subsection II A), then the harmonics used to expand the LHSF (subsection II B) and the full wavefunction expansion on the LHSF (subsection II C). Subsection II D provides technical details on our implementation of the Ne 4 problem. In section III, we first present the obtained vibrational spectra and we discuss their accuracy by comparison with the reference data of ref. [12]. Then, we describe the physical insight provided by the hyperspherical method : we present the LSHF energy curves as a function of ρ (subsection III B), the LHSF themselves (subsection III C) and the bound states (subsection III D).

II. METHOD

A. Coordinates and kinetic energy

The ROHC formalism used in the present paper have been described previously [15,20,27]. We summarize here important information and refer the reader to these papers for further details. Let ρ sf be the Jacobi matrix, whose columns are the three cartesian coordinates in an arbitrary space-fixed frame of the three mass-scaled Jacobi vectors (r (1) , r (2) , r (3) ) associated to a given clustering scheme of the 4 particles which we choose as follows : attributing arbitrarily numbers 1-4 to the 4 identical neon atoms, r (1) connects atoms 1 and 2, r (2) atom 3 to the center of mass of atoms 1 and 2, r (3) atom 4 to the center of mass of the 3 others. The ROHC χ, ρ, a, δ associated to this clustering scheme are defined implicitely by :

ρ sf = (-1) χ R(a)ρN(θ, φ) R(δ) (1)
χ is the chirality variable, a refers collectively to the three Euler angles (a, b, c) of the principal axes of inertia frame and ρ is the hyperradius which, together with the angles θ, φ, parametrizes the three principal moments of inertia of the system and defines the kinematic invariant space S KI . Three additional angles (δ 1 , δ 2 , δ 3 ) collectively denoted by δ parametrize the kinematic rotation space S KR . The matrix R(a) is the transpose of the rotation matrix defined by [28] :

R(a) = M 1 (c)M 2 (b)M 1 (a) (2)
where :

M 1 (ω) =      cos ω sin ω 0 -sin ω cos ω 0 0 0 1      (3)
and :

M 2 (ω) =      cos ω 0 -sin ω 0 1 0 sin ω 0 cos ω      (4)
These definitions are also valid for R(δ 1 , δ 2 , δ 3 ). The matrix N(θ, φ) is a diagonal matrix, with diagonal elements :

N 11 (θ, φ) = sin θ cos φ, N 22 (θ, φ) = sin θ sin φ, N 33 (θ, φ) = cos θ (5)
Equation 1 is called the row-orthonormal form of the Jacobi matrix ρ sf because matrix R(δ)

is row-orthonormal. In the present 4-particle case, the matrix R(δ) is square and thus also column-orthonormal.

It was previously shown [20] that a one to one correspondence can be achieved between physical configurations and coordinate space by restricting the hyperspherical angles to :

0 ≤ a, c < 2π, 0 ≤ b ≤ π (6) 0 ≤ δ 1 , δ 3 < π, 0 ≤ δ 2 ≤ π (7) 0 ≤ φ ≤ π/4, 0 ≤ θ ≤ θ M = arcsin[1/(1 + cos 2 φ) 1/2 ] ≤ arcsin(2/3) 1/2 54.7 o (8)
and χ = 0 or 1.

In terms of these ROHC, the kinetic energy operator associated to the motion of the 4 nuclei is given by [27] :

T = - 2 2µ 1 ρ 8 ∂ ∂ρ ρ 8 ∂ ∂ρ + Λ2 (a, θ, φ, δ) 2µρ 2 (9)
where µ = m/4

B. Hyperspherical harmonics

Since the 6 operators Λ2 , Ĵ2 , Ĵz (the component of Ĵ on the space-fixed z axis), L2 , L3

(the 3rd body-fixed type component of L) and Î (inversion of the Jacobi vectors r (1) , r (2) , r (3) to their opposites, which, according to eq. 1, means that the chirality coordinate changes according to : χ ↔ χ ± 1, while the other 9 coordinates are unchanged) commute with each other, we can find functions

F nJ M J L M L D d
, the hyperspherical harmonics, which are simultaneous eigenfunctions of those 6 operators :

Λ2 F nJ M J L M L D d (χ, a, θ, φ, δ) = n(n + 7) 2 F nJ M J L M L D d (χ, a, θ, φ, δ) (10) Ĵ2 F nJ M J L M L D d (χ, a, θ, φ, δ) = J(J + 1) 2 F nJ M J L M L D d (χ, a, θ, φ, δ) (11)
Ĵz

F nJ M J L M L λ D d (χ, a, θ, φ, δ) = M J F nJ M J L M L λ D d (χ, a, θ, φ, δ) (12)
L2

F nJ M J L M L λ D d (χ, a, θ, φ, δ) = L(L + 1) 2 F nJ M J L M L D d (χ, a, θ, φ, δ) (13)
L3

F nJ M J L M L D d (χ, a, θ, φ, δ) = M L F nJ M J L M L D d (χ, a, θ, φ, δ) (14)
ÎF nJ

M J L M L D d (χ, a, θ, φ, δ) = (-1) Π F nJ M J L M L D d (χ, a, θ, φ, δ) (15)
The quantum numbers n, J, M J , L, M L , Π appearing in these expressions are all integers, satisfying the constraints :

n ≥ 0 0 ≤ J, L ≤ n (16) -J ≤ M J ≤ J -L ≤ M L ≤ L (17)
Π = 0 or 1 (18)
The inversion quantum number Π is implicitely contained in the principal quantum number n through the relation : (-1) Π = (-1) n . The sub-and super-scripts d and D indicate that the F functions can be degenerate, i.e. for a given set of quantum numbers n, J, M J , L, M L , there are D linearly independent F functions. The subscript d indicates which of those degenerate functions is being considered.

These functions can be expanded on symmetrized Wigner functions defined by :

D J,η M J Ω J (a) = 1 2(1 + δ Ω J 0 ) D J M J Ω J (a) + ηD J M J -Ω J (a) (19)
D L,η Ω L M L (δ) = 1 2(1 + δ Ω L 0 ) D L Ω L M L (δ) + ηD L -Ω L M L (δ) (20)
where D J M J Ω J (a) and D L Ω L M L (δ) are the usual Wigner rotation matrices (ref. [START_REF] Davydov | Quantum Mechanics[END_REF], p. 154-156) and η = ±1. This expansion can be written as :

F nJ M J L M L D d (χ, a, θ, φ, δ) = (-1) Πχ J Ω J =0 L Ω L =0 η=±1 D J,η J η M J Ω J (a)D L,η L η Ω L M L (δ)G nJ Ω J L Ω L D d η (θ, φ) (21)
with η J = (-1) J and η L = (-1) L . The G(θ, φ) functions, which appear in the expansion of the F functions, are two-angle principal-axes-of-inertia hyperspherical harmonics. In the present implementation, they are obtained by expansion on simple symmetry adapted trigonometric primitive basis functions and by diagonalization of the Λ2 operator matrix representation in this basis, as described in ref. [27].

These harmonics can be further symmetrized with respect to permutations of the nuclei.

This requires the knowledge of the effects of these permutations on the coordinates of the system. These are given in ref. 20 (eq. 3.37-3.47) for the 6 elements of the permutation group S 3 of 3 particles. This group contains the identity operation, 3 permutations of 2 particles ( P12 , P13 and P23 for the permutation of particles 1 and 2, 1 and 3, 2 and 3, respectively) and 2 cyclic permutations. It has 2 one-dimensional irreducible representations (IRREP) and one two-dimensional one, which are labeled Γ = A

(3)

1 , A (3)
2 and E (3) . We add the unusual superscript (3) to these labels to distinguish them from those of the permutation group S 4 of 4 particles which are : 2 one-dimensional representations A 2 . More information on S 3 and S 4 , and in particular expressions of the matrices of their IRREPs, is given in ref. [START_REF] Hamermesh | Group theory and its applications to physical problems[END_REF], p. 224-226. S 4 is the full permutation group of the Ne 4 system and must be used to classify its vibrational states. The effets of the permutations on the ROHC are known for the 6 elements of S 4 which also belong to S 3 , but not for its 18 remaining other elements. Hyperspherical harmonics fully symmetrized with respect to S 4 are thus not available. Instead, we expand in the present work the Ne 4 vibrational states on the harmonics symmetrized with respect to S 3 and we show how to classify the resulting bound states according to the IRREPs of S 4 .

Using eq. 3.30, 3.37-3.47 of ref. 20 which give the effects of the permutations of S 3 on the coordinates of the system, we obtain :

P12 F nJ M J L M L D d (χ, a, θ, φ, δ) = (-1) n (-1) L F nJ M J L -M L D d (χ, a, θ, φ, δ) P13 F nJ M J L M L D d (χ, a, θ, φ, δ) = (-1) n (-1) L e i4π M L 3 F nJ M J L -M L D d (χ, a, θ, φ, δ) (22)
P23

F nJ M J L M L D d (χ, a, θ, φ, δ) = (-1) n (-1) L e i2π M L 3 F nJ M J L -M L D d (χ, a, θ, φ, δ)
The combinations :

F n η J M J L M L D d (χ, a, θ, φ, δ) = F nJ M J L M L D d (χ, a, θ, φ, δ) + ηF nJ M J L -M L D d (χ, a, θ, φ, δ) (23)
are even (η = (-1) n (-1) L ) or odd (η = (-1) n (-1) L+1 ) with respect to P12 . Furthermore, when M L is a multiple of 3, they are similarly even and odd with respect to P13 and P23 .

Thus, the harmonics

F n η J M J L M L D d belong to A (3)
1 when η = (-1) n (-1) L or to A

2 when η = (-1) n (-1) L+1 . When M L is not a multiple of 3, the pair of harmonics

F n +1 J M J L M L D d and F n -1 J M J L M L D d form
a basis for the 2 dimensional IRREP E (3) . In the following, we consider only the harmonics

F n η J M J L M L D d
with M L multiple of 3 and with η = (-1) L , consequently, we compute bound states which belong exclusively to A

C. Bound state computation

The method used here to obtain 4-body bound states is an extension of the 3-body roworthonormal hyperspherical harmonics expansion method already used to obtain Ar 3 [START_REF] Karlický | [END_REF], H + 3 [32], He + 3 [33] and Ar + 3 [34] rovibrational spectra. The Hamiltonian of the system is obtained by adding the 6 dimensional interaction potential of the system V (ρ, θ, φ, δ) to the kinetic energy operator given by eq. 9. We wish to compute a set of

N b eigen- states Ψ JM J ΠΓ k (ρ, χ, a, θ, φ, δ) (k = 1 -N b )
of this Hamiltonian with appropriate symmetry properties. These bound states are obtained by expansion on a basis of N s LHSF

Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ) (i = 1 -N s )
, which depend parametrically on ρ and which are defined as eigenstates of a partial fixed-ρ Hamiltonian :

Λ2 2µρ 2 + V (ρ, θ, φ, δ) Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ) = JΠΓ i (ρ)Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ) (24)
The LHSF potential energy curves JΠΓ i (ρ) are a useful tool to gain insight on the physical properties of the system [8,9] (see section III B below).

We consider a regular hyperradial grid ρ m of N ρ points and the associated d j (ρ) [35,36]. These are localized functions obtained from a basis of sine functions such that : d j (ρ m ) ∝ δ jm (δ jm being here the Kronecker symbol). We now form the N ρ ×N s product functions

(j = 1, N ρ ) Discrete Variable Representation (DVR) functions
d j (ρ)Φ JM ΠΓ i (χ, a, θ, φ, δ; ρ j ),
these are non-direct products since each LHSF is considered at the grid point ρ j associated to the DVR function d j (ρ). These product functions form a basis set on which to expand the bound states :

Ψ JM J ΠΓ k (ρ, χ, a, θ, φ, δ) = 1 ρ 4 i,j a k JΠΓ ij d j (ρ)Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ j ) (25)
The coefficients a k JΠΓ ij are obtained by diagonalizing the matrix representation of the Hamiltonian in this product basis. It is given after renormalization of the Hamiltonian by :

< d j (ρ)Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ j )|ρ 4 Ĥρ -4 |d j (ρ)Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ j ) >= - 2 2µ < d j (ρ)| d 2 dρ 2 |d j (ρ) > O JM J ΠΓ ij,i j + 6 2 µρ 2 j + JΠΓ i (ρ j ) δ ii δ jj (26)
δ ii and δ jj are usual Kronecker symbols, < d j (ρ)| d 2 dρ 2 |d j (ρ) > is the usual DVR matrix representation of the second derivative operator, as given for instance by appendix A of ref. 37. The overlap matrix elements are given by :

O JM J ΠΓ ij,i j =< Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ j )|Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ j ) >.
They replace the usual first derivative couplings which would appear if an adiabatic basis (see eq. 29 below) had been used to expand the wavefunction.

The surface functions Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ) are obtained from eq. 24 by expansion on a basis of hyperspherical harmonics

F n η J M J L M L D d (χ, a, θ, φ, δ) with quantum numbers selected according to the symmetry JM J ΠΓ considered. Numerical quadratures on a 5-dimensional grid N θ × N φ × N δ 1 × N δ 2 × N δ 3 are
performed for each hyperradial grid point ρ j to obtain the matrix elements of the potential in the harmonic basis. This task represents the most time consuming part of the computation, but it can be performed fully in parallel on different processors for different hyperradii. The matrix elements of the kinetic energy operator are obtained analytically from eq. 10.

In summary, we compute the vibrational states of Ne 4 with a basis expansion method which involves 3 successive diagonalization-contraction steps : hyperspherical harmonics are obtained from trigonometric functions, LHSF from harmonics, bound states from LHSF.

D. Implementation

We use the same interaction potential V (ρ, θ, φ, δ) as in ref. 12, this allows for a quantitative comparison between both sets of results. It is a sum of Lennard-Jones pair potentials v(r) which can be expressed as a function of the interatomic distance r by :

v(r) = 4 σ r 12 - σ r 6 ( 27
)
where =24.743 267 cm -1 is the 2 body binding potential energy and σ=5.195 au is a parameter related to the interatomic equilibrium distance r 0 by : r 0 = 2 1 6 σ. In the following, all data concerning energy will be given in units of . The 6 interatomic distances needed in the Ne 4 full potential expression are obtained from the Jacobi vectors, which are related to the ROHC by eq. 1. The Neon mass is, as in ref. [12], 19.99244 u.

The hyperspherical harmonics belong to an IRREP Γ of the subgroup S 3 of the full 4-body permutation group in our present implementation. We will present in the following results obtained for Γ = A

1 (Π = 0) and for Γ = A

2 (Π = 1). The total angular momentum is restricted to J=0. An irreducible representation of the group S 4 is a representation which may be reducible for the subgroup S 3 and we have the correlations :

A (4) 1 = A (3) 1 , A (4) 2 = A (3) 2 , T (4) 1 = A (3) 2 ⊕ E (3) , T (4) 2 = A (3) 1 ⊕ E (3) , (28)
For completeness, we recall the correlation : 3) . S 4 is isomorphic to T d and the naming conventions of the T d IRREPs can be used for S 4 . T

E (4) = E (
1 and T

2 are defined such that the character of the permutations of 2 particles i and j Pij (i, j = 1-4) is -1 for T (4)

1 and +1 for T (4)
2 . This convention is the one found in most textbooks (see for instance ref. [START_REF] Tinkham | Group theory and quantum mechanics[END_REF]

Γ = A (3) 1 (Γ = A (3) 2 ) must belong to the A (4) 1 or T (4) 2 (A (4) 2 or T (4) 1
respectively) IRREPs of S 4 . Indeed, the numerical diagonalization process generates automatically eigenstates which belong to IRREPs of the full symmetry group of the operator which is diagonalized, even if these states are expanded on a basis which belongs to a subgroup IRREP.

However, the full group IRREPs to which each eigenstate belongs automatically remains to be found. Usually, it is obtained by repeating the bound state calculations for all the subgroup IRREPs and by analyzing the degeneracies between the resulting bound states. This method was used in ref. 12 to obtain the full group IRREPs which the Ne 4 rovibrational states belong to from Jacobi coordinate calculations done for the IRREPs of the S 2 × S 2 subgroup, which contains 4 elements. In the present case, we could have generated the bound states for all S 3 IRREPs, not only Γ = A

(3) 1 which contains the physically allowed states, but also for Γ = A

(3) 2 and Γ = E (3) . We would know from eq. 28 that

A (3) 1 (A (3)
2 ) bound states which are degenerate with E (3) states belong to the full group IRREP T A distinctive advantage of the hyperspherical method is that the full group IRREPs which each bound state belongs to can be identified without extra calculations. Because the LHSF are obtained by diagonalization of the fixed ρ-Hamiltonian (eq. 24) which contains, similarly to the full Hamiltonian, all symmetries of the system, they also belong automatically to IRREPs of the full symmetry group, and not only of the subgroup associated to the hyperspherical harmonics on which they are expanded. Due to the symmetry invariance of the hyperradius ρ, a vibrational state belonging to a given (but at this time still unknown) full group IRREP of S 4 has non zero expansion coefficients a k JΠΓ ij (eq. 25) on a subset of the LHSF, the one which contains LHSF belonging to the same full group IRREP. By inspection of the expansion coefficient values (0 or not), we can separate the set of LHSF resulting from a single, say A Our convergence goal is to obtain bound state energies within 1 % of the reference data given by [12] and which turn out to be fully converged. This objective was achieved by choosing a truncated harmonic basis set such that n ≤ n max =50 and L ≤ L max =30 for Π = 0, Γ = A (3) 1 . A basis set of 20289 symmetry adapted harmonics was thus generated on a grid having the following size :

N θ = 51, N φ = 51, N δ 1 = N δ 2 = N δ 3 = 37.
The θ, φ, δ 2 grids are of Gauss-Legendre type while the δ 1 , δ 2 grids ones are regular Fourier type ones. N S =100 surface functions are obtained by diagonalization of the fixed-ρ Hamiltonian in the harmonic basis. A grid of N ρ = 36 DVR functions is used in the ρ interval [8.4,15] au. Finally, the full Hamiltonian representation in the product basis (eq. 25) is a matrix of size N ρ × N S =3600 which is diagonalized to provide the desired vibrational states and energies. Similar convergence parameters were used for the other symmetry considered

Π = 1, Γ = A (3) 2 .

III. RESULTS

A. Energy spectra

The energy levels resulting from our computations are shown on Table I. Our spectra are in good agreement with those of ref. 12, the relative difference being within 1 %. Our results have been obtained using an harmonic basis without any grid optimization. The harmonics which analytically regularize all the singularities of the kinetic energy operator [27] are fully delocalized in the configuration space. When the potential-induced localization of the LHSF in configuration space increases, the convergence rate of their expansion into harmonics decreases. Convergence is thus slower in the large ρ domain and for excited states which extend in this region. Grid optimization would be useful to speed up convergence and limit the CPU time spent in potential quadratures in our present implementation.

The states are labeled on Table I with respect to the IRREPs of S 4 . These labels are obtained by inspection of the weights a k JΠΓ ij in the expansions of the vibrational states into LHSF (eq. 25), according to the procedure described in subsection II D. Fig. 1 shows the weights integrated over ρ : j |a k JΠΓ ij | 2 , for different vibrational states (labelled by k), as a function of the surface function number i. There is a striking oscillatory dependence of the integrated weights as a function of i, for a given k. This is a symmetry effect : a weight a k JΠΓ ij is non zero only when the vibrational state Ψ JM J ΠΓ k and the LHSF Φ JM J ΠΓ i belong to the same full group IRREP. However, as the symmetry character of a given LHSF depends on ρ (there are avoided crossings between LHSF of different S 4 symmetries), it is somewhat blurred by integration over ρ and the alternation seen on Fig. 1 for integrated weights is not perfect. It is nevertheless sufficient to assign the vibrational state symmetries. Surface functions and vibrational states can be separated into 2 groups, the group which contains the surface function # 1 (ground state) is associated to the A 2 . Alternatively, one may prefer to look at the expansion coefficients at a chosen ρ instead integrating over ρ. Because the symmetry character of the surface function at fixed ρ is unambiguous, the alternation equivalent to that of Fig. 1 would be in this case perfect.

However, according to the ρ value chosen, the weight of a given LHSF on vibrational states may be very small, and the correlation between bound states and LHSF symmetries may be more difficult to identify.

The hyperspherical method has 2 other distinctive features which are useful to gain some physical insight. One is to provide surface function energy diagrams JΠΓ i (ρ) (eq. 24) which allow to follow the dynamics as a function of the hyperradius. The other is the separation of the 6 dimensional configuration space into S KI and S KR . We now show how to take advantage of these features to gain a better qualitative understanding of the Ne 4 spectra.

B. Potential and surface function energy curves

Fig. 1 shows that the expansion of a given vibrational state given by eq. 25 converges for a small number of surface functions. It can be seen on the figure that a single surface function is enough to account for 90 % (80 %) of the ground state of A (4)

1 (T (4)
2 ) symmetry. This justifies a posteriori the use of the hyperspherical adiabatic approximation to compute low energy vibrational states of rare gas clusters, as for instance in ref. [9]. Using now the LHSF basis as an adiabatic basis on which to expand the full wavefunction, eq. 25 turns into :

Ψ JM J ΠΓ k (ρ, χ, a, θ, φ, δ) = 1 ρ 4 i h JM J ΠΓ i (ρ)Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ) (29)
where the h Interestingly, the expansion of the excited vibrational states converges also with a small number of surface functions. For instance, 2 surfaces functions are sufficient to account for 90 % (85 %) of the first

A (4) 1 (T (4)
2 ) excited state at -2.961 (-2.885 ). The two-color segments drawn on Fig. 2 Asymptotically, the surface functions converge toward the Ne+Ne 3 dissociation limit and their energies are those of the trimer, -1.71 (Table II in ref. [START_REF] Yang | [END_REF], which corresponds to 0.57 of binding energy for each the 3 identical Ne-Ne bonds. The system can also dissociate into two dimers, Ne 2 +Ne 2 . This configuration corresponds however to a higher energy, -1.13360 [40], than for the trimer. At the largest ρ value shown on Fig. 2, ρ = 14.5 au, the energy of the lowest surface function is nearly -2.2 , which is significantly smaller than the asymptote -1.71 , the difference being related to the attractive nature of the long range part -r -6 of the Lennard-Jones potential. The range of ρ shown on Fig. 2 can be separated in 3 domains defined by : P SF : ρ 11, R SF : 11 ρ 13, T SF : 13 ρ (in au).

The superscript SF refers to the fact that these domains are defined with respect to the characteristics of the Surface Function energy curves i . Similar domains were defined on the ground surface function energy curve computed in ref. 9 and called "plateau-like regions".

The definition of these domains can be extended for higher energy surface function curves, but their boundaries must be shifted to increasingly larger ρ.

In the P SF domain, Π = 0 and Π = 1 surface function energies are degenerate by pairs. The potential curves have minima, near ρ = 9.6 au for the ground state, and near increasingly larger ρ for excited states. States of the same parity Π are non degenerate in this domain.

On the contrary, in the R SF domain, the Π = 0 -Π = 1 degeneracy is not present, but several states of the same parity are doubly degenerate. In this domain, energy curves have no minima but are convex. In the T SF domain, the energy curves are also convex and some are triply degenerate. II which can be formed from the Pyramid. 5 bonds provide the cohesion of the rhombus, its potential energy is minimum and is equal to nearly -5 when its side length is r 0 , which occurs when ρ = 2σ=10.390 au ∈ R min (see Fig. 3).

As the hyperradius increases further, the rhombus size increases, its energy also, but this increase is minimized by a slight distortion which consists in a decrease from π/3 of the vertex angle α (see Fig. 3). At some point, however, a Trigonal structure becomes energetically more favorable : 3 Ne atoms form an equilateral triangle from 3 identical bonds, a fourth bond is formed between one of these 3 atoms (say atom 3) and the 4th one, such that the (34) bond bisects bonds ( 13) and (23). We label this structure 4 -3 -12. There are 12 such equivalent trigonal structures (see Table III). The energy of these structures is nearly -4 when the common bond length is r 0 , which occurs when : ρ = 4 + √ 3σ=12.438 au ∈ T min (see Fig. 3). At some larger ρ, a structure more stable than this symmetric trigonal one can obtained by increasing the (34) bond length, but leaving the other 3 bonds unchanged, leading eventually to the asymptotic structure Ne+Ne 3 .

The abrupt transitions in the geometries of these constrainted fixed-ρ minimum potential energy structures as the hyperradius increases define 3 domains, P min , R min and T min , which have different bounds from those, P SF , R SF and T SF , observed on the surface function energy curves. We now analyze the surface functions themselves to understand how they reflect these structural changes as a function of ρ. This will allow us to understand why the ρ domain bounds are different for the surface function energies and for the minimum potential energy.

C. Surface function analysis

Analyzing the LHSF themselves is obviously difficult because they belong to the 5dimensional configuration space spanned by the 5 angles (θ, φ, δ 1 , δ 2 , δ 3 ). We can take advantage however of the separation of this space into the S KI and S KR subspaces.

1. In S KI S KI contains θ and φ and its structure is depicted in ref. [15], Fig. Fig. 4 shows ground and excited LHSF in S KI for ρ=9.6 au. These functions are all localized in S KI near the P point (tetrahedron), as expected since ρ=9.6 au ∈ P SF (and incidentally also ∈ P min ). The plots are performed somewhat arbitrarily for : δ 1 = δ 2 = δ 3 = π/2 since the wavefunction value is strictly independent of δ 1 , δ 2 , δ 3 at the P point, and weakly dependent in its vicinity. The localization of the wavefunctions near P results from the presence of a potential barrier near coplanarity (φ = 0). This barrier prevents the system from changing chirality during its vibration as it would if it could pass through the coplanar geometry. As a result, Π = 0 wavefunctions are degenerate with those for Π = 1 and the plots for Π = 1 would be identical to those of Fig. Fig. 5 and6 show the changes of the ground LHSF characteristics as a function of ρ. The 3 ρ values are selected such that : ρ=10.6 au ∈ P SF , ρ=11.6 au ∈ R SF and ρ=13.6 au ∈ T SF , so that the gradual change in S KI (left column) from pyramid to rhombus and trigonal configurations is visible on the figures. At ρ=10.6 au (Fig. 5 top left), the wavefunction is still mainly localized in the vicinity of P , as for ρ=9.6 au (Fig. 4), however a transition toward coplanarity starts to appear by an extension of the wavefunction with significant amplitude in this direction (φ = 0). At ρ=11.6 au (Fig. 5 bottom left), the surface function already has a significant amplitude in the coplanarity region as a result, the degeneracy between Π = 0 and Π = 1 pair of states disappears.

It is interesting to notice that the wavefunction at ρ = 10.6 au is not localized in the minimum potential region (φ = 0 line = coplanar configurations), but in the low kinetic energy region. In fact, the wavefunction remains localized far from coplanarity until the potential near P becomes larger than the surface function energy. This is shown on Fig. 3 : the transition between P SF and R SF occurs near ρ =11 au when the potential at P (dashed line) becomes larger than the lowest surface function energy. The P point region then becomes energetically unreachable for the wavefunction, which is pushed toward the coplanar region by the potential as ρ increases. The fact that P SF extends up to 11. au, as compared to P min which extends up to 10. au only, is thus related to the fact that the wavefunction localization is governed by the turning point region, rather than by the minimum potential region. A similar phenomenon occurs for the transition from R SF to T SF , between rhombus and trigonal configurations. The transition occurs when the potential associated to the rhombus configuration (dashed line, α=60 • ) crosses the lowest surface function energy curves near ρ=13. au (Fig. 3) so that the rhombus region becomes energetically unreachable, and not near ρ=12. au (the boundary between R min and T min ) when the potential minimum switches from rhombus to trigonal geometries. Because ρ = 11.6 au ∈ R SF and simultaneously ∈ R min , the corresponding wavefunction (Fig. 5, left column, bottom row) is localized near coplanarity, which is also the region of minimum potential. A similar observation can be made in the trigonal region for the wavefunction at ρ = 13.6 au ∈ T SF and simultaneously ∈ T min (Fig. 6, left column, bottom row). In summary, we can say that the LHSF localization is dynamically (turning points driven) rather than statically (potential driven) controlled.

In S KR

Plots of the LHSF in S KR provides complementary useful information (Fig. 5, central and right columns, Fig. 6, right and Fig. 7, right). Coplanar configurations of special interest are the 6 rhombus and 12 trigonal ones already mentioned, their coordinates have been obtained with the method described in appendix A and are listed in Tables II andIII.

The plots are presented in specific planes of S KR selected such that they contain several of these important configurations and such that the surface function amplitude is significant.

The ground surface function is shown on Fig. 5 and6. The most striking feature is the localization of the functions near the 6 points R 1 -R 6 at ρ = 11.6 au and near the 12 points T 1 -T 12 at ρ = 13.6 au. This localization was expected since ρ=11.6 au ∈ R SF and ρ=13.6 au ∈ T SF . It is important to realize that this localization is the reason why surface functions are degenerate in the R SF and T SF domains (Fig. 2 and3). Indeed, let us call Φ Rm (m = 1 -6) (Φ Tm (m = 1 -12)) wavefunctions which are localized exclusively in the vicinity one of these 6 (12) points. The symmetrized LHSF shown on the figures are linear combinations of these. Indeed, the 6 Φ Rm are degenerate in R SF , the 12 Φ Tm in T SF , and each of these is, like the LHSF Φ JM J ΠΓ , an eigenfunction of the fixed ρ-Hamiltonian, but without the appropriate permutation symmetry properties. We show in appendix B how to use group representation theory to form the appropriate linear combinations of these localized functions which restore the LHSF with the expected permutation symmetries. We also show that this S KR localization explains the appearance of degenerate pairs (triplets) in the R SF (T SF ) domain. As surface function localization decreases as ρ decreases (compare for instance the surface functions in S KR at ρ = 10.6 au, Fig. 5 top row, vs ρ =11.6 au, Fig. 5 bottom row), S KR -induced degeneracies disappear in the P SF domain (but S KI -induced ones appear, see section III C 1). Fig. 7 also shows excited surface functions in S KR . Surface functions at the same ρ may have different energies but similar shapes in S KR , as if they had similar excitation in this subspace. For instance, the surface function # 3 at ρ =11.6 au (Fig. 7, top-right) is similar in S KR to function # 1 (Fig. 5, bottom right), as well as to functions #2 and #4 (not shown). In this case, excitation is visible in S KI (similar to that shown on Fig. 4 at ρ =9.6 au). In other cases, excitation also appears in S KR , as shown on Fig. 7 (bottom right) for state # 5 which has a δ 1 excitation. We show in appendix B that the state of excitation has impact on the symmetry properties and on the degeneracies of the surface functions, for instance states #1-4 without excitation in S KR are doubly degenerate in R SF whereas state # 5 with δ 1 -excitation is non degenerate.

D. Bound state analysis

As the bound states are linear combinations of the LHSF (eq. 29) which involve a few termes only (Fig. 1), the properties of the bound states are directly related to those of the surface functions. Moreover, at least for the low energy states, the coefficients of this expansion are localized in ρ so that the properties of the bound state reflect those of the surface functions in a limited number of domains. For instance, bound state # 1 (for both parities) is localized in the P SF domain (see Fig. 2). As a result, the bound state wavefunction in localized near P in S KI (see Fig. 4) and is associated to the tetrahedral geometry.

States from #2 to #10 have energies below -2.55 (-2.50 ) for Π = 0 (Π = 1), see Table I. These states are localized in P SF + R SF , the corresponding motion is one of a pyramid experiencing large amplitude rotations around its edges (say the one connecting atoms 1 and 2) so that the equilateral triangles ( 123) and (124) can explore the coplanar geometries and in particular the rhombus one R 1 = 3 -12 -4. The motion extends symmetrically (Π = 0) or antisymmetrically (Π = 1) on the other side of the rhombus plane and the tetraedron geometry with opposite chirality can be explored. This vibrational motion between tetrahedron and rhombus configurations corresponds to an oscillation motion in ρ between P SF and R SF described by the h JM J ΠΓ i (ρ) functions (eq. 29). The same rotational motion can occur around the 6 edges of the tetraedron, which corresponds to exploring the 6 rhombus configurations R 1 -R 6 (Table II) with identical, opposite or different amplitudes, according to the possible symmetries of the wavefunction in S KR . As long as the bound states are in P SF , they do not explore the coplanarity region and tunneling splitting is very small. This splitting increases however for higher energy bound states which extend in the R SF region which corresponds to coplanarity. At even higher energies (state #11 and above) the motion can be even more delocalized in P SF + R SF + T SF , which corresponds to an even larger amplitude motion in ρ. In this case, an in-plane motion, where the rhombus is distorted into the trigonal geometry, adds to the out of plane motion between tetrahedron and rhombus configurations. Tunneling splitting remains large.

A splitting of the Ne 4 vibrational spectrum in 3 regions is proposed in ref. [12]. It is based on the analysis of the increase of the tunneling splittings and of the density of vibrational states (DOS) with energy. The low energy region contains a few low lying states (states # 1 of both parities on Table I) which are "highly localized". The intermediate energy region corresponding approximately to the interval [-3.4 , -2.6 ] (from states # 2 to #8 on Table I) is associated to an intermediate DOS of "fluxional" states with intermediate tunneling splittings. The high energy region (above -2.6 ) corresponds to a high DOS of completely delocalized vibrational states. The present hyperspherical analysis provides a complementary view on this spectrum splitting. From inspection of Fig. 2 and3 in relation to the energy bounds of these 3 regions, it is clear that the low energy region corresponds to the states which extend in P SF , the intermediate energy region to the states which extend in P SF +R SF and finally the high energy region to the states which extend in P SF +R SF +T SF .

The increase with energy of the vibrational DOS mentioned in ref. [12] is thus related to an increase of the configuration space volume occupied by vibrational states of increasing energy. The increase of the tunneling splitting is related to the accessibility of the coplanar region in R SF and T SF .

IV. SUMMARY AND CONCLUSIONS

We have presented in the present paper the first implementation of the row-orthonormal hyperspherical coordinate formalism on a tetratomic system. The method provides not only spectra with reasonable accuracy, but also physical insight into the vibrational dynamics of the system. In particular, we have related features of the spectra, like permutation and inversion symmetry degeneracies, inversion tunneling splittings and densities of states, to the way the wavefunction is localized in configuration space. There is however an obvious shortcoming in our approach which should be remedied in the future. We have used the subgroup S 3 of the full permutation group S 4 of the system. As a result, our computed spectra are superpositions of spectra associated to different irreducible representations of S 4 . Although we have shown how to disentangle simply mixed symmetries within the hyperspherical formalism, performing separate computations for the different irreducible representations of S 4 would reduce the computational cost of the calculations. This requires first to identify the effect of all the 24 elements of S 4 on the coordinates of the system and to symmetrize accordingly the Wigner rotation matrices on which the hyperspherical harmonics are expanded. without sign change of the amplitude. Using Table II, we obtain for instance:

P12 Φ R 1 = Φ R 1 , P23 Φ R 1 = Φ R 3 , P34 Φ R 1 = Φ R 1 (B1)
and similarly for the other Φ Rm functions. The effect of the other permutations of S 4 can be obtained by combination of these 3 generators of S 4 . We thus obtain the 6 dimensional reducible representations of the 3 generators of S 4 in the Φ Rm basis :

Γ R ( P12 ) =
1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 Γ R ( P23 ) = 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 Γ R ( P34 ) = 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1

The IRREPs of S 4 , 1-dimensional ones A

1 and A

2 , the 2-dimensional one E (4) and the 3-dimensional ones T (4) 1 and T (4) 2 , are well known (ref. 30, p. 224-226). Using these and the orthogonality relation for characters (eq. 3.20, p. 30 in ref. [START_REF] Tinkham | Group theory and quantum mechanics[END_REF], we obtain the following

decomposition of Γ R : Γ R = A (4) 1 ⊕ E (4) ⊕ T (4)
2 . Because of the known correlation between S 3 and S 4 IRREPs, we know that only 1 of the 3 basis functions belonging to T We now consider a function with a δ 1 excitation, like the state # 5 shown on Fig. 7 (ρ = 11.6 au). Although it is like states # 2 and 4 a T and is non degenerate.

In the T SF domain, the surface functions are localized in the vicinity of the trigonal points T 1 -T 12 (see for example Fig. 6, right column) and they can be obtained as linear combinations of the localized functions Φ Tm (m = 1 -12). These form a 12-dimensional representation Γ T of S 4 which can be reduced into its IRREPs using the same techniques as in the R SF domain. We obtain :

Γ T = A (4) 1 ⊕ E (4) ⊕ 2T (4) 2 ⊕ T (4)
1 . As only A

1 , we expect a triple degeneracy for the lowest energy states in the T SF , as observed on Fig. 3, which consists of a single A and is non degenerate in the T SF domain.

Notice finally that the degeneracies observed in the R SF and T SF domains result from surface function localizations in S KR . In the P SF domain, surface functions are delocalized in S KR and these degeneracies disappear. Instead, new degeneracies between Π = 0 and Π = 1 surface functions appear in this domain and they are induced by localization of the surface functions near P in S KI (see section III C 1). The present results are compared with those of ref. [12] IRREPs is different from the one used in ref. [12] (see section II D). d and the kinematic invariant angles by : θ R = α/2 and φ R = 0 (coplanar configuration).

A (3) 1 , Π=0 A (3) 2 , Π=1 IRREP E (present) E ref (ref. [12]) E -E ref IRREP E (present) E ref (ref. [12]) E -E ref 1 A (
In the Table : δ R = arcsin 1 √ 3 ≈ 35 • . The 4 configurations i -jk -l, l -jk -i, i -kj -l and l -kj -i correspond to the same set of coordinates (ρ, θ, φ, δ 1 , δ 2 , δ 3 ). Appendix A explains how these coordinates have been obtained. and the kinematic invariant angles θ and φ are given by : θ T = arcsin

Rhombus # Configuration δ 1 , δ 2 , δ 3 R 1 3 -12 -4 π 2 , π -δ R , π 2 R 2 1 -23 -4 π 2 , δ R , 5π 6 R 3 2 -13 -4 π 2 , δ R , π 6 R 4 2 -14 -3 π -δ R , π 2 , π 3 R 5 1 -24 -3 δ R , π 2 , 2π 3 R 6 1 -34 -2 δ R , π 2 , 0
1 √ 4+ √ 3
≈ 25 • and φ T = 0 (coplanar configuration).

The values of the angles used in the Table are:

δ 1 T = π 4 -δ R 2 = arcsin 1 2 -1 2 √ 3 ≈ 27 • , δ 2 T = 3δ 1 T = arccos 1 2 -5 6 √ 3 ≈ 82 • , δ 3 T = arccos -30 37 + 2 √ 3 37 ≈ 162 • , δ 4 T = arccos 1 6 + 1 6 √ 3 ≈ 59 • , δ 5 T = arctan -1 + 1 √ 3 ≈ 122 • .
δ R is defined in the caption of Table II. The 2 configurations i -j -kl and i -j -lk correspond to the same set of coordinates ρ, θ, φ, δ 1 , δ 2 , δ 3 . Appendix A explains how these coordinates have been obtained.

1 permutation symmetry. Right : Π = 0, T II for the definition of δ R . The point labelled P corresponds to the tetraedral (Pyramid) configuration, the point S to coplanar symmetric top configurations and the line φ = 0 to all coplanar configurations. Right : results in the (δ 1 , δ 3 ) plane, θ = π 6 , φ = 0, δ 2 = π 2 .

  two-dimensional representation E(4) and 2 three-dimensional ones T

  that those which are not belong to A

  . Such an IRREP identification method thus usually requires the computation of extra bound states which are not of direct interest.

2 .

 2 into 2 subsets, one for each full group IRREP, say A The subset which corresponds to A (4) 1 is the one which contains the lowest energy LHSF. The full group IRREP to which a given bound state belongs is now known : it is the one of the LHSF subset on which it has non zero expansion coefficients. Neon is a boson with zero nuclear spin and the only physically relevant bound states belong to the full group IRREP A(4) 1 and are obtained from a single Γ = A (3) 1 calculation. Extra calculations were performed for Γ = A

  physical insight and to allow for a larger comparison with the results of ref.12.

  JM J ΠΓ i (ρ) functions are the coefficients of the expansion in the adiabatic representation. A single product h JM J ΠΓ i (ρ)Φ JM J ΠΓ i (χ, a, θ, φ, δ; ρ) is enough to represent the lowest energy bound state of each symmetry and in this case the hyperradial function h JM J ΠΓ i (ρ) is the one-dimensional bound state supported by the potential energy curve JΠΓ i (ρ) associated to Φ JM J ΠΓ i . Fig. 2 shows a set of such energy curves JΠΓ i (ρ) as a function of the hyperradius, as well as a few low-lying bound state energies.

  are bounded by pairs of potentials JΠΓ i (ρ) of the same two colors. These potentials are associated to the pair of surface functions involved in the expansions of these vibrational states. Fig. 1 also shows that a few surface functions are sufficient to converge each of the more excited vibrational states. This has the important implication that most of the qualitative characteristics of a given vibrational state are embedded in a few surface functions, those which are involved in the expansion of this vibrational state. This is why it is interesting to consider further the properties of the surface function energy diagram shown on Fig. 2.

R min : 10 ρ 12 ,

 12 figure also shows the potential energy minimum curve. At fixed ρ, minimization of the potential is performed over the 5 angles θ, φ, δ. Minima and curvature effects similar to those associated to i can be observed on this minimum potential energy curve, but the bounds of the corresponding domains are shifted in ρ. We now have : P min : ρ 10, R min : 10 ρ 12, T min : 12 ρ (in au). In each domain, the potential minimum corresponds to a different geometry. P min is the domain where the potential minimum corresponds to a tetrahedron (a triangular Pyramid) geometry. Its 6 edges are formed from

1 . 1 √ 1+cos 2 φ

 112 It contains a point P given by θ = arcsin2 3 ≈ 55 • , φ = π 4 which corresponds to the same tetrahedral configuration (symmetric top, Pyramid) for all δ 1 , δ 2 , δ 3 values. The point θ = 0 (arbitrary φ) corresponds to the collinear configurations, the point S (θ = π/4, φ = 0) to the coplanar symmetric top configurations, different configurations of each type corresponding to different δ 1 , δ 2 , δ 3 , i.e. different points in S KR . One example of such coplanar symmetric top is the square formed by 4 Neon atoms, another one is the equilateral triangle formed by 3 of them, with the 4th one placed at the center of mass of the triangle. The segment which connects P to S and defined by : θ = arcsin contains all oblate symmetric top configurations, like the triangular pyramids derived from the tetrahedron by approaching one apex Ne atom perpendicularly to its equilateral triangular base.

  4. Fig. 4 also shows the nodal lines associated to vibrational excitations. The P S line is an antinodal (nodal) line for the A . A motion along the P S line corresponds to vibrations of the Ne atoms which preserve the symmetric top character. The nodal lines which are approximately perpendicular to the P S line in S KI (bottom row of Fig. 4 for both permutation symmetries) are other traces of vibrational excitation.

1 = 1 √ 6 m=1- 6 Φ 2 = 1 √ 6 m=1- 3 Φ

 11662163 pairs (functions #1-2 and #3-4) observed on Fig. 3 in R SF . Using the projector operator technique (eq. 3.36, p. 41, in ref. 38), these functions are expressed in terms of the basis function of the reducible representation as : Φ A (4) Rm and Φ T (4) Rm -m=4-6 Φ Rm . This explains why the amplitude of the A (4) 1surface function of a degenerate pair in R SF is the same for all rhombus points R 1 -R 6 , while the amplitude near R 4 -R 6 is opposite to the one near R 1 -R 3 for the other surface function of T

1 ,

 1 has significant and similar amplitudes only in the vicinity of R 4 -R 6 and negligible ones near R 1 -R 3 . Also, there is no corresponding A (4) 1 state to form a degenerate pair with. These observations can also be accounted for by group theory considerations. This surface function with δ 1 excitation is obtained by linear combination of new Φ Rm (m = 1 -6) basis functions, which are localized similarly to the Φ Rm ones but differ from them by a δ 1 excitation similar to that of the surface function. A function Φ Rm associated to the rhombus m = i -jk -l (see TableII) behave like Φ Rm under the action of P12 , P23 and P34 , except for the central ones : P12 Φ R 1 = -Φ R 1 ... We thus obtain the following new 6 dimensional reducible representations Γ R of the 3 generators of S 4 in the Φ Rm basis : this explains why state # 5 in R SF belongs to T (4) 2

  belonging to A

  #4 has a δ 1 excitation in the T SF domain. We use localized functions Φ Tm (m = 1 -12) with appropriate excitations such that they are unchanged by permutation of 2 atoms except the pair of end atoms kl in the trigonal structure i -j -kl (TableIII): P12 Φ T 1 = -Φ T 1 ...The 12 dimensional reducible representation Γ T of S 4 in the Φ Tm is likewise reduced into its IRREPs : Γ T = A

Trigonal geom. # Configuration δ 1 , δ 2 , δ 3 T 1 4 -FIG. 1 :

 1341 FIG. 1: Integrated weight of the surface function #i (Φ JM J ΠΓ i

1 FIG. 2 : 3 ) 1 , 3 ) 2 ,

 123132 FIG. 2: Surface function energies (unit : , the Ne 2 pair-potential well depth) as a function of the hyperradius (atomic unit). The 0-energy corresponds to 4-body dissociation. They result from 2 calculations, Γ = A (3) 1 , Π = 0 (full lines) and Γ = A(3) 2 , Π = 1 (symbols). S 4 symmetry labels are shown for the four lowest potentials. The 4 lowest bound state energies (Π = 0) are shown as horizontal segments. The lowest energy bound state (for each A (4) 1 and T (4) 2 symmetry) is well represented by the adiabatic approximation and is supported almost entirely by a single surface function potential, the color of segment is the one of the corresponding potential. As at least 2 surface functions are necessary in the expansion of the first excited states (in each symmetry), their energies are represented by segments with 2 colors identical to those of the corresponding potentials.

FIG. 3 : 7 FIG. 4 : 2 .FIG. 5 :

 37425 FIG.3: Potential energy minimum (blue continuous line) as a function of ρ. 0-energy corresponds to 4-body dissociation. The dashed lines correspond to the potential energies for selected geometries of the system : tetrahedron, rhombus (vertex angle fixed : α = π3 or optimized to minimize potential), trigonal (bond angles with respect to the central atom : 5π 6 , 5π 6 , π 3 ). These geometries correspond to the potential minimum in the vicinity of ρ = 8.998 au, ρ = 10.390 au and ρ = 12.438 au, respectively. The lowest Π = 0 A (4) 1 and T (4) 2 surface function energies are also shown as black continuous lines. The excitation numbers of the surface functions plotted on Fig. 4-7 are also indicated.

  .

3 FIG. 6 :FIG. 7 :

 367 FIG.6: Surface function amplitude color plots for the ground state surface function for Π = 0 (A (4) 1 permutation symmetry) at ρ = 13.6 au. The color scale gives the amplitude from blue (negligible amplitude) to red (max amplitude, given by the value of M on the graph). The white lines are contour plots of the potential with labels in units of .Left : plot in the (θ, φ) plane for fixed values of the other coordinates :δ 1 = π 2 , δ 2 = δ 1 T , δ 3 = π 2 .The point labelled P corresponds to the tetraedral (Pyramid) configuration, the point S to coplanar symmetric top configurations and the line φ = 0 to all coplanar configurations. Right : plot in the (δ 2 , δ 3 ) plane for δ 1 = π 2 , θ = θ T , φ = 0. See the caption of TableIIIfor the definitions of δ 1T and θ T . A thick white bullet in the (θ, φ) plot (left frame) indicates the values of θ and φ chosen for the plots in the (δ 2 , δ 3 ) plane. Similarly, a thick white bullet in the (δ 2 , δ 3 ) plot (right frame) indicates the values of δ 2 and δ 3 used in the corresponding (θ, φ) plot. The labels T 1 -T 6 correspond to the coordinates listed in TableIII.

  , p. 329), albeit opposite to that of ref. 12 following ref. 11. Eq. 28 implies that a Ne 4 state computed for

TABLE I :

 I Ne 4 energy levels in unit of for the Γ = A

  (Table VII). The IRREPs of S 4 to which each state belongs are also given. For each of the T Table VII), only the lowest of both is given in the present Table. The present naming convention of the T

		(4) 1	or T 2 (4)	states, 2 almost equal energies are
	given in ref. [12] ((4) 1	(4) 2 and T

TABLE II :

 II Row orthonormal hyperspherical coordinates (δ 1 , δ 2 , δ 3 ) for R 1 to R 6 rhombuses which are labelled i -jk -l, where i, j, k, l are Ne atom numbers. Atoms i and l form the pair of opposite vertices connected by the longest diagonal, j and k by the shortest one. The angle α at vertices i and l is thus such that α ≤ π/2. Calling d the side of the rhombus, the hyperradius is given by

	ρ = 2	5 6
			4) 1	-3.454	-3.4494	-0.005	A	(4) 2	-3.454	-3.4494	-0.005
	2	T 2 (4)	-2.984	-2.9844	0.000	T 1 (4)	-2.983	-2.9830	0.000
	3	A	(4) 1	-2.961	-2.9638	0.003	A	(4) 2	-2.957	-2.9586	0.002
	4	T 2 (4)	-2.885	-2.9073	0.022	T 1 (4)	-2.847	-2.8591	0.012
	5	A	(4) 1	-2.867	-2.8890	0.022	A	(4) 2	-2.808	-2.8224	0.014
	6	T 2 (4)	-2.778	-2.7913	0.013	T 1 (4)	-2.694	-2.7093	0.015
	7	A	(4) 1	-2.733	-2.7450	0.012	A	(4) 2	-2.677	-2.6840	0.007
	8	A	(4) 1	-2.652	-2.6573	0.005	A	(4) 2	-2.604	-2.6158	0.012
	9	T 2 (4)	-2.605	-2.6171	0.012	T 1 (4)	-2.539	-2.5490	0.010
	10 A	(4) 1	-2.571	-2.5836	0.012	A	(4) 2	-2.515	-2.5249	0.010
	11 T 2 (4)	-2.529	-2.5459	0.016	T 1 (4)	-2.494	-2.5065	0.012
	12 A	(4) 1	-2.510	-2.5243	0.014	T 1 (4)	-2.472	-2.4862	0.014
	13 T 2 (4)	-2.494	-2.5162	0.022	T 1 (4)	-2.431	-2.4564	0.025
	14 T 2 (4)	-2.472	-2.5090	0.037	A	(4) 2	-2.415	-2.4291	0.014
	15 T 2 (4)	-2.458	-2.4794	0.021	A	(4) 2	-2.404	-2.4230	0.019
	16 A	(4) 2	-2.425	-2.4473	0.022	T 1 (4)	-2.395	-2.4221	0.027

TABLE III :

 III Row orthonormal hyperspherical coordinates (δ 1 , δ 2 , δ 3 ) for T 1 -T 12 trigonal geometries. These configurations are labelled i -j -kl, where i, j, k, l are Ne atom numbers. Atoms jkl form an equilateral configuration and atom i is bonded to the central atom j. Bond distances ji, jk, jl and kl are the same (= d), however these trigonal configurations are not fully symmetric as ji, jk and jl, ji bond angles are 150 • whereas the jk, jl bond angle is 60 • . For all trigonal configurations of the Table, the hyperradius is given by ρ = d

	(4+ √	3)	1 2
	2	1 6

Acknowledgments

We thank prof. Christoph Meier for stimulating discussions.

Data availability

The data and computer programs that support the findings of this study are available from the author upon reasonable request.

This is a natural continuation of the present work.

Appendix A: The hyperspherical coordinates for rhombus and trigonal geometries Finding the hyperspherical coordinates for a given geometry of the system requires solving non linear equations to obtain them from the mass-scaled Jacobi coordinates in the body frame : ρ bf = ρN(θ, φ) R(δ 1 , δ 2 , δ 3 ). Solving these equations is in general a difficult task but for specific geometries simple solutions exist. For the rhombus geometries, we first consider the configurations where atoms 1 and 2 are located at opposite vertices, namely R 1 : 3-12-4 and R 6 : 1-34-2. In these cases, the Jacobi vectors r (2) and r (3) are collinear and both perpendicular to r (1) , this special geometries make the extraction of the hyperspherical coordinates easier. The hyperspherical coordinates for all other rhombus configurations could be obtained from a single one, say R 1 , if the effects on coordinates of all atom permutations belonging to S 4 were known. Unfortunately, simple expressions are available only for the transformations which belong to the subgroup S 3 of S 4 . The permutations of S 3 must be applied on 2 configurations, say R 1 and R 6 , to cover all 6 configurations. The hyperspherical coordinates for R 2 are obtained from those for R 1 by applying the coordinate transformation associated to the permutation of atoms 1 and 3 and given by eq. 3.40 of ref. [20].

Similarly, the hyperspherical coordinates for R 3 are obtained from those for R 1 by applying the coordinate transformation associated to the permutation of atoms 2 and 3 and given by eq. 3.38 of ref. [20]. A similar procedure can be used to obtain the coordinates associated to R 4 and R 5 starting from R 6 . The results for the 6 rhombus configurations are collected on Table II.

A similar procedure can be used to obtain the hyperspherical coordinates of the 12 trigonal configurations. T 1 , T 4 and T 12 are special configurations for which the hyperspherical coordinates can be obtained easily from the Jacobi vectors. T 2 and T 3 coordinates are obtained from T 1 by permutations, T 5 and T 6 from T 4 ; the 5 configurations T 7 to T 11 are similarly obtained from T 12 . The results for the 12 trigonal configurations are collected on Table III.

Appendix B: Group theory, degeneracy and symmetry of surface functions

We discuss the degeneracies of the surface functions observed on Fig. 2 and3 : double degeneracies appear in the R SF domain (rhombus) and triple degeneracies in the T SF domain (trigonal). We relate these degeneracies to the symmetry properties of the surface functions. In the T SF domain near ρ=13.6 au, the 3 lowest surface functions are degenerate (see Fig. 3) and among them, 2 have the T We now show how these degeneracy and symmetry properties can be accounted for using group theory. We take advantage of the fact that the surface functions are linear combinations of disconnected wavefunctions localized in the (δ 1 , δ 2 , δ 3 ) subspace near the rhombus R 1 -R 6 or trigonal T 1 -T 12 points (see for example Fig. 5, central and right columns, Fig. 6, right column, as well as Fig. 7, right column). We call Φ Rm (m = 1 -6) and Φ Tm (m = 1 -12) these localized wavefunctions. They are degenerate solutions of the fixed-ρ Hamiltonian, but on the contrary of the surface functions Φ JM J ΠΓ i which are also solutions of eq. 24, they are not basis functions of an IRREP of the permutation group S 4 . However, a 6-dimensional (12-dimensional) reducible representation of this group can be obtained from the Φ R i (i = 1 -6) basis functions (Φ T i (i = 1 -12)) which we call Γ R (Γ T ) and which can be decomposed into the IRREPs of S 4 using the usual projection techniques ( [START_REF] Tinkham | Group theory and quantum mechanics[END_REF], p. 29-43).

We first assume that the Φ Rm functions are nodeless in the (δ 1 , δ 2 , δ 3 ) subspace. This is the case for instance for the functions # 1 to # 4 at ρ=11.6 au (function # 1 is shown on Fig. 5 and #3 on Fig. 7). In this case, applying Pij on Φ Rm changes it into another Φ R m


O Adrot
email: oadrot@ensem.u-nancy.fr

D Maquin
email: dmaquin@ensem.u-nancy.fr

J Ragot
email: jragot@ensem.u-nancy.fr

FAULT DETECTION WITH MODEL PARAMETER STRUCTURED UNCERTAINTIES

Keywords: Fault detection, Parity space, Structured uncertainties, Intervals, Enclosures

This paper focuses on an original fault detection method, able to take into account parameter uncertainties described by bounded variables. Intervals offer an attractive alternative for uncertain modelling and allow a straightforward generation of adaptive thresholds. This method is based on the parity space approach. The original point treated herein concerns the choice of the parity matrix, according to uncertain matrices of state representations. After expressing external forms of parity relations, interval arithmetic is used to evaluate internal forms and to deduce residual enclosures delimiting normal operation field. Then, we stress on the problems caused by dependencies between interval variables.

INTRODUCTION

Fault detection schemes often use the concept of analytical redundancy. They are based on consistency tests where sensory observations of a physical system are confronted with the information embodied in its model. Then residuals are generated and a normal operation procedure is usually made. An inconsistency is revealed when at least one residual becomes non-zero or one estimated parameter abnormally deviates. Reviews of model-based fault detection and isolation (F.D.I.) have been published for twenty years: [START_REF] Willsky | A survey of design methods for failure detection in dynamic systems[END_REF], [START_REF] Frank | Fault Diagnosis in Dynamical Systems Using Analytical and Knowledge based Redundancy -a Survey and Some New Results[END_REF], [START_REF] Gertler | Analytical redundancy methods in fault detection and isolation[END_REF] for state estimation; [START_REF] Isermann | Process fault detection based on modeling and estimation methods -a survey[END_REF] for parameter identification.

A problem met in the field of F.D.I. schemes lies in the fact that a model only defines an approximate behaviour of a physical system. This is caused by modelling errors when a model is made linear or when some physical phenomena are not taken into account. However, a modelling error has not to be identified with a fault. Some methods generate robust residuals using uncertainties de-coupling techniques. These are considered as disturbances whose distribution matrices are well-known. Some investigators have applied the unknown input observer (U.I.O.) or have built observers with the assignment of eigenstructure scheme [START_REF] Patton | Fault Diagnosis in Dynamic Systems[END_REF]. The major drawback is that diagnosis is only made on the certain part of the model, without taking into account the information embodied in the eliminated uncertainties. In some other methods, uncertainties receive a stochastic description and are described as an additive random noise. If statistical hypotheses about parameters are made (multiplicative uncertainties), the problem becomes insoluble because of operations on probability density functions.

This paper focuses on a fault detection method taking into account the structure and the amplitude of uncertainties. Each uncertain parameter is described by a bounded variable (or interval variable). This representation is used in parameter identification [START_REF] Milanese | Bounding approches to system identification[END_REF] but is not widely-known in diagnosis [START_REF] Chang | Model-Based Approach for Fault Diagnosis. 2. Extension to Interval Systems[END_REF]. A major drawback of interval arithmetic being its explosive nature in the case of recursive systems (because of dependencies between interval variables), a parity space approach has been chosen. In this way, the dynamic equations of a model can be formulated in the form of algebraic relations. Sensory observations are stacked on a finite time horizon and a static representation is obtained. A procedure allows to optimise the parity matrix W to obtain a total or partial de-coupling of residuals from state variables.

Then, the enclosures of residuals are built with the help of interval arithmetic. These take into account the ranges of uncertain parameters and define the bounds of the normal operation field. This approach is comparable with stochastic methods, where thresholds depend on the no-detection and false alarm rates as well as the statistic characteristics of the measured signals. Nevertheless, interval variables are interesting because they allow to straightforwardly describe parameter uncertainties in a model and then thresholds are built more naturally.

Thus, this paper has the following structure. The general expressions of the static forms are defined in the next section. In section 3, a procedure allowing to choose the parity matrix W is presented. Internal forms using this result are given with respect to the uncertain matrix of the state representation in section 4. At last, interval arithmetic is presented to evaluate previous internal forms and build the enclosures of residuals.

STATIC FORMS

2.a. First static form

Structured uncertain models allow to take into account the lack of knowledge on a physical system by indicating which parameters are uncertain. We only consider dynamic systems given by linear discrete state equations. The structured uncertainties are described by additive terms, which modify the matrices A, B, C of the next model:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) X k A X k B U k X k Y k C X k Y k U k A B n C m q + = + ∈ = ∈ ∈      1 θ θ θ r r r , , (1)
where: ( )

{ } G G G G A B C G G i i θ θ = + ∈ ∑ 0 , , , , .
U(k), X(k) and Y(k) respectively define the state, actuator input and sensor output vectors at time k. The certain part of G is described by G 0 and the time invariant scalar θ G,i is the i th uncertain parameter associated with the matrix G i . Its value is unknown but its bounds are well-known, so, interval formalism is used to represent this notion [START_REF]Methods and Applications of Interval Analysis[END_REF], [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF]. A real interval is a closed, connected and bounded subset of r . If x is a real and scalar variable, then the associated interval is defined by a letter in italics x:

{ } [ ] x = ∈ ≤ ≤ = x x x x x x r / , ,
where x and x respectively denote the lower and upper bounds of x. To boil down to the case all interval variables have the same bounds -α, α and the same midpoint 0, the next modification is made:

( ) [ ] ( ) x x x = + - mid width α α α , 2 ,
where mid and width are the midpoint and the width of x:

( )

mid x x x = + 2 ,
( )

width x x x = -.
That is the reason why the case all intervals have the same bounds -α and α, remains general.

From (1), we can easily express the equality which links up the state vector X(k+j) (j∈n * ) to X(k) and inputs:

( ) ( ) ( ) ( ) ( ) ( ) X k j A X k A B U k i j A j i A B i j + = + + -- = - ∑ θ θ θ 1 0 1 . (2)
After multiplying (2) by C(θ C ) and piling measures on an observation window [k,k+s], the first static form is deduced:

( ) ( ) ( ) ( ) ( ) Y C H U k s X k k s s A C s A B C , , , , , = + - θ θ θ θ θ 1 , (3) ( ) ( ) ( )
[ ]

Y k, s = + Y k Y k s T T T L , ( ) Y , k s s t ∈r , ( )
s s m t = + 1 , ( ) ( ) ( ) [ ] U k, s -1 = + - U k U k s T T T L 1 , ( )
U , k s sq -∈ 1 r , ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) C s A C C T C A T C s A T T C C A C A θ θ θ θ θ θ θ , =       L , ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) H s A B C C B C s A B C B C B C A θ θ θ , , =               - 0 0 0 1 L M O M L , C s s n t ∈ × r , H s s sq t ∈ × r .
Generalised parity space methods are based on this static form (3) obtained for a certain model [START_REF] Chow | Analytical Redundancy and the Design of Robust Failure Detection System[END_REF], [START_REF] Lou | Optimally Robust Redundancy Relations for Failure Detection in Uncertain Systems[END_REF]. Then, external and internal forms of parity relations are defined and allow to respectively calculate residuals and their enclosures. Since the first only contains measured variables, uncertain parameters and terms which depend on the state vector, are put together in the second form.

To obtain a parity relation p(k+s), ( 3) is multiplied by a suitable row vector Ω T (whose choice will be later specified):

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) p k s k s k s X k k s T s s s T s A C residual terms T s A B C uncertain terms + = - - = = + - Ω Ω Ω Y H U H H C H U , , , , , , , , , , , 0 0 1
0 0 0 1 , θ θ θ θ θ θ 1 2 444 3 444 1 2 44444 3 44444 . (4)
Remarks. Since the external form is independent of uncertainties, it is valid whatever the uncertain matrices of the model. H s,0 and H s,θ respectively denote the certain and uncertain part of ( )

H s A B C θ θ θ , ,
. Under ideal circumstances (no noise and no fault), residuals are non-zero since they depend on uncertain parameters. Residual terms describe the coupling between residuals and X(k), which has to be minimised by the choice of Ω T .

2.b. Second static form

Instead of expressing X(k+j) according to the initial state vector X(k) as in (2), it is deduced from X(k+j-1):

( ) [ ] ( ) ( ) ( ) ( ) A I X k j X k j B U k j A B θ θ - + - +       = - + - 1 1 .
From this relation, state and input vectors are stacked on a time horizon [k,k+s] which leads to the next equation:

( ) ( ) ( ) ( ) M k s N k s A B θ θ X U , , = -1 , with: (5) ( ) ( ) ( )
[ ]

X k, s = + X k X k s T T T L , ( ) X , k s s n ∈r , ( )
s s n n = +1 , ( ) ( ) ( ) M A I A I A A n A n θ θ θ = - -           0 0 0 0 0 0 O O , M sn s n ∈ × r , ( ) ( ) ( )
N diag B B s B θ θ = - , N sn sq ∈ × r
, where diag s () is a bloc diagonal matrix and I n the n-order identity matrix.

Using the relationship between measures and state variables:

( ) ( ) ( ) ( ) ( ) ( ) Y X k s L k s diag C L C s C s s t n , , , = = ∈ + × θ θ θ , L C 1 r
, and the relation [START_REF] Isermann | Process fault detection based on modeling and estimation methods -a survey[END_REF], a new static form is finally deduced:

( ) ( ) ( ) ( ) ( ) ~, , ~, , C X H U Y s A C s B k s k s k s θ θ θ = -       1 , (6)
with:

( ) ( ) ( ) ~, C s A C C A L M θ θ θ θ =         , ( ) ( ) H s B st B I N θ θ =       0 0 , { } C s sn s s t n ∈ + × r , { } { } H s sn s sq s t t ∈ + × + r .
The expression ( 6) is general and valid whatever the uncertain matrices of the model. As the first static form, a parity relation ( ) Both equalities [START_REF] Massoumnia | Generating Parity Relations for Detecting and Identifying Control System Component Failures[END_REF] respectively define external and internal forms of the second static form. ΩT is chosen in order to minimise the influence of residual terms.

p
) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ~~~, , ~, , ~~, , , , , p k s k s k s k s k s k s T s s s T s A C residual t T s B uncertain + = -       = = - -       Ω Ω Ω H U Y H H C X H U Y 0 0 1 0

CHOICE OF THE PARITY MATRIX

3.a. Optimisation procedure

The aim of the method treated herein is the search of a row vector Ω T orthogonal to C s . Then, there are two possibilities: the matrix C s may be well-known or depend on parameter uncertainties. If the actuator matrix is uncertain, the classic generalised parity space method is used to determinate Ω T [2], [START_REF] Lou | Optimally Robust Redundancy Relations for Failure Detection in Uncertain Systems[END_REF], [START_REF] Massoumnia | Generating Parity Relations for Detecting and Identifying Control System Component Failures[END_REF]. Under present conditions, a total de-coupling with respect to state variables is made. Nevertheless, if the state or observation matrices are uncertain, the quantity Ω T c s (θ) has to be minimised with respect to the vector θ containing uncertain parameters θ i , i∈{1,...,r}. Therefore, the row vector Ω T is determined by minimising the criterion:

( ) ( ) ( ) J d d d d S T s r T s s T r = =         - - - - ∫ ∫ ∫ ∫ L L L L 1 2 44444 4 3 444444 α α α α α α α α θ θ θ θ θ θ θ Ω Ω Ω C C C 2 1 1 , (8)
under the constraint Ω T Ω = 1 (in order to eliminate the trivial solution). The orthogonality between Ω T and c s (θ) is measured by this criterion J, which describes the influence of the state vector X(k) on residuals. It is minimised by taking Ω equal to the eigenvector associated with the smallest eigenvalue of S. The value of the criterion J is well-known for a chosen eigenvalue λ and its eigenvector ω:

J T = = λω ω λ .
But, the analytical expression of S must be known to use the previous criterion. To do it, two cases must be studied according to the uncertain matrix.

3.b. Uncertain observation matrix C(θ θ c )

Let us assume the matrix C(θ c ) depends on r uncertain parameters θ i , then, it is written as:

( ) C C C c i i i r θ θ = + = ∑ 0 1 , (9)
where C 0 is the nominal (certain) part of C(θ c ).

The matrix C s (θ c ) defined in (3) is straightforwardly expressed from the equality (9): [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF] Using the properties:

( ) C C C s c s i s i i r θ θ = + = ∑ , , 0 1 , ( )
C s i i T i s T C C A , =       L T .
θ θ α α α i i d 2 3 2 3 = - ∫ and θ θ α α i i d = - ∫ 0 ,
the expression of the matrix S defined in (8) can be written as:

( ) ( ) S r s s T s s T s r s r T = + + +         2 3 0 0 2 1 1 α α C C C C C C , , , , , , L .
Since the analytical expression of S is known when the observation matrix C(θ c ) is uncertain, Ω T can be determined by the minimisation of the criterion J.

3.c. Uncertain state matrix A(θ θ A )

Some difficulties exist if the previous method is used when the state matrix A(θ A ) is uncertain:

( ) A A A A i i i r θ θ = + = ∑ 0 1 ( 11
)
where A 0 is the nominal part of the matrix A(θ A ).

Since the structure of the matrix C s (θ A ) (3) is not additive as [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF] and is highly non linear in respect to θ A , it is very difficult to determine the analytical expression of S. That is why the second static form ( 6) is used to solve this problem.

The structure of the matrix ( )

Cs A θ is deduced from (11):
( )

M M M A i i i r θ θ = + = ∑ 0
( ) ~~, , C C C s A s i i r s i θ θ = + = ∑ 0 1 , (12) with:
[ ]

~, C s T T L M 0 0 = T , [ ] ~, C s i i T T M = 0 .
So, the most orthogonal row vector ΩT to the uncertain matrix ( )

Cs A θ is obtained when the quantity ( ) ~Ω T C s A θ is
minimised with respect to the range of θ A . ΩT is obtained by minimising the following criterion:

( ) ( ) ( ) ~~~~~~J d d d d S T s A r T s A s A T r = =         - - - - ∫ ∫ ∫ ∫ L L L L 1 2 444444 4 3 4444444 α α α α α θ θ θ Ω Ω Ω C C C 2 1 1 , under the constraint ~Ω Ω T = 1 .
Since the analytical expressions of ( )

Cs
= + + +         2 3 0 0 2 1 1 α α C C C C C C L .
Ω is logically given by the eigenvector associated with the smallest eigenvalue of the matrix S .

INTERNAL FORMS OF PARITY RELATIONS

In this section, internal forms (4) or [START_REF] Massoumnia | Generating Parity Relations for Detecting and Identifying Control System Component Failures[END_REF] and the procedure allowing to calculate the parity matrix are expressed with respect to the uncertain matrix of the state representation.

4.a. Uncertain actuator matrix B(θ θ B )

If the actuator matrix is uncertain:

( ) B B B B i i i r θ θ = + = ∑ 0 1
, where B 0 is the nominal part of B(θ B ), then the matrix H s (θ B ) is defined by the additive expression:

( )

H H H s B s i s i i r θ θ = + = ∑ , , 0 1 , with: H s
i i s i i CB CA B CB , =             - 0 0 0 1 L L M O M L , i∈{0 ,1,...,r}.
The parity relation (4) becomes:

( ) ( ) ( ) ( ) ( ) p k s k s k s X k k s T s T s i s i i r ( ) , , , , , + = -
- = + -       = ∑ Ω Ω Y H U C H U 0 1 1 1 θ . ( 13
)
Since the matrix C s is certain, a total de-coupling with respect to state variables can be obtained. Generalised parity space method is applied to determine a vector Ω T orthogonal to the matrix c s : Ω T c s = 0. Therefore, residual terms are equal to zero and the previous internal form (13) becomes:

( ) ( ) p k s k s i i r T s i ( ) , , + = - = ∑ θ 1 1 Ω H U .

4.b. Uncertain observation matrix C(θ θ C )

From ( 9), the additive structure of H s (θ C ) is deduced:

( )

H H H s c s i s i i r θ θ = + = ∑ , , 0 1
, with:

H s i i i s i C B C A B C B , =             - 0 0 0 1 L L M O M L , i∈{0,1,...,r}.
Using (4) and the expression of C s (θ C ) [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF], the following parity relation is obtained:

( ) ( ) ( ) ( ) ( ) p k s k s k s X k k s T s T s i s i i r i s i i r ( ) , , , , , , , +
= - - = +       + -           = = ∑ ∑ Ω Ω Y H U C C H U 0 0 1 1 1 1 θ θ .
The calculus of Ω T using the criterion J ( 8) is interesting because if the chosen eigenvalue is equal to zero, then the residual is de-coupled from unknown state variables.

The term Ω T C s (θ C ) is equal to zero for all values of the r uncertain parameters θ i if the following properties hold: Ω T C S,0 = 0 and Ω T C S,i = 0, i ∈ {1,...,r}, θ i ∈ [-α,α].

(14) If the criterion J is equal to zero (as the associated eigenvalue), then the following equality holds:

( ) Ω Ω Ω T r s T s T s r T S =         = 2 3 3 0 0 1 2 α α α C C C , , , L . ( 15
)
If ( 14) is verified, so does (15), and the reverse is always true. Then, a minimal observation window size always exists for which the matrix C S (θ C ) defined in ( 10) is full row rank, which allows a total de-coupling. To show that, the Cayley-Hamilton theorem is applied on the state matrix A. Thus, it is proved there is at least one non-zero row vector Ω T which satisfies the following property:

( ) { } Ω Ω T j T j n T T T n j C C A L L       = = ∀ ∈ C , , 0 j 0, , r .
For an observation window size n, the product on the left hand size of C n (θ C ) by Ω T is equal to zero and so does Ω T S. These relations hold for an integer greater than n. In conclusion, several row vectors Ω T orthogonal to C S (θ C ) are obtained for an observation window size s greater or equal to n. The vectors Ω are defined by the eigenvectors associated with the eigenvalues of S equal to zero. The observation window size s has to be sufficient to obtain rank(C) eigenvalues equal to zero or smaller than a fixed tolerance.

4.c. Uncertain state matrix A(θ θ A )

Only the second form static is interesting because the determination of Ω T is otherwise very problematic. From [START_REF] Massoumnia | Generating Parity Relations for Detecting and Identifying Control System Component Failures[END_REF], the expression of the parity relation becomes:

( ) ( ) ( ) ( ) ( ) ~~~, , ~~, p k s k s k s k s T s T s A + = -       = Ω Ω H U Y C X 1 θ ,
where the matrix ( ) [START_REF] Ragot | A propos de l'algèbre des intervalles, Application à la validation de données[END_REF]. Since the matrix Hs is certain, uncertain terms are thus equal to zero. Like the previous paragraph, the eigenvectors associated with the eigenvalues of S equal to zero allow a total de-coupling.

Cs A θ is defined in
Nevertheless, this situation is not always possible. Then, the internal form becomes:

( ) ( ) ( ) ~~~, ~~, , , p k s k s k s T s i i r T s i + = + = ∑ Ω Ω C X C X 0 1
θ .

Since it depends on the state vector expressed on an observation window [k,k+s], the evaluation of its enclosure is problematic when the observation matrix C is not full column rank. Otherwise, the state variables are estimated:

( ) ( ) ( ) $ , , X Y k s diag C C C k s s T T =       + - 1 1
.

ENCLOSURES EVALUATION

5.a. Interval arithmetic

Uncertain parameters are defined by bounded variables intervening in the previous internal forms. To evaluate their bounds at any time, interval formalism is used [START_REF]Methods and Applications of Interval Analysis[END_REF], [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF], [START_REF] Ragot | A propos de l'algèbre des intervalles, Application à la validation de données[END_REF]. This procedure determines enclosures which represent the normal operation field. A fault is detected if a residual calculated from its external form, goes out its enclosure.

The arithmetic operations (+, -, ×, /) on real variables can be reformulated in the case of independent interval variables.

Operation

Interval obtained

Addition

[ ]

x y , + = + + x y x y Subtraction [ ] x y , -= - - x y x y Product ( ) ( ) [ ]

Table 1. Interval arithmetic operations

An elementary function g is a real function, continuous on every closed interval on which it is defined (absolute value, square, square root, exponential, logarithm, sine, cosine, arc tangent,...). Then, it can be reformulated in the case of interval variables: ( ) ( )

{ } g x x = ∈ g x x / .
An interval extension f of a real function f is an interval function with the following property: f(x) = f(x), x ∈ x. If this extension is inclusion monotonic [START_REF]Methods and Applications of Interval Analysis[END_REF], it satisfies: if x ⊂ y then f(x) ⊂ f(y).

For any given real function f, there are infinitely many extensions. One of them, called united extension, satisfies:

( ) ( ) f x x = ∈ U x f x
The global extrema of an elementary function g are often known on a range x, then the bounds of the interval function ( ) g x (and the united extension of g), are defined by: ( ) ( ) ( )

g x x x = ∈ ∈       inf x g x x g x , sup
.

But, for more complicated functions, this would require solving two global optimisation problems, which often exceeds time consuming. Then, if f results from the composition of real operators or elementary functions, it is impossible to compute the united extension. So, natural extensions f are used. They are obtained by replacing, in the expression of f, its real argument x by x, and the elementary functions or operators by the associated united extensions. An inclusion monotonic and natural extension always includes the united extension but they are different in general.

5.b. Dependence between interval variables

In the case of dependence, where some interval variables occur repeatedly in an expression, interval arithmetic leads to an overestimation of the minimal range associated with an united extension. To show that, the following interval extension is considered: f(x) = xx. From definition, the united extension is given by: ( ) { }

~/ x f x = - ∈ x x x
. So, for every value of the real variable x, this expression is equal to {0}. Nevertheless, when the formula giving the subtraction is used, ( ) [ ]

f x = - - x x x x , is deduced. If x is not reduced to a
real number, f leads to an overestimation which depends on the width of x. In fact, the dependence is not taken into account, because the following operation is made: xz, where z is an independent interval variable with the same bounds as x. This problem leads to the subdistributivity property:

x(y ± z) ⊆ xy ± xz or (z ± y) x ⊆ zx ± yx.
So, the natural extensions of two equivalent functions in the arithmetic of real numbers are not necessary equivalent in the interval arithmetic. Nevertheless, if each interval variable appears at most once in the expression of an natural extension, then its evaluation by interval arithmetic leads to the united extension. Otherwise, for determining it without searching the extrema of the real function by non linear computing (Newton-Raphson, simplex,...), some recursive algorithms based on Mean Value or Centered Forms as the Monoticity Test Form can be used [START_REF]Methods and Applications of Interval Analysis[END_REF], [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF].

EXAMPLE

Let us consider the following system where the state matrix depends on an uncertain parameter θ 1 ∈[-0.5,0.5]:

A A B C 0 1 0 8 1 0 0 2 0 1 0 0 0 0 1 1 0 0 1 =       =       =       =       . . , . , , .
The study of C 1 (θ 1 ) (θ 1 is time invariant):

( ) ( )

C 1 1 1 1 1 0 0 1 0 8 0 1 1 0 0 2 θ θ θ =       = +             C CA . . .
, shows this matrix is always of rank 2 for every possible value of θ 1 . So, the observation window size is chosen equal to 1. Since ( ) C1 1 θ is uncertain, the matrix S is built. It is not full rank and one of its eigenvalues is equal to zero. The others are equal to: 2.30.10 -4 , 1, 1, 2.02 and 3.67. In case the chosen eigenvalue is equal to zero, the enclosure of the residual is reduced to 0. In the other cases, interval arithmetic is used to evaluate the internal form, which is very simple in that example:

( ) ~, p k a b k k θ θ 1 1 = + ,
where the parameters a k and b k depend on measures. Since there is no dependency between interval variables in a same residual, its bounds are easily deduced: The model is simulated during one second. The input is piece wise constant with a magnitude comprised between -1 and 1. Its value changes every 0.2s. The measures y 1 and y 2 obtained when θ 1 is a constant and is successively equal to -0.5 and 0.5 are represented below. The residual and its enclosure obtained for λ = 2.30.10 -4 are presented below. Both situations where θ 1 is equal to -α or α are examined. At time t = 0.5s, a bias equal to 0.625 occurs on the first sensor. It represents 10 per cent of the maximal magnitude of the measure y 1 obtained for an input equal to 1. During a normal operation, the residual stays inside its enclosure (t < 5s). But, when the fault is present, the detection depends on the operating point and the value of the uncertain parameter.

CONCLUSION

This paper deals with a fault detection method suited to structured uncertain models. According to the uncertain matrix of the state representation, different techniques are proposed to choose the parity matrix. In the case of an uncertain observation matrix, the choice of eigenvalues equal to zero during the minimisation of the criterion J allows to solve the problem of state estimation. Otherwise, the second static form allows to treat the case of an uncertain state matrix. Once the external forms of parity relations were expressed, the internal forms are evaluated by means of interval arithmetic. Then the enclosures, which define the normal operation field, are obtained. If a residual leaves its enclosure, a fault is detected.

(

Figure 1 .

 1 Figure 1. Input and sensory observations

Figure 2 .

 2 Figure 2. Bias on the first sensor