# Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach

Fabien Gatti, Christophe Iung

► To cite this version:

Fabien Gatti, Christophe Iung. Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach. Physics Reports, 2009, 484, pp.10.1016/j.physrep.2009.05.003. hal-00412650

HAL Id: hal-00412650

https://hal.science/hal-00412650v1

Submitted on 18 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Exact and constrained kinetic energy operators for polyatomic
molecules: The polyspherical approach

Fabien Gatti and Christophe Iung*

CTMM, Institut Charles Gerhardt, UMR 5253,

CC 1501, Université Montpellier II,

F - 34095 Montpellier, Cedex 05, France

(Dated: May 12, 2009)

Abstract

A review of the *polyspherical approach* to the kinetic operators for polyatomic molecules is given.
This approach provides general and correct forms of the kinetic energy operator (KEO) expressed
in terms of curvilinear coordinates. These forms are well adapted to the physical description of
molecular systems and to the numerical methods used to solve the Schrödinger equation. The
approach derives its name, *polyspherical*, from the fact that the operators are expressed in terms
of spherical coordinates eventually. These kinetic energy operators can be exploited to treat a
large variety of problems such as the calculation of infrared or photo-absorption spectra or the
study of reactive scattering systems. Special emphasis is placed on concrete examples.

**Keywords:** quantum molecular dynamics, curvilinear coordinates, polyspherical coordinates, kinetic energy operator.

<footnote>*E-mail:gatti@univ-montp2.fr,iung@univ-montp2.fr</footnote>
# Contents

## I. Introduction 4

## II. The polyspherical approach: Theory 8

### A. Theoretical background 8

### B. Vector parametrization and properties of angular momenta 15

#### 1. Examples 16

#### 2. General formulation 20

### C. A general expression of the KEO in standard polyspherical coordinates 25

#### 1. Expression in angular momenta 26

#### 2. Matrix expression of the KEO 41

### D. Separation into subsystems 46

### E. Checking the correctness of the operators 56

## III. Explicit expressions for particular cases and applications 58

### A. Flexible systems 59

#### 1. Tri- and tetra-atomic flexible systems 59

#### 2. Weakly bound molecular trimers 62

#### 3. Protonated methane in Radau coordinates 63

### B. Semi-rigid systems 65

#### 1. Tetra-atomic molecules 65

#### 2. Penta-atomic molecules 70

### C. Separation into subsystems 73

#### 1. Protonated water dimer (Zundel cation) 73

#### 2. Protonated water tetramer (Eigen cation) 81

## IV. Description of large systems with active polyspherical coordinates 84

### A. Kinetic energy operators for rigidly and adiabatically constrained molecular systems 86

#### 1. Rigid-constraint approach 88

#### 2. Adiabatic-constraint approach 92

### B. Combination with the polyspherical approach 95
<table>
  <tbody>
    <tr>
      <td>C. Examples of applications</td>
      <td>97</td>
    </tr>
    <tr>
      <td>1. Rigid- and Adiabatic-constraint approaches: the water dimer</td>
      <td>97</td>
    </tr>
    <tr>
      <td>2. The photo-induced non-adiabatic dynamics of ethene</td>
      <td>104</td>
    </tr>
    <tr>
      <td>3. General strategies for very large systems</td>
      <td>107</td>
    </tr>
    <tr>
      <td>V. Conclusion and outlook</td>
      <td>118</td>
    </tr>
    <tr>
      <td>Acknowledgements</td>
      <td>119</td>
    </tr>
    <tr>
      <td>A. Properties of the partial angular momenta</td>
      <td>120</td>
    </tr>
    <tr>
      <td>B. General expression of the KEO in angular momenta for non-orthogonal coordinates</td>
      <td>125</td>
    </tr>
    <tr>
      <td>C. General expression of the KEO in conjugate momenta</td>
      <td>129</td>
    </tr>
    <tr>
      <td>D. KEO for HCF<sub>3</sub> in valence coordinates</td>
      <td>134</td>
    </tr>
    <tr>
      <td>References</td>
      <td>138</td>
    </tr>
  </tbody>
</table>
# I. INTRODUCTION

During the past years, considerable progress has been made in the field of molecular quantum dynamics. The domain of applications encompasses broad areas in physics or chemistry: photo-dissociations or excitations, dissociations of an adsorbate on a metal surface, intramolecular vibrational energy redistribution and predissociation, infrared spectroscopy, inelastic surface scattering, molecular reactive scattering, evolution of a molecular system excited by an ultra-short laser pulse. On the experimental front, decisive progress, in particular in the field of femto-spectroscopy [1-3] or infra-red/microwave spectroscopy [4-7], allows the scientists to probe chemical phenomena at an atomic time scale or to obtain fully resolved spectra of highly excited polyatomic systems. These ultra-sophisticated experiments require new theoretical tools to interpret, predict and, in other words, accompany these experimental works.

In addition, it should be emphasized that these molecular processes are generally impacted to a significant extent by nuclear quantum mechanical effects [8-10] such as zero-point energy effects or tunnelling of light atoms through barriers (for instance in the case of electron or proton transfers), transitions due to strong vibronic couplings such as conical intersections, which seem to play a crucial role in many organic or biological systems. Several algorithms for solving the Schrödinger equation have thus been developed which have wide applicability and allow to treat larger systems than in the past. In particular, one can cite the Multi-Configuration Time Dependent Hartree (MCTDH) [11-20] approach, the Multimode code [21-25] based on a time independent Vibrational Self-Consistent Field approach coupled to a Configuration Interaction (CI-VSCF) procedure (see also [26-32] for other CI-VSCF methods), or the Wave Operator Sorting Algorithm (WOSA) based on the extraction of an active space within the framework of the Bloch formalism, [33-36].

When solving the Schrödinger equation, the choice of the set of coordinates used to describe the physical systems is of great importance. First, for a given system treated in full dimensionality, the convergence of the algorithms that solve the Schrödinger equation depends on the correlation between the coordinates. More precisely, a set of coordinates that minimizes the correlation will improve convergence, while a coordinate set that
introduces strong artificial correlations, i. e. correlations that are entirely due to an
unsuitable choice of coordinates, will slow down convergence. Let us consider an example.
For small amplitude motions around a well defined equilibrium geometry, the vibrations are
often rather harmonic. The well-known normal-mode rectilinear coordinates then make the
Hamiltonian operator almost separable – it is exactly separable for infinitesimal distortions
– and the use of these coordinates will be optimal in a low-energy domain. However, for
more floppy systems exhibiting two or several minima, or at higher excitation energies,
the vibrational amplitudes become larger and the rectilinear normal coordinates cease to
describe the motion in a natural way. This introduces strong artificial correlations. In such
situations, in particular when studying scattering, dissociation, or isomerization problems,
the use of appropriate internal coordinates becomes crucial. In general, curvilinear coordi-
nates, involving angles, are the natural choice, as they usually lead to a more separable and
hence less artificially correlated Hamiltonian operator.

Second, for very large systems, the correspondingly large number of internal coordinates
makes approximations unavoidable. This is true not only for the dynamics, but also,
and above all, for the calculation of the potential operator, i.e. the generation of global
potential energy surfaces, which is restricted to systems with only a few degrees of
freedom. Fortunately, chemical processes can often be described by a limited number of
active degrees of freedom. Under these conditions it is reasonable to invoke simplifying
approximations for the other degrees of freedom. This permits the construction of simple
model Hamiltonians [37–41]. Several general models have been formulated, such as the
reaction surface Hamiltonian [42–45], which makes use of a harmonic approximation locally
defined along one or several reaction coordinates, and the multistate-multimode vibronic
coupling model [46–50] to calculate photoabsorption spectra involving conical intersections.
The choice of the pertinent coordinates is crucial and determined by physical intuition.
It may happen that curvilinear coordinates are the only coordinates that can correctly
describe these chemical processes with a limited number of degrees of freedom. For instance,
several works have proved (i) that the reaction path Hamiltonian may sometimes be far
better described in terms of curvilinear coordinates than of rectilinear coordinates [51, 52];
(ii) that the crossing seam of a conical intersection is, in general, curved [53, 54] and, thus,
that curvilinear coordinates are essential to describe the behavior of the extended seam
[55, 56].

Unfortunately, the use of curvilinear coordinates often leads to very complicated expressions of the kinetic energy operator (KEO) [57, 58], which need to be derived for a particular system [59–66] and are not easy to generalize. This is in contrast to the rectilinear coordinates, which simplify the mathematical formulation of the same operator in a systematic way. The problem is not primarily to derive a formula for the KEO. An algorithmic program such as Mathematica [67] can be used to calculate the operators analytically. A numerical computation of the action of the kinetic operator is also possible and several contributions have been put forth in this direction (see, e. g., Refs. [68–71]). Consequently, the crucial point is to find a general form of the KEO that is (i) as compact as possible and (ii) well adapted to the numerical methods used in dynamics.

The *polyspherical approach*, presented in a series of papers [72–83], accomplishes both of these requirements. It is a general formulation of the exact KEO of an *N*-atom system whose configuration, in the center-of-mass frame, is described by *N*−1 relative position vectors. The approach derives its name *polyspherical* from the fact that the operators are expressed in terms of spherical coordinates that parametrize the *N*−1 relative position vectors eventually.

The polyspherical approach is characterized by the five following properties: (i) it explicitly provides rather compact and general expressions of the exact kinetic energy operator including rotation and Coriolis coupling and avoids the use of differential calculus when deriving these operators. (ii) Within this approach it is very easy to find a primitive basis set (e. g. a basis set of spherical harmonics) that discards all singularities that may occur in the KEO. (iii) General expressions for the KEO are explicitly provided. (iv) There is much freedom in choosing the underlying set of vectors: they may be of Jacobi, Radau, Valence, satellite type, or a combination of these. (v) When polyspherical coordinates are used, the KEO always is separable: i.e. it can be written as a sum of products of monomodal operators.

The purpose of the review is to provide an overview of the polyspherical approach. In Part II, the polyspherical method is introduced. The background theory is dealt with briefly,
before the polyspherical approach is looked at in detail. In Part III, concrete applications
on systems from 4 to 13 atoms are presented: the exact KEOs of these systems involving
all the degrees of freedom derived from the general expressions are given explicitly. Part IV
is devoted to the study of large systems: starting from a limited set of active curvilinear
coordinates, we present a general strategy to obtain the correct KEOs. Several examples
are given for systems such as organometallic complexes or organic molecules involved in
photochemical reactivity.
## II. THE POLYSPHERICAL APPROACH: THEORY

### A. Theoretical background

(i) Kinetic energy operators in curvilinear coordinates:

Systems of $N$ particles possess $3N$ degrees of freedom, of which three can be eliminated by translational invariance when there is no external field. The reduced system thus obtained has $3N-3$ degrees of freedom, and its configuration can be specified by $N-1$ relative position vectors. Hereafter, $G$ will denote the center of mass of the molecular system and the Space-Fixed (SF) frame will denote the center of mass frame whose origin coincides with $G$ and whose axes coincide with the axes of the Laboratory-Fixed (LF) frame. The convention $\hbar = 1$ is used everywhere in the present review. Let the configuration of an $N$-particle system be described by $3N-3$ generalized coordinates $q_i$ ($i = 1, \dots, 3N-3$), and

$$x_{\alpha} = x_{\alpha}(q_1, \dots, q_i, \dots, q_{3N-3}) \quad (\alpha = 1, \dots, 3N-3), \qquad (1)$$

be the expression of the SF mass-weighted Cartesian coordinates associated with the $N-1$ vectors describing the system in terms of the chosen generalized coordinates. (The mass-weighted coordinates are obtained by multiplying the standard Cartesian coordinates by $\sqrt{m}$, $m$ being the mass of the particle). Coordinates that may be expressed as linear combinations of Cartesian coordinates are called rectilinear, and those that are nonlinear functions of Cartesian coordinates are called curvilinear. As shown by Podolsky, the expression of the kinetic energy operator in curvilinear coordinates $q_i$ ($i = 1, \dots, 3N-3$) is given by the expression $(-\frac{1}{2})$ times the Laplacian) [57, 58, 84]:

$$\hat{T} = -\frac{1}{2} \sum_{i,j}^{3N-3} J_e^{-1} \frac{\partial}{\partial q_i} J_e G_{ij} \frac{\partial}{\partial q_j}, \qquad (2)$$

with

$$J_e = |Det(\frac{\partial x_{\alpha}}{\partial q_i})|, \qquad (3)$$

the absolute value of the Jacobian (in the present review, 'the Jacobian' stands for the Jacobian determinant, i.e. the determinant of the Jacobian matrix, and not for the matrix itself. Det denotes the determinant of a matrix) corresponding to the coordinate transformation
from $q_i$ to $x_\alpha$. We assume that $J_e > 0$ except, possibly, on a null set. $G_{ij}$ ($i, j = 1, \dots, 3N-3$) in Eq. (2) is given by

$$G_{ij} = \sum_{\alpha=1}^{3N-3} \frac{\partial q_i}{\partial x_\alpha} \frac{\partial q_j}{\partial x_\alpha} \quad (i, j = 1, \dots, 3N-3). \qquad (4)$$

In addition, the inverse of $G$,

$$g = [g_{ij}] = G^{-1} = [G_{ij}]^{-1}, \qquad (5)$$

satisfies

$$g_{ij} = \sum_{\alpha=1}^{3N-3} \frac{\partial x_{\alpha}}{\partial q_{i}} \frac{\partial x_{\alpha}}{\partial q_{j}} \quad (i, j = 1, \dots, 3N-3). \qquad (6)$$

Eqs. (5,6) come from the fact that $G_{ij}$ and $g_{ij}$ are the contravariant and covariant components, respectively, of the metric tensor used to measure distances and angles in terms of curvilinear coordinates $q_i$ in the Euclidean space spanned by the mass-weighted Cartesian coordinates $x_\alpha$. Note, however, that tensor notation is avoided in the present review in order to simplify the equations. In other words, the matrix elements $G_{ij}$ and $g_{ij}$ correspond to $g^{ij}$ and $g_{ij}$, respectively, in tensor notation. It is then straightforward to rewrite equation (2) as

$$2\hat{T} = \sum_{i,j=1}^{3N-3} \hat{P}_{q_i}^{\dagger} G_{ij} \hat{P}_{q_j}, \qquad (7)$$

with the momentum operators:

$$\hat{P}_{q_i} = \frac{1}{i} \frac{\partial}{\partial q_i}, \qquad (8)$$

and their adjoints [58]

$$\hat{P}_{q_i}^{\dagger} = \frac{1}{i} J_e^{-1} \frac{\partial}{\partial q_i} J_e, \qquad (9)$$

which can be rewritten as

$$\begin{aligned} \hat{P}_{q_i}^{\dagger} &= \hat{P}_{q_i} + \Lambda_i, \\ \Lambda_i &= (J_e^{-1} \hat{P}_{q_i} J_e), \end{aligned} \qquad (10)$$

with $\Lambda_i$, a purely multiplicative operator. The parentheses in Eq. (10) indicate that $\hat{P}_{q_i}$ does not operate beyond the parentheses. In other words,

$$(J_e^{-1} \hat{P}_{q_i} J_e) = -i J_e^{-1} \frac{\partial(J_e)}{\partial q_i}. \qquad (11)$$
Note also that

$$J_e = (|Det(G_{ij})|)^{-\frac{1}{2}} = g^{-\frac{1}{2}}. \qquad (12)$$

(ii) *Change in the convention of normalization:*

The kinetic energy operator of Eqs. (2,7) is correct in case the wave functions are normalized as follows:

$$\int \Psi^* \Psi d\tau = 1, \qquad (13)$$

where

$$d\tau = J_e dq_1 \dots dq_{3N-3} \qquad (14)$$

is the usual Euclidean volume element.

Now, if another volume element is used in normalizing the wave-functions, i.e.

$$d\tau' = \rho dq_1 \dots dq_{3N-3}, \qquad (15)$$

where $\rho$ is an arbitrary weight function (we only assume that $\rho \ge 0$ and that, if $\rho = 0$, it is only on a null set), the wave function $\Psi$ and the kinetic energy operator are to be replaced respectively by [58]

$$\Psi' = J_e^{\frac{1}{2}} \rho^{-\frac{1}{2}} \Psi, \qquad (16)$$

and

$$\hat{T}' = J_e^{\frac{1}{2}} \rho^{-\frac{1}{2}} \hat{T} \rho^{\frac{1}{2}} J_e^{-\frac{1}{2}}, \qquad (17)$$

so as to preserve the normalization of the wave-functions. The adjoint of the conjugate momentum $\hat{P}_{q_i}$ becomes
$$ \begin{aligned} \hat{P}_{q_i}^{\dagger'} &= \hat{P}_{q_i} + \Lambda_i', \\ \Lambda_i' &= (\rho^{-1} \hat{P}_{q_i} \rho). \end{aligned} \qquad (18) $$

Let us consider an important particular case: the $\rho = 1$ case. If $\rho = 1$, Eq. (17) expands into [58, 85–87]

$$ 2\hat{T}' = \sum_{i,j}^{3N-3} \hat{P}_{q_i} \mathcal{G}_{ij} \hat{P}_{q_j} + 2V', \qquad (19) $$

with $V'_D$ a purely multiplicative operator called 'extra potential term', which reads

$$ 2V' = \left( \sum_{i,j}^{3N-3} J_e^{-\frac{1}{2}} \hat{P}_{q_i} J_e \mathcal{G}_{ij} \hat{P}_{q_j} J_e^{-\frac{1}{2}} \right). \qquad (20) $$

(Again, the parentheses in Eq. (20) indicate that $\hat{P}_{q_i}$ and $\hat{P}_{q_j}$ do not operate beyond those).

### (iii) Introduction of quasi-momenta:

Now, let us introduce a new set of linear combinations of the $\hat{P}_{q_i}$ ($n \ge 3N - 3$):

$$ \hat{P}_K = \sum_{i=1}^{3N-3} A_K^i \hat{P}_{q_i} (K = 1, \dots, n). \qquad (21) $$

The following cases may occur:

(1) there exists a coordinate transformation from $q_i$ to some new coordinates $Q_K$ such that $A_K^i$ is the ordinary Jacobian matrix:

$$ A_K^i = \frac{\partial q_i}{\partial Q_K}, \qquad (22) $$

so that $\hat{P}_K$ is the momentum conjugate to the curvilinear coordinate $Q_K$.

(2) no such coordinate transformation exists and the $\hat{P}_K$ that cannot be written as momentum conjugates are usually said to be quasi-momenta [58, 81]. The usual angular momentum component operators introduced in quantum mechanics are a particular case of such operators.

The KEO can be recast in terms of these new momenta:

$$ 2\hat{T} = \sum_{K,L=1}^{n} \hat{P}_{K}^{\dagger} \mathcal{G}'_{KL} \hat{P}_{L}, \qquad (23) $$
with

$$ \hat{P}_K^\dagger = \sum_{i=1}^{3N-3} \hat{P}_{q_i}^\dagger A_K^i, \qquad (24) $$

the adjoint of $\hat{P}_K$. The elements of matrix $\mathcal{G}'$ are functions of the elements of the matrices $\mathcal{G}$ and $\mathbf{A}$:

$$ \mathcal{G} = \mathbf{A}^T \mathcal{G}' \mathbf{A}, \qquad (25) $$

$\mathbf{A}^T$ denotes the transpose of $\mathbf{A}$.

Let us consider an important particular case: the three BF-components of $\hat{\mathcal{J}}$, the total angular momentum of a molecular system. In order to separate the overall rotation from the internal deformation, let us introduce a Body-Fixed (BF) frame, $\{G; \vec{e}_{xBF}, \vec{e}_{yBF}, \vec{e}_{zBF}\}$. The BF frame is linked to the system, and its axes $\vec{e}_{xBF}$, $\vec{e}_{yBF}$ and $\vec{e}_{zBF}$ rotate in a conventional manner when the particles move. The orientation of the BF frame with respect to the SF frame is determined by three Euler angles: $\alpha, \beta, \gamma$, and, after definition of the three Euler angles, the shape and the size of the molecules are described by $3N-6$ internal BF coordinates. We have $3N-3$ coordinates such as $\{q_i, i = 1, \dots, 3N-6\}$ correspond to the $3N-6$ internal BF coordinates and $q_{3N-5} = \alpha$, $q_{3N-4} = \beta$, and $q_{3N-3} = \gamma$. This separation greatly facilitates the construction of the irreducible representations of the rotation symmetry group, thus reducing the size of the calculations needed to solve the Schrödinger equation. The three BF-components of $\hat{\mathcal{J}}$ can be expressed in terms of the three momentum operators associated with the Euler angles as

$$ \begin{bmatrix} \hat{J}_{xBF} \\ \hat{J}_{yBF} \\ \hat{J}_{zBF} \end{bmatrix} = \mathbf{A}_J \begin{bmatrix} \frac{1}{i} \frac{\partial}{\partial \alpha} \\ \\ \frac{1}{i} \frac{\partial}{\partial \beta} \\ \\ \frac{1}{i} \frac{\partial}{\partial \gamma} \end{bmatrix}. \qquad (26) $$

(see Eq. (77), in Section II C 1 below, for an explicit expression of $\mathbf{A}_J$). The KEO can be rewritten in terms of the $3N-6$ momentum operators associated with the $3N-6$ internal
BF coordinates and of the three quasi-momenta, $\hat{J}_{xBF}$, $\hat{J}_{yBF}$, and $\hat{J}_{zBF}$:

$$ \hat{T} = \sum_{l,m=1}^{3N-6} \frac{\hat{P}_{ql}^{\dagger} \sum_{q_l q_m} \hat{P}_{q_m}}{2} + \sum_{l=1}^{3N-6} \sum_{\alpha=x,y,z} \frac{\hat{P}_{ql}^{\dagger} \sigma_{q_l \alpha} \hat{J}_{\alpha BF} + \hat{J}_{\alpha BF} \sigma_{\alpha q_l} \hat{P}_{ql}}{2} \\ + \sum_{\alpha=x,y,z} \sum_{\beta=x,y,z} \frac{\hat{J}_{\alpha BF} \Gamma_{\alpha \beta} \hat{J}_{\beta BF} + \hat{J}_{\beta BF} \Gamma_{\beta \alpha} \hat{J}_{\alpha BF}}{2}. \qquad (27) $$

The matrices $\Sigma$, $\sigma$ and $\Gamma$ parametrize the vibrational ($\hat{T}_{vib}$), Coriolis ($\hat{T}_{cor}$) and rotational ($\hat{T}_{rot}$) parts of the KEO, respectively. (To be rigorous, it is never possible to fully isolate the different parts of the KEO corresponding to the vibration, the rotation and the Coriolis coupling. This depends on the coordinate system and the definition of the BF frame. Hence, there is also some Coriolis coupling in $\hat{T}_{vib}$ and $\hat{T}_{rot}$. However, we keep this notation for simplicity). The elements of the matrices $\sigma$ and $\Gamma$ depend on the elements of the matrices $\mathcal{G}$ and $A_J$ of Eqs. (7) and (26), respectively. Indeed, we have

$$ \mathcal{G} = \begin{bmatrix} 1 & 0 \\ 0 & A_J^T \end{bmatrix} \begin{bmatrix} \Sigma & \sigma^T \\ \sigma & \Gamma \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & A_J \end{bmatrix} = \begin{bmatrix} \Sigma & \sigma^T A_J \\ A_J^T \sigma & A_J^T \Gamma A_J \end{bmatrix}. \qquad (28) $$

The elements of the matrix $\Sigma$ are identical to the elements of the matrix $\mathcal{G}$ corresponding to the $3N-6$ internal coordinates. In other word, $\Sigma_{q_l q_j} = \mathcal{G}_{ij}$, $(i, j = 1, \dots, 3N - 6)$. The matrix $\Sigma$ does not depend on the definition of the BF frame.

Let us point out that it is possible to introduce other quasi-momenta, for instance, the BF-components of the angular momenta associated with the $N-1$ relative position vectors that will appear, in the polyspherical approach, in Section II B.

Another important 'historical' case of quasi-momenta are the so-called vibrational momentum components, $\hat{\pi}_{\alpha}$ ($\alpha = x, y, z$), introduced by Watson [88, 89]. They read

$$ \hat{\pi}_{\alpha} = -i \sum_{k,l=1}^{3N-6} \zeta_{kl}^{\alpha} Q_k \frac{\partial}{\partial Q_l}, \qquad (29) $$

where $Q_k$ are the normal coordinates and $\zeta_{km}^{\alpha}$ are Coriolis parameters defined in Ref. [88]. It goes without saying that they are different from the angular momenta introduced in the polyspherical approach.
Finally, let us point out that it is also possible to introduce hermitian quasi-momentum operators by the standard procedure:

$$ \hat{P}_K^H = \hat{P}_K + \hat{P}_K^\dagger. \qquad (30) $$

It can be shown (see Chapter 3 in Ref. [58]) that

$$ \hat{P}_K^H = J^{-\frac{1}{2}} \hat{P}_K J^{\frac{1}{2}} + \frac{1}{2} \lambda_K, \qquad (31) $$

with

$$ \lambda_K = \sum_{i=1}^{3N-3} (\hat{P}_{q_i} A_K^i), \qquad (32) $$

and that

$$ 2\hat{T} = \sum_{K,L=1}^{n} g^{\frac{1}{4}}(\hat{P}_{K}^{H} + \frac{1}{2}\lambda_{K}) \mathcal{G}_{KL}^{\prime} g^{-\frac{1}{2}}(\hat{P}_{L}^{H} - \frac{1}{2}\lambda_{L})g^{\frac{1}{4}}. \qquad (33) $$

The latter expression of the KEO is similar to the general expression given by Watson in Ref. [88, 89]: see, for instance, Eq. (A.23) in Ref. [89]. As explained by Watson the terms in $\lambda_K$ or $\lambda_L$ are closely connected with the non-commutativity of the quasi-momenta.

### (iv) Derivation of kinetic energy operator:

The derivation of the KEO of Eqs. (7,23,27) is a difficult task. In concrete terms, it is possible to follow one of two routes: either to obtain the classical kinetic energy expression and to quantize it [59, 73, 74, 87] or to use the chain rule to transform the space-fixed Cartesian KEO to derive the new KEO [60, 61, 63–65, 90–93]. Note that the polyspherical approach uses the first route and this is why we shall start from the classical expression of the kinetic energy in Sections IIB 1 and IIB 2 and quantize it in a second step. In this context, Lukka greatly simplified the derivation of KEOs in bond coordinates by suggesting the use of infinitesimal rotation angles [94]. In the same manner, Pesonen proposed to use geometric algebra to easily calculate the KEOs [95–97]. However, with the development of algorithmic programs such as Mathematica, the main problem is no longer the derivation of the KEO itself but rather the availability of a general expression of the KEO in curvilinear coordinates that is not too complicated and well adapted to the numerical methods used in
dynamics. To our knowledge, there are not many available general expressions of the KEO in curvilinear coordinates: the most famous one is given by the Wilson G matrix formulation [98, 99] (see also Ref. [100–103] for applications). One can also cite the formulation of Gu *at al.* [104, 105] and, finally, the polyspherical approach, which is the subject of the present review. They all correspond to different parametrizations of the N-body system. Regarding the polyspherical approach, Chapuisat and Iung [73] were perhaps the first to adopt this family of coordinates in quantum dynamics in a systematic way. The exact quantum mechanical formulation was given in a series of papers a little bit later [72, 74–77]. It should also be pointed out that Mladenović has given a very clear presentation of the equations of the polyspherical formulation in angular momenta [79, 106]. After this brief presentation of the general background, let us describe the polyspherical approach in detail.

## B. Vector parametrization and properties of angular momenta

As aforementioned, in the absence of external fields, it is possible to separate the translational degrees of freedom from the rest of the kinetic energy and to describe the configuration of an *N*-particle system in the SF frame, by *N*−1 relative position vectors. In all that follows, we no longer use mass-weighted Cartesian coordinates but, instead, the standard Cartesian coordinates. In addition, we confine ourselves to studying isolated molecular systems. Note, however, that if an external field is present or if the molecule is located in a noble gas matrix or on metal surface, all the results subsequently presented remain true. One simply has (i) to add one vector more: $\vec{OG}$, O being the center of the Laboratory-Fixed (LF) frame and G the center of mass of the molecular system, and (ii) to add the kinetic energy of the center of mass of the molecular system to the final expression of the KEO.

As said in Sec. IIA, in the polyspherical approach, we start from the classical expression of the kinetic energy and quantize it afterwards. More precisely, we start with the description of molecular systems by *N*−1 vectors. To illustrate such vector parametrization, we consider two examples in Sec. IIB1. Thereafter, we introduce the angular momenta associated with these vectors and give a general expression of the KEO in terms of these angular momenta. The latter are quasi-momenta (see the definition of the quasi-momenta in Section IIA (iii)), and they may be looked upon as 'intermediates' in the polyspherical approach.
FIG. 1: Definition of the three Jacobi vectors for an ABCD system. Here, $g_1$ is the center of mass of the (AB) subsystem, $g_2$ is the center of mass of the (CD) subsystem.

<figure><img src="image_2.png" /></figure>

The introduction of these intermediates allows (i) to obtain a very compact expression of the KEO whatever the set of vectors. For orthogonal coordinates, this compact expression is given in Eq. (65) in Sec. II B 2. (ii) It allows to derive a general expression of the KEO for the family of standard polyspherical coordinates (Eq. (92) in Sec. IIC1 for orthogonal coordinates) without using differential calculus. It is also the introduction of the angular momenta that avoids the use of differential calculus. (iii) It allows to find straightforwardly a primitive basis set (e. g. the basis set of spherical harmonics associated with the angular momenta) discarding all singularities that may occur in the KEO. To illustrate our ideas we will consider in Section IIB 1 a four-atom system, ABCD, first using Jacobi vectors and, second, using valence vectors.

## 1. Examples

### (i) Jacobi vectors:

Let $\vec{r}_A$, $\vec{r}_B$, $\vec{r}_C$, and $\vec{r}_D$ be the position vectors from the center of mass G of the atoms A,B,C and D (see Fig. 1). Jacobi vectors point from one atomic group center of mass to another one. In general, the choice of the Jacobi vectors is not unique [107] and several 'clustering' schemes are possible, one of which is shown in Fig. 1. Thus, we will use the Jacobi vectors as shown in Figure 1:
$$ \vec{R}_1 = \vec{r}_B - \vec{r}_A, \quad \vec{R}_2 = \vec{r}_D - \vec{r}_C, \qquad (34) $$

$$ \vec{R}_3 = \frac{m_C}{m_C + m_D} \vec{r}_C + \frac{m_D}{m_C + m_D} \vec{r}_D - \frac{m_A}{m_A + m_B} \vec{r}_A - \frac{m_B}{m_A + m_B} \vec{r}_B. $$

Let {$G, \vec{e}_{xSF}, \vec{e}_{ySF}, \vec{e}_{zSF}$} be an inertial Space-Fixed (SF) frame, where $\vec{e}_{\lambda SF}$ ($\lambda = x, y, z$) are fixed unit vectors. The position vector $\vec{r}_A$, when referred to the SF frame, is given by

$$ \vec{r}_A = \sum_{\lambda=x,y,z} r_{A\lambda^{SF}} \vec{e}_{\lambda^{SF}}, \qquad (35) $$

where $r_{A\lambda^{SF}}$ ($\lambda = x, y, z$) are the SF Cartesian coordinates. The velocity vector of the atom A, with respect to the SF frame, is given by

$$ \dot{\vec{r}}_A = \sum_{\lambda=x,y,z} \dot{r}_{A\lambda^{SF}} \vec{e}_{\lambda^{SF}}, \qquad (36) $$

where the dots denote time-derivatives. Similar relations hold of course for the atoms B,C and D and the Jacobi vectors. The kinetic energy of the system relative to the center of mass G, i.e. relative to the SF frame, is given by

$$ 2T^{SF} = m_A \dot{\vec{r}}_A \dot{\vec{r}}_A + m_B \dot{\vec{r}}_B \dot{\vec{r}}_B + m_C \dot{\vec{r}}_C \dot{\vec{r}}_C + m_D \dot{\vec{r}}_D \dot{\vec{r}}_D, \qquad (37) $$

and which, in view of Eq. (34), can be rewritten as

$$ 2T^{SF} = \mu_1 \dot{\vec{R}}_1 \dot{\vec{R}}_1 + \mu_2 \dot{\vec{R}}_2 \dot{\vec{R}}_2 + \mu_3 \dot{\vec{R}}_3 \dot{\vec{R}}_3, \qquad (38) $$

with the corresponding 'reduced masses':

$$ \mu_1 = \frac{m_A m_B}{m_A + m_B}, \quad \mu_2 = \frac{m_C m_D}{m_C + m_D}, \qquad (39) $$

$$ \mu_3 = \frac{(m_C + m_D)(m_A + m_B)}{M_T}, $$

and with $M_T = m_A + m_B + m_C + m_D$. Equation (38) shows that the kinetic energy has a diagonal form, and this is why the Jacobi coordinates are said to be orthogonal coordinates [108]. In other words, the kinetic energy of a four-atom system, relative to its center of mass,
can be described as the kinetic energy of three fictitious particles of masses $\mu_1$, $\mu_2$, and $\mu_3$ and of the position vectors $\vec{R}_1$, $\vec{R}_2$, and $\vec{R}_3$. To each vector $\vec{R}_i$ ($i = 1, 2, 3$), it is possible to associate a conjugate momentum

$$P_{i\lambda^{SF}} = \frac{\partial \mathcal{L}}{\partial \dot{R}_{i\lambda^{SF}}} = \frac{\partial T^{SF}}{\partial \dot{R}_{i\lambda^{SF}}}, \lambda = x, y, z; \qquad (40)$$

where $\mathcal{L}$ is the Lagrangian of the system. It is also possible to define a partial angular momentum corresponding to each Jacobi vector:

$$\vec{L}_{i}^{SF} = \vec{R}_{i} \times \vec{P}_{i}^{SF}. \qquad (41)$$

This leads to

$$2T^{SF} = \sum_{i=1}^{3} \frac{(\vec{P}_{i}^{SF})^{2}}{\mu_{i}} = \sum_{i=1}^{3} \frac{P_{R_{i}}^{2}}{\mu_{i}} + \frac{(\vec{L}_{i}^{SF})^{2}}{\mu_{i}R_{i}^{2}}. \qquad (42)$$

Here,

$$R_i = ||\vec{R}_i|| \qquad (43)$$

denotes a vector length and

$$P_{R_i} = \vec{P}_i \cdot \vec{e}_{R_i} \qquad (44)$$

is a radial momentum

$$(P_{R_i} = \frac{\partial L}{\partial \dot{R}_i} = \frac{\partial T^{SF}}{\partial \dot{R}_i}). \qquad (45)$$

The momentum $\vec{P}_i$ is given by its radial and angular parts as

$$\vec{P}_i = P_{R_i} \vec{e}_{R_i} - \frac{\vec{e}_{R_i} \times \vec{L}_i}{R_i}, \qquad (46)$$

where

$$\vec{e}_{R_i} = \vec{R}_i / R_i \qquad (47)$$

denotes a unit vector along $\vec{R}_i$.

The total angular momentum of the molecular system ABCD is

$$\vec{J} = \vec{r}_A \times m_A \dot{\vec{r}}_A + \vec{r}_B \times m_B \dot{\vec{r}}_B + \vec{r}_C \times m_C \dot{\vec{r}}_C + \vec{r}_D \times m_D \dot{\vec{r}}_D. \qquad (48)$$
<figure><img src="image_3.png" /><figcaption>FIG. 2: Definition of the three valence vectors for an ABCD system.</figcaption></figure>

It is straightforward, but tedious to show that

$$ \vec{J} = \sum_{i=1,2,3} \vec{L}_i^{SF}. \qquad (49) $$

Eq. (49) is true whatever the set of N-1 vectors used to describe the system.

### (ii) Valence vectors:

The valence vectors are vectors which join one atom to another one. The choice of the valence vectors is not unique and we will use the following definitions (see Figure 2) :

$$ \vec{R}_1 = \vec{r}_B - \vec{r}_A, \quad \vec{R}_2 = \vec{r}_C - \vec{r}_A, \quad \vec{R}_3 = \vec{r}_D - \vec{r}_A. \qquad (50) $$

Inserting Eq. (50) in Eq. (37) leads to a non-diagonal expression of the kinetic energy :

$$ 2T^{SF} = \sum_{i,j=1}^{3} \mu_{ij} \dot{\vec{R}}_i \dot{\vec{R}}_j. \qquad (51) $$

Since the expression of $T^{SF}$ is non-diagonal, the valence coordinates are said to be non-orthogonal. Again, to each vector $\dot{\vec{R}}_i$ ($i = 1, 2, 3$), it is possible to associate a conjugate momentum (Eq. (40)) and to define a corresponding partial angular momentum (Eq.
(41)). Finally, after having inserted Eq. (50) into Eq. (37) and used Eq. (40), we obtain a non-diagonal expression of the kinetic energy in terms of the momenta

$$2T^{SF} = \sum_{i,j=1}^{3} \vec{P}_{i}^{SF} M_{ij} \vec{P}_{j}^{SF}, \qquad (52)$$

with the matrix M, which reads

$$M = \begin{bmatrix} \frac{1}{m_A} + \frac{1}{m_B} & \frac{1}{m_A} & \frac{1}{m_A} \\ \frac{1}{m_A} & \frac{1}{m_A} + \frac{1}{m_C} & \frac{1}{m_A} \\ \frac{1}{m_A} & \frac{1}{m_A} & \frac{1}{m_A} + \frac{1}{m_D} \end{bmatrix}. \qquad (53)$$

## 2. General formulation

(i) Defining a set of $N-1$ vectors and the corresponding classical kinetic energy: Let $\vec{R}_1, \vec{R}_2, \dots, \vec{R}_{N-1}$ be any set of vectors chosen for the description of a molecular system. $\vec{R}_1, \vec{R}_2, \dots, \vec{R}_{N-1}$ can be orthogonal vectors such as Jacobi vectors (i.e. vectors pointing from one atomic group center of mass to another one) or Radau vectors (see Sec. III A 3 for the definition of the Radau vectors) or non-orthogonal vectors such as valence vectors joining two atoms. The classical SF kinetic energy can be written as

$$2T^{SF} = \sum_{i,j=1}^{N-1} \vec{P}_{i}^{SF} M_{ij} \vec{P}_{j}^{SF} = \sum_{i,j=1,...,N-1} \sum_{\lambda=x,y,z} P_{i\lambda SF} M_{ij} P_{j\lambda SF}, \qquad (54)$$

where $M_{ij}$ are the elements of the mass matrix M. We often use, mainly for orthogonal coordinates, the notation

$$\mu_i = \frac{1}{M_{ii}} \qquad (55)$$

for the diagonal elements of the matrix M. $\mu_i$ is called the 'reduced mass' associated with the vector $\vec{R}_i$. To determine M, we introduce a (redundant) set of position vectors, $\vec{r}_1, \vec{r}_2, \dots, \vec{r}_N$, where $\vec{r}_i$ points from the center of mass G to the position of the i-th atom. The kinetic energy, relative to the center of mass, is simply given by

$$2T^{SF} = \sum_{i=1,...,N} \sum_{\lambda=x,y,z} \frac{(p_{i\lambda SF})^2}{m_i}, \qquad (56)$$
where $\vec{p}_i$ is the conjugate momentum associated with $\vec{r}_i$, and $m_i$ the mass of the i-th atom. Let us now return to a set of $(N-1)$ internal vectors, $\vec{R}_1, \vec{R}_2, \dots, \vec{R}_{N-1}$, chosen to describe the same N-atom system. Then there exists an $(N-1) \times N$ matrix $\mathbf{A}$, which connects the two sets of coordinates (see Eqs. (34) and (50) for two particular cases of an ABCD system):

$$ \begin{bmatrix} \vec{R}_1 \\ \vec{R}_2 \\ \dots \\ \vec{R}_{N-1} \end{bmatrix} = \mathbf{A} \begin{bmatrix} \vec{r}_1 \\ \vec{r}_2 \\ \dots \\ \vec{r}_N \end{bmatrix} . \qquad (57) $$

The symmetric mass-dependent constant matrix $\mathbf{M}$ for the internal set of $N-1$ vectors appearing in the KEO is given by $\mathbf{M} = \mathbf{A}\mathbf{m}^{-1}\mathbf{A}^T$ where $\mathbf{m}$ denotes the diagonal matrix of particle masses. In this formulation, it is thus very easy to switch from one set of vectors to another, just by changing the matrix $\mathbf{M}$.

For valence and Jacobi vectors, $\mathbf{M}$ can be obtained directly by using the following general property [73]:

(1) The diagonal element $M_{ii}$ is the inverse of the reduced mass associated with the two groups of particles whose centers of mass are joined by $\vec{R}_i$.

For instance, in the case of the ABCD system in the Jacobi vectors of Figure 1, we get

$$ \begin{aligned} M_{11} &= \frac{1}{m_A} + \frac{1}{m_B} = \frac{1}{\mu_1}, \\ M_{22} &= \frac{1}{m_C} + \frac{1}{m_D} = \frac{1}{\mu_2}, \\ M_{33} &= \frac{1}{m_A + m_B} + \frac{1}{m_C + m_D} = \frac{1}{\mu_3}. \end{aligned} \qquad (58) $$

In contrast, using the valence vectors of Figure 2, we obtain

$$ \begin{aligned} M_{11} &= \frac{1}{m_A} + \frac{1}{m_B}, \\ M_{22} &= \frac{1}{m_A} + \frac{1}{m_C}, \\ M_{33} &= \frac{1}{m_A} + \frac{1}{m_D}. \end{aligned} \qquad (59) $$

(2) For the off-diagonal elements, $M_{ij}$ is zero for Jacobi vectors and for valence vectors with no common atom. Otherwise, $|M_{ij}|$ is the inverse of the mass of the atom common to both $\vec{R}_i$ and $\vec{R}_j$. The sign of $M_{ij}$ is + if the two vectors originate from or point towards the same atom; it is - in the other case.
<figure><img src="image_4.png" /><figcaption>FIG. 3: (AB)CB(EF) system parametrised by five valence vectors.</figcaption></figure>

For instance, in the case of the ABCD system in Jacobi vectors of Figure 1, all the off-diagonal elements are equal to zero, while for the valence vectors described in Figure 2, we obtain $M_{12} = M_{13} = M_{23} = + \frac{1}{m_A}$.

Let us give a final example: in the case of an (AB)CB(EF) system parametrized by five valence vectors as depicted in Figure 3, the matrix M can be obtained using the previous result without any calculation:

$$ M = \begin{bmatrix} \frac{1}{m_A} + \frac{1}{m_C} & + \frac{1}{m_C} & 0 & 0 & + \frac{1}{m_C} \\ + \frac{1}{m_C} & \frac{1}{m_B} + \frac{1}{m_C} & 0 & 0 & + \frac{1}{m_C} \\ 0 & 0 & \frac{1}{m_D} + \frac{1}{m_E} & + \frac{1}{m_D} & - \frac{1}{m_D} \\ 0 & 0 & + \frac{1}{m_D} & \frac{1}{m_D} + \frac{1}{m_F} & - \frac{1}{m_D} \\ + \frac{1}{m_C} & + \frac{1}{m_C} & - \frac{1}{m_D} & - \frac{1}{m_D} & \frac{1}{m_C} + \frac{1}{m_D} \end{bmatrix} . \quad (60) $$

In the general case, one should calculate the matrix A and derive the matrix M from it.

### (ii) Introduction of the Body-Fixed frame and quantization:

In order to separate the overall rotation from the internal deformation, let us now introduce
a Body-Fixed (BF) frame, {$G$; $\vec{e}_{xBF}$, $\vec{e}_{yBF}$, $\vec{e}_{zBF}$}. Its orientation with respect to the SF frame is determined by three Euler angles: $\alpha$, $\beta$, $\gamma$. After definition of the three Euler angles, the shape and the size of the molecules are described by $3N-6$ internal BF coordinates. At this level, the exact definition of the BF frame is not yet specified (and not necessary), but the classical kinetic energy can already be recast as follows:

$$2T^{SF} = \sum_{i,j=1,...,N-1} \sum_{\lambda=x,y,z} P_{i\lambda BF} M_{ij} P_{j\lambda BF} , \qquad (61)$$

where $P_{i\lambda BF}$ are the Cartesian components of the SF (and not BF! see below) conjugate momenta in the BF frame [81], or in other words, the projections of the SF conjugate momenta onto the BF-axes:

$$\vec{P}_i = P_{ixBF} \vec{e}_{xBF} + P_{iyBF} \vec{e}_{yBF} + P_{izBF} \vec{e}_{zBF} . \qquad (62)$$

According to Refs. [58, 81], the exact quantum mechanical counterpart of Eq. (61) is given by:

$$2\hat{T} = \sum_{i,j=1,...,N-1} \sum_{\lambda=x,y,z} \hat{P}_{i\lambda BF}^{\dagger} M_{ij} \hat{P}_{j\lambda BF} . \qquad (63)$$

where the adjoint $\hat{P}_{i\lambda BF}^{\dagger}$ is given by Eq. (9).

### (iii) Introduction of the Body-Fixed projections of the angular momenta associated with the N-1 vectors

It should be clear that the kinetic energy operator will always be the Space-Fixed one, and not the Body-Fixed kinetic energy, even after introduction of the BF frame. Furthermore, the quantum mechanical counterpart of Eq. (46) is given by

$$\hat{P}_i = \hat{P}_{R_i} \vec{e}_i - \frac{\vec{e}_i \times \hat{L}_i}{R_i} . \qquad (64)$$

After the introduction of a BF frame, projections of the angular momenta onto BF axes will be used. For instance, in the special case of orthogonal vectors, Eq. (63) yields (here, $\mu_i = \frac{1}{M_{ii}}$):

$$2\hat{T} = \sum_{i=1}^{N-1} \frac{\hat{P}_{R_i}^{\dagger} \hat{P}_{R_i}}{\mu_i} + \frac{(\hat{L}_i^{\dagger} \cdot \hat{L}_i)_{BF}}{\mu_i R_i^2} . \qquad (65)$$

If the set of vectors is non-orthogonal, the structure of the kinetic energy operator is more complex (with crossing terms between the radial and angular momenta proportional to
the masses $M_{ij}$, with $i \neq j$: see Appendix B). At this level, we have already a general formulation of the kinetic energy operator for any set of vectors expressed in terms of the BF projections of the momenta associated with these vectors. The $3N-3$ BF-projections of the angular momenta are quasi-momenta and Eq. (65) appears as a particular case of Eq. (23) of section II A.

Let us explain why it is necessary to introduce the BF projections of the momenta when making a change of coordinates to separate the overall rotation from the internal vibrations. The action of SF-components of the angular momenta on the primitive basis functions expressed in BF coordinates is very complicated, whereas the action of the BF-components onto the primitive basis functions expressed in BF coordinates is rather simple, see, e. g., the action of the KEO onto a basis set of spherical harmonics of BF spherical coordinates discussed in Sec. IIC 1 below. As BF coordinates are needed to separate the overall rotation from the internal vibrations, it becomes necessary to express the KEO in terms of BF-components of momenta and angular momenta as in Eqs. (63, 65), rather than in terms of the SF-components.

It is to be emphasized again that angular momenta that will appear in the kinetic energy operators (such as in Eq. (65)) are all computed in the SF frame, but are projected onto the axes of several frames (e. g. the BF frame in Eq. (65)):

$$ \begin{bmatrix} \hat{L}_{ix^{BF}} \\ \hat{L}_{iy^{BF}} \\ \hat{L}_{iz^{BF}} \end{bmatrix} = \mathcal{R}^T(\alpha, \beta, \gamma) \begin{bmatrix} \hat{L}_{ix^{SF}} \\ \hat{L}_{iy^{SF}} \\ \hat{L}_{iz^{SF}} \end{bmatrix}, \qquad (66) $$

where $\mathcal{R}(\alpha, \beta, \gamma)$ is the Euler rotation matrix [109] (see Eq. (70) below for the explicit expression of this matrix). As aforementioned, the introduction of these projections is necessary when using the BF coordinates, but raises a new technical problem. The projections of the angular momenta onto the SF-axes satisfy the usual commutation relations, and their action onto a basis set of spherical harmonics in terms of the SF spherical coordinates is well known (see for instance Ref. [109]). However, the projections of the same angular momenta onto the axes of a moving frame may satisfy anomalous, non-definite commutation relations [81]. Luckily, this problem does not occur for all angular momenta. If a vector is not involved in the definition of a frame F (see Section IIC 1 (i) to understand what we exactly
mean by 'a vector involved in the definition of a frame'), the expression for the projection of the corresponding angular momentum onto the F-axes expressed in the coordinates in this frame is identical to the usual one in an SF frame [58]. For instance, if one vector $\vec{R}_j$ is not involved in the definition of the BF frame, the projections of $\vec{L}_j$ onto the BF-axes and expressed in terms of the BF-Cartesian coordinates are given by:

$$ \begin{bmatrix} \hat{L}_{j x^{BF}} \\ \hat{L}_{j y^{BF}} \\ \hat{L}_{j z^{BF}} \end{bmatrix} = \begin{bmatrix} R_{j x^{BF}} \\ R_{j y^{BF}} \\ R_{j z^{BF}} \end{bmatrix} \times \begin{bmatrix} \frac{1}{i} \frac{\partial}{\partial R_{j x^{BF}}} \\ \\ \frac{1}{i} \frac{\partial}{\partial R_{j y^{BF}}} \\ \\ \frac{1}{i} \frac{\partial}{\partial R_{j z^{BF}}} \end{bmatrix} . \qquad (67) $$

This expression is identical to the expression of the projections of an angular momentum onto the SF-axes in terms of SF coordinates. The origin of this property is explained in Appendix A. It is easy to verify that the usual commutation relations,

$$ [\hat{L}_{j \lambda^{BF}}, \hat{L}_{j \nu^{BF}}] = i \sum_{\rho} \epsilon_{\lambda \nu \rho} \hat{L}_{j \rho^{BF}} , \qquad (68) $$

are satisfied, where $\epsilon$ denotes the well-known totally antisymmetric tensor, i.e. $\epsilon_{\lambda\nu\rho} = 1$ ($\epsilon_{\lambda\nu\rho} = -1$) if $\{\lambda\nu\rho\}$ is an even (odd) permutation of $\{xyz\}$, and zero otherwise. This property (especially Eq. (67)) is remarkable and very helpful, since combined with the previous vector parametrization of N-atom systems, it will allow us to derive kinetic energy operators in a compact and general form with no need to use differential calculus.

### C. A general expression of the KEO in standard polyspherical coordinates

Until now, we have described an N-particle system by a set of N-1 vectors and provided the expression of the KEO in terms of the conjugate momenta associated with these vectors. We have explained in Sec. II B 2 (i) how to switch from one set of vectors to another one. Thereafter, we have introduced the BF projections of the angular momenta associated with the N-1 vectors in Sec. II B 2 (ii) and brought out some general properties regarding these projections in Sec. II B 2 (iii). All of this constitutes a framework that makes possible to derive a general expression of the KEO in terms of the spherical coordinates that parametrize the vectors in the BF frame whose precise definition must be given now.
1. Expression in angular momenta

<figure><img src="image_5.png" /><figcaption>FIG. 4: Definition of the $E_2$ Frame and of the first two Euler angles.</figcaption></figure>

(i) Definition of the BF frame: Fig. 4, 5, and 6.

The orientation of the BF frame with respect to the SF frame is determined by the Euler angles $\alpha, \beta, \gamma$. There are numerous ways of specifying Euler angles. The one chosen here is that of Biedenharn and Louck [110] where $\alpha$ and $\beta$ are the ordinary spherical coordinates of the $\vec{e}_z^{BF}$ vector in the SF system and $\gamma$ is an angle measuring a clockwise rotation about $\vec{e}_z^{BF}$. To be very precise, it means that, for a given vector $\vec{R}$, the relation between SF and BF components of this vector is given by

$$ \begin{bmatrix} R_{x^{SF}} \\ R_{y^{SF}} \\ R_{z^{SF}} \end{bmatrix} = \mathcal{R}(\alpha, \beta, \gamma) \begin{bmatrix} R_{x^{BF}} \\ R_{y^{BF}} \\ R_{z^{BF}} \end{bmatrix}, \qquad (69) $$

with $\mathcal{R}(\alpha, \beta, \gamma)$ the Euler rotation matrix, which takes the following form in the convention of Biedenharn and Louck [110] (see also Ref. [109]):
<figure><img src="image_6.png" /><figcaption>FIG. 5: Definition of the BF frame and of the third angle.</figcaption></figure>

$$ \begin{aligned} \mathcal{R}(\alpha, \beta, \gamma) &= D_z(\alpha) \times D_y(\beta) \times D_z(\gamma) \\ &= \begin{bmatrix} \cos \alpha & -\sin \alpha & 0 \\ \sin \alpha & \cos \alpha & 0 \\ 0 & 0 & 1 \end{bmatrix} \times \begin{bmatrix} \cos \beta & 0 & \sin \beta \\ 0 & 1 & 0 \\ -\sin \beta & 0 & \cos \beta \end{bmatrix} \times \begin{bmatrix} \cos \gamma & -\sin \gamma & 0 \\ \sin \gamma & \cos \gamma & 0 \\ 0 & 0 & 1 \end{bmatrix}. \end{aligned} \quad (70) $$

Throughout this Section, a very particular definition of the BF frame is used:

(1) Let {$G; \vec{e}_{xE_1}, \vec{e}_{yE_1}, \vec{e}_{zE_1}$} be the $E_1$ frame resulting from the first Euler rotation, i.e. whose origin coincides with G, where $\vec{e}_{zE_1} = \vec{e}_{zSF}$ and where $\vec{e}_{xE_1}$ and $\vec{e}_{yE_1}$ are obtained by a rotation through the angle $\alpha$ about $\vec{e}_{zE_1}$.

(2) Let {$G; \vec{e}_{xE_2}, \vec{e}_{yE_2}, \vec{e}_{zE_2}$} be the $E_2$ frame resulting from the first two Euler rotations, where $\vec{e}_{yE_2} = \vec{e}_{yE_1}$ and where $\vec{e}_{xE_2}$ and $\vec{e}_{zE_2}$ are obtained by a rotation through the angle $\beta$ about $\vec{e}_{yE_2}$ (see Fig. 4 and 6).

(3) Let {$G; \vec{e}_{xE_3}, \vec{e}_{yE_3}, \vec{e}_{zE_3}$} = {$G; \vec{e}_{xBF}, \vec{e}_{yBF}, \vec{e}_{zBF}$} be the $E_3$ or BF frame, i.e. the frame resulting from the three Euler rotations, where $\vec{e}_{zE_3} = \vec{e}_{zE_2}$, and where $\vec{e}_{xE_3}$ and $\vec{e}_{yE_3}$ are obtained by a rotation through the angle $\gamma$ about $\vec{e}_{zE_2}$.
(4) The Euler angles are chosen such that the BF or $E_3$ frame is oriented in a way that $\vec{e}_{zBF} = \vec{e}_{zE_2}$ is parallel to the vector $\vec{R}_{N-1}$ and that $\vec{R}_{N-2}$, $\vec{e}_{xBF}$ and $\vec{e}_{zBF}$ lie in the same half-plane (see Figs. 5 and 6).
FIG. 6: Definition of the BF frame and of the three Euler angles. The $E_2$ Frame is the frame obtained after the first two Euler rotations.

<figure><img src="image_7.png" /></figure>

We have thus chosen a very specific definition of the BF frame (note, however, that the choice of the set of vectors remains free). It is important to notice that this definition of the BF frame has an important advantage: *only two vectors* are involved in this definition of the BF frame. This property will greatly simplify the quantum mechanical expression of the kinetic energy operator. Indeed, as outlined in Sec. II B 2 (iii), the expressions for the angular momenta associated with the $N-3$ remaining vectors are the usual ones when projected onto the BF-axes and expressed in terms of BF-coordinates. It should be pointed out that this is very different from the Eckart frame [111], which minimizes the Coriolis coupling. The definition of this latter frame, however, depends on *all* vectors rather than on two. Consequently, using the Eckart frame in general leads to a much more complex expression of the KEO. This is particularly true for the part of the KEO containing the Coriolis coupling whose expression directly depends upon the definition of the BF frame. On the other hand, the KEO for $J=0$ is independent of the definition of the BF frame (see II A (iii)).

Generally speaking, it is advantageous to choose the vector $\vec{R}_{N-1}$ (that defines the $z^{BF}$
axis) such that its 'reduced mass',

$$ \mu_{N-1} = \frac{1}{M_{N-1N-1}} \qquad (71) $$

in Eqs. (54,63,65), is larger than the reduced masses associated with the other vectors ($\mu_i$, $i = 1, \dots, N-2$). Indeed, we shall see later that the part of the KEO containing the Coriolis coupling is proportional to

$$ \frac{1}{\mu_{N-1}} = M_{N-1N-1}. \qquad (72) $$

If it is possible to find such a vector, one reduces the Coriolis coupling without resorting to a very complex expression of the KEO. It is when such a vector cannot be found that changing the definition of the BF frame might be appropriate. In the present review, we do not want to tackle the subject of how to change this definition, and we simply refer the reader to the article of Wang and Carrington [112], which deals with this issue (see also Ref. [113] and Sec. IV C 3 (ii) for examples).

### (ii) Polyspherical parametrization.

Until now, we have not defined the 3N-6 internal coordinates that describe the deformations of the molecular systems. Let us now parametrize the vectors by spherical coordinates called the 'polyspherical coordinates'. The spherical coordinates in a given frame F are defined as

$$ R_{i x^F} = R_i \sin \theta_i^F \cos \varphi_i^F, \quad R_{i y^F} = R_i \sin \theta_i^F \sin \varphi_i^F, \quad R_{i z^F} = R_i \cos \theta_i^F, $$

with $0 \le \theta_i^F \le \pi$ and $0 \le \varphi_i^F < 2\pi$. The parametrization of the standard polyspherical type consists of three Euler angles for the overall rotation of the BF frame and 3N-6 internal coordinates. With our definition of the BF frame, it is clear that the first two Euler angles simply are the two spherical angles of $\vec{R}_{N-1}$ in the SF frame, $\alpha = \varphi_{N-1}^{SF}$ and $\beta = \theta_{N-1}^{SF}$. The third Euler angle is given by $\gamma = \varphi_{N-2}^{E2}$. The other 3N-6 coordinates are the BF spherical coordinates, i. e. the $N-1$ vector lengths $R_i \in [0, \infty)$, $N-2$ (BF) planar angles $\theta_i^{BF} \in [0, \pi]$ between the vectors $\vec{R}_{N-1}$ and $\vec{R}_i$, and $N-3$ (BF) dihedral angles $\varphi_i^{BF} \in [0, 2\pi)$ between the two vectors $\vec{R}_i$, $\vec{R}_{N-2}$ around the vector $\vec{R}_{N-1}$. They are depicted in Figure 7. If an Euclidean convention of normalization is used [58], the volume element reads

$$ d\tau = R_{N-1}^2 \sin \beta \, dR_{N-1} d\alpha \, d\beta R_{N-2}^2 \sin \theta_{N-2}^{BF} \, dR_{N-2} d\gamma \, d\theta_{N-2}^{BF} \\ \times \prod_{i=1}^{N-3} R_i^2 \sin \theta_i^{BF} \, dR_i \, d\varphi_i^{BF} \, d\theta_i^{BF}, \qquad (73) $$
i.e. it is the usual Euclidean volume element:

$$d\tau = J_e \, dR_{N-1} \, d\alpha \, d\beta \, dR_{N-2} \, d\gamma \, d\theta_{N-2}^{BF} \prod_{i=1}^{N-3} dR_i \, d\varphi_i^{BF} \, d\theta_i^{BF}, \qquad (74)$$

where $J_e$ is the absolute value of the Jacobian (see Eq. (3)) corresponding to the coordinate transformation from the Cartesian coordinates to the polyspherical coordinates.

<figure><img src="image_8.png" /><figcaption>FIG. 7: Definition of the standard polyspherical coordinates: the $N-1$ vector lengths $R_i \in [0, \infty)$, $N-2$ (BF) planar angles $\theta_i^{BF} \in [0, \pi]$, and $N-3$ (BF) dihedral angles $\varphi_i^{BF} \in [0, 2\pi)$.</figcaption></figure>

### (iii) Properties of the BF-projections of the angular momenta.

In the expression of the kinetic energy operator in Eq. (65), $N-1$ angular momenta appear (the same is true for non-orthogonal coordinates). Since only $\vec{R}_{N-1}$ and $\vec{R}_{N-2}$ are involved in the definition of the BF frame, the projections of all the $N-3$ angular momenta $\vec{L}_1, \dots,$
$\hat{\vec{L}}_{N-3}$ onto the BF axes are regular. They are characterized by the usual formula [109, 110]

$$ \begin{bmatrix} \hat{\vec{L}}_{ix^{BF}} \\ \hat{\vec{L}}_{iy^{BF}} \\ \hat{\vec{L}}_{iz^{BF}} \end{bmatrix} = \begin{bmatrix} -\sin \varphi_i^{BF} & -\cos \varphi_i^{BF} \cot \theta_i^{BF} \\ \cos \varphi_i^{BF} & -\sin \varphi_i^{BF} \cot \theta_i^{BF} \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \frac{1}{i} \frac{\partial}{\partial \theta_i^{BF}} \\ \frac{1}{i} \frac{\partial}{\partial \varphi_i^{BF}} \end{bmatrix}. \quad (75) $$

Hence, they obey the usual commutation relations and are self-adjoint.

Now, what about the two last angular momenta, $\hat{\vec{L}}_{N-1}$ and $\hat{\vec{L}}_{N-2}$? First, it is trivial to eliminate $\hat{\vec{L}}_{N-1}$ by substituting $\hat{\vec{L}}_{N-1}$ with $\hat{\vec{J}} - \sum_{i=1}^{N-2} \hat{\vec{L}}_i$ where $\hat{\vec{J}}$ is the total angular momentum, which is a constant of motion when no external field is present. Considering the particular case of orthogonal vectors, one can recast Eq. (65) as

$$ 2\hat{T} = \sum_{i=1}^{N-1} \frac{\hat{P}_{Ri}^{\dagger} \hat{P}_{Ri}}{\mu_i} + \sum_{i=1}^{N-2} \frac{(\hat{\vec{L}}_i^{\dagger} \cdot \hat{\vec{L}}_i)_{BF}}{\mu_i R_i^2} + \frac{((\hat{\vec{J}}^{\dagger} - \sum_{i=1}^{N-2} \hat{\vec{L}}_i^{\dagger}) \cdot (\hat{\vec{J}} - \sum_{i=1}^{N-2} \hat{\vec{L}}_i))_{BF}}{\mu_{N-1} R_{N-1}^2}. \quad (76) $$

The expression of the BF-projections of the total angular momentum in terms of the Euler angles is well-known [109, 110, 114]:

$$ \begin{bmatrix} \hat{J}_{x^{BF}} \\ \hat{J}_{y^{BF}} \\ \hat{J}_{z^{BF}} \end{bmatrix} = \begin{bmatrix} -\frac{\cos \gamma}{\sin \beta} & \sin \gamma & \cot \beta \cot \gamma \\ \frac{\sin \gamma}{\sin \beta} & \cos \gamma & -\cot \beta \sin \gamma \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \frac{1}{i} \frac{\partial}{\partial \alpha} \\ \frac{1}{i} \frac{\partial}{\partial \beta} \\ \frac{1}{i} \frac{\partial}{\partial \gamma} \end{bmatrix}, \quad (77) $$

and we obtain $[\hat{J}_{\lambda^{BF}}, \hat{J}_{\nu^{BF}}] = -i \sum_{\rho} \epsilon_{\lambda\nu\rho} \hat{J}_{\rho^{BF}}$, which are the well-known anomalous commutation relations [114]. It is also worth noting that all the components are self-adjoint.

Now, what about the angular momentum $\hat{\vec{L}}_{N-2}$ associated with $\vec{R}_{N-2}$? Since $\vec{R}_{N-2}$ is involved in the definition of the BF frame, nothing can be said *a priori* about the properties of the projections of $\hat{\vec{L}}_{N-2}$ onto the BF-axes. When such a case occurs, one proceeds as follows. One starts with the projections of the same angular momentum but onto another frame. Here, it is judicious to start with the projections of $\hat{\vec{L}}_{N-2}$ onto the $E_2$-axes, because $\vec{R}_{N-2}$ is not used when defining the $E_2$ frame. (We recall that the $E_2$ frame is the frame resulting from the first two Euler rotations only, see Fig. 4). Indeed, $\vec{R}_{N-2}$ is needed only
for defining the third Euler angle. Consequently, the $E_2$-components of $\hat{\vec{L}}_{N-2}$ have a regular form

$$ \begin{bmatrix} \hat{L}_{N-2x^{E_2}} \\ \hat{L}_{N-2y^{E_2}} \\ \hat{L}_{N-2z^{E_2}} \end{bmatrix} = \begin{bmatrix} -\sin \varphi_{N-2}^{E_2} & -\cos \varphi_{N-2}^{E_2} \cot \theta_{N-2}^{E_2} \\ \cos \varphi_{N-2}^{E_2} & -\sin \varphi_{N-2}^{E_2} \cot \theta_{N-2}^{E_2} \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \frac{1}{i} \frac{\partial}{\partial \theta_{N-2}^{E_2}} \\ \frac{1}{i} \frac{\partial}{\partial \varphi_{N-2}^{E_2}} \end{bmatrix}. \quad (78) $$

These components obey the usual commutation relations and are self-adjoint. The BF-projections of $\hat{\vec{L}}_{N-2}$ are obtained by rotating the $E_2$-projections of $\hat{\vec{L}}_{N-2}$ through the angle $\varphi_{N-2}^{E_2} = \gamma$ around the $z^{E_2}$ axis

$$ \begin{bmatrix} \hat{L}_{N-2x^{BF}} \\ \hat{L}_{N-2y^{BF}} \\ \hat{L}_{N-2z^{BF}} \end{bmatrix} = \begin{bmatrix} \cos \varphi_{N-2}^{E_2} & \sin \varphi_{N-2}^{E_2} & 0 \\ -\sin \varphi_{N-2}^{E_2} & \cos \varphi_{N-2}^{E_2} & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \hat{L}_{N-2x^{E_2}} \\ \hat{L}_{N-2y^{E_2}} \\ \hat{L}_{N-2z^{E_2}} \end{bmatrix}. \quad (79) $$

A final change of coordinates is then performed

$$ \begin{align} \gamma &= \varphi_{N-2}^{E_2} , \tag{80} \\ \theta_{N-2}^{BF} &= \theta_{N-2}^{E_2} , \nonumber \\ \theta_i^{BF} &= \theta_i^{E_2}, \quad i = 1, \dots, N-3 , \nonumber \\ \varphi_i^{BF} &= \varphi_i^{E_2} - \varphi_{N-2}^{E_2}, \quad i = 1, \dots, N-3 . \nonumber \end{align} $$

Applying the chain rule to $\frac{\partial}{\partial \theta_{N-2}^{E_2}}$ and $\frac{\partial}{\partial \varphi_{N-2}^{E_2}}$, gives respectively

$$ \begin{align} \frac{\partial}{\partial \theta_{N-2}^{E_2}} &= \frac{\partial}{\partial \theta_{N-2}^{BF}} , \tag{81} \\ \frac{\partial}{\partial \varphi_{N-2}^{E_2}} &= \frac{\partial}{\partial \gamma} - \sum_{i=1}^{N-3} \frac{\partial}{\partial \varphi_i^{BF}} . \nonumber \end{align} $$
Inserting these relations into Eq. (79), we obtain

$$ 
\begin{bmatrix} \hat{L}_{N-2} x^{BF} \\ \hat{L}_{N-2} y^{BF} \\ \hat{L}_{N-2} z^{BF} \end{bmatrix} = \begin{bmatrix} -\cot \theta_{N-2}^{BF} & 0 & \cot \theta_{N-2}^{BF} \\ 0 & 1 & 0 \\ 1 & 0 & -1 \end{bmatrix} \begin{bmatrix} \frac{1}{i} \frac{\partial}{\partial \gamma} \\ \frac{1}{i} \frac{\partial}{\partial \theta_{N-2}^{BF}} \\ \sum_{j=1}^{N-3} \frac{1}{i} \frac{\partial}{\partial \varphi_j^{BF}} \end{bmatrix}, \quad (82) $$

i. e.

$$ 
\begin{bmatrix} \hat{L}_{N-2} x^{BF} \\ \hat{L}_{N-2} y^{BF} \\ \hat{L}_{N-2} z^{BF} \end{bmatrix} = \begin{bmatrix} -\cot \theta_{N-2}^{BF} (\hat{J}_{zBF} - \sum_{i=1}^{N-3} \hat{L}_{izBF}) \\ -i \frac{\partial}{\partial \theta_{N-2}^{BF}} \\ \hat{J}_{zBF} - \sum_{i=1}^{N-3} \hat{L}_{izBF} \end{bmatrix}, \quad (83) $$

which satisfy non-definite commutation relations

$$ 
\begin{align} 
[\hat{L}_{N-2} x^{BF}, \hat{L}_{N-2} y^{BF}] &= \frac{1}{(\sin \theta_{N-2}^{BF})^2} \left( \frac{\partial}{\partial \gamma} - \sum_{i=1}^{N-3} \frac{\partial}{\partial \varphi_i^{BF}} \right) , \nonumber \\ 
[\hat{L}_{N-2} y^{BF}, \hat{L}_{N-2} z^{BF}] &= 0 , \tag{84} \\ 
[\hat{L}_{N-2} z^{BF}, \hat{L}_{N-2} x^{BF}] &= 0 . \nonumber 
\end{align} $$

The y-BF component of $\hat{\vec{L}}_{N-2}$ is not Hermitian and its adjoint, $\hat{\vec{L}}_{N-2}^{\dagger}$, is given by the following equation

$$ (\hat{\vec{L}}_{N-2})^{\dagger} = \hat{\vec{L}}_{N-2} + \begin{pmatrix} 0 \\ i \cot \theta_{N-2}^{BF} \\ 0 \end{pmatrix} . \quad (85) $$

We have used Eqs. (10,73) to obtain Eq. (85).

Two points must be addressed here: (1) the BF components of $\hat{\vec{L}}_{N-2}$ do not commute with those of the other angular momenta (see Eq. (83): $\hat{L}_{N-2} x^{BF}$ and $\hat{L}_{N-2} z^{BF}$ depend upon $\hat{J}_{zBF}$ and $\hat{L}_{izBF}$ ($i = 1, \dots, N-3$) so that they do not commute with $\hat{J}_{\lambda BF}$, $\hat{L}_{i\lambda BF}$ ($i = 1, \dots, N-3$), $\lambda = x, y$), (2) the fact that we have used an intermediate frame
(the $E_2$ frame here) to determine the projections of $\hat{\vec{L}}_{N-2}$ onto the BF-axes is of general character and can be used to determine the projections of other angular momenta with a non-regular behavior.

(iv) General expression of the KEO in polyspherical coordinates.

If we keep in mind that $\hat{\vec{L}}_{N-2}$ has a particular behavior, we can establish that [74, 76]:

$$ \begin{aligned} \hat{\vec{L}}_{N-1}^{\dagger} \hat{\vec{L}}_{N-1} &= \left( (\hat{\vec{J}}^{\dagger} - \sum_{i=1}^{N-2} \hat{\vec{L}}_{i}^{\dagger}) \cdot (\hat{\vec{J}} - \sum_{i=1}^{N-2} \hat{\vec{L}}_{i}) \right) \\ &= \left( \hat{\vec{J}}^{\dagger} \hat{\vec{J}} + (\sum_{i=1}^{N-2} \hat{\vec{L}}_{i})(\sum_{j=1}^{N-2} \hat{\vec{L}}_{j}) - 2(\sum_{i=1}^{N-2} \hat{\vec{L}}_{i}) \hat{\vec{J}} \right)_{E_2} \\ &= \left( \hat{\vec{J}}^2 + (\sum_{i=1}^{N-2} \hat{\vec{L}}_{i}^{\dagger})(\sum_{j=1}^{N-2} \hat{\vec{L}}_{j}) - 2\hat{\vec{J}}(\sum_{i=1}^{N-2} \hat{\vec{L}}_{i}) \right)_{BF} \\ &= \left( \hat{\vec{J}}^2 + \sum_{i=1}^{N-2} \hat{\vec{L}}_{i}^{\dagger} \hat{\vec{L}}_{i} + \sum_{i=1}^{N-2} \sum_{j=1, j \neq i}^{N-2} \hat{\vec{L}}_{i}^{\dagger} \hat{\vec{L}}_{j} - 2\hat{\vec{J}}(\sum_{i=1}^{N-2} \hat{\vec{L}}_{i}) \right)_{BF} \\ &= \left( \hat{\vec{J}}^2 + \sum_{i=1}^{N-2} \hat{\vec{L}}_{i}^{\dagger} \hat{\vec{L}}_{i} + 2 \sum_{i=1}^{N-3} \sum_{j=1, j>i}^{N-2} \hat{\vec{L}}_{i} \hat{\vec{L}}_{j} - 2\hat{\vec{J}}(\sum_{i=1}^{N-2} \hat{\vec{L}}_{i}) \right)_{BF} . \end{aligned} \quad (86) $$

Note that the order of the operators in the scalar products is strictly fixed! The key point is the fact that the BF-components of $\hat{\vec{L}}_{N-2}$ must remain on the right hand side of all the scalar products. The first line of Eq. (86) is correct whatever the frame used to calculate the scalar product. The parentheses $(\cdot)_{F}$ indicate that all the scalar products inside the parentheses are calculated using the F-components of the angular momenta. The non-trivial results of Eq. (86) are obtained from Eqs. (75,77,82). Now, since

$$ \hat{P}_{Ri} = \frac{1}{i} \frac{\partial}{\partial R_i} \qquad (87) $$

and, in view of Eqs. (73,10),

$$ \hat{P}_{Ri}^{\dagger} = \frac{1}{i} \frac{\partial}{\partial R_i} + \frac{2}{i R_i}, \qquad (88) $$

and since

$$ \hat{P}_{Ri}^{\dagger} \hat{P}_{Ri} = - \frac{1}{R_i^2} \frac{\partial}{\partial R_i} R_i^2 \frac{\partial}{\partial R_i} = - \frac{1}{R_i} \frac{\partial^2}{\partial R_i^2} R_i, \qquad (89) $$

Eq. (76) for orthogonal coordinates can be written as [74, 76]:
$$ \hat{T} = \sum_{i=1}^{N-1} \left( -\frac{1}{2\mu_i} \frac{1}{R_i} \frac{\partial^2}{\partial R_i^2} R_i \right) + \sum_{i=1}^{N-2} \left( \frac{1}{2\mu_{N-1}R_{N-1}^2} + \frac{1}{2\mu_i R_i^2} \right) \hat{\vec{L}}_i^\dagger \hat{\vec{L}}_i + \sum_{i,j=1}^{N-2} \frac{(\hat{\vec{L}}_i \hat{\vec{L}}_j)_{BF}}{\mu_{N-1}R_{N-1}^2} - \sum_{i=1}^{N-2} \frac{(\hat{\vec{J}}_i \hat{\vec{L}}_i)_{BF}}{\mu_{N-1}R_{N-1}^2} + \frac{\hat{\vec{J}}_{BF}^2}{2\mu_{N-1}R_{N-1}^2}. \quad (90) $$

We recall that if we do not indicate the frame in which the scalar product is computed, like for the term $\hat{\vec{L}}_i^\dagger \hat{\vec{L}}_i$ in Eq. (90), it means that the scalar product can be calculated in any frame. If we now introduce the operators $\hat{L}_+ = \hat{L}_x + i\hat{L}_y$ and $\hat{L}_- = \hat{L}_x - i\hat{L}_y$ for all the angular momenta and take into account the fact that

$$ \hat{L}_{N-2z^{BF}} = \hat{J}_{z^{BF}} - \sum_{i=1}^{N-3} \hat{L}_{iz^{BF}} \quad (91) $$

(see Eq. (83)), Eq. (90) yields

$$ \begin{aligned} \hat{T} = & \sum_{i=1}^{N-1} \left( -\frac{1}{2\mu_i} \frac{1}{R_i} \frac{\partial^2}{\partial R_i^2} R_i \right) + \sum_{i=1}^{N-2} \left( \frac{1}{2\mu_{N-1}R_{N-1}^2} + \frac{1}{2\mu_i R_i^2} \right) \hat{\vec{L}}_i^\dagger \hat{\vec{L}}_i \\ & + \frac{\hat{\vec{J}}_{BF}^2}{2\mu_{N-1}R_{N-1}^2} - \frac{2\hat{J}_{z^{BF}}(\hat{J}_{z^{BF}} - \sum_{i=1}^{N-3} \hat{L}_{iz^{BF}}) + \sum_{i=1}^{N-3} \hat{L}_{iz^{BF}}(2\hat{L}_{iz^{BF}} + \sum_{j=1, j \neq i}^{N-3} \hat{L}_{jz^{BF}})}{2\mu_{N-1}R_{N-1}^2} \\ & + \frac{\sum_{i=1}^{N-3} \sum_{j=1, j>i}^{N-2} (\hat{L}_{i+BF} \hat{L}_{j-BF} + \hat{L}_{i-BF} \hat{L}_{j+BF}) - \sum_{j=1}^{N-2} (\hat{J}_{+BF} \hat{L}_{j-BF} + \hat{J}_{-BF} \hat{L}_{j+BF})}{2\mu_{N-1}R_{N-1}^2}, \end{aligned} \quad (92) $$

where, again, the order of the operators in the scalar products is strictly fixed: the key point is the fact that the BF-components of $\hat{\vec{L}}_{N-2}$ must remain on the right hand side of all scalar products. All the explicit expressions of the projections onto the BF-axes are given by Eqs. (75,77,82). Consequently, we have explicitly a general expression of the KEO whatever the number of atoms and whatever the set of $N-1$ orthogonal vectors for the family of polyspherical coordinates.

Several points must be emphasized here:

(1) except in Eq. (81), the derivation of the results is free of differential calculus. We have simply exploited the properties of the angular momenta, which are highlighted in Sec. II B 2 (iii).

(2) The physical meaning of the different parts of the KEO is quite transparent since it can
be interpreted in terms of the couplings between the momenta associated with the different vectors.

(3) As already noted, the Coriolis coupling, which is mainly included in the term $-\sum_{i=1}^{N-2} \frac{(\hat{J}\hat{L}_i)_{BF}}{\mu_{N-1} R_{N-1}^2}$, is proportional to $\frac{1}{\mu_{N-1}}$, the higher the mass $\mu_{N-1}$, the smaller the Coriolis coupling.

(4) In view of Eq. (83), it clearly appears that the expressions of $\hat{L}_{N-2+BF}$ and $\hat{L}_{N-2-BF}$ are not the usual ones (compare for instance with Eq. (75)). This is due to $\hat{R}_{N-2}$ being involved in the definition of the BF frame. However, we shall see later that their action onto an appropriate primitive basis set of spherical harmonics in terms of the BF angles is the normal one.

(5) The generalization to non-orthogonal coordinates is straightforward and given in Appendix B.

(6) Finally, if an external field is present or if the molecular system is located in a noble gas matrix or on metal surfaces, one simply has to add the kinetic energy $\hat{T}_G$ of the center of mass of the system, G, in Eq. (92). In other words, let $X_G$, $Y_G$, and $Z_G$ be the Cartesian coordinates determining the position of G in the Laboratory-Fixed (LF) frame. The KEO for the center-of-mass motion then reads

$$ \hat{T}_G = -\frac{1}{2M_T} \left( \frac{\partial^2}{\partial X_G^2} + \frac{\partial^2}{\partial Y_G^2} + \frac{\partial^2}{\partial Z_G^2} \right), \qquad (93) $$

where $M_T$ is the total mass of the molecular system and the final KEO reads $\hat{T} + \hat{T}_G$.

(7) All the BF-projections of the angular momenta appearing in Eqs. (90,92) are quasi-momenta and Eqs. (90,92) are particular cases of Eq. (23) of Section II A.

### v) Introduction of a primitive basis set of spherical harmonics.

As often when using curvilinear coordinates, several configurations of the molecular systems correspond to a singular situation for which the KEO in Eq. (92) has infinite contributions. Let us consider an example. In the first line of Eq. (92) the following terms do appear: $(-\frac{1}{2\mu_i} \frac{1}{R_i} \frac{\partial^2}{\partial R_i^2} R_i) + (\frac{1}{2\mu_{N-1} R_{N-1}^2} + \frac{1}{2\mu_i R_i^2}) \hat{L}_i^\dagger \hat{L}_i$, with $i = 1, \dots, N-3$. In view of Eq. (75), these terms can be recast as
$$ \left(-\frac{1}{2\mu_i}\frac{1}{R_i}\frac{\partial^2}{\partial R_i^2}R_i\right) + \left(\frac{1}{2\mu_{N-1}R_{N-1}^2} + \frac{1}{2\mu_i R_i^2}\right) \hat{L}_i^\dagger \hat{L}_i = \left(-\frac{1}{2\mu_i}\frac{1}{R_i}\frac{\partial^2}{\partial R_i^2}R_i\right) - \left(\frac{1}{2\mu_{N-1}R_{N-1}^2} + \frac{1}{2\mu_i R_i^2}\right) \left[\frac{1}{\sin\theta_i^{BF}}\frac{\partial}{\partial\theta_i^{BF}}\sin\theta_i^{BF}\frac{\partial}{\partial\theta_i^{BF}} + \frac{1}{\sin^2\theta_i^{BF}}\frac{\partial^2}{\partial\varphi_i^{BF}^2}\right]. \quad (94) $$

In Eq. (94), if $\theta_i^{BF}$ is equal to 0 or $\pi$ or $R_i$ is equal to zero, the terms including $\frac{1}{\sin^2\theta_i^{BF}}$ and $\frac{1}{R_i^2}$ diverge and the the volume element

$$ d\tau \propto R_i^2 \sin \theta_i^{BF} dR_i d\varphi_i^{BF} d\theta_i^{BF} \quad (95) $$

is equal to zero. In the first case ($\theta_i^{BF} = 0$ or $\pi$) $\varphi_i^{BF}$ is no longer defined, in the second case ($R_i = 0$) $\varphi_i^{BF}$ and $\theta_i^{BF}$ are no longer defined. When such situations arise, it is necessary to find an appropriate primitive basis set that removes all the singularities. Thanks to the introduction of the angular momenta in Sec. II B 2 (iii), the appropriate angular basis set appears naturally: it is the basis set of spherical harmonics that is introduced in the present section.

For orthogonal coordinates, it is also straightforward to find the radial basis set that discards the singularities for $R_i$ equal to zero: it is a basis set of Bessel functions along with an appropriate pseudo-spectral approach [115] (see Section II A in Ref. [77]). However, since they correspond to geometries that rarely occur, the topic of the radial singularities will not be addressed here.

Let us now introduce the appropriate angular basis set. The overall rotation of the molecule can be described by a set of (normalized) Wigner rotation matrix elements [109, 110]:

$$ \langle \alpha, \beta, \gamma | J, M, \Omega \rangle = D_{M\Omega}^{J*}(\alpha, \beta, \gamma), \quad (96) $$

where the star denotes complex conjugate and M is the projection of the total angular momentum onto $z^{SF}$ and $\Omega = \sum_{i=1}^{N-2} \Omega_i$ is the projection of the total angular momentum onto $z^{BF}$, $\Omega_i$ being the projection of $\vec{L}_i$ onto the same axis (note that $\Omega_{N-1} = 0$, since $\vec{R}_{N-1}$ is parallel to $z^{BF}$).

In the particular case of an isolated molecule, SF is isotropic, i.e. the orientation of $z^{SF}$ is arbitrary, and any observable must be $\alpha$-independent. The overall rotation of the molecule can thus be described by the following basis set [109, 110]:

$$ \langle (\alpha), \beta, \gamma | J, 0, \Omega \rangle = Y_J^\Omega(\beta, \gamma)(-1)^\Omega, \quad (97) $$
$$ \hat{L}_{(N-2)\pm BF} < angles | \dots > = c_{\pm}(\ell_{N-2}, \Omega_{N-2}) \times \\ < angles | \ell_1, \Omega_1, \dots, \ell_{N-3}, \Omega_{N-3}, \ell_{N-2}, \Omega_{N-2} \pm 1; J, \Omega \pm 1 > $$

$$ \hat{L}_{(N-2)zBF} < angles | \dots > = (\Omega_{N-2}) \times < angles | \dots > $$

$$ \hat{J}_{\mp BF} \hat{L}_{N-2 \pm BF} < angles | \dots; J, \Omega > = c_{\pm}(J, \Omega) c_{\pm}(\ell_{N-2}, \Omega_{N-2}) \times \\ < angles | \ell_1, \Omega_1, \dots, \ell_{N-3}, \Omega_{N-3}, \ell_{N-2}, \Omega_{N-2} \pm 1; J, \Omega \pm 1 >, \qquad (99) $$

with

$$ c_{\pm}(J, \Omega) = \sqrt{J(J+1) - \Omega(\Omega \pm 1)}. \qquad (100) $$

(We remind that the BF-projections of $\hat{J}$ have an 'anomalous' behavior [114], i.e. $\hat{J}_{\pm} | J, \Omega > = c_{\mp}(J, \Omega) | J, \Omega \mp 1 >$.) The symbol $< angles | \dots >_J$ denotes the current element of the angular basis set. Even though these results are not obvious they are not surprising as all angular momenta are calculated in the SF frame and several changes of coordinates and projections onto adequate axes is all that has been performed. Eqs. (92) and (99) lead straightforwardly, for orthogonal coordinates, to

$$ \begin{aligned} \hat{T} & | \ell_1, \Omega_1, \dots, \ell_{N-3}, \Omega_{N-3}, \ell_{N-2}; J, \Omega > \\ &= \left[ \left( -\frac{1}{2\mu_i} \frac{1}{R_i} \frac{\partial^2}{\partial R_i^2} R_i \right) + \sum_{i=1}^{N-2} \left( \frac{1}{2\mu_{N-1} R_{N-1}^2} + \frac{1}{2\mu_i R_i^2} \right) \ell_i(\ell_i + 1) \right. \\ & \qquad \left. + \frac{J(J+1)}{2\mu_{N-1} R_{N-1}^2} - \frac{2\Omega(\Omega - \sum_{i=1}^{N-3} \Omega_i) + \sum_{i=1}^{N-3} \Omega_i (2\Omega_i + \sum_{j=1, j \neq i}^{N-3} \Omega_j)}{2\mu_{N-1} R_{N-1}^2} \right] | \dots; J, \Omega > \\ & \qquad + \frac{\sum_{i=1}^{N-3} \sum_{j=1, j>i}^{N-2} (c_+(\ell_i, \Omega_i) c_-(\ell_j, \Omega_j) | \dots \Omega_i + 1 \dots \Omega_j - 1, \dots; J, \Omega >)}{2\mu_{N-1} R_{N-1}^2} \\ & \qquad + \frac{\sum_{i=1}^{N-3} \sum_{j=1, j>i}^{N-2} (c_-(\ell_i, \Omega_i) c_+(\ell_j, \Omega_j) | \dots \Omega_i - 1 \dots \Omega_j + 1, \dots; J, \Omega >)}{2\mu_{N-1} R_{N-1}^2} \\ & \qquad \left. - \frac{\sum_{j=1}^{N-2} (c_-(J, \Omega) c_-(\ell_j, \Omega_j) | \dots \Omega_j - 1, \dots; J, \Omega - 1 >)}{2\mu_{N-1} R_{N-1}^2} \right. \\ & \qquad \left. - \frac{\sum_{j=1}^{N-2} (c_+(J, \Omega) c_+(\ell_j, \Omega_j) | \dots \Omega_j + 1, \dots; J, \Omega + 1 >)}{2\mu_{N-1} R_{N-1}^2} \right]. \end{aligned} \qquad (101) $$
Consequently, the basis set of Eq. (98) along with an appropriate pseudo-spectral approach [116–118] eliminates all singularities appearing in the KEO. The same is true for non-orthogonal coordinates. In other words, equations such as Eq. (99) can be applied to evaluate the action of the BF-components appearing in Eqs. (B5,B6,B7). In addition, for orthogonal coordinates, the number of terms in the KEO always remains moderate. The only disadvantage of this formulation is that, in the grid representation, one has to deal with complicated multi-dimensional Discrete Variable Representations (DVRs) for the angles [118–120].

## 2. Matrix expression of the KEO

As above-mentioned, the conjugate momenta are only intermediates in the polyspherical approach. There is no obligation to utilize them, except in the case of the presence of singularities. If the potential prevents the system from reaching such singular geometries, one can express the KEO explicitly in differential operators, $\frac{\partial}{\partial q_i}$, where $q_i$ denote the $3N-6$ polyspherical coordinates as in the original formulation of Podolsky (Eq. (2)). As we will see, the number of terms in the KEO increases, but only simple one-dimensional DVRs are required in a numerical calculation. The present section is devoted to the derivation of the general expression of the KEO in terms of differential operators for the family of standard polyspherical coordinates.

In Section II C 1, we have established a general expression of the KEO in polyspherical coordinates. The operator is available for orthogonal vectors (Eq. 92) and for non-orthogonal vectors (Eqs. (B5,B6,B7)). As we also know the explicit expressions of the projections onto the BF-axes of all the angular momenta (Eqs. (75,77,82)), we are thus in a position to provide a general expression of the KEO in the following form, similar to Eq. (27) of Section II A, whatever the set of vectors (orthogonal or non-orthogonal):

$$ \begin{aligned} \hat{T} = & \sum_{l,m=1}^{3N-6} \frac{\hat{P}_{ql}^{\dagger} \sum_{q_l q_m} \hat{P}_{q_m}}{2} + \sum_{l=1}^{3N-6} \sum_{\alpha=x,y,z} \frac{\hat{P}_{ql}^{\dagger} \sigma_{q_l \alpha} \hat{J}_{\alpha BF} + \hat{J}_{\alpha BF} \sigma_{\alpha q_l} \hat{P}_{ql}}{2} \\ & + \sum_{\alpha=x,y,z} \sum_{\beta=x,y,z} \frac{\hat{J}_{\alpha BF} \Gamma_{\alpha \beta} \hat{J}_{\beta BF} + \hat{J}_{\beta BF} \Gamma_{\beta \alpha} \hat{J}_{\alpha BF}}{2}. \end{aligned} \quad (102) $$
In contrast with the matrix M, the matrices Σ, σ and Γ do not depend on the masses of the atoms only. They also depend on the 3N-6 internal polyspherical coordinates. In addition, they are all symmetrical, i.e. Σ$_{q_l q_m}$ = Σ$_{q_m q_l}$, Γ$_{λλ'}$ = Γ$_{λ'λ}$, and σ$_{λq_l}$ = σ$_{q_lλ}$. Ĵ$_{x}$ BF , Ĵ$_{y}$ BF , and Ĵ$_{z}$ BF are the self-adjoint BF-projections of the total angular momentum, the expressions of which are given by Eq. (77). P̂$_{q_i}$ = 1 $_{i}$∂ ∂q$_{i}$ are the (3N-6) momentum operators associated with the (3N-6) spherical coordinates { q$_{i}$ ; i = 1 , . . . , 3 N − 6 } = { ( R$_{i}$ ; i = 1 , . . . , N − 1) , ( θ BF i ; i = 1 , . . . , N − 2) , ( ϕ BF i ; i = 1 , . . . , N − 3) } . Thus, q collectively denotes the (3N-6) spherical coordinates.

(i) New convention of normalization of the wave-function:

Before providing all the matrix elements, let us point out that all the operators introduced so far are correct only for the case in which the wave functions are normalized with an Euclidean convention of normalization, i.e. with the volume element given by Eq. (73). Next, in view of Eqs. (73) and (10), it appears that P̂ θ BF i ( i = 1 , . . . , N − 2) and P̂ R$_{j}$ ( j = 1 , . . . , N − 3) are not self-adjoint, i.e. P̂ † q$_{i}$ ̸ = P̂ q$_{i}$ for these coordinates. For the sake of simplicity, we shall prefer to work with momentum operators that are Hermitian. First, we replace the coordinates θ$_{i}$ by u$_{i}$ = cos( θ BF i ); ( i = 1 , . . . , N − 2). Indeed, since

$$
du_i = - \sin \theta_i^{BF} d\theta_i^{BF}, \tag{103}
$$

after this transformation, the Euclidean volume element now reads [121]

$$
d\tau = R_{N-1}^2 \sin\beta \, dR_{N-1} \, d\alpha \, d\beta \, R_{N-2}^2 \, dR_{N-2} \, d\gamma \, du_{N-2} \prod_{i=1}^{N-3} R_i^2 \, dR_i \, d\varphi_i^{BF} \, du_i, \quad (104)
$$

so that, in view of (10), P̂ † u$_{i}$ = P̂ u$_{i}$ with ( i = 1 , . . . , N − 2). Eq. (104) can be formally written as

$$
d\tau = J'_{e} dR_{N-1} d\alpha d\beta dR_{N-2} d\gamma du_{N-2} \prod_{i=1}^{N-3} du_{i}, \quad (105)
$$

with J ′ e = sin β ∏ N − 1 i =1 R 2 $_{i}$, the absolute value of the Jacobian (see Eq. (3)) correspond- ing to the coordinate transformation from the Cartesian coordinates to the polyspherical coordinates, taking into account the fact that we use u$_{i}$ = cos( θ BF i ) instead of θ BF i for
(i = 1, ..., N − 2). The KEO now reads

$$
\begin{equation}
\begin{split}
\hat{T} = \sum_{l,m=1}^{3N-6} \frac{\hat{P}_{ql}^{\dagger} G_{qlq_m} \hat{P}_{q_m}}{2} &+ \sum_{l=1}^{3N-6} \sum_{\alpha=x,y,z} \frac{\hat{P}_{ql}^{\dagger} C_{ql\alpha} \hat{J}_{\alpha BF} + \hat{J}_{\alpha BF} C_{\alpha q_l} \hat{P}_{q_l}}{2} \\
&+ \sum_{\alpha=x,y,z} \sum_{\beta=x,y,z} \frac{\hat{J}_{\alpha BF} \Gamma_{\alpha\beta} \hat{J}_{\beta BF} + \hat{J}_{\beta BF} \Gamma_{\beta\alpha} \hat{J}_{\alpha BF}}{2}.
\end{split}
\tag{106}
\end{equation}
$$

The relation between Σ, and G, on one hand, and between σ and C, on the other hand, is
rather simple. When q_m = θ_i^{BF} in Eq. (106), one has to multiply the corresponding matrix
element by − sin θ_i^{BF}. Let us give two examples:

$$
\begin{align}
G_{u_1 R_2} &= -\sin \theta_1^{BF} \times \Sigma_{\theta_1^{BF} R_2}, \nonumber \\
G_{u_1 u_2} &= \sin \theta_1^{BF} \times \Sigma_{\theta_1^{BF} \theta_2^{BF}} \times \sin \theta_2^{BF}. \tag{107}
\end{align}
$$

Second, we also change the convention of normalization of the wave-function, i.e. we now
use the following volume element in normalizing wave-functions:

$$
d\tau' = \sin\beta \, dR_{N-1} \, d\alpha \, d\beta \, dR_{N-2} \, d\gamma \, du_{N-2} \prod_{i=1}^{N-3} dR_i \, d\varphi_i^{BF} \, du_i. \quad (108)
$$

Let us also take

$$
d\tau' = \rho dR_{N-1} d\alpha d\beta dR_{N-2} d\gamma du_{N-2} \prod_{i=1}^{N-3} dR_i d\varphi_i^{BF} du_i, \quad (109)
$$

with $\rho = \sin \beta$. This new convention of normalization entails that the operators $\hat{P}_{R_i}$ ($i = 1, \dots, N-1$) are now self-adjoint (see Eq. (18)). However, the KEO, $\hat{T}$, has to be replaced by

$$
\hat{T}' = J_e^{\frac{1}{2}} \rho^{-\frac{1}{2}} \hat{T} \rho^{\frac{1}{2}} J_e^{-\frac{1}{2}} = \prod_i^{N-1} (R_i) \hat{T} \prod_i^{N-1} (R_i^{-1}), \quad (110)
$$

so as to preserve the normalization of the wave-functions (see Eq. (17)). The explicit
calculation yields:
$$ \begin{aligned} \hat{T}' = & \sum_{l,m=1}^{3N-6} \frac{\hat{P}_{ql} G_{qlqm} \hat{P}_{qm}}{2} + \sum_{l=1}^{3N-6} \sum_{\alpha=x,y,z} \frac{\hat{P}_{ql} C_{ql\alpha} \hat{J}_{\alpha BF} + \hat{J}_{\alpha BF} C_{\alpha q\alpha} \hat{P}_{ql}}{2} \\ & + \sum_{\alpha=x,y,z} \sum_{\beta=x,y,z} \frac{\hat{J}_{\alpha BF} \Gamma_{\alpha\beta} \hat{J}_{\beta BF} + \hat{J}_{\beta BF} \Gamma_{\beta\alpha} \hat{J}_{\alpha BF}}{2} + V_{extra}, \end{aligned} \quad (111) $$

with $V_{extra}$ a purely multiplicative operator called 'extra potential term' (see Eqs. (19,20)), which, in our particular case, reads

$$ V_{extra} = - \sum_{i \neq j}^{N-1} \frac{G_{R_i R_j}}{2 R_i R_j}. \quad (112) $$

Interestingly enough, in the case of orthogonal coordinates $V_{extra} = 0$. For non-orthogonal coordinates, $V_{extra} \neq 0$ (see Appendix B: in Eq. (B7), there is a coupling between $R_i$ and $R_j$. Consequently, $G_{R_i R_j} \neq 0$).

Inserting Eqs. (75,77,82) into Eq. (92) provides all the matrix elements appearing in **G**, **C** and **Γ**, which are given in Appendix C. We have used the notation **G** for the vibrational part of the KEO to keep the notation of Wilson *et al* [98] but it should be emphasized that the **G**-matrix elements in Appendix C are different from those in Ref. [98] since they correspond to a different family of coordinates.

Let us consider only one term in Eq. (111), the first appearing for $G_{\varphi_i \varphi_j}$ in Eq. (C4):

$$ G_{\varphi_i \varphi_j} = M_{ij} \frac{\cos(\varphi_i - \varphi_j)}{R_i R_j \sin \theta_i \sin \theta_j}, \quad (113) $$

Hence, the following terms crop up in the KEO

$$ \begin{aligned} & M_{ij} \left( \hat{P}_{\varphi_i} \frac{\cos(\varphi_i - \varphi_j)}{R_i R_j \sin \theta_i \sin \theta_j} \hat{P}_{\varphi_j} + \hat{P}_{\varphi_j} \frac{\cos(\varphi_i - \varphi_j)}{R_i R_j \sin \theta_i \sin \theta_j} \hat{P}_{\varphi_i} \right) \\ & = M_{ij} \left( \hat{P}_{\varphi_i} \frac{\cos \varphi_i \cos \varphi_j}{R_i R_j \sin \theta_i \sin \theta_j} \hat{P}_{\varphi_j} + \hat{P}_{\varphi_j} \frac{\cos \varphi_i \cos \varphi_j}{R_i R_j \sin \theta_i \sin \theta_j} \hat{P}_{\varphi_i} \right) \\ & + M_{ij} \left( \hat{P}_{\varphi_i} \frac{\sin \varphi_i \sin \varphi_j}{R_i R_j \sin \theta_i \sin \theta_j} \hat{P}_{\varphi_j} + \hat{P}_{\varphi_j} \frac{\sin \varphi_i \sin \varphi_j}{R_i R_j \sin \theta_i \sin \theta_j} \hat{P}_{\varphi_i} \right). \end{aligned} \quad (114) $$

Two points must be addressed here: (1) as long as the singular geometries (for instance for $\theta_i^{BF} = 0$ or $\pi$) are not physically accessible, the action of all the terms in the KEO
(such as those in Eq. (114)) can be evaluated with a simple one-dimensional primitive basis set (DVR or Finite Basis Representation (FBR)) for each degree of freedom. (2) The last two lines of Eq. (114) clearly show that the terms are products of one-dimensional functions. This property is true for all the terms in the KEO in Eq. (111) and can be very useful for many numerical methods used in dynamics to solve the Schrödinger equation.

(ii) Connection between polyspherical coordinates and normal modes:

The previous formulation allows a simple correspondence between normal and polyspherical coordinates. This may turn out to be very helpful since it permits to analyze the results of the dynamics in the normal mode picture. The method is as follows [122]: we define curvilinear normal modes from a zero order harmonic Hamiltonian, $\hat{H}^o$, expressed as

$$ \hat{H}^o = \frac{1}{2} \sum_{l,m=1}^{3N-6} (q_l - q_l^{eq}) F_{q_l q_m} (q_m - q_m^{eq}) + \hat{P}_{q_l} G_{q_l q_m}^o \hat{P}_{q_m}, \quad (115) $$

where $G^o$ represents the $G$ matrix of Eq. (111) but evaluated at the equilibrium geometry $q_{eq}$. Here, $q_l$ and $q_m$ denote the $3N-6$ polyspherical coordinates and the $F$ matrix corresponds to the harmonic approximation for the potential (including the extra potential term): $F_{q_l q_m} = \partial^2 V'/\partial q_l \partial q_m|_{q_{eq}}$; $V'$ is the potential $V$ plus $V_{extra}$. One can then proceed with the Wilson $G$ matrix formulation [98, 123], and define curvilinear normal modes $\{Q_\alpha\}$ in terms of the polyspherical coordinates. In this, one diagonalizes the matrix $FG^o$

$$ F G^o L = L \omega^2, \quad (116) $$

where $\omega^2$ denotes the diagonal eigenvalue matrix with diagonal elements $\omega_\alpha^2$ and $L$ is the eigenvector matrix subject to the normalization

$$ L^T G^o L = 1. \quad (117) $$

$A^T$ denotes the transpose of $A$, $A$ being a matrix or a vector. The mass and frequency weighted normal coordinates $Q_\alpha$ are related to the polyspherical coordinates as

$$ Q_{\alpha} = \omega_{\alpha}^{1/2} \sum_{l=1}^{3N-6} (q_l - q_l^{eq}) L_{q_l \alpha} . \quad (118) $$

These dimensionless normal coordinates let $\hat{H}^o$ become separable,

$$ \hat{H}^o = \sum_{\alpha=1}^{3N-6} \hat{h}_{\alpha}^o, \quad (119) $$
and the one-dimensional operators $\hat{h}_{\alpha}^{o}$ read

$$ \hat{h}_{\alpha}^{o} = \frac{\omega_{\alpha}}{2} (Q_{\alpha}^{2} - \frac{\partial^{2}}{\partial Q_{\alpha}^{2}}). \qquad (120) $$

Note that these operators depend on all $3N-6$ polyspherical coordinates because of Eq. (118).

Finally, it should be emphasized that the curvilinear normal mode frequencies, $\omega_{\alpha}$, are identical to the usual normal mode frequencies.

## D. Separation into subsystems

The general approach outlined in sections IIB and IIC can be generalized to derive kinetic energy operators in terms of other kinds of coordinates than the standard polyspherical coordinates. The latter are defined with respect to two vectors only: $\vec{R}_{N-1}$ and $\vec{R}_{N-2}$. For large systems, this can introduce artificial couplings between motions that are physically decoupled, e. g. between atoms that are located in very different parts of the molecule. A simple solution is to integrate atoms into subsystems and to define for each subsystem an intermediate BF frame. In each subsystem, the angles are defined with respect to the intermediate frame, which decouples the motions between atoms belonging to different subsystems.

### (i) Vector parametrization and separation into subsystems:

Let an $N$-atom system be separated into $n$ subsystems. Let $N^{(j)}$ be the number of atoms of the $j$-th subsystem, ($j = 1, ..., n$). The system thus obtained is parametrized by $\sum_{j=1}^{n}(N^{(j)} - 1) + n - 1 = N - 1$ vectors. The vectors are divided into two groups: first, the $\sum_{j=1}^{n}(N^{(j)} - 1)$ vectors $\vec{R}_{k}^{(j)}$ with ($k = 1, ..., N^{(j)} - 1$) parametrizing the subsystems and, second, the $n-1$ 'remaining' vectors $\vec{R}_{i}$ with ($i = 1, ..., n - 1$) that parametrize the relative positions of the subsystems. Let $\hat{J}^{(j)}$ be the total angular momentum of the $j$-th subsystem. This definition entails that

$$ \hat{J}^{(j)} = \sum_{k=1,...,N^{(j)}-1} \hat{L}_{k}^{(j)}, \qquad (121) $$

$\hat{L}_{k}^{(j)}$ being the angular momentum associated with $\vec{R}_{k}^{(j)}$.
FIG. 8: Definition of the eight Jacobi vectors for the water trimer. Here, $G_j$ is the center of mass of the $j$-th monomer. $\vec{R}_1^{(j)}$ is the vector joining the two H atoms, $\vec{R}_2^{(j)}$ is the vector joining the center of mass of $H_2$ to the O atom. $\vec{R}_2$ joins $G_3$ to $G_2$ and $\vec{R}_2$ joins the center of mass of $(G_3G_2)$ to $G_1$.

<figure><img src="image_9.png" /></figure>

To illustrate our ideas, let us consider an example: the water trimer using eight Jacobi vectors, which are shown in Figure 8. In the present case, $n = 3$, $N^{(j)} - 1 = 2$ with $j = 1, 2, 3$, and each subsystem is described by two Jacobi vectors. $\vec{R}_1$ and $\vec{R}_2$ are the two 'remaining' vectors.

For the definition of the BF frame, we choose $\vec{R}_{n-1}$ as the vector that defines the $z^{BF}$ axis and, if $n > 2$, $\vec{R}_{n-2}$ as the vector that defines the $((x^{BF}, z^{BF}); x^{BF} > 0)$ half-plane (as in Figure 6). If $n = 2$, i.e. for a separation into two subsystems, or if, for any physical reasons, it is not judicious to choose one vector $\vec{R}_j$ to define the $((x^{BF}, z^{BF}); x^{BF} > 0)$ half-plane, this half-plane must be defined by a vector belonging to one of the subsystems. The latter case is slightly more complicated and is addressed in Section III C 1 (i).

In the standard polyspherical parametrization, the system is described by $3N-3$ coordinates: (1) the three Euler angles $\alpha, \beta$, and $\gamma$;
(2) N-1 vector lengths {$R_i, (i = 1, ..., n-1); R_k^{(j)}, (k = 1, ..., N^{(j)} - 1; j = 1, ..., n)$};

(3) N-2 planar angles {$\theta_i^{BF}, (i = 1, ..., n-2); \theta_k^{(j,BF)}, (k = 1, ..., N^{(j)} - 1; j = 1, ..., n)$};

(4) N-3 dihedral angles {$\varphi_i^{BF}, (i = 1, ..., n-3); \varphi_k^{(j,BF)}, (k = 1, ..., N^{(j)} - 1; j = 1, ..., n)$}.

(ii) New polyspherical parametrization:

Let us now perform a new coordinate transformation from the 3N-3 standard polyspherical coordinates to the new 3N-3 coordinates defined as follows.

Let BF<sup>(j)</sup> be an intermediate frame associated with the j-th subsystem such as $\vec{e}_{zBF(j)}$ is parallel to the vector $\vec{R}_{N(j)-1}^{(j)}$ and $\vec{R}_{N(j)-2}^{(j)}, \vec{e}_{xBF(j)}$ and $\vec{e}_{zBF(j)}$ lie in the same half-plane (as in Fig. 6). The orientation of the BF<sup>(j)</sup> frame with respect to the BF frame is determined by three Euler angles: $\alpha^{(j)}, \beta^{(j)}, \gamma^{(j)}$. With our definition of the BF<sup>(j)</sup> frame, it is clear that the first two Euler angles simply are the two spherical angles of $\vec{R}_{N(j)-1}^{(j)}$ in the BF frame, $\alpha^{(j)} = \varphi_{N(j)-1}^{(j,BF)}$ and $\beta^{(j)} = \theta_{N(j)-1}^{(j,BF)}$. The third Euler angle is given by $\gamma^{(j)} = \varphi_{N(j)-2}^{(j,E_2^{(j)})}$, where the $E_2^{(j)}$ frame is the frame resulting from the first two Euler rotations $\alpha^{(j)}, \beta^{(j)}$ only (starting from the BF frame). Next, the $N^{(j)} - 1$ vectors are parametrized by their spherical angles in the BF<sup>(j)</sup> frame: $\theta_k^{BF(j)}$ and $\varphi_k^{BF(j)}$.

The new polyspherical parametrization thus obtained consists of 3N-3 coordinates that are:

(1) the three Euler angles $\alpha, \beta$, and $\gamma$;

(2) N-1 vector lengths {$R_i, (i = 1, ..., n-1); R_k^{(j)}, (k = 1, ..., N^{(j)} - 1; j = 1, ..., n)$};

(3) n-2 planar angles {$\theta_i^{BF}, (i = 1, ..., n-2)$};

(4) n-3 dihedral angles {$\varphi_i^{BF}, (i = 1, ..., n-3)$};

(5) and for each subsystem : three Euler angles $\alpha^{(j)}, \beta^{(j)}, \gamma^{(j)}$;

$N^{(j)} - 2$ planar angles $\theta_k^{BF(j)} (k = 1, ..., N^{(j)} - 2)$;

$N^{(j)} - 3$ dihedral angles $\varphi_k^{BF(j)} (k = 1, ..., N^{(j)} - 3)$.

(iii) Expression of the KEO:

Using the properties of the angular momenta highlighted in Sec. II B 2 (iii) allows the expression of the KEO to be derived without resorting to differential calculus.

Indeed, suffice it to note for our present purpose that:

(1) since the vectors in the different subsystems are not involved in the definition of the
BF frame, the expression of the BF-components of their angular momenta in terms of the BF coordinates is identical to the expression of the SF-projections of any arbitrary angular momentum in terms of SF coordinates.

(2) For each subsystem, the coordinate transformation from $\{\theta_k^{(j,BF)}, (k=1, \dots, N^{(j)}-1); \varphi_k^{(j,BF)}, (k=1, \dots, N^{(j)}-1)\}$ to $\{\alpha^{(j)}, \beta^{(j)}, \gamma^{(j)}; \theta_k^{BF(j)}(k=1, \dots, N^{(j)}-2); \varphi_k^{BF(j)}(k=1, \dots, N^{(j)}-3)\}$ only affects the coordinates of the j-th subsystem and not the other coordinates.

From this, we can deduce that the previous coordinate transformation for the j-th subsystem alone is mathematically identical to the coordinate transformation from the SF to the BF coordinates thoroughly studied in Sec. IIB 2 and IIC 1. Therefore, the BF(j)-components of the angular momenta can be readily obtained from the formulae derived in Sec. IIC 1.

In view of Eq. (77)

$$ \begin{bmatrix} \hat{J}_{x^{BF(j)}} \\ \hat{J}_{y^{BF(j)}} \\ \hat{J}_{z^{BF(j)}} \end{bmatrix} = \begin{bmatrix} -\frac{\cos \gamma^{(j)}}{\sin \beta^{(j)}} & \sin \gamma^{(j)} & \cot \beta^{(j)} \cot \gamma^{(j)} \\ \frac{\sin \gamma^{(j)}}{\sin \beta^{(j)}} & \cos \gamma^{(j)} & -\cot \beta^{(j)} \sin \gamma^{(j)} \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \frac{1}{i} \frac{\partial}{\partial \alpha^{(j)}} \\ \frac{1}{i} \frac{\partial}{\partial \beta^{(j)}} \\ \frac{1}{i} \frac{\partial}{\partial \gamma^{(j)}} \end{bmatrix}, \quad (122) $$

or still in view of (75) for $k < N^{(j)} - 2$

$$ \begin{bmatrix} \hat{L}_k x^{BF(j)} \\ \hat{L}_k y^{BF(j)} \\ \hat{L}_k z^{BF(j)} \end{bmatrix} = \begin{bmatrix} -\sin \varphi_k^{BF(j)} & -\cos \varphi_k^{BF(j)} \cot \theta_k^{BF(j)} \\ \cos \varphi_k^{BF(j)} & -\sin \varphi_k^{BF(j)} \cot \theta_k^{BF(j)} \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \frac{1}{i} \frac{\partial}{\partial \theta_k^{BF(j)}} \\ \frac{1}{i} \frac{\partial}{\partial \varphi_k^{BF(j)}} \end{bmatrix}, \quad (123) $$

and, from (83)
$$
\begin{bmatrix}
\hat{L}_{(N^{(j)}-2) x^{BF^{(j)}}} \\
\hat{L}_{(N^{(j)}-2) y^{BF^{(j)}}} \\
\hat{L}_{(N^{(j)}-2) z^{BF^{(j)}}}
\end{bmatrix}
=
\begin{bmatrix}
-\cot \theta_{N^{(j)}-2}^{BF^{(j)}}(\hat{J}_{z^{BF^{(j)}}}^{(j)} - \sum_{k=1}^{N^{(j)}-3} \hat{L}_{k z^{BF^{(j)}}}) \\
-i \frac{\partial}{\partial \theta_{N^{(j)}-2}^{BF^{(j)}}} \\
\hat{J}_{z^{BF^{(j)}}}^{(j)} - \sum_{k=1}^{N^{(j)}-3} \hat{L}_{k z^{BF^{(j)}}}
\end{bmatrix}
. \quad (124)
$$

In order to obtain the KEO in terms of the new set of polyspherical coordinates, let us first
confine ourselves to studying the case of orthogonal vectors. Eq. (76) can be recast as

$$
2\hat{T} = \sum_{i=1}^{n-1} \frac{\hat{P}_{R_i}^{\dagger} \hat{P}_{R_i}}{\mu_i} + \sum_{i=1}^{n-2} \frac{\hat{L}_i^{\dagger} \hat{L}_i}{\mu_i R_i^2} + 2 \sum_{j=1}^{n} \hat{T}^{(j)} \\
+ \frac{(\hat{J}^{\dagger} - (\sum_{i=1}^{n-2} \hat{L}_i^{\dagger} + \sum_{j=1}^{n} \hat{J}^{(j)})^{\dagger} \cdot (\hat{J} - (\sum_{i=1}^{n-2} \hat{L}_i + \sum_{j=1}^{n} \hat{J}^{(j)}))}{\mu_{n-1} R_{n-1}^2}, \quad (125)
$$

with

$$
2\hat{T}^{(j)} = \sum_{k=1}^{N^{(j)}-1} \frac{\hat{P}_{R_k^{(j)}}^{\dagger} \hat{P}_{R_k^{(j)}}}{\mu_k^{(j)}} + \sum_{k=1}^{N^{(j)}-1} \frac{\hat{L}_k^{(j)\dagger} \hat{L}_k^{(j)}}{\mu_k^{(j)} R_k^{(j)2}}, \quad (126)
$$

where $\mu_i$ denotes the reduced mass associated with $\vec{R}_i$ and $\mu_k^{(j)}$, the reduced mass associated
with $\vec{R}_k^{(j)}$. Taking into account that $\hat{L}_{n-2}$ has a particular behavior, we obtain from Eq.
(125) (the key point is the fact that the BF-components of $\hat{L}_{n-2}$ must remain on the right
hand side of all the scalar products):

$$
\begin{align}
\hat{T} ={}& \sum_{i=1}^{n-1} \left( -\frac{1}{2\mu_i} \frac{1}{R_i} \frac{\partial^2}{\partial R_i^2} R_i \right) + \sum_{i=1}^{n-2} \left( \frac{1}{2\mu_{n-1} R_{n-1}^2} + \frac{1}{2\mu_i R_i^2} \right) \hat{L}_i^{\dagger} \hat{L}_i + \sum_{j=1}^{n} \hat{T}^{(j)} \nonumber \\
& + \sum_{i,j=1}^{n-2} \frac{(\hat{L}_i \hat{L}_j)_{BF}}{\mu_{n-1} R_{n-1}^2} - \sum_{i=1}^{n-2} \frac{(\hat{J} \hat{L}_i)_{BF}}{\mu_{n-1} R_{n-1}^2} + \frac{\hat{J}_{BF}^2}{2\mu_{n-1} R_{n-1}^2} \nonumber \\
& + \sum_{i=1, \dots, n-2}^{j=1, \dots, n} \frac{(\hat{J}^{(j)} \hat{L}_i)_{BF}}{\mu_{n-1} R_{n-1}^2} + \sum_{l,m=1}^{n} \frac{(\hat{J}^{(l)} \hat{J}^{(m)})_{BF}}{\mu_{n-1} R_{n-1}^2} - \sum_{j=1}^{n} \frac{(\hat{J} \hat{J}^{(j)})_{BF}}{\mu_{n-1} R_{n-1}^2}. \tag{127}
\end{align}
$$

Eq. (127) is strictly equivalent to Eq. (90) except that (1) N has been replaced by n,
(2) $\hat{J}$ has been replaced by $\hat{J} - \sum_{j=1,...,n} \hat{J}^{(j)}$, and (3) $\sum_{j=1}^{n} \hat{T}^{(j)}$ has been added to the
KEO. The BF-components of $\hat{J}$ and $\hat{L}_i$, ($i = 1, \dots, n-1$) are given by Eq. (77) and Eq. (75), respectively. In the last two lines of Eq. (127), only the BF-projections of $\hat{J}^{(j)}$ with ($j = 1, \dots, n$) appear. They can be obtained from Eqs. (122,69,70) (they are similar to the SF-projections of a total angular momentum [109]):

$$ \begin{bmatrix} \hat{J}_{x^{BF}}^{(j)} \\ \hat{J}_{y^{BF}}^{(j)} \\ \hat{J}_{z^{BF}}^{(j)} \end{bmatrix} = \begin{bmatrix} -\cos \alpha^{(j)} \cot \beta^{(j)} & -\sin \alpha^{(j)} & \frac{\cos \alpha^{(j)}}{\sin \beta^{(j)}} \\ -\sin \alpha^{(j)} \cot \beta^{(j)} & \cos \alpha^{(j)} & \frac{\sin \alpha^{(j)}}{\sin \beta^{(j)}} \\ 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} \frac{1}{i} \frac{\partial}{\partial \alpha^{(j)}} \\ \frac{1}{i} \frac{\partial}{\partial \beta^{(j)}} \\ \frac{1}{i} \frac{\partial}{\partial \gamma^{(j)}} \end{bmatrix} . \quad (128) $$

Finally, Eq. (92) or Eq. (76) can be applied to express $\hat{T}^{(j)}$:

$$ \begin{aligned} \hat{T}^{(j)} = & \sum_{k=1}^{N^{(j)}-1} \left( -\frac{1}{2\mu_k^{(j)}} \frac{1}{R_k^{(j)}} \frac{\partial^2}{\partial R_k^{(j)^2}} R_k^{(j)} \right) + \sum_{k=1}^{N^{(j)}-2} \left( \frac{1}{2\mu_{N^{(j)}-1}^{(j)} (R_{N^{(j)}-1}^{(j)})^2} + \frac{1}{2\mu_k^{(j)} R_k^{(j)^2}} \right) \hat{L}_k^{(j)\dagger} \hat{L}_k^{(j)} \\ & - \frac{2\hat{J}_{z^{BF}(j)}^{(j)} (\hat{J}_{z^{BF}(j)}^{(j)} - \sum_{k=1}^{N^{(j)}-3} \hat{L}_{kz^{BF}(j)}^{(j)}) + \sum_{k=1}^{N^{(j)}-3} \hat{L}_{kz^{BF}(j)}^{(j)} (2\hat{L}_{kz^{BF}(j)}^{(j)} + \sum_{k'=1, k' \neq k}^{N^{(j)}-3} \hat{L}_{k'z^{BF}(j)}^{(j)})}{2\mu_{N^{(j)}-1}^{(j)} (R_{N^{(j)}-1}^{(j)})^2} \\ & + \frac{\sum_{k=1}^{N^{(j)}-3} \sum_{k'=1, k'>k}^{N^{(j)}-2} (\hat{L}_{k+^{BF}(j)}^{(j)} \hat{L}_{k'-^{BF}(j)}^{(j)} + \hat{L}_{k-^{BF}(j)}^{(j)} \hat{L}_{k'+^{BF}(j)}^{(j)})}{2\mu_{N^{(j)}-1}^{(j)} (R_{N^{(j)}-1}^{(j)})^2} \\ & - \frac{\sum_{k=1}^{N^{(j)}-2} (\hat{J}_{+^{BF}(j)}^{(j)} \hat{L}_{k-^{BF}(j)}^{(j)} + \hat{J}_{-^{BF}(j)}^{(j)} \hat{L}_{k+^{BF}(j)}^{(j)})}{2\mu_{N^{(j)}-1}^{(j)} (R_{N^{(j)}-1}^{(j)})^2} + \frac{\hat{J}_{BF(j)}^{(j)}{}^2}{2\mu_{N^{(j)}-1}^{(j)} (R_{N^{(j)}-1}^{(j)})^2}. \end{aligned} \quad (129) $$

To be more specific, let us consider various parametrizations of the water trimer. In the case of the set of orthogonal vectors depicted in Figure 8, the matrix $M$ of Eq. (54) becomes
$$ M = \begin{bmatrix} \frac{2}{m_H} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{2m_H+m_O}{2m_H m_O} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{2}{m_H} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{2m_H+m_O}{2m_H m_O} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{2}{m_H} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{2m_H+m_O}{2m_H m_O} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{m_H+m_O} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{3}{4m_H+2m_O} \end{bmatrix} . \quad (130) $$

In the matrix of Eq. (130), eight reduced masses appear. They are associated with $\vec{R}_1^{(1)}$, $\vec{R}_2^{(1)}$, $\vec{R}_1^{(2)}$, $\vec{R}_2^{(2)}$, $\vec{R}_1^{(3)}$, $\vec{R}_2^{(3)}$, $\vec{R}_1$, $\vec{R}_2$ (see Figure 8). The KEO is given by Eqs. (127) and (129) with $n=3$, $N^{(j)} = 2$, ($j=1, 2, 3$).

Now, what happens if non-orthogonal coordinates are introduced?

Several cases must be distinguished depending on the structure of the matrix M from Eq. (54).

(1) If the subsystems are described by non-orthogonal vectors and if there is no coupling in the matrix M between one subsystem and the rest of the molecular system, then Eq. (111) or Eqs. (B5,B6,B7) can be used for $\hat{T}^{(j)}$. In Eqs. (111) or (B5,B6,B7) N must be replaced by $N^{(j)}$ and $\hat{J}$ by $\hat{J}^{(j)}$.

Let us give a new example: the mixed Jacobi/valence parametrization for the water trimer depicted in Figure 9. As far as the matrix M is concerned, one obtains
<figure><img src="image_10.png" /><figcaption>FIG. 9: Definition of the second set of vectors for the water trimer. $\vec{R}_1^{(j)}$ and $\vec{R}_2^{(j)}$ are the two valence vectors joining one H atom to the O atom.</figcaption></figure>

$$ M = \begin{bmatrix} \frac{m_H+m_O}{m_H m_O} & \frac{1}{m_O} & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{m_O} & \frac{m_H+m_O}{m_H m_O} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{m_H+m_O}{m_H m_O} & \frac{1}{m_O} & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{m_O} & \frac{1}{m_H} + \frac{1}{m_O} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{m_H+m_O}{m_H m_O} & \frac{1}{m_O} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{m_O} & \frac{m_H+m_O}{m_H m_O} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{2m_H+m_O} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{3}{4m_H+2m_O} \end{bmatrix} . \quad (131) $$

The KEO is given by Eq. (127), $\hat{T}^{(j)}$ is given by Eq. (111) or Eqs. (B5,B6,B7) with $n=3$, $N^{(j)} = 2$, ($j=1, 2, 3$). In Eqs. (111) or (B5,B6,B7) $N$ must be replaced by $N^{(j)} = 2$ and $\hat{J}$ by $\hat{J}^{(j)}$.

(2) If the $n-1$ vectors $\vec{R}_i$ are non-orthogonal vectors and if there is no coupling in the matrix $M$ between the $n-1$ vectors $\vec{R}_i$ and all the subsystems, then Eq. (111) or Eqs.
<figure><img src="image_11.png" /><figcaption>FIG. 10: Definition of the third set of vectors for the water trimer. Here, $G_j$ is the center of mass of the j-th monomer.</figcaption></figure>

(B5,B6,B7) can be used for $\hat{T} - \sum_{i=1,...,n} \hat{T}^{(j)}$ with two changes : $N$ must be replaced by $n$ and $\hat{J}$ by $\hat{J} - \sum_{j=1,...,n} \hat{J}^{(j)}$.

Let us consider the third parametrization of the water trimer depicted in Figure 10. The new matrix $M$ reads (note its block-diagonal form)

$$ M = \begin{bmatrix} \frac{m_H+m_O}{m_H m_O} & \frac{1}{m_O} & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{m_O} & \frac{m_H+m_O}{m_H m_O} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{m_H+m_O}{m_H m_O} & \frac{1}{m_O} & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{m_O} & \frac{m_H+m_O}{m_H m_O} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{m_H+m_O}{m_H m_O} & \frac{1}{m_O} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{m_O} & \frac{m_H+m_O}{m_H m_O} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{2m_H+m_O} & -\frac{1}{2m_H+m_O} \\ 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{2m_H+m_O} & \frac{2}{2m_H+m_O} \end{bmatrix} . \quad (132) $$

The KEO is given by Eq. (111) or Eqs. (B5,B6,B7) for $\hat{T} - \sum_{j=1,...,n} \hat{T}^{(j)}$ (with $N^{(j)}$ instead of $N$ and $\hat{J}^{(j)}$ instead of $\hat{J}$) and $\hat{T}^{(j)}$ is given by Eq. (111) or Eqs. (B5,B6,B7) with
<figure><img src="image_12.png" /><figcaption>FIG. 11: Definition of the fourth set of vectors for the water trimer. Here, only valence vectors are used.</figcaption></figure>

$n=3, N^{(j)} = 2, (j = 1, 2, 3).$

(3) The situation becomes more complex otherwise.

Let the water trimer be parametrized by the valence vectors depicted in Figure 11. Here, the matrix $M$ becomes (note that the matrix is no longer block-diagonal)

$$ M = \begin{bmatrix} \frac{m_H+m_O}{m_H m_O} & \frac{1}{m_O} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{m_O} \\ \frac{1}{m_O} & \frac{m_H+m_O}{m_H m_O} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{m_O} \\ 0 & 0 & \frac{m_H+m_O}{m_H m_O} & \frac{1}{m_O} & 0 & 0 & -\frac{1}{m_O} & \frac{1}{m_O} \\ 0 & 0 & \frac{1}{m_O} & \frac{m_H+m_O}{m_H m_O} & 0 & 0 & -\frac{1}{m_O} & \frac{1}{m_O} \\ 0 & 0 & 0 & 0 & \frac{m_H+m_O}{m_H m_O} & \frac{1}{m_O} & \frac{1}{m_O} & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{m_O} & \frac{m_H+m_O}{m_H m_O} & \frac{1}{m_O} & 0 \\ 0 & 0 & -\frac{1}{m_O} & -\frac{1}{m_O} & \frac{1}{m_O} & \frac{1}{m_O} & \frac{2}{m_O} & -\frac{1}{m_O} \\ -\frac{1}{m_O} & -\frac{1}{m_O} & \frac{1}{m_O} & \frac{1}{m_O} & 0 & 0 & -\frac{1}{m_O} & \frac{2}{m_O} \end{bmatrix}. \quad (133) $$

The situation becomes more intricate since scalar products such as $\frac{1}{m_O} \vec{L}_2 \hat{L}_1^{(1)}$ appear in the
KEO (see the second line of Eq. (176) in Section III C 1 (ii) and Eq. (133) to understand the origin of this term). These scalar products raise a new technical problem as, on one hand, the expression of the BF-projections of $\hat{L}_1^{(1)}$ are very complicated and, on the other hand, the $BF^{(1)}$-projections of $\vec{L}_2$ are complicated. Obviously, one of these two frames has to be chosen to calculate the scalar product. If the scalar product is calculated in the BF frame for instance, Eq. (75) provides the BF-projections of $\vec{L}_2$. Unfortunately, the BF-projections of $\hat{L}_1^{(1)}$ are more complicated. They can be obtained from (123) combined with Eqs. (69) and (70): the result is more complicated than in the case of standard polyspherical coordinates because of the presence of the three Euler angles $\alpha^{(1)}, \beta^{(1)}, \gamma^{(1)}$. One example is given in Section III C 1 (ii).

However, it should be emphasized that, whatever the situation, all the terms can be calculated in a systematic way. Therefore, it is possible, in principle, to derive a general expression of the KEO whatever the set of vectors, whatever the number of atoms, *and* whatever the choice of the subsystems. It is also possible to separate each subsystem into a new group of subsystems: see Sec. IV C 3 (i) for an example.

E. Checking the correctness of the operators

It is not enough to derive the kinetic energy operator, it is also important to check its correctness. For this purpose, two strategies have been adopted. First, for several particular cases, we have compared the operators obtained from Eqs. (111,C4,C5,C6) with those obtained using the expression of Podolsky from Eq. (2) derived with Mathematica or Maple [67]. We have checked, for instance, that Eq. (C4) gives the same matrix elements as those previously calculated with Maple and published in Ref. [124] for an (AB)CD(EF) system in valence coordinates (Figure 3). We have also verified that Eqs. (111,C4,C5,C6) can reproduce the KEOs (including the rotation and the Coriolis coupling) published in Ref. [59]. In Ref. [87], the KEOs for tetra-atomic systems parametrized by Jacobi, valence, and satellite coordinates were given. They were first calculated 'by hand' and later checked using Mathematica [125].

Second, for several systems, the correctness of the derivation of the KEO was checked
by comparing it with numerical results provided by the program TNUM of Lauvergnat and Nauts [68]. TNUM computes **G**, **Γ**, and **C** of Eqs. (111,C4) numerically. We have verified that the numerical values of all the functions $G_{ij}(\mathbf{q})$ at several non-symmetrical grid points **q** agree with those provided by the program TNUM for systems such as $C_2H_4$ [126], $HCF_3$ in valence coordinates [127] or the Zundel cation, $H_5O_2^+$, in Jacobi coordinates [128] or in mixed Jacobi/valence coordinates (in all cases for $J = 0$). The latter test does not only guarantee that the KEO is correct but also that the operator is correctly implemented in the code used to solve the Schrödinger equation, in this case the Heidelberg MCTDH package [129].
### III. EXPLICIT EXPRESSIONS FOR PARTICULAR CASES AND APPLICATIONS

In Part II, we have described the configuration of an *N*-particle system, in the SF frame, by *N*-1 relative position vectors. We have considered only isolated molecular systems but it is straightforward to generalize the expressions presented above to molecules that are not isolated. We have presented several vector parametrizations of an *N*-particle system, for instance, parametrizations based on Jacobi or valence vectors (see, e.g., Figures 1,2,3). We have expressed the kinetic energy in terms of quasi-momenta (see Section II A (iii) for a definition of the quasi-momenta): the angular momenta associated with the vectors. We have explained how to switch from one set of vectors to another one and how to calculate the matrix *M* that appears in the kinetic energy. In Section II B 2 (ii), we have shown that the quantization of the kinetic energy is straightforward. In order to separate the overall rotation from the internal vibration, we have introduced, in Section II B 2 (iii), the Body-Fixed (BF) frame of the molecular system and highlighted some general properties regarding the projections of the angular momenta onto the BF axes. In Section II C 1 (i), we have chosen a very particular definition of the BF frame (see Figure 6), and, in Section II C 1 (ii), we have parametrized the vectors, in turn, by (poly)spherical coordinates. We have then provided the general expression of the KEO for orthogonal vectors (Section II C 1 (iv)) and for non-orthogonal vectors (Appendix B) in terms of the angular momenta. From these general expressions, we have deduced, in Section II C 2, another general expression of the KEO in terms of the conjugate momenta of the 3*N*-6 internal polyspherical coordinates. Finally, we have explained how to combine the polyspherical approach with the separation into subsystems. The present Part of this review is devoted to the presentation of different applications treated in full dimensionality, using the results obtained in Part II.

For concrete applications, two cases must be distinguished:

(1) Geometries corresponding to angular singularities in the KEO are accessible during the physical process. This is usually true for so-called *flexible systems*, i.e. systems presenting motions of large amplitude such as scattering systems. In this case, one should use the KEO in terms of angular momenta, see Eqs. (92,B5,B6,B7) together with the basis
set of spherical harmonics of Eq. (98). The action of the KEO is rather simple in this basis set and the singularities are removed analytically. In a grid representation, however, one has to deal with complicated multi-dimensional DVRs.

(2) The potential prevents the system from reaching singular geometries. Such a system is called semi-rigid, although it may exhibit large amplitude motions. For semi-rigid systems one should utilize the expression of the KEO in derivative operators $\frac{\partial}{\partial q_i}$ of Eqs. (111,C4,C5,C6). The number of terms in the KEO increases, but only simple one-dimensional primitive functions (DVRs or FBRs) are required.

## A. Flexible systems

### 1. Tri- and tetra-atomic flexible systems

For the explicit expression of the KEOs in angular momenta for tri-atomic systems with the conventions of Section II C 1, we refer the reader to Ref. [130] for orthogonal vectors and to Ref. [131] for non-orthogonal vectors. Let us focus on tetra-atomic systems and consider an example: the $H_2+H_2$ inelastic collision studied in Ref. [132].

In order to describe this process, it is natural to parametrize the system by three Jacobi vectors as shown in Figure 1. The BF frame is oriented such that $\vec{e}_{zBF}$ is parallel to the vector $\vec{R}_3$ and that $\vec{R}_2$, $\vec{e}_{xBF}$ and $\vec{e}_{zBF}$ lie in the same half-plane (or, in other words, $\vec{R}_2$ is parallel to the $((x^{BF}Gz^{BF}), x^{BF} > 0)$ half-plane). Before and after the collision, the two molecules can freely rotate, and therefore the singularities corresponding to $\theta_i^{BF} = 0$ or $\pi$ are accessible. Consequently, it is the expression of the KEO in terms of angular momenta of Eq. (92) that must be applied.

The reduced masses of Eq. (39) become

$$ \{\mu_1 = \mu_2 = \frac{m_H}{2}, \mu_3 = m_H\}; \qquad (134) $$

and the KEO of Eq. (76) becomes
$$2\hat{T} = \frac{(\hat{J}^{\dagger}\hat{J} - (\hat{L}_1^{\dagger} + \hat{L}_2^{\dagger})\hat{J} - \hat{J}^{\dagger}(\hat{L}_1 + \hat{L}_2) + \hat{L}_1^{\dagger}\hat{L}_1 + \hat{L}_2^{\dagger}\hat{L}_2 + \hat{L}_2^{\dagger}\hat{L}_1 + \hat{L}_1^{\dagger}\hat{L}_2)_{BF}}{\mu_3 R_3^2} + \sum_{i=1}^{3} \frac{\hat{P}_{R_i}^{\dagger} \hat{P}_{R_i}}{\mu_i} + \frac{(\hat{L}_1^{\dagger}\hat{L}_1)_{BF}}{\mu_1 R_1^2} + \frac{(\hat{L}_2^{\dagger}\hat{L}_2)_{BF}}{\mu_2 R_2^2}. \qquad (135)$$

Furthermore, in view of (86),

$$(\hat{L}_1^\dagger + \hat{L}_2^\dagger)\hat{J} + \hat{J}^\dagger(\hat{L}_1 + \hat{L}_2) = (2\hat{J}(\hat{L}_1 + \hat{L}_2))_{BF}, \qquad (136)$$

and

$$\hat{L}_2^\dagger \hat{L}_1 + \hat{L}_1^\dagger \hat{L}_2 = 2(\hat{L}_1 \hat{L}_2)_{BF}. \qquad (137)$$

In Eqs. (136,137), the order of the operators in the scalar products is again strictly fixed. The BF-components of $\hat{L}_2$, the angular momentum associated with $\vec{R}_2$, the vector that is involved in the definition of the BF frame, must remain on the right hand side of the scalar products. Eq. (92) becomes

$$\begin{aligned} \hat{T} = & \sum_{i=1}^{3} \left( -\frac{1}{2\mu_i} \frac{1}{R_i} \frac{\partial^2}{\partial R_i^2} R_i \right) + \sum_{i=1}^{2} \left( \frac{1}{2\mu_3 R_3^2} + \frac{1}{2\mu_i R_i^2} \right) \hat{\vec{L}}_i^\dagger \hat{\vec{L}}_i \\ & + \frac{\hat{J}_{BF}^2}{2\mu_3 R_3^2} - \frac{\hat{J}_{z^{BF}}(\hat{J}_{z^{BF}} - \hat{L}_{1z^{BF}}) + \hat{L}_{1z^{BF}} \hat{L}_{1z^{BF}}}{\mu_3 R_3^2} \\ & + \frac{(\hat{L}_{1+z^{BF}} \hat{L}_{2-z^{BF}} + \hat{L}_{1-z^{BF}} \hat{L}_{2+z^{BF}}) - \sum_{j=1}^{2} (\hat{J}_{+z^{BF}} \hat{L}_{j-z^{BF}} + \hat{J}_{-z^{BF}} \hat{L}_{j+z^{BF}})}{2\mu_3 R_3^2}. \end{aligned} \qquad (138)$$

The angular basis of Eq. (98) becomes

$$\langle \text{angles} | J, \Omega, \ell_1, \Omega_1, \ell_2, \Omega_2 \rangle = (2\pi)^{1/2} Y_J^{\Omega}(\beta, \gamma) (-1)^{\Omega} Y_{\ell_2}^{\Omega_2}(\theta_2^{BF}, 0) Y_{\ell_1}^{\Omega_1}(\theta_1^{BF}, \varphi_1^{BF}), \qquad (139)$$

with $\Omega = \Omega_1 + \Omega_2$. The action of the projections of angular momenta appearing in Eq. (138) onto the angular basis functions of Eq. (139) is the usual one and is given by Eq. (101). In particular, if we use the volume element

$$\sin \beta \, d\beta \, d\gamma \, dR_1 \, dR_2 \, dR_3 \sin \theta_1^{BF} d\theta_1^{BF} \sin \theta_2^{BF} d\theta_2^{BF} d\varphi_1^{BF}, \qquad (140)$$
Eq. (138) becomes

$$
\begin{equation}
\begin{split}
\hat{T} = \sum_{i=1}^{3} \left( -\frac{1}{2\mu_i} \frac{\partial^2}{\partial R_i^2} \right) &+ \sum_{i=1}^{2} \left( \frac{1}{2\mu_3 R_3^2} + \frac{1}{2\mu_i R_i^2} \right) \hat{\vec{L}}_i^\dagger \hat{\vec{L}}_i \\
& + \frac{\hat{J}_{BF}^2}{2\mu_3 R_3^2} - \frac{\hat{J}_{z^{BF}}(\hat{J}_{z^{BF}} - \hat{L}_{1z^{BF}}) + \hat{L}_{1z^{BF}} \hat{L}_{1z^{BF}}}{\mu_3 R_3^2} \\
& + \frac{(\hat{L}_{1+z^{BF}} \hat{L}_{2-z^{BF}} + \hat{L}_{1-z^{BF}} \hat{L}_{2+z^{BF}}) - \sum_{j=1}^{2} (\hat{J}_{+z^{BF}} \hat{L}_{j-z^{BF}} + \hat{J}_{-z^{BF}} \hat{L}_{j+z^{BF}})}{2\mu_3 R_3^2}.
\end{split}
\tag{141}
\end{equation}
$$

In addition, in view of Eq. (101), we obtain

$$
\begin{align}
\hat{\vec{L}}_1^\dagger \hat{\vec{L}}_1 Y_{\ell_1}^{\Omega_1}(\theta_1^{BF}, \varphi_1^{BF}) Y_J^\Omega(\beta, \gamma) &= \ell_1(\ell_1 + 1) Y_{\ell_1}^{\Omega_1}(\theta_1^{BF}, \varphi_1^{BF}) Y_J^\Omega(\beta, \gamma), \nonumber \\
\hat{\vec{L}}_2^\dagger \hat{\vec{L}}_2 Y_{\ell_2}^{\Omega_2}(\theta_2^{BF}, 0) Y_J^\Omega(\beta, \gamma) &= \ell_2(\ell_2 + 1) Y_{\ell_2}^{\Omega_2}(\theta_2^{BF}, 0) Y_J^\Omega(\beta, \gamma), \nonumber \\
\hat{\vec{L}}_{1z^{BF}} Y_{\ell_1}^{\Omega_1}(\theta_1^{BF}, \varphi_1^{BF}) Y_J^\Omega(\beta, \gamma) &= \Omega_1 Y_{\ell_1}^{\Omega_1}(\theta_1^{BF}, \varphi_1^{BF}) Y_J^\Omega(\beta, \gamma), \nonumber \\
\hat{\vec{L}}_{1z^{BF}} Y_{\ell_2}^{\Omega_2}(\theta_2^{BF}, 0) Y_J^\Omega(\beta, \gamma) &= \Omega_2 Y_{\ell_2}^{\Omega_2}(\theta_2^{BF}, 0) Y_J^\Omega(\beta, \gamma), \nonumber \\
\hat{L}_{1+z^{BF}} \hat{L}_{2-z^{BF}} Y_{\ell_1}^{\Omega_1}(\theta_1^{BF}, \varphi_1^{BF}) Y_{\ell_2}^{\Omega_2}(\theta_2^{BF}, 0) Y_J^\Omega(\beta, \gamma) &= \nonumber \\
&\quad C_+(\ell_1, \Omega_1) C_-(\ell_2, \Omega_2) Y_{\ell_1}^{\Omega_1+1}(\theta_1^{BF}, \varphi_1^{BF}) Y_{\ell_2}^{\Omega_2-1}(\theta_2^{BF}, 0) Y_J^\Omega(\beta, \gamma), \nonumber \\
\hat{L}_{1-z^{BF}} \hat{L}_{2+z^{BF}} Y_{\ell_1}^{\Omega_1}(\theta_1^{BF}, \varphi_1^{BF}) Y_{\ell_2}^{\Omega_2}(\theta_2^{BF}, 0) Y_J^\Omega(\beta, \gamma) &= \nonumber \\
&\quad C_-(\ell_1, \Omega_1) C_+(\ell_2, \Omega_2) Y_{\ell_1}^{\Omega_1-1}(\theta_1^{BF}, \varphi_1^{BF}) Y_{\ell_2}^{\Omega_2+1}(\theta_2^{BF}, 0) Y_J^\Omega(\beta, \gamma), \nonumber \\
\hat{J}_{\mp z^{BF}} \hat{\vec{L}}_{1\pm z^{BF}} Y_{\ell_1}^{\Omega_1}(\theta_1^{BF}, \varphi_1^{BF}) Y_J^\Omega(\beta, \gamma) (-1)^\Omega &= \nonumber \\
&\quad C_\pm(J, \Omega) C_\pm(\ell_1, \Omega_1) Y_{\ell_1}^{\Omega_1\pm 1}(\theta_1^{BF}, \varphi_1^{BF}) Y_J^{(\Omega\pm 1)}(\beta, \gamma) (-1)^{(\Omega\pm 1)}, \nonumber \\
\hat{J}_{\mp z^{BF}} \hat{\vec{L}}_{2\pm z^{BF}} Y_{\ell_2}^{\Omega_2}(\theta_2^{BF}, 0) Y_J^\Omega(\beta, \gamma) (-1)^\Omega &= \nonumber \\
&\quad C_\pm(J, \Omega) C_\pm(\ell_2, \Omega_2) Y_{\ell_2}^{\Omega_2\pm 1}(\theta_2^{BF}, 0) Y_J^{(\Omega\pm 1)}(\beta, \gamma) (-1)^{(\Omega\pm 1)}. \tag{142}
\end{align}
$$

The KEO of Eq. (138) was also applied to the study of the photo-dissociation dynamics
of the Ar₂HBr van der Waals molecule with MCTDH [133]. Obviously, the reduced masses
are different for Ar₂HBr (see Ref. [133] for their actual values). In addition, it goes without
saying that the KEO of Eq. (138) for standard diatom-diatom Jacobi coordinates is not a
new result. Some operators similar to Eq. (138) have already been applied in the past, often
with different conventions (see for instance [134–137]).
FIG. 12: Definition of the five Jacobi vectors coordinates for the $(H_2)_3$ trimer. Here, $g_i$ is the center of mass of the i-th monomer $H_2$, $G_1$ is the center of mass of the $(H_2)_2$ subsystem.

<figure><img src="image_13.png" /></figure>

## 2. Weakly bound molecular trimers

Weakly bound molecular systems feature several minimum energy structures and large amplitude motions between regions connected by low energy barriers. Consequently, they can easily reach singular geometries. Let us consider an example in molecular spectroscopy, the calculation of the vibrational energy levels of the hydrogen trimer $(H_2)_3$ achieved by Costa and Clary [138] and Yu [139]. The system is described by five Jacobi vectors as shown in Figure 12. The five reduced masses read

$$ \{\mu_1 = \mu_2 = \mu_3 = \frac{m_H}{2}, \mu_4 = m_H, \mu_5 = \frac{4}{3}m_H\}; \qquad (143) $$

and Eq. (92) straightforwardly leads to
$$
\begin{equation}
\begin{split}
\hat{T} = \sum_{i=1}^{5} \left( -\frac{1}{2\mu_i} \frac{1}{R_i} \frac{\partial^2}{\partial R_i^2} R_i \right) &+ \sum_{i=1}^{4} \left( \frac{1}{2\mu_4 R_4^2} + \frac{1}{2\mu_i R_i^2} \right) \hat{L}_i^\dagger \hat{L}_i \\
&+ \frac{\hat{J}_{BF}^2}{2\mu_5 R_5^2} - \frac{2\hat{J}_{zBF}(\hat{J}_{zBF} - \sum_{i=1}^3 \hat{L}_{izBF}) + \sum_{i=1}^3 \hat{L}_{izBF}(2\hat{L}_{izBF} + \sum_{j=1, j \neq i}^3 \hat{L}_{jzBF})}{2\mu_5 R_5^2} \\
&+ \frac{\sum_{i=1}^3 \sum_{j=1, j>i}^4 (\hat{L}_{i+BF} \hat{L}_{j-BF} + \hat{L}_{i-BF} \hat{L}_{j+BF}) - \sum_{j=1}^4 (\hat{J}_{+BF} \hat{L}_{j-BF} + \hat{J}_{-BF} \hat{L}_{j+BF})}{2\mu_5 R_5^2}.
\end{split}
\tag{144}
\end{equation}
$$

The angular basis of Eq. (98) becomes

$$
\begin{equation}
\begin{split}
\langle \text{angles} | J, \Omega, \ell_1, \Omega_1, \ell_2, \Omega_2, \ell_3, \Omega_3, \ell_4, \Omega_4 \rangle = & \\
& (2\pi)^{1/2} Y_J^{\Omega_1}(\beta, \gamma) (-1)^{\Omega_1} Y_{\ell_4}^{\Omega_4}(\theta_4^{BF}, 0) Y_{\ell_3}^{\Omega_3}(\theta_3^{BF}, \varphi_3^{BF}) Y_{\ell_2}^{\Omega_2}(\theta_2^{BF}, \varphi_2^{BF}) Y_{\ell_1}^{\Omega_1}(\theta_1^{BF}, \varphi_1^{BF}),
\end{split}
\tag{145}
\end{equation}
$$

with $\Omega = \Omega_1 + \Omega_2 + \Omega_3 + \Omega_4$. In Ref. [138], as well as in Eq. (38) of Ref. [79], appears the 'sign of (k-K)' in the matrix elements (K-k corresponds to our $\Omega_{N-2}$). This sign does not appear in Eq. (101). This difference can be traced back to a slightly different convention in the angular basis set. Indeed, Ref. [138] and Ref. [79] use $P_{\ell}^{|\Omega_4|}(\cos(\theta_4^{BF}))$ instead of our

$$
\tilde{P}_{\ell}^{\Omega_4}(\cos(\theta_4^{BF})) = (-1)^{\Omega_4} P_{\ell}^{\Omega_4}(\cos(\theta_4^{BF})) = (2\pi)^{1/2} Y_{\ell_4}^{\Omega_4}(\theta_4^{BF}, 0) \quad (146)
$$

in Eq. (145) (we recall that $P_{\ell}^{\Omega_4}(\cos(\theta_4^{BF}))$ is a normalized associated Legendre function).

3. Protonated methane in Radau coordinates

Isolated CH<sub>5</sub><sup>+</sup> is a stable but very floppy system, the protons continually hopping between
120 minima [140]. Very recently, Wang and Carrington [140] reported the first vibrational
energy levels including all the 12 internal degrees of freedom. They adopted the polyspherical
approach based on a parametrization of five Radau vectors, as shown in Figure 13. For the
reader who is not familiar with Radau coordinates, let us add that these are often introduced
to replace valence coordinates in describing the vibrational motions of molecules in which
the central atom is much heavier than the end atoms [141]. For this type of molecules, both
systems of coordinates are very similar, the former having the advantage of being orthogonal.

The Radau vectors, $\vec{R}_i$, are the vectors $\vec{B}\vec{H}_i$, B being the canonical point of CH<sub>5</sub><sup>+</sup> defined
as the geometric mean of CG and GG<sub>H5</sub>. Here, C is the C atom, G is the center of mass
FIG. 13: Definition of the five Radau vectors coordinates for the CH5$^+$ system. Here, B is the canonical point of CH$_5^+$, G is the center of mass of CH$_5^+$, G$_{H5}$ is the center of mass of the H$_5$ subsystem.

<figure><img src="image_14.png" /></figure>

of CH$_5^+$, G$_{H5}$ is the center of mass of the H$_5$ subsystem. The canonical point is defined by $(BG_{H5})^2 = GG_{H5} \times CG_{H5}$ [141]. Let us also put $\vec{r}_i = \vec{G} \vec{H}_i$. It entails that

$$ 
\begin{aligned}
\vec{R}_1 &= (1 - \alpha m_H)\vec{r}_1 - \alpha m_H \vec{r}_2 - \alpha m_H \vec{r}_3 - \alpha m_H \vec{r}_4 - \alpha m_H \vec{r}_5, \\
\vec{R}_2 &= -\alpha m_H \vec{r}_1 + (1 - \alpha m_H)\vec{r}_2 - \alpha m_H \vec{r}_3 - \alpha m_H \vec{r}_4 - \alpha m_H \vec{r}_5, \\
\vec{R}_3 &= -\alpha m_H \vec{r}_1 - \alpha m_H \vec{r}_2 + (1 - \alpha m_H)\vec{r}_3 - \alpha m_H \vec{r}_4 - \alpha m_H \vec{r}_5, \\
\vec{R}_4 &= -\alpha m_H \vec{r}_1 - \alpha m_H \vec{r}_2 - \alpha m_H \vec{r}_3 + (1 - \alpha m_H)\vec{r}_4 - \alpha m_H \vec{r}_5, \\
\vec{R}_5 &= -\alpha m_H \vec{r}_1 - \alpha m_H \vec{r}_2 - \alpha m_H \vec{r}_3 - \alpha m_H \vec{r}_4 + (1 - \alpha m_H)\vec{r}_5,
\end{aligned}
\qquad (147) $$

with

$$ \alpha = (1 - \sqrt{\frac{M}{m_C}})/m_{H_5}, \qquad (148) $$
and $M=m_C+5m_H$, $m_{H_5} = 5m_H$. The KEO is formally identical to the operator of Eq. (144) except that $\mu_i = m_H$, with $(i = 1, \dots, 5)$.

## B. Semi-rigid systems

### 1. Tetra-atomic molecules

#### (i) Tetra-atomic molecules in orthogonal coordinates:

If the potential prevents the system from reaching singular geometries, the KEO in derivative operators of Eqs. (111,C4,C5,C6) can be used. Let us detail two examples to explain how one proceeds. The two examples are the HOOH molecule in Jacobi vectors as shown in Figure 1, and the NHD$_2$ molecule in Radau vectors as shown in Figure 14. If one aims at calculating the infra-red spectrum of HOOH as in Ref. [65] or the tunneling dynamics of the NH chromophore in NHD2 after coherent infrared excitation as in Ref. [142], no singularity is accessible. The KEO for $J = 0$ is given by Eqs. (111,C4) for $N = 4$ and with the following matrix $M$:

$$ M = \begin{bmatrix} \frac{1}{\mu_1} & 0 & 0 \\ 0 & \frac{1}{\mu_2} & 0 \\ 0 & 0 & \frac{1}{\mu_3} \end{bmatrix}, \qquad (149) $$

with

$$ \{\mu_1 = \mu_2 = \frac{m_O m_H}{m_O + m_H}, \mu_3 = \frac{m_O + m_H}{2}\} \qquad (150) $$

for HOOH in Jacobi coordinates and $\{\mu_1 = \mu_2 = m_D, \mu_3 = m_H\}$ for NHD$_2$ in Radau coordinates. The three vectors are in turn parametrized by six internal polyspherical coordinates: $R_1, R_2, R_3, u_1, u_2, \varphi_1^{BF}$. Eq. (C7) yields (we use $\varphi_2^{BF} = 0$)
<figure><img src="image_15.png" /><figcaption>FIG. 14: Definition of the three Radau vectors for the NHD₂ molecule. Here, B is the canonical point of NHD₂, G is the center of mass of the molecule, G<sub>HD2</sub> is the center of mass of the HD₂</figcaption></figure>

$$G_{R_i R_i} = \frac{1}{\mu_i}; \quad i = 1, 2, 3,$$

$$G_{u_i u_i} = (\sin \theta_i^{BF})^2 (\frac{1}{\mu_i R_i^2} + \frac{1}{\mu_3 R_3^2}); \quad i = 1, 2,$$

$$G_{u_1 u_2} = \frac{\sin \theta_1^{BF} \sin \theta_2^{BF} \cos \varphi_1^{BF}}{\mu_3 R_3^2},$$

$$G_{u_1 \varphi_1^{BF}} = \frac{\sin \theta_1^{BF} \sin \varphi_1^{BF} \cot \theta_2^{BF}}{\mu_3 R_3^2},$$

$$G_{u_2 \varphi_1^{BF}} = \frac{\sin \theta_2^{BF} \sin \varphi_1^{BF} \cot \theta_1^{BF}}{\mu_3 R_3^2},$$

$$G_{\varphi_1^{BF} \varphi_1^{BF}} = \frac{1}{\mu_3 R_3^2} ((\cot \theta_1^{BF})^2 + (\cot \theta_2^{BF})^2) + \frac{1}{\mu_1 (\sin \theta_1^{BF})^2 R_1^2} + \frac{1}{\mu_2 (\sin \theta_2^{BF})^2 R_2^2} \\ - \frac{2}{\mu_3 R_3^2} \cot \theta_2^{BF} \cos \varphi_1^{BF} \cot \theta_1^{BF}. \qquad (151)$$

All the other matrix elements (for instance $G_{R_i u_j}$ or $G_{R_i \varphi_j^{BF}}$) are equal to zero. Now, the choice of the following volume element in normalizing wave-functions
$$d\tau' = dR_3 dR_2 dR_1 du_2 du_1 d\varphi_1^{BF}, \qquad (152)$$

entails that the KEO reads

$$ \begin{aligned} \hat{T} = & - \sum_{i=1,2,3} \frac{1}{2\mu_i} \frac{\partial^2}{\partial R_i^2} - \sum_{i=1}^2 \left( \frac{1}{2\mu_i R_i^2} + \frac{1}{2\mu_3 R_3^2} \right) \frac{\partial}{\partial u_i} v_i^2 \frac{\partial}{\partial u_i} \\ & - \frac{\cos \varphi_1^{BF}}{2\mu_3 R_3^2} \left( v_1 \frac{\partial^2}{\partial u_1 \partial u_2} v_2 + v_2 \frac{\partial^2}{\partial u_1 \partial u_2} v_1 \right) \\ & - \frac{u_2}{2\mu_3 R_3^2 v_2} \left( \sin \varphi_1^{BF} \frac{\partial^2}{\partial u_1 \partial \varphi_1^{BF}} v_1 + v_1 \frac{\partial^2}{\partial u_1 \partial \varphi_1^{BF}} \sin \varphi_1^{BF} \right) \\ & - \frac{u_1}{2\mu_3 R_3^2 v_1} \left( \sin \varphi_1^{BF} \frac{\partial^2}{\partial u_2 \partial \varphi_1^{BF}} v_2 + v_2 \frac{\partial^2}{\partial u_2 \partial \varphi_1^{BF}} \sin \varphi_1^{BF} \right) \\ & - \sum_{i=1}^2 \left( \frac{1}{2\mu_i R_i^2 v_i^2} + \frac{u_i^2}{2R_3^2 \mu_3 v_i^2} \right) \frac{\partial^2}{\partial \varphi_1^{BF 2}} \\ & + \frac{u_2 u_1}{\mu_3 v_2 v_1 R_3^2} \frac{\partial}{\partial \varphi_1^{BF}} \cos \varphi_1^{BF} \frac{\partial}{\partial \varphi_1^{BF}}, \end{aligned} \qquad (153) $$

Here, and in the following, we use the notation:

$$v_i = \sqrt{1 - u_i^2} = \sin \theta_i^{BF}. \qquad (154)$$

Note that this KEO applies to both HOOH and NHD₂, despite the two molecules are described by different coordinates systems. In fact, Eq. (153) applies to all 4-atom systems described by orthogonal coordinates, only the masses, $\mu_i$, have different values. Furthermore, the KEO of Eq. (153) is identical to the operator given in Eq. (A.5) of Ref. [65] by Bramley and Carrington (they used Mathematica to calculate the operator). Note, however, that these authors do not use the same conventions for the coordinates, for instance their 'r2' is our $R_3$. In addition, they use $\theta_i^{BF}$ instead of $u_i = \cos \theta_i^{BF}$: this is why there is a difference of sign for the terms $\frac{\partial^2}{\partial \theta_i \partial \varphi_1^{BF}}$ (see Eq. (107) to understand the origin of the change of the sign).

(ii) Tetra-atomic molecules in valence coordinates:
FIG. 15: Definition of the three valence vectors for the HONO molecule. Here, $\theta_2^{new} = \pi - \theta_2^{BF}$.

<figure><img src="image_16.png" /></figure>

(1) A-B-C-D structure:

Let us now consider a non-orthogonal parametrization of tetra-atomic molecules. For valence vectors, only two 'clustering' schemes are possible, denoted as A-B-C-D and A(CBD) structures. They are shown in Figures 15 and 16, respectively.

In Ref. [143-145], the trans-cis Isomerization of HONO was studied with MCTDH. A parametrization of three valence vectors as shown in Figure 15 was adopted. This results in the following matrix $\mathbf{M}$:

$$ \mathbf{M} = \begin{bmatrix} \frac{1}{m_H} + \frac{1}{m_O} & 0 & \frac{1}{m_O} \\ 0 & \frac{1}{m_O} + \frac{1}{m_N} & -\frac{1}{m_N} \\ \frac{1}{m_O} & -\frac{1}{m_N} & \frac{1}{m_O} + \frac{1}{m_N} \end{bmatrix}. \qquad (155) $$

In order to be consistent with the conventions used for the potential energy surface calculated in Ref. [143], we also slightly changed the convention for the angle $\theta_2^{BF} \rightarrow \theta_2^{new} = \pi - \theta_2^{BF}$ (see Figure 15). Because of this, there appears the phase factor $S_i$ in Eq. (A.2) of Ref. [144]. The internal motions are described by six coordinates: $R_1, R_2, R_3, u_1 = \cos \theta_1^{BF}$, $u_2 = \cos \theta_2^{new}$, $\phi_1^{BF}$. If the following volume element,
$$d\tau' = dR_3 dR_2 dR_1 du_2 du_1 d\varphi_1^{BF}, \qquad (156)$$

is used in normalizing wave-functions, the KEO is obtained from Eqs. (111,C4) and its explicit expression is given by Eq. (A.2) of Ref. [144].

### (2) A(BCD) structure:

Let us consider the second 'clustering' scheme for a tetra-atomic molecule parametrized by valence vectors. In Ref. [146-148], the Intra-molecular energy redistribution (IVR) in HFCO and DFCO with MCTDH was studied. The valence parametrization shown in Figure 16 was adopted.

The matrix M reads

$$M = \begin{bmatrix} \frac{1}{m_H} + \frac{1}{m_C} & \frac{1}{m_C} & \frac{1}{m_C} \\ \frac{1}{m_C} & \frac{1}{m_F} + \frac{1}{m_C} & \frac{1}{m_C} \\ \frac{1}{m_C} & \frac{1}{m_C} & \frac{1}{m_O} + \frac{1}{m_C} \end{bmatrix}. \qquad (157)$$

If we use the volume element of Eq. (156) in normalizing wave-functions, the KEO can be obtained from Eqs. (111,C4) and its explicit expression is given by Eq. (15) of Ref. [146].

The KEO of Eq. HFCO shows more terms than those of HONO and NHD₂. The two last-mentioned operators can be derived from the KEO of HFCO just by replacing the matrix elements of M appearing in Eq. (157) by those of Eq. (149) and (155), respectively. It is worth noting that the KEO used in Ref. [149] to study the selective population of the vibrational levels of H₂CS with MCTDH is formally identical to the operator of HFCO. The only difference is that the matrix M, for H₂CS, reads

$$M = \begin{bmatrix} \frac{1}{m_H} + \frac{1}{m_C} & \frac{1}{m_C} & \frac{1}{m_C} \\ \frac{1}{m_C} & \frac{1}{m_H} + \frac{1}{m_C} & \frac{1}{m_C} \\ \frac{1}{m_C} & \frac{1}{m_C} & \frac{1}{m_C} + \frac{1}{m_S} \end{bmatrix}. \qquad (158)$$
<figure><img src="image_17.png" /><figcaption>FIG. 16: Valence BF polyspherical coordinates for the HFCO system. $\vec{R}_2$ belongs to the xz plane.</figcaption></figure>

## 2. Penta-atomic molecules

### (i) Orthogonal coordinates:

In Ref. [150–152], the computation of vibrational levels of methane was reported. A set of Radau polyspherical coordinates was used. The Radau vectors, for a molecule such as methane, are identical to those depicted in Figure 13 except that there are only four Radau vectors for $CH_4$ instead of five for $CH_5^+$. The internal motions are described by nine coordinates: $R_1, R_2, R_3, R_4, u_1 = \cos \theta_1^{BF}, u_2 = \cos \theta_2^{BF}, u_3 = \cos \theta_3^{BF}, \varphi_1^{BF}$, and $\varphi_2^{BF}$. If the following non-Euclidean volume element,
$$d\tau' = dR_4 dR_3 dR_2 dR_1 du_3 du_2 du_1 d\varphi_2^{BF} d\varphi_1^{BF}, \qquad (159)$$

is used, the KEO of Eqs. (111,C4) becomes (we use again $\varphi_i$ instead of $\varphi_i^{BF}$ to simplify the notation)

$$
\begin{aligned}
\hat{T} = & -\sum_{i=1}^{4} \frac{M_{ii}}{2} \frac{\partial^2}{\partial R_i^2} - \sum_{i=1}^{3} \left( \frac{1}{2\mu_i R_i^2} + \frac{1}{2\mu_4 R_4^2} \right) \frac{\partial}{\partial u_i} v_i^2 \frac{\partial}{\partial u_i} \\
& - \frac{\cos \varphi_1 \cos \varphi_2 + \sin \varphi_1 \sin \varphi_2}{2\mu_4 R_4^2} \left( v_2 \frac{\partial}{\partial u_2} \frac{\partial}{\partial u_1} v_1 + v_1 \frac{\partial}{\partial u_1} \frac{\partial}{\partial u_2} v_2 \right) \\
& - \frac{\cos \varphi_1}{2\mu_4 R_4^2} \left( v_3 \frac{\partial}{\partial u_3} \frac{\partial}{\partial u_1} v_1 + v_1 \frac{\partial}{\partial u_1} \frac{\partial}{\partial u_3} v_3 \right) - \frac{\cos \varphi_2}{2\mu_4 R_4^2} \left( v_3 \frac{\partial}{\partial u_3} \frac{\partial}{\partial u_2} v_2 + v_2 \frac{\partial}{\partial u_2} \frac{\partial}{\partial u_3} v_3 \right) \\
& - \frac{u_3}{2v_3\mu_4 R_4^2} \left[ \frac{\partial}{\partial u_1} v_1 \sin \varphi_1 \frac{\partial}{\partial \varphi_1} + v_1 \frac{\partial}{\partial u_1} \frac{\partial}{\partial \varphi_1} \sin \varphi_1 + \left( \frac{\partial}{\partial u_2} v_2 \sin \varphi_2 \frac{\partial}{\partial \varphi_2} + v_2 \frac{\partial}{\partial u_2} \frac{\partial}{\partial \varphi_2} \sin \varphi_2 \right) \right] \\
& + \frac{1}{2\mu_4 R_4^2} \left[ \frac{\partial}{\partial u_1} v_1 \left( \frac{u_2}{v_2} (\sin \varphi_1 \cos \varphi_2 - \cos \varphi_1 \sin \varphi_2) - \frac{u_3}{v_3} \sin \varphi_1 \right) \frac{\partial}{\partial \varphi_2} \right. \\
& \left. + v_1 \frac{\partial}{\partial u_1} \frac{\partial}{\partial \varphi_2} \left( \frac{u_2}{v_2} (\sin \varphi_1 \cos \varphi_2 - \cos \varphi_1 \sin \varphi_2) - \frac{u_3}{v_3} \sin \varphi_1 \right) \right] \\
& + \frac{1}{2\mu_4 R_4^2} \left[ \frac{\partial}{\partial u_2} v_2 \left( \frac{u_1}{v_1} (\sin \varphi_2 \cos \varphi_1 - \cos \varphi_2 \sin \varphi_1) - \frac{u_3}{v_3} \sin \varphi_2 \right) \frac{\partial}{\partial \varphi_1} \right. \\
& \left. + v_2 \frac{\partial}{\partial u_2} \frac{\partial}{\partial \varphi_1} \left( \frac{u_1}{v_1} (\sin \varphi_2 \cos \varphi_1 - \cos \varphi_2 \sin \varphi_1) - \frac{u_3}{v_3} \sin \varphi_2 \right) \right] \\
& - \sum_{i=1}^{2} \left\{ \frac{\partial}{\partial u_3} v_3 \frac{\sin \varphi_i \cot \theta_i}{2\mu_4 R_4^2} \frac{\partial}{\partial \varphi_i} + \frac{\partial}{\partial \varphi_i} v_3 \frac{\sin \varphi_i \cot \theta_i}{2\mu_4 R_4^2} \frac{\partial}{\partial u_3} \right\} \\
& - \frac{1}{2\mu_4 R_4^2} \left[ \left( \frac{u_1^2}{v_1^2} + \frac{u_3^2}{v_3^2} \right) \frac{\partial^2}{\partial \varphi_1^2} - 2 \frac{u_3 u_1}{v_3 v_1} \frac{\partial}{\partial \varphi_1} \cos \varphi_1 \frac{\partial}{\partial \varphi_1} \right] \\
& - \frac{1}{2\mu_4 R_4^2} \left[ \left( \frac{u_2^2}{v_2^2} + \frac{u_3^2}{v_3^2} \right) \frac{\partial^2}{\partial \varphi_2^2} - 2 \frac{u_3 u_2}{v_3 v_2} \frac{\partial}{\partial \varphi_2} \cos \varphi_2 \frac{\partial}{\partial \varphi_2} \right] \\
& - \sum_{i=1,2} \left( \frac{1}{2\mu_i R_i^2 v_i^2} + \frac{1}{2\mu_3 R_3^2 v_3^2} \right) \frac{\partial^2}{\partial \varphi_i^2} - \frac{u_3^2}{\mu_4 R_4^2 v_3^2} \frac{\partial^2}{\partial \varphi_1 \partial \varphi_2} - \frac{1}{\mu_3 R_3^2 v_3^2} \frac{\partial^2}{\partial \varphi_1 \partial \varphi_2} \\
& - \frac{1}{2\mu_4 R_4^2} \frac{\partial}{\partial \varphi_1} \left[ \frac{u_1 u_2}{v_1 v_2} (\cos \varphi_1 \cos \varphi_2 + \sin \varphi_1 \sin \varphi_2) - \frac{u_3}{v_3} \left( \cos \varphi_1 \frac{u_1}{v_1} + \cos \varphi_2 \frac{u_2}{v_2} \right) \right] \frac{\partial}{\partial \varphi_2} \\
& - \frac{\partial}{\partial \varphi_2} \left[ \frac{u_1 u_2}{v_1 v_2} (\cos \varphi_1 \cos \varphi_2 + \sin \varphi_1 \sin \varphi_2) - \frac{u_3}{v_3} \left( \cos \varphi_1 \frac{u_1}{v_1} + \cos \varphi_2 \frac{u_2}{v_2} \right) \right] \frac{\partial}{\partial \varphi_1},
\end{aligned}
\qquad (160)
$$

with $\mu_i = m_H$, and ($i = 1, \dots, 4$). Of course, the operator of Eq. (160) could be applied to any set of orthogonal vectors for a penta-atomic system. One simply has to change the four
<figure><img src="image_18.png" /><figcaption>FIG. 17: Definition of the four valence vectors for the CHF$_3$ molecule.</figcaption></figure>

reduced masses.

(ii) Non-orthogonal coordinates:

For the study of the vibrations of a molecule such as fluoroform, CHF$_3$, the previous Radau parametrization is no longer adapted since the F atoms are rather heavy. Hence, the canonical point, used to define the Radau vectors, is located far from the C atom. The resulting Radau vectors are therefore very different from the valence vectors and cannot correctly describe the vibrations in the normal or local mode domains. It is more pertinent to parametrize CHF$_3$ by four valence vectors, which are shown in Figure 17.

This results in a full matrix $M$ that reads

$$ M = \begin{bmatrix} \frac{1}{m_C} + \frac{1}{m_H} & \frac{1}{m_C} & \frac{1}{m_C} & \frac{1}{m_C} \\ \frac{1}{m_C} & \frac{1}{m_C} + \frac{1}{m_F} & \frac{1}{m_C} & \frac{1}{m_C} \\ \frac{1}{m_C} & \frac{1}{m_C} & \frac{1}{m_C} + \frac{1}{m_F} & \frac{1}{m_C} \\ \frac{1}{m_C} & \frac{1}{m_C} & \frac{1}{m_C} & \frac{1}{m_C} + \frac{1}{m_F} \end{bmatrix}. \quad (161) $$

If one uses the volume element of Eq. (159), the KEO is given by Eqs. (111,C4), as usual.
Its explicit expression can be found in Appendix D. Let us only point out that we have recently implemented this nine-dimensional operator in the Heidelberg MCTDH package [11–15, 129]. In addition, we have used the program TNUM of Lauvergnat and Nauts [68] to check the correctness of the operator [127]. As already mentioned, the test does not only guarantee that the KEO is correct but also that the operator is correctly implemented in the code.

### C. Separation into subsystems

Until now, we have confined ourselves to considering only particular cases of the standard polyspherical approach. For clusters involving molecules of water, the separation into subsystems is more appropriate than the standard parametrization since the separation into subsystems allows to have purely intramolecular angles instead of angles mixing the intra-monomer and inter-monomer motions. The present section aims at giving the explicit expression of the KEO for several protonated water clusters.

#### 1. Protonated water dimer (Zundel cation)

##### (i) Jacobi description:

In Ref. [153, 154], the simulation of the IR linear absorption spectrum of $H_5O_2^+$ was reported. A set of six Jacobi vectors was chosen [128], as depicted in Figure 18. The system is separated into three subsystems: the two molecules of water and the proton. Each molecule of water is parametrized by two Jacobi vectors: $\vec{R}_1^{(j)}$, the vector connecting the center of mass of the H₂ subsystem to the O atom, and $\vec{R}_2^{(j)}$, the vector joining one H atom to the other one (here, j = 1, 2). The two 'remaining' vectors are the vector $\vec{R}_2$, which connects the center of mass of one monomer to the center of mass of the other monomer, and $\vec{R}_1$, the Jacobi vector from the center of mass of the water dimer to the proton. Let $\hat{\vec{J}}^{(j)}$ be the total angular momentum of the j-th molecule of water. This definition entails that $\hat{\vec{J}}^{(j)} = \sum_{k=1,2} \hat{\vec{L}}_k^{(j)}, \hat{\vec{L}}_k^{(j)}$ being the angular momentum associated with $\vec{R}_k^{(j)}$.

According to the general strategy adopted in Section II D, the $z^{BF}$ axis should be parallel to $\vec{R}_2$ and $\vec{R}_1$ should be parallel to the $((x^{BF}, z^{BF}); x^{BF} > 0)$ half-plane. However, for the
FIG. 18: Jacobi description of the $H_5O_2^+$ system. $G_1$ is the center of mass of the first monomer, $G_2$ is the center of mass of the second monomer. $G_{H4O_2}$ is the center of the dimer, $H_4O_2$.

<figure><img src="image_19.png" /></figure>

definition of the $((x^{BF}, z^{BF}); x^{BF} > 0)$ half-plane in $H_5O_2^+$, it is more natural to choose a vector belonging to one monomer, let say $\vec{R}_2^{(2)}$ (it is, indeed, more appropriate to link the overall rotation of the Zundel cation to the molecules of water than to the proton alone since the mass of the latter is much smaller). Therefore, as already mentioned in Section IID (i), we have to face a particular case which is slightly more complicated than the general case thoroughly treated in Section IID. Using the properties of the angular momenta highlighted in Sec. IIB 2 (iii), we can already predict that the BF-projections of the total angular momenta of the second monomer have a non-regular behavior. As in Section IIC 1 (iii) for the BF-components of $\vec{L}_{N-2}$, one needs to introduce an intermediate frame, the $E_2$ frame.

But, prior to the introduction of the $E_2$ frame, let us define a BF frame for each molecule of water, as in Section IID (ii). Let $BF^{(j)}$ be the intermediate frame associated with the $j$-th molecule of water ($j = 1$ or 2) such as $\vec{e}_{zBF^{(j)}}$ is parallel to the vector $\vec{R}_2^{(j)}$ and $\vec{R}_1^{(j)}$, $\vec{e}_{xBF^{(j)}}$ and $\vec{e}_{zBF^{(j)}}$ lie in the same half-plane (as in Fig. 6). Clearly, for the molecule of water, the vector $\vec{R}_1^{(j)}$ would have been a better choice for the definition of the $z^{BF^{(j)}}$ axis since the reduced mass of $\vec{R}_1^{(j)}$ is larger than that of $\vec{R}_2^{(j)}$. However, in Ref. [128], we preferred to use $\vec{R}_2^{(j)}$ in order to avoid the singularity which appears in the KEO when $\vec{R}_2$
and the $z^{BF(j)}$ axis are parallel. The orientation of the BF(j) frame with respect to the E₂ frame is determined by three Euler angles: $\alpha^{(j)}$, $\beta^{(j)}$, $\gamma^{(j)}$. However, $\alpha^{(2)} = 0$ since the $((x^{BF}, z^{BF}); x^{BF} > 0)$ half-plane is parallel to $\vec{R}_2^{(2)}$.

For the KEO, one obtains in view of Eqs. (125,126)

$$2\hat{T} = \sum_{i=1}^{2} \frac{\hat{P}_{R_i}^{\dagger} \hat{P}_{R_i}}{\mu_i} + \frac{\hat{L}_1^{\dagger} \hat{L}_1}{\mu_1 R_1^2} + 2 \sum_{j=1}^{2} \hat{T}^{(j)} + \frac{(\hat{J}^{\dagger} - (\hat{L}_1^{\dagger} + \sum_{j=1}^{2} \hat{J}^{(j)})^{\dagger}) \cdot (\hat{J} - (\hat{L}_1 + \sum_{j=1}^{2} \hat{J}^{(j)}))}{\mu_2 R_2^2}, \quad (162)$$

with

$$2\hat{T}^{(j)} = \sum_{k=1}^{2} \frac{\hat{P}_{R_k^{(j)}}^{\dagger} \hat{P}_{R_k^{(j)}}}{\mu_k^{(j)}} + \sum_{k=1}^{2} \frac{\hat{L}_k^{(j)\dagger} \hat{L}_k^{(j)}}{\mu_k^{(j)} R_k^{(j)2}}, \quad (163)$$

where $\mu_i$ denotes the reduced mass associated with $\vec{R}_i$ and $\mu_k^{(j)}$, the reduced mass associated with $\vec{R}_k^{(j)}$. They read

$$\begin{aligned}
\mu_1^{(j)} &= \frac{2m_H m_O}{2m_H + m_O}, \quad j = 1, 2; \\
\mu_2^{(j)} &= m_H/2, \quad j = 1, 2; \\
\mu_1 &= \frac{m_H(4m_H + 2m_O)}{5m_H + 2m_O}; \\
\mu_2 &= \frac{2m_H + m_O}{2}. \quad (164)
\end{aligned}$$

Next, since it is $\vec{R}_2^{(2)}$ that plays a role in the definition of the BF frame, Eq. (127) is no longer valid. In fact, it is easy to derive an expression similar to Eq. (127) for our particular case: we would obtain something similar to Eq. (127) except that it is $\hat{L}_2^{(2)}$ which would appear on the right hand side of the scalar products instead of $\hat{L}_{n-2} (= \hat{L}_1$ for the Zundel cation). But, here, let us simply put $\hat{J} = 0$ and recast Eq. (162) as
$$2\hat{T} = \sum_{i=1}^{2} \frac{\hat{P}_{R_i}^{\dagger} \hat{P}_{R_i}}{\mu_i} + \frac{(\hat{\vec{L}}_1^{\dagger} \hat{\vec{L}}_1)_{BF}}{\mu_1 R_1^2} \\ + \sum_{j=1}^{2} \left[ \sum_{k=1}^{2} \frac{\hat{P}_{R_k^{(j)}}^{\dagger} \hat{P}_{R_k^{(j)}}}{\mu_k^{(j)}} + \frac{((\hat{\vec{J}}^{(j)})^{\dagger} - \hat{\vec{L}}_1^{(j)}) (\hat{\vec{J}}^{(j)} - \hat{\vec{L}}_1^{(j)}))_{BF}^{(j)}}{\mu_k^{(j)} R_k^{(j)} R_k^{(j)} R_k^{(j)}} + \frac{(\hat{\vec{L}}_1^{(j)})^{\dagger} (\hat{\vec{L}}_1^{(j)})_{BF}^{(j)}}{\mu_k^{(j)} R_k^{(j)} R_k^{(j)} R_k^{(j)}} \right] \\ + \frac{((\hat{\vec{L}}_1^{\dagger} + \sum_{j=1}^{2} \hat{\vec{J}}^{(j)})^{\dagger} \cdot (\hat{\vec{L}}_1 + \sum_{j=1}^{2} \hat{\vec{J}}^{(j)}))_{BF}}{\mu_2 R_2^2}. \quad (165)$$

The BF-components of $\hat{\vec{L}}_1$ are given by Eq. (75) and $BF^{BF(j)}$-components of $\hat{\vec{L}}_1^{(j)}$, with $(j=1, 2)$, and $\hat{\vec{J}}^{(1)}$ by Eqs. (122) and (124), respectively. But what about the BF- and $BF_{BF(2)}$-components of $\hat{\vec{J}}^{(2)}$? As in Section IIC 1 (iii), we start from another frame, the $E_2$ frame (we recall that the $E_2$ frame is the frame resulting from the first two Euler rotations only). Since the proton and the two molecules of water are not involved in the definition of the $E_2$ frame, the expressions of the $E_2$-components in terms of the $E_2$ coordinates of all the angular momenta appearing in Eq. (165) have a regular form. In particular, one obtains

$$\begin{bmatrix} \hat{J}_{xE_2}^{(2)} \\ \hat{J}_{yE_2}^{(2)} \\ \hat{J}_{zE_2}^{(2)} \end{bmatrix} = \begin{bmatrix} -\cos \alpha_{E_2}^{(2)} \cot \beta_{E_2}^{(2)} & -\sin \alpha_{E_2}^{(2)} & \frac{\cos \alpha_{E_2}^{(2)}}{\sin \beta_{E_2}^{(2)}} \\ -\sin \alpha_{E_2}^{(2)} \cot \beta_{E_2}^{(2)} & \cos \alpha_{E_2}^{(2)} & \frac{\sin \alpha_{E_2}^{(2)}}{\sin \beta_{E_2}^{(2)}} \\ 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} \frac{1}{i} \frac{\partial}{\partial \alpha_{E_2}^{(2)}} \\ \frac{1}{i} \frac{\partial}{\partial \beta_{E_2}^{(2)}} \\ \frac{1}{i} \frac{\partial}{\partial \gamma_{E_2}^{(2)}} \end{bmatrix}, \quad (166)$$

where $\alpha_{E_2}^{(2)}$, $\beta_{E_2}^{(2)}$, and $\gamma_{E_2}^{(2)}$ are the three Euler angles of the second molecule of water in the $E_2$ frame.

The BF-projections of $\hat{\vec{J}}^{(2)}$ are obtained by rotating the $E_2$-projections of $\hat{\vec{J}}^{(2)}$ through the angle $\alpha_{E_2}^{(2)} = \gamma$ around the $z^{E_2}$ axis

$$\begin{bmatrix} \hat{J}_{x^{BF}}^{(2)} \\ \hat{J}_{y^{BF}}^{(2)} \\ \hat{J}_{z^{BF}}^{(2)} \end{bmatrix} = \begin{bmatrix} \cos \alpha_{E_2}^{(2)} & \sin \alpha_{E_2}^{(2)} & 0 \\ -\sin \alpha_{E_2}^{(2)} & \cos \alpha_{E_2}^{(2)} & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \hat{J}_{x^{E_2}}^{(2)} \\ \hat{J}_{y^{E_2}}^{(2)} \\ \hat{J}_{z^{E_2}}^{(2)} \end{bmatrix}. \quad (167)$$
A final change of coordinates from the $E_2$ to the BF coordinates is then performed

$$
\begin{align}
\gamma &= \alpha_{E_2}^{(2)} , \tag{168} \\
\theta_1^{BF} &= \theta_1^{E_2} , \nonumber \\
\varphi_1^{BF} &= \varphi_1^{E_2} - \alpha_{E_2}^{(2)} , \nonumber \\
\alpha_{BF}^{(1)} &= \alpha_{E_2}^{(1)} - \alpha_{E_2}^{(2)} , \nonumber \\
\beta_{BF}^{(j)} &= \beta_{E_2}^{(j)}, \quad (j=1,2) , \nonumber \\
\gamma_{BF}^{(j)} &= \gamma_{E_2}^{(j)}, \quad (j=1,2) , \nonumber \\
\theta_1^{BF(j)} &= \theta_1^{BF(j)}, \quad (j=1,2) . \nonumber
\end{align}
$$

Applying the chain rule to the partial derivatives, gives

$$
\begin{align}
\frac{\partial}{\partial \alpha_{E_2}^{(2)}} &= \frac{\partial}{\partial \gamma} - \frac{\partial}{\partial \alpha_{BF}^{(1)}} - \frac{\partial}{\partial \varphi_1^{BF}} , \tag{169} \\
\frac{\partial}{\partial \theta_1^{E_2}} &= \frac{\partial}{\partial \theta_1^{BF}} , \nonumber \\
\frac{\partial}{\partial \varphi_1^{E_2}} &= \frac{\partial}{\partial \varphi_1^{BF}} , \nonumber \\
\frac{\partial}{\partial \alpha_{E_2}^{(1)}} &= \frac{\partial}{\partial \alpha_{BF}^{(1)}} , \nonumber \\
\frac{\partial}{\partial \beta_{E_2}^{(j)}} &= \frac{\partial}{\partial \beta_{BF}^{(j)}} , \quad (j=1,2) , \nonumber \\
\frac{\partial}{\partial \gamma_{E_2}^{(j)}} &= \frac{\partial}{\partial \gamma_{BF}^{(j)}} , \quad (j=1,2) , \nonumber \\
\frac{\partial}{\partial \theta_1^{BF(j)}} &= \frac{\partial}{\partial \theta_1^{BF(j)}} . \nonumber
\end{align}
$$

It is important not to forget the term $-\frac{\partial}{\partial \varphi_1^{BF}}$ in the first line of Eq. (169). Inserting Eq. (169) into Eq. (167), we obtain
$$ \begin{bmatrix} \hat{J}_{x^{BF}}^{(2)} \\ \hat{J}_{y^{BF}}^{(2)} \\ \hat{J}_{z^{BF}}^{(2)} \end{bmatrix} = \begin{bmatrix} -\cot \beta_{BF}^{(2)} & 0 & \frac{1}{\sin \beta_{BF}^{(2)}} \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} \frac{1}{i} \left( \frac{\partial}{\partial \gamma} - \frac{\partial}{\partial \alpha_{BF}^{(1)}} - \frac{\partial}{\partial \varphi_{1}^{BF}} \right) \\ \frac{1}{i} \frac{\partial}{\partial \beta_{BF}^{(2)}} \\ \frac{1}{i} \frac{\partial}{\partial \gamma_{BF}^{(2)}} \end{bmatrix}. \quad (170) $$

In addition, in view of Eqs. (69,70)

$$ \begin{bmatrix} \hat{J}_{x^{BF(2)}}^{(2)} \\ \hat{J}_{y^{BF(2)}}^{(2)} \\ \hat{J}_{z^{BF(2)}}^{(2)} \end{bmatrix} = \begin{bmatrix} -\frac{\cos \gamma^{(2)}}{\sin \beta^{(2)}} & \sin \gamma^{(2)} & \cot \beta^{(j)} & \cot \gamma^{(j)} \\ \frac{\sin \gamma^{(2)}}{\sin \beta^{(2)}} & \cos \gamma^{(2)} & -\cot \beta^{(2)} & \sin \gamma^{(2)} \\ 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} \frac{1}{i} \left( \frac{\partial}{\partial \gamma} - \frac{\partial}{\partial \alpha_{BF}^{(1)}} - \frac{\partial}{\partial \varphi_{1}^{BF}} \right) \\ \frac{1}{i} \frac{\partial}{\partial \beta^{(2)}} \\ \frac{1}{i} \frac{\partial}{\partial \gamma^{(2)}} \end{bmatrix}. \quad (171) $$

Consequently, we explicitly have all the components of the angular momenta appearing in Eq. (165). It results in a KEO with 95 terms whose explicit expression is given in the Appendix of Ref. [128]. As already explained in Section II E, the correctness of the derivation of the KEO was checked by comparing it with numerical results provided by the program TNUM of Lauvergnat and Nauts [68].

### (ii) Mixed Jacobi/valence description:

In Ref. [155], the computation of the Infra-red spectra of the isotopically substituted forms of the Zundel cation $D(D_2O)_2^+$, $H(D_2O)_2^+$, and $D(H_2O)_2^+$ was reported. For technical reasons (the cluster expansion of the potential failed to converge), it was necessary to switch from the previous Jacobi parametrization to a mixed valence/Jacobi parametrization which is depicted in Figure 19. The vectors $\vec{R}_2$ is no longer the Jacobi vector joining the center of mass of one monomer to the center of mass of the second monomer but a valence vector joining the two O atoms. The vector $\vec{R}_1$ joins the center of mass of the $O_2$ subsystem to the proton. The matrix $M$ is no longer diagonal and now reads
FIG. 19: Mixed Jacobi/valence description of the H$_5$O$_2^+$ system. G$_{O2}$ is the center of mass of the O$_2$ subsystem, or, in other words, the middle of $\vec{R}_2$, the valence vector joining one O atom to the other one.

<figure><img src="image_20.png" /></figure>

$$ M = \begin{bmatrix} \frac{2m_H+m_O}{2m_O m_H} & 0 & 0 & 0 & -\frac{1}{2m_O} & +\frac{1}{m_O} \\ 0 & \frac{2}{m_H} & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{2m_H+m_O}{2m_O m_H} & 0 & -\frac{1}{2m_O} & -\frac{1}{m_O} \\ 0 & 0 & 0 & \frac{2}{m_H} & 0 & 0 \\ -\frac{1}{2m_O} & 0 & -\frac{1}{2m_O} & 0 & \frac{m_H+2m_O}{2m_O m_H} & 0 \\ +\frac{1}{m_O} & 0 & -\frac{1}{m_O} & 0 & 0 & \frac{2}{m_O} \end{bmatrix}. \quad (172) $$

The part of the KEO coming from the diagonal part of the matrix M of Eq. (172) is identical to the operator for the Jacobi description given in the Appendix of Ref. [128] except that the reduced masses are different. With the new description, four more terms crop up in the KEO. They read
$$
\begin{align}
\hat{T}_{\vec{R}_1^{(1)} \vec{R}_2} &= -\frac{1}{m_o} \frac{(\hat{P}_{\vec{R}_1^{(1)}}^{\dagger} \cdot \hat{P}_{\vec{R}_2} + \hat{P}_{\vec{R}_2}^{\dagger} \cdot \hat{P}_{\vec{R}_1^{(1)}})}{2}, \nonumber \\
\hat{T}_{\vec{R}_1^{(2)} \vec{R}_2} &= +\frac{1}{m_o} \frac{(\hat{P}_{\vec{R}_1^{(1)}}^{\dagger} \cdot \hat{P}_{\vec{R}_2} + \hat{P}_{\vec{R}_2}^{\dagger} \cdot \hat{P}_{\vec{R}_1^{(2)}})}{2}, \nonumber \\
\hat{T}_{\vec{R}_1^{(1)} \vec{R}_1} &= -\frac{1}{2m_o} \frac{(\hat{P}_{\vec{R}_1^{(1)}}^{\dagger} \cdot \hat{P}_{\vec{R}_1} + \hat{P}_{\vec{R}_1}^{\dagger} \cdot \hat{P}_{\vec{R}_1^{(1)}})}{2}, \nonumber \\
\hat{T}_{\vec{R}_1^{(2)} \vec{R}_1} &= -\frac{1}{2m_o} \frac{(\hat{P}_{\vec{R}_1^{(2)}}^{\dagger} \cdot \hat{P}_{\vec{R}_1} + \hat{P}_{\vec{R}_1}^{\dagger} \cdot \hat{P}_{\vec{R}_1^{(2)}})}{2}. \tag{173}
\end{align}
$$

Let us recast only the first term. Since the conjugate momentum can be written as:

$$
\hat{P}_i = \hat{P}_{R_i} \vec{e}_i - \frac{\vec{e}_i \times \hat{\vec{L}}_i}{R_i}, \qquad (174)
$$

with $\vec{e}_i = \vec{R}_i/R_i$ a unit vector along $\vec{R}_i$; $\hat{P}_{R_i}$, the radial momentum and $\hat{\vec{L}}_i$, the angular momentum; and its adjoint as

$$
(\vec{P}_i)^{\dagger} = \hat{P}_{R_i}^{\dagger} \vec{e}_i + \frac{\hat{\vec{L}}_i^{\dagger} \times \vec{e}_i}{R_i}, \qquad (175)
$$

$\hat{T}_{\vec{R}_1^{(1)} \vec{R}_2}$ can be rewritten as

$$
\begin{equation}
\begin{split}
2\hat{T}_{\vec{R}_1^{(1)} \vec{R}_2} ={}& \frac{1}{m_o} [\vec{e}_1^{(1)} \cdot \vec{e}_2 (\hat{P}_{\vec{R}_1^{(1)}}^{\dagger} \hat{P}_{\vec{R}_2} + \hat{P}_{\vec{R}_2}^{\dagger} \hat{P}_{\vec{R}_1^{(1)}}) \\
& \qquad - \frac{\hat{\vec{L}}_1^{(1)}^{\dagger} \times \vec{e}_1^{(1)}}{R_1^{(1)}} \cdot \frac{\vec{e}_2 \times \hat{\vec{L}}_2}{R_2} - \frac{\hat{\vec{L}}_2^{\dagger} \times \vec{e}_2}{R_2} \cdot \frac{\vec{e}_1^{(1)} \times \hat{\vec{L}}_1^{(1)}}{R_1^{(1)}} \\
& \qquad + \hat{P}_{R_2} \frac{\hat{\vec{L}}_1^{(1)}^{\dagger} \times \vec{e}_1^{(1)}}{R_1^{(1)}} \cdot \vec{e}_2 - \hat{P}_{R_2}^{\dagger} \vec{e}_2 \cdot \frac{\vec{e}_1^{(1)} \times \hat{\vec{L}}_1^{(1)}}{R_1^{(1)}} \\
& \qquad - \hat{P}_{R_1^{(1)}}^{\dagger} \vec{e}_1^{(1)} \cdot \frac{\vec{e}_2 \times \hat{\vec{L}}_2}{R_2} + \hat{P}_{R_1^{(1)}} \frac{\hat{\vec{L}}_2^{\dagger} \times \vec{e}_2}{R_2} \cdot \vec{e}_1^{(1)}].
\end{split}
\tag{176}
\end{equation}
$$

As already mentioned at the end of Section II D, these terms raise a new technical problem
since the expression of the BF projections of $\hat{\vec{L}}_1^{(1)}$ and of the BF$^{(1)}$ projections of $\vec{L}_2$ are
complicated. Obviously, one of these two frames must be chosen to calculate the scalar
products appearing in Eq. (176), $-\frac{\hat{L}_1^{(1)\dagger} \times e_1^{(1)}}{R_1^{(1)}} \cdot \frac{e_2 \times \hat{L}_2}{R_2}$, for instance. If one chooses the BF
frame, equations similar (with $\alpha^{(1)}$, $\beta^{(1)}$, and $\gamma^{(1)}$ instead of $\alpha$, $\beta$, and $\gamma$) to Eqs. (69,70)
along with (123) must be used to obtain the BF components of $\hat{L}_1^{(1)}$. The appearance of
the three Euler matrices of $\alpha^{(1)}$, $\beta^{(1)}$, and $\gamma^{(1)}$ complicates the final expression. However, all
these terms can be calculated.

## 2. Protonated water tetramer (Eigen cation)

For H<sub>9</sub>O<sub>4</sub><sup>+</sup>, the Eigen cation, a judicious vector parametrization is the mixed Jacobi/Radau one as shown in Figure 20. In Figure 20, the system is separated into four subsystems: the three molecules of water and the H<sub>3</sub>O<sup>+</sup> ion. Each molecule of water is parametrized by two Jacobi vectors: R&#773;<sub>1</sub><sup>(j)</sup>, the vector connecting the center of mass of the H<sub>2</sub> subsystem to the O atom, and R&#773;<sub>2</sub><sup>(j)</sup>, the vector joining one H atom to the other one (here, j = 1, 2, 3). The H<sub>3</sub>O<sup>+</sup> system is parametrized by three Radau vectors, joining B<sub>OH3</sub>, the canonical point of H<sub>3</sub>O<sup>+</sup>, to the H atoms. These three Radau are denoted R<sub>i</sub><sup>(4)</sup>, with (i = 1, 2, 3), in Figure 20.

By definition, $(B_{OH3}G_{H3})^2 = G_{OH3}G_{H3} \times OG_{OH3}$, with $G_{H3}$, the center of mass of the H<sub>3</sub> subsystem; $G_{OH3}$, the center of mass of H<sub>3</sub>O<sup>+</sup>. Now, the three 'remaining' atoms are the three 'Radau' vectors $\overrightarrow{BG_j}$, with (j = 1, 2, 3). Here, B is the canonical point of the Eigen cation and $G_j$ is the center of mass of the j-th molecule of water. By definition, $(BG_{H6O3})^2 = GG_{H6O3} \times G_{H3O}G_{H6O3}$, with G, the center of mass of the Eigen cation; $G_{H6O3}$, the center of mass of the water trimer, i.e. the Eigen cation without the H<sub>3</sub>O<sup>+</sup> ion. The KEO can be straightforwardly obtained from the formulation of Section II D with n=4, $N^{(j)}$ = 3 with j = 1, 2, and 3 and $N^{(4)}$ = 4. To be more specific, in view of Eqs. (127,129) this reads
<figure><img src="image_21.png" /><figcaption>FIG. 20: Mixed Jacobi/Radau description of the H$_9$O$_4^+$ system (see the text).</figcaption></figure>

$$
\begin{equation}
\begin{aligned}
\hat{T} = & \sum_{i=1}^{3} \left( -\frac{1}{2\mu_i} \frac{1}{R_i} \frac{\partial^2}{\partial R_i^2} R_i \right) + \sum_{i=1}^{2} \left( \frac{1}{2\mu_3 R_3^2} + \frac{1}{2\mu_i R_i^2} \right) \hat{\bar{L}}_i^\dagger \hat{\bar{L}}_i + \sum_{j=1}^{4} \hat{T}^{(j)} \\
& + \sum_{i,j=1}^{2} \frac{(\hat{\bar{L}}_i \hat{\bar{L}}_j)_{BF}}{\mu_3 R_3^2} - \sum_{i=1}^{2} \frac{(\hat{\bar{J}} \hat{\bar{L}}_i)_{BF}}{\mu_3 R_3^2} + \frac{\hat{\bar{J}}_{BF}^2}{2\mu_3 R_3^2} \\
& + \sum_{i=1,...,2}^{j=1,...,4} \frac{(\hat{\bar{J}}^{(j)} \hat{\bar{L}}_i)_{BF}}{\mu_3 R_3^2} + \sum_{l,m=1}^{4} \frac{(\hat{\bar{J}}^{(l)} \hat{\bar{J}}^{(m)})_{BF}}{\mu_3 R_3^2} - \sum_{j=1}^{4} \frac{(\hat{\bar{J}} \hat{\bar{J}}^{(j)})_{BF}}{\mu_3 R_3^2},
\end{aligned}
\tag{177}
\end{equation}
$$

and
$$
\begin{equation}
\begin{split}
\hat{T}^{(j)} = \sum_{k=1}^{2} & \left( -\frac{1}{2\mu_k^{(j)}} \frac{1}{R_k^{(j)}} \frac{\partial^2}{\partial R_k^{(j)^2}} R_k^{(j)} \right) + \left( \frac{1}{2\mu_2^{(j)}} \frac{1}{R_2^{(j)^2}} + \frac{1}{2\mu_1^{(j)}} \frac{1}{R_1^{(j)^2}} \right) \hat{L}_1^{(j)\dagger} \hat{L}_1^{(j)} \right. \\
& \left. - \frac{\hat{J}_{z^{BF}(j)}^{(j)2}}{\mu_2^{(j)} R_2^{(j)^2}} - \frac{(\hat{J}_{+^{BF}(j)}^{(j)} \hat{L}_{1-^{BF}(j)}^{(j)} + \hat{J}_{-^{BF}(j)}^{(j)} \hat{L}_{1+^{BF}(j)}^{(j)})}{2\mu_2^{(j)} R_2^{(j)^2}} + \frac{\hat{J}_{BF(j)}^{(j)2}}{2\mu_2^{(j)} R_2^{(j)^2}},
\end{split}
\tag{178}
\end{equation}
$$

and

$$
\begin{equation}
\begin{split}
\hat{T}^{(4)} = \sum_{k=1}^{3} & \left( -\frac{1}{2\mu_k^{(4)}} \frac{1}{R_k^{(4)}} \frac{\partial^2}{\partial R_k^{(4)^2}} R_k^{(4)} \right) + \sum_{k=1}^{2} \left( \frac{1}{2\mu_3^{(4)} R_3^{(4)^2}} + \frac{1}{2\mu_k^{(4)} R_k^{(4)^2}} \right) \hat{L}_k^{(4)\dagger} \hat{L}_k^{(4)} \\
& - \frac{\hat{J}_{z^{BF}(4)}^{(4)} (\hat{J}_{z^{BF}(4)}^{(4)} - \hat{L}_{1z^{BF}(4)}^{(4)}) + (\hat{L}_{1z^{BF}(4)}^{(4)})^2}{\mu_3^{(4)} R_3^{(4)^2}} \\
& + \frac{(\hat{L}_{1+^{BF}(4)}^{(4)} \hat{L}_{2-^{BF}(4)}^{(4)} + \hat{L}_{1-^{BF}(4)}^{(4)} \hat{L}_{2+^{BF}(4)}^{(4)})}{2\mu_3^{(4)} R_3^{(4)^2}} \\
& - \frac{\sum_{k=1}^{2} (\hat{J}_{+^{BF}(4)}^{(4)} \hat{L}_{k-^{BF}(4)}^{(4)} + \hat{J}_{-^{BF}(4)}^{(4)} \hat{L}_{k+^{BF}(4)}^{(4)})}{2\mu_3^{(4)} R_3^{(4)^2}} + \frac{\hat{J}_{BF(4)}^{(4)2}}{2\mu_3^{(4)} R_3^{(4)^2}},
\end{split}
\tag{179}
\end{equation}
$$

with the corresponding reduced masses:

$$
\begin{align}
\mu_1^{(j)} &= \frac{2m_H m_O}{2m_H + m_O}, \quad j = 1, 2, 3; \notag \\
\mu_2^{(j)} &= m_H/2, \quad j = 1, 2, 3; \notag \\
\mu_i^{(4)} &= m_H, \quad i = 1, 2, 3; \notag \\
\mu_i &= 2m_H + m_O, \quad i = 1, 2, 3. \tag{180}
\end{align}
$$
IV. DESCRIPTION OF LARGE SYSTEMS WITH ACTIVE POLYSPHERICAL
COORDINATES

The rigorous quantum mechanical description of the dynamics of polyatomic molecular
systems remains a formidable challenge for the theoretical chemists [156, 157]. Indeed, the
cost of quantum dynamics simulations grows very rapidly with the number of degrees of
freedom involved. For the numerical solution of the Schrödinger equation of the nuclei,
very important developments have been recently introduced to describe the motion of
molecular systems with 10-30 degrees of freedom: for instance, the MCTDH, CI-VSCF,
WOSA approaches already mentioned in the introduction. More recently, the multilayer
formulation of the MCTDH algorithm [11, 12, 158–161] seems promising for the treatment
of hundreds of particles. However, the solution of the Schrödinger equation of the nuclei will
still be a challenging task, in particular due to the complexity of the Hamiltonian operator.

Another bottleneck in the future treatment of large systems will be the need for potential
energy surfaces (the potential energy surfaces of the ground and of the excited states).
In principle, since the potential energy surfaces for an *N*-atom system depend on 3*N*-6
coordinates and about 10 points are needed on the surfaces for each degree of freedom in
order to characterize them, about 10<sup>3*N*−6</sup> *ab initio* quantum chemistry calculations should
be required for each electronic state. However, important progress have been made in the
surface-fitting procedures as well, which allows to describe PESs with a limited number of
*ab initio* calculations. For instance, we mention the Sheppard interpolation scheme and
the fitting with redundant coordinates [162–171]. However, the generation of an accurate
multi-dimensional PES remains a difficult task.

Fortunately, dynamical processes in complex molecular systems often evolve around
a 'system' or 'reactive core' composed of relatively few 'active' degrees of freedom that
carry the dominant part of the dynamics. The latter are generally coupled to a 'bath'
composed of a much larger set of degrees of freedom that act as an environment. Under
these conditions, it is reasonable to invoke simplifying approximations for the other degrees
of freedom when calculating the potential and simulating dynamics. These approximations
permit the construction of simple model Hamiltonians with a relatively small number of
calculations for the electronic energy. A widely used approach is the rigid-constraint one, which consists in freezing some bond lengths, angles or entire atomic groups: the dynamics simulations are performed within a subspace of 'active' coordinates only and, for each set of active coordinates, the potential is calculated by relaxing all remaining coordinates according to a local potential minimization (see Refs. [68, 69], and references therein for applications).

An improvement of the previous model can also be envisioned by means of the reaction path Hamiltonian approach [42, 43] that makes use of a harmonic approximation locally defined along one or several reaction or 'active' coordinates (see also [21, 69, 172] for other models related to reaction path Hamiltonian). In its simplest version, this approach becomes a model composed of one or several active coordinates simply coupled to a bath [43]. With the reaction path Hamiltonian approach, many degrees of freedom can be taken into account in the dynamics, using MCTDH or its multilayer formulation for instance, whereas the number of *ab initio* quantum chemistry calculations does not dramatically increase (see for instance [173-175] for applications).

Another fruitful model has been formulated [46-48] to describe multimode vibronic coupling effects. This model provides the first terms in a Taylor expansion of a diabatic representation of the potential surfaces around the Franck-Condon point and can correctly reproduce the nonadiabatic effects in its vicinity. Thanks to this model, even though the number of degrees of freedom in the dynamics can be relatively large, the number of *ab initio* quantum chemistry calculations and the number of parameters in the analytical expression of the potential are rather limited [176-178]. More recently [179, 180], the model has been generalized to study the quantum dynamics through a conical intersection in a macro-system. The system is also split into a 'system' that collects the most important 'active' coordinates and a 'bath'. The bath modes are transformed to three 'effective' modes only, and it has been shown that the effective modes reproduce the short-time dynamics of the macro-system correctly.

These two models, the reaction path Hamiltonian and the vibronic coupling model, demonstrate that it is possible to correctly reproduce elementary processes in chemical dynamics with a relatively modest number of calculations for the electronic energy and for
dynamics. Obviously, the choice of the active coordinates is crucial. More precisely, only a judicious choice of active coordinates allows to carry out the dominant part of the dynamics using a small number of degrees of freedom.

Now, rectilinear coordinates, such as the Cartesian or normal-mode coordinates, may be
unsuitable since a huge number of them is necessary to describe chemical processes when
motions of large amplitude are involved. Consequently, *curvilinear coordinates are sometimes the only coordinates that can correctly describe these chemical processes with a limited number of degrees of freedom*. The importance of curvilinear coordinates to describe elementary processes in chemistry is well known and has already been highlighted by numerous publications: see, for instance, Ref. [51, 52] for the calculation of thermal rate constants of chemical reactions or Ref. [55, 56] for photochemical reactivity involving conical intersections. The number of quantum mechanical simulations with curvilinear coordinates for these processes is however rather limited. This difficulty probably lies in obtaining the correct kinetic energy operator in a form that is not too involved and well-adapted to the numerical methods used to solve the Schrödinger equation. The goal of the present Part of the present review is (i) to present a rigorous theoretical background to derive kinetic energy operators when only a limited set of 'active' coordinates is considered to simulate the dynamics and, (ii) to show that the polyspherical approach provides an ideal framework to obtain these operators in a form that is suitable to numerical applications.

**A. Kinetic energy operators for rigidly and adiabatically constrained molecular**
**systems**

Let the *N*-atom system be subject to constraints, so that its motion is no longer
free. The purpose of the additional constraints is to reduce the dimensionality of the
problem, from 3*N*-6 (free system) internal degrees of freedom to *n* < 3*N*-6 internal degrees
of freedom by means of *m* = 3*N*-*n*-6 additional constraints. Two cases must be distinguished:

(1) The rigid-constraint approach: it consists in freezing some bond lengths or angles.

(2) The adiabatic-constraint approach: it consists in adjusting the variations of the inactive coordinates to those of the active coordinates.
The rigorous derivation of the corresponding constrained kinetic energy operators is a rather delicate task [68–70, 72, 98, 181–187] and is reviewed in the present Section. The results presented below (for instance Eqs. (190,211)) do not depend on the choice of the internal coordinates, and, consequently, they are not just valid for polyspherical coordinates.

Let the 3N-6 internal degrees of freedom {$q_i; i = 1, \dots, 3N-6$} be collectively denoted $\mathbf{q}$. The starting point in applying constraints to a molecular system to constraints consists in partitioning the set of internal degrees of freedom into two subsets: $\mathbf{q}'$ and $\mathbf{q}''$. $\mathbf{q}'$ corresponds to the set of 'active' coordinates, {$q'_i; i = 1, \dots, n$}, and $\mathbf{q}''$ corresponds to the set of 'inactive' coordinates, {$q''_i; i = 1, \dots, m$}, with $n+m = 3N-6$. In addition, $\hat{\mathbf{P}}'$, $\hat{\mathbf{P}}''$, $\hat{\mathbf{J}}$ correspond to {$\hat{P}_{q'_i} = \frac{1}{i} \frac{\partial}{\partial q'_i}; i = 1, \dots, n$}, {$\hat{P}_{q''_i} = \frac{1}{i} \frac{\partial}{\partial q''_i}; i = 1, \dots, m$}, and {$\hat{J}_{x^{BF}}, \hat{J}_{y^{BF}}, \hat{J}_{z^{BF}}$}, respectively (see Eq. (77) for the expression of $\hat{J}_{x^{BF}}, \hat{J}_{y^{BF}}$, and $\hat{J}_{z^{BF}}$ in terms of the conjugate momenta associated with the Euler angles). After this separation into active and inactive coordinates, the KEO can be written as follows

$$ 
\begin{aligned}
2\hat{T} &= \left[ \hat{\mathbf{P}}^{\dagger} \hat{\mathbf{P}}^{\dagger} \hat{\mathbf{J}} \right] \Xi \left[ \begin{array}{c} \hat{\mathbf{P}}^{T} \\ \hat{\mathbf{P}}^{\prime\prime T} \\ \hat{\mathbf{J}}^{T} \end{array} \right] \\
&= \left[ \hat{\mathbf{P}}^{\dagger} \hat{\mathbf{P}}^{\prime\prime} \hat{\mathbf{J}} \right] \left[ \begin{array}{ccc} \Sigma' & \Omega^T & \sigma'^T \\ \Omega & \Sigma'' & \sigma''^T \\ \sigma' & \sigma'' & \Gamma \end{array} \right] \left[ \begin{array}{c} \hat{\mathbf{P}}^{T} \\ \hat{\mathbf{P}}^{\prime\prime T} \\ \hat{\mathbf{J}}^{T} \end{array} \right].
\end{aligned} 
\tag{181} $$

The matrix $\Omega$, appearing in Eq. (181), corresponds to the coupling between active and inactive coordinates. In the particular case of standard polyspherical coordinates, all the matrix elements of Eq. (181) could be obtained from the matrices $G$, $\Gamma$, and $C$, which are given by Eqs. (C4,C5,C6), along with equations similar to Eq. (107). However, it should be made clear that the KEO in Eq. (181) can be expressed in terms of any set of coordinates and the results presented below are valid whatever the set of coordinates.
1. Rigid-constraint approach

If the inactive coordinates are fixed once and for all, the $m = 3N - 6 - n$ rigid constraints
can be expressed as

$$
q''_{i} = q''_{i}|_{o}; \quad i = 1, \dots, m; \tag{182}
$$

or equivalently

$$
\dot{q}_i'' = 0; \quad i = 1, \dots, m; \tag{183}
$$

where the dots denote the time-derivatives.

Two important points must be addressed at this level:

(1) Eq. (183) is not equivalent to

$$
\hat{P}_{q_i''} = 0; \quad i = 1, \dots, m; \tag{184}
$$

In other words, starting from the constraint-free KEO of Eq. (181) and applying the rules
($\hat{P}_{q_i''} = 0; i = 1, \dots, m$) does not lead to the correct constrained KEO.

(2) Regarding the overall rotation of the system, two cases must be distinguished. These
two cases correspond to two different kinetic energy operators and to two different experi-
mental conditions:

(i) the total angular momentum $J$ is constant. It is the subject of the present section to
derive the constrained KEOs for this case.

(ii) The molecule is in a fixed orientation or, more precisely, the BF frame is in a fixed orientation, that corresponds to

$$
\dot{\alpha} = \dot{\beta} = \dot{\gamma} = 0, \tag{185}
$$

with $\alpha, \beta, \gamma$, the three Euler angles introduced in Section IIB 2. Here also Eq. (185) is not equivalent to $\hat{P}_{\alpha} = \hat{P}_{\beta} = \hat{P}_{\gamma} = 0$, and thus, is not equivalent to $J = 0$. The derivation of the constrained KEOs corresponding to this second case is not explicitly treated in this section. However, using the results presented below (Eqs. (190,211)), its derivation is almost trivial. Indeed, if the case of a molecule oriented in space is considered, the three Euler angles play the role of inactive coordinates. The corresponding KEO could be obtained simply by
recasting Eq. (181) in terms of the three conjugate momenta associated with the Euler angles ($\hat{P}_{\alpha} = \frac{1}{i} \frac{\partial}{\partial \alpha}$, $\hat{P}_{\beta} = \frac{1}{i} \frac{\partial}{\partial \beta}$, and $\hat{P}_{\gamma} = \frac{1}{i} \frac{\partial}{\partial \gamma}$) instead of {$\hat{J}_{xBF}$, $\hat{J}_{yBF}$, $\hat{J}_{zBF}$}. The transformation is straightforward and given by Eq. (77) of Section II C 1. After this transformation, the three Euler angles can then simply be added to the group of the inactive coordinates $q''$, which now contains $m+3$ coordinates and the new KEO reads

$$2\hat{T} = \left[ \hat{P}'^{\dagger} \hat{P}''^{\dagger} \right] \left[ \begin{array}{cc} \Sigma' & \Omega^T \\ \Omega & \Sigma'' \end{array} \right] \left[ \begin{array}{c} \hat{P}'^T \\ \hat{P}''^T \end{array} \right]. \qquad (186)$$

From now on, we restrict ourselves to considering the case (i) only, i.e. the Euler angles never play the role of inactive coordinates and we start from the KEO given by Eq. (181). In order to derive the constrained kinetic energy operator, it is necessary to come back to the general definition of the conjugate momenta and angular momentum components without constraints [98]:

$$\left[ \begin{array}{c} \dot{\mathbf{q}}'^T \\ \dot{\mathbf{q}}''^T \\ \omega^T \end{array} \right] = \left[ \begin{array}{ccc} \Sigma' & \Omega^T & \sigma'^T \\ \Omega & \Sigma'' & \sigma''^T \\ \sigma' & \sigma'' & \Gamma \end{array} \right] \left[ \begin{array}{c} \hat{\mathbf{P}}'^T \\ \hat{\mathbf{P}}''^T \\ \hat{\mathbf{J}}^T \end{array} \right], \qquad (187)$$

where $\omega$ denotes {$\omega_{xBF}$, $\omega_{yBF}$, $\omega_{zBF}$}, i.e. the angular velocity vector that specifies the angular speed and axis about which a framework is rotating. Let $\hat{P}'|_o$, $\hat{P}''|_o$, $\hat{J}|_o$ be the (conjugate and angular) momenta with the extra constraints of Eq. (182): the mark $|_o$ indicates that the system being subject to rigid constraints. It is important to point out that, in general, $\hat{P}'|_o \neq \hat{P}'$ and $\hat{P}''|_o \neq 0$. Imposing the constraints of Eq. (182), i.e. $\dot{\mathbf{q}}'' = \mathbf{0}$, in Eq. (187) results in [181]

$$\Omega|_o \hat{P}'|_o^T + \Sigma''|_o \hat{P}''|_o^T + \sigma''|_o^T \hat{J}|_o^T = 0, \qquad (188)$$

where for instance: $\Sigma''|_o = \Sigma''(\mathbf{q}'; \mathbf{q}''|_o) = \Sigma''(q'_i, i = 1, \dots, n; q''_j = q''_j|_o, j = 1, \dots, m)$. In other words, the inactive coordinates are fixed in $\Sigma''|_o$. The same notation holds, of course, for all the matrices. Eq. (188) yields

$$\hat{P}''|_o^T = -\Sigma''|_o^{-1}(\Omega|_o \hat{P}'|_o^T + \sigma''|_o^T \hat{J}|_o^T). \qquad (189)$$
In view of Eq. (189), it appears that Eq. (183) is not equivalent to Eq. (184). Finally, inserting Eq. (189) into Eq. (181) gives

$$2\hat{T}|_o = [\hat{P}'|_o^\dagger \hat{J}|_o] \xi|_o \begin{bmatrix} \hat{P}'|_o^T \\ \hat{J}|_o^T \end{bmatrix} = [\hat{P}'|_o^\dagger \hat{J}|_o] \begin{bmatrix} (\Sigma'|_o - \Omega|_o^T \Sigma''|_o^{-1} \Omega|_o) & (\sigma'|_o^T - \Omega|_o^T \Sigma''|_o^{-1} \sigma''|_o^T) \\ (\sigma'|_o - \sigma''|_o \Sigma''|_o^{-1} \Omega|_o) & (\Gamma|_o - \sigma''|_o \Sigma''|_o^{-1} \sigma''|_o^T) \end{bmatrix} \begin{bmatrix} \hat{P}'|_o^T \\ \hat{J}|_o^T \end{bmatrix}. \quad (190)$$

In order to complete the derivation of the quantum-mechanical kinetic energy operator $\hat{T}|_o$, we also have to calculate the Jacobian of the transformation from the constrained coordinates $\{q'; q''|_o\}$ to the Cartesian coordinates. The Jacobian is equal by construction, to $|Det[\xi|_o]|^{-\frac{1}{2}}$, where $\xi|_o$ is the matrix appearing in Eq. (190). It is useful to express this Jacobian in terms of the Jacobian of the 'free' system:

$$J_e^{free} = |Det[\Xi]|^{-\frac{1}{2}}, \quad (191)$$

i.e. the Jacobian of the transformation from the coordinates $\{q'; q''\}$ to the Cartesian co-ordinates ($\Xi$ is the matrix appearing in Eq. (181)). It can be shown (see Section 2 in Ref. [72]) that

$$Det[\Xi] = Det[\Sigma''] Det[\xi], \quad (192)$$

where

$$\xi = \begin{bmatrix} (\Sigma' - \Omega^T \Sigma''^{-1} \Omega) & (\sigma'^T - \Omega^T \Sigma''^{-1} \sigma''^T) \\ (\sigma' - \sigma'' \Sigma''^{-1} \Omega) & (\Gamma - \sigma'' \Sigma''^{-1} \sigma''^T) \end{bmatrix}. \quad (193)$$

By introducing the constraints $(q''_i = q''_i|_o; i = 1, \dots, m)$, one obtains

$$Det[\xi|_o] = Det[\Xi|_o] Det[\Sigma''|_o]^{-1}, \quad (194)$$

so that, for the Jacobian

$$J_e^{con}|_o(q') = J_e^{free}|_o(q') \sqrt{|Det[\Sigma''|_o]|}. \quad (195)$$
In other words, the kinetic energy operator of Eq. (190) is correct only for the case in which the wave functions are normalized as follows:

$$ \int \Psi^* \Psi d\tau = 1, \qquad (196) $$

where

$$ d\tau = J_e^{\text{con}}|_o(\mathbf{q}') d\mathbf{q}' d\alpha \sin\beta d\beta d\gamma, \qquad (197) $$

with $d\mathbf{q}' = dq'_1 \dots dq'_n$ and $\alpha, \beta, \text{and } \gamma$ denote the three Euler angles. However, since $\sqrt{|Det[\Sigma''|_o]|}$ is, in general, complicated, it is often profitable to normalize the wave functions as follows:

$$ \int \Psi^* \Psi d\tau' = 1, \qquad (198) $$

where

$$ d\tau' = d\mathbf{q}' d\alpha \sin\beta d\beta d\gamma, \qquad (199) $$

then an extra-potential (i.e. purely multiplicative) term pops up (see Eq. (17)):

$$ V^{\text{con}} = (J_e^{\text{con}}|_o^{\frac{1}{2}} \hat{T}|_o J_e^{\text{con}}|_o^{-\frac{1}{2}}). \qquad (200) $$

In Eq. (200) $\hat{T}|_o$ does not operate beyond the parentheses. In other words, the KEO $\hat{T}|_o + V^{\text{con}}$, with $\hat{T}|_o$ coming from Eq. (190), $V^{\text{con}}$ coming from Eq. (200), must be used along with Eqs. (198,199). Note also that, with the volume element of Eq. (199), $\hat{P}'|_o^\dagger = \hat{P}'|_o$.
## 2. Adiabatic-constraint approach

An improvement on the rigid-constraint approach consists in adjusting the inactive co-ordinates to the active coordinates, i.e. Eq. (182) is replaced by

$$q''_{i} = q''_{i}|_{ad}(q'); i = 1, \dots, m; \qquad (201)$$

which can be collectively rewritten as $q'' = q''|_{ad}(q')$ (see Eq. (236) for example). Let us first perform a new change of coordinates from $\{q'; q''\}$ to $\{q'; f\}$, with

$$f = q'' - q''|_{ad}(q'). \qquad (202)$$

It is clear that Eq. (201) is equivalent to

$$f = 0, \qquad (203)$$

i.e. the adiabatic constraints of Eq. (201) are equivalent to rigid constraints on the 'new' set of coordinates, $\{q'; f\}$: we can thus apply the results of Section IV A 1. However, we first have to perform the change of coordinates from $\{q'; q''\}$ to $\{q'; f\}$. Let $\{\hat{P}_{fi} = \frac{1}{i} \frac{\partial}{\partial f_i}; i = 1, \dots, m\}$ be collectively denoted $\hat{P}_f$. The KEO of Eq. (181) can be rewritten as

$$2\hat{T} = \left[ \hat{P}'^\dagger \hat{P}_f^\dagger \hat{J} \right] \Xi' \begin{bmatrix} \hat{P}'^T \\ \hat{P}_f^T \\ \hat{J}^T \end{bmatrix} \\ = \left[ \hat{P}'^\dagger \hat{P}_f^\dagger \hat{J} \right] \begin{bmatrix} \Sigma_{new}' & \Omega_{new}^T & \sigma_{new}'^T \\ \Omega_{new} & \Sigma_{new}'' & \sigma_{new}''^T \\ \sigma_{new}' & \sigma_{new}'' & \Gamma_{new} \end{bmatrix} \begin{bmatrix} \hat{P}'^T \\ \hat{P}_f^T \\ \hat{J}^T \end{bmatrix}, \qquad (204)$$

with
$$ 
\begin{bmatrix}
\Sigma'_{new} & \Omega_{new}^T & \sigma'_{new}^T \\
\Omega_{new} & \Sigma''_{new} & \sigma''_{new}^T \\
\sigma'_{new} & \sigma''_{new} & \Gamma_{new}
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & 0 \\
D & 1 & 0 \\
0 & 0 & 1
\end{bmatrix}
\times
\begin{bmatrix}
\Sigma' & \Omega^T & \sigma'^T \\
\Omega & \Sigma'' & \sigma''^T \\
\sigma' & \sigma'' & \Gamma
\end{bmatrix}
\begin{bmatrix}
1 & D^T & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{bmatrix},
\quad (205)
$$

where 1 and 0 denote the identity matrix and the zero matrix, respectively. Moreover,

$$ D_{ij} = - \frac{\partial q''_j|_{eq}(\boldsymbol{q}')}{\partial q'_i}. \quad (206) $$

In other words, the matrix

$$ d = \begin{bmatrix} 1 & 0 & 0 \\ D & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad (207) $$

is the Jacobian matrix of the transformation from {$q'; q''; \alpha, \beta, \gamma$} to {$q'; f; \alpha, \beta, \gamma$}. If the motion of the inactive coordinates is decoupled from the active coordinates, all the terms $-\frac{\partial q''_j|_{ad}(q')}{\partial q'_i}$ must be very small. This is why, we proposed, in Ref. [185], to assume that, in Eq. (207),

$$ -\frac{\partial q''_j|_{ad}(\boldsymbol{q}')}{\partial q'_i} = 0; \quad (208) $$

i.e.

$$ d = 1, \quad (209) $$

or

$$ \Xi' = \Xi, \quad (210) $$

instead of $\Xi' = d \Xi d^T$. After this simplification, in order to derive the KEO, we simply have to apply the results of Section IV A 1. It follows that
$$2\hat{T}|_{ad} = [\hat{\mathbf{P}}'_{|ad} \hat{\mathbf{J}}|_{ad}] \xi|_{ad} \begin{bmatrix} \hat{\mathbf{P}}'_{|ad} \\ \hat{\mathbf{J}}'_{|ad} \end{bmatrix} \\ = [\hat{\mathbf{P}}'_{|ad} \hat{\mathbf{J}}|_{ad}] \begin{bmatrix} (\Sigma'_{|ad} - \Omega|_{ad}^T \Sigma''_{|ad}^{-1} \Omega|_{ad}) & (\sigma'_{|ad}^T - \Omega|_{ad}^T \Sigma''_{|ad}^{-1} \sigma''_{|ad}^T) \\ (\sigma'_{|ad} - \sigma''_{|ad} \Sigma''_{|ad}^{-1} \Omega|_{ad}) & (\Gamma|_{ad} - \sigma''_{|ad} \Sigma''_{|ad}^{-1} \sigma''_{|ad}^T) \end{bmatrix} \begin{bmatrix} \hat{\mathbf{P}}'_{|ad} \\ \hat{\mathbf{J}}'_{|ad} \end{bmatrix}, \qquad (211)$$

where for instance:

$$\Sigma''|_{ad} = \Sigma''(\mathbf{q}'; \mathbf{q}''|_{ad}(\mathbf{q}')) = \Sigma''(q'_i, i = 1, \dots, n; q''_j = q''_j|_{ad}(\mathbf{q}')), j = 1, \dots, m). \quad (212)$$

The same notation holds, of course, for all the matrices. In Eq. (211), we do not impose $-\frac{\partial q''_j|_{ad}(\mathbf{q}')}{\partial q'_i} = 0$, since it would lead us back to the rigid model. Paradoxically, this small (numerical) additional approximation (Eq. (210)) allows to keep a much simpler (analytical) form of the KEO. In this adiabatic model, we do not have to add some new terms to the rigid KEO (compare Eq. (190) with Eq. (211)): this adiabatic formulation allows to recast the problem into an equivalent 'flexible' pseudo-rigid case. We repeat that, from a physical point of view, this new approximation (Eq. (210)) is justified since it is well-known [188–191] that the adiabatic approximation is valid provided that there is a clear decoupling or hierarchy between the inactive coordinates and the active coordinates, i.e. provided that the functions $-\frac{\partial q''_j|_{ad}(\mathbf{q}')}{\partial q'_i}$ are close to zero.

The general equations without the additional approximation of Eq. (210) are given in Section 3 of Ref. [72]: this 'exact' adiabatic model leads to numerical results that are slightly better but the corresponding KEOs are considerably more complicated.

Let us add that the operator $\hat{T}|_{ad}$ of Eq. (211) is correct only in the case that the wave functions are normalized with an Euclidean convention of normalization, i.e. with the volume element given by

$$d\tau = J_e^{\text{con}}|_{ad}(\mathbf{q}') d\mathbf{q}' d\alpha \sin\beta d\beta d\gamma, \qquad (213)$$

with
$$J_e^{\text{con}}|_{\text{ad}}(\mathbf{q}') = J_e^{\text{free}}|_{\text{ad}}(\mathbf{q}') \sqrt{|\det[\Sigma''|_{\text{ad}}]|}. \quad (214)$$

(This result is not affected by the use of the additional approximation of Eq. (210)). However, if the wave functions are normalized according to

$$\int \Psi^* \Psi d\tau' = 1, \quad (215)$$

where

$$d\tau' = d\mathbf{q}' d\alpha \sin\beta d\beta d\gamma, \quad (216)$$

then an extra-potential (i.e. purely multiplicative) term defined as (see Eq. (17))

$$V^{\text{con}} = (J_e^{\text{con}}|_{\text{ad}}^{\frac{1}{2}} \hat{T}|_{\text{ad}} J_e^{\text{con}}|_{\text{ad}}^{-\frac{1}{2}}) \quad (217)$$

must be added to $\hat{T}|_{\text{ad}}$ of Eq. (211).

Finally, it is worth noting that the program TNUM of Lauvergnat and Nauts [68], already mentioned in Section II E, numerically computes all the functions appearing in the constrained operators. To be more specific, (i) for the rigid-constraint approach, TNUM numerically provides the matrix elements in Eq. (190) and the extra-potential term of Eq. (200); (ii) for the adiabatic-constraint approach, TNUM numerically provides the extra-potential term of Eq. (217) and the matrix elements appearing in the adiabatic operator with and without the approximation of Eq. (210).

## B. Combination with the polyspherical approach

The constrained operators of Eqs. (190,211) are a sum of two terms. The first term, called the 'main' term, reads
$$2\hat{T}_{\text{main}|o} = [\hat{\mathbf{P}}'_{|o}{}^{\dagger} \hat{\mathbf{J}}_{|o}] \begin{bmatrix} \Sigma'_{|o} & \sigma'_{|o}{}^{T} \\ \sigma'_{|o} & \Gamma_{|o} \end{bmatrix} \begin{bmatrix} \hat{\mathbf{P}}'_{|o}{}^{T} \\ \hat{\mathbf{J}}_{|o}^{T} \end{bmatrix}, \qquad (218)$$

in the rigid model and

$$2\hat{T}_{\text{main}|ad} = [\hat{\mathbf{P}}'_{|ad}{}^{\dagger} \hat{\mathbf{J}}_{|ad}] \begin{bmatrix} \Sigma'_{|ad} & \sigma'_{|ad}{}^{T} \\ \sigma'_{|ad} & \Gamma_{|ad} \end{bmatrix} \begin{bmatrix} \hat{\mathbf{P}}'_{|ad}{}^{T} \\ \hat{\mathbf{J}}_{|ad}^{T} \end{bmatrix}, \qquad (219)$$

in the adiabatic model. This 'main' term is directly obtained from the operator without constraints (Eq. (181)) simply by removing all terms involving conjugate momenta associated with the inactive coordinates. If polyspherical coordinates are used, this term can always be written as a sum of products of mono-mode operators and is directly obtained from Eq. (C4, C5, C6, 107).

The second term

$$2\hat{T}_{\text{corr}|o} = -[\hat{\mathbf{P}}'_{|o}{}^{\dagger} \hat{\mathbf{J}}_{|o}] \begin{bmatrix} \Omega_{|o}{}^{T} \Sigma''_{|o}{}^{-1} \Omega_{|o} & \Omega_{|o}{}^{T} \Sigma''_{|o}{}^{-1} \sigma''_{|o}{}^{T} \\ \sigma''_{|o} \Sigma''_{|o}{}^{-1} \Omega_{|o} & \sigma''_{|o} \Sigma''_{|o}{}^{-1} \sigma''_{|o}{}^{T} \end{bmatrix} \begin{bmatrix} \hat{\mathbf{P}}'_{|o}{}^{T} \\ \hat{\mathbf{J}}_{|o}^{T} \end{bmatrix}, \qquad (220)$$

for the rigid model, and

$$2\hat{T}_{\text{corr}|ad} = -[\hat{\mathbf{P}}'_{|ad}{}^{\dagger} \hat{\mathbf{J}}_{|ad}] \begin{bmatrix} \Omega_{|ad}{}^{T} \Sigma''_{|ad}{}^{-1} \Omega_{|ad} & \Omega_{|ad}{}^{T} \Sigma''_{|ad}{}^{-1} \sigma''_{|ad}{}^{T} \\ \sigma''_{|ad} \Sigma''_{|ad}{}^{-1} \Omega_{|ad} & \sigma''_{|ad} \Sigma''_{|ad}{}^{-1} \sigma''_{|ad}{}^{T} \end{bmatrix} \begin{bmatrix} \hat{\mathbf{P}}'_{|ad}{}^{T} \\ \hat{\mathbf{J}}_{|ad}^{T} \end{bmatrix}, \qquad (221)$$

for the adiabatic model, is a 'correction'. This term vanishes if there is no coupling between the active and the inactive coordinates, i.e. if $\sigma''$ and $\Omega$ in Eq. (181) are equal to zero. In general, if the choice of the 'active' coordinates is relevant, this 'correction' is numerically smaller than the 'main' term but may play a non-negligible role in the dynamics, as aforementioned. According to our experience, this correction is very important for accurate
ro-vibrational spectroscopy (see [4, 183, 192] for instance) but of minor importance for
the calculation of the photo-absorption spectra involving conical intersections (see [126] for
instance), which, in general, are not of high resolution.

The combination of the polyspherical approach and the previously developed rigidly and adiabatically constrained models presents several advantages. First, the polyspherical approach straightforwardly provides the 'main' term in a form adapted to the methods used to solve the Schrödinger equation. Second, as far as the 'correction' is concerned, the situation is a little bit more complex: this correction can be very complicated, and above all, this correction, even when using polyspherical coordinates, cannot be written as a sum of products of mono-mode operators, because the inverse of the matrices $\Sigma''$ will in general depend on the coordinates in a complicated non-separable way (see Eqs. (234,240) of Section IVC for examples). It is all the more important to exploit the polyspherical coordinates. Indeed, if one uses the polyspherical coordinates in orthogonal coordinates along with the separation into subsystems of Section IID, it is easy to almost separate the active coordinates from the inactive coordinates in the free KEO (see Eq. (228) of Section IVC for an example). In other words, almost all the terms in the matrices $\sigma''$ and $\Omega$ in Eq. (181) are equal to zero so that the correction almost vanishes. If the correction is neglected, the effect of the latter approximation is drastically reduced. If it is not neglected, for instance to calculate accurate infra-red spectra, it is, in general, much easier to derive $\hat{T}_{corr}$ (in polyspherical coordinates than in another set of curvilinear coordinates): see Section IVC for illustrations.

### C. Examples of applications

This Section is devoted to the presentation of different applications using the results obtained in Section IVA combined with the polyspherical approach.

#### 1. *Rigid- and Adiabatic-constraint approaches: the water dimer*

##### (i) Exact operator:

For van der Waals polymers, a decoupling between the intra- and intermolecular modes
<figure><img src="image_22.png" /><figcaption>FIG. 21: Jacobi description of the water dimer.</figcaption></figure>

is often very relevant since the intramolecular coordinates move much faster than the intermolecular ones. In Ref. [192], the determination of the water dimer potential energy surface via direct inversion of spectroscopic data was reported. An adiabatic model was adopted. A set of five Jacobi vectors was chosen. This set is depicted in Figure 21.

The dimer is split into two subsystems that are simply the two molecules of water. As far as the exact KEO is concerned, one obtains in view of Eqs. (125,126) of Section IID

$$2\hat{T} = \frac{\hat{P}_{R_1}^{\dagger} \hat{P}_{R_1}}{\mu_1} + 2 \sum_{j=1}^{2} \hat{T}^{(j)} + \frac{(\hat{\vec{J}}^{\dagger} - (\sum_{j=1}^{2} \hat{\vec{J}}^{(j)})^{\dagger}) \cdot (\hat{\vec{J}} - (\sum_{j=1}^{2} \hat{\vec{J}}^{(j)}))}{\mu_1 R_1^2}, \quad (222)$$

with

$$2\hat{T}^{(j)} = \sum_{k=1}^{2} \frac{\hat{P}_{R_k^{(j)}}^{\dagger} \hat{P}_{R_k^{(j)}}}{\mu_k^{(j)}} + \frac{((\hat{\vec{J}}^{(j)})^{\dagger} - \hat{\vec{L}}_1^{(j)}) (\hat{\vec{J}}^{(j)} - \hat{\vec{L}}_1^{(j)})_{BF^{(j)}}}{\mu_2^{(j)} R_2^{(j)2}} + \frac{(\hat{\vec{L}}_1^{(j)})^{\dagger} (\hat{\vec{L}}_1^{(j)})_{BF^{(j)}}}{\mu_1^{(j)} R_1^{(j)2}}, \quad (223)$$

where $\mu_i$ denotes the reduced mass associated with $\vec{R}_i$ and $\mu_k^{(j)}$, the reduced mass associated with $\vec{R}_k^{(j)}$. They read
$$ 
\begin{aligned}
\mu_1^{(j)} &= \frac{2m_H m_O}{2m_H + m_O}, \quad j = 1, 2; \\
\mu_2^{(j)} &= m_H/2, \quad j = 1, 2; \\
\mu_1 &= \frac{2m_H + m_O}{2}. \qquad (224)
\end{aligned}
$$

Furthermore, we have (see Eq. (124) of Section IID)

$$ 
\begin{bmatrix} \hat{L}_1 x^{BF(j)} \\ \hat{L}_1 y^{BF(j)} \\ \hat{L}_1 z^{BF(j)} \end{bmatrix} = \begin{bmatrix} -\cot \theta_1^{BF(j)} \hat{J}_{z^{BF(j)}}^{(j)} \\ -i \frac{\partial}{\partial \theta_1^{BF(j)}} \\ \hat{J}_{z^{BF(j)}}^{(j)} \end{bmatrix} . \qquad (225)
$$

The six intramolecular coordinates are the four lengths $R_1^{(1)}$, $R_2^{(1)}$, $R_1^{(2)}$, $R_2^{(2)}$, and the two angles $\theta_1^{BF(1)}$ and $\theta_1^{BF(2)}$ (see Section IID (ii) for their precise definition). The six intermolecular coordinates are $R_1$, and the Euler angles $\alpha^{(1)}, \beta^{(1)}, \gamma^{(1)}, \beta^{(2)}, \gamma^{(2)}$. $\alpha^{(2)} = \gamma$, i.e. $\alpha^{(2)}$ defines the $((x^{BF}, z^{BF}); x^{BF} > 0)$ half-plane.

(ii) Rigidly constrained operator (rigid rotor):

Before adopting an adiabatic model, as in Ref. [192], the water dimer was studied by means of a rigid monomer description for each molecule of water (see Ref. [4, 193], for instance), following the general formulation of Brocks et al. [60] for van der Waals dimers. In order to derive the corresponding rigid operator, let us go back to the exact formulation. Since we have used a polyspherical description of the water dimer with orthogonal vectors along with the separation into subsystems, it results in an exact KEO, Eqs. (222,223), in which there is almost no coupling between the active coordinates (the six intermolecular coordinates) and the inactive coordinates (the six intramolecular coordinates). More precisely, in view of Eqs. (222,223), we see that $R_1^{(1)}, R_2^{(1)}, R_1^{(2)}, R_2^{(2)}$ are not coupled with the other coordinates; that $\theta_1^{BF(1)}$ is coupled with $\hat{J}_{x^{BF(1)}}^{(1)}, \hat{J}_{y^{BF(1)}}^{(1)}, \hat{J}_{z^{BF(1)}}^{(1)}$, only; and that $\theta_1^{BF(2)}$ is coupled with $\hat{J}_{x^{BF(2)}}^{(2)}, \hat{J}_{y^{BF(2)}}^{(2)}, \hat{J}_{z^{BF(2)}}^{(2)}$, only. It means that freezing the intramolecular co-
ordinates comes down to freezing the internal coordinates of each molecule of water taken separately. For each monomer, the KEO can be written as (see Eq. (181))

$$2\hat{T}^{(j)} = [\hat{P}^{(j)\dagger} \hat{J}^{(j)}] \begin{bmatrix} \Sigma''^{(j)} & \sigma''^{(j)T} \\ \sigma''^{(j)} & \Gamma^{(j)} \end{bmatrix} \begin{bmatrix} \hat{P}^{(j)T} \\ \hat{J}^{(j)T} \end{bmatrix}, \qquad (226)$$

with $\hat{P}''^{(j)} = \{\hat{P}_{R_1^{(j)}}, \hat{P}_{R_2^{(j)}}, \hat{P}_{\theta_1^{BF(j)}}\}$ and with (we use Eqs. (223,225))

$$\Sigma''^{(j)} = \begin{bmatrix} \frac{1}{\mu_1^{(j)}} & 0 & 0 \\ 0 & \frac{1}{\mu_2^{(j)}} & 0 \\ 0 & 0 & \frac{1}{\mu_1^{(j)} R_1^{(j)2} + \mu_2^{(j)} R_2^{(j)2}} \end{bmatrix}, \qquad (227)$$

$$\sigma''^{(j)} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & \frac{1}{\mu_2^{(j)} R_2^{(j)2}} & 0 \end{bmatrix}, \qquad (228)$$

and

$$\Gamma^{(j)} = \begin{bmatrix} \frac{1}{\mu_2^{(j)} R_2^{(j)2}} & 0 & \frac{\cot \theta_1^{BF(j)}}{\mu_2^{(j)} R_2^{(j)2}} \\ 0 & \frac{1}{\mu_2^{(j)} R_2^{(j)2}} & 0 \\ \frac{\cot \theta_1^{BF(j)}}{\mu_2^{(j)} R_2^{(j)2}} & 0 & \frac{1}{\sin^2 \theta_1^{BF(j)} \mu_1^{(j)} R_1^{(j)2}} + \frac{\cot^2 \theta_1^{BF(j)}}{\mu_2^{(j)} R_2^{(j)2}} \end{bmatrix}. \qquad (229)$$

If we now subject each molecule of water ($j = 1, 2$) to the following rigid constraints:

$$R_1^{(j)} = R_1^{(j)}|_{eq},$$

$$R_2^{(j)} = R_2^{(j)}|_{eq},$$

$$\theta_1^{BF(j)} = \theta_1^{BF(j)}|_{eq} = \frac{\pi}{2}; \qquad (230)$$

where $|_{eq}$ means that the coordinate is fixed at its value at the equilibrium geometry. The 'main' term of the rigidly constrained operator reads (see Eqs. (218,230))
$$2\hat{T}_{\text{main}}^{(j)}|_o = [\hat{\mathbf{J}}|_o^{(j)}] \begin{bmatrix} \frac{1}{\mu_2^{(j)} R_2^{(j)2}}|_o & 0 & 0 \\ 0 & \frac{1}{\mu_2^{(j)} R_2^{(j)2}}|_o & 0 \\ 0 & 0 & \frac{1}{\mu_1^{(j)} R_1^{(j)2}}|_o \end{bmatrix} [\hat{\mathbf{J}}|_o^{(j)T}], \quad (231)$$

and it is worth noting that the resulting coefficients are not the principal values of the inertia tensor of the water molecule. The 'correction' reads (see Eq. (220) and Eqs. (227,228))

$$\begin{align}
2\hat{T}_{\text{corr}}^{(j)}|_o &= -[\hat{\mathbf{J}}|_o^{(j)}] \sigma''^{(j)}|_o \Sigma''^{(j)}|_o^{-1} \sigma''^{(j)}|_o^T [\hat{\mathbf{J}}|_o^{(j)T}] \nonumber \\
&= -[\hat{\mathbf{J}}|_o^{(j)}] \begin{bmatrix} 0 & 0 & 0 \\ 0 & \frac{\mu_1^{(j)} R_1^{(j)2}}{\mu_2^{(j)} R_2^{(j)2} (\mu_1^{(j)} R_1^{(j)2} + \mu_2^{(j)} R_2^{(j)2})}|_o & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} [\hat{\mathbf{J}}|_o^{(j)T}]. \tag{232}
\end{align}$$

The final rigidly constrained operator reads

$$2\hat{T}|_o = \frac{\hat{P}_{R_1}^\dagger \hat{P}_{R_1}}{\mu_1} + 2 \sum_{j=1}^2 \hat{T}^{(j)}|_o \\ + \frac{(\hat{J}^\dagger - (\sum_{j=1}^2 \hat{J}^{(j)})^\dagger) \cdot (\hat{J} - (\sum_{j=1}^2 \hat{J}^{(j)}))}{\mu_1 R_1^2}, \quad (233)$$

with

$$\begin{align}
2\hat{T}^{(j)}|_o &= 2\hat{T}_{\text{main}}^{(j)}|_o + 2\hat{T}_{\text{corr}}^{(j)}|_o \nonumber \\
&= [\hat{\mathbf{J}}|_o^{(j)}] \begin{bmatrix} \frac{1}{\mu_2^{(j)} R_2^{(j)2}}|_o & 0 & 0 \\ 0 & \frac{1}{\mu_1^{(j)} R_1^{(j)2} + \mu_2^{(j)} R_2^{(j)2}}|_o & 0 \\ 0 & 0 & \frac{1}{\mu_1^{(j)} R_1^{(j)2}}|_o \end{bmatrix} [\hat{\mathbf{J}}|_o^{(j)T}]. \tag{234}
\end{align}$$

The operator in Eq. (234) is the rigid rotor of the water molecule, i.e. the resulting coefficients are the principal values of the inertia tensor of the water molecule. The result
given by Eqs. (233,234) is a particular case of the general formulation of Brocks et al. [60] for van der Waals dimers.

(iii) Adiabatically constrained operator (flexible rotor):

In Ref. [192], it was shown that the flexibility of each water molecule plays a crucial role in the description of the hydrogen bond between the two molecules. A very relevant adiabatic formulation was adopted for the water dimer. Indeed, for the water dimer, a single quantum of excitation energy in one intramolecular mode of vibration is sufficient to break the hydrogen bond. This implies that only the ground vibrational state of each monomer has to be considered to study the rovibrational bound states of the dimer. Another possibility is to take into account the first excited vibrational state of each monomer alone if the dissociation of the dimer is to be considered. Such a situation is reminiscent of the Born-Oppenheimer approximation that allows one to confine the description of a bound molecule to one electronic state. Consequently, it turns out to be very relevant to define the adiabatic dependence of the inactive coordinates as follows. We first work out the vibrational equation for the six inactive coordinates (see Eq. (181) for the definition of $\Sigma''$):

$$ (\frac{1}{2} \dot{\mathbf{P}}''^{\dagger} \mathbf{\Sigma}'' \dot{\mathbf{P}}'' + V(\mathbf{q}''; \mathbf{q}')) \Phi_n(\mathbf{q}''; \mathbf{q}') = E_n(\mathbf{q}') \Phi_n(\mathbf{q}''; \mathbf{q}') , \quad (235) $$

with $\mathbf{q}''$ that collectively denotes $\{R_1^{(1)}, R_2^{(1)}, R_1^{(2)}, R_2^{(2)}, \theta_1^{BF(1)}, \theta_1^{BF(2)}\}$ and with $\mathbf{q}'$ that collectively denotes $\{R_1, \alpha^{(1)}, \beta^{(1)}, \gamma^{(1)}, \beta^{(2)}, \gamma^{(2)}\}$. The inactive coordinates are then adapted to the active ones depending on the intramolecular state to be considered (for instance, the ground state corresponds to n=0). In other words, Eq. (201) becomes

$$ q''_i = q''_{i|ad}(\mathbf{q}') = \langle \Phi_n(\mathbf{q}''; \mathbf{q}') | q''_i | \Phi_n(\mathbf{q}''; \mathbf{q}') \rangle ; i = 1, 2, 3 . \quad (236) $$

$q''_{i|ad}(\mathbf{q}')$ is not known analytically but numerically on a grid. For each monomer, the KEO is given by Eq. (226) and if we now subject each molecule of water ($j = 1, 2$) to the adiabatic constraints of Eq. (236), the 'main' term of the adiabatically constrained operator for each water molecule ($j = 1, 2$) reads (see Eq. (219))