# Study of vibrational stationary states of H/Pd(111)

Oussama Bindech, Fabien Gatti, Souvik Mandal, Roberto Marquardt, Jean Christophe Tremblay

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

## To cite this version:

Oussama Bindech, Fabien Gatti, Souvik Mandal, Roberto Marquardt, Jean Christophe Tremblay. Study of vibrational stationary states of H/Pd(111). Molecular Physics, 2025, 10.1080/00268976.2025.2473661. hal-05368502

HAL Id: hal-05368502

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

Submitted on 17 Nov 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Study of vibrational stationary states of H/Pd(111)

Oussama Bindech<sup>a</sup>, Fabien Gatti<sup>b</sup>, Souvik Mandal<sup>c</sup>, Roberto Marquardt<sup>a</sup>, and
Jean Christophe Tremblay<sup>c</sup>

<sup>a</sup> Laboratoire de Chimie Quantique, Institut de Chimie, UMR 7177 CNRS/Unistra
Université de Strasbourg, 4, Rue Blaise Pascal, CS 90032, 67081 Strasbourg Cedex, France;

<sup>b</sup> Institut de Sciences Moléculaires d'Orsay, UMR 8214 CNRS/Université Paris-Saclay, Bât
520, Rue André Rivière, 91405 Orsay Cedex, France; Laboratoire Interdisciplinaire Carnot de
Bourgogne, UMR 6303 CNRS/Université de Bourgogne Europe, BP 47870, 21078 Dijon,
France

<sup>c</sup> Laboratoire de Physique et Chimie Théoriques, UMR 7019 CNRS/Université de Lorraine, 1
Blvd. Arago, 57070 Metz, France

## ARTICLE HISTORY

Compiled February 12, 2025

## Abstract

We study the stationary states and energies of the hydrogen atoms when they are adsorbed on the palladium surface. Our calculations are performed with the improved relaxation method based on the Multi-Configuration Time-Dependent Hartree (MCTDH) ansatz for the solution of the Schrödinger equation. In order to simulate an infinite surface in our calculations, a periodic approach is adopted. It includes an elementary cell given by a slab of five Pd layers where each one contains 3×3 Pd atoms. Vibrational eigenstates reveal several significant quantum mechanical phenomena. In particular, we were able to show in this study the presence of tunneling effects and of a Fermi resonance coupling the localized vibrational modes of the system.

## KEYWORDS

adsorption eigenstates, hydrogen adsorption, Fermi resonance

## 1. Introduction

Phenomena involving molecule/atom-surface interactions are ubiquitous in everyday life. One may think, for instance, of the corrosion of metallic objects. In addition, metal surfaces are used as catalysts in many chemical reactions of industrial importance.

corresponding author: F. Gatti, fabien.gatti@u-bourgogne.fr

corresponding author: R. Marquardt, roberto.marquardt@unistra.fr
Consequently, a detailed understanding of the properties and interactions of atoms and molecules on the surfaces of solids is crucial. Generally, researchers in this field employ numerical modeling based on classical molecular dynamics simulations since a full quantum treatment is always challenging when many atoms are involved in the process. However, it is now well-established, among others by the work of Z. Bačić [1–4], that many chemical processes are impacted by inherent quantum mechanical phenomena involving the nuclei [5, 6, and many references cited therein].

Adopting a quantum approach is all the more essential that, on the experimental front, techniques such as Scanning Tunneling Microscopy (STM) [7], High-Resolution Electron Energy Loss Spectroscopy (HREELS) [8], neutron or helium atom scattering [9] have been developed. In principle, it is possible, with these techniques, to describe the phenomena occurring when an atom or a molecule is moving on a surface at their elementary level, including even the most subtle quantum mechanical phenomena. In this context, we mention the pioneering study of the symmetry of the frustrated translational motion of CO/Cu(100) by Bahel and Bačić [10].

The quantum description of the interaction of atoms or molecules with surfaces can be tackled either within the framework of a time-independent picture, i.e. by the assessment of the stationary states, the eigenstates of the system, or by studying the time evolution of the system starting from specific initial conditions. In this paper, we focus on the first approach. More precisely, we calculate energy eigenstates of a single hydrogen atom adsorbed on a surface of palladium atoms, here the Pd(111) surface. The study of hydrogen atoms or molecules on such a surface is important for several reasons. Pd(111) is an excellent catalyst for various chemical reactions and the description of the quantum behavior of hydrogen atoms on this surface is crucial for applications such as hydrogen storage [11]. Due to its importance, this system has given rise to theoretical studies such as in Ref. 12 based on the potential energy surface (PES) of Ref. 13 or in Ref. 14 using a variant of the PES from Ref. 15. In addition, experimental data is available from measurements using the High Resolution Electron Energy Loss Spectroscopy (HREELS) method. These measurements were carried out by Conrad *et al.* [8].

In this paper, we thus present a detailed analysis of the energies and the wave-
functions of stationary states of a hydrogen atom adsorbed on the Pd(111) surface.
The dynamics of this H/Pd(111) system involves vibrational motions of the atom in
the highly anharmonic and periodic PES of the system’s electronic ground state, and
the study focuses on the set of eigenstates containing up to two quanta of vibrational
energy in certain vibrational modes. In the present approach, palladium atoms are
considered to be rigidly constrained to their equilibrium positions in the crystal, and
the Born Oppenheimer approximation is tacitly assumed to be valid. We use the PES
developed by W. Dong et al. [15, 16] based on density functional theory (DFT) and define the position of the hydrogen atom in terms of its three Cartesian coordinates. We use the Multi-Configuration Time-dependent Hartree (MCTDH) method [17–20], more precisely its *improved relaxation* [21, 22] variant, to calculate the vibrational eigenstates on the given PES.

As shown below, the main difficulty is to converge eigenstates on large grids, some with a large degeneracy, others separated by small splittings due to different tunneling processes. Thus, the calculations must be accurate enough to clearly distinguish the different components of the diversely split tunnel multiplets generated by these effects. We will see that these energy splittings can be smaller than 0.1 meV.

Approaching the dynamics of the adsorbates while keeping the substrate atoms rigid in space is a rather strong approximation, as clearly shown by Bahel and Bačić [10]. Additional degrees of freedom related to the motion of the substrate atoms may couple to the adsorbate’s ones and potentially influence the energy position of their stationary states as well as increase the local density of states in the region of these states, leading hence to a homogeneous broadening of the spectra. Recently, the multi-layer variant of the MCTDH method was used to include the motion of 70 carbon atoms in the quantum dynamical simulation of hydrogen scattering at a graphene sheet [23]. In the present work, we intentionally concentrate on the stationary states of the adsorbates in the rigid substrate approximation, in which the dynamics proves to be already quite complex. We leave the investigation of the role of the substrate’s degrees of freedom for future studies. The MCTDH program was used in scattering studies of atoms [24] and molecules [25] at surfaces, as well as for the study of stationary states [26].

The outline of this paper is as follows. In the first section, the theoretical background of the models used as well as the computational details are discussed. In the second section, we discuss thoroughly the nature of the calculated eigenstates. The paper concludes presenting outlooks for the future.

## 2. Theoretical background

### 2.1. Potential energy surface and coordinates

As mentioned in the introduction, we use the PES developed by W. Dong et al. [15, 16] based on density functional theory (DFT). Palladium is a metal that crystallizes in the fcc (face-centred cubic) structure. In order to simulate an infinite surface in our calculations, a periodic approach is adopted, which includes an elementary cell given by a slab of five Pd layers where each one contains 3×3 Pd atoms. The vacuum space corresponds to five Pd(111) layers. The elementary cell used in our calculations is the same as the one used by the authors to develop the PES. Figure 1 depicts a scheme of
this grid showing only the surface plane (solid line) and the one below (dashed line).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** (3×3)-grid used for periodic calculations. hcp (filled blue circles) and fcc (red circumferences) denote the hexagonal closest packed and face-centered cubic sites, respectively. The larger circles correspond to the Pd atoms. *d* is the distance between two neighbour atoms, *d* = a/√2 = 275.114 pm, with *a* the lattice parameter and *a* = 389 pm. Panel (a) gives a general overview of the grid while indicating the first sublayer of palladium atoms and the discrimination of hcp and fcc sites. Panels (b) and (c) indicate relevant geometrical parameters as well as the twisted (b) and Cartesian (c) coordinate system used to represent the grid. The fcc to hcp distance is 2h = d/√3 ≈ 159 pm.</figcaption></figure>

As often with a hexagonal structure, we use "twisted" coordinates $x_t, y_t$ and $z_t$ instead of the usual Cartesian coordinates $x_c, y_c$ and $z_c$, as shown in Fig. 1(b), to carry out the MCTDH calculations. The origin, and thus $z_t = 0$, is located in the center of the Pd atom located in the middle of the topmost layer of the elementary cell.

The transformation between the two sets of coordinates is given by :

$$ 
\left.
\begin{aligned}
x_t &= x_c + \frac{1}{\sqrt{3}} y_c \\
y_t &= \frac{2}{\sqrt{3}} y_c \\
z_t &= z_c
\end{aligned}
\right\} \qquad (1)
$$
This grid contains 9 equivalent hcp sites and 9 equivalent fcc sites. The latter cor-
respond to the most stable adsorption sites at the surface level (the electronic energy
in this site is 0 eV with our conventions). The nine fcc geometries correspond to the
values $x_t = d/3 + k \times d, y_t = -d/3 + k \times d (k = -1,0,1)$ and $z_t = 90$ pm. The
nine hcp geometries correspond to the values $x_t = -d/3 + k \times d, y_t = d/3 + k \times d,$
$(k = -1,0,1)$ and $z_t = 91.6$ pm.

The PES is originally represented in the twisted coordinate system and as such ex-
actly translational symmetry adapted. Cartesian coordinates, Fig. 1(c), are not symme-
try adapted. Fig. 2 shows two and one dimensional sections of the PES with energies
relevant for the present study. In particular the two dimensional section shows the
region that stationary states will predominantly occupy. These contour lines will be
repeated in subsequent figures in gray shaded tones to recall precisely this region.
The Cartesian coordinates of the fcc site indicated by 0 meV are $x_c = 137.6$ pm,
$y_c = -79.4$ pm, and $z_c = 90.0$ pm. The value of the energy in the hcp sites is about
15 meV above that. The analysis of the PES shows that all potential barriers that
separate the potential wells are equivalent with a height of 185 meV at $z_t = 90$ pm,
which reduces to 120 meV, when the z-coordinate is relaxed, see Fig. 2(b).

The DFT calculations underlying the PES were reported to yield a dihydrogen
associative desorption energy of 554 meV [15], which is about 61% of the experimental
value from isosteric heat measurements [27]. The latter value potentially needs to be
reviewed, to our knowledge there has not been any recent experimental reevaluation of
the adsorption energy of hydrogen or dihydrogen on palladium. Our evaluation of the
PES yields a desorption energy for a single hydrogen atom of about 3.5 eV. We consider
the energy domain beyond 1 to 2 eV above the minimum of the PES to be less well
represented by the expressions from Ref. 15, although a satisfactory agreement between
ab initio molecular dynamics using this PES and experimental sticking probabilities of
hydrogen atoms and molecules as a function of the coverage was reported [28] . The
quality of the PES in the low energy domain, such as depicted in Fig. 2, in particular
in relation to experimental and theoretical values of vibrational energies, was assessed
in Ref. 14, and will be further elucidated below in the discussion section. We should
mention that the analytical representation of the PES developed so far allows us to
straightforwardly include couplings between the motion of the substrate and adsorbate
atoms in the harmonic approximation.

## 2.2. The Multi-Configuration Time-dependent Hartree (MCTDH) method

The Multi-Configuration Time-dependent Hartree (MCTDH) [17–20] method, more
precisely the Heidelberg MCTDH package[29], is used. It was originally devised for the
numerical resolution of the time-dependent Schrödinger equation. The program allows
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The left panel shows a two dimensional representation of the potential energy surface in Cartesian coordinates for $z_c = 90$ pm with energy cut-off at 300 meV above the fcc minimum. The right panel shows a section of the PES along the $y_c$ coordinate for $x_c = d/2$ and $z = z_c = z_t = 90$ pm (solid lines) and relaxed $z$ (dashed lines). Upon relaxation, the $z$ coordinate stretches from about 90 to 106 pm at the $y_c = 0$ barrier.</figcaption></figure>

one to treat a wide range of problems in quantum molecular dynamics, in particular those with large dimension [23, 30]. The power of this method comes from the multi-configuration compact treatment of the wave function using a set of a time-dependent basis functions that are optimized by the Dirac-Frenkel variational principle[17–20]. All the calculations presented in this paper are thus based on the MCTDH representation of the nuclear wave function that is written as a sum of products of so-called single particle functions SPF, $\phi(q_k, t)$, where $q_k$ is one nuclear degree of freedom (DOF). The ansatz for the MCTDH wave function reads :

$$ \Psi(q_1, \dots, q_f, t) = \sum_{j_1}^{n_1} \dots \sum_{j_f}^{n_f} A_{j_1, \dots, j_f}(t) \prod_{\kappa=1}^{f} \phi_{j_\kappa}^{(\kappa)}(q_\kappa, t) \quad (2) $$

where $f$ denotes the number of degrees of freedom. There are $n_\kappa$ SPF for the $\kappa$th degree of freedom. For a given set of indices $j_1, \dots, j_f$ the multidimensional Hartree product $\prod_{\kappa=1}^{f} \phi_{j_\kappa}^{(\kappa)}(q_\kappa, t)$ in Eq. (2) defines a configuration. The number of config-
urations is given by the product $n_1...n_f$. The SPF, $\varphi_{j_\kappa}^{(\kappa)}(q_\kappa, t)$, are expressed in a time-independent basis set over one degree of freedom

$$ \varphi_{j_\kappa}^{(\kappa)}(q_\kappa, t) = \sum_{i_\kappa=1}^{N_\kappa} c_{i_\kappa}^{(\kappa, j_\kappa)}(t) \chi_{i_\kappa}^{(\kappa)}(q_\kappa), \quad (3) $$

where $\chi_{i_\kappa}^{(\kappa)}(q_\kappa)$ is a primitive basis functions, in general a Discrete variable representation (DVR) function [2], that depends on coordinate $q_\kappa$.

Large primitive bases are employed in order to describe the long-range interactions between the hydrogen atoms. Here, we used an exponential DVR for $x_t$ and $y_t$ with 61 DVR grid points and a sine DVR with 31 grid points for $z_t$. In comparison with most of the conventional MCTDH calculations, this primitive basis is quite large, albeit very large primitive basis sets have already employed in the past for similar problems [25]. We need a very large primitive basis to accurately describe the delocalization of one hydrogen atom in one cell that is treated periodically. In a separate work, we investigated the situation where two hydrogen atoms are included in the cell. The paper is in preparation and we shall show that the numerical treatment is possible. Generally speaking, MCTDH is not very sensitive to the size of the primitive basis set, the bottleneck is rather the size of the SPF basis set.

The exponential DVR is the appropriate choice for the DOF parallel to the substrate, as it satisfies the periodical boundary conditions of the dynamics. A sine DVR is reasonable for the perpendicular DOF, as wave functions vanish at the boundary of the grid, at which energies are much higher than those of interest in this study.

The Heidelberg MCTDH package can also be used to solve the time-independent Schrödinger equation[21, 22]. This is a convenient approach to obtain eigenstates with high accuracy. For this purpose, the package combines iteratively a propagation in imaginary time, a relaxation, of the SPF with an extraction of the low-lying eigenstates using the Davidson algorithm to determine the $A$-coefficients. The procedure converges to one or several eigenstate(s) that are expanded as linear combinations of products of SPF (configurations). Both the SPF and the coefficients in the expansion are optimized with a time-independent variational procedure. This iterative approach is called *improved relaxation* [21, 22] and this is the approach that is used in the present work. More precisely, we use the block improved relaxation method[22] which allows one to propagate simultaneously a set of initial vectors collectively called vector block to the system's eigenstates. In our simulations, we used up to 40, 40 and 8 SPF for $x_t$, $y_t$ and $z_t$, respectively. The natural populations of the weakest populated SPF generated in each final relaxation are less than $10^{-6}$. The elementary cell depicted in Fig. 1 is considered to be periodically extendable along $x_t$ and $y_t$. The primitive basis was also very well converged with an maximal error of 0.02% on the $z_t$ DVR for the highest
block of eigenstates. The main challenge in this paper, more specifically for the use of
the improved relaxation method, comes from the very high density of eigenstates in
the system under interest as well as their high degeneracy.

The use of MCTDH requires a structure of the Hamiltonian given by a sum of
products of low-dimensional operators to avoid large integrals. The kinetic operator
is in most cases in this form, but an *ab initio* potential energy surface is generally
not provided in this form. The Heidelberg MCTDH package includes a program called
potential fit (potfit) which, using a numerical procedure, rewrites the potential in the
required form. We will not describe the potfit algorithm here: we refer the reader to
Refs. 31–33 for a thorough description of this algorithm. We just briefly explain that
the potential is expressed in terms of (here one-dimensional) Single-Particle Potentials
(SPP) and that a “contraction” over one degree of freedom can be performed. This
contraction reduces the number of terms since only one function is used for this con-
tracted degree of freedom (DOF). We used 30 SPP for both $x_t$ and $y_t$, and $z_t$ was
the contracted mode. The root-mean square (rms) error of the refitting with potfit is
smaller than 0.04 meV in the physical region.

The fine grid mesh used on the twisted coordinates allows us to obtain numerically
exact the translational symmetry of the PES. The use of the twisted coordinates is
essential as they allow us to easily take into account that symmetry. The kinetic energy
operator contains coupling terms between these coordinates that are of the sum-of-
product form and thus straightforwardly implemented in MCTDH [14].

**3. Results and discussion**

**3.1. General considerations and conventions**

In this section, we present a thorough analysis of the first eigenstates of the hydrogen
adsorbed on the palladium Pd(111) surface. Eigenstates are superposition states of
zeroth order vibrational modes parallel or perpendicular to the substrate. Motivated
by the local $C_{3v}$ symmetry of the potential that holds at both the hcp and fcc sites, as
shown in Fig. 1, we adopt the following convention: $n^v$ labels states in a zeroth order
vibrational mode, where *n* is the number of the mode and *v* is the vibrational quantum
number of that mode. This notation was used in Ref. 34 and follows conventions for
the numbering of modes [35] and labelling of levels [36]. The latter builds on an early
suggestion by Callomon and Innes [37] who used the notation $n_{v''}$ to refer to a vibronic
transition from a vibrational state with $v''$ quanta of vibration in mode number *n* in
the lower electronic state to a vibrational state with *v*' quanta of vibration in the same
mode in the upper electronic state. This notation is reported in Ref. 38 to have first
been used by Watson in his Ph.D. thesis.
In the present case $n \in \{1, 2\}$, where 1 and 2 refer to vibrational modes perpendicular and parallel to the surface, respectively. Indeed, the perpendicular mode, associated with the adsorption or desorption motion, is locally of A₁ type, while the parallel mode, associated with a frustrated translation along the substrate, is two-fold degenerate of type E in a region close to the fcc (A) and hcp (B) equilibrium adsorption sites. These sites are energetically different and we therefore add subscripts A and B to mode 2 in order to refer to vibrational states pertaining to these specific adsorption sites on the Pd(111) substrate.

To zeroth order, a vibrational state is a configuration state that we might label $1^{v_1} 2^{v_2}$, where $S = A, B$ and $v_i = 0, 1, 2, \dots$ ($i = 1, 2$). For the sake of clarity, let us consider several examples. $1^0 2_B^1$ refers to a state located on an hcp potential well having a single vibrational quantum excitation parallel to the substrate, while it is in the ground vibrational state of mode 1. Of course, the eigenstates can be delocalized on both hcp and fcc sites. In this case, the eigenstates become superposition states of localized zeroth order states, which then become strongly mixed. To indicate a strong superposition or mixing of zeroth order states, we combine the labels of zeroth order states with the summation symbol "+". Note that this symbol is generic and refers to any possible relative phase between the mixed zeroth order states. For instance, $1^0 2_A^1 + 1^0 2_B^1$ represents a state with one quantum of vibration along mode 2, parallel to the substrate, which is delocalized on the hcp and fcc sites. Factorization is possible: $1^0 2_B^1 + 1^0 2_A^1 = 1^0 (2_A^1 + 2_B^1)$.

Table 1 gives the lowest energy vibrational eigenstates of H/Pd(111). The table lists levels (1st column) and sub-levels (2nd column), which are defined by specified sequences of states indicated by state number ranges (3rd column). As eigenstates are normally superpositions of many zeroth order states, the labels used in the table (4th column) indicate the zeroth order state having the strongest weight in the superposition, as determined by visual inspection of the wave functions. The eigenstate is then said to carry the signature of the corresponding zeroth order state. As explained above, when several zeroth order states contribute strongly to an eigenstate, they become mixed, which is indicated by the generic summation symbol "+". The listed energies (5th column) are degenerate at the given decimal precision. The degeneracies $g$ are given in the sixth column. Some energies are given in bold face, which shall be addressed below.
<table><caption>**Table 1.** Vibrational eigenenergies of H/Pd(111) in the (3×3)-grid. The last column of the table gives the near-degeneracy *g* associated with each level or sub-level. Bold faced energies refer to values obtained both on the 3 × 3 grid of Pd atoms as well as with 1 × 1-calculations (see text). The zero point energy at the lowest fcc site is around 160 meV. Notes: <sup>a)</sup> states contain some mixing of states of the A<sub>1</sub> type; <sup>b)</sup> states contain some mixing of states of the E type.</caption><thead><tr><th>Level</th><th>Sub-level</th><th>State number</th><th>Label</th><th>E<sub>n</sub>/meV</th><th>g</th></tr></thead><tbody><tr><td>1</td><td>-</td><td>1 ... 9</td><td>1<sup>0</sup> 2<sup>0</sup><sub>A</sub></td><td><strong>0.0</strong></td><td>9</td></tr><tr><td>2</td><td>-</td><td>10 ... 18</td><td>1<sup>0</sup> 2<sup>0</sup><sub>B</sub></td><td><strong>6.1</strong></td><td>9</td></tr><tr><td rowspan="3">3</td><td>3.a</td><td>19 ... 28</td><td>1<sup>0</sup> (2<sup>1</sup><sub>A</sub> + 2<sup>1</sup><sub>B</sub>)</td><td><strong>73.9</strong></td><td>10</td></tr><tr><td>3.b</td><td>29 ... 34</td><td>1<sup>0</sup> (2<sup>1</sup><sub>A</sub> + 2<sup>1</sup><sub>B</sub>)</td><td>74.5</td><td>6</td></tr><tr><td>3.c</td><td>35 ... 36</td><td>1<sup>0</sup> 2<sup>1</sup><sub>A</sub></td><td>74.9</td><td>2</td></tr><tr><td rowspan="3">4</td><td>4.c</td><td>37 ... 38</td><td>1<sup>0</sup> 2<sup>1</sup><sub>B</sub></td><td>78.0</td><td>2</td></tr><tr><td>4.b</td><td>39 ... 44</td><td>1<sup>0</sup> (2<sup>1</sup><sub>B</sub> + 2<sup>1</sup><sub>A</sub>)</td><td>78.6</td><td>6</td></tr><tr><td>4.a</td><td>45 ... 54</td><td>1<sup>0</sup> (2<sup>1</sup><sub>B</sub> + 2<sup>1</sup><sub>A</sub>)</td><td><strong>79.7</strong></td><td>10</td></tr><tr><td rowspan="3">5</td><td>5.a</td><td>55</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(A<sub>1</sub>) + 2<sup>2</sup><sub>B</sub>(A<sub>1</sub>)) + 1<sup>1</sup> (2<sup>0</sup><sub>A</sub> + 2<sup>0</sup><sub>B</sub>)</td><td><strong>99.5</strong></td><td>1</td></tr><tr><td>5.b</td><td>56 ... 61</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(A<sub>1</sub>) + 2<sup>2</sup><sub>B</sub>(A<sub>1</sub>)) + 1<sup>1</sup> (2<sup>0</sup><sub>A</sub> + 2<sup>0</sup><sub>B</sub>)</td><td>102.1</td><td>6</td></tr><tr><td>5.c</td><td>62 ... 63</td><td>1<sup>0</sup> 2<sup>2</sup><sub>A</sub>(A<sub>1</sub>) + 1<sup>1</sup> 2<sup>0</sup><sub>A</sub></td><td>105.4</td><td>2</td></tr><tr><td rowspan="3">6</td><td>6.c</td><td>64 ... 65</td><td>1<sup>0</sup> 2<sup>2</sup><sub>B</sub>(A<sub>1</sub>) + 1<sup>1</sup> 2<sup>0</sup><sub>B</sub></td><td>106.5</td><td>2</td></tr><tr><td>6.b</td><td>66 ... 71</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(A<sub>1</sub>) + 2<sup>2</sup><sub>B</sub>(A<sub>1</sub>)) + 1<sup>1</sup> (2<sup>0</sup><sub>A</sub> + 2<sup>0</sup><sub>B</sub>)</td><td>110.2</td><td>6</td></tr><tr><td>6.a</td><td>72</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(A<sub>1</sub>) + 2<sup>2</sup><sub>B</sub>(A<sub>1</sub>)) + 1<sup>1</sup> (2<sup>0</sup><sub>A</sub> + 2<sup>0</sup><sub>B</sub>)</td><td><strong>114.1</strong></td><td>1</td></tr><tr><td rowspan="3">7</td><td>7.a</td><td>73 ... 82</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(E) + 2<sup>2</sup><sub>B</sub>(E))</td><td><strong>125.5</strong></td><td>10</td></tr><tr><td>7.b</td><td>83 ... 88</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(E) + 2<sup>2</sup><sub>B</sub>(E))<sup>a)</sup></td><td>127.8</td><td>6</td></tr><tr><td>7.c</td><td>89 ... 90</td><td>1<sup>0</sup> 2<sup>2</sup><sub>A</sub>(E)<sup>a)</sup></td><td>129.3</td><td>2</td></tr><tr><td rowspan="3">8</td><td>8.c</td><td>91 ... 92</td><td>1<sup>0</sup> 2<sup>2</sup><sub>A</sub>(A<sub>1</sub>) + 1<sup>1</sup> 2<sup>0</sup><sub>A</sub><sup>b)</sup></td><td>130.2</td><td>2</td></tr><tr><td>8.b</td><td>93 ... 98</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(A<sub>1</sub>) + 2<sup>2</sup><sub>B</sub>(A<sub>1</sub>)) + 1<sup>1</sup> (2<sup>0</sup><sub>A</sub> + 2<sup>0</sup><sub>B</sub>)</td><td>130.9</td><td>6</td></tr><tr><td>8.a</td><td>99</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(A<sub>1</sub>) + 2<sup>2</sup><sub>B</sub>(A<sub>1</sub>)) + 1<sup>1</sup> (2<sup>0</sup><sub>A</sub> + 2<sup>0</sup><sub>B</sub>)</td><td><strong>131.4</strong></td><td>1</td></tr><tr><td rowspan="3">9</td><td>9.a</td><td>100</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(A<sub>1</sub>) + 2<sup>2</sup><sub>B</sub>(A<sub>1</sub>)) + 1<sup>1</sup> (2<sup>0</sup><sub>A</sub> + 2<sup>0</sup><sub>B</sub>)</td><td><strong>149.1</strong></td><td>1</td></tr><tr><td>9.b</td><td>101 ... 106</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(A<sub>1</sub>) + 2<sup>2</sup><sub>B</sub>(A<sub>1</sub>)) + 1<sup>1</sup> (2<sup>0</sup><sub>A</sub> + 2<sup>0</sup><sub>B</sub>)</td><td>150.6</td><td>6</td></tr><tr><td>9.c</td><td>107 ... 108</td><td>1<sup>0</sup> 2<sup>2</sup><sub>B</sub>(A<sub>1</sub>) + 1<sup>1</sup> 2<sup>0</sup><sub>B</sub><sup>b)</sup></td><td>152.4</td><td>2</td></tr><tr><td rowspan="3">10</td><td>10.c</td><td>109 ... 110</td><td>1<sup>0</sup> 2<sup>2</sup><sub>B</sub>(E)<sup>a)</sup></td><td>153.3</td><td>2</td></tr><tr><td>10.b</td><td>111 ... 116</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(E) + 2<sup>2</sup><sub>B</sub>(E))</td><td>155.0</td><td>6</td></tr><tr><td>10.a</td><td>117 ... 126</td><td>1<sup>0</sup> (2<sup>2</sup><sub>A</sub>(E) + 2<sup>2</sup><sub>B</sub>(E))</td><td><strong>156.3</strong></td><td>10</td></tr></tbody></table>
As aforementioned, the elementary cell shown in Fig. 1 displays 9 global equilibrium fcc sites and 9 local equilibrium hcp sites. This should lead to a degeneracy for each level of energy by multiples of 9. However, we will see that several tunneling effects split the energies in all the blocks of eigenvalues in multiple ways. The relaxation calculations have some numerical error stemming from the cumulative effect of the convergence of the relaxation process, the primitive basis size and the potfit program. It is estimated to be of the order of ±0.2 meV. For infinite large cells, sub-levels are members of the same vibrational sub-band and will likely show some dispersive behavior. In practice, for finite cells, they are therefore not perfectly degenerate, and the estimated level of accuracy is sufficient for the physical interpretations of the eigenstates. The labels proposed in Tab. 1 will be discussed in detail in the following. The separation of the states into pertaining to vibrational modes parallel and perpendicular to the substrate and the degeneracy linked to the number of sites included in the elementary cell can be seen from the labels that induce a "zero-order" description of our system. However, we shall see that the "zero-order" description is not satisfactory for a correct description of the eigenstates. In particular, we will see that the system is strongly impacted by several tunneling effects and a low-lying Fermi resonance.

### 3.2. Levels 1 and 2

In Tab. 1, the levels called 1 and 2 display a degeneracy factor of $g = 9$ at the numerical accuracy of the calculations. For the analysis of the eigenstates, two strategies have been adopted. First, we calculate the reduced probability densities (RPD) corresponding to the eigenstates. This will allow us to visualize properties such as the localization and the nodal structure of the eigenfunctions, which is crucial for their physical interpretation. Graphical representations of the RPD of all states documented in Tab. 1 are reproduced in the supplementary material. Some selected states will be analyzed in more detail in the following.

In addition, in order to investigate patterns that may arise from the local shape of the PES at the equilibrium sites we considered using a reduced elementary cell, which contains a single fcc and a single hcp site. We will refer to these calculations as the "1 × 1-calculations". More precisely, the definition of this 1 × 1 elementary cell is given by $x_t \in [-d/2, d/2]$ and $y_t \in [-d/2, d/2]$. Bold faced values are obtained likewise from the calculations on the 3 × 3 grid as well as from the 1 × 1-calculations. In the following discussion, they will prove to be useful reference energies.

In the case of the levels 1 and 2, the same values are obtained with the 1 × 1- and 3 × 3-calculations: 0.0 meV for the $1^0 2_A^0$ and 6.1 meV for the $1^0 2_B^0$ levels. In the 1 × 1-calculation, level degeneracy is $g = 1$, of course. Interestingly enough, the difference of 6.1 meV is smaller than the 15 meV between the two potential well minima see Fig. 2.
This is simply due to the difference in the local curvature of the stiffer fcc and more
widen hcp potential wells.

Figure 3 shows the RPD of two selected eigenstates in the two blocks of degenerate states, more precisely the states 7 on the left and the state 18 on the right hand side. The eigenstates are delocalized along the surface. By visual inspection, we observe that state 7 corresponds to hydrogen atoms localized on the fcc sites and state 18 to hydrogen atoms localized on the hcp sites. The node-less probability densities indicate that levels 1 and 2 correspond to ground states localized on the fcc and hcp sites, respectively.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Two-dimensional representation of the reduced probability densities (RPD) of two selected eigenstates in levels 1 (state 7, left hand side) and 2 (state 18, right hand side). Here and in the remaining figures showing RPD, the *x*, *y* and *z* coordinates refer to Cartesian coordinates *x*<sub>c</sub>, *y*<sub>c</sub> and *z*<sub>c</sub>, respectively, and 1 Å ≅ 100 pm. At the top, the probability densities pertaining to the eigenfunctions are integrated over *y*, at the bottom they are integrated over *z*. RPD are superimposed on equipotentials ranging from 0 to 300 meV represented by gray shaded contour lines (z = 90 pm= 0.9 Å); see also Fig. 2 for more details on the equipotential contour lines.</figcaption></figure>

3.3. Levels 3 and 4

Level 3 is a group of 18 states gathered in three degenerate sub-levels of 10, 6 and 2
states denoted as 3.a , 3.b and 3.c, respectively. Inspection of the RPD of all states
belonging to this level reveal nodal structures in the xy-plane. The RPD of state 28
shown in Fig. 4 is a typical example. The corresponding vibrational motion therefore
occurs parallel to the substrate. In a localized well, this motion is two-fold degenerate
which yields, together with the nine-fold configurational degeneracy of each stable well
in the 3 × 3 elementary cell, the 18 states composing this level. These 18 states split
in the three sub-levels.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Two-dimensional representation of the RPD of eigenstates 28 (left hand side) and 54 (right hand side), lying in sub-levels 3.a and 4.a, respectively. See also the caption of Fig. 3.</figcaption></figure>

The energy of the sub-level 3.a is also obtained from the 1 × 1-calculation. This is the reference energy of level 3. The energy of the sub-level 3.c is the highest in this cluster of sub-levels, about 1 meV above the energy of the sub-level 3.a. The energy difference between 3.c and 3.a is the width of the energy band pertaining to a continuum of Bloch waves that would be obtained when the grid was enlarged from the 1 × 1 elementary cell to the infinitely large cell. The band width can be readily obtained from a calculation using a 2 × 2 elementary cell, but also from the calculation on the basis of the 3 × 3 elementary cell. The splitting of level 3 into sub-levels results from tunneling, as will be further elucidated below. The 10+6+2 partitioning into three sub-levels can be rationalised from a group theoretical investigation based on the translational symmetry along the substrate which is not further developed here.

Level 4 includes also 18 states that split into three sub-levels of 2, 6 and 10 eigenstates called 4.c, 4.b and 4.a, respectively. The energy difference between 4.c and 4.a corresponds to the band width that would be obtained if the size of the cell tended to infinity. We note that the ordering is the exact opposite of that of level 3. The energy of sub-level 4.a is the reference energy obtained with the 1 × 1-calculation. By visual inspections of the RPD of the eigenstates one confirms that there is a correspondence between the subgroups of levels 3 and 4. Let us consider one example: Fig. 4 depicts the RPD of two typical eigenstates of the subgroups 3.a and 4.a. More precisely, we chose the eigenstates 28 and 54. We clearly see the node of a single mode at the equilibrium sites in the *xy*-plane parallel to the surface, while the RPD remain node-less along *z*. By visual inspection of Fig. 4 (left hand side), we see that state 26 has non-zero probability density at the saddle point linking the hcp and fcc sites. On the other hand,
state 54 (Fig. 4, right hand side) has an additional node at this point.

Based on the RPD obtained from the 1 × 1-calculation, we can deduce that the lower state of the group 3, i.e. state 26 in 3.a and the higher state in the group 4, i.e. state 54 in 4.a are two mutually corresponding components of the tunneling doublet through the barrier between the hcp and fcc sites. We shall refer this pair of states as constituting a "AB"-tunneling doublet in the following. Due to the fact that the energy of the fcc and hcp sites is not the same, the weights of the two sites constituting the doublet is not symmetrically distributed and we always observe a certain preponderance on one site depending on the state.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Two-dimensional representation of the RPD of eigenstates 32 (left hand side) and 43 (right hand side), lying in sub-levels 3.b and 4.b, respectively. See also the caption of Fig. 3.</figcaption></figure>

The splitting of levels 3 and 4 into three sub-levels of degenerate states is further analyzed as follows. The sub-levels 3.b, 4.b, 3.c and 4.c do not appear in the 1 × 1-calculation. Fig. 5 depicts the RPD of the eigenstates 32 and 43 that belong to the sub-levels 3.b and 4.b, respectively. The shape of the wavefunctions leads to the same conclusion as for sub-levels 3.a and 4.a: the eigenstates correspond to two mutually corresponding components of the AB-tunneling doublet. The main difference with respect to sub-levels 3.a and 4.a comes from the fact that the preponderance of the densities on the fcp and hcp sites is more pronounced. The same is true for the eigenstates 36 and 38 belonging to the sub-levels 3.c and 4.c except that the situation is more extreme and we postulate that, at the level of the accuracy of the figure, there is hardly any mixing of A and B sites: see Fig. 6.

Moreover, the extreme delocalization of the eigenstates and the fact that the eigenstates of the different sub-levels do not have the same energy leads to the following conclusion. The splitting of the levels 3 and 4 into sub-levels comes from a different
tunneling effect: from an A to an A site via a B site, and inversely. We will refer to this
process as an "ABA" or a "BAB" tunneling. The eigenstates of level 3 may be viewed
as being tunneling states of the ABA type, while the eigenstates of the level 4 are
tunneling states of the BAB type. In contrast, eigenstates of level 3, on one hand, and
level 4, on the other hand, form groups of mutually corresponding tunneling states of
the AB type. The detailed partitioning of the sub-level structure can be inferred from
a group theoretical account based on the translational symmetry along the substrate
and depends on the size of the elementary cell. For the 1 × 1 and 3 × 3 cells discussed
here, the partitioning is discrete with 1 and 3 sublevels, respectively. For infinitely large
cells, levels become continuous bands of Bloch states with band widths given by the
extreme sublevels: 3.c and 3.a, for level (band) 3, 4.a and 4.c, for level (band) 4.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Two-dimensional representation of the RPD of eigenstates 36 (left hand side) and 38 (right hand side), lying in sub-levels 3.c and 4.c, respectively. See also the caption of Fig. 3.</figcaption></figure>

To conclude with, eigenstates in levels 3 and 4 are characterized by the presence of
one quantum of vibrational motion parallel to the surface. Inspired by the local C<sub>3v</sub>
symmetry, modes with this property are labelled "2" and eigenstates result correspond-
ingly from configuration states 1<sup>0</sup> 2<sup>1</sup>. The eigenstates of level 3 are mainly located on
the fcc (A) site while the eigenstates of level 4 are mainly located on the hcp (B) site.
The mixing of probability density between the two sites is due to tunneling through
the barrier separating the A and B site, dubbed AB-tunneling. In addition, each level
spreads itself in energy due to another tunneling toward sites of the same type via the
site of the opposite type, called ABA or BAB-tunneling.
### 3.4. Levels 5, 6, 8 and 9

In this section we analyze the levels 5, 6, 8 and 9 that are all related together. Each level contains a group of 9 eigenstates with three sub-levels with different energies and degeneracies 1, 2 and 6. The 1 × 1-calculations provide the values of the reference energies in sub-levels L.a, with $L = 5, 6, 8$ and $9$. They correspond to the eigenstates 55, 72, 99 and 100 shown in Fig. 7: they are all delocalized.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Two-dimensional representation of the RPD of eigenstates 55 and 99 (left hand side), as well as 72 and 100 (right hand side). These are the states composing levels L.a, with $L = 5, 6, 8$ and $9$. See also the caption of Fig. 3.</figcaption></figure>

Compared to the states in levels 3 and 4 these eigenstates display an additional nodal structure along the z-axis. Thus, they correspond to states having one quantum in vibration perpendicular to the substrate. In local $C_{3v}$ symmetry, this motion has an irreducible representation of type $A_1$. Modes with this property are labeled "1". However, the situation is more complex as the RPD also spread in a highly symmetric
way in the *xy*-plane both at the fcc and at the hcp sites, which indicates that these states also involve vibrational excitations parallel to the substrate. Indeed, the second excited vibrational state of mode 2, in local *C*<sub>3v</sub> symmetry, reduces to one state of type *A*<sub>1</sub> and one doubly degenerate level of states of type E so that the 1<sup>0</sup> 2<sup>2</sup>(*A*<sub>1</sub>) configuration state may thus couple with configuration state 1<sup>1</sup> 2<sup>0</sup> and form a mixed state. This quantum phenomenon is known as Fermi resonance [39] and the corresponding doublet of eigenstates formed from the mixing of the 1<sup>1</sup> 2<sup>0</sup> and 1<sup>0</sup> 2<sup>2</sup>(*A*<sub>1</sub>) configuration states is called a Fermi pair.

Consequently, the eigenstates are characterized by (i) a global spreading in the *xy*-plane related to the tunneling between the different equilibrium sites and (ii) a local mixing of motions perpendicular and parallel to the substrate due to a Fermi resonance between states characterised by having two quanta in the parallel vibration and one quantum in the perpendicular vibration to the substrate. Configuration states are strongly coupled and the situation is too involved to pick up any residual feature of these states by a mere inspection of the RPD of the eigenstates. In order to disentangle the role of the different physical phenomena, we discuss in the appendix the interplay between Fermi resonance and tunneling coupling. These two quantum phenomena are well described in the literature [34, 40–43], but little is known of systems including both effects simultaneously.

In the appendix we define a local potential model based on the PES for H/Pd(111) and show that tunneling can reduce or enhance the effective Fermi resonance coupling strength, and vice-versa. In a supplementary material we propose a second discussion based on several model potentials to unravel further the interplay between Fermi resonance and tunneling coupling.

The main conclusion of the study carried out in the appendix, and supported by the studies in the supplementary material, is that there is the formation of a Fermi pair between two vibrational configuration states in the local potential. The ensuing eigenstates have energies of 102.7 meV and 129.6 meV (see Tab. A-2). Consequently, we can conclude firmly that state 55 is the lower component of the Fermi resonance doublet and that state 99 is its upper component. In the same manner, states 72 and 100 are the lower and upper energy components of a second Fermi resonance pair, respectively. From the energies, states 55 and 99 are more likely related to the fcc (A) site, while states 72 and 100 are more likely related to the hcp (B) site. However, RPD are delocalized over both sites. By visual inspection of Fig. 7, we notice that states 55 and 99 have non-vanishing probability densities at the saddle points between the fcc and hcp sites, whereas states 72 and 100 have a nodal structure at these sites. Consequently, we can conclude that state 99 is simultaneously the upper energy Fermi resonance partner of state 55 and the lower energy tunneling partner of state 100, while
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Two-dimensional representation of the RPD of eigenstates 59 and 97 (left hand side), as well as 69 and 104 (right hand side), as typical states composing levels L.b, with L as in Fig. 7. See also the caption of Fig. 3.</figcaption></figure>

state 72 is simultaneously the lower energy Fermi resonance partner of state 100 and the upper energy tunneling partner of state 55.

Armed with the conclusions of the appendix and the supplementary material, we can also say that the eigenstates in levels 5 and 8 may be considered to be nearly equal weight mixtures of $1^1 2^0$ and $1^0 2^2$ zeroth order configuration states, while the states in level 6 seem to have a slight preponderance of $1^1 2^0$ configuration states, and in level $9 1^0 2^2$ configuration states preponderate. In addition, since level 6 also has some preponderance of zeroth order states localized on the fcc site, and level 9 of states localized on the hcp site, one might say that the Fermi resonance coupling between zeroth order states located in the hcp site is slightly weaker than at the fcc site. The reader is referred to the supplementary material for the visual inspection of the complete set of graphical representations of RPD.
The situation is even more involved with the sub-levels 5.α, 6.α, 8.α and 9.α, where α = b and c since these sub-levels are degenerate. Fig. 8 shows RPD of typical states in the sub-levels 5.b, 6.b, 8.b and 9.b. Due to the degeneracy the exact attribution is very tricky by visual inspection of the RPD. However, by careful comparison of RPD within a sub-level, we can nevertheless conclude, with the help of the analysis carried out in the appendix, that levels 5.b with 8.b on one hand, and 6.b with 9.b on the other hand contain respectively lower and upper energy components of a Fermi pair, while 5.b with 6.b on one hand, and 8.b with 9.b on the other hand contain respectively lower and upper energy components of tunneling doublets. Similar conclusions can be drawn for α = c. The latter case is even more complex since with increasing distance from the most stable adsorption sites, which exhibit a local $C_{3v}$ symmetry, i.e. with increasing energy the local symmetry ceases to hold, and hence some mixing of E and $A_1$ type states become apparent in the densities. The remaining global symmetry is then given by the translational group. The supplementary material contains the graphical representations of the RPD for α = c.

As for levels 3 and 4, the lifting of configurational degeneracy in all four levels is
due to ABA or BAB tunneling. The details in the remaining quasi degeneracies may
be deduced from further considerations based on the translational symmetry.

To conclude the analysis of this section, level 5 is an ABA, level 8 is a BAB tun-
neling level; they respectively correspond to lower and upper energy components of
Fermi resonance pairs of zeroth order vibrational states containing one quantum of
perpendicular and two quanta of parallel vibration; levels 6 and 9 are related one to
the other in the same manner. Levels 5 and 6 on one hand, and levels 8 and 9 on the
other hand are pairs of AB-tunneling levels.

3.5. Levels 7 and 10

Let us now discuss the last two levels: levels 7 and 10. These states are essentially composed of the E-type of zeroth order states with two vibrational quanta in parallel direction, i.e. of configuration states of the kind 1º 2º(E). They are also impacted by tunneling. The RPS of typical states of the reference sub-levels 7.a and 10.a are depicted in Fig. 9. They contain some components of zero order vibrational states of the $A_1$ type. These mixtures become stronger with increasing energy due to the gradual breakdown of the $C_{3v}$ symmetry as aforementioned. States in the other sub-levels are less delocalized.

We can however clearly see that states in level 7 correspond to the lower energy
components of AB tunneling pairs of zeroth order vibrational states containing two
quanta of parallel vibration, while level 10 correspond to the upper energy components.
The relation between these levels is thus similar to that between levels 3 and 4, but
both their average energy gap as well as their internal ABA or BAB tunneling band
spread are larger than in the case of the lower energy levels. Interestingly enough,
the average energy of level 7 is of the same order than the barrier of ~ 120 meV to
AB-tunneling on the electronic potential, see Fig. 2. The average energy of level 10 is
even beyond. However, this does not mean that the hydrogen atoms can move more
freely over the barrier, since there will still be a zero point energy from the motion
orthogonal to the tunneling path that is added to the electronic barrier, by which the
effective barrier is increased. Indeed, it is proven that, in multidimensional systems,
tunneling can occur "over the barrier" [34, 42], and can even be inhibited far above
the barrier [41].

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Two-dimensional representation of the RPD of eigenstates 74 (left hand side) and 118 (right hand side), lying in sub-levels 7.a and 10.a, respectively. See also the caption of Fig. 3.</figcaption></figure>

3.6. *Comparison with other work*

Based on the PES derived in Ref. 13, Saalfrank and Tremblay in Ref. 12 obtained from
calculations on a 1 × 1-grid roughly 89 and 114 meV for the corresponding levels 3 and
5 in table 1. The values are larger than the values from the present work, but they
do see the Fermi resonance between zeroth order modes of the type 1¹ and 2² as in
the present work, although this finding remained unpublished. The present vibrational
energies are also smaller than those published earlier on the same PES but using a 2 × 2
elementary cell [14]. The reason for this mismatch is a scaling factor mistakenly used
in the evaluation of the energy levels obtained from the MCTDH program in Ref. 14.

Moreover, experimental data are available from measurements using the High Res-
olution Electron Energy Loss Spectroscopy (HREELS) method. These measurements
were carried out by Conrad *et al.* [8] and do not dispose of sufficient resolution to
extract tunneling splittings. Nevertheless, these measurements enabled the extraction of two peaks at 96 meV and 124 meV, which were associated to states with parallel and perpendicular excitation, respectively, in Ref. 8. However, additional aspects hamper the direct comparison between these experimental and the theoretical data investigated in this work.

First, a detailed theory for the intensity distribution of the HREELS spectrum as a function of the loss energy is highly complicated [44, 45]. Selection rules are not clear or difficult to be assessed. Secondly, the experimental approach followed to extract the two transition energies does not allow us to know whether these energies correspond to that of an isolated hydrogen atom located on an fcc, or on an hcp site, or whether the latter is isolated or surrounded by other hydrogen atoms adsorbed on the surface. The set of vibrational levels of dissociatively adsorbed H₂ molecules with mutually interacting H atoms differs quite significantly from that of isolated adsorbed H atoms [46].

The calculation of the dynamical structure factor (DSF) of diffracted low energy electrons might enable us to directly link experimental observations in the HREELS experiments to a simulated HREELS spectrum. Work is in progress to simulate this spectrum and preliminary results seem to confirm the quality of the underlying PES.

## 4. Conclusion

In this work, we calculate the vibrational states of H/Pd(111) up to doubly excited vibrations parallel to the palladium substrate represented by a periodic 3 × 3-grid. We thus obtain the 126 lowest lying vibrational eigenstates of this system. The calculation presents a real challenge for numerical computation. The presence of long-range interactions between the hydrogen atoms mean that large periodic calculations cells are required to describe the system correctly, necessitating large primitive bases. The presence of small splittings necessitates a high control of the accuracy of the calculations. To this end, the improved relaxation approach based on the MCTDH method is well adapted as it allows us to reach the necessary level of accuracy. The analysis of the vibrational eigenstates shows that the physics is complex and rich, with as highlight the manifestation of a Fermi resonance between local vibrational modes of single hydrogen atoms where a state with a single quantum of vibrational excitation perpendicular to the substrate strongly couples with a state carrying two quanta of vibrational excitation parallel to the substrate. In addition, the presence of tunneling affects all excited vibrational states. This tunneling effect manifests itself in a variety of ways across different potential barriers, leading to multiple splittings of several energy levels into sub-levels with relative gaps ranging from 1 meV to 8 meV. This work can be seen as a first step to rationalize the interaction between hydrogen and a surface
of palladium. Calculations of the vibrational stationary states of H₂/Pd(111) and of the quantum dynamics of hydrogen and dihydrogen adsorbed on Pd(111) to quantify the diffusion of these species will be presented in future papers, where particular attention will be given to discern states pertaining to the ortho and para modifications of dissociatively adsorbed dihydrogen.

## Acknowledgement

We acknowledge financial support from ANR through the QDDA project (ANR-19-CE29-0007) as well as from Université de Strasbourg, Université de Lorraine, Université de Paris-Saclay, and CNRS for computational resources. The authors also thank the CNRS International Research Network (IRN) "MCTDH" for funding support.

## References

[1] Z. Bačić, R. Whitnell, D. Brown and J. Light, *Comp. Phys. Comm.* **51** (1), 35–47 (1988).

[2] Z. Bačić and J.C. Light, *Annu. Rev. Phys. Chem.* **40**, 469 (1986).

[3] M. Von Dirke, Z. Bačić, D.H. Zhang and J.Z.H. Zhang, *J. Chem. Phys.* **102** (11), 4382–4389 (1995).

[4] Y. Qiu and Z. Bačić, in *Advances in Molecular Vibrations and Collision Dynamics, Vol. III Molecular Clusters*, edited by J. Bowman and Z. Bačić (JAI Press, London, 1998), pp. 183–204.

[5] F. Gatti, editor, *Molecular Quantum Dynamics* (Springer, Heidelberg, 2014).

[6] R. Marquardt and M. Quack, editors, *Molecular Spectroscopy and Quantum Dynamics* (Elsevier, St. Louis, Missouri, 2020).

[7] A.D. Jewell, G. Peng, M.F. Mattera, E.A. Lewis, C.J. Murphy, G. Kyriakou, M. Mavrikakis and E.C.H. Sykes, *ACS nano* **6** (11), 10115–10121 (2012).

[8] H. Conrad, M.E. Kordesch, R. Scala and W. Stenzel, *J. Electr. Spectr. Rel. Phen.* **38**, 289–298 (1986).

[9] A. Jardine, H. Hedgeland, G. Alexandrowicz, W. Allison and J. Ellis, *Progr. Surf. Sci.* **84** (11-12), 323–379 (2009).

[10] A. Bahel and Z. Bačić, *J. Chem. Phys.* **111** (24), 11164–11176 (1999).

[11] J. Reilly, *Hydrogen: Its Technology and Implications* **1**, 13–48 (2018).

[12] J.C. Tremblay and P. Saalfrank, *J. Chem. Phys.* **131** (8), 084716 (2009).

[13] N. Ozawa, T.A. Roman, H. Nakanishi, H. Kasai, J. Nelson B. Arboleda and W.A. Dino, *J. Appl. Phys.* **101** (12), 123530 (2007).

[14] T. Firmino, R. Marquardt, F. Gatti and W. Dong, *J. Phys. Chem. Lett.* **5** (24), 4270–4274 (2014).

[15] Y. Xiao, W. Dong and H. Busnengo, *J. Chem. Phys.* **132** (1) (2010).

[16] W. Dong, G. Kresse and J. Hafner, *J. Mol. Cat. A: Chemical* **119** (1-3), 69–76 (1997).
[17] H.D. Meyer, U. Manthe and L.S. Cederbaum, Chem. Phys. Lett. **165**, 73–78 (1990).

[18] U. Manthe, H. Köppel and L.S. Cederbaum, J. Chem. Phys. **95**, 1708 (1991).

[19] M.H. Beck, A. Jäckle, G.A. Worth and H.D. Meyer, Phys. Rep **324**, 1–105 (2000).

[20] H.D. Meyer, F. Gatti and G.A. Worth, *Multidimensional quantum dynamics: MCTDH theory and applications* (Wiley, New York, 2009).

[21] H.D. Meyer, F. Le Quéré, C. Léonard and F. Gatti, Chem. Phys. **329**, 179–192 (2006).

[22] L.J. Doriol, F. Gatti, C. Iung and H.D. Meyer, J. Chem. Phys. **129**, 224109 (2008).

[23] L. Shi, M. Schröder, H.D. Meyer, D. Peláez, A.M. Wodtke, K. Golibrzuch, A.K. A.-M. Schönemann and F. Gatti, J. Chem. Phys. **159**, 194102 (2023).

[24] M. Nest and R. Kosloff, J. Chem. Phys. **127** (13), 134711 (2007).

[25] Q. Meng and H.D. Meyer, J. Chem. Phys. **146** (18), 184305 (2017).

[26] Q. Meng and H.D. Meyer, J. Chem. Phys. **139** (16), 164709 (2013).

[27] H. Conrad, G. Ertl and E. Latta, Surf. Sci. **41** (2), 435 – 446 (1974).

[28] Y. Xiao and W. Dong, Phys. Rev. B **83** (12), 125418 (2011).

[29] G.A. Worth, M.H. Beck, A. Jäckle and H.D. Meyer, The MCTDH Package, Version 8.2, (2000). H.-D. Meyer, Version 8.3 (2002), Version 8.4 (2007). See http://mctdh.uni-hd.de.

[30] M. Schröder, F. Gatti, D. Lauvergnat, H.D. Meyer and O. Vendrell, Nat. Com. **13**, 618 (2023).

[31] A. Jäckle and H.D. Meyer, J. Chem. Phys. **109**, 2614 (1998).

[32] A. Jäckle and H.D. Meyer, J. Chem. Phys. **109**, 3772 (1998).

[33] F. Gatti and H.D. Meyer, Chem. Phys. **304**, 3–15 (2004).

[34] R. Marquardt, M. Sanrey, F. Gatti and F. Le Quéré, J. Chem. Phys. **133**, 174302 (2010).

[35] R.S. Mulliken, J. Chem. Phys. **23** (11), 1997–2011 (1955), Joint Commission for Spectroscopy from IUPAC and IUPAP.

[36] E.R. Cohen, T. Cvitas, J.G. Frey, B. Holmström, K. Kuchitsu, R. Marquardt, I. Mills, F. Pavese, M. Quack, J. Stohner, H.L. Strauss, M. Takami and A.J. Thor, *Quantities, Units and Symbols in Physical Chemistry ('The Green Book')*, 3rd ed. (RSC, London, 2007).

[37] J. Callomon and K. Innes, Journal of Molecular Spectroscopy **10** (1–6), 166 – 181 (1963).

[38] J. Brand, J. Callomon, K. Innes, J. Jortner, S. Leach, D. Levy, A. Merer, I. Mills, C. Moore, C. Parmenter, D. Ramsay, K. Narahari Rao, E. Schlag, J. Watson and R. Zare, Journal of Molecular Spectroscopy **99** (2), 482–483 (1983).

[39] E. Fermi, Z. Phys. **71** (3-4), 250–259 (1931).

[40] M. Razavy, *Quantum theory of tunneling* (World Scientific, Singapur, 2013).

[41] B. Fehrensen, D. Luckhaus and M. Quack, Chem. Phys. Lett. **300**, 312–320 (1999).

[42] M. Quack and G. Seyfang, in *Molecular Spectroscopy and Quantum Dynamics*, edited by R. Marquardt and M. Quack, Chap. 7 (Elsevier, St. Louis, Missouri, 2020).

[43] R. Marquardt and M. Quack, in *Encyclopedia of Chemical Physics and Physical Chemistry*, edited by Moore, J. and Spencer, N., Vol. I, Chap. A3.13 (IOP, Bristol, 2001), pp. 897–936.

[44] H. Conrad, M. Kordesch, W. Stenzel, M. Sunjić and B. Trninic-Radja, Surf. Sci. **178** (1), 578–588 (1986).
[45] A.M. Bradshaw and N.V. Richardson, *Pure Appl. Chem.* **68**, 457–467 (1996).

[46] O. Bindech, *Quantum dynamics of the diffusion of adsorbates*, Ph.D. thesis, Université de Strasbourg, 2023.

[47] P.F. Bernath, *Spectra of atoms and molecules* (Oxford University Press, Oxford, 2005).

[48] R. Marquardt, *J. Math. Chem.* **50** (3), 577–587 (2012).
# Appendix: Local potential study of the H/Pd(111) system

The time (in)-dependent calculations we have carried out on the H/Pd(111) system are based on a periodic approach with an elementary cell containing 18 potential wells, 9 of which are of the fcc type and 9 of the hcp type. The latter are separated by a potential barrier of finite size, allowing tunneling through it to take place. The presence of tunneling has an impact not only on the surface dynamics of the adsorbates, but also on the structure of their vibrational stationary states. It is interesting to study the local properties of the vibrational states of the hydrogen atom in the vicinity of a single potential well in the absence of tunneling. To this end, an analytical three-dimensional potential was constructed to reproduce locally the shape of the global potential from ref. [15] at an fcc site. We will refer to these two potentials as the local PES and global PES in the following.

The local PES is designed to fit the global PES at low energies around the fcc potential well. As energy increases, the global PES tends to fall off in certain directions to form a potential barrier, while the local PES continues to increase. The common vibrational states of the two systems will then begin to show differences as they approach the potential barrier. All the local properties of the global PES have been integrated into this potential to preserve the local $C_{3v}$ symmetry of the vibrational states, as well as the different correlations between the DOF of the system.

In the following, we first give the technical details of the construction of this local PES and its implementation in MCTDH, then discuss the results from calculations of the obtained local vibrational hydrogen states at the fcc site.

## A-1. Details about the construction of the local PES

The analytical expression of the function used to reproduce the PES contains elements chosen to match the shape of the potential and its symmetry. In the following, we will work in the Cartesian coordinate system $(x, y, z)$ to locate the position of the hydrogen atom. The potential model used is of the form :

$$V(x, y, z) = V_{\perp}(z) + V_{\parallel}(x, y)f_{\text{damp}}(z) \quad (A-1)$$

Where

$$V_{\perp}(z) = \frac{1}{2} f_{\perp} y_{\perp}^2(z) \quad (A-2)$$
$$ 
\begin{aligned}
V_{\parallel}(x, y) = & \frac{1}{2} f_{\parallel} R^2(x, y) \\
& \times \exp[-a_{\parallel 2} R^2(x, y) - a_{\parallel 3} R^3(x, y) \sin(3\varphi(x, y))] \\
& + a_{\parallel 6a} R^6(x, y) \sin^2(3\varphi(x, y)) + a_{\parallel 6b} R^6(x, y) \sin^2(3\varphi(x, y))
\end{aligned} 
\quad (A-3) 
$$

With

$$ 
\left. 
\begin{aligned}
y_{\perp}(z) &= \frac{e^{-a_{\perp}(z - z_{\text{eq}})}}{a_{\perp}} \\
R^2(x, y) &= (x - x_{\text{eq}})^2 + (y - y_{\text{eq}})^2 \\
\varphi(x, y) &= \arctan\left(\frac{y - y_{\text{eq}}}{x - x_{\text{eq}}}\right)
\end{aligned}
\right\} 
\quad (A-4) 
$$

$x_{eq}$, $y_{eq}$ and $z_{eq}$ are the equilibrium coordinates of a fcc site. The function $f_{\text{damp}}$ is given by

$$ f_{\text{damp}}(z) = e^{-a_{\text{damp}}(e^{b_{\text{damp}}(z - z_{\text{eq}})} - 1))} \quad (A-5) $$

and has a double role: first, in the case that the parameters $a_{\text{damp}}$ and $b_{\text{damp}}$ are positive, it leads to a damping of the parallel potential $V_{\parallel}$ for large distances from the substrate, i.e. the parallel potential becomes a flat potential for large values of $z$. Secondly, because of the low order expansion :

$$ V(x, y, z) \approx \frac{1}{2} f_{\perp} \Delta z^2 + \frac{1}{2} f_{\parallel} (\Delta x^2 + \Delta y^2) - \frac{1}{2} f_{\parallel} a_{\text{damp}} b_{\text{damp}} (\Delta x^2 + \Delta y^2) \Delta z \quad (A-6) $$

where $\Delta\kappa = \kappa - \kappa_{\text{eq}}$ for $\kappa = \{x, y, z\}$. The term $\propto f_{\parallel}a_{\text{damp}}b_{\text{damp}}$ is potentially generating a Fermi resonance coupling between the first overtone of the parallel and the fundamental of the perpendicular mode.

By judiciously setting $f_{\parallel} = e^{l_{f_{\parallel}}}$, $f_{\perp} = e^{l_{f_{\perp}}}$, $a_{\text{damp}} = e^{l_{a_{\text{damp}}}}$ and $b_{\text{damp}} = e^{l_{b_{\text{damp}}}}$, and letting the logarithmic parameters $l_{f_{\parallel}}$, $l_{f_{\perp}}$, $l_{a_{\text{damp}}}$ and $l_{b_{\text{damp}}}$ be actual fitting parameters, one naturally ensures positive definiteness of the potential and the said damping behavior.

The terms

$$ 
\begin{aligned}
R^3 \sin(3\varphi) &= y^3 - 3xy^2 \\
R^3 \cos(3\varphi) &= x^3 - 3yx^2
\end{aligned} 
\quad (A-7) 
$$

are respectively $A_1$ and $A_2$ labeled cubic terms of $C_{3v}$-symmetry [47]. The factor $\exp(-a_{\parallel 2}R^2...)$ in equation A-1 is hence a lowest order $C_{3v}$-symmetric modulation
<table><caption>**Table A-1.** Optimal parameter values for the model potential equation A-1.</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>f<sub>&perp;</sub> &times; &Aring;<sup>2</sup>/eV</td><td>2.9027</td></tr><tr><td>f<sub>&parallel;</sub> &times; &Aring;<sup>2</sup>/eV</td><td>2.9496</td></tr><tr><td>a<sub>&parallel;2</sub> &times; &Aring;<sup>2</sup></td><td>2.22</td></tr><tr><td>a<sub>&parallel;3</sub> &times; &Aring;<sup>3</sup></td><td>9.5</td></tr><tr><td>a<sub>&parallel;6a</sub> &times; &Aring;<sup>6</sup></td><td>-4.2</td></tr><tr><td>a<sub>&parallel;6b</sub> &times; &Aring;<sup>6</sup></td><td>4.0</td></tr><tr><td>a<sub>damp</sub></td><td>3.56</td></tr><tr><td>b<sub>damp</sub> &times; &Aring;</td><td>1.34</td></tr></tbody></table>

of the otherwise cylindrical potential. By setting $a_{\parallel 2} = e^{l_{\parallel 2}}$, one additional enforces a damping of the modulation with larger displacements from the stable equilibrium site, avoiding in this way side minima. The various parameters given in this model are obtained by adjustment to the global PES around one of the fcc wells on the surface. We have taken the region of space defined by:

$$ \begin{array}{l} 0.9756 \AA \le x \le 1.8756 \AA \\ -1.3756 \AA \le y \le 0.0 \AA \\ 0.6 \AA \le z \le 1.5 \AA \end{array} \qquad (A-8) $$

The model is fitted to the 600 data points with the simple version of the Marquardt-Levenberg algorithm [48] by which the optimum set of parameters is obtained and given in Tab. A-1.

The root mean square deviation of this fit is 34.3 meV over all points. The total energy reached in the spatial volume used for the fit is 6.2 eV. The energy interval in which the energy difference between the two potentials must be minimal is below that of the potential barrier, i.e. around 186 meV (majoring). The vast majority of points showing significant deviation are all above this value. A two-dimensional representation of the local PES (solid line) superimposed on the global PES (dashed line) is given in figure A-1.

## A-2. Results and discussion

The vibrational stationary states of this system are calculated with MCTDH using the block improved relaxation program. The analytical model of the local PES does not allow decomposition as a sum of product a one-dimensional operator. We therefore used POTFIT to represent the potential in a form suitable for MCTDH calculations. The same primitive basis used for the global potential is used for these calculations.
<figure><img src="image_12.png" /><figcaption>**Figure A-1.** Two dimensional section of the PES along the substrate coordinates (z = 90Å). The local PES (solid line) is superimposed on the global PES (dashed line).</figcaption></figure>

The rms error on the fit is totally negligible for this system.

We are only interested in the vibrational states corresponding to those calculated in the global potential in the core of the paper. That is ground state, the levels of excited states with a single excitation, then part of the levels of doubly-excited states. This corresponds to a total of 7 vibrational states. Table A-2 and figure A-2 show the vibrational eigenenergies and two-dimensional representations of the reduced probability densities of the local vibrational states.

<table><caption>**Table A-2.** Vibrational eigenenergies of H/Pd(111) in the local PES.</caption><thead><tr><td>Level</td><td>State number</td><td>State label</td><td>E<sub>n</sub>/meV</td><td>Degeneracy g</td></tr></thead><tbody><tr><td>1</td><td>1</td><td>1<sup>0</sup>2<sup>0</sup></td><td>0.0</td><td>1</td></tr><tr><td>2</td><td>2,3</td><td>1<sup>0</sup>2<sup>1</sup></td><td>80.0</td><td>2</td></tr><tr><td>3</td><td>4</td><td>1<sup>0</sup>2<sup>2</sup>(A<sub>1</sub>) + 1<sup>1</sup>2<sup>0</sup></td><td>102.7</td><td>1</td></tr><tr><td>4</td><td>5</td><td>1<sup>0</sup>2<sup>2</sup>(A<sub>1</sub>) + 1<sup>1</sup>2<sup>0</sup></td><td>129.6</td><td>1</td></tr><tr><td>5</td><td>6,7</td><td>1<sup>0</sup>2<sup>2</sup>(E)</td><td>132.4</td><td>2</td></tr></tbody></table>

To label the vibrational states of the system, we have adopted the nomenclature introduced in the main text for configuration states. The symbol $n^v$ designates a state having $v$ quanta of energy in mode $n$, where $n \in \{1, 2\}$. Vibrations perpendicular and parallel to the substrate are denoted by 1 and 2 respectively. The ground state is denoted by $1^02^0$. For example, $2^1$ refers to a state with a single vibrational excitation parallel to the surface. In this simplified nomenclature, no distinction is made between vibrational modes parallel to the substrate separately. These two modes are degenerate because of the local $C_{3v}$ symmetry of the potential, so we speak only of excitation parallel or perpendicular to the substrate.

The structure of the reduced probability density of state 1 is node-less, indicating
a non-excited state. This state is considered to be the ground state of the system. Its energy is around 155 meV (with respect to the potential minimum), and is taken as a reference in table A-2.

The first level of excited states contains two vibrational states excited in the parallel modes. Figure A-2 shows that the wave function of these two states has a node localized at the center of the potential well, with two anti-node structures (states 2 and 3). Given the $C_{3v}$ symmetry of the potential, partial or total degeneracy results in excited modes parallel to the substrate. Here, states 2 and 3 are perfectly degenerate, with an energy of 80 meV.

This degeneracy appears partially in the doubly excited state level $1^02^2$, which contains two degenerate states and a third one with a slightly lower energy. These three states differ in their $C_{3v}$ symmetry labels. The two degenerate states of mode 2 have E label and are denoted by $2^2(E)$ (states 6 and 7), while the third one has an $A_1$ label denoted by $2^2(A_1)$.

The latter appears, together with state $1^1$ of mode 1, in states 4 and 5. The reduced probability density of these states in the $xz$-plane also shows the presence of a node at this plane. These states therefore simultaneously have one quantum of vibrational excitation in mode 1 and two quanta of vibrational excitation in mode 2. In a harmonic approximation, this would correspond to the coupling between the fundamental of the perpendicular mode and the first overtone of parallel mode. This is akin to a Fermi resonance between these two modes. The latter belong to the same symmetry representation group ($A_1$) and may thus couple. Given that the potential is more contracted in the perpendicular mode, a state with one quantum of excitation in mode 1 could have an energy close to that with two quanta of excitation in mode 2. They are therefore susceptible to couple resonantly. The resulting vibrational states, known as Fermi pairs, are given by the combination of these two modes. We use this image here for the nomenclature of these states. States 4 and 5 are therefore both labeled by $1^0 2^2(A_1)+1^1 2^0$.
<figure><img src="image_13.png" /><figcaption>**Figure A-2.** Two-dimensional representation of the reduced probability densities of the vibrational eigenstates of the local PES. The latter are superimposed on equipotentials represented by dashed red lines.</figcaption></figure>

If we now consider a second neighboring potential well located at the hcp site with vibrational states that possess the same local properties as those of the fcc well, then the presence of the finite-size potential barrier between the two wells can lead certain
vibrational states at the fcc site to couple with those located at the hcp site and form tunnel doublets. If states 4 and 5 are also affected by tunneling, then each of them will form a tunnel doublet with a corresponding state in the hcp well. We thus have 4 vibrational states forming two Fermi pairs and two tunnel doublets simultaneously. Figure A-3 shows schematically the impact of successively introducing the tunneling effect between the two potential wells, then Fermi resonance between the $1^02^2(A_1)$ and $1^12^0$ states along the modes localized at each well. In the following, we refer to the fcc and hcp potential wells as A and B respectively. Well labels are attributed to mode 2, i.e. $2_A$ and $2_B$, so that state labels become $1^02_A^2(A_1)$ and $1^12_A^0$, or $1^02_B^2(A_1)$ and $1^12_B^0$. For the sake of simplicity, we shall use the short-cut notations $1^12_A^0 \equiv 1_A^1$, $1^12_B^0 \equiv 1_B^1$, $1^02_A^2 \equiv 2_A^2$, and $1^02_B^2 \equiv 2_B^2$, where the 2 mode symmetry label will be omitted, see Figure A-3.

<figure><img src="image_14.png" /><figcaption>**Figure A-3.** Schematic representation showing the tunnel doublets and Fermi pairs that can be formed between a quadruplet of states.</figcaption></figure>

Figure A-3a represents a system where the vibrational states at each well are assumed to be fully decoupled. Introducing tunneling between the parallel modes located at A and B shifts their energies and leads to the formation of a first tunnel doublet (Figure A-3b). If we also introduce a local coupling that induces Fermi resonance between modes 1 and $2^2$ (Figure A-3c), this will lead to the formation of two Fermi pairs localized at each well. The $1_A^1$ and $1_B^1$ modes are now also impacted by tunneling indirectly through coupling with the $2^2$ modes, forming a second tunnel doublet. The four vibrational states thus obtained contain two Fermi pairs (blue) and two tunnel doublets (red), as illustrated in Figure A-3d. If the system contains 9 equivalent potential wells of each type, as in the example of the H/Pd(111) system studied in the main text, then each of these states will be 'in principle' 9 fold times degenerate. Coupling will take place between 4 levels, each containing 9 vibrational states. We refer to these levels as tunnel doublet and Fermi pair levels.

Tunneling is generally followed by a change in the wavefunction of the involved
state, as well as its energy. If this state is part of a Fermi pair, tunneling can weaken or strengthen the intensity of the local resonance. Inversely, the Fermi resonance can modify considerably the local structure of the wavefunction, leading thus to a weakening or enhancement of the tunneling effect. Depending on the relative strength of the tunneling and Fermi resonance coupling, the resulting vibrational state will be dominated by the 'parent' mode that constitutes it, and by looking at the structure of the reduced probability densities of these states, we can trace them back to the 'parent' mode of the vibrational state under consideration. When both the tunneling and Fermi resonance coupling are strong, i.e. at close to optimal resonance conditions, it becomes difficult to extract the dominant character of one coupling scheme, given that the nodes are perfectly shared between the two coupled modes.

This desire to distinguish between these coupled states is of no particular interest when studying the spectrum of this local PES, as it contains only one easily identifiable Fermi pair. But in the case of the global PES, where the latter are mixed with tunneling doublets, it becomes complicated to find the different conjugate pairs formed within the spectrum of the system. The structure of vibrational states that have been less affected by the Fermi resonance informs us about the primary nature of the levels to which they belong, and this makes it easier to identify the different tunneling doublets and Fermi pairs in the spectrum.

Finally, by making the comparison with the vibrational states of the global PES given in Table 1 of the main text, we identify the following correspondence between the local and global PES levels (the first number refers to the local, the second to the global PES): 1 ↔ 1, 2 ↔ 3, 3 ↔ 5, and 4 ↔ 8. Energies are not exactly comparable due to the approximate character of the local PES, which departs from the global PES at higher energies, and tunneling, which is absent in the local PES model. It is therefore difficult to extract the exact contribution of tunneling to the energy gap between the two levels based on these elements only.